# Body forces in static equilibrium (FEM issue)

• skrat
In summary, the conversation discusses the issue of understanding body forces in FEM simulations using freefem. A simple unit cube is used to demonstrate the problem and the script used for the simulations is provided. The results look fine, but when trying to compute the body forces, the average norm is of order 10^7, which is unexpected. The expert suggests that the medium is not in static equilibrium and an external force needs to be applied to resist the medium's desire to return to its equilibrium state. The expert also points out an error in the definition of the stress tensor and suggests using the weak formulation of the hookean law to compute the body forces.
skrat
TL;DR Summary
According to Zienkiewicz, the body forces in static equilibrium should equal zero, but in my simulations they don't.
Hi,

I have an issue with understanding the body forces in the context of FEM simulations. I am using freefem to do the simulations.

So I took a simple unit cube, fixed the displacements at side ##x=0## and prescribed a displacement ##u_x = 0.025## at ##x=1##. If required, the entire script is given below.

Cube example:
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 );

Now the results look completely fine to me:

However, the problem.

The problems start when I try to compute the body forces. Now, according to O.C. Zienkiewicz, ... J.Z. Zhu, in The Finite Element Method: its Basis and Fundamentals (Seventh Edition), 2013 (snapshot of a chapter given here) for any element in the static equilibrium we can write
$$\nabla \cdot \sigma + b = 0$$
where ##\sigma## is the stress tensor and ##b## is a vector of body forces.

Now in my case, there are no external body forces (e.g. gravity), so ##b=0## which in turn means ##\nabla \cdot \sigma = 0##.

So I computed the divergence of the stress tensor, where the divergence of a tensor ##T## is defined as

$$(\nabla \cdot T)_i = \frac{\partial T_{ij}}{\partial x_j}$$

assuming the Einstein summation rule.

Here is my script:

Stress divergence computation.:
// 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;

But I was extremely surprised to find out that the average norm of body forces is of order ##10^7##, which is not even close to 0.

Question

I have two questions.

1. I am not expecting that any of you is familiar with the FreeFem++. If you are, even better, you can help me understand what I am doing wrong in my script.
2. Am I making wrong assumptions? Should I not be surprised that the divergence of the stress tensor is not zero? Is this not even a problem but rather an expected result? If yes, could you help me understand why?

I think the conclusion is that your medium is not in static equilibrium.

You are essentially pulling the medium along the $x$-axis in order to achieve the displacement. If you want the result to be in static equilibrium, you need to apply an external force to resist the medium's desire to return to its equilibrium ($u \equiv 0$) state.

I'm also not convinced that your definition of $\sigma$ is consistent with your definition of $\epsilon$. Generally in a linear theory you would have $$\epsilon_{ij} = \frac12 \left( \frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right)$$ and $$\sigma_{ij} = 2\mu \epsilon_{ij} + \lambda \epsilon_{kk} \delta_{ij}.$$

Yes, you are correct, my definition of ##\sigma## is a wrong by factor 2 on all off-diagonal terms.

What if, excuse me for brainstorming, I wanted to take the last term from the weak formulation of the hookean law
$$\int _\Omega \nabla\cdot \sigma v \mathrm{d} \Omega$$
and decided to play with a little. For example use the calculus identity
$$\nabla \cdot (\sigma v) = (\nabla \cdot \sigma)v + \sigma \nabla v.$$
This would allow me to rewrite the term to
$$\int _\Omega \nabla\cdot \sigma v \mathrm{d} \Omega = \int _\Omega \sigma \nabla v \mathrm{d} \Omega + \int _\Omega \nabla\cdot (\sigma v) \mathrm{d} \Omega$$
and using the Divergence theorem the last integral can be over ##\partial \Omega## finally resulting in
$$\int _\Omega \nabla\cdot \sigma v \mathrm{d} \Omega = \int _\Omega \sigma \nabla v \mathrm{d} \Omega+ \int _{\partial \Omega} \hat{n} \cdot \sigma v \mathrm{d} \partial \Omega.$$

I may be shooting blanks with this but the second interval could be zero as ##v## can always be set to 0 on the boundary ##\partial \Omega##. If that is true than computing the first integral is all that remains.

Does that make any sense? Would this help at all? If nothing else it should improve the numerical stability where the cube is fixed as computation of derivatives of stress tensor components there is no longer required.

## 1. What is static equilibrium in relation to body forces?

Static equilibrium refers to a state in which all forces acting on a body are balanced, resulting in no movement or acceleration. In the context of FEM, it is important to consider the effects of body forces, such as gravity or electromagnetic forces, on the equilibrium of a structure.

## 2. How do body forces affect the analysis of a structure using FEM?

Body forces can significantly impact the behavior and stability of a structure, and therefore must be taken into account in FEM analysis. These forces can induce stresses and deformations in a structure, which can affect its overall performance and safety.

## 3. What types of body forces are commonly encountered in FEM analysis?

The most common types of body forces in FEM analysis include gravitational forces, electromagnetic forces, and thermal forces. Other types of body forces may also be present, depending on the specific application.

## 4. How are body forces represented in FEM models?

In FEM, body forces are typically represented as distributed loads, which are applied over a specific area or volume of a structure. These loads are then incorporated into the FEM model to accurately simulate the effects of body forces on the structure.

## 5. How can body forces be minimized in FEM analysis?

To minimize the effects of body forces in FEM analysis, it is important to carefully consider the design and material selection of a structure. For example, using lightweight materials or incorporating structural supports can help reduce the impact of body forces on a structure.