Model displacement in unconstrained thermoelasticity problems

General
1

Hello everyone, I am solving for the unconstrained thermal expansion displacements of six 2D circular shapes (all at 1000℃), but the results are quite strange: as shown in the figure, the displacement of the bottom-right circle is very large. I cannot determine whether the model as a whole has undergone translation/rotation, or if only the bottom-right circle has translated/rotated. If it’s the latter case, it implies that the relative positions between the bottom-right circle and the others have changed, which is obviously unreasonable. How can I solve this issue? I hope to get some help.
2d0012a1a787122e6718bbd95c4ce27
Below are my code and geo file:

from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import io
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
from dolfinx.fem import (functionspace)
from dolfinx.fem.petsc import LinearProblem

mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=2)

kappa = 2.8e-3
E = 2e5
nu = 0.26
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
T = 1000       #Temperature

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(2) + 2.0 * mu * eps(v)

VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim,)))
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, T), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)

problem = LinearProblem(aM, LM, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Uf = problem.solve()

r1 = 7.06;
r2 = 7.875;
dh = 27.713;
fh = 16;
he = 1.5;

Circle(1) = {r2+he+dh*Sqrt(3)/2-fh, dh, 0, r1, 0, 2*Pi};
Circle(2) = {r2+he+fh, dh*1.5, 0, r1, 0, 2*Pi};
Circle(3) = {r2+he+fh, dh*0.5, 0, r1, 0, 2*Pi};
Circle(4) = {r2+he+dh*Sqrt(3)-fh, dh*1.5, 0, r1, 0, 2*Pi};
Circle(5) = {r2+he+dh*Sqrt(3)-fh, dh*0.5, 0, r1, 0, 2*Pi};
Circle(6) = {r2+he+dh*Sqrt(3)/2+fh, dh, 0, r1, 0, 2*Pi};

Curve Loop(1) = {1};
Plane Surface(1) = {1};
Curve Loop(2) = {2};
Plane Surface(2) = {2};
Curve Loop(3) = {3};
Plane Surface(3) = {3};
Curve Loop(4) = {4};
Plane Surface(4) = {4};
Curve Loop(5) = {5};
Plane Surface(5) = {5};
Curve Loop(6) = {6};
Plane Surface(6) = {6};

Physical Curve("1", 7) = {1};
Physical Curve("2", 8) = {2};
Physical Curve("3", 9) = {3};
Physical Curve("4", 10) = {4};
Physical Curve("5", 11) = {5};
Physical Curve("6", 12) = {6};
Physical Surface("1", 13) = {1};
Physical Surface("2", 14) = {2};
Physical Surface("3", 15) = {3};
Physical Surface("4", 16) = {4};
Physical Surface("5", 17) = {5};
Physical Surface("6", 18) = {6};
Mesh 2;
RecombineMesh;
Mesh.SubdivisionAlgorithm = 1;