How to set the thermal expansion coefficient in a certain direction (z) to 0?

dolfinx
1

It is well known that the definitions of the three-dimensional thermoelastic constitutive equations and variational formulations are as follows. Based on this, can the thermal expansion coefficient in a certain direction be set to 0? For example, let alpha = [a, a, 0] transpose, which indicates that the thermal expansion coefficients in the xy-plane are a, while there is no expansion in the z-direction due to temperature changes. If this is achievable, how should the current code be modified? I look forward to your help. Thank you.

# Parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
Delta_T = 1000

def eps(v):
    return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
    return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
#Function space
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
# Variational formula
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, Delta_T), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
2

Hello, yes the thermal expansion strain is in general:

\boldsymbol{\varepsilon_{th}} = \boldsymbol{\alpha}\Delta T

where \boldsymbol{\alpha} is a 2nd rank tensor

You can modify your stress computation as follows:

alpha = as_matrix([[1e-5, 0, 0], [0, 1e-5, 0], [0, 0, 0]])
def sigma(v, dT):
    eps_th = alpha*dT
    eps_el = eps(v) - eps_th
    return lmbda*tr(eps_el)*Identity(3)+2.0*mu*eps_el
3

Thank you very much, your method is very useful! It successfully eliminates the thermal expansion strain in the z direction. I apologize for wanting to ask further: Is it possible to eliminate all strains in the z direction (thermal strain + elastic strain)?

4

Why not just do everything in plane strain 2D then ?

5

I would like to see the result of a 3D model considering only the strains in the x and y directions (thermal strain + elastic strain), even though it is equivalent to 2D plane strain. In other words, I want to visualize a 3D cylinder with strains only in the x and y directions. If this can be achieved by modifying the constitutive equations, I hope you can assist me; otherwise, I might consider giving up on this idea. Thank you.

6

Then just define a 2D function space on a 3D mesh and change the total strain
to

def eps(v):
    grad = as_matrix([[u[0].dx(0), u[0].dx(1), 0], 
                      [u[1].dx(0), u[1].dx(1), 0], 
                      [0, 0, 0]])
    return sym(grad)
7

Hi, sorry for replying so late. I tried your suggestions in the code, but it seems they are not easy to implement:

# The thermal elastic calculation eliminates the strain in the z-direction
from mpi4py import MPI
from dolfinx import io
from dolfinx.io import gmshio
from petsc4py import PETSc
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, as_matrix, sym)
from dolfinx.fem import (functionspace, Function, form)
from dolfinx.fem.petsc import (LinearProblem, assemble_matrix, assemble_vector)
# Mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
# Parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
# Function space
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim - 1, )))   #
U = Function(VU)
def eps(v):   #
    grad = as_matrix([[U[0].dx(0), U[0].dx(1), 0],
                      [U[1].dx(0), U[1].dx(1), 0],
                      [0, 0, 0]])
    return sym(grad)
def sigma(v, dT):
    return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
# Variational formula
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
# Assemble matrix
A = assemble_matrix(form(aM))
A.assemble()
b = assemble_vector(form(LM))
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Set solver options
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_mat_solver_type"] = "mumps"
#Create PETSc Krylov solver and turn convergence monitoring on
solver = PETSc.KSP().create(mesh.comm)
solver.setFromOptions()
solver.setOperators(A)
# Solve
solver.solve(b, U.vector)

And I encountered the following error:

Traceback (most recent call last):
  File "/home/dyfluid/Desktop/data_test2/test1.py", line 34, in <module>
    A = assemble_matrix(form(aM))
        ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/functools.py", line 909, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/petsc.py", line 424, in assemble_matrix
    A = _cpp.fem.petsc.create_matrix(a._cpp_object)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.
8

You must use v in eps not U

9

Indeed, as you said, thank you very much for your patient explanation!