<figure><img src="image_1.png" /></figure>

HAL
open science

# Nonlocal damage based failure models, extraction of crack opening and transition to fracture

Frédéric Dufour, Gilles Pijaudier-Cabot, Grégory Legrain

► To cite this version:

Frédéric Dufour, Gilles Pijaudier-Cabot, Grégory Legrain. Nonlocal damage based failure models, extraction of crack opening and transition to fracture. EURO-C 2010 - Computational Modelling of Concrete Structures, Mar 2010, Schladming, Austria. pp.301-308. hal-00867866

HAL Id: hal-00867866

https://hal.science/hal-00867866

Submitted on 29 Oct 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Nonlocal damage based failure models, extraction of crack opening and
transition to fracture

Frédéric Dufour

3S-R Lab. (Soils, Solids, Structures - Risks), Grenoble Institute of Technology, France

Gilles Pijaudier-Cabot

Laboratoire des Fluides Complexes, Univ. de Pau et des Pays de l'Adour, France

Grégory Legrain

GeM Research Institute, Ecole Centrale de Nantes, France

ABSTRACT: Damage models are capable to represent initiation and somehow crack propagation in a contin-
uum framework. Thus crack openings are not explicitly described. However for concrete structures durability
analysis, crack opening through transfer properties is a key issue. Therefore, in this contribution we present a
new approach that is able from a continuum modelling to locate a crack from internal variable field and then
to estimate crack opening along its path. Results compared to experimental measures for a three point bending
test are in a good agreement with an error lower than 10% for widely opened crack (40µm).

# 1 INTRODUCTION

For many concrete structures, crack opening is a
key parameter needed in order to estimate durabil-
ity. Cracks are preferential paths along which flu-
ids or corrosive chemical species may penetrate in-
side concrete structural elements. For structures such
as confinement vessels, reservoirs or nuclear waste
disposals for instance, tightness to gas or liquids is
a major serviceability criterion that is governed by
Darcy's relation in which permeability of the mate-
rial is involved. Hence, the prediction of the durabil-
ity of structural components requires models that de-
scribe failure, crack locations and crack openings in
the present example too when damage has localised.

Enhanced continuum and integral damage mod-

els are capable of representing diffuse damage, crack
initiation and possibly crack propagation (Pijaudier-
Cabot and Bažant 1987; Peerlings et al. 1996). They
regard cracking as an ultimate consequence of a grad-
ual loss of material integrity. These models, however,
do not predict crack opening as they rely on a contin-
uum approach to fracture.

Ideally, the prediction of durability that involves
inception of failure, crack location, propagation and
crack opening would require to merge the continuum
damage approach and the discrete crack approach
(for instance, the cohesive crack model (Hillerborg
et al. 1976)) into a single. Bridges between damage
and fracture have been devised in the literature (see
e.g. (Mazars and Pijaudier-Cabot 1996)). They rely
on the equivalence between the dissipation of energy due to damage and the energy dissipated in order to propagate a crack. Given the energy dissipated in the damage process, the equivalent crack length is computed, knowing the fracture energy. Generally, the entire energy that is dissipated in the fracture process zone is "converted" into a crack length (Mazars and Pijaudier-Cabot 1996). Some part of this energy may be dissipated in the process zone outside from the crack and it follows that the crack length and opening are probably overestimated.

The strong discontinuity approach initiated by (Simo et al. 1993) offers the possibility of merging in the same formulation a continuous damage model for the bulk response and a cohesive model for the discontinuous part of the kinematics. It is certainly a combination of continuum - discrete modelling that is sound from a theoretical point of view and appealing from the point of view of the physics of fracture. The issue in combining the continuum based model for crack initiation and then a discrete crack model for propagation is, however, the threshold upon which one switches from one analysis to the other. Usually, it is considered that the discontinuity appears when damage, stresses or strain energy reach a certain threshold fixed beforehand, which remains arbitrary. Besides as damage and fracture models do not rely on the same material description and thus on the same internal variables, jumps in time are observed on variables of interest (strain and stress) at the switch time.

As we will see further, one of the outcome of the present paper is to provide an indicator on the basis of which the appearance of a discontinuity during a damage process can be defined, with an indication of accuracy. Instead of trying to combine continuum and discrete models in computational analyses, it would be attractive to derive from the continuum approach an estimate of crack opening, without considering the explicit description of a discontinuous displacement field in the computational model. This derivation could be based on some post-processing of the distribution of strain and damage in the considered structure. The main purpose of this paper is to present such an estimate of crack opening derived from a continuum model description.

First, we recall the continuum approach that will be considered: the (integral) nonlocal damage model. The location of the crack in the computational domain and the estimate of its opening are discussed in the second part in which we propose an improvement of an existing approach (Dufour et al. 2008). Finally, we compare our numerical procedure with experimental results obtained on a 3 point bending test on plain concrete beam.

## 2 NUMERICAL MODELLING

### 2.1 Nonlocal damage approach

The scalar isotropic damage model (Mazars and Pijaudier-Cabot 1989) will be used in the FE computations for representing the progressive failure. This constitutive relation exhibits strain softening. Thus a regularization technique shall be considered in order to avoid mesh dependency and ill-posedness of the governing equations of equilibrium. In this model the tensorial stress $\sigma$ - strain $\varepsilon$ relationship is expressed as follows:

$$ \sigma = (1 - D)C : \varepsilon \quad (1) $$

where D is the damage scalar variable and C is the elastic stiffness tensor of the sound material. Damage is a combination of two components: $D_t$ and $D_c$ which are damages due to tension and compression based loads respectively:

$$ D = \alpha_t D_t + \alpha_c D_c \quad (2) $$

$\alpha_t$ and $\alpha_c$ depend on both strain and stress tensors. Damage evolution laws for both traction and compression components read:

$$ D_{t,c} = 1 - \frac{Y_{D0}(1 - A_{t,c})}{\bar{Y}} - \frac{A_{t,c}}{e^{[B_{t,c}(\bar{Y} - Y_{D0})]}} \quad (3) $$

where $A_t$, $A_c$, $B_t$, $B_c$ and $Y_{D0}$ are model parameters and $\bar{Y}$ is defined by:

$$ \bar{Y} = \max(\bar{Y}, \bar{\varepsilon}_{eq}) \quad (4) $$

with $\bar{Y} = Y_{D0}$ initially. The nonlocal equivalent strain $\bar{\varepsilon}_{eq}$ (Pijaudier-Cabot and Bažant 1987) is defined as a
weighted average of the local equivalent strain $\varepsilon_{eq}$:

$$ \bar{\varepsilon}_{eq}(\boldsymbol{x}) = \frac{\int_{\Omega} \phi(\boldsymbol{x} - \boldsymbol{s}) \varepsilon_{eq}(\boldsymbol{s}) d\boldsymbol{s}}{\int_{\Omega} \phi(\boldsymbol{x} - \boldsymbol{s}) d\boldsymbol{s}} \quad (5) $$

Several weight functions exist in the literature, we choose the most used, i.e. the Gaussian function:

$$ \phi(\boldsymbol{x} - \boldsymbol{s}) = \exp\left(-\left(\frac{2\|\boldsymbol{x} - \boldsymbol{s}\|}{l_c}\right)^2\right) \quad (6) $$

where $l_c$ is the internal length of the model. Finally the local equivalent strain is defined according to Mazars criterion:

$$ \varepsilon_{eq} = \sqrt{\sum_{i=1}^{3} \langle \varepsilon_i \rangle_{+}^2} \quad (7) $$

$\langle \rangle_+$ denotes the positive part of the principal strain $\varepsilon_i$.

## 2.2 Location of a crack

In (Dufour et al. 2008), the extraction procedure of the crack opening supposed a-priori known the crack position and the computational domain was reduced to 1D. In order to extend this approach in a more general context (2D and 3D with unknown crack position), it is necessary to locate an idealized crack from the non-local computation.

Some approaches have already been proposed in the field of damage/fracture transition in order to update the crack position during the propagation. In (Comi et al. 2007), the authors proposed to fit a fourth-order polynomial on the damage field, then to propagate the crack in a direction that is perpendicular to the maximum curvature of the polynomial at the crack-tip. The main drawback is that when the damage profile does not exhibit a clear peak but a region with a small curvature, the fitting may be obtained with a degraded accuracy. Moreover, the accuracy of the fitting may not be sufficient at the crack-tip, which should lead to extra difficulties for the estimation of the crack direction.

In (Mariani and Perego 2003), the authors proposed a similar procedure, but working on the stress-field in a half disc centered at the crack-tip. Finally, the

crack is introduced perpendicular to the fitted maximal principal stress. This approach allows to work on a “sharper” mechanical field, but the influence of the degree of the polynomial fitting was not discussed. However, the authors reported that a third order polynomial fit was not sufficient, and that in the proposed examples a fourth order one provided consistent results.

A last approach was proposed by (Oliver and Huespe 2004a), called “Global tracking algorithm”. This approach was first used within the strong discontinuity approach (SDA) (Oliver and Huespe 2004b) to evaluate the crack propagation direction. The resolution of a heat-conduction like problem leads to a scalar function whose iso-values represent all the possible directions of propagation. The selection of the iso-value emanating from the crack-tip makes possible its propagation. The approach has been modified (Feist and Hofstetter 2007) in order to restrict the heat-conduction problem on a subset of elements already or potentially crossed by the crack.

Here, we propose to use this approach in order not to propagate the crack since we use a continuum modeling but rather to locate it from mechanical variables at hand.

**Global tracking algorithm** According to this approach, the evaluation of the propagation direction is obtained as a separate problem (linked to the mechanical one). The crack is assumed to be located along a surface (or a line in 2D) which is tangent to a vector field $\vec{T}(\boldsymbol{x})$ (with unit norm). The construction of the envelopes of $\vec{T}(\boldsymbol{x})$ supplies all the possible discontinuity lines at time $t$ (see Figure 1).

The envelopes of $\vec{T}(\boldsymbol{x})$ are described by a function $\theta(\boldsymbol{x})$ whose level contours ($\theta(\boldsymbol{x}) = \text{constant}$) define all the possible discontinuity lines, as described in Figure 1. The gradient of this function must be normal to $\vec{T}(\boldsymbol{x})$ in each point:

$$ \vec{T}(\boldsymbol{x}) \cdot \text{Grad} \theta = \frac{\partial \theta}{\partial \vec{T}} = 0 \quad (8) $$

This condition can be formulated as the following linear boundary value problem (Oliver and Huespe
<figure><img src="image_2.png" /><figcaption>Figure 1: Global tracking algorithm: envelopes of the vector field T(x), possible crack path and real crack.</figcaption></figure>

2004a):

$$ \mathrm{div}(\mathbf{K} \mathrm{Grad} \theta) = 0 \quad \text{in } (\Omega) $$

$$ (-\mathbf{K} \mathrm{Grad} \theta) \cdot \vec{n} = 0 \quad \text{on } \partial_q \Omega \quad (9) $$

$$ \theta = \theta_d \quad \text{on } \partial_{\theta} \Omega $$

where (Ω) is the domain occupied by the solid, $\vec{n}$ is the unit vector normal to $\partial_q \Omega$, $\theta_d$ is a prescribed value for the Dirichlet boundary condition and $\mathbf{K}$ is a second order tensor defined as:

$$ \mathbf{K}(\mathbf{x}) = \vec{T}(\mathbf{x}) \otimes \vec{T}(\mathbf{x}) \quad (10) $$

The $\theta$ field can be assimilated as a temperature field, $-\mathbf{K} \mathrm{Grad} \theta$ as a heat flux, and $\mathbf{K}$ as an anisotropic conductivity tensor. If the Dirichlet boundary conditions are compatible with Equation (8), then a solution satisfying:

$$ \theta(\mathbf{x}) \neq \text{constant} ; \quad \frac{\partial \theta}{\partial \vec{T}} = 0 \quad (11) $$

is solution of the boundary value problem presented in Equation (9). In practice, the temperature is arbitrarily defined at a single node in order to avoid loss of uniqueness of the thermal problem. In order to overcome the singularity of the problem ($\mathbf{K}$ is rank one), the conductivity tensor is modified as (Oliver and Huespe 2004a):

$$ \mathbf{K}(\mathbf{x}) = \vec{T}(\mathbf{x}) \otimes \vec{T}(\mathbf{x}) + \epsilon \mathbf{I} \quad (12) $$

<figure><img src="image_3.png" /><figcaption>Figure 2: Algorithm for the location of the idealized crack</figcaption></figure>

where $\epsilon$ is a small isotropic algorithmic conductivity, and $\mathbf{I}$ is the second order identity tensor. Once the problem is solved, the crack can be propagated along the path defined by the iso-value of $\theta$ that passes at the crack-tip.

**Location of the crack using the global tracking algorithm** To apply this approach to the problem at hand, two main ingredients have to be adapted: (1) the definition of the $\vec{T}$ field, and (2) the location of one point of the crack. We make here the hypothesis that the idealized crack is perpendicular to the principal direction associated to the maximum principal strain $\epsilon_{max}$ which represents the opening direction in mode I dominated loading. The $\vec{T}$ field is thus taken perpendicular to the principal direction associated to $\epsilon_{max}$. The knowledge of this field in the body makes possible to solve the boundary value problem and obtain the $\theta$ field. The last operation consists in selecting the right iso-value. We make here a last hypothesis by considering that the crack passes at the Gauss point where $\epsilon_{max}$ is maximal on the body. The algorithm that summarizes the process is presented in Figure 2. In practice, the thermal-like problem is not solved on the full structure, but only on the damaged zone. This allows to speed-up the process and decrease the computer requirements.
<figure><img src="image_4.png" /><figcaption>Figure 3: 1D profiles (blue) generated from the crack mesh (red)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: 1D profiles (blue) generated from the crack mesh (red)</figcaption></figure>

## 2.3 Reduction to a 1D problem

At the end of the tracking process, a mesh of the crack is built using the iso-temperature defining the crack. The second step consists now in evaluating the opening across this idealized crack. In this contribution, it is proposed to re-use the approach that was presented in (Dufour et al. 2008). Once the crack is meshed, it is possible to apply the 1D approach on lines perpendicular to the elements (segments in 2D, triangles in 3D) defining the crack surface. In 2D, for example, a set of lines is generated from the middle of each segment of the crack (see Figure 3).

Once these profiles are defined, the component of the local strain field along the 1D profile $\varepsilon_N = \vec{N} \cdot \varepsilon$. $\vec{N}$ is first computed (see Figure 4). Then, this axial strain field is projected on the 1D profile as an input for the 1D crack opening procedure.

## 2.4 Estimation of a crack opening

We summarize in this part the key idea developed by (Dufour et al. 2008) to estimate the crack opening in a 1D structure that we use along perpendicular profiles to the idealized crack. If we assume a bar upon failure, the displacement field is a step with a jump $\{U\}$ at the

crack location $x_0$. The derivation of the displacement field gives a Dirac function for the local strain and a nonlocal strain with an amplitude of $\{U\}$ and the same shape than the averaging function $\psi$ used in the convolution product. This nonlocal strain is denoted as $\bar{\varepsilon}_{sd}$.

*Remark:* For any regularized damage models $\psi$ can be defined independently of the mechanical model. However, since with the nonlocal model used in the present work we already have defined a weighting function $\phi$, we keep it, thus $\psi = \phi$ (see Equation (6)).

With this procedure a nonlocal measure of strain is analytically obtained assuming a strong discontinuity kinematical field upon failure. The key point of (Dufour et al. 2008) is to compare this function with the nonlocal strain obtained by the FE mechanical computations. Several possibilities do exist in order to compare this strain profile to the nonlocal strain obtained from the FE calculation. In the original paper, only the strong link were developed, i.e. the crack opening is computed so that both profiles are equal at their maximum $x = x_0$, this yields:

$$ \{U\} = \frac{\bar{\varepsilon}_{eq}(x_0) \int_{\Gamma} \psi(x_0 - s) ds}{\psi(0)} \quad (13) $$

Furthermore, the distance between both profiles gives an indicator of the quality of the solution obtained by the FE computation using nonlocal damage model with respect to an analytical strong discontinuity approach:

$$ \varepsilon_I(x) = \frac{\int_{\Gamma} |\bar{\varepsilon}_{sd} - \bar{\varepsilon}_{eq}| d\Gamma}{\int_{\Gamma} \bar{\varepsilon}_{eq} d\Gamma} \quad (14) $$

Thus, it is not an error on the crack opening itself but on the capacity of nonlocal damage to reproduce local kinematical field across the crack as in the strong discontinuity approach. The strain profile width is related to the internal length of the model, the crack opening corresponds to some integration over this profile and it depends on the internal length (Giry et al. 2010).

A new profile comparison technique, named weak form, is proposed in the present work by equating the
integral of both profiles, i.e.:

$$ \int_{\Gamma} \bar{\epsilon}_{eq} d\Gamma = \int_{\Gamma} \bar{\epsilon}_{sd} d\Gamma \quad (15) $$

Thus it gives a different value for the crack opening. Since this error measure gives only a quality estimation of the model, we have performed experimental test in order to estimate by comparison the error on the crack opening itself.

## 2.5 Comparison vs experimental results

In order to quantify our approach against experimental measurements, three point bending test were performed on a notched beam. The beam dimensions are 400 mm span, 100 mm high, 50 mm thick and the notch is 20 mm high. The test is driven by the CMOD measure at the notch mouth.

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Initial picture in undeformed state, (b) picture during the crack propagation, (c) horizontal displacement field and its 3D view (d).</figcaption></figure>

In order to measure the crack length and opening,

we use a Digital Image Correlation technique (see Figure 5). For practical reasons, the picture frame is limited to 55 mm high from the notch tip. The crack is assumed to be vertical and thus the crack opening is estimated as the horizontal displacement jump. In order to get the crack opening evolution along the crack, 30 horizontal profiles are drawn and the displacement jumps are estimated along those profiles. For particular values of CMOD 20 (corresponding to peak load), 30, 40, 50, 60, 80, 100, 150 and 200 µm, a linear curve is fitted through the 30 measurements (see Figure 6). The fitted solid lines are prolonged in dashed line up to the CMOD value.

<figure><img src="image_7.png" /><figcaption>Figure 6: Crack opening at different stage of the loading process.</figcaption></figure>

For the numerical simulation, we use the nonlocal version of Mazars’ damage model described in 2.1. Model parameters are fitted (before any crack opening estimation) on the experimental global response, i.e. force vs CMOD. Comparison between experimental and numerical curves is shown in Figure 7. A good fit is obtained for material parameters summarized in Table 1.

<table><caption>Table 1: Parameter fitting using nonlocal damage model</caption><thead><tr><th>E [GPa]</th><th>ν</th><th>A<sub>t</sub></th><th>B<sub>t</sub></th><th>Y<sub>D0</sub></th><th>l<sub>c</sub> [mm]</th></tr></thead><tbody><tr><td>30</td><td>0.2</td><td>0.9</td><td>4 000</td><td>410<sup>-5</sup></td><td>8</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 7: Experimental and numerical Force vs CMOD responses.</figcaption></figure>

For a given CMOD, the crack shape is compared between experimental measurements and numerical estimation both using the strong and the weak link approaches. A relative error is computed between experimental crack opening and its numerical counterpart. Just after the peak (CMOD = 50 µm), the two numerical approaches are quite similar (see Figure 8.a) and slightly underestimate the measured crack opening. However for large CMOD (200 µm) the strong approach yield a large error (see Figure 8.b and d) and the weak approach always provides a better estimation of the measured crack opening.

The strong approach relies only on the regularized equivalent strain at one given point that may be affected by boundary effect for instance (Pijaudier-Cabot et al. 2009) and is thus more sensitive to numerical perturbations.

The larger the crack opening, the better the estimation. This is a rather important result since the transfer properties for a structure are naturally dominated by large crack openings.

The numerical approach systematically underestimates the experimental crack opening, at least for a 3 point bending test. Although it is not on the safe side for an engineering use, it can be clearly explained from crack propagation considerations and by recalling that experimental crack opening are measured on the surface whereas the numerical one is performed on a 2D plane stress simulation:

* The stress state is close to a plane stress condi-

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison between strong and weak approaches vs experimental crack opening for COD = 50 µm (a) and 200 µm (b). Corresponding errors between numerical and experimental crack openings for COD = 50 µm (c) and 200 µm (d).</figcaption></figure>

tion at the beam free surfaces whereas it is close to a plane strain condition in the bulk of the beam that reduces the crack propagation velocity due to confinement.

* Due to casting process the material contains less aggregate close to the boundaries and is thus weaker in the sense that aggregates are obstacles for cracking.

For these two reasons, on the surface the crack is more developed in length and opening than in the core of the beam. It is clearly proved if one looks carefully at the experimental measurements of the crack opening for CMOD of 200 µm (see Figure 8-b). The extension of the plot gives a zero opening above the top of the beam, i.e. the neutral axis is out of the beam. For a bending test it means that the applied load is null. However in Fig 7 one can see that for CMOD of 200 µm the bearing capacity is not yet zero and thus outside the surface the crack has not yet propagated
to the beam limits. Besides due to 2D assumption the numerical modeling gives an average crack geometry between the surface and the core of the beam. Note also that part of the inaccuracy in the crack opening estimation is due to the spreading of strain profile that occurs in nonlocal damage models.

## 3 CONCLUSIONS

We have presented a complete procedure used in a post-treatment analysis to get crack location, crack opening and an estimation of the error done for tensile failure. The tracking is perform solving a conduction-like FE problem based on mechanical variables. The crack opening is estimated by nonlocal strain profile comparisons with those analytically obtained from the strong discontinuity approach. Results are in good agreement with crack opening measured on a 3 point bending test by DIC technique.

## REFERENCES

* Comi, C., S. Mariani, and U. Perego (2007). An extended fe strategy for transition from continuum damage to mode i cohesive crack propagation. *International Journal for Numerical and Analytical Methods in Geomechanics* 31(2), 213–238.

* Dufour, F., G. Pijaudier-Cabot, M. Choinska, and A. Huerta (2008). Extraction of a crack opening from a continuous approach using regularized damage models. *Computers and Concrete* 5(4), 375–388.

* Feist, C. and G. Hofstetter (2007). Three-dimensional fracture simulations based on the sda. *International Journal for Numerical and Analytical Methods in Geomechanics* 31(2), 189–212.

* Giry, C., F. Dufour, J. Mazars, and P. Kotronis (2010). Stress state influence on nonlocal interactions in damage modelling. In *Computational Modelling of Concrete Structures*. Rohrmoos-Shladming, Austria. Euro-C.

* Hillerborg, A., M. Modeer, and P. E. Pertersson (1976). Analysis of crack formation and crack

growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research* 6, 773–782.

* Mariani, S. and U. Perego (2003). Extended finite element method for quasi-brittle fracture. *International Journal for Numerical Methods in Engineering* 58(1), 103–128.

* Mazars, J. and G. Pijaudier-Cabot (1989). Continuum damage theory - application to concrete. *Journal of Engineering Mechanics* 115(2), 345–365.

* Mazars, J. and G. Pijaudier-Cabot (1996). From damage to fracture mechanics and conversely: a combined approach. *International Journal of Solids and Structures* 33, 3327–3342.

* Oliver, J. and A. E. Huespe (2004a). Continuum approach to material failure in strong discontinuity settings. *Computer Methods in Applied Mechanics and Engineering* 193(30-32), 3195–3220.

* Oliver, J. and A. E. Huespe (2004b). Theoretical and computational issues in modelling material failure in strong discontinuity scenarios. *Computer Methods in Applied Mechanics and Engineering* 193(27-29), 2987–3014.

* Peerlings, R. H. J., R. de Borst, W. A. M. Brekelmans, and J. H. P. de Vree (1996). Gradient enhanced damage for quasi-brittle materials. *International Journal for Numerical Methods in Engineering* 39, 937–953.

* Pijaudier-Cabot, G. and Z. Bažant (1987). Nonlocal damage theory. *Journal of Engineering Mechanics* 113, 1512–1533.

* Pijaudier-Cabot, G., A. Krayani, and F. Dufour (2009). Boundary effect on weight function in nonlocal damage model. *Engineering Fracture Mechanics* 76(14), 2217–2231.

* Simo, J. C., J. Oliver, and F. Armero (1993). An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. *Computational Mechanics* 12, 277–296.