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:
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:
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}\).
Recommended References
For an exposition of nonlinear elasticity and of the underlying linear and tensor algebra, see [OGDEN1984]. For an advanced mathematical analysis of the Finite Element Method, see [RAVIART1998].
A Neo-Hookean Compressible Material
Constitutive Theory and Tangent Stress Measures
The strain energy density function is given by:
(see [HORGAN2004], formula (12)).
The corresponding 2nd Piola-Kirchoff stress tensor is given by:
The Kirchhoff stress, then, is:
The tangent Kirchhoff stress tensor at \(\bF_{n}\) acting on \(\delta \bF_{n+1}\) is, consequently:
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:
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:
by solving the weak formulation:
where we have taken:
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:
will be implemented as the macros \(<TanK11>\), \(<TanK12>\), …
The individual components of the tensor quantities:
and
will be implemented as the macros:
respectively.
In the above notation, the tangent Kirchhoff stress term becomes
while the weak BVP formulation acquires the form:
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]);