# TL;DR

According to the analytical solution (link) for a deflection of a cantilever beam, the divergence of the stress field should equal 0, but my FEM simulations do not result in zero.

# Problem

A simple two-dimensional cantilever beam is taken. Displacements on one side of the beam are fixed using the analytical solution, similar applies to the other side of the beam where traction boundary is assigned using the analytical solution.

# FreeFEM solution

I have used this example to get the freefem code. The only changes I did is that I wrote the full expressions for the ux and uy.

`````` // Parameters
real E = 3.0e7, nu = 0.3;
real P = -1000;
real L = 48, D = 12;
real I = D^3/12;

// Mesh
mesh Th = square(40,10,[L*x,-D/2+D*y]);

// Macros
macro u [ux,uy] // displacement
macro v [vx,vy] // test function
macro div(u) (dx(u)+dy(u)) // divergence
macro eM(u) [dx(u), dy(u), sqrt(2)*(dx(u) + dy(u))/2] // strain tensor

// Finite element space
fespace Vh(Th,[P2,P2]);
Vh u, v;

// Boundary conditions
func uxb = -P*y/(6.*E*I)*(3.*(x*x - L*L) - (2.+nu)*y*y + 6.*(1.+nu)*D*D/4.);
func uyb = P/(6.*E*I)*(3. * nu * x * y * y + x*x*x - 3.*L*L*x + 2.*L*L*L);
func ty = -P/(2*I)*(D^2/4-y^2);

// Convert E and nu to plane stress
E = E/(1-nu^2); nu = nu/(1-nu);

// Lame parameters
real mu = E/(2*(1+nu));
real lambda = E*nu/((1+nu)*(1-2*nu));

// Solve problem
solve cantileverBeam(u,v) =
int2d(Th)(lambda*div(u)*div(v) + 2*mu*(eM(u)'*eM(v)))
- int1d(Th,4)([0,-ty]'*v)
+ on(2,ux=uxb,uy=uyb);

// Plot solution
real coef = 1000;
Th = movemesh(Th,[x+ux*coef,y+uy*coef]);
plot(Th,wait=1);
``````

This results in a deflected beam that seems to be ok.

# Modifications to the code

So to compute the divergence of the stress field I have to be able to compute the stresses sxx, syy and sxy. Therefore, I added a bit of macro, including the definition of the divergence of the tensor field.

`````` fespace Nh(Th, P2);
Nh sxx, sxy, syy;
Nh rx, ry, r;

macro e(u)
[
dx(u),
0.5 * (dy(u) + dx(u)),
0.5 * (dx(u) + dy(u)),
dy(u)
] // EOM
macro sigma(u)
[
(lambda + 2. * mu) * e(u) + lambda * e(u),
2. * mu * e(u),
lambda * e(u) + (lambda + 2. * mu) * e(u)
] // EOM
macro stressdiv(Sxx, Syy, Sxy) [dx(Sxx) + dy(Sxy),
dx(Sxy) + dy(Syy)] // EOM

sxx = sigma(u);
sxy = sigma(u);
syy = sigma(u);

rx = stressdiv(sxx, syy, sxy);
ry = stressdiv(sxx, syy, sxy);

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

plot(sxx, wait=1, cmm="sxx");
plot(syy, wait=1, cmm="syy");
plot(sxy, wait=1, cmm="sxy");
plot([rx, ry], wait=1, cmm="residuum");
``````

# Results

However, the results are far from the expected. The analytical solution claims the divergence of the stress field should be zero, but the numerical simulations using the code above result in non-zero values:

I am also not sure if these numbers are “small” or “negliable”. I think the divergence equals to the density of the internal forces. If that is the case, I don’t think these numbers are small unless it is safe to compare them to the loading force P.

# Question

So, of course, this makes little to no sense to me. I would like to understand if there is a bug in my code or understand the reason why the divergence is not zero?
Since the deformation is elastic, I think the divergence equals to density of internal forces which in stationary state should most likely be zero, but again, simulations show otherwise.