FEM stress analysis by considering thermal strain

Following the question in this post, I am trying to do a stress analysis by including thermal strain. @user21 already showed how to do a pre-stressed analysis. Analogically, I modified the code and included a thermal strain term. Let's consider a very simple case that the thermal strain [Epsilon]t is uniform in the body.

\[Epsilon]t = 0.0001.
pst = {Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)))}} .
Inactive[Plus][
Inactive[Grad][u[x, y], {x, y}], {\[Epsilon]t, 0}], {x, y}] +
Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-((Y*(1 - \[Nu]))/(2*(1 \- \[Nu]^2))), 0}} .
Inactive[Plus][
Inactive[Grad][v[x, y], {x, y}], {0, \[Epsilon]t}], {x, y}],
Inactive[
Div][{{0, -((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)))}, \
{-((Y*\[Nu])/(1 - \[Nu]^2)), 0}} .
Inactive[Plus][
Inactive[Grad][u[x, y], {x, y}], {\[Epsilon]t, 0}], {x, y}] +
Inactive[
Div][{{-((Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))),
0}, {0, -(Y/(1 - \[Nu]^2))}} .
Inactive[Plus][
Inactive[Grad][v[x, y], {x, y}], {0, \[Epsilon]t}], {x, y}]};


Let's consider the same geometry, but with the beam bottom and left surfaces restrained and no external forces applied:

Needs["NDSolveFEM"];
L = 1;
h = 0.125;
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^3, \[Nu] -> 33/100};

{uif, vif} =
NDSolveValue[{pst == {0, 0},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, y == -h]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];

dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


I think the problem is straightforward, but I cannot make this code run. Can anybody help?

• Where you got the system of equations? Usually the thermal stress is defined as a volumetric force proportional to $\nabla T$. Commented May 3, 2020 at 9:08
• @AlexTrounev, thank you! Sorry, I am new to FEM stress analysis. Can you please further explain why we should use et=T0 (x+y) ? I actually need to do an analysis with nonuniform thermal strain in the body, so your explain would help me solve my final problem Commented May 3, 2020 at 18:24
• We can use any function et[x,y]. I used linear function as an example. Practically we should include equation for temperature in the system and solve it together with elasticity equations. Commented May 3, 2020 at 18:52
• @AlexTrounev Thanks. How about modeling a uniform thermal strain for verification purpose? A uniform thermal strain in the body can be et=T0, right? Then this becomes my original code. However, the displacements are 0 if I use et=T0, which does not make sense (there should be displacements). Commented May 3, 2020 at 19:11
• We have et under Div[] operator. So with et=T0 we have zero effect. Linear function it is the simplest case of thermal effect. Also we can put this on a border as a DirichletCondition[] or NeumannValue[]. But you put both zero. So what effect you expect? Commented May 3, 2020 at 19:29

We have replaced $$\epsilon t\rightarrow et, \nu\rightarrow nu$$ and put et=T0 (x+y). Also we used Activate in NDSolve. Code:

et = T0 (x + y);
pst = {Inactive[
Div][{{-(Y/(1 - nu^2)),
0}, {0, -((Y*(1 - nu))/(2*(1 - nu^2)))}}.Inactive[Plus][
Inactive[Grad][u[x, y], {x, y}], {et, 0}], {x, y}] +
Inactive[
Div][{{0, -((Y*nu)/(1 -
nu^2))}, {-((Y*(1 - nu))/(2*(1 - nu^2))), 0}}.Inactive[
Plus][Inactive[Grad][v[x, y], {x, y}], {0, et}], {x, y}],
Inactive[
Div][{{0, -((Y*(1 - nu))/(2*(1 - nu^2)))}, {-((Y*nu)/(1 -
nu^2)), 0}}.Inactive[Plus][
Inactive[Grad][u[x, y], {x, y}], {et, 0}], {x, y}] +
Inactive[
Div][{{-((Y*(1 - nu))/(2*(1 - nu^2))),
0}, {0, -(Y/(1 - nu^2))}}.Inactive[Plus][
Inactive[Grad][v[x, y], {x, y}], {0, et}], {x, y}]};

Needs["NDSolveFEM"];
L = 1;
h = 0.125;
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^3, nu -> 33/100, T0 -> 1/10};

{uif, vif} =
NDSolveValue[{Activate[pst == {0, 0} /. materialParameters],
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, y == -h]}, {u,
v}, {x, y} \[Element] mesh];


Visualisation

dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

{DensityPlot[uif[x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Bottom], PlotLabel -> "u",
AspectRatio -> Automatic],
DensityPlot[vif[x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Bottom], PlotLabel -> "v",
AspectRatio -> Automatic]}


There is a PDE model collection in the documentation. Among other examples there is a Thermal Load on a Beam application example. You can base your work off of that.

• thank you! Very useful! Commented May 4, 2020 at 6:48

I have also been struggling to understand the implementation of thermal strains, and here is an alternative attempt that may also help others. These strains are essentially the same as "eigenstrains" or "stress-free" strains. The idea being that an element without constraint will have a certain strain due to swelling, thermal expansion/contraction, phase transformation etc. It is still unclear to me whether it is possible to directly compute the role of thermal strains (like @Alex Trounev) has done, or whether one needs to use a time dependent model. As @XIJUN SHI suggested in the comments simulations of a uniform thermal strain (i.e. a mono material with a constant temperature) should be possible, however if one does this by for example inserting et = T0 ; in @Alex Trounev's code, as this is a constant it disappears and the final strain/deformation is zero.

@Alex Trounev suggested coupling this with a model for temperature changes. We can do this by using some of the code from the Thermal Load on a Beam example mentioned by @user21. The idea is to have a temperature load that is time dependent, and we include this as a body force and solve as a function of time.

To keep the code consistent we use the same as above:

Needs["NDSolveFEM"];
L = 1;
h = 0.125;
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^3, nu -> 33/100, \[Rho] -> 0.001, \[Alpha] -> 0.001};


Rho is a density term , alpha is our thermal expansion coefficient.

We define a stress operator the same as done in the example

pst={\[Rho] D[u[t,x,y],{t,2}]+
Inactive[Div][{{-(Y/(1 - nu^2)), 0}, {0, -((Y*(1 - nu))/(2*(1 - nu^2)))}} . Inactive[Grad][u[t, x, y], {x, y}], {x, y}] +
Inactive[Div][{{0, -((Y*nu)/(1 - nu^2))}, {-((Y*(1 - nu))/(2*(1 - nu^2))), 0}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}],
\[Rho] D[v[t,x,y],{t,2}]+Inactive[Div][{{0, -((Y*(1 - nu))/(2*(1 - nu^2)))}, {-((Y*nu)/(1 - nu^2)), 0}} . Inactive[Grad][u[t, x, y], {x, y}], {x,y}] + Inactive[Div][{{-((Y*(1 - nu))/(2*(1 - nu^2))), 0}, {0, -(Y/(1 - nu^2))}} .Inactive[Grad][v[t, x, y], {x, y}], {x, y}]} /.materialParameters;


We define a temperature function that goes from T = 0 to T = 1 and is uniform over the mesh.

Tfun[t, x, y] = t;


We now introduce the thermal load as a body force.

fx = Inactive[Div][{-(Y \[Alpha]/(1 - nu)) Tfun[t, x, y], 0}, {x, y}] /.materialParameters;
fy = Inactive[Div][{0, -(Y \[Alpha]/(1 - nu)) Tfun[t, x, y]}, {x, y}] /. materialParameters;


Add in boundary conditions and initial conditions:

gd = {DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, x == 0 && y == 0], DirichletCondition[{u[t, x, y] == 0.}, x == 0]};

ic = {u[0, x, y] == 0, v[0, x, y] == 0, (D[u[t, x, y], t] == 0) /. t -> 0, (D[v[t, x, y], t] == 0) /. t -> 0};

pde = {pst == {fx, fy}, gd, ic};


We can now solve:

{uif, vif} = NDSolveValue[pde, {u, v}, {t, 0, 1}, {x, y} \[Element] mesh];


and plot now the strain in the x-direction

Plot[D[uif[1, x, 0.0], x] // Evaluate, {x, 0, 1},PlotRange -> {{0, 1.1}, {0, 0.002}}]


The strain is now 0.001, (same as the thermal strain). A similar graph is achieved for the strain in the y-direction.