Modeling damage behavior of concrete subjected to cyclic and multiaxial loading conditions

Concrete exhibits a different deformation behavior in tension and compression when subjected to uniaxial and biaxial loads. It also reveals crack‐opening and crack‐closing effects under cyclic loads. This paper discusses a gradient‐enhanced continuum damage model to investigate the deformation and damage behavior of concrete subjected to cyclic and multiaxial loading conditions. An activating variable introduced in the model distinguishes the multiaxial stress states. A single damage criterion that relates the history variable to the equivalent strain limits the elastic domain. The evolution of internal variables confirms their thermodynamic consistency. The gradient enrichment of equivalent strains leads to proper localization of deformation. Comparison of numerical predictions with relevant experimental data and some previous models validates the capability of the model to describe softening/crack‐closure behavior under various loading. Application of the model to fracture problems demonstrates the ability of the model to describe damage initiation and propagation, and its localization.


| INTRODUCTION
The failure of many civil engineering structures is caused by the failure of materials used in construction. Concrete is a widely used construction material. It exhibits brittle behavior under tension due to its reduced tensile strength but exhibits ductile behavior under compression. The brittleness of concrete results in the opening of tensile cracks. The propagation and coalescence of microcracks consequently lead to failure of structural components. Therefore, a proper understanding of the mechanical deformation behavior of concrete is essential to use continuum damage mechanics to predict its material behavior under various loading like cyclic/reverse loading and multiaxial loading. Several authors employed theories of plasticity to simulate concrete behavior, 1 although the inelastic behavior of concrete is dominated by microcracking. On the other hand, elastic-damage models, 2,3 coupled plastic-damage models, 4,5 freeze-thaw damage model, 6 and lattice model 7 have also been adopted to describe the fracture in concrete under various loading and environmental conditions.
The asymmetric behavior of concrete in tension and compression is characterized by means of a single damage Discussion on this paper must be submitted within two months of the print publication. The discussion will then be published in print, along with the authors' closure, if any, approximately nine months after the print publication. variable 8 or additive or multiplicative coupling of damage variables. 4,[9][10][11] The growth of damage variables is associated with either a single loading surface 10,12 or two different (tensile/compressive) loading surfaces. 8,13 The history variable governing the evolution of damage is usually related to a local state of deformation. The local deformation can be measured in different ways. 3,14 Most of the elasticity based models adopt the equivalent strain definition namely Mazars criterion, 8 modified von-Mises strain, 3,15 modified criteria of Mazars or Drucker Prager, 16 and cracking and crushing strains, 13 respectively. Nevertheless, these definitions do not completely capture the experimentally shown limiting domains of elasticity and ultimate stress. While comparing with experimental data, numerical results show considerable differences, especially in bicompression and complex regions. 16 Likewise, the numerical predictions 5,13 are not smooth enough in the complex region, though these models offer good results with few differences near bisector of bicompression region. Moreover, the direct use of unified damage model 12 would lead to an overestimation of the damage initiation surface in the complex regions where the compression is dominant.
On the other hand, the implementation of damage models of strain-softening materials like concrete yields mesh-sensitive and physically objectionable solutions upon mesh refinement and thereby causes ill-posedness of the boundary value problem. Thus, local continuum damage approach remaining a mathematically meaningful concept does not describe the microcracking pattern in strain-softening materials. Therefore, several methods such as viscous regularization, Cosserat continua, and nonlocal continuum methods (integral methods or implicit/explicit gradient methods) that adopt internal length scales are suggested. Although the nonlocal continuum models have become popular in describing microstructural material interaction, the use of nonlocal integral models is not efficient due to its global averaging method and difficulty in linearization from a computational point of view. In contrast, the linearization in case of gradient continuum models is straightforward and the introduction of gradients in softening material models is now well established. 17,18 Hence, we apply a gradient enhancement of the continuum model herein.
The main objective of the present paper is twofold. Firstly, it has already been shown that the gradientenhanced unified damage model proposed by authors 12,19 is capable of describing softening, asymmetric/cyclic/reverse cyclic behavior of concrete using a unified damage criterion and also results in mesh-insensitive solutions in a onedimensional setting, but suitability of the model is verified herein under multiaxial loading conditions using a few model parameters identified by monotonic-uniaxial tests. To enable the model to capture appropriate envelopes of damage initiation and maximum stresses, the model defines an activating variable that appears in the function of history variables to identify tension or compression dominant regions adequately. To verify thermodynamic consistency, the non-negativity of internal variables is proven. Secondly, to validate the modified model and to demonstrate its ability in describing damage initiation, propagation, and localization, several proportional/non-proportional biaxial tests or cyclic tests and benchmark fracture problems like a single-edge notched specimen under tension and shear and three-point bending tests are simulated.

| GRADIENT-ENHANCED UNIFIED DAMAGE MODEL
The constitutive material law for quasi-brittle materials like concrete showing linear elastic behavior can be written as follows: where, σ and ε are the Cauchy stress tensor and the Green-Lagrangian strain tensor, respectively; C 0 is the fourth-order elasticity tensor of the undamaged material. The principle of energy equivalence is assumed. The scalar variable D ranging from 0 to 1 refers to an active isotropic damage governed by a history variable κ such that: where, κ 0 , β 1 , and β 2 are model parameters controlling the exponential softening behavior of the material either in tension or in compression. 14 As inelastic strain evolution is not considered here, the constitutive law (1) describes the nonlinear deformation behavior of concrete as a pure damage dissipation. The averaging of damage variable leads to spurious residual stresses and an expansion of the softening zone across the bar. 20,21 Therefore, the damage growth is considered to be driven by the local equivalent strain ε eq and the nonlocal averaging of the equivalent strain ε eq at each material point is considered in the gradient-enhanced modeling. Consequently, κ must be related to nonlocal or gradient-enhanced equivalent strain b ε eq . Hence, analogous to yield or failure surfaces of the multiaxial material behavior, an associated single damage surface criterion f is conveniently adopted as to limit the elastic domain and also to allow the growth of damage. Thus, the nonlocal variable b ε eq is defined over its surrounding volume domain V as where, ξ is the coordinate referring the distance from the point x to a point in its vicinity and g ξ ð Þ is a certain weight function. As per, Refs. [3] and [18] the implicitgradient method that requires only a C 0 continuous displacement approximation is adopted herein. Thus, the partial differential approximation of the above integral (4) is expressed as follows: where, l c and r 2 are the internal length scale and the Laplacian operator, respectively. While solving the system of equations of the boundary value problem, in addition to imposing an essential boundary condition on displacements, the partial differential Equation (5) requires an additional boundary condition at every point of the boundary.
Since imposing a boundary condition on nonlocal strain b ε eq has physically no meaning, the nonlocal strain b ε eq is left free. However, the use of a boundary condition on the gradient of nonlocal strain b ε eq is naturally satisfied with the weak form. Hence, the natural boundary condition is adopted herein as follows: such that the average of nonlocal variable b ε eq over the entire domain equals that of its local variable ε eq . 3 The evolution of the history variable can be defined by the Kuhn-Tucker loading/unloading relations 2,15 as follows: The consistency condition _ f ¼ 0 must always be fulfilled during the damage loading process. If _ κ ¼ 0, the damage remains constant, that is, _ D ¼ 0. On the other hand, unlike various criteria determining the local deformation such as Brekelmans criterion based on principal stresses 2 or Mazars criterion based on principal strains, 10 or μ-model criterion 13 and de Vree criterion 15 based on strain invariants, the local counterpart ε eq is defined here in terms of stress invariants as a unified measure: for characterizing both tensile cracking and crushing failure. Herein, σ p is the elastically predicted stress tensor; I 1 is the first invariant of σ p ; J 2 is the second invariant of the deviatoric part of σ p ; σ max is the maximum principal stress of σ p ; and a Heaviside function H is equal to 1 for σ max > 0 and 0 for other values. According to, Refs. [9] and [22] the dimensionless constants α and β are defined depending on the maximum strength properties of material in uniaxial tension f t and in equibiaxial and uniaxial compression f bc and f c , respectively.
To account for the crack-opening/closure effects under reverse cyclic loading, the variable κ is written as which describes a maximum deformation of the material occurred during tensile or compression load history. Herein, two independent history variables κ t and κ c as a monotonically increasing parameters during loading/unloading processes are mathematically expressed as where, κ 0t and κ 0c are the initial threshold values at which damage initiates either in tension or compression. Thus, the expression (9) along with the unified measure (8) must be sufficient to describe the asymmetric or direct/reverse cyclic behavior of concrete, since the local measure (8) itself recognizes whether the loading is in tension or in compression. In a similar fashion, the model parameters κ 0 , β 1 , and β 2 are subsequently updated as follows: using the independent model parameters under tension (κ 0t , β it ) or under compression (κ 0c , β ic ). Figure 1 displays the uniaxial tension and compression responses and biaxial damage initiation surfaces predicted by previous models from literature and the damage model using the unified damage criterion (9). Comparing to experimental biaxial surface, the numerical predictions of biaxial damage initiation envelopes due to various damage criteria show considerable differences, especially in bicompression and complex regions, even though the models offer good results under uniaxial tension or compression. The differences are appeared due to the use of different local damage measures. It can also be observed that the damage criterion (9) overestimates the damage initiation or the elastic domain when complex load states occur, especially in the compression predominant regions.

| MODIFICATION FOR MULTIAXIAL STRESS STATES
To overcome the limitations discussed in Section 2, the damage criterion (9) has to be modified to capture the elastic domain and failure surfaces appropriately to validate with the experimental data. This necessitates further introduction of a variable in the model that distinguishes between tension dominant and compression dominant areas. Thus, the additional variable δ r is embedded into the criterion (9) to deactivate tensile domination where the compression is predominant. Thus, the expression (9) and corresponding initial threshold (14) are modified as F I G U R E 1 Uniaxial responses and biaxial damage initiation surfaces for various models The activation variable δ r is defined as where, r is the triaxial weight factor 9 and its range lies 0 ≤ r ≤ 1. For pure tension, r ¼ 1 and, for pure compression, r ¼ 0. Similarly, the model parameters used in (2) are written for various load paths as follows: Thus, one can observe that the original model is resumed when δ r ¼ 0, that is, Equation (13) is reduced to Equation (9).

| VERIFICATION OF THERMODYNAMICS
Any process of undergoing deformations is generally thermodynamically irreversible. 23 Therefore, any constitutive model must be consistent with the thermodynamic laws. According to thermodynamics, a generalized form of the Clausius-Duhem inequality for a process of an allowable dissipation is given by where, ψ, ρ, η, T, and q refer to the Helmholtz free energy per unit volume, the material density, the entropy density per unit mass, the absolute temperature, and the heat flux, respectively. For any isothermal and adiabatic mechanical dissipation process ( _ T ¼ 0 and q ¼ 0), the Clausius-Duhem inequality (18) reduces to Within this context, the Helmholtz free energy ψ can be expressed in terms of internal state variables characterizing the nonlinear deformation behavior of concrete both in tension and compression based on energy equivalence such that: After taking time derivative of Equation (20), subsequent substitution of the resulting equation for the rate of change of free energy _ ψ into Equation (19) This inequality (21) should be satisfied for any admissible process of energy dissipation, and hence the following state law becomes valid Consequently, by defining the damage energy release rate Y that describes the thermodynamic conjugate force to damage variable D and substituting into Equation (21), the inequality reduces to Herein, D denotes the intrinsic mechanical dissipation potential function of the material per unit volume expressed in terms of associated thermodynamic forces and the evolution of internal variables. During the dissipation process, Y remains always non-negative (Y ≥ 0), as _ D ≥ 0 and also _ κ t ≥ 0 and _ κ c ≥ 0 are valid. Moreover, the derivative ∂D ∂κ t=c is non-negative, as long as 0 ≤ D ≤ 1 holds. Thus, the non-negativity of internal variables, such as Y , D, κ t , and κ c proves that the damage process respects the Clausius-Duhem inequality (24).
To verify the validity of the admissible thermodynamic inequality, the evolution of internal state variables is taken into account under uniaxial reverse cyclic loading. Figure 2 displays the evolution of external/internal state variables during the loading-unloading and reverse loading processes. As seen in Figure 2 (left), ε 1 is a signed variable depending on tensile/compression loading, which may be increasing while loading/reloading or decreasing while unloading and changes negative sign when loading in compression. But the expression (8) of ε eq yields always positive during the entire load history. Consequently, b ε eq also becomes a positive function. Figure 2 (right) clearly discloses that the damage energy release rate Y is always a non-negative function, regardless of various stages of loading (oab, oc), unloading (bo, co, do), and reloading (od, oef) shown in the left. This is only possible, when the following two statements are valid as observed in Figure 2 (right) during the loading/ unloading processes: • κ t or κ c must offer a monotonic increase independently and • consequently, the active damage D reveals a monotonic increase in the respective load phases.
Besides the non-negativity of the internal state variables, such as Y , D, κ t , and κ c , the derivative term that appears in the inequality would keep the dissipation function (24) non-negative, because of bounding limits of D [0, 1] and the monotonic increase of history variables. Thus, the positivity of the Clausius-Duhem inequality (24) is numerically proven.

| VALIDATION OF THE MODEL
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 an in-house finite element program, which is also used for solving coupled problems like soil-structure-problems. 24 The validation of predicted results under uniaxial direct cyclic tension and compression and reverse cyclic tests are already reported, see. 12 The performance of the model is investigated under several biaxial tests such as several proportional load paths and non-proportional loading. The numerical results under the biaxial tests are presented in this paper for a comparison with the experimental results 25 26 is adopted for the numerical analyses using 8 Gauss integration points. Boundary conditions and prescribed displacements at the boundaries as imposed loading are shown in Figure 3. Table 1 provides the material parameters from the work 9 adopted in the validation process. The model parameters β 1 and β 2 are simply calibrated using uniaxial test results. The implicit-Euler backward method is adopted for time discretization. Knowing the strain increment at time t nþ1 , the values of stress and internal variables at t n , the updated values are estimated at the current step at t nþ1 . Use of smaller time steps yields accurate results in damage evolution. To investigate the solution procedure without localization, the material length scale l c is assumed as 200 [mm]. Linear shape functions are used for the nonlocal equivalent strain b ε eq . The system of resulting nonlinear equations has been solved with the Newton-Raphson method.
F I G U R E 2 Evolution of variables ε 1 , κ, Y, and D

| Monotonic proportional loading
Several biaxial load tests along various prescribed load paths are conducted to analyze the model. The present numerical results are plotted against the experimental data 26 and results of other models such as elastic-damage model 13 and plastic-damage model 9,22 in Figure 4 for bicompression tests, respectively. In all cases, the present model yields very good results by capturing maximum strength as well as sufficient softening; it also offers satisfactorily good results for the strain evolution in loaded directions, except for extension while comparing to experiments; yet the present results agree with the results of the model. 13 As reported from the experimental investigations, the larger extensional strains in the load free directions are due to stronger damage observed in that direction before fracture. This phenomenon has not been modeled by the present model. On the other hand, the other model 9 that considered volumetric dilatancy during the inelastic evolution offers a little different results. Hence, this present output becomes possible that the effect of isotropic damage is not considered on the Poisson's ratio and the volumetric dilatancy is not taken into account as well. Moreover, we underline that the parameters from Table 1 are kept same for all the test cases and the experimental data are only available until collapse of the specimen due to load control.
To illustrate the model performance in capturing biaxial elastic or failure envelopes, the elastic domains and failure surfaces from the present results are compared with the numerical results available from various other researchers 5,10,13,16 and the experimental results 26 in Figure 5. As can be seen from Figure 5(top), the envelope of damage surface (D ¼ 0), which limits elasticity domain obtained from the present model is compared with other numerical and experimental data. The numerical data 16 shown are obtained for the identical uniaxial tension and compression responses as mentioned in. 16 Nevertheless, when comparing to experiments, the bicompression region of elastic domains is either underestimated by both criteria of Mazars 10 and de Vree 16 or overestimated by all the combined or modified criteria, 16 except the μ model 13 with few differences and the present model with good agreement. Admittedly, the differences in other regions are due to assumption of damage initiation after the peak tensile strength. On the other hand, the present model captures the envelope of ultimate stresses (i.e., failure surface) well in bi-tension, bi-compression, and compression-tension regions, as pictured in Figure 5 (bottom). In fact, the other models 5,13 exhibit few differences in bi-compression region, yet better than the Mazars criterion. 10 In addition, the model 5 offers different results in the compression-tension region. Given the above, the present model provides quite good results in all the regions either at the damage initiation or at the failure during the damage process.

| Cyclic proportional loading
Biaxial cyclic compression with a stress ratio σ 2 =σ 1 ¼ 0:333 were conducted by Beam et al. 25 under stress control, which are extracted from the paper 29 to compare with the present results obtained under displacement control using the imposed loading as displayed in Figure 3 (right). For such purpose, the material parameters are chosen from the work 29 and the corresponding model parameters are calibrated using a monotonic compression test by 28 reported in. 29 Therefore, the comparison between the experimental data 28 and model results under monotonic compression for calibration (top) and the experimental data 25 and model simulations of the previous model 29 and the present model under biaxial cyclic compression (bottom) are presented in Figure 6, respectively. The agreement between the experimental and the previous model results is quite good; values of strains ε 2 both in the experiment and in the models are nearly zero; the few differences are due to omission of inelastic strains and no consideration of the effect of isotropic damage on the Poisson's ratio in the present model; however, the agreement between the results of experiments and the present model is similar only in the qualitative sense but not quantitatively.

| Non-proportional loading
Biaxial experiments with non-proportional loading were conducted by Maekawa and Okamura, 27 which were stress-controlled tests. But the numerical simulations are performed under displacement control. Therefore, the compression loading is first applied monotonically up to a certain level and then the tensile loading is superimposed in the other principal direction as shown in Figure 3 (right). Both the experimental results 27 and the numerical predictions using material parameters are illustrated in Figure 7. To compare the numerical results with experimental results, the calculated principal stresses are normalized by uniaxial compressive strength f c and tensile strength f t and the calculated principal strains are also normalized by the compressive strain ε c ¼ À25 Â 10 À3 ð Þ corresponding to f c . Additionally, the results of another model 29 are also compared against the experimental results 27 in Figure 7 (left). As can be observed from Figure 7 (right), the numerically predicted tensile stresses and corresponding strains behave nearly linear for compressive strains, say ε 2 =ε c < 0:5. The nonlinear behavior appears with increasing compressive strains. Subsequently, the increase in the tensile stress accelerates the deformation of concrete. Thus, the present predictions describe favorably the biaxial deformation behavior F I G U R E 4 Comparison between numerical curves and experimental data for biaxial compression loading observed in the experiments. In fact, the present model offers better results on the deformation behavior in the secondly loaded direction than the previous model. 29

| NUMERICAL SIMULATIONS
The capability of the proposed model in describing fracture processes is demonstrated by means of numerical analyses of well-studied benchmark problems such as a single-edge notched specimen subjected to tension, shear and a symmetric three-point bending of notched beam often found in literature 30,31 and a rectangular prism subjected to uniaxial reverse cyclic loading. 32 Prior to loading, the mesh is refined in the regions where crack propagation is expected to occur to capture the crack pattern properly. All the simulations are performed under displacement control. The material and model parameters are adopted from Table 1. For details on mesh sensitivity and the influence of the internal length scales on the simulation results of the present damage model, the work 12,19 may be referred.

| Tension test on single-edge notched specimen
A single-edge notched specimen is investigated using the proposed model under tension. It has a straight vertical notch running until its mid-height, located at a half distance of the left or right edge. The internal length is chosen as l c ¼ 0:2 mm ½ . The geometry and boundary conditions of the test are depicted in Figure 8 (left). As seen, horizontal displacements are prescribed at the complete right edge. The specimen is discretized with 3933 solid elements of 27-nodes. The stress-displacement curves obtained at the point A on the right edge of the specimen are shown in Figure 8 (right). The growth of damage occurs along the crack path normal to the loading direction and in turn the decrease of the stresses occurs in the post-peak region. Figure 9 displays the growth of nonlocal strains and corresponding damage propagation at several load stages in the post-peak region. As observed, the growth of nonlocal strains and damage run until the complete failure, as the deformation is well localized.

| Shear test on single-edge notched specimen
The same notched specimen used in the previous example is further investigated under shear loading. The geometry and boundary conditions of the test are depicted in Figure 10 (left). As shown, the test is simulated by keeping the right edge subjected to downward vertical displacements. The specimen is discretized with 4046 solid elements of 27-nodes. The internal length is chosen as l c ¼ 0:2 mm ½ . The stress response curves obtained at the point A are plotted against the applied displacements u y in Figure 10 (right). The initial decrease of stresses at A occurs due to the strain-softening at B during the growth of damage. Thereafter, further decrease of stresses at B upon subsequent load steps causes consecutive stress drops at A. However, the bumps that appear in the responses are due to the fact that the model also allows the damage development in the compressed region near the vicinity of notch tip. Figure 11 illustrates the results of nonlocal strains and damage propagation at several loading stages in the post-peak regime. As noticed, the nonlocal strains and damage grow until the damage reaches unity.

| Three-point symmetric bending test on notched beam
An another classical benchmark three-point bending test on a simply supported notched beam yielding a symmetric F I G U R E 9 Nonlocal strains (top) and damage (bottom) at u x displacements F I G U R E 1 0 Shear test: Geometry, boundary conditions, and stress response bending is investigated. The geometry and boundary conditions are illustrated in Figure 12 (left). The specimen is discretized with 1744 solid elements of 27-nodes. The internal length is set to l c ¼ 0:03 mm ½ . Figure 12 (right) illustrates the stress response at the point A where downward displacements are imposed. The distributions of nonlocal equivalent strains and damage evolution which are provided in Figure 13 are obtained at various loading stages after the strain-softening begins.
Thus, the growth of nonlocal strains clones the crack propagation and the stress-displacement behaviors are physically meaningful, although any relevant test data is not available. The present results of the proposed damage model for concrete fracture yield qualitatively a very good agreement with the results obtained by alternative computational modeling methods of brittle fracture. However, the selection of a proper internal length scale l c plays an essential role in fracture phenomena to achieve physically acceptable solutions.

| Reverse cyclic test on prism
In order to test the model further in the case of damage development under uniaxial reverse cyclic loading, the test conducted by Mazars et al. 32 on the concrete specimen 80 mm ½ Â40 mm ½ Â160 mm ½ subjected to uniaxial tension-compression loading is simulated. The geometric setup, boundary conditions, and imposed loading history are depicted in Figure 14. The material and model parameters are adopted from Table 1. The specimen is discretized with 400 solid elements of 20-nodes. The length scale l c is assumed as 160 [mm] to avoid localization such that the simulation is carried out to produce homogeneous deformation in the specimen. The actual stress-strain response is shown in Figure 14(c) and normalized stress-strain behavior with the respective maximum uniaxial tensile/compression strengths (f c and f t ) and maximum tensile strain ε c up to which material behaves linearly is plotted against the corresponding  Figure 14(d). The discrepancies between the numerical and experimental results are due to the fact that the model does not include inelastic strain evolution to capture the permanent strains. Nevertheless, the stiffness degradation and crack-opening and closure effects with full/partial stiffness recovery are completely described.

| SPECIAL FEATURES
The presented model provides a better understanding of concrete material behavior under cyclic and biaxial loading. Unlike other models, it is evidenced that: • the model could describe the nonlinear behavior of the material by an effective scalar damage itself in tension as well as in compression, • a single damage criterion could be easily adopted to capture the elastic domain, damage imitation and failure surfaces, appropriately, • the model could distinguish between tensile/compression stress states by the two independent history variables based on a Heaviside function and consequently, the damage growth in tension and in compression could be described using a single damage variable, • the various test scenarios could be simulated with the same set of a few model parameters that are calibrated from uniaxial tension and compression tests, and • the model could be conveniently be incorporated by a nonlocal method to avoid numerical difficulties in the FEM simulations.
Thus, the performance of the present model is adequately demonstrated by means of the numerical simulations of • tension, shear, and bending tests on notched specimen and • reverse cyclic test on prism like specimen.

| CONCLUSION
The gradient-enhanced unified constitutive model for concrete has been improved to take its asymmetric deformation behavior under multiaxial loading conditions into account by introducing an activating variable. The cyclic effects like crack-opening/closure of microcracks are also efficiently described by Heaviside function in terms of the two independent history variables for tension and compression. The thermodynamic principles are satisfied within the model formulation. The validation of numerical responses with experimental data proved that the model is capable of simulating the nonlinear deformation behavior of concrete under several biaxial load cases such as proportional or non-proportional and monotonic or cyclic loading conditions. Moreover, the agreement of numerical predictions with experimental results demonstrated the capability of the model to capture the material behaviors like softening, stiffness degradation, and stiffness recovery, adequately. Finally, the ability of the model in solving fracture problems has been illustrated by numerical simulations of benchmark problems such as a single-edge notched specimen subjected to tension, shear, and also symmetric three-point bending. The results evidenced the importance of such a regularized damage model to achieve damage initiation, propagation, and localization appropriately so as to understand the fracture behavior of concrete.