Concrete failure based on a localizing gradient‐enhanced damage formulation coupled with plasticity

Investigations of the mechanical behavior of concrete, including plastic deformability and failure mechanisms, have always been of vital importance. In this work, a modified yield criterion with three surfaces is proposed to capture the plastic yielding behavior under different loading conditions. The failure, that is, the gradual degradation and the final loss of material integrity, can be described through a damage concept. For the purpose of robust finite element (FE) implementations as well as the failure description at an appropriate scale, a localizing gradient‐enhanced damage model has been adopted. As such, one obtains mesh‐insensitive material responses, while ruling out the unphysical damage broadening observed in conventional gradient‐enhanced damage models with a constant regularization length. Following a consistent derivation, the plasticity‐damage coupled model is implemented into an in‐house FE framework, based on which an experiment‐inspired example is simulated to illustrate the capability of the proposed description in concrete modeling, especially the mixed‐mode failure. Lastly, the conclusion provides some final insights.

plasticity [4], Mohr-Coulomb plasticity [5], Drucker-Prager plasticity [6], amongst others.Before the eventual material failure takes place, many engineering materials firstly show plastic hardening behavior to a certain extent.Apart from the isotropic hardening that describes the isotropic expansion of the initial yield surface, kinematic hardening often serves to model the Bauschinger's effect in metallic materials when undergoing cyclic loading, see for example [7].Plastic softening manages to describe the softening behavior governed by the shrinkage of the yield surface with or without considering initial hardening behavior, see for example [8,9].However, during complex loadings, for example, a loading-unloading-reloading condition, it fails to capture the degradation in material stiffness.Therefore, it is necessary to incorporate certain criteria to account for the irreversible material deterioration.
For many materials, the major failure mechanism is a damage process -the formation, propagation and further coalescence of microscopic cracks.To overcome the numerical issues of the local damage description in the context of finite element (FE) modeling, that is, the loss of well-posedness of the partial differential equation (PDE) system, the nonlocal damage concept has been introduced.In the original nonlocal concept, it is interpreted such that the response of the material at a certain point depends not only on the interested material point itself, but also on its neighborhood.Therefore, the interested local internal variable is replaced by a weighted volume average of itself by means of volume integration, see [10].To further facilitate nonlocal damage modeling in the framework of FE, an alternative formulation is proposed in [11], by using the Taylor expansion of interested local variables, which in general would require shape functions of high order.This has motivated a different formulation -the so-called implicit gradient damage model in [12], where the mathematical reformulation is achieved by calculating the second order gradient of the nonlocal quantity, and further neglecting higher order gradient terms of the local quantity.The implicit gradient damage has been extended to numerous models, including anisotropic damage in for example [13], and microplane damage in for example [14,15].Despite the straightforward implementation, the implicit gradient damage model quite often leads to a spuriously broadening damage profile.In the attempt to address the broadening issue, it has been assumed in [16,17], that the regularization length in the implicit gradient damage model should be a variable.Even with different physical interpretations, the results in [16,17] both render thin damage profiles.The so-called localizing gradient-enhanced approach in [17] has attracted much research interest and has been investigated and extended to inelastic characteristics, for example, ductile failure [18][19][20].Another approach -micromorphic damage in [21,22], could also render a thin damage profile that even resembles fracture.The energetic formulation, however, poses a great challenge for the further employment in many constitutive material models.
In this work, a modified cap plasticity model has been proposed by a new definition of the yield function.The modification enables a more clear definition of the evolution law for plastic quantities as well as a simpler derivation for the FE implementation, while maintaining the advantages of the original model.Subsequently, the localizing gradient-enhanced damage model is coupled to the modified plasticity model.As such, the coupled model succeeds in the numerical simulation of a concrete structures failure characterized by multiple damage zones.For more detailed insights, the interested read is referred to [23].

A SMOOTH THREE-SURFACE CAP YIELD CRITERION
For the modeling of the complex behavior of concrete under different loading conditions, several plasticity models have been developed, including models with "caps".Generally, a nonlinear yield criterion is formulated, by partially replacing linear terms with nonlinear terms so as to render complex material behavior and/or stable numerical solutions.However, most of the existing cap yield criteria have proven to be insufficient due to the lack of  1 − continuity at certain points of the yield surface.Another issue that often arises with the cap models, is triggered by the piece-wise definition of the yield criteria.Potential oscillation of the solutions at Gauss point level could take place, when the resultant stress state is close enough to the intersection of different pieces of the yield function.The cap plasticity model proposed in [14] is capable of addressing these issues by a continuous definition of the yield surface.In their formulation, the nonlinear segments of the yield function, that is, the caps, are recovered by considering the multiplication of the linear part with additional terms.However, the nature of their yield criterion, that is, a function of stress quantities to the power of two, leads to the plastic strain with a unit of stress -a rather confusing concept.Therefore, an alternative yield function is proposed here, reading with As in [14], the definitions of  1 ,   and   remain unchanged.  (  , ) and   (  , ) are the components to be activated for the tension cap and compression cap. 0 and   are the yield stress and hardening modulus for the classic Drucker-Prager yield criterion. is another material constant that represents the tension-compression asymmetry of the yield criterion.(•) is the Heaviside function with subscripts  and  for the tension and compression cap, respectively.The threshold value    /   indicates the intersection of the linear segment and the tension/compression cap.Denominators of   and   in Equation ( 2), together with the   in  1 , allow the expansion of the yield surface with the plastic flow evolution.The material constant   /  serves to represent the ratio, in terms of the thermodynamic force evolution of the plastic hardening variable , between the tension/compression cap and the linear Drucker-Prager segment.As such, the advantages from the work [14] are still observed, while the derivation is somewhat simpler.As for the plastic strain evolution, it is defined to be split into two parts, that is, Meanwhile, the evolution of the plastic hardening variable is assumed to be The evolution of the plastic flow is summarized by the Kuhn-Tucker condition that reads λ ≥ 0, ℱ ≤ 0 and λ ℱ = 0. (5)

A LOCALIZING GRADIENT-ENHANCED DAMAGE MODEL
In the context of continuum damage mechanics, the implicit gradient-enhanced approach is probably the most popular one.For its straightforwardness in FE implementation, the implicit gradient approach has been intensively investigated and widely employed for the modeling of damage or other softening behavior.
The classic implicit gradient nonlocal formulation reads Here, η is the nonlocal counterpart of the local variable .The constant  0 serves as the regularization parameter, and is often linked to a micromechanics-based length-scale parameter via  0 =  2 . is the outward normal at the continuum boundary.With the implicit model in Equation ( 6), the regularization parameter  0 should be sufficiently large, often several times of the element size in FE implementations, in order to obtain a mesh-insensitive result with a satisfactory convergence rate of the numerical solution.As observed and discussed in several publications, however, the implicit model would lead to excessive damage evolution in the form of unphysical broadening of the damage zone.
In this work, the localizing gradient-enhanced damage formulation serves as the remedy for this unwanted effect.Specifically, the gradient effect decreases with the damage evolution.In other words, the model behaves initially as a nonlocal model, and gradually evolves to a somewhat local one.The strong form now reads The dependence of the regularization parameter  on the damage variable  provides a straightforward coupling of the nonlocal effect and the damage mechanism.The regularization degradation function () is adopted from [17], reading It satisfies the conditions (0) = 1 and (1) = , together with a negative first order derivative, that is,    < 0. The residual of the degradation  has a small yet non-negative value, in order to preserve the nonlocality of the model at a damaged state.The parameter  ≠ 0 controls the decreasing rate of the degradation.

PLASTICITY-DAMAGE COUPLING
When the damage variable  acts as the local variable  to be regularized by Equation ( 7), an insufficient evolution for the nonlocal damage variable would occur.Because the gradient enhancement would render a smeared distribution as well as a smaller magnitude for the interested variable.Even though the localizing gradient enhancement might alleviate this issue to some extent with  ≈ 0, potential numerical convergence difficulties might occur.Therefore, the solution of many research works is to take an intermediate variable for the regularization.Afterwards, the damage is to evolve based on the regularized variable.
Here, the plastic equivalent strain   eq acts as the local variable for regularization.To cover the difference in damage evolution rate between tension and compression, the definition is adopted.Here,  1 = tr(  ) is the trace of plastic strain, and the  2 quantity is defined as  2 = The material parameter  defines the damage evolution rate, and  0 is the threshold for damage initiation.Further, it is worth to note, that damage evolution is usually considered as an irreversible process.Therefore, the maximal nonlocal quantity during the entire loading history is considered to drive damage.Given that the damage is considered to be isotropic, the resultant stress is written as further leading to the strong form of the displacement field that reads Characteristics of the notched block test: (A) adopted geometry for the simulation, (B) experimental result [24], (C)-(D) damage evolution in the simulation and (E) representative load-displacement relationship.

NUMERICAL SIMULATION
The cement-based concrete material shows different behaviors under tension and compression.The cap plasticity model could partially capture such characteristics, while the failure mechanism also needs to be distinguished by the parameter , allowing different damage evolution rates from tensile and compressive loadings.The experimental investigations in [24] show that the failure mechanism of a SHCC structure with notches, when undergoing a shear loading, could be rather complex.Given the symmetry of geometry, boundary conditions as well as loading, the simplified setup is adopted for the simulation, see Figure 1A.With the area between two notches confined for both upper and lower sides by the experimental equipment, the specimen undergoes a loading from the remaining upper edges, rendering a complex failure pattern, see Figure 1B.
The evolution of failure consists of early tensile damage zones and a subsequent shear damage zone through the specimen.As shown in Figure 1C, the tensile damage initiates in both notches, from where the structure tends to be torn apart.This process occurs at an early loading stage, showing the brittleness of the material.Moreover, the lower damage zone appears slightly later than the upper one due to the boundary conditions as well as the loading condition.The shear damage starts to develop after a significant tensile damage evolution.Instead of the notch edges, the shear damage initiates from the interior region between the two notches, and subsequently cuts through the specimen rather rapidly, leading to the complete loss of load-bearing capacity of the structure, see Figure 1D.By the comparison between the experimental result (Figure 1B) and the simulation result (Figure 1D), one could conclude that the proposed model has successfully reproduced the failure pattern.The modeling of mixed-mode failure, that is, the coexistence of tensile failure and shear failure, could benefit from the appropriate definitions of separate criteria for tensile and shear failure evolution.This is quite often not so straightforward in brittle damage modeling.However, the three-surface cap plasticity model coupled with damage naturally provides this possibility with its different thresholds for plasticity initiation in different loading scenarios.
Further, the two parameters of the localizing model, that is, the residual percentage  and the deceasing rate  of the length-scale, are investigated in the context of this example.Figure 2A-D show the damage profile with different percentages for the residual length-scale, that is, 100%, 50%, 20%, 10%. Figure 2A represents the case with a constant length-scale parameter.The sufficiently large gradient effect for the PDE well posedness would in turn prevent the initial tensile damage to propagate much further.Instead, the damage in these regions becomes wider.The shear damage appears later, and also suffers from the widening issue.With decreasing residual percentage, one could observe more brittle results, see Figure 2B-D as well as the load-displacement relations in Figure 2E.During the first half of the simulation when only tensile damage evolves, there is no softening in the load-displacement relationship, meaning that the load bearing capacity remains still high.This is why the further evolution of tensile damage (with a smaller ) leads to a later appearance of shear damage and a higher peak load.The necessity of adopting the localizing gradient-enhanced approach is clearly explained by the comparison between Figure 2A and Figure 2B-D.The subsequent comparison between different values of the decreasing order parameter  shows, that a positive  leads to a large decreasing rate at an early stage of damage, while a negative  would have an opposite effect.It can be easily understood from Figure 3, that a negative  would preserve  a high level of gradient effect until the damage becomes significant; therefore, it leads to the broadening of the damage zone -a similar effect from a large .A positive , like a small , leads to a more brittle behavior, see also Figure 3E.

SUMMARY
In this work, a modified smooth three-surface cap model for plasticity has been coupled with a localizing gradientenhanced damage formulation for the modeling of concrete materials.The cap plasticity model serves to address the difference in plastic evolution under tension and compression.To avoid the unwanted damage zone expansion at a completely damaged state and to preserve the PDE well-posedness, the damage behavior is modeled by a localizing gradient-enhanced approach.Therein, the nonlocal counterpart of the plastic equivalent strain is to govern to the damage evolution.The work at hand is consistently formulated, derived and further implemented into the FE framework.The example illustrates the capability of the model in capturing the complex failure pattern of concrete materials.

A C K N O W L E D G M E N T S
The authors would like to acknowledge the financial support of the China Scholarship Council (No. 202008080328), and of Ansys Germany GmbH.Open access funding enabled and organized by Projekt DEAL.
[Correction added on 5 September 2023, after first online publication: Projekt DEAL funding statement has been added.] is Poisson's ratio.The magnitude of  indicates the level of the aforementioned difference, and  = 1 leads to a classical von-Mises definition.Based on the nonlocal quantity η calculated from   eq , the damage can be defined by  = 1 − exp(−( −  0 )), where  = max ≤ η().
U R E 2 (A)-(D) damage profiles (E) displacement-reaction force relations, with  = 1 and different residual parameter  for ().