Hello everyone, when using the following code to calculate the thermal expansion displacement in a 2D plane, the theoretical value for the displacement of the right side edge’s position, denoted as x, in the rectangular region should be 0.72. However, the simulation result is 0.94. On the other hand, when I used the same code to calculate a 3D model, the simulation result matched the theoretical value perfectly. I am unsure of the problem and would appreciate your guidance. Thank you! Below is my code and the 2D model is a rectangle, while the 3D model is a cuboid：
（The boundary condition for the 2D case is that the displacement of the left side boundary is fully fixed. For the 3D case, the boundary condition is that the x and y displacements at z=0 are both zero.）
codes
from dolfin import *
import matplotlib.pyplot as plt
mesh = Mesh(‘test.xml’)
boundary = MeshFunction(“size_t”, mesh, “test_facet_region.xml”)
subdomains = MeshFunction(“size_t”, mesh, “test_physical_region.xml”)
ds = Measure(‘ds’, domain=mesh, subdomain_data=boundary) # dispensable
dx = Measure(‘dx’, domain=mesh, subdomain_data=subdomains)
E = Constant(2e5) # MPa
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
alpha = Constant(1.2e-5) # linear
Delta_T = 600
f = Constant((0, 0))
def eps(v):
return sym(grad(v))
def sigma(v, dT):
return (lmbdatr(eps(v))- alpha(3lmbda+2mu)dT)Identity(2) + 2.0mueps(v)
Vu = VectorFunctionSpace(mesh, ‘CG’, 2)
du = TrialFunction(Vu)
u_ = TestFunction(Vu)
Wint = inner(sigma(du, Delta_T), eps(u_))*dx
aM = lhs(Wint)
LM = rhs(Wint) + inner(f, u_)*dx
A = assemble(aM)
b = assemble(LM)
bcu = [DirichletBC(Vu, Constant((0, 0)), boundary, 7)]
u = Function(Vu, name=“Displacement”)
solve(aM == LM, u, bcu)
vtkfile = File(‘displacement/u.pvd’)
vtkfile << u
2D-model:
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {100, 0, 0, 1.0};
//+
Point(3) = {100, 5, 0, 1.0};
//+
Point(4) = {0, 5, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
Transfinite Line {1, 3} = 11;
Transfinite Line {2, 4} = 6;
//+
Physical Curve(“bottom”, 5) = {1};
//+
Physical Curve(“top”, 6) = {3};
//+
Physical Curve(“left”, 7) = {4};
//+
Physical Curve(“right”, 8) = {2};
//+
Physical Surface(“2D”, 9) = {1};
3D-model:
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {5, 0, 0, 1.0};
//+
Point(3) = {5, 5, 0, 1.0};
//+
Point(4) = {0, 5, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
Transfinite Line {1, 3} = 11;
Transfinite Line {2, 4} = 6;
Extrude {0, 0, 100} { Surface{1}; Layers {10}; }
//+
Physical Surface(“bottom”, 5) = {1};
//+
Physical Surface(“top”, 27) = {26};
//+
Physical Surface(“sidese”, 28) = {17, 21, 25, 13};
//+
Physical Volume(“body”, 29) = {1};