Error in Newton-Raphson Incremental Iteration

General Discussion
1

Hi,

I am trying to use FreeFem++ to solve plasticity mechanics problems using an incremental iterative approach, with the specific code documentation being:

test.edp (15.6 KB)

tension.zip (307.4 KB) (it’s the mesh file)

, but I have discovered an issue:

When calculating stress, I use

Vh [u, v], [u1, v1]；
Eh[int] sigma(3)；
sigma[0] = E/(1-nu^2)*dx(u) + nu*E/(1-nu^2)*dy(v);
sigma[1] = nu*E/(1-nu^2)*dx(u) + E/(1-nu^2)*dy(v);
sigma[2] = 0.5*E/(1+nu)*(dx(v)+dy(u));

where u and v represent displacements, and sigma[3] represents stress.

The process of assigning values computed through dx(u) at integration points to the nodes of the sigma variable appears to introduce errors.

To verify this conjecture, I designed an extremely simple numerical example of linear elasticity

test_simple.edp (3.1 KB)

tension.zip (307.4 KB) (it’s the mesh file)

that should fully converge after the first iteration step. However, when I run the code:

solve A11([u,v],[u1,v1]) = int2d(Th)(epsilon(u,v)'*Copen(E,nu)*epsilon(u1,v1))
                                -int2d(Th)(epsilon(ubefore,vbefore)'*Copen(E,nu)*epsilon(u1,v1))               
                                -int2d(Th)([sigma[0],sigma[1],sigma[2]]'*epsilon(u1,v1))
                                +on(1,u=0,v=-displacement)+on(2,v=displacement);

, it fails to converge. But when I run:

solve A11([u,v],[u1,v1]) = int2d(Th)(epsilon(u,v)'*Copen(E,nu)*epsilon(u1,v1))
                                -int2d(Th)(epsilon(ubefore,vbefore)'*Copen(E,nu)*epsilon(u1,v1))               
                                -int2d(Th)(epsilon(ubefore,vbefore)'*Copen(E,nu)*epsilon(u1,v1))   
                                +on(1,u=0,v=-displacement)+on(2,v=displacement);

it converges completely in the first step, just as expected.

It is strange that these two codes are essentially identical. I am unsure whether this discrepancy is caused by errors introduced when interpolating values from integration points to nodes. I would appreciate any help from friends here to resolve this issue.

2

Alexander,
you are right, the difference in the results obtained by the two codes is due to the interpolation of dx(u) into the space Eh (=P1 with your definitions). If you take Eh=P0 there is no more discrepancy, the code with sigma converges also. Indeed in this case since u is P1,
dx(u) is P0 thus there is no interpolation error.

3

Thank you for your prompt reply, Professor. It’s really helpful to me.