Stabilized mixed finite elements with embedded strong discontinuities for shear band modelling

A stabilized mixed finite element with elemental embedded strong discontinuities for shear band modelling is presented. The discrete constitutive model, representing the cohesive forces acting across the shear band, is derived from a rate-independent J2 plastic continuum material model with strain softening, by using a projectiontype procedure determined by the Continuum-Strong Discontinuity Approach. The numerical examples emphasize the increase of the numerical solution accuracy obtained with the present strategy as compared with alternative procedures using linear triangles.


Introduction
Shear bands in plastic solids arise as a typical deformation mode due to a strain localization phenomenon, during the inelastic deformation processes, when the material becomes unstable. Numerical modelling of shear bands using discontinuous velocity fields was previously proposed by several authors (see for example Armero et al. [1], Regueiro et al. [2] and Samaniego et al. [3]).
In a number of problems, a large stable process of irreversible isochoric plastic deformation precedes the inception of a strain localized mode. In these cases, and from a computational point of view, it should be considered the deficient response provided by standard finite elements when kinematics incompressibility constraints are present. This particular aspect of the numerical approach is a classical, and extensively studied, issue in computational mechanics (see Zienkiewicz et al. [4], Hughes [5]).
In this paper we present a stabilized mixed finite element formulation, which has been recently developed for J 2 plasticity [6,7]. The kinematics is enriched with the addition of an embedded strong discontinuity mode with elemental support, like that proposed in Oliver [8,9], for capturing the characteristic shear band type deformation mechanisms. The idea of using a well behaved finite element for plasticity in conjunction with an embedded strong discontinuity kinematics is not new in shear band modelling. Armero et al. [1] have used a triangular MINI element and Regueiro et al. [2] the classical quadrilateral BBAR element, both of them enriched with an embedded strong discontinuity. Nevertheless, the authors understand that the problem remains open since, in their opinion, the linear triangle has a number of advantages which make it particularly suitable to be enriched with embedded discontinuities. The stabilized element here presented is a linear triangle.
The Continuum-Strong Discontinuity Approach [10] adopted in this work, determines the shear strain rate-traction rate separation law of the shear band. A characteristic of this procedure is that the resulting discrete law governing the shear band evolution, i.e. the cohesive force acting across the shear band surface, is a projection onto the discontinuity surface of the bulk material constitutive model. In this work, the non-localized (bulk) material behavior follows a rate-independent J 2 elastoplastic law with strain softening response.
Alternative models for simulating shear bands have been numerous in the past. Recently, Cervera et al. [11] have presented a model that uses the same stabilized mixed finite element shown here, but without introducing the embedded strong discontinuity mode into the finite element. However, the authors think that the additional features provided by the CSDA deserve to be studied too.
The paper proceeds as follows: in "Section 2", we present the enriched kinematics with the strong discontinuity mode and the discrete constitutive model governing the shear band evolution. "Section 3" presents the finite element formulation with the stabilization procedure and "Section 4" its numerical implementation. In the present work, we are interested in the analysis of the numerical stabilization effect, its influence on the shear-band capturing and the subsequent post-critical response, particularly when embedded strong discontinuities are used. This analysis is presented in "Section 5" by means of two numerical applications. In the first case, a slope instability problem, we compare the numerical response obtained by different finite element implementations, including standard and stabilized mixed linear triangles with and without embedded strong discontinuities, quadrilaterals, etc. Also, we analyze the convergence rate of the solution with the finite element mesh size. In the second example, the near incompressibility constraint is imposed already at the beginning of elastic regime. In this context, again we study the ability of the model to capture the shear band and the obtained peak load is compared with an analytical solution taken from the literature. Finally, the conclusions are presented.

Strong discontinuity kinematics
Let Ω be a body which experiences a shear band failure mode. The material surface S, with normal n intersecting the body Ω, represents the zone with localized strain rate, as it is shown in " Fig. 1". The appropriate kinematics describing this phenomenon should account for a discontinuous velocity field across S, such as the following one: whereu(x, t) represents a continuum field, H S (x) is the Heaviside's step function shifted to S (H S (x) = 1 ∀x ∈ Ω + and H S (x) = 0 ∀x ∈ Ω − ), that, multiplied by the velocity jump vectoṙ β, introduces the discontinuity term into the velocity field. The infinitesimal strain rate that is compatible with this velocity field, is a generalized func- composed of a regular termε = ∇ symu + H S (∇ symβ ) and a singular one, given by the Dirac's delta function (δ S ) shifted to S. The boundary value problem (BVP) of a quasi-static elastoplastic body showing a strong discontinuity kinematics, such as a shear band, is described (in rate form) by the following equations: where the Cauchy's equation (3), relating the stress rateσ with the rate of volumetric forces ρḃ, and ignoring the inertial effects, is defined in the regular part of the body (Ω/S), i.e. the points in Ω excluding those in S and where no strain rate localization effects are observed. The boundary condition in velocitiesu * and rate of tractionsṫ * are imposed on Γ u and Γ σ ("Eq. (4)" and "Eq. (5)") respectively. Furthermore, the equilibrium condition across the discontinuity surface S requires that:ṫ + =σ + · n =σ − · n =ṫ − ∀x ∈ S whereṫ + (ṫ − ) is the traction vector applied to the body part Ω + (or Ω − ) on the boundary S. If cohesive tractions (t S ) are considered in the shear band interface, the equilibrium condition in rates also requires:ṫ S =σ S · n =σ + · n =ṫ + ∀x ∈ S where, and consistently with the Continuum-Strong Discontinuity Approach, a fundamental hypothesis has been adopted: a stress state σ S exists into the discontinuity zone S (where singular strain rates are present), which is defined by a regularized version of the constitutive model that describes the regular part of the body Ω/S, see [8,12,13].

Continuum constitutive model and discrete cohesive law
We assume for the Ω/S domain a rate-independent J 2 elastoplastic material model with strain softening described by the equations: x 1 x G Figure 1: Strong discontinuity problem.
where C is the fourth order elastic constitutive tensor depending on the Lamé's parameters (λ and µ), with I 1 and I being the second and fourth order unit tensors respectively,ε p is the plastic strain rate tensor, q and α are scalar internal variables and φ is the yield function describing the elastic domain evolution depending on the deviatoric stress tensor S = dev(σ) (through the second invariant J 2 ) and the yield strength σ y . We denote M the plastic deviatoric strain rate direction (being tr(M) = 0) and γ the plastic multiplier. From "Eq. (10-12)", α is identified as the total equivalent plastic strain. Special attention should be paid, in the present setting, to the softening modulus H (H < 0), which plays a main role in the localization condition.
In the Continuum-Strong Discontinuity Approach, followed in the present work, it is assumed that the stress σ S is determined by a regularized version of the model given by "Eq. (8)(9)(10)(11)(12)(13)". This stress state, which due to equilibrium conditions must be a bounded tensor, defines the cohesive behavior of the interface S.
Following Simo et al. [12] and Oliver [10], and considering the regularized sequence of functions δ S = lim h→0 µs h (where µ S (x ∈ S) = 1, µ S (x / ∈ S) = 0), it can be shown that (variables with subindex (·) S are referred to their evaluation at the domain S): is a bounded term whenever: Condition (15) can be verified by introducing a singular measures for the plastic multiplier γ and the inverse of the softening modulus H: whereH is an intrinsic softening modulus, determined by the material fracture energy G f . Therefore, from "Eq. (11)" and "Eq. (16)",q becomes a regular term, even whenε S is singular: and replacing "Eq. (9)" into "Eq. (14)", yields: which has been termed the "strong discontinuity equation" [10]. Recalling that: and given the particular structure of tensor (n ⊗β) sym , "Eq. (18)" imposes the strong discontinuity condition on σ S , which establishes that σ S is only characterized by the traction vector t S , see " Fig. 2". Additional details on this aspect can be found in [10].
Then, the yield function and the consistency equation, in a loading process, can be written as follows: We remark that t dev S is normally termed the Schmidt resolved shear stresses for the slip plane S, see [14].
In a loading state (γ > 0), "Eq. (18)" and "Eq. (20-b)" determine the velocity jumpβ: Equation (21) is consistent with a classical constitutive assumption on the slip phenomenon in single-crystal plasticity: the shear rate (ζ = β ), in a slip system, depends on the stresses only through the Schimdt resolved shear stress (t dev S ). Implicit in "Eq. (21)" is the fact that the velocity jumpβ is compatible with a slip line mode (β · n = 0) and that: where C ep is the perfectly-plastic constitutive tensor and Q ep is the so-called "localization tensor". The degenerated (projected) cohesive model, traction-separation law, derived from the continuum model and induced by a strong discontinuity kinematics, is displayed in Box 1.

Stabilized mixed variational formulation using embedded strong discontinuities
Decomposing the stress rate into its deviatoricṠ, and spherical −ṗ I 1 (ṗ = − 1 3 tr(σ)), parts: and considering from the constitutive model thatṗ = −κ∇ ·u, where κ is the volumetric modulus, the BVP can be set within a classical variational mixed (velocity, pressure) format: findu ∈ V u and p ∈ Q such that: The admissible functional space for q is Q ≡ L 2 (Ω/S) . We define the space of admissible functions for velocities V u by assuming the existence of non-smooth terms representing the velocity jumps developed in the shear band zone. These terms are included via the embedded strong discontinuity technique. Let the velocity space V u be defined by: where M S is the so-called elemental unit jump function [9], whose support is a given domain Ω h that includes S. The ϕ(x) term can be taken as an arbitrary smooth function such that: ϕ(x ∈ Ω + ) = 1 and ϕ(x ∈ Ω − ) = 0. Alsoβ ∈ IR dim , with dim standing for the space dimension, is the velocity jump vector. The virtual (kinematically admissible) velocities lie on the space: It should be mentioned thatu andη are smooth functions (V u ⊂ H 1 ). Introducing the spaces (29) and (31) into "Eq. (28)", the governing equations can be alternatively written as follows: where P (ext) u is the virtual power of the body forces and external loads. Recalling that ∇M S = (δ S n − ∇ϕ), then "Eq. (32-c)" imposes a weak traction continuity condition on the discontinuity surface and it can be rewritten as: where we have identified the matrixG S with the operator ∇M S applied to stresses. A widely used variational non-symmetric (not-consistent) formulation, redefines the weak traction continuity "Eq. (33)" by exchanging ∇ϕ by n, in the second left hand side term, and computing the mean values of the traction continuity [15]: where l S is the length of discontinuity S intersecting the finite element, see " Fig. 4". In the numerical examples, we present solutions considering both procedures. We use the term "symmetric formulation" when "Eq. (33)" is implemented and "non-symmetric formulation" if condition "Eq. (34)" governs the traction continuity.

Stabilization
It is well known that mixed formulations like "Eq. (28)" suffer from numerical instability issues [4,5]. The instability problem becomes particularly serious when piecewise linear polynomial functions of continuity C 0 are chosen for interpolation of both spacesV u and Q, because in that case the so-called Ladyzhenskya-Babuska-Brezzi condition (or simply LBB) is not satisfied [16]. A remedy for this unwanted effect has been the introduction of stabilization terms S st into the variational principle "Eq. (28)". Particularly, this term is added to the left hand side of "Eq.
The stabilization term used in this work has been introduced by Codina [17] in the fluid mechanics context and extended by Cervera et al. [6] to J 2 -plasticity problems. It has been termed the orthogonal sub-scale method, PGP, and is defined by: whereΠ(∈V u ) is the projection-L 2 (Ω/S) of the discrete pressure rate gradient (∇ṗ) on the regular finite element approximation space (V u ), see " Fig. 3": This procedure considers the term S st proportional to a stabilization factor τ , depending on the shear modulus µ and a characteristic finite element size h (we have adopted h to be the square root of the finite element area): where the scalar coefficient c is a constant parameter (c ≈ O(1)).
Introducing the stabilization term (35) into the variational equation (28), and considering that "Eq. (36)" shall be included as an additional restriction, it is possible to rewrite the variational principle, in (u,β,ṗ,Π), as follows: "Equations (38)-(a-d)" at time t, can be alternatively written in terms of the total stresses and displacements.

Numerical implementation
Considering Ω ∈ IR 2 , simplicial finite elements (linear triangles) with C 0 piecewise linear interpolation polynomials for pressure and regular displacement fields have been chosen for the present implementation.

Displacement field approximation
The continuous part of the displacementū = {ū x ,ū y } T is interpolated in the standard way by using piecewise linear shape functions N e u (x) (supra index (·) e refers to element e). The elemental unit jump function M S e (x) = H S e (x) − (N e u ) node+ (x) is built by using the linear shape function (N e u ) node+ corresponding to that nodes belonging to the Ω + region, see " Fig. 4". The support of M S e is, therefore, one element: where(.) refers to nodal values. The strains, in a vectorial format (ε e = ε e x , ε e y , ε e xy T ), can be written as follows: where B e = (∇N e u ) sym is the strain-displacement matrix and G e is the matrix given by:

Interpolations of the pressure and L 2 -projected pressure gradient fields
The L 2 -projected pressure gradient field (Π e ) is interpolated by using identical shape functions to those chosen for the velocity approximation. In the same way, the pressure is also interpolated by means of C 0 piecewise linear functions: where N e p are, again, the classical linear shape functions.

Discrete equations. Internal force evaluation
The discrete version of the variational principle "Eq. (38)" can be formulated as follows: findû, p,Π and β such that they verify the essential boundary condition "Eq. (4)" and the following system of equations: where the internal F (int) and external F (ext) generalized forces are defined as: A being the finite element assembling operator, and matrices G 0 , M p , M u , L, H and Q are computed as follows: Implicitly, "Eq. (43)-(45)" introduce the strategy of assuming the uncoupling of the fieldΠ. Its value at the end of step n (Π (n) ), that is determined by using "Eq (38)-c" with the previously known variablep (n) : is used for solving the system "Eq. (43)-(45)" at step n + 1. This strategy has been previously utilized by Codina et al. and Chiumenti et al. [17][18][19], allowing for a more efficient computational treatment of the problem.
Following the same integration procedure presented in Oliver [9], one additional Gauss point is considered for evaluation of strains and stresses at S. Thus, integrals on S in "Eq. (44)", are referred to terms evaluated in those additional Gauss point multiplied by an adequate weight.

Tangent matrix.
The use of the Newton-Raphson scheme for solving "Eq. (43)-(45)" requires the evaluation of the system jacobian matrix J . Considering that X = [ûp β] T is the independent variable vector, J can be evaluated as follows:  where submatrices K i,j result:

Numerical simulations
The numerical response of the present model is analyzed by means of two bidimensional problems. Particularly, we are addressing our study to determine the ability of the numerical model for capturing the strain localization mode and the structural peak load. Also, we analyze other fundamental aspects in failure mechanics analysis under softening regime, such as the objectivity of the numerical results with independence of the mesh size and orientation.
The mathematical verification and consistency of the model is studied by comparing alternative finite element formulations, which are denoted using the nomenclature in " Table 1". As it can be seen there, the set of elements that we use for this comparison belongs either to the generalized displacement finite element formulation (the constant strain triangle STDSD in the Table and the BBAR quadrilateral element taken from Simo et al. [20]) or to the mixed (pressure-velocity) formulations including the PGP stabilization scheme. All of them, excepting the first one, are enriched with an embedded strong discontinuity kinematics with elemental support. The traction continuity condition is implemented using both procedures: the symmetric element type given by "Eq. (33)" and the non-symmetric element type given by "Eq. (34)".  In the PGP formulation without embedded strong discontinuities (denoted "smooth velocity kinematics" in the Table 1), solutions have been obtained by regularization of the softening modulus H, redefining it in accordance with : where h is the characteristic size of the element andH the intrinsic softening modulus computed as in "Eq. (16)". A comparison of the relative computational cost between PGPSD and BBARSD elements is also reported. For this purpose, it must be considered that the examples have been run in a PC equipped with a single Pentium 4 -3.0 GHz, 512 MB Ram -processor.
For all cases, a stability factor "c" near to unity (see "Eq. 37") was adopted to perform the numerical tests.

2D Slope stability problem
When undrained loading conditions are assumed, the constitutive behavior of saturated cohesive soils can be approximately modelled by an associative deviatoric plastic flow law. In this context, we use a J2 model to simulate a typical plane strain geotechnical slope stability problem and its corresponding shear band failure mode. A similar example was presented in Regueiro et al. [2] and in Oliver et al. [21] where a BBAR element with embedded strong discontinuities was used. Due to the lack, at least up to the author's knowledge, of an analytical or exact solution for this problem, the above mentioned strategies (denoted as BBARSD-N in " Table 1"), will be used as a reference solution to compare quantitative results.
The effects of including, or not, the strong discontinuity mode are particularly remarked in the present analysis. Also, the numerical stabilization influence on the solution, which is contrasted with similar formulations that do not use such strategy, is studied. respectively. Notice the particular mesh configuration that has been intentionally generated against to the expected strain localization path. This situation represents a challenge for the linear triangle kinematics. A fourth mesh, of quadrilaterals (M4 in " Fig. 6-(d)"), with element size similar to M3, is used to obtain the BBARSD-N reference solution.
" Fig. 7" shows, in gray color, the evolution of those PGPSD-N elements that are subjected to plastic loading conditions in four different stages, as the process advances along the time. It is clear from this figure how the strain localization phenomenon is developed, inducing the shear band mode.
In " Fig. 8" again we show, in gray color, those elements in the meshes M1, M2 and M3 and using the PGPSD-N approach, that are post bifurcation regime (were the strong discontinuity is active) at the end of the simulated process. We can observe that the three meshes display a qualitative agreement respect to the shear band trajectory, with a clear tendency to converge with the mesh refinement, toward a well defined curve which compares well with that reported by Regueiro et al. [2]). " Fig. 9-(a)" and " Fig. 9-(b)" compare the deformed mesh solutions obtained using the mixed stabilized formulations either without embedded strong discontinuity (PGP) or with it (PGPSD-N). In the first case, it is observed that the zone of strain localization has a pronounced trend to follow the mesh direction (mesh bias). Furthermore, the solution of the PGP procedure  presents a more diffuse deformation pattern respect to that shown by the PGPSD-N. Both effects determine a noticeable difference in the structural response (see " Fig. 9-(c)"), mainly in the limit load prediction. Next, we report the structural response in terms of load versus the vertical displacement δu curves (point A). " Fig. 10-(a)" shows these results, which correspond to the PGPSD-N element and different meshes. The M3 solution compares well with those obtained using the BBARSD-N strategy. In " Fig. 10-(b)" we plot the same results corresponding to the M3 mesh, but using different finite element formulations. The two responses obtained with the STDSD procedure reveals a locking (spurious) effect produced by the isochoric deformation constraints, overestimating the dissipated energy and peak load.  It must be observed that the PGPSD-N scheme shows a good prediction of the limit load P u , as compared with the reference solution, and also in terms of the dissipated energy during the localization process. To quantify both features, we plot in " Fig. 11" the convergence analysis of the PGPSD-N and STDSD-N solutions. " Fig. 11-(a)" displays in a logarithmic plot the linear regression curve of the dissipated energy error ( e L 2 ) as a function of the mesh size h. The relative error ( e L 2 ) of every solution S (M i) , where S (M i) is the load vs. displacement (δu) curve of the mesh M i (i = 1, 2, 3), is computed in terms of a L 2 norm as follow: where the integration parameter τ corresponds to the vertical displacement (δu) and τ max = max(δu) is the same for all cases, while is the reference solution. Similarly, " Fig. 11-(b)" displays the linear regression curve of the limit load prediction error as a function of the size mesh h. The relative error of the peak load solution is determined by means of: . From " Fig. 11" it is clearly observed a higher accuracy and convergence rate, either in limit load prediction as also in the dissipated energy, of the PGPSD-N model if compared with the STDSD-N element. Finally, the comparative computational cost for PGPSD-N element, relative to BBARSD-N formulation, is outlined in "   The analytical peak load solution for a problem with sharp crack type and considering perfect elastoplasticity, is available from Limit Analysis Theory (LAT) [22]. " Fig. 13-(a)" shows the equilibrium curves (vertical displacement vs. resultant force P) obtained for both meshes, using the PGPSD-N formulation, and also the PGP without discontinuous enriching modes. Again, in the first case, an adequate convergence can be observed with mesh refinement toward the BBARSD-N solution, and a reasonable accuracy respect to the analytical peak load solution. It must be reported that the standard triangle (STDSD-N) fails, dramatically, in simulating this near incompressible test.

Center cracked panel
In addition, the deformed mesh configuration of the PGPSD-N model, see " Fig. 13-(b)", displays the predicted collapse mechanism.

Conclusions
The main contribution of the present work is the presentation of a new simplicial finite element, called PGPSD, which appears as an improvement respect to previous models, leading to robust and accurate simulations for shear band problems induced by strain softening in plastic materials models. It has been developed within the context of the Continuum-Strong Discontinuity Approach. The proposed formulation is based on a consistent coupling of two techniques: (i) The Pressure Gradient Projection stabilization scheme (PGP).
(ii) The elemental embedded strong discontinuity kinematics. From the study reported in the above examples, we can extract the following conclusions in reference with the PGPSD model in general, and in particular with the PGPSD-N implementation: • The PGPSD-N element has shown an adequate performance when the strain localization phenomenon happens in a dominant quasi incompressible regime; • The numerical behavior of the PGPSD-N element has been proven quantitatively through a classical convergence study based on a structured mesh refinement from a reference solution taken from different sources ( [2,23]); • The PGPSD element shows an improvement in the convergence rates and diminution in the relative error magnitude in comparison with the standard (non-stabilized) enriched element (STDSD) and also with respect to the non-enriched PGP strategy. In addition, it compares very well with BBARSD formulation.
• The computational cost, for the two bidimensional cases presented in this work, seems to be reasonable (1.3 to 1.4 times greater than that obtained with the BBARSD procedure) considering that both set of d.o.f.'s, Π and β, can be decoupled and statically condensed in the numerical implementation respectively.
However, we have observed some troubles that must be remarked. The linear kinematics of the simplicial elements (triangle), before the activation of enriched modes, seems to be fairly stiff, which produces a noticeable effect on the bifurcation conditions, delaying the activation of the shear-band. This effect induces, in certain pathological mesh orientations, a serious kinematical locking. This unsolved limitation and the extension to 3D context, motivate future research works.