FreeFEM Documentation on GitHub

stars - forks

Compressible Neo-Hookean materials

Author: Alex Sadovsky

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)
  3 macro grad11(u1, u2) (dx(u1)) //
  4 macro grad21(u1, u2) (dy(u1)) //
  5 macro grad12(u1, u2) (dx(u2)) //
  6 macro grad22(u1, u2) (dy(u2)) //
  7 
  8 //Macros for the deformation gradient
  9 macro F11(u1, u2) (1.0 + grad11(u1, u2)) //
 10 macro F12(u1, u2) (0.0 + grad12(u1, u2)) //
 11 macro F21(u1, u2) (0.0 + grad21(u1, u2)) //
 12 macro F22(u1, u2) (1.0 + grad22(u1, u2)) //
 13 
 14 //Macros for the incremental deformation gradient
 15 macro dF11(varu1, varu2) (grad11(varu1, varu2)) //
 16 macro dF12(varu1, varu2) (grad12(varu1, varu2)) //
 17 macro dF21(varu1, varu2) (grad21(varu1, varu2)) //
 18 macro dF22(varu1, varu2) (grad22(varu1, varu2)) //
 19 
 20 //Macro for the determinant of the deformation gradient
 21 macro 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
 27 macro Finv11 (u1, u2) (
 28       F22(u1, u2) / J(u1, u2)
 29 ) //
 30 macro Finv22 (u1, u2) (
 31       F11(u1, u2) / J(u1, u2)
 32 ) //
 33 macro Finv12 (u1, u2) (
 34     - F12(u1, u2) / J(u1, u2)
 35 ) //
 36 macro Finv21 (u1, u2) (
 37     - F21(u1, u2) / J(u1, u2)
 38 ) //
 39 
 40 //Macros for the square of the inverse of the deformation gradient
 41 macro FFinv11 (u1, u2) (
 42       Finv11(u1, u2)^2
 43     + Finv12(u1, u2)*Finv21(u1, u2)
 44 ) //
 45 
 46 macro FFinv12 (u1, u2) (
 47       Finv12(u1, u2)*(Finv11(u1, u2)
 48     + Finv22(u1, u2))
 49 ) //
 50 
 51 macro FFinv21 (u1, u2) (
 52       Finv21(u1, u2)*(Finv11(u1, u2)
 53     + Finv22(u1, u2))
 54 ) //
 55 
 56 macro 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
 62 macro FinvT11(u1, u2) (Finv11(u1, u2)) //
 63 macro FinvT12(u1, u2) (Finv21(u1, u2)) //
 64 macro FinvT21(u1, u2) (Finv12(u1, u2)) //
 65 macro FinvT22(u1, u2) (Finv22(u1, u2)) //
 66 
 67 //The left Cauchy-Green strain tensor
 68 macro B11(u1, u2) (
 69       F11(u1, u2)^2 + F12(u1, u2)^2
 70 )//
 71 
 72 macro B12(u1, u2) (
 73       F11(u1, u2)*F21(u1, u2)
 74     + F12(u1, u2)*F22(u1, u2)
 75 )//
 76 
 77 macro B21(u1, u2) (
 78       F11(u1, u2)*F21(u1, u2)
 79     + F12(u1, u2)*F22(u1, u2)
 80 )//
 81 
 82 macro 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
 88 macro d0Aux11 (u1, u2, varu1, varu2) (
 89       dF11(varu1, varu2) * F11(u1, u2)
 90     + dF12(varu1, varu2) * F12(u1, u2)
 91 )//
 92 
 93 macro d0Aux12 (u1, u2, varu1, varu2) (
 94       dF21(varu1, varu2) * F11(u1, u2)
 95     + dF22(varu1, varu2) * F12(u1, u2)
 96 )//
 97 
 98 macro d0Aux21 (u1, u2, varu1, varu2) (
 99       dF11(varu1, varu2) * F21(u1, u2)
100     + dF12(varu1, varu2) * F22(u1, u2)
101 )//
102 
103 macro 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
109 macro d1Aux11 (u1, u2, varu1, varu2) (
110       2.0 * d0Aux11 (u1, u2, varu1, varu2)
111 )//
112 
113 macro d1Aux12 (u1, u2, varu1, varu2) (
114       d0Aux12 (u1, u2, varu1, varu2)
115     + d0Aux21 (u1, u2, varu1, varu2)
116 )//
117 
118 macro d1Aux21 (u1, u2, varu1, varu2) (
119       d1Aux12 (u1, u2, varu1, varu2)
120 )//
121 
122 macro 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}
127 macro d2Aux11 (u1, u2, varu1, varu2) (
128       dF11(varu1, varu2) * FinvT11(u1, u2)
129     + dF21(varu1, varu2) * FinvT12(u1, u2)
130 )//
131 
132 macro d2Aux12 (u1, u2, varu1, varu2) (
133       dF12(varu1, varu2) * FinvT11(u1, u2)
134     + dF22(varu1, varu2) * FinvT12(u1, u2)
135 )//
136 
137 macro d2Aux21 (u1, u2, varu1, varu2) (
138       dF11(varu1, varu2) * FinvT21(u1, u2)
139     + dF21(varu1, varu2) * FinvT22(u1, u2)
140 )//
141 
142 macro 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}
148 macro 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 
155 macro 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 
162 macro 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 
169 macro 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
177 macro d4Aux11 (w1, w2, u1, u2) (
178       Finv11(u1, u2)*grad11(w1, w2)
179     + Finv21(u1, u2)*grad12(w1, w2)
180 )//
181 
182 macro d4Aux12 (w1, w2, u1, u2) (
183       Finv12(u1, u2)*grad11(w1, w2)
184     + Finv22(u1, u2)*grad12(w1, w2)
185 )//
186 
187 macro d4Aux21 (w1, w2, u1, u2) (
188       Finv11(u1, u2)*grad21(w1, w2)
189     + Finv21(u1, u2)*grad22(w1, w2)
190 )//
191 
192 macro d4Aux22 (w1, w2, u1, u2) (
193       Finv12(u1, u2)*grad21(w1, w2)
194     + Finv22(u1, u2)*grad22(w1, w2)
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
200 macro StressK11(u1, u2) (
201       mu * (B11(u1, u2) - 1.0)
202 )//
203 
204 //The Kirchhoff stress tensor
205 macro StressK12(u1, u2) (
206       mu * B12(u1, u2)
207 )//
208 
209 //The Kirchhoff stress tensor
210 macro StressK21(u1, u2) (
211       mu * B21(u1, u2)
212 )//
213 
214 //The Kirchhoff stress tensor
215 macro StressK22(u1, u2) (
216       mu * (B22(u1, u2) - 1.0)
217 )//
218 
219 //The tangent Kirchhoff stress tensor
220 macro TanK11(u1, u2, varu1, varu2) (
221       mu * d1Aux11(u1, u2, varu1, varu2)
222 )//
223 
224 macro TanK12(u1, u2, varu1, varu2) (
225       mu * d1Aux12(u1, u2, varu1, varu2)
226 )//
227 
228 macro TanK21(u1, u2, varu1, varu2) (
229       mu * d1Aux21(u1, u2, varu1, varu2)
230 )//
231 
232 macro 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
238 real mu = 5.e2; //Elastic coefficients (kg/cm^2)
239 real D = 1.e3; //(1 / compressibility)
240 real Pa = -3.e2; //Stress loads
241 
242 real InnerRadius = 1.e0; //The wound radius
243 real OuterRadius = 4.e0; //The outer (truncated) radius
244 real tol = 1.e-4; //Tolerance (L^2)
245 real InnerEllipseExtension = 1.e0; //Extension of the inner ellipse ((major axis) - (minor axis))
246 
247 int m = 40, n = 20;
248 
249 // Mesh
250 border InnerEdge(t=0, 2.*pi){x=(1.0 + InnerEllipseExtension)*InnerRadius*cos(t); y=InnerRadius*sin(t); label=1;}
251 border OuterEdge(t=0, 2.*pi){x=(1.0 + 0.0*InnerEllipseExtension)*OuterRadius*cos(t); y=OuterRadius*sin(t); label=2;}
252 mesh Th = buildmesh(InnerEdge(-m) + OuterEdge(n));
253 int bottom = 1, right = 2, upper = 3, left = 4;
254 plot(Th);
255 
256 // Fespace
257 fespace Wh(Th, P1dc);
258 fespace Vh(Th, [P1, P1]);
259 Vh [w1, w2], [u1n,u2n], [varu1, varu2];
260 Vh [ehat1x, ehat1y], [ehat2x, ehat2y];
261 Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress
262 Vh [ef1, ef2];
263 
264 fespace Sh(Th, P1);
265 Sh p, ppp;
266 Sh StrK11, StrK12, StrK21, StrK22;
267 Sh u1, u2;
268 
269 // Problem
270 varf vfMass1D(p, q) = int2d(Th)(p*q);
271 matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);
272 
273 p[] = 1;
274 ppp[] = Mass1D * p[];
275 
276 real DomainMass = ppp[].sum;
277 cout << "DomainMass = " << DomainMass << endl;
278 
279 varf vmass ([u1, u2], [v1, v2], solver=CG)
280     = int2d(Th)( (u1*v1 + u2*v2) / DomainMass );
281 
282 matrix Mass = vmass(Vh, Vh);
283 
284 matrix 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 
290 real ContParam, dContParam;
291 
292 problem 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
324 matrix auxMat;
325 
326 // Newton's method
327 ContParam = 0.;
328 dContParam = 0.01;
329 
330 //Initialization:
331 [varu1, varu2] = [0., 0.];
332 [u1n, u2n] = [0., 0.];
333 real res = 2. * tol;
334 real eforceres;
335 int loopcount = 0;
336 int loopmax = 45;
337 
338 // Iterations
339 while (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
363 plot(Th, wait=true);
364 
365 // Movemesh
366 mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);
367 
368 // Plot
369 plot(Th1, wait=true);
370 plot([u1n,u2n]);
Table of content