Assume I am solving the thermal expansion of a cylinder at 100°C (with the axial direction as the z-axis and the height being 10mm), where the mechanical boundary conditions restrict the z-direction displacement at the top and bottom faces to be 0. Clearly, this is a problem that is symmetric about z = 5mm. If I only model half of the cylinder (with a height of 5mm) and apply a displacement restriction of 0 in the z-direction at one of the bottom faces, how can I apply symmetric boundary conditions on the symmetry plane to obtain the same results as the original model?
Here is part of my code:
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import gmshio
from dolfinx import default_scalar_type
from ufl import (TestFunction, TrialFunction, dx, grad, tr, Identity, inner, lhs, rhs)
from dolfinx.fem import (functionspace, Function, dirichletbc, locate_dofs_topological, form, Constant)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector, apply_lifting, set_bc)
import nullspace1
# mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("sym.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
Delta_T = 100
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)
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim,)))
#The z-direction displacement at one bottom face is 0.
U_D = Constant(mesh, default_scalar_type(0))
fdim = mesh.topology.dim - 1
facets_bottom = facet_markers.find(4)
dofs_bottom_z = locate_dofs_topological(VU.sub(2), fdim, facets_bottom)
bcU1 = dirichletbc(U_D, dofs_bottom_z, VU.sub(2))
bcU = [bcU1]
#symmetry boundary condition:
facets_sym = facet_markers.find(5) #The symmetry plane
"...Pending"
#Variational problem
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, Delta_T), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
A = assemble_matrix(form(aM), bcs = bcU)
A.assemble()
b = assemble_vector(form(LM))
apply_lifting(b, [form(aM)], bcs=[bcU])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcU)
#Create Krylov solver
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
#Set the null space
null_space = nullspace1.null_space(VU)
A.setNearNullSpace(null_space)
#solve
U = Function(VU)
solver.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
solver.solve(b, U.vector)
and below is the mesh for the half-model：
N = 31;
Circle(1) = {0, 0, 0, 7, 0, 2*Pi};
Transfinite Curve {1} = N Using Progression 1;
Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 5} {Surface{1};Layers {5};}
Physical Surface("bottom", 4) = {1};
Physical Surface("sym", 5) = {3};
Physical Surface("cylinder", 6) = {2};
Physical Volume("v", 7) = {1};