Issue 
Mechanics & Industry
Volume 20, Number 1, 2019



Article Number  105  
Number of page(s)  16  
DOI  https://doi.org/10.1051/meca/2018043  
Published online  02 April 2019 
Regular Article
A micromechanical inspired model for the coupled to damage elastoplastic behavior of geomaterials under compression
^{1}
Laboratoire de Mécanique des Solides, École Polytechnique, 91128 Palaiseau, France
^{2}
Institut des Sciences de la Mécanique et Applications Industrielles, UMR EDFCNRSCEAENSTA 9219, Université Paris Saclay, 91120 Palaiseau, France
^{*} email: kyrylo.kazymyrenko@edf.fr
Received:
15
December
2017
Accepted:
13
November
2018
We propose an elastoplastic model coupled with damage for the behavior of geomaterials in compression. The model is based on the properties, shown in [S. Andrieux, et al., Un modèle de matériau microfissuré pour les bétons et les roches, J. Mécanique Théorique Appliquée 5 (1986) 471?513], of microcracked materials when the microcracks are closed with a friction between their lips. That leads to a macroscopic model coupling damage and plasticity where the plasticity yield criterion is of the Drucker–Prager type with kinematical hardening. Adopting an associative flow rule for the plasticity and a standard energetic criterion for damage, the properties of such a model are illustrated in a triaxial test with a fixed confining pressure.
Key words: Plasticity / damage / compression test / geomaterials / dilatancy
© J.J. Marigo and K. Kazymyrenko, published by EDP Sciences 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Many works have been devoted to the modeling of the mechanical behavior of geomaterials like concrete, rocks, and soils and many phenomenological models have been proposed to give an account of the main aspects of the observed phenomena, namely, dilatancy and stress softening, while those materials are submitted to triaxial compressions. The majority of these approaches is based on elastoplastic formalism with some of them having resort to damage coupling. For instance, in soil or rock mechanics, one generally uses elastoplastic or more generally (thermo)elastoviscoplastic models without introducing damage variables [1–3], whereas it is generally admitted that one cannot reproduce the experimental results observed for concrete without considering a damage evolution law [4–10]. In all the cases the elastoplastic models are based on a restricted choice of yield criteria (like Drucker–Prager one or some variants like Hoek–Brown criterion [1,2] or the cap model [11]), with common feature being the hypothesis of a nonassociative flow rule. The reason generally invoked to put aside the normality rule for the plasticity evolution is that otherwise a too important dilatancy effects would be obtained in the evolution of the volumetric strain once the normality rule is used, for instance, for a standard Drucker–Pragertype criterion.
This, evoked earlier, macroscopic elastoplastic coupled behavior is interpreted on a lower scale as the consequence of the presence of microcracks inside the material. Because of the compression, a part or all of those cracks are closed and a friction between the lips of the cracks can prevent a free sliding of the lips. This impossibility to relax to the initial strain state results at the macrolevel in the presence of residual plastic strains. Moreover, if one assumes that the sliding with friction between the lips follows Coulomb law, then a Drucker–Pragertype criterion for the plasticity law is naturally obtained. Finally, the growth of the microcracks is in straight relation to the observed progressive loss of rigidity of the samples and can be reproduced at the macrolevel by using a damage model. But if those mechanisms related to the microcracking are usually invoked to build macroscopic models, very few works have been devoted to an authentic micromechanical approach. An exception is reference [12] where such an attempt is made by using Coulomb friction law for the sliding of the crack lips and Griffith law for the crack propagation; see also a more recent paper [13] where a full 3D approach is proposed. Even if the analysis in reference [12] is conducted in a simplified twodimensional setting where the uncracked material is assumed isotropic, elastic, and homogeneous whereas the microcracks are assumed straight and small enough to neglect their interactions, some general properties are obtained and can be used to construct more general macroscopic models. In particular, the authors show that the elastic energy stored in the microcracked material reads as the following function of the total strain and the plastic strain :(1)
In equation (1), one uses the notations precised at the end of this section. In particular, denotes the stiffness tensor of the sound material, whereas represents a loss of stiffness due to the microcracks. Therefore, the micromechanical approach leads to a term in the elastic energy which corresponds to the energy blocked by the contact with friction of the lips of the cracks. Let us note that this latter term involves only the plastic strain (and eventually the damage state via dependence) not the total strain. The damage state is supposed not to enter in the first term of the energy equation (1) and, consequently, even if the stress–strain relation still reads asthe stress which appears in the plasticity yield criterion deduced from Coulomb friction law is not but the tensor given by
In other words, the yield criterion is of the Drucker–Prager type with a kinematical hardening where the back stress depends linearly on the plastic strain but also nonlinearly on the damage state. These two particularities (the presence of a kinematical hardening in the plasticity yield criterion and the fact that damage appears only in the blocked elastic energy) can generally not be seen in the various models proposed in the literature; see however references [10,14] where kinematical hardening is also introduced. It turns out that they have fundamental and rather unexpected consequences on the response of the material under triaxial compression tests. In particular, one can account for reasonable contractance and dilatancy effects without leaving an associate flow rule. The goal of the present paper is to study the specificity of the responses of this type of models.
The paper is organized as follows. In Section 2 we con sider the elastoplastic model without damage. Presenting first the general ingredients of an associative plasticity law, we progressively refine the model by introducing the particularities coming from the micromechanical considerations. We finally study the response of a volume element submitted to a triaxial test at fixed confining pressure. The damage is introduced in Section 3, still on the basis of the micromechanical considerations. Using a standard law for the damage evolution law, we obtain the complete model of elastoplasticity coupled with damage. After establishing some general properties, we then consider a family of particular models to finally calculate the response of the volume element submitted to the triaxial test at fixed confining pressure.
Throughout the paper, we use the following notations (see also Tab. 1). The summation convention on repeated indices is implicitly adopted. The vectors and secondorder tensors are indicated by boldface letters, like and for the unit normal vector and the stress tensor. Their components are denoted by italic letters, like n _{ i } and σ _{ ij }. The fourthorder tensors as well as their components are indicated by a sans serif letter, like or for the stiffness tensor. Such tensors are considered as linear maps applying on vectors or secondorder tensors and the application is denoted without dots, like whose ijcomponent is . The inner product between two vectors or two tensors of the same order is indicated by a dot, like which stands for a _{ i } b _{ i } or for σ _{ ij } ε _{ ij }.
Notation used throughout the paper for main mechanical quantities.
2 The elastoplastic behavior at fixed damage state
This section describes a pure elastoplastic behavior law. Even so the damage is not explicitly present in this formulation, the latter could be considered as the model at fixed damage state. The idea is to identify all the relevant coupling terms which influence dilatancy and softening behavior.
2.1 General considerations
An associative model of elastoplasticity is defined by two potentials:

The volume free energy function of the state variables.

The dissipation potential or equivalently the convex set in which the thermodynamical forces associated with the internal variables must lie.
Here we will only consider isothermal processes and a plastic behavior with linear kinematical hardening. Accordingly, the state variables are constituted by the strain tensor and the plastic tensor without other internal variables (the damage variable will be introduced in Section 3). Therefore, the free energy density ψ is given by(2)where the state function is assumed to be convex and at least continuously differentiable. By differentiation, one obtains the stress tensor and the tensor of the thermodynamical forces associated with the plastic strain tensor:(3)
From the mathematical point of view, , , , and are order 2 symmetrical tensors which can be identified with 3 × 3 symmetrical matrices after the choice of an orthonormal basis. In other words, those tensors will be considered as elements of .
The plasticity (or yield) criterion is defined by giving the (closed with nonempty interior) convex set of where the thermodynamical forces must lie. Here we make the strong assumption that this domain is fixed, that is independent of time and of the plasticity evolution. (However, we will see in the next section that the image of this domain in the stress space evolves with the plastic strain.) In practice, the convex is characterized by a convex function so that the plasticity criterion reads as(4)
The interior points of correspond to the forces such that and the points on the boundary of to the forces such that .
Since we only consider an associative model, the evolution of the plastic strain follows the normality rule. When the convex set has a smooth boundary without angular points, the function f is differentiable and the flow rule can read as(5)where denotes the plastic multiplier. If the boundary of is not smooth and the normal is not defined at some points, then the normality rule is extended by using the Hill maximal work principle. In such a case, the flow rule is given by(6)
The inequality (6) says that must belong to the cone of the outer normals to at the point of the boundary. That inequality leads to equation (5) when the normal is well defined at .
From the energetic point of view, the energy which is dissipated by plasticity can be defined in terms of the support function of the convex . Specifically, let us introduce the support function of :(7)
By construction, is convex and positively homogeneous of degree 1,
Moreover, vanishes at , is nonnegative provided that , and can take the value +∞ when is not bounded. Its interpretation as the volume dissipated power by plasticity comes from Clausius–Duhem inequality. Indeed, in isothermal condition, the dissipated power (by volume unit) is defined byand the second principle of thermodynamics requires that be nonnegative. Owing to equations (2) and (3), reads as
By virtue of Hill maximal work principle (6) and the definition (7) of the support function, one gets(8)which means that plays also the role of the dissipation potential. Furthermore, it is sufficient to have for Clausius–Duhem inequality to be satisfied, i.e., .
2.2 Choice of the form of the free energy state function
Guided by the micromechanical considerations presented in reference [12] for microcracked materials, we consider that the plastic strain is due to the friction between the lips of the (closed) microcracks. However, we do not assume that and hence abandon the usual plastic incompressibility condition. Consequently, the plastic strain is a full symmetric secondorder tensor and no more a pure deviator. Moreover, the spherical part of which comes from the normal displacements between the lips of the microcracks will be also governed by an irreversible evolution law. Accordingly, we choose for the free energy the following quadratic function of :(9)
In equation (9) the first term represents the elastic energy and A _{0} denotes the (fourthorder) stiffness tensor of the sound material. The second term in equation (9) represents the blocked energy by friction and the (fourthorder) stiffness tensor A _{1} represents the loss of stiffness due to the presence of microcracks. As one will see later, this loss of stiffness is visible only when the sliding or the opening between the lips of the cracks are active, i.e., when .
The tensor A _{0} is positive definite and the tensor A _{1} is assumed to be nonnegative. In fact, A _{1} is also positive definite as soon as the material is microcracked, as it is shown in reference [12]. But, we will also consider the case when A _{1} vanishes to emphasize the strong influence of the blocked energy on the mechanical behavior of the microcracked material. Consequently, is a (strictly) convex function of .
To simplify the presentation, we will assume that the material remains isotropic (even when the damage grows). Therefore, both tensors A _{0} and A _{1} have only two independent moduli. Decomposing the strain and the plastic strain tensors into their spherical and deviatoric parts(10)the free energy finally reads as(11)where K_{0} > 0 and μ _{0} > 0 are the compressibility and shear moduli of the sound material, whereas K_{1} ≥ 0 and μ _{1} ≥ 0 will be called the kinematical hardening moduli. The elastic moduli K_{0} and μ _{0} are related to the Young modulus E_{0} and the Poisson ratio ν _{0} of the sound material by(12)
By differentiation of equation (11), the stress–strain relation is linear and reads as(13)
Decomposing the stress tensor into its spherical and deviatoric parts,(14)the relation (13) becomes(15)
Still by differentiation, the thermodynamical forces associated with the plastic strains read as(16)and they differ in general from the stresses when A _{1} ≠ 0. Decomposing the forces X into their spherical and deviatoric parts,(17)the relation (16) becomes(18)
2.3 Choice of the plasticity criterion
Drucker–Prager criterion (or some variants like Mohr–Coulomb or Hoek–Brown criteria) is widely used in the modeling of the behavior of geomaterials under compression, see references [1,2,5]. In general, those criteria are used in a nonstandard context, i.e., with a nonassociative flow rule for the plasticity evolution, like in reference [3]. The main usual argument which is evoked to use a nonassociative flaw rule is that one overestimates the dilatancy effect with the normality rule associated with Drucker–Prager criterion. On the contrary, we propose in this paper to keep an associative flaw rule but with a Drucker–Prager criterion which contains a kinematical hardening. One of the main goals is to show how a kinematical hardening can lead to relevant dilatancy effects.
Specifically, in the present context, the Drucker–Prager criterion reads in terms of the thermodynamical forces X as(19)
In equation (19) k ∈ (0, 1) is a dimensionless coefficient associated with the internal friction, τ _{ c } ≥ 0 represents a critical stress which vanishes in absence of cohesion, and ∥· ∥ denotes the euclidean norm of a secondorder tensor (whereas the factor is introduced for convenience in order to simplify forthcoming expressions),
Consequently, the elastic domainis really a closed convex subset of (with nonempty interior). In fact, is a convex cone which has an angular point at and which is unbounded in the direction of hydrostatic compressions. Furthermore, even if the elastic domain is fixed in the space of the thermodynamical forces, it varies with the plastic strain in the stress space. Indeed, by virtue of equation (16), the elastic domain expressed in term of becomes the set given by
Thus, is subjected to a linear translation when the plastic strain evolves:
Starting from equation (7), the support function of the convex is obtained after tedious calculations which are not reproduced here and finally reads as(20)
Since we adopt the normality rule, the evolution of the plastic strain is directly deduced from equation (19). Considering first a regular point of the boundary of , the flow rule reads as follows:(21)
Let us note that at such a point one has . At the angular point, the rate of the plastic strain tensor must belong to the cone of the outer normal to . That leads to(22)
Consequently, the normality rule forces the trace of the plastic strain to only increase with time, a fundamental property to account for the dilatancy effects.
2.4 Response under a triaxial test with a confining pressure
The elastoplastic behavior of a material governed by the associative Drucker–Prager law presented in the previous sections is illustrated by considering the response of a volume element during a triaxial test with a confining pressure. Assuming that the volume element starts from a natural reference configuration without plastic strain, the test is divided into the following three stages:

First, the volume element is submitted to a hydrostatic compression where the pressure is progressively increased from 0 to a final value p _{0 }> 0.

Then, maintaining the lateral pressure to the value p _{0}, one compresses in the axial direction z by prescribing the axial strain ε _{ z } with . Accordingly, the response will be purely elastic as long as ε _{ z } is small enough so that the plastic yield criterion is not reached.

Finally, wheε _{ z } is larger than a critical value, the volume element plastifies.
Let us determine the evolution of the strains, the plastic strains, and the stresses during those different stages with respect to ε _{ z } (which could be considered as the loading parameter).^{1}
2.4.1 Confining stage
During this stage the strain and stress tensors are purely spherical and the plastic strain remains equal to 0. At the end of this stage, the stresses and the strains are given by(23)
Since p = 0, one has and DuckerPrager criterion gives
The axial strain ϵ _{ z } decreases from 0 to −p _{0}/3K_{0}, whereas the volumetric strain decreases from 0 to −p _{0}/K_{0}. During all this stage, one hasand hence a contractance due to the hydrostatic compression.
2.4.2 Elastic stage at fixed confining pressure
During this stage, the plastic strain remains equal to 0:
By symmetry, the stress and the strain tensors are of the following form^{2}:where p _{0} and ε _{ z } are prescribed. The deviatoric parts read aswith σ _{ z } and ε _{ T } remain to be determined. They are given by equation (15) with the help (12):(24)
One deduces in particular the evolution of the volumetric strain:(25)
This stage stops when the stresses reach the yield criterion, i.e., when . Sinceand since σ _{ z } + p _{0} < 0, the yield surface is reached when(26)
The corresponding value of ϵ _{ z } is deduced from equation (24):(27)
2.4.3 Plastification stage
Then, if the confining pressure is constant and if one prescribes , the material will plastify. This plastification stage strongly depends on the presence or not of the kinematical hardening. Therefore, we will distinguish between the two situations by considering first the case where the hardening moduli μ _{1} and K_{1} vanish, i.e., the case without hardening.

Case without kinematical hardening. In such a case, one gets and since the elastic domain is fixed in the stress space, the stresses will remain blocked to the value reached at the end of the elastic stage. In other words, the lateral pressure is equal to p _{0} and the axial stress σ _{ z } is given by equation (26). Since , one getsand hence p takes the form(28)
Since σ _{ z } + p _{0} < 0, the flow rule (21) gives
Therefore, the volumetric strain rate is given by(29)
Thus, the normality rule predicts a dilatancy of the material during the plastification stage. In general, this predicted dilatancy is too large by comparison with experimental observations and that leads people to use a nonassociative flow rule. We will see just below that the kinematical hardening can strongly modify these results.

Case with kinematical hardening. In such a case, the stresses do not remain blocked. Their evolution is obtained by assuming that the strain tensor remains of the form (28) and that the yield criterion is always satisfied at a nonangular point, i.e.,
Let us note that this assumption is not restrictive because we can use the general result of uniqueness of the response in the case of an associative law with kinematical hardening. Indeed, if we find a solution, then we are sure that is the good one. (Note that such a uniqueness result is no more ensured in the case of nonassociative flow rule.)
Differentiating these relations with respect to time and taking into account equation (18) lead to
The flow rule (21) gives(30)
Combining the four previous relations gives the following relationship between and :(31)
Differentiating the stress–strain relations (13) gives
Eliminating and using the flow rule give two other relations between , , , and :(32)
Using equation (31) allows us to determine , , and in terms of . We finally get(33) (34)
Therefore, equation (33) shows that the kinematical hardening tends to decrease and even equation (34) shows that the sign of depends on the intensity of the hardening. Specifically, if the hardening moduli are small enough (by comparison with the elastic moduli), then and there is dilatancy. But, if the hardening moduli are large enough, then and there is con tractance. Moreover, let us note that if μ _{1} = K_{1} = 0, then we recover that and the relation (29). But, if µ _{1} and K_{1} are very large by comparison to µ _{0} and K_{0}, then one gets
which is nothing but the relation (25) obtained in the elastic stage for the volumetric strain change. In conclusion, according to the ratio between the hardening moduli and the elastic moduli, the volumetric strain evolution can lead to a strong dilatancy effect corresponding to perfect plasticity as well as a contractance effect like in elasticity. This last conclusion leads to a foretaste consideration (developed in detail in the next section) for a relevant damage coupling − once the damage affects A _{1} tensor, it directly impacts the dilatancy strength of the given law.
3 The elastoplastic behavior coupled with damage
In this section we introduce the damage into the elastoplastic model. As mentioned at the end of Section 2, the coupling could be introduced purely from macroscopic considerations. We made a choice for more physical approach starting from a set of hypothesis which will finally lead to the same coupling. Let us first present a list of assumptions used in that construction before particularizing the model.
3.1 Main assumptions on the damage dependence
Although there exist an infinite number of ways to couple damage with plasticity [5–9], we will make the simplest choices. Specifically, the following are the main assumptions:

The damage state is characterized by a scalar variable α which can grow from 0 to 1, α = 0 corresponding to a sound material, and α = 1 corresponding to a full damaged material.

The damage state enters in the expression of the free energy, but not in the plasticity criterion. In other words, the convex domain is independent of α.

The dissipated energy by damage is only function of the damage state.
By virtue of these assumptions, the state of the material point is characterized by the triple and the free energy is still a function of state which now be read as(35)
The function is still assumed to be continuously differentiable from which one deduces the stresses and the thermodynamical forces X and Y associated with the plastic strain and the damage. Specifically, one sets(36)
In equation (36), Y can be interpreted as the free energy release rate associated with a growth of the damage at constant strain. Since it is natural to assume that the release of energy is really nonnegative, one adds the following condition of positivity for Y:(37)condition which restricts the dependence of the free energy on α.
Since the plasticity criterion is assumed to be independent of α (hypothesis (2)), the plastic dissipation potential still reads aswhere is a closed convex subset of . Finally, by hypothesis (3), the damage dissipated energy d reads as(38)where is a nonnegative differentiable function of α which vanishes at α = 0:
This assumption is also motivated by the micromechanical approach of reference [12]: indeed, if one considers microcracked materials whose microcracking evolution is governed by Griffith's criterion, then dissipated energy by damage corresponds to the Griffith surface energy.
Therefore, the part of the dissipated power due to the evolution of damage simply reads in terms of the damage state and the damage rate as
whereas the total dissipated power due to the evolution of both the plastic strain and the damage is given by
It remains to formulate the damage evolution law. As for the elastoplastic behavior, we adopt a standard law in the sense of the concept of Generalized Standard Materials proposed in reference [15] or [16] in a general context and particularized by [17,18] for damaging materials. Specifically, the evolution of the damage consists in the following three items:

The irreversibility condition: Damage can only grow and hence one requires that at any time.

The damage criterion: The free energy release rate must remain less or equal to the critical value given by the damage dissipated power:(39)

The consistency relation: Damage can only evolve when the free energy release rate is equal to its critical value, condition which can read as
3.2 The associative elastoplastic Drucker–Prager model with kinematical hardening coupled with damage
Let us particularize now the model that we will use for geo materials in compression.
3.2.1 Choice of the free energy and of the dissipated energy by damage
Let us first introduce the damage variable into the expression (9) of the free energy. The micromechanical approach developed in reference [12] for a microcracked material with contact and friction between the lips of the cracks shows that only the hardening tensor A _{1} depends on the crack state, the tensor A _{0} representing the stiffness of the sound material. Therefore, we will assume that A _{1} only depends on the damage variable and hence the expression of the free energy becomes(41)
Moreover, if we assume that the material is isotropic, then one gets(42)
The positivity of the stiffness tensor A _{0} and of the hardening tensor A _{1} requires that K_{0}, µ _{0}, K_{1}(α), and μ _{1}(α) satisfy the following inequalities:
The energy release rate reads as(43)and hence it depends only on the plastic strain and the damage, and not on the total strain. In order that Y be nonnegative, refer to the condition (37), it is necessary and sufficient that
The damage criterion becomes then
Since the Drucher–Prager yield criterion is assumed to be independent of damage, the model will be complete once the three functions α ↦ K_{1} (α), α ↦ μ _{1}(α), and are given for α ∈ [0, 1]. In fact since the choice of the damage variable is arbitrary, it is always possible, after a suitable change of variable, to fix one of the three functions. Here we choose to fix the function by setting(44) where D_{1} is a positive constant which has the dimension of a stress (or equivalently an energy by volume unit). Accordingly, the dimensionless damage variable α can be interpreted as the fraction of the volume energy dissipated by damage at the current damage state by comparison to D_{1} which represents the volume dissipated energy by damage when the volume element is totally damaged.
3.2.2 The complete model
Finally, the associative elastoplastic Drucker–Prager model with kinematical hardening coupled with damage consists in the following set of definitions and relations:
Remark 3.1 This model can be written in a variational form following the general presentation proposed in reference [19] for rateindependent evolution laws. After the choice of the expression of the total energy, sum of the free energy, and the dissipated energy, as a function of the state variables, the evolution of the internal variables is deduced from three general physical principles: (i) an irreversibility principle, (ii) a stability criterion, and (iii) an energy balance. The interested reader for such an approach could refer to references [20–22] where its application to elastoplasticity and damage models is developed. That approach can even be used to construct nonassociative plasticity model, see reference [11] for an example. Note also that the model has a stresssoftening character induced by the damage evolution. That means that the present local model must be regularized (for instance, by introducing gradient of damage terms in the energy) in order that the damage localization in a body be limited and controlled. For the construction of the regularized model, the variational approach should be really useful and even necessary.
3.3 Application to the triaxial test with containing pressure
3.3.1 Some general properties
Let us reconsider the triaxial test by including now the possibility of damage of the volume element. The results obtained in Section 2.4 remain valid as long as the plasticity yield criterion is not reached. In other words, the confining stage and the elastic stage remain unchanged. What happens after depends in particular on the choice of the functions μ _{1}(α) and K_{1}(α). However we can obtain some general properties without specifying those functions provided that they satisfy the following conditions:(45) (46)
Those conditions are inspired by the micromechanical model presented in reference [12]. In particular, the fact that the hardening moduli are infinite when α = 0 ensures that the elastic moduli K_{0} and µ _{0} are those of the undamaged material. In the same manner, the fact that the hardening moduli vanish when α = 1 corresponds to the total loss of stiffness when the material is fully damaged. The hypothesis that the second derivatives of K_{1} and µ _{1} are positive plays an important role as we will see below. Let us now establish some general properties based on the above assumptions.

Damage evolves only with plasticity. Let us first show that damage can grow only when plasticity evolves. Indeed, by virtue of the damage criterion, damage evolves only when Y = D_{1}, i.e., when
Differentiating this relation with respect to t leads to(47)from which one deduces that if , which is the desired result. Let us note that this result is essentially due to the fact that the stiffness tensor A _{0} is assumed to be damage independent.

Damage starts at the same time as plasticity. Indeed, by virtue of the hypothesis , the damage criterion at the onset of damage givesand hence can be satisfied only if p = 0.

Form of the strain, plastic strain, and stress tensors during the plasticity stage with damage. The previous properties suggest to search a response such that damage and plasticity evolve simultaneously. During this damage with plasticity stage, the strain and stress tensors are assumed to be of the form^{3} where p _{0} remains constant and ϵ _{ z } is a given function of time which is decreasing from its critical value given in reference (27).

Relationship between and α. Inserting the form of the different tensors into the plasticity flow rule (which does not involve the damage state) gives(48)
Reporting this relation into the damage criterion gives in function of α:(49) where R (α) is the following combination of the hardening moduli:(50)
Let us note that, since the second derivatives of K_{1} and µ _{1} are positive, is a strictly increasing function of α. That means that the required inequality is satisfied provided that .

Relationship between σ _{ z } and α. The plasticity criterion with the definition of the plasticity thermodynamical forces give the following three equations:
Eliminating X_{T} , X_{z} , X_{m} and using equation (49) allow us to obtain σ _{ z } as the following function of α:(51)where
One then checks that σ _{ z } is always negative. Its absolute value will be either increasing or decreasing when α grows according to the second derivative of the compliance function S (α) which is negative or positive. The former case, when , corresponds to a stresshardening response, whereas the second one, when , corresponds to a stresssoftening response. Note that the sign of the second derivative of S(α) is not given by the sign of the second derivative of R(α) and hence constitutes an additional choice for the model.

Relationships between ϵ _{ z }, ϵ _{ v }, and α. The stress–strain relations allow us to express ϵ _{ z } and ϵ _{ v } in terms of α. Specifically, taking into account equations (49) and (51), one gets(52) (53)

The relation between the axial stress σ _{ z } and the axial strain ϵ _{ z }. The relations (51) and (52) give the relation between ϵ _{ z } and σ _{ z } under the form of a curve parametrized by α. Let us note that we are not ensured that ϵ _{ z } is a monotonic function of α. Specifically, when (stresshardening behavior), then when and hence damage grows when ϵ _{ z } (which is negative) is decreasing. But, when , then the sign of is not guaranteed when , because the term in σ _{ z } is increasing, whereas the term in is decreasing. Therefore, a snapback is possible. If such a snapback exists, the control of the axial strain ϵ _{ z } necessarily implies a discontinuous evolution of the damage. That can even lead to a brutal rupture of the volume element with a sudden jump of α to 1 when ϵ _{ z } reaches its limit value.

The relation between the volumetric strain ϵ _{ v } and the axial strain ϵ _{ z }. The relations (51)–(53) give the relation between ϵ _{ v } and σ _{ z } under the form of a curve parametrized by α. Here also one can have a competition between the term in σ _{ z } and the term in . During the hardening phases, one has and . Consequently, according to the respective weight of the two terms, one will observe either a contractance or a dilatancy. On the contrary, during the softening phases, since and , one will always observe a dilatancy.
These different properties are illustrated in the next section by particularizing the model. Let us note that in the triaxial test with a fixed confining pressure, the functions μ _{1}(α) and K_{1}(α) appear only by their com bination R(α) (which involves also the internal friction coefficient k). Therefore, it is sufficient to know the function R(α) and the constants µ _{0}, K_{0}, D_{1}, k, and τ _{ c } to characterize the model, the confining pressure p _{0} playing the role of a parameter.
3.3.2 Study of a family of models
Let us consider the following family of functions R(α):(54)
This choice of function is governed mainly by our intention to control the asymptotic behavior both at the initiation of damage (α ≈ 0 power law of n) and at the end of microcracking (α ≈ 1 power law of m). One easily checks that such a function α ↦ R(α) is really compatible with the hypotheses (45) and (46). Indeed, R (α) is strictly decreasing from +∞ to 0 when α grows from 0 to 1. Its first and second derivatives are given by
It is also not difficult to verify that for every α ∈ (0, 1) because m > 1 > n > 0. The first derivative R^{′}(α) is strictly increasing from −∞ to 0 when α grows from 0 to 1. By virtue of equation (49), one deduces that grows from 0 to +∞ with the evolution of damage.
The compliance function S (α) and its first and second derivatives are given by See equation above.
The study of the sign of S^{′′}(α) shows that S^{′′}(α) is negative in the interval (0, α _{0}), positive in the interval (α _{0}, 1) with α _{0} given by(55)
Therefore, by virtue of equation (51), the absolute value of the axial stress σ _{ z } increases and one observes a stresshardening behavior when α grows from 0 to α _{0}, then σ _{ z } decreases and one observes a stresssoftening behavior when α grows from α _{0} to 1. Moreover, since S^{′} (0) = S^{′}(1) =+ ∞, σ _{ z } starts from the value given by equation (26) and corresponding to the end of the elastic stage, then increases (in absolute value) until its maximal value corresponding to the time at which the damage reaches the value α _{0}, and finally decreasing (in absolute value) to come back to the value (26) reached at the end of the elastic stage. The absolute value of the overstress Δσ _{ z } associated with the stresshardening phase is given by(56)
Let us first study the variations of the axial strain ϵ _{ z } in function of the damage α and let us analyze under which condition a snapback occurs. Starting from equation (52) and using the fact that ϵ _{ z } is negative, one deduces that the derivative of ϵ _{ z } with respect to α is always negative. Therefore, a snapback exists if and only if the following condition is satisfied:(57)
Since R (α) is of the form , the lefthand side above can read as(58)where the function α ↦ φ (α) depends only on the exponents m and n of the model. In the neighborhood of α = 1, the different terms behave as follows:
Therefore, since φ is continuous, negative in the interval (0, α _{0}), positive in the interval (α _{0}, 1), and since φ (α _{0}) = φ (1) =0, the function φ reaches its upper bound in the interval (α _{0}, 1), the upper bound depending only on m and n. Hence, the condition of nonsnapback (57) can read as(59)
In other words, there is no snapback provided that the hardening modulus R_{1} is small enough. Let us note that this condition depends on the parameters E_{0}, K, m, and n of the model, but is independent of D_{1}, τ _{ c }, and the confining pressure p _{0}.
Let us now study the variations of the volumetric strain ϵ _{ v } as a function of the damage α and let us analyze when one observes a contractance or a dilatancy. We assume that R_{1} is small enough and satisfies equation (59) so that there is no snapback in the response σ _{ z } − ϵ _{ z } and hence . By virtue of equation (53), one has and hence con tractance if and only if
With the help of equations (49) and (51), that condition can be expressed in terms of φ (α) defined by equation (58). Specifically, one observes a contractance at the time when the damage reaches the value α if the following condition is satisfied:(60)
That can happen only during the phase where φ is negative, that is only during the stresshardening phases. During the stresssoftening phases, since φ > 0, one observes a dilatancy. Therefore, a contractance can happen only at the beginning of the test, as long as α < α _{0}. Since the behavior of the first and the second derivatives of R and S in the neighborhood of α = 0 is given byone gets
Consequently, since φ (0) =−∞, one obtains a contractance when α is small enough. On the contrary, as soon as α becomes larger than a critical value (which is necessarily less than α _{0}), one observes a dilatancy until the complete rupture of the volume element. Let us note that since tends to 0 when α tends to 1, one gets by virtue of equations (52) and (53)
when α is close to 1. That is in conformity with the result (29) obtained in the case of a perfect plasticity behavior without kinematical hardening, situation to which one tends when α becomes close to 1.
Let us finally study the relation between the axial strain ϵ _{ z } and the axial stress σ _{ z } in the case when the condition (59) of nonsnapback is satisfied. Using the behavior of R^{′} and S^{′} in the neighborhood of α = 0, one deduces from equations (49) and (51) the following behaviors for and σ _{ z }:
Thus, the plastic strain is on the order of α ^{ (1+n)/2} while the variation of the elastic strain is on the order of α ^{ (1−n)/2}. Hence, at the beginning of the plasticity stage, the plastic strain is negligible by comparison with the variation of the elastic strain. Inserting these estimates into equation (52), one deduces that the relation between ϵ _{ z } and σ _{ z } is continuous and continuously differentiable at the transition between the elastic stage and the plasticity stage, the slope being equal to E_{0}. After the elastic stage, if the axial strain ε _{ z } is increased until +∞, then σ _{ z } begins by increase, and then passes by a maximum before to decrease to finally take the value that it had at the end of the elastic stage. All those variations are in fact independent of the confining pressure.
3.3.3 Graphs of the response according to the values of the parameters m, n, and R_{1}
In all the pictures presented in this section, one uses the following values for the Poisson ratio ν _{0} and the parameters τ _{ c } and k of the Drucker–Prager criterion:
The material constant σ _{ c } fixes the scale of the stresses, whereas the material constant ϵ _{ c } fixes the scale of the strains. The responses depend on the two exponents m and n, and on the two ratios R_{1}/E_{0} and p _{0}/σ _{ c }.

Case where m = 2, n = 1/2. The value α _{0} of the damage from which there is stresssoftening is given by equation (55) and hence one finds α _{0} = 0.0883. By virtue of equation (56), the overstress due to the stresshardening phase depends on R_{1}/E_{0} (but not on the confining pressure p _{0}) and is given by(62)
In order that there is no snapback in the response σ _{ z } − ε _{ z }, it is necessary that the condition (59) be satisfied. For the considered values of m and n, one gets and hence the nonsnapback condition requires that R_{1} ≤ 0.275E_{0}.
Figure 1 shows the graphs of σ _{ z } and ϵ _{ v } as functions of ϵ _{ z } when there is no confining pressure (p _{0} = 0) and in the case where R_{1} = 0.1E_{0} (hence there is no snapback). In that case the value of the overstress is .
In Figure 2 are plotted the same graphs for different values of R_{1}/E_{0}. One can note that the larger the R_{1}, the larger the overstress, which conforms to equation (62), but also the more rapid is the decrease of the axial stress after the peak. A similar behavior can be seen for the volumetric strain: the larger the R_{1}, the greater the maximal contraction, but also the more rapid the growing of the dilatancy after the peak.
When R_{1} is greater than 0.275E_{0}, the graph of σ _{ z } contains a snapback which induces, during a test where the axial strain is controlled and continuously decreasing, a discontinuity of the evolution of the axial stress and the damage (cf. Fig. 3).
The graph of the axial stress vs. the axial strain for a given confining pressure during the damageplasticity phase is simply a translation of that corresponding to the uniaxial compression test, the translation being given by the values of the axial strain and the axial stress at the end of the elastic stage (cf. Fig. 4). This property is due to the linear character of Drucker–Prager criterion. That would not be true anymore if one used a more general plasticity criterion like Hoek–Brown criterion [1,2].

Influence of the exponents m and n. The shape of the graphs giving the axial stress in terms of the axial strain (during an uniaxial compression test without confining pressure) is always the same whatever the values of m and n provided that they remain inside the admissible intervals, i.e., m > 1 > n > 0. One can see the dependence of the response on m an n in Figure 5.
We note that as expected the m parameter influence mainly the postpic stress behavior (for m = 1 the ultimate damage state achieved even for the finite value of ϵ _{ z }), while n parameter alternate mainly the initial hardening state leaving the large ε _{ z } asymptote intact. Both parameters (m and n) contribute (Eqs. (55) and (56)) to the overstress Δσ _{ z } estimation that reads as
where ϕ (m, n) is only dependent on m and n. At given n, ϕ (m, n) is a decreasing function of m, whereas at given m, ϕ (m, n) is first a decreasing function, then an increasing function of n when n varies from 0 to 1 (cf. Fig. 6).

Possible parameter fitting strategy. Even so our model depends just on three “damage” parameters (R_{1}, n, m), the common task of parameter choice for the given material could be a rather laborious exercise. Nevertheless, based on the obtained earlier parameter influence properties, one could develop straight experimental data fitting strategy. From the experimental data set characterizing geomaterials, one has usually access to dilatancy () and softening curves (σ _{ z } (loading)). The softening postpic behavior would allow us to fit first the m value. Once the m parameter is fixed, one could proceed to R_{1} value choice which allow adjusting the strength of dilatancy curve. Finally, the experimental overstress threshold could be fixed by adapting the last n parameter.
4 Conclusion and perspectives
Starting from lowscale mechanical properties of microcracked elastic materials and taking into account the closed cracks friction between its lips, we propose a phenomenological model with an associative elastoplastic behavior coupled with damage. The main feature of the model which is inspired by the micromechanical considerations is that the free energy contains not only the usual elastic energy but also a stored energy which is due at the microlevel to the friction between the lips of the microcracks. This blocked energy is assumed to depend at the macrolevel on the damage state and the plastic strain and induces a kinematical hardening in the plasticity yield criterion. Considering a Drucker–Prager yield criterion and adopting the normality rule for the plasticity flow rule, we have applied such a model to a triaxial test with a fixed confining pressure. It turns out that the model is able to account for the main observed properties of geomaterials in such a situation, like contractance or dilatancy effect on the evolution of volumetric strain and stresshardening or stresssoftening effect for the axial stress.
The next step will be to identify from experimental results the parameters of the proposed model (which can be enriched by considering more complex yield criterion than the Drucker–Prager one). Its validity will be also tested by comparing its predictions with other families of experimental tests like the oedometric test or cyclic tests. An interesting issue is also to use such a model to calculate the response of a true threedimensional sample (and no more of the volume element). Since the damage evolution is accompanied by stresssoftening effect, one can expect that the response of the sample is no more homogeneous (in space). That could require to enhance the damage model by introducing nonlocal terms like in references [23,24].
Fig. 1 Graphs of the axial stress σ _{ z } (left) and of the volumetric strain ϵ _{ v } (right) versus the axial strain ϵ _{ z } when R_{1} = 0.1E_{0} and p _{0} = 0 (uniaxial compression test). 
Fig. 2 Influence of R_{1}/E_{0} on the graphs of the axial stress and the volumetric strain versus the axial strain when there is no confining pressure (uniaxial compression test). The value of R_{1}/E_{0} associated with each graph is indicated in red on the graph. 
Fig. 3 When R_{1}/E_{0} is large enough, the graph of σ _{ z } vs. ϵ _{ z } contains a snapback (here R_{1} = 0.4E_{0}). Consequently, when the axial strain reaches the value corresponding to the limit point A, the response jumps from A to B. 
Fig. 4 Graph of the axial stress vs. the axial strain for different values (in red on the graphs) of the confining pressure (in fact of the ratio p _{0}/σ _{ c }) when R_{1} = 0.1E_{0}. The first part (the linear part) of the graphs correspond to the hydrostatic and elastic stages. 
Fig. 5 Influence of the exponent m (left) and n (right) on the graph of the axial stress vs. the axial strain in a uniaxial compression test (p _{0} = 0). The values of m and n are indicated in red on the graph. In the left figure where m varies, one takes n = 1/2, whereas on the right figure where n varies, one takes m = 2. In all cases R_{1} = 0.1E_{0}. 
Fig. 6 Dependence of the overstress Δσ _{ z } on the exponents m and n of the model. Left: dependence on m when n = 1/2. Right: dependence on n when m = 2. In all cases R_{1} = 0.1E_{0} and k = 0.2. 
References
 E. Hoek, E.T. Brown, The Hoek–Brown failure criterion: a 1988 update, in: J. Curran (Ed.), Proceedings of the 15th Canadian Rock Mechanics Symposium, University of Toronto, Dept. of Civil Engineering, Toronto, Canada, 1988, pp. 31–38 [Google Scholar]
 E. Eberhardt, The HoekBrown failure criterion, Rock Mech. Rock Eng. 45 (2012) 981–988 [CrossRef] [Google Scholar]
 S. Cuvilliez, I. Djouadi, S. Raude, R. Fernandes, An elastoviscoplastic constitutive model for geomaterials: application to hydromechanical modelling of claystone response to drift excavation, Comput. Geotech. 85 (2017) 321–340 [CrossRef] [Google Scholar]
 J. Lemaître, Coupled elastoplasticity and damage constitutive equations, Comput. Methods Appl. Mech. Eng. 51 (1985), 31–49 [CrossRef] [Google Scholar]
 J. Lemaître, J.L. Chaboche, Mécanique des Matériaux Solides, Dunod, Paris, 1985 [Google Scholar]
 J. Mazars, A model of a unilateral elastic damageable material and its application to concrete, in: F.H. Wittmann (Ed.), Fracture toughness and fracture energy of concrete, Elsevier, New York, 1986 [Google Scholar]
 F. Ragueneau, C. La Borderie, J. Mazars, Damage model for concretelike materials coupling cracking and friction, contribution towards structural damping: first uniaxial applications, Mech. Cohesive Frict. Mater. 5 (2000) 607–625 [CrossRef] [Google Scholar]
 C. Comi, U. Perego, Fracture energy based bidissipative damage model for concrete, Int. J. Solids Struct. 38 (2001) 6427–6454 [CrossRef] [Google Scholar]
 C. Pontiroli, A. Rouquand, J. Mazars, Predicting concrete behaviour from quasistatic loading to hypervelocity impact: an overview of the PRM model, Eur. J. Environ. Civil Eng. 14 (2010) 703–727 [CrossRef] [Google Scholar]
 B. Richard, F. Ragueneau, C. Cremona, L. Adelaide, Isotropic continuum damage mechanics for concrete under cyclic loading: stiffness recovery, inelastic strains and frictional sliding, Eng. Fract. Mech. 77 (2010) 1203–1223 [CrossRef] [Google Scholar]
 J.F. Babadjian, G.A. Francfort, M.G. Mora, Quasistatic evolution in nonassociative plasticity: the cap model, SIAM J. Math. Anal. 44 (2012) 245–292 [CrossRef] [MathSciNet] [Google Scholar]
 S. Andrieux, Y. Bamberger, J.J. Marigo, Un modèle de matériau microfissuré pour les bétons et les roches, J. Mécanique Théorique Appliquée 5 (1986) 471–513 [Google Scholar]
 Q.Z. Zhu, J.F. Shao, D. Kondo, A micromechanicsbased thermodynamic formulation of isotropic damage with unilateral and friction effects, Eur. J. Mech. A 30 (2011) 316–325 [CrossRef] [Google Scholar]
 E. Lanoye, F. Cormery, D. Kondo, J.F. Shao, An isotropic unilateral damage model coupled with frictional sliding for quasibrittle materials, Mech. Res. Commun. 53 (2013) 31–35 [CrossRef] [Google Scholar]
 B. Halphen, Q.S. Nguyen, Sur les matériaux standards généralisés, J. Mécanique 14 (1975) 39–63 [Google Scholar]
 P. Germain, Q.S. Nguyen, P. Suquet, Continuum thermodynamics, J. Appl. Mech. 50 (1983), 1010–1020 [CrossRef] [Google Scholar]
 J.J. Marigo, Formulation d'une loi d'endommagement d'un matériau élastique, C. R. Acad. Sci. Sér II 292 (1981) 1309–1312 [Google Scholar]
 J.J. Marigo, Constitutive relations in plasticity, damage and fracture mechanics based on a work property, Nucl. Eng. Des. 114, (1989) 249–272 [CrossRef] [Google Scholar]
 A. Mielke, Evolution of rateindependent systems, in: Handbook of Differential Equations: Evolutionary Equations Vol. 2, Elsevier, Amsterdam, pp. 461–559 [Google Scholar]
 R. Alessi, J.J. Marigo, S. Vidoli, Gradient damage models coupled with plasticity: variational formulation and main properties, Mech. Mater. 80 (2015) 351–367 [CrossRef] [Google Scholar]
 G.A. Francfort, A. Giacomini, J.J. Marigo, The taming of plastic slips in von Mises elastoplasticity, Interfaces Free Boundaries 17 (2015) 497–516 [CrossRef] [Google Scholar]
 G.A. Francfort, A. Giacomini, J.J. Marigo, The elastoplastic exquisite corpse: a Suquet legacy, J. Mech. Phys. Solids 97 (2016) 125–139 [CrossRef] [Google Scholar]
 C. Comi, A nonlocal model with tension and compression damage mechanisms, Eur. J. Mech. A 20 (2001) 1–22 [CrossRef] [Google Scholar]
 K. Pham, J.J. Marigo, Approche variationnelle de l'endommagement: II. Les modèles à gradient, C. R. Mécanique 338 (2010), 199–206 [Google Scholar]
Cite this article as: J.J. Marigo, K. Kazymyrenko, A micromechanical inspired model for the coupled to damage elastoplastic behavior of geomaterials under compression, Mechanics & Industry 20, 105 (2019)
All Tables
All Figures
Fig. 1 Graphs of the axial stress σ _{ z } (left) and of the volumetric strain ϵ _{ v } (right) versus the axial strain ϵ _{ z } when R_{1} = 0.1E_{0} and p _{0} = 0 (uniaxial compression test). 

In the text 
Fig. 2 Influence of R_{1}/E_{0} on the graphs of the axial stress and the volumetric strain versus the axial strain when there is no confining pressure (uniaxial compression test). The value of R_{1}/E_{0} associated with each graph is indicated in red on the graph. 

In the text 
Fig. 3 When R_{1}/E_{0} is large enough, the graph of σ _{ z } vs. ϵ _{ z } contains a snapback (here R_{1} = 0.4E_{0}). Consequently, when the axial strain reaches the value corresponding to the limit point A, the response jumps from A to B. 

In the text 
Fig. 4 Graph of the axial stress vs. the axial strain for different values (in red on the graphs) of the confining pressure (in fact of the ratio p _{0}/σ _{ c }) when R_{1} = 0.1E_{0}. The first part (the linear part) of the graphs correspond to the hydrostatic and elastic stages. 

In the text 
Fig. 5 Influence of the exponent m (left) and n (right) on the graph of the axial stress vs. the axial strain in a uniaxial compression test (p _{0} = 0). The values of m and n are indicated in red on the graph. In the left figure where m varies, one takes n = 1/2, whereas on the right figure where n varies, one takes m = 2. In all cases R_{1} = 0.1E_{0}. 

In the text 
Fig. 6 Dependence of the overstress Δσ _{ z } on the exponents m and n of the model. Left: dependence on m when n = 1/2. Right: dependence on n when m = 2. In all cases R_{1} = 0.1E_{0} and k = 0.2. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.