# 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
  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 = " <