A Coupled Plastic Damage Model for Concrete considering the Effect of Damage on Plastic Flow

A coupled plastic damagemodel with two damage scalars is proposed to describe the nonlinear features of concrete.The constitutive formulations are developed by assuming that damage can be represented effectively in the material compliance tensor. Damage evolution law and plastic damage coupling are described using the framework of irreversible thermodynamics.The plasticity part is developed without using the effective stress concept. A plastic yield function based on the true stress is adopted with two hardening functions, one for tensile loading history and the other for compressive loading history. To couple the damage to the plasticity, the damage parameters are introduced into the plastic yield function by considering a reduction of the plastic hardening rate. The specific reduction factor is then deduced from the compliance tensor of the damaged material. Finally, the proposed model is applied to plain concrete. Comparison between the experimental data and the numerical simulations shows that the proposed model is able to describe the main features of the mechanical performances observed in concrete material under uniaxial, biaxial, and cyclic loadings.


Introduction
The mechanical behavior of concrete is unique, due to the influence of micromechanisms involved in the nucleation and growth of microcracks and plastic flow.This behavior is characterized by several features as follows: the tensile strength, which is different from the compressive strength; the irreversible plastic deformation, which is found as the pressure on the concrete exceeds a certain value; beyond the peak stress, the stress-strain curve displays an unstable region, accompanied by stiffness degradation and strength softening.These features have to be considered in constitutive models of concrete materials.To model these features, several mechanics theories have been used.In general, the damage mechanics theories can be used to model the nucleation and growth of microcracks [1][2][3], whereas the plasticity theories can be used to model the plastic flow component of the deformation [4].
Because damage mechanics can be used to describe the progressive weakening of solids due to the development of microcracks, and because plasticity theories have been successfully applied to the modeling of slip-type processes, several plastic damage models of concrete have been developed through combining damage mechanics and plasticity theories [5][6][7][8][9][10][11][12][13].These coupled plastic damage models (CPDMs) could be formulated in the irreversible thermodynamics framework and can be easily applied to describe the essential nonlinear performances of concrete including the strain softening and the stiffness degradation.One type of CPDMs is based on the concept of effective stress, which was initially proposed by Kachanov [14] for metal creep failure.In these models, the plastic yield function is defined in the effective configuration pertaining to the stresses in the undamaged material [5,[7][8][9].Many authors used this approach to couple damage to plasticity.Some of these models, developed by the use of two damage scalars and damage energy release rate-based damage criteria, show excellent performance in reproducing the typical nonlinear behavior of concrete materials under different monotonic and cyclic load conditions.The other type of CPDMs is opposite of the above.In these models, the true stress appears in the plastic process, which clearly couples plasticity to damage [6,10,11,[15][16][17].However, considering the influence of damage on plastic flow, it is difficult to develop a plastic yield function that can be used to describe the plastic deformation and damage growth of concrete simultaneously.Taqieddin et al. [17] suggested that a damage-sensitive plastic yield function should be introduced in terms of true stress configuration.Even if several types of expressions for the plastic yield function written in terms of the effective stress have been successfully applied to model some of the typical nonlinearities of concrete (such as the volumetric dilation and strength increase under multidimensional compression), they cannot be directly used in the true stress space.A common approach to solve this problem is to introduce a reduction factor of plastic hardening rate into the plastic yield function.The plastic yield function is usually expressed by a function of the stress tensor and plastic hardening function, so the damage parameters are included in the plastic yield function with the introduction of the reduction factor in the plastic hardening function.Several authors applied this approach to develop the CPDM with a true stress space plastic yield function.Shao et al. [15] proposed a single scalar damage framework based on the elastoplastic damage release rate.The influence of damage on the plastic flow is calculated by considering a reduction of the plastic hardening rate.The reduction factor is introduced into the plastic yield function for the coupling between the damage evolution and the plastic flow.Nguyen and Houlsby [10] proposed a double scalar damage-plasticity model for concrete based on thermodynamic principles.In this model, the plasticity part is based on the true stress using a yield function with two hardening functions, one for the tensile loading history and the other for the compressive loading history.Two reduction factors are introduced into the tensile and the compressive hardening functions to consider the influence of the tensile and compressive damage mechanisms on the plastic flow, respectively.Taqieddin et al. [17] proposed a double scalar damage-plasticity model for concrete based on the hypothesis of strain energy equivalence.The model incorporates a plastic yield function written in terms of the true stress.Damage variables are introduced all over the plastic yield function.
For the models mentioned above, the reduction factor of the plastic hardening rate is equal to the damage variable defined in these models, and there is no distinction in the aspect of choosing a reduction factor between uniaxial and biaxial loadings.The plastic hardening rate is generally defined by the equivalent plastic strain, so the sum of the principal damage values (in three directions: , , and ) seems reasonable to use as the reduction factor of the plastic hardening rate under multiaxial stresses.By introducing such a reduction factor, the widely used plastic yield functions, such as those applied by Lee and Fenves [7], Wu et al. [9], and Lubliner et al. [18], can be used in the true stress space.
In the present work, a coupled plastic damage model is formulated in the framework of thermodynamics.The constitutive formulations are developed by considering an increment in the concrete compliance due to microcrack propagation.The plasticity part is based on the true stress using a yield function with tensile and compressive hardening functions.The plasticity yield function widely used in effective stress space is modified to be applied in this study by considering a reduction in the plastic hardening rate.Finally, the proposed model is applied to plain concrete.Comparison between the experimental data and the numerical simulations shows that the proposed model is able to represent the main features of mechanical performances observed in concrete material under uniaxial, biaxial, and cyclic loadings.

General Framework of Coupled Plastic Damage Model
The model presented in this work is thermodynamically consistent and is developed using internal variables to represent the material damage state.The assumption of small strains will be adopted in this work.In the isothermal conditions, the state variables are composed of the total strain tensor , scalar damage variables  + and  − , plastic strain   , and internal variables for plastic hardening  + and  − .The total strain tensor is decomposed into an elastic part   and plastic part   : To establish the constitutive law, a thermodynamic potential for the damaged elastoplastic material should be introduced as a function of the internal state variables.Considering that the damage process is coupled with plastic deformation and plastic hardening, the thermodynamic potential can be expressed by  =   ( ± ,   ) +   ( ± ,  ± ) , where   and   are the elastic and plastic energy for plastic hardening of damaged material, respectively.S is the fourthorder stiffness tensor of the damaged material.
According to the second principle of thermodynamics, any arbitrary irreversible process satisfies the Clausius-Duhem inequality as Substituting the time derivative of the thermodynamic potential into the inequality yields Because the inequality must hold for any value of ε  , ε  , ḋ+ , ḋ − , κ + , and κ − , the constitutive equality can be obtained as follows: The constitutive relation between the stress and strain tensors can also be expressed utilizing the material's compliance tensor as where C is the fourth-order compliance tensor of the damaged material, and (⋅) −1 is the matrix inverse function.It is assumed in (6) that the damage can be represented effectively in the material compliance tensor.This assumption is in line with many published papers.Budiansky and O' connell [19] and Horii and Nemat-Nasser [20] indicated that damage induces degradation of the elastic properties by affecting the compliance tensor of the material.Based on this consideration, several researchers proposed an elastic degradation model in which the material stiffness or compliance was adopted as the damage variable.The elastic degradation theory was then extended to the multisurface elastic-damage model [21] and to the combined plastic damage model [6,15,16,22].
Considering that an added compliance tensor is induced by the microcrack propagation, the fourth-order compliance tensor C is decomposed as follows [23]: where C 0 defines the compliance tensor of the undamaged material and C  defines the total added compliance tensor due to microcracks.Microcracks in concrete are induced in two modes: the splitting mode and compressive mode [6,23].In general, the splitting mode is dominant in tension, whereas the compressive mode is dominant in compression.
To incorporate these modes into the formulation, the stress tensor is decomposed into a positive part  + and negative part  − : where the eigenvalues of  + are nonnegative, and the eigenvalues of  − are all nonpositive.According to Faria et al. [8], Wu et al. [9], and Pelà et al. [24], each part of the stress tensor  can be expressed as where the fourth-order projection tensors P + and P − are expressed as in which I is the fourth-order identity tensor, σ and p  are the th eigenvalue and the normalized eigenvector of the stress tensor , respectively, and (⋅) is the Heaviside function (() = 0 for  < 0 and () = 1 for  ≥ 0).By introducing the positive/negative projection operators P ± , the total added compliance tensor C  is decomposed as [25] C  = P + : C +  : where C +  and C −  represent the added compliance tensors due to microcracks induced in tension and compression, respectively.
To progress further, evolutionary relations are required for the added compliance tensors C +  and C −  .Based on the previous work of Yazdani and Karnawat [6], Ortiz [23], and Wu and Xu [25], the added compliance tensors are expressed in terms of response tensors M + and M − such that where the response tensors M + and M − determine the evolution directions of the added compliance tensors C +  and C −  , respectively.The thermodynamic forces conjugated to the corresponding damage variables are given by where the terms S/ + and S/ − are the derivatives of the stiffness tensor of the damaged material with respect to the damage variables  + and  − , respectively.These terms can be calculated utilizing the material's compliance tensor as Substituting ( 14) into (13) and calling for (5), one obtains Finally, the rate form of the constitutive equation can be given as The derivation of the detailed rate equation from ( 16) requires determining the evolution laws of the damage variables and plastic strains.

Damage Characterization. The damage criteria are defined as [15]
where  + / − (tension/compression) represents the threshold of the damage energy release rate at a given value of damage.
On the basis of the normality structure in the continuum mechanics, the evolution law for the damage variables can be derived by in which λ +  ≥ 0 and λ −  ≥ 0 denote the tensile and compressive damage multipliers, respectively.Under the loadingunloading condition, λ +  and λ −  are expressed according to the Kuhn-Tucker relations: In the special case of elastic damage loading without plastic flow ( λ  = 0), the damage multipliers λ +  and λ −  can be determined by calling for the damage consistency condition under static loading:

Plastic Characterization.
Irreversible plastic deformation will be formed during the deformation process of concrete.The plastic strain rate can be determined using a plastic yield function with multiple isotropic hardening criteria and a plasticity flow rule.The yield function determines under what conditions the concrete begins to yield and how the yielding of the material evolves as the irreversible deformation accumulates [26].According to Shao et al. [15] and Salari et al. [27], the plastic yield criterion   can be expressed by a function of the stress tensor  and the thermodynamic forces  + ( + ,  + ) and  − ( − ,  − ), respectively, conjugated with the internal hardening variables  + and  − : The plastic hardening functions  + and  − can be deduced by the derivative of the thermodynamic potential: The plastic potential   (,  ± ,  ± ) is also a function of the stress tensor, the scalar damage variables, and the internal hardening variables.Note that   is the yield function in the associated flow rule.At the given function of the plastic potential, the evolution law of the plastic strain is expressed as follows: (,  ± ) = 0, λ ≥ 0,   (,  ± ,  ± ) λ = 0. (24) In the special case of plastic loading without damage evolution ( λ +  = 0 and λ −  = 0), the plastic multiplier λ can be determined from the plastic consistency condition: 2.3.Computational Procedure.The numerical algorithm of the proposed constitutive model is implemented in a finite element code.The elastic predictor-damage and plastic corrector integrating algorithm proposed by Crisfield [28] and then used by Nguyen and Houlsby [10] is used here for calculating the coupled plastic damage model.In this algorithm, the damage and plastic corrector is along the normal at the elastic trial point, which avoids considering the intersection between the predicting increments of elastic stress and the damage surfaces.Furthermore, this algorithm ensures that the plastic and damage consistent conditions are fulfilled at any stage of the loading process.For details on the numerical scheme and the associated algorithmic steps, refer to the published reports [10,28].
The following notational convention is used in this section: (1)  trial  represents  at the trial point in time step .
The elastic trial stress is defined as follows: By solving ( 16), (17), and ( 21) in terms of the trial stress, the increments of the equivalent plastic strains Δ + and Δ − , plastic strain Δ  , and damage variables Δ + and Δ − can be obtained: The expression for the increments of the plastic strain Δ  can be obtained by substituting (27b) into (23): Substituting ( 28) into (27a) yields the relationship between the stress and the damage variables, written in incremental form as follows: where the fourth-order tensor C st is the constitutive matrix, which is the tangent with respect to the plasticity and the secant with respect to damage Substituting ( 29) into (27c) and (27d) yields the increments of the damage variables Δ + and Δ − .Back substituting the increments of the damage variables into (29) results in the stress increment Δ.

Determination of Specific Functions For Concrete
Specific functions for concrete are proposed based on the general framework of the coupled plastic damage model given in the previous section.
Referring to the definition of the isotropic compliance tensor, the evolution directions of added compliance tensors can be defined as [25] where  0 is the initial Young's modulus, and V 0 is the initial Poisson's ratio.The symmetrized outer product "⊗" is defined as Generally, the parameters  + and  − are the functions of the damage scalars, respectively.Their mathematical expressions can be logarithmic, power, or exponential functions.The specific expressions of  + and  − are given as where  + 0 and  − 0 are, respectively, the initial damage energy release threshold under tension and compression,  + and  + are the parameters controlling the damage evolution rate under tension, and  − and  − are the parameters controlling the damage evolution rate under compression.
The plastic yield function adopted in this work was first introduced in the Barcelona model by Lubliner et al. [18] and was later modified by Lee and Fenves [7], Wu et al. [9], and Voyiadjis et al. [11].It has proven to be excellent in modeling the biaxial strength of concrete.However, the yield function is usually used in effective stress.To apply this function in the CPDM written in terms of true stress, Taqieddin et al. [17] modified the function by introducing damage parameters and achieved better results.Taqieddin et al. showed that using the same model and material parameters, fully coupled plastic damage yield criterion showed a more pronounced effect of damage on the mechanical behavior of the material than did the weakly coupled plastic damage model yield criterion.Considering that damage induces a reduction in the plastic hardening rate, the damage parameters are then introduced into the yield function by the plastic hardening functions Mathematical Problems in Engineering of damaged material [ + ( + ,  + ) and  − ( − ,  − )].The yield function used in the present work is expressed as follows: where  1 =  11 +  22 +  33 is the first invariant of the stress tensor . 2 =  : /2 is the second invariant of the nominal deviatoric stress   =   −     /3.The parameter  is a dimensionless constant given by Lubliner et al. [18] as follows: where  − 0 and  − 0 are the initial equibiaxial and uniaxial compressive yield stresses, respectively.Based on the experimental results, the ratio  − 0 / − 0 lies between 1.10 and 1.20; then, as derived by Wu et al. [9], the parameter  is between 0.08 and 0.14.In the present work,  is chosen as 0.12.The parameter  is a parameter expressed from the tensile and compressive plastic hardening functions: Because the associative flow rule is adopted in the present model, the plastic yield function   is also used as the plastic potential   to obtain the plastic strain.
According to (22), a specific expression of the plastic energy   should be given to obtain the plastic hardening functions.Considering, as previously, a reduction in the plastic hardening rate due to damage, the plastic energy is expressed as follows: where the reduction factors  ±  are functions of the damage scalars  ± , and   0 is the plastic hardening energy for undamaged material, which is defined by in which  + 0 and  − 0 are the initial uniaxial tensile and compressive yield stresses, respectively.ℎ + and ℎ − are the tensile and compressive plastic hardening moduli, respectively.The hardening parameters  + and  − are the equivalent plastic strains under tension and compression, respectively, defined as [1,17] where κ + and κ − are the tensile and compressive equivalent plastic strain rates, respectively.
According to Wu et al. [9], the evolution laws for κ + may be postulated as follows: where ̂ε  max and ̂ε  min are the maximum and minimum eigenvalues of the plastic strain rate tensor ε  , respectively.σmax and σmin are the maximum and minimum eigenvalues of the stress tensor . is a weight factor, expressed as in which the symbol ⟨⋅⟩ is the Macaulay bracket, defined as ⟨⟩ = ( + ||)/2.Note that  is equal to one if all of the eigenstresses σ are positive, and it is equal to zero if all of the eigenstresses σ are negative.
The plastic tensile and compressive hardening functions of the damaged material are then specified as It can be observed from ( 33) and (41) that the damage variables are introduced into the plastic yield function.
The last expressions to be specified are the reduction factors  ±  .Defined as the sum of the principal damage values in three directions (, , and ), the reduction factor is obtained from the definition of a particular thermodynamic energy function in the present paper.Based on the previous works of Lemaitre and Desmorat [29], the following particular form of the elastic thermodynamic energy function is used: where   is the hydrostatic stress, and   =  1 + 2 + 3 , and  1 ,  2 , and  3 are the damage values in , , and  directions, respectively.The effective damage tensor H is given by Differentiation of the elastic thermodynamic energy yields the strain-stress relations where (⋅)  is the component of the deviatoric tensor.
Substituting (44) into (6), the relation between the second-order damage tensor and the compliance tensor of the damaged material can be given by By solving (45), the reduction factor can be obtained.For example, considering that concrete is subjected to biaxial tension loading ( 2 / 1 =  and  3 = 0), the reduction factor  +  can be obtained Similarly,  −  , corresponding to compression loading, can be calculated In the special case of uniaxial loadings or combined tension-compression loading ( = 0), the reduction factors can be expressed as

Numerical Examples and Discussion of Results
In this section, the procedure for the determination of the model parameters is first presented.Then, several numerical examples are provided to investigate the capability of the proposed model in capturing material behavior in both tension and compression under uniaxial and biaxial loadings.The examples are taken from [7], with corresponding experimental data provided by Kupfer et al. [30], Karson and Jirsa [31], Gopalaratnam and Shah [32], and Taylor [33].To comply with the results of the studies [7,9], these numerical examples were analyzed using the single quadrilateral finite element shown in Figure 1.

Identification of Model
Parameters.The proposed model contains 12 parameters: four elastic constants for the undamaged material ( 0 , V 0 ,  + 0 , and  − 0 ), two parameters for plasticity characterization (ℎ + and ℎ − ), and six parameters for damage characterization ( + 0 ,  − 0 ,  ± , and  ± ).All of the parameters can be identified from the uniaxial tension and compression tests.
The initial elastic constants ( 0 , V 0 ,  + 0 , and  − 0 ) are determined from the linear part of the typical stress-strain curve in Figure 2.
As mentioned before, the damage evolution is coupled with the plastic flow.According to Shao, it is reasonable to assume that the damage initiation occurs at the same time as the plastic deformation.They are identified as the end of the linear part of the stress-strain curves (point  in Figure 2).In this case, considering  ± = 0 and  ± = 0 in the damage criterion and the plastic yield function, the initial damage thresholds  + 0 and  − 0 can be determined as follows: where the stress ratio is expressed as  =  2 / 1 .
The two parameters ( + and  + ) characterizing the failure surface and the plastic hardening parameter (ℎ + ) can be determined by fitting both curves of stress-strain and stress-plastic strain obtained in the uniaxial tension test.Similarly, the parameters ( − ,  − , and ℎ − ) are determined by fitting both curves of stress-strain and stress-plastic strain obtained in the uniaxial compression test.
The effects of the model parameters ( + ,  + , and ℎ + ) on the stress-strain response, the stress-plastic strain response, and the damage evolution under uniaxial tension are shown in Figure 3.Each model parameter in turn is varied, while the others are kept fixed, to show the corresponding effect on the behavior of the concrete material.A good agreement exists between the experimental data [33] and the numerical simulations obtained with the associated material parameters given in Table 1.It can also be observed in Figure 3 that not only is the stress-strain curve governed by all of the parameters ( + ,  + , and ℎ + ) of the model but also the plastic strain and damage evolution are controlled by all of the parameters.This result is in agreement with that described by the proposed damage energy release rate  +  and the plastic yield criterion   in which both the damage scalar  + and the plastic hardening rate  + are introduced.These responses observed in Figure 3 show the coupled effect of damage and plasticity on the predicted response.
To illustrate the effects of the model parameters ( − ,  − , and ℎ − ) on the behavior of the proposed model, the stress-strain curve, the stress-plastic strain curve, and the damage evolution under uniaxial compression are plotted in Figure 4 in which the results obtained from varying each parameter while keeping the others fixed are also presented.The model parameters obtained from the experimental data [31] are listed in Table 1.It can be observed in Figure 4 that the numerical predictions, including not only the hardening and softening regimes of the stress-strain curve but also the stress-plastic strain curve, are in good agreement with the experimental data [31].It can also be observed that all of the responses of the model, including the stress-strain curve, the stress-plastic strain curve, and the damage evolution, are controlled by all of the parameters ( − ,  − , and ℎ − ).This finding corresponds to the expressions of the proposed damage energy release rate  −  and the plastic yield criterion   in which both the damage scalar  − and the plastic hardening rate  − are introduced.These responses observed in Figure 4 again indicate the coupled effect of damage and plasticity on the predicted behavior.

Biaxial Loading.
In this section, the behavior of the proposed model in combined loadings (biaxial compression and biaxial tension-compression) is investigated.The experimental data of Kupfer et al. [30] are used to compare with the numerical predictions of the model.The model parameters used in this series of tests are determined from the uniaxial test ( 2 / 1 = 0/−1).Then, these parameters are used to model the behavior of the concrete under different stress ratios ( 2 / 1 = −0.52/− 1,  2 / 1 = −1/−1, and  2 / 1 = 0.052/ − 1).The model parameters obtained for the concrete are listed in Table 2.Note that after being given the parameters  ± 0 and  ( =  2 / 1 ), the initial damage threshold  ± 0 can be calculated by (49).In addition, as the experimental data are available only in the prepeak regime under the uniaxial tension test, some of the model parameters can be assumed to yield good fit in biaxial tension-compression, such as ℎ + = 31000 MPa,  + = 0.6, and  + = 1.4.By using these parameters, simulations of biaxial compression-compression and compression-tension tests with different stress ratios have been performed.As Figure 5 shows, the numerical results obtained from the model are observed to be in good agreement with the experimental data.

Cyclic Uniaxial
Loading.The cyclic uniaxial tensile test of Taylor [33] and the cyclic compressive test of Karson and Jirsa [31] are also compared with the numerical results.The properties and model parameters for Taylors's and Karsan and Jirsa's simulation are listed in Table 1.As shown in Figure 6, the strength softening and stiffness degrading, as well as the irreversible strains upon unloading, can be clearly seen under both cyclic uniaxial tension and compression.This demonstrates the coupling between damage and plasticity, as well as the capability of the model in reproducing mechanical features of concrete.

Conclusions
A coupled plastic damage model written in terms of the true stress has been proposed in this paper to describe the nonlinear features of concrete in uniaxial and biaxial loadings.Based on the theoretical development of the model formulation, several conclusions have been summarized as follows.
The constitutive formulations are developed by considering an added flexibility due to microcrack growth.It is assumed that damage can be represented effectively in the material compliance tensor.The framework of irreversible thermodynamics is adopted to describe the damage evolution and plasticity damage coupling.
The plasticity part is based on the true stress using a yield function with two hardening variables, one for the tensile loading history and the other for the compressive loading history.To couple the damage to the plasticity, the damage parameters are introduced into the plastic yield function by considering a reduction in the plastic hardening rate.The specific reduction factor is defined as the sum of the principal values of a second-order damage tensor, which is deduced from the compliance tensor of the damaged material.
The model contains 12 parameters, which can be obtained by fitting both curves of the stress-strain and the stress-plastic strain in the uniaxial tension and compression tests.The numerical results suggest that the proposed model is able to describe the main features of the mechanical behavior for concrete under uniaxial, biaxial, and cyclic loadings.

Table 2 :
Model parameters for concrete material under biaxial loadings.

Figure 6 :
Figure 6: Comparison of the model predictions with the experimental results in (a) cyclic uniaxial tension and (b) cyclic uniaxial compression.