A convex anisotropic damage model based on the compliance tensor

Abstract

This work is devoted to anisotropic continuum-damage mechanics in the quasi-static, isothermal, small-strain setting. We propose a framework for anisotropic damage evolution based on the compliance tensor as primary damage variable, in the context of generalized standard models for dissipative solids. Based on the observation that the Hookean strain energy density of linear elasticity is jointly convex in the strain and the compliance tensor, we design thermodynamically consistent anisotropic damage models that satisfy Wulfinghoff’s damage-growth criterion and feature a convex free energy. The latter property permits obtaining mesh-independent results on component scale without the necessity of introducing gradients of the damage field. We introduce the concepts of stress-extraction tensors and damage-hardening functions, implicitly describing a rigorous damage-analogue of yield surfaces in elastoplasticity. These damage surfaces may be combined in a modular fashion and give rise to complex damage-degradation behavior. We discuss how to efficiently integrate Biot’s equation implicitly, and show how to design specific stress-extraction tensors and damage-hardening functions based on Puck’s anisotropic failure criteria. Last but not least we demonstrate the versatility of our proposed model and the efficiency of the integration procedure for a variety of examples of interest.

Introduction

State of the art

Damage mechanics describes the progressive degradation of the elastic stiffness of materials upon loading, and is typically attributed to growing voids or cracks on a lower length scale (Lemaitre, 1996), see Figure 1. There are two predominant approaches to continuum damage-mechanics (Krajcinovic, 1984; Lemaitre and Chaboche, 1990). The first approach accounts for the origin of damage on a lower length scale in terms to micromechanics (Fitoussi et al., 1996; Guo et al., 1997), see also Section 3 in Krajcinovic (1989) for an early account. With qualitative predictions in mind, the second strategy is of phenomenological nature. After selecting a suitable damage variable (or a collection thereof), suitable kinetic laws are postulated taking continuum thermodynamics into account, Section 4 in Krajcinovic (1989).

Figure 1.

Schematics of microstructures with growing microscopic cracks, passing from state ① to state ②, similar to Fassin et al. (2019). Growing microcracks induce a reduction of the effective stiffness. (a) Polycrystalline microstructure. (b) Microstructure composed of fiber bundles.

The micromechanics-based approach to damage mechanics takes the damage mechanisms on a lower scale into account and is still subject of current research, for instance concerning mesh-size objective modeling (Liang et al., 2018), a coupling to model-order reduction (Bhattacharyya et al., 2020) or accounting for micro-computed tomography data (Luo et al., 2020). Micromechanics-informed damage models permit taking the stochastics on the microscale into account naturally, e. g., for progressive fiber breakage in fiber-reinforced composites (Ju and Wu, 2016; Wu and Ju, 2017), interfacial transition-zone effects (Chen et al., 2018), uncertainty in the elastic moduli of fiber-reinforced concrete (Liu et al., 2020), localized microcracks (Li et al., 2020) or random loading in fatigue processes (Franko et al., 2017). Another advantage concerns modeling the unilateral character of brittle damage, i. e., a different damaging behavior under tension compared to compression (Goidescu et al., 2015; Zhang et al., 2019), and accounting for interface debonding (Pupurs and Varna, 2017; Schemmann et al., 2018b; Yang et al., 2020). However, care has to be taken as homogenization and localization are incompatible (Gitman et al., 2007), in general, i.e., upon localization, the volume elements considered will not be representative for the effective mechanical behavior (Drugan and Willis, 1996; Hill, 1963; Kanit et al., 2003).

As an alternative to micromechanics-type strategies, phenomenological approaches to continuum-damage mechanics may be pursued. In a first step, a (scalar- or tensor-valued) damage variable is selected which describes the reduction of the effective cross-section of a typical material sample undergoing material degradation (Gurson, 1977; Voyiadjis, 2015). Then, suitable kinetic laws are postulated on the basis of continuum thermodynamics (Hansen and Schreyer, 1994; Simo and Ju, 1987).

The tensor order of the damage variable naturally distinguishes different phenomenological damage models. Even today, the classical scalar isotropic damage variable serves as a reliable workhorse with numerous applications including cast steel with pores (Yan et al., 2020), concrete (Li and Wu, 2018), rocks (Liu et al., 2018; Xu et al., 2018), framed structures (Yang et al., 2017), unidirectional glass fiber-reinforced plastic composite plies (Sharma and Daggumati, 2020), fibrous composite laminae (Abu-Farsakh and Asfa, 2020), notched epoxy resin specimens (Rahimi et al., 2020) and steel-fiber reinforced concrete (Moradi et al., 2020).

Damage variables with higher tensor order permit modeling an emerging anisotropy of damage. As working with second-order tensors comes naturally to disciples of continuum mechanics, it is not surprising that second-order damage-tensors (Murakami and Ohno, 1981) are used frequently in continuum damage-mechanics. Recent applications include concrete (Desmorat, 2016; Wardeh and Toutanji, 2017), metal-forming processes (Nasab and Mashayekhi, 2019), rock materials (Wang and Xu, 2020), composite fabrics and laminated panels (Wei et al., 2020) and composite laminates (Okabe et al., 2018; Onodera and Okabe, 2020). Second-order damage-tensors are always orthotropic w. r. t. their eigenbasis, limiting their degree of generality. More often than not, such a limitation is interpreted as a feature, and specific orthotropic damage models are developed, for instance for brittle materials (Kim et al., 2016), in elastoplastic and finite-strain damage coupling (Ganjiani, 2018; Reese et al., 2021), or for ceramic-matrix composites (Alabdullah and Ghoniem, 2020).

As continuum damage-models primarily seek to describe a loss of stiffness due to emerging defects in solids, using a fourth-order damage-tensor (Chaboche, 1981), the same tensor order as the stiffness tensor, appears reasonable. In Section 4.3.4, Krajcinovic (1989) even notes that “an appropriate description of damage […] must involve at least a fourth-rank tensor.” This idea was pursued for the stiffness or compliance tensors as the primary damage variable (Dougill, 1976; Ortiz, 1985; Ortiz and Popov, 1982), also coupled to plasticity (Ju, 1989; Simo and Ju, 1987; Yazdani and Schreyer, 1990). We refer to Zhang and Cai (2010) for a modern account of anisotropic damage mechanics. However, some care has to be taken when working with tensor-valued damage variables due to possible inconsistencies arising for complex non-radial loading-unloading scenarios, see Simon et al. (2017).

The unilateral character of pores and cracks (see Figure 1) often leads to a tension-compression asymmetry of the material behavior upon damage loading, see Chaboche (1993) for a discussion. To incorporate the latter effect in phenomenological models, one may introduce different damage variables for the tensile and the compressive regime (Cicekli et al., 2007; Ramtani et al., 1992). For three-dimensional stress states, spectral decompositions of either the strain or the stress tensor may form the basis of continuum damage models that differentiate between damage evolution due to tension and compression (Ladeveze and Lemaitre, 1984; Ortiz, 1985).

Whenever damage models exhibit a softening behavior, their use in a continuum formulation leads to an ill-posed problem due to localization effects (Lemaitre, 1986), which is reflected by strongly mesh-dependent results in numerical simulations (De Borst, 1996). Countermeasures in the framework of local damage models were investigated (Becker et al., 1988; Beremin et al., 1983; Tvergaard, 1982). Non-local formulations (Bažant, 1991; Belytschko et al., 1986) prevent the localization responsible for the ill-posedness, and may be realized by an explicit convolution with a tapering function (Pijaudier-Cabot and Bažant, 1987), by augmenting the damage evolution equation by an elliptic differential operator (Aifantis, 1984) or by employing a gradient-enhanced formulation (Abu Al-Rub and Voyiadjis, 2009; Brünig and Ricci, 2005; Germain et al., 2007), which may also be coupled to Hamilton’s least-action principle (Junker et al., 2019, 2021). As long as the softening is not too pronounced, existence of results for non-local damage models (Thomas and Mielke, 2010) may be established. However, except for specific models (Roubíček, 2009; Susu, 2017), uniqueness (and, thus, well-posedness) cannot be ensured. For a review on ill-posedness due to localization problems and appropriate regularization methods, the reader is referred to Forest et al. (2004). Also, for a general overview on continuum damage-mechanics and further literature, the reader may consult the books of Murakami (2012) and of Voyiadjis (2015).

Oftentimes, the ill-posedness of local damage models is taken for granted, and appropriate countermeasures are taken. A charming strategy takes a conventional local damage model with softening (but sufficient growth at infinity), and applies relaxation techniques (Balzani and Ortiz, 2012; Schmidt and Balzani, 2016; Schwarz et al., 2021), which are typically used for studying solids with emerging microstructure. When describing stable damage processes, these countermeasures should not be neccesary, however. Indeed, for a moderate degree of loading, localization is excluded, and manifests only at a specific turning point in loading level. For component-scale simulations, this loading level is not readily apparent, and depends on the specimen geometry via solving the equations of continuum mechanics. To sum up and loosely speaking, we know that local damage models are perfectly reasonable up to a specific level of loading, but we do not know this level in advance. Thus, interest arose to design damage models which give rise to a meaningful response for the entire range of loading, and which are intended to be complemented by a classical failure criterion.

Contributions and organization of this article

We contribute to phenomenological continuum damage-mechanics with a tensorial damage variable. We advocate using the full compliance tensor as a rather natural and observable damage variable, liberating the engineer of the burden of selecting the appropriate damage variable in advance, permitting her to focus the attention on appropriate kinetic laws. Thus, when it comes to continuum damage-mechanics of phenomenological type, the proposed framework is as ab-initio as possible, since only the evolution of the damage surface in stress space needs to be identified.

The compliance tensor has been used as the primary damage variable before (Baranger, 2018; Ladevèze, 1983, 2002; Ladevèze et al., 2014). Yet, this approach has not yet entered the main stream of damage-modeling frameworks. Our theoretical contributions to compliance-based damage models are actually twofold. For a start, we point out that the standard Hookean strain energy density, regarded as a function of the strain tensor and the full compliance tensor, is de facto jointly convex in both arguments. This result is surprising, and we are not aware of an account in the literature (although we sincerely believe that others have presumably noticed this fact before without stating it explicitly, see Thomas and Mielke (2010) for a special case).

Based on the compliance tensor, we develop a simple, modular framework for anisotropic damage mechanics. The framework provides the working engineer with a number of options which we believe to be of advantage. Indeed, due to the convexity property of the Hookean elastic energy, it is possible to develop a purely hardening damage-mechanics modeling-framework, where localization does not become an issue. Very much, there are materials which show a purely damage-hardening material response prior to sudden and brutal failure, e.g., Sheet Molding Compound (SMC) composites (Anagnostou et al., 2018; Fitoussi et al., 1996, 1998) comprising an unsaturated polyester-polyurethane hybrid (UPPH) resin (Kehrer et al., 2018; Schemmann et al., 2018a; Trauth et al., 2017) reinforced by glass fibers (Görthofer et al., 2019; Meraghni and Benzeggagh, 1995; Schemmann et al., 2018c).

Of course, the modeling framework is not restricted to damage-hardening, but may be adapted to softening in a straightforward manner. However, the latter scenario is rather classical in continuum damage-mechanics, and we decided to work out the details of a hardening framework in the paper at hand, essentially due to our desire to model SMC materials.

To highlight the simplicity of our proposed compliance-type damage modeling framework, we present a first-principles development in the context of generalized standard models (GSMs) for dissipative solids (Halphen and Nguyen, 1975) and discuss the efficient resolution of the evolution equations in a predictor-corrector framework.

Our second contribution concerns a design methodology for the damage surfaces which draws upon similar approaches in (associative) elastoplasticity (Bertram, 2011; Chaboche, 2008; McDowell, 2008), but takes failure criteria and multiple damage surfaces (Bakhshan et al., 2018; Jin and Arson, 2018; Khayyam Rayeni et al., 2020) into account. More precisely, building upon Puck’s anisotropic failure criteria developed for continuously reinforced polymers (Knops, 2008; Puck and Schürmann, 2002), we design specific damage-extraction tensors and damage-activation functions which present a flexible arsenal of tools that, taking the individual damaging mechanisms into consideration, permit building up an accurate and fully anisotropic continuum damage model.

For anisotropic damage models not to be judged as purely academic, it is of utmost importance to establish links to experimental data and to compare it to (dis)similar modeling approaches. After conducting computational investigations which clarify the influence of the different model parameters on the damage evolution and expose the developing elastic anisotropy upon loading, we study a plain-weave mesostructure of a woven carbon-fiber reinforced thermoset investigated by Simon et al. (2017). We show that the convex modeling framework permits reproducing the effective mechanical behavior of the individual tows and the composite quantitatively within the loading range of interest.

Notation

We follow a direct tensor notation throughout the text, representing vectors and tensors by their components or using matrix representations (in an orthonormal basis) only when necessary. Vectors and second-order tensors are denoted by lower case and upper case bold letters, respectively (e.g., a and A ). Fourth-order tensors are denoted by, e.g., $\mathbb{A},\mathbb{B}$ . Scalars and arrays of quantities are represented by non-bold letters (e.g., H, w or z). The transposition of a vector and second-order tensor reads ${\mathit{a}}^{\mathrm{T}}$ and ${\mathit{A}}^{\mathrm{T}}$ , respectively. The principal transposition of a fourth-order tensor is denoted via ${\mathbb{A}}^{{\mathrm{T}}_{\mathrm{M}}}$ and the left and right transpositions are ${\mathbb{A}}^{{\mathrm{T}}_{\mathrm{L}}}$ and ${\mathbb{A}}^{{\mathrm{T}}_{\mathrm{R}}}$ . The linear mappings induced by second-order and fourth-order tensors are written as $\mathit{a}=\mathit{Cb}$ and $\mathit{A}=ℂ[\mathit{B}]$ , respectively. The composition of two second-order or two fourth-order tensors is denoted by $\mathit{AB}$ and $\mathbb{A}\mathbb{B}$ . The Frobenius inner product is denoted by $\mathit{A}·\mathit{B}=\text{tr}(\mathit{A}{\mathit{B}}^{\mathrm{T}})$ . The tensor product is symbolized by ⊗. Its symmetrized version ${\otimes }_{\mathrm{S}}$ is defined via $\mathit{a}{\otimes }_{\mathrm{S}}\mathit{b}=\left(\mathit{a}\otimes \mathit{b}+\mathit{b}\otimes \mathit{a}\right)/2$ . We introduce the abbreviation ${\mathit{a}}^{\otimes n}=\mathit{a}\otimes \mathit{a}\dots \otimes \mathit{a}$ (n repetitions). The material time derivative of a quantity w is expressed as $\dot{w}=\text{d}w/\text{d}t$ . We denote by $\text{Sym}(d)$ the space of symmetric second-order tensors on ${\mathbb{R}}^d$ . The unit sphere in ${\mathbb{R}}^3$ reads $S^2$ . The vector space of fourth-order tensors with minor symmetries ( $\mathbb{A}={\mathbb{A}}^{{\mathrm{T}}_{\mathrm{L}}},\mathbb{A}={\mathbb{A}}^{{\mathrm{T}}_{\mathrm{R}}}$ ) is written as $L(\text{Sym}(d))$ , whereas $\text{Sym}(\text{Sym}(d))$ denotes those-fourth order tensors that have minor and major symmetries ( $\mathbb{A}={\mathbb{A}}^{{\mathrm{T}}_{\mathrm{M}}}$ ). In general, details on further spaces of interest, domains of definition and corresponding explicit expressions are given upon their first appearance.

A compliance-based anisotropic damage model

A convex standard model for anisotropic damage

We will describe the damage model, in a small-strain and isothermal setting, as a generalized standard model (GSM) (Halphen and Nguyen, 1975), whose framework we briefly recall. In addition to the symmetric $d\times d$ infinitesimal strain tensor $\mathit{\epsilon }\in \text{Sym}(d)$ , where d = 2, 3 denotes the dimension of the ambient space, a (Banach) space $\mathcal{Z}$ of internal variables is postulated. Furthermore, a free energy (density)

$$w:\text{Sym}(d)\times \mathcal{Z}\to \mathbb{R},(\mathit{\epsilon },z)\hspace{0.17em}\mapsto \hspace{0.17em}w(\mathit{\epsilon },z),$$

(2.1)

a continuously differentiable function of the strain tensor $\mathit{\epsilon }$ and the internal variables $z\in \mathcal{Z}$ , and a force potential ${\Phi }^{*}:\mathcal{Z}\prime \to [0,\infty ]$ , a lower semicontinuous, non-negative and convex function on the continuous dual space $\mathcal{Z}\prime$ satisfying ${\Phi }^{*}(0)=0$ , are introduced. To ensure thermodynamic consistency, the Clausius-Duhem inequality (CDI), see Chapter 13 in Haupt (2000),

$$0\stackrel{!}{\le }\mathit{\sigma }·\dot{\mathit{\epsilon }}-\frac{\text{d}}{\text{d}t}\left[w(\mathit{\epsilon },z)\right]\equiv \left[\mathit{\sigma }-\frac{\partial w}{\partial \mathit{\epsilon }}(\mathit{\epsilon },z)\right]·\dot{\mathit{\epsilon }}-\frac{\partial w}{\partial z}(\mathit{\epsilon },z)·\dot{z},$$

(2.2)

where $\dot{z}\equiv \text{d}z/\text{d}t$ denotes the material time derivative of the internal variables, needs to be satisfied. Associated to a current equilibrium state $(\mathit{\epsilon },z)$ of a hyperelastic material, the Cauchy stress tensor $\mathit{\sigma }\in \text{Sym}(d)$ is defined (Halphen and Nguyen, 1975; Miehe, 2002) by

$$\mathit{\sigma }=\frac{\partial w}{\partial \mathit{\epsilon }}(\mathit{\epsilon },z).$$

(2.3)

For a prescribed loading path $\mathit{\epsilon }:[0,T]\to \text{Sym}(d)$ on a given interval of time and the initial condition $z(0)=z_0$ for some $z_0\in \mathcal{Z}$ , the evolution of the internal variables is governed by Biot’s (dual) equation

$$\dot{z}\in \partial {\Phi }^{*}\left(-\frac{\partial w}{\partial z}(\mathit{\epsilon },z)\right),$$

(2.4)

where $\partial {\Phi }^{*}$ denotes the subdifferential of the convex function ${\Phi }^{*}$

$$\partial {\Phi }^{*}(\xi )=\left\{z\in \mathcal{Z}|{\Phi }^{*}\left(\tilde{\xi }\right)-{\Phi }^{*}(\xi )\ge \left(\tilde{\xi }-\xi \right)·z\text{for all}\tilde{\xi }\in \mathcal{Z}\prime \right\},$$

(2.5)

a subset of $\mathcal{Z}\prime$ , see Borwein and Lewis (2006) for details. Due to these definitions, see Halphen and Nguyen (1975), generalized standard materials are automatically thermodynamically consistent. Indeed, by Biot’s (dual) equation (2.4),

$${\Phi }^{*}(0)-{\Phi }^{*}(\xi )\ge (0-\xi )·\dot{z}\text{ for }\xi =-\frac{\partial w}{\partial z}(\mathit{\epsilon },z)$$

(2.6)

holds. Using ${\Phi }^{*}(0)=0$ , rearranging the latter inequality yields

$$-\frac{\partial w}{\partial z}(\mathit{\epsilon },z)·\dot{z}\equiv \xi ·\dot{z}\ge {\Phi }^{*}(\xi )\ge 0,$$

(2.7)

i.e., the Clausius-Duhem inequality (2.2) holds in view of the definition of stress (2.3). As it drives the evolution of the internal variables, the quantity $\xi \equiv -\frac{\partial w}{\partial z}(\mathit{\epsilon },z)$ is called driving force.

As internal variables z of our proposed continuum damage-mechanics model, we consider an elastic compliance tensor

$$\mathbb{S}\in {\mathcal{S}}_d=\left\{\mathbb{S}\in \text{Sym}(\text{Sym}(d))|\mathit{\tau }·\mathbb{S}\left[\mathit{\tau }\right]>0\text{for all}\mathit{\tau }\in \text{Sym}(d)\backslash \{0\}\right\},$$

(2.8)

and a general variable $q\in \mathcal{Q}$ which describes the shape and size of the damage surface, s. t.

$$z=\left(\mathbb{S},q\right)\in {\mathcal{S}}_d\times \mathcal{Q}.$$

(2.9)

Notice that the set ${\mathcal{S}}_d$ of (positive definite) compliance tensors is not a linear space. Instead, it is an open, convex subset of the linear space of fourth-order tensors $\text{Sym}(\text{Sym}(d))$ with minor and major symmetries.

For a GSM, the CDI (2.2) will always be satisfied. However, we need to ensure that the (dual) Biot’s equation (2.4) guarantees that $\mathbb{S}$ remains an element of ${\mathcal{S}}_d$ , i.e., that the compliance tensor $\mathbb{S}$ remains positive definite. In contrast, the damage-surface variables we consider live in a linear space $\mathcal{Q}$ (which we deliberately keep abstract). For the specific models presented in the next section q is just a finite collection of scalar values. However, our arguments cover the more general case, accounting for vector- or tensor-valued damage-surface variables in a natural way.

The free energy (density) we consider is defined by

$$w:\text{Sym}(d)\times {\mathcal{S}}_d\times \mathcal{Q}\to \mathbb{R},(\mathit{\epsilon },\mathbb{S},q)\hspace{0.17em}\mapsto \hspace{0.17em}{w}_e(\mathit{\epsilon },\mathbb{S})+h(q),$$

(2.10)

involving the Hookean elastic energy (density)

$$w_e:\text{Sym}(d)\times {\mathcal{S}}_d\to [0,\infty ),(\mathit{\epsilon },\mathbb{S})\hspace{0.17em}\mapsto \hspace{0.17em}\frac{1}{2}\mathit{\epsilon }·{\mathbb{S}}^{-1}\left[\mathit{\epsilon }\right],$$

(2.11)

and an energy (density) related to the progressive degradation of the material,

$$h:\mathcal{Q}\to \mathbb{R},q\hspace{0.17em}\mapsto \hspace{0.17em}h(q),$$

(2.12)

which we assume to be convex and continuously differentiable. Notice that the Hookean elastic energy w_e (2.11) is jointly convex in both variables and infinitely often differentiable. The latter property is immediate ¹ , as w_e depends on $\mathit{\epsilon }$ quadratically and the Neumann-series representation

$${(\mathbb{S}+\mathbb{L})}^{-1}={\displaystyle \sum _{k=0}^{\infty }{\left(-{\mathbb{S}}^{-1}\mathbb{L}\right)}^k}{\mathbb{S}}^{-1},\mathbb{S}\in {\mathcal{S}}_d,\mathbb{L}\in \text{Sym}(\text{Sym}(d)),$$

(2.13)

valid for sufficiently small $\mathbb{L}$ , shows that w_e is even analytic. For the convexity, recall that a twice differentiable function defined on an open convex set is convex if and only if its Hessian is positive semidefinite everywhere, see Theorem 3.1.11 in Borwein and Lewis (2006). A direct computation shows that the Hessian admits the representation

$$D^2{w}_e(\mathit{\epsilon },\mathbb{S})[\mathit{\xi },\mathbb{L}]=\frac{1}{2}\left(\mathit{\xi }-\mathbb{L}{\mathbb{S}}^{-1}\left[\mathit{\epsilon }\right]\right)·{\mathbb{S}}^{-1}\left[\mathit{\xi }-\mathbb{L}{\mathbb{S}}^{-1}\left[\mathit{\epsilon }\right]\right],$$

(2.14)

for $(\mathit{\epsilon },\mathbb{S})\in \text{Sym}(d)\times {\mathcal{S}}_d$ and $(\mathit{\xi },\mathbb{L})\in \text{Sym}(d)\times \text{Sym}(\text{Sym}(d))$ , see Appendix B. Any $\mathbb{S}\in {\mathcal{S}}_d$ is positive definite, and thus, the Hessian in equation (2.14) is non-negative. Consequently, the elastic energy is convex (but not strictly convex). As we assumed the energy h (2.12) to be continuously differentiable and convex, the smoothness and convexity properties of w_e imply that the free energy w (2.10) is continuously differentiable and convex, as well. Furthermore, the formula for the Cauchy stress (2.3) becomes

$$\mathit{\sigma }=\frac{\partial w}{\partial \mathit{\epsilon }}(\mathit{\epsilon },\mathbb{S},q)\equiv {\mathbb{S}}^{-1}\left[\mathit{\epsilon }\right],$$

(2.15)

i.e., for a fixed compliance tensor $\mathbb{S}$ , the stress-strain relationship reduces to Hooke’s law. To conclude this paragraph, several remarks are in order.

1. Using the framework of generalized standard materials for phenomenological modeling of damage is classical. For instance, Hansen and Schreyer (1994) study a general tensor-valued damage variable coupled to plasticity in such a framework, apparently unaware of the connection.

2. In phenomenological continuum damage-mechanics, choosing the damage variable typically comes first, and the damage kinetics needs to be set up based on the resulting driving forces. Our approach frees the reader of an a priori selection of damage variable, and permits her to focus on the kinetics in terms of the quite natural stress-based driving force.

3. The compliance tensor has the attractive characteristic that it is a physical quantity which can be determined experimentally. Of course, determining all 21 independent parameters of a stiffness tensor in three spatial dimensions is a daunting task from an experimental perspective. Still, observability of the damage variable is not ensured for purely phenomenological damage vectors and tensors.

4. The compliance tensor has been used as a damage variable before (Ladevèze, 1983, 2002). However, its use seemed restricted to specific situations, e.g., damage modeling in ceramic-matrix composites (Baranger, 2018; Ladevèze et al., 2014). In this work, we advocate using the compliance tensor as the damage variable of choice in greater generality.

5. It is more than well-known that the Hookean energy (2.11) is convex in the strain tensor. It appears much less known that the Hookean energy is jointly convex in the strain and the compliance tensor. When coupled to an energy h which makes the condensed incremental potential strictly convex, the resulting framework produces an anisotropic damage model which does not permit localization. In particular, associated finite-element computations are not affected by mesh sensitivity induced by softening behavior. We do not want to argue against damage localization. Rather, we wish to add a powerful weapon to the arsenal of continuum damage-mechanics when it comes to modeling stable anisotropic damage phenomena.

6. In classical small-strain elasto(visco)plasticity the (visco)plastic strain ${\mathit{\epsilon }}_p\in \text{Sym}(d)$ serves as an internal variable. The corresponding stored energy (density)

   $$(\mathit{\epsilon },{\mathit{\epsilon }}_p)\hspace{0.17em}\mapsto \hspace{0.17em}\frac{1}{2}\left(\mathit{\epsilon }-{\mathit{\epsilon }}_p\right)·ℂ\left[\mathit{\epsilon }-{\mathit{\epsilon }}_p\right]$$

   (2.16)

with a fixed stiffness tensor $ℂ={\mathbb{S}}^{-1}$ is smooth and jointly convex in both arguments, but not strictly convex. The Hookean elastic stored energy function w_e (2.11) may be considered as a damage-analog of the elastic stored energy in classical elasto(visco)plasticity (2.16). The combined energy taking into account damage (2.11) and elasto(visco)plasticity (2.16) is jointly convex in all variables. If plasticity is neglected ( ${\mathit{\epsilon }}_p\equiv 0$ ) we recover the damage case and for a constant stiffness ( $ℂ={\mathbb{S}}^{-1}\equiv \text{const}.$ ) we recover classical elasto(visco)plasticity. Such a model differs from the classical presentation, which is typically based on either strain or energy equivalence (Hansen and Schreyer, 1994; Sections 3.2.1 and 3.2.2).

7. If we regard the Hookean elastic stored energy function $(\mathit{\epsilon },ℂ)\mapsto \frac{1}{2}\mathit{\epsilon }·ℂ\left[\mathit{\epsilon }\right]$ as a function of the stiffness tensor $ℂ$ , it will not be convex. Indeed, its Hessian at $(\mathit{\epsilon },ℂ)$ computes as

$$(\mathit{\xi },\mathbb{L})\mapsto 2\mathit{\xi }·\mathbb{L}\left[\mathit{\epsilon }\right]+\mathit{\xi }·ℂ\left[\mathit{\xi }\right],$$

(2.17)

which may become negative (take, for instance $\mathit{\xi }=\mathit{\epsilon }$ and $\mathbb{L}=-ℂ$ ). This lack of convexity is the reason why it is so difficult to design convex damage models for stable damage processes. Using the compliance tensor eradicates these issues with the help of a nonlinear transformation.

8. For the thermodynamics considerations at the beginning of this section to be valid, the “interfacial” energy (2.12) need not be convex, see, for example, Govindjee et al. (1995). In particular, softening behavior can be modeled in the compliance-tensor framework, as well. In that case, for obtaining a well-defined boundary-value problem, damage localization has to be overcome, for instance by adding gradient terms of the variable q to the energy (2.12).

9. The presented model cannot distinguish tensile and compressive loading. Indeed, the driving force $\mathbb{T}$ (2.18) for the compliance evolution computes as

$$\mathbb{T}=-\frac{\partial w(\mathit{\epsilon },\mathbb{S},q)}{\partial \mathbb{S}}=-\frac{\partial w_e(\mathit{\epsilon },\mathbb{S})}{\partial \mathbb{S}}=\frac{1}{2}\mathit{\sigma }\otimes \mathit{\sigma }\in \text{Sym}(\text{Sym}(d)),$$

(2.18)

which is insensitive to the involution $\mathit{\sigma }\mapsto -\mathit{\sigma }$ . Thus, in order to extend our model to account for tension-compression asymmetry, the free energy w requires a modification, see Ladevèze and co-workers (Ladevèze, 1983, 2002; Ladevèze et al., 2014).

To finish presenting the two-potential model, a suitable force potential ${\Phi }^{*}$ needs to be provided, entering the evolution equation of the internal variables, see equation (2.4),

$$(\dot{\mathbb{S}},\dot{q})\in \partial {\Phi }^{*}\left(\mathbb{T},\beta \right).$$

(2.19)

In the quasi-static setting targeting a rate-independent damage model, we describe the force potential ${\Phi }^{*}$ in terms of M continuously differentiable and convex damage functions ${\varphi }_i:\text{Sym}(\text{Sym}(d))\times \mathcal{Q}\prime \to \mathbb{R}$ , i.e.,

$${\Phi }^{*}(\mathbb{T},\beta )=\{\begin{array}{ll}0,& {\varphi }_i(\mathbb{T},\beta )\le 0\text{for all}i=1,\dots ,M,\\ +\infty ,& \text{otherwise}.\end{array}$$

(2.20)

Such a force potential gives rise to a quasi-static damage evolution in terms of an elastic domain defined by the functions ${\varphi }_i$ , in strict analogy to associated elastoplasticity at small strains, see Chapter 5 in Simo and Hughes (1998). A schematic of the admissible region based on the force potential (2.20) with corresponding driving forces $\mathbb{T}$ and $\beta$ is shown in Figure 2.

Figure 2.

Schematic of the admissible elastic region in ( $\mathbb{T},\beta$ )-space.

However, some care has to be exercised, as the elastic domain is defined in terms of the compliance driving-force $\mathbb{T}$ , which takes the form $\mathbb{T}=\frac{1}{2}\mathit{\sigma }\otimes \mathit{\sigma }$ for the free energy w (2.10), in contrast to elastoplasticity, where the stress tensor $\mathit{\sigma }$ (or a shifted version thereof) serves as the driving force. For the specific force potential ${\Phi }^{*}$ (2.20), Biot’s (dual) equation (2.4) becomes

$$\dot{\mathbb{S}}={\displaystyle \sum _{i=1}^M{\dot{\mu }}_i}\frac{\partial {\varphi }_i(\mathbb{T},\beta )}{\partial \mathbb{T}}\text{ and }\dot{q}={\displaystyle \sum _{i=1}^M{\dot{\mu }}_i}\frac{\partial {\varphi }_i(\mathbb{T},\beta )}{\partial \beta },$$

(2.21)

involving the driving forces

$$\beta =-\frac{\partial w(\mathit{\epsilon },\mathbb{S},q)}{\partial q}\equiv -\frac{\partial h(q)}{\partial q}\in \mathcal{Q}\prime$$

(2.22)

for the evolution of the damage-surface variables q and consistency parameters ${\mu }_1,\dots ,{\mu }_M$ which obey the Karush-Kuhn-Tucker (KKT) conditions

$${\dot{\mu }}_i\ge 0,{\varphi }_i(\mathbb{T},\beta )\le 0,{\dot{\mu }}_i{\varphi }_i(\mathbb{T},\beta )=0,i=1,\dots ,M.$$

(2.23)

To ensure that $\mathbb{S}$ remains in the set ${\mathcal{S}}_d$ of positive definite compliance tensors, a condition of the form

$$\frac{\partial {\varphi }_i(\mathbb{T},\beta )}{\partial \mathbb{T}}\ge 0\text{for all}i=1,\dots ,M,$$

(2.24)

on the damage functions ${\varphi }_i$ is sufficient. The latter condition was established by Wulfinghoff et al. (2017) as a criterion any physically meaningful vectorial or tensorial continuum-damage model should satisfy. In our context, the compliance tensor $\mathbb{S}$ serves as the damage variable, and Wulfinghoff’s criterion becomes “ $\dot{\mathbb{S}}\ge 0$ ”, i.e., $\dot{\mathbb{S}}$ is positive semidefinite.

To complete describing our model, we restrict the space of damage variables to $\mathcal{Q}={\mathbb{R}}^M$ , i.e., one scalar damage variable per damage-activation function ${\varphi }_i$ . We define the damage-activation function to be

$${\varphi }_i:\text{Sym}(\text{Sym}(d))\times \mathbb{R}\to \mathbb{R},(\mathbb{T},{\beta }_i)\hspace{0.17em}\mapsto \hspace{0.17em}2\mathbb{T}·{\mathbb{B}}_i^2-{\sigma }_{0,i}^2+H_i{\beta }_i,i=1,\dots ,M,$$

(2.25)

involving a (fourth-order, dimension-free) damage-extraction tensor ${\mathbb{B}}_i\in L(\text{Sym}(d))$ with minor and major symmetries, a damage-activation threshold ${\sigma }_{0,i}$ (analogous to the yield stress in elastoplasticity), and a positive parameter H_i with the dimensions of stress.

In principle, the damage-extraction tensor ${\mathbb{B}}_i$ need not have the major symmetry for equation (2.29) to make sense. In this non-symmetric case, the term ${\mathbb{B}}_i^2$ in equation (2.25) needs to be replaced by ${\mathbb{B}}_i^{{\mathrm{T}}_{\mathrm{M}}}{\mathbb{B}}_i$ in terms of the transpose ${\mathbb{B}}_i^{{\mathrm{T}}_{\mathrm{M}}}$ of the damage-extraction tensor ${\mathbb{B}}_i$ . However, the framework (2.25) may be recovered by defining ${\tilde{\mathbb{B}}}_i=\sqrt{{\mathbb{B}}_i^{{\mathrm{T}}_{\mathrm{M}}}{\mathbb{B}}_i}$ . Thus, by restricting the damage-extraction tensor to have major symmetries we do not lose generality. Furthermore, as we consider the variable q_i to be dimensionless, the associated driving force β_i has dimensions of stress and the parameter H_i is necessary for dimensional reasons.

In any case, for the damage function (2.25) the condition (2.24) to fulfill Wulfinghoff’s damage growth criterion, is automatically satisfied. Indeed, for any $i=1,\dots ,M$ , we obtain,

$$\mathit{\tau }·\frac{\partial {\varphi }_i(\mathbb{T},{\beta }_i)}{\partial \mathbb{T}}\left[\mathit{\tau }\right]=\mathit{\tau }·{\mathbb{B}}_i^2\left[\mathit{\tau }\right]={\mathbb{B}}_i\left[\mathit{\tau }\right]·{\mathbb{B}}_i\left[\mathit{\tau }\right]=\left|\right|{\mathbb{B}}_i\left[\mathit{\tau }\right]|{|}^2\ge 0\text{for all}\mathit{\tau }\in \text{Sym}(d).$$

(2.26)

In addition to the damage functions, we assume a hardening-type damage-surface potential of power-law type

$$h(q)={\displaystyle \sum _{i=1}^M\frac{G_i}{m_i+1}}{q}_i^{m_i+1},m_i>0,$$

(2.27)

involving a positive, dimension-free power-law exponent $m_i$ and a positive hardening parameter $G_i$ with dimensions of stress. Thus, according to (2.22), the damage-driving forces compute as

$${\beta }_i=-G_i{q}_i^{m_i},i=1,\dots ,M.$$

(2.28)

In view of the force potential ${\Phi }^{*}$ (2.20) and the driving forces $\mathbb{T}$ (2.18) and $\beta$ (2.28), there is an elastic domain in the (extended) stress space, described by the conditions

$$\left|\right|{\mathbb{B}}_i\left[\mathit{\sigma }\right]|{|}^2\le {\sigma }_{0,i}^2+G_i{H}_i{q}_i^{m_i},i=1,\dots ,M,$$

(2.29)

where no damage occurs. As defined in equations $(2.21)$ , the evolutions of the compliance and the damage-surface variables are governed by

$$\dot{\mathbb{S}}=2{\displaystyle \sum _{i=1}^M{\dot{\mu }}_i}{\mathbb{B}}_i^2\text{ and }{\dot{q}}_i={\dot{\mu }}_i{H}_i,i=1,\dots ,M,$$

(2.30)

in case of an active damage system at index i – otherwise, ${\dot{\mu }}_i=0$ holds.

Several simplifications are in order. First, notice that the parameters G_i and H_i only enter (2.29) as the product $G_i{H}_i$ . As we may redefine $\widetilde{G_i}=\widetilde{H_i}=\sqrt{G_i{H}_i}$ without changing the elastic domain (2.29), we assume $G_i=H_i$ . Secondly, we may eliminate the consistency parameter from the evolution of the compliance $(2.30)$ and integrate to get

$$\mathbb{S}(t)={\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i(t)}{H_i}}{\mathbb{B}}_i^2,$$

(2.31)

where ${\mathbb{S}}_0=\mathbb{S}(0)$ is the initial compliance. Thirdly, in three spatial dimensions d = 3, the compliance tensor $\mathbb{S}$ is described by 21 independent parameters. The latter formula (2.31) permits us to express the current compliance tensor $\mathbb{S}$ in terms of the internal variables q. Thus, with an eye towards an efficient implementation, we may a posteriori eliminate the compliance tensor $\mathbb{S}$ from the model. Furthermore, as stated above, notice that the equations (2.30) permit us to eliminate the parameters $\dot{\mu }$ completely. Last but not least, in view of the elastic domain (2.29), we may work with the damage-activation functions f_i

$$f_i:\text{Sym}(d)\times \mathbb{R}\to \mathbb{R},(\mathit{\sigma },q_i)\hspace{0.17em}\mapsto \hspace{0.17em}\left|\right|{\mathbb{B}}_i\left[\mathit{\sigma }\right]|{|}^2-{\sigma }_{0,i}^2-H_i^2{q}_i^{m_i},i=1,\dots ,M,$$

(2.32)

instead of the original functions ${\varphi }_i$ (2.25). For the convenience of the reader, we summarize the key aspects of the model in the following box.

Summary of compliance-based convex damage model (primal formulation)

Input Initial compliance tensor ${\mathbb{S}}_0$ , extraction tensors ${\mathbb{B}}_i$ , hardening moduli $H_i>0$ , damage thresholds ${\sigma }_{0,i}>0$ , power-law exponents $m_i>0 (i=1,\dots ,M)$ .

Evolution equations For given strain path $\mathit{\epsilon }:[0,T]\to \text{Sym}(d)$ , find damage-hardening variables $q:[0,T]\to {\mathbb{R}}^M$ and a stress path $\mathit{\sigma }:[0,T]\to \text{Sym}(d)$ , s. t.

$$f_i(\mathit{\sigma },q_i)\le 0,{\dot{q}}_i\ge 0,{\dot{q}}_i{f}_i(\mathit{\sigma },q_i)=0,i=1,\dots ,M,$$

(2.33)

holds, with initial conditions $q(0)=0$ , and where

$$f_i(\mathit{\sigma },q_i)=\left|\right|{\mathbb{B}}_i\left[\mathit{\sigma }\right]|{|}^2-{\sigma }_{0,i}^2-H_i^2{q}_i^{m_i}\text{ and }\mathit{\epsilon }=\left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i}{H_i}}{\mathbb{B}}_i^2\right)\left[\mathit{\sigma }\right].$$

(2.34)

Computational predictor-corrector framework

In this section, upon an implicit Euler discretization in time, we discuss a predictor-corrector solution strategy for the model introduced in the previous section in strict analogy to associative elastoplasticity, see Chapter 2 in Simo and Hughes (1998). Suppose that a number of discrete time steps $0=t_0<t_1<\dots <t_{N-1}<t_N=T$ is given, together with prescribed strain tensors ${\mathit{\epsilon }}_0,{\mathit{\epsilon }}_1,\dots {\mathit{\epsilon }}_N$ , an initial compliance tensor ${\mathbb{S}}_0$ and the initial damage-hardening variable $q_0\equiv 0\in {\mathbb{R}}^M$ . For any $n=0,\dots ,N-1$ , dropping the subscript n + 1 for simplicity of notation, we seek $(\mathit{\sigma },q)\in \text{Sym}(d)\times {\mathbb{R}}^M$ solving the system of equations

$$\begin{array}{l}\mathit{\epsilon }=\left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i}{H_i}}{\mathbb{B}}_i^2\right)\left[\mathit{\sigma }\right]\hfill \\ f_i(\mathit{\sigma },q_i)\le 0,q_i-q_{i,n}\ge 0,(q_i-q_{i,n})f_i(\mathit{\sigma },q_i)=0,i=1,\dots ,M,\hfill \end{array}$$

(2.35)

with the damage functions f_i ${(2.34)}_1$ . With a computational resolution in mind, we rewrite the system (2.35) in terms of active sets. For any $(\mathit{\sigma },q)\in \text{Sym}(d)\times {\mathbb{R}}^M$ , the active set $\mathcal{A}(\mathit{\sigma },q)$ is defined as

$$\mathcal{A}(\mathit{\sigma },q)=\left\{i\in \left\{1,2,\dots ,M\right\}|\hspace{0.17em}{f}_i(\mathit{\sigma },q_i)\ge 0\right\},$$

(2.36)

collecting all indices of inequality constraints that are either violated or satisfied exactly. Then, as a consequence of the complementarity condition in the system (2.35), $(\mathit{\sigma },q)\in \text{Sym}(d)\times {\mathbb{R}}^M$ solves the system (2.35) precisely if it satisfies $q_i\ge q_{i,n} (i=1,\dots ,M)$ and solves

$$\begin{array}{l}\left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i}{H_i}}{\mathbb{B}}_i^2\right)\left[\mathit{\sigma }\right]=\mathit{\epsilon }\hfill \\ f_i(\mathit{\sigma },q_i)=0\text{for all}i\in \mathcal{A}(\mathit{\sigma },q).\hfill \end{array}$$

(2.37)

We solve the latter problem by an active set strategy (Bergounioux et al., 1999, 2000), i.e., by solving the system (2.37) with a Newton method, updating the currently active set at each Newton iteration and accounting for the constraints $q_i\ge q_{i,n} (i=1,\dots ,M)$ via backtracking. The details comprise Alg. 1, where $\gamma \in (0,1)$ is a backtracking factor. We use a backtracking factor of $\gamma =0.9$ in our presented examples.

As long as the damage constraints are linearly independent, due to the established connections of active set strategies to semi-smooth Newton methods, see Hintermüller et al. (2002), a locally superlinear convergence behavior can be expected. A schematic of the predictor-corrector strategy is shown in Figure 3 with

$$\text{residual}(\mathit{\sigma },q)=\frac{\sqrt{{\Vert \mathit{\sigma }-{\left({\mathbb{S}}_0+2{\displaystyle {\sum }_{i=1}^M{q}_i}{\mathbb{B}}_i^2/H_i\right)}^{-1}\left[\mathit{\epsilon }\right]\Vert }^2+{\displaystyle {\sum }_{i=1}^M\max }\{0,f_i(\mathit{\sigma },q_i)\}}}{\left|\right|\mathit{\sigma }\left|\right|}$$

(2.38)

for measuring convergence. Whenever the trial stress fails to be contained in the elastic region, an iterative process is initiated which ensures that the final stress state again lies on the boundary of the elastic domain. For the latter, both the elastic region may grow – as a result of the damage-hardening – and the stress may decrease – owing to increasing compliance.

For solving problem (2.37), we assemble the Newton system for the active set $\mathcal{A}(\mathit{\sigma },q)$

$$\begin{array}{c}\left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i}{H_i}}{\mathbb{B}}_i^2\right)\left[△\mathit{\sigma }\right]+{\displaystyle \sum _{i=1}^M\frac{2}{H_i}}{\mathbb{B}}_i^2\left[\mathit{\sigma }\right]△q_i=\mathit{\epsilon }-\left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M\frac{q_i}{H_i}}{\mathbb{B}}_i^2\right)\left[\mathit{\sigma }\right]\\ \left(\frac{2}{H_i}{\mathbb{B}}_i^2\left[\mathit{\sigma }\right]\right)·△\mathit{\sigma }-m_i{H}_i{q}_i^{m_i-1}△q_i=-\frac{1}{H_i}\left(\left|\right|{\mathbb{B}}_i\left[\mathit{\sigma }\right]|{|}^2-{\sigma }_{0,i}^2-H_i^2{q}_i^{m_i}\right)\end{array}$$

(2.39)

for all $i\in \mathcal{A}(\mathit{\sigma },q)$ , where we divided the second line by $H_i$ to ensure a symmetric Newton system.

Figure 3.

Evolution of the elastic region upon loading within a predictor-corrector framework.

Algorithm 1 Predictor-corrector strategy $\left(\mathit{\epsilon },q_{i,n}\right)$ with model parameters $\left({\mathbb{S}}_0,{\mathbb{B}}_i,H_i,{\sigma }_{0,i},m_i\right)$ and algorithm parameters $\left(\text{maxit},\text{tol},\gamma \right)$

1: Elastic predictor

2: $q\leftarrow q_n$

3: $\mathit{\sigma }\leftarrow \left({\mathbb{S}}_0+2{\displaystyle \sum _{i=1}^M{q}_i}{\mathbb{B}}_i^2/H_i\right)\left[\mathit{\epsilon }\right]$

4: if all $f_i\le 0$ then

5:  no damage evolution, elastic predictor step correct

6: else

7:  Damage corrector

8:   $k\leftarrow 1$ ▷ Iteration counter

9:  Update residual (2.38)

10:  while $k<\text{maxit}\mathbf{and}\text{residual}>\text{tol}$ do

11:    $\mathcal{A}\leftarrow \mathcal{A}(\mathit{\sigma },q)$

12:   assemble and solve the Newton system (2.39) for $(△\mathit{\sigma },△q)$ ▷ $△q_i:=0$ for $i\ne \mathcal{A}$

13:    $s\leftarrow 1$   ▷Full step size

14:    $(\mathit{\sigma },q)\leftarrow (\mathit{\sigma }+s△\mathit{\sigma },q+s\Delta q)$

15:    $j\leftarrow 0$ ▷ Counts backtracking steps

16:    ${\text{residual}}_{\text{old}}\leftarrow \text{residual}$

17:   Update residual (2.38)

18:   while $\text{residual}>{\text{residual}}_{\text{old}}\mathbf{or}{q}_i<q_{i,n}\text{ for some }i$ do▷ Backtracking, typically $\gamma =0.9$

19:      $(\mathit{\sigma },q)\leftarrow (\mathit{\sigma }+(\gamma s-s)△\mathit{\sigma },q+(\gamma s-s)△q)$

20:      $s\leftarrow \gamma s$ ▷ Reduce current step size

21:     Update residual (2.38)

22:      $j\leftarrow j+1$

23:   end while

24:    $k\leftarrow k+1$

25:  end while

26: end if

27: compute ${ℂ}_{\text{algo}}$

28: Output

29: return $\mathit{\sigma }$ , q, ${ℂ}_{\text{algo}}$

Damage models with Puck-type extraction tensors

Basic idea

Puck (Knops, 2008; Puck and Schürmann, 2002) introduced strength-estimation models for composites reinforced by continuous fibers based on specfic failure scenarios that are commonly observed in post-critical investigations of failed specimens. For the current article at hand, we will use these so-called Puck cases as primary drivers of the anisotropic damage evolution presented in the previous section. More precisely, we will investigate the Puck cases individually, and determine proper extraction tensors ( ${\mathbb{B}}_{\text{I}}$ - ${\mathbb{B}}_{\text{IV}}$ ). These Puck-type extraction tensors are motivated by the stress and corresponding damage states present in the fiber bundle mesostructure of Sheet Molding Compound (SMC) composites (Dumont et al., 2007; Görthofer et al., 2019).

We introduce a local Cartesian coordinate system $\{{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3\}$ , s. t. the fibers are aligned to the ${\mathit{e}}_1$ -direction, see Figure 4(a). Then, the stress state $\mathit{\sigma }$ may be decomposed into blocks

$$\mathit{\sigma }=\left(\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{12}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{13}& {\sigma }_{23}& {\sigma }_{33}\end{array}\right),$$

(3.1)

where ${\sigma }_{11}$ is the stress in fiber direction, the lower right block describes the stresses in the plane orthogonal to the fiber direction, and $({\sigma }_{12},{\sigma }_{13})$ collects the remaining shear stresses. Adopting ideas of Puck (Knops, 2008; Puck and Schürmann, 2002), we distinguish four basic cases which drive the damage evolution in a fiber bundle, for instance.

(I)	Normal loading in fiber direction	${\sigma }_{11}$	$▷$ Figure 4(a) and (b)
(II)	Normal loading perpendicular to fiber direction	${\sigma }_{22},{\sigma }_{33}$	$▷$ Figure 4(c) and (d)
(III)	Shear loading perpendicular to fiber direction	${\sigma }_{23}$	$▷$ Figure 4(e)
(IV)	Shear loading in fiber direction	${\sigma }_{12},{\sigma }_{13}$	$▷$ Figure 4(f)

Figure 4.

Regions of major damage (blue) resulting from different loading scenarios in a cell with aligned fibers (dark green). (a) Extension in fiber direction. (b) Compression in fiber direction. (c) Extension ⊥ to fiber direction. (d) Compression ⊥ to fiber direction. (e) Shearing ⊥ to fiber direction. (f) Shearing in fiber direction.

The loading scenarios shown in Figure 4 are only examples, e.g., loadings perpendicular to the fiber direction need not necessarily follow direction ${\mathit{e}}_2$ . Instead, any other direction in the ${\mathit{e}}_2$ - ${\mathit{e}}_3$ -plane could be used, as well. Nevertheless, we may regard a general loading scenario as a superposition of the four introduced cases. In the following sections, we will derive appropriate extraction tensors ( ${\mathbb{B}}_{\text{I}}$ - ${\mathbb{B}}_{\text{IV}}$ ) corresponding to each of the four presented cases based on averaged stress conditions. The presented model cannot distinguish between tensile and compressive loading, as the driving force $\mathbb{T}$ (2.18) is quadratic in the stress $\mathit{\sigma }$ . Consequently, the six sketched loading scenarios in Figure 4 reduce to the mentioned four cases, as the scenarios Figure 4(a) and (b), as well as Figure 4(c) and (d) coincide for our model.

Case I: Normal loading in fiber direction

The first damage case is governed by loading in fiber or bundle direction, respectively, and thus solely concerns the stress ${\sigma }_{11}$ . For fiber direction ${\mathit{e}}_1$ , the fourth-order extraction tensor ${\mathbb{B}}_{\text{I}}$ extracting the stress in bundle direction ${\sigma }_{11}$ from an arbitrary stress state $\mathit{\sigma }$ is given by

$${\mathbb{B}}_{\text{I}}={\mathit{e}}_1^{\otimes 4}.$$

(3.2)

The associated damage function ${(2.34)}_1$ with case I extraction tensor ${\mathbb{B}}_{\text{I}}$ (3.2) reads

$$f_{\text{I}}(\mathit{\sigma },q)={\sigma }_{11}^2-H^2{q}^m-{\sigma }_0^2,$$

(3.3)

and will solely induce a decrease in the Young’s modulus in ${\mathit{e}}_1$ -direction. As a general note, although the damage parameters like ${\sigma }_0$ , H, m may differ for the considered cases I to IV, we do not include additional subscripts for the sake of readability.

Case II: Normal loading perpendicular to fiber direction

To quantify damaging due to normal loading in any direction $S^2\mathrm{\ni }\mathit{n}\perp {\mathit{e}}_1$ perpendicular to the fiber direction (a unit vector with $S^2=\{\mathit{x}\in {\mathbb{R}}^3|\left|\right|\mathit{x}\left|\right|=1\}$ ), for a general stress state $\mathit{\sigma }$ , we measure the average normal stress perpendicular to the fiber direction

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\mathit{n}}{(\theta )}^{\otimes 4}\left[\mathit{\sigma }\right]\text{d}\theta \text{ with }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta ).$$

(3.4)

The latter average may be represented in the form

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\mathit{n}}{(\theta )}^{\otimes 4}\left[\mathit{\sigma }\right]\text{d}\theta \stackrel{!}{=}{\mathbb{B}}_{\text{II}}\left[\mathit{\sigma }\right]$$

(3.5)

in terms of the extraction tensor ${\mathbb{B}}_{\text{II}}$

$${\mathbb{B}}_{\text{II}}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\mathit{n}}{(\theta )}^{\otimes 4}\text{d}\theta \text{ with }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta ).$$

(3.6)

Evaluating the integral explicitly, see Appendix C for details, yields

$${\mathbb{B}}_{\text{II}}=\frac{1}{4}{({\mathit{e}}_2^{\otimes 2}+{\mathit{e}}_3^{\otimes 2})}^{\otimes 2}+\frac{1}{8}{({\mathit{e}}_2^{\otimes 2}-{\mathit{e}}_3^{\otimes 2})}^{\otimes 2}+\frac{1}{2}{({\mathit{e}}_2{\otimes }_{\mathrm{S}}{\mathit{e}}_3)}^{\otimes 2}.$$

(3.7)

This extraction tensor ${\mathbb{B}}_{\text{II}}$ is identical to the fourth-order fiber-orientation tensor for a planar isotropic orientation, see Advani and Tucker (1987). The damage function ${(2.34)}_1$ involving the case II extraction tensor ${\mathbb{B}}_{\text{II}}$ (3.7) reads

$$f_{\text{II}}(\mathit{\sigma },q)=\frac{1}{32}\left(5{\sigma }_{22}^2+5{\sigma }_{33}^2+6{\sigma }_{22}{\sigma }_{33}+4{\sigma }_{23}^2\right)-H^2{q}^m-{\sigma }_0^2.$$

(3.8)

Case III: Shear loading perpendicular to fiber direction

In addition to damage caused by normal loading, we also want to account for shear-loading induced damage. For a general stress state $\mathit{\sigma }$ , the average shear stress transverse to the fiber direction is given by

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left(\mathit{n}(\theta ){\otimes }_{\mathrm{S}}\mathit{m}(\theta )\right)}^{\otimes 2}}\left[\mathit{\sigma }\right]\text{d}\theta$$

(3.9)

$$\text{ with }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta )\text{ and }\mathit{m}(\theta )\stackrel{\wedge}{=}(0,-\sin \theta ,\cos \theta ).$$

We may rewrite this expression

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left(\mathit{n}(\theta ){\otimes }_{\mathrm{S}}\mathit{m}(\theta )\right)}^{\otimes 2}}\left[\mathit{\sigma }\right]\text{d}\theta \stackrel{!}{=}{\mathbb{B}}_{\text{III}}\left[\mathit{\sigma }\right]$$

(3.10)

in terms of the extraction tensor ${\mathbb{B}}_{\text{III}}$ as

$${\mathbb{B}}_{\text{III}}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left(\mathit{n}(\theta ){\otimes }_{\mathrm{S}}\mathit{m}(\theta )\right)}^{\otimes 2}}\text{d}\theta$$

(3.11)

$$\text{ with }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta )\text{ and }\mathit{m}(\theta )\stackrel{\wedge}{=}(0,-\sin \theta ,\cos \theta ).$$

Explicitly evaluating the integral (3.11) yields

$${\mathbb{B}}_{\text{III}}=\frac{1}{8}{\left({\mathit{e}}_2^{\otimes 2}-{\mathit{e}}_3^{\otimes 2}\right)}^{\otimes 2}+\frac{1}{2}{({\mathit{e}}_2{\otimes }_{\mathrm{S}}{\mathit{e}}_3)}^{\otimes 2}.$$

(3.12)

The damage function ${(2.34)}_1$ for the case III extraction tensor ${\mathbb{B}}_{\text{III}}$ (3.12) reads

$$f_{\text{III}}(\mathit{\sigma },q)=\frac{1}{32}\left({\sigma }_{22}^2+{\sigma }_{33}^2-2{\sigma }_{22}{\sigma }_{33}+4{\sigma }_{23}^2\right)-H^2{q}^m-{\sigma }_0^2.$$

(3.13)

Comparing the extraction tensors ${\mathbb{B}}_{\text{II}}$ (3.7) and ${\mathbb{B}}_{\text{III}}$ (3.12), some similarities between these tensors become apparent. In fact, these similarities reflect the relationship between normal loadings (especially compression) and shear loadings perpendicular to the fiber direction, that are familiar from undergraduate engineering mechanics (Hibbeler, 2001), i.e., Mohr’s circle (Mohr, 1900).

Case IV: Shear loading in fiber direction

To evaluate damage induced by shearing in fiber direction ${\mathit{e}}_1$ , we evaluate

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left(\mathit{n}(\theta ){\otimes }_{\mathrm{S}}{\mathit{e}}_1\right)}^{\otimes 2}}\left[\mathit{\sigma }\right]\text{d}\theta \text{ with }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta ),$$

(3.14)

which we represent in the form

$$\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{\left(\mathit{n}(\theta ){\otimes }_{\mathrm{S}}{\mathit{e}}_1\right)}^{\otimes 2}}\left[\mathit{\sigma }\right]\text{d}\theta \stackrel{!}{=}{\mathbb{B}}_{\text{IV}}\left[\mathit{\sigma }\right]$$

(3.15)

with extraction tensor ${\mathbb{B}}_{\text{IV}}$

$${\mathbb{B}}_{\text{IV}}=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }(\mathit{n}(}\theta ){\otimes }_{\mathrm{S}}{\mathit{e}}_1{)}^{\otimes 2}\text{d}\theta \text{ and }\mathit{n}(\theta )\stackrel{\wedge}{=}(0,\cos \theta ,\sin \theta ).$$

(3.16)

The resulting extraction tensor for case IV may be expressed via

$${\mathbb{B}}_{\text{IV}}=\frac{1}{2}{({\mathit{e}}_1{\otimes }_{\mathrm{S}}{\mathit{e}}_2)}^{\otimes 2}+\frac{1}{2}{({\mathit{e}}_1{\otimes }_{\mathrm{S}}{\mathit{e}}_3)}^{\otimes 2}.$$

(3.17)

The damage function ${(2.34)}_1$ with case IV extraction tensor ${\mathbb{B}}_{\text{IV}}$ (3.17) becomes

$$f_{\text{IV}}(\mathit{\sigma },q)=\frac{1}{8}\left({\sigma }_{12}^2+{\sigma }_{13}^2\right)-H^2{q}^m-{\sigma }_0^2.$$

(3.18)

Computational investigations

Setup

We integrated the proposed damage model as a user-defined subroutine into an in-house OpenMP-parallel FFT-based computational homogenization code written in Python 3.7 with Cython extensions (Behnel et al., 2011) and FFTW (Frigo and Johnson, 2005) bindings, as described, e.g., by Schneider (2019). The balance of linear momentum was discretized on a staggered grid (Schneider et al., 2016b) and the ensuing nonlinear systems of equations were solved by a Newton-CG scheme (Gélébart and Mondon-Cancel, 2013; Kabel et al., 2014; Wicht et al., 2020). The computations ran on 6 - 12 threads on a desktop computer with 32 GB RAM and an Intel i7-8700K CPU with 6 cores and a clock rate of 3.7 GHz. The plain-weave presented at the end of this section was computed on a workstation with two AMD EPYC 7642 and 48 physical cores each, enabled SMT, and 1024 GB of DRAM.

For the following studies, we use the isotropic elastic parameters of unsaturated polyester-polyurethane hybrid (UPPH) resin and E-glass fibers, see Kehrer et al. (2018), respectively, if not specified otherwise. Furthermore, we use the damage parameters ${\sigma }_0$ , H and m listed in Table 1 as a point of departure for parameter variation and the different introduced damage cases. Due to the small-strain setting, we limit the strain axes to $5\%$ in magnitude.

Table 1.

Standard material parameters (Kehrer et al., 2018) and reference damage parameters, serving as point of departure depending on the corresponding damage case.

UPPH matrix	E-glass fibers	Damage parameters
$E_{\text{M}}=3.4 \text{GPa}$	$E_{\text{F}}=72 \text{GPa}$	${\sigma }_0\in \left[5,30\right] \text{MPa}$ , $m=1$
${\nu }_{\text{M}}=0.385$	${\nu }_{\text{F}}=0.22$	$H\in \left[30,80\right] \text{MPa}$

Numerical studies on integration-point level

Parameter study for Puck-type extraction tensor I

The first study concerns the effects of the damage parameters ${\sigma }_0$ , H and m on the stress and damage evolution. For this purpose, we investigate the model behavior for one active damage function and the Puck-type extraction tensor I only. We vary the parameters and evaluate the stress-strain curves, as well as the normalized Young’s modulus $E_{11}/E_{11}^0$ for uniaxial extension in ${\mathit{e}}_1$ -direction and a prescribed strain ${\epsilon }_{11}$ .

For the case at hand, we extract the current Young’s modulus in ${\mathit{e}}_1$ -direction from equation (2.31) by

$$E_{11}=H{E}_{11}^0/(H+2q{E}_{11}^0),$$

(4.1)

where $E_{11}^0$ stands for the initial Young’s modulus in ${\mathit{e}}_1$ -direction. Please note that $E_{11}^{-1}$ may be regarded as a component of a (fourth-order) tensor, as it is the 1111 component of the corresponding compliance $\mathbb{S}$ . For the sake of simplicity, we will use such shorthand index notation for certain Young’s moduli throughout the remainder of this work.

The influence of the damage-activation threshold ${\sigma }_0$ on the stress and normalized stiffness in ${\mathit{e}}_1$ -direction is shown in Figure 5. The higher the damage-activation threshold ${\sigma }_0$ , the later the damage evolution initiates w. r. t. the applied strain ${\epsilon }_{11}$ . In Figure 5(a), we observe damage to initiate as soon as the stress ${\sigma }_{11}$ equals the damage-activation threshold ${\sigma }_0$ , which is expected. The convex hardening nature of our model gives rise to a decreasing slope of the stress-strain curve. This slope tends to zero at infinity, but remains non-negative.

With increasing damage-activation thresholds ${\sigma }_0$ , the stress-strain behavior approaches a plateau beyond the elastic region, in which an increase of strain does not induce a further increase of stress. The reduced stiffness $E_{11}$ in ${\mathit{e}}_1$ -direction equals the slope of the stress-strain curve during unloading (which returns to the origin in our pure elastic-damage framework). Whereas the plateaus are more pronounced for higher damage-activation thresholds ${\sigma }_0$ , the increase in damage and the (normalized) stiffness reduction are less pronounced, see Figure 5(b).

Due to the thermodynamically consistent GSM framework of our proposed model, an upper bound for the damage variables is ensured, governing the asymptotic behavior of the (normalized) stiffness reduction, see Figure 5(b). Evaluating the CDI (2.2) for the considered case at hand, we find the upper bound for the damage variable w. r. t. Puck case I to be $q_{\text{I}}\le \sqrt[m]{{\sigma }_{11}^2/H^2}$ .

The effect of changing the hardening parameter H on the stress and damage evolution is shown in Figure 6. As the damage-activation threshold ${\sigma }_0$ remains unchanged for this study, the elastic regions ( ${\sigma }_{11}$ below 30 MPa) coincide for all stress-strain curves, see Figure 6(a).

In the damage region, the slope of the stress-strain curve increases with the hardening parameter H. Indeed, the hardening parameter H describes the hardening capacity of the model. For $H\to 0$ , the slope tends to zero and approaches the plateau already observed in Figure 5(a). For $H\to \infty$ , an active damage function f (2.34)₁ will become inactive for small values of the damage variable $q\to 0$ , resulting in a damage region that can hardly be distinguished from a purely elastic material behavior.

Independent of the hardening parameter, damage evolution initiates at the same strain level (and therefore the same stress level), see Figure 6(b). The stiffness reduction is inversely proportional to the hardening parameter H.

In Figure 7, we fix the damage-activation threshold ${\sigma }_0$ as well as the hardening parameter H and vary the power-law exponent m. In contrast to the influence of the hardening parameter H (see Figure 6), the slope of the stress-strain curve is inversely proportional to the exponent m in the damage region, see Figure 7(a).

For increasing exponents m, the stress-strain curves approach the plateau-like behavior. For small values of m, after exceeding the damage-activation threshold ${\sigma }_0$ , the stress-strain curves remain approximately linear and only develop a significant amount of damage for higher loading levels.

Figure 7(b) shows the damage evolution to be inversely proportional to the exponent m, leading to a lower remaining (normalized) stiffness component $E_{11}$ for higher exponents m.

For representing the stiffness tensor $ℂ={\mathbb{S}}^{-1}$ graphically, we use the Young’s modulus surface (YMS) plots introduced by Böhlke and Brüggemann (2001), i.e., for fixed compliance tensor $\mathbb{S}\in {\mathcal{S}}_d$ , the image of the nonlinear mapping

$$S^2\to {\mathbb{R}}^3,\mathit{n}\hspace{0.17em}\mapsto \hspace{0.17em}E(\mathbb{S},\mathit{n})\mathit{n},$$

(4.2)

where the Young’s modulus $E(\mathbb{S},\mathit{n})$ in direction $\mathit{n}$ is determined by

$$E(\mathbb{S},\mathit{n})=\frac{1}{\left(\mathit{n}\otimes \mathit{n}\right)·\mathbb{S}\left[\mathit{n}\otimes \mathit{n}\right]}.$$

(4.3)

Figure 5.

Varying the initial stress ${\sigma }_0$ for the proposed model. (a) Stress-strain curve. (b) Normalized Young’s modulus vs. applied strain.

Figure 6.

Varying the hardening parameter $H$ for the proposed model. (a) Stress-strain curve. (b) Normalized Young’s modulus vs. applied strain.

Figure 7.

Varying the power-law exponent $m$ for the proposed model. (a) Stress-strain curve. (b) Normalized Young’s modulus vs. applied strain.

Asymmetry properties of the stiffness tensor $ℂ$ become apparent in the corresponding YMS plot.

Examples of such YMS plots are shown in Figure 8. The initially isotropic stiffness tensor with UPPH material parameters (see Table 1) has a spherical shape, as shown in Figure 8(a). As discussed above, the induced damage model based on Puck case I leads to a reduction of the stiffness $E_{11}$ in the ${\mathit{e}}_1$ -direction. The YMS plot is contracted in the direction of loading, whereas the directional Young’s moduli in the orthogonal plane remain unaffected, see Figure 8(b).

Figure 8.

YMS plots (see Böhlke and Brüggemann, 2001) showing the reduction of the stiffness tensor based on Puck case I. (a) Initial isotropic state. (b) Final anisotropic state.

Stiffness reduction for different Puck-type extraction tensors

In the following, we will discuss possible damage evolutions and corresponding stiffness reductions based on Puck-type extraction tensors II, III and IV. For these cases, the influence of the damage-hardening parameters ${\sigma }_0$ , H and m is similar to case I, which we discussed in the previous section. We consider an initially isotropic stiffness tensor with UPPH material parameters. The corresponding YMS plot is shown in Figure 8(a). Specific loadings will evoke a damage evolution due to the Puck cases II, III and IV. We apply a normal strain ${\epsilon }_{22}$ for case II, a shear strain ${\epsilon }_{23}$ for case III and a shear strain ${\epsilon }_{13}$ for case IV, forcing the respective complementary stress components to zero. The resulting YMS plots are shown in Figure 9. For all three cases the Young’s modulus in ${\mathit{e}}_1$ -direction remains unchanged. As shortly discussed in the previous chapter, Puck cases II and III are interlinked due to similar effects of averaged normal loadings and shear loadings perpendicular to the fiber direction. Inspecting Figure 9(a) and (b), corresponding to Puck cases II and III, we observe a reduction of the Young’s moduli within the ${\mathit{e}}_2$ - ${\mathit{e}}_3$ -plane for both cases, but with different characteristics. For the loading scenarios considered here, the stiffness reduction in directions ${\mathit{e}}_2$ and ${\mathit{e}}_3$ is more pronounced for Puck case II compared to Puck case III. Figure 9(c) shows that damage based on Puck case IV does not affect the Young’s moduli in the ${\mathit{e}}_2$ - ${\mathit{e}}_3$ -plane. Young’s moduli in the ${\mathit{e}}_1$ - ${\mathit{e}}_3$ -plane and ${\mathit{e}}_1$ - ${\mathit{e}}_2$ -plane are reduced equally, leading to a transversely isotropic stiffness with the fiber direction ${\mathit{e}}_1$ as the axis of symmetry.

Figure 9.

YMS plots (see Böhlke and Brüggemann, 2001) illustrating the reduction of the stiffness tensor based on the Puck cases II, III and IV. (a) Case II. (b) Case III. (c) Case IV.

Non-monotonic loading

To show the capabilities of our model in general, we perform loading-unloading experiments for different loading directions in a successive fashion. To mimic uniaxial normal loadings and corresponding shear loadings, we subsequently apply normal and shear strains ${\epsilon }_{11},{\epsilon }_{22},{\epsilon }_{33},{\epsilon }_{23},{\epsilon }_{13} \text{ and } {\epsilon }_{12}$ , each with mixed boundary conditions permitting solely the corresponding stresses ${\sigma }_{11},{\sigma }_{22},{\sigma }_{33},{\sigma }_{23},{\sigma }_{13} \text{ and } {\sigma }_{12}$ , to be non-zero, see Figure 10(a) and (b). Lateral contraction is permitted. Between each of these six loading steps we unload to zero strain and stress. The resulting evolution of the stress components is plotted in Figure 10(b). We see a linear elastic region for each individual loading step and a damage region for all but the ${\sigma }_{33}$ and ${\sigma }_{12}$ cases. On these occasions, the threshold for damage initiation is not triggered. As all normal loadings in the ${\mathit{e}}_2\text{-}{\mathit{e}}_3$ -plane evoke damage due to case II, ${\sigma }_{33}$ cannot induce damage beyond the level previously induced. The stress-strain curves for the shear stresses ${\sigma }_{13}$ and ${\sigma }_{12}$ exhibit a similar behavior.

Figure 10.

Complex loading history addressing each stress-tensor component separately. (a) Strain vs. time. (b) Stress vs. time. (c) Damage evolution over time.

Figure 11.

Evolution of an initially isotropic stiffness upon complex loading, see Figure 10(a), visualized via YMS plots. (a) After loading step 1 (ε₁₁). (b) After loading step 2 (ε₁₁ and ε₂₂). (c) After loading step 3 (ε₁₁, ε₂₂ and ε₃₃). (d) After loading step 4 (ε₁₁, ε₂₂, ε₃₃ and ε₂₃). (e) After loading step 5 (ε₁₁, ε₂₂, ε₃₃, ε₂₃ and ε₁₃). (f) After loading step 6 (ε₁₁, ε₂₂, ε₃₃, ε₂₃, ε₁₃ and ε₁₂).

These observations are also reflected in the damage evolution, see Figure 10(c). Due to the applied stress ${\sigma }_{22}$ , the damage variable representing case II increases. For increasing stress ${\sigma }_{33}$ , the damage variables remain unaffected up to the level ${\sigma }_{33}={\sigma }_{22}$ . The same line of argument applies to ${\sigma }_{13}$ and ${\sigma }_{12}$ . Again, we observe the connection between Puck cases II and III, i. e., each of the damage variables increases whenever a loading scenario is applied which progresses the other case. The damaged stiffness tensors corresponding to each loading step are visualized via their YMS plots in Figure 11. The different colors of the plots represent the different loading steps shown in Figure 10. Note that the ranges of the axes are adjusted accordingly and therefore vary from plot to plot as the damage increases from step to step.

Comparing Figure 11(b) and (c), as well as Figure 11(e) and (f), we see that, in accordance with Figure 10(c), no further damage is induced between these loading steps. The presented YMS plots demonstrate the capability of our model to evolve the stiffness tensor in a complex and anisotropic way. The model is capable of handling any initial stiffness, not restricted by a specific symmetry class, i.e., transversely isotropic or orthotropic. Furthermore, the stiffness tensor may also develop anisotropy - within the permissible set ${\mathcal{S}}_d$ (2.8) - as a result of a damaging process, owing to the introduced damage functions.

Multiaxial loading with increasing loading level

In a study combining the six loading steps from the previous section to one superposed loading case, we apply a three-dimensional strain state s. t. all strains are active simultaneously. To investigate the model behavior, the predicted stiffness degradation, as well as the evolution of the damage variables more closely, we gradually increase the strain levels in five steps from 0% to 5% via 1% increments with intermediate unloading, see Figure 12(a). As in the previous section, we analyze the stress and the damage evolution, as well as the stiffness reduction. Figure 12(b) shows the evolution of the individual stress components during the combined loading. We see that after each loading-unloading step the level of damage increases. This is also reflected in the evolution of the damage variables, see Figure 12(c). After each loading-unloading step, the damage variables continue to increase whenever the maximum stresses of the previous step are exceeded. The YMS plots corresponding to $2\%$ and $5\%$ loading, as well as the initially isotropic YMS plot, are shown in Figure 13. The stiffness tensor gradually reduces based on all four introduced Puck cases simultaneously. Apparently, the multiaxial loading evokes a stiffness reduction in all directions, as we observe a superposition of the individual loading scenarios investigated in the previous sections. Please note, that, similar to these previous section, the ranges of the axes are adjusted accordingly and therefore vary from plot to plot.

Figure 12.

Step wise increase of multiaxial loading to evoke all stresses and damage functions simultaneously. (a) Strain vs. time. (b) Stress vs. time. (c) Damage evolution over time.

Figure 13.

Evolution of an initially isotropic stiffness during multiaxial loading steps, visualized via YMS plots. (a) Strain 0%. (b) Strain 2%. (c) Strain 5%.

Cyclic tensile loading with increasing loading level

We conduct an analysis of the model response upon cyclic loading via uniaxial extension. We successively apply a normal strain ${\epsilon }_{11}$ to induce cyclic uniaxial loading in the ${\mathit{e}}_1$ -direction. We apply the strain ${\epsilon }_{11}$ in five cycles from $0\% \text{to} 5\%$ with an increasing magnitude of $1\%$ per cycle. For this analysis, we restrict to Puck case I. The resulting stress-strain curves and damage-strain curves are shown in Figure 14.

Figure 14.

Cyclic tensile loading with increasing loading level for Puck case I. (a) Stress-strain curve. (b) Damage-strain curve.

Upon loading, the material behaves linear elastically until a specific critical stress threshold or the maximum stress of the previous cycle is reached, see Figure 14(a). During the in-between unloading to ${\epsilon }_{11}=0\%$ and reloading, the damage variables do not evolve further. Besides, the pure damaging character of the model is highlighted, as no remaining residual strains occur. The step wise evolution of the damage variable $q$ is shown in Figure 14(b).

Note that the presented model is capable to predict damage onset upon both, tensile and compressive loading. Due to the definition of the damage functions the quadratic nature of the driving force $\mathbb{T}$ (2.18), the damage evolution is driven in similar fashions by both, tensile and compressive stresses. Considering a single damage variable $q$ for, both, the tensile and the compressive regime, would lead to a combined damage evolution. Accounting for tension-compression asymmetry requires an extension of the model at hand, see the conclusion.

Model response for a continuous-fiber microstructure

Application of separate loading cases with Puck-type extraction tensors

After discussing the model response for homogeneous stress states, we account for heterogeneous stress states in two ways to show the basic feasibilities of our damage model. First, we shall investigate a microstructure with a continuous fiber reinforcement. In the next section, we will turn our attention to a mesoscale simulation.

To account for damage evolution in the matrix, we introduce two extraction tensor corresponding to the spherical and deviatoric projectors of fourth-order

$${\mathbb{B}}_{\text{sph}}={\mathbb{P}}_1=\frac{1}{3}\mathit{I}\otimes \mathit{I},{\mathbb{B}}_{\text{dev}}={\mathbb{P}}_2={\mathbb{I}}^{\text{S}}-{\mathbb{P}}_1,$$

(4.4)

that allow describing a damage evolution in response to dilatation and distortion. We use damage-activation functions based on these two extraction tensors (4.4) and corresponding damage parameters ${\sigma }_0=30 \text{MPa},H=130 \text{MPa}$ and $m=1$ for both cases. Furthermore, the matrix of the fiber-reinforced microstructure is endowed with the elastic properties of UPPH, as defined in Table 1. The fibers are modeled in a purely elastic fashion using the elastic moduli of E-glass, see Table 1. The continuous-fiber reinforced microstructure is geometrically modeled by $50$ inclusions with a diameter of 13 µm and a total volume fraction of about $40\%$ . We generated the microstructure by the adaptive shrinking-cell algorithm of Torquato and Jiao (2010). The setup represents aligned fibers in a UPPH matrix as present in SMC composite bundles (Dumont et al., 2007; Kim et al., 2011; Meyer et al., 2020), where the fiber direction ${\mathit{e}}_1$ coincides with the primary direction of the Puck cases.

In three different loading scenarios, we apply three different macroscopic strains via mixed boundary conditions, see Kabel et al. (2016) for details. For each scenario, we analyze the induced damage fields of the associated variables $q_{\text{sph}}$ and $q_{\text{dev}}$ . In scenario 1, we apply the macroscopic normal strain ${\overline{\epsilon }}_{22}$ perpendicular to the fiber direction (in horizontal direction). In scenarios 2 and 3, we apply macroscopic shear strains ${\overline{\epsilon }}_{12}$ and ${\overline{\epsilon }}_{23}$ , in fiber direction and transverse to the fiber direction, respectively.

The average runtime for a resolution of $256\times 256$ pixels and $50$ time steps was about $100 \text{s}$ on 12 threads. An accompanying resolution study is discussed in the following section. Figure 15 shows the damage fields for $q_{\text{sph}} \text{ and } q_{\text{dev}}$ for the introduced cases and corresponding to the different loading scenarios. Figure 15(a) shows that normal loading perpendicular to the fiber direction leads to a dilatation-triggered damage evolution in the respective direction, as a result of stress concentrations at the inclusion boundaries. As a consequence of the complexity of the induced stress state, damage due to distortion is initiated, as well, see Figure 15(a). In general, the stress level is higher for regions with more closely packed inclusions, inducing a significantly higher level of damage in those regions. Damage initiates at the inclusion boundaries and evolves in the loading direction, deflected by other inclusions. Shear loading in fiber direction leads to an associated damage evolution due to distortion, as shown in Figure 15(f). As the applied shear is oriented in fiber direction, the deformation is not hindered by these fibers and spherical stresses do not occur. Hence, damage due to dilatation does not evolve, see Figure 15(b). Applying a macroscopic shear perpendicular to fiber direction evokes both, the evolution of damage due to spherical stresses, see Figure 15(c), as well as deviatoric stresses, see Figure 15(g).

Figure 15.

Model response for a continuous-fiber reinforced microstructure based on spherical and deviatoric damage and three different loading cases. (a) Loading ${\overline{\epsilon }}_{22},q_{\text{sph}}$ . (b) Loading ${\overline{\epsilon }}_{12},q_{\text{sph}}$ . (c) Loading ${\overline{\epsilon }}_{23},q_{\text{sph}}$ . (d) Loading ${\overline{\epsilon }}_{22},q_{\text{dev}}$ . (e) Loading ${\overline{\epsilon }}_{12},q_{\text{dev}}$ . (f) Loading ${\overline{\epsilon }}_{23},q_{\text{dev}}$ .

Resolution study

The presented resolution study demonstrates that the proposed damage model leads to mesh-independent results even without gradient enhancement. This does not come by surprise, as we specifically designed such a hardening-type damage model. Still, even in the case of hardening, a resolution study is imperative to ensure mesh-independent results. In particular, we will justify the resolution employed in the previous section.

Figure 16 shows a continuous-fiber reinforced microstructure with the same properties as for Figure 15. We vary the resolution from $64\times 64$ to $1024\times 1024$ pixels. Similar to scenario 1, see Figure 15(a), we apply a macroscopic strain ${\overline{\epsilon }}_{22}$ for all resolutions, so that both damage cases, i.e., dilatation and distortion, are being activated. The strain is successively increased from $0\%$ to $5\%$ within $50$ equidistant loading steps.

Figure 16.

Model response for a continuous-fiber reinforced microstructure evaluated at five different resolutions. (a) 64². (b) 128². (c) 256². (d) 512². (e) 1024².

The resulting distribution of the predominant damage variable $q_{\text{sph}}$ for damage due to dilatation is shown in Figure 16. We observe that areas of low and high damage level are captured also for low resolution, but there are slight differences in the achieved damage level. Also, as expected, localization behavior is not evident. To get a more qualitative insight, the macroscopic stress-strain curves are shown in Figure 17(a). For a resolution of $64\times 64$ pixels, the computed stresses are overestimated. For higher resolutions with $128\times 128$ to $1024\times 1024$ pixels, the differences are small. Investigating the relative deviations $\left({\overline{\sigma }}_{22}^{\text{reso}}-{\overline{\sigma }}_{22}^{1024}\right)/{\overline{\sigma }}_{22}^{1024}$ of the computed effective stress ${\overline{\sigma }}_{22}$ relative to the stress at a resolution of $1024\times 1024$ , see Figure 17(b), we observe that, for resolutions with $256\times 256$ pixels and higher, the deviations are below $1\%$ . The iteration counts and timings are collected in Table 2. The total outer iterations (including Newton and CG iterations) for all loading steps vary in a narrow window around approximately $2950$ for all resolutions considered. To give a comparable standard for the inner (material) iterations, we compute the average number of inner iterations over all voxels and subsequently take the maximum value over all loading steps. A value of about $2.8$ inner iterations per voxel, irrespective of the resolution, indicates quadratic convergence of the Newton method. Both, the overall timing for computing all inner iterations, as well as the total timing, increase roughly with the degrees of freedom.

Figure 17.

Resolution study for the continuous-fiber reinforced microstructure. (a) Stress-strain curve. (b) Relative stress deviation w. r. t. ${\overline{\sigma }}_{22} \text{ for } {1024}^2$ pixels.

Table 2.

Iterations and timings for the conducted resolution study.

Resolution	Iterations		Timings
	Total outer	Max. average inner	$t_{\text{inner}}$	$t_{\text{total}}$
${64}^2$	$2953$	$2.774$	$2\text{\hspace{0.17em}s}$	$9\text{\hspace{0.17em}s}$
${128}^2$	$2897$	$2.814$	$6\text{\hspace{0.17em}s}$	$26\text{\hspace{0.17em}s}$
${256}^2$	$3012$	$2.843$	$22\text{\hspace{0.17em}s}$	$108\text{\hspace{0.17em}s}$
${512}^2$	$3023$	$2.853$	$88\text{\hspace{0.17em}s}$	$524\text{\hspace{0.17em}s}$
${1024}^2$	$2987$	$2.858$	$352\text{\hspace{0.17em}s}$	$1793\text{\hspace{0.17em}s}$

Based on this resolution study, a resolution of $256\times 256$ pixels represents a compromise between an accurate prediction and a short runtime.

A Plain-weave composite under shear loading

Last but not least, we demonstrate the utility of our model framework for modeling anisotropic damage evolution in a woven fiber-reinforced composite. Simon et al. (2017) investigated the mechanical behavior of a plain-weave composite manufactured from continuous carbon fibers reinforcing an epoxy matrix resin, see Figure 18(a). The carbon fibers are aligned unidirectionally in fiber tows that are regularly interwoven. As each of these tows consists of thousands of carbon fibers, it is customary to work with a multiscale scheme that considers three different scales: the macroscopic scale is large compared to the woven unit cell, see Figure 18(a), which constitutes the mesoscale. Within the latter, the tows are considered homogeneous and anisotropic. On the microscale, the tows get resolved in terms of continuous carbon fibers in an epoxy resin.

Figure 18.

Microstructure and predicted relative reduction of the shear modulus $G_{12}$ in a plain weave composite. (a) Voxelized weave microstructure with four tows. (b) Relative reduction (4.8) of shear modulus $G_{12}$ .

The linear elastic moduli of the considered materials are listed in Table 3. These comprise the isotropic epoxy matrix and the transversely isotropic carbon fibers. The transversely isotropic engineering constants for the tows were obtained by linear elastic homogenization. Please note that the subscript “L” and “T” refer to longitudinal and transverse, respectively.

Table 3.

Elastic moduli of matrix, fibers and tows (Simon et al., 2017, Tables 1 and 2).

Constituent	Young’s modulus	Shear modulus	Poisson’s ratio
	in GPa	in GPa
Epoxy	$E$ $=$ $3$	$G$ $=$ $1.09$	$\nu$ $=$ $0.38$
Carbon fibers	$E_{\text{L}}$ $=$ $290$	$G_{\text{LT}}$ $=$ $20$	${\nu }_{\text{LT}}$ $=$ $0.2$
	$E_{\text{T}}$ $=$ $20$	$G_{\text{TT}}$ $=$ $9$	${\nu }_{\text{TT}}$ $=$ $0.11$
Tows	$E_{\text{L}}$ $=$ $144$	$G_{\text{LT}}$ $=$ $2.58$	${\nu }_{\text{LT}}$ $=$ $0.29$
	$E_{\text{T}}$ $=$ $7.84$	$G_{\text{TT}}$ $=$ $1.91$	${\nu }_{\text{TT}}$ $=$ $0.39$

Based on earlier work (Bednarcyk et al., 2015, 2014; Stier et al., 2015), Simon et al. (2017) presented a regularized orthotropic continuum damage-model based on the framework developed by Barbero and co-workers (Barbero and Lonetti, 2002; Lonetti et al., 2003, 2004) and concisely summarized in his book (Barbero, 2013). More precisely, their strategy takes the orthotropic engineering constants as the point of departure, and models their degradation on an individual basis in terms of associated scalar damage variables. Based on the associated driving forces, damage surfaces are defined, together with appropriate kinetic laws.

Expressing the dependence of the stiffness tensor on the orthotropic engineering constants is most easily realized in terms of the compliance tensor, the approach of Simon et al. (2017) appears superficially similar to our approach. However, we do not fix the damage variables a priori. Rather, they emerge naturally in our framework based on the chosen extraction tensors and damage-activation functions.

In this paragraph, we demonstrate that our model is capable of reproducing the damage behavior upon quasi-static loading of the weave composite. The protocol we present is straightforward and proceeds step by step. As a first step, we introduce a number of extraction tensors which capture elementary damage cases evoked by pure normal- and shear-loading scenarios. The tensors extract the associated normal and shear stress components from the applied stress state $\mathit{\sigma }$ . These extraction tensors are uncoupled. Hence, damage in a certain direction is solely driven by the associated loading case, e. g., normal damage in ${\mathit{e}}_1$ -direction due to normal loading in ${\mathit{e}}_1$ -direction,

$${\mathbb{B}}_{11}={\mathit{e}}_1^{\otimes 4},{\mathbb{B}}_{22}={\mathit{e}}_2^{\otimes 4},{\mathbb{B}}_{33}={\mathit{e}}_3^{\otimes 4},$$

(4.5)

$${\mathbb{B}}_{23}={\left({\mathit{e}}_2{\otimes }_{\mathrm{S}}{\mathit{e}}_3\right)}^{\otimes 2},{\mathbb{B}}_{13}={\left({\mathit{e}}_1{\otimes }_{\mathrm{S}}{\mathit{e}}_3\right)}^{\otimes 2},{\mathbb{B}}_{12}={\left({\mathit{e}}_1{\otimes }_{\mathrm{S}}{\mathit{e}}_2\right)}^{\otimes 2}.$$

(4.6)

Combining suitable damage-activation functions based on these extraction tensors permit modeling a wide range of damage-evolution predictions. In particular, they enable us to describe the stiffness reduction for the scenario considered by Simon et al. (2017).

We first capture the damage evolution in a neat epoxy sample under non-monotonic uniaxial loading and choose an extraction tensor of type (4.5) ${\hspace{0.17em}}_1$ . Subsequently, we account for the damage onset due to shear loading by using a second damage-activation function in combination with an extraction tensor of type (4.6) ${\hspace{0.17em}}_3$ . The identified parameters for the epoxy damage-model are summarized in Table 4, which were chosen to reproduce the results of Simon et al. (2017) best.

Table 4.

Extraction tensors and identified damage parameters to capture the mechanical behavior of epoxy and tow.

Constituent	Extraction tensor	${\sigma }_0$ in MPa	$H$ in MPa	$m$
Epoxy	(4.5) ${\hspace{0.17em}}_1$	$0.5$	$180$	$0.47$
	(4.6) ${\hspace{0.17em}}_3$	$101$	$41$	$0.97$
Tow	(4.6) ${\hspace{0.17em}}_1$	$8$	$275$	$0.3$
	(4.6) ${\hspace{0.17em}}_3$	$10$	$245$	$0.3$
	(4.5) ${\hspace{0.17em}}_2$	$79$	$30$	$0.9$
	(4.7)	$200$	$630$	$0.4$

Furthermore, we employ a number of damage-activation functions and suitable extraction tensors to capture the damage evolution in the fiber tows. We fix the longitudinal tow direction to correspond to the local ${\mathit{e}}_1$ -direction. The response to shear loading in longitudinal and transverse directions is best described with extraction tensors of the forms (4.6) ${\hspace{0.17em}}_1$ and (4.6) ${\hspace{0.17em}}_3$ . As the reduction of the two orthotropic Young’s moduli in the transverse plane (and hence the associated damage evolutions) is not identical, we introduce an additional extraction tensor

$$\mathbb{B}=4{\mathit{e}}_2^{\otimes 4}+{\mathit{e}}_3^{\otimes 4},$$

(4.7)

which is supplemented by a fourth damage-activation function with an extraction tensor that drives damage in ${\mathit{e}}_2$ -direction (4.5) ${\hspace{0.17em}}_2$ only. Table 4 comprises a complete list of the identified damage parameters. The listed extraction tensors and damage parameters at hand allow us to reproduce the structural behavior of, both, the neat epoxy and the tows, the latter in terms of stress-strain curves and the reduction of the orthotropic engineering constants. The corresponding results are shown in Figure 19.

With the introduced extraction tensors and damage parameters at hand, we are able to reproduce the experimental results obtained for the neat epoxy resin, as well es the predictions computed by Simon et al. (2017) quite accurately, see Figure 19(a). The decrease in the individual orthotropic engineering-constants are shown in Figure 19(b) to (d), where dashed lines correspond to our model and solid lines refer to the references Simon et al. (2017) and Bednarcyk et al. (2015). For all considered loading cases, our proposed modeling framework makes it simple to account for those engineering constants which remain unaffected during the loading. Figure 19(b) shows that the reduction of Young’s modulus $E_{33}$ is predicted correctly up to an applied strain of $10\%$ . The shear moduli are predicted quite well up to a strain level of $5\%$ , see Figure 19(c) and (d). Only the evolution of Young’s modulus $E_{22}$ shows some deviations beyond a strain level of $3\%$ , as our model is unable to capture the slope of reduction accurately enough for this case. To sum up, the proposed modeling framework permits reproducing the reduction of all affected engineering constants accurately for small deformations of up to $3\%$ . For the identified parameter set, we investigate the plain-weave composite that was also studied by Simon et al. (2017), see Figure 18(a). The mesostructure, with a required tow-volume fraction of approximately $75\%$ , was generated by a level-set approach developed by Sonon and co-workers (Sonon et al., 2012; Sonon and Massart, 2013; Wintiba et al., 2017) and discretized on a regular grid with $512\times 512\times 56$ voxels. Just as Simon et al. (2017), we investigate the longitudinal shear response of the plain-weave cell. For this purpose, we analyze the effective stress-strain curve, as well as the reduction of the shear modulus $G_{12}$ measured in terms of the relative error

$$△G_{12}=G_{12}^0-G_{12}$$

(4.8)

w. r. t. the initial, undamaged shear modulus $G_{12}^0$ . The predicted full-field distribution of this relative reduction is shown in Figure 18(b), and coincides well with the results presented by (Simon et al., 2017, Figure 9). Moreover, the effective stress-strain curves of the plain-weave composite subjected to longitudinal shear match well for the entire loading regime, see Figure 20.

Figure 19.

Comparison of predicted stress-strain curve and reductions of the orthotropic engineering constants based on introduced extraction tensors (see Table 4). Our model predictions are dashed in (b)–(d). (a) Neat epoxy behavior under normal loading. (b) Tow behavior under transverse normal loading. (c) Tow behavior under transverse shear loading. (d) Tow behavior under longitudinal shear loading.

Figure 20.

Comparison of effective stress-strain curves for plain weave under shear loading.

Summary and conclusions

In this article, a generalized standard material (GSM) model for anisotropic damage evolution based on the compliance tensor as the primary damage variable was developed. Based on the insight that the Hookean elastic energy density, considered as a function of the elastic strain and the compliance tensor, is a convex function of both arguments, a convex framework for quasi-static damage evolution was established, preventing damage localization intrinsically. Indeed, by choosing the energy (density) related to the progressive degradation of the material appropriately, the condensed incremental potential (Miehe, 2002) is strictly convex and of superlinear growth, which prevents localization for such a model. Of course, working with a softening-type energy for the damage-surface variables is also possible, and should be studied more closely in subsequent work.

The second section is organized to emphasize the modular fashion that the compliance-based damage model is built up. The model might be extended in subsequent work, for instance accounting for strain-rate sensitivity within the model. For an overview of the assumptions leading to specific specializations of the evolution of internal variables, we refer to the overview in Appendix A. The modeling framework is general enough to incorporate coupling to other inelastic models, such as plasticity or viscoplasticity (McDowell, 2008; Rousselier, 1979), entirely within the proposed framework. Also, due to its inherent stability, an extension to fatigue damage, as observed for certain fiber-reinforced polymers (Bartkowiak et al., 2019, 2020; Sauer and Richardson, 1980), appears promising (Magino et al., 2021).

The modular character of the model was exemplified by specific damage-extraction tensors motivated by Puck’s anisotropic failure criteria (Knops, 2008; Puck and Schürmann, 2002). With these ingredients at hand, we demonstrated the model’s capabilities of developing complex anisotropic stiffness states, not restricted a priori by a specific degree of (an)isotropy of the stiffness tensor, emphasizing that the model is capable of handling any initial stiffness. We also showed the model’s capabilities on meso and volume-element scale, based upon a straightforward numerical treatment. With these achievements at hand, accounting for additional failure criteria (Bouhala et al., 2013; Fritzsche et al., 2008; Kaddour et al., 2004) or coupling the model to phase-field fracture models (Gerasimov and De Lorenzis, 2019; Miehe et al., 2010; Schneider et al., 2016a) appears possible.

Returning to our original motivation, i.e., modeling anisotropic damage of SMC composite materials, requires incorporating the presented modeling framework into a three-scale homogenization scheme (Anagnostou et al., 2018). The underlying fiber bundle mesostructure (Dumont et al., 2007; Meyer et al., 2020; Schöttl et al., 2019) has to be accounted for, and the model parameters have to be fitted to experimental data. For the latter purpose, a convenient experimental program is necessary (Schemmann et al., 2018c).

From a mathematical perspective, a thorough mathematical analysis of our model is desirable, whereas an extension to tension-compression asymmetry appears imperative in order to model load reversals.