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

Abstract

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.

1. 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 non-smooth displacements, or even singular (unbounded), caused by the discontinuous displacements. These strain (or weak) and displacement (or strong) discontinuities are usually regarded as prognostics of localized failure in frictional-cohesive materials and resulting catastrophic collapse of engineering structures. Therefore, it is of vital importance to predict the occurrence of strain localization and quantify its adverse influence on the overall response.

In the theoretical context, strain localization and the induced localized failure in frictional-cohesive materials can be characterized by either generalized continuum models or nonlinear crack/fracture models. In the former approach, the effects of strain/displacement discontinuities are smoothed or smeared. Accordingly, the overall nonlinear behavior of the weakened material can be described by a tensorial constitutive law in terms of stress vs. strain. Among many alternatives, plasticity [1] and damage mechanics [2] or their combination [3,4,5] are frequently employed to develop appropriate inelastic constitutive laws, usually based on the irreversible thermodynamics with internal variables. Comparatively, in nonlinear crack/fracture models, displacement jumps are explicitly accounted for by embedding the discontinuities into a solid matrix (bulk) along preferred orientations. It is in general assumed that the discontinuity in which energy dissipation localizes is characterized by a vectorial traction-based frictional-cohesive zone model while the bulk remains elastic, between which the traction continuity condition is imposed. Similarly, depending on the recoverable/irreversible properties of the discontinuities, frictional-cohesive zone models of either the plastic [6], damage [7,8,9] or combined plastic-damage [10,11] type can be established. Note that in both approaches, the introduction of fracture energy per discontinuity surface, usually regarded as a material property, is indispensable to guarantee the objectivity of energy dissipation during the whole failure process [12,13].

Though both the continuum and crack/fracture models are able to quantify the overall effects of localized failure in frictional-cohesive materials, their theoretical fundamentals and physical motivations are rather distinct from each other. It is of great significance to investigate and clarify the interrelations between them. On the one hand, the constitutive relations of a frictional-cohesive zone model, in terms of tractions versus displacement jumps (or the regularized inelastic deformation vector), are generally established in a heuristic or ad hoc manner. This heavily restrains developing a physically-sound and theoretically-rational cohesive zone model for mixed-mode failure. On the other hand, the smeared crack model can be regarded as a particular generalized continuum model in between both families. More specifically, it can be recovered by combining of the elastic bulk and the nonlinear cohesive discontinuity upon the assumption of a continuous stress field overall the whole solid; see [14,15]. Therefore, it would be rather enlightening if the above procedure could be performed inversely, i.e., developing a frictional-cohesive zone model in the context of discrete crack methods [16,17] from its continuum counterpart with softening material laws.

Employing the strong discontinuity approach [18,19], Oliver et al. [20,21,22,23,24] derived frictional-cohesive zone models by projecting inelastic material laws onto the discontinuity orientation. Similar work has also been done in [25,26]. However, only very simple continuum models, e.g., the classical isotropic damage model [22,23,24], the Rankine and plane strain $J_2$ (von Mises) plasticity models [20,21], can be considered, whereas more general material constitutive laws cannot be sufficiently accounted for [21]: “obtaining such explicit forms of the discrete constitutive equations is not so straight-forward for other families of elastoplastic models”. Furthermore, as the discontinuity orientation is determined from the discontinuous bifurcation condition [27,28,29,30,31,32] together with the null softening modulus, strong (regularized) discontinuities cannot form in general cases. Consequently, some kinematic mismatches are observed [24,33], resulting in spurious stress locking [34,35].

To overcome the issue of mispredicted discontinuity orientation, Cervera et al. [34] suggested directly using Maxwell’s kinematic conditions of strong (or regularized) discontinuities to determine the discontinuity orientation, so that the stress locking-free property can be guaranteed for a fully-softened discontinuity. The closed-form results of the localization angles for the $J_2$ plasticity model were validated by numerical simulations in the cases of plane stress and plane strain. More recently, the authors [15,36] successfully extended the above strain localization analysis to a unified stress-based 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 $\lambda$); vectors and second-order tensors are signified by italic boldface minuscule and majuscule letters like $\mathit{a}$ and $\mathit{A}$, respectively. Fourth-order tensors are identified by blackboard-bold majuscule characters (e.g., $\mathbb{A}$). Symbols $\mathit{I}$ and $\mathbb{I}$ represent the second-order and symmetric fourth-order identity tensors, respectively. Superscripts ‘${}^{\mathrm{T}}$’ and ‘${}^{\mathrm{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 ${\displaystyle \underline{\overline{\otimes }}}$ are defined as:

$${\displaystyle {\left(\mathit{A}\otimes \mathit{B}\right)}_{ijkl}=A_{ij}{B}_{kl},{\left(\mathit{A}\underline{\overline{\otimes }}\mathit{B}\right)}_{ijkl}=\frac{1}{2}\left(A_{ik}{B}_{jl}+A_{il}{B}_{jk}\right)}$$

2. Strain Localization in Damaging Solids

2.1. Continuum Damage Models

For a damaging solid, the constitutive relation between the second-order stress $\sigma$ and (infinitesimal) strain $ϵ$ is expressed as:

$$\begin{array}{c}\hfill \sigma =\mathbb{E}:ϵ,ϵ=\mathbb{C}:\sigma \end{array}$$

(1)

where $\mathbb{E}$ and $\mathbb{C}:={\mathbb{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:

$$\begin{array}{c}\hfill ϵ={ϵ}^{\mathrm{e}}+{ϵ}^{\mathrm{d}}\end{array}$$

(2)

where the elastic strain ${ϵ}^{\mathrm{e}}$ and the damage one ${ϵ}^{\mathrm{d}}$, both being recoverable upon unloading, are given by:

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{e}}={\mathbb{C}}_0:\sigma ,{ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma \end{array}$$

(3)

for the damage compliance tensor ${\mathbb{C}}^{\mathrm{d}}:=\mathbb{C}-{\mathbb{C}}_0$ defined as the difference of the total secant one $\mathbb{C}$ with respect to the undamaged elasticity one ${\mathbb{C}}_0$. The above kinematic decomposition is depicted in Figure 1 for the 1D case.

In the above constitutive relations, the damage compliance ${\mathbb{C}}^{\mathrm{d}}$ (or, equivalently, the secant compliance $\mathbb{C}$ or stiffness $\mathbb{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 ${\mathbb{C}}^{\mathrm{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\in [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:

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{d}}={\displaystyle \frac{d}{1-d}}{\overline{\mathbb{C}}}^{\mathrm{d}}=\omega {\overline{\mathbb{C}}}^{\mathrm{d}}⟹{ϵ}^{\mathrm{d}}={\overline{\mathbb{C}}}^{\mathrm{d}}:\sigma =\omega \Lambda \end{array}$$

(4)

where $\omega :=d/(1-d)$ is an alternative damage variable in the range $[0,\infty ]$; the fourth-order tensor ${\overline{\mathbb{C}}}^{\mathrm{d}}$ and second-order one $\Lambda :={\overline{\mathbb{C}}}^{\mathrm{d}}:\sigma$, both being model-dependent and bounded, characterize the damage compliance ${\mathbb{C}}^{\mathrm{d}}$ and strain ${ϵ}^{\mathrm{d}}$, respectively.

With the above definitions, only the evolution law for the damage scalar d, rather than for the complicate tensor ${\mathbb{C}}^{\mathrm{d}}$, needs to be postulated. Without loss of generality, this can be expressed as:

$$\begin{array}{c}\hfill d\left(r\right)=1-\frac{q}{r},q=\widehat{q}\left(r\right)=\left(1-d\right)r\end{array}$$

(5)

where the stress-like internal variable $\widehat{q}\left(r\right)$ is a monotonically-decreasing function of the history variable (threshold) r, with identical initial value $q_0=r_0$. For instance, the linear function:

$$\begin{array}{c}\hfill \widehat{q}\left(r\right)=\left\{\begin{array}{cc}\left[1+H·\left({\displaystyle \frac{r-r_0}{r_0}}\right)\right]q_0\hfill & r_0\le r\le r_u\hfill \\ 0\hfill & r\ge r_u\hfill \end{array}\right.\end{array}$$

(6a)

or the exponential one:

$$\begin{array}{c}\hfill \widehat{q}\left(r\right)=q_{0}exp\left[H·\left({\displaystyle \frac{r-r_0}{r_0}}\right)\right]\end{array}$$

(6b)

is frequently adopted, where the parameter $H<0$ controls the softening function $\widehat{q}\left(r\right)$; see Figure 2.

Let us introduce the secant slope $H_s$ of the q versus $r-r_0$ softening curve, so that:

$$\begin{array}{c}\hfill H_s:={\displaystyle \frac{q-q_0}{r-r_0}}={\displaystyle \frac{\mathcal{H}}{\mathcal{H}+1}}<0⟺r-q={\displaystyle \frac{1}{\mathcal{H}}}\left(q-q_0\right)\end{array}$$

(7)

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 $\omega$ can be expressed as:

$$\begin{array}{c}\hfill \omega :={\displaystyle \frac{d}{1-d}}={\displaystyle \frac{r-q}{q}}={\displaystyle \frac{1}{\mathcal{H}}}\left(1-{\displaystyle \frac{q_0}{q}}\right)\end{array}$$

(8)

which results in the following damage strain tensor ${ϵ}^{\mathrm{d}}$:

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma =\omega \Lambda ={\displaystyle \frac{1}{\mathcal{H}}}\left(1-{\displaystyle \frac{q_0}{q}}\right)\Lambda \end{array}$$

(9)

Note that the entity $(1-q_0/q)$, related to the stress-like internal variable q, and the characteristic tensor $\Lambda$, dependent on the stress $\sigma$, are both bounded.

2.2. 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 $\Omega$ containing a strong discontinuity $\mathcal{S}$ as depicted in Figure 3a. The orientation of the discontinuity $\mathcal{S}$ is characterized by the normal vector $\mathit{n}$, as well as two tangential vectors $\mathit{m}$ and $\mathit{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 $\mathcal{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 $〚\mathit{u}〛=\mathit{w}$ across the discontinuity $\mathcal{S}$ is expressed as:

$$\begin{array}{c}\hfill ϵ={ϵ}^{\mathrm{e}}+{\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}{\delta }_{{}_{\mathcal{S}}}\end{array}$$

(10)

where ${ϵ}^{\mathrm{e}}$ represents the bulk strain field, being elastic but not necessarily continuous across the discontinuity $\mathcal{S}$; the Dirac-delta ${\delta }_{{}_{\mathcal{S}}}$ is introduced to characterize the singular strain field at the interface $\mathcal{S}$; see Figure 3b.

Recalling the aforementioned equivalence, for strain localization in damaging solids to occur, it follows from Maxwell’s kinematic (10) that:

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma =\omega \Lambda ={\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}{\delta }_{{}_{\mathcal{S}}}\end{array}$$

(11)

As the stress $\sigma$ and the resulting characteristic tensor $\Lambda$ are both bounded, upon strain localization, the damage compliance ${\mathbb{C}}^{\mathrm{d}}$ and the variable $\omega$ in Equation (9) have to be singular, i.e.,

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{d}}={\overline{\mathbb{C}}}^{\mathrm{d}}{\delta }_{{}_{\mathcal{S}}},\omega =\overline{\omega }{\delta }_{{}_{\mathcal{S}}}\end{array}$$

(12)

with ${\overline{\mathbb{C}}}^{\mathrm{d}}$ and $\overline{\omega }$ both being bounded.

Calling for the relation (8), the singularity of $\omega$ implies the existence of a bounded softening parameter $\overline{\mathcal{H}}<0$, such that:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\mathcal{H}}}={\displaystyle \frac{1}{\overline{\mathcal{H}}}}{\delta }_{{}_{\mathcal{S}}},\overline{\omega }={\displaystyle \frac{1}{\overline{\mathcal{H}}}}\left(1-{\displaystyle \frac{q_0}{q}}\right)\end{array}$$

(13)

It then follows that:

$$\begin{array}{c}\hfill \omega \Lambda =\overline{\omega }\Lambda {\delta }_{{}_{\mathcal{S}}}={\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}{\delta }_{{}_{\mathcal{S}}}⟹\overline{\omega }\Lambda ={\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}\end{array}$$

(14)

Therefore, the displacement jump $\mathit{w}$ can be solved in terms of the characteristic tensor $\Lambda$ as:

$$\begin{array}{c}\hfill \mathit{w}=w_n\mathit{n}+w_m\mathit{m}+w_p\mathit{p}=\overline{\omega }\left(2\mathit{n}·\Lambda -\mathit{n}{\Lambda }_{nn}\right)\end{array}$$

(15)

where the local components $(w_n,w_m,w_p)$ of the displacement jump vector $\mathit{w}$ in the orthogonal system $(\mathit{n},\mathit{m},\mathit{p})$ of the discontinuity $\mathcal{S}$ are given by:

$$\begin{array}{cc}\hfill & w_n:=\mathit{w}·\mathit{n}=\overline{\omega }{\Lambda }_{nn},w_m:=\mathit{w}·\mathit{m}=2\overline{\omega }{\Lambda }_{nm}\hfill \end{array}$$

(16a)

$$\begin{array}{c}{w}_p:=\mathit{w}·\mathit{p}=2\overline{\omega }{\Lambda }_{np}\hfill \end{array}$$

(16b)

Substitution of the above results into the relation (14) yields:

$$\begin{array}{cc}\hfill & {\Lambda }_{mm}\left({\theta }^{\mathrm{cr}}\right)=0\hfill \end{array}$$

(17a)

$$\begin{array}{c}{\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)=0,{\Lambda }_{mp}\left({\theta }^{\mathrm{cr}}\right)=0\hfill \end{array}$$

(17b)

for the discontinuity angles ${\theta }^{\mathrm{cr}}$ upon which the kinematic conditions (14) are satisfied. Calling for the relations of the local components $({\Lambda }_{mm},{\Lambda }_{pp},{\Lambda }_{mp})$ between the principal values ${\Lambda }_i(i=1,2,3)$ and a set of characteristic angles $\theta$, the kinematic constraints (17) yield a system of nonlinear equations, such that the discontinuity angles ${\theta }^{\mathrm{cr}}$ can be determined.

In particular, let us consider the 2D plane stress and plane strain conditions (${\sigma }_{np}={\sigma }_{mp}=0$). In such cases, the discontinuity orientation can be characterized by the inclination angle (anti-clockwise) ${\theta }^{\mathrm{cr}}\in [-\pi /2,\pi /2]$ between the normal vector $\mathit{n}$ of the discontinuity and the principal vector ${\mathit{v}}_1$ of the tensor $\Lambda$; see Figure 4. Furthermore, the in-plane components $({\Lambda }_{nn},{\Lambda }_{mm},{\Lambda }_{nm})$ are given by:

$$\begin{array}{c}{\Lambda }_{nn}={\Lambda }_1{cos}^2\theta +{\Lambda }_2{sin}^2\theta =\frac{{\Lambda }_1+{\Lambda }_2}{2}+\frac{{\Lambda }_1-{\Lambda }_2}{2}cos\left(2\theta \right)\hfill \end{array}$$

(18a)

$$\begin{array}{c}{\Lambda }_{mm}={\Lambda }_1{sin}^2\theta +{\Lambda }_2{cos}^2\theta =\frac{{\Lambda }_1+{\Lambda }_2}{2}-\frac{{\Lambda }_1-{\Lambda }_2}{2}cos\left(2\theta \right)\hfill \end{array}$$

(18b)

$$\begin{array}{c}{\Lambda }_{nm}=\left({\Lambda }_1-{\Lambda }_2\right)cos\theta sin\theta =\frac{{\Lambda }_1-{\Lambda }_2}{2}sin\left(2\theta \right)\hfill \end{array}$$

(18c)

for the in-plane principal values ${\Lambda }_1$ and ${\Lambda }_2$. As the out-of-plane constraint ${\Lambda }_{mp}\left({\theta }^{\mathrm{cr}}\right)=0$ is automatically fulfilled, it follows from the condition ${\Lambda }_{mm}\left({\theta }^{\mathrm{cr}}\right)=0$ that:

$$\begin{array}{c}\hfill {sin}^2{\theta }^{\mathrm{cr}}=-\frac{{\Lambda }_2}{{\Lambda }_1-{\Lambda }_2},{cos}^2{\theta }^{\mathrm{cr}}=\frac{{\Lambda }_1}{{\Lambda }_1-{\Lambda }_2}\end{array}$$

(19)

where the in-plane principle values ${\Lambda }_1$ and ${\Lambda }_2$ of the tensor $\Lambda$, satisfying ${\Lambda }_1\ge 0$ and ${\Lambda }_2\le 0$, are constrained by the conditions:

$$\begin{array}{cc}\hfill & {\sigma }_3\left({\theta }^{\mathrm{cr}}\right)={\sigma }_{pp}\left({\theta }^{\mathrm{cr}}\right)=0\mathrm{Plane}\mathrm{stress}\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill & {\Lambda }_3\left({\theta }^{\mathrm{cr}}\right)={\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)=0\mathrm{Plane}\mathrm{strain}\hfill \end{array}$$

(20b)

Note that in the plane stress state (${\sigma }_{pp}=0$), the constraint ${\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)=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 $\mathcal{B}$ with finite bandwidth $b\nrightarrow 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 ${\delta }_{{}_{\mathcal{S}}}$ is regularized by ${\delta }_{{}_{\mathcal{S}}}\left(\mathit{x}\right)\approx \Xi \left(\mathit{x}\right)/b$, where the collocation function $\Xi \left(\mathit{x}\right)$ is defined as $\Xi \left(\mathit{x}\right)=1$ for $\mathit{x}\in \mathcal{B}$ and $\Xi \left(\mathit{x}\right)=0$ otherwise; see [15,36].

Remark 2.

To gain further insight into the above strain localization analysis, let us consider a regularized discontinuity $\mathcal{B}$ in the solid Ω with the following internal virtual work:

$$\begin{array}{c}\hfill \delta W={\int }_{\Omega }\sigma :\delta {ϵ}^{\mathrm{e}}\mathrm{d}V+{\int }_{\mathcal{B}}\sigma :\delta {ϵ}^{\mathrm{d}}\mathrm{d}V\end{array}$$

(21)

For general 3D and plane strain cases, provided that the discontinuity angles ${\theta }^{\mathrm{cr}}$ satisfying the constraints (17) exist, the inelastic internal work localized in the discontinuity band $\mathcal{B}$ is expressed as:

$$\begin{array}{c}\hfill {\int }_{\mathcal{B}}\sigma :\delta {ϵ}^{d}d\mathcal{B}={\int }_{\mathcal{B}}\left({\sigma }_{nn}\delta {ϵ}_{nn}^d+2{\sigma }_{nm}\delta {ϵ}_{mm}^d+2{\sigma }_{np}\delta {ϵ}_{nm}^d\right)d\mathcal{B}={\int }_{\mathcal{S}}\left(t_n\delta w_n+t_m\delta w_m+t_m\delta w_p\right)d\mathcal{S}\end{array}$$

(22)

where the results (16) with the regularization relation $\omega =\overline{\omega }/b$ have been recalled. This is exactly the result for the elastic solid containing a frictional-cohesive crack. In the plane stress state (${\sigma }_{pp}=0$), only the condition (17a) is necessary for guaranteeing the identity:

$$\begin{array}{c}\hfill {\int }_{\mathcal{B}}\sigma :\delta {ϵ}^{d}d\mathcal{B}={\int }_{\mathcal{B}}\left({\sigma }_{nn}\delta {ϵ}_{nn}^d+2{\sigma }_{nm}\delta {ϵ}_{nm}^d\right)d\mathcal{B}={\int }_{\mathcal{S}}\left(t_n\delta w_n+t_m\delta w_m\right)d\mathcal{S}\end{array}$$

(23)

with the constraint ${\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)=0$ being irrelevant. That is, regardless of the constraint ${\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)=0$, strain localization also occurs in the plane stress condition, with the resulting discontinuity characterized by the cohesive relations (16a).

3. 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.

3.1. $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.

3.1.1. Constitutive Relations

In the $J_2$ damage model, the stress and strain tensors are decomposed as:

$$\begin{array}{c}\hfill \sigma =p\mathit{I}+\mathit{s},ϵ=\frac{1}{3}{ϵ}_v\mathit{I}+\mathit{e}\end{array}$$

(24)

where $p=\frac{1}{3}\mathrm{tr}\sigma$ and $\mathit{s}=\sigma -p\mathit{I}$ are the volumetric and deviatoric parts of the stress tensor $\sigma$, respectively; ${ϵ}_v=\mathrm{tr}ϵ$ and $\mathit{e}:=ϵ-\frac{1}{3}{ϵ}_v\mathit{I}$ represent the trace and the deviatoric part of the strain tensor $ϵ$, respectively.

While the volumetric behavior remains elastic, a scalar variable $d\in [0,1]$ is introduced to characterize the deviatoric behavior. It then follows that:

$$\begin{array}{c}\hfill p=K_0{ϵ}_v,\mathit{s}=\left(1-d\right)\overline{\mathit{s}}=\left(1-d\right)2G_0\mathit{e}\end{array}$$

(25)

where $\overline{\mathit{s}}=2G_0\mathit{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:

$$\begin{array}{c}\hfill ϵ=\mathbb{C}:\sigma ,\mathbb{C}=\frac{1}{3K_0}{\mathbb{I}}_{\mathrm{vol}}+\frac{1}{2\left(1-d\right)G_0}{\mathbb{I}}_{\mathrm{dev}}\end{array}$$

(26)

for the volumetric and deviatoric projection operators:

$$\begin{array}{c}\hfill {\mathbb{I}}_{\mathrm{vol}}=\frac{1}{3}\mathit{I}\otimes \mathit{I},{\mathbb{I}}_{\mathrm{dev}}=\mathbb{I}-{\mathbb{I}}_{\mathrm{vol}}=\mathbb{I}-\frac{1}{3}\mathit{I}\otimes \mathit{I}\end{array}$$

(27)

The damage compliance tensor ${\mathbb{C}}^{\mathrm{d}}$ and the resulting inelastic strain ${ϵ}^{\mathrm{d}}$ are then given by:

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{d}}=\mathbb{C}-{\mathbb{C}}_0={\displaystyle \frac{\omega }{2G_0}}{\mathbb{I}}_{\mathrm{dev}},{ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma ={\displaystyle \frac{\omega }{2G_0}}\mathit{s}=\omega \Lambda \end{array}$$

(28)

for the variable $\omega :=d/(1-d)$ and the characteristic tensor $\Lambda =\mathit{s}/\left(2G_0\right)$.

The above $J_2$ damage model corresponds to the following Helmholtz free energy $\psi$:

$$\begin{array}{c}\hfill \psi ={\displaystyle \frac{1}{2}}{K}_0{ϵ}_v^2+{\displaystyle \frac{1}{2}}\left(1-d\right)2G_0\mathit{e}:\mathit{e}\end{array}$$

(29)

together with the energy dissipation inequality:

$$\begin{array}{c}\hfill \mathcal{D}=Y·\dot{d}\ge 0Y=-{\displaystyle \frac{\partial \psi }{\partial d}}=G_0\mathit{e}:\mathit{e}=\frac{1}{4G_0}\overline{\mathit{s}}:\overline{\mathit{s}}\end{array}$$

(30)

for the energy release rate Y, conjugate to the damage variable d.

Additionally, a strain-based (or effective stress-based) damage criterion is introduced:

$$\begin{array}{c}\hfill F_{ϵ}(ϵ,r)={ϵ}_{\mathrm{eq}}\left(ϵ\right)-r\le 0\end{array}$$

(31)

in terms of the equivalent strain ${ϵ}_{\mathrm{eq}}\left(ϵ\right)$ and the threshold r:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{eq}}\left(ϵ\right)=\sqrt{\frac{1}{2}\overline{\mathit{s}}:\overline{\mathit{s}}}=\sqrt{{\overline{J}}_2},r=\underset{t\in [0,T]}{max}\left(r_0,{ϵ}_{\mathrm{eq}}^t\right)\end{array}$$

(32)

for the second invariant ${\overline{J}}_2=\frac{1}{2}\overline{\mathit{s}}:\overline{\mathit{s}}$ of the effective deviatoric stress tensor $\overline{\mathit{s}}=2G_0\mathit{e}$. Upon damage loading, it follows from the consistency condition ${\dot{F}}_{ϵ}(ϵ,r)=0$ that the threshold r, with the initial value $r_0$, represents the maximum value of ${ϵ}_{\mathrm{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:

$$\begin{array}{c}\hfill F_{\sigma }(\sigma ,q)={\sigma }_{\mathrm{eq}}\left(\sigma \right)-q\le 0\end{array}$$

(33)

in terms of the following equivalent stress ${\sigma }_{\mathrm{eq}}$ and the corresponding softening variable $q\left(r\right)$:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{eq}}=\sqrt{\frac{1}{2}\mathit{s}:\mathit{s}}=\sqrt{J_2}=\left(1-d\right){ϵ}_{\mathrm{eq}},q=\left(1-d\right)r\end{array}$$

(34)

where $J_2=\frac{1}{2}\mathit{s}:\mathit{s}$ represents the second invariant of the deviatoric stress tensor $\mathit{s}$.

3.1.2. Orientation of the Discontinuity

For the $J_2$ damage model with the characteristic tensor $\Lambda =\mathit{s}/\left(2G_0\right)$, the conditions (17) upon strain localization become:

$$\begin{array}{c}\hfill s_{mm}\left({\theta }^{\mathrm{cr}}\right)=s_{pp}\left({\theta }^{\mathrm{cr}}\right)=s_{mp}\left({\theta }^{\mathrm{cr}}\right)=0\end{array}$$

(35)

for the deviatoric stress components $(s_{mm},s_{pp},s_{mp})$ on the discontinuity surface. Accordingly, the orientation angles ${\theta }^{\mathrm{cr}}$ can be solved in terms of the principal values $s_i(i=1,2,3)$ of the deviatoric stress tensor $\mathit{s}$.

In particular, for the 2D cases of plane stress and plane strain (${\sigma }_{np}={\sigma }_{mp}=0$), the condition $s_{mm}\left({\theta }^{\mathrm{cr}}\right)=0$ gives:

$$\begin{array}{c}\hfill {sin}^2{\theta }^{\mathrm{cr}}=-\frac{s_2}{s_1-s_2},{cos}^2{\theta }^{\mathrm{cr}}=\frac{s_1}{s_1-s_2}\end{array}$$

(36)

Furthermore, it follows from the constraint $s_3=s_{pp}=0$ that:

$$\begin{array}{c}\hfill s_1=-s_2⟹{\theta }^{\mathrm{cr}}=\pm {45}^{°}\end{array}$$

(37)

That is, in 2D conditions strain localization of a $J_2$ damaging material with stress continuity can always occur along the discontinuity angle ${\theta }^{\mathrm{cr}}=\pm {45}^{°}$. Note that the above analytical result coincides with the numerical simulations [46].

3.1.3. Localized Damage Model

Provided that the solution to the constraints (35) exists, it follows from the vanishing trace $\mathrm{tr}\mathit{s}=0$ that $s_{nn}=-\left(s_{mm}+s_{pp}\right)=0$ holds on the discontinuity surface with the normal $\mathit{n}\left({\theta }^{\mathrm{cr}}\right)$. Accordingly, the relation (14) becomes:

$$\begin{array}{c}\hfill \overline{\omega }{\displaystyle \frac{\mathit{s}}{2G_0}}={\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}=w_m{\left(\mathit{n}\otimes \mathit{m}\right)}^{\mathrm{sym}}+w_p{\left(\mathit{n}\otimes \mathit{p}\right)}^{\mathrm{sym}}\end{array}$$

(38a)

or, equivalently,

$$\begin{array}{c}\hfill \mathit{s}={\displaystyle \frac{2G_0}{\overline{\omega }}}\left[w_m{\left(\mathit{n}\otimes \mathit{m}\right)}^{\mathrm{sym}}+w_p{\left(\mathit{n}\otimes \mathit{p}\right)}^{\mathrm{sym}}\right]\end{array}$$

(38b)

The above frictional-cohesive relations can also be given straightforwardly from Equation (16):

$$\begin{array}{c}\hfill w_n=\overline{\omega }{\displaystyle \frac{s_{nn}}{2G_0}}=0,w_m=\overline{\omega }{\displaystyle \frac{s_{nm}}{G_0}}=\overline{\omega }{\displaystyle \frac{t_m}{G_0}},w_p=\overline{\omega }{\displaystyle \frac{s_{np}}{G_0}}=\overline{\omega }{\displaystyle \frac{t_p}{G_0}}\end{array}$$

(39)

where the identities $s_{nm}={\sigma }_{nm}=t_m$ and $s_{np}={\sigma }_{np}=t_p$ have been considered.

With the frictional-cohesive relations (39) upon strain localization, the stress-based damage criterion (33) becomes:

$$\begin{array}{c}\hfill f(\mathit{t},q)=F_{\sigma }(\sigma ,q)={\displaystyle \frac{G_0}{\overline{\omega }}}{w}_{\mathrm{eq}}-q=t_{\mathrm{eq}}-q\le 0\end{array}$$

(40)

or the alternative displacement jump-based one:

$$\begin{array}{c}\hfill g(\mathit{w},\kappa )=w_{\mathrm{eq}}-\kappa \le 0,\kappa ={\displaystyle \frac{\overline{\omega }}{G_0}}q={\displaystyle \frac{1}{G_0}}{\displaystyle \frac{q-q_0}{\overline{\mathcal{H}}}}\end{array}$$

(41)

where $w_{\mathrm{eq}}$ and $t_{\mathrm{eq}}$ denote the equivalent displacement jump and traction, defined as:

$$\begin{array}{c}\hfill w_{\mathrm{eq}}=\sqrt{w_m^2+w_p^2},t_{\mathrm{eq}}=\sqrt{t_m^2+t_p^2}\end{array}$$

(42)

Similarly to the strain-based threshold r, the displacement jump-like internal variable $\kappa$ represents the maximum value of the equivalent displacement jump $w_{\mathrm{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.

3.2. 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.

3.2.1. Constitutive Relations

The Simó and Ju [41] damage model is characterized by the following well-defined Helmholtz free energy:

$$\begin{array}{c}\hfill \psi =\frac{1}{2}\left(1-d\right)ϵ:{\mathbb{E}}_0:ϵ\end{array}$$

(43)

where ${\mathbb{E}}_0={\lambda }_0\mathit{I}\otimes \mathit{I}+2{\mu }_0\mathbb{I}$ represents the fourth-order material elasticity tensor, with ${\lambda }_0$ and ${\mu }_0$ being the Lamé constants. Standard arguments then yield:

$$\begin{array}{c}\hfill \sigma =\frac{\partial \psi }{\partial ϵ}=\left(1-d\right)\overline{\sigma },\overline{\sigma }={\mathbb{E}}_0:ϵ\end{array}$$

(44)

or, equivalently,

$$\begin{array}{c}\hfill ϵ=\mathbb{C}:\sigma ,\mathbb{C}=\frac{1}{1-d}{\mathbb{C}}_0\end{array}$$

(45)

Similarly, the damage compliance tensor ${\mathbb{C}}^{\mathrm{d}}$ and the resulting damage strain ${ϵ}^{\mathrm{d}}$ are characterized in the form (4):

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{d}}=\mathbb{C}-{\mathbb{C}}_0=\omega {\mathbb{C}}_0,{ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma =\omega {\mathbb{C}}_0:\sigma =\omega \Lambda \end{array}$$

(46)

for the variable $\omega :=d/(1-d)$ and the characteristic tensor $\Lambda ={\mathbb{C}}_0:\sigma$.

The above constitutive relations are constrained by the damage energy dissipation inequality:

$$\begin{array}{c}\hfill \mathcal{D}=Y·\dot{d}\ge 0,Y=-\frac{\partial \psi }{\partial d}=\frac{1}{2}ϵ:{\mathbb{E}}_0:ϵ={\displaystyle \frac{1}{2}}\overline{\sigma }:{\mathbb{C}}_0:\overline{\sigma }\end{array}$$

(47)

It is then possible to introduce the Beltrami strain-based damage criterion:

$$\begin{array}{c}\hfill F_{ϵ}(ϵ,r)={ϵ}_{\mathrm{eq}}\left(ϵ\right)-r\le 0\end{array}$$

(48)

where the equivalent strain ${ϵ}_{\mathrm{eq}}$ and the threshold r are defined as:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{eq}}\left(ϵ\right)=\sqrt{M_0\left(\overline{\sigma }:{\mathbb{C}}_0:\overline{\sigma }\right)},r=\underset{t\in [0,T]}{max}\left(r_0,{ϵ}_{\mathrm{eq}}^t\right)\end{array}$$

(49)

for the longitudinal or constrained modulus $M_0={\lambda }_0+2{\mu }_0$. Alternatively, for the damage evolution law (5), a stress-based failure criterion can be expressed as:

$$\begin{array}{c}\hfill F_{\sigma }(\sigma ,q)={\sigma }_{\mathrm{eq}}\left(\sigma \right)-q\le 0\end{array}$$

(50)

for the equivalent stress ${\sigma }_{\mathrm{eq}}$ and the corresponding stress-like internal variable q:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{eq}}\left(\sigma \right)=\sqrt{M_0\left(\sigma :{\mathbb{C}}_0:\sigma \right)}=\left(1-d\right){ϵ}_{\mathrm{eq}},q=\left(1-d\right)r\end{array}$$

(51)

Once the softening function $q\left(r\right)$ is postulated, the damage variable d can be determined as shown later.

3.2.2. Orientation of the Discontinuity

For the Simó and Ju [41] damage model, the characteristic tensor $\Lambda$ is expressed as:

$$\begin{array}{c}\hfill \Lambda ={\mathbb{C}}_0:\sigma =\frac{1}{E_0}\left[\left(1+{\nu }_0\right)\sigma -{\nu }_0\mathrm{tr}\sigma \mathit{I}\right]\end{array}$$

(52)

so that the constraints (17) give:

$$\begin{array}{c}\hfill {\sigma }_{mm}-{\nu }_0\left({\sigma }_{nn}+{\sigma }_{pp}\right)=0,{\sigma }_{pp}-{\nu }_0\left({\sigma }_{nn}+{\sigma }_{mm}\right)=0,{\sigma }_{mp}=0\end{array}$$

(53)

With the stress components expressed in terms of the principal stresses ${\sigma }_i(i=1,2,3)$ and a set of characteristic angles $\theta$, the discontinuity angles ${\theta }^{\mathrm{cr}}$ can be determined.

In particular, let us focus on the 2D plane stress and plane strain cases in which the identity ${\sigma }_{np}={\sigma }_{mp}=0$ holds. For the plane stress state, the conditions (19) and (20a) give:

$$\begin{array}{c}\hfill {\sigma }_{mm}={\nu }_0{\sigma }_{nn}⟹cos\left(2{\theta }^{\mathrm{cr}}\right)={\displaystyle \frac{1-{\nu }_0}{1+{\nu }_0}}·\frac{{\sigma }_1+{\sigma }_2}{{\sigma }_1-{\sigma }_2}\end{array}$$

(54a)

For the plane strain case (i.e., ${\sigma }_3={\sigma }_{pp}\ne 0$), it follows from the conditions (19) and (20a) that:

$$\begin{array}{c}\hfill {\sigma }_{mm}={\sigma }_{pp}=\frac{{\nu }_0}{1-{\nu }_0}{\sigma }_{nn}⟹cos\left(2{\theta }^{\mathrm{cr}}\right)=\left(1-2{\nu }_0\right)\frac{{\sigma }_1+{\sigma }_2}{{\sigma }_1-{\sigma }_2}\end{array}$$

(54b)

Note that the convention ${\sigma }_1>{\sigma }_2$ is assumed as usual. As can be seen, both results depend on Poisson’s ratio ${\nu }_0$ and coincide for a vanishing one.

3.2.3. Localized Damage Model

Provided that the solution to the conditions (53) exists, it follows from the relation (14) that:

$$\begin{array}{c}\hfill \sigma ={\displaystyle \frac{1}{\overline{\omega }}}{\mathbb{E}}_0:{\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}\end{array}$$

(55)

Accordingly, the frictional-cohesive relations are given by:

$$\begin{array}{c}\hfill \mathit{t}=\mathit{n}·\sigma ={\displaystyle \frac{1}{\overline{\omega }}}{\mathit{Q}}_0·\mathit{w},\mathit{w}=\overline{\omega }{\mathit{Q}}_0^{-1}·\mathit{t}\end{array}$$

(56a)

or, in the component form,

$$\begin{array}{c}\hfill w_n=\overline{\omega }{\Lambda }_{nn}={\displaystyle \frac{\overline{\omega }}{M_0}}{\sigma }_{nn},w_m=\overline{\omega }{\Lambda }_{nm}={\displaystyle \frac{\overline{\omega }}{G_0}}{\sigma }_{nm},w_p=\overline{\omega }{\Lambda }_{np}={\displaystyle \frac{\overline{\omega }}{G_0}}{\sigma }_{np}\end{array}$$

(56b)

where the second-order elastic acoustic tensor ${\mathit{Q}}_0$ and its inverse ${\mathit{Q}}_0^{-1}$ are expressed as:

$$\begin{array}{cc}\hfill {\mathit{Q}}_0:& =\mathit{n}·{\mathbb{E}}_0·\mathit{n}=M_0\mathit{n}\otimes \mathit{n}+G_0\left(\mathit{m}\otimes \mathit{m}+\mathit{p}\otimes \mathit{p}\right)\hfill \end{array}$$

(57a)

$$\begin{array}{cc}\hfill {\mathit{Q}}_0^{-1}& ={\displaystyle \frac{1}{M_0}}\mathit{n}\otimes \mathit{n}+{\displaystyle \frac{1}{G_0}}\left(\mathit{m}\otimes \mathit{m}+\mathit{p}\otimes \mathit{p}\right)\hfill \end{array}$$

(57b)

in the local $(n,m,p)$ coordinate system.

With the relation (55), upon strain localization, the stress-based damage criterion (50) becomes:

$$\begin{array}{c}\hfill f(\mathit{t},q)=F_{\sigma }(\sigma ,q)={\displaystyle \frac{M_0}{\overline{\omega }}}{w}_{\mathrm{eq}}-q=t_{\mathrm{eq}}-q\le 0\end{array}$$

(58)

or, the alternative form,

$$\begin{array}{c}\hfill g(\mathit{w},\kappa )=w_{\mathrm{eq}}-\kappa \le 0,\kappa ={\displaystyle \frac{\overline{\omega }}{M_0}}q={\displaystyle \frac{1}{M_0}}{\displaystyle \frac{q-q_0}{\overline{\mathcal{H}}}}\end{array}$$

(59)

where the equivalent displacement jump $w_{\mathrm{eq}}$ and traction $t_{\mathrm{eq}}$ are expressed as:

$$\begin{array}{c}{w}_{\mathrm{eq}}=\sqrt{\frac{1}{M_0}\left(\mathit{w}·{\mathit{Q}}_0·\mathit{w}\right)}=\sqrt{w_n^2+{\beta }_0^2\left(w_m^2+w_p^2\right)}\hfill \end{array}$$

(60a)

$$\begin{array}{cc}\hfill & t_{\mathrm{eq}}=\sqrt{M_0\left(\mathit{t}·{\mathit{Q}}_0^{-1}·\mathit{t}\right)}=\sqrt{t_n^2+{\beta }_0^{-2}\left(t_m^2+t_p^2\right)}\hfill \end{array}$$

(60b)

for the coefficient ${\beta }_0=\sqrt{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.

3.3. 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.

3.3.1. Constitutive Relations

In the modified damage model, the constitutive relations read:

$$\begin{array}{c}\hfill ϵ=\mathbb{C}:\sigma ,\mathbb{C}={\mathbb{C}}_0+{\mathbb{C}}^{\mathrm{d}}\end{array}$$

(61)

which gives the following damage strain tensor ${ϵ}^{\mathrm{d}}$:

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{d}}={\mathbb{C}}^{\mathrm{d}}:\sigma =\omega \Lambda ,\Lambda ={\overline{\mathbb{C}}}_0:\sigma \end{array}$$

(62)

where the compliance tensor ${\mathbb{C}}^{\mathrm{d}}$ is expressed as the form (4):

$$\begin{array}{c}\hfill {\mathbb{C}}^{\mathrm{d}}=\mathbb{C}-{\mathbb{C}}_0=\omega {\overline{\mathbb{C}}}_0,{\overline{\mathbb{C}}}_0=\mathrm{DIAG}\left[{\displaystyle \frac{1}{E_0}},{\displaystyle \frac{1}{E_0}},{\displaystyle \frac{1}{E_0}},{\displaystyle \frac{1}{G_0}},{\displaystyle \frac{1}{G_0}},{\displaystyle \frac{1}{G_0}}\right]\end{array}$$

(63)

for the variable $\omega :=d/(1-d)$. Here, the diagonal matrix ${\overline{\mathbb{C}}}_0$ is the Voigt representation of a fourth-order reference compliance tensor, which is obtained from the elastic one ${\mathbb{C}}_0$ by ignoring the lateral deformations caused by Poisson’s ratio. The explicit expressions of the resulting secant compliance $\mathbb{C}$ and stiffness $\mathbb{E}$ are given in Appendix A.

As these constitutive tensors are of major symmetry, the corresponding energy function $\psi$ is defined as:

$$\begin{array}{c}\hfill \psi (ϵ,d)={\displaystyle \frac{1}{2}}ϵ:\mathbb{E}:ϵ={\displaystyle \frac{1}{2}}\sigma :\mathbb{C}:\sigma ={\displaystyle \frac{1}{2}}\sigma :{\mathbb{C}}_0:\sigma +{\displaystyle \frac{d}{2\left(1-d\right)}}\sigma :{\overline{\mathbb{C}}}_0:\sigma \end{array}$$

(64)

Standard arguments give the constitutive relations (61) and the following energy dissipation inequality:

$$\begin{array}{c}\hfill \mathcal{D}=Y·\dot{d}\ge 0,Y=-{\displaystyle \frac{\partial \psi }{\partial d}}={\displaystyle \frac{\sigma :{\overline{\mathbb{C}}}_0:\sigma }{2{\left(1-d\right)}^2}}={\displaystyle \frac{1}{2}}\overline{\sigma }:{\overline{\mathbb{C}}}_0:\overline{\sigma }\end{array}$$

(65)

for the effective stress tensor $\overline{\sigma }=\sigma /(1-d)$. Accordingly, the following strain-based damage criterion can be considered:

$$\begin{array}{c}\hfill F_{ϵ}(ϵ,r)={ϵ}_{\mathrm{eq}}\left(ϵ\right)-r\le 0\end{array}$$

(66)

where the equivalent strain ${ϵ}_{\mathrm{eq}}$ and the threshold r are given by:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{eq}}\left(ϵ\right)=\sqrt{E_0\left(\overline{\sigma }:{\overline{\mathbb{C}}}_0:\overline{\sigma }\right)},r=\underset{t\in [0,T]}{max}\left(r_0,{ϵ}_{\mathrm{eq}}^t\right)\end{array}$$

(67)

Alternatively, the stress-based damage criterion is expressed as:

$$\begin{array}{c}\hfill F_{\sigma }(\sigma ,q)={\sigma }_{\mathrm{eq}}\left(\sigma \right)-q\le 0\end{array}$$

(68)

for the equivalent stress ${\sigma }_{\mathrm{eq}}$ and internal variable q:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{eq}}=\sqrt{E_0\left(\sigma :{\overline{\mathbb{C}}}_0:\sigma \right)},q=\left(1-d\right)r\end{array}$$

(69)

Note that the above equivalent strain and stress are subtly different from their counterparts in the Simó and Ju [41] damage model.

3.3.2. Orientation of the Discontinuity

For the damage strain (62), the constraints (17) become:

$$\begin{array}{c}\hfill {\Lambda }_{mm}\left({\theta }^{\mathrm{cr}}\right)={\displaystyle \frac{{\sigma }_{mm}}{E_0}}=0,{\Lambda }_{pp}\left({\theta }^{\mathrm{cr}}\right)={\displaystyle \frac{{\sigma }_{pp}}{G_0}}=0,{\Lambda }_{mp}\left({\theta }^{\mathrm{cr}}\right)={\displaystyle \frac{{\sigma }_{mp}}{G_0}}=0\end{array}$$

(70)

With the stress components expressed in terms of the principal stresses ${\sigma }_i(i=1,2,3)$ and a set of characteristic angles $\theta$, the discontinuity angles ${\theta }^{\mathrm{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 ${\Lambda }_{mm}\left({\theta }^{\mathrm{cr}}\right)=0$ that:

$$\begin{array}{c}\hfill {\sigma }_{mm}\left({\theta }^{\mathrm{cr}}\right)=0⟹cos\left(2{\theta }^{\mathrm{cr}}\right)=\frac{{\sigma }_1+{\sigma }_2}{{\sigma }_1-{\sigma }_2}\end{array}$$

(71)

The above result applies for the biaxial tension-compression quadrant, i.e., ${\sigma }_1>0>{\sigma }_2$. For the case ${\sigma }_1>{\sigma }_2>0$, the discontinuity angle is given from $cos\left(2{\theta }^{\mathrm{cr}}\right)=1$, i.e., ${\theta }^{\mathrm{cr}}=0$. That is, the discontinuity orientation $\mathit{n}$ coincides with the major principle vector ${\mathit{v}}_1$ of the stress tensor $\sigma$. 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].

3.3.3. Localized Damage Model

Provided the solution to the conditions (53) exists, it follows from the relation (14) that:

$$\begin{array}{c}\hfill \sigma ={\displaystyle \frac{1}{\overline{\omega }}}{\overline{\mathbb{E}}}_0:{\left(\mathit{w}\otimes \mathit{n}\right)}^{\mathrm{sym}}\end{array}$$

(72)

for the displacement jumps $\mathit{w}:={\left\{w_n,w_m,w_p\right\}}^{\mathrm{T}}$. Accordingly, the frictional-cohesive tractions $\mathit{t}:={\left\{t_n,t_m,t_p\right\}}^{\mathrm{T}}$ are given by:

$$\begin{array}{c}\hfill \mathit{t}=\mathit{n}·\sigma ={\overline{\mathit{E}}}_{{}_{{}_{\mathcal{S}}}}·\mathit{w},{\overline{\mathit{E}}}_{{}_{{}_{\mathcal{S}}}}={\displaystyle \frac{1}{\overline{\omega }}}{\overline{\mathit{E}}}_0\end{array}$$

(73a)

or, inversely,

$$\begin{array}{c}\hfill \mathit{w}={\overline{\mathit{C}}}_{{}_{{}_{\mathcal{S}}}}·\mathit{t},{\overline{\mathit{C}}}_{{}_{{}_{\mathcal{S}}}}={\mathit{E}}_{{}_{{}_{\mathcal{S}}}}^{-1}=\overline{\omega }{\overline{\mathit{C}}}_0\end{array}$$

(73b)

for the reference stiffness ${\overline{\mathit{E}}}_0$ and compliance ${\overline{\mathit{C}}}_0$ of the discontinuity:

$$\begin{array}{c}{\overline{\mathit{E}}}_0=\mathit{n}·{\overline{\mathbb{E}}}_0·\mathit{n}=E_0\mathit{n}\otimes \mathit{n}+G_0\left(\mathit{m}\otimes \mathit{m}+\mathit{p}\otimes \mathit{p}\right)\hfill \end{array}$$

(74a)

$$\begin{array}{c}{\overline{\mathit{C}}}_0={\overline{\mathit{E}}}_0^{-1}={\displaystyle \frac{1}{E_0}}\mathit{n}\otimes \mathit{n}+{\displaystyle \frac{1}{G_0}}\left(\mathit{m}\otimes \mathit{m}+\mathit{p}\otimes \mathit{p}\right)\hfill \end{array}$$

(74b)

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:

$$\begin{array}{c}\hfill f(\mathit{t},q)=F_{\sigma }(\sigma ,q)={\displaystyle \frac{E_0}{\overline{\omega }}}{w}_{\mathrm{eq}}-q=t_{\mathrm{eq}}-q\le 0\end{array}$$

(75)

or, equivalently,

$$\begin{array}{c}\hfill g(\mathit{w},\kappa )=w_{\mathrm{eq}}-\kappa \le 0,\kappa ={\displaystyle \frac{\overline{\omega }}{E_0}}q={\displaystyle \frac{1}{E_0}}{\displaystyle \frac{q-q_0}{\overline{\mathcal{H}}}}\end{array}$$

(76)

where the equivalent displacement jump $w_{\mathrm{eq}}$ and traction $t_{\mathrm{eq}}$ are expressed as:

$$\begin{array}{c}\hfill w_{\mathrm{eq}}=\sqrt{w_n^2+{\beta }_0^2\left(w_m^2+w_p^2\right)},t_{\mathrm{eq}}=\sqrt{t_n^2+{\beta }_0^{-2}\left(t_m^2+t_p^2\right)}\end{array}$$

(77)

for the material property ${\beta }_0=\sqrt{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 ${ϵ}^{\mathrm{e}}\left(\mathit{x}\right)$ 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.

4. 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 $ϵ-\mathit{u}$ finite element method, recently developed by Cervera et al. [39] for elasticity and extended to isotropic damage materials [40,49], is considered.

4.1. 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 $\Omega$ is described by the compatibility of deformations and the equilibrium of body forces:

$$\begin{array}{c}-ϵ+{\nabla }^{\mathrm{sym}}\mathit{u}=\mathit{0}\hfill \end{array}$$

(78a)

$$\begin{array}{c}\nabla ·\sigma +{\mathit{b}}^{*}=\mathit{0}\hfill \end{array}$$

(78b)

where $\mathit{u}$ is the displacement vector, $ϵ$ is the strain tensor, $\sigma$ represents the stress tensor, ${\nabla }^{\mathrm{sym}}$ and $\nabla ·(·)$ are the adjoint symmetric gradient and the divergence operators, respectively; ${\mathit{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 $\partial \Omega$.

After pre-multiplying the secant stiffness tensor $\mathbb{E}$, the strong form reads:

$$\begin{array}{c}-\mathbb{E}:ϵ+\mathbb{E}:{\nabla }^{\mathrm{sym}}\mathit{u}=\mathit{0}\hfill \end{array}$$

(79a)

$$\begin{array}{c}\nabla ·\left(\mathbb{E}:ϵ\right)+{\mathit{b}}^{*}=\mathit{0}\hfill \end{array}$$

(79b)

for the unknowns fields of total strains $ϵ$ and displacements $\mathit{u}$. Introducing the test function $\gamma \in \mathcal{G}\subset {\left[L^2\left(\Omega \right)\right]}^{\dim }$ for strains and $\mathit{v}\in \mathcal{V}\subset {\left[H^1\left(\Omega \right)\right]}^{\dim }$ for displacements and applying Gauss’s divergence theorem to the strong form yields the following weak form of the mixed problem:

$$\begin{array}{c}-{\int }_{\Omega }\gamma :\mathbb{E}:ϵ\mathrm{d}V+{\int }_{\Omega }\gamma :\mathbb{E}:{\nabla }^{\mathrm{sym}}\mathit{u}\mathrm{d}V=\mathit{0}\forall \gamma \in \mathcal{G}\hfill \end{array}$$

(80a)

$$\begin{array}{c}{\int }_{\Omega }{\nabla }^{\mathrm{sym}}\mathit{v}:\left(\mathbb{E}:ϵ\right)\mathrm{d}V=F\left(\mathit{v}\right)\forall \mathit{v}\in \mathcal{V}\hfill \end{array}$$

(80b)

where $F\left(\mathit{v}\right)$ represents the work done by boundary tractions ${\mathit{t}}^{*}$ and body forces ${\mathit{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:

$$\begin{array}{c}\hfill ϵ\to {ϵ}_h=\sum _{i=1}^{n_{pts}}{\gamma }_h^{\left(i\right)}{ϵ}_h^{\left(i\right)}{\gamma }_h\in {\mathcal{G}}_h;\mathit{u}\to {\mathit{u}}_h=\sum _{i=1}^{n_{pts}}{\mathit{v}}_h^{\left(i\right)}{\mathit{u}}_h^{\left(i\right)}{\mathit{v}}_h\in {\mathcal{V}}_h\end{array}$$

(81)

where ${ϵ}_h$ and ${\mathit{u}}_h$ are the nodal degrees of freedom; ${\gamma }_h$ and ${\mathit{v}}_h$ are the discrete test functions for the strains and displacements pertaining to the spaces ${\mathcal{G}}_h$ and ${\mathcal{V}}_h$, respectively, the discrete counterparts of $\mathcal{G}$ and $\mathcal{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:

$$\begin{array}{c}-\left(1-{\tau }_{ϵ}\right){\int }_{\Omega }{\gamma }_h:\mathbb{E}:\left({ϵ}_h-{\nabla }^{\mathrm{sym}}{\mathit{u}}_h\right)\mathrm{d}V=0\forall {\gamma }_h\in {\mathcal{G}}_h\hfill \end{array}$$

(82a)

$$\begin{array}{c}{\int }_{\Omega }{\nabla }^{\mathrm{sym}}{\mathit{v}}_h:\mathbb{E}:{\widetilde{ϵ}}_h\mathrm{d}V=F\left({\mathit{v}}_h\right)\forall {\mathit{v}}_h\in {\mathcal{V}}_h\hfill \end{array}$$

(82b)

where the stabilized discrete strain field ${\widetilde{ϵ}}_h$ is a blending of the continuous and discontinuous strain fields $({ϵ}_h,{\nabla }^{\mathrm{sym}}{\mathit{u}}_h)$ weighted by the stabilization parameter ${\tau }_{ϵ}$:

$$\begin{array}{c}\hfill {\widetilde{ϵ}}_h=\left(1-{\tau }_{ϵ}\right){ϵ}_h+{\tau }_{ϵ}{\nabla }^{\mathrm{sym}}{\mathit{u}}_h\mathrm{with}{\tau }_{ϵ}=c_{ϵ}{\displaystyle \frac{h}{L_0}}\end{array}$$

(83)

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 ${\mathit{u}}_h$, the contribution of the stabilization terms (those multiplied by ${\tau }_{ϵ}$) vanishes.

4.2. 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 $\times 1$ m (width × height × thickness). An imperfection is introduced with a slanted perforation of diameter $D=1$ m such that symmetric solutions are excluded.

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:

$$\begin{array}{c}\hfill H=-{\displaystyle \frac{l_{\mathrm{ch}}}{L_{\mathrm{H}}-0.5l_{\mathrm{ch}}}}\end{array}$$

(84)

where $l_{\mathrm{ch}}$ and $L_{\mathrm{H}}:=E_0{G}_{\mathrm{F}}/q_0^2$ represent the localization and Griffith’s characteristic lengths, respectively. In the mixed formulation, $l_{\mathrm{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 ${\nu }_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].

4.2.1. Simó and Ju [41] Damage Model

Let us first discuss the Simó and Ju [41] damage model. For the considered uniaxial stress state (${\sigma }_1>0$ and ${\sigma }_2=0$), it follows from the closed-form results (54) that:

$$\begin{array}{c}\hfill cos\left(2{\theta }^{\mathrm{cr}}\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1-{\nu }_0}{1+{\nu }_0}}\hfill & \mathrm{Plane}\mathrm{stress}\hfill \\ 1-2{\nu }_0\hfill & \mathrm{Plane}\mathrm{strain}\hfill \end{array}\right.\end{array}$$

(85)

The resulting discontinuity angles ${\theta }^{\mathrm{cr}}$ for various Poisson’s ratio ${\nu }_0$ are summarized in

Table 1. Comparison between the analytically- and numerically-obtained discontinuity angles ${\theta }^{\mathrm{cr}}$ for the Simó and Ju [41] damage model.

Stress State	${\mathit{\nu }}_0=0.00$		${\mathit{\nu }}_0=0.15$		${\mathit{\nu }}_0=0.30$		${\mathit{\nu }}_0=0.45$
	Analytical	Numerical	Analytical	Numerical	Analytical	Numerical	Analytical	Numerical
Plane stress	$0.{00}^{°}$	$0.{00}^{°}$	$21.{17}^{°}$	$20.{75}^{°}$	$28.{71}^{°}$	$28.{69}^{°}$	$33.{85}^{°}$	$33.{13}^{°}$
Plane strain	$0.{00}^{°}$	$0.{00}^{°}$	$22.{79}^{°}$	$22.{31}^{°}$	$33.{21}^{°}$	$32.{27}^{°}$	$42.{13}^{°}$	$41.{48}^{°}$

.

Figure 6 and Figure 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 ${\theta }^{\mathrm{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 ${\nu }_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 ${\nu }_0=0.0,0.15,0.30$ and 0.45 are given in Figure 9.

4.2.2. Modified Simó and Ju Damage Model

For the modified Simó and Ju damage model, the discontinuity angle ${\theta }^{\mathrm{cr}}$ given by the analytical results (71) is independent of Poisson’s ratio ${\nu }_0$. In particular, for the uniaxial stress state, it follows that:

$$\begin{array}{c}\hfill cos{\theta }^{\mathrm{cr}}=1⟹{\theta }^{\mathrm{cr}}=0\end{array}$$

(86)

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., ${\theta }^{\mathrm{cr}}=0$. The novel strain localization analysis and the mixed stabilized finite element method are verified again.

5. 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.