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

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 // Macro
  2 macro EL(u, v) [dx(u), (dx(v)+dy(u)), dy(v)] //is [epsilon_11, 2epsilon_12, epsilon_22]
  3 
  4 macro 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 
 10 macro 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 
 16 macro E(u, v) (EL(u,v) + ENL(u,v)) //is [E_11, 2E_12, E_22]
 17 macro dE(u, v, uu, vv) (EL(uu, vv) + dENL(u, v, uu, vv)) //
 18 macro ddE(u, v, uu, vv, uuu, vvv) dENL(uuu, vvv, uu, vv) //
 19 macro F2(u, v) (E(u, v)'*A*E(u, v)) //
 20 macro dF2(u, v, uu, vv) (E(u, v)'*A*dE(u, v, uu, vv)*2.) //
 21 macro 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 
 26 macro f(u) ((u)*(u)*0.25) //
 27 macro df(u) ((u)*0.5) //
 28 macro ddf(u) (0.5) //
 29 
 30 // Parameters
 31 real mu = 0.012e5; //kg/cm^2
 32 real lambda = 0.4e5; //kg/cm^2
 33 real 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 
 47 real a11 = 2*mu + lambda;
 48 real a22 = mu; //because [0, 2*t12, 0]' A [0, 2*s12,0] = 2*mu*(t12*s12 + t21*s21) = 4*mu*t12*s12
 49 real a33 = 2*mu + lambda;
 50 real a12 = 0;
 51 real a13 = lambda;
 52 real a23 = 0;
 53 // symetric part
 54 real a21 = a12;
 55 real a31 = a13;
 56 real a32 = a23;
 57 
 58 //the matrix A
 59 func A = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]];
 60 
 61 // Mesh
 62 int n = 30;
 63 int m = 10;
 64 mesh Th = square(n, m, [x, .3*y]); //label: 1 bottom, 2 right, 3 up, 4 left;
 65 int bottom = 1, right = 2, upper = 3, left = 4;
 66 plot(Th);
 67 
 68 // Fespace
 69 fespace Wh(Th, P1dc);
 70 Wh e2, fe2, dfe2, ddfe2;
 71 
 72 fespace Vh(Th, [P1, P1]);
 73 Vh [uu, vv] = [0, 0], [w, s], [un, vn] = [0, 0];
 74 
 75 fespace Sh(Th, P1);
 76 Sh u1, v1;
 77 
 78 // Problem
 79 varf vmass ([uu, vv], [w, s], solver=CG) = int2d(Th)(uu*w + vv*s);
 80 matrix M = vmass(Vh, Vh);
 81 problem 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
 96 for (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
132 plot([un, vn], wait=true);
133 
134 // Movemesh
135 mesh th1 = movemesh(Th, [x+un, y+vn]);
136 
137 // Plot
138 plot(th1, wait=true);
Table of content