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):
    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()
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