Mechanically consistent continuum damage model for anisotropic composites including damage deactivation

Abstract

Due to crack bridging effects, ceramic matrix composites (CMCs) have outstanding properties that combine a quasi-ductile material behaviour with the high-temperature properties of ceramics. Combined with their high specific strength, this makes them perfectly suitable for high temperature safety relevant components. In view of the design process of CMC components elaborated continuum damage models are required that most importantly consider their anisotropy and damage deactivation effects in a mechanically and mathematically consistent manner. With respect to their damage effect, most of the existing anisotropic models fail with regard to the damage growth criterion leading to the unphysical effect of an increasing stiffness due to damage. Motivated by the modelling process of initially anisotropic composite materials like CMCs, this paper presents the systematic formulation and validation of a mechanically consistent damage effect model together with crack closure effects.

Introduction

The present paper investigates the formulation and validation of a continuum damage model for the stiffness degradation of anisotropic brittle materials due to microcracks. Even though this model is applicable to anisotropic brittle composite materials in general, it is motivated by ceramic matrix composites (CMC). CMCs gain their interesting properties by the fibre crack bridging (Beyerle et al., 1992) that prevent a component failure due to microcracks and lead to a quasi-ductile behaviour with increased toughness, making this material a perfect candidate for safety relevant high temperature components in aviation turbines (Zawada et al., 2017).

The concept of continuum damage modelling was initially introduced and motivated by creep damage effects of metals by Kachanov (1958). The main idea is to include the effect of microcracks, or other evolving microdefects, into a continuum model by a damage variable $\mathcal{D}$ without explicitly modelling the microdefects. This clearly differentiate the field of continuum damage mechanics from the field of fracture mechanics which aims at explicitly modelling the effect and the evolution of cracks. Mostly the choice of the appropriate method is a matter of the related scales. If the effect of a single, sufficient large crack, is of interest, fracture mechanics is the method of choice. If the effect of multiply, sufficiently small cracks is of interest, a fracture mechanics approach becomes inaccurate or expensive, and a continuum mechanics approach is chosen instead. In any case, there is the possibility that a continuum damage mechanics approach is devolved into a fracture mechanics approach by the localization of microcracks into a single macroscopic crack (Geers et al., 1998).

With a focus on brittle homogeneous materials, the effect of microcracks on the effective material properties has been intensively investigated in the literature both experimentally and theoretically. From a theoretical point of view, effective properties are derived by homogenization of interacting and non-interacting cracks in (Kachanov, 1994). Effects of crack closure were for instance investigated in (Horii and Nemat-Nasser, 1983) using homogenization. The effect of friction of closed crack faces were analysed by homogenization in (Telega, 1990). Based on these results, different authors have proposed continuum damage models that are built upon homogenization approaches. Andrieux et al. (1986), Pensée et al. (2002), and Pensee and Kondo (2003) present damage mechanics models for initially isotropic materials like concrete considering crack closure effects. In Zhu et al. (2008) and Zhu et al. (2009) homogenization based models are presented including friction and the interaction of microcracks respectively.

Whereas all previously investigations mentioned focus on initially isotropic media and the anisotropic effect is only included by the damage, (Mauge and Kachanov, 1994) investigates in the damage effect of initially anisotropic media. Based on analytical homogenization results assuming a dilute crack distribution, a damage effect model for initially anisotropic media including crack closure effects was presented in Goidescu et al., 2013. Based thereon the model was supplemented with a damage evolution and implemented into a finite element framework in Goidescu et al., 2015. However these models are restricted to the 2-dimensional case, probably because no closed form solution exists for the general three dimensional anisotropic case, as pointed out by Mauge and Kachanov (1994).

Not relying on homogenization results, a phenomenological damage model is an alternative. This approach is especially relevant for composite materials where the microstructure of the composite materials is expected to strongly interact with evolving microcracks (Baste and Aristégui, 1998; Bernachy-Barbe et al., 2015a, 2015b) and a closed analytical solution is not available. Using a phenomenological model, the model complexity and accuracy can be increased by the rank of the damage variable $\mathcal{D}$ (Lubarda and Krajcinovic, 1993). Scalar damage variables d are typically used to describe damage with an a priori fixed direction (Besson et al., 2010; Maire and Chaboche, 1997) or as an approximation of higher rank tensors as in (Baranger, 2013; Marcin et al., 2011). Second rank tensors D are typically used to describe damage with a variable, initially unknown, direction (Besson et al., 2010; Maire and Chaboche, 1997). Fourth rank damage tensors $\mathbb{D}$ offer some additional modelling capabilities for non-proportional loading (Baranger, 2018; Carol et al., 2001), however, are rarely used due to the high model complexity.

To model the damage effect, i.e. the way the damage affects the stiffness degradation, one possibility is to refer to the concept of strain (Lemaitre and Chaboche, 1978) or energy (Cordebois and Sidoroff, 1979) equivalence. By this, modeling the damage effect is reduced to choosing the 4th rank damage effect tensor $\mathbf{M}\left(\mathbf{D}\right)$ as a function of the, mostly second rank, damage tensor D. However, as pointed out in (Lu and Chow, 1990), the strain equivalence leads to a non-symmetric degraded stiffness tensor such that this approach is questionable even if only the symmetric part is used as in (Cordebois and Sidoroff, 1982). Moreover it is unclear how the damage effect tensor $\mathbf{M}\left(\mathbf{D}\right)$ should be chosen to result in a mechanically consistent and accurate model. This aspect is also seen by the diversity of damage models that can build upon the very general hypothesis of strain or energy equivalence such as (Cordebois and Sidoroff, 1982; Lu and Chow, 1990; Murakami and Kamiya, 1997), to name only a few. Finally, as pointed out Wulfinghoff et al., 2017, and will further be discussed in this paper, most of the aforementioned models are mechanically inconsistent as they do not satisfy the micromechanically motivated damage growth criterion. This implies that an increase in damage can actually lead to an increase in stiffness, which contradicts the notion of damage modelling.

A second key challenge in the phenomenological modelling of damage is the description of damage deactivation that takes place as microcracks are closing under a compressive state. Beside the question of a mechanically feasible closure criterion (Cormery and Welemane, 2010) the key problem is rather mathematically, i.e. the formulation of a temporarily damage deactivation without violating the necessary smoothness of the strain energy density function (Chaboche et al., 1995).

In this paper a model for the stiffness degradation for anisotropic (brittle) materials is derived and validated that is consistent with the damage growth criterion. In a second step the model is furthermore supplemented with a damage deactivation satisfying all of the necessary smoothness requirements. The derivation is done in a systematically and rigorous way by making use of the representation theory of invariant tensor functions. This is a common way (Betten, 1992, 1998; Zheng et al., 1992; Zheng and Betten, 1995) to formulate material models relying on the very well established theory of invariants (Zheng, 1994). In this way the symmetry of the undamaged material is a priori included in the model formulation. The approach followed here shares some similarities with the derivation of the model in (Fassin et al., 2019), with the important difference that the latter is limited to isotropic undamaged materials which precludes its application to composite materials.

In the paper standard symbolic tensor notation is used and only the dyadic products ${\otimes }_{-}$ and ${\otimes }^{-\hspace{0.17em}}$ are further specified in the following to prevent confusion with other definitions sometimes used in the literature

$$\begin{array}{c}{(\mathbf{A}{\otimes }_{-}\mathbf{B})}_{\mathit{ijkl}}=A_{ik}{B}_{lj},{(\mathbf{A}{\otimes }^{-\hspace{0.17em}}\mathbf{B})}_{\mathit{ijkl}}=A_{il}{B}_{kj}.\end{array}$$

Model formulation for active damage

General framework

As a first assumption, the damage model defined here is restricted to the use of a second rank damage tensor D. This simplification is motivated by the fact that crack directions, and hence the brittle damage in composite materials, especially in CMCs, is strongly influenced by the loading direction (Bernachy-Barbe et al., 2015) and can not be predicted a priori. A second rank tensor is the minimum damage descriptor to describe such a load induced anisotropy (Besson et al., 2010) where the damage direction is given by the eigenvectors ${\mathbf{t}}_i$ of D and its principal values d_i. On the other side, using a second rank tensor, the model formulation and implementation is still practicable, especially since the corresponding representation theory is developed to a large extend (Zheng, 1994). An extension to a fourth rank tensor already becomes cumbersome especially since the corresponding representation theory for fourth rank tensors is not jet fully developed (Betten and Helisch, 1992; Zheng, 1994) and becomes rather unhandy.

The general modelling framework used as a starting point for this paper, is briefly summarized in the following. Let

$$\phi \left(\mathit{\epsilon },\mathbf{D}\right)={\phi }^0\left(\mathit{\epsilon }\right)+{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)=\frac{1}{2}\mathit{\epsilon }:{\mathbb{C}}^0:\mathit{\epsilon }+{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)$$

(1)

be the strain energy density function with ${\phi }^0$ and ${\phi }^D$ being the undamaged elastic and the damaged elastic part respectively. For brittle materials, the former is restricted to the usual anisotropic linear elasticity described by the 4^th rank stiffness tensor ${\mathbb{C}}^0$ . Introducing the damage power-conjugate forces

$$\mathbf{Y}:=-{\partial }_{\mathbf{D}}\phi =-{\partial }_{\mathbf{D}}{\phi }^D,$$

(2)

evaluation of the second law of thermodynamics leads to the stress-strain relation (assuming ${\phi }^D$ is at most quadratic in D)

$$\mathit{\sigma }\left(\mathit{\epsilon },\mathbf{D}\right)={\mathbb{C}}^0:\mathit{\epsilon }+{\partial }_{\mathit{\epsilon }}{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)$$

(3)

$$=\left({\mathbb{C}}^0+{\partial }_{\mathit{\epsilon }\mathit{\epsilon }}^2{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)\right):\mathit{\epsilon }$$

(4)

$$=:\left({\mathbb{C}}^0+\Delta \mathbb{C}\right):\mathit{\epsilon },$$

(5)

as well as the dissipation inequality

$$\mathbf{Y}:{\partial }_t\mathbf{D}\ge 0.$$

(6)

Following the maximum dissipation principle, the dissipation inequality is automatically satisfied if the damage evolution is derived from a convex dissipation potential $f\left(\mathit{\epsilon },\mathbf{D}\right)$ as

$${\partial }_t\mathbf{D}={\partial }_t\lambda {\partial }_{\mathbf{Y}}f\left(\mathit{\epsilon },\mathbf{D}\right),$$

(7)

whereas the damage multiplier λ and the dissipation potential must satisfy the Kuhn-Tucker relations

$${\partial }_t\lambda \ge 0,f\left(\mathit{\epsilon },\mathbf{D}\right)\le 0,{\partial }_t\lambda f\left(\mathit{\epsilon },\mathbf{D}\right)=0.$$

(8)

Symmetry considerations

In contrast to pure linear elasticity, the strain energy density (1) is defined by two tensor quantities, the strain tensor $\mathit{\epsilon }$ and the damage tensor D, both having principal values and principal directions. This gives rise to different symmetry definitions illustrated in Figure 1. Differentiation of these is important for the modelling process and will also help in the subsequent model interpretation.

Figure 1.

Illustration of the different symmetry definition, warp and weft yarns are indicated by gray lines, symmetry axes of the undammaged material are assumed to be aligned with the yarn directions which are identical with the basis vectors ${\mathbf{e}}_i$ of the global coordinate system. (a) Damaged material symmetry, (b) Undamaged material symmetry and (c) Symmetry of damage.

Symmetry of the damaged material

In the following, a material is defined as symmetric with respect to the damaged state if

$$\phi \left(\mathbf{Q}*\mathit{\epsilon },\mathbf{D}\right)=\phi \left(\mathit{\epsilon },\mathbf{D}\right),\forall \mathbf{Q}\in \chi ,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}$$

(9)

where $\chi \subseteq \mathrm{Orth}$ is the corresponding symmetry group, the star operator $(*)$ is used as an abbreviation for the rotation and reflection transformation, i.e. $\mathbf{Q}*\mathbf{A}={\mathbf{Q}}^T\mathbf{A}\mathbf{Q}$ and $\mathbf{Q}*\mathbb{B}={\left(\mathbf{Q}{\otimes }_{-}\mathbf{Q}\right)}^T:\mathbb{B}:\left(\mathbf{Q}{\otimes }_{-}\mathbf{Q}\right)$ for a second rank tensor A, and a fourth rank tensor $\mathbb{B}$ respectively. Again assuming the strain energy density $\phi$ is only quadratic in $\mathit{\epsilon }$ , this implies the following for the stress $\mathit{\sigma }$ and the stiffness tensor $\mathbb{C}$

$$\begin{array}{l}\mathbf{Q}*\mathit{\sigma }\left(\mathit{\epsilon },\mathbf{D}\right)=\mathit{\sigma }\left(\mathbf{Q}*\mathit{\epsilon },\mathbf{D}\right)=\mathbb{C}\left(\mathbf{D}\right):\mathbf{Q}*\mathit{\epsilon },\\ \mathbf{Q}*{\mathbb{C}}^0+\mathbf{Q}*\Delta \mathbb{C}\left(\mathbf{D}\right)={\mathbb{C}}^0+\Delta \mathbb{C}\left(\mathbf{D}\right),\\ \forall \mathbf{Q}\in \chi ,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}.\end{array}$$

(10)

The mechanical interpretation of this condition is as follows: If, in an experiment, the damaged material is transformed (e.g. rotated) by Q then the result of this experiment (e.g. the stress reaction) is the same as if the material were rotated. Strictly speaking, testing a material for this symmetry condition would require testing a material in different material orientations for all (relevant) possible damage states. Whereas ultrasonic time-of-flight measurements (Baste, 2001; Baste and Aristégui, 1998; Grippon et al., 2012) are perfectly suitable for the first part, testing all possible damage states is practically impossible.

Furthermore, it has been argued in the literature, that a certain symmetry property is very likely to be destroyed by damage (see for instance the discussions in Mauge and Kachanov, 1994; Goidescu et al., 2013). The same holds for the model derived later, in which Condition (9) is only satisfied in some special cases. However, it is helpful for the discussion of the model to look at the symmetry properties of the damage parts itself, i.e. the condition

$$\mathbf{Q}*\Delta {\mathbb{C}}_i\left(\mathbf{D}\right)=\Delta {\mathbb{C}}_i\left(\mathbf{D}\right),$$

(11)

where $\Delta {\mathbf{C}}_i$ are the individual damage effect tensors derived later (see equation (54)).

Symmetry of the undamaged material

If a material satisfies

$$\phi \left(\mathit{\epsilon },\mathbf{D}\right)=\phi \left(\mathbf{Q}*\mathit{\epsilon },\mathbf{Q}*\mathbf{D}\right),\forall \mathbf{Q}\in \chi ,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}$$

(12)

it is defined to be symmetric with respect to the undamaged state. The physical interpretation is as follows: if both, the strain and the damage are transformed by Q, the energy and the stress reaction remain unchanged. This is equivalent to transforming the undamaged material by ${\mathbf{Q}}^T$ (see Figure 1(b) for a graphical illustration). As the name already suggests, this property is solely defined by the undamaged material and hence can be analysed by for example ultrasonic measurements on the undamaged material or simply by geometrical arguments of the microstructure, without considering any damage state as in condition (9). In this context, it is also important to note that this condition does not, in any sense, restrict the symmetry properties of the damaged material.

Similar as before, this condition implies the following properties of the stress-strain relation and the stiffness tensor

$$\begin{array}{l}\mathbf{Q}*\mathit{\sigma }\left(\mathit{\epsilon },\mathbf{D}\right)=\mathbb{C}\left(\mathbf{Q}*\mathbf{D}\right):\mathbf{Q}*\mathit{\epsilon },\\ \mathbf{Q}*{\mathbb{C}}^0+\mathbf{Q}*\Delta \mathbb{C}\left(\mathbf{Q}*\mathbf{D}\right)={\mathbb{C}}^0+\Delta \mathbb{C}\left(\mathbf{D}\right),\\ \forall \mathbf{Q}\in \chi ,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym},\end{array}$$

(13)

Symmetry of damage

For the sake of completeness, the symmetry of damage is defined by the requirement

$$\phi \left(\mathit{\epsilon },\mathbf{D}\right)=\phi \left(\mathit{\epsilon },\mathbf{Q}*\mathbf{D}\right),\forall \mathbf{Q}\in \chi ,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}$$

(14)

From a mechanical point of view, this requires the behaviour to remain the same if the damage state is transformed by Q. Experimentally this condition is hard to verify since it would require samples with a known damage state having the same magnitude but different directions. Such a sample preparation is merely impossible to realize, especially for materials that are anisotropic in the undamaged state since damage evolution is largely influenced by the microstructure. However, from the mechanical interpretation it immediately follows that it should be satisfied for (a) a reflection transformation on a plane perpendicular to the eigenvectors ${\mathbf{t}}_i$ of D and (b) all orthogonal transformations Q if the damage is isotropic, i.e. $\mathbf{D}\sim \mathbf{1}$ .

Model formulation based on symmetry considerations

Based on these preliminaries the damage model is now defined by making use of the representation theory for tensor functions. This enables us, considering the symmetry of the undamaged material, to formulate the strain energy density in terms of elements of a complete function basis. According to (Zheng, 1994) a function basis is the set of isotropic invariants $\left\{I_1,\dots ,I_n\right\}$ such that every function $\phi \left(\mathit{\epsilon },\mathbf{D}\right)$ is expressible in terms of $I_1,\dots ,I_n$ only

$$\phi \left(\mathit{\epsilon },\mathbf{D}\right)=\tilde{\phi }\left(I_1,\dots ,I_n\right)$$

(15)

The function basis is called irreducible if none of the I_i is expressible by a single-valued function of the remainder. In this context, isotropy always refers to a certain symmetry group χ as for example transversal-isotropy or orthotropy.

The focus of this paper lays on the damage modelling of brittle composite materials where the individual ply is made from unidirectional or bi-directional (woven) reinforced materials. According to this, the assumption of an orthotropic undamaged (12) material is applied here. The symmetry group characterizing this type of symmetry is the orthotropic group

$$\chi =\left\{\mathbf{R}\left(\pi {\mathbf{e}}_3\right),{\mathbf{R}}_1,{\mathbf{R}}_2,-\mathbf{1}\right\},$$

(16)

with $\mathbf{R}\left(\pi {\mathbf{e}}_3\right)$ being a $180°$ rotation around the out-of-plane axis ${\mathbf{e}}_3$ and ${\mathbf{R}}_i$ being a reflection along a plane perpendicular to i.

A complete and irreducible function basis for an orthotropic $\phi \left(\mathit{\epsilon },\mathbf{D}\right)$ can be derived as a special case from the general case of a strain energy density $\mathit{\phi }$ dependent on n tensor arguments (Zheng, 1993) and is given by the set

$$\begin{array}{c}\mathcal{F}=\{\text{tr}\left(\mathbf{D}\right),\text{tr}\left({\mathbf{D}}^2\right),\text{tr}\left({\mathbf{D}}^3\right),\text{tr}\left(\mathbf{MD}\right),\text{tr}\left(\mathbf{M}{\mathbf{D}}^2\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\text{tr}\left({\mathbf{M}}^2{\mathbf{D}}^2\right),\text{tr}\left(\mathbf{M}{\mathbf{D}}^2\right),\text{tr}\left(\mathit{\epsilon }\right),\text{tr}\left({\mathit{\epsilon }}^2\right),\text{tr}\left({\mathit{\epsilon }}^3\right),\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right),\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right),\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right),\text{tr}\left(\mathit{\epsilon }\mathbf{D}\right),\text{tr}\left(\mathit{\epsilon }{\mathbf{D}}^2\right)\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right),\text{tr}\left(\mathbf{M}\mathit{\epsilon }\mathbf{D}\right)\}.\end{array}$$

(17)

The symmetry axes ${\mathbf{e}}_i$ enter the function basis by the structural tensor

$$\mathbf{M}={\mathbf{e}}_2\otimes {\mathbf{e}}_2-{\mathbf{e}}_3\otimes {\mathbf{e}}_3$$

(18)

which uniquely defines the symmetry group of interest as

$$\chi =\left\{\mathbf{Q}\in \mathrm{Orth}:\mathbf{Q}*\mathbf{M}=\mathbf{M}\right\},$$

(19)

Based on the function basis (17), a general prototype of the damage effect function can be specified. For this it is assumed that the damage effect ${\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)$ in (1) is polynomial. For brittle materials like CMCs, numerous experimental results are available showing that the stress-strain behaviour remains fully linear for a fixed damage D. Therefore, the polynomial degree with respect to strain is restricted to be exactly two. A linear part in $\mathit{\epsilon }$ would involve an irreversible plastic strain which is a secondary effect in CMCs and not the focus of the current investigation.

To further reduce the complexity, the polynomial order, in the damage variable D, is restricted to be of the order one. This simplifies the damage effect to

$${\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)=\underset{(1)}{\underbrace{{\xi }^0{I}^g}}+\underset{(9)}{\underbrace{{\xi }_{ij}^1{I}_i^a{I}_j^f}}+\underset{(18)}{\underbrace{{\xi }_{\mathit{ijk}}^2{I}_i^a{I}_j^c{I}_k^c}}+\underset{(6)}{\underbrace{{\xi }_{ij}^3{I}_i^c{I}_j^d}}$$

(20)

whereas the Invariants $I^a,\dots ,I^g$ are defined in Table 1. However, even this simple linear damage form includes 34 material parameters ${\xi }^0,{\xi }_{ij}^1,{\xi }_{\mathit{ijk}}^2$ and ${\xi }_{ij}^3$ which is still not manageable for most practical applications and therefore the stated damage effect function need to be further reduced in the following.

Table 1.

Orthotopic function basis used for the formulation of the damage effect function.

Order of $\mathit{\epsilon }$	Order of D	Invariants	Notation
0	1	trD, trMD, trM2D	Ia
0	2	trD2, trMD2, trM2D2	Ib
0	3	trD3	–
1	0	trε, trMε, trM2ε	Ic
1	1	trεD, trMεD	Id
1	2	trεD2	Ie
2	0	trε2, trMε2, trM2ε2	If
2	1	trε2D	Ig
3	0	trε3	–

The gray elements refer to the redundant (reducible) invariants for the special case of a planar damage investigated later.

Reduction for in-plane damage

In the first step the function basis (17) is reduced by further specifying the damage state. This is motivated by separating the damage of composite materials according to the scale-separation. The damage of composite materials naturally consists of two damage mechanisms, the inter-laminar and the intra-laminar damage. The former is naturally occurring in between the plies where the damage has a preferential out-of-plane direction. Furthermore, due to the macroscopic, structural character, the separation of scales is not applicable here and inter-laminar damage is commonly explicitly modelled by cohesive zone models (Kumar and Welsh, 2012). The intra-laminar damage in turn has a microscopic character and is caused by the formation of micro-cracks. The crack normal n normally has a preferential in-plane direction. As long as these microcracks do not coalesce to a macro-crack, this damage mechanism is modelled by a continuum damage model with a damage tensor restricted to the in-plane components:

$$\mathbf{D}=d_1{\mathbf{t}}_1\otimes {\mathbf{t}}_1+d_2{\mathbf{t}}_2\otimes {\mathbf{t}}_2,$$

(21)

whereas ${\mathbf{t}}_1$ and ${\mathbf{t}}_2$ are the corresponding eigenvectors laying in the composite (or ply) plane. By this requirement the number of independent components in the damage tensor D has been reduced. This in turn can be used to reduce the number of irreducible invariants in (17).

For the reduction the isotropization theorem is applied. The isotropization theorem (Itzkov, 2015) is the fundamental result that helps in deriving a function basis for anisotropic functions. It basically states that an anisotropic function, depending on a certain number of tensor valued arguments, say $\phi \left(\mathit{\epsilon },\mathbf{D}\right)$ , is expressible as an isotropic function $\tilde{\varphi }$ adding the structural tensor M as an additional argument

$$\phi \left(\mathit{\epsilon },\mathbf{D}\right)=\tilde{\phi }\left(\mathit{\epsilon },\mathbf{D},\mathbf{M}\right).$$

(22)

The function basis for an anisotropic function can thus be derived in the same manner as the function basis of an isotropic function adding the fixed structural tensor M as an additional entity. Therefore the following proof follows the steps in deriving a complete and irreducible function basis for isotropic functions with an arbitrary number of symmetric, and skew-symmetric tensors and vector arguments as presented in (Smith, 1971). However the proof here is concise, since according to (22) only three symmetric tensor arguments need to be considered, one of them being constant. Furthermore, since only the scalar valued strain energy density is of interest and the tensor valued stress $\mathit{\sigma }\left(\mathit{\epsilon },\mathbf{D}\right)$ relation as well as the conjugate forces Y are derived by differentiation from (2) and (5), the derivation is restricted to the scalar invariants and the corresponding tensor-valued forminvariants are not derived.

The first step is to derive a complete set of orthotropic invariants and show its irreducibility in the second step. For the first step a set $I_1,\dots ,I_n$ is assembled such that every function (22) is expressible as

$$\varphi \left(\mathit{\epsilon },\mathbf{D},\mathbf{M}\right)=\overline{\varphi }\left(I_1,\dots ,I_n\right),$$

(23)

with the entities given by

$$\begin{array}{l}\mathbf{D}=d_{11}{\mathbf{e}}_1\otimes {\mathbf{e}}_1+d_{12}\left({\mathbf{e}}_1\otimes {\mathbf{e}}_2+{\mathbf{e}}_2\otimes {\mathbf{e}}_1\right)+d_{22}{\mathbf{e}}_2\otimes {\mathbf{e}}_2,\\ \mathit{\epsilon }={\epsilon }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j,\text{with}{\epsilon }_{ij}={\epsilon }_{ji},\\ \mathbf{M}={\mathbf{e}}_2\otimes {\mathbf{e}}_2-{\mathbf{e}}_3\otimes {\mathbf{e}}_3.\end{array}$$

(24)

Thus we simply have to find a set of invariants such that the components of $\mathbf{D},\mathit{\epsilon },\mathbf{M}$ can be uniquely determined from this set, which then proofs the equality (23). For this the principal axis coordinate system ${\mathbf{e}}_i$ is chosen for all entities since here components of M are known a priori. Indeed it would also be possible to determine the components of D in the principal system ${\mathbf{t}}_i$ however, then the relation between ${\mathbf{t}}_i$ and ${\mathbf{e}}_i$ is additionally unknown.

Since M has three distinct principal values, the components d₁₁, d₂₂ and $d_{12}^2$ can be derived from

$$\text{tr}\left(\mathbf{D}\right)=d_{11}+d_{22},$$

(25)

$$\text{tr}\left(\mathbf{MD}\right)=d_{22},$$

(26)

$$\text{tr}\left({\mathbf{D}}^2\right)=d_{11}^2+2d_{12}^2+d_{22}^2.$$

(27)

Since $M_{33}^2\ne M_{33}$ , the diagonal strain components ${\epsilon }_{ii}$ are obtained from

$$\text{tr}\left(\mathit{\epsilon }\right)={\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33},$$

(28)

$$\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)=M_{22}{\epsilon }_{22}+M_{33}{\epsilon }_{33},$$

(29)

$$\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)=M_{22}^2{\epsilon }_{22}+M_{33}^2{\epsilon }_{33}.$$

(30)

The squared off diagonal components ${\epsilon }_{ij}^2$ can be determined from the the following system of equations

$$\text{tr}\left({\mathit{\epsilon }}^2\right)={\epsilon }_{11}^2+2{\epsilon }_{12}^2+2{\epsilon }_{23}^2+2{\epsilon }_{13}^2+{\epsilon }_{22}^2+{\epsilon }_{33}^2,$$

(31)

$$\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right)=M_{22}\left({\epsilon }_{12}^2+{\epsilon }_{22}^2+{\epsilon }_{23}^2\right)+M_{33}\left({\epsilon }_{13}^2+{\epsilon }_{23}^2+{\epsilon }_{33}^2\right),$$

(32)

$$\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)=M_{22}^2\left({\epsilon }_{12}^2+{\epsilon }_{22}^2+{\epsilon }_{23}^2\right)+M_{33}^2\left({\epsilon }_{13}^2+{\epsilon }_{23}^2+{\epsilon }_{33}^2\right).$$

(33)

Next, the product ${\epsilon }_{12}{\epsilon }_{13}{\epsilon }_{23}$ is obtained from

$$\begin{array}{l}\text{tr}\left({\mathit{\epsilon }}^3\right)={\epsilon }_{11}^3+{\epsilon }_{22}^3+{\epsilon }_{33}^3+3{\epsilon }_{11}\left({\epsilon }_{12}^2+{\epsilon }_{13}^2\right)+3{\epsilon }_{22}\left({\epsilon }_{12}^2+{\epsilon }_{23}^2\right)\\ +3{\epsilon }_{33}\left({\epsilon }_{13}^2+{\epsilon }_{23}^2\right)+6{\epsilon }_{12}{\epsilon }_{13}{\epsilon }_{23},\end{array}$$

(34)

which can be used to determine the sign of the off diagonals ${\epsilon }_{ij}$ by investigating the following cases

1. Let all off-diagonals be non-zero ( ${\epsilon }_{ij}\ne 0,\forall i\ne j$ ), then, without affecting the validity of the previous results, the direction of ${\mathbf{e}}_1$ and ${\mathbf{e}}_3$ is chosen such that ${\epsilon }_{12}$ and ${\epsilon }_{13}$ are positive. By this, (34) gives the sign of ${\epsilon }_{23}$ . The off diagonal d₁₂ is defined by

   $$\text{tr}\left(\epsilon \mathbf{D}\right)={\epsilon }_{11}{d}_{11}+2{\epsilon }_{12}{d}_{12}+{\epsilon }_{22}{d}_{22}$$

   (35)

In that case $\text{tr}\left({\mathbf{D}}^2\right)$ from (27) is not required.

2. Exactly one off-diagonals from $\mathit{\epsilon }$ is zero, suppose ${\epsilon }_{12}=0$ . The direction of ${\mathbf{e}}_3$ and ${\mathbf{e}}_2$ is chosen such that ${\epsilon }_{12}$ and ${\epsilon }_{13}$ are positive. In contrast to the first case

$$\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right)=\left({\epsilon }_{11}^2+{\epsilon }_{12}^2+{\epsilon }_{13}^2\right)d_{11}+\left({\epsilon }_{12}^2+{\epsilon }_{22}^2+{\epsilon }_{23}^2\right)d_{22}+2\left({\epsilon }_{11}{\epsilon }_{12}+{\epsilon }_{12}{\epsilon }_{22}+{\epsilon }_{13}{\epsilon }_{23}\right)d_{12}$$

(36)

is now required to determine the sign of d₁₂ since the coefficient ${\epsilon }_{12}$ is zero in (35). If ${\epsilon }_{23}$ or ${\epsilon }_{13}$ are the only off-diagonals being zero, the procedure is analogous, although in that situation (35) would be sufficient to determine d₁₂.

3. The strain $\mathit{\epsilon }$ does not satisfy the previous two cases, but D and $\mathit{\epsilon }$ have different non-zero off diagonals. Since $d_{13}=d_{23}=0$ by definition, this requires $d_{12}\ne 0$ and

$${\epsilon }_{12}=0,{\epsilon }_{23}=0,{\epsilon }_{13}\ne 0$$

(37)

$$\text{or}{\epsilon }_{12}=0,{\epsilon }_{23}\ne 0,{\epsilon }_{13}=0$$

(38)

In the first case the direction of ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ is defined such that $d_{12}>0$ and ${\epsilon }_{13}>0$ , the value of d₁₂ is then given by (27). Case (38) can be proceeded in the same way.

4. It remains to investigate the two cases

$$d_{12}\ne 0,{\epsilon }_{23}={\epsilon }_{13}={\epsilon }_{12}=0\text{\hspace{1em}or\hspace{1em}}{\epsilon }_{12}\ne 0$$

(39)

$${\epsilon }_{12}\ne 0\text{\hspace{1em}and\hspace{1em}}{d}_{12}=0\text{\hspace{1em}or\hspace{1em}}{d}_{12}\ne 0$$

(40)

In the former case, the direction of ${\mathbf{e}}_2$ is chosen such that $d_{12}>0$ , whereas its value is given by (27). ${\epsilon }_{12}$ is then calculated from (35). In the latter case, the direction of ${\mathbf{e}}_2$ is chosen such that ${\epsilon }_{12}>0$ and d₁₂ is calculated from (35).

In total, restricting the damage state to in-plane damage, leads to the complete function basis

$$\begin{array}{l}\text{tr}\left(\mathbf{D}\right),\text{tr}\left({\mathbf{D}}^2\right),\text{tr}\left(\mathbf{MD}\right),\text{tr}\left(\mathit{\epsilon }\right),\text{tr}\left({\mathit{\epsilon }}^2\right),\text{tr}\left({\mathit{\epsilon }}^3\right),\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right),\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right),\\ \text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right),\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right),\text{tr}\left(\mathit{\epsilon }\mathbf{D}\right),\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right)\end{array}$$

(41)

having only 12 elements instead of 18 in (17). This in turn implies a reduction from 34 to 22 material parameters in (20).

The second part of the proof, i.e. the irreducibility of (41) is given in the Appendix.

Reduction by the damage growth criterion

The damage growth criterion derived and further discussed in (Wulfinghoff et al., 2017) is motivated by a micromechanical investigation of damage. Assuming a separation of scales, i.e. excluding macro cracks and assuming a displacement field with sufficiently small fluctuations, the usual micro-macro relations following from the Hill-Principle (Hill, 1963) can be applied. In this context the macroscopic strain energy density function is derived from a minimizer of the attached micro-problem Ω _m

$$\phi \left(\mathit{\epsilon }\right)=\underset{u_m\in \Theta \left(\mathit{\epsilon }\right)}{\min }〈{\phi }_m\left({\mathit{\epsilon }}_m\right)〉,$$

(42)

with $\Theta \left(\mathit{\epsilon }\right)$ being the space of admissible micro-displacement fields, which, restricting ourselves to the case where cracks are the only cause of damage, is expressible as

$$\Theta \left(\mathit{\epsilon }\right)=\left\{{\mathbf{u}}_m:{\mathbf{u}}_m=\mathit{\epsilon }{\mathbf{x}}_m,\forall {\mathbf{x}}_m\in \partial {\Omega }_m,[[\mathbf{u}]]·\mathbf{n}\ge 0\text{\hspace{0.17em}\hspace{0.17em}on\hspace{0.17em}\hspace{0.17em}}\Gamma \right\},$$

(43)

with Γ being the crack surface. Next, consider a second (primed) state with additional, or grown, cracks $\Gamma \prime$ with

$$\Gamma \subset \Gamma \prime \text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}\text{hence}\Theta \left(\mathit{\epsilon }\right)\subset \Theta \prime \left(\mathit{\epsilon }\right).$$

(44)

Since the space of admissible deformation increases by the growth of microcracks, the following inequality holds for the macroscopic strain energy density $\phi \prime$

$$\phi \prime \left(\mathit{\epsilon }\right)=\underset{{\mathbf{u}}_m\in {\Theta }^{\prime }\left(\mathit{\epsilon }\right)}{\min }〈{\phi }_m\left({\mathit{\epsilon }}_m\right)〉\le \underset{{\mathbf{u}}_m\in \Theta \left(\mathit{\epsilon }\right)}{\min }〈{\phi }_m\left({\mathit{\epsilon }}_m\right)〉=\phi \left(\mathit{\epsilon }\right).$$

(45)

Importantly, this relation holds for arbitrary admissible macroscopic strain fields $\mathit{\epsilon }$ . If the damage state is alternatively parametrized by a positive semi-definite damage tensor D, the inequality (45) can be written as

$$\phi \left(\mathit{\epsilon },\mathbf{D}+\text{d}\mathbf{D}\right)\le \phi \left(\mathit{\epsilon },\mathbf{D}\right),\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}.$$

(46)

Writing this in the equivalent form

$$-{\partial }_{\mathbf{D}}\phi \left(\mathit{\epsilon },\mathbf{D}\right):\text{d}\mathbf{D}=\mathbf{Y}:\text{d}\mathbf{D}\ge 0,\forall \left(\mathbf{D},\mathit{\epsilon }\right)\in \mathrm{psym}\times \mathrm{sym}$$

(47)

shows a strong analogy to the dissipation inequality (6), whereas $\text{d}\mathbf{D}$ could easily be substituted by ${\partial }_t\mathbf{D}$ considering a parametrization by time. The crucial difference now being that the dissipation inequality (6) is only required for the actual deformation-damage history $\left\{\mathit{\epsilon },\mathbf{D}\right\}$ which also involves the dissipation potential. The damage growth criterion in turn is required for all possible deformation-damage states. In that sense, the damage growth criterion can be interpreted as a generalization of the second law.

The damage growth criterion is next used for a further reduction of the formulated damage effect function. Applying (47) to (83) or (20) leads

$$\mathbf{Y}:\text{d}\mathbf{D}=-{\Sigma }_i{\partial }_{\mathbf{D}}{\phi }_i^D\left(\mathit{\epsilon },\mathbf{D}\right):\text{d}\mathbf{D}\ge 0$$

(48)

where ${\partial }_{\mathbf{D}}{\phi }_i^D$ are the derivatives with respect to each term in (83) or (20). Since this criterion is hard to handle mathematically, the more severe set of sufficient criteria

$$-{\partial }_{\mathbf{D}}{\phi }_i^D\left(\mathit{\epsilon },\mathbf{D}\right):\text{d}\mathbf{D}\ge 0,\forall i=1,\dots ,N$$

(49)

are used instead. The damage D is always positive definite by definition, thus the damage growth criteria (49) are satisfied if all parts of the conjugate forces ${\mathbf{Y}}_i=-{\partial }_{\mathbf{D}}{\phi }_i^D\left(\mathit{\epsilon },\mathbf{D}\right)$ are positive definite. This statement is trivial if $\text{d}\mathbf{D}$ and Y are colinear, but also holds if their eigensystems do not coincide. The latter can be seen by using the component representation of Y and $\text{d}\mathbf{D}$ in the corresponding principal frames ${\mathbf{a}}_i$ and ${\mathbf{b}}_i$ respectively

$$\begin{array}{l}\mathbf{Y}:\text{d}\mathbf{D}=Y_i{\mathbf{a}}_i\otimes {\mathbf{a}}_i:\text{d}{D}_j{\mathbf{b}}_j\otimes {\mathbf{b}}_j=\left(Y_1{\left({\mathbf{a}}_1·{\mathbf{b}}_1\right)}^2+Y_2{\left({\mathbf{a}}_2·{\mathbf{b}}_1\right)}^2\right)\text{d}{D}_1\\ +\left(Y_1{\left({\mathbf{a}}_1·{\mathbf{b}}_2\right)}^2+Y_2{\left({\mathbf{a}}_2·{\mathbf{b}}_2\right)}^2\right)\text{d}{D}_2\ge 0\end{array}$$

(50)

with $\text{d}{D}_i$ being positive but otherwise arbitrary, it follows that $Y_i\ge 0$ ensures a positive damage growth rate.

Applying the criterion (49) to the individual terms in the linear damage effect (20) helps to further drastically reduce the complexity of the damage effect. The individual terms, together with its conjugate forces ${\mathbf{Y}}_i\ge 0$ are summarized in Table 2, herein only the conjugate forces printed in black are positive definite if the related material parameter ξ is negative.

Table 2.

Terms of the damage effect function ${\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)$ under consideration of a planar damage state according to (21).

ξi	${\phi }_i^D\left(\mathit{\epsilon },\mathbf{D}\right)$	$\mathbf{Y}=-{\partial }_{\mathbf{D}}\phi \left(\mathit{\epsilon },\mathbf{D}\right)$
${\xi }^0{I}^g$
${\xi }^0$	$\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right)$	$-\mathit{\epsilon }$
${\xi }_{ij}^1{I}_i^a{I}_j^f$
${\xi }_{11}^1$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left({\mathit{\epsilon }}^2\right)$	$-\text{tr}\left({\mathit{\epsilon }}^2\right)\mathbf{1}$
${\xi }_{12}^1$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right)$	$-\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right)\mathbf{1}$
${\xi }_{13}^1$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)$	−tr(M2ε2)1
${\xi }_{11}^1$	tr (MD) tr(ε2)	$-\text{tr}\left({\mathit{\epsilon }}^2\right)\mathbf{M}$
${\xi }_{12}^1$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right)$	$-\text{tr}\left(\mathbf{M}{\mathit{\epsilon }}^2\right)\mathbf{M}$
${\xi }_{13}^1$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)$	$-\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)\mathbf{M}$
${\xi }_{\mathit{ijk}}^3{I}_i^a{I}_j^c{I}_k^c$
${\xi }_{111}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}{\left(\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left(\mathit{\epsilon }\right)}^2\mathbf{1}$
${\xi }_{112}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\mathbf{1}$
${\xi }_{113}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\mathbf{1}$
${\xi }_{122}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2\mathbf{1}$
${\xi }_{123}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\mathbf{1}$
${\xi }_{123}^2$	$\text{tr}\left(\mathbf{D}\right)\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2\mathbf{1}$
${\xi }_{211}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}{\left(\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left(\mathit{\epsilon }\right)}^2\mathbf{M}$
${\xi }_{212}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\mathbf{M}$
${\xi }_{213}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\mathbf{M}$
${\xi }_{222}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2\mathbf{M}$
${\xi }_{223}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)$	$-\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\mathbf{M}$
${\xi }_{223}^2$	$\text{tr}\left(\mathbf{MD}\right)\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2$	$-\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2\mathbf{M}$
${\xi }_{ij}^7{I}_i^c{I}_j^d$
${\xi }_{11}^3$	$\text{tr}\left(\mathit{\epsilon }\right)\text{tr}\left(\mathit{\epsilon }\mathbf{D}\right)$	$-\text{tr}\left(\mathit{\epsilon }\right)\mathit{\epsilon }$
${\xi }_{21}^3$	$\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\text{tr}\left(\mathit{\epsilon }\mathbf{D}\right)$	$-\text{tr}\left(\mathbf{M}\mathit{\epsilon }\right)\mathit{\epsilon }$
${\xi }_{31}^3$	$\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\text{tr}\left(\mathit{\epsilon }\mathbf{D}\right)$	$-\text{tr}\left({\mathbf{M}}^2\mathit{\epsilon }\right)\mathit{\epsilon }$

The gray and black printed conjugate forces are either negative definite or indefinite and positive definite respectively.

Summarizing this result and renumbering the material parameters to ξ_i with $i=1,\dots ,6$ leads the compact linear damage effect function

$${\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)={\xi }_1\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right)+\text{tr}\left(\mathbf{D}\right)[{\xi }_2\text{tr}\left({\mathit{\epsilon }}^2\right)+{\xi }_3\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)+{\xi }_4\text{tr}{\left(\mathit{\epsilon }\right)}^2+{\xi }_5\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2+{\xi }_6\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2]$$

(51)

which is consistent with the damage growth criterion if

$${\xi }_i\le 0$$

(52)

By this, the degraded stiffness tensor can be derived from (5) as

$$\mathbb{C}\left(\mathbf{D}\right)={\partial }_{\mathit{\epsilon }}^2\phi \left(\mathit{\epsilon },\mathbf{D}\right)={\mathbb{C}}^0+{\Sigma }_i{\xi }_i\Delta {\mathbb{C}}_i\left(\mathbf{D}\right)$$

(53)

with

$$\begin{array}{l}\Delta {\mathbb{C}}_1=\frac{1}{2}\left[\mathbf{1}{\otimes }_{-}\mathbf{D}+\mathbf{1}{\otimes }^{-\hspace{0.17em}}\mathbf{D}+\mathbf{D}{\otimes }_{-}\mathbf{1}+\mathbf{D}{\otimes }^{-\hspace{0.17em}}\mathbf{1}\right],\Delta {\mathbb{C}}_2=\text{tr}\left(\mathbf{D}\right)\left[\mathbf{1}{\otimes }_{-}\mathbf{1}+\mathbf{1}{\otimes }^{-\hspace{0.17em}}\mathbf{1}\right],\Delta {\mathbb{C}}_3=\text{tr}\left(\mathbf{D}\right)\frac{1}{2}\left[\mathbf{1}{\otimes }_{-}{\mathbf{M}}^2+\mathbf{1}{\otimes }^{-\hspace{0.17em}}{\mathbf{M}}^2+{\mathbf{M}}^2{\otimes }_{-}\mathbf{1}+{\mathbf{M}}^2{\otimes }^{-\hspace{0.17em}}\mathbf{1}\right],\Delta {\mathbb{C}}_4=\text{tr}\left(\mathbf{D}\right)2\left[\mathbf{1}\otimes \mathbf{1}\right],\Delta {\mathbb{C}}_5=\text{tr}\left(\mathbf{D}\right)2\left[\mathbf{M}\otimes \mathbf{M}\right],\Delta {\mathbb{C}}_6=\text{tr}\left(\mathbf{D}\right)2\left[{\mathbf{M}}^2\otimes {\mathbf{M}}^2\right].\end{array}$$

(54)

Model and parameter interpretation

For a deeper insight into the derived damage model and an interpretation of the material parameters it is helpful to look at the elongation modulus $L\left(\mathbf{m}\right)$ and the bulk modulus $\kappa \left(\mathbf{m}\right)$ that fully characterize the anisotropic material behaviour (Böhlke and Brüggemann, 2001). From a pure elongation test in the direction m, i.e. $\mathit{\epsilon }={\epsilon }_m\mathbf{m}\otimes \mathbf{m}$ , the elongation- and bulk modulus are defined by

$${\sigma }_m=L\left(\mathbf{m}\right){\epsilon }_m\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}\text{tr}\left(\mathit{\sigma }\right)=\kappa \left(\mathbf{m}\right){\epsilon }_m$$

(55)

It follows that

$$L\left(\mathbf{m}\right)=\mathbb{C}::\mathbf{m}\otimes \mathbf{m}\otimes \mathbf{m}\otimes \mathbf{m},$$

(56)

$$\kappa \left(\mathbf{m}\right)=\text{tr}\left(\mathbb{C}:\mathbf{m}\otimes \mathbf{m}\right).$$

(57)

For the following discussion, the individual contributions to the elongation- and bulk-modulus according to the stiffness terms in (54) are summarized in Table 3 for in-plane strain states with $\mathbf{m}=\left(\cos \omega ,\sin \omega ,0\right)$ .

Table 3.

In-plane damage in terms of the elongation- $\Delta L_i\left(\omega \right)$ and bulk modulus $\Delta K_i\left(\omega \right)$ .

i	$\Delta L_i$	$\Delta K_i$
1	$2\left({\cos }^2\left(\omega \right)d_{11}+\sin \left(2\omega \right)d_{12}+{\sin }^2\left(\omega \right)d_{22}\right)$	$\Delta L_1$
2	$2\text{tr}\left(\mathbf{D}\right)$	$2\text{tr}\left(\mathbf{D}\right)$
3	$2\text{tr}\left(\mathbf{D}\right)$	$2\text{tr}\left(\mathbf{D}\right)$
4	$2\text{tr}\left(\mathbf{D}\right)$	$6\text{tr}\left(\mathbf{D}\right)$
5	$2{\cos }^2\left(2\omega \right)\text{tr}\left(\mathbf{D}\right)$	0
6	$2\text{tr}\left(\mathbf{D}\right)$	$4\text{tr}\left(\mathbf{D}\right)$

Inspecting Table 3 it is readily seen that $\Delta {\mathbb{C}}_2,\Delta {\mathbb{C}}_3,\Delta {\mathbb{C}}_4$ and $\Delta {\mathbb{C}}_6$ are in-plane isotropic with respect to the damaged material (see Condition (9)). It follows that the in-plane symmetry of the undamaged material, defined by ${\mathbb{C}}_0$ , is not affected by these terms.

Next, Table 3 shows that $\Delta {\mathbb{C}}_2$ is redundant to $\Delta {\mathbb{C}}_4$ and $\Delta {\mathbb{C}}_6$ is redundant to $\Delta {\mathbb{C}}_4$ if either of $\Delta {\mathbb{C}}_2$ or $\Delta {\mathbb{C}}_3$ is included. This motivates the formulation of the following simplified version of the proposed damage model

$$\begin{array}{c}{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)={\xi }_1\text{tr}\left({\mathit{\epsilon }}^2\mathbf{D}\right)+\text{tr}\left(\mathbf{D}\right)[\hspace{0.17em}{\xi }_3\text{tr}\left({\mathbf{M}}^2{\mathit{\epsilon }}^2\right)+{\xi }_5\text{tr}{\left(\mathbf{M}\mathit{\epsilon }\right)}^2+{\xi }_6\text{tr}{\left({\mathbf{M}}^2\mathit{\epsilon }\right)}^2\hspace{0.17em}],\end{array}$$

(58)

which can be used if the focus lays on in-plane loading scenarios.

The terms $\Delta {\mathbb{C}}_1$ and $\Delta {\mathbb{C}}_5$ both possess a damaged orthotropic symmetry, whereas the latter is even tetragonal for in-plane loadings. To illustrate this and gain further insights, their corresponding elongation- and bulk modulus damage are plotted in Figure 2 for different damage states d_i and principal directions, denoted herein as α.

Figure 2.

Illustration of the symmetry properties of $\Delta {\mathbb{C}}_1$ and $\Delta {\mathbb{C}}_5$ investigated by their in-plane elongation damage moduli $\Delta L_1\left(\omega \right)$ (a) and $\Delta L_5\left(\omega \right)$ (b). For $\Delta L_1$ eigenvectors of D are exemplarily illustrated for the case of an uniaxial damage with a principal direction of $\alpha =22.5$ . For $\Delta L_5$ no further differentiation with respect to the damage state is necessary.

For $\Delta {\mathbb{C}}_1$ Figure 2(a) clearly reveals the dependency on the principal direction of D: The damage effect is orthotropic with respect to the damaged material whereas the symmetry axes are given by the eigenvectors ${\mathbf{t}}_i$ . Thus, if the undamaged material ${\mathbb{C}}^0$ is orthotropic, the symmetry of the undamaged material is fully disturbed by the damage effect $\Delta {\mathbb{C}}_1$ only if ${\mathbf{t}}_i$ , the eigenvectors of D, are not aligned with ${\mathbf{e}}_i$ , the symmetry axes of the undamaged material. In the case where ${\mathbb{C}}^0$ has a higher symmetry (e.g. tetragonality, transversal-isotropy or isotropy) the effect of $\Delta {\mathbb{C}}_1$ is similar with the exception that the symmetry of the undamaged material is at least degraded to orthotropy by damage.

Furthermore Figure 2(a) demonstrates the two interesting special cases of an uniaxial and isotropic damage with $\mathbf{D}=d_1{\mathbf{t}}_1\otimes {\mathbf{t}}_1$ and $\mathbf{D}=d\mathbf{1}$ respectively. In the case of an uniaxial damage, the transverse damage effect, i.e. damage in the direction ${\mathbf{t}}_2$ , fully vanishes. In case of an isotropic damage the damage effect also becomes fully isotropic, which is expected since in this case every vector is an eigenvector.

For $\Delta {\mathbb{C}}_5$ , being tetragonal for in-plane loadings, the picture is completely different compared to $\Delta {\mathbb{C}}_1$ . As is seen from Figure 2(b) the symmetry axes of the damage effect are now fixed and identical to the symmetry axes of the undamaged material. Therefore the damage effect is independent of the direction of D, since the damage only enters the formulation by the isotropic invariant $\text{tr}\left(\mathbf{D}\right)$ . But in contrast to $\Delta {\mathbb{C}}_2,\Delta {\mathbb{C}}_3,\Delta {\mathbb{C}}_4,\Delta {\mathbb{C}}_6$ the basistensor $\mathbf{M}\otimes \mathbf{M}$ (see equation (54)) is not isotropic. Thus $\Delta {\mathbb{C}}_5$ disturbs the symmetry of the undamaged material to a damaged tetragonal symmetry if ${\mathbf{C}}_0$ inherits a symmetry higher than tetragonality.

We next turn to the out-of-plane behavior: Excluding the case of an isotropic damage with $\mathbf{D}=d\mathbf{1},\Delta {\mathbb{C}}_2$ and $\Delta {\mathbb{C}}_4$ can be shown to be the only damage effects that are fully isotropic with respect to the damaged material. In the same sense, the remaining terms $\Delta {\mathbb{C}}_3,\Delta {\mathbb{C}}_5$ and $\Delta {\mathbb{C}}_6$ are orthotropic, tetragonal and orthotropic, with fixed symmetry axes ${\mathbf{e}}_i$ , respectively. Investigating the out-of-plane damage effect for the non-isotropic terms, one further finds the interesting property

$$\Delta L_i\left({\mathbf{e}}_3\right)=\Delta {\kappa }_i\left({\mathbf{e}}_3\right),\forall i=1,3,5,6$$

(59)

Together with a more strict interpretation of the damage mode separation, this can be used to formulate a simplified version of the proposed damage model where the out-of-plane damage effect due to microcracks is fully deactivated by setting ${\xi }_2={\xi }_4=0$ leading to the same simplified version as already proposed in (58).

To finish the interpretation of the proposed model, the symmetry of damage (Condition (14)) is investigated. Here it is directly concluded that $\Delta {\mathbb{C}}_2,\dots ,\Delta {\mathbb{C}}_6$ are isotropic, meaning the direction of damage D does not influence the damage effect at all. $\Delta {\mathbb{C}}_1$ is orthotropic with respect to damage. The symmetry axes are given by the eigenvectors of D, which is, as discussed above, also the minimal requirement that has to be satisfied.

To sum up this discussion, we can conclude that $\Delta {\mathbb{C}}_1$ is the most important part to describe the anisotropic damage effect. Beside the discussed special cases, it is also the only part of the damage effect that can fully degrade a certain symmetry if the principal directions of damage do not coincide with a principal direction of the undamaged material. This aspect will be further illustrated in the next section by investigating in two example materials.

Example of a ceramic matrix composite

To demonstrate the accuracy of the proposed model and compare it with an well established pre-existing model, an analytical homogenization procedure is used to provide virtual test data. In contrast to data derived from laboratory experiments, this has the advantage that the stiffness degradation is accessible for a broad range of damage states and loading directions that are not accessible in an experimental procedure. This gives more options for a validation and is sufficient to demonstrate the capabilities of the proposed model and compare it with a pre-existing model.

For the homogenization we resort to results from Goidescu et al., 2013. In there analytical closed form solutions are given assuming non-interacting (dilute concentration) microcracks in an orthotropic and homogeneous matrix. As boundary conditions either constant stress or strain boundary conditions can be used, which does not influence the results due to the assumption of a dilute concentration. Furthermore, the homogenization is limited to the 2D case, which does not significantly limit the validation process, since out-of-plane damage is assumed to play a minor role. Finally, as in the modeling process itself, we restrict the homogenization to the special case with no crack closure.

Let $\tilde{L}\left(\omega ,\mathbf{D}\right),L\left(\omega ,\mathbf{D}\right)$ be the in-plane elongation and let $\tilde{\kappa }\left(\omega ,\mathbf{D}\right),\kappa \left(\omega ,\mathbf{D}\right)$ be the bulk modulus in the direction ω with damage state D as predicted from the analytical model and the proposed damage effect model respectively. The material parameters $\mathit{\xi }$ can then be determined by minimizing the average error as

$$\mathit{\xi }=\arg \underset{{\xi }_i<0}{\min }\frac{1}{IJ}{\displaystyle \sum _{i,j=1}^{I,J}\left(\frac{\left|\tilde{L}\left({\omega }_i,{\mathbf{D}}_j\right)-L\left({\omega }_i,{\mathbf{D}}_j\right)\right|}{\tilde{L}\left({\omega }_i,{\mathbf{D}}_j\right)}+\frac{\left|\tilde{\kappa }\left({\omega }_i,{\mathbf{D}}_j\right)-\kappa \left({\omega }_i,{\mathbf{D}}_j\right)\right|}{\tilde{\kappa }\left({\omega }_i,{\mathbf{D}}_j\right)}\right)}.$$

(60)

In the optimization procedure the virtual test data is sampled by considering the four uni-axial damage states ${\mathbf{D}}_j=d_j{\mathbf{t}}_j\otimes {\mathbf{t}}_j$ , whereas the principal directions are $\alpha =0,\pi /2,\pi /4$ and $3/4\pi$ respectively and considering 50 loading directions ω_i, equally distributed from 0 to $2\pi$ . At this point it is important to note that both, the proposed damage model, and the homogenization procedure are linear in the damage state D. Due to this linearity, neither the combination of the above damage states, the case of isotropic damage with $\mathbf{D}=d\mathbf{1}$ , nor uni-axial damage states with the same principal direction, but different eigenvalues, contribute additional information and therefore does not need to be considered.

In what follows the fitting procedure is applied using different simplifications of the proposed model. These are (see also Table 5):

1. The full proposed model from (51) with the restriction ${\xi }_i\le 0$

2. The simplified in-plane model from (58) with ${\xi }_2={\xi }_4=0$

3. A model with the isotropic part simplified ${\xi }_2={\xi }_3={\xi }_4={\xi }_6\le 0$ and ${\xi }_5=0$

4. A pure anisotropic model with ${\xi }_2={\xi }_3={\xi }_4={\xi }_6=0$ and ${\xi }_1\le 0$

Additionally, the damage effect model from Maire and Chaboche (1997) and Chaboche and Maire (2001) are fitted. This fitting procedure is repeated for two examples of CMC materials, one with an initial undamaged tetragonal symmetry and one with an initial undamaged orthotropic symmetry. The corresponding material parameters are given in Table 4, the resulting modeling errors are summarized in Table 5.

Table 4.

Undamaged stiffness components in GPa used for the model validation, components are given in the (normalized) Kelvin-Mandel Notation.

	C 11	C 22	C 33	C 44	C 55	C 66	C 12	C 13	C 23
tetragonal example	144	C 11	73	C 55	58	72	38	16	C 13
orthotropic example (Mouchtachi et al., 2004)	217	278	84	81	76	135	59	13	34

Table 5.

Modeling errors for the two investigated examples together with different model simplifications.

		Tetragonal	Orthotropic
No.	Model	Example	Example
1	${\xi }_i<0$	4.4%	5.2%
2	${\xi }_2={\xi }_4=0$	4.4%	5.1%
3	${\xi }_2={\xi }_3={\xi }_4={\xi }_6,{\xi }_5=0$	5.1%	5.5%
4	${\xi }_2={\xi }_3={\xi }_4={\xi }_5={\xi }_6=0$	9.0%	8.5%
5	Reference model (Chaboche and Maire, 2001)	1.9%	2.3%

The first example is a balanced woven SiC-SiC ceramic matrix composite produced by Schunk Carbon-Technology in Heuchelheim, Germany. The corresponding elastic material constants, i.e. the irreducible material parameters of ${\mathbb{C}}^0$ , have been measured by an ultrasonic immersion technique at the Fraunhofer IWM and are summarized in Table 4. In this example the 0° and 90° direction are identical, leading to a tetragonal symmetry of the undamaged material. Figure 3 shows the undamaged properties in terms of the in-plane elongation and bulk modulus and the damaged state obtained from the homogenization approach. Furthermore, the latter is compared to a fitting using the two proposed extreme cases of the full model (Number 1 in Table 5) and the pure anisotropic simplification (Number 4 in Table 5) as well as the reference model.

Figure 3.

Damaged and undamaged elongation- (top) and bulk- (bottom) modulus for the two damage states with uni-axial damage with a principal direction of 0° (left) and 45° (right). Undamaged properties are plotted in green, damaged properties are presented as a result from the homogenization approach (dotted black) and compared to the proposed full model (red), the pure anisotropic model (dashed red) and the reference model from (Chaboche and Maire, 2001) (blue).

First of all, it is observed that both, the prediction from the homogenization model and all of the investigated damage effect models disturb the initial tetragonal symmetry and predict a orthotropic symmetry with symmetry axes that coincide with the principal direction ${\mathbf{t}}_i$ of the damage tensor D. Noting that the 0°, 45° and 90° directions are symmetry planes for tetragonal materials, this behaviour is expected from the discussion above and is attributed to the symmetry properties of $\Delta {\mathbb{C}}_1$ . However, it is important to note that this is also consistent with the analytical homogenization results.

Having a closer look at the different proposed model variants in terms of the modelling error (Table 5) helps to get further insight. First of all, the modelling error of the second model simplification (the in-plane model) is of the same order of magnitude as the full model which confirms the previous findings that $\Delta {\mathbb{C}}_2$ and $\Delta {\mathbb{C}}_4$ are indeed redundant for the in-plane damage behaviour. Compared to these two, the pure anisotropic model, where only $\Delta {\mathbb{C}}_1$ is active, leads to a higher modelling error. Figure 3 reveals the reason for this: Even though the accuracy of the anisotropic model is of the same order (or even better) for the elongation modulus, it completely fails to describe the bulk modulus compared to the full model. Finally, as intended, the model with a simplified isotropic part (Model 3 in Table 5) yields a very good compromise between complexity (2 material parameters) and accuracy.

Before investigating the reference model in detail, a second example is presented. This is a non-balanced SiC-SiC material with properties varying from warp to weft direction and thus present an initially orthotropic material. Undamaged material parameters are given in Table 4. Again, results from the homogenization together with the fitting of the two proposed extreme models are presented in Figure 4. In contrast to the previous case with an initially tetragonal symmetry, it is now necessary to investigate in all three damage directions with ${\mathbf{t}}_1$ having an angle of 0°, 45° and 90°. However, for the reasons of clarity and since the $90°$ direction does not reveal any additional findings, only the first two directions are presented in Figure 4. With reference to the previous example, it is now observed from Figure 4(b) that the symmetry can be degraded to a triclinic symmetry in the case of an uniaxial damage in the $45°$ direction (i.e $\mathbf{Q}=-\mathbf{1}$ being the only element in the symmetry group χ). Again this is attributed to the principal damage direction that does not coincide with any symmetry axis of the undamaged material and most importantly is consistent with the analytical homogenization approach assuming a dilute crack concentration in an homogeneous (matrix) material. Beside this aspect, even though the modelling error is slightly higher in that case (Table 5), this second example confirms the previous finding both qualitatively and quantitatively.

Figure 4.

The same as in Figure 4 but for the second investigated orthotropic Material: Damaged and undamaged elongation- (top) and bulk- (bottom) modulus for the two damage states with uni-axial damage with a principal direction of $0°$ (left) and $45°$ (right). Undamaged properties are plotted in green, damaged properties are presented as a result from the homogenization approach (dotted black) and compared to the proposed full model (red), the pure anisotropic model (dashed red) and the reference model from (Chaboche and Maire, 2001) (blue).

We finally turn to the comparison with the referece model presented in Maire and Chaboche (1997) and Chaboche and Maire, (2001). The damage effect of this model is described by the equation

$$\mathbb{C}\left(\mathbf{D}\right)=\frac{1}{2}\left(\left[11-\mathbb{D}\left(\mathbf{D}\right)\right]:{\mathbb{C}}^0+{\mathbb{C}}^0:\left[11-\mathbb{D}\left(\mathbf{D}\right)\right]\right).$$

(61)

with D being the second rank damage tensor entering the formulation by

$$\mathbb{D}\left(\mathbf{D}\right)=\frac{\xi }{2}\left(\mathbf{1}\otimes \mathbf{D}+\mathbf{D}\otimes \mathbf{1}\right)+\frac{1-\xi }{2}\left(\mathbf{1}\underset{¯}{\otimes }\mathbf{D}+\mathbf{1}\underset{¯}{\otimes }\mathbf{D}+\mathbf{D}\underset{¯}{\otimes }\mathbf{1}+\mathbf{D}\underset{¯}{\otimes }\mathbf{1}\right).$$

(62)

The material parameter ξ is supposed to be in the range $\left[0,1\right]$ . As can be seen from the above references, this tensorial damage model is a promising candidate to model the damage behavior of CMCs if it is supplemented by additional scalar damage variables. This is also verified by the present investigation where the fitting of equation (61) shown in Figures 3 and 4 is of similar quality as the proposed model, and this without considering the additional scalar variables. Moreover, the modelling error is even slightly lower compared to the proposed models (see Table 5).

However beside this high accuracy, the Chaboche Model inherents a fundamental deficit: It does not fullfill the damage growth criterion. This can can easily be proved for the simple isotropic case, but also holds for the investigated anisotropic cases. To illustrate this, Figure 5 presents the value of the first two eigenvalues of the conjugate Forces Y. It is cleary seen that these become negative for a broad range of in-plane strain states with the out-of-plane strain set to its equilibrium state for a uni-axial tension, i.e. ${\epsilon }_{33}=-{\nu }_{23}{\epsilon }_{22}$ .

Figure 5.

Numerical evaluation of the damage growth criterion for the damage effect tensor given in (61) in terms of the first eigenvalues y₁ (top) and the second eigenvalue y₂ (bottom) for the two investigated materials with tetragonal undamaged symmetry (left) and orthotropic undamaged undamaged symmetry (right) showing the region with a negative damage growth rate.

To demonstrate the significance of the positive damage growth rate for the critical strain regions, the damage dependent modulus $E_2\left(\mathbf{D}\right)$ of an uni-axial tension experiment in the ${\mathbf{e}}_2$ -directions, i.e.

$$E_2\left(\mathbf{D}\right)=C_{22}\left(\mathbf{D}\right)-{\nu }_{21}{C}_{21}\left(\mathbf{D}\right)-{\nu }_{23}\left(\mathbf{D}\right)C_{23}\left(\mathbf{D}\right)$$

(63)

is calculated explicitly for a bi-axial damage state with principal directions ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ . For the first example material the modulus is given by

$$E_2\left(\mathbf{D}\right)=C_{22}^0+5.5d_1-129d_2,$$

(64)

whereas for the second example, the orthotropic material, it is given by

$$E_2\left(\mathbf{D}\right)=C_{22}^0+10.9d_1-251d_2.$$

(65)

In both cases it is clearly seen how a damage in the ${\mathbf{e}}_1$ directions increases the modulus in the ${\mathbf{e}}_2$ direction. This is not necessarily observable in the simulation of an uni-axial tension experiment, since the damage is expected to be mainly orientated in the direction of the loading. However it can become relevant if bi-axial tension-internal pressure experiments are investigated. These experiments are often chosen to characterize the anisotropic damage effect and damage evolution (Bernachy-Barbe et al., 2015; Chaboche and Maire, 2001). Therefore, the same situation is observable from Figure 2 in Chaboche and Maire (2001) where tension/internal pressure experiments are used to fit a model based on equation (61) and the model clearly predicts an increasing modulus E₂ while the tube is loaded with an axial stress σ₁. Bernachy-Barbe (a) and (b).

Clearly, it can be argued that the increasing transverse Young’s modulus due to the violation of the damage growth criterion is of minor importance. However, the fact that it plays a minor role is related to the small amount of transverse damage effect in the present cases. Therefore the violation of the damage growth criterion will become more important for materials with a larger transverse damage effect. On the other side, the increasing stiffness can also become relevant in the simulation of components with strain states being in the critical regions highlighted in Figure 5. As an additional aspect, fulfilling the damage growth criterion helps in formulating the damage growth criterion since the second law is automatically satisfied even if the dissipation potential becomes non-convex. The latter is often satisfied by formulating the dissipation potential as a function of the positive definite part of the conjugate forces only (e.g. Maire and Chaboche, 1997; Chaboche and Maire, 2001; Chaboche and Maire, 2002). Realising such a positive definite decomposition is however quite complex since it involves a spectral decomposition as well as the heaviside distribution, which are both not easy to handle numerically. To prevent these difficulties, the proposed model is highly advantageous since it a priori satisfies the damage growth criterion irrespective of the chosen material parameters, loading scenarios and damage evolution.

Damage deactivation

As a final aspect, it is next demonstrated that the proposed model can systematically be supplemented with a damage deactivation for “compression”. The challenge in modelling damage deactivation is to alter the strain energy density for compressive stress or strain states in a physically reasonable manner such that both, the resulting stress-strain and the conjugate forces do not posses any discontinuity. I.e. the strain energy density must possess at least a C¹ continuity with respect to its arguments $\mathit{\epsilon }$ and D.

To include damage deactivation, we first of all note, that the damage growth Criterion (45) already involves the case of closed microcracks by the unilateral condition $[[\mathbf{u}]]·\mathbf{n}\ge 0,\forall x\in \Gamma$ in equation (43) and is thus also expected to be valid for closed microcracks.

To include the damage deactivation, a deactivation (or crack closure) criterion is introduced. With the damage tensor D in its spectral representation

$$\mathbf{D}=d_1{\mathbf{t}}_1\otimes {\mathbf{t}}_1+d_2{\mathbf{t}}_2\otimes {\mathbf{t}}_2,$$

(66)

it is proposed to use the two deactivation criteria

$$g_1\left(\mathit{\epsilon },\mathbf{D}\right)={\mathbf{t}}_1\mathit{\epsilon }{\mathbf{t}}_1,$$

(67)

$$g_2\left(\mathit{\epsilon },\mathbf{D}\right)={\mathbf{t}}_2\mathit{\epsilon }{\mathbf{t}}_2,$$

(68)

where the first (second) crack, associated to d₁ (d₂) is closed if $g_1<0$ ( $g_2<0$ ), i.e. if there is a compressive normal strain to the crack face. Otherwise the damage is assumed to be fully active. To simplify the following derivation, only the deactivation with respect to the first eigenvalue d₁ is considered subsequently. The generalization with respect to the second eigenvalue d₂ will be given thereafter. Associated with the deactivation criterion (67), the strain energy density (1) is modified as follows

$$\begin{array}{c}\phi \left(\mathit{\epsilon },\mathbf{D}\right)=\{\begin{array}{ll}{\phi }^0\left(\mathit{\epsilon }\right)+{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right),\hfill & g_1\left(\mathit{\epsilon },\mathbf{D}\right)>0\hfill \\ {\phi }^0\left(\mathit{\epsilon }\right)+{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)+{\phi }_1^{-}\left(\mathit{\epsilon },\mathbf{D}\right),\hfill & g_1\left(\mathit{\epsilon },\mathbf{D}\right)\le 0\hfill \end{array}=:\{\begin{array}{cc}{\phi }_I& g_1\left(\mathit{\epsilon },\mathbf{D}\right)>0\\ {\phi }_{II}& g_1\left(\mathit{\epsilon },\mathbf{D}\right)\le 0\end{array}\end{array}$$

(69)

where ${\phi }^0$ and ${\phi }^D$ are the previously derived undamaged and damaged parts for a fully active damage and ${\phi }_1^{-}\left(\mathit{\epsilon }\right)$ is introduced to describe the damage deactivation. Following similar arguments as in Curnier et al., (1995), Welemane (2002), and Cormery and Welemane (2010) the strain energy density (69) is continuously differentiable across the deactivation surface $\mathcal{A}=\left\{\mathit{\epsilon },\mathbf{D}:g_1\left(\mathit{\epsilon },\mathbf{D}\right)=0\right\}$ if $\exists \left(\mathit{\epsilon },\mathbf{D}\right)\in \mathcal{A}$ such that

$${\phi }_I\left(\mathit{\epsilon },\mathbf{D}\right)={\phi }_{II}\left(\mathit{\epsilon },\mathbf{D}\right),$$

(70)

$${\partial }_{\mathit{\epsilon }}{\phi }_I\left(\mathit{\epsilon },\mathbf{D}\right)={\partial }_{\mathit{\epsilon }}{\phi }_{II}\left(\mathit{\epsilon },\mathbf{D}\right),$$

(71)

$${\partial }_{\mathbf{D}}{\phi }_I\left(\mathit{\epsilon },\mathbf{D}\right)={\partial }_{\mathbf{D}}{\phi }_{II}\left(\mathit{\epsilon },\mathbf{D}\right),$$

(72)

and $\forall \left(\mathit{\epsilon },\mathbf{D}\right)\in \mathcal{A}$ we have:

$$-{\partial }_{\mathbf{D}\mathit{\epsilon }}^2{\phi }_1^{-}=s_1\left(\mathit{\epsilon },\mathbf{D}\right){\partial }_{\mathbf{D}}{g}_1\otimes {\partial }_{\mathit{\epsilon }}{g}_1,$$

(73)

$$-{\partial }_{\mathit{\epsilon }\mathbf{D}}^2{\phi }_1^{-}=s_1\left(\mathit{\epsilon },\mathbf{D}\right){\partial }_{\mathit{\epsilon }}{g}_1\otimes {\partial }_{\mathbf{D}}{g}_1,$$

(74)

$$-{\partial }_{\mathbf{DD}}^2{\phi }_1^{-}=s_1\left(\mathit{\epsilon },\mathbf{D}\right){\partial }_{\mathbf{D}}{g}_1\otimes {\partial }_{\mathbf{D}}{g}_1,$$

(75)

$$-{\partial }_{\mathit{\epsilon }\mathit{\epsilon }}^2{\phi }_1^{-}=s_1\left(\mathit{\epsilon },\mathbf{D}\right){\partial }_{\mathit{\epsilon }}{g}_1\otimes {\partial }_{\mathit{\epsilon }}{g}_1.$$

(76)

In here $s\_1\left(\mathit{\epsilon },\mathbf{D}\right)$ is an arbitrary scalar valued C⁰ function. Conditions (70) to (72) are satisfied for $\mathit{\epsilon }=\mathbf{D}=\mathbf{0}$ if the principles applied to ${\phi }_I$ are transferred to ${\phi }_{II}$ , i.e. ${\phi }_{II}$ is linear in D and quadratic in $\mathit{\epsilon }$ . Based on these properties it can further be shown that Condition (76) is sufficient to satisfy Conditions (73) to (75). Evaluating (76) for (69) together with (67) leads

$${\partial }_{\mathit{\epsilon }\mathit{\epsilon }}^2{\phi }^{-}=-s_1\left(\mathit{\epsilon },\mathbf{D}\right){\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1.$$

(77)

However, since ${\phi }^{-}$ is assumed to be only quadratic in $\mathit{\epsilon }$ we have $s\_1=s\_1\left(\mathbf{D}\right)$ and integration leads

$${\phi }^{-}\left(\mathit{\epsilon },\mathbf{D}\right)=-s_1\left(\mathbf{D}\right)\text{tr}\left(\mathit{\epsilon }·{\mathbf{t}}_1\otimes {\mathbf{t}}_1\right)$$

(78)

By this, the choice of the damage deactivation function has been successfully reduced to the choice of $s\_1\left(\mathbf{D}\right)$ . To go further we can next resort to experimental and analytical results for the damage behaviour of ceramic composites and require that the damage effect normal to the crack face is fully deactivated if the crack is under compression:

$${\partial }_{\mathit{\epsilon }\mathit{\epsilon }}^2\phi \left(\mathit{\epsilon },\mathbf{D}\right)::{\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1={\mathbb{C}}^0::{\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1\otimes {\mathbf{t}}_1$$

(79)

This is satisfied if and only if

$$s_1\left(\mathbf{D}\right)=2{\xi }_1{d}_1+2\text{tr}\left(\mathbf{D}\right){\Sigma }_{j=2}^6{\xi }_j$$

(80)

which complements the formulation of the damage deactivation. Following an analogous procedure for the second eigenvalue d₂ the proposed model can be supplemented by damage deactivation as follows

$$\begin{array}{l}\phi \left(\mathit{\epsilon },\mathbf{D}\right)={\phi }^0\left(\mathit{\epsilon }\right)+{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)-{\displaystyle \sum _{i=1,2}H}\left(-g_i\left(\mathit{\epsilon },\mathbf{D}\right)\right)s_i\left(\mathbf{D}\right)\text{tr}\left(\mathit{\epsilon }·{\mathbf{t}}_i\otimes {\mathbf{t}}_i\right)\end{array}$$

(81)

where H is the heaviside distribution and

$$s_i\left(\mathbf{D}\right)=2{\xi }_1{d}_i+2\text{tr}\left(\mathbf{D}\right){\Sigma }_{j=2}^6{\xi }_j.$$

(82)

Finally, according to the spectral representation (66) two special cases must be discussed in the following. If the eigenvalues are the same d₁ = d₂, no unique spectral representation exists and an arbitrary, orthonormal $\left\{{\mathbf{t}}_1,{\mathbf{t}}_2\right\}$ may be chosen. The mechanical interpretation of this case is, that there is no preferential crack orientation. In that case it is suggested to keep the direction from a previous time $t-ϵ$ , if it exists. If, instead, D is isotropic from the beginning, a feasible choice would be to chose ${\mathbf{t}}_i$ as the principal material axes ${\mathbf{e}}_i$ , to have a coordinate system independent behaviour.

To illustrate the damage deactivation behaviour, Figure 6 shows both the elongation and the bulk modulus for fully active ( $g_1>0,g_2>0$ ), partially deactivated ( $g_1<0,g_2>0$ ) and fully deactivated ( $g_1<0,g_2<0$ ) damage for two principal directions $\alpha =0°,45°$ . First of all it is clearly seen, that in case of the partial deactivation the damaged elongation modulus in the direction of the closed microcrack is fully restored to the undamaged value. Transversal to the closed crack, the damage effect is in turn fully active if the transversal direction is active. As expected, the elongation modulus is fully restored in both principal directions if both damage components are deactivated. However, in between the principal damage directions, some amount of damage remains. The latter can be attributed to the shear components and the fact that crack faces are allowed to slip in the present formulation even if the crack is under compression. The behaviour of the bulk modulus is very similar with the difference that the bulk modulus is not fully restored in any direction.

Figure 6.

Illustration of the damage deactivation by the orthotropic example material showing the undamed state (black), fully active damage (red), partially deactivated (blue) and fully deactivated (green) states in terms of the elongation modulus $L\left(\omega \right)$ (top) and bulk modulus (bottom) for a damage state with $d_1=2d_2=0.05$ and a principal direction of $\alpha =0°$ (left) and $\alpha =45°$ (right).

Summary and discussion

This paper presents a systematic way to derive a damage effect model with a special focus on the damage behaviour of initially anisotropic composite materials. The result is a compact formulation of the damage effect, which is, by a thorough interpretation of the model parameters, further simplified for special cases. A validation using results from a homogenization approach and a comparison with a well established reference model from the literature shows the good capabilities of the proposed model, even though the reference model results in a slightly smaller modelling error. The latter however comes with the drawback of not being mechanically consistent with the damage growth criterion, which can lead to the unphysical behaviour of an increasing stiffness caused by damage. On the other hand this can lead to problems in the formulation of the damage growth equation by a convex and thermodynamical consistent dissipation potential. This for example seemed to be a challenge in the model from Maire and Chaboche, 1997, where using the positive definite part of the conjugate forces in the dissipation potential requires a spectral decomposition. Such a positive definite decomposition is numerically difficult to handle since it involves the heaviside distribution and the calculation of eigensystems and their derivatives which are furthermore not uniquely defined in some special cases. These difficulties are a priori prevented by the present formulation. Finally, as an add-on, the model has been extended with a mathematically consistent damage deactivation describing the crack closure effect under compressive states.

Four main ingredients are used for the formulation: The symmetry class of the initial undamaged material, a specification of the damage state to in-plane damage, the mechanically consistency in terms of the damage growth criterion and a plausible limitation of the order of the relevant variables. With respect to theses assumptions, the question arises: what could be done differently to further improve the model? First of all, the chosen set of invariants (17) and subsequently also (41) is certainly not unique. As was also pointed out in (Zheng, 1994), other choices are possible. In fact it is easy to construct an alternative set by substituting any I_i and I_j with $I_i+I_j$ and $I_i-I_j$ , giving many alternative options. This fact would be irrelevant if condition (48) would be applied directly. However, since only the sufficient condition (49) is applied here, choosing other sets of invariants would result in a different model. This directly leads to the next aspect: applying the necessary condition (48) instead of the sufficient condition (49), would help in improving the modelling accuracy, however would also lead to a much more complex model with more complex constraints on the material parameters.

Next, most obviously, the polynomial order in the damage variable D could be increased. Instead of (20) a quadratic damage effect would result in the much more complex damage effect function

$$\begin{array}{l}{\phi }^D\left(\mathit{\epsilon },\mathbf{D}\right)=\xi I^g+{\xi }_{ij}{I}_i^a{I}_j^f+{\xi }_{\mathit{ijk}}{I}_i^a{I}_j^a{I}_k^f+{\xi }_{\mathit{ijk}}{I}_i^a{I}_j^c{I}_k^c\\ +{\xi }_{\mathit{ijk}}{I}_i^a{I}_j^d{I}_k^c+{\xi }_{ij}{I}_i^b{I}_j^f+{\xi }_{\mathit{ijk}}{I}_i^b{I}_j^c{I}_k^c+{\xi }_{ij}{I}_i^c{I}_j^d\\ +{\xi }_{\mathit{ijk}}{I}_i^c{I}_j^c{I}_k^e+{\xi }_{\mathit{ijkl}}{I}_i^a{I}_j^a{I}_k^c{I}_l^c\end{array}$$

(83)

involving 139 material parameters. Beside this non-practical dimension, it would have helped only a little bit to decrease the modelling error in the presented examples. The reason being that the applied homogenization, assuming a dilute crack distribution in an homogeneous matrix, also results in a linear damage effect. Therefore, increasing the order of the approximating model might help in approximating the used data, however would lead to a large error for other cases with other damage states, say $\lambda \mathbf{D}$ which scale linear in the homogenization result, but includes quadratic parts ${\lambda }^2{\mathbf{D}}^2$ in the approximating model leading to large modelling errors. From this aspect the conclusion is drawn that the first improvement, i.e. generalizing the way the damage growth criterion is handled, should be tried first.

This conclusion is only valid under the assumption of a homogeneous matrix, whereas matrix here refers to the initial undamaged material and not the matrix phase of the composite material. This assumption is also valid for a cracked composite material as long as the scales are sufficiently separated, i.e. the crack length must be much smaller than a characteristic length of the microstructure, e.g. the fibre diameter. This however will only be valid for a very limited damage phase at the very beginning. As soon as the cracks reach the size of the microstructure, interactions between the cracks and the microstructure prohibits the assumption of an homogeneous matrix material. It is then also expected that the damage behaviour becomes more complex where for instance quadratic interaction terms like $d_1{d}_2$ might indeed help to increase the modelling accuracy. With regard to this aspect, computational homogenization on cracked inhomogeneous composite materials might give some insights and would help in further validating and improving the proposed model.

Finally the proposed damage deactivation, even though being plausible so far, needs to be further validated. Special focus should be paid to the behaviour under combined compressive-shear loads. As shown experimentally (see Chaboche and Maire (2002)) a shear deformation superimposed on a compressive load can lead to stiffening effects related to the friction of the crack faces, this in turn can also be associated with a small amount of plastic deformation. Introducing a so-called stored strain energy as in Chaboche and Maire (2002) which is related to the crack face friction and is released upon crack opening, seems to be a consistent way, even though the corresponding model implementation becomes quite complex.