1<br> 2 3<i>This program was contributed by Jean-Paul Pelteret and Andrew McBride. 4<br> 5This material is based upon work supported by the German Science Foundation (Deutsche 6Forschungsgemeinschaft, DFG), grant STE 544/39-1, and the National Research Foundation of South Africa. 7</i> 8 9@dealiiTutorialDOI{10.5281/zenodo.439772,https://zenodo.org/badge/DOI/10.5281/zenodo.439772.svg} 10 11<a name="Intro"></a> 12<h1>Introduction</h1> 13 14The subject of this tutorial is nonlinear solid mechanics. 15Classical single-field approaches (see e.g. step-18) can not correctly describe the response of quasi-incompressible materials. 16The response is overly stiff; a phenomenon known as locking. 17Locking problems can be circumvented using a variety of alternative strategies. 18One such strategy is the three-field formulation. 19It is used here to model the three-dimensional, fully-nonlinear (geometrical and material) response of an isotropic continuum body. 20The material response is approximated as hyperelastic. 21Additionally, the three-field formulation employed is valid for quasi-incompressible as well as compressible materials. 22 23The objective of this presentation is to provide a basis for using deal.II for problems in nonlinear solid mechanics. 24The linear problem was addressed in step-8. 25A non-standard, hypoelastic-type form of the geometrically nonlinear problem was partially considered in step-18: a rate form of the linearised constitutive relations is used and the problem domain evolves with the motion. 26Important concepts surrounding the nonlinear kinematics are absent in the theory and implementation. 27Step-18 does, however, describe many of the key concepts to implement elasticity within the framework of deal.II. 28 29We begin with a crash-course in nonlinear kinematics. 30For the sake of simplicity, we restrict our attention to the quasi-static problem. 31Thereafter, various key stress measures are introduced and the constitutive model described. 32We then describe the three-field formulation in detail prior to explaining the structure of the class used to manage the material. 33The setup of the example problem is then presented. 34 35@note This tutorial has been developed (and is described in the introduction) for the problem of elasticity in three dimensions. 36 While the space dimension could be changed in the main() routine, care needs to be taken. 37 Two-dimensional elasticity problems, in general, exist only as idealizations of three-dimensional ones. 38 That is, they are either plane strain or plane stress. 39 The assumptions that follow either of these choices needs to be consistently imposed. 40 For more information see the note in step-8. 41 42<h3>List of references</h3> 43 44The three-field formulation implemented here was pioneered by Simo et al. (1985) and is known as the mixed Jacobian-pressure formulation. 45Important related contributions include those by Simo and Taylor (1991), and Miehe (1994). 46The notation adopted here draws heavily on the excellent overview of the theoretical aspects of nonlinear solid mechanics by Holzapfel (2001). 47A nice overview of issues pertaining to incompressible elasticity (at small strains) is given in Hughes (2000). 48 49<ol> 50 <li> J.C. Simo, R.L. Taylor and K.S. Pister (1985), 51 Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, 52 <em> Computer Methods in Applied Mechanics and Engineering </em>, 53 <strong> 51 </strong>, 1-3, 54 177-208. 55 DOI: <a href="http://doi.org/10.1016/0045-7825(85)90033-7">10.1016/0045-7825(85)90033-7</a>; 56 <li> J.C. Simo and R.L. Taylor (1991), 57 Quasi-incompressible finite elasticity in principal stretches. Continuum 58 basis and numerical algorithms, 59 <em> Computer Methods in Applied Mechanics and Engineering </em>, 60 <strong> 85 </strong>, 3, 61 273-310. 62 DOI: <a href="http://doi.org/10.1016/0045-7825(91)90100-K">10.1016/0045-7825(91)90100-K</a>; 63 <li> C. Miehe (1994), 64 Aspects of the formulation and finite element implementation of large strain isotropic elasticity 65 <em> International Journal for Numerical Methods in Engineering </em> 66 <strong> 37 </strong>, 12, 67 1981-2004. 68 DOI: <a href="http://doi.org/10.1002/nme.1620371202">10.1002/nme.1620371202</a>; 69 <li> G.A. Holzapfel (2001), 70 Nonlinear Solid Mechanics. A Continuum Approach for Engineering, 71 John Wiley & Sons. 72 ISBN: 0-471-82304-X; 73 <li> T.J.R. Hughes (2000), 74 The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, 75 Dover. 76 ISBN: 978-0486411811 77</ol> 78 79An example where this three-field formulation is used in a coupled problem is documented in 80<ol> 81 <li> J-P. V. Pelteret, D. Davydov, A. McBride, D. K. Vu, and P. Steinmann (2016), 82 Computational electro- and magneto-elasticity for quasi-incompressible media immersed in free space, 83 <em> International Journal for Numerical Methods in Engineering </em>. 84 DOI: <a href="http://doi.org/10.1002/nme.5254">10.1002/nme.5254</a> 85</ol> 86 87<h3> Notation </h3> 88 89One can think of fourth-order tensors as linear operators mapping second-order 90tensors (matrices) onto themselves in much the same way as matrices map 91vectors onto vectors. 92There are various fourth-order unit tensors that will be required in the forthcoming presentation. 93The fourth-order unit tensors $\mathcal{I}$ and $\overline{\mathcal{I}}$ are defined by 94@f[ 95 \mathbf{A} = \mathcal{I}:\mathbf{A} 96 \qquad \text{and} \qquad 97 \mathbf{A}^T = \overline{\mathcal{I}}:\mathbf{A} \, . 98@f] 99Note $\mathcal{I} \neq \overline{\mathcal{I}}^T$. 100Furthermore, we define the symmetric and skew-symmetric fourth-order unit tensors by 101@f[ 102 \mathcal{S} \dealcoloneq \dfrac{1}{2}[\mathcal{I} + \overline{\mathcal{I}}] 103 \qquad \text{and} \qquad 104 \mathcal{W} \dealcoloneq \dfrac{1}{2}[\mathcal{I} - \overline{\mathcal{I}}] \, , 105@f] 106such that 107@f[ 108 \dfrac{1}{2}[\mathbf{A} + \mathbf{A}^T] = \mathcal{S}:\mathbf{A} 109 \qquad \text{and} \qquad 110 \dfrac{1}{2}[\mathbf{A} - \mathbf{A}^T] = \mathcal{W}:\mathbf{A} \, . 111@f] 112The fourth-order <code>SymmetricTensor</code> returned by identity_tensor() is $\mathcal{S}$. 113 114 115<h3>Kinematics</h3> 116 117Let the time domain be denoted $\mathbb{T} = [0,T_{\textrm{end}}]$, where $t \in \mathbb{T}$ and $T_{\textrm{end}}$ is the total problem duration. 118Consider a continuum body that occupies the reference configuration $\Omega_0$ at time $t=0$. 119%Particles in the reference configuration are identified by the position vector $\mathbf{X}$. 120The configuration of the body at a later time $t>0$ is termed the current configuration, denoted $\Omega$, with particles identified by the vector $\mathbf{x}$. 121The nonlinear map between the reference and current configurations, denoted $\boldsymbol{\varphi}$, acts as follows: 122@f[ 123 \mathbf{x} = \boldsymbol{\varphi}(\mathbf{X},t) \, . 124@f] 125The material description of the displacement of a particle is defined by 126@f[ 127 \mathbf{U}(\mathbf{X},t) = \mathbf{x}(\mathbf{X},t) - \mathbf{X} \, . 128@f] 129 130The deformation gradient $\mathbf{F}$ is defined as the material gradient of the motion: 131@f[ 132 \mathbf{F}(\mathbf{X},t) 133 \dealcoloneq \dfrac{\partial \boldsymbol{\varphi}(\mathbf{X},t)}{\partial \mathbf{X}} 134 = \textrm{Grad}\ \mathbf{x}(\mathbf{X},t) 135 = \mathbf{I} + \textrm{Grad}\ \mathbf{U} \, . 136@f] 137The determinant of the of the deformation gradient 138$J(\mathbf{X},t) \dealcoloneq \textrm{det}\ \mathbf{F}(\mathbf{X},t) > 0$ 139maps corresponding volume elements in the reference and current configurations, denoted 140$\textrm{d}V$ and $\textrm{d}v$, 141respectively, as 142@f[ 143 \textrm{d}v = J(\mathbf{X},t)\; \textrm{d}V \, . 144@f] 145 146Two important measures of the deformation in terms of the spatial and material coordinates are the left and right Cauchy-Green tensors, respectively, 147and denoted $\mathbf{b} \dealcoloneq \mathbf{F}\mathbf{F}^T$ and $\mathbf{C} \dealcoloneq \mathbf{F}^T\mathbf{F}$. 148They are both symmetric and positive definite. 149 150The Green-Lagrange strain tensor is defined by 151@f[ 152 \mathbf{E} \dealcoloneq \frac{1}{2}[\mathbf{C} - \mathbf{I} ] 153 = \underbrace{\frac{1}{2}[\textrm{Grad}^T \mathbf{U} + \textrm{Grad}\mathbf{U}]}_{\boldsymbol{\varepsilon}} 154 + \frac{1}{2}[\textrm{Grad}^T\ \mathbf{U}][\textrm{Grad}\ \mathbf{U}] \, . 155@f] 156If the assumption of infinitesimal deformations is made, then the second term 157on the right can be neglected, and $\boldsymbol{\varepsilon}$ (the linearised 158strain tensor) is the only component of the strain tensor. 159This assumption is, looking at the setup of the problem, not valid in step-18, 160making the use of the linearized $\boldsymbol{\varepsilon}$ as the strain 161measure in that tutorial program questionable. 162 163In order to handle the different response that materials exhibit when subjected to bulk and shear type deformations we consider the following decomposition of the deformation gradient $\mathbf{F}$ and the left Cauchy-Green tensor $\mathbf{b}$ into volume-changing (volumetric) and volume-preserving (isochoric) parts: 164@f[ 165 \mathbf{F} 166 = (J^{1/3}\mathbf{I})\overline{\mathbf{F}} 167 \qquad \text{and} \qquad 168 \mathbf{b} 169 = (J^{2/3}\mathbf{I})\overline{\mathbf{F}}\,\overline{\mathbf{F}}^T 170 = (J^{2/3}\mathbf{I})\overline{\mathbf{b}} \, . 171@f] 172Clearly, $\textrm{det}\ \mathbf{F} = \textrm{det}\ (J^{1/3}\mathbf{I}) = J$. 173 174The spatial velocity field is denoted $\mathbf{v}(\mathbf{x},t)$. 175The derivative of the spatial velocity field with respect to the spatial coordinates gives the spatial velocity gradient $\mathbf{l}(\mathbf{x},t)$, that is 176@f[ 177 \mathbf{l}(\mathbf{x},t) 178 \dealcoloneq \dfrac{\partial \mathbf{v}(\mathbf{x},t)}{\partial \mathbf{x}} 179 = \textrm{grad}\ \mathbf{v}(\mathbf{x},t) \, , 180@f] 181where $\textrm{grad} \{\bullet \} 182= \frac{\partial \{ \bullet \} }{ \partial \mathbf{x}} 183= \frac{\partial \{ \bullet \} }{ \partial \mathbf{X}}\frac{\partial \mathbf{X} }{ \partial \mathbf{x}} 184= \textrm{Grad} \{ \bullet \} \mathbf{F}^{-1}$. 185 186 187<h3>Kinetics</h3> 188 189Cauchy's stress theorem equates the Cauchy traction $\mathbf{t}$ acting on an infinitesimal surface element in the current configuration $\mathrm{d}a$ to the product of the Cauchy stress tensor $\boldsymbol{\sigma}$ (a spatial quantity) and the outward unit normal to the surface $\mathbf{n}$ as 190@f[ 191 \mathbf{t}(\mathbf{x},t, \mathbf{n}) = \boldsymbol{\sigma}\mathbf{n} \, . 192@f] 193The Cauchy stress is symmetric. 194Similarly, the first Piola-Kirchhoff traction $\mathbf{T}$ which acts on an infinitesimal surface element in the reference configuration $\mathrm{d}A$ is the product of the first Piola-Kirchhoff stress tensor $\mathbf{P}$ (a two-point tensor) and the outward unit normal to the surface $\mathbf{N}$ as 195@f[ 196 \mathbf{T}(\mathbf{X},t, \mathbf{N}) = \mathbf{P}\mathbf{N} \, . 197@f] 198The Cauchy traction $\mathbf{t}$ and the first Piola-Kirchhoff traction $\mathbf{T}$ are related as 199@f[ 200 \mathbf{t}\mathrm{d}a = \mathbf{T}\mathrm{d}A \, . 201@f] 202This can be demonstrated using <a href="http://en.wikipedia.org/wiki/Finite_strain_theory">Nanson's formula</a>. 203 204The first Piola-Kirchhoff stress tensor is related to the Cauchy stress as 205@f[ 206 \mathbf{P} = J \boldsymbol{\sigma}\mathbf{F}^{-T} \, . 207@f] 208Further important stress measures are the (spatial) Kirchhoff stress $\boldsymbol{\tau} = J \boldsymbol{\sigma}$ 209and the (referential) second Piola-Kirchhoff stress 210$\mathbf{S} = {\mathbf{F}}^{-1} \boldsymbol{\tau} {\mathbf{F}}^{-T}$. 211 212 213<h3> Push-forward and pull-back operators </h3> 214 215Push-forward and pull-back operators allow one to transform various measures between the material and spatial settings. 216The stress measures used here are contravariant, while the strain measures are covariant. 217 218The push-forward and-pull back operations for second-order covariant tensors $(\bullet)^{\text{cov}}$ are respectively given by: 219@f[ 220 \chi_{*}(\bullet)^{\text{cov}} \dealcoloneq \mathbf{F}^{-T} (\bullet)^{\text{cov}} \mathbf{F}^{-1} 221 \qquad \text{and} \qquad 222 \chi^{-1}_{*}(\bullet)^{\text{cov}} \dealcoloneq \mathbf{F}^{T} (\bullet)^{\text{cov}} \mathbf{F} \, . 223@f] 224 225The push-forward and pull back operations for second-order contravariant tensors $(\bullet)^{\text{con}}$ are respectively given by: 226@f[ 227 \chi_{*}(\bullet)^{\text{con}} \dealcoloneq \mathbf{F} (\bullet)^{\text{con}} \mathbf{F}^T 228 \qquad \text{and} \qquad 229 \chi^{-1}_{*}(\bullet)^{\text{con}} \dealcoloneq \mathbf{F}^{-1} (\bullet)^{\text{con}} \mathbf{F}^{-T} \, . 230@f] 231For example $\boldsymbol{\tau} = \chi_{*}(\mathbf{S})$. 232 233 234<h3>Hyperelastic materials</h3> 235 236A hyperelastic material response is governed by a Helmholtz free energy function $\Psi = \Psi(\mathbf{F}) = \Psi(\mathbf{C}) = \Psi(\mathbf{b})$ which serves as a potential for the stress. 237For example, if the Helmholtz free energy depends on the right Cauchy-Green tensor $\mathbf{C}$ then the isotropic hyperelastic response is 238@f[ 239 \mathbf{S} 240 = 2 \dfrac{\partial \Psi(\mathbf{C})}{\partial \mathbf{C}} \, . 241@f] 242If the Helmholtz free energy depends on the left Cauchy-Green tensor $\mathbf{b}$ then the isotropic hyperelastic response is 243@f[ 244 \boldsymbol{\tau} 245 = 2 \dfrac{\partial \Psi(\mathbf{b})}{\partial \mathbf{b}} \mathbf{b} 246 = 2 \mathbf{b} \dfrac{\partial \Psi(\mathbf{b})}{\partial \mathbf{b}} \, . 247@f] 248 249Following the multiplicative decomposition of the deformation gradient, the Helmholtz free energy can be decomposed as 250@f[ 251 \Psi(\mathbf{b}) = \Psi_{\text{vol}}(J) + \Psi_{\text{iso}}(\overline{\mathbf{b}}) \, . 252@f] 253Similarly, the Kirchhoff stress can be decomposed into volumetric and isochoric parts as $\boldsymbol{\tau} = \boldsymbol{\tau}_{\text{vol}} + \boldsymbol{\tau}_{\text{iso}}$ where: 254@f{align*} 255 \boldsymbol{\tau}_{\text{vol}} &= 256 2 \mathbf{b} \dfrac{\partial \Psi_{\textrm{vol}}(J)}{\partial \mathbf{b}} 257 \\ 258 &= p J\mathbf{I} \, , 259 \\ 260 \boldsymbol{\tau}_{\text{iso}} &= 261 2 \mathbf{b} \dfrac{\partial \Psi_{\textrm{iso}} (\overline{\mathbf{b}})}{\partial \mathbf{b}} 262 \\ 263 &= \underbrace{( \mathcal{I} - \dfrac{1}{3} \mathbf{I} \otimes \mathbf{I})}_{\mathbb{P}} : \overline{\boldsymbol{\tau}} \, , 264@f} 265where 266$p \dealcoloneq \dfrac{\partial \Psi_{\text{vol}}(J)}{\partial J}$ is the pressure response. 267$\mathbb{P}$ is the projection tensor which provides the deviatoric operator in the Eulerian setting. 268The fictitious Kirchhoff stress tensor $\overline{\boldsymbol{\tau}}$ is defined by 269@f[ 270 \overline{\boldsymbol{\tau}} 271 \dealcoloneq 2 \overline{\mathbf{b}} \dfrac{\partial \Psi_{\textrm{iso}}(\overline{\mathbf{b}})}{\partial \overline{\mathbf{b}}} \, . 272@f] 273 274 275@note The pressure response as defined above differs from the widely-used definition of the 276pressure in solid mechanics as 277$p = - 1/3 \textrm{tr} \boldsymbol{\sigma} = - 1/3 J^{-1} \textrm{tr} \boldsymbol{\tau}$. 278Here $p$ is the hydrostatic pressure. 279We make use of the pressure response throughout this tutorial (although we refer to it as the pressure). 280 281<h4> Neo-Hookean materials </h4> 282 283The Helmholtz free energy corresponding to a compressible <a href="http://en.wikipedia.org/wiki/Neo-Hookean_solid">neo-Hookean material</a> is given by 284@f[ 285 \Psi \equiv 286 \underbrace{\kappa [ \mathcal{G}(J) ] }_{\Psi_{\textrm{vol}}(J)} 287 + \underbrace{\bigl[c_1 [ \overline{I}_1 - 3] \bigr]}_{\Psi_{\text{iso}}(\overline{\mathbf{b}})} \, , 288@f] 289where $\kappa \dealcoloneq \lambda + 2/3 \mu$ is the bulk modulus ($\lambda$ and $\mu$ are the Lamé parameters) 290and $\overline{I}_1 \dealcoloneq \textrm{tr}\ \overline{\mathbf{b}}$. 291The function $\mathcal{G}(J)$ is required to be strictly convex and satisfy the condition $\mathcal{G}(1) = 0$, 292among others, see Holzapfel (2001) for further details. 293In this work $\mathcal{G} \dealcoloneq \frac{1}{4} [ J^2 - 1 - 2\textrm{ln}J ]$. 294 295Incompressibility imposes the isochoric constraint that $J=1$ for all motions $\boldsymbol{\varphi}$. 296The Helmholtz free energy corresponding to an incompressible neo-Hookean material is given by 297@f[ 298 \Psi \equiv 299 \underbrace{\bigl[ c_1 [ I_1 - 3] \bigr] }_{\Psi_{\textrm{iso}}(\mathbf{b})} \, , 300@f] 301where $ I_1 \dealcoloneq \textrm{tr}\mathbf{b} $. 302Thus, the incompressible response is obtained by removing the volumetric component from the compressible free energy and enforcing $J=1$. 303 304 305<h3>Elasticity tensors</h3> 306 307We will use a Newton-Raphson strategy to solve the nonlinear boundary value problem. 308Thus, we will need to linearise the constitutive relations. 309 310The fourth-order elasticity tensor in the material description is defined by 311@f[ 312 \mathfrak{C} 313 = 2\dfrac{\partial \mathbf{S}(\mathbf{C})}{\partial \mathbf{C}} 314 = 4\dfrac{\partial^2 \Psi(\mathbf{C})}{\partial \mathbf{C} \partial \mathbf{C}} \, . 315@f] 316The fourth-order elasticity tensor in the spatial description $\mathfrak{c}$ is obtained from the push-forward of $\mathfrak{C}$ as 317@f[ 318 \mathfrak{c} = J^{-1} \chi_{*}(\mathfrak{C}) 319 \qquad \text{and thus} \qquad 320 J\mathfrak{c} = 4 \mathbf{b} \dfrac{\partial^2 \Psi(\mathbf{b})} {\partial \mathbf{b} \partial \mathbf{b}} \mathbf{b} \, . 321@f] 322The fourth-order elasticity tensors (for hyperelastic materials) possess both major and minor symmetries. 323 324The fourth-order spatial elasticity tensor can be written in the following decoupled form: 325@f[ 326 \mathfrak{c} = \mathfrak{c}_{\text{vol}} + \mathfrak{c}_{\text{iso}} \, , 327@f] 328where 329@f{align*} 330 J \mathfrak{c}_{\text{vol}} 331 &= 4 \mathbf{b} \dfrac{\partial^2 \Psi_{\text{vol}}(J)} {\partial \mathbf{b} \partial \mathbf{b}} \mathbf{b} 332 \\ 333 &= J[\widehat{p}\, \mathbf{I} \otimes \mathbf{I} - 2p \mathcal{I}] 334 \qquad \text{where} \qquad 335 \widehat{p} \dealcoloneq p + \dfrac{\textrm{d} p}{\textrm{d}J} \, , 336 \\ 337 J \mathfrak{c}_{\text{iso}} 338 &= 4 \mathbf{b} \dfrac{\partial^2 \Psi_{\text{iso}}(\overline{\mathbf{b}})} {\partial \mathbf{b} \partial \mathbf{b}} \mathbf{b} 339 \\ 340 &= \mathbb{P} : \mathfrak{\overline{c}} : \mathbb{P} 341 + \dfrac{2}{3}[\overline{\boldsymbol{\tau}}:\mathbf{I}]\mathbb{P} 342 - \dfrac{2}{3}[ \mathbf{I}\otimes\boldsymbol{\tau}_{\text{iso}} 343 + \boldsymbol{\tau}_{\text{iso}} \otimes \mathbf{I} ] \, , 344@f} 345where the fictitious elasticity tensor $\overline{\mathfrak{c}}$ in the spatial description is defined by 346@f[ 347 \overline{\mathfrak{c}} 348 = 4 \overline{\mathbf{b}} \dfrac{ \partial^2 \Psi_{\textrm{iso}}(\overline{\mathbf{b}})} {\partial \overline{\mathbf{b}} \partial \overline{\mathbf{b}}} \overline{\mathbf{b}} \, . 349@f] 350 351<h3>Principle of stationary potential energy and the three-field formulation</h3> 352 353The total potential energy of the system $\Pi$ is the sum of the internal and external potential energies, denoted $\Pi_{\textrm{int}}$ and $\Pi_{\textrm{ext}}$, respectively. 354We wish to find the equilibrium configuration by minimising the potential energy. 355 356As mentioned above, we adopt a three-field formulation. 357We denote the set of primary unknowns by 358$\mathbf{\Xi} \dealcoloneq \{ \mathbf{u}, \widetilde{p}, \widetilde{J} \}$. 359The independent kinematic variable $\widetilde{J}$ enters the formulation as a constraint on $J$ enforced by the Lagrange multiplier $\widetilde{p}$ (the pressure, as we shall see). 360 361The three-field variational principle used here is given by 362@f[ 363 \Pi(\mathbf{\Xi}) \dealcoloneq \int_\Omega \bigl[ 364 \Psi_{\textrm{vol}}(\widetilde{J}) 365 + \widetilde{p}\,[J(\mathbf{u}) - \widetilde{J}] 366 + \Psi_{\textrm{iso}}(\overline{\mathbf{b}}(\mathbf{u})) 367 \bigr] \textrm{d}v 368 + \Pi_{\textrm{ext}} \, , 369@f] 370where the external potential is defined by 371@f[ 372 \Pi_{\textrm{ext}} 373 = - \int_\Omega \mathbf{b}^\text{p} \cdot \mathbf{u}~\textrm{d}v 374 - \int_{\partial \Omega_{\sigma}} \mathbf{t}^\text{p} \cdot \mathbf{u}~\textrm{d}a \, . 375@f] 376The boundary of the current configuration $\partial \Omega$ is composed into two parts as 377$\partial \Omega = \partial \Omega_{\mathbf{u}} \cup \partial \Omega_{\sigma}$, 378where 379$\partial \Omega_{\mathbf{u}} \cap \partial \Omega_{\boldsymbol{\sigma}} = \emptyset$. 380The prescribed Cauchy traction, denoted $\mathbf{t}^\text{p}$, is applied to $ \partial \Omega_{\boldsymbol{\sigma}}$ while the motion is prescribed on the remaining portion of the boundary $\partial \Omega_{\mathbf{u}}$. 381The body force per unit current volume is denoted $\mathbf{b}^\text{p}$. 382 383 384 385The stationarity of the potential follows as 386@f{align*} 387 R(\mathbf\Xi;\delta \mathbf{\Xi}) 388 &= D_{\delta \mathbf{\Xi}}\Pi(\mathbf{\Xi}) 389 \\ 390 &= \dfrac{\partial \Pi(\mathbf{\Xi})}{\partial \mathbf{u}} \cdot \delta \mathbf{u} 391 + \dfrac{\partial \Pi(\mathbf{\Xi})}{\partial \widetilde{p}} \delta \widetilde{p} 392 + \dfrac{\partial \Pi(\mathbf{\Xi})}{\partial \widetilde{J}} \delta \tilde{J} 393 \\ 394 &= \int_{\Omega_0} \left[ 395 \textrm{grad}\ \delta\mathbf{u} : [ \underbrace{[\widetilde{p} J \mathbf{I}]}_{\equiv \boldsymbol{\tau}_{\textrm{vol}}} 396 + \boldsymbol{\tau}_{\textrm{iso}}] 397 + \delta \widetilde{p}\, [ J(\mathbf{u}) - \widetilde{J}] 398 + \delta \widetilde{J}\left[ \dfrac{\textrm{d} \Psi_{\textrm{vol}}(\widetilde{J})}{\textrm{d} \widetilde{J}} 399 -\widetilde{p}\right] 400 \right]~\textrm{d}V 401 \\ 402 &\quad - \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\textrm{d}V 403 - \int_{\partial \Omega_{0,\boldsymbol{\sigma}}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\textrm{d}A 404 \\ 405 &=0 \, , 406@f} 407for all virtual displacements $\delta \mathbf{u} \in H^1(\Omega)$ subject to the constraint that $\delta \mathbf{u} = \mathbf{0}$ on $\partial \Omega_{\mathbf{u}}$, and all virtual pressures $\delta \widetilde{p} \in L^2(\Omega)$ and virtual dilatations $\delta \widetilde{J} \in L^2(\Omega)$. 408 409One should note that the definitions of the volumetric Kirchhoff stress in the three field formulation 410$\boldsymbol{\tau}_{\textrm{vol}} \equiv \widetilde{p} J \mathbf{I}$ 411 and the subsequent volumetric tangent differs slightly from the general form given in the section on hyperelastic materials where 412$\boldsymbol{\tau}_{\textrm{vol}} \equiv p J\mathbf{I}$. 413This is because the pressure $\widetilde{p}$ is now a primary field as opposed to a constitutively derived quantity. 414One needs to carefully distinguish between the primary fields and those obtained from the constitutive relations. 415 416@note Although the variables are all expressed in terms of spatial quantities, the domain of integration is the initial configuration. 417This approach is called a <em> total-Lagrangian formulation </em>. 418The approach given in step-18, where the domain of integration is the current configuration, could be called an <em> updated Lagrangian formulation </em>. 419The various merits of these two approaches are discussed widely in the literature. 420It should be noted however that they are equivalent. 421 422 423The Euler-Lagrange equations corresponding to the residual are: 424@f{align*} 425 &\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p} = \mathbf{0} && \textrm{[equilibrium]} 426 \\ 427 &J(\mathbf{u}) = \widetilde{J} && \textrm{[dilatation]} 428 \\ 429 &\widetilde{p} = \dfrac{\textrm{d} \Psi_{\textrm{vol}}(\widetilde{J})}{\textrm{d} \widetilde{J}} && \textrm{[pressure]} \, . 430@f} 431The first equation is the (quasi-static) equilibrium equation in the spatial setting. 432The second is the constraint that $J(\mathbf{u}) = \widetilde{J}$. 433The third is the definition of the pressure $\widetilde{p}$. 434 435@note The simplified single-field derivation ($\mathbf{u}$ is the only primary variable) below makes it clear how we transform the limits of integration to the reference domain: 436@f{align*} 437\int_{\Omega}\delta \mathbf{u} \cdot [\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p}]~\mathrm{d}v 438&= 439\int_{\Omega} [-\mathrm{grad}\delta \mathbf{u}:\boldsymbol{\sigma} + \delta \mathbf{u} \cdot\mathbf{b}^\text{p}]~\mathrm{d}v 440 + \int_{\partial \Omega} \delta \mathbf{u} \cdot \mathbf{t}^\text{p}~\mathrm{d}a \\ 441&= 442- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V 443+ \int_{\Omega_0} \delta \mathbf{u} \cdot J\mathbf{b}^\text{p}~\mathrm{d}V 444 + \int_{\partial \Omega_0} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\ 445&= 446- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V 447+ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V 448 + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\ 449&= 450- \int_{\Omega_0} [\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} :\boldsymbol{\tau}~\mathrm{d}V 451+ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V 452 + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \, , 453@f} 454where 455$[\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} = 1/2[ \mathrm{grad}\delta\mathbf{u} + [\mathrm{grad}\delta\mathbf{u}]^T] $. 456 457We will use an iterative Newton-Raphson method to solve the nonlinear residual equation $R$. 458For the sake of simplicity we assume dead loading, i.e. the loading does not change due to the deformation. 459 460The change in a quantity between the known state at $t_{\textrm{n}-1}$ 461and the currently unknown state at $t_{\textrm{n}}$ is denoted 462$\varDelta \{ \bullet \} = { \{ \bullet \} }^{\textrm{n}} - { \{ \bullet \} }^{\textrm{n-1}}$. 463The value of a quantity at the current iteration $\textrm{i}$ is denoted 464${ \{ \bullet \} }^{\textrm{n}}_{\textrm{i}} = { \{ \bullet \} }_{\textrm{i}}$. 465The incremental change between iterations $\textrm{i}$ and $\textrm{i}+1$ is denoted 466$d \{ \bullet \} \dealcoloneq \{ \bullet \}_{\textrm{i}+1} - \{ \bullet \}_{\textrm{i}}$. 467 468Assume that the state of the system is known for some iteration $\textrm{i}$. 469The linearised approximation to nonlinear governing equations to be solved using the Newton-Raphson method is: 470Find $d \mathbf{\Xi}$ such that 471@f[ 472 R(\mathbf{\Xi}_{\mathsf{i}+1}) = 473 R(\mathbf{\Xi}_{\mathsf{i}}) 474 + D^2_{d \mathbf{\Xi}, \delta \mathbf{\Xi}} \Pi(\mathbf{\Xi_{\mathsf{i}}}) \cdot d \mathbf{\Xi} \equiv 0 \, , 475@f] 476then set 477$\mathbf{\Xi}_{\textrm{i}+1} = \mathbf{\Xi}_{\textrm{i}} 478+ d \mathbf{\Xi}$. 479The tangent is given by 480 481@f[ 482 D^2_{d \mathbf{\Xi}, \delta \mathbf{\Xi}} \Pi( \mathbf{\Xi}_{\mathsf{i}} ) 483 = D_{d \mathbf{\Xi}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) 484 =: K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi}) \, . 485@f] 486Thus, 487@f{align*} 488 K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi}) 489 &= 490 D_{d \mathbf{u}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) \cdot d \mathbf{u} 491 \\ 492 &\quad + 493 D_{d \widetilde{p}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) d \widetilde{p} 494 \\ 495 &\quad + 496 D_{d \widetilde{J}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) d \widetilde{J} \, , 497@f} 498where 499@f{align*} 500 D_{d \mathbf{u}} R( \mathbf{\Xi}; \delta \mathbf{\Xi}) 501 &= 502 \int_{\Omega_0} \bigl[ \textrm{grad}\ \delta \mathbf{u} : 503 \textrm{grad}\ d \mathbf{u} [\boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}] 504 + \textrm{grad}\ \delta \mathbf{u} :[ 505 \underbrace{[\widetilde{p}J[\mathbf{I}\otimes\mathbf{I} - 2 \mathcal{I}]}_{\equiv J\mathfrak{c}_{\textrm{vol}}} + 506 J\mathfrak{c}_{\textrm{iso}}] :\textrm{grad} d \mathbf{u} 507 \bigr]~\textrm{d}V \, , 508 \\ 509 &\quad + \int_{\Omega_0} \delta \widetilde{p} J \mathbf{I} : \textrm{grad}\ d \mathbf{u} ~\textrm{d}V 510 \\ 511 D_{d \widetilde{p}} R( \mathbf{\Xi}; \delta \mathbf{\Xi}) 512 &= 513 \int_{\Omega_0} \textrm{grad}\ \delta \mathbf{u} : J \mathbf{I} d \widetilde{p} ~\textrm{d}V 514 - \int_{\Omega_0} \delta \widetilde{J} d \widetilde{p} ~\textrm{d}V \, , 515 \\ 516 D_{d \widetilde{J}} R( \mathbf{\Xi}; \delta \mathbf{\Xi}) 517 &= -\int_{\Omega_0} \delta \widetilde{p} d \widetilde{J}~\textrm{d}V 518 + \int_{\Omega_0} \delta \widetilde{J} \dfrac{\textrm{d}^2 \Psi_{\textrm{vol}}(\widetilde{J})}{\textrm{d} \widetilde{J}\textrm{d}\widetilde{J}} d \widetilde{J} ~\textrm{d}V \, . 519@f} 520 521Note that the following terms are termed the geometrical stress and the material contributions to the tangent matrix: 522@f{align*} 523& \int_{\Omega_0} \textrm{grad}\ \delta \mathbf{u} : 524 \textrm{grad}\ d \mathbf{u} [\boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}]~\textrm{d}V 525 && \quad {[\textrm{Geometrical stress}]} \, , 526 \\ 527& \int_{\Omega_0} \textrm{grad} \delta \mathbf{u} : 528 [J\mathfrak{c}_{\textrm{vol}} + J\mathfrak{c}_{\textrm{iso}}] :\textrm{grad}\ d \mathbf{u} 529 ~\textrm{d}V 530 && \quad {[\textrm{Material}]} \, . 531@f} 532 533 534<h3> Discretization of governing equations </h3> 535 536The three-field formulation used here is effective for quasi-incompressible materials, 537that is where $\nu \rightarrow 0.5$ (where $\nu$ is <a 538href="http://en.wikipedia.org/wiki/Poisson's_ratio">Poisson's ratio</a>), subject to a good choice of the interpolation fields 539for $\mathbf{u},~\widetilde{p}$ and $\widetilde{J}$. 540Typically a choice of $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ is made. 541Here $DGPM$ is the FE_DGPMonomial class. 542A popular choice is $Q_1 \times DGPM_0 \times DGPM_0$ which is known as the mean dilatation method (see Hughes (2000) for an intuitive discussion). 543This code can accommodate a $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ formulation. 544The discontinuous approximation 545allows $\widetilde{p}$ and $\widetilde{J}$ to be condensed out 546and a classical displacement based method is recovered. 547 548For fully-incompressible materials $\nu = 0.5$ and the three-field formulation will still exhibit 549locking behaviour. 550This can be overcome by introducing an additional constraint into the free energy of the form 551$\int_{\Omega_0} \Lambda [ \widetilde{J} - 1]~\textrm{d}V$. 552Here $\Lambda$ is a Lagrange multiplier to enforce the isochoric constraint. 553For further details see Miehe (1994). 554 555The linearised problem can be written as 556@f[ 557 \mathbf{\mathsf{K}}( \mathbf{\Xi}_{\textrm{i}}) d\mathbf{\Xi} 558 = 559 \mathbf{ \mathsf{F}}(\mathbf{\Xi}_{\textrm{i}}) 560@f] 561where 562@f{align*} 563 \underbrace{\begin{bmatrix} 564 \mathbf{\mathsf{K}}_{uu} & \mathbf{\mathsf{K}}_{u\widetilde{p}} & \mathbf{0} 565 \\ 566 \mathbf{\mathsf{K}}_{\widetilde{p}u} & \mathbf{0} & \mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}} 567 \\ 568 \mathbf{0} & \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}} & \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}} 569 \end{bmatrix}}_{\mathbf{\mathsf{K}}(\mathbf{\Xi}_{\textrm{i}})} 570 \underbrace{\begin{bmatrix} 571 d \mathbf{\mathsf{u}}\\ 572 d \widetilde{\mathbf{\mathsf{p}}} \\ 573 d \widetilde{\mathbf{\mathsf{J}}} 574 \end{bmatrix}}_{d \mathbf{\Xi}} 575 = 576 \underbrace{\begin{bmatrix} 577 -\mathbf{\mathsf{R}}_{u}(\mathbf{u}_{\textrm{i}}) \\ 578 -\mathbf{\mathsf{R}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}}) \\ 579 -\mathbf{\mathsf{R}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}}) 580 \end{bmatrix}}_{ -\mathbf{\mathsf{R}}(\mathbf{\Xi}_{\textrm{i}}) } 581= 582 \underbrace{\begin{bmatrix} 583 \mathbf{\mathsf{F}}_{u}(\mathbf{u}_{\textrm{i}}) \\ 584 \mathbf{\mathsf{F}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}}) \\ 585 \mathbf{\mathsf{F}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}}) 586 \end{bmatrix}}_{ \mathbf{\mathsf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, . 587@f} 588 589There are no derivatives of the pressure and dilatation (primary) variables present in the formulation. 590Thus the discontinuous finite element interpolation of the pressure and dilatation yields a block 591diagonal matrix for 592$\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}$, 593$\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}$ and 594$\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}}$. 595Therefore we can easily express the fields $\widetilde{p}$ and $\widetilde{J}$ on each cell simply 596by inverting a local matrix and multiplying it by the local right hand 597side. We can then insert the result into the remaining equations and recover 598a classical displacement-based method. 599In order to condense out the pressure and dilatation contributions at the element level we need the following results: 600@f{align*} 601 d \widetilde{\mathbf{\mathsf{p}}} 602 & = \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \bigl[ 603 \mathbf{\mathsf{F}}_{\widetilde{J}} 604 - \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}} d \widetilde{\mathbf{\mathsf{J}}} \bigr] 605 \\ 606 d \widetilde{\mathbf{\mathsf{J}}} 607 & = \mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}^{-1} \bigl[ 608 \mathbf{\mathsf{F}}_{\widetilde{p}} 609 - \mathbf{\mathsf{K}}_{\widetilde{p}u} d \mathbf{\mathsf{u}} 610 \bigr] 611 \\ 612 \Rightarrow d \widetilde{\mathbf{\mathsf{p}}} 613 &= \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \mathbf{\mathsf{F}}_{\widetilde{J}} 614 - \underbrace{\bigl[\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}} 615 \mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}^{-1}\bigr]}_{\overline{\mathbf{\mathsf{K}}}}\bigl[ \mathbf{\mathsf{F}}_{\widetilde{p}} 616 - \mathbf{\mathsf{K}}_{\widetilde{p}u} d \mathbf{\mathsf{u}} \bigr] 617@f} 618and thus 619@f[ 620 \underbrace{\bigl[ \mathbf{\mathsf{K}}_{uu} + \overline{\overline{\mathbf{\mathsf{K}}}}~ \bigr] 621 }_{\mathbf{\mathsf{K}}_{\textrm{con}}} d \mathbf{\mathsf{u}} 622 = 623 \underbrace{ 624 \Bigl[ 625 \mathbf{\mathsf{F}}_{u} 626 - \mathbf{\mathsf{K}}_{u\widetilde{p}} \bigl[ \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \mathbf{\mathsf{F}}_{\widetilde{J}} 627 - \overline{\mathbf{\mathsf{K}}}\mathbf{\mathsf{F}}_{\widetilde{p}} \bigr] 628 \Bigr]}_{\mathbf{\mathsf{F}}_{\textrm{con}}} 629@f] 630where 631@f[ 632 \overline{\overline{\mathbf{\mathsf{K}}}} \dealcoloneq 633 \mathbf{\mathsf{K}}_{u\widetilde{p}} \overline{\mathbf{\mathsf{K}}} \mathbf{\mathsf{K}}_{\widetilde{p}u} \, . 634@f] 635Note that due to the choice of $\widetilde{p}$ and $\widetilde{J}$ as discontinuous at the element level, all matrices that need to be inverted are defined at the element level. 636 637The procedure to construct the various contributions is as follows: 638- Construct $\mathbf{\mathsf{K}}$. 639- Form $\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}^{-1}$ for element and store where $\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}$ was stored in $\mathbf{\mathsf{K}}$. 640- Form $\overline{\overline{\mathbf{\mathsf{K}}}}$ and add to $\mathbf{\mathsf{K}}_{uu}$ to get $\mathbf{\mathsf{K}}_{\textrm{con}}$ 641- The modified system matrix is called ${\mathbf{\mathsf{K}}}_{\textrm{store}}$. 642 That is 643 @f[ 644 \mathbf{\mathsf{K}}_{\textrm{store}} 645\dealcoloneq 646 \begin{bmatrix} 647 \mathbf{\mathsf{K}}_{\textrm{con}} & \mathbf{\mathsf{K}}_{u\widetilde{p}} & \mathbf{0} 648 \\ 649 \mathbf{\mathsf{K}}_{\widetilde{p}u} & \mathbf{0} & \mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}^{-1} 650 \\ 651 \mathbf{0} & \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}} & \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}} 652 \end{bmatrix} \, . 653 @f] 654 655 656<h3> The material class </h3> 657 658A good object-oriented design of a Material class would facilitate the extension of this tutorial to a wide range of material types. 659In this tutorial we simply have one Material class named Material_Compressible_Neo_Hook_Three_Field. 660Ideally this class would derive from a class HyperelasticMaterial which would derive from the base class Material. 661The three-field nature of the formulation used here also complicates the matter. 662 663The Helmholtz free energy function for the three field formulation is $\Psi = \Psi_\text{vol}(\widetilde{J}) + \Psi_\text{iso}(\overline{\mathbf{b}})$. 664The isochoric part of the Kirchhoff stress ${\boldsymbol{\tau}}_{\text{iso}}(\overline{\mathbf{b}})$ is identical to that obtained using a one-field formulation for a hyperelastic material. 665However, the volumetric part of the free energy is now a function of the primary variable $\widetilde{J}$. 666Thus, for a three field formulation the constitutive response for the volumetric part of the Kirchhoff stress ${\boldsymbol{\tau}}_{\text{vol}}$ (and the tangent) is not given by the hyperelastic constitutive law as in a one-field formulation. 667One can label the term 668$\boldsymbol{\tau}_{\textrm{vol}} \equiv \widetilde{p} J \mathbf{I}$ 669as the volumetric Kirchhoff stress, but the pressure $\widetilde{p}$ is not derived from the free energy; it is a primary field. 670 671In order to have a flexible approach, it was decided that the Material_Compressible_Neo_Hook_Three_Field would still be able to calculate and return a volumetric Kirchhoff stress and tangent. 672In order to do this, we choose to store the interpolated primary fields $\widetilde{p}$ and $\widetilde{J}$ in the Material_Compressible_Neo_Hook_Three_Field class associated with the quadrature point. 673This decision should be revisited at a later stage when the tutorial is extended to account for other materials. 674 675 676<h3> Numerical example </h3> 677 678The numerical example considered here is a nearly-incompressible block under compression. 679This benchmark problem is taken from 680- S. Reese, P. Wriggers, B.D. Reddy (2000), 681 A new locking-free brick element technique for large deformation problems in elasticity, 682 <em> Computers and Structures </em>, 683 <strong> 75 </strong>, 684 291-304. 685 DOI: <a href="http://doi.org/10.1016/S0045-7949(99)00137-6">10.1016/S0045-7949(99)00137-6</a> 686 687 <img src="https://www.dealii.org/images/steps/developer/step-44.setup.png" alt=""> 688 689The material is quasi-incompressible neo-Hookean with <a href="http://en.wikipedia.org/wiki/Shear_modulus">shear modulus</a> $\mu = 80.194e6$ and $\nu = 0.4999$. 690For such a choice of material properties a conventional single-field $Q_1$ approach would lock. 691That is, the response would be overly stiff. 692The initial and final configurations are shown in the image above. 693Using symmetry, we solve for only one quarter of the geometry (i.e. a cube with dimension $0.001$). 694The inner-quarter of the upper surface of the domain is subject to a load of $p_0$. 695 696 697