Plasto-damage modelling for semi-brittle geomaterials

This paper presents an elastoplastic damage model for constitutive modelling of semi-brittle geomaterials showing two irreversible mechanisms. On one hand, the model deals with the plastic behaviour of a porous medium by a new variant of Barcelona Basic Model. On the other hand, the model combines the micromechanical definition of damage and phenomenological concepts in the framework of Continuum Damage Mechanics (CDM) for damage modelling. A second order tensorial damage variable is adopted for the model. Damaged effective stress variables are employed for formulation of elastoplastic behaviour laws and the plastic yield surface is a damage dependent one. The model has been validated by comparing the numerical results with experimental results of argillites.


Introduction
In the framework of radioactive waste disposal, argillaceous formations are under widespread consideration as host rock because of their favourable properties.Accordingly, constitutive modelling of clay formations as geological barrier is also an active topic for the researchers of this field.
Most of the currently existing geotechnical plastodamage models, mainly for the sake of simplicity, have made use of isotropic definition of damage (e.g.[1,2]).Contrary to this assumption, the experimental results on clay stones show an important plastic flow being coupled with an induced anisotropic damage [3,4].As a result, an elastoplastic damage model might be an appropriate choice for constitutive modelling of semibrittle geomaterials like argillites.In view of this, an anisotropic damage mechanism is formulated for the current study using a second order damage variable.
Although a number of plasto-damage models are also introduced in the literature, most of these models are based on the use of Bishop's like effective stress (e.g.[5,6]).
As pointed out by Sheng et al. [7], the complex stress variables (e.g.Bishop's effective stress) are dependent on material states and therefore are not straightforwardly controllable in conventional laboratory testing procedures.In the models which make use of the Bishop's effective stress with the parameter χ=S r , the inherent hydro-mechanical coupling is greatly dependent on the water retention relation used in the model.This entails utilizing sophisticated water retention models.Consequently the water retention properties of unsaturated medium should be specifically determined and it may necessitate the use of hysteretic water retention models in addition to mechanical formulation.
On the other hand the experimental results on the retention properties of damaged geomaterials are scarce in the literature.Therefore independent stress state variables including net stress and suction are adopted as stress variables of the proposed elastoplastic damage model.
This paper first presents the damage-related formulation of the model in the framework of Continuum Damage Mechanics (CDM).Afterwards, elastoplastic formulation is proposed based on Barcelona Basic Model concepts incorporating damage dependency of yield surface.At last, mechanical aspects of the model are validated by simulation of some triaxial tests and the results of the proposed model are compared with the experimental ones.

Formulation of the damage 2.1 Definition of an anisotropic damage variable
The damage variable is typically defined as a macroscopic measure of the microscopic degradation of a Representative Elementary Volume [8,9].To date, different kinds of damage variables have been used in the context of Continuum Damage Mechanics to represent the density and orientation of micro-cracks.One or more scalar damage variables for isotropic damage modelling and vectorial damage variables, second and fourth rank damage tensors or even eighth order tensors [10] for anisotropic damage models.Induced damage in brittle geomaterials is strongly anisotropic in nature.Only second or higher even-order tensors can be adopted as damage variable to characterize the damage induced anisotropy [11].As pointed out by [12], owing to the difficulties which would be encountered in postulating appropriate evolution laws for higher order tensors, the second-order tensor is preferred in the modelling of damage induced anisotropy.Having considered the above description, the second order variable  is adopted as the damage variable at the mesoscopic scale: This definition is based on spectral decomposition.From the physical point of view the approximately parallel damages developed in the Representative Elementary Volume (REV) are considered as three orthogonal equivalent meso-cracks.As it can be illustrated for a two dimensional case, each bunch of micro-defects with approximately parallel direction is homogenized in an equivalent meso-scale crack.The equivalent meso-cracks are characterized by unit normal vector of crack plane, n j , and a crack density parameter d j .

Damaged stress variables
As the meso-cracks forms within the medium, the material starts losing effective load carrying surface.As a consequence, stress increases in the undamaged part of the material.Mechanically, the new stress field of the damaged medium can be obtained by removal of the defects and replacement of nominal stresses by effective stresses.For the case of isotropic damage (that is, the scalar damage variable d), one can write the following in multiaxial stress space: However, in the case of the anisotropic damage, a fourth-order operator, M should be used for the definition of the damaged effective stress.
This is appropriate for a single phase damaged media, but for an unsaturated porous media both stress variables must be transformed to the effective state.Several mathematical forms have been proposed for the operator in the CDM context; however the operators which ensure symmetrisation of the damaged effective stress are preferred [13].Therefore, the operator of Cordebois and Sidoroff [14] is extended to the unsaturated state variables for the definition of the damaged net stress, σ ˆ  , and the damaged suction, s ˆ.
where δ is Kronecker delta.Having obtained the damaged stress variables (suction and net stress), the damaged medium behaviour can be formulated in the effective stress and actual strain space.This treatment can markedly simplify the computational implementation.

Damage criterion and damage evolution law
According to the experimental results reported in the literature, there is an obvious correlation between total tensile strains and the induced damage in brittle geomaterials.We can, therefore, assume that the damage growth is associated with the development of generalized tensile strain and the damage start developing within the medium as the strains reach a threshold which in itself is dependent on the current damaged state.Thus the Eq. 5 would be used as damaging yield function or damage criterion.It should be noted that this equation is inspired from [15].
C 0 , C 1 and are material parameters; wherein C 0 designates the initial damage threshold, and C 1 characterizes the degradation of the medium materials with the cumulated damage.This means that under the multiaxial stress states, the REV still continues damaging in tensile mode when the combination of major tensile strains attain damage driving threshold In addition, g characterizes the resistance of the material to cracking of the medium and propagation of damage.It should be noted that the term  ε g is part of damageconjugated stress or 'damage affinity' that the product of it by damage is an energy release rate.

E 2016 -
The sign convention of solid mechanics (compression negative) is adopted here.Accordingly, the positive part of total strain, ε + , represents the extensional components of deformations of the medium.Thus the spectral decomposition of the total strains can be used as follows: Where k  is k th eigenvalue of the strain tensor which corresponds to k th principal strain, k W is eigenvector (principal vector) associated to is the Heaviside function of the principal strain, H(x) = 1 for x ≥ 0, H(x) = 0 for x < 0.
Assuming an associated flow rule (or normal dissipation scheme) for damage, one can write the damage evolution law as follows: ) .( 2) , ( This treatment can guarantee that, the evolution rate of damage variable and the tensor of tensile strains are developing in the same direction, as it was desirable from the beginning.The evolution of positive scalar damage multiplier, d   , can be determined by the consistency condition of the damage criterion The loading-unloading conditions known as so-called Kuhn-Tucker conditions should be satisfied as follows: 3 Characterization of elasto-plastic behaviour

Degradation of elastic properties
Degradation of elastic properties is one of main features of mechanical behaviour of argillites which has often been attributed to the induced damage by micro-cracks.Different forms of relations are introduced in the literature to reproduce such effect.Equation ( 10) is employed here to modify the bulk modulus and shear modulus of undamaged material undergoing damage propagation.Non-interacting micro-cracks and moderate damage assumption are made to obtain these relations [16].

Formulation of plastic behaviour
Barcelona Basic Model (BBM) [17] and its variants are one of well known categories for constitutive modelling of clayey media.In this research, another elasto-plastic model based on the extension of BBM to damage effects is presented.The yield surface is modified so as to incorporate damage effects and new model is referred to as Damage Barcelona Basic Model (DBBM).It is assumed that the change of the yield function can be reproduced by variation of maximum apparent preconsolidation stress 0 p and cohesion related variable s p ˆ. Thus one can express yield function as a function of damaged suction, accumulated volumetric thermoplastic strains p v  as well as damage itself: In which: c p ˆis the apparent preconsolidation stress for the saturated state and 0 p is its counterpart for the unsaturated state.It is evident that the yield function is written based on the damaged stress components p q , ˆ.In fact, the damaged effective stress variables play the role of coupling variables for the current elastoplastic damage model which can readily couple the damage influence with the elastic criterion and the plastic formulation.The   The slope of the normal compression line (NCL) or virgin compressibility index can be defined by following relation [17]: Where r is a parameter controlling the maximum stiffness of the porous medium for an infinite suction, , and β is a parameter which characterizes the rate of increase of soil stiffness with suction.
With the assumption of small strains, the total strain rate can be decomposed into an elastic part and a plastic part as follows: Plastic strain increments are deduced from an associated flow-rule (F p =G p ): In which the rate of plastic multiplier p d is derived from the consistency condition given by: This leads to: Substitution of the plastic multiplier into flow rule yields the following expression for incremental constitutive equation: The terms of the equation are given in Appendix.

Numerical modelling and verification
To verify the adequacy of the model, triaxial testes have been modelled hereafter.The tests are conducted on a clay stone from the East region of France called 'Argilites de l'Est' [3,26,27].The specimens were 40 mm in height and 20 mm in radius and the tests were carried out in isothermal state.In order to validate model for different loading paths, triaxial lateral extension as well as triaxial compression test with unloading -reloading cycles are simulated.Having considered the geometry and the loading condition of the tests, an axisymmetric configuration is adopted for the model.
Preliminary parameters for elastoplastic and damage behaviour were mainly taken from argillite parameters introduced in the literature (e.g.[25][26][27]).First, a triaxial lateral extension test was back-analyzed to adjust the preliminary parameters so as to obtain parameters that yielded the best fit to the experimental result.The Final model parameters obtained are summarized in Table 1 and the result of the back analysis is shown in Fig. 3.  To validate the proposed model in complex loading conditions, a triaxial compression tests with unloadingreloading cycles from Chiarelli's thesis is simulated.Comparison of numerical results and experimental data are illustrated in Fig. 4. As can be seen, the model provided a reasonable result for the test with unloading cycles and it could reproduce degradation of elastic properties properly.Using the determined parameters, the authors have also simulated some triaxial compression tests with the proposed plastic damage model [28] and satisfactory results have been obtained although, for the sake of brevity of the text, these results are not included here.

Conclusions
This paper has dealt with the formulation of a poroplasto-damage model for semi-brittle geomaterials.A second-order tensor is adopted as damage variable, characterizing crack density in principal directions.Damage is assumed to grow with tensile strains of medium and constitutive equations are formulated in terms of effective stress and actual strain.To this end, damaged effective terms of stress state variables including the net stress and suction are defined according to the Continuum Damage Mechanics principles.
A new elastoplastic formulation is presented based on Barcelona Basic Model (BBM) concepts in which damage dependency of the yield surface and the potential function have been taken into account by introducing new relationships.The proposed formulation is implemented in Θ-STOCK Code being a finite element code for modelling of coupled problems of multiphase media.
To investigate capability of the new model, it is verified against the experimental results from the literature.In order to validate model for different loading paths, triaxial lateral extension as well as triaxial compression tests with unloading -reloading cycles are simulated.Verification results indicated that the developed poro-plasto-damage model can provide an adequate framework for the modelling of the behaviour of the semi-brittle materials like argillite.

Fig. 1 (
Fig. 1 (left) illustrates a cylindrical damaged Representative Elementary Volume (REV) submitted to a uniaxial tensile stress, σ.Fig. 1 (right) shows the undamaged counterpart of REV as a fictitious equivalent state in which the cross-sectional area of the medium is decreased due to the removal of the cracks and the nominal stress is replaced by the intensified stress, that is to say uniaxial effective damaged stress d   1 ˆ  .

Figure 1 .
Figure 1.Representation of the effective damaged stress.

2 
are material parameters characterizing the damage effects on the elastic properties of the geomaterial.

sp
ˆ represents the variation of cohesion with the corresponding suction or damage changes.k, ξ, γ H are material parameters through which the above mentioned variables are defined.ref p is the reference stress for the plastic yield surface. term the undamaged variables p s and p c to reflect variation of plastic yield locus with damage.

Fig. 2 Figure 2 .
Fig.2illustrates the yield loci of DBBM for the undamaged and damaged states.As can be seen, damage can shrink or swell the elastic domain.Adoption of positive value for ξ can enlarge the elastic domain while a negative value of ξ has an opposite effect resulting in shrinkage of the extent of the yield locus with the increase of damage intensity.

Figure 3 .
Figure 3.The results of back-analysis of triaxial lateral extension test.

Figure 4 .
Figure 4. Comparison of model results and experimental data for unloadingreloading cycles.
Matrices of incremental constitutive equation are as follows:

Table 1 .
Parameters of the Model.