How to compute body forces

General Discussion
1

Hook law in weak formulation is

image
image859×201 10 KB

where f are body forces (gravity for example). If there are no body forces, then the last integral simply equals zero.

So I made an example, where one side of the cube as displacement prescribed.

load "medit"
include "cube.idp"
int[int] numberOfNodes = [4, 4, 4]; // Nodes per side.
real [int, int] boundaryRange = [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]; // Side lengths.
int [int, int] labels = [[1, 2], [3, 4], [5, 6]]; // Labels.
mesh3 Th = Cube(numberOfNodes, boundaryRange, labels);

real E = 2230e6;
real sigma = 0.4;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));

fespace Vh(Th,P2);
Vh u1,u2,u3, v1,v2,v3;

real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3)  [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

solve Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
//   - int3d(Th) (gravity*v3)
  + on(1,u1=0,u2=0,u3=0)
  + on(2,u1=0.025)
  ;

real dmax = u3[].max;
cout << " max deplacement = " << dmax << endl;
real coef= 0.1/dmax;
int[int] ref2=[1,0,2,0];
mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2);
Thm=change(Thm,label=ref2);
plot(Th,Thm, wait=1,cmm="coef  amplification = "+coef );

The results look fine:

image
image1339×983 23.1 KB

Problem

According to O.C. Zienkiewicz, … J.Z. Zhu, in The Finite Element Method: its Basis and Fundamentals (Seventh Edition), 2013

image
image1176×240 22.6 KB

meaning the sum of stress tensor divergence div(sigma) and body forces b should be zero. Now in my case, there are NO body forces, because displacement is prescribed, meaning the divergence of the stress tensor should equal zero.

But it does not.

Here is how I computed the stress divergence:

// UPDATE
macro getTensorDivX(sxx, syy, szz, syz, sxz, sxy) (dx(sxx) + dy(sxy) + dz(sxz)) // EOM
macro getTensorDivY(sxx, syy, szz, syz, sxz, sxy) (dx(sxy) + dy(syy) + dz(syz)) // EOM
macro getTensorDivZ(sxx, syy, szz, syz, sxz, sxy) (dx(sxz) + dy(syz) + dz(szz)) // EOM
macro sigma(u1, u2, u3) [(lambda + 2 * mu) * dx(u1) + lambda * dy(u2) + lambda * dz(u3),
                          lambda * dx(u1) + (lambda + 2 * mu) * dy(u2) + lambda * dz(u3),
                          lambda + dx(u1) + lambda * dy(u2) + (lambda + 2 * mu) * dz(u3),
                          2 * mu * (dz(u2)+dy(u3)),
                          2 * mu * (dz(u1)+dx(u3)),
                          2 * mu * (dy(u1)+dx(u2))] // EOM

fespace Nh(Th, P1);
Nh sxx, syy, szz, syz, sxz, sxy;
Nh rx, ry, rz, r;

// Compute stresses.
sxx = sigma(u1, u2, u3)[0];
syy = sigma(u1, u2, u3)[1];
szz = sigma(u1, u2, u3)[2];
syz = sigma(u1, u2, u3)[3];
sxz = sigma(u1, u2, u3)[4];
sxy = sigma(u1, u2, u3)[5];

// Compute residuums.
rx = getTensorDivX(sxx, syy, szz, syz, sxz, sxy);
ry = getTensorDivY(sxx, syy, szz, syz, sxz, sxy);
rz = getTensorDivZ(sxx, syy, szz, syz, sxz, sxy);

for (int i = 0; i < r.n; i++) {
  r[][i] = dist(rx[][i], ry[][i], rz[][i]);
}

cout << "Average body force (stress divergence) = " << r[].sum / r.n << endl;

The average residuum value is of order 10^7, which is not even close to zero. Does anybody know how to compute body forces in freefem?

2

Your first equation is Navier equation, Hooke’s law is the linear constitutive relation between stress and strain tensors. The body forces are prescribed, normally you don’t ‘compute’ them. Now I understand you want to estimate to what degree the equilibrium equation (2) is satisfied by the FEM solution. Note that you use a weak displacement formulation, so you can’t expect it to hold exactly.
That the average residual value you compute is of the order of 10^7 does not mean that there is necessarily a problem. You are using E=2.23 GPa=223 10^7 Pa to set the Lamé constants lambda and mu, so you should compare stresses to a similar scale.
Maybe you can observe the values of the different terms of the stress divergence to check if they actually are larger than 10^7 and cancel each other when you sum them up?

Now, you can also simplify your experiment by imposing displacement along only one axis, similar to a compression-expansion experiment:

solve Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
//   - int3d(Th) (gravity*v3)
  + on(1,u1=0)
  + on(2,u1=0.025)
  ;

The result has very small average residual value. Clearly, this way you avoid the fast deformation around the clamped side resulting from on(1,u1=0,u2=0,u3=0).

3

Yes, normally you do not compute body forces. But this is only a small part of an iterative process, where computation of body forces is the main thing.

I am not sure I understand what you mean by this. Could you be more specific?

Yes, sure, but that’s a different problem. I may have better idea but not enough freefem knowledge. I had an idea to rewrite the body force term from the Navier equation using some calculus and divergence theorem, see MFEM - Finite Element Discretization Library chapter Weak divergence. This mathematically seems correct, yet I am lacking some freefem knowledge to implement that.