# Two-dimensional elasticity problem using variational definition

**URL:** https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177
**Category:** General Discussion
**Created:** [August 25, 2021, 7:08am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177 "2021-08-25T07:08:26Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![distractor](https://avatars.discourse-cdn.com/v4/letter/d/57b2e6/32.png) [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [August 25, 2021, 7:08am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177/1 "2021-08-25T07:08:27Z")

</div>

# 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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/7cd99c359dd9e29c0d542b63d98ee21afe2336b4.png)

# 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:

```auto
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

```auto
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

```auto
//max x: 0.00247184
//max y: 0.00852047

```

while my reference (Abaqus) solution results in

```auto
//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

```auto
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

```auto
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;

```

---

<div class="post-metadata">

### Author: ![fivanci](https://avatars.discourse-cdn.com/v4/letter/f/e95f7d/32.png) [@fivanci](https://community.freefem.org/u/fivanci)
#### Post date: [August 25, 2021, 8:59am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177/2 "2021-08-25T08:59:27Z")

</div>

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

> [@distractor](#):
>
> Attempt at solution

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

```auto
-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

```auto
+ 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

---

<div class="post-metadata">

### Author: ![distractor](https://avatars.discourse-cdn.com/v4/letter/d/57b2e6/32.png) [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [August 26, 2021, 9:14am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177/3 "2021-08-26T09:14:44Z")

</div>

> [@fivanci](#):
>
> I don’t think it’s a `freefem` problem, but, rather, “problem setup” issue

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.

---

<div class="post-metadata">

### Author: ![fivanci](https://avatars.discourse-cdn.com/v4/letter/f/e95f7d/32.png) [@fivanci](https://community.freefem.org/u/fivanci)
#### Post date: [August 27, 2021, 1:48am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177/4 "2021-08-27T01:48:15Z")

</div>

Not sure what you mean by

> [@distractor](#):
>
> a mathematical way to describe the tension in a way FreeFem would understand it.

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

---

<div class="post-metadata">

### Author: ![distractor](https://avatars.discourse-cdn.com/v4/letter/d/57b2e6/32.png) [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [August 27, 2021, 6:11am UTC](https://community.freefem.org/t/two-dimensional-elasticity-problem-using-variational-definition/1177/5 "2021-08-27T06:11:35Z")

</div>

# Solution

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

```auto
// 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:

```auto
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.