I have posted a previous question at this post. Thanks @AlexTrounex and @user21 for providing very useful comments. The more rigorous 2D plane stress stress operator by considering both position and time dependent temperature fields is: (a=thermal expansion coefficient, ρρ=density)
pst = {Inactive[
Div][{{-(Y/(1 - nu^2)),
0}, {0, -((Y*(1 - nu))/(2*(1 - nu^2)))}} .
Inactive[Grad][u[x, y, t], {x, y}], {x, y}] +
Inactive[
Div][{{0, -((Y*nu)/(1 -
nu^2))}, {-((Y*(1 - nu))/(2*(1 - nu^2))), 0}} .
Inactive[Grad][v[x, y, t], {x, y}], {x, y}] +
(a*Y*Inactive[D][T[x, y, t], {x}])/(1 - nu) - ρρ*
Inactive[D][u[x, y, t], {t}],
Inactive[
Div][{{0, -((Y*(1 - nu))/(2*(1 - nu^2)))}, {-((Y*nu)/(1 -
nu^2)), 0}} . Inactive[Grad][u[x, y, t], {x, y}], {x,
y}] +
Inactive[
Div][{{-((Y*(1 - nu))/(2*(1 - nu^2))),
0}, {0, -(Y/(1 - nu^2))}} .
Inactive[Grad][v[x, y, t], {x, y}], {x,
y}] + (a*Y*Inactive[D][T[x, y, t], {y}])/(1 -
nu) - ρρ*Inactive[D][v[x, y, t], {t}]};
By defining displacement and temperature fields as a function of time, I initially thought it would naturally handle a non-steady problem. However, I still have troubles getting reasonable results. Let's consider a simple example. The temperature field is just a function of time: T[t]=10t. The boundary conditions are the left and bottom surfaces of the body are fixed. The code for this problem is:
Needs["NDSolve`FEM`"];
L = 1;
h = 1;
reg = Rectangle[{0, 0}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^10, nu -> 33/100,
a -> -0.001, ρρ -> 2000};
T = 10 t;
pst = {Inactive[
Div][{{-(Y/(1 - nu^2)),
0}, {0, -((Y*(1 - nu))/(2*(1 - nu^2)))}} .
Inactive[Grad][u[x, y, t], {x, y}], {x, y}] +
Inactive[
Div][{{0, -((Y*nu)/(1 -
nu^2))}, {-((Y*(1 - nu))/(2*(1 - nu^2))), 0}} .
Inactive[Grad][v[x, y, t], {x, y}], {x, y}] +
(a*Y*Inactive[D][T, {x}])/(1 - nu) - ρρ*
Inactive[D][u[x, y, t], {t,2}],
Inactive[
Div][{{0, -((Y*(1 - nu))/(2*(1 - nu^2)))}, {-((Y*nu)/(1 -
nu^2)), 0}} . Inactive[Grad][u[x, y, t], {x, y}], {x,
y}] +
Inactive[
Div][{{-((Y*(1 - nu))/(2*(1 - nu^2))),
0}, {0, -(Y/(1 - nu^2))}} .
Inactive[Grad][v[x, y, t], {x, y}], {x,
y}] + (a*Y*Inactive[D][T, {y}])/(1 - nu) - ρρ*
Inactive[D][v[x, y, t], {t,2}]};
ic = {u[x, y, 0] == 0,
v[x, y, 0] == 0, (D[u[x, y, t], t] /. t -> 0) ==
0, (D[v[x, y, t], t] /. t -> 0) == 0};
{uif, vif} =
NDSolveValue[{Activate[pst == {0, 0} /. materialParameters], ic,
DirichletCondition[u[x, y, t] == 0, x == 0],
DirichletCondition[v[x, y, t] == 0, y == 0]}, {u, v}, {t, 0,
3}, {x, y} ∈ mesh];
{DensityPlot[uif[x, y, 3], {x, y} ∈ mesh,
ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Bottom], PlotLabel -> "u",
AspectRatio -> Automatic],
DensityPlot[vif[x, y, 3], {x, y} ∈ mesh,
ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Bottom], PlotLabel -> "v",
AspectRatio -> Automatic]}
However, by running this code, all the solved displacements are 0, which obviously does not make any sense.
Looks like that the system equation lacks a term which can account for the time change of the temperature field. Can anyone please help me?
