# Trouble with sigma formulation in a simple thermoelastic problem in 3D

I’m pretty new to fenicsx so apologies if this is a dumb question. I’m just trying to simulate thermal expansion with a constant delta_T. I’ve visited old tutorials in 2D and they add: - delta_T * (alpha*(2mu + 3lambda)) to the stress expression. My case should be even easier as my delta_T is a constant, but as soon as I add the term, I get solver errors. What am I missing?

This is the code:

``````import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, io, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx import default_scalar_type
import ufl

# Define the box domain dimensions
Lx, Ly, Lz = 1.0, 1.0, 1.0  # Length in x, y, z directions
nx, ny, nz = 10, 10, 10  # Number of elements in x, y, z directions

# Create the box mesh
domain = mesh.create_box(MPI.COMM_WORLD,
[np.array([0.0, 0.0, 0.0]), np.array([Lx, Ly, Lz])],
[nx, ny, nz],
mesh.CellType.tetrahedron)

# Define the function space
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

# Define boundary condition
def clamped_boundary(x):
return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

# Define material properties
E = 10.0  # Young's modulus
nu = 0.3  # Poisson's ratio
alpha = 1e-5  # Coefficient of thermal expansion
delta_T = 100  # Temperature change

lambda_ = (E * nu / ((1 + nu) * (1 - 2 * nu)))
mu = (E / (2 * (1 + nu)))
kappa = delta_T * (alpha*(2*mu + 3*lambda_))

print(f"lambda = {lambda_}, mu = {mu}, kappa = {kappa}")

# Define the stress and strain tensors with thermal expansion
def epsilon(u):

def sigma_elastic(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Stress including thermal expansion: NOT WORKING
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u) - kappa * ufl.Identity(len(u))

# Define the variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ufl.Measure("ds", domain=domain)

# Solve the problem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
``````
1 Like

The final added does not depend on u [*]. If that is the term you are really trying to implement, then it must go in `L` on the right-hand side.

[*] I do understand that you use `u` in the expression, but that is not a true dependence on the solution u. You could have used `domain.geometry.dim` instead of `len(u)`.

For future posts, please make sure to include the error message you get.

Thanks! you are totally correct.
I did not attach the error because it was incredibly long, but next time I will. Cheers