# Thermal body force definition
**URL:** https://community.freefem.org/t/thermal-body-force-definition/1470
**Category:** General Discussion
**Created:** [February 4, 2022, 7:40am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470 "2022-02-04T07:40:10Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [February 4, 2022, 7:40am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/1 "2022-02-04T07:40:10Z")
Hi,
I am trying to use FreeFem to solve thermal expansion, but my results make no sense.
With this problem definition
```auto
solve ThermoPlasticity([tdux, tduy, tduz], [vx, vy, vz], solver = CG) =
int3d(Th)(
lam * div(tdux, tduy, tduz) * div(vx, vy, vz)
+ 2. * G * (epsilon(tdux, tduy, tduz)' * epsilon(vx, vy, vz))
)
- int3d(Th)(beta * grad(T)' * [vx, vy, vz])
+ on(bottomLabel, tdux = 0.0, tduy = 0.0, tduz = 0.0)
;
```
Here `T` is the temperature field and the results are completely different from Abaqus, for example.
Does the problem definition look OK to you or do you see an obvious mistake?
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [February 4, 2022, 9:12am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/2 "2022-02-04T09:12:05Z")
Are you sure that the CG converge?
---
### Author:  [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [February 4, 2022, 9:14am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/3 "2022-02-04T09:14:09Z")
In fact, I am not. Is there a better solver?
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [February 4, 2022, 9:17am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/4 "2022-02-04T09:17:58Z")
You can remove solver=CG to see if UMFPACK solver works. (this the default ones).
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [February 4, 2022, 9:24am UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/5 "2022-02-04T09:24:17Z")
Without the full model is hard to say some thing,
I think in Abaqus documentation you have the full model,
Where the temperature field come
I found so information on the web
> **[Simulation of the thermomechanical and metallurgical behavior of steels by...](https://www.sciencedirect.com/science/article/abs/pii/S0927025612005824)**
>
> Simulation is a very helpful and valuable work tool in the field of heat treatment of steels. It allows behavior laws and algorithms to be learned and…
The model is coupled so I think you make a erreur in your modélisation.
---
### Author:  [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [February 4, 2022, 1:03pm UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/6 "2022-02-04T13:03:30Z")
Hi,
I am not so sure my freefem code is correct. Here is my code:
```auto
load "msh3"
load "medit"
load "iovtk"
// geometry
int N = 20; // number of nodes in a single dimension
real L = 1; // cube side length
// mesh
mesh3 Th = cube(N, N, N, [x * L - L * 0.5, y * L - L * 0.5, z * L - L * 0.5]);
// print mesh & geometry statistics
cout << "Volume = " << Th.measure << ", border area = " << Th.bordermeasure << endl;
cout << "Number of vertices = " << Th.nv << endl;
// The finite element space defined over mesh
fespace Vh(Th, P1);
fespace Nh(Th, [P1, P1, P1]);
// Variables
real sqrt2 = sqrt(2.);
real Tinit = 20.; // initial temperature
real Thot = 50.;
real totalTime = 60.; // total simulation time
real dt = 1.0; // time step
// matrial properties
real lambda = 1.; // thermal conducitvity
real cp = 1.; // specific heat
real rho = 1.; // material density
real E = 2230; // elasticity (Young) modulus
real poisson = 0.4; // poisson ratio
// Lame parameters
real mu = E / 2. / (1. + poisson);
real lam = E * poisson / (1. - 2. * poisson) / (1. + poisson);
real G = E / 2. / (1. + poisson);
real beta = 0.00014;
// FEM variables
Vh Told, T = Tinit, v;
Nh [u1,u2,u3], [v1,v2,v3];
// Macros
macro stress [Sxx, Syy, Szz, Sxy, Sxz, Syz] // stress tensor
macro Grad(u) [dx(u), dy(u), dz(u)] //EOM
macro Div(u1, u2, u3) (dx(u1) + dy(u2) + dz(u3)) // EOM
macro Epsilon(u1, u2, u3) [dx(u1), dy(u2), dz(u3), (dy(u1) + dx(u2)) / sqrt2, (dz(u1) + dx(u3)) / sqrt2, (dz(u2) + dy(u3)) / sqrt2] // EOM
macro e(u1, u2, u3) [dx(u1), dy(u2), dz(u3), dy(u1) + dx(u2), dz(u1) + dx(u3), dz(u2) + dy(u3)] // EOM
macro matE [[(lam + 2. * mu), lam, lam, 0., 0., 0.],
[lam, (lam + 2. * mu), lam, 0., 0., 0.],
[lam, lam, (lam + 2. * mu), 0., 0., 0.],
[0., 0., 0., mu, 0., 0.],
[0., 0., 0., 0., mu, 0.],
[0., 0., 0., 0., 0., mu]] // EOM plane strain stress-strain matrix
// Heat problem
problem Heat(T,v) = int3d(Th)(T * v / dt)
+int3d(Th)((lambda / rho /cp) * Grad(T)' * Grad(v))
-int3d(Th)(Told * v / dt)
+on(6, T = Thot)
+on(5, T = Tinit);
;
// Navier problem
problem Navier([u1, u2, u3], [v1, v2, v3]) = int3d(Th)(lam * Div(u1, u2, u3) * Div(v1, v2, v3))
+int3d(Th)(2. * mu * (Epsilon(u1, u2, u3)' * Epsilon(v1, v2, v3)))
-int3d(Th)(beta * Grad(T)' * [v1, v2, v3])
+on(5, u1 = 0, u2 = 0, u3 = 0) // fixed bottom plate
;
// Solve Heat problem
cout << "Solving body temperature ..." << endl;
for (real t = 0.; t < totalTime; t += dt) {
cout << "Solving for t = " << t << " seconds." << endl;
// Heat
Told = T;
Heat;
}
// Solve Navier
cout << "Solving Navier from given temperature field ..." << endl;
Navier;
// ************
string filename = "../results/mj.vtk";
cout << "*" << endl;
cout << "*" << endl;
cout << "Saving results to: " << filename << endl;
int[int] Order = [1, 1, 1];
string DataName = "u temp";
savevtk(filename, Th, [u1, u2, u3], T, dataname=DataName, order=Order);
```
The displacements due to expansion are e-8, which doesn’t make sense to me.
Seems to me like the gradient is not computed correctly or I am missing something else.
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [February 4, 2022, 3:52pm UTC](https://community.freefem.org/t/thermal-body-force-definition/1470/7 "2022-02-04T15:52:01Z")
I have try you code,
just two remark , in fact the temperature affine T = (Tinit-Thot)/2 + z/(Tinit-Thot);
so
grad(T) = [0,0,1/(Tinit-Thot)]
with your constant de displacement is of order 1e-7 so a add 3 lignes to see the deforme mesh
```auto
real cc=0.2/u1[].linty;
mesh3 Thm = movemesh(Th,[x+cc*u1,y+cc*u2,z+cc*u3]);
plot(Thm,T,wait=1);
```