Hello everyone, I am solving an unconstrained thermal expansion problem in a 2D circular domain(The temperature is 1000K). I am using the method from this link to eliminate translations and rotations in the model:https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/3?u=french_fries
but the visualization results are not correct. Displacements caused by translations and rotations still exist. How should I proceed?
Here are my code:
from mpi4py import MPI
from dolfinx.fem import form
from dolfinx.fem.petsc import create_vector,assemble_matrix,assemble_vector
from petsc4py import PETSc
from dolfinx.io import gmshio
from dolfinx.fem import (functionspace, Function)
import materials
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
mesh, cell_markers, facet_markers = gmshio.read_from_msh("f1.msh", MPI.COMM_WORLD, gdim=2)
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
uh = Function(VU)
def eps(v):
return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
return (materials.lmbda_f * tr(eps(v)) - materials.alpha_f * (3 * materials.lmbda_f + 2 * materials.mu_f) * dT) * Identity(2) + 2.0 * materials.mu_f * eps(v)
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
# Assemble system
A = assemble_matrix(form(aM))
A.assemble()
b=create_vector(form(LM))
with b.localForm() as b_loc:
b_loc.set(0)
assemble_vector(b,form(LM))
# Create Krylov solver
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
A.setNullSpace(nullspace)
nullspace.remove(b)
solver.solve(b,uh.vector)
and geo files:
Point(1) = {0, 0, 0, 1.0};
Point(2) = {0, 7, 0, 1.0};
Point(3) = {0, -7, 0, 1.0};
Circle(1) = {2, 1, 3};
Circle(2) = {3, 1, 2};
Line(3) = {2, 3};
Curve Loop(1) = {2, 3};
Plane Surface(1) = {1};
Curve Loop(2) = {1, -3};
Plane Surface(2) = {2};
Physical Curve("f_s", 4) = {2, 1};
Physical Surface("f", 5) = {1, 2};
Mesh 2;
RecombineMesh;
Mesh.SubdivisionAlgorithm = 1;
The visualization results are shown in the figure:
