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)
  3macro grad11(u1, u2) (dx(u1)) //
  4macro grad21(u1, u2) (dy(u1)) //
  5macro grad12(u1, u2) (dx(u2)) //
  6macro grad22(u1, u2) (dy(u2)) //
  7
  8//Macros for the deformation gradient
  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) (
178      Finv11(u1, u2)*grad11(w1, w2)
179    + Finv21(u1, u2)*grad12(w1, w2)
180)//
181
182macro d4Aux12 (w1, w2, u1, u2) (
183      Finv12(u1, u2)*grad11(w1, w2)
184    + Finv22(u1, u2)*grad12(w1, w2)
185)//
186
187macro d4Aux21 (w1, w2, u1, u2) (
188      Finv11(u1, u2)*grad21(w1, w2)
189    + Finv21(u1, u2)*grad22(w1, w2)
190)//
191
192macro 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
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
242real InnerRadius = 1.e0; //The wound radius
243real OuterRadius = 4.e0; //The outer (truncated) radius
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
250border InnerEdge(t=0, 2.*pi){x=(1.0 + InnerEllipseExtension)*InnerRadius*cos(t); y=InnerRadius*sin(t); label=1;}
251border OuterEdge(t=0, 2.*pi){x=(1.0 + 0.0*InnerEllipseExtension)*OuterRadius*cos(t); y=OuterRadius*sin(t); label=2;}
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]);
Table of content