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 macro EL(u, v) [dx(u), (dx(v)+dy(u)), dy(v)] //is [epsilon_11, 2epsilon_12, epsilon_22]
 2
 3 macro ENL(u, v) [
 4     (dx(u)*dx(u) + dx(v)*dx(v))*0.5,
 5     (dx(u)*dy(u) + dx(v)*dy(v)),
 6     (dy(u)*dy(u) + dy(v)*dy(v))*0.5
 7 ] //
 8
 9 macro dENL(u, v, uu, vv) [
10     (dx(u)*dx(uu) + dx(v)*dx(vv)),
11     (dx(u)*dy(uu) + dx(v)*dy(vv) + dx(uu)*dy(u) + dx(vv)*dy(v)),
12     (dy(u)*dy(uu) + dy(v)*dy(vv))
13 ] //
14
15 macro E(u, v) (EL(u,v) + ENL(u,v)) //is [E_11, 2E_12, E_22]
16 macro dE(u, v, uu, vv) (EL(uu, vv) + dENL(u, v, uu, vv)) //
17 macro ddE(u, v, uu, vv, uuu, vvv) dENL(uuu, vvv, uu, vv) //
18 macro F2(u, v) (E(u, v)'*A*E(u, v)) //
19 macro dF2(u, v, uu, vv) (E(u, v)'*A*dE(u, v, uu, vv)*2.) //
20 macro ddF2(u, v, uu, vv, uuu, vvv) (
21     (dE(u, v, uu, vv)'*A*dE(u, v, uuu, vvv))*2.
22     + (E(u, v)'*A*ddE(u, v, uu, vv, uuu, vvv))*2.
23 ) //

The Newton Method is:

choose \(n=0\),and \(u_0,v_0\) the initial displacement

  • loop:

    • find \((du,dv)\) : solution of

      \[D^2J(u_n,v_n)((w,s),(du,dv)) = DJ(u_n,v_n)(w,s) , \quad \forall w,s\]
    • \(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.

../_images/NonLinearElasticity_Mesh1.png

Fig. 182 The deformed domain

  1// Macro
  2macro EL(u, v) [dx(u), (dx(v)+dy(u)), dy(v)] //is [epsilon_11, 2epsilon_12, epsilon_22]
  3
  4macro ENL(u, v) [
  5    (dx(u)*dx(u) + dx(v)*dx(v))*0.5,
  6    (dx(u)*dy(u) + dx(v)*dy(v)),
  7    (dy(u)*dy(u) + dy(v)*dy(v))*0.5
  8    ] //
  9
 10macro dENL(u, v, uu, vv) [
 11    (dx(u)*dx(uu) + dx(v)*dx(vv)),
 12    (dx(u)*dy(uu) + dx(v)*dy(vv) + dx(uu)*dy(u) + dx(vv)*dy(v)),
 13    (dy(u)*dy(uu) + dy(v)*dy(vv))
 14    ] //
 15
 16macro E(u, v) (EL(u,v) + ENL(u,v)) //is [E_11, 2E_12, E_22]
 17macro dE(u, v, uu, vv) (EL(uu, vv) + dENL(u, v, uu, vv)) //
 18macro ddE(u, v, uu, vv, uuu, vvv) dENL(uuu, vvv, uu, vv) //
 19macro F2(u, v) (E(u, v)'*A*E(u, v)) //
 20macro dF2(u, v, uu, vv) (E(u, v)'*A*dE(u, v, uu, vv)*2.) //
 21macro ddF2(u, v, uu, vv, uuu, vvv) (
 22      (dE(u, v, uu, vv)'*A*dE(u, v, uuu, vvv))*2.
 23    + (E(u, v)'*A*ddE(u, v, uu, vv, uuu, vvv))*2.
 24    ) //
 25
 26macro f(u) ((u)*(u)*0.25) //
 27macro df(u) ((u)*0.5) //
 28macro ddf(u) (0.5) //
 29
 30// Parameters
 31real mu = 0.012e5; //kg/cm^2
 32real lambda = 0.4e5; //kg/cm^2
 33real Pa = 1e2;
 34
 35// sigma = 2 mu E + lambda tr(E) Id
 36// A(u,v) = sigma(u):E(v)
 37//
 38// ( a b )
 39// ( b c )
 40//
 41// tr*Id : (a,b,c) -> (a+c,0,a+c)
 42// so the associed matrix is:
 43// ( 1 0 1 )
 44// ( 0 0 0 )
 45// ( 1 0 1 )
 46
 47real a11 = 2*mu + lambda;
 48real a22 = mu; //because [0, 2*t12, 0]' A [0, 2*s12,0] = 2*mu*(t12*s12 + t21*s21) = 4*mu*t12*s12
 49real a33 = 2*mu + lambda;
 50real a12 = 0;
 51real a13 = lambda;
 52real a23 = 0;
 53// symetric part
 54real a21 = a12;
 55real a31 = a13;
 56real a32 = a23;
 57
 58//the matrix A
 59func A = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]];
 60
 61// Mesh
 62int n = 30;
 63int m = 10;
 64mesh Th = square(n, m, [x, .3*y]); //label: 1 bottom, 2 right, 3 up, 4 left;
 65int bottom = 1, right = 2, upper = 3, left = 4;
 66plot(Th);
 67
 68// Fespace
 69fespace Wh(Th, P1dc);
 70Wh e2, fe2, dfe2, ddfe2;
 71
 72fespace Vh(Th, [P1, P1]);
 73Vh [uu, vv] = [0, 0], [w, s], [un, vn] = [0, 0];
 74
 75fespace Sh(Th, P1);
 76Sh u1, v1;
 77
 78// Problem
 79varf vmass ([uu, vv], [w, s], solver=CG) = int2d(Th)(uu*w + vv*s);
 80matrix M = vmass(Vh, Vh);
 81problem NonLin([uu, vv], [w, s], solver=LU)
 82    = int2d(Th, qforder=1)( //(D^2 J(un))
 83           dF2(un, vn, uu, vv)*dF2(un, vn, w, s)*ddfe2
 84        + ddF2(un, vn, uu, vv, w, s)*ddfe2
 85    )
 86    - int1d(Th, upper)(
 87          Pa*s
 88    )
 89    - int2d(Th, qforder=1)( //(D J(un))
 90          dF2(un, vn, w, s)*dfe2
 91    )
 92    + on(right, left, uu=0, vv=0)
 93    ;
 94
 95// Newton's method
 96for (int i = 0; i < 10; i++){
 97    cout << "Loop " << i << endl;
 98
 99    // Update
100    e2 = F2(un, vn);
101    dfe2 = df(e2) ;
102    ddfe2 = ddf(e2);
103    cout << "e2 max = " <<e2[].max << ", min = " << e2[].min << endl;
104    cout << "de2 max = "<< dfe2[].max << ", min = " << dfe2[].min << endl;
105    cout << "dde2 max = "<< ddfe2[].max << ", min = " << ddfe2[].min << endl;
106
107    // Solve
108    NonLin;
109    w[]  = M*uu[];
110
111    // Residual
112    real res = sqrt(w[]' * uu[]); //L^2 norm of [uu, vv]
113    cout << " L^2 residual = " << res << endl;
114
115    // Update
116    v1 = vv;
117    u1 = uu;
118    cout << "u1 min = " <<u1[].min << ", u1 max = " << u1[].max << endl;
119    cout << "v1 min = " <<v1[].min << ", v2 max = " << v1[].max << endl;
120
121    // Plot
122    plot([uu, vv], wait=true, cmm="uu, vv");
123
124    // Update
125    un[] -= uu[];
126    plot([un, vn], wait=true, cmm="displacement");
127
128    if (res < 1e-5) break;
129}
130
131// Plot
132plot([un, vn], wait=true);
133
134// Movemesh
135mesh th1 = movemesh(Th, [x+un, y+vn]);
136
137// Plot
138plot(th1, wait=true);
Table of content