Strain Localization of Elastic-Damaging Frictional-Cohesive Materials: Analytical Results and Numerical Verification

Damage-induced strain softening is of vital importance for the modeling of localized failure in frictional-cohesive materials. This paper addresses strain localization of damaging solids and the resulting consistent frictional-cohesive crack models. As a supplement to the framework recently established for stress-based continuum material models in rate form (Wu and Cervera 2015, 2016), several classical strain-based damage models, expressed usually in total and secant format, are considered. Upon strain localization of such damaging solids, Maxwell’s kinematics of a strong (or regularized) discontinuity has to be reproduced by the inelastic damage strains, which are defined by a bounded characteristic tensor and an unbounded scalar related to the damage variable. This kinematic constraint yields a set of nonlinear equations from which the discontinuity orientation and damage-type localized cohesive relations can be derived. It is found that for the “Simó and Ju 1987” isotropic damage model, the localization angles and the resulting cohesive model heavily depend on lateral deformations usually ignored in classical crack models for quasi-brittle solids. To remedy this inconsistency, a modified damage model is proposed. Its strain localization analysis naturally results in a consistent frictional-cohesive crack model of damage type, which can be regularized as a classical smeared crack model. The analytical results are numerically verified by the recently-proposed mixed stabilized finite element method, regarding a singly-perforated plate under uniaxial tension. Remarkably, for all of the damage models discussed in this work, the numerically-obtained localization angles agree almost exactly with the closed-form results. This agreement, on the one hand, consolidates the strain localization analysis based on Maxwell’s kinematics and, on the other hand, illustrates versatility of the mixed stabilized finite element method.


Introduction
It is well known that under certain circumstances, frictional-cohesive materials with a softening regime exhibit strain localization prior to the occurrence of macroscopic failure. For instance, in the so-called frictional J 2 (von Mises) materials, the shear (or slip) strains usually concentrate, leading to the formation of shear bands with a small, but finite width or slip lines with a vanishing width. Similarly, cohesive crack bands or surfaces are often observed in frictional-cohesive geomaterials like concrete and rocks. Once such strain localization occurs, inside and outside these domains with highly localized deformations, the strain fields are either discontinuous due to the continuous, but plastic-damage model with general (e.g., Rankine, von Mises, Mohr-Coulomb, Drucker-Prager and more complex elliptic, parabolic, hyperbolic, etc.) failure criteria. Both general 3D and 2D (plane stress and plane strain) cases were considered. Not only the discontinuity orientation, but also the corresponding cohesive zone model, i.e., constitutive relations, evolution equations, traction-based failure criterion, softening functions, etc., are determined consistently from the given stress-based counterpart. Furthermore, the bi-directional connections and in particular the equivalence conditions between two complementary methodologies for the modeling of localized failure in quasi-brittle solids, i.e., traction-based discontinuities localized in an elastic solid and strain localization of a stress-based inelastic softening solid, have also been fully established. Numerical results obtained from the stress accurate stabilized mixed elements [34,37] coincide with the theoretical predictions, validating the developed framework.
On the one hand, so far, Maxwell's kinematics-based strain localization has been employed mainly to analyze stress-based elastoplastic damage models with the inelastic strain expressed in rate form. Consequently, in general load scenarios, the normal strains not acting on the discontinuity surface, caused by Poisson's lateral effects, can vanish in such a way that Maxwell's kinematics of strong (or regularized) discontinuities is accommodated upon strain localization. On the other hand, the continuum damage model, usually expressed in the strain-based total form, has also been widely accepted as an alternative to deal with complex material behavior. The application includes, but is not limited to, creep, fatigue and other nonlinear behavior of frictional-cohesive materials; see [38] for a review. The reason for its popularity is as much the intrinsic simplicity and versatility of the theory, as well as its consistency based on the irreversible thermodynamics with the internal variable. However, it is not so straightforward to extend the aforesaid Maxwell's kinematics-based strain localization to strain-based continuum models, since the resulting inelastic (damage) strain depends heavily on the lateral deformations induced by Poisson's ratio. This fact is inconsistent with the frictional-cohesive zone models, which generally neglect the strain and stress triaxiality. As will be shown, owing to this conceptual contradiction, strain localization of a general strain-based isotropic damage model cannot always be guaranteed as of stress-based ones.
This paper addresses the above challenging issue. Its objectives are four-fold: (i) to present strain localization analysis of strain-based isotropic damage models in a unified manner; (ii) to analyze strain localization of two classical damage models and, in particular, to derive closed-form results of the discontinuity orientation and corresponding localized frictional-cohesive model; (iii) to advocate a new modified damage model for frictional-cohesive materials, with its localized counterpart consistent with the assumption of classical discrete and smeared crack models; and (iv) to numerically validate the analytical results by the mixed stabilized finite element method [39,40].
The structure of this paper is outlined as follows. After this Introduction, Maxwell's kinematics-based strain localization analyses of strain-based isotropic damage models are presented in Section 2 in a unified approach. In Section 3, the classical J 2 and Simó and Ju [41] damage models, as well as a new modified Simó and Ju [41] model are analytically investigated. In particular, the closed-form results of the localization angles and the corresponding localized frictional-cohesive models are derived consistently. Section 4 addresses the numerical verification of the proposed strain localization analysis, regarding numerical simulations of a singly-perforated plate under uniaxial tension. The effects of Poisson's ration on the localization angles are investigated. The most relevant conclusions are drawn in Section 5. Finally, two appendices are attached to close this paper. Notation: Compact tensor notation is used in this paper. As a general rule, scalars are denoted by italic light-face Greek or Latin letters (e.g., a or λ); vectors and second-order tensors are signified by italic boldface minuscule and majuscule letters like a and A, respectively. Fourth-order tensors are identified by blackboard-bold majuscule characters (e.g., A). Symbols I and I represent the second-order and symmetric fourth-order identity tensors, respectively. Superscripts ' T ' and ' sym ' indicate the transposition and symmetrization operations, respectively. The inner products with single and double contractions are denoted by '·' and ':', respectively. The dyadic product ' ' and the symmetrized Kronecker product are defined as:

Continuum Damage Models
For a damaging solid, the constitutive relation between the second-order stress σ and (infinitesimal) strain is expressed as: where E and C := E −1 represent the fourth-order (variable) secant stiffness and compliance tensors, respectively. As depicted in Figure, the strain tensor can be kinematically decomposed into the following additive form: where the elastic strain e and the damage one d , both being recoverable upon unloading, are given by: for the damage compliance tensor C d := C − C 0 defined as the difference of the total secant one C with respect to the undamaged elasticity one C 0 . The above kinematic decomposition is depicted in Figure 1 for the 1D case. In the above constitutive relations, the damage compliance C d (or, equivalently, the secant compliance C or stiffness E) is an internal variable. Based on different approximations, either isotropic, orthotropic or fully-anisotropic evolution laws can be postulated for the fourth-order damage compliance C d ; see [3,11,36,[42][43][44] for more details. Among them, due to the conceptual simplicity, the strain-based 1 − d damage model, with d ∈ [0, 1] representing the single damage scalar index, has been widely adopted in the literature. For such a model, as will be shown later, the damage compliance can be expressed in the following total format: where ω := d/(1 − d) is an alternative damage variable in the range [0, ∞]; the fourth-order tensor C d and second-order one Λ :=C d : σ, both being model-dependent and bounded, characterize the damage compliance C d and strain d , respectively. With the above definitions, only the evolution law for the damage scalar d, rather than for the complicate tensor C d , needs to be postulated. Without loss of generality, this can be expressed as: where the stress-like internal variableq(r) is a monotonically-decreasing function of the history variable (threshold) r, with identical initial value q 0 = r 0 . For instance, the linear function: or the exponential one:q is frequently adopted, where the parameter H < 0 controls the softening functionq(r); see Figure 2. Let us introduce the secant slope H s of the q versus r − r 0 softening curve, so that: Materials 2017, 10, 434 6 of 27 Note that the secant value H s < 0 is not coincident with the softening parameter H introduced in Equation (6); only for the linear softening case, the relation H = H s holds; see Figure 2 for the illustration.
Accordingly, the alternative damage variable ω can be expressed as: which results in the following damage strain tensor d : Note that the entity (1 − q 0 /q), related to the stress-like internal variable q, and the characteristic tensor Λ, dependent on the stress σ, are both bounded.

Strain Localization and Localized Damage Models
For strain localization to occur in a softening solid and to develop eventually into a fully-softened discontinuity at the final stage of the deformation process, material points inside the strain localization band undergo inelastic loading, while those outside it unload elastically [20,21,34]. As recently clarified by the authors [45], this standpoint is equivalent to strong discontinuities of the displacement field localized in an elastic solid. Based on this equivalence, novel strain localization analysis [15,36] based on the satisfaction of Maxwell's kinematics [32] and of stress continuity has been proposed. Particularly, for a stress-based material model with a softening regime, not only the discontinuity (band) orientation, but also the corresponding traction-based cohesive model can be consistently derived. In this section, the aforesaid method is extended to the strain-based damage models addressed in Section 2.1.  Let us first consider an elastic solid Ω containing a strong discontinuity S as depicted in Figure 3a. The orientation of the discontinuity S is characterized by the normal vector n, as well as two tangential vectors m and p, constituting a local coordinate system (n, m, p). Upon strain localization, the standard kinematics of a continuum medium are replaced by Maxwell's kinematics [32], which allows for jumps in the derivatives of the displacement field with respect to the direction normal to the discontinuity S, but not in the derivatives with respect to the direction tangential to it. According to the above Maxwell kinematics, the jump in the strain field between the inside and the outside of the localization band may be expressed exclusively in terms of the unit normal vector and a deformation vector. Therefore, the singular strain field caused by the displacement jump u = w across the discontinuity S is expressed as: where e represents the bulk strain field, being elastic but not necessarily continuous across the discontinuity S; the Dirac-delta δ S is introduced to characterize the singular strain field at the interface S; see Figure 3b.
Recalling the aforementioned equivalence, for strain localization in damaging solids to occur, it follows from Maxwell's kinematic (10) that: As the stress σ and the resulting characteristic tensor Λ are both bounded, upon strain localization, the damage compliance C d and the variable ω in Equation (9) have to be singular, i.e., withC d andω both being bounded. Calling for the relation (8), the singularity of ω implies the existence of a bounded softening parameterH < 0, such that: It then follows that: Therefore, the displacement jump w can be solved in terms of the characteristic tensor Λ as: where the local components (w n , w m , w p ) of the displacement jump vector w in the orthogonal system (n, m, p) of the discontinuity S are given by: Substitution of the above results into the relation (14) yields: for the discontinuity angles θ cr upon which the kinematic conditions (14) are satisfied. Calling for the relations of the local components (Λ mm , Λ pp , Λ mp ) between the principal values Λ i (i = 1, 2, 3) and a set of characteristic angles θ, the kinematic constraints (17) yield a system of nonlinear equations, such that the discontinuity angles θ cr can be determined.
In particular, let us consider the 2D plane stress and plane strain conditions (σ np = σ mp = 0). In such cases, the discontinuity orientation can be characterized by the inclination angle (anti-clockwise) θ cr ∈ [−π/2, π/2] between the normal vector n of the discontinuity and the principal vector v 1 of the tensor Λ; see Figure 4. Furthermore, the in-plane components (Λ nn , Λ mm , Λ nm ) are given by: for the in-plane principal values Λ 1 and Λ 2 . As the out-of-plane constraint Λ mp (θ cr ) = 0 is automatically fulfilled, it follows from the condition Λ mm (θ cr ) = 0 that: where the in-plane principle values Λ 1 and Λ 2 of the tensor Λ, satisfying Λ 1 ≥ 0 and Λ 2 ≤ 0, are constrained by the conditions: Note that in the plane stress state (σ pp = 0), the constraint Λ pp (θ cr ) = 0 makes no sense, and the cohesive relations (16a) always hold; see Remark 2 for more discussion.

Remark 1.
The above arguments also apply to a regularized discontinuity B with finite bandwidth b 0. Note that the bandwidth b is a regularized parameter, which can be taken as small as desired. In such a case, the Kronecker-delta δ S is regularized by [15,36].

Remark 2.
To gain further insight into the above strain localization analysis, let us consider a regularized discontinuity B in the solid Ω with the following internal virtual work: Materials 2017, 10, 434

of 27
For general 3D and plane strain cases, provided that the discontinuity angles θ cr satisfying the constraints (17) exist, the inelastic internal work localized in the discontinuity band B is expressed as: where the results (16) with the regularization relation ω =ω/b have been recalled. This is exactly the result for the elastic solid containing a frictional-cohesive crack. In the plane stress state (σ pp = 0), only the condition (17a) is necessary for guaranteeing the identity: with the constraint Λ pp (θ cr ) = 0 being irrelevant. That is, regardless of the constraint Λ pp (θ cr ) = 0, strain localization also occurs in the plane stress condition, with the resulting discontinuity characterized by the cohesive relations (16a).

Closed-Form Results
In this section, the strain localization analysis presented in Section 2.2 is applied to several classical 1 − d damage models. Both the discontinuity orientation and the corresponding localized constitutive laws are given.

J 2 Damage Model
Let us first discuss the J 2 damage model. It is usually employed for the modeling of shear bands or slip lines in J 2 softening materials; see [46] for the details.

Constitutive Relations
In the J 2 damage model, the stress and strain tensors are decomposed as: where p = 1 3 trσ and s = σ − pI are the volumetric and deviatoric parts of the stress tensor σ, respectively; v = tr and e := − 1 3 v I represent the trace and the deviatoric part of the strain tensor , respectively.
While the volumetric behavior remains elastic, a scalar variable d ∈ [0, 1] is introduced to characterize the deviatoric behavior. It then follows that: wheres = 2G 0 e is the effective deviatoric stress tensor; K 0 and G 0 are the elastic bulk and shear moduli of the material, respectively. In this case, the constitutive relation (2) is particularized as: for the volumetric and deviatoric projection operators: The damage compliance tensor C d and the resulting inelastic strain d are then given by: for the variable ω := d/(1 − d) and the characteristic tensor Λ = s/(2G 0 ). The above J 2 damage model corresponds to the following Helmholtz free energy ψ: together with the energy dissipation inequality: for the energy release rate Y, conjugate to the damage variable d.
Additionally, a strain-based (or effective stress-based) damage criterion is introduced: in terms of the equivalent strain eq ( ) and the threshold r: for the second invariantJ 2 = 1 2s :s of the effective deviatoric stress tensors = 2G 0 e. Upon damage loading, it follows from the consistency conditionḞ ( , r) = 0 that the threshold r, with the initial value r 0 , represents the maximum value of eq ever reached.
In order to postulate the damage evolution law in the format (5), it is convenient to rewrite the above strain-based damage criterion as: in terms of the following equivalent stress σ eq and the corresponding softening variable q(r): where J 2 = 1 2 s : s represents the second invariant of the deviatoric stress tensor s.

Orientation of the Discontinuity
For the J 2 damage model with the characteristic tensor Λ = s/(2G 0 ), the conditions (17) upon strain localization become: for the deviatoric stress components (s mm , s pp , s mp ) on the discontinuity surface. Accordingly, the orientation angles θ cr can be solved in terms of the principal values s i (i = 1, 2, 3) of the deviatoric stress tensor s.
In particular, for the 2D cases of plane stress and plane strain (σ np = σ mp = 0), the condition s mm (θ cr ) = 0 gives: Furthermore, it follows from the constraint s 3 = s pp = 0 that: That is, in 2D conditions strain localization of a J 2 damaging material with stress continuity can always occur along the discontinuity angle θ cr = ±45 • . Note that the above analytical result coincides with the numerical simulations [46].

Localized Damage Model
Provided that the solution to the constraints (35) exists, it follows from the vanishing trace tr s = 0 that s nn = − s mm + s pp = 0 holds on the discontinuity surface with the normal n(θ cr ). Accordingly, the relation (14) becomes: or, equivalently, The above frictional-cohesive relations can also be given straightforwardly from Equation (16): where the identities s nm = σ nm = t m and s np = σ np = t p have been considered.
With the frictional-cohesive relations (39) upon strain localization, the stress-based damage criterion (33) becomes: or the alternative displacement jump-based one: where w eq and t eq denote the equivalent displacement jump and traction, defined as: Similarly to the strain-based threshold r, the displacement jump-like internal variable κ represents the maximum value of the equivalent displacement jump w eq ever reached.
As expected, upon strain localization, the J 2 damage model is characterized by displacement jumps and tractions. That is, a pure mode-II discontinuity in the sense of fracture mechanics occurs in the J 2 damaging material.

Simó and Ju [41] Damage Model
Next, let us consider the Simó and Ju [41] damage model, which has been widely adopted in the literature. As will be shown, strain localization with stress continuity cannot always occur in any case.

Constitutive Relations
The Simó and Ju [41] damage model is characterized by the following well-defined Helmholtz free energy: where E 0 = λ 0 I I + 2µ 0 I represents the fourth-order material elasticity tensor, with λ 0 and µ 0 being the Lamé constants. Standard arguments then yield: or, equivalently, Similarly, the damage compliance tensor C d and the resulting damage strain d are characterized in the form (4): for the variable ω := d/(1 − d) and the characteristic tensor Λ = C 0 : σ.
The above constitutive relations are constrained by the damage energy dissipation inequality: It is then possible to introduce the Beltrami strain-based damage criterion: where the equivalent strain eq and the threshold r are defined as: for the longitudinal or constrained modulus M 0 = λ 0 + 2µ 0 . Alternatively, for the damage evolution law (5), a stress-based failure criterion can be expressed as: for the equivalent stress σ eq and the corresponding stress-like internal variable q: Once the softening function q(r) is postulated, the damage variable d can be determined as shown later.

Orientation of the Discontinuity
For the Simó and Ju [41] damage model, the characteristic tensor Λ is expressed as: so that the constraints (17) give: With the stress components expressed in terms of the principal stresses σ i (i = 1, 2, 3) and a set of characteristic angles θ, the discontinuity angles θ cr can be determined.
In particular, let us focus on the 2D plane stress and plane strain cases in which the identity σ np = σ mp = 0 holds. For the plane stress state, the conditions (19) and (20a) give: For the plane strain case (i.e., σ 3 = σ pp = 0), it follows from the conditions (19) and (20a) that: Note that the convention σ 1 > σ 2 is assumed as usual. As can be seen, both results depend on Poisson's ratio ν 0 and coincide for a vanishing one.

Localized Damage Model
Provided that the solution to the conditions (53) exists, it follows from the relation (14) that: Accordingly, the frictional-cohesive relations are given by: or, in the component form, where the second-order elastic acoustic tensor Q 0 and its inverse Q −1 0 are expressed as: in the local (n, m, p) coordinate system. With the relation (55), upon strain localization, the stress-based damage criterion (50) becomes: or, the alternative form, where the equivalent displacement jump w eq and traction t eq are expressed as: for the coefficient β 0 = √ G 0 /M 0 . The above localized cohesive law for the Simó and Ju [41] damage model was first derived in [22,24] where the discontinuity orientation was determined from the classical discontinuous bifurcation analysis rather than the Maxwell kinematic condition considered in this work.

Modified Damage Model
As can be seen from the cohesive relation (56b) 1 , the normal behavior upon strain localization is affected by the lateral deformation due to the non-vanishing Poisson's ratio. This fact is inconsistent with the concept of a cohesive zone model in which the lateral deformations on the discontinuity surface are not considered. To reconcile this inconsistency, let us consider a modified Simó and Ju damage model.

Constitutive Relations
In the modified damage model, the constitutive relations read: which gives the following damage strain tensor d : where the compliance tensor C d is expressed as the form (4): for the variable ω := d/(1 − d). Here, the diagonal matrixC 0 is the Voigt representation of a fourth-order reference compliance tensor, which is obtained from the elastic one C 0 by ignoring the lateral deformations caused by Poisson's ratio. The explicit expressions of the resulting secant compliance C and stiffness E are given in Appendix A.
As these constitutive tensors are of major symmetry, the corresponding energy function ψ is defined as: Standard arguments give the constitutive relations (61) and the following energy dissipation inequality: for the effective stress tensorσ = σ/(1 − d). Accordingly, the following strain-based damage criterion can be considered: where the equivalent strain eq and the threshold r are given by: Alternatively, the stress-based damage criterion is expressed as: for the equivalent stress σ eq and internal variable q: Note that the above equivalent strain and stress are subtly different from their counterparts in the Simó and Ju [41] damage model.

Orientation of the Discontinuity
For the damage strain (62), the constraints (17) become: With the stress components expressed in terms of the principal stresses σ i (i = 1, 2, 3) and a set of characteristic angles θ, the discontinuity angles θ cr can thus be determined.
For the plane stress and plane stress conditions, both Constraints (20a) and (20b) coincide, and the last two conditions in Equation (70) are automatically fulfilled. It then follows from the remaining one Λ mm (θ cr ) = 0 that: The above result applies for the biaxial tension-compression quadrant, i.e., σ 1 > 0 > σ 2 . For the case σ 1 > σ 2 > 0, the discontinuity angle is given from cos(2θ cr ) = 1, i.e., θ cr = 0. That is, the discontinuity orientation n coincides with the major principle vector v 1 of the stress tensor σ. Compared to the Simó and Ju [41] model, the discontinuity angle (71) does not depend on the elastic Poisson's ratio. This is consistent with the conclusion drawn for stress-based material models [15,36].

Localized Damage Model
Provided the solution to the conditions (53) exists, it follows from the relation (14) that: for the displacement jumps w := w n , w m , w p T . Accordingly, the frictional-cohesive tractions t := t n , t m , t p T are given by: or, inversely, for the reference stiffnessĒ 0 and complianceC 0 of the discontinuity: Note that in 2D plane stress and plane strain conditions, the above frictional-cohesive relations also apply with a vanishing out-of-plane jump w p = 0 (or, equivalently, t p = 0).
After calling for the relation (72) and performing some straightforward manipulations, upon strain localization, the stress-based damage criterion (50) becomes: or, equivalently, where the equivalent displacement jump w eq and traction t eq are expressed as: for the material property β 0 = √ G 0 /E 0 . Compared to the frictional-cohesive relations (56) resulting from the Simó and Ju [41] damage model, the longitudinal modulus M 0 in Equation (57) is replaced here by Young's modulus E 0 . This modification is consistent with the concept of a frictional-cohesive crack model in which only the strain components acting on the discontinuity surface are accounted for [47]. Moreover, upon the assumption of a continuous elastic strain field e (x) across the discontinuity, it can be regularized as the smeared crack model discussed in [38,48]; see Appendix B for more details. The above facts justify the above modified Simó and Ju damage model, which will be addressed in the numerical context elsewhere.

Numerical Verification
In this section, the analytical results presented in Section 3 are numerically verified. Due to the poor resolution upon strain localization, the irreversible displacement-based finite element method is not sufficient for this purpose. In order to circumvent this difficulty, the mixed stabilized strain/displacement − u finite element method, recently developed by Cervera et al. [39] for elasticity and extended to isotropic damage materials [40,49], is considered.

Mixed Stabilized Strain/Displacement Element
The aforesaid mixed stabilized strain/displacement element is briefly recalled in a secant format appropriate for the continuum damage models discussed in this work. The mechanical behavior of the solid body Ω is described by the compatibility of deformations and the equilibrium of body forces: where u is the displacement vector, is the strain tensor, σ represents the stress tensor, ∇ sym and ∇ · (·) are the adjoint symmetric gradient and the divergence operators, respectively; b * is the body force vector. The strain and stress tensors are linked by the secant constitutive relations (1) with appropriate damage evolution laws. The strong form of the boundary value problem is completed by imposing proper traction and displacement boundary conditions on ∂Ω.
After pre-multiplying the secant stiffness tensor E, the strong form reads: for the unknowns fields of total strains and displacements u. Introducing the test function γ ∈ G ⊂ L 2 (Ω) dim for strains and v ∈ V ⊂ H 1 (Ω) dim for displacements and applying Gauss's divergence theorem to the strong form yields the following weak form of the mixed problem: where F(v) represents the work done by boundary tractions t * and body forces b * . Note that this weak form is symmetric. The discrete version of the mixed weak form is obtained by substituting the unknown fields with their finite element interpolation counterparts: where h and u h are the nodal degrees of freedom; γ h and v h are the discrete test functions for the strains and displacements pertaining to the spaces G h and V h , respectively, the discrete counterparts of G and V . As implied by the Inf-Supcondition [50], equal interpolations for strains and displacements are bound to be unstable. To overcome this issue, a stabilization procedure needs be introduced by modifying the discrete variational form with numerical stability while maintaining consistency. For the stabilized problem [39,40] derived the following system of equations: where the stabilized discrete strain field˜ h is a blending of the continuous and discontinuous strain fields ( h , ∇ sym u h ) weighted by the stabilization parameter τ : for the arbitrary positive parameter c , the representative mesh size h and the characteristic length L 0 of the problem. Note that the above stabilized formulation is consistent with the original discrete weak form since, with converging values of the unknowns h and u h , the contribution of the stabilization terms (those multiplied by τ ) vanishes.

Numerical Results
The above mixed stabilized strain-displacement element is then applied to analyses of strain localization. As the J 2 damage model has already been thoroughly verified in [46], only the Simó and Ju [41] damage model and the modified one are presented here.
The example is a 2D singly perforated strip loaded in uniaxial tension by stretching via imposed vertical displacements at the top and bottom ends; horizontal movement is not restrained. Both plane strain and plane stress conditions are considered and compared. Figure 5a depicts the geometry of the problem with dimensions 20 m × 40 m ×1 m (width × height × thickness). An imperfection is introduced with a slanted perforation of diameter D = 1 m such that symmetric solutions are excluded.

Singly perforated strip
Uniaxial stretch In all simulations, the exponential softening curve (6b) is considered and regularized according to the crack band theory [12], with the softening parameter H determined as: where l ch and L H := E 0 G F /q 2 0 represent the localization and Griffith's characteristic lengths, respectively. In the mixed formulation, l ch = 2h, with h being the finite element mesh size.
The following material properties are assumed: Young's modulus E 0 = 10 MPa, uniaxial failure stress q 0 = E/1000 = 10 KPa and fracture energy is G F = 500 J/m 2 . Poisson's ratio ν 0 takes different values for comparison purposes.
The mesh used in the analyses consists of 45,750 mixed P1P1triangles (23,183 nodes) as shown in Figure 5b, with an average mesh size of h = 0.20 m. This level of refinement ensures convergence of the numerical simulation and allows comparison with the analytical solution. The mesh is fully unstructured. All analyses are performed on the same mesh.
Loading is applied by imposed vertical displacements at both ends of the strip. The Newton-Raphson method is used to solve the nonlinear system of equations arising from the spatial and temporal discretization of the problem. An automatic procedure is used to decide the step size, and about 200 steps are necessary to complete the analyses. Convergence of a time step is attained when the ratio between the norms of the residual and the total forces is lower than 10 −3 .
Calculations are performed with an enhanced version of the finite element program COMET [51], developed by the authors at the International Center for Numerical Methods in Engineering (CIMNE). Pre-and post-processing is done with GiD, also developed at CIMNE [52].

Simó and Ju [41] Damage Model
Let us first discuss the Simó and Ju [41] damage model. For the considered uniaxial stress state (σ 1 > 0 and σ 2 = 0), it follows from the closed-form results (54) that: The resulting discontinuity angles θ cr for various Poisson's ratio ν 0 are summarized in Table 1. Figures 6 and 7 show the numerical contours of vertical displacements with respect to Simó and Ju [41] damage model in the plane stress and plane strain conditions. Strain localization is manifested by the discontinuous displacement field. The numerical discontinuity angles θ cr are then determined as in Table 1   As can be seen, the analytical results agree almost exactly with the numerical ones, verifying the validity of both the novel strain localization analysis presented in this work and the mixed stabilized finite element method proposed in [39,40]. Slight differences are attributed to the unavoidable perturbation of the hole and to the boundary effects.
For illustration, the evolution of damage and vertical displacements for the Simó and Ju [41] model in the case of plane strain with Poisson's ratio ν 0 = 0.30 is shown in Figure 8. Furthermore, load-displacement curves for the Simó and Ju [41] model in the conditions of plane strain and plane stress with various Poisson's ratio ν 0 = 0.0, 0.15, 0.30 and 0.45 are given in Figure 9.

Modified Simó and Ju Damage Model
For the modified Simó and Ju damage model, the discontinuity angle θ cr given by the analytical results (71) is independent of Poisson's ratio ν 0 . In particular, for the uniaxial stress state, it follows that: for both plane stress and plane strain states. Figure 10 shows the numerical contours of vertical displacements obtained from the modified Simó and Ju [41] damage model in the plane stress and plane strain conditions. As predicted from the analytical result (86), the displacement field exhibits a jump across the horizontal discontinuity surface, normal to the principal major stress, i.e., θ cr = 0. The novel strain localization analysis and the mixed stabilized finite element method are verified again.

Conclusions
This paper presents Maxwell's kinematics-based strain localization analysis of frictional-cohesive materials characterized by strain-based damage models of total form, extending our recent work for stress-based ones of rate type. In such models, the (inelastic) damage strains are characterized by a bounded characteristic tensor and an unbounded variable related to the damage scalar. For strain localization to occur, Maxwell's kinematics of a strong (or regularized) discontinuity has to be reproduced by the inelastic damage strains. This kinematic constraint yields a set of nonlinear equations from which the discontinuity orientation and localized frictional-cohesive relations of damage type are derived consistently.
Two classical isotropic damage models, i.e., the J 2 model and the Simó and Ju [41] one, are then analyzed. In particular, upon strain localization, the J 2 damaging solid is characterized by a pure mode-II discontinuity in which only the tangential tractions (or displacement jumps) are activated. Comparatively, a mixed-mode discontinuity occurs in the Simó and Ju [41] isotropic damaging solid. Furthermore, the localization angles and the resulting cohesive model depend explicitly on Poisson's ratio. That is, the lateral deformations, not acting on the discontinuity surface, affect heavily strain localization, inconsistent with classical frictional-crack models in which the strain and stress triaxiality is not accounted for.
To remove the above inconsistency, we further proposed a modified Simó and Ju damage model. Its strain localization analysis naturally yields a consistent damage-type cohesive crack model, with the discontinuity orientation independent of the material elastic property (i.e., Poisson's ratio). It is found that, on the one hand, the discontinuous version of the cohesive crack model is formally identical, but not equivalent to that first proposed in [7] and widely applied in the literature. On the other hand, the regularized version recovers the classical smeared crack model [47,48,53] upon the assumption of a continuous stress field across the discontinuity band.
The analytical results are numerically-verified by the mixed stabilized finite element method, regarding a singly perforated plate under uniaxial tension. Remarkably, for all of the damage models discussed in this work, the numerically-obtained localization angles agree almost exactly with the closed-form results. This agreement, on the one hand, consolidates the strain localization analysis based on Maxwell's kinematics and, on the other hand, illustrates versatility of the mixed stabilized finite element method.
Finally, so far, only the material models with associated evolution laws have been considered in Maxwell's kinematics-based strain localization analysis. The extension to non-associated cases will be the forthcoming work. b 0; see Figure B1a. For such a regularized discontinuity, an inelastic deformation vector e := w/b can be introduced such that: or, equivalently, e = ωC 0 · t =C S · t,C S = ωC 0 (A4b) for the damage variable ω := d/(1 − d) =ω/b and the complianceC S of the regularized discontinuity. In both the frictional-cohesive relations and the regularized counterparts, the kinematic unknowns related to the discontinuity, i.e., the displacement jump w or the inelastic deformation vector e, are determined through the statement of the traction continuity condition, i.e., t = σ + S · n = σ − S · n (A5) where σ + S and σ − S represent the bulk stresses "ahead of" and "behind" the discontinuity S. Alternatively, as depicted in Figure B1b, if the stresses are assumed to be continuous across the discontinuity, i.e., σ + S = σ − S = σ, traction continuity t = σ · n follows automatically, and no specific kinematic unknowns for the discontinuity are necessary. In this case, it follows from the regularized cohesive relations (A4) that: d = e n sym = C S · t n sym = C d : σ (A6) for the fourth-order damage compliance C d : The resulting stress versus strain relation become: where the compliance tensor C is expressed as: for the second-order geometric tensor N = n n. In matrix notation, the compliance C of the smeared crack model is given by: which applies to general 3D and plane strain conditions. In the case of plane stress, it follows that: As can be seen, upon strain localization, the damage affects only the behavior normal to the discontinuity surface, while the one parallel to it remains elastic. This idea dates back to the classical smeared crack models [53][54][55][56]; see also the reviews [38,48].