A modified cap plasticity description coupled with a localizing gradient-enhanced approach for concrete failure modeling

The modeling of concrete has always been of great interest. In this work, a yield criterion with three surfaces is proposed to capture the plastic yielding behavior under different loading conditions, including uniaxial tension/compression, biaxial tension/compression and tension-compression. Material failure, characterized by the degradation and finally the complete loss of material integrity, can be modeled by the continuum damage approach. Within finite element methods, the nonlocal enhancement of integral-type or gradient-type is often required for the well-posedness of the partial differential equation system. In this work, a localizing gradient damage model has been adopted to obtain mesh-insensitive material responses, while ruling out the unphysical broadening of the damage zone often observed in constant length scale gradient damage models. Following a consistent derivation, the plasticity-damage coupled model has been implemented into an in-house finite element framework. Several representative and demonstrative examples serve to illustrate the capability of the proposed description in concrete modeling, followed by the conclusion, where some final insights are provided.


Introduction
Concrete has been one of the most important materials for engineering constructions.As a typical composite material, concrete is usually composed of aggregates of different sizes and types, that are jointed together by the hardened cement base.In terms of mechanical properties, it is renowned for its high compressive strength.Therefore, it is the ideal material for large scale infrastructures, where heavy loading is to be expected.However, the fact that its tensile strength is significantly lower makes it vulnerable to failure, e.g., when employed in a bending beam structure.In the process of the usually brittle failure development, complex mechanisms take place, along with observable cracks localized at certain critical regions.
B Bo Yin liam.yin@ansys.comB Michael Kaliske michael.kaliske@tu-dresden.de 1 Institute for Structural Analysis, Technische Universität Dresden, 01062 Dresden, Germany 2 Ansys Germany GmbH, 99423 Weimar, Germany Different strategies of reinforcement have been taken to alleviate the issues stemming from the aforementioned disadvantage.First possibility is through the adoption of long reinforcement bars made of steel, carbon, amongst many others.In this case, the reinforcement would keep bearing significant load, even after the failure of the concrete matrix.As a result, its load bearing capacity greatly increases, when the bonding between concrete matrix and reinforcement is preserved to some extent.The cracks, which would otherwise be highly localized, are more distributed and of smaller size.Textile reinforcement, which has become popular in recent years, also makes the overall performance of concrete structures much more ductile and durable when undergoing excessive loadings.The geometry of the reinforcement makes it especially suitable for structures like slabs; see the investigations in e.g.[16].Another category is the reinforcement by short fibers, e.g.Polyvinyl Alcohol (PVA) fibers and steel fibers.As an example, the strain hardening cement-based composite (SHCC) behaves as an approximately homogeneous and isotropic material when investigated on a macroscopic level.Meanwhile, the material has become rather ductile and its capability of enduring extensive deformation has greatly increased; see [2,18].The failure of SHCC is accompanied by multiple micro-cracking, fiber debonding, fiber rupture, and, finally, the formation of a macroscopic failure zone.
The deformation of concrete material, especially with reinforcement, shows a certain level of irreversibility.This, in part, could be explained by the plasticity mechanism.The canonical plasticity theory is characterized by the plastic flow evolution, governed by a certain yield criterion, through internal variables.Classical criteria for yield surfaces include Tresca plasticity [19], von-Mises plasticity [9], Mohr-Coulomb plasticity [1], Drucker-Prager plasticity [3], amongst others.With the evolution of the plastic flow, the response could be classified into perfect plasticity, plastic hardening and plastic softening.Even though, there would be material failure eventually, many of the engineering materials firstly show hardening behavior to a certain extent.Isotropic hardening is quite often employed, in order to describe the isotropic expansion of the initial yield surface.Another type of plastic hardening -kinematic hardening -serves to model the Bauschinger's effect in metallic materials when undergoing cyclic loading; see e.g.[27].Plastic softening manages to describe the softening behavior governed by the shrinkage of the yield surface with or without considering initial hardening behavior; see e.g.[24,28].It fails, however, to capture the degradation in the material stiffness during complex loadings, e.g. a loading-unloading-reloading condition.Therefore, this has necessitated the modeling of failure phenomenon by a different methodology.
For many materials, the major failure mechanism is the formation, propagation and further coalescence of microscopic cracks.The phenomenon of the loss of material integrity is termed as a damage process.The progressive degradation of the material stiffness is governed by an internal damage variable d, which is a load-history-dependent function.To overcome the numerical issues of the local damage description in the context of finite element (FE) modeling, i.e., the loss of well-posedness of the partial differential equation (PDE) system, the Nonlocal Damage concept has been introduced.The original nonlocal concept is interpreted such that the response of 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 [12].To further facilitate nonlocal damage modeling in the framework of FE, an alternative formulation is proposed in [10], by using the Taylor expansion of interested local variable, which in general would require shape functions of high order in FE method.This disadvantage has motivated a different formulation -the so-called Implicit Gradient Damage Model in [11], where 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.In spite of the straightforwardness in FE implementation, the implicit gradient damage model quite often leads to a spurious broadening damage profile.For this, the assumption has been made in [5,13], that the regularization parameter in the implicit gradient damage model should be non-constant.Even that the regularization parameter is defined in quite different, even somewhat opposite, manners, the results in these two works both render thin damage profiles.The so-called Localizing Gradient-enhanced Approach in [13] has been proven to be a promising method in computational damage modeling, and has been investigated and extended to inelastic characteristics, e.g.ductile failure [14,22,23].Another approach -Micromorphic Damage in [4,25], could also render a thin damage profile that even resembles fracture behavior.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 compared to [29].While maintaining the advantages of the original model, the modification enables a more clear definition of the evolution law for plastic quantities as well as a simpler derivation for the FE implementation.Moreover, the localizing gradient-enhanced damage model is coupled to the modified plasticity model.As such, the plasticity-damage coupled model succeeds in the numerical simulation of the failure in different concrete structures by means of capturing both the failure patterns and the load-displacement relationships in general.
The work at hand is structured as follows.In Sect.2, the governing equations of a three-surface plasticity model is firstly introduced, followed by the concept of non-constant implicit gradient-enhanced damage formulations, and the consideration of the plasticity-damage coupling.Section 3 describes the aspects regarding the numerical implementation, including the solution scheme on the material level, the stress and the algorithmic tangent, and the further derivations for FE computations.Section 4 follows with several demonstrative examples, where the model is firstly investigated from a pure numerical perspective and further evaluated through the comparison with experimental results.Section 5 summarizes the findings and closes the contribution with potential perspectives.
which is further split into elastic strain and plastic strain, reading The stress tensor is defined as where C e is the linear elastic material tangent that reads C e = κ1 ⊗ 1 + 2μP.Here, 1 is the second order identity tensor, and P is the fourth order deviatoric projection tensor defined as P = I − 1 3 1 ⊗ 1.I is the symmetric fourth order identity tensor.κ is the bulk modulus and μ is the shear modulus.The scalar-valued volumetric quantity of a second order tensor is defined as where The deviatoric quantity is defined as a second order tensor, reading Therefore, the stress can be split into volumetric and deviatoric components, reading

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 cap models.Generally, the yield function from a classic plasticity criterion, which is usually a linear function of stress quantities, are partially replaced by nonlinear terms so as to render complex material behavior and/or stable numerical solutions.However, most of the proposed nonlinear yield criteria are proven to be insufficient, due to the lack of C 1 continuity at the intersection of linear and nonlinear segments of the yield surface.Another issue often arises with the cap models, that the local numerical iteration for the evolution of plastic quantities might lead to the oscillation of the solutions that belong to two segments of the complete yield surface, since most of the aforementioned yield functions are defined as piece-wise ones.The cap plasticity model proposed in [29] is capable of addressing these issues by a continuous definition of the yield surface, reading with With such a formulation, the nonlinear segments of the yield function, i.e. the caps, are recovered by considering the multiplication of the linear part Here, σ 0 and H p are the yield stress and hardening modulus for the classic Drucker-Prager yield criterion.η is another material constant that represents the tensioncompression asymmetry of the yield criterion.H (•) is the Heaviside function with subscripts t and c for the tension and compression cap, respectively.The threshold value σ t V /σ c V indicates the intersection of the linear segment and the tension/compression cap.Denominators of f t and f c in Eq. ( 8), together with the H p in f 1 , allow the expansion of yield surface with the plastic flow evolution.The material constant R t /R c 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.
In [29], the evolution of plastic strain is derived based on the yield function from Eq. ( 7), reading The very nature of the yield criterion in Eq. ( 7), i.e. 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 in this work, reading 123 where the definitions of f 1 , f t and f c remain unchanged, i.e. as in Eq. ( 8).As such, the advantages from the work [29] are still observed, while the derivation, which is to be elaborated in the following sections, is somewhat simpler.As for the plastic strain evolution in Eq. ( 9), it is further split into two parts, i.e., 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

A localizing gradient-enhanced damage model
In the context of continuum damage mechanics, especially with respect to numerical investigations, the nonlocal concept of the continuum body often needs to be considered.Typical nonlocal approaches include the integral-type nonlocal approach, the explicit gradient-enhanced approach and the implicit gradient-enhanced approach.For its straightforwardness in finite element 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 c 0 serves as the regularization parameter, and is often linked to a micromechanics-based length-scale parameter via c 0 = l 2 .n is the outward normal at the continuum boundary.With the implicit model in Eq. ( 14), the regularization parameter c 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.One of the remedies for this unwanted effect is to consider a non-constant regularization parameter that evolves with certain local quantities.In [5,15], it is proposed that the regularization parameter in the governing equation Eq. ( 14) 1 is to be multiplied with another factor ξ , leading to where The factor ξ is defined as a function of the equivalent strain ε eq , reading The positively defined parameters ζ and ε crit eq lead to a monotonically increasing factor ξ before reaching its maximum 1.This so-called Transient Gradient Damage Model could help, at least in some cases, to avoid the aforementioned unphysical broadening of damage zone, even though the increasing gradient effect may seem counterintuitive.This is due to the fact that the model would initially behave somewhat like a local damage model when the factor ξ is still small.If damage occurs (d ≈ 1), when ξ 1, yet ξ is large enough to provide sufficient regularization, it is quite possible that one would obtain a thin damage zone from numerically convergent solutions.It is claimed by some researchers, however, that with sufficient loading after the initial evolution of the damage zone, the broadening might reoccur due to the fully evolved factor ξ , i.e., ξ = 1.In this case, the transient gradient damage model only serves to postpone the unphysical damage broadening.Furthermore, what still remains questionable is the physical interpretation of the increasing gradient effect.
The alternative, which might be easier to comprehend, is a decreasing gradient effect -quite the opposite to Eq. (17).It has been firstly proposed in [13], that 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.Analogously to Eq. ( 16), the regularization parameter is defined as with the strong form The dependence of the regularization parameter c on the damage variable d, instead of the equivalent strain ε eq , provides a more clear and straightforward coupling of the nonlocal effect and the damage mechanism.In [13], the regularization degradation function g (d) is defined as It satisfies the conditions g (0) = 1 and g (1) = R, together with a negative first order derivative, i.e., ∂ d g < 0. The residual of the degradation R has a small yet non-negative value, in order to preserve the nonlocality of the model at a damaged state.The nonzero order parameter n controls the decreasing rate of the degradation.Considering its aforementioned advantages over the model with increasing gradient effect, the gradient-enhanced model with decreasing effect, termed in [13] the Localizing Gradient Damage Model, is employed in this work.

Plasticity-damage coupling
The term Gradient-Enhanced Damage Model could lead to two interpretations.The first consideration is that the damage variable d acts as the local variable η to be regularized by Eq. ( 14).This would, however, lead to an insufficient evolution for the nonlocal damage variable, due to the fact that the gradient enhancement would render, not only a smeared distribution, but also a smaller magnitude for the interested variable.Even though the localizing gradient enhancement might alleviate this issue to some extent with R ≈ 0, potential numerical convergence difficulties might occur.As for the second interpretation, it is an intermediate variable, quite often a strain-based scalar quantity, that one considers to regularize.Afterwards, the damage evolution is considered as a function of the resultant quantity from the regularization.In this work, the second concept is adopted, with the plastic equivalent strain ε p eq as the local variable η.To cover the difference in damage evolution rate between tension and compression, one adopts the definition ν is Poisson's ratio.The magnitude of k indicates the level of the aforementioned difference, and k = 1 leads to a classical von-Mises definition.Based on the nonlocal quantity η calculated from ε p eq in Eq. ( 21), the damage can be defined by Here, β is a material parameter that defines the damage evolution rate, and ε 0 is the threshold for damage initiation.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 as the driving quantity of damage.Since isotropic damage is considered in this work, the resultant stress is therefore rewritten as 3 Numerical aspects

Return mapping algorithm for plasticity evolution
This subsection deals with the solution scheme for the evolution of plastic variables at the Gauss point level by a return mapping algorithm.When investigating into the variable evolution from the previous pseudo time step t n to the current one t n+1 with t n+1 = t n + t, two states are considered, i.e., the elastic trial state and the actual state.In the trial step, it is assumed that no increment occurs for the plastic variables, and the increment is therefore of pure elastic nature.Under this condition, the variables for the trial state read The subscript n + 1 in the subsequent equations is neglected for conciseness.Further, one obtains the stress quantities in the trial state and further splits it into volumetric and deviatoric components, reading σ tr = C e : ε e,tr ⇒ σ tr V = σ tr : V and σ tr D = σ tr : P.
Further, the yield function at trial state is evaluated.The consideration of intact stress σ , instead of the damaged stress σ = (1 − d) σ , is quite conventional in the damage community and leads to a less complex derivation.The condition F tr ≤ 0 indicates that the trial state represents the actual state of the material at the interested loading increment, whereas F tr > 0 requires an update of the so far neglected evolution of the plastic variables.In the latter case, the conditions F = 0 and λ > 0 need to be fulfilled by adopting an established Newton-Raphson iterative solution scheme at the material level.By including Eq. ( 25) into the stress at a plastic state, one obtains which can be further split into with the consideration of the plastic flow direction Therefore, the equation system with the unknowns { λ, σ V } that needs to be solved reads The consistent linearization reads with Therefore, the unknown increments are with which the unknowns are updated as for each internal iteration step.

Stress and algorithmic tangent
With the plastic evolution resolved through a local Newtonraphson scheme, the stress and algorithmic material tangent need to be subsequently calculated.From Eqs. ( 27) and ( 29), the stress reads Therefore, the algorithmic tangent is derived as where (37) To solve the remaining two unknowns, i.e. ∂ λ ∂ε and ∂σ V ∂ε , the unconditional equilibria are considered.The unknowns are therefore solved by (40)

Finite element implementation
For the implementation into the context of the finite element method, the global fields considered are the displacement field u and the nonlocal field η.The strong forms for these global fields read where σ = (1 − d) σ .Subsequently, the weak forms are obtained by the multiplication with test functions, i.e. δu and δ η, respectively, and the volume integration over the continuum domain .They read The finite element discretization is then achieved by where the shape functions for the displacement field N I u and for the nonlocal field N I η are considered to be identical for simplicity, i.e.N I u = N I η = N I .With Eq. ( 44) and the weak together with the element stiffness tangent where In Eq. (47) 3 , the derivative is nonzero only in the case of loading, i.e., when the current value of η is also the maximum in loading history.It is also worth to note, that the resultant stiffness tangent is unsymmetric in general.Therefore, an unsymmetric solver is required for the numerical solution.

Failure of an L-shaped structure
Inspired by [6,21,26], the first example studies the damage evolution of an L-shaped structure with the intention to illustrate the advantage of the localizing gradient-enhanced damage model (LGDM) over the conventional gradientenhanced damage model (CGDM), as well as over the local damage model (LDM).The geometry of the structure is shown in Fig. 1a along with the boundary constraints at the bottom edge and the displacement-type loading.
To simply elaborate LGDM, the influence of plasticity and plasticity-damage coupling are excluded from this example.Therefore, the bulk material is considered linear elastic, before damage comes into play.Meanwhile, instead of the plastic strain ε p , the total strain ε is considered for the calculation of the equivalent strain in Eq. ( 21).For detailed derivations, the interested reader is referred to [26].
Using hexahedral elements, the three dimensional boundary value problem is discretized into three different meshes, with 6867 elements for mesh1 (M1), 8640 elements for mesh2 (M2) and 11371 elements for mesh3 (M3), see Fig. 1b-d, respectively.As shown in Fig. 2, the damage initiates at the inner corner, where the shape change in geometry lies.The damage firstly evolves in an inclined manner, followed by a slight change in the direction of further evolution.For better visualization, only the damaged zones are presented in Fig. 3, with the parameters in Table 1.As can be seen in the figure, both LGDM and CGDM yield mesh-insensitive damage profiles.In Fig. 5, the load-displacement relationships for LGDM and CGDM also show reasonably mesh-insensitive results.In terms of the damage zone width, LGDM manages to preserve it to some extent, while CGDM suffers from an unphysical broadening effect.In Fig. 3g-i, the LDM result shows an unwanted zigzag damage profile during the first half of the damage evolution, see also Fig. 4 for magnified view.Moreover, damage evolution with LDM is usually strongly affected or even confined by the element boundaries.This dependency on the element geometric properties renders a difference in damage evolution rate, further leading to a difference in the position of damage front (e.g. in height) at the end of the first half of damage evolution.This is especially problematic when it comes to a direction change in the damage evolution, e.g. in this example, when small deviations of damage evolution in the early stage could lead to a greater difference in the later evolution of damage profile.The mesh sensitivity of LDM can also be clearly observed in Fig. 5, where the load-displacement characteristics show stronger deviations.Worse still, the numerical convergence is quite

Biaxial test of a cuboid concrete block
In [30], the biaxial strength envelop obtained from a microplane damage model is compared with that from the experimental data in [8].Motivated by their investigations, this numerical example focuses on a similar biaxial loading setup, as illustrated in Fig. 6a.The block is of the size 200 mm×200 mm×60 mm.For the reason of simplicity, however, the effect of damage is neglected for this example and the yield envelop is therefore presented.
With the parameters identified in Table 2 and f uc = 39 MPa for the uniaxial compressive strength of the considered concrete material, one obtains the yield envelop in Fig. 6b, together with the experimental data.It is noteworthy, that the yield envelop is plotted with the reaction forces at the first step of plastic yielding, the effect of the plastic hardening is therefore negligible and the corresponding parameters are not given in Table 2.
As shown in Fig. 6b, the results from the numerical simulations succeed to reproduce the material response measured Fig. 5 Load-displacement relationships with the localizing model in the experiment in general.The relatively large strength in biaxial compression, the small strength in biaxial tension, as well as the smooth transition in tension-compression loading are well captured.

Failure of a notched block
As already mentioned, the cement-based concrete material shows different behavior under tension and compression.With the cap plasticity model partially capturing the characteristics, the failure mechanism also needs to be distinguished by considering the plastic equivalent strain defined in Eq. ( 21).This implies that the plastic strain from tensile and compressive loadings result in different rates of damage evolution.
It has been observed in the experimental investigations in [7], that a complex failure mechanism exists when a SHCC structure with imperfections, e.g.notches, undergoes a shear loading.In their setup, the block SHCC specimen is notched as shown in Fig. 7a.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 Fig. 8b.In the current numerical example, one attempts to reproduce the failure pattern under quasi-static load regime, with the consideration of the right half of the structure due to symmetry of geometry, boundary conditions as well as loading.The simplified numerical setup is shown in Fig. 7b.
In general, the evolution of failure consists of firstly developed tensile damage zones and a subsequent shear damage zone through the specimen.As in Fig. 8c, the tensile damage initiates in both notches, from where the structure tends to be torn apart.This process appears to be rather brittle, occurring at an early loading stage.Moreover, due to the considered boundary conditions and displacement loading from the upper side, the lower damage zone appears slightly later than the upper one.The shear damage starts to grow once the tensile damage evolution stabilizes, see Fig. 8d.The shear damage initiates, from the interior region between the two notches instead of the notch edges, and subsequently goes through the specimen, leading to a sharp decrease of the reaction force and then a structural failure.It appears that the proposed model has successfully reproduced the failure pattern in principle, with the parameters in Table 3.
Further, the two parameters of the localizing model, i.e. the length-scale residual percentage R and the length-scale     deceasing rate n, are investigated in the context of the notched block failure.Figure 10 shows the damage profile with different percentages for the residual length-scale.Figure 10a represents the case with a constant length-scale parameter.
The sufficiently large gradient effect, intended to preserve the well posedness of the PDE system, would prevent the initial tensile damage to propagate much further.Instead, the damage in these regions becomes wider in distribution but smaller in value.The shear damage that appears later would also suffer from the widening in distribution, even though the constant length scale is relatively small in this case, i.e. l = √ c ≈ 0.63 mm -approximately three times of the element size h e = 0.2 mm.With decreasing residual percentage, more brittle results are observed, see Fig. 10b-d as well as the load-displacement characteristics in Fig. 12a.There is no softening in the load-displacement relationship while the tensile damage is still evolving, which means that the load bearing capacity remains still at a high level.Therefore, the further evolution of tensile damage, with a smaller R, would lead to a later appearance of shear damage and, thus, a higher peak load.The comparison between Fig. 10a,  b-d clearly explains the necessity of adopting the localizing gradient-enhanced approach.Next, the comparison between different values of the decreasing order parameter n is conducted, with the residual R = 20%.In general, a positive n leads to a large decreasing rate at an early stage of damage, while a negative n would have an opposite effect.It can be easily understood from Fig. 11, that a negative n would preserve a high level of gradient effect until the damage becomes significant, therefore, the damage zone is likely to widen.A positive n, again, leads to a more brittle behavior, see also Fig. 12b.
It is worth to mention that several research works have already discussed similar setups where tensile failure and shear failure both have taken place.For example, in [17] a brittle damage model has been proposed with the consideration of the modified Ottosen equivalent strain and further employed to model the damage evolution in a Double Edge Notch (DEN) specimen.With an initial evolution of the tensile damage, the subsequent shear damage localizes at the exterior boundary where the stress concentration lies, see Fig. 8d in [17].However, the further evolution of the shear damage is not provided in their work, and, therefore, it is not clear regarding the capability of their model in properly capturing the final failure pattern.In general, the modeling of mixed-mode failure, i.e. the coexistence of tensile fail-Fig.13 Geometric parameters (in mm) and loading conditions of the four point bending test of the RC beam ure and shear failure, could benefit from the appropriate definitions of separate criteria for tensile and shear failure evolution.While this is not so straightforward in brittle damage modeling, the three-surface cap plasticity model coupled with damage naturally provides this possibility with its different thresholds for plasticity initiation in different loading scenarios, including tension-tension, tension-compression and compression-compression as elaborated in the previous numerical example.In this model, the damage initiation is always expected to take place where plastic deformation has already occurred.This is especially meaningful when the J 2 component is considered in the plastic yield function, since this could lead to a shear band regardless of the effects of the boundary, e.g.notches.Such a consideration would lead to a final shear damage zone, which initiates inside the geometry and propagates toward the exterior boundary in a numerically stable manner.
For the simulation results shown in Fig. 8, quadratic numerical convergence has been observed for all the loading steps, see Fig. 9 for the three representative phases indicated in Fig. 8a.It is noteworthy, that for shear damage evolution accompanied by the drastic decrease of reaction force, i.e. from phase 2 to phase 3 in Fig. 8a, it is necessary to reduce the load increment in order to obtain satisfactory convergence rates.

Four point bending of a reinforced concrete (RC) beam
The final example shows the capability of the proposed model in describing the failure of a reinforced concrete beam in a four point bending test.The matrix material is represented by the proposed approach with 3D hexahedral elements and the reinforcement is modeled by considering linear elastic truss elements, whose nodes are perfectly tied to the block elements.In other words, the debonding failure of the reinforcement-matrix interface as well as the nonlinear material properties of the reinforcement are neglected for the reason of simplicity.Moreover, a special treatment is required for the truss elements so that they would only bear tensile load, given the fact that the rebar reinforcement, especially the hanger bars, can easily buckle when under compression.The geometry of the simply supported beam is illustrated in Fig. 13.The displacement-type loading is applied at two positions on the upper edge.Given the symmetry in geometry, boundary condition as well as loading conditions in the thickness direction, only half of the complete structure is considered, rendering 7711 elements in the FE simulation.
With parameters identified in Table 4 for the concrete matrix, and E b = 105 GPa and d b = 10 mm for the Young's modulus and diameter of the reinforcement respectively, the results are obtained as shown in Fig. 14a.After an early initiation of damage at the supports, subsequent damage zones occur at the lower edge of the beam, almost coinciding with the loading positions (in longitudinal direction).With further loading, two more damage zones appear close to mid-span, where the largest deflection lies, in an almost symmetric manner.Together with the further evolution of the existing damage zones, more minor damage zones follow, finally rendering the damage distribution shown in Fig. 14d.By comparison of the simulation result with the experimental investigations shown in Fig. 14b, it could be concluded that the failure pattern is successfully captured in general.As for the gradient damage model with a constant length-scale parameter for regularization, the resultant damage distribution from the work [29] is shown in Fig. 14c.The critical region of the beam, i.e. the region with a relatively large deformation, is reasonably well captured.However, the failure pattern, i.e. the multiple cracking, is lost due to the very nature of regularization by the gradient approach.
Additionally, the load-deflection relationship helps to investigate the four point bending test from a different perspective.As shown in Fig. 15, the reaction force at the loaded positions initially increases rapidly, representing the quasi-elastic phase.It it after the plastic yielding point is reached, that the damage evolution begins at the lower edge of the beam, and significant damage evolution is subsequently observed in the structure.In Fig. 15, the four representative phases previously discussed in Fig. 14a are marked by dark green points.It could be clearly seen that, when the concrete matrix surrounding the reinforcement is sufficiently degraded, the otherwise expected softening phase on the structural level is eliminated, because the reinforcement starts to bear more load together with the matrix failure.However, the slope of the load-deflection characteristic at this stage is much smaller than the initial one, given the smaller portion of the reinforcement compared to the bulk material.This would in general increase the deformation capacity of the structure, and also introduce more damage zones.Apart from the aforementioned failure pattern investigation, the load-deflection property from the simulation also agrees well with the experimental data.

Conclusion
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 in general.The difference in plastic evolution under tension and compression is addressed by the cap plasticity model, which would further lead to different points of dam-age initiation.In order to avoid the spurious damage zone expansion at a completely failed state of the structure and to preserve the well-posedness of the partial differential equation system that describes the plasticity-damage coupled problem, the damage behavior is modeled by a localizing gradient-enhanced approach, where the nonlocal counterpart of the plastic equivalent strain is considered to govern to the damage evolution.The work at hand is consistently formulated, derived and further numerically implemented into the framework of the finite element method.Several representative examples have subsequently served to illustrate the capability of the coupled model, in rendering the damage zone of physical size, in capturing the complex behavior of concrete materials in different scenarios of application, and in the enhancement of the stability of numerical solutions.Certainly, further investigations would be required for the parameter calibration, so as to better describe the experimental observations.Moreover, the extension to finite strain kinematics is one of the next priorities, which would allow the damage modeling at large deformations that are quite often observed in the experimental investigations of reinforced concrete materials including SHCC.

P
− N tr ⊗ N tr .

Fig. 1 Fig. 2
Fig. 1 Setup for the test of the L-shaped structure, a geometry and loading conditions and b-d finite element meshes

Fig. 3
Fig. 3 Damage profiles of the L-shaped structure: a-c with the localizing gradient damage model, d-f with the conventional gradient damage model and g-i with the local damage model

Fig. 4
Fig. 4 Magnified damage profiles of the L-shaped structure with the local damage model

Fig. 6
Fig.6 Biaxial test: a test setup and b the yield envelop from simulation compared to experimental results in[8]

Fig. 7 Fig. 8
Fig. 7 Setup for failure test of the block with notches: a experimental geometry [7] and b simplification for the simulations

Fig. 10 Fig. 11 Fig. 12
Fig. 10 Damage profiles with n = 1 and different residual parameter R for g(d)

Fig. 14 Fig. 15
Fig. 14 Failure in the RC beam: a damage evolution, b experiment [20], c reference simulation [29] and d our simulation result

Table 1
Material parameters for the L-shaped structure test

Table 2
Material parameters for

Table 3
Material parameters for the failure of a notched block

Table 4
parameters for the four point bending test