Stress divergence after elastic deformation does not equal 0, but it should

General Discussion
1

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[0])+dy(u[1])) // divergence
 macro eM(u) [dx(u[0]), dy(u[1]), sqrt(2)*(dx(u[1]) + dy(u[0]))/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.

image
image2445×984 198 KB

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]),
		0.5 * (dy(u[0]) + dx(u[1])),
		0.5 * (dx(u[1]) + dy(u[0])), 
		dy(u[1])
	] // EOM
 macro sigma(u)
	[
		(lambda + 2. * mu) * e(u)[0] + lambda * e(u)[3],
		2. * mu * e(u)[2],
		lambda * e(u)[0] + (lambda + 2. * mu) * e(u)[3]
	] // EOM
macro stressdiv(Sxx, Syy, Sxy) [dx(Sxx) + dy(Sxy), 
                                dx(Sxy) + dy(Syy)] // EOM


sxx = sigma(u)[0];
sxy = sigma(u)[1];
syy = sigma(u)[2];

rx = stressdiv(sxx, syy, sxy)[0];
ry = stressdiv(sxx, syy, sxy)[1];

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:

image
image1262×591 39.9 KB

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.