Thermal body force definition

General Discussion
1

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?

2

Are you sure that the CG converge?

3

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

4

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

5

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.

6

Hi,

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

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.

7

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

real cc=0.2/u1[].linty;
mesh3 Thm = movemesh(Th,[x+cc*u1,y+cc*u2,z+cc*u3]);
plot(Thm,T,wait=1);