Two-dimensional elasticity problem using variational definition

General Discussion
1

Mission

Consider a plane (a*b rectangular). On the left side, the edge is attached to a wall but can move freely along a vertical axis. On the right side, a pressure p is applied that stretches the plate.

A sketch:

image
image926×717 15.7 KB

Attempt at solution

So I took a look at a couple of examples and I came up with the following code relevant to the problem definition:

varf elas([u1,u2],[v1,v2]) =
  int2d(Th)(  
	    lam*div(u1,u2)*div(v1,v2)	
	    +2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) )
	      )
  + on(4, u1 = 0)
  ;

varf l([u1,u2],[v1,v2]) = int1d(Th, 2)(p * v1);

matrix A = elas(Nh, Nh); // Stiffness matrix.
real[int] b = l(0, Nh); // External forces.

real[int] U(A.n);
U = A^-1*b;

So displacement in x direction on the left side is fixed to 0 and pressure is applied to the right side. I did not fix the displacement along the y axis.

However, this does not produce correct results, as

real[int] U1(A.n/2),U2(A.n/2);
for(int i=0; i<A.n/2; i++){
    U1(i) = U(2*i);
    U2(i) = U(2*i+1);
}

cout << "max x: " << U1.max << endl; 
cout << "max y: " << U2.max << endl;

prints out

//max x: 0.00247184
//max y: 0.00852047

while my reference (Abaqus) solution results in

//max x: 2.774e-3
//max y: 2.735e-4

Question

I checked the code several times but I can not find any mistakes. I don’t know what is wrong, but I am quite sure the results are NOT ok. Please help.

I am guessing this boundary condition

varf l([u1,u2],[v1,v2]) = int1d(Th, 2)(p * v1);

is not correctly written, but I am not sure how to fix it.

Entire code

load "medit"

// Generate mesh.
real x0 = 0;
real x1 = 1;
real y0 = 0;
real y1 = 0.5;
int n = 20;
real m = 10;
mesh Th = square(n, m, [x0+(x1-x0)*x, y0+(y1-y0)*y]);
plot(Th, wait=1, cmm="Mesh");

// **********
// Variables.
// **********
real sqrt2 = sqrt(2.);
real p = 2e8;

// Material properties.
real youngModul = 72.1e9; // Elasticity (Young) modulus.
real poisson = 0.33; // Poisson ratio.

// Print constants.
if (1) {
    cout << "=========================" << endl;
    cout << "E = " << youngModul << endl;
    cout << "poissonRatio = " << poisson << endl;
    cout << "=========================" << endl;
}

// Lame parameters
real mu = youngModul / 2. / (1. + poisson);
real lam = youngModul * poisson / (1. - 2. * poisson) / (1. + poisson);
real G = youngModul / 2. / (1. + poisson);

// Print derived constants.
if (1) {
    cout << "=========================" << endl;
    cout << "mu = " << mu << endl;
    cout << "lam = " << lam << endl;
    cout << "G = " << G << endl;
    cout << "=========================" << endl;
}

// *********************
// Finite element space.
// *********************
// The finite element space defined over mesh
fespace Vh(Th, P1);
fespace Nh(Th, [P1, P1]);

// FEM variables.
Vh v;
Nh [u1,u2], [v1,v2];

// *******
// Macros.
// *******
macro grad(u) [dx(u), dy(u)] //EOM
macro div(u1, u2) (dx(u1) + dy(u2)) // EOM
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2)) / sqrt2] // EOM
macro e(u1, u2) [dx(u1), dy(u2), dz(u3), dy(u1) + dx(u2)] // EOM

macro matE [[(lam + 2. * mu), lam, 0.],
            [lam, (lam + 2. * mu), 0.],
            [0., 0., mu]] // EOM 

// ***********************
// Variational definition.
// ***********************
varf elas([u1,u2],[v1,v2]) =
  int2d(Th)(  
	    lam*div(u1,u2)*div(v1,v2)	
	    +2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) )
	      )
  + on(4, u1 = 0)
  ;

varf l([u1,u2],[v1,v2]) = int1d(Th, 2)(p * v1);

matrix A = elas(Nh, Nh); // Stiffness matrix.
real[int] b = l(0, Nh); // External forces.

real[int] U(A.n);
U = A^-1*b;

real[int] U1(A.n/2),U2(A.n/2);

for(int i=0; i<A.n/2; i++){
    U1(i) = U(2*i);
    U2(i) = U(2*i+1);
}

cout << "max x: " << U1.max << endl;
cout << "max y: " << U2.max << endl;
2

I don’t think it’s a freefem problem, but, rather, “problem setup” issue. From your variational formulation given in

it would seem to me your problem is not well posed. While I didn’t read the whole code, it seems to me that you are solving some simple elasticity problem of form

-div(\sigma)=0
+ boundary conditions

with boundary conditions of mixed type. Now, most of the boundary is Neumann with prescirbed stress (zero or p*[N.x,N.y] on border 2) and Dirichlet bc on border 4 but only for variable u1. For u2 you only have Neumann b.c. (something like “no stress” from what I can see) so u2 is defined only up to a constant. Try to plot [u1,u2] which you get from your system, probably won’t be what you expect. Also, try to add

+ on(4, u1 = 0, u2=0)

This will make your system determined. So, from what I can see, it would seem your mathematical formulation is ill-posed.

if

3

I couldn’t agree more! It is for sure not a FreeFem problem, but rather a problem of me being unable to set the boundary conditions correctly. That is exactly why I came to this forum, to figure out the right way to set up the boundary conditions.

It is quite easy and straight forward for Dirichlet BC, but so far I was unable to solve this problem where the right side us under tension. I can’t figure about a mathematical way to describe the tension in a way FreeFem would understand it.

4

Not sure what you mean by

From what I can see in your first post, your actual mathematical problem is ill-posed in y-displacement, i.e. u2 variable (unless you made a mistake in implementation which, I assume, you would have already noticed). So, if you check your mathematical formulation of the problem I believe you have (if I am following correctly)

  1. Dirichlet and Neumann b.c. for ux, and
  2. only Neumann b.c. for uy.

Since there are only Neumann b.c. for uy field, your actual mathematical problem is ill-posed (without any other assumptions). That being said, there are couple of ways to “fix” this.

  1. The easy way is, what I mentioned in the first post, specify uy on part of the boundary (impose Dirichlet b.c.).
  2. Alternatively, you may specify uy in just a single node - this practically means you are fixing the constant for your uy which is now “unique up to a constant” (something similar what you would do for pressure in Stokes system for incompressible fluids)
  3. Similar as in (2), you may fix this constant by throwing Lagrange multipliers into the mix (but then you will have to deal with saddle point problem)
  4. Ensure well posedness of your elasticity problem. What you are dealing with is basically a singular Neumann problem in linear elasticity. If I am not mistaken, this problem is not well posed in general, but it is if your external forces satisfy some compatibility conditions; net force and net torque on the body should be zero (net force = sum of body and surface forces).
  5. What you might want to look into, are the Robin boundary conditions, which would be a kind of a trade-off if you wish to keep your system as close to as what it is now.

if

5

Solution

Ok, so I figured out the solution. The biggest trick was in setting the plane stress correction:

// Lame parameters
real mu = youngModul / 2. / (1. + poisson);
real lam = youngModul * poisson / (1. - 2. * poisson) / (1. + poisson);
lam = 2 * mu * lam / (2 * mu + lam); // Plane stress correction!

Then the problem can be defined as:

problem elas([u1,u2],[v1,v2]) =
  int2d(Th)(  
	    lam*div(u1,u2)*div(v1,v2)	
	    +2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) )
	      )
  - int1d(Th, 2)(p * v1)
  + on(4, u1 = 0)
  + on(1, u2 = 0)
  ;

Which results in displacements identical to the displacements computed in two other softwares.