FEM stress analysis by considering thermal strain

Question

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["NDSolve`FEM`"];
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?

Answer 1

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["NDSolve`FEM`"];
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]}

Answer 2

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.

Answer 3

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["NDSolve`FEM`"];
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.

I hope this is helpful.