FreeFEM Documentation on GitHub

stars - forks

Non-linear elasticity

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

\[\min J(u_{1},u_{2}) = \int_{\Omega} f(F2) - \int_{\Gamma_{p}} P_{a} \, u_{2}\]

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:

\[E_{ij} = 0.5 \big( ( \p_i u_j + \p_j u_i ) + \sum_k \p_i u_k {\times} \p_j u_k \big)\]

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:

\[DJ(\mathbf{u})(\mathbf{v}) = \int DF2(\mathbf{u})(\mathbf{v}) \;f'(F2(\mathbf{u}))) - \int_{\Gamma_{p}} P_{a} v_{2}\]

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{split}\begin{array}{rcl} 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{array}\end{split}\]

where:

\[D^2F2(\mathbf{u})(\mathbf{v},\mathbf{w}) = 2 \; A(\;D^2E[\mathbf{u}](\mathbf{v},\mathbf{w})\;,\;E[\mathbf{u}]\;) + 2 \; A(\;DE[\mathbf{u}](\mathbf{v})\;,DE[\mathbf{u}](\mathbf{w})\;) .\]

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

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.

../_images/NonLinearElasticity_Mesh1.png

Fig. 180 The deformed domain

  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);
Table of content