A Thermodynamically Consistent Model of Quasibrittle Elastic Damaged Materials Based on a Novel Helmholtz Potential and Dissipation Function

Abstract

In the paper, a thermodynamically consistent model of elastic damaged material in the framework of small strain theory is formulated, describing the process of deterioration in quasibrittle materials, concrete in particular. The main goal is to appropriately depict the distinction between material responses in tension and compression. A novel Helmholtz energy and a dissipation potential including three damage parameters are introduced. The Helmholtz function has a continuous first derivative with respect to strain tensor. Based on the assumed functions, the strain–stress relationship, the damage condition, the evolution laws, and the tangent stiffness tensor are derived. The model’s predictions for uniaxial tension, uniaxial compression, uniaxial cyclic compression–tension, and pure shear tests are calculated using Wolfram Mathematica in order to identify the main features of the model and to grasp the physical meaning of an isotropic damage parameter, a tensile damage parameter, and a compressive damage parameter. Their values can be directly bound to changes of secant stiffness and generalized Poisson’s ratio. An interpretation of damage parameters in association with three mechanisms of damage is given. The considered dissipation potential allows a flexible choice of a damage condition. The influence of material parameters included in dissipation function on damage mode interaction is discussed.

1. Introduction

Many microscopic cavities, microcracks, flaws, and so forth can be observed on the microscale in solid materials. Those defects owe their origins to disturbances in the interatomic or intermolecular bonds. The technological casting process or loading leads to the development of existing micro- and mesoscale structural irregularities as well as to the initiation of new breakages in material structure. The evolution of those defects causes fracture of a material together with deterioration of its mechanical properties called a damage phenomenon. Damage can occur on the atomic scale, microscale, and mesoscale, leading to macroscale effects. The physical mechanisms of damage include cleavage (decohesion), growth, and coalescence of microvoids, glide plane decohesion, and grain boundary decohesion. Deterioration of mechanical properties mainly depends on the type of material and loading conditions. Ductile damage can be observed in metallic materials, while brittle damage appears in rocks, concrete, ceramics, and composite materials. At an elevated temperature, polycrystalline metals exhibit creep damage. Majority of materials deteriorate due to cyclic loading. Damage and fracture processes ongoing in a material treated as a macroscopic material continuum are successfully described in the framework of continuum damage mechanics [1,2,3,4].

The growing complexity of structures demands improvements of existing constitutive models of quasibrittle materials, concrete in particular. An accurate description of their behavior is a subject of an ongoing research [5,6,7,8,9,10,11,12,13,14,15,16]. Designing complex and innovative structures requires a sound understanding of the mechanical behavior of concrete. One of the crucial characteristics of concrete (and many other quasibrittle materials) is its low tensile strength, which allows tensile cracking at very low stresses compared with stresses in the compression range. This tensile cracking significantly reduces the stiffness of structural elements. Nucleation, growth, and coalescence of microcracks present in concrete lead to degradation of stiffness and irreversible deformations. Therefore, advanced and accurate constitutive modeling of concrete under complex loading conditions is of keen interest from the point of view of rational prediction of damage and failure.

In phenomenological constitutive modeling, two main aspects are of utmost importance: the thermodynamic consistency of the model and its compatibility with experimental data [6,13,17,18,19,20,21,22]. In the article, we are predominately focused on the first aspect, although we reference to basic experiments as well. In the framework of small strain, we propose a model of an elastic damaged material based on novel Helmholtz and dissipation functions that both satisfy the laws of thermodynamics and capture the crucial difference between tensile and compressive behavior of materials. The fundamentals of thermodynamically consistent constitutive modeling based on two potentials were laid in the last four decades and have well consolidated since [1,2,7,14,17,18,19,20]. The particulars of the framework for geotechnical materials have been improved upon [19,20], but the actual usage is limited due to mathematical difficulties, especially performing the Fenchel–Legendre transformation, although some examples are present in the appropriate literature [6,7,8,9,16,21]. In the paper, we introduce two potentials: a Helmholtz free energy and a dissipation and use the discussed above approach to derive relations that can be directly applied in numeric calculations and do not diverge from the structure of commonly used models.

Quasibrittle materials are characterized by two features: failure is caused mainly by fracture rather than plastic yield, and “the fracture front is surrounded by a large fracture-process zone in which progressive distributed cracking or other damage takes place” [5]. We focus on concrete as the main representative of quasibrittle materials’ class extensively used in civil engineering, but the regarded setting also includes rocks [19] and some alloys [23,24]. In the rough approximation, the degradation of its properties can be seen as isotropic and thus described by one damage parameter, $d$, which reflects the decrease in Young’s modulus (or an actual decrease in a cross section carrying loading). As a result, the Helmholtz function is the strain energy of an isotropic body multiplied by factor $1-d$ [1,2,3]. However, concrete subjected to a sufficiently high level of loading exhibits major differences between the behavior in tension and compression due to complex interactions between its components [25,26,27,28,29,30,31,32]. Therefore, at least two different damage parameters are needed, connected to positive and negative parts of the strain (or stress) tensor [33]. This idea is widely used (see, for example, [1,2,4,11,33,34,35,36]) and well grounded in experiments, but it is cumbersome to employ it in the numerical calculations due to the discontinuous relation between the stress tensor and the strain tensor for the uniaxial tests. Among other approaches, there are concepts to use the fully anisotropic models [10,13,15,16,21,22,37], to bind the Helmholtz function with the direction of emerging cracks [9], or to connect damage parameters to the split of energy into bulk and distortional parts [4]. We make use of the last idea but introduce three damage parameters instead of two: one describing isotropic damage ($d$) and two other, bound to the volumetric part of the Helmholtz function, separating the degradation into tension ($d_T$) and compression ($d_C$). This assumption allows for describing the basic features of concrete’s behavior.

Evolution laws of damage parameters are often making degradation growth dependent on parts of strain energy [1,2], while damage conditions are similar or, in models including plasticity, the same as yield conditions [13,22,34,35,36,38,39]. We propose a dissipation function that allows for including those key features. The resultant damage condition conjugated to the introduced dissipation is the one presented in [40], inspired by the Ottosen’s shape function [41].

After this introduction, the paper is organized as follows. As a starting point, a brief summary of the framework of constitutive modeling of thermodynamically consistent rate-independent materials based on [7,14,19,20] is given in Section 2. In the main part of the paper (Section 3), we introduce two potentials: a dissipation and a novel Helmholtz function. We employ three scalar damage parameters, which improves the model’s flexibility and allows for overcoming differentiability issues. Next, using the adopted framework, we derive all the relations necessary for the numerical implementation. All material parameters included in the damage condition are given by closed analytical formulae. Then, in Section 4 we show the results of calculations for simple tests: the uniaxial compression, uniaxial tension, uniaxial cyclic tension–compression, and simple shear, which allow for understanding the meaning of the introduced damage parameters and parameters involved in the dissipation potential. The presented outcomes are obtained by direct numerical solving of sets of the model differential equations: calculations are performed using Wolfram Mathematica. In Section 5, we interpret three damage parameters in reference to the physical course of deterioration. The influence of material parameters, included in the definition of dissipation function, on damage mode interaction is discussed. The most relevant research outcomes and conclusions are drawn in Section 6.

2. Framework for Thermodynamically Consistent Constitutive Modeling of Elastic Damaged Materials Using Helmholtz Energy and Dissipation Function

A thermodynamically consistent material model is a model that ensures the satisfaction of the first and the second law of thermodynamics in the form of the Clausius–Duhem inequality [18]. Basically, it states that dissipation in any real process, reversible or not, has to be non-negative. In order to fulfill this condition, using an approach described in [7,14,19,20], we have to obey the below-described algorithm of consecutive steps in the design of a material’s model.

The employed framework for elastic-plastic materials is described in [19,20] with extension to damaged materials in [7,14]. It is quite general and allows any of four energy potentials (Helmholtz or Gibbs or internal energy or enthalpy) and either dissipation or damage condition to be the initial points of the formulation. It admits tensorial internal damage variables. The presented version of the framework is reduced to suit the considered problem: the equations are narrowed down to the case of elasticity and damage described by a Helmholtz potential and a dissipation function dependent on multiple scalar damage parameters. Based on this kinematic description of damage, the goal is to obtain the following equations: a strain–stress relationship, a damage condition, evolution equations, and a definition of the fourth-order stiffness tensor of an elastic damaged material.

The first and the most difficult step of the scheme is to assume two appropriate functions: a Helmholtz potential and a dissipation. This requires careful consideration regarding both mathematical properties, such as convexity, and the ability to describe the real behavior of material.

The Helmholtz energy $H\left(\mathsf{\epsilon },d_i\right)$ must be a function that is strictly convex with respect to strain tensor $\mathsf{\epsilon }$, convex with respect to damage parameters $d_i$ ($i=1,2,..,N$, where $N$ is the number of damage parameters in a considered model), and strictly non-negative with respect to $\mathsf{\epsilon }$, meaning that it is null for $\mathsf{\epsilon }=0$ and positive elsewhere. On the other hand, dissipation $D\left({\dot{d}}_i\right)$ is a convex, strictly non-negative, and, for rate-independent materials, a homogenous of degree 1 function of rates, ${\dot{d}}_i\left(t\right)$, where $t$ is time. It can also have non-rate-type arguments; for example, when hardening is involved, it can additionally depend on the current values of damage parameters $D\left({\dot{d}}_i;d_i\right)$ (compare in [7,14,19,20]).

Now, the Cauchy stress tensor $\mathsf{\sigma }$ and the generalized stresses ${\sigma }_{di}$ can be calculated, differentiating the Helmholtz potential:

$$\mathsf{\sigma }=\frac{\partial H}{\partial \mathsf{\epsilon }},\hspace{1em}\hspace{1em}{\sigma }_{di}=-\frac{\partial H}{\partial d_i}$$

(1)

and next, the dissipative stresses related to the dissipation function can be obtained using the following formula:

$${\overline{\sigma }}_{di}=\frac{\partial D}{\partial {\dot{d}}_i}.$$

(2)

Dissipation, being a homogenous function of degree 1, satisfies the suitable Euler’s theorem of the form:

$$D\left({\dot{d}}_i\right)={\displaystyle \sum _{i=1}^N\frac{\partial H}{\partial {\dot{d}}_i}{\dot{d}}_i}={\displaystyle \sum _{i=1}^N{\overline{\sigma }}_{di}{\dot{d}}_i}.$$

(3)

As a result, the following Fenchel–Legendre transformation binds the dissipation and the conjugated damage function $\overline{F}\left({\overline{\sigma }}_{di}\right)$:

$$D\left({\dot{d}}_i\right)+\lambda \hspace{0.17em}\overline{F}\left({\overline{\sigma }}_{di}\right)={\displaystyle \sum _{i=1}^N{\overline{\sigma }}_{di}{\dot{d}}_i}.$$

(4)

Like in the plastic flow theory, this is equivalent to:

$$\lambda \hspace{0.17em}\overline{F}\left({\overline{\sigma }}_{di}\right)=0\hspace{1em}\mathrm{and}\hspace{1em}\lambda \ge 0\hspace{1em}\mathrm{and}\hspace{1em}\overline{F}\left({\overline{\sigma }}_{di}\right)\le 0,$$

(5)

$\lambda$ being a damage multiplier. In turn, if degradation is evolving, the following damage condition occurs:

$$\overline{F}\left({\overline{\sigma }}_{di}\right)=0$$

(6)

and $\lambda >0$ during the purely elastic process, $\lambda =0$, and $\overline{F}\left({\overline{\sigma }}_{di}\right)<0$. The damage condition (6) needs to be complemented by the consistency condition:

$$\dot{\overline{F}}\left({\overline{\sigma }}_{di}\right)=0,$$

(7)

where the upper dot denotes derivative with respect to time. The laws of evolution of damage parameters are associated with Equation (6), that is,

$${\dot{d}}_i=\lambda \frac{\partial \hspace{0.17em}\overline{F}}{\partial \hspace{0.17em}{\overline{\sigma }}_{di}}.$$

(8)

It needs to be stressed that, in general, performing transformation described by Equation (4) is not a trivial task. In this paper, we will use a quite simple dissipation function, but getting the conjugated function, still, is not an “automatic” process.

At this point, we have not established a connection between the two potentials yet. So far, we have defined two independent sets of the generalized stresses in Equation (1) resulting from the Helmholtz energy and the dissipative stresses in Equation (2) based on the dissipation function. The role of the connector is played by the orthogonality principle of Ziegler [17,19,20]. The principle is based on equalization of the dissipation due to the rate of the Helmholtz energy and the proposed dissipation potential. Consequently, it states that the generalized and the respective dissipative stresses are equal:

$${\overline{\sigma }}_{di}={\sigma }_{di}.$$

(9)

Equation (9) is a categorizing hypothesis: it suits a wide class of engineering materials (compare in [19]). The presented framework for constitutive modelling gives enough flexibility in the selection of appropriate Helmholtz potential and dissipation function.

To derive the incremental form of the constitutive relationships and calibrate the model, it is convenient to transform the above obtained equations to the Cauchy stress space. As Equation (9) is in force, we can express the damage condition (Equation (6)), the consistency condition (Equation (7)), and the evolution laws (Equation (8)) via the generalized stresses and, in turn through Equation (1), via the Cauchy stress tensor and damage parameters, that is,

$$F\left(\mathsf{\sigma }\right)=0,\hspace{1em}\hspace{1em}\dot{F}\left(\mathsf{\sigma }\right)=0,\hspace{1em}{\dot{d}}_i=\lambda \hspace{0.17em}{f}_i\left(\mathsf{\sigma }\right).$$

(10)

In general, case functions $f_i$ are not bound directly to $F$ (are not its derivatives), since in the space of generalized stresses, the evolution laws do not have to be associated.

Hardening can be easily incorporated by making material constants dependent on the current values of damage parameters, which will be discussed briefly in the subsequent sections.

For numerical implementation, it is a vital part of the formulation to provide components of the fourth-order tensor of tangent stiffness. The procedure is classic: from the consistency condition included in Equation (10), we obtain multiplier $\lambda$ and incorporate it in the formula for stress increment. The details are given in [18,20].

3. Formulation of Elastic Damaged Material Model

In this section, a model of an elastic damaged material is presented. We propose two novel governing functions—a Helmholtz potential and a dissipation—and obtain all the useful relations following the steps described in Section 2.

3.1. Helmholtz Energy

Let us define the cylindrical invariants of the strain tensor:

$$p=\frac{1}{\sqrt{3}}\mathrm{tr}\mathsf{\epsilon },\hspace{1em}q=\sqrt{\mathrm{tr}{e}^2},$$

(11)

where $\mathrm{tr}$ denotes the trace of a second-order tensor and $e=\mathsf{\epsilon }-\frac{1}{3}\left(\mathrm{tr}\mathsf{\epsilon }\right)I$ is the strain tensor’s deviator, $I$ being the second-order unit tensor. Now we can introduce the following Helmholtz function:

$$H\left(\mathsf{\epsilon },d,d_T,d_C\right)=\left(1-d\right)\left[\frac{3}{2}{K}_M\left(d_T,d_C\right)p^2+\frac{3}{2}{K}_D\left(d_T,d_C\right)p\sqrt{p^2}+G\hspace{0.17em}{q}^2\right].$$

(12)

$0\le d<1$ is called an isotropic damage parameter, $G$ denotes the initial shear modulus (the shear modulus of an isotropic undamaged material), while $K_M\left(d_T,d_C\right)$ and $K_D\left(d_T,d_C\right)$ are partial bulk moduli dependent on the tensile damage parameter $0\le d_T<1$ and the compressive damage parameter $0\le d_C<1$:

$$\begin{gathered} G=\frac{E}{2\left(1+\nu \right)},\hspace{1em}\hspace{1em}{K}_M\left(d_T,d_C\right)=\frac{E}{6}\left[\frac{1-d_T}{1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)}+\frac{1-d_C}{1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)}\right], \\ K_D\left(d_T,d_C\right)=\frac{E\left(d_C-d_T\right)}{2\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}. \end{gathered}$$

(13)

The tensile damage parameter $d_T$ is linked to the uniaxial tension, while the compressive damage parameter $d_C$ is related to the uniaxial compression. $E>0$ is the initial Young’s modulus of an isotropic undamaged material, and $\nu$ is the respective Poisson’s ratio.

Let us emphasize that the proposed Helmholtz energy is an isotropic and positive homogeneous function of degree 2 with respect to the strain tensor: $H\left(\beta \mathsf{\epsilon }\right)={\beta }^{2}H\left(\mathsf{\epsilon }\right)$ for any $\beta >0$. In general, this property will result in a nonlinear constitutive relationship; in the particular case of the usually used quadratic form of $H$, we obtain a linear relationship. Two first terms in Equation (12) describe bulk energy, while the last one constitutes the distortional part. Potential expressed by Equation (12) satisfies all the requirements described in Section 2, that is, convexity and non-negativity. The introduction of the term $p\sqrt{p^2}$ allows for obtaining a continuous first derivative with respect to $\mathsf{\epsilon }$, which implicates the existence of a continuous relation between the stress and the strain tensors. Denoting:

$$\begin{gathered} K_T\left(d_T\right)=K_M+K_D=\frac{E\left(1-d_T\right)}{3\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]}>0, \\ {K}_C\left(d_C\right)=K_M-K_D=\frac{E\left(1-d_C\right)}{3\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}>0, \end{gathered}$$

(14)

the Helmholtz function’s values can be calculated as:

$$H=\left\{\begin{array}{ll}\left(1-d\right)\left(\frac{3K_T}{2}{p}^2+G\hspace{0.17em}{q}^2\right)\hfill & \mathrm{for}\hspace{1em}p\ge 0,\hfill \\ \left(1-d\right)\left(\frac{3K_C}{2}{p}^2+G\hspace{0.17em}{q}^2\right)\hfill & \mathrm{for}\hspace{1em}p<0.\hfill \end{array}\right.$$

(15)

Equation (15) clearly shows that there is a difference between the stiffnesses for the extending strain states ($p>0$, compatible with tensile stress states, $\xi >0$) and the contracting strain states ($p<0$, compatible with compressive stress states, $\xi <0$), which will be further investigated. For $p=0$, the second derivative of the considered Helmholtz potential with respect to $\mathsf{\epsilon }$ is not determined. $K_T\left(d_T\right)$ and $K_C\left(d_C\right)$ are bulk moduli for tension and compression, respectively. They are always both positive, irrespective of the level of degradation, which is shown in Figure 1a. The difference between $K_T$ and $K_C$ is imprinted on the shape of the contour lines of the free energy potential (see Figure 1b). Initially, for an isotropic undamaged material, the regarded contour is elliptic (yellow line), but due to the growth of $d_T$ or $d_C$ (decrease of the bulk moduli), it can stretch unequally in the $p$ invariant direction (dashed red or green lines), respectively. When $d_T=d_C=0$ and $d\ne 0$, the contour stretches uniformly (blue line) with respect to the undamaged material. If all three parameters increase equally, the contour stretches to a new elliptic curve (pink line).

3.2. Dissipation Function

For concrete, and for quasibrittle materials in general, the form of the dissipation potential is rarely given in the literature. Since most of the elastic damage material models described in the literature are formulated by the damage condition dependent on the stress tensor, see [13,22,34,35], among others. However, the following general clues should be taken into consideration when designing the dissipation function’s form:

• The conjugate damage condition resembles yield conditions, and the damage surface in the Haigh–Westergaard coordinate system has curved meridian, open in the direction of compression and approximately triangular deviatoric cross section [40,41];
• Evolution laws connect damage parameters’ rates to some portions of the strain energy.

Adopting those suggestions, we propose the following dissipation potential:

$$D\left(\dot{d},{\dot{d}}_T,{\dot{d}}_C;\mathsf{\sigma },\mathsf{\epsilon },d,d_T,d_C\right)={\kappa }_D\left(\mathsf{\sigma },\mathsf{\epsilon },d,d_T,d_C\right)\sqrt{\frac{{\dot{d}}_{}^2}{Q_{}}+\frac{{\dot{d}}_T^2}{Q_T}+\frac{{\dot{d}}_C^2}{Q_C}},$$

(16)

where $Q>0$, $Q_T>0$, and $Q_C>0$ are material parameters, while ${\kappa }_D$ is:

$$\begin{gathered} {\kappa }_D\left(\mathsf{\sigma },\mathsf{\epsilon },d,d_T,d_C\right)=\sqrt{Q\hspace{0.17em}{h}_1\left(\mathsf{\epsilon },d_T,d_C\right)+Q_T{h}_2\left(\mathsf{\epsilon },d,d_T\right)+Q_C{h}_3\left(\mathsf{\epsilon },d,d_C\right)-F\left(\mathsf{\sigma }\right)}, \\ {h}_1\left(\mathsf{\epsilon },d_T,d_C\right)=\frac{3K_M\left(d_T,d_C\right)}{2}{p}^2+\frac{3K_D\left(d_T,d_C\right)}{2}\hspace{0.17em}p\sqrt{p^2}+G\hspace{0.17em}{q}^2, \\ {h}_2\left(\mathsf{\epsilon },d,d_T\right)=\frac{3\left(1-d\right)E\hspace{0.17em}p\left(p+\sqrt{p^2}\right)}{4{\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]}^2},\hspace{1em}{h}_3\left(\mathsf{\epsilon },d,d_C\right)=\frac{3\left(1-d\right)E\hspace{0.17em}p\left(p-\sqrt{p^2}\right)}{4{\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}^2}. \end{gathered}$$

(17)

$F\left(\mathsf{\sigma }\right)$ is a convex function such that for any stress tensor, $F\left(\mathsf{\sigma }\right)\le 0$. Those conditions guarantee the required convexity and non-negativity of the considered potential. Dissipation expressed by Equation (16) is homogeneous of the degree 1 function of the rates $\dot{d}$, ${\dot{d}}_T$, and ${\dot{d}}_C$. Interpretation of the material parameters $Q$, $Q_T$, and $Q_C$ will be given in Section 4.

At this point, $F\left(\mathsf{\sigma }\right)$ remains undefined to maintain the flexibility of the model. Function $F\left(\mathsf{\sigma }\right)$ defines the target damage condition via the equation $F\left(\mathsf{\sigma }\right)=0$. It is supposed to properly describe the real material behavior. We will give the details about it after some further derivations in Section 3.6.

The definition of functions in Equation (17) means that $D$ is a function of both rates and nonrates ($\mathsf{\sigma }$, $\mathsf{\epsilon }$, $d$, $d_T$, and $d_C$). The successive differentiations concern only $\dot{d}$, ${\dot{d}}_T$, and ${\dot{d}}_C$; the rest of the arguments are somewhat secondary. In the literature, this distinction is often denoted by a semicolon, that is, $D=D\left(\dot{d},{\dot{d}}_T,{\dot{d}}_C;\mathsf{\sigma },\mathsf{\epsilon },d,d_T,d_C\right)$.

3.3. Cauchy Stress and Relationships between Invariants of Stress and Strain Tensors

Using the first formula in Equation (1) for potential described by Equation (12), the following is derived:

$$\begin{array}{l}\mathsf{\sigma }=\frac{\partial \hspace{0.17em}H}{\partial \hspace{0.17em}\mathsf{\epsilon }}=\left(1-d\right)\left[\sqrt{3}\left(K_{M}p+K_D\sqrt{p^2}\right)\hspace{0.17em}I+2G\hspace{0.17em}e\right]=\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}=\left\{\begin{array}{ll}\left(1-d\right)\left(\sqrt{3}{K}_{T}p\hspace{0.17em}I+G\hspace{0.17em}e\right)\hfill & \mathrm{for}\hspace{1em}p\ge 0,\hfill \\ \left(1-d\right)\left(\sqrt{3}{K}_{C}p\hspace{0.17em}I+G\hspace{0.17em}e\right)\hfill & \mathrm{for}\hspace{1em}p<0.\hfill \end{array}\right.\end{array}$$

(18)

In an undamaged state (null values of damage parameters), Equation (18) reduces to the linear Hooke’s law of an isotropic material:

$$K_C=K_T=\frac{E}{3\left(1-2\nu \right)}\hspace{1em}\Rightarrow \hspace{1em}\mathsf{\sigma }=\frac{E}{3\left(1-2\nu \right)}\left(\mathrm{tr}\mathsf{\epsilon }\right)\hspace{0.17em}I+\frac{E}{1+\nu }e.$$

(19)

It is convenient to use the cylindrical invariants of the stress tensor $\mathsf{\sigma }$, that is:

$$\xi =\frac{1}{\sqrt{3}}\mathrm{tr}\mathsf{\sigma },\hspace{1em}r=\sqrt{\mathrm{tr}{s}^2},\hspace{1em}\mathrm{and}\hspace{1em}\theta =\frac{1}{3}\arccos \frac{\sqrt{6}\hspace{0.17em}\mathrm{tr}{s}^3}{\sqrt{{\mathrm{tr}}^3{s}^2}},$$

(20)

where $s=\mathsf{\sigma }-\frac{1}{3}\left(\mathrm{tr}\mathsf{\sigma }\right)I$ is the stress tensor’s deviator. Splitting stress tensor described by Equation (18) into the hydrostatic and the deviatoric part, we obtain the following constitutive relationships between invariants specified by Equations (11) and (20):

$$\begin{gathered} \xi =3\left(1-d\right)\left(K_{M}p+K_D\hspace{0.17em}\sqrt{p^2}\right)=\left\{\begin{array}{ll}3\left(1-d\right)K_{T}p\hfill & \mathrm{for}\hspace{1em}p\ge 0,\hfill \\ 3\left(1-d\right)K_{C}p\hfill & \mathrm{for}\hspace{1em}p<0,\hfill \end{array}\right. \\ s=2G\hspace{0.17em}\left(1-d\right)e=\hspace{0.17em}\left(1-d\right)\frac{E}{1+\nu }e, r=2G\hspace{0.17em}\left(1-d\right)q=\left(1-d\right)\frac{E}{1+\nu }q. \end{gathered}$$

(21)

Equation (21) can be easily inverted to get the following useful relations:

$$\begin{gathered} p=\frac{1}{3\left(1-d\right)}\left(\frac{K_M}{K_M^2-K_D^2}\xi -\frac{K_D}{K_M^2-K_D^2}\sqrt{{\xi }^2}\right)=\left\{\begin{array}{ll}\frac{\xi }{3\left(1-d\right)K_T}\hfill & \mathrm{for}\hspace{1em}\xi \ge 0,\hfill \\ \frac{\xi }{3\left(1-d\right)K_C}\hfill & \mathrm{for}\hspace{1em}\xi <0,\hfill \end{array}\right. \\ e=\frac{1}{2G\hspace{0.17em}\left(1-d\right)}s=\frac{1+\nu }{E\hspace{0.17em}\left(1-d\right)}s, q=\frac{1}{2G\hspace{0.17em}\left(1-d\right)}r=\frac{1+\nu }{E\hspace{0.17em}\left(1-d\right)}r \end{gathered}$$

(22)

and as a result:

$$\mathsf{\epsilon }=\frac{1}{1-d}\left[\frac{K_M\xi -K_D\sqrt{{\xi }^2}}{3\sqrt{3}\left(K_M^2-K_D^2\right)}\hspace{0.17em}I+\frac{1}{2G}\hspace{0.17em}s\right]=\left\{\begin{array}{ll}\frac{1}{1-d}\left(\frac{\xi }{3\sqrt{3}{K}_T}\hspace{0.17em}I+\frac{1}{2G}\hspace{0.17em}s\right)\hfill & \mathrm{for}\hspace{1em}\xi \ge 0,\hfill \\ \frac{1}{1-d}\left(\frac{\xi }{3\sqrt{3}{K}_C}\hspace{0.17em}I+\frac{1}{2G}s\right)\hfill & \mathrm{for}\hspace{1em}\xi <0.\hfill \end{array}\right.$$

(23)

Derived dual constitutive relationships (Equations (18) and (23)) are homogeneous functions of degree 1 with respect to the strain and stress tensors, respectively. For fixed values of damage parameters, they define nonlinear constitutive relationships of an isotropic elastic material. The jump of stiffness regards the volumetric part of energy. It is clear that passing through $p=0$ ($\xi =0$), the slope of $\xi -p$ curve changes discontinuously for a damaged material according to Equations (21) and (22) (see Figure 2). Graphical representations of the constitutive relations for the selected values of the damage parameters are associated with the contours of the Helmholtz energy given in Figure 1b. Due to an increase in the $d$ parameter, the damaged material responds linearly with reduced bulk stiffness, so the parameter describes isotropic damage. If $d_T\ne 0$ or $d_C\ne 0$, a nonlinear response of material is described. An increase in $d_T$ results in a reduction of the bulk stiffness for tension, while an increase in $d_C$ describes a decrease in bulk stiffness for compression. So far, those responses are not limited by any damage (strength) criterion.

3.4. Generalized Stresses

Differentiating potential specified by Equation (12) according to the second formula in Equation (1), the following relations are obtained:

$$\begin{gathered} {\sigma }_d=-\frac{\partial \hspace{0.17em}H}{\partial \hspace{0.17em}d}=\frac{H}{1-d}=\frac{3}{2}{K}_M{p}^2+\frac{3}{2}{K}_D\hspace{0.17em}p\sqrt{p^2}+G\hspace{0.17em}{q}^2, \\ {\sigma }_{dT}=-\frac{\partial \hspace{0.17em}H}{\partial \hspace{0.17em}{d}_T}=\frac{3\left(1-d\right)E\hspace{0.17em}p\left(p+\sqrt{p^2}\right)}{4{\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]}^2}=\left\{\begin{array}{ll}\frac{3\left(1-d\right)E\hspace{0.17em}}{2{\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]}^2}{p}^2\hfill & \mathrm{for}\hspace{1em}p\ge 0,\hfill \\ 0\hfill & \mathrm{for}\hspace{1em}p<0,\hfill \end{array}\right. \\ {\sigma }_{dC}=-\frac{\partial \hspace{0.17em}H}{\partial \hspace{0.17em}{d}_C}=\frac{3\left(1-d\right)E\hspace{0.17em}p\left(p-\sqrt{p^2}\right)}{4{\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}^2}=\left\{\begin{array}{ll}0\hfill & \mathrm{for}\hspace{1em}p\ge 0,\hfill \\ \frac{3\left(1-d\right)E\hspace{0.17em}}{2{\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}^2}{p}^2\hfill & \mathrm{for}\hspace{1em}p<0.\hfill \end{array}\right. \end{gathered}$$

(24)

Let us notice that, in agreement with indications given in Section 3.2, the generalized stresses have units of energy. ${\sigma }_d$ is the Helmholtz energy divided by the factor $1-d$, so it does not take into account isotropic damage, but it registers the degradation of material connected to the volumetric portion of the energy through moduli $K_M$ and $K_D$, which in turn depend on $d_T$ and $d_C$. The remaining two generalized stresses are associated with the scaled volumetric part $H_K$ of the considered Helmholtz energy in Equation (12), that is:

$$\begin{gathered} \frac{{\sigma }_{dT}}{H_K}=\frac{{\sigma }_{dT}}{\frac{3}{2}\left(1-d\right)\left(K_M{p}^2+K_{D}p\sqrt{p^2}\right)}=\frac{3}{\left(1-d_T\right)\left[1-2\nu +2\hspace{0.17em}{d}_T\left(1+\nu \right)\right]}\hspace{1em}\mathrm{for}\hspace{1em}p>0 \\ ,\frac{{\sigma }_{dC}}{H_K}=\frac{{\sigma }_{dC}}{\frac{3}{2}\left(1-d\right)\left(K_M{p}^2+K_{D}p\sqrt{p^2}\right)}=\frac{3}{\left(1-d_C\right)\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}\hspace{1em}\mathrm{for}\hspace{1em}p<0. \end{gathered}$$

(25)

The generalized stresses may be expressed via the stress tensor’s invariants. Using Equation (22) in Equation (24) results in relationships:

$$\begin{gathered} {\sigma }_d=\left\{\begin{array}{ll}\frac{1}{2{\left(1-d\right)}^{2}E}\left[\frac{1-2\nu +2\left(1+\nu \right)d_T}{1-d_T}{\xi }^2+\left(1+\nu \right)r^2\right]\hfill & \mathrm{for}\hspace{1em}\xi \ge 0,\hfill \\ \frac{1}{2{\left(1-d\right)}^{2}E}\left[\frac{1-2\nu +2\left(1+\nu \right)d_C}{1-d_C}{\xi }^2+\left(1+\nu \right)r^2\right]\hfill & \mathrm{for}\hspace{1em}\xi <0,\hfill \end{array}\right. \\ {\sigma }_{dT}=\left\{\begin{array}{ll}\frac{3\hspace{0.17em}{\xi }^2\hspace{0.17em}}{2\left(1-d\right){\left(1-d_T\right)}^{2}E}\hfill & \mathrm{for}\hspace{1em}\xi \ge 0,\hfill \\ 0\hfill & \mathrm{for}\hspace{1em}\xi <0,\hfill \end{array}\right.{\sigma }_{dC}=\left\{\begin{array}{ll}0\hfill & \mathrm{for}\hspace{1em}\xi \ge 0,\hfill \\ \frac{3\hspace{0.17em}{\xi }^2\hspace{0.17em}}{2\left(1-d\right){\left(1-d_C\right)}^{2}E}\hfill & \mathrm{for}\hspace{1em}\xi <0.\hfill \end{array}\right. \end{gathered}$$

(26)

3.5. Dissipative Stresses, Damage Condition, and Evolution Laws

The dissipative stresses described by Equation (2) are calculated as the following derivatives of the proposed dissipation in Equation (16):

$${\overline{\sigma }}_d=\frac{\partial \hspace{0.17em}D}{\partial \hspace{0.17em}\dot{d}}=\frac{{\kappa }_D^2{\dot{d}}_C^{}}{Q\hspace{0.17em}D},\hspace{1em}{\overline{\sigma }}_{dT}=\frac{\partial \hspace{0.17em}D}{\partial \hspace{0.17em}{\dot{d}}_T}=\frac{{\kappa }_D^2{\dot{d}}_T^{}}{Q_{T}D},\hspace{1em}{\overline{\sigma }}_{dC}=\frac{\partial \hspace{0.17em}D}{\partial \hspace{0.17em}{\dot{d}}_C}=\frac{{\kappa }_D^2{\dot{d}}_C^{}}{Q_{C}D}.$$

(27)

In turn, the damage condition (Equation (6)) becomes:

$$\overline{F}\left({\overline{\sigma }}_d,{\overline{\sigma }}_{dT},{\overline{\sigma }}_{dC}\right)=Q\hspace{0.17em}{\overline{\sigma }}_d^2+Q_T{\overline{\sigma }}_{dT}^2+Q_C{\overline{\sigma }}_{dC}^2-{\kappa }_D^2=0,$$

(28)

and the evolution laws (Equation (8)) are:

$$\dot{d}=\lambda \frac{\partial \overline{F}}{{\overline{\sigma }}_d}=2\lambda Q\hspace{0.17em}{\overline{\sigma }}_d^{},\hspace{1em}{\dot{d}}_C=\lambda \frac{\partial \overline{F}}{{\overline{\sigma }}_{dC}}=2\lambda Q_C{\overline{\sigma }}_{dC}^{},\hspace{1em}{\dot{d}}_T=\lambda \frac{\partial ;\overline{F}}{{\overline{\sigma }}_{dT}}=2\lambda Q_T{\overline{\sigma }}_{dT}^{}.$$

(29)

Equations (28) and (29) are complemented by the consistency condition (Equation (7)).

3.6. Orthogonality Principle, Damage Condition, and Evolution Laws Expressed via Cauchy Stress

In the considered case, the orthogonality principle (Equation (9)) establishes three equalities:

$${\overline{\sigma }}_d={\sigma }_d,\hspace{1em}{\overline{\sigma }}_{dT}={\sigma }_{dT},\hspace{1em}\mathrm{and}\hspace{1em}{\overline{\sigma }}_{dC}={\sigma }_{dC},$$

(30)

which bind spaces of the generalized and the dissipative stresses. Now, using the definition of ${\kappa }_D^{}$ according to Equation (17) and of the generalized stresses (Equation (26)), Equation (28) can be written as:

$$F\left(\mathsf{\sigma }\right)=0$$

(31)

and the evolution laws (Equation (29)) are in the form:

$$\dot{d}=2\lambda Q\hspace{0.17em}{\sigma }_d^{},\hspace{1em}{\dot{d}}_C=2\lambda Q_C{\sigma }_{dC}^{},\hspace{1em}{\dot{d}}_T=2\lambda Q_T{\sigma }_{dT}^{}.$$

(32)

Consistency condition expressed via Cauchy stress tensor reads:

$$\dot{F}\left(\mathsf{\sigma }\right)=0.$$

(33)

In order to conduct further derivations, we need to specify the function $F\left(\mathsf{\sigma }\right)$, which properly captures a real material behavior. For concrete and rocks, we may assume that:

$$\begin{gathered} F\left(\xi ,r,\theta \right)=3\hspace{0.17em}C\left(\xi \right)r^2-2\sqrt{\left(1-{\gamma }^2\right){\left(C{\left(\xi \right)}^2-r^2\right)}^3}-\hspace{0.17em}\left(2-{\gamma }^2\right)C{\left(\xi \right)}^3-2\hspace{0.17em}\gamma \hspace{0.17em}{r}^3\cos 3\theta =0, \\ \mathrm{with}\hspace{1em}C\left(\xi \right)=\sqrt{{\left(\frac{c-\xi }{B}\right)}^2-b^2}, \end{gathered}$$

(34)

where $c>0$, $b>0$, $B>0$, and $-1<\gamma <1$ are material parameters. Function $F$ is convex for any $\xi$, $r$, and $\theta$. The damage surface $F=0$ possesses the desired features mentioned in Section 3.2; see the cross sections of the damage surface in Figure 3, where the shapes of tensile ($\theta =0$) and compressive ($\theta =\pi /3$) meridians are compared with experimental data, and possible shapes of the deviatoric cross sections are shown.

Finding the four parameters $c$, $b$, $B$, and $\gamma$ involved in the damage condition requires knowing the damage limits (strengths) for four independent tests. For concrete, the basic experiment is the uniaxial compression. Additionally, the uniaxial tension, biaxial equilateral compression, and triaxial compression tests are representative and well grounded in the literature. Let ${\sigma }_T$, ${\sigma }_C$, ${\sigma }_{BC}$, and $\eta \hspace{0.17em}{\sigma }_{TC}$ with $\eta >1$ be the damage limits for the uniaxial tension, uniaxial compression, biaxial equilateral compression, and triaxial compression test with $\mathsf{\sigma }\to \mathrm{diag}[-\eta {\sigma }_{TC},-{\sigma }_{TC},-{\sigma }_{TC}]$ (the stress tensor is represented by the shown diagonal matrix), respectively. The tests are located on either the tensile meridian ($\theta =0$) or the compressive meridian ($\theta =\pi /3$), and for those cases, damage condition expressed by Equation (34) reduces to quite simple expressions (for particulars, see [40]). This useful feature allows for solving the set of four equations—$F\left({\xi }_T,r_T,0\right)=0$, $F\left({\xi }_{BC},r_{BC},0\right)=0$, $F\left({\xi }_C,r_C,\pi /3\right)=0$, and $F\left({\xi }_{TC},r_{TC},\pi /3\right)=0$—for the unknown parameters $c$, $b$, $B$, and $\gamma$. The values of the cylindrical invariants, defined by Equation (20), for the specified cases are:

$$\begin{gathered} {\xi }_T=\frac{{\sigma }_T}{\sqrt{3}},\hspace{1em}{\xi }_C=-\frac{{\sigma }_C}{\sqrt{3}},\hspace{1em}{\xi }_{BC}=-\frac{2{\sigma }_{BC}}{\sqrt{3}},\hspace{1em}{\xi }_C=-\frac{\left(\eta +2\right){\sigma }_C}{\sqrt{3}}, \\ {r}_T=\sqrt{\frac{2}{3}}\hspace{0.17em}{\sigma }_T,\hspace{1em}{r}_C=\sqrt{\frac{2}{3}}\hspace{0.17em}{\sigma }_C,\hspace{1em}{r}_{BC}=\sqrt{\frac{2}{3}}\hspace{0.17em}{\sigma }_{BC},\hspace{1em}{r}_{TC}=\sqrt{\frac{2}{3}}\hspace{0.17em}\left(\eta +1\right){\sigma }_{TC}. \end{gathered}$$

(35)

Introducing notations:

$$\tilde{a}=\frac{{\xi }_{BC}{r}_T^2-{\xi }_T{r}_{BC}^2}{r_{BC}^2-r_T^2}, \tilde{b}=\frac{{\xi }_{TC}{r}_C^2-{\xi }_C{r}_{TC}^2}{r_{TC}^2-r_C^2}, \tilde{c}=\frac{{\xi }_{BC}^2{r}_T^2-{\xi }_T^2{r}_{BC}^2}{r_{BC}^2-r_T^2}, \tilde{d}=\frac{{\xi }_{TC}^2{r}_C^2-{\xi }_C^2{r}_{TC}^2}{r_{TC}^2-r_C^2}$$

(36)

the following formulae for parameters involved in the damage condition are derived:

$$\begin{gathered} c=\frac{\tilde{c}-\tilde{d}}{2(\tilde{a}-\tilde{b})}, t=\frac{r_C}{r_T}\sqrt{\frac{\left(\tilde{b}-\tilde{a}\right){\xi }_T^2+\left(\tilde{c}-\tilde{d}\right){\xi }_T-\tilde{a} \tilde{d}+\tilde{b} \tilde{c}}{\left(\tilde{b}-\tilde{a}\right){\xi }_C^2+\left(\tilde{c}-\tilde{d}\right){\xi }_C-\tilde{a} \tilde{d}+\tilde{b} \tilde{c}}}, \gamma =\frac{3\sqrt{3}\hspace{0.33em}t\hspace{0.17em}(1-t)}{2\sqrt{{\left(t^2-t+1\right)}^3}}, \\ b=\frac{r_T}{2\cos \left(\frac{1}{3}\arccos \gamma \right)}\sqrt{\frac{{\left(\tilde{c}-\tilde{d}\right)}^2+4\left(\tilde{a} \tilde{d}-\tilde{b} \tilde{c}\right)\left(\tilde{b}-\tilde{a}\right)}{\left(\tilde{b}-\tilde{a}\right)\left[\left(\tilde{b}-\tilde{a}\right){\xi }_T^2+\left(\tilde{c}-\tilde{d}\right){\xi }_T-\tilde{a} \tilde{d}+\tilde{b} \tilde{c}\right]}}, \\ B=\frac{1}{b}\sqrt{\frac{1}{4}{\left(\frac{\tilde{c}-\tilde{d}}{\tilde{b}-\tilde{a}}\right)}^2+\frac{\tilde{a} \tilde{d}-\tilde{b} \tilde{c}}{\tilde{b}-\tilde{a}}}. \end{gathered}$$

(37)

A detailed description of the function’s $F$ properties with suitable graphs and particulars of the calibration procedure can be found in [40].

Obtaining the closed formulae for material parameters is useful, especially with regard to hardening. Hardening can be taken into account if we make the limits (strengths) ${\sigma }_T$, ${\sigma }_C$, ${\sigma }_{BC}$, and $\eta \hspace{0.17em}{\sigma }_{TC}$ dependent on the hardening variables, which can be the damage parameters, meaning that the limits, instead of being constant, are assumed to be the functions ${\sigma }_T={\sigma }_T\left(d,d_T,d_C\right)$, and so forth. Then Equation (37) allow for directly finding $c$, $b$, $B$, and $\gamma$ as functions of the damage parameters. In this way, the damage function $F$ becomes dependent on $\sigma$, $d$, $d_T$, and $d_C$. The subsequent derivations will be based on the assumption that hardening is dependent on those parameters. If any of the parameters $c$, $b$, $B$, and $\gamma$ were sought by curve fitting, the procedure would be quite complicated. The existence of the closed calibration formulae simplifies the issue. Still, the damage condition is fairly complex and quite general.

3.7. Tensor of Tangent Stiffness of Elastic Damaged Material

Let us calculate the increment of $\mathsf{\sigma }\left(\mathsf{\epsilon },d,d_T,d_C\right)$ based on Equation (18):

$$\begin{array}{c}\dot{\mathsf{\sigma }}=\frac{\partial \hspace{0.17em}\mathsf{\sigma }}{\partial \hspace{0.17em}\mathsf{\epsilon }}\mathbf{·}\dot{\mathsf{\epsilon }}+\frac{\partial \hspace{0.17em}\mathsf{\sigma }}{\partial \hspace{0.17em}d}\dot{d}+\frac{\partial \hspace{0.17em}\mathsf{\sigma }}{\partial \hspace{0.17em}{d}_T}{\dot{d}}_T+\frac{\partial \hspace{0.17em}\mathsf{\sigma }}{\partial \hspace{0.17em}{d}_C}{\dot{d}}_C=\\ =\left(1-d\right)\left[\left(K_M+K_D\frac{p}{\sqrt{p^2}}-\frac{2}{3}G\right)\hspace{0.17em}I\otimes I+2G\hspace{0.17em}1\right]\mathbf{·}\dot{\mathsf{\epsilon }}+\\ -\dot{d}\left[\sqrt{3}\left(K_{M}p+K_D\sqrt{p^2}\right)\hspace{0.17em}I+2G\hspace{0.17em}e\right]+\sqrt{3}\left(1-d\right)\left\{{\dot{d}}_T{g}_T\left(p+\sqrt{p^2}\right)+{\dot{d}}_C{g}_C\left(p-\sqrt{p^2}\right)\right\}I,\end{array}$$

(38)

where $\mathbf{·}$ denotes the full contraction of tensors, $\otimes$ denotes the tensor product, $1\hspace{0.17em}=\frac{1}{2}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)b_i\otimes b_j\otimes b_k\otimes b_l$ is for $i,j,k,l=1,2,3$ using the summation convention, $\left\{b_1,b_2,b_3\right\}$ is an orthonormal basis, and ${\delta }_{ik}$ is the Kronecker delta, and:

$$\begin{gathered} g_T\left(d_T\right)=\frac{\partial K_M}{\partial d_T}=\frac{\partial K_D}{\partial d_T}=-\frac{E}{2{(1-2\nu +2d_T(1+\nu ))}^2}, \\ {g}_C\left(d_C\right)=\frac{\partial K_M}{\partial d_C}=-\frac{\partial K_D}{\partial d_C}=-\frac{E}{2{(1-2\nu +2d_C(1+\nu ))}^2}. \end{gathered}$$

(39)

Now, we use the consistency condition described by Equation (33):

$$\dot{F}=\left(\frac{\partial F}{\partial \xi }\frac{\partial \xi }{\partial \mathsf{\sigma }}+\frac{\partial F}{\partial r}\frac{\partial r}{\partial \mathsf{\sigma }}+\frac{\partial F}{\partial \cos 3\theta }\frac{\partial \cos 3\theta }{\partial \mathsf{\sigma }}\right)\mathbf{·}\dot{\sigma }+\frac{\partial F}{\partial d}\dot{d}+\frac{\partial F}{\partial d_T}{\dot{d}}_T+\frac{\partial F}{\partial d_C}{\dot{d}}_C=0.$$

(40)

Three last components appear as a result of hardening (or softening) dependent on the current values of damage parameters. Introducing:

$$\begin{gathered} f_1=\frac{\partial F}{\partial C}=3\left[r^2-C^2(2-{\gamma }^2)-2C\sqrt{(1-{\gamma }^2)(C^2-r^2)}\right], \\ {C}_{\xi }=\frac{\partial C}{\partial \xi }=-\frac{c-\xi }{B\sqrt{{\left(c-\xi \right)}^2-b^2{B}^2}}, \\ {f}_2=\frac{\partial F}{\partial r}=6\hspace{0.17em}r\left(C+\sqrt{\left(1-{\gamma }^2\right)\left(C^2-r^2\right)}-\gamma \hspace{0.17em}r\hspace{0.17em}\cos 3\theta \right), f_3=\frac{\partial F}{\partial \cos 3\theta }=-2\hspace{0.17em}\gamma \hspace{0.17em}{r}^3, \\ {g}_1=\frac{\partial \xi }{\partial \mathsf{\sigma }}=\frac{1}{\sqrt{3}}{\rm I}, g_2=\frac{\partial r}{\partial \mathsf{\sigma }}=\frac{1}{r}s, g_3=\frac{\partial \cos 3\theta }{\partial \mathsf{\sigma }}=\frac{3\sqrt{6}}{r^3}\left(s^2-\frac{r\hspace{0.17em}\cos 3\theta }{\sqrt{6}}s-\frac{r^2}{3}I\right), \\ {K}_0=K_M+K_D\frac{p}{\sqrt{p^2}}\hspace{1em}\mathrm{and}\hspace{1em}{S}_0=Q_T{\sigma }_{dT}^{}{g}_T\left(p+\sqrt{p^2}\right)+Q_C{\sigma }_{dC}^{}{g}_C\left(p-\sqrt{p^2}\right) \end{gathered}$$

(41)

and taking into account the evolution laws (Equation (32)), the consistency condition (Equation (40)) and the stress increment (Equation (38)) become:

$$\dot{F}=\left(f_1{C}_{\xi }{g}_1+f_2\hspace{0.17em}{g}_2+f_3\hspace{0.17em}{g}_2\right)\mathbf{·}\dot{\mathsf{\sigma }}+2\lambda \left(\frac{\partial F}{\partial d}Q{\sigma }_d^{}+\frac{\partial F}{\partial d_T}{Q}_T{\sigma }_{dT}^{}+\frac{\partial F}{\partial d_C}{Q}_C{\sigma }_{dC}^{}\right)=0,$$

(42)

$$\dot{\mathsf{\sigma }}=\left(1-d\right)\left[\left(K_0-\frac{2}{3}G\right)I\otimes I+2G\hspace{0.17em}1\right]\mathbf{·}\dot{\mathsf{\epsilon }}-2\lambda \left[Q{\sigma }_d^{}\left(\sqrt{3}{K}_{0}p\hspace{0.17em}I+2G\hspace{0.17em}e\right)+\sqrt{3}\left(1-d\right)S_{0}I\right].$$

(43)

Using the above formula for the stress increment, we solve Equation (42) to obtain the damage multiplier:

$$\begin{gathered} \lambda =\frac{1}{2R_0}\left(1-d\right)\left[\sqrt{3}{f}_1{C}_{\xi }{K}_{0}I+2G\left(f_2{g}_2+f_3{g}_3\right)\right]\mathbf{·}\dot{\mathsf{\epsilon }},\hspace{1em}\mathrm{where} \\ \begin{array}{c}{R}_0=Q{\sigma }_d^{}\left(3f_1{C}_{\xi }{K}_{0}p+\hspace{0.17em}\frac{f_{2}r}{\hspace{0.17em}1-d}\right)-\sqrt{3}\left(1-d\right)f_1{C}_{\xi }{S}_0+\\ -\left(Q{\sigma }_d^{}\frac{\partial F}{\partial d}+Q_T{\sigma }_{dT}^{}\frac{\partial F}{\partial d_T}\hspace{0.33em}+Q_C{\sigma }_{dC}^{}\frac{\partial F}{\partial d_C}\right).\end{array} \end{gathered}$$

(44)

The final step in the derivation of the incremental form of the constitutive relationship is to apply Equation (44) in Equation (43), which results in:

$$\dot{\mathsf{\sigma }}=C^{\tan }·\dot{\mathsf{\epsilon }},$$

(45)

where $C^{\tan }$ is the sought tensor of the tangent stiffness of the considered elastic damaged material:

$$C^{\tan }\hspace{0.33em}=\left(1-d\right)\left(C_0-C_D\right)$$

(46)

with:

$$\begin{gathered} C_0=\left(K_0-\frac{2}{3}G\right)\hspace{0.17em}I\otimes I+2G\hspace{0.17em}1, C_D=\frac{1}{R_0}{c}_1\otimes c_2, \\ {c}_1=\left(\sqrt{3}{Q}_{}{\sigma }_d{K}_{0}p-\left(1-d\right)S_0\right)I+2G\hspace{0.17em}Q{\sigma }_d\hspace{0.17em}e, \\ {c}_2=\sqrt{3}{f}_1{C}_{\xi }{K}_{0}I+2G\left(f_2{g}_2+f_3{g}_3\right). \end{gathered}$$

(47)

Let us emphasize that all the proper symmetries, that is, $C_{ijkl}^{\tan }=\hspace{0.17em}{C}_{klij}^{\tan }=C_{jikl}^{\tan }=C_{ijlk}^{\tan }$, occur. The tensor is determined for $p\in R\backslash \left\{0\right\}$ ($\xi \in R\backslash \left\{0\right\}$).

We can also derive the fourth-order tensor of the secant stiffness $C^{\sec }$:

$$\mathsf{\sigma }=C^{\sec }\mathbf{·}\mathsf{\epsilon },$$

(48)

merely by rearranging Equation (18). The following is obtained:

$$C^{\sec }=\left(1-d\right)\left[\left(K_0-\frac{2}{3}G\hspace{0.17em}\right)I\otimes I+2G\hspace{0.17em}1\right]=\left(1-d\right)C_0.$$

(49)

4. Model’s Predictions for Uniaxial Stress States and Pure Shear

Below we present the numerical results for the uniaxial tension, compression, and cyclic tension–compression and pure shear omitting hardening. The aim of the calculations is to show the basic features of the model, to assess the quality of results rather than their quantity. We seek for physical (experimental) interpretation of the damage parameters $d$, $d_T$, and $d_C$ and investigate the influence of the parameters $Q$, $Q_T$, and $Q_C$ on the model predictions. Including hardening is quite simple, but it blurs the equations, so for the sake of clarity, it is left out in most cases.

For the chosen tests, the equations are simple enough to use basic solvers available in Wolfram Mathematica; there is no necessity for creating a special finite element code, although for more complex problems, it would be inevitable.

4.1. Uniaxial Tension

For the uniaxial tension, the strain and the stress states are:

$$\mathsf{\epsilon }={\mathsf{\epsilon }}_{11}{b}_1\otimes b_1+{\mathsf{\epsilon }}_{22}{b}_2\otimes b_2+{\mathsf{\epsilon }}_{33}{b}_3\otimes b_3,\hspace{1em}\hspace{1em}\mathsf{\sigma }={\sigma }_{11}\hspace{0.17em}{b}_1\otimes b_1,\hspace{1em}\mathrm{and}\hspace{1em}{\sigma }_{11}\ge 0.$$

(50)

As invariants defined by Equation (20) for this case are:

$$\xi =\frac{1}{\sqrt{3}}{\sigma }_{11},\hspace{1em}r=\sqrt{\frac{2}{3}}\hspace{0.17em}{\sigma }_{11},\hspace{1em}\theta =0,$$

(51)

the invariants of $\mathsf{\epsilon }$ can be found using constitutive relations for the volumetric and distortional parts according to Equation (22), become:

$$p=\frac{1-2\nu +2d_T\left(1+\nu \right)}{\sqrt{3}\left(1-d\right)\left(1-d_T\right)E}{\sigma }_{11} \mathrm{and} q=\sqrt{\frac{2}{3}}\frac{1+\nu }{\hspace{0.17em}E\hspace{0.17em}\left(1-d\right)}{\sigma }_{11},$$

(52)

so the strain tensor is:

$$\mathsf{\epsilon }=\frac{1}{\left(1-d\right)E}\left[\frac{1-2\nu +2d_T\left(1+\nu \right)}{3\left(1-d_T\right)}{\sigma }_{11}\hspace{0.17em}I+\left(1+\nu \right)\hspace{0.17em}s\right]\to \left[\begin{array}{ccc}{\mathsf{\epsilon }}_{11}& 0& 0\\ 0& {\mathsf{\epsilon }}_{22}& 0\\ 0& 0& {\mathsf{\epsilon }}_{33}\end{array}\right].$$

(53)

The above constitutive equation is equivalent to the following relation between the stress and strain components:

$${\mathsf{\epsilon }}_{11}=\frac{{\sigma }_{11}}{\left(1-d\right)\left(1-d_T\right)E},\hspace{1em}{\mathsf{\epsilon }}_{22}={\mathsf{\epsilon }}_{33}=-\left[\nu -d_T\left(1+\nu \right)\right]{\mathsf{\epsilon }}_{11},$$

(54)

which, for an undamaged state of material, reduce to the well-known Hooke’s law:

$${\mathsf{\epsilon }}_{11}=\frac{{\sigma }_{11}}{E},\hspace{1em}{\mathsf{\epsilon }}_{22}={\mathsf{\epsilon }}_{33}=-\nu \hspace{0.17em}{\mathsf{\epsilon }}_{11}.$$

(55)

Let ${\mathsf{\epsilon }}_{11}$ be dependent on time $t$ and prescribed by the formula:

$${\mathsf{\epsilon }}_{11}\left(t\right)=a\hspace{0.17em}t,\hspace{1em}a>0.$$

(56)

Then using Equation (54), we obtain:

$${\sigma }_{11}\left(t\right)=\left[1-d\left(t\right)\right]\left[1-d_T\left(t\right)\right]E\hspace{0.17em}a\hspace{0.17em}t,\hspace{1em}{\mathsf{\epsilon }}_{22}\left(t\right)={\mathsf{\epsilon }}_{33}\left(t\right)=-\left[\nu -d_T\left(t\right)\left(1+\nu \right)\right]a\hspace{0.17em}t.$$

(57)

The loading process starting from a natural state, that is, $\mathsf{\sigma }=0$, $\mathsf{\epsilon }=0$, $d=0$, $d_T=0$, and $d_C=0$, is purely elastic until the damage limit for the uniaxial tension ${\sigma }_T$ is attained. The following equations remain true:

$${\sigma }_{11}\left(t\right)=Eat,\hspace{1em}{\mathsf{\epsilon }}_{22}\left(t\right)={\mathsf{\epsilon }}_{33}\left(t\right)=-\nu \hspace{0.17em}a\hspace{0.17em}t,\hspace{1em}d\left(t\right)=0,\hspace{1em}{d}_T\left(t\right)=0,\hspace{1em}{d}_C\left(t\right)=0$$

(58)

for $t<t_0={\sigma }_T/\left(E\hspace{0.17em}a\right)$. Further increase in strain according to the program prescribed by Equation (56) leads to a development of damage. The damage condition defined by Equation (34) reduces to:

$${\sigma }_{11}-{\sigma }_T=0.$$

(59)

The generalized stresses defined by Equation (24) become:

$${\sigma }_d^{}=\frac{E\left(1-d_T\right)}{2}{\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\sigma }_{dT}^{}=\frac{E\left(1-d\right)}{2}{\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\sigma }_{dC}^{}=0.$$

(60)

In turn, the evolution laws (Equation (32)) take the form:

$$\dot{d}=\lambda QE\left(1-d_T\right){\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\dot{d}}_T=\lambda Q_{T}E\left(1-d\right){\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\dot{d}}_C=0.$$

(61)

As a result, the compressive damage parameter stays null throughout the process, that is, $d_C\left(t\right)=0$, and for $t\ge t_0$, the other quantities must meet the following conditions:

$$\left\{\begin{array}{c}\left(1-d\left(t\right)\right)\left(1-d_T\left(t\right)\right)=\frac{t_0}{t},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \dot{d}\left(t\right)\left(1-d\left(t\right)\right)Q_T={\dot{d}}_T\left(t\right)\left(1-d_T\left(t\right)\right)Q,\\ d_T\left(t_0\right)=0,\hspace{1em}\hspace{1em}d\left(t_0\right)=0.\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(62)

The above differential equations are solved numerically using Wolfram Mathematica for the following values of material constants: $E=20\hspace{0.17em}\mathrm{GPa}$, $\nu =0.2$, ${\sigma }_T=2\hspace{0.17em}\mathrm{MPa}$, and $a=0.0001\hspace{0.17em}{\mathrm{s}}^{-1}$, $Q=1$, $Q_T=2$. For these data, the damage limit is reached at $t_0=1\hspace{0.17em}\mathrm{s}$.

Initially, the material is elastic with stiffness equal to the initial Young’s modulus $E$. After reaching the damage limit, the stress becomes constant due to the assumed lack of hardening, as presented in Figure 4a. The secant stiffness becomes $\left(1-d\right)\left(1-d_T\right)E$. The degradation of the initial stiffness (see Figure 5a) mirrors the decrease in the current slope of a hypothetic unloading curve, ${\sigma }_{11}-{\mathsf{\epsilon }}_{11}$. As shown in Figure 4b, the lateral strain components decrease to the limit point, but after the degradation process origins, they start to increase. As a result, the trace of the strain tensor, that is, the relative change of volume, increases approximately piecewise linearly (for ${\sigma }_{11}={\sigma }_T$, it is nonlinear). A more dynamic growth of the volumetric strain is observed as the degradation enters. As $d\to 1$ or $d_T\to 1$ (Figure 5b), the ability of material to carry any loading vanishes, and total failure occurs. It is reasonable to establish arbitrary limit values of the damage parameters $d_0<1$ and $d_{T0}<1$, marking a partial but severe-enough level of degradation.

Knowing the experimental curves of the axial and lateral strains, we can easily find a function of the tensile damage parameter based on Equation (54):

$$d_T=-\frac{1}{1+\nu }\left(\frac{{\mathsf{\epsilon }}_{22}}{{\mathsf{\epsilon }}_{11}}+\nu \right).$$

(63)

Thus, the progress of $d_T$ can be found directly in experimental results. In the uniaxial tension, factor $\left(1+\nu \right)d_T$ lowers the initial Poisson’s ratio $\nu$.

Figure 6 depicts the evolution of damage parameters for various ratios, $Q/Q_T$. We set $Q=1$ and control the model prediction by changing $Q_T$. It is apparent that by increasing $Q/Q_T$, the isotropic damage parameter prevails. On the contrary, a decrease in the quotient lets the tensile damage parameter be the leading function in the material degradation process. For $Q/Q_T<1$, the degradation due to the change of volume is bigger than the isotropic degradation ($d/d_T<1$), and vice versa for $Q/Q_T>1$, we have $d/d_T>1$.

4.2. Uniaxial Compression

The strain and stress tensors for the uniaxial compression are the same as for the uniaxial tension excluding signs of the components, that is:

$$\mathsf{\epsilon }={\mathsf{\epsilon }}_{11}{b}_1\otimes b_1+{\mathsf{\epsilon }}_{22}{b}_2\otimes b_2+{\mathsf{\epsilon }}_{33}{b}_3\otimes b_3,\hspace{1em}\hspace{1em}\mathsf{\sigma }={\sigma }_{11}\hspace{0.17em}{b}_1\otimes b_1,\hspace{1em}\mathrm{and}\hspace{1em}{\sigma }_{11}\le 0.$$

(64)

As a result, the cylindrical invariants defined by Equation (20) become:

$$\xi =\frac{1}{\sqrt{3}}{\sigma }_{11},\hspace{1em}r=-\sqrt{\frac{2}{3}}\hspace{0.17em}{\sigma }_{11},\hspace{1em}\theta =\pi /3.$$

(65)

Next, using the constitutive relationship between invariants (Equation (22)), we obtain:

$$p=\frac{1-2\nu +2d_C\left(1+\nu \right)}{\sqrt{3}\left(1-d\right)\left(1-d_C\right)E}{\sigma }_{11},\hspace{1em}q=-\frac{1}{\sqrt{6}\hspace{0.17em}G\hspace{0.17em}\left(1-d\right)}{\sigma }_{11}$$

(66)

and:

$$\mathsf{\epsilon }=\frac{1}{\left(1-d\right)E}\left[\frac{1-2\nu +2d_C\left(1+\nu \right)}{3\left(1-d_C\right)}{\sigma }_{11}\hspace{0.17em}I+\left(1+\nu \right)\hspace{0.17em}s\right]\to \left[\begin{array}{ccc}{\mathsf{\epsilon }}_{11}& 0& 0\\ 0& {\mathsf{\epsilon }}_{22}& 0\\ 0& 0& {\mathsf{\epsilon }}_{33}\end{array}\right].$$

(67)

The above constitutive equation is equivalent to the following relationships between the stress and strain components:

$${\mathsf{\epsilon }}_{11}=\frac{{\sigma }_{11}}{\left(1-d\right)\left(1-d_C\right)E},\hspace{1em}{\mathsf{\epsilon }}_{22}={\mathsf{\epsilon }}_{33}=-\left[\nu -d_C\left(1+\nu \right)\right]{\mathsf{\epsilon }}_{11}.$$

(68)

As in Section 4.1, Equation (68) for the undamaged material reduce to Hooke’s law (Equation (55)).

Again, we assume that ${\mathsf{\epsilon }}_{11}$ is dependent on time as follows:

$${\mathsf{\epsilon }}_{11}\left(t\right)=-a\hspace{0.17em}t,\hspace{1em}a>0.$$

(69)

Then using Equation (68), we obtain:

$${\sigma }_{11}\left(t\right)=-\left[1-d\left(t\right)\right]\left[1-d_C\left(t\right)\right]E\hspace{0.17em}a\hspace{0.17em}t,\hspace{1em}{\mathsf{\epsilon }}_{22}\left(t\right)={\mathsf{\epsilon }}_{33}\left(t\right)=\left[\nu -d_C\left(t\right)\left(1+\nu \right)\right]a\hspace{0.17em}t.$$

(70)

For the elastic range, that is, as long as the axial stress does not reach the damage limit for the uniaxial compression ${\sigma }_C$ or $t<t_0={\sigma }_C/\left(E\hspace{0.17em}a\right)$, the governing equations are:

$${\sigma }_{11}\left(t\right)=-Ea\hspace{0.17em}t,\hspace{1em}{\mathsf{\epsilon }}_{22}\left(t\right)={\mathsf{\epsilon }}_{33}\left(t\right)=\nu \hspace{0.17em}a\hspace{0.17em}t,\hspace{1em}d\left(t\right)=0,\hspace{1em}{d}_T\left(t\right)=0,\hspace{1em}{d}_C\left(t\right)=0.$$

(71)

The subsequent loading leads to deterioration of material. The damage condition expressed by Equation (34) takes the form:

$${\sigma }_{11}+{\sigma }_C=0.$$

(72)

The generalized stresses (Equation (24)) are:

$${\sigma }_d^{}=\frac{E\left(1-d_C\right)}{2}{\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\sigma }_{dT}^{}=0,\hspace{1em}{\sigma }_{dC}^{}=\frac{E\left(1-d\right)}{2}{\mathsf{\epsilon }}_{11}^2,$$

(73)

so the evolution laws (Equation (32)) reduce to:

$$\dot{d}=\lambda QE\left(1-d_C\right){\mathsf{\epsilon }}_{11}^2,\hspace{1em}{\dot{d}}_T=0,\hspace{1em}{\dot{d}}_C=\lambda Q_{C}E\left(1-d\right){\mathsf{\epsilon }}_{11}^2.$$

(74)

Contrary to the uniaxial tension, $d_T\left(t\right)=0$, while $d_C$ increases. For $t\ge t_0$, the following system governs the process:

$$\left\{\begin{array}{c}\left(1-d\left(t\right)\right)\left(1-d_C\left(t\right)\right)=\frac{t_0}{t},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \dot{d}\left(t\right)\left(1-d\left(t\right)\right)Q_C={\dot{d}}_C\left(t\right)\left(1-d_C\left(t\right)\right)Q,\\ d_C\left(t_0\right)=0,\hspace{1em}\hspace{1em}d\left(t_0\right)=0.\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(75)

The above set of equations is solved numerically using the NDSolve procedure of Wolfram Mathematica for $E=20\hspace{0.17em}\mathrm{GPa}$, $\nu =0.2$, ${\sigma }_C=20\hspace{0.17em}\mathrm{MPa}$, and $a=0.0001\hspace{0.17em}{\mathrm{s}}^{-1}$, $Q=1$ and $Q_C=5$. The damage limit is attained at $t_0=10\hspace{0.17em}\mathrm{s}$.

As shown in Figure 7 and Figure 8, the characteristics of the obtained function are similar to that for the uniaxial tension. When the degradation process is active, the secant stiffness is lower compared with the initial value and is scaled by factor $\left(1-d\right)\left(1-d_T\right)$. As damage develops factor $\left(1-d\right)\left(1-d_C\right)$ tends to zero, which leads to complete failure of material.

Analogous to Equation (63), the compressive damage parameter can be determined by knowing the ratio of the strain tensor components (the generalized Poisson’s ratio function) and the initial Poisson’s ratio according to:

$$d_C=-\frac{1}{1+\nu }\left(\frac{{\mathsf{\epsilon }}_{22}}{{\mathsf{\epsilon }}_{11}}+\nu \right).$$

(76)

This time factor $\left(1+\nu \right)d_C$ is subtracted from $\nu$ to obtain the current ratio of strains.

The ratio $Q/Q_C$ governs the proportions of the isotropic and compressive damage. Based on Figure 9, if $Q/Q_C>1$, then throughout the process, $d/d_C>1$, and for $Q/Q_T<1$, $d/d_C<1$. In the case of $Q=Q_C$, a balanced $d\left(t\right)=d_C\left(t\right)$ exists as a result of the symmetry of system defined by Equation (75) with reference to the sought functions.

During the loading process, the relative change of volume (the trace of the strain tensor) decreases, initially slowly, then rapidly. As shown above, for both the uniaxial tension and the uniaxial compression, the generalized Poisson’s ratio $-{\mathsf{\epsilon }}_{22}/{\mathsf{\epsilon }}_{11}$ can be either negative or positive (Figure 10a), but $\mathrm{tr}\mathsf{\epsilon }/{\mathsf{\epsilon }}_{11}$ is always non-negative (Figure 10b). As a result, the phenomenon of dilation, present in concrete and rocks, cannot be described properly only by this model.

As mentioned in Section 3, the generalized stresses are some parts of the strain energy, so they can be interpreted as proper areas (triangles $QA\prime A″$) connected to $\mathsf{\sigma }-\mathsf{\epsilon }$ curves, which is depicted in Figure 11. They can be associated to the strain energies of material in the undamaged state.

The actual behavior of material can include hardening/softening. For such a case, a preliminary model prediction is compared with experimental data for concrete [25,42,43]. Equations (65)–(74) remain true, but the system described by Equation (75) changes as follows:

$$\left\{\begin{array}{c}\left(1-d\left(t\right)\right)\left(1-d_C\left(t\right)\right)=\frac{t_0}{t}\frac{{\sigma }_C\left(t\right)}{{\sigma }_{C0}},\hspace{1em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\\ \dot{d}\left(t\right)\left(1-d\left(t\right)\right)Q_C={\dot{d}}_C\left(t\right)\left(1-d_C\left(t\right)\right)Q,\\ d_C\left(t_0\right)=0,\hspace{1em}\hspace{1em}d\left(t_0\right)=0.\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(77)

${\sigma }_{C0}$ denotes an initial damage limit in the uniaxial compression, and $t_0={\sigma }_{C0}/Ea$. The current damage limit ${\sigma }_C$ is assumed as the following function:

$${\sigma }_C\left(t\right)={\sigma }_{C0}-\left(f_C-{\sigma }_{C0}\right)\left[\frac{d_C\left(t\right)}{d_0}-2\right]\frac{d_C\left(t\right)}{d_0}.$$

(78)

The maximum value of the damage limit $f_C$ is attained for $d_0$. The considered system is solved numerically using the NDSolve procedure of Wolfram Mathematica for $E=31\hspace{0.17em}\mathrm{GPa}$, $\nu =0.2$, ${\sigma }_{C0}=1\hspace{0.17em}\mathrm{MPa}$, $f_C=32\hspace{0.17em}\mathrm{MPa}$, $d_0=0.08$, and $a=0.0001\hspace{0.17em}{\mathrm{s}}^{-1}$, $Q=5$, $Q_C=1$. The damage limit is attained at $t_0=0.323\hspace{0.17em}\mathrm{s}$. Figure 12a shows the ${\sigma }_{11}-{\mathsf{\epsilon }}_{11}$ and ${\sigma }_{11}-{\mathsf{\epsilon }}_{22}$ curves in comparison with experimental results [25]. The material constants and the function describing ${\sigma }_C\left(t\right)$ in Equation (78) were assumed on the basis of the relationship between the axial strain and stress, that is, ${\sigma }_{11}-{\mathsf{\epsilon }}_{11}$.

The process of deterioration begins at ${\sigma }_{C0}=0.03\hspace{0.17em}{f}_C$. The experimental results taken from the literature [25,42,43] indicate that both bulk and shear moduli start to decrease for fairly low levels of loading. Let us introduce the following dimensionless secant moduli:

$$\begin{array}{l}{k}^{\sec }=\frac{K_C}{K}=\frac{\frac{E\left(1-d_C\right)}{3\left[1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)\right]}}{\frac{E}{3\left(1-2\nu \right)}}=\frac{\left(1-2\hspace{0.17em}\nu \right)\left(1-d_C\right)}{1-2\nu +2\hspace{0.17em}{d}_C\left(1+\nu \right)},\\ g^{\sec }=\frac{G^{\sec }}{G}=\frac{\frac{E\left(1-d\right)}{2\left(1+\nu \right)}}{\frac{E}{2\left(1+\nu \right)}}=1-d.\end{array}$$

(79)

A comparison of the scaled current shear moduli $g^{\sec }$ is presented in Figure 12b. In experiments, the shear modulus decreases with growing strain [42], which is reflected in numerical results. In Figure 13, we compare the scaled moduli with experimental regression curves obtained in [43]. In Figure 13a, values of $g^{\sec }$ are presented, and for the considered strain range, the compatibility of results is good. The agreement of predictions of the scaled bulk moduli $k^{\sec }$ and experiments seems satisfactory (see Figure 13b). Both functions, $k^{\sec }$ and $g^{\sec }$, diverge from the results adopted from the literature for higher levels of loading. This is due to the fact that after reaching some level of loading, plastic strains strongly influence the behavior of the material and interact with the advancing damage, which is not included in the considered model. The lateral strains (Figure 12a) obtained numerically are in good agreement with experimental results, until plastic dilation starts to play a major role in the process.

4.3. Uniaxial Cyclic Compression–Tension

Now, we investigate the material’s behavior in a uniaxial cyclic test. The primary focus is the transitions from tension to compression and from compression to tension and the changes of stiffness connected to them. All the governing equations for the uniaxial tests have already been given in Section 4.1 and Section 4.2 . Therefore, we do not describe the process in detail. Instead, we concentrate on the graphical representations of the obtained results.

The following cycles and phases are considered (see also Figure 14, Figure 15 and Figure 16, where phases are denoted as ①–⑨):

• Cycle I: the uniaxial compression including: phase 1 (elastic loading for ${\sigma }_{11}>-{\sigma }_C$), phase 2 (loading with active damage process, ${\sigma }_{11}=-{\sigma }_C$), phase 3 (elastic unloading up to ${\sigma }_{11}=0$);
• Cycle II: the uniaxial tension including: phase 4 (elastic loading for ${\sigma }_{11}<{\sigma }_T$), phase 5 (loading with active damage process, ${\sigma }_{11}={\sigma }_T$), phase 6 (elastic unloading to ${\sigma }_{11}=0$);
• Cycle III: the uniaxial compression including: phase 7 (elastic loading for ${\sigma }_{11}>-{\sigma }_C$), phase 8 (loading with active damage process, ${\sigma }_{11}=-{\sigma }_C$), phase 9 (elastic unloading to ${\sigma }_{11}=0$).

According to the assumed strain driven loading program, ${\sigma }_{11}$ stress can change its sign. The strain and stress tensors are diagonal in accordance with Equations (50) and (64). A piecewise linear program of axial strain ${\mathsf{\epsilon }}_{11}$ is assumed to be similar to Equations (56) and (69). The direct strain ${\mathsf{\epsilon }}_{11}$, lateral strains ${\mathsf{\epsilon }}_{22}={\mathsf{\epsilon }}_{33}$, and axial stress ${\sigma }_{11}$ are bound by either Equation (54) or Equation (68). Functions of the strains and the stress with respect to time are shown in Figure 14. The results are delivered using Wolfram Mathematica for $E=20\hspace{0.17em}\mathrm{GPa}$, $\nu =0.2$, ${\sigma }_C=20\hspace{0.17em}\mathrm{MPa}$, ${\sigma }_T=5\hspace{0.17em}\mathrm{MPa}$, and $Q=1$, $Q_T=Q_C=6$.

The axial stress and strain change piecewise linearly throughout the process, while the lateral components develop nonlinearly when the degradation advances (phases 2, 5, and 8). The volume change ($\mathrm{tr}\mathsf{\epsilon }$) is always positive (expansion) for the tension (cycle II) and always negative (contraction) for the compression (cycles I and III).

The compressive damage parameter $d_C$ increases in phases 2 and 8, while the tensile damage parameter $d_T$ evolves only in phase 5. The isotropic damage parameter $d$ grows whenever the degradation process is active (phases 2, 5, 8) (compare in Figure 15). This affects the slope of ${\sigma }_{11}-{\mathsf{\epsilon }}_{11}$ (Figure 16a) for the elastic loading and unloading. That is, the secant stiffnesses change as follows:

$$\begin{gathered} \frac{E^{\sec ,T}\left(t\right)}{E}=\frac{{\sigma }_{11}\left(t\right)}{E\hspace{0.17em}{\mathsf{\epsilon }}_{11}\left(t\right)}=\left(1-d_T\left(t\right)\right)\left(1-d\left(t\right)\right) \mathrm{for} \mathrm{uniaxial} \mathrm{tension}, \\ \frac{E^{\sec ,C}\left(t\right)}{E}=\frac{{\sigma }_{11}\left(t\right)}{E\hspace{0.17em}{\mathsf{\epsilon }}_{11}\left(t\right)}=\left(1-d_C\left(t\right)\right)\left(1-d\left(t\right)\right) \mathrm{for} \mathrm{uniaxial} \mathrm{compression}, \end{gathered}$$

(80)

which is shown in Figure 16b. Let us highlight that the abrupt changes of the stiffness occur when passing through ${\mathsf{\epsilon }}_{11}=0$ ($p=0$ and $\xi =0$), that is, for transitions from phase 3 to phase 4 (compression to tension) and from phase 6 to phase 7 (tension to compression). The discontinuity of stiffness in the graph is caused by the switch between $d_T$ and $d_C$. In an arbitrary test, they cannot increase simultaneously in the regarded point, which ensures the clear distinction between tension and compression. The isotropic damage parameter $d$ connects those two states: the tensile secant stiffness $E^{\sec ,T}$ decreases during the uniaxial tension due to the lowering of $\left(1-d\right)$, and analogously, $E^{\sec ,C}$ drops during the uniaxial compression. Therefore, the damage parameter $d$ captures the history of damage that occurred before a certain instant of time but still affects the actual state of material; even the sign of loading is switched. In Figure 16b, there is an apparent difference between the secant stiffness for the first unloading in the uniaxial compression (phase 3) and the elastic loading in the next compression cycle (phase 7). It is due to the change of $d$ during tension (cycle II). The absolute value of volume change ($\mathrm{tr}\mathsf{\epsilon }$) always grows during loading phases, as pictured in Figure 14 and Figure 17.

4.4. Pure Shear

For the pure shear, the strain and stress tensors and their invariants are as follows:

$$\mathsf{\epsilon }={\mathsf{\epsilon }}_{12}\left(b_1\otimes b_2+b_2\otimes b_1\right), \mathsf{\sigma }={\sigma }_{12}\left(b_1\otimes b_2+b_2\otimes b_1\right) \mathrm{for} {\sigma }_{12}\ge 0,$$

(81)

$$\xi =0, r=\sqrt{2}\hspace{0.17em}{\sigma }_{12}, \theta =\pi /6,$$

(82)

$$p=0, q=\sqrt{2}\hspace{0.17em}{\mathsf{\epsilon }}_{12}=\frac{1+\nu }{\hspace{0.17em}E\hspace{0.17em}\left(1-d\right)}r$$

(83)

and as a result:

$$\mathsf{\epsilon }=\frac{1+\nu }{E\left(1-d\right)}s \mathrm{and} {\mathsf{\epsilon }}_{12}=\frac{1+\nu }{E\left(1-d\right)}{\sigma }_{12}.$$

(84)

We define the strain-driven loading program:

$${\mathsf{\epsilon }}_{12}\left(t\right)=a\hspace{0.17em}t, a>0.$$

(85)

Until the damage limit for the pure shear ${\sigma }_S$ is reached, the material’s response to loading is purely elastic in the form:

$${\sigma }_{12}\left(t\right)=\frac{E\left(1-d\right)}{1+\nu }a\hspace{0.17em}t, d\left(t\right)=0, d_T\left(t\right)=0, d_C\left(t\right)=0 \mathrm{for} t<t_0=\frac{{\sigma }_S\left(1+\nu \right)}{E\hspace{0.17em}a\hspace{0.17em}}.$$

(86)

Further loading leads to degradation. The damage condition (Equation (34)) reduces to:

$${\sigma }_{12}-{\sigma }_S=0.$$

(87)

Using Equation (24), the generalized stresses are obtained:

$${\sigma }_d^{}=\frac{E}{1+\nu }{\mathsf{\epsilon }}_{12}^2, {\sigma }_{dT}^{}=0, {\sigma }_{dC}^{}=0,$$

(88)

and subsequently, the evolution laws (Equation (32)) take the form:

$$\dot{d}=2\lambda \hspace{0.17em}Q\frac{E}{1+\nu }{\mathsf{\epsilon }}_{12}^2, {\dot{d}}_T=0, {\dot{d}}_C=0.$$

(89)

As a result, $d_T\left(t\right)=0$ and $d_C\left(t\right)=0$; thus only the isotropic damage parameter evolves. To get $d\left(t\right)$ for $t\ge t_0$, it suffices to solve the first relationship of Equation (86), taking into consideration the damage condition from Equation (87):

$$d\left(t\right)=1-\frac{\left(1+\nu \right){\sigma }_S}{Ea\hspace{0.17em}t}.$$

(90)

The graphs are obtained for $E=20\hspace{0.17em}\mathrm{GPa}$, $\nu =0.2$, ${\sigma }_S=5\hspace{0.17em}\mathrm{MPa}$, and $a=0.00001\hspace{0.17em}{\mathrm{s}}^{-1}$. Let us notice that we do not need the value of $Q$ for the investigated case. The ratios $Q/Q_T$ and $Q/Q_C$ govern the proportions of suitable damage parameters ($d/d_T$ and $d/d_C$) for tests with coinciding evolutions of two damage parameters. The absolute values of $Q$, $Q_T$, and $Q_C$ do not have a meaning of their own, so it is reasonable to set $Q=1$ and control the processes only through the ratios.

As depicted in Figure 18a, after the damage onset, the secant Kirchhoff modulus decreases compared with the initial value. The reduction is by factor $\left(1-d\right)$, so $d$ can be directly interpreted in the pure shear experiment. When $d\to 1$, the failure occurs—the secant stiffness tends toward zero (compare in Figure 18b). Additionally, for the pure shear the generalized stress ${\sigma }_d$ is the doubled strain energy of the virgin material (see Figure 19).

5. Discussion

In the proposed model, three damage parameters, $d$, $d_T$, and $d_C$, are used. In any active damage process, the isotropic damage parameter changes; that is, $\dot{d}\ne 0$. The isotropic damage parameter $d$ is introduced in a standard way (see in [2,3]), multiplying the Helmholtz energy (strain energy) by $\left(1-d\right)$, which can be directly caught in the pure shear test as a scaling factor for the stiffness of a damaged material. It can also be interpreted as a percent of a specimen’s area that underwent complete failure and thus does not carry any loading. The parameter $d$ affects both the volumetric (bulk) and the distortional part of energy. Its evolution is associated with all stress states that satisfy the assumed damage condition.

On the contrary, two remaining parameters describe damage separately—evolution of $d_T$ excludes the simultaneous growth of $d_C$ and vice versa. They apply only to the strain (stress) states that have nonzero isotropic parts. The reason for including the damage parameters $d_C$ and $d_T$ into the model is to split the material’s response in two cases, so the stiffnesses in compression ($p<0$, $\xi <0$) and tension ($p>0$, $\xi >0$) can be governed individually. The parameters are introduced in such a way that they have clear interpretations in the uniaxial tests: their values control the current stiffnesses of the damaged material and establish unambiguous relationships between the axial and the lateral components of the strain tensor.

Why use three parameters instead of one or two? Using one parameter, $d$ allows for describing only isotropic damage. This is not the realistic case for concrete and rocks, which behave in a more complex way. There is a crucial difference between tension and compression. For those materials, failure in tension is brittle and sudden, and experimental results show a narrow zone of softening [13,15,30]. For compression, the deterioration develops slowly, distinct yielding and hardening are observed, and the failure occurs through crushing [13,30]. The damage (and yield) limits for the uniaxial tension and compression are markedly different. Thus, we need at least two variables to describe the damage process. Models with two damage parameters, in particular the most used concrete damaged plasticity of Abaqus [36], are suitable for describing this dissimilarity but usually use, directly or indirectly, a split of stress or strain tensors into negative and positive parts, which entails differentiability problems. Using three damage parameters reduces this problem: the first derivative of the Helmholtz energy with respect to strain tensor is continuous (the relation between $\mathsf{\sigma }$ and $\mathsf{\epsilon }$ is always unambiguous). The second derivative (the tensor of tangent stiffness) is not determined for $p=0$, but it can be improved by a simple regularization.

Based on the analyses carried out in Section 4, let us make an attempt to associate the variables $d$, $d_C$, and $d_T$ with physical fracture phenomena. For the considered model, there are three damage modes (mechanisms) shown schematically in Figure 20 (compared with two modes described by Ortiz [33]). Mode I, connected to parameter $d$, is the isotropic damage as represented in Figure 20a. Simultaneously, depending on the sign of $\mathrm{tr}\mathsf{\epsilon }$ ($\mathrm{tr}\mathsf{\sigma }$), one of two remaining modes can happen within a material point. In uniaxial tension, the fracture usually occurs in a plane perpendicular to the direction of subjected loading, while for the uniaxial compression, a specimen breaks either parallel or at a slight angle to the loading direction [27,28,31]. Mode II, associated with $d_T$, mirrors the evolution of cracks resulting from the extension coincident with the tensile loading’s direction. An increase in $d_T$ is represented by a black ellipse, while a simultaneous increase in $d$ is depicted by a red circle within the ellipse in Figure 20b. During unloading with material stiffness $\left(1-d\right)\left(1-d_T\right)E$, those cracks partially close, which can be interpreted as regaining the stiffness up to $\left(1-d\right)E$ when the loading is changed to compressive. In Figure 20b, this process is represented by a black ellipse (mode II—opening) contracting to a line (mode II—closure). The red circle is a symbol of isotropic degradation (mode I) that transfers to compressive states. Mode III, associated with the evolution of $d_C$, represents crushing, which is a situation where the compressive loading causes cracks to open in planes parallel to the loading’s direction. Additionally, in this case, the damage can be partly “undone” by switching to tension as a result of rearrangement (wedging) of grains, so mode III is accompanied by mode I, as shown in Figure 20c. In general, the isotropic damage parameter $d$ representing mode I plays the role of a restitution parameter, compare with the concrete damaged plasticity model in Abaqus [12,35,36]. It preserves the memory of the former damage when the sign of loading is changed, which is indispensable in the analyses of cyclic loadings.

Any of the three damage parameters, $d$, $d_T$, and $d_C$, can be easily removed from the model. It suffices to eliminate in Equation (16) the unwanted term connected to one of the coefficients, $Q$, $Q_T$, or $Q_C$, and exclude the damage parameter from energy (Equation (12)) related to it, but their presence seems to be reasoned when comparing with the experimental evidence, at least for concrete. It is operative to set $Q=1$ and establish the ratios $Q/Q_T$ and $Q/Q_C$. This should be performed by comparing the results of numerical calculations for the model with experimental data. That is, having ${\mathsf{\epsilon }}_{11}$, ${\mathsf{\epsilon }}_{22}={\mathsf{\epsilon }}_{33}$, and ${\sigma }_{11}$ from tests, we can use relationships ins Equations (63), (76) and (90) to estimate $d\left(t\right)$, $d_T\left(t\right)$, and $d_C\left(t\right)$ and compare them with the functions received using the considered model. This allows for determining the appropriate levels of the ratios $Q/Q_T$ and $Q/Q_C$.

6. Conclusions

In the paper, a thermodynamically consistent model of an elastic damaged material is proposed. The emphasis is placed on the Helmholtz energy, which is not a quadratic form of the elastic strain tensor. It is required that the function enables a sharp distinction between tensile and compressive states, typical for quasibrittle materials. Compared with the functions present in the literature, the unquestioned advantage of the considered Helmholtz function is its differentiability, which results from omitting the split of the strain (stress) tensor into positive and negative parts. Still, introduction of three damage parameters allows for describing different modes of cracking development. Unfortunately, the assumed energy has an indeterminate second derivative with respect to $\mathsf{\epsilon }$ when $\mathrm{tr}\mathsf{\epsilon }=0$. This problematic flaw in the finite element calculations can be dealt with effectively using a regularization parameter in the tensor of tangent stiffness described by Equation (46).

The presented model is reasonably flexible. The introduction of three parameters in the definition of the dissipation function allows for effectively maneuvering between modes of damage evolution. Due to the proposed dissipation potential, it is possible to obtain any convex damage surface. A selection of damage function would influence the consistency condition, the damage multiplier, and the form of the tensor of tangent stiffness. However, the damage condition given in Equation (34) used in this paper possesses all the features required for concrete. Undoubtedly, the existence of the closed-form calibration formulae (Equations (35)–(37)) is an asset to the model. The evolution laws are constructed so as to guarantee the proportionality of damage parameters’ rates and portions of energy, which is thermodynamically consistent and broadly accepted.

Throughout the text, we refer to the quasibrittle materials, especially concrete, but do not offer extensive comparison with experimental data. This is due to the fact that a full description of such materials requires including plasticity. Elastic damage models not only neglect permanent deformation but also fail to depict dilation. The shown Helmholtz energy and dissipation potential constitute a model of their own, but they can be incorporated in a more complex description of an elastic-plastic damaged material.