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.

  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