A macroscopic gradient-enhanced damage model for deformation behavior of concrete under cyclic loadings

Concrete materials show a very complex macroscopic deformation behavior under tension and compression, accompanied by crack opening and crack closing phenomena under cyclic loading. The continuum damage mechanics offers a promising framework for the description of the damage deformation behavior. This paper proposes a continuum damage model, which is formulated based on energy equivalence using a unified equivalent strain. The evolution of isotropic damage is governed by two independent history variables to describe the crack opening and closing behavior, i.e., unilateral behavior, of concrete. The evolution of damage and inelastic strains are described by a single damage function and a modified failure surface, respectively. Moreover, the implicit gradient method is applied to the equivalent strains to achieve proper localization of deformation. The stiffness recovery and crack opening/closing mechanisms are simulated considering the thermodynamically consistent framework. Validation of numerical results with experimental data and the previous models demonstrates the efficiency of the model to simulate concrete behavior under monotonic and cyclic/reverse loadings.


Introduction
Concrete remains as a broadly used construction material in the field of civil engineering. Due to its main characteristics of low tensile strength compared to compressive strength, it exhibits ductile behavior under compression and brittle behavior under tension. Thus, the brittleness of concrete causes the damage of the material in the form of opening of tensile cracks and consequently may lead to failure of structural components. In addition, the material exhibits a gradual decrease in stress with the increase in strains in the post-peak region, i.e., strain softening. Therefore, it is very essential to understand the mechanical behavior of concrete under various loadings such as monotonic, cyclic, and reverse loadings.
In the past, several researchers have developed various models in the framework of continuous media to describe the deformation behavior of concrete. Although the irreversible deformation of concrete is not fully understood, several models apply theories of plasticity to simulate concrete behavior [2,8,34,41]. On the other hand, models based on linear elastic fracture mechanics [3] yield unrealistic results. Unlike the ductile fracture of metals, the fracture process zone of concrete is inhomogeneous and governed by multiple microcracks. Therefore, a fictitious crack model or cohesive crack model (a line crack model), which is based on nonlinear fracture mechanics characterizing cracks by a stress-displacement rather than stress-strain relation, was developed by Hillerborg [17]. Later, it is found that the material loses stability, as the stiffness ceases to be positive definite once the strain softening begins. As a consequence, the strain softening tends to localize strains within a zone of vanishing volume [4]. The disadvantages of the line crack model that separates the crack process zone are avoided by smeared crack models or Crack Band Models [6], which are characterized by three parameters, namely the tensile strength, the fracture energy, and the cracking front width.
On the contrary, many continuum models have been developed based on the work of Kachanov [18,19] who first introduced a scalar internal variable to give a physical meaning of damage. Subsequently, several continuum damage models [9,25,27,28,36] are developed for concrete by introducing a set of internal state variables solely contributing to the post-peak nonlinear softening behavior. On the other hand, there are several other models, namely coupled plastic damage models [23,24,40,42] and multi-surface plastic damage models [7,11,33], that couple plasticity and either isotropic or anisotropic damage theories for concrete.
Most of the elastic damage or plastic damage models introduce two independent internal damage variables to represent the distinct behavior of concrete in tension as well as in compression. The evolution of these damage variables is based on a damage surface. In another approach, a total damage is a weighted sum of these damages in tension and compression related to two different loading surfaces, respectively [28]. Occasionally, the total damage is estimated from two damage variables using a multiplicative relation [11,23] to consider elastic stiffness recovery during loading/unloading processes. Likewise, the positive (tension) and negative (compression) parts of effective stresses, which are obtained using spectral decomposition, are sometimes linked to these damage variables [10,42]. However, an effective damage could also be considered as a function of two independent history deformation variables associated with two different loading surfaces [27]. Thereby, these history variables are related to damage equivalent strains in tension and compression, respectively. Furthermore, the history variable governing the evolution of damage is usually related to a scalar measure of local deformation via the loading criterion function. It is also noteworthy that there are various definitions available in the literature for the consideration of local measure of deformation [29,36], which can be determined by using damage equivalent strains in terms of either principal stresses [9] or principal strains [27,28], damage energy release rate [10,42], effective stress [26], equivalent effective stress [14] and inelastic strains [7] or sometimes inelastic multipliers [11]. Similarly, there are also different alternatives to express the damage equivalent strain, namely Mazars criterion [28], modified von Mises strain [13,36], and two separate expressions, for cracking and crushing, respectively [27]. Even though the equivalent strain is modified in order to account for predicting both the tensile and compressive behaviors, these modifications are not yet fully capturing the initial elastic domain and ultimate stress domain. Numerical results show considerable differences, especially in bicompression and complex regions, while comparing with experimental data [38]. Moreover, the numerical predictions of Mazars's unilateral model [27] are not smooth enough in the complex region, though this model offers good results with few differences near bisector of bicompression region.
Moreover, it is witnessed from the experimental studies on concrete [16,20,26,32] that concrete exhibits inelastic/permanent deformation to some extent under direct cyclic loading or reverse cyclic loading even though this inelastic deformation is not as large as that observed in ductile materials like steel. Thence, the incorporation of inelastic deformation behavior is unavoidable to model the deformation behavior of concrete correctly.
When these local damage models are implemented into a finite element program, the localization occurring in fracture process zone leads to the loss of well-posedness of boundary value problems and mesh objective results. Therefore, an internal length scale can be incorporated by an implicit gradient method of nonlocal enrichment to avoid these numerical difficulties, as the implicit gradient formulation remains as a one of the most efficient approaches in regularization. The nonlocal enrichment should be applied to the variable that drives the damage growth. Nevertheless, standard gradient damage models [12,15] based on strain equivalence concept may show broadening effect upon increase in loading. This kind of continued broadening is caused, as the local equivalent strains only depending on total strains go on increasing upon loading.
The main objective of this paper is to formulate a continuum damage model to describe crack opening/closing behavior exhibited by concrete using a unified equivalent strain based on energy equivalence principle. To address the important aspects of concrete behavior under tension and compression, the proposed model defines a single loading criterion for damage evolution. To capture inelastic strains, a slightly modified version of extended Lubliner-Lee failure criterion [22,23] is used. To take stiffness recovery/unilateral effects due to crack closure during cyclic loading into account, the model uses a history variable and inelastic multiplier as a function of two independent history variables and inelastic multipliers, respectively. To avoid ill-posedness of problem and achieve a proper localization, the present model applies implicit gradient method to equivalent strains, which act as the driving force for the damage growth. To validate the present model with experimental tests [16,20,26], some numerical examples using the model are analyzed under uniaxial tensile and compression and direct cyclic/reverse loading conditions. Finally, the localization phenomena are also illustrated using a numerical example.

Constitutive modeling of elastic damage
The stress-strain relation for elastic damage material behavior can be written using the principle of energy equivalence as follows: where σ and ε are the second-order Cauchy stress tensor and the linear strain tensor, respectively, and H is the fourth-order elasticity tensor of the undamaged material. The damage variable D is a scalar measure of effective isotropic damage. Inelastic strain evolution is not considered. Thus, the strain softening behavior of concrete is described by the constitutive law (1).

Damage evolution law
The damage variable D is generally considered to be driven by a history variable κ, which is described as a maximum deformation of the material occurred during tension or compression loading path. Therefore, κ is a monotonically increasing parameter during loading/unloading processes and the accumulation of damage. Hence, the damage evolution law [29] that exhibits exponential softening is explicitly expressed as a function of κ as follows: where κ 0 is an initial threshold to set the initial elastic domain below which there is no damage. β 1 and β 2 are model parameters characterizing the softening behavior of the material after the peak stress up to a complete loss of stiffness. β 1 controls the initial growth of damage, and β 2 ensures the later damage growth and its large extent. Thus, D is limited to the range 0 ≤ D ≤ 1 at any case.

Damage equivalent strain
In the framework of continuum damage mechanics, it is obvious that κ is often related to a scalar measure of deformation, a damage equivalent strain at a point of continuum considered locally. There are several alternatives for determining the local measure , as discussed in Sect. 1. In this work, is defined as a unified damage equivalent strain for describing tensile cracking as well as crushing failure, which is given by where I 1 and J 2 are the first invariant of elastically predicted stress tensor σ p and second invariant of the deviatoric part of σ p , respectively; σ max is the maximum principal stress of σ p ; and H is a Heaviside function, whereby H = 1 if σ max > 0, else H = 0. The parameters α L and β L are defined as follows [23,24]: where f bc and f c are the characteristic strength values in equibiaxial and uniaxial compression, respectively. The experiments confirm the ratio f bc / f c0 ranging from 1.10 to 1.16. For more details, the reader may refer to [23,24]. The expression (3) is inspired from the failure criterion of Lubliner [23,24,42] and slightly transformed herein in order to express the damage criterion in strain space, which will be discussed in the following section.

Damage surface criterion
As similar to a yield or failure surface of the material, a function is generalized here to limit the initial elastic domain and also to allow the growth of damage. The damage criterion f is conveniently used to relate κ to the unified damage equivalent strain . Therefore, the associated single loading surface is written as with κ as a maximum threshold value reached during tensile/compression loading path. Thus, it is defined by where κ 0 is the initial threshold value at which damage initiates. The evolution of the history variable κ can mathematically be expressed by the Kuhn-Tucker loading/unloading relations [9,13] as follows: where (˙) represents the derivative of any variable with respect to time t. The function f < 0 means that the material response remains linear elastic, as there is no damage development. The damage initiates when f = 0, since this is the merely dissipative mechanism considered in the model for the brittle failure of concrete in the absence of inelastic strains. During the damage evolution, the consistency conditionḟ = 0 must always be valid. Moreover, the damage remains constant ifκ = 0.

Tension-compression behavior
Two independent history variables such as κ t and κ c can be defined for tension and compression, respectively, as follows: Subsequently, the expression for κ (6) can be expressed in terms of κ t and κ c using the Heaviside function H as to account for the crack opening/closure effects (i.e., unilateral behavior). Consequently, the model parameters can also be determined from the respective parameters under tension as well as compression Remarks In the case of the model including crack opening/closure behavior, one needs to understand correctly how the evolution of the internal variables occurs. Figure 1 shows the evolution of internal variables and strain measures during tension-compression loading/unloading and reverse loading processes. As can be seen, in the first tension phase, ε increases when t < t 2 and decreases when t 3 > t > t 2 depending on whether loading or unloading occurs. But the adopted expression (3) of yields always positive during the entire tensioncompression loading processes, but increases during loading (t < t 2 or t 3 < t < t 4 ) and decreases during unloading (t 3 > t > t 2 or t 5 < t < t 6 ). It is clearly shown that κ t or κ c , which describes the maximum threshold value that the material experienced during its tensile or compression load history, respectively, always shows monotonic increase. Therefore, either the history variable κ or the effective damage D refers to the active stress state, because the effective damage D becomes zero if the subsequent loading occurs in compression direction (t 3 < t < t 4 ). By contrast, if the loading in tensile direction (t 6 < t) occurs, the effective damage D recovers from the previous tensile loading, as the cracks open up. Thus, bothκ andḊ remain positive for all positive or negative strain states.

Coupling of inelasticity
The Lubliner-Lee [22][23][24] yield criterion is found to be successful in modeling the nonlinear behavior of concrete under both monotonic and cyclic loads. This model uses a loading surface that couples plasticity and isotropic damage by effective plastic strain. The evolution of the isotropic damage variable is described by means of two damage variables in multiplicative manner in order to take the crack opening/closing effects into account. However, the following section describes a slightly modified version of this criterion considering the evolution of inelastic strains that is coupled with the unilateral description of isotropic damage. Isotropic linear hardening is considered in tension as well as compression [42].

Inelastic deformation
In the theory of inelasticity under small deformations, the total strain tensor ε is decomposed into the elastic part of strain ε el and the inelastic part of strain ε in such that: Subsequently, by taking time derivative of Eq. (12), the total strain rateε can be written aṡ whereε el andε in refer to the rates of the elastic and inelastic strain tensors, respectively.

Constitutive law
After coupling of the inelastic strains into the isotropic damage model based on the energy equivalence, the generalized Hooke's law (1) becomes as where H is the fourth-order undamaged elastic stiffness tensor. Here, σ , ε, ε el and ε in are considered as macroscopic homogeneous variables. Thereafter, the elastic strain ε el follows Subsequently, the effective stress and effective strain tensors of the damaged material can be written as In the case of inelastic consideration, the same damage evolution law (2) models the proper softening behavior of concrete under tension, but yields an underestimated hardening-softening behavior under compression. It means that a lower peak stress value from which the softening begins is predicted due to the incorporation of inelastic deformation. Thus, the damage evolution law (2) with the existing model parameters is not sufficient to model the behavior of concrete material perfectly. As a result, it entails a slight modification in the damage evolution so as to capture an appropriate behavior in compression as well. Thus, the softening law (2) is further modified with an additional model parameter β 3 as follows: Note that the influence of β 3 on the damage evolution or on the related stress behavior is considered only for κ > κ 0 within this modification. Refer to Sect. 3.5.

Inelastic failure surface
It is familiar that concrete exhibits distinct behavior in tension and compression. The failure criterion initially developed by Lubliner et al. [24] and later extended by Lee and Fenves [23] is adopted here. The extended version of the model, which was based on effective stress concept or strain equivalence concept, was expressed in effective stress space. As we use the energy equivalence principle of damage description in this work, both equivalent stress and limiting strength need to be expressed in the effective stress space to recover the original model behavior. Accordingly, a slight modification is embodied in this work by including a reduction of (1− D) to the limiting strength as follows: whereĨ 1 andJ 2 are the first invariant of the effective stressσ and the second invariant of the deviatoric part of the effective stressσ , respectively, H is the Heaviside step function (H = 1 forσ max > 0 and H = 0 for σ max < 0), andσ max is the maximum principal stress. The dimensionless constants α and β are defined as follows [24]: where f bc and f c0 are the initial flow stresses in equibiaxial and uniaxial compression, respectively. The experiments confirm the ratio f bc / f c0 ranging from 1.10 to 1. 16.
In order to consider isotropic hardening, the hardening function is adopted in tension as well as compression by assuming linear functions as below: The evolution functions of tensile and compressive hardening can be defined aṡ where H t and H c are the hardening modulus in tension and compression, respectively, and f t0 is the initial flow stress in uniaxial tension. The hardening variables h + in tension and h − in compression are defined in rate form [23,24] as follows: where the dimensionless parameter w, which is a triaxial weight factor obtained from the effective stress σ = 0, is defined according to [23] as where x = x+|x| 2 refers to Macaulay bracket and principal directions k = 1, 2, 3.ε in max =ε in 1 andε in min =ε in 3 are the maximum and minimum inelastic strain rates of the inelastic strain rate tensorε in such that ifε in 1 > ε in 2 >ε in 3 . Thus, Eq. (22) can be equivalently written in tensor form as ⎛ ⎜ ⎝ḣ where w + = w and w − = −(1 − w).

Non-associative flow rule
In order to obtain the volumetric expansion or dilatancy of frictional materials like concrete properly under given loading, the definition of a non-associative flow rule is important. As soon as the failure surface criterion F is satisfied, the material starts to deform inelastically upon further increase in the loading. The inelastic deformation is measured by the evolution of inelastic strains ε in , which is evaluated by the flow rule as follows:ε whereλ is a nonnegative multiplier as a consistency parameter that can be obtained from the consistency conditionḞ = 0 for inelastic flow in such a way that the following loading/unloading conditions must be satisfied during the inelastic deformation process. The scalar inelastic potential function Q in (25) is usually assumed different from the failure surface function F to govern non-associated flow rule, and thereby, the direction ∂ Q ∂σ of the inelastic strains is not normal to the failure function F. This can be fulfilled by adopting the Drucker-Prager surface as Q expressed in the effective stress space such that: where α p is the dilation constant. Thus, the flow rule (25) relates the failure loading surface and the constitutive law (14).
According to the experimental investigations on crack opening/closure effects in concrete [32], the inelastic strains progressively vanish and the material recovers its initial stiffness during the first compression loading phase, when load changes from tension to compression. During the crack closure, the evolution of inelastic strains is partially explained by the friction developed by the discontinuity lips. Moreover, a crack closing stress σ cr at which the accumulated damage in tension has no effect on the stiffness of concrete once cracks are completely closed can be taken as 10% of f c [1]. Therefore, the present model is enabled to predict crack closure mechanisms such as stiffness recovery and vanishing inelastic strains by allowing σ cr in the compression region. In this work, we adopted a value of about 2.5 (MPa) for the selected material parameters as σ cr . Thus, the damaged stiffness remains constant as long as the stress in compression is lesser than σ cr . As soon as the compressive stress goes beyond σ cr , the initial stiffness is recovered at the corresponding state of strain. In addition, the multiplier λ appearing in the flow rule (25) is introduced with two separate multipliers for tension and compression as follows:λ Similar to the definition of history variable (9),λ t andλ c deactivate the effect of inelastic strains developed due to tensile loading on the compression phase. Hence, the evolution of inelastic strains due to tension does not have any influence on the compressive behavior and vice versa. Meanwhile, the loading/unloading conditions (26) are still valid, butλ t orλ c is itself a monotonically increasing multiplier.

Influence of model parameters
In order to study the influence of the model parameters on the softening behavior of the material, tensile loading under 1D-situation is considered to cause homogeneous deformation. In such case, κ = = ε. The damage evolution and the corresponding stress-strain curves are obtained for varying the values of model parameters β 1 , β 2 , and β 3 . Figure 2 shows the stress-strain responses for increasing values of β 1 , whereas β 2 is kept constant as 0.18. It is observed that a higher value of β 1 makes the initial damage growth faster and thus leads to a more brittle response of the stress-strain behavior. Moreover, it is worth mentioning that the same evolution law (2) can On the other hand, Fig. 3 shows the stress-strain responses for increasing values of β 2 , whereas β 1 is kept constant as 0.75. It is noticed that a higher value of β 2 leads to a faster damage growth and results in zero residual stresses in the post-peak regime of the curve.
Additionally, the influence of model parameter β 3 on the stress-strain behavior is studied for β 1 = 0, β 2 = 0.105 with inelastic strain evolution and corresponding stress responses are displayed in Fig. 4. Having β 3 = 1.0 reduces the damage law (17) to its initial form (2), which means that there is no effect of β 3 on the stress response. Therefore, β 3 ≥ 1.0. The higher values of β 3 yield higher peak stresses by delaying the damage growth.

Consistency of thermodynamics
The process of undergoing deformations of any kind is generally a thermodynamic irreversible one [31]. Therefore, any constitutive modeling of material behavior must be consistent with the thermodynamic principles. According to the thermodynamic laws, a generalized form of the Clausius-Duhem inequality for any allowable dissipation process is given by where ρ, η, T , and q refer to the material density, the entropy density per unit mass, the absolute temperature, and the heat flux, respectively. Within this context, the Helmholtz free energy ψ per unit volume can be written in terms of internal state variables characterizing the nonlinear deformation behavior of concrete both in tension and in compression such that: As similar to the split of total strains (12), the Helmholtz free energy ψ can also be given as a decomposition of elastic ψ el and inelastic ψ in parts such that: The above decomposition is valid when the effect of damage is considered on the elastic properties and implicitly on the inelastic properties by means of inelastic strains and a reduction in limiting strength, as the proposed model is intended for brittle materials like concrete. The elastic free energy ψ el based on the principle of energy equivalence is given by Herein, the elastic free energy (32) uses a quadratic degrading function due to the assumption of energy equivalence principle in the material law (1). This kind of quadratic degrading function is usually adopted for phase-field variables in phase-field models [12] that lead to a constant broadening zone of localization. A detailed comparison between the gradient damage models and phase-field approaches discussed can be found in [21]. Nevertheless, the model presented in this paper does not fall into the category of gradient damage models discussed in their comparison, as the present model will use a nonlocal equivalent strain by following the work [36], which is explained in Sect. 5. Furthermore, the assumption of linear isotropic hardening functions (20) and knowing the dependency of ψ in on hardening variables h ± enable one to postulate the inelastic free energy ψ in as in the following form which is similar to the form of inelastic dissipation used in the Lubliner-Lee criterion [23]. Herein, g ± is the dissipated energy density during the fracture process [23,24]. Taking time derivative of Eq. (31), the rate of change of free energy can be written aṡ For any isothermal and adiabatic mechanical dissipation process (Ṫ = 0 and q = 0), the Clausius-Duhem inequality (29) reduces to σ :ε −ψ ≥ 0.

Substitution of Eqs. (13) and (34) into Eq. (35) yields
where Y , ∂ψ in ∂h + , and ∂ψ in ∂h − describe the thermodynamic conjugate force to damage (often termed as damage energy release rate) and thermodynamic conjugate forces to inelastic strain evolution under tension and compression, respectively. The damage energy release rate Y is defined as For any admissible process of energy dissipation that satisfies Eq. (36), it follows that: and the intrinsic mechanical dissipation potential function D can be expressed in terms of associated thermodynamic forces and the evolution of internal variables as The process of damage and inelastic deformation can occur without any influence of each other, and therefore, the dissipation potential function D can be decoupled into dissipation potentials due to damage D d and inelasticity D in such that where (or) To satisfy the dissipation function (40), the following relations must be valid: remains positive because of the convexity of the selected inelastic potential function (27) and Eq. (44b) remains negative, asḣ ± ≥ 0. Thus, the non-negativity of Eqs. (41)-(43) proves that both the damage and inelastic processes respect the Clausius-Duhem inequality (39) or (40).

Gradient enhancement
In this section, how the gradient enhancement of the local continuum model from a general concept of nonlocal continua is achieved is presented. Any nonlocal variablek at each material point x can be written as the weighted average of its local counterpart k over its surrounding volume domain V [5,37]: where g(ξ ) is a certain weight function and the vectorial variable ξ represents the distance from the point x to a point in its vicinity. As the averaging of damage variable leads to spurious residual stresses and an expansion of the softening zone across the bar [5], and also the damage growth in the present model is driven by the equivalent strain , the averaging procedure is therefore applied to the equivalent strain . Consequently, the nonlocal quantity¯ is defined asˆ The integral (46) is approximated by a partial differential equation according to [30,35], which is given bŷ where l c is the characteristic length scale and ∇ 2 is the Laplacian operator. An appropriate boundary condition (either a Dirichlet boundary or a Neumann boundary) must be used to solve the partial differential equation (47). Nevertheless, the present work adopts a natural boundary condition at every point of the boundary as provided below: The above boundary condition (48) ensures that the average of nonlocal equivalent strainˆ over the entire domain equals that of local equivalent strain [36]. Hence, κ is related to the distributed nonlocal equivalent strainˆ in the gradient-enhanced damage model. Subsequently, the damage criterion (5) is rewritten as follows: The Kuhn-Tucker loading/unloading relations (7) even so remain valid. Assuming the length scale l c → 0, the local continuum model is resumed asˆ = .
0.6 f c

Validation of the model equations
In order to check the ability of the proposed damage model in predicting the material behavior, the numerical algorithm to solve the model equations has been implemented into a finite element program (an in-house code [39]). The performance of the model is investigated under several loadings such as monotonic uniaxial tensile and compression, direct cyclic tension and compression as well as reverse cyclic loading. In order to prove the model equations, a single 8-noded brick element of size 200 (mm) × 200 (mm) × 50 (mm) is adopted for the numerical analyses using 8 Gauss integration points. The material parameters are taken from the work [23] so as to compare the numerical results of the present model with the respective experimental results [16,20,26]. The model parameters β 1 and β 2 have been chosen accordingly so as to reproduce the material behavior under respective monotonic loadings. Both the material and model parameters used in this work are provided in Table 1.
In order to investigate the solution procedure without localization, the material length scale l c is considered as 200 (mm). Linear shape functions have been used for the nonlocal equivalent strainˆ . The nonlinear system of resulting equations has been solved with the Newton-Raphson method. Displacements are imposed as prescribed boundary conditions. Numerical results are obtained and plotted against the experimental results in Figs. 6, 7, 8, 9, whereby the quality of the numerical simulation depends on the selected model parameters.

Uniaxial tension and compression
The single element is subjected to uniaxial tensile and compression. The numerical responses obtained using the model are compared with experimental results in Figs. 6 and 5.
As shown in Fig. 6, the simulated stress-strain behavior agrees well with the experimental results [16]. It can be noted that the model exhibits a gradual decrease in stress with increasing strains, i.e., strain softening and ensuring the residual stresses to a large extent. It means that damage increases monotonically.
Similarly, the numerically predicted compressive stress-strain responses depict very good agreement with the experimental compressive curve [20], which is noticed in Fig. 5. Moreover, the proposed model can predict the hardening as well as the post-peak softening behavior well.
In addition, it is obvious that the damage initiation in the case of compressive loading is later than that in case of tension, which is physically reasonable. Furthermore, the nonlinear behavior in the case of coupled inelastic damage model is provided by both the dissipative mechanisms of damage and inelasticity, even though hardening functions are assumed as linear functions.

Uniaxial cyclic tension and compression
The ability of the model is checked under uniaxial cyclic tensile as well as compressive loadings. Figures 8  and 7 depict the simulated tensile and compressive stress-strain behavior, respectively. The numerical results are compared with the experimental investigations according to Gopalaratnam and Shah [16] and Karsan and Jirsa [20] and also compared with previously proposed models of Lee and Fenves [23], which is the extended model version of Lubliner et al. [24] and Zhang and Li [42]. The model predictions agree with the phenomena observed by experiments: (a) the softening behavior under tension predicted by the model is fairly comparable with the experimental results, as shown in Fig. 8, and the results are also close to the previously proposed models; (b) As shown in Fig. 7, the nonlinear hardening behavior under compression is close to the experiments and previous models, but the softening behavior reasonably agrees with the test data with few differences due to the calibration of model parameters using uniaxial monotonic tests; Nonetheless, the capture of the inelastic strains gives good approximation of the material behavior; (c) the history variable κ t or κ c grows monotonically, as every reloading follows the path of the previous unloading, and thus, thermodynamic consistency is verified; (d) the tensile damage or compressive damage is monotonically increasing, which can be observed from the decrease in the slopes of the unloading behavior; (e) the history variable and damage remain constant during unloading regions;  Thus, the model is able to simulate nonlinear deformation behavior of concrete satisfactorily, but only the coupled model with inelasticity and damage is able to capture the residual permanent strains during the unloading cycles.
The hysteresis behavior of the experimental investigations is not included in the model, since unloading and reloading are described with constant stiffnesses. An extension of the model with respect to the hysteresis behavior is possible, if the inelastic strains are modified with respect to the strain velocity.

Tension-compression loading
In order to investigate the effect of damage under uniaxial cyclic tension-compression loading, the displacement history is used to be cyclic with an increasing magnitude alternatively tension and compression as an imposed loading as performed in [26]. The numerical results are illustrated in Fig. 9 with the loading path (O-A-B-O-B-D'-O-C-E). Additionally, loading and unloading cycles are labeled as 1 -4 . During cycles 1 and 2 , the first and the second tensile loading and unloading occur, respectively. During cycle 3, the loading changes from tension to first compression phase and further continues in cycle 4 . As observed in Fig. 9a, κ = κ t or κ = κ c is monotonically increasing during tension-compression loading/unloading processes, respectively. Consequently, a monotonic increase in the effective damage D during the respective load stages is depicted in Fig. 9b. It is also observed that D becomes zero after the tensile loadings (cycle 1 and 2 ), and thus, the accumulated effective damage in the previous tensile loadings does not affect the stiffness of the material in the first phase of compression (cycle 3 ). It means that the initial stiffness of the material is recovered for the first compression phase. Further loading in the compression phase (cycle 4 ) causes continuous growth The corresponding stress-strain responses using elastic damage and inelastic damage models, which show unilateral behavior, are compared with the previous model [23], as illustrated in Fig. 9e, f. The degradation of stiffness is shown by softening, as damage develops due to the opening of micro-cracks during the first tensile loading path (O-A-B). Then, the tensile unloading that occurs follows the path (B-O), and thus, the unloading slope refers to the damaged stiffness. Afterward, the tensile reloading (O-B-D') that occurs recovers the damaged stiffness caused by the previous tensile loading (Ḋ = 0 asκ = 0). Subsequent degradation of stiffness occurs, as cracks open up during the cycle 2 . Subsequent unloading that follows (D'-O) keeps the damaged stiffness unaltered temporarily. When the sign of loading changes, i.e., the first compressive loading (O-C) occurs, the initial stiffness is completely recovered due to the closure of opened micro-cracks. At this level, there is no damage development, as observed in Fig. 9e, since the damage criterion is not fulfilled. During the cycles 3 and 4 , if the loading occurs along the path (C-E) and further continues, damage grows due to compression, and consequently, degradation in stiffness is resulted which is indicated by the nonlinear increase in stress and subsequent softening after the peak stress. But there are certain discrepancies due to the fact that the permanent inelastic strains are not taken into account in the elastic damage model.
Nevertheless, as can be seen in Fig. 9f, the evolution of the history variable and damage during the loading/unloading and reloading cycles demonstrates that the thermodynamic laws are respected by showing the monotonic increase in the history variable and corresponding damage either in tension or in compression. Moreover, the main discrepancy of the elastic damage model is overcome, as the inelastic strain evolution (OO' and O'O") is observed after every unloading cycle. The initial stiffness recovery also occurs during the first compression phase, since the inelastic strains from the previous tensile load histories become inactive at the pointŌ. Thus, the present model is able to reproduce similar results of the previous model and the present crack closing behavior shows an improved response with the crack closing feature such that the inelastic strains are vanished compared to the response of the previous model, because the previous model recovers the initial stiffness as soon as the compressive stresses come into picture.
Further, to verify the crack closing behavior observed by the model, the numerical tension-compression stresses are normalized and compared with the experimental results [26,32] in Figs. 10 and 11. In Fig. 10, the tensile and compressive stresses are normalized by uniaxial tensile strength f t and uniaxial compressive strength f c , respectively. Similarly, the experimental response [26] is also normalized by respective material data from their work to compare with the numerical results. On other hand, Fig. 11 shows both the tensile and compressive stresses of numerical and experimental data [32] are normalized by the uniaxial tensile strength f t from the respective data. From both the figures, it can be noticed that the normalized stress-strain responses agree fairly with the experimental data [26,32] by showing adequate crack opening/closing effects such as damaged stiffness, capturing residual strains, stiffness recovery, and vanishing of inelastic strains on the compression behavior.

Localization phenomena
In order to visualize the localization phenomena of the proposed gradient-enhanced damage model, a simple one-dimensional problem is analyzed using a bar of size 100 (mm) × 10 (mm) × 10 (mm) subjected to a pure tensile loading, as shown in Fig. 12. The left end of the bar is fixed. A pure tensile loading is applied at the free end of the bar by prescribing an imposed displacement. Model parameters and material parameters (Data 2) used are from Table 1. To trigger the damage localization, the initial damage threshold κ 0 is adopted from Table 1 in a 10-mm-wide zone in the middle of the bar, whereas κ 0 = is set to 1 × 10 −3 in the rest of the portion. Thus, the reduction in κ 0 in the middle zone is about 13%. It is a well-known fact that the order of shape functions for the displacements should be at least one order higher than that of shape functions used for nonlocal variables so as to avoid stress oscillations. Therefore, linear shape functions have been used for the nonlocal equivalent strainˆ . The bar is discretized into 40, 80, 160 equal-sized solid (Hex27) elements along the loading direction (x-axis), respectively. The three different meshes are shown in Fig. 13. The computations are performed using different internal length scales, namely l c = {5, 7.5, 10} (mm). Figures 14, 15 show the nonlocal results obtained during the numerical analyses. Figure 14a displays the applied displacement during analysis. The stress-displacement responses generated for the three meshes using various internal length scales are depicted in Fig. 14b. It is clearly observed that the softening responses of every mesh are almost identical and converging faster. It is obvious that a finite energy dissipation has been obtained due to mesh-independent solutions. In addition, Fig. 14b illustrates the influence of characteristic internal length scale l c on the softening behavior. It is also apparent that the brittleness of softening responses is influenced by the internal length. The model responses become brittle; the model loses its stiffness and the peak stress value reduces for decreasing internal lengths, and vice versa. It is due to the fact that the larger internal length makes the localization zone wider. Moreover, the softening responses are also physically meaningful. Figure 15a, b displays the distributions of nonlocal equivalent strains and resulting damage over the localization zone that are obtained using the different length scales for the fine mesh (n = 160) at the loading step t 3 . Indeed, a larger value of l c allows wider spreading of the nonlocal strains, and accordingly, the damage becomes wider. In other words, a smaller value of l c yields strong (narrowing down) localization of nonlocal strains as well as damage at the center of triggered zone.
Furthermore, Fig. 15c compares the evolution of local and nonlocal equivalent strains corresponding to the three different loading steps (t 1 , t 2 , t 3 ) in the post-peak regime, which are obtained for the fine mesh (n = 160) with l c = 5 (mm). As can be observed, the local equivalent strains (dashed lines) are concentrated near the vicinity of the triggered region and further loading increases the narrowing of zone leading to propagation of macro-crack. On the other hand, the nonlocal equivalent strains (thick lines) are spread over the localization zone, physically meaning the occurrence of microstructural interaction in the neighborhood.
The respective evolution of damage is plotted over the length of the bar in Fig. 15d. The continued loading leads to the subsequent development of damage. The damage formed in a broader area. Nonetheless, the narrow region of strains that obtained replicates formation of the finite width of localization. Thus, the gradient enhancement of the model describes the physical meaning of crack initiation, growth of micro-cracks, and coalescence of these into a macro-crack, as observed in [35].

Concluding remarks
A simplified constitutive model for concrete has been developed based on energy equivalence concept to take its distinct deformation behavior under monotonic and cyclic loading conditions into account. A unified equivalent strain is considered as the local measure of deformation driving the evolution of isotropic damage. The crack opening and closure of micro-cracks under cyclic and reverse loadings are described by the evolution of two independent history variables. The implicit gradient method is incorporated to enrich the local equivalent strain. The efficiency of the damage models is demonstrated by validating the numerical responses with the experimental results under several load tests. Thus, the model can simulate the unilateral behavior of concrete well by showing adequate crack opening/closing effects such as damaged stiffness, inelastic residual strains, and stiffness recovery and vanishing of inelastic strains for the compression phase. The investigations of the one-dimensional bar problem prove that the ill-posedness of the boundary value problem is avoided by the incorporation of the implicit gradient method. Moreover, the ability of the present model in producing a constant width of damage zone is illustrated qualitatively.