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*(2*mu + 3*lambda)) 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):
return ufl.sym(ufl.grad(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()
```