# Thermal body force definition

Hi,

I am trying to use FreeFem to solve thermal expansion, but my results make no sense.

With this problem definition

``````            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?

Are you sure that the CG converge?

In fact, I am not. Is there a better solver?

You can remove solver=CG to see if UMFPACK solver works. (this the default ones).

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

The model is coupled so I think you make a erreur in your modélisation.

Hi,

I am not so sure my freefem code is correct. Here is my code:

``````load "msh3"

// 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)(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.

I have try you code,

just two remark , in fact the temperature affine T = (Tinit-Thot)/2 + z/(Tinit-Thot);
so
``````real cc=0.2/u1[].linty;