Gradient-Enhanced Continuum Damage Model for Residually Stressed Fibre-Reinforced Materials at Finite Strains

–$$\,d]$$-damage formulation, where only the anisotropic elastic part is assumed to damage. The resulting coupled, highly non-linear system of equations is symmetric and can conveniently be solved by standard incremental-iterative Newton-Raphson schemes or arc-length-based solution methods without any need for advanced and computationally expensive solution methods such as a global active-set-search as applied by Liebe et al. [11]. Furthermore, the approach proves to be robust, even at high levels of degradation, and allows to incorporate any suitable scalar-valued damage formulation.


Apart from the highly nonlinear elastic response of soft biological tissues, it is well-known that the structural design of arteries is characterised by a fibre-reinforced multi-layered composite subjected to pronounced residual stresses. The complex interaction of material properties with these residual stress effects and the geometry guarantees the optimal support under different blood pressures within the vessel. As a further key aspect of this contribution, we therefore incorporate residual stresses by means of a multiplicative composition of the deformation gradient.

This article is structured as follows: In Sect. 2, we summarise relevant kinematic relations for the geometrically non-linear case and the balance equations of the coupled boundary value problem in weak form. In Sect. 3, we specify the underlying constitutive equations, containing the isotropic and anisotropic non-linear elastic and gradient-enhanced free energies, as well as the continuum damage formulation. In Sect. 4, we discretise the governing weak forms by means of the finite element method resulting in a coupled non-linear system of equations. Last, in Sect. 6, we apply the model to illustrative three-dimensional inhomogeneous deformation problems. In order to show the capabilities of the approach with regard to biomechanics-related problems, we study a force-driven finite element example by means of an anisotropic artery-like tube subjected to internal pressure and under consideration of residual stresses. We conclude with a summary and future perspectives in Sect. 7.



2 Gradient Enhancement of a Continuum Damage Formulation


This section summarises the essential kinematic relations and presents the governing coupled balance equations of the boundary value problem on the basis of the principle of minimum total potential energy. The related variational form provides the basis for the finite element discretisation described in Sect. 4.


2.1 Basic Kinematics


Let $${\varvec{x}} = {\varvec{\varphi }}({\varvec{X,t}}) $$ describe the deformation of the body, which transforms referential placements $${\varvec{X}} \in \mathcal {B}_0$$ to their spatial counterparts $${\varvec{x}}\in \mathcal {B}_t$$. Based on this, the deformation gradient is defined as $${\varvec{F}}=\nabla _{\!{\varvec{X}}} {\varvec{\varphi }}$$ which transforms infinitesimal referential line elements $$\mathrm {d}{\varvec{X}}$$ onto their spatial counterparts $$\mathrm {d}{\varvec{x}}$$. Furthermore, let $$\mathrm {d}V$$ and $$\mathrm {d}v$$ denote an infinitesimal volume element in referential and spatial setting. Accordingly, the Jacobian $$ J=\mathrm {d}v/\mathrm {d}V=\det ({\varvec{F}})>0 $$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq11.gif”></SPAN> is the ratio of the deformed to the undeformed volume. Finally, let <SPAN id=IEq12 class=InlineEquation><IMG alt= and $$\mathrm {d}{\varvec{a}}={\varvec{n}}\,\mathrm {d}a$$ define the referential and spatial area normals. Then, Nanson’s formula $$ {\varvec{n}} \,\mathrm {d}a=\mathrm {cof}({\varvec{F}})\cdot {\varvec{N}}\,\mathrm {d}A $$ describes the transformation of infinitesimal area elements between the reference and the spatial configuration with the co-factor of $${\varvec{F}}$$ defined as $$\mathrm {cof}({\varvec{F}})=J{\varvec{F}}^{\,-\mathrm {t}}$$. Fibre-reinforcement of the material is incorporated by assuming two families of fibres to be embedded in the continuum. Their orientation is characterised by referential unit vectors $${\varvec{a}}_{0\,i}$$, $$i=1,2$$ with $$||{\varvec{a}}_{0\,i}||=1$$.


2.2 General Gradient-Enhanced Format of the Free Energy


We assume the local free energy


$$\begin{aligned} \varPsi _\mathrm {loc}({\varvec{F}},\varvec{a}_{0\,i},\kappa )=\varPsi _\mathrm {iso}({\varvec{F}}\,)+f_\mathrm {d}(\kappa )\,\varPsi _\mathrm {ani}({\varvec{F}},\varvec{a}_{0\,i}) \end{aligned}$$

(1)
to account for anisotropic non-linear elastic material response under the influence of scalar damage. Basically, we additively compose the effective free energy of the undamaged material of an isotropic contribution $$\varPsi _\mathrm {iso}$$ representing the ground substance and of an anisotropic contribution $$\varPsi _\mathrm {ani}$$ associated with $$i$$ fibre families. We assume only the anisotropic part to be affected by the damage, whereas the isotropic matrix material remains elastic. In Eq. (1), $$\kappa \in [0,\infty )$$ is a scalar internal damage variable, characterising a material stiffness loss of the fibres, while $$f_\mathrm {d}(\kappa )=1-d\in (0,1]$$ represents an appropriate damage function that is at least twice differentiable and satisfies $$f_\mathrm {d}(0)=1$$ and $$f_\mathrm {d}(\kappa \rightarrow \infty )\rightarrow 0$$. This ensures purely elastic behaviour of the undamaged material and complete loss of the related material stiffness for $$\kappa \rightarrow \infty $$. Conceptually following the approach by Dimitrijević and Hackl [2], we introduce a gradient-enhanced non-local free energy $$\varPsi _\mathrm {nloc}$$ as


$$\begin{aligned} \varPsi _\mathrm {nloc}(\phi ,\nabla _{\!{\varvec{X}}}\phi ,\kappa ;{\varvec{F}})=\varPsi _\mathrm {grd}(\nabla _{\!{\varvec{X}}}\phi ;{\varvec{F}})+\varPsi _\mathrm {plty}(\phi ,\kappa )\, . \end{aligned}$$

(2)
Here, $$\varPsi _\mathrm {grd}$$ contains the referential gradient of the non-local damage field variable $$\phi $$ while $$\varPsi _\mathrm {plty}$$ incorporates a penalisation term which links the non-local damage variable $$\phi $$ to the local damage variable $$\kappa $$. Consequently, we obtain an enhanced free energy as


$$\begin{aligned} \varPsi _\mathrm {int}({\varvec{F}},\phi ,\nabla _{\!{\varvec{X}}}\phi ,{\varvec{a}}_{0\,i},\kappa )=\varPsi _\mathrm {loc}({\varvec{F}},\varvec{a}_{0\,i},\kappa )+\varPsi _\mathrm {nloc}(\phi ,\nabla _{\!{\varvec{X}}}\phi ,\kappa ;{\varvec{F}})\, . \end{aligned}$$

(3)
Provided that the external load can be derived from a potential, we can specify the local external energy function as $$\varPsi _\mathrm {ext}({\varvec{\varphi }})$$. In summary, the total potential energy function is additively composed of the internal and external contribution so  that its local form reads


$$\begin{aligned} \varPsi ({\varvec{\varphi }},{\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ,{\varvec{a}}_{0\,i},\kappa )=\varPsi _\mathrm {int}({\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ,{\varvec{a}}_{0\,i},\kappa )+\varPsi _\mathrm {ext}({\varvec{\varphi }})\, . \end{aligned}$$

(4)


2.3 Total Potential Energy


The total potential energy of a system additively combines the internal contribution $$\varPi _\mathrm {int}$$, reflecting the action of internal forces, and an external contribution $$\varPi _\mathrm {ext}= \varPi _\mathrm {vol}+ \varPi _\mathrm {sur}$$ due to volume and surface forces, i.e.


$$\begin{aligned} \varPi ({\varvec{\varphi }},{\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ;{\varvec{a}}_{0\,i},\kappa )=\varPi _\mathrm {int}+\varPi _\mathrm {ext}\, . \end{aligned}$$

(5)
The internal energy contribution can be written as


$$\begin{aligned} \varPi _\mathrm {int}({\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ;{\varvec{a}}_{0\,i},\kappa )=\int \limits _{\mathcal {B}_0}\varPsi _\mathrm {int}\,\mathrm {d}V\, . \end{aligned}$$

(6)
while the external contributions, assuming ‘dead’ loads, are provided by


$$\begin{aligned} \varPi _\mathrm {vol}({\varvec{\varphi }})&=\int \limits _{\mathcal {B}_0}\varPsi _\mathrm {vol}\,\mathrm {d}V=-\int \limits _{\mathcal {B}_0}\bar{{\varvec{B}}}\cdot {\varvec{\varphi }}\,\mathrm {d}V\, ,\end{aligned}$$

(7)



$$\begin{aligned} \varPi _\mathrm {sur}({\varvec{\varphi }})&=\int \limits _{\partial \mathcal {B}_0}\varPsi _\mathrm {sur}\,\mathrm {d}A=-\int \limits _{\partial \mathcal {B}_0}\bar{{\varvec{T}}}\cdot {\varvec{\varphi }}\,\mathrm {d}A\, , \end{aligned}$$

(8)
where $$\bar{{\varvec{B}}}$$ denotes the body force vector per unit reference volume and $$\bar{{\varvec{T}}}$$ characterises the traction vector per unit reference surface area. In this regard, see, for instance, Waffenschmidt and Menzel [15] where a double-layered thick-walled cylindrical tube subjected to internal pressure is analysed on the basis of a total potential.


2.4 Variational Form


The boundary value problem is governed by the principle of minimum potential energy


$$\begin{aligned} \min \limits _{\varphi ,\,\phi }{\varPi ({\varvec{\varphi }},{\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ;{\varvec{a}}_{0\,i},\kappa )}\, , \end{aligned}$$

(9)
which requires the first variation of the total potential energy with respect to $${\varvec{\varphi }}$$ and $$\phi $$ to vanish, i.e.


$$\begin{aligned} \delta _{\varvec{\varphi }}\varPi ({\varvec{\varphi }},{\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ;{\varvec{a}}_{0\,i},\kappa )&\doteq {\varvec{0}}\qquad \forall \quad \delta {\varvec{\varphi }}\, ,\end{aligned}$$

(10)



$$\begin{aligned} \delta _{\phi }\varPi ({\varvec{\varphi }},{\varvec{F}},\phi ,\nabla _{\!\varvec{X}}\phi ; {\varvec{a}}_{0\,i},\kappa )&\doteq 0\qquad \forall \quad \delta \phi \, . \end{aligned}$$

(11)
Taking into account that $$\delta _{\varvec{\varphi }}W_\mathrm {int}=\delta _{\varvec{\varphi }}\varPi _\mathrm {int}$$ and $$\delta _{\varvec{\varphi }}W_\mathrm {ext}=-\,\delta _{\varvec{\varphi }}\varPi _\mathrm {ext}$$, we obtain a coupled system of variational equations, i.e.


$$\begin{aligned} \delta _{\varvec{\varphi }}W=\delta _{\varvec{\varphi }}W_\mathrm {int}-\delta _{\varvec{\varphi }}W_\mathrm {ext}&={\varvec{0}} \qquad \forall \quad \delta {\varvec{\varphi }}\, ,\end{aligned}$$

(12)



$$\begin{aligned} \delta _{\phi }W=\delta _{\phi }W_\mathrm {int}-\delta _{\phi }W_\mathrm {ext}&=0 \qquad \forall \quad \delta \phi \, , \end{aligned}$$

(13)
where the internal and external contributions are given in spatial form as


$$\begin{aligned}&\delta _{\varvec{\varphi }}W_\mathrm {int}=\int \limits _{\mathcal {B}_t}{\varvec{\sigma }}:\nabla _{\!\varvec{x}}\delta {\varvec{\varphi }}\,\mathrm {d}v\, , \quad \quad \delta _{\varvec{\varphi }}W_\mathrm {ext}=\int \limits _{\mathcal {B}_t}\bar{{\varvec{b}}}\cdot \delta {\varvec{\varphi }}\,\mathrm {d}v+\int \limits _{\partial \mathcal {B}_t}\bar{{\varvec{t}}}\cdot \delta {\varvec{\varphi }}\,\mathrm {d}a\, ,\end{aligned}$$

(14)



$$\begin{aligned}&\delta _{\phi }W_\mathrm {int}=\int \limits _{\mathcal {B}_t}{\varvec{y}}\cdot \nabla _{\!\varvec{x}}\delta \phi \,\mathrm {d}v\, , \quad \quad \delta _{\phi }W_\mathrm {ext}=\int \limits _{\mathcal {B}_t}y\,\delta \phi \,\mathrm {d}v\, . \end{aligned}$$

(15)
Here the Cauchy stress $${\varvec{\sigma }}$$ and the vectorial damage quantity $${\varvec{y}}$$ are related to flux terms, whereas the body force $$\bar{{\varvec{b}}}$$ and the scalar damage quantity $$y$$ are associated to source terms. They are defined as


$$\begin{aligned}&{\varvec{\sigma }}=\partial _{\varvec{F}}\varPsi \cdot \mathrm {cof}{({\varvec{F}}^{\,-1})}\, , \qquad \quad \bar{{\varvec{b}}}= -\,J^{\,-1}\partial _{\varvec{\varphi }}\varPsi _\mathrm {vol}\, , \end{aligned}$$

(16)



$$\begin{aligned}&{\varvec{y}} = \partial _{\nabla _{\!\varvec{X}}\phi }\varPsi \cdot \mathrm {cof}{({\varvec{F}}^{\,-1})}\, , \qquad y =-\,J^{\,-1} \partial _{\phi }\varPsi \, . \end{aligned}$$

(17)
Relations (14) and (15) provide the basis for the finite element discretisation in Sect. 4.


3 Constitutive Relations


In this section, we first review the hyperelastic constitutive equations adopted on the basis of an isotropic neo-Hookean relation and an anisotropic exponential part. These relations characterise the elastic anisotropic response of the fibre-reinforced material. Secondly, we specify the gradient-enhanced, non-local contribution to the free energy, followed by the continuum damage formulation.


3.1 Hyperelastic Part of the Free Energy


From Sect. 2.2, we recall the local free energy density $$\varPsi _\mathrm {loc}$$, Eq. (1), to be additively composed of an isotropic part $$\varPsi _\mathrm {iso}$$, representing the isotropic matrix material, and of an anisotropic part $$\varPsi _\mathrm {ani}$$, representing the individual families of fibres. In the following, we assume the isotropic part to be specified by a compressible neo-Hookean format


$$\begin{aligned} \varPsi _\mathrm {iso}=\frac{\mu _\mathrm {e}}{2}\left[ I_1-3\right] -\mu _\mathrm {e}\ln (J)+\frac{\lambda _\mathrm {e}}{2}\,[\ln (J)]^2\, , \end{aligned}$$

(18)
with $$I_1={\varvec{F}}:\varvec{F}$$ denoting the first invariant. The elastic constants are represented by the Lamé-parameters $$\mu _\mathrm {e}$$ and $$\lambda _\mathrm {e}=\kappa _\mathrm {e}-2/3\,\mu _\mathrm {e}$$ in terms of the shear modulus $$\mu _\mathrm {e}$$ and the bulk modulus $$\kappa _\mathrm {e}$$. The anisotropic contribution of the local free energy (1) is based on an orthotropic exponential model with two families of fibres including fibre dispersion according to Gasser et al. [5] or Menzel et al. [12], i.e.


$$\begin{aligned} \varPsi _\mathrm {ani}=\frac{k_1}{2\,k_2}\sum _{i=1}^N\left[ \,\exp \left( k_2\left\langle E_i\right\rangle ^2\right) -1\,\right] \, , \end{aligned}$$

(19)
with the strain-like quantity $$E_i=\varkappa \, I_1+[1-3\varkappa ]\,I_{4\,i}-1$$ and the invariant $$I_{4\,i}={\varvec{a}}_{0\,i}\,\cdot \,{\varvec{F}}^\mathrm {t}\cdot {\varvec{F}} \,\cdot \,{\varvec{a}}_{0\,i}$$ for $$N=2$$ fibre families. The term $$\left\langle E_i\right\rangle $$, where $$\left\langle \bullet \right\rangle =\left[ \left| \bullet \right| +\bullet \right] /2$$ is the Macaulay bracket, reflects the assumption that fibres can support tension only. Consequently, $$\varPsi _\mathrm {ani}>0$$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq60.gif”></SPAN> only if the fibre-related strain is positive, i.e. <SPAN id=IEq61 class=InlineEquation><IMG alt=0$$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq61.gif”>. Fibre dispersion is introduced by means of the parameter $$\varkappa \in [0,\,1/3]$$, where $$\varkappa =0$$ corresponds to no dispersion, i.e. transverse isotropy, and where $$\varkappa =1/3$$ renders an isotropic distribution. Table 1 summarises the structural and elastic material quantities included in constitutive Eqs. (18) and (19) together with their units. It is important to note that the fibre orientations may be defined arbitrarily, but the present formulation uses only one non-local damage variable so that both fibre families undergo identical degradation. This is physically meaningful as long as both families of fibers possess one and the same stretch history, otherwise a second non-local damage variable should be included in the formulation.


Table 1
Constitutive parameters as used in Sects. 3.1–3.3


























































Type
Symbol
Description
Unit

Structural
$${\varvec{a}}_{0\,i}$$

Fibre orientation vectors
–

$$\varkappa $$

Dispersion parameter
–

Elastic
$$\mu _\mathrm {e}$$

Shear modulus
$$\text {kPa}$$

$$\kappa _\mathrm {e}$$

Bulk modulus
$$\text {kPa}$$

$$k_1$$

Elastic constant
$$\text {kPa}$$

$$k_2$$

Elastic constant
–

Regularisation
$$c_\mathrm {d}$$

Degree of regularisation
$$\text {kPa}^{-1}\text {mm}^2$$

$$\beta _\mathrm {d}$$

Penalty parameter
$$\text {kPa}^{-1}$$

Damage
$$\eta _\mathrm {d}$$

Saturation parameter
$$\text {kPa}^{-1}$$

$$\kappa _\mathrm {d}$$

Damage threshold
$$\text {kPa}$$


3.2 Gradient-Enhanced Part of the Free Energy


According to Eq. (2), we assume the non-local part of the free energy to be additively composed of a gradient-related term $$\varPsi _\mathrm {grd}$$ and of a penalty term $$\varPsi _\mathrm {plty}$$ and specify these terms as


$$\begin{aligned} \varPsi _\mathrm {grd}(\nabla _{\!\varvec{X}}\phi ; {\varvec{F}}) =\frac{c_\mathrm {d}}{2}\,||\nabla _{\!\varvec{x}}\phi \,||^2 \quad \text {and} \quad \varPsi _\mathrm {plty}(\phi ,\kappa )=\frac{\beta _\mathrm {d}}{2}\,[\phi -\kappa ]^2\, . \end{aligned}$$

(20)
The energy-related penalty parameter $$\beta _\mathrm {d}$$ approximately enforces the local damage field $$\kappa $$ and the non-local field $$\phi $$ to coincide. Furthermore, the gradient parameter $$c_\mathrm {d}$$ controls the quasi-non-local character of the formulation and characterises the degree of gradient regularisation: $$c_\mathrm {d}=0$$ results in a local model, while $$c_\mathrm {d}> 0$$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq89.gif”></SPAN> leads to the regularised gradient-enhanced model. The damage-related parameters included in constitutive Eq. (<SPAN class=InternalRef><A href=20) together with their units are summarised in Table 1.


3.3 Gradient-Enhanced Damage Model


In order to obtain the stress-like thermodynamic forces driving the local dissipative damage process, we follow the standard Coleman-Noll procedure. Differentiation of the general format of the free energy (3) with respect to time and application of the Clausius-Planck inequality yields, amongst others, a contribution including the thermodynamic force $$q=q_\mathrm {loc}+q_\mathrm {nloc}=-\,\partial _d\varPsi $$ conjugate to the damage variable $$d$$, i.e.


$$\begin{aligned} q_\mathrm {loc}=\varPsi _\mathrm {ani}\quad \text {and} \quad q_\mathrm {nloc}=\beta _\mathrm {d}\,[\phi -\kappa ]\,\partial _d\kappa \, . \end{aligned}$$

(21)
Furthermore, we adopt the damage condition


$$\begin{aligned} \varPhi _\mathrm {d}=q-\kappa \le 0 \end{aligned}$$

(22)
where $$\varPhi _\mathrm {d}<0$$ refers to the purely elastic loading and $$\varPhi _\mathrm {d}=0$$ to damage evolution. Based on the postulate of maximum dissipation, we construct a constrained optimisation problem involving the Lagrange multiplier $$\lambda $$. This results in the following associated evolution equation for the damage variable


$$\begin{aligned} \dot{\kappa }=\lambda \,\frac{\partial \varPhi _\mathrm {d}}{\partial q}=\lambda \quad \text{ with } \quad \kappa \vert _{t=0}=\kappa _\mathrm {d}\, . \end{aligned}$$

(23)
where initiation and termination of damage are governed by the Karush-Kuhn-Tucker complementary conditions


$$\begin{aligned} \lambda \ge 0\, ,\quad \varPhi _\mathrm {d}\le 0\, ,\quad \lambda \,\varPhi _\mathrm {d}=0\, . \end{aligned}$$

(24)
We assume an exponential behaviour for the damage function


$$\begin{aligned} f_\mathrm {d}(\kappa )=1-d=\exp (\eta _\mathrm {d}\,[\kappa _\mathrm {d}-\kappa ])\, , \end{aligned}$$

(25)
with $$\eta _\mathrm {d}>0$$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq95.gif”></SPAN> so that <SPAN id=IEq96 class=InlineEquation><IMG alt=0$$” src=”/wp-content/uploads/2016/10/A324942_1_En_2_Chapter_IEq96.gif”> and introduce an initial damage threshold $$\kappa _\mathrm {d}$$, which must be exceeded in order to activate damage evolution. Furthermore, we include a saturation parameter $$\eta _\mathrm {d}$$. It becomes apparent that larger values of $$\eta _\mathrm {d}$$ accelerate the damage process, whereas larger values of $$\kappa _\mathrm {d}$$ lead to a delay of the damage initiation. Note, that for the limiting case $$\kappa _\mathrm {d}=0$$, damage is initiated from the very beginning of the loading process, whereas damage does not evolve for $$\kappa _\mathrm {d}\rightarrow \infty $$.


4 Finite Element Discretisation


This section deals with the spatial finite element discretisation of the underlying coupled system of non-linear equations. This includes a combination of tri-quadratic serendipity interpolation functions with respect to the displacement field, and tri-linear interpolation functions with respect to the non-local damage field variable where we outline an efficient and compact FE-implementation using a common Voigt-notation-based vector-matrix-format.


4.1 Discretisation


We discretise the domain $$\mathcal {B}_0$$ by $$n_\mathrm {el}$$ finite elements, so that $$\mathcal {B}_0\approx \mathcal {B}_0^\mathrm {h}=\bigcup _{e=1}^{n_\mathrm {el}} \, \mathcal {B}_{0\,e}$$, where every finite element $$\mathcal {B}_{0\,e}$$ is characterised by $$n_\mathrm {en}^\varphi $$ placement-nodes and $$n_\mathrm {en}^\phi $$ non-local-damage-nodes. According to the isoparametric concept, we interpolate the field variables $${\varvec{\varphi }}$$ as well as the geometry $${\varvec{X}}$$ by the same shape functions $$N^\varphi $$
Tags:
Oct 21, 2016 | Posted by in BIOCHEMISTRY | Comments Off on Gradient-Enhanced Continuum Damage Model for Residually Stressed Fibre-Reinforced Materials at Finite Strains

Full access? Get Clinical Tree
Get Clinical Tree app for offline access