Strain Localization of Orthotropic Elasto–Plastic Cohesive–Frictional Materials: Analytical Results and Numerical Verification

Strain localization analysis for orthotropic-associated plasticity in cohesive–frictional materials is addressed in this work. Specifically, the localization condition is derived from Maxwell’s kinematics, the plastic flow rule and the boundedness of stress rates. The analysis is applicable to strong and regularized discontinuity settings. Expanding on previous works, the quadratic orthotropic Hoffman and Tsai–Wu models are investigated and compared to pressure insensitive and sensitive models such as von Mises, Hill and Drucker–Prager. Analytical localization angles are obtained in uniaxial tension and compression under plane stress and plane strain conditions. These are only dependent on the plastic potential adopted; ensuing, a geometrical interpretation in the stress space is offered. The analytical results are then validated by independent numerical simulations. The B-bar finite element is used to deal with the limiting incompressibility in the purely isochoric plastic flow. For a strip under vertical stretching in plane stress and plane strain as well as Prandtl’s problem of indentation by a flat rigid die in plane strain, numerical results are presented for both isotropic and orthotropic plasticity models with or without tilting angle between the material axes and the applied loading. The influence of frictional behavior is studied. In all the investigated cases, the numerical results provide compelling support to the analytical prognosis.


Introduction
Orthotropic Materials such as wood and masonry have been traditionally used in construction and are very much used today. Other frequently used materials, such as rolled metals, are orthotropic because of their manufacturing process. This is also very much the case of metallic and polymeric materials and components produced layer-by-layer using modern additive manufacturing (AM) techniques, now increasingly used. In the field of geological engineering, the analysis of orthotropic materials is of interest in ground excavation, tunnel construction and landslides prevention.
Hill [1][2][3], a pioneer in the mathematical research of plasticity, proposed several constitutive orthotropic plasticity models for sheet metals and investigated strain localization and failure of orthotropic plastic materials. Based on Hill's works, many isotropic and orthotropic plastic criteria have been later proposed, such as the Drucker-Prager model [4][5][6][7], Hoffman model [8], Tsai-Wu model [9], and many more [10][11][12]. Purely cohesive models that are insensitive to pressure and yield an isochoric plastic flow, such as the von Mises and the Hill models, are appropriate for metallic materials. Associated elasto-plastic cohesive-frictional models such as the Drucker-Prager, Hoffman and Tsai-Wu models are suitable for simulating polymeric materials such as PVC H100, H250 and carbon fiber not depend on the elastic properties or on the yield surface. The localization angles can be analytically predicted from the inelastic flow tensor alone. This paper addresses the analytical determination of the orientation of slip lines in orthotropic elasto-plastic cohesive-frictional materials by extending the strain localization analysis developed in previous works. The objectives are fourfold: (i) to extend the strain localization analysis to orthotropic elasto-plastic cohesive-frictional materials; (ii) to derive analytically localization angles in plane stress and plane strain conditions for these models; (iii) to verify these analytical results via independent numerical simulations; and (iv) to investigate the influence of plastic material properties on strain localization in orthotropic cohesive-frictional materials.
The paper is structured as follows. Section 2 briefly presents the analytical framework: constitutive relations, kinematics for strong and weak discontinuities, and strain localization conditions. Section 3 introduces orthotropic plasticity and develops the analytical results for the localization angles in plane stress and plane strain conditions with some examples. In Section 4, numerical verification of the analytical results using B-bar finite elements is offered. Section 5 closes the paper with some conclusions.

Strain Localization in Elasto-Plastic Solids
In this section, the mechanics of strain localization in elasto-plastic media is addressed. Using Maxwell's kinematics and assuming boundedness of the stress rates, the necessary condition for strain localization in elasto-plastic materials is obtained. The results hold both for strong (displacement) discontinuities and for regularized strain localization bands limited by weak (strain) discontinuities.
Let Ω ⊂ R n dim (n dim = 1, 2, 3) be an elasto-plastic solid domain, with the reference position vector x ⊂ R n dim . The outer boundary is denoted by Γ ⊂ R n dim −1 , with the outward unit normal vector n * . Deformations of the solid are characterized by the displacement field u(x) and the infinitesimal strain field (x) = ∇ sym u(x), where ∇ sym ( · ) is the symmetric gradient operator.

Elasto-Plasticity Model
In the following, tensorial notation is used. The inner products with single and double contractions are denoted by '·' and ':', respectively, while the dyadic operator is signified by ' '.
For the elasto-plastic model, the constitutive equation is expressed in total form as = e + p , σ = E 0 : e = E 0 : ( − p ) (1) where the second-order strain tensor is decomposed into its elastic and plastic parts, e and p . The second-order stress tensor σ is proportional to the elastic strain tensor e , through the fourth-order elasticity tensor E 0 . All the tensors involved are symmetric. The elastic properties may be orthotropic.
The admissible stress domain is determined by the yield criterion Φ(σ, ζ) = φ(σ) − q(ζ) ≤ 0, defined in terms of the equivalent stress φ(σ) and a stress-like internal variable q(ζ), which determine the shape and size of the domain, respectively. Yield criteria for orthotropic elasto-plasticity are discussed in Section 3.
The plastic strain is defined in rate form, its direction is derived from a plastic potential. In associated plasticity, the plastic potential is equal to the yield surface, so that where . λ ≥ 0 denotes the plastic multiplier; . is the time derivative and the plastic flow tensor Λ = ∂φ /∂σ is normal to the yield surface Φ = 0. Similarly, the evolution of the size of the yield surface is determined by The constitutive equation in rate form follows from Equation (1), where the fourth-order elasto-plasticity tangent tensor E ep is obtained from the Kuhn-Tucker and consistency conditions as where H = ∂q /∂ζ is the hardening or softening modulus. For perfect plasticity, q = q 0 , and H = 0. Note that in associated plasticity, the elasto-plastic tangent tensor is symmetric.

Kinematics of Strong and Regularized Discontinuities
In the early stages of the loading and deformation process of an elasto-plastic solid, standard kinematics applies and both the displacement rate and strain rate fields are continuous. However, in softening and associated perfect plasticity, and even in hardening non-associated plasticity, slip lines (in 2D) or slip surfaces (in 3D) may form. Across these, the deformation can grow unbounded, displacement and/or strain discontinuities may appear and Maxwell's compatibility condition needs to be considered. Figure 1a shows the elasto-plastic solid domain Ω split by a displacement discontinuity S (the slip line or slip surface) into two parts Ω + and Ω − . The orientation of the discontinuity is denoted with the unit normal vector n with direction from Ω − to Ω + . Let L be a characteristic size of the domain.
The constitutive equation in rate form follows from Equation (1), where the fourth-order elasto-plasticity tangent tensor is obtained from the Kuhn-Tucker and consistency conditions as = = − : ⨂ : + : : where = ∂ / ∂ζ is the hardening or softening modulus. For perfect plasticity, = , and = 0. Note that in associated plasticity, the elasto-plastic tangent tensor is symmetric.

Kinematics of Strong and Regularized Discontinuities
In the early stages of the loading and deformation process of an elasto-plastic solid, standard kinematics applies and both the displacement rate and strain rate fields are continuous. However, in softening and associated perfect plasticity, and even in hardening non-associated plasticity, slip lines (in 2D) or slip surfaces (in 3D) may form. Across these, the deformation can grow unbounded, displacement and/or strain discontinuities may appear and Maxwell's compatibility condition needs to be considered. Figure 1a shows the elasto-plastic solid domain Ω split by a displacement discontinuity (the slip line or slip surface) into two parts Ω and Ω . The orientation of the discontinuity is denoted with the unit normal vector with direction from Ω to Ω . Let L be a characteristic size of the domain.  Figure 2a shows the corresponding kinematics: the velocity and strain rate fields are not regular. There is a discontinuity of the displacement rate at of value ; correspondingly, the strain rate at is where denotes the Dirac delta function. Note that this strain rate is unbounded and has a very definite structure determined by the direction of the discontinuity surface, as it allows for unbounded strain rate components at due to the discontinuity of the displacement in the normal direction , but not in those directions tangential to .  where δ S denotes the Dirac delta function. Note that this strain rate is unbounded and has a very definite structure determined by the direction of the discontinuity surface, as it allows for unbounded strain rate components at S due to the discontinuity of the displacement in the normal direction n, but not in those directions tangential to S. For the analysis of strain localization in the continuum setting and also for its numerical verification using FEM, it is convenient to consider a regularized discontinuity, as shown in Figure 1b. Here, subdomains Ω and Ω are separated by a regularized discontinuity band of finite width b, as the distance between surfaces and ; these are weak (strain) discontinuities. The bandwidth b is small compared to the characteristic size of the domain L, so that / << 1. Figure 2b shows the corresponding regularized kinematics. Note that the strain localizes in the regularized band . The deformation rate vector in the strain localization band is defined as the (apparent) jump of displacement rate across the regularized discontinuity band divided by the band width, = / .
Let be a characteristic displacement in domain and the jump be of the same order. Deformations outside the localization band are of the order = / , while inside the band they are of order = / . As / << 1, / << 1 even for a finite, small bandwidth. Denoting by and the strain rates inside and outside of the localization band, respectively, and being the corresponding strain rate jump, Maxwell's compatibility condition [20] is now expressed as Equation (7) is the regularized counterpart of Equation (6). Note that for the band width b → 0, the strain rate in the regularized discontinuity band tends to the strain rate in the strong discontinuity .

Strain Localization and Stress Boundedness
Upon strain localization inside the band, and ongoing deformation, the deformation vector rate in the band, = / , the strain rate jump, ( ⨂ ) , and the plastic strain rate in the band, , will grow much larger than the total strain rate outside the band, , and the corresponding plastic strain rate, , will either vanish (on elastic unloading) or remain small (on plastic loading); this ensures boundedness of the stress rate outside the band, . The terms that will grow upon strain localization, inversely proportional to b, are underlined in the following. For the analysis of strain localization in the continuum setting and also for its numerical verification using FEM, it is convenient to consider a regularized discontinuity, as shown in Figure 1b. Here, subdomains Ω + and Ω − are separated by a regularized discontinuity band B of finite width b, as the distance between surfaces S + and S − ; these are weak (strain) discontinuities. The bandwidth b is small compared to the characteristic size of the domain L, so that b/L << 1. Figure 2b shows the corresponding regularized kinematics. Note that the strain localizes in the regularized band B. The deformation rate vector in the strain localization band . e is defined as the (apparent) jump of displacement rate . w across the regularized discontinuity band divided by the band width, . int the strain rates inside and outside of the localization band, respectively, and being 〚 . 〛 the corresponding strain rate jump, Maxwell's compatibility condition [20] is now expressed as Equation (7) is the regularized counterpart of Equation (6). Note that for the band width b → 0, the strain rate in the regularized discontinuity band B tends to the strain rate in the strong discontinuity S.

Strain Localization and Stress Boundedness
Upon strain localization inside the band, and ongoing deformation, the deformation vector rate in the band, . p int , will grow much larger than the total strain rate outside the band, . ext , and the corresponding plastic strain rate, . p ext , will either vanish (on elastic unloading) or remain small (on plastic loading); this ensures boundedness of the stress rate outside the band, . σ ext . The terms that will grow upon strain localization, inversely proportional to b, are underlined in the following.
From the constitutive relation of the elasto-plastic solids, the stress rates inside and outside of the localization band are given by Note that plastic behavior is considered inside and outside the localization band. The jump of stress rate 〚 . σ〛 is expressed as where the compatibility Condition (7) has been used and the jump of plastic strain rate is Equations (8)-(10) are derived from the constitutive behavior and the compatibility conditions across the weak discontinuities S + and S − ; as strain localization has not been invoked, all the terms involved are bounded.
Inside the localization band, elasto-plastic behavior and satisfaction of the yield criterion ensure that the stress rate needs to remain bounded even if the strain rate is not. Consequently, the jump of the stress rate in Equation (9) may not be null, but it is bounded; therefore, stress rate boundedness requires that The entire jump of the strain rate is due to the plastic strain rate inside the band. This a necessary condition for strain localization to occur. Some Remarks are in order.

Remark 1.
This condition holds for small finite bandwidths b, as in regularized discontinuities and standard FEM simulations. The condition for strong discontinuities follows for the limit case of vanishing bandwidth b → 0.

Remark 2.
This condition does not necessarily occur upon plastic yielding or strain bifurcation. Therefore, a transition stage may be necessary in most situations during which plastic behavior happens without strain localization. Only when the localization condition is fulfilled, might true strain localization happen.

Remark 3.
Only kinematic conditions depending on the plastic flow rule are implied; therefore, the condition may be extended to non-associated plasticity.

Remark 4.
For the same reason, the condition is independent of the elastic properties. Application to rigid-plastic materials can be implied from this independence. This is not the case for classical conditions related to strain bifurcation.

Remark 5.
Stress rate continuity upon strain localization follows from Equation (9) if unloading occurs outside the band, that is, . λ ext = 0. In this case, 〚 . 〛 = 〚 . p 〛 and 〚 . σ〛 = 0. This is usually the case when softening plasticity is considered.

Strain Localization Plastic Flow Vector and Tensor
In the following, the subscript () int will be omitted for the sake of simplicity, as all quantities refer to points inside the localization band. From Equation (10), a plastic flow localization vector, γ, can be defined so that the deformation rate vector . e and the plastic flow tensor Λ are written as . e = . λ γ, Λ = γ n sym (12) where n is the unit vector normal to the discontinuity S. Note that Λ is a second order tensor, while . e, γ and n are vectors. Let m and p be unit vectors on the plane of the discontinuity S such that (n, m, p) is a basis of orthonormal vectors. Then, the plastic flow localization vector, γ, can be equivalently defined so that The components of the plastic flow localization vector γ = (γ n , γ m , γ p ) are determined so that Accordingly, the other components of the strain localization plastic flow tensor are zero: From these equations the orientation of the slip surface may be derived.

Application to Orthotropic Cohesive-Frictional Plastic Materials
In this section, the above results for strain localization in elasto-plastic materials are purposedly applied to orthotropic cohesive-frictional plastic materials. A general form of the considered yield criteria is given that allows closed-form solutions for the orientation of the slip lines in 2D plane strain and plane stress conditions.

Orthotropic Cohesive-Frictional Plasticity
Orthotropic cohesive-frictional yield criteria of the form Φ(σ, ζ) = φ(σ) − q(ζ) ≤ 0 are now considered. Let (1,2,3) be the material orthotropy axes and σ T = [σ 11 , σ 22 , σ 33 , σ 12 , σ 13 , σ 23 ] Voigt's representation of the symmetric second-order stress tensor is used in those axes. Voigt's notation will be used in the following for symmetric second-order tensors. The equivalent stress φ(σ) is expressed as The generalized orthotropic matrix P and Q vectors read where the material parameters F, G, H, F, G, H, L, M, N, I, J and K are given in terms of the material strengths (superscripts c and t denote compression and tension, respectively): Unless otherwise stated: The initial stress threshold is defined as Different well-known quadratic isotropic and orthotropic yield criteria are obtained by appropriately selecting the material parameters: Von Mises criterion: Parabolic Drucker-Prager (DP) criterion: Hill criterion: Hoffman criterion: Tsai-Wu criterion: Remark 6. The effective stress in Equation (18) defines a quadratic yield surface, with a quadratic dependence of the friction-angle on pressure. Alternatively, an effective stress defined as withQ i = √ Q i allowing for a yield surface with straight meridians; the isotropic criterion would the more conventional DP cone.
Orthotropic criteria cannot be represented graphically in the Haigh-Westergaard (HW) stress space because they depend on the six stress components. A partial graphical representation can be obtained by considering them projected into the HW space when the principal stresses act on the material axis, that is, no shear stress appears on the material system. Such representation, generally as an elliptic paraboloid, is offered in Figure 3. All strengths are scaled to 1. Figure 3a shows an orthotropic Hill cylinder, with f 1 / f 2 = f 1 / f 3 = 1.5, tensile and compresive strenth are equal. Figure 3b show the isotropic parabolic Drucker-Prager for compressive to tensile strength ratio κ = f c / f t = 1.5. Figure 3d,c show the orthotropic Hoffman and Tsai-Wu criteria, respectively, for ratios all tensile strengths are taken equal to 1.

Orthotropic Plastic Flow
From the effective stress in Equation (18), the components of the plastic flow tensor plastic flow are obtained: So that The identity tr = + + = ( + + ) holds.

Strain Localization Angle
In this section, the orientation of the slip lines is analytically obtained for orthotropic and pressure-dependent plastic solids subjected to plane strain and plane stress conditions. The strain localization angle is measured counter-clockwise ∈ [− , ] as the angle between the vector normal to the discontinuity and the material axis 1; see Figure  4.

Orthotropic Plastic Flow
From the effective stress in Equation (18), the components of the plastic flow tensor plastic flow are obtained: So that The identity trΛ = Λ 11

Strain Localization Angle
In this section, the orientation of the slip lines is analytically obtained for orthotropic and pressure-dependent plastic solids subjected to plane strain and plane stress conditions. The strain localization angle is measured counter-clockwise θ cr ∈ [− π 2 , π 2 ] as the angle between the vector n normal to the discontinuity and the material axis 1; see Figure 4. Let ( , , ) be the basis formed by the orthonormal vectors normal and tangential to the discontinuity , such that vectors and are respectively normal and tangential to the trace of S in the reference plane xy and vector points in the out-of-plane z direction, as shown in Figure 4.
The strain localization Equation (16) requires the flow tensor in Equation (31) to be written in this system. Let be the angle between the material system (1, 2, 3) and the ( , , ) system. Then The strain localization angle is obtained from the kinematic constraints in Equation (16), that is, equating these components to zero. Solving Λ ( ) = 0 for tan : As can be seen, the strain localization angle depends on the stress state upon strain localization. The condition ( ) = = 0 (38) needs to be imposed in plane stress and strain conditions.  Let (n, m, p) be the basis formed by the orthonormal vectors normal and tangential to the discontinuity S, such that vectors n and m are respectively normal and tangential to the trace of S in the reference plane xy and vector p points in the out-of-plane z direction, as shown in Figure 4.
The strain localization Equation (16) requires the flow tensor in Equation (31) to be written in this system. Let θ cr be the angle between the material system (1, 2, 3) and the (n, m, p) system. Then The strain localization angle θ cr is obtained from the kinematic constraints in Equation (16), that is, equating these components to zero. Solving Λ mm (θ cr ) = 0 for tan θ cr : As can be seen, the strain localization angle θ cr depends on the stress state upon strain localization. The condition Λ pp (θ cr ) = Λ 33 = 0 (38) needs to be imposed in plane stress and strain conditions.

Remark 7.
For the case of Λ 12 = 0, where the no shear stress acts on the material axes, Equation (39) simplifies to Remark 8. The kinematic constraints produce two alternative strain localization angles, see also Figure 3. The in-between angles that follow from Equation (40) are Remark 9. In purely isochoric models (von Mises, Hill), Λ 11 = −Λ 22 , and tan θ cr Remark 10. The angle of the slip lines (counter-clockwise from 1-axis) is θ slip = π 2 − θ cr : Remark 11. The above expressions are obtained for the stress expressed in the material system. These are obtained from the stresses in the global (x,y,z) system by standard transformation. For instance, in plane strain conditions where α is the tilt angle between the global axis x and the material local axis 1 measured counter-clockwise.

Plane Stress
In plane stress, σ 33 = σ pp = 0. In this case, the non-zero plastic flow components Λ ij in Equation (40) are These components can be substituted in Equation (39).

Plane Strain
In this case, the non-zero plastic flow components Λ ij in Equation (40) are considered with Equation (38).

Geometrical Interpretation of the Strain Localization Angle in the Stress Space
In the following, a geometrical interpretation of the strain localization angles obtained analytically is offered. As explained, Figure 3 gives a partial graphical representation of the orthotropic yield criteria projected into the HW space when the principal stresses act on the material axis, that is, no shear stress appears on the material system.
In Figure 5, a cross section of those paraboloids by a horizontal plane is given. For plane stress, the plane σ 33 = 0 is used; for plane strain, the plane σ 33 = [2( Gσ 11 + Hσ 22 ) − K]/ 2(G + H), from Equation (46), is used. The isotropic Drucker-Prager and the orthotropic Hoffmann and Tsai-Wu criteria are depicted for ratios κ = f c / f t = 1.5 and 3.0 and f c 2 / f t 2 = f c 3 / f t 3 = 1; all tensile strengths are taken equal to 1. For plane stress, the intersected quadratic curves are ellipses; for plane strain, they are parabolas. They are more stretched for higher ratios κ. tained analytically is offered. As explained, Figure 3 gives a partial graphical representation of the orthotropic yield criteria projected into the HW space when the principal stresses act on the material axis, that is, no shear stress appears on the material system.
In Figure 5, a cross section of those paraboloids by a horizontal plane is given. For plane stress, the plane = 0 is used; for plane strain, the plane = 2 + − / 2( + ), from Equation (46), is used. The isotropic Drucker-Prager and the orthotropic Hoffmann and Tsai-Wu criteria are depicted for ratios κ = / = 1. 5 3.0 and / = / = 1; all tensile strengths are taken equal to 1. For plane stress, the intersected quadratic curves are ellipses; for plane strain, they are parabolas. They are more stretched for higher ratios κ. In Figure 5, the projection of the plastic flow vector, normal to the yield surface, for uniaxial tension and compression in the 2-direction, is also plotted. See the next section for the analytical values. In Figure 5, the projection of the plastic flow vector, normal to the yield surface, for uniaxial tension and compression in the 2-direction, is also plotted. See the next section for the analytical values.

Remark 12.
The angle θ between this projected flow vector and the 2-axis is related to the strain localization angle θ cr , because

Uniaxial Tension and Compression: Analytical Strain Localization Angles
In the following, the analytical values of the strain localization angle are obtained for the uniaxial tension and compression cases illustrated in Figure 5. Material strengths are those indicated in the previous section; with those, the coefficients for matrix P and vector Q are computed for the three different criteria (Drucker-Prager, Hoffmann and Tsai-Wu) and listed in Table 1.
Therefore, the kinematical condition Λ 33 = 0 (38) needs not be considered, and no extra constraint imposes on the stress state upon strain localization. Therefore, once the initial yield surface, Φ(σ, ζ) = 0, is reached, strain localization occurs at the same instant, with the orientation determined from the corresponding flow tensor.
As Λ 12 = 0, The obtained values for θ slip are given in Table 2. Results corresponding to uniaxial compression are also given in the Table. Note that the localization angles under tension and compression are very different for the various yield criteria, as depicted graphically in Figure 5.
This happens if the compressive strength is sufficiently larger than the tensile strength; for instance, it happens for the ratio κ = 2 for the Drucker-Prager criterion. For larger ratios, there is no real value for the localization angles. In compression, for σ 22 = σ = −1, this happens reciprocally, that is, if the compressive strength is sufficiently smaller than the tensile strength.

Plane Strain
In the plane strain case, the kinematical condition Λ 33 = 0 (38) needs to be enforced. From this, A stress σ 22 = σ is found so that the point [0, σ, σ 33 ] is on the corresponding yield surface, see Figure 5a. Then, The obtained values for θ slip are given in Table 3. Results corresponding to uniaxial compression are also given. Note that the angles under tension and compression are distinct, and they are also different for the various yield criteria, as shown in Figure 5.

Numerical Verification
In this section, FEM analyses are performed to numerically verify the analytical results obtained in Section 3 and derived from the strain localization analysis in Section 2.
It is emphasized that the numerical verification is totally independent of the analytical results. That is, the numerical analyses follow the standard procedure for solving the nonlinear mechanical problem and plastic behavior appears and evolves into the formation of slip lines; and the analytical results are not by any means used.
In previous works [45], it has been demonstrated that the strain localization angle is independent of the elastic properties. Therefore, the argument is not pursued here. Similarly, the localization angle does not depend on the softening behavior [47], so perfect plasticity is assumed in the following.
Although pressure-dependent plasticity models are to be investigated, they are compared to the isochoric von Mises model. To avoid volumetric locking in nearly incompressible situations, B-bar finite elements [49] are used in 2D and 3D.

B-Bar Finite Element
The B-bar element is a particular implementation of the mixed displacement/pressure Q1P0 element in which the constant pressure has been eliminated at element level at the expense of renouncing the incompressible limit. This is accomplished by evaluating the constant mean stress in terms of the mean volumetric strain, the latter computed from the nodal displacements.
The standard discrete strain-displacement B matrix, computed at each integration point from the Cartesian derivatives of the nodal shape functions, is split into its volumetric and deviatoric parts A mean volumetric sub-matrix B vol is computed as where n g is the number of integration points in the element. The B-bar discrete strain-displacement matrix is obtained as The B-bar element has some zero-energy modes that may show as spurious hourglassing in some instances. This may be avoided by using For τ = 1, then B stab = B vol + B dev = B is identical to the B-bar formulation. For τ = 0, then B stab = B vol + B dev is identical to the standard formulation.

Uniaxial Tension and Compression: Numerical Verification
In this section, the above B-bar finite element is used to perform benchmark verifications in strain localization analysis. The benchmark example is a strip loaded in uniaxial tension and compression via imposed vertical displacements at the top and bottom ends; the horizontal movement is not restrained. As shown in Figure 6, the strip has dimensions 10 m × 20 m (width × height). A sharp horizontal slit (2 m) is inserted in the center of strip to introduce the perturbation necessary to trigger strain localization. where n is the number of integration points in the element. The B-bar discrete strain-displacement matrix is obtained as

= +
The B-bar element has some zero-energy modes that may show as spuriou glassing in some instances. This may be avoided by using For τ = 1, then = + = is identical to the B-bar formulat τ = 0, then = + is identical to the standard formulation.

Uniaxial Tension and Compression: Numerical Verification
In this section, the above B-bar finite element is used to perform benchmark tions in strain localization analysis. The benchmark example is a strip loaded in tension and compression via imposed vertical displacements at the top and botto the horizontal movement is not restrained. As shown in Figure 6, the strip has dim 10 m × 20 m (width × height). A sharp horizontal slit (2 m) is inserted in the cente to introduce the perturbation necessary to trigger strain localization. In this problem the stress field is known a priori. Plane strain and plane stres tions are investigated. In both cases, the far field stress state corresponds exactly assumed for the analytical results in Section 3: The sharp horizontal slit causes a stress concentration that triggers the onset o behavior and strain localization; subsequently, straight slip lines stem from th In this problem the stress field is known a priori. Plane strain and plane stress conditions are investigated. In both cases, the far field stress state corresponds exactly to those assumed for the analytical results in Section 3: The sharp horizontal slit causes a stress concentration that triggers the onset of plastic behavior and strain localization; subsequently, straight slip lines stem from these and cross the strip at well-defined slopes that must follow the angles analytically predicted in Sections 2 and 3. The numerical results obtained in the FE analysis are used to validate the strain localization analysis in Section 2 and the analytical results in Section 3 that follow from it.
The following material properties are used: Young's modulus E = 1.0 × 10 7 MPa, Poisson's ratio ν = 0.2. Several orthotropic elasto-plastic criteria are compared; the different plastic yield strengths along the material axes are detailed for each case. Perfect plasticity is assumed.
Structured meshes of regular quadrilateral are employed. Square elements (0.05 m × 0.05 m) are arranged 200 horizontally and 400 vertically, with a total of 80,000 elements used for plane strain 2D simulations. Plane stress cases are simulated in 3D, with as many hexahedral elements arranged in a mesh 1 element thick. In all cases, 500 time steps are performed to complete the analyses. The constitutive laws and finite elements used have been implemented in the COMET finite element program, developed by the authors at the International Center for Numerical Methods in Engineering (CIMNE). Preand post-processing are done with GiD, also developed at CIMNE.

Isotropic Incompressible and Cohesive-Frictional Models
In this subsection, strain localization is first investigated for isotropic incompressible and pressure sensitive models.
Isotropic The following material properties are used: Young's modulus = 1.0 × 10 Poisson's ratio ν = 0.2. Several orthotropic elasto-plastic criteria are compared; the d ent plastic yield strengths along the material axes are detailed for each case. Perfect ticity is assumed.
Structured meshes of regular quadrilateral are employed. Square elements (0.0 0.05 m) are arranged 200 horizontally and 400 vertically, with a total of 80,000 elem used for plane strain 2D simulations. Plane stress cases are simulated in 3D, with as hexahedral elements arranged in a mesh 1 element thick. In all cases, 500 time step performed to complete the analyses. The constitutive laws and finite elements used been implemented in the COMET finite element program, developed by the auth the International Center for Numerical Methods in Engineering (CIMNE). Pre-and processing are done with GiD, also developed at CIMNE.

Isotropic Incompressible and Cohesive-Frictional Models
In this subsection, strain localization is first investigated for isotropic incompre and pressure sensitive models.
Isotropic von Mises J2 plasticity with yield strength = 1.0 × 10 MPa is us reference case. Insensitive to pressure, under plane strain, tensile and compression show the same localization angles (±45°), while under plane stress, the localization a are ±35.26° from the horizontal axis, measured in a counter-clockwise manner.
Isotropic Parabolic Drucker-Prager models are also considered. A tensile str = 1.0 × 10 MPa and different compressive strengths according to the ratio / are used, κ = 1.25, 1.50 for tension, κ = 2, 3 for compression ; the isotropic strength is      In the figures in this and the following sections, contour fills of the equivalent plastic strain are depicted to show the orientation of the slip lines and the corresponding failure mechanisms. The resolution of the mesh and the color pattern are selected so that these can be easily perceived. The red to blue color range indicates the largest to smallest magnitude of the equivalent plastic strain. For the numerical results, the numerical Lode angle is measured at the point in the slip line located 1 m to the right from the right end of the slit. The angle of the slide slip line is measured counter-clockwise from the x-axis.
For all cases, the numerical results are coincident with the analytical results. Correct angles of the slip lines are depicted in Figures 7-10. Also, the coincidence between analytical and numerical results is shown in Tables 4-7. The strain localization angle decreases with increasing ratios κ in the tensile tests (Figures 7 and 9), while the strain localization angle increases with κ in the compressive tests (Figures 8 and 10). The coincidence of the analytical and numerical Lode angles in the plane strain cases indicates that the kinematical constraint imposed by the Λ 33 = Λ zz = 0 condition is verified.

Isotropic and Orthotropic Cohesive-Frictional Models
In this subsection, the formation of slip lines is now investigated for orthotropic Hoffman and Tsai-Wu pressure-sensitive models and compared to the isotropic counterpart. The orthotropy material axes (1,2,3) are coincident with the global axes (x,y,z); relative tilting is investigated in Appendix B.
For the comparison, a ratio of compressive to tensile strengths κ = 1.5 is taken for the tension tests and κ = 3.0 for the compression tests. For the orthotropic models, all the yield strengths are taken as f = 1.0 × 10 4 MPa, except the compressive f c x , which is taken to the κ ratio; shear strength f xy = f c x f c y / 3. Some of the corresponding analytical results are given in Section 3.5. In the following, the analytical and numerical results are presented for comparison in Figures 11-14 and Tables 8-11.       Table 11. Analytical and numerical Lode and strain localization angles for frictional-cohesive models under plane strain compression, κ = 3.0.  As previously, for all cases, the numerical and analytical results are coincident. Note that Hoffman and Tsai-Wu models produce different outcomes for the same material properties, as they use different F, G and H parameters. Lode angles in plane stress are 0 • under tensile loading and 60 • under compressive loading; they vary in plane strain.

Prandtl's Punch Test
The second example is Prandtl's punch test, a 2D plane problem in which a flat rigid die punches into an elasto-plastic semi-infinite medium. Classical solutions to this problem for rigid-plastic materials are well-known.
As shown in Figure 15, the computational domain of the elasto-plastic medium is 10 m × 3 m (width × height). Boundary conditions consist of a fixed bottom edge, left and right edges horizontally restrained. Punching is applied by imposing an increasing vertical displacement at the base of the rigid die; the horizontal movement is restrained at this base. As shown in Figure 15, the computational domain of the elasto-plastic medium is 10 m × 3 m (width × height). Boundary conditions consist of a fixed bottom edge, left and right edges horizontally restrained. Punching is applied by imposing an increasing vertical displacement at the base of the rigid die; the horizontal movement is restrained at this base. The mechanics of the failure are as follows. Plastic yielding starts at the singular points at the extreme ends of the rigid die. From here, two slip lines dig into the supporting elasto-plastic medium at diverging angles. Further loading causes the formation of a collapse mechanism in which the triangular wedge of material immediately under the punch moves vertically, causing the outward lateral displacement of adjoining material and the upwards displacement of the material located in the triangular wedges close to the surface and next to the flat punch. Figure 16 shows the numerically obtained failure mechanisms for the four different cases studied, depending on the plastic criterion used in each one: (a) isotropic pressureindependent von Mises; (b) isotropic pressure-dependent Parabolic Drucker-Prager, κ = 3.0; (c) orthotropic Hoffman, κ = 3.0 in the horizontal direction; and (d) orthotropic Tsai-Wu, κ = 3.0, in the horizontal direction. Associated perfect plasticity is used, so that the The mechanics of the failure are as follows. Plastic yielding starts at the singular points at the extreme ends of the rigid die. From here, two slip lines dig into the supporting elasto-plastic medium at diverging angles. Further loading causes the formation of a collapse mechanism in which the triangular wedge of material immediately under the punch moves vertically, causing the outward lateral displacement of adjoining material and the upwards displacement of the material located in the triangular wedges close to the surface and next to the flat punch. Figure 16 shows the numerically obtained failure mechanisms for the four different cases studied, depending on the plastic criterion used in each one: (a) isotropic pressure-independent von Mises; (b) isotropic pressure-dependent Parabolic Drucker-Prager, κ = 3.0; (c) orthotropic Hoffman, κ = 3.0 in the horizontal direction; and (d) orthotropic Tsai-Wu, κ = 3.0, in the horizontal direction. Associated perfect plasticity is used, so that the plastic potential coincides with the described yield criteria. The plots are zoomed in the region of interest, with identical magnification. cases, 1000 time steps are performed to complete the analyses.
The mechanics of the failure are as follows. Plastic yielding starts at the singular points at the extreme ends of the rigid die. From here, two slip lines dig into the supporting elasto-plastic medium at diverging angles. Further loading causes the formation of a collapse mechanism in which the triangular wedge of material immediately under the punch moves vertically, causing the outward lateral displacement of adjoining material and the upwards displacement of the material located in the triangular wedges close to the surface and next to the flat punch. Figure 16 shows the numerically obtained failure mechanisms for the four different cases studied, depending on the plastic criterion used in each one: (a) isotropic pressureindependent von Mises; (b) isotropic pressure-dependent Parabolic Drucker-Prager, κ = 3.0; (c) orthotropic Hoffman, κ = 3.0 in the horizontal direction; and (d) orthotropic Tsai-Wu, κ = 3.0, in the horizontal direction. Associated perfect plasticity is used, so that the plastic potential coincides with the described yield criteria. The plots are zoomed in the region of interest, with identical magnification. As can be observed, similar but notably different failure mechanisms form, depending on the plastic potential that applies. Although the process of the formation of the slip lines and the failure mechanism is analogous in all cases, the observed discrepancies in As can be observed, similar but notably different failure mechanisms form, depending on the plastic potential that applies. Although the process of the formation of the slip lines and the failure mechanism is analogous in all cases, the observed discrepancies in the slopes of the slip lines, and the corresponding amounts of mobilized material, depend on the plastic material properties.
Contrariwise to the case studied in the previous section, here the stress field is known a priori. Furthermore, substantial stress redistribution happens in the transition between the initial elastic stage and the final elasto-plastic state in which the failure mechanism is completely formed and yielding. This can be observed in Figure 17, where the distribution of the principal stresses in the elastic (initial) and plastic (stationary) states in the region below the punch are compared for the Drucker-Prager case (b). It can be seen that the stress state in the elastic range consists mainly of vertical stress σ yy and the corresponding out of plane σ zz (not shown in the figure), due to Poisson's ratio ν = 0.2. In contrast, in the stationary plastic stress state, in-plane horizontal σ xx have noticeably developed. Contrariwise to the case studied in the previous section, here the stress field is known a priori. Furthermore, substantial stress redistribution happens in the transition between the initial elastic stage and the final elasto-plastic state in which the failure mechanism is completely formed and yielding. This can be observed in Figure 17, where the distribution of the principal stresses in the elastic (initial) and plastic (stationary) states in the region below the punch are compared for the Drucker-Prager case (b). It can be seen that the stress state in the elastic range consists mainly of vertical stress and the corresponding out of plane (not shown in the figure), due to Poisson's ratio ν = 0.2. In contrast, in the stationary plastic stress state, in-plane horizontal have noticeably developed. The extension and nature of this stress redistribution is further investigated in Figure  18, where the evolution of the normal stress components , , against the vertical displacement of the die is plotted for the four cases. The stresses are sampled at the corresponding crossing point of the slip lines, in the symmetry axis below the center of the punch. Figure 19 further summarizes the comparison of stress evolution by plotting the evolution of the stress Invariant and the Lode angle . The extension and nature of this stress redistribution is further investigated in Figure 18, where the evolution of the normal stress components σ xx , σ yy , σ zz against the vertical displacement of the die is plotted for the four cases. The stresses are sampled at the corresponding crossing point of the slip lines, in the symmetry axis below the center of the punch. Figure 19 further summarizes the comparison of stress evolution by plotting the evolution of the stress Invariant I 1 and the Lode angle ϑ. The extension and nature of this stress redistribution is further investigated in Figure  18, where the evolution of the normal stress components , , against the vertical displacement of the die is plotted for the four cases. The stresses are sampled at the corresponding crossing point of the slip lines, in the symmetry axis below the center of the punch. Figure 19 further summarizes the comparison of stress evolution by plotting the evolution of the stress Invariant and the Lode angle . Several remarks are in order: (a) the extension of the transition phase largely differs from one case to the other; it is shorter for von Mises and longer for Drucker-Prager; (b) due to increasing vertical loading, the out of plane develops due to the plane strain constraint; (c) concurrently, the in-plane horizontal also develops, very much depending on the yield criterion; this later development shows the frictional and/or orthotropic nature of the plastic behavior.
It can be verified that this stress redistribution during the formation of the slip lines occurs precisely as dictated by the strain localization condition. This is done in Table 12 by comparing the value , measured directly from Figure 16, with the value , obtained by applying the analytical condition (Sections 3.2.2 and 3.5.2) to the numerically obtained values for the stresses. The correspondence between both value is remarkable. Several remarks are in order: (a) the extension of the transition phase largely differs from one case to the other; it is shorter for von Mises and longer for Drucker-Prager; (b) due to increasing vertical loading, the out of plane develops due to the plane strain constraint; (c) concurrently, the in-plane horizontal also develops, very much depending on the yield criterion; this later development shows the frictional and/or orthotropic nature of the plastic behavior.
It can be verified that this stress redistribution during the formation of the slip lines occurs precisely as dictated by the strain localization condition. This is done in Table 12 by comparing the value , measured directly from Figure 16, with the value , obtained by applying the analytical condition (Sections 3.2.2 and 3.5.2) to the numerically obtained values for the stresses. The correspondence between both value is remarkable. Several remarks are in order: (a) the extension of the transition phase largely differs from one case to the other; it is shorter for von Mises and longer for Drucker-Prager; (b) due to increasing vertical loading, the out of plane σ zz develops due to the plane strain constraint; (c) concurrently, the in-plane horizontal σ xx also develops, very much depending on the yield criterion; this later development shows the frictional and/or orthotropic nature of the plastic behavior.
It can be verified that this stress redistribution during the formation of the slip lines occurs precisely as dictated by the strain localization condition. This is done in Table 12 by comparing the value θ slip num , measured directly from Figure 16, with the value θ slip ana , obtained by applying the analytical condition (Sections 3.3.2 and 3.5.2) to the numerically obtained values for the stresses. The correspondence between both value is remarkable.

Conclusions
In this work, the strain localization analysis of cohesive-frictional elasto-plastic materials is addressed, which applies to both strong and regularized slip lines and surfaces. Maxwell kinematics, stress boundedness and plastic consistency are invoked to derive the necessary strain localization conditions. Contrariwise to the usually studied conditions for strain bifurcation, these proffer requirements that do not depend on the elastic properties of the medium, but only on the plastic flow provided by the adopted plastic potential.
Expanding on previous works, application of the above localization conditions to isotropic and orthotropic cohesive-frictional plastic models derives analytical solutions for the strain localization angle and the slopes of the ensuing slip lines. The distinct effects of compressive and tensile loading are also evaluated.
The analytical results are validated independently by 2D plane stress and plane strain FE simulations using the B-bar element; namely, a strip under vertical tension and compression tests and Prandtl's punch problem are investigated. In the first problem, the far field stress state is known, and the analytical results can be verified directly from the numerical simulations. In the second problem, once the failure mechanism and the corresponding stress field are computationally evaluated, these are shown to conform precisely with those anticipated by the strain localization condition.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.
In this Appendix, the effect of the tilting of the orthotropy material axes with respect the global axes and the orientation of the loading is demonstrated. The tilt angle α is measured counter-clockwise between the global x and the material 1 axes. The rotation transformation matrix was introduced in Remark 3.6. Figure A1 and Table A1 show the results for the strip under vertical plane strain tension (Parabolic Drucker-Prager, κ = 1.5) and different tilting; as the model is isotropic, the results are insensitive to the rotation of the material axes.   Figure A2 and Table A2 show corresponding results for the orthotropic Hoffman model, κ = 1.5. Here, the effect of the tilting of the material axes is evident. Analytical and numerical results coincide for all slip lines.   Figure A2 and Table A2 show corresponding results for the orthotropic Hoffman model, κ = 1.5. Here, the effect of the tilting of the material axes is evident. Analytical and numerical results coincide for all slip lines.  Finally, Figure A3 shows the effect of the tilting of the material axes in Prandtl's punch test (Plane Strain, Hoffman model, κ = 3.0).    Figure A3 shows the effect of the tilting of the material axes in Prandtl's punch test (Plane Strain, Hoffman model, κ = 3.0).  Finally, Figure A3 shows the effect of the tilting of the material axes in Prandtl's punch test (Plane Strain, Hoffman model, κ = 3.0).