Non-linear elasticity

The nonlinear elasticity problem is: find the displacement (u_{1},u_{2}) minimizing J

where F2(u_{1},u_{2}) = A(E[u_{1},u_{2}],E[u_{1},u_{2}]) and A(X,Y) is bilinear symmetric positive form with respect two matrix X,Y.

where f is a given \mathcal{C}^2 function, and E[u_{1},u_{2}] = (E_{ij})_{i=1,2,\,j=1,2} is the Green-Saint Venant deformation tensor defined with:

Denote \mathbf{u}=(u_{1},u_{2}), \mathbf{v}=(v_{1},v_{2}), \mathbf{w}=(w_{1},w_{2}). So, the differential of J is

where DF2(\mathbf{u})(\mathbf{v}) = 2 \; A(DE[\mathbf{u}](\mathbf{v}),E[\mathbf{u}]) and DE is the first differential of E.

The second order differential is

\begin{eqnarray*} D^2 J(\mathbf{u})((\mathbf{v}),(\mathbf{w})) &=& \displaystyle\int DF2(\mathbf{u})(\mathbf{v}) \; DF2(\mathbf{u})(\mathbf{w}) \; f''(F2(\mathbf{u}))) \\ &+& \displaystyle\int \; D^2F2(\mathbf{u})(\mathbf{v},\mathbf{w}) \; f'(F2(\mathbf{u}))) \end{eqnarray*}

where

and D^{2}E is the second differential of E.

So all notations can be define with macro:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
macro EL(u, v) [dx(u), (dx(v)+dy(u)), dy(v)] //is [epsilon_11, 2epsilon_12, epsilon_22]

macro ENL(u, v) [
    (dx(u)*dx(u) + dx(v)*dx(v))*0.5,
    (dx(u)*dy(u) + dx(v)*dy(v)),
    (dy(u)*dy(u) + dy(v)*dy(v))*0.5
    ] //

macro dENL(u, v, uu, vv) [
    (dx(u)*dx(uu) + dx(v)*dx(vv)),
    (dx(u)*dy(uu) + dx(v)*dy(vv) + dx(uu)*dy(u) + dx(vv)*dy(v)),
    (dy(u)*dy(uu) + dy(v)*dy(vv))
    ] //

macro E(u, v) (EL(u,v) + ENL(u,v)) //is [E_11, 2E_12, E_22]
macro dE(u, v, uu, vv) (EL(uu, vv) + dENL(u, v, uu, vv)) //
macro ddE(u, v, uu, vv, uuu, vvv) dENL(uuu, vvv, uu, vv) //
macro F2(u, v) (E(u, v)'*A*E(u, v)) //
macro dF2(u, v, uu, vv) (E(u, v)'*A*dE(u, v, uu, vv)*2.) //
macro ddF2(u, v, uu, vv, uuu, vvv) (
    (dE(u, v, uu, vv)'*A*dE(u, v, uuu, vvv))*2.
    + (E(u, v)'*A*ddE(u, v, uu, vv, uuu, vvv))*2.
    ) //

The Newton Method is

choose $ n=0$,and u_0,v_0 the initial displacement

  • loop:
    • find (du,dv) : solution of
    • un = un - du,\quad vn =vn - dv
    • until (du,dv) small is enough

The way to implement this algorithm in FreeFem++ is use a macro tool to implement A and F2, f, f',f''.

A macro is like in ccp preprocessor of C++, but this begin by macro and the end of the macro definition is before the comment //. In this case the macro is very useful because the type of parameter can be change. And it is easy to make automatic differentiation.

Fig. 36: The deformed domain
Mesh1
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
// Macro
macro EL(u, v) [dx(u), (dx(v)+dy(u)), dy(v)] //is [epsilon_11, 2epsilon_12, epsilon_22]

macro ENL(u, v) [
    (dx(u)*dx(u) + dx(v)*dx(v))*0.5,
    (dx(u)*dy(u) + dx(v)*dy(v)),
    (dy(u)*dy(u) + dy(v)*dy(v))*0.5
    ] //

macro dENL(u, v, uu, vv) [
    (dx(u)*dx(uu) + dx(v)*dx(vv)),
    (dx(u)*dy(uu) + dx(v)*dy(vv) + dx(uu)*dy(u) + dx(vv)*dy(v)),
    (dy(u)*dy(uu) + dy(v)*dy(vv))
    ] //

macro E(u, v) (EL(u,v) + ENL(u,v)) //is [E_11, 2E_12, E_22]
macro dE(u, v, uu, vv) (EL(uu, vv) + dENL(u, v, uu, vv)) //
macro ddE(u, v, uu, vv, uuu, vvv) dENL(uuu, vvv, uu, vv) //
macro F2(u, v) (E(u, v)'*A*E(u, v)) //
macro dF2(u, v, uu, vv) (E(u, v)'*A*dE(u, v, uu, vv)*2.) //
macro ddF2(u, v, uu, vv, uuu, vvv) (
      (dE(u, v, uu, vv)'*A*dE(u, v, uuu, vvv))*2.
    + (E(u, v)'*A*ddE(u, v, uu, vv, uuu, vvv))*2.
    ) //

macro f(u) ((u)*(u)*0.25) //
macro df(u) ((u)*0.5) //
macro ddf(u) (0.5) //

// Parameters
real mu = 0.012e5; //kg/cm^2
real lambda = 0.4e5; //kg/cm^2
real Pa = 1e2;

// sigma = 2 mu E + lambda tr(E) Id
// A(u,v) = sigma(u):E(v)
//
// ( a b )
// ( b c )
//
// tr*Id : (a,b,c) -> (a+c,0,a+c)
// so the associed matrix is:
// ( 1 0 1 )
// ( 0 0 0 )
// ( 1 0 1 )

real a11 = 2*mu + lambda;
real a22 = mu; //because [0, 2*t12, 0]' A [0, 2*s12,0] = 2*mu*(t12*s12 + t21*s21) = 4*mu*t12*s12
real a33 = 2*mu + lambda;
real a12 = 0;
real a13 = lambda;
real a23 = 0;
// symetric part
real a21 = a12;
real a31 = a13;
real a32 = a23;

//the matrix A
func A = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]];

// Mesh
int n = 30;
int m = 10;
mesh Th = square(n, m, [x, .3*y]); //label: 1 bottom, 2 right, 3 up, 4 left;
int bottom = 1, right = 2, upper = 3, left = 4;
plot(Th);

// Fespace
fespace Wh(Th, P1dc);
Wh e2, fe2, dfe2, ddfe2;

fespace Vh(Th, [P1, P1]);
Vh [uu, vv] = [0, 0], [w, s], [un, vn] = [0, 0];

fespace Sh(Th, P1);
Sh u1, v1;

// Problem
varf vmass ([uu, vv], [w, s], solver=CG) = int2d(Th)(uu*w + vv*s);
matrix M = vmass(Vh, Vh);
problem NonLin([uu, vv], [w, s], solver=LU)
    = int2d(Th, qforder=1)( //(D^2 J(un))
           dF2(un, vn, uu, vv)*dF2(un, vn, w, s)*ddfe2
        + ddF2(un, vn, uu, vv, w, s)*ddfe2
    )
    - int1d(Th, upper)(
          Pa*s
    )
    - int2d(Th, qforder=1)( //(D J(un))
          dF2(un, vn, w, s)*dfe2
    )
    + on(right, left, uu=0, vv=0)
    ;

// Newton's method
for (int i = 0; i < 10; i++){
    cout << "Loop " << i << endl;

    // Update
    e2 = F2(un, vn);
    dfe2 = df(e2) ;
    ddfe2 = ddf(e2);
    cout << "e2 max = " <<e2[].max << ", min = " << e2[].min << endl;
    cout << "de2 max = "<< dfe2[].max << ", min = " << dfe2[].min << endl;
    cout << "dde2 max = "<< ddfe2[].max << ", min = " << ddfe2[].min << endl;

    // Solve
    NonLin;
    w[]  = M*uu[];

    // Residual
    real res = sqrt(w[]' * uu[]); //L^2 norm of [uu, vv]
    cout << " L^2 residual = " << res << endl;

    // Update
    v1 = vv;
    u1 = uu;
    cout << "u1 min = " <<u1[].min << ", u1 max = " << u1[].max << endl;
    cout << "v1 min = " <<v1[].min << ", v2 max = " << v1[].max << endl;

    // Plot
    plot([uu, vv], wait=true, cmm="uu, vv");

    // Update
    un[] -= uu[];
    plot([un, vn], wait=true, cmm="displacement");

    if (res < 1e-5) break;
}

// Plot
plot([un, vn], wait=true);

// Movemesh
mesh th1 = movemesh(Th, [x+un, y+vn]);

// Plot
plot(th1, wait=true);