Hyperelastic material

Yeoh strain energy density function:

$\Psi = \sum_{i=1}^{3} {c_i}\left(\bar{\mathbf{I}}_{1e} - 3\right)^{1}$

../_images/instable.png — Example of used material parameters and possibile instable combinations created interpolating the coefficient (see myLoading.m):

Hyperelasticity background

Warning

The following is partially taken from my Master thesis.

We consider a continuous body in the reference configuration $\Omega_0$ that could move in space toward the spatial configuration $\Omega$ . The map between reference position $\mathbf{X}$ and the spatial position $\mathbf{x}$ of the material particle is represented by the motion $\mathbf{x}(t)=\chi(\boldsymbol{X},t)$ . In finite deformations, no assumptions are made on the magnitude of displacement $\mathbf{u}=\mathbf{x}-\mathbf{X}$ , that as no small displacements are required.

At the time $t$ the motion in the neighborhood of a material point is represented by the deformation gradient:

$\mathbf{F}(\mathbf{X}, t)=\operatorname{Grad} \chi(\mathbf{X}, t)={\partial \chi(\mathbf{X}, t) \over \partial \mathbf{X}}.$

It is the primary measure of deformation and allows to relate the relative position of neighboring particles after and before the deformation in terms of their relative material position.

The local motion could be also decomposed into a rigid body motion and deformation using a polar decomposition of the deformation gradient as

$\mathbf{F}=\mathbf{R U}=\mathbf { vR }.$

$\mathbf{U}$ and $\mathbf{v}$ represent unique, positive definite, symmetric tensors called respectively as right (or material) stretch tensor and left (or spatial) stretch tensor. They measure local stretching along mutually orthogonal eigenvectors. $\mathbf{R}$ is an orthogonal tensor called as rotation tensor representing the rigid body rotation.

Strain

At finite deformations, non-linear strain measures are used to correctly represent the differences between the reference and the spatial configuration.

The right Cauchy-Green strain is formulated with respect to the reference configuration and could be expressed as:

$\mathbf{C}=\mathbf{U}^2=\mathbf{F}^{\mathrm{T}} \mathbf{F}.$

An additional deformation tensor is the Green-Lagrange strain tensor, reading:

$\mathbf{E}=\frac{1}{2}\left(\mathbf{F}^{\mathrm{T}} \mathbf{F}-\mathbf{I}\right).$

The latter is used to represent the change in the squared lengths $\lambda_a^2 -1$ , where $\lambda_a$ is the linear stretch and represents how the line element is deformed from the original direction $\mathbf{a}_0$ to the spatial configuration $\mathbf{a}=\mathbf{F} \mathbf{a}_0$ (with $|\mathbf{a}_0|=1$ ). This is:

$\lambda_{\mathbf{a}}=\frac{\mathrm{d} l}{\mathrm{~d} L}=\frac{\left|\mathbf{F} \mathbf{a}_0\right|}{\left|\mathbf{a}_0\right|}=\sqrt{\mathbf{a}_0 \mathbf{F}^{\mathrm{T}} \mathbf{F} \mathbf{a}_0}=|\mathbf{a}|.$

A line element is extended for $\lambda_a>1$ , unstretched for $\lambda_a=1$ or compressed for $\lambda_a<1$ .

The Cauchy-Green strain tensor can be also represented through eigenvalue representation, with the principle stretches $\lambda_i$ and the principal stretch direction $\widehat{\mathbf{N}}_i$ , as:

$\mathbf{C}=\lambda_i^2 \widehat{\mathbf{N}}_i \otimes \widehat{\mathbf{N}}_i.$

Stress

To represent the load on the material particle a stress measure is required. The Cauchy stress $\boldsymbol{{\sigma}}$ is a second-order tensor used to represent the stress state containing both the normal stress and the shear stress components acting respectively perpendicular or in parallel to the local area $ds$ as shown in cref{fig:cauchy_stress}.

However, at finite deformation, where the current configuration is different from the reference one, other stress tensors may be defined. While the Cauchy stress is represented in the current configuration, the first Piola-Kirchoff relates the force in the spatial configuration to the area element in the reference configuration. Through the backward Piola transform it could be related to the Cauchy stress as:

$\mathbf{P}=J \boldsymbol{\sigma} \mathbf{F}^{-\mathrm{T}},$

where $J=\operatorname{det}(\mathbf{F})$ and $\mathbf{F}$ is the deformation gradient.

The second Piola-Kirchhoff stress pull-back the force in the spatial configuration to the reference one and so it is entirely expressed in the reference configuration $\Omega_0$ . It could be expressed as:

$\mathbf{S}=\mathbf{F}^{-1} \mathbf{P}.$

Furthermore, the Cauchy stress can be represented in a coordinate system within which the shear stress disappears. This stress state is called the principal stress state and the principal stress directions $\widehat{\mathbf{n}}_i$ represent the eigenvectors of the stress state. Introducing the principal stresses $\sigma_i$ , the Cauchy stress in a principal coordinate system reads:

$\boldsymbol{\sigma}=\sigma_i \widehat{\mathbf{n}}_i \otimes \widehat{\mathbf{n}}_i.$

Vascular vessels are cylinder-like structures that can be approximated as a typical engineering problem named as pressurized tube. Inflated vessels may be characterized by a mutual tri-orthogonal axial system representing cylindrical coordinates. For this simplified setup, the three principal directions are circumferential, axial and radial one with the related stresses shown in cref{fig:cylindricalStress}. Due to the boundary conditions of the problem the radial stress has to balance the internal pressure inside and the external one outside while the circumferential stresses are generated to resist the resulting bursting effect. Since the internal pressure is higher than the external one, radial stresses are generally higher inside and gradually decrease across the thickness. Circumferential stresses follow a similar trend depending on the internal pressure but also on wall thickness. For a thick-wall cylinder, they could be calculated using the Lam'e equations with the resulting trend in cref{fig:cylindricalStress}.

Governing laws

The motion of a material particle of density $\rho$ is governed by the linear momentum balance equation:

$\rho \frac{\mathrm{D} \mathbf{v}}{\mathrm{D} t}=\operatorname{div} \boldsymbol{\sigma}+\mathbf{b}_{\mathrm{f}},$

where $\boldsymbol{\sigma}$ is the Cauchy stress tensor, $\mathbf{v}$ is the particle velocity and $\mathbf{b}_{\mathrm{f}}$ is the body force per unit volume.

The balance of angular momentum is a priori satisfied by the symmetry of the Cauchy stress tensor: $\boldsymbol{\sigma}=\boldsymbol{\sigma}^{\mathrm{T}}$ .

Boundary Value Problem

In continuum mechanics, the primary unknown variable is the displacement field $\mathbf{u}$ and this can be obtained from the solution of the Boundary Value Problem (BVP) derived from the balance of linear momentum and constitutive relations.

The body is characterized by the boundary $\partial \Omega$ , split into $\partial \Omega_{\mathrm{u}}$ and $\partial \Omega_{\mathrm{t}}$ . At $\partial \Omega_{\mathrm{u}}$ Dirichlet or essential boundary conditions are prescribed in order to impose the displacement, while at $\partial \Omega_{\mathrm{t}}$ Neumann or natural boundary conditions prescribe the traction through the boundary normal.

For a transient problem, the mathematical description should be closed by the initial conditions on the velocity and displacement field. hence, the BVP is described by the following set of equations:

$\left\{ \begin{array}{rlrl} \operatorname{div} \boldsymbol{\sigma}+\mathbf{b}_{\mathrm{f}} & =\rho {D^2{\mathbf{u}}\over dt^2} && \text {in } \Omega, \\ \mathbf{u} & =\overline{\mathbf{u}} & & \text {at } \partial \Omega_{\mathrm{u}}, \\ \mathbf{t}=\sigma \mathbf{n} & =\overline{\mathbf{t}} & & \text {at } \partial \Omega_{\mathrm{t}}, \\ \frac{\mathrm{D} \mathbf{u}}{\mathrm{D} t}=\mathbf{v}_0 ; \mathbf{u} & =\mathbf{u}_0 & & \text {in } \Omega \text {, at } t=t_0 \end{array}\right. \label{eq:BVP}$

Finite Element Method

The Finite Element Method is a widely used method to numerically solve partial differential equations (PDE). It is a procedure based on the representation of the continuum with a finite number of points, called nodes, within the continuum body in order to discretize the full body (of arbitrary geometry) by means of a finite set of sub-regions of standard geometries called finite elements. Over every region, both the geometry and the variables of the governing equations are approximated in order to solve complex systems of equations over elaborate domains.

The governing equations are usually reformulated in the weak integral form and the finite element approximation leads to a finite set of non-linear algebraic equations in the primary variable. Then, these equations could be solved with an iterative technique such as the Netwton-Rapson method.

Spatial description

The domain is approximated as a discretized domain segmented into $n$ subdomains:

$\Omega \simeq \Omega_h=\sum_{e=1}^{n} \Omega_e.$

This discretization allows us to focus only on a single finite element $\Omega_e$ for the governing equations and material description.

The primary variables are approximated through the shape functions $N_i$ . The latter interpolate the primary variables in space over the finite element $\Omega_e$ through their values at nodes as:

$\mathbf{u}(\xi)=\sum_{i=1}^{n_{ {i}}} N_i(\boldsymbol{\xi}) \mathbf{u}_i,$

where $\xi$ represents the coordinates in the single element domain with in $n_{i}$ nodes.

By storing all the independent variables $\mathbf{u}_i \:( i=1, \ldots, n_{\mathrm{n}_{dim}})$ into a vector $\mathbf{h}$ of dimension $\mathrm{n}_{\mathrm{dof}}=\mathrm{n}_{\mathrm{i}} \:\mathrm{n}_{\mathrm{dim}}$ the finite element interpolation can be reformulated as

$\mathbf{u}=\mathbf{N h}, \label{eq:interpolation}$

with the interpolation matrix $\mathbf{N}$ collecting shape functions $N_i$ .

By using an isoparametric finite element formulation, both the essential variables and the spatial coordinates $\mathbf{x}$ use the same shape function, so:

$\mathbf{x}(\xi)= \sum_{i=1}^{n_{\text {n }}} N_i(\boldsymbol{\xi}) \mathbf{x}_i.$

The shape functions could be defined in the natural coordinate of the finite element $\boldsymbol{\xi}$ and linked to the global coordinate $\mathbf{x}$ through the Jacobian $\mathbf{J}=\partial \mathbf{x} / \partial \xi$ . It allows us to define shape function and further numerical quadrature only for each specific typology of finite element, and then map it to the current geometries.

Due to the interpolation in cref{eq:interpolation} also the symmetric spatial gradient of the essential variables could be expressed through an interpolation matrix as:

$\operatorname{grad}_{\mathrm{s}} \mathbf{u}=\mathbf{B h},$

where $\mathbf{h}$ expresses the $n_s$ independent components of the gradients as a single vector with the convention of $\left[(\bullet)_{11}(\bullet)_{22}(\bullet)_{33}(\bullet)_{23}(\bullet)_{12}(\bullet)_{12}\right]^{\mathrm{T}}$ . $\mathbf{B}$ represents the gradient interpolation matrix and stores the gradient of the shape functions. This gradient requires the Jacobian matrix because the shape functions are expressed with respect to the natural coordinates $\boldsymbol{\xi}$ but they are derived with respect to the global coordinate $\mathbf{x}$ . Each row of $\mathbf{B}$ is expressed as:

$\frac{\partial N_i(\boldsymbol{\xi})}{\partial \mathbf{x}}=\frac{\partial N_i(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}}: \frac{\partial \boldsymbol{\xi}}{\partial \mathbf{x}}=\frac{\partial N_i(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}}: \mathbf{J}^{-1}.$

Calculus of variations

Physical problems are often described by local governing equations with the description in the so-called strong form of the problem. However, FEM requires the weak or integral form in order to discretize and then solve the problem. A general method to derive the weak form is the Calculus of Variations. Multiplying the governing equations with a test function, and then integrating over the body domain, leads to the weak form of the problem.

For a solid mechanics problem, the body occupies the initial configuration $\Omega$ and the problem is formulated with respect to the essential variable, the displacement $\mathbf{u}$ . The state of motion is described by the displacement and its time derivatives, velocity $\dot{\mathbf{u}}$ and acceleration $\ddot{\mathbf{u}}$ . They have to satisfy the equilibrium equation. Furthermore, a Dirichlet boundary $\Omega_u$ allows for prescribed displacement condition $\mathbf{u}=\overline{\mathbf{u}}$ and a Neumann boundary $\partial \Omega_t$ for traction condition oriented across the outward-normal $\mathbf{t}=\sigma \mathbf{n}=\overline{\mathbf{t}}$ . This leads to the BVP as expressed in cref{eq:BVP}.

The equilibrium equation is then multiplied with an admissible variation $\delta \mathbf{u}$ and it is integrated over the body domain.

$\int_{\Omega} \delta \mathbf{u} \cdot\left(\operatorname{div} \boldsymbol{\sigma}+\mathbf{b}_{\mathrm{f}}-\rho \ddot{\mathbf{u}}\right) \mathrm{d} v=0.$

The first term can be split following the integration by parts:

$\int_{\Omega} \delta \mathbf{u} \cdot \operatorname{div} \boldsymbol{\sigma} \mathrm{d} v=\int_{\Omega} \operatorname{div}(\boldsymbol{\sigma} \delta \mathbf{u}) \mathrm{d} v-\int_{\Omega}(\boldsymbol{\sigma}: \operatorname{grad} \delta \mathbf{u}) \mathrm{d} v.$

Furthermore, the divergence theorem can be applied:

$\int_{\Omega} \operatorname{div}(\boldsymbol{\sigma} \delta \mathbf{u}) \mathrm{d} v=\int_{\partial \Omega} \delta \mathbf{u} \cdot \boldsymbol{\sigma} \mathbf{n} \mathrm{d} s.$

Then, the weak form can be expressed as:

$\underbrace{\int_{\partial \Omega_{\mathrm{t}}} \overline{\mathbf{t}} \cdot \delta \mathbf{u d} s+\int_{\Omega} \delta \mathbf{u} \cdot\left(\mathbf{b}_{\mathrm{f}}-\rho \ddot{\mathbf{u}}\right) \mathrm{d} v}_{\delta \Pi_{\mathrm{ext}}}-\underbrace{\int_{\Omega} \sigma: \operatorname{grad}_{\mathrm{s}} \delta \mathbf{u} \mathrm{d} v}_{\delta \Pi_{\text {int }}}=0, \label{eq:weakform}$

where the two $\delta \Pi$ terms represent respectively the external and internal work upon the admissible variation.

It is also interesting to observe that cref{eq:weakform} correspond to the principle of virtual work when the admissible variation $\delta \mathbf{u}$ is regarded as corresponding to an admissible variation of the engineering strain $\delta \boldsymbol{\varepsilon}=\operatorname{grad}_{\mathrm{s}} \delta \mathbf{u}$ .

The weak form is irreducible and involves only the first-order derivative of the essential variable. So, the derivative order of the displacement field appears to be reduced with respect to the strong from and shared with the derivative of the test function.

In non-linear mechanical problems, where large deformations are involved together with non-linear constitutive models, it is convenient to express the internal work term with respect to the reference configuration through the second Piola-Kirchhoff stress and the variation of the Green-Lagrange strain, as:

$\label{eq:deltaPi} \delta \Pi_{\mathrm{int}}=\int_{\Omega} \boldsymbol{\sigma}:\delta \boldsymbol{\varepsilon}\mathrm{d} v=\int_{\Omega_0} \mathbf{S} :\delta \mathbf{E} \mathrm{d} V,$

where $\delta \mathbf{E}$ is the admissible variation of the Green-Lagrange strain.

The balance equation obtained from Eq. eqref{eq:deltaPi} is nonlinear and, to be solved, has to be linearized along the increment of the displacement, that is:

$D_{\mathbf{u}} \delta \Pi_{\text {int}} = \delta \Pi_{\text {int}}(\mathbf{u}+\Delta \mathbf{u})- \delta \Pi_{\text {int}}(\mathbf{u})$

With the direction derivative $D_{\mathbf{u}}$ along the increment of displacement $\Delta \mathbf{u}$ :

$D_{\mathbf{u}} \delta \Pi_{\text {int }}=D_{\mathbf{u}}\left(\int_{\Omega_0} \mathbf{S}: \delta \mathbf{E} \mathrm{d} V\right)=\int_{\Omega_0}\left(\mathbf{S}: D_{\mathbf{u}} \delta \mathbf{E}+\delta \mathbf{E}: \mathbb{C}: D_{\mathbf{u}} \mathbf{E}\right) \mathrm{d} V$

where $\mathbb{C}(\mathbf{E})=\frac{\partial \mathbf{S}}{\partial \mathbf{E}}$ is the non-linear material stiffness.

By pushing it forward to the spatial configuration:

$D_{\mathbf{u}} \delta \Pi_{\mathrm{int}}=\int_{\Omega}\left(\operatorname{grad}_{\mathrm{s}} \delta \mathbf{u}: \operatorname{grad}_{\mathrm{s}} \Delta \mathbf{u} \sigma+\operatorname{grad}_{\mathrm{s}} \delta \mathbf{u}: \mathbb{C}: \operatorname{grad}_{\mathrm{s}} \Delta \mathbf{u}\right) \mathrm{d} v \label{eq:variation}$

There are two different contributions to the linearized internal work, that is a geometric contribution:

$D_{\mathbf{u}} \delta \Pi_{\text {int geo }}=\int_{\Omega}\left(\operatorname{grad}_{\mathrm{s}} \delta \mathbf{u}: \operatorname{grad}_{\mathrm{s}} \Delta \mathbf{u} \sigma\right) \mathrm{d} v,$

and a material contribution:

$D_{\mathbf{u}} \delta \Pi_{\text {int mat }}=\int_{\Omega}\left(\operatorname{grad}_{\mathrm{s}} \delta \mathbf{u}: \mathbb{C}: \operatorname{grad}_{\mathrm{s}} \Delta \mathbf{u}\right) \mathrm{d} v.$

Both terms depend on the choice of the constitutive model to be determined.

Discretization

Galerkin’s approach uses the same interpolation for essential variable and their test function. In other words, the nodal displacements, stored in $\mathbf{h}$ , and their admissible variations, stored in $\delta \mathbf{h}$ , are interpolated with the same interpolation functions, collected in $\mathbf{N}$ . Then the discretized weak form for the $e-th$ finite element reads:

$\delta h_i \left[\underbrace{\int_{\Omega_{\mathrm{e}}} N_{a i} \rho N_{a j} \mathrm{~d} v}_{M_{i j}} \ddot{h}_j+{K_{i j}} h_j -\underbrace{\left(\int_{\Omega_{\mathrm{e}}} b_{\mathrm{f} a} N_{a i} \mathrm{~d} v+\int_{\partial \Omega_{\mathrm{et}}} \bar{t}_a N_{a i} \mathrm{~d} s\right)}_{f_i}\right]=0.$

The $K_{ij}$ term contains two different contributions. The material contribution leads to the material stiffness matrix of the finite element:

$K_{\text {mat } i j}=\int_{\Omega_e} B_{a i} C_{a b} B_{b j} \mathrm{~d} v,$

where $C_{ab}$ denotes the material stiffness linked to the constitutive model.

The geometrical contribution denotes the geometric stiffness matrix:

$K_{\mathrm{geo}\: i j}=\int_{\Omega_e} \frac{\partial N_{a i}}{\partial x_c} \sigma_{a b} \frac{\partial N_{c j}}{\partial x_b} \mathrm{~d} v,$

where $\sigma_{ab}$ are the component of the Cauchy stress.

These integrals are commonly computed through numerical quadrature formulas.

Given $n_n$ arbitrary admissible variation $\delta \mathbf{h}$ , the discretized weak form yields to a system of algebraic equations:

$\mathbf{M} \ddot{\mathbf{h}}+\mathbf{K h}-\mathbf{f}=\mathbf{0}$

This system of equations can be solved with numerical iterative methods. If the velocity and acceleration terms are different from zero also a time discretization is required. \

A non-linear stationary FEM problem can be expressed as:

$\mathbf{K}(\mathbf{h}) \Delta \mathbf{h}-\mathbf{f}=\mathbf{r}(\mathbf{h})$

and iteratively solved through a numerical procedures that minimize the residuum $\mathbf{r}(\mathbf{h})$ , ideally until $|\mathbf{r}|\to 0$ . Adopting a Newton-Rapson iterative scheme, the current solution at step $i+1$ is computed from the previous step $i$ as:

$\mathbf{h}_{i+1}=\mathbf{h}_i+\left.\Delta \mathbf{h}\right|_{\mathbf{h}=\mathbf{h}_i}\:,$

with the increment $\Delta \mathbf{h}$ :

$\Delta \mathbf{h}=-\left(\partial \mathbf{r} \over \partial \mathbf{h}\right)^{-1} \mathbf{r}=-\left( \mathbf{K}(\mathbf{h}) \right)^{-1} \mathbf{r}(\mathbf{h}).$

Both the material and the geometrical stiffness matrixes are functions of the essential variable and need to be evaluated in each step of the iterative solver. Furthermore, also a load increment approach could be used. The external boundary load is incremented with a load factor and an equivalent system is solved for each load step, as two nested loops. newline

Further technical details for implementation with the Comsol finite element solver are described in the appendix.

Kinematic decomposition

../_images/decomposition.png — Representation of the tissue growth and remodeling through a kinematic description. The stress-free intermediate (eventually incompatible) configuration Ω0 separates the growth-related and the non-growth-related motions. The deformation gradient is multiplicatively decomposed and the final deformation is due to both the deformation gradient and the growth tensor.

The total deformation gradient reads:

$\mathbf{F}=\mathbf{F}_{\mathrm{e}} \mathbf{F}_{\mathrm{g}},$

where $\mathbf{F}^{\mathrm{e}}$ is the elastic deformation gradient and $\mathbf{F}^{\mathrm{g}}$ is the growth tensor.

The elastic deformation gradient can be exprssed as:

$\mathbf{F}_{\mathrm{e}} =\mathbf{F} \mathbf{F}_{\mathrm{g}}^{-1}$