FreeFEM Documentation on GitHub

stars - forks

# Compressible Neo-Hookean materials

## Notation

In what follows, the symbols $$\mathbf{u}, \bF, \bB, \bC, \stress$$ denote, respectively, the displacement field, the deformation gradient, the left Cauchy-Green strain tensor $$\bB = \bF \bF^T$$, the right Cauchy-Green strain tensor $$\bC =\bF^T \bF$$, and the Cauchy stress tensor.

We also introduce the symbols $$I_1 := \tr \bC$$ and $$J := \det\bF$$. Use will be made of the identity:

${\p J \over \p \bC} = J \bC^{-1}$

The symbol $$\Id$$ denotes the identity tensor. The symbol $$\Omega_{0}$$ denotes the reference configuration of the body to be deformed. The unit volume in the reference (resp., deformed) configuration is denoted $$dV$$ (resp., $$dV_{0}$$); these two are related by:

$dV = J dV_{0},$

which allows an integral over $$\Omega$$ involving the Cauchy stress $$\bT$$ to be rewritten as an integral of the Kirchhoff stress $$\kappa = J \bT$$ over $$\Omega_{0}$$.

## A Neo-Hookean Compressible Material

Constitutive Theory and Tangent Stress Measures

The strain energy density function is given by:

$W = {\mu \over 2}(I_1 - \tr \Id - 2 \ln J)$

(see [HORGAN2004], formula (12)).

The corresponding 2nd Piola-Kirchoff stress tensor is given by:

$\bS_{n} := {\p W \over \p\bE} (\bF_{n}) = \mu (\Id - \bC^{-1})$

The Kirchhoff stress, then, is:

$\kappa = \bF \bS \bF^{T} = \mu (\bB - \Id)$

The tangent Kirchhoff stress tensor at $$\bF_{n}$$ acting on $$\delta \bF_{n+1}$$ is, consequently:

${\p \kappa \over \p \bF} (\bF_{n}) \delta \bF_{n+1} = \mu \left[ \bF_{n} (\delta \bF_{n+1})^T + \delta \bF_{n+1} (\bF_{n})^T \right]$

The Weak Form of the BVP in the Absence of Body (External) Forces

The $$\Omega_0$$ we are considering is an elliptical annulus, whose boundary consists of two concentric ellipses (each allowed to be a circle as a special case), with the major axes parallel. Let $$P$$ denote the dead stress load (traction) on a portion $$\partial \Omega_0^{t}$$ (= the inner ellipse) of the boundary $$\partial \Omega_0$$. On the rest of the boundary, we prescribe zero displacement.

The weak formulation of the boundary value problem is:

$\begin{split}\arr{lll} 0 & = & \int_{\Omega_0} \kappa[\bF] \: : \: \left\{ (\Grad \otimes \mathbf{w}) (\bF)^{-1} \right\}\\ & - & \int_{\p \Omega_0^{t}} P \cdot \hat{N}_0\\ \rra\end{split}$

For brevity, in the rest of this section we assume $$P = 0$$. The provided FreeFEM code, however, does not rely on this assumption and allows for a general value and direction of $$P$$.

Given a Newton approximation $$\mathbf{u}_n$$ of the displacement field $$\mathbf{u}$$ satisfying the BVP, we seek the correction $$\delta \mathbf{u}_{n+1}$$ to obtain a better approximation:

$\mathbf{u}_{n+1} = \mathbf{u}_{n} + \delta \mathbf{u}_{n+1}$

by solving the weak formulation:

$\begin{split}\arr{lll} 0 &=& \int_{\Omega_0}\kappa[\bF_{n} + \delta \bF_{n+1}]\: :\: \left\{(\Grad \otimes \mathbf{w}) (\bF_{n} + \delta\bF_{n+1})^{-1}\right\}- \int_{\p \Omega_0} P \cdot \hat{N}_0\\ &=& \int_{\Omega_0}\left\{\kappa[\bF_{n}] + {\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w})(\bF_{n} + \delta \bF_{n+1})^{-1}\right\}\\ &=& \int_{\Omega_0}\left\{\kappa[\bF_{n}] + {\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w}) (\bF_{n}^{-1} + \bF_{n}^{-2} \delta \bF_{n+1})\right\}\\ \\ &=& \int_{\Omega_0}\kappa[\bF_{n}]\: :\: \left\{(\Grad \otimes \mathbf{w})\bF_{n}^{-1}\right\}\\ &-& \int_{\Omega_0}\kappa[\bF_{n}]\: :\: \left\{(\Grad \otimes \mathbf{w})(\bF_{n}^{-2} \delta \bF_{n+1})\right\}\\ &+& \int_{\Omega_0}\left\{{\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w}) \bF_{n}^{-1} \right\} \\ \rra \quad \mbox{for all test functions} \mathbf{w},\end{split}$

where we have taken:

$\delta \bF_{n+1} = \Grad \otimes \delta \mathbf{u}_{n+1}$

Note

Contrary to standard notational use, the symbol $$\delta$$ here bears no variational context. By $$\delta$$ we mean simply an increment in the sense of Newton’s Method. The role of a variational virtual displacement here is played by $$\mathbf{w}$$.

## An Approach to Implementation in FreeFEM

Introducing the code-like notation, where a string in $$< >$$’s is to be read as one symbol, the individual components of the tensor:

$<TanK> := {\p \kappa \over \p \bF}[\bF_{n}] \delta \bF_{n+1}$

will be implemented as the macros $$<TanK11>$$, $$<TanK12>$$, …

The individual components of the tensor quantities:

$\bD_{1} := \bF_{n} (\delta \bF_{n+1})^T + \delta \bF_{n+1} (\bF_{n})^T,$
$\bD_{2} := \bF_{n}^{-T} \delta \bF_{n+1},$
$\bD_{3} := (\Grad \otimes \mathbf{w}) \bF_{n}^{-2} \delta \bF_{n+1},$

and

$\bD_{4} := (\Grad \otimes \mathbf{w}) \bF_{n}^{-1},$

will be implemented as the macros:

$\begin{split}\arr{l} <d1Aux11>, <d1Aux12>, \quad \ldots \quad, <d1Aux22>,\\ <d2Aux11>, <d2Aux12>, \quad \ldots \quad, <d2Aux22>\\ <d3Aux11>, <d3Aux12>, \quad \ldots \quad, <d3Aux22>\\ <d4Aux11>, <d4Aux12>, \quad \ldots \quad, <d4Aux22>\\ \rra,\end{split}$

respectively.

In the above notation, the tangent Kirchhoff stress term becomes

${\p \kappa \over \p \bF} (\bF_{n}) \: \delta \bF_{n+1} = \mu \: \bD_{1}$

while the weak BVP formulation acquires the form:

$\begin{split}\arr{lll} 0 & = & \int_{\Omega_0} \kappa[\bF_{n}] \: : \: \bD_{4} \\ &-& \int_{\Omega_0} \kappa[\bF_{n}] \: : \: \bD_{3} \\ &+& \int_{\Omega_0} \left\{ {\p \kappa \over \p \bF}[\bF_{n}] \delta \bF_{n+1} \right\} \: : \: \bD_{4} \\ \rra \quad \mbox{for all test functions} \mathbf{w}\end{split}$
  1// Macro
2//Macros for the gradient of a vector field (u1, u2)
7
9macro F11(u1, u2) (1.0 + grad11(u1, u2)) //
10macro F12(u1, u2) (0.0 + grad12(u1, u2)) //
11macro F21(u1, u2) (0.0 + grad21(u1, u2)) //
12macro F22(u1, u2) (1.0 + grad22(u1, u2)) //
13
14//Macros for the incremental deformation gradient
15macro dF11(varu1, varu2) (grad11(varu1, varu2)) //
16macro dF12(varu1, varu2) (grad12(varu1, varu2)) //
17macro dF21(varu1, varu2) (grad21(varu1, varu2)) //
18macro dF22(varu1, varu2) (grad22(varu1, varu2)) //
19
20//Macro for the determinant of the deformation gradient
21macro J(u1, u2) (
22      F11(u1, u2)*F22(u1, u2)
23    - F12(u1, u2)*F21(u1, u2)
24) //
25
26//Macros for the inverse of the deformation gradient
27macro Finv11 (u1, u2) (
28      F22(u1, u2) / J(u1, u2)
29) //
30macro Finv22 (u1, u2) (
31      F11(u1, u2) / J(u1, u2)
32) //
33macro Finv12 (u1, u2) (
34    - F12(u1, u2) / J(u1, u2)
35) //
36macro Finv21 (u1, u2) (
37    - F21(u1, u2) / J(u1, u2)
38) //
39
40//Macros for the square of the inverse of the deformation gradient
41macro FFinv11 (u1, u2) (
42      Finv11(u1, u2)^2
43    + Finv12(u1, u2)*Finv21(u1, u2)
44) //
45
46macro FFinv12 (u1, u2) (
47      Finv12(u1, u2)*(Finv11(u1, u2)
48    + Finv22(u1, u2))
49) //
50
51macro FFinv21 (u1, u2) (
52      Finv21(u1, u2)*(Finv11(u1, u2)
53    + Finv22(u1, u2))
54) //
55
56macro FFinv22 (u1, u2) (
57      Finv12(u1, u2)*Finv21(u1, u2)
58    + Finv22(u1, u2)^2
59) //
60
61//Macros for the inverse of the transpose of the deformation gradient
62macro FinvT11(u1, u2) (Finv11(u1, u2)) //
63macro FinvT12(u1, u2) (Finv21(u1, u2)) //
64macro FinvT21(u1, u2) (Finv12(u1, u2)) //
65macro FinvT22(u1, u2) (Finv22(u1, u2)) //
66
67//The left Cauchy-Green strain tensor
68macro B11(u1, u2) (
69      F11(u1, u2)^2 + F12(u1, u2)^2
70)//
71
72macro B12(u1, u2) (
73      F11(u1, u2)*F21(u1, u2)
74    + F12(u1, u2)*F22(u1, u2)
75)//
76
77macro B21(u1, u2) (
78      F11(u1, u2)*F21(u1, u2)
79    + F12(u1, u2)*F22(u1, u2)
80)//
81
82macro B22(u1, u2)(
83      F21(u1, u2)^2 + F22(u1, u2)^2
84)//
85
86//The macros for the auxiliary tensors (D0, D1, D2, ...): Begin
87//The tensor quantity D0 = F{n} (delta F{n+1})^T
88macro d0Aux11 (u1, u2, varu1, varu2) (
89      dF11(varu1, varu2) * F11(u1, u2)
90    + dF12(varu1, varu2) * F12(u1, u2)
91)//
92
93macro d0Aux12 (u1, u2, varu1, varu2) (
94      dF21(varu1, varu2) * F11(u1, u2)
95    + dF22(varu1, varu2) * F12(u1, u2)
96)//
97
98macro d0Aux21 (u1, u2, varu1, varu2) (
99      dF11(varu1, varu2) * F21(u1, u2)
100    + dF12(varu1, varu2) * F22(u1, u2)
101)//
102
103macro d0Aux22 (u1, u2, varu1, varu2) (
104      dF21(varu1, varu2) * F21(u1, u2)
105    + dF22(varu1, varu2) * F22(u1, u2)
106)//
107
108///The tensor quantity D1 = D0 + D0^T
109macro d1Aux11 (u1, u2, varu1, varu2) (
110      2.0 * d0Aux11 (u1, u2, varu1, varu2)
111)//
112
113macro d1Aux12 (u1, u2, varu1, varu2) (
114      d0Aux12 (u1, u2, varu1, varu2)
115    + d0Aux21 (u1, u2, varu1, varu2)
116)//
117
118macro d1Aux21 (u1, u2, varu1, varu2) (
119      d1Aux12 (u1, u2, varu1, varu2)
120)//
121
122macro d1Aux22 (u1, u2, varu1, varu2)(
123      2.0 * d0Aux22 (u1, u2, varu1, varu2)
124)//
125
126///The tensor quantity D2 = F^{-T}_{n} dF_{n+1}
127macro d2Aux11 (u1, u2, varu1, varu2) (
128      dF11(varu1, varu2) * FinvT11(u1, u2)
129    + dF21(varu1, varu2) * FinvT12(u1, u2)
130)//
131
132macro d2Aux12 (u1, u2, varu1, varu2) (
133      dF12(varu1, varu2) * FinvT11(u1, u2)
134    + dF22(varu1, varu2) * FinvT12(u1, u2)
135)//
136
137macro d2Aux21 (u1, u2, varu1, varu2) (
138      dF11(varu1, varu2) * FinvT21(u1, u2)
139    + dF21(varu1, varu2) * FinvT22(u1, u2)
140)//
141
142macro d2Aux22 (u1, u2, varu1, varu2) (
143      dF12(varu1, varu2) * FinvT21(u1, u2)
144    + dF22(varu1, varu2) * FinvT22(u1, u2)
145)//
146
147///The tensor quantity D3 = F^{-2}_{n} dF_{n+1}
148macro d3Aux11 (u1, u2, varu1, varu2, w1, w2) (
149      dF11(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
150    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
151    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
152    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
153)//
154
155macro d3Aux12 (u1, u2, varu1, varu2, w1, w2) (
156      dF12(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
157    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
158    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
159    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
160)//
161
162macro d3Aux21 (u1, u2, varu1, varu2, w1, w2) (
163      dF11(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
164    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
165    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
166    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
167)//
168
169macro d3Aux22 (u1, u2, varu1, varu2, w1, w2) (
170      dF12(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
171    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
172    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
173    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
174)//
175
176///The tensor quantity D4 = (grad w) * Finv
177macro d4Aux11 (w1, w2, u1, u2) (
180)//
181
182macro d4Aux12 (w1, w2, u1, u2) (
185)//
186
187macro d4Aux21 (w1, w2, u1, u2) (
190)//
191
192macro d4Aux22 (w1, w2, u1, u2) (
195)//
196//The macros for the auxiliary tensors (D0, D1, D2, ...): End
197
198//The macros for the various stress measures: BEGIN
199//The Kirchhoff stress tensor
200macro StressK11(u1, u2) (
201      mu * (B11(u1, u2) - 1.0)
202)//
203
204//The Kirchhoff stress tensor
205macro StressK12(u1, u2) (
206      mu * B12(u1, u2)
207)//
208
209//The Kirchhoff stress tensor
210macro StressK21(u1, u2) (
211      mu * B21(u1, u2)
212)//
213
214//The Kirchhoff stress tensor
215macro StressK22(u1, u2) (
216      mu * (B22(u1, u2) - 1.0)
217)//
218
219//The tangent Kirchhoff stress tensor
220macro TanK11(u1, u2, varu1, varu2) (
221      mu * d1Aux11(u1, u2, varu1, varu2)
222)//
223
224macro TanK12(u1, u2, varu1, varu2) (
225      mu * d1Aux12(u1, u2, varu1, varu2)
226)//
227
228macro TanK21(u1, u2, varu1, varu2) (
229      mu * d1Aux21(u1, u2, varu1, varu2)
230)//
231
232macro TanK22(u1, u2, varu1, varu2) (
233      mu * d1Aux22(u1, u2, varu1, varu2)
234)//
235//The macros for the stress tensor components: END
236
237// Parameters
238real mu = 5.e2; //Elastic coefficients (kg/cm^2)
239real D = 1.e3; //(1 / compressibility)
240real Pa = -3.e2; //Stress loads
241
244real tol = 1.e-4; //Tolerance (L^2)
245real InnerEllipseExtension = 1.e0; //Extension of the inner ellipse ((major axis) - (minor axis))
246
247int m = 40, n = 20;
248
249// Mesh
252mesh Th = buildmesh(InnerEdge(-m) + OuterEdge(n));
253int bottom = 1, right = 2, upper = 3, left = 4;
254plot(Th);
255
256// Fespace
257fespace Wh(Th, P1dc);
258fespace Vh(Th, [P1, P1]);
259Vh [w1, w2], [u1n,u2n], [varu1, varu2];
260Vh [ehat1x, ehat1y], [ehat2x, ehat2y];
261Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress
262Vh [ef1, ef2];
263
264fespace Sh(Th, P1);
265Sh p, ppp;
266Sh StrK11, StrK12, StrK21, StrK22;
267Sh u1, u2;
268
269// Problem
270varf vfMass1D(p, q) = int2d(Th)(p*q);
271matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);
272
273p[] = 1;
274ppp[] = Mass1D * p[];
275
276real DomainMass = ppp[].sum;
277cout << "DomainMass = " << DomainMass << endl;
278
279varf vmass ([u1, u2], [v1, v2], solver=CG)
280    = int2d(Th)( (u1*v1 + u2*v2) / DomainMass );
281
282matrix Mass = vmass(Vh, Vh);
283
284matrix Id = vmass(Vh, Vh);
285
286//Define the standard Euclidean basis functions
287[ehat1x, ehat1y] = [1.0, 0.0];
288[ehat2x, ehat2y] = [0.0, 1.0];
289
290real ContParam, dContParam;
291
292problem neoHookeanInc ([varu1, varu2], [w1, w2], solver=LU)
293    = int2d(Th, qforder=1)(
294        - (
295              StressK11 (u1n, u2n) * d3Aux11(u1n, u2n, varu1, varu2, w1, w2)
296            + StressK12 (u1n, u2n) * d3Aux12(u1n, u2n, varu1, varu2, w1, w2)
297            + StressK21 (u1n, u2n) * d3Aux21(u1n, u2n, varu1, varu2, w1, w2)
298            + StressK22 (u1n, u2n) * d3Aux22(u1n, u2n, varu1, varu2, w1, w2)
299        )
300        + TanK11 (u1n, u2n, varu1, varu2) * d4Aux11(w1, w2, u1n, u2n)
301        + TanK12 (u1n, u2n, varu1, varu2) * d4Aux12(w1, w2, u1n, u2n)
302        + TanK21 (u1n, u2n, varu1, varu2) * d4Aux21(w1, w2, u1n, u2n)
303        + TanK22 (u1n, u2n, varu1, varu2) * d4Aux22(w1, w2, u1n, u2n)
304    )
305    + int2d(Th, qforder=1)(
306          StressK11 (u1n, u2n) * d4Aux11(w1, w2, u1n, u2n)
307        + StressK12 (u1n, u2n) * d4Aux12(w1, w2, u1n, u2n)
308        + StressK21 (u1n, u2n) * d4Aux21(w1, w2, u1n, u2n)
309        + StressK22 (u1n, u2n) * d4Aux22(w1, w2, u1n, u2n)
310    )
311    //Choose one of the following two boundary conditions involving Pa:
312    // Load vectors normal to the boundary:
313    - int1d(Th, 1)(
314          Pa * (w1*N.x + w2*N.y)
315    )
316    //Load vectors tangential to the boundary:
317    //- int1d(Th, 1)(
318    //    Pa * (w1*N.y - w2*N.x)
319    //)
320    + on(2, varu1=0, varu2=0)
321    ;
322
323//Auxiliary variables
324matrix auxMat;
325
326// Newton's method
327ContParam = 0.;
328dContParam = 0.01;
329
330//Initialization:
331[varu1, varu2] = [0., 0.];
332[u1n, u2n] = [0., 0.];
333real res = 2. * tol;
334real eforceres;
335int loopcount = 0;
336int loopmax = 45;
337
338// Iterations
339while (loopcount <= loopmax && res >= tol){
340    loopcount ++;
341    cout << "Loop " << loopcount << endl;
342
343    // Solve
344    neoHookeanInc;
345
346    // Update
347    u1 = varu1;
348    u2 = varu2;
349
350    // Residual
351    w1[] = Mass*varu1[];
352    res = sqrt(w1[]' * varu1[]); //L^2 norm of [varu1, varu2]
353    cout << " L^2 residual = " << res << endl;
354
355    // Newton
356    u1n[] += varu1[];
357
358    // Plot
359    plot([u1n,u2n], cmm="displacement");
360}
361
362// Plot
363plot(Th, wait=true);
364
365// Movemesh
366mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);
367
368// Plot
369plot(Th1, wait=true);
370plot([u1n,u2n]);