The paper presents a new tension failure criterion which generalizes the so-called Galileo–Rankine formulation. The criterion can be used as a component of the so-called perfectly no-tension model for masonry and cements as well as for establishing a tension cut-off in complex constitutive models for soils, granular materials and powders. The criterion is described by means of a very concise equation based on the third invariant of the stress tensor, approximating the boundaries of the compressive octant of the principal stress space. This sheds new light on the physical significance of the third invariant of the stress tensor. The new criterion has been validated against two known analytical solutions for no-tension materials and also effectively applied for solving two geotechnical and structural engineering problems. The proposed formulation allows for an efficient implementation in finite-element programmes, removing some of the numerical difficulties associated with the Galileo–Rankine criterion.
Accounting for the limited ability of materials to resist tensile stresses introduces considerable difficulties when performing numerical analyses of boundary value problems. This is predominantly relevant for bonded and unbonded particulate aggregates (e.g. powders, snow, soils, weak rocks) which are characterized by a variable capability to resist tensile stresses, building materials for which the tensile strength is disregarded (e.g. masonries and cement) and materials for which the tensile strength is considered in calculations (e.g. rocks and ice).
Natural and artificial unbonded aggregates of particles (e.g. powders, seeds, clean sands, clean gravels, etc.) cannot resist tensile stresses. In strain hardening/softening plasticity, this can be accounted for by adopting a yield surface that stems from the origin of the principal stress space and projects into the compressive octant. However, yield surfaces usually cross the coordinate planes which delimit the compressive octant, thus including negative stress states within the admissible elastic domain. This is clearly shown in figures 1 and 2 for two classical yield surfaces, namely that of the modified Cam–clay model  and that proposed in . Figure 1 presents a section of the two surfaces with the axis-symmetric stress plane σ2=σ3 (usually indicated as triaxial plane), plotted in terms of the and p=I1/3, where p is the mean pressure, I1 is the first invariant of the stress tensor and J2D is the second invariant of the deviatoric stress tensor. The σ1 and the σ2=σ3 principal stress axes are also indicated for clarity with dotted lines. It clearly appears that if the tangent to the yield curve at the origin of the triaxial plane is steeper than , the curve is also defined for negative values of σ2=σ3. This implies that smooth yield surfaces, which are characterized by a tangent , necessarily include negative stress states within the elastic domain. As smoothness is usually considered a desirable feature (e.g. [3,4]), most yield surfaces are affected by this problem. This is also clearly shown in figure 2, where a three-dimensional representation of the two surfaces is given. Constitutive models with such yield surfaces require tension cut-off to confine admissible stresses within the compressive octant. It should be noted, however, that a consistent elasto-plastic cut-off should be defined, so that plastic deformations develop when a stress predictor violates the no-tension condition.
Bonded particulate materials, on the other hand, exhibit tensile strength. These include, among others, compressed ceramic and metal powders and snow. Soils too can resist tensile stresses. When no ‘bonds’ are present at the grain contacts, as occurs in freshly sedimented deposits and clean sands and gravels, the soil cannot withstand any tension. However, some natural soils exhibit a variable, but usually limited, amount of tensile strength as a result of bonding forces between grains which are associated to the presence of ‘cement’ at particles contacts. This occurs in diagenized and in residual soils. In the former case, sedimentary soils undergo physical, chemical and/or biological processes which bind together the particles, while in the latter case, a parent rock is weathered into a weaker and porous material. On this regard, diagenesis and weathering can be viewed as opposite phenomena that can progress in nature to various degrees and therefore the amount of ‘bonding’ in diagenized and residual soils can vary from that of an unbounded particle aggregate to that of a rock (e.g. ). A well-documented example of diagenesis inducing tensile strength in a bioclastic sediment can be found in . As tensile strength is beneficial for the stability of many geotechnical structures, both natural (e.g. slopes) and man-made (e.g. tunnel excavation), it can also be artificially generated on soils by means of cement injection, freezing of pore water and, more recently, by stimulating bacteria activity. Constitutive models for bonded particulate materials are characterized by a yield surface which starts at a negative mean pressure p (e.g. ) and provides an isotropic tensile strength. Also in this case, however, a cut-off is necessary so that tensile stresses larger than the uni-axial tensile strength cannot be attained.
In Structural Mechanics, a simplified approach based on perfect elasto-plasticity is very often adopted to account for the inability of materials to resist tensile stresses. As failure in masonry is rarely caused by compressive forces (e.g. ), the material behaviour can be reproduced by means of the so-called perfectly no-tension constitutive model. This assumes that the stress tensor is always positive semidefinite (in this manuscript the Soil Mechanics sign convention is adopted, therefore compressive stresses are considered positive) and that the stress–strain behaviour is hyper-elastic. Failure occurs if one of the principal stresses becomes nil, and therefore the yield and the plastic potential surfaces of the model are given by the Galileo–Rankine failure criterion (e.g. [9,10]), the tensile and compressive strengths being set as nil and as infinite, respectively (figure 3). With these parameters, the failure criterion plots as the three coordinate planes which delimit the compressive octant of the principal stress space. When one of the principal stresses becomes nil, plastic deformations develop in the same direction, and therefore crack development is simulated by means of a continuous strain and displacement field, rather than a discontinuous one. Despite its simplicity, it has been shown that the model is capable of predicting experimentally observed crack patterns (e.g. ). However, its numerical implementation, even in plane stress conditions, is not trivial and convergence is slow and difficult. These problems are associated with the activation of zero energy modes, as described for example in . Extensive literature can be found on this subject showing how various methods, often very complex, have been suggested to try to eliminate the problems associated with the implementation of the no-tension model. However, numerical difficulties are also related to the non-smoothness of the Galileo–Rankine surface and, as it will shown later on, to the lack of full convexity when formulated in terms of invariants.
All applications discussed so far indicate that the coordinate planes of the principal stress space play a fundamental role as a tension cut-off. However, the Galileo–Rankine failure criterion provides a mathematical description of each of those planes separately and, in addition, severe numerical difficulties arise when implemented in a numerical program (see §5). In this paper, a very simple equation is presented for approximating the three coordinate planes which delimit the compressive octant of the principal stress space. The advantage over the original criterion being that a smooth surface is now provided, hence eliminating the discontinuities between each pair of coordinate planes. The formulation is suitable to be used both for defining the yield and plastic potential surfaces of a perfectly no-tension model and within more complex models such as those adopted in Soil Mechanics to prevent attaining negative stresses. This formulation also sheds light on the physical interpretation of the third invariant of the characteristic equation of the stress tensor.
2. Proposed formulation for a no-tension surface
The coordinate planes which delimit the compressive octant of the principal stress space represent the natural tension cut-off surface to be adopted within a constitutive model. A mathematical representation of such a surface is provided by the Galileo–Rankine failure criterion (e.g. [9,10]) when a specific set of parameters is selected. The criterion is suitable for brittle materials and states that failure is attained when the maximum or minimum principal stresses reach the limit values σft and σfc in tension and compression, respectively. These limit values can be experimentally determined by means of mono-axial tension and compression tests. The Galileo–Rankine criterion is therefore formulated as 2.1 where σ1, σ2 and σ3 are the principal stresses. It should be noted that the Galileo–Rankine criterion plots for a general set of parameters as a cube, its faces parallel to the coordinate planes. However, for a material which can withstand zero tensile stress and infinite compressive stresses, the criterion simply describes the required coordinate planes.
The Galielo–Rankine criterion is frequently adopted for the definition of the no-tension model used in structural mechanics and, more seldom, as a tension cut-off within more complex elasto-plastic constitutive models (e.g. ). However, the mathematical formulation of the Galileo–Rankine criterion is not ideal for being used within a constitutive model, as the complete tension cut-off surface is described by three different equations, and, more importantly, discontinuities occur at the intersection between each couple of coordinate planes. The implementation of no-tension models based on the GR failure criterion, has in fact proved to be very problematic, and there is general agreement in the scientific community that this is related to the activation of zero energy modes (e.g. ).
An alternative formulation that approximates the boundary of the positive octant of the stress space is here proposed 2.2 where σ is the principal stress tensor and k is a small constant. It will be noted that as the third invariant of the characteristic equation of the stress tensor is I3=σ1σ2σ3, the surface can be rewritten as 2.3 which is a more convenient formulation as it does not require the evaluation of principal stresses. This expression can be considered as a smooth approximation of the Galileo–Rankine criterion for brittle no-tension materials.
Figure 4a shows that equations (2.2) or (2.3) are indeed capable of approximating the coordinate planes which delimit the positive octant of the principal stress space. The parameter k controls how closely the surface coincides with the coordinate planes and while the surface in the figure has been drawn with a relatively large value of k for the sake of clarity, a very good delimitation of the octant can be achieved by selecting a small value. It should be noted that k1/3 measures the distance of the tip of the surface from the origin of the stress space along the mean pressure axis and that the smaller its value the closer the surface approximates the boundaries of the compressive octant. However, selecting too small a value may affect the rate of convergence of numerical analyses without improving the accuracy and its usefulness should be evaluated case by case. For instance, in the finite-element (FE) analyses presented in the following sections, very good results were achieved with k1/3 ranging between 10−3p0 to 10−6p0, where p0 is the atmospheric pressure.
Equation (2.3) is very attractive for its conciseness, however, as shown in figure 4b, it plots in four of the eight octants of the principal stress space and therefore there are three more branches (hereafter indicated as negative branches) other than that in the positive octant. This is a frequently occurring problem which affects many classical criteria and a source of numerical difficulties. As shown in , one such example is the well-known Matsuoka–Nakai  failure criterion, for which an exact mathematical reformulation has been recently proposed in , that eliminates the undesired branches and plots as a convex surface. Similarly, in the next section a mathematical approach is proposed which makes equation (2.3) convex by removing all branches other than that of the compressive octant.
3. Convexification of multi-branch surfaces
The equation proposed in the previous section for approximating the boundaries of the positive octant of the stress space, plots as a multi-branch surface which spans also in negative octants (figure 4). As shown in , there are numerous problems associated with the implementation within a FE programme of constitutive models which adopt a multi-branch surface for reproducing the yield and/or plastic potential surfaces. The most obvious inconvenience is that the branches, which plot in octants other than the positive, define false elastic domains. If, for instance, the constitutive law is numerically integrated by means of a full backward Euler implicit scheme stress states resulting from an elastic predictor can happen to lie within such a false elastic domain. Clearly, this eventuality is more likely to occur for Gauss points with an initial stress state close to the origin. However, more complex scenarios can be envisaged, if one considers that the return algorithm relies on the normal to the plastic potential to bring the stress state back to a consistent yield surface. In this section, a general approach is proposed which effectively eliminates the three undesired branches of equation (2.3) and is likely to be effective also for other surfaces with similar characteristics.
Failure and yield surfaces bound admissible and elastic stress states and are generally expressed as 3.1 where σ is the stress tensor and pk is a vector of parameters. It should be noted that in the case of strain hardening/softening plasticity some of the parameters of pk are ‘hidden variables’ which are defined by evolution laws and describe how the surface expands or contracts.
A stress state σ* will be inside, outside or on the yield surface, depending on whether f(σ*,pk) evaluates to <0, >0 or =0, respectively. Moreover, an important requirement is that for stress states outside the surface, the larger the distance the larger the value of f>0. Clearly, this is not the case for expressions which plot as a multi-branch surface and, indeed f does not grow monotonically with the distance from the surface (e.g. ).
However, an inspection of figure 4 suggests that the positive branch of the surface given by equation (2.3) can span the whole stress space if it is translated along the trisectrix of the positive octant. The amount of translation, T, can be measured along the trisectrix, as the difference in mean pressure between the intersections with that axis of the original and of the translated positive branches, as indicated in figure 5. For any given stress state σ*, a well-defined translation exists, T(σ*,pk), which moves the positive branch until it includes that stress state. T(σ*,pk) will be positive, negative or nil, depending on whether the given stress state σ* lies outside, inside or on the original branch of the positive octant. This observation is fundamental for the following developments, as it suggests that T(σ*,pk) can actually be adopted as an alternative to the evaluation of the yield function f(σ*,pk).
In mathematical terms, the equation of the surface translated along the trisectrix of the positive octant can be easily obtained by modifying equation (2.2) as 3.2 where T is a positive quantity when the translation is in the negative direction of the mean pressure axis. After trivial algebraic manipulations, equation (3.2) can be rewritten as 3.3 where the coefficients of the equation Ii, i=1,2,3, are the invariants of the characteristic equation of the stress tensor σ (e.g. ) 3.4 and principal stresses reduce to 3.5
The cubic equation (3.3) can be solved for T, yielding in general three different solutions T(σ,k). However, the positive branch is the furthest away from any negative stress state, and therefore we are only interested in the solution which yields the major translation T.
The cubic equation (3.3) is of the form 3.6 the coefficients being 3.7 By applying Tschirnhaus transformation 3.8 the depressed cubic equation is obtained 3.9 where 3.10 The latter transformation resulting from the relationship between the invariants of the stress tensor defined by equation (3.4) and those of the deviatoric stress tensors, J2D and J3D 3.11 where 3.12 and the deviatoric stress tensor is defined as , I indicating the second-order identity tensor.
By adopting the trigonometric approach to the solution of a cubic equation, the three solutions of equation (3.6) are 3.13 and x0≥x1≥x2. It should be noted that equation (3.13) has real terms only if u<0 and the argument of the is bound between [−1,1], which requires that 3.14 While the former condition in our case is always satisfied as u=−J2D and J2D is by definition always positive, the latter condition is not always satisfied as equation (3.14) is equivalent to 3.15 Therefore, we need to include also the only real solution of equation (3.6) when u<0 and 4u3+27v2>0 3.16 Finally, an additional solution needs to be included for the case of u=0. Recalling the definition of Lode's angle 3.17 and setting 3.18 the complete solution of equation (3.3) which provides the largest value of the translation T is therefore 3.19 The third of the solution indicated in (3.19) being valid for the case of u=0 , i.e. J2D=0, and can be easily retrieved from equation (3.2) by imposing an isotropic stress state. If the translation T(σ,pk) in equations (3.19) is set to zero 3.20 a mathematical description of the only branch of figure 4 which lies in the positive octant is obtained 3.21 T(σ,k)=0 can therefore be used for defining the yield surface f=0 of a constitutive model. Moreover, as T(σ,k)≠0 provides a translated instance of the same surface, not only f=0 but also f≠0 is convex, which considerably improves the efficiency of the numerical integration of the constitutive law.
It should be noted that the same result can be obtained by writing the invariant I3, in equation (2.3), as a function of the invariants I1, J2D and J3D, solving the resulting equation for I1 and selecting the solution which provides the largest value of I1. This approach will be adopted in the plane stress version of the tension cut-off.
Finally, the failure criterion can be easily extended to materials with a uni-axial tensile strength σt, by replacing I1 with in the previous equation.
4. Plane stress version of the no-tension criterion
A version of the proposed tension criterion for analyses in plane stress conditions can be easily formulated. If we assume that the principal stress component σ3 is nil, then a tension cut-off curve in the σ1−σ2 plane is given by 4.1
As shown in figure 6, equation (4.1) plots as two symmetric curves in the first, positive, and in the third, negative, quadrants. Convexification of the curve can be achieved by following an approach similar to that described in the previous section and translating the positive curve along the diagonal of the first quadrant. However, a slightly different approach is described here, based on the definition of the two stress variables which are often adopted in Soil Mechanics, and which represent the σ coordinate of the centre and the radius of Mohr's circle for plane stress conditions 4.2 By inverting equation (4.2), the principal stresses can then be expressed as 4.3 and equation (4.1) rewritten as 4.4 However, as we are only interested in the curve of the positive quadrant, the solution to our problem is given by 4.5 and, equivalently, the no-tension curve in plane stress conditions is given by 4.6 and its graphical representation is shown in figure 6b.
It should be noted that the plane stress version of the proposed tension criterion can be generalized to materials characterized by a uni-axial tensile strength σt≥0, by reformulating equation (4.1) as 4.7 which can be made convex by following the same procedure which lead to equation (4.5), resulting in 4.8 or, equivalently, 4.9 A representation of equation (4.9) is shown in figure 7.
For implementation purposes, it is convenient to rewrite equation (4.9) in terms of invariants, I1 and , resulting in 4.10 while an explicit formulation can be obtained by solving the previous equation for I1, yielding 4.11
Equations (4.10) and (4.11) are not completely equivalent, and the latter only should be used for implementation purposes. The reader is referred to appendix A where a detailed discussion of the problems associated to equation (4.10) is presented.
5. Validation against analytical solutions
The difficulty of analysing boundary value problems involving no-tension materials has been discussed by many authors (e.g. [12,16,17]) and numerous approaches have been proposed, often very complex, some of them based on ad hoc formulations. For instance, in  a method is proposed which is based on the complementary energy. The advantage over other similar solutions, which are usually slow to converge and unstable, is related to the adoption of augmented Lagrangian techniques to regularize the functionals. However, as reported by the authors, the method is computationally expensive. Another example of an ad hoc approach can be found in , where a modified procedure is suggested for the construction of the global tangent stiffness matrix.
In this section, the possibility is investigated of handling the numerical difficulties reported in the literature at the constitutive level, hence allowing for the use of a standard displacement-based FE method. For this purpose, a linear elastic-perfect plastic constitutive model is adopted, in which the yield and the plastic potential surfaces are described by the proposed failure criterion. The model has been implemented in two different versions, for general three-dimensional and plane stress conditions, respectively. While in the former case, a standard implicit implementation has been adopted (e.g. ), in the latter the mathematical development described in appendix B has been used. Both versions of the model have been implemented in Abaqus, a commercial FE programme, and employed to analyse no-tension problems often used in the literature as benchmark tests (e.g. [12,16,17,19,20]).
These benchmark tests are particularly severe as they evaluate the performance of the analysis involving a no-tension material with continuous displacement and stress fields, with discontinuous displacement and stress fields and, finally, with continuous stress and discontinuous displacement fields. More importantly, for these benchmark tests, the exact analytical solutions are available  and therefore the performance of the analyses can be rigorously evaluated.
(a) Hollow cylinder with prescribed internal and external pressure
Numerical analyses of an infinitely long, thick-walled cylinder loaded with internal and external compressive pressures is considered first. The same geometry, boundary conditions and model parameters adopted in [16,17] have been used in the present analysis and therefore results can be directly compared. The analysis is carried out under the assumption of plane strain conditions and due to the double symmetry of the problem only one quarter of the cylinder has been analysed. The mesh, shown in figure 8, includes 288 eight-noded, plane strain elements, with quadratic shape functions and reduced integration (CPE8R in ABAQUS). The inner and outer radii of the hollow cylinder are ri=0.2 m and re=1. m, respectively, while the elastic behaviour is characterized by Young's modulus E=100×106 Pa and Poisson's ratio ν=0.1. The one-dimensional tensile strength given by the parameter σt of the model was set to 1×10−4, while the rounding parameter was k=1×10−3 Pa3. It should be noted that, as a plane strain analysis is here carried out, the three-dimensional version of the proposed failure criterion is here employed. Therefore, no constrains exists on the value of the rounding parameter k. The applied internal and external compressive pressures are pi=1000 Pa and pe=220 Pa. A comparison between the analytical and numerical distributions of the radial and hoop stresses as a function of the radius of the tube is shown in figure 9a,b, respectively. The figures show that the numerical analyses were capable of accurately reproducing the analytical solution. As expected radial stresses are compressive throughout the thickness of the tube, while hoop stresses are zero up to a radius r0=0.641 m which coincides with the analytical value, and thereafter compressive. No numerical difficulties were encountered during the FE analyses nor problems of slow convergence, as usually reported when the Galileo–Rankine no-tension model is adopted.
(b) Plane stress panel
The performance of the proposed failure criterion has then been tested to simulate discontinuous displacement and/or stress fields. This benchmark test was originally devised in  and includes the numerical analysis of two different isotropic panels (figure 10). Both panels are considered to have zero Poissons' ratio, so that no coupling exists between normal strains.
Panel A is made of a uniform material, with Young's modulus E=10 000 kPa, while the vertical pressure applied at the top horizontal boundary is 200 and 100 kPa for the left and right half, respectively. In this case, both vertical stress and vertical displacement fields are discontinuous.
Panel B is non-homogeneous as Young's moduli of the left and right half are E=10 000 and 20 000 kPa, respectively. The panel is loaded by a uniform vertical pressure distribution equal to 200 kPa.
Also in this case, the same geometry and boundary conditions of the orginal study  have been used, while the panels were discretized into 1600 plane stress linear elements (CPS4 in ABAQUS). The plane stress version of the proposed model was considered, the tensile strength and rounding parameters being σ0=2×10−3 and k=1.2×10−5 kPa2, respectively.
Figure 11 shows the comparison between the analytical and numerical results for panel A. The exact solution is characterized by an independent behaviour of each part of the panel, with a discontinuity in both vertical displacements and vertical stresses at the geometrical symmetry axis, which is well reproduced by the numerical analysis.
Figure 12 shows the same comparison for panel B. In this case, the analytical solution describes a discontinuous vertical displacement field, while the vertical stress remains homogeneous. Also in this case, the numerical results are capable of reproducing the theoretical behaviour, although due to the displacement-based formulation of FE method, a vertical discontinuity in the stress field at the restrained bottom of the panel is observed.
6. Engineering applications
(a) Circular shallow footing on elastic no-tension soil
One of the major applications of the proposed no-tension criterion is in numerical analyses of boundary value problems involving soils and, more generally, uncemented particulate materials where it can be adopted for defining a consistent, three-dimensional, elasto-plastic tension cut-off. In fact, although such materials cannot sustain tensile stresses, yield surfaces of complex elasto-plastic models adopted for simulating their behaviour usually expand beyond the compressive octant of the stress space, thus enabling tensile stress states. Depending on the engineering problem, numerical analyses can therefore present large volumes of material under tension, making their validity questionable.
An example of such an application is described in this section, where the three-dimensional version of the proposed criterion has been adopted for preventing tensile stresses in a linear elastic soil loaded by a shallow circular footing. Owing to symmetry only half of the soil was considered, and the model, shown in figure 13a, consists of a 50 m wide by 50 m deep domain. This was discretized, as shown in figure 13b, in 457 eight-nodes isoparametric and axis-symmetric quadratic elements with reduced Gaussian integration (CAX8R in Abaqus). A rough interface was simulated between the soil and the 2 m diameter circular footing by applying vertical displacements and restraining the horizontal ones. The elastic behaviour of the soil was characterized by Young's modulus E=100 000 kPa and Poisson's ratio ν=0.2, while the no-tension parameters were σt=0.1 kPa and k=0.005 kPa3, respectively.
Figure 14a shows contours of plastic deformations resulting from the numerical analysis performed with the proposed tension cut-off. It appears that a large volume of soil adjacent to the footing satisfies the criterion, indicating that one or more principal stresses have reached the boundary of the compressive octant of the principal stress space and were prevented from attaining tensile stresses, which on the contrary characterize the truly elastic solution. Despite the different stress state in the domain for the two analyses, its influence on the behaviour of the footing is negligible, as it appears from figure 14b. This figure shows the evolution of the bearing pressure as a function of the normalized settlement for the elastic and for the elastic no-tension soil. As expected, the load-deflection curve for the latter material is less stiff, but the differences appear to be very limited. In general, considering the widespread use of elastic solutions in geotechnical engineering practice, the proposed criterion can provide a useful tool for their validation under more realistic no-tension conditions.
(b) Plane stress wall with openings
The performance of the plane stress version of the proposed tension criterion has been compared against numerical results of a masonry wall with opening first analysed in  and then in . The material behaviour was again reproduced by means of an elastic-perfect plastic model, and, as in the original investigations, Young's modulus and Poissons' ratio were set to E=10×109 Pa and ν=0.2, respectively. The geometry and the loading conditions are indicated in figure 15, while the adopted mesh, consisting of 368 quadratic elements, is shown in figure 16. The tensile strength adopted in the proposed failure criterion was set to σt=20 Pa and the rounding parameter to k=110 Pa2, thus satisfying the constrain given by equation (A 2). Referring to figure 15, two loading steps were applied. In the first, all loads with the exception of the parametrized horizontal load αq0 were applied, while during the second the parameter α was increased up to failure. Figure 17a shows contours of the accumulated modulus of the plastic deformation tensor, which indicates that tensile failure is achieved at the base of the masonry and along an inclined shear band, before and after the opening, respectively. The load factor α versus the horizontal displacement of the top right node is shown in figure 17b and compared with the results of . This indicates that by adopting the proposed failure criterion, the panel behaviour can be successfully reproduced by means of a standard displacement-based FE approach. This is in contrast to the conclusions attained in , which suggest that ordinary FE approaches are unsuitable for analysing such a no-tension problem and that an ad hoc approach, based on a specific construction of the global tangent stiffness matrix is required.
7. Discussion of the numerical performance
The failure criterion has been adopted within an elastic-perfect plastic constitutive model in the framework of rate independent plasticity. When no-tension materials are involved, the consistent Jacobian matrix is ill-conditioned and this may cause numerical problems during the equilibrium iteration of the global Newton loop.
However, with the implemented constitutive model, no severe numerical problems were encountered when the three boundary value problems described in the previous section were analysed, although convergence of the global equilibrium equations was not particularly fast. This can be possibly explained with the small, but non-nil value of uni-axial tensile strength, σt, adopted in the proposed failure criterion to reproduce material behaviour.
An alternative approach to the evaluation of the consistent Jacobian matrix is to operate with an elastic global stiffness matrix, i.e. adopting a quasi-Newton approach to solve the equilibrium equations. Unfortunately, this does not guarantee quadratic convergence. As a result, analyses performed with this approach did not differ significantly from the previous one.
It should be noted that the problem of preventing a stress state from trespassing the boundaries of the compressive octant bears similarities with a contact problem. Based on this observation in , an augmented Lagrangian approach has been proposed to regularize the adopted hybrid stress and displacement-based functional. As an alternative, contact problems are often numerically solved using a penalty method, and this can be applied to no-tension problems operating at the constitutive level via rate-dependent plasticity. The Perzyna model  has been here adopted to regularize the equilibrium equations. This requires very limited changes to the implementation of the constitutive model described in the appendix. This only requires changing the consistency condition f(Δε,Δλ)=0 with ( f(Δε,Δλ)/τ0)nΔt=Δλ, where τ0 and n are two material parameters which control the time-dependency response. As in visco-plasticity, the consistency condition can be violated, so that f>0 is acceptable, this results in a penalty method.
The plane stress masonry wall, previously described, has also been analysed using this method. An appreciable improvement in terms of rate of convergence was achieved, while the numerical results were unchanged.
A tension failure criterion has been presented generalizing the so-called Galileo–Rankine criterion within the elastic-perfectly plastic model often adopted for reproducing the behaviour of no-tension materials. The formulation can also be used for establishing a tension cut-off in complex constitutive models, such as those used for reproducing the behaviour of soils, granular materials and powders, to prevent tensile stresses.
The criterion is described by means of a very concise equation based on the third invariant of the stress tensor. It can be used as a mathematical approximation of the coordinate planes which delimit the compressive octant of the principal stress space. This sheds new light on the physical significance of the third invariant of the stress tensor, which can then be considered to provide a quantitative measure of how close the stress state is to the boundary of the first octant. Therefore, the third invariant provides a direct indication of the negative semi-definiteness of the stress state.
The failure surface is defined by two parameters. The first parameter is the the uniaxial tensile strength allowing for the criterion to be used for materials characterized by finite tensile strength. The second parameter sets the roundness of the surface, thus ruling how close the coordinate planes are approximated. As the formulation plots as a multi-branch surface, a technique is presented which removes the undesired branches, producing a smooth convex surface.
The paper also presents a two-dimensional version of the failure criterion, which is suitable for plane stress conditions. It has been shown that although different equations can be obtained via algebraic manipulations, only one is suitable for achieving a very fast convergence when the constitutive law is numerically integrated by means of a full backward Euler scheme.
The new criterion, used as a component of a linear elastic-perfectly plastic no-tension model, has been numerically validated against analytical solutions of two boundary value problems involving no-tension materials, which are often adopted in the literature as benchmark tests. Two additional engineering applications have demonstrated its effectiveness in solving geotechnical and structural boundary value problems.
At a constitutive level, the model based on the proposed formulation appeared to be extremely fast and convergence is achieved in few iterations even when large strain increments are applied. The problem of preventing a stress state from trespassing the boundaries of the compressive octant bears similarities with a contact problem. Therefore, a penalty approach can be adopted for avoiding ill-conditioning of the global stiffness matrix. As no-tension is here simulated by means of a plastic model, a viscous regularization has been adopted based on the Perzyna model , which improved the rate of convergence of the global equilibrium equations.
The authors wish to thank Prof. Rosati and Prof. Valoroso for kindly providing their data.
Appendix A. Plane stress criterion: a note of warning on the formulations in terms of invariants
Equations (4.10) and (4.11) are not completely equivalent. If plane stress conditions are satisfied, the argument of the square root in equation (4.10) is always positive, while in equation (4.11) this occurs only if the value of the rounding parameter k and of the tensile strength σt satisfy the condition A 1 or, equivalently, A 2
For no-tension materials with σt=0, the choice of the rounding parameter k is unrestricted, and, as shown in figure 18a, for any given value both equations (4.10) (full line) and (4.11) (dashed line) produce the same failure curve. However, if a non-nil tensile strength σt>0 is set, so that the rounding parameter k does not satisfy the condition of equation (A 2), the two curves are translated in the negative stress plane but the domain of definition of equation (4.11) is a subset of that of equation (4.10) and the resulting curve is incomplete. Setting k=σ2t/4 in equation (4.11) results in a curve defined in the whole relevant domain, however a more rounded criterion is obtained (figure 18a). The increase in roundness as the tensile strength σt increases, is shown in figure 18b, when k=σ2t/4.
It should also be noted that for equation (4.10) any f≠0 contour with zero tensile strength is equivalent to a f=0 contour drawn for an appropriate value of σt≠0, as shown in figure 19a,b (full lines). This does not occur with equation (4.11), as σt also appears in the argument of the square root. As a consequence, f=0 contours for different values of σt and f≠0 contours for the same σt are not identical parallel curves. However, it can be easily shown that for both equations (4.10) and (4.11) f≠0 contours constitute an homothetic expansion of the f=0 curve, which is a convenient feature for achieving high rate of convergence in the integration of the constitutive law.
Equation (4.11) is preferable to equation (4.10) for implementation within a constitutive model. In fact, an explicit functional dependence of I1 on q allows for a very efficient numerical integration of the constitutive law within a backward Euler scheme. After some algebraic manipulation, the Newton's loop of the stress point algorithm can be made to iterate on a single equation only, which results in an extremely fast convergence (see appendix B).
More importantly, the argument of the square root in equation (4.10) is always positive no matter the value of k only when plane stress conditions are satisfied. In this case, however, the out of plane deformation ε3 is an unknown rather than a control variable and during the iterations, before convergence is achieved, the stress state is in general not plane. Hence, the argument of the square root in equation (4.10) is not necessarily positive which can result in numerical problems. By contrast, equation (4.11), which requires condition (A 2) for the positivity of the argument of the square root independently of the stress state, is not affected by that problem.
Finally, the representation of the plane stress failure criterion in the q versus I1 plane is convex only when equation (4.11) is adopted. Figure 20a shows that for no-tension materials with σt=0, the failure curve, f=0, is convex for both equations (4.10) (full lines) and (4.11) (dashed lines). However, f≠0 contours of the implicit equation are not convex, as indicated in figure 20b (full line). The explicit formulation of the failure criterion given by equation (4.11), on the other hand, is always convex, provided that the condition given by equation (A 2) is satisfied.
This indicates that the original formulation of the Galileo–Rankine criterion in plane stress conditions, which is given by equation (4.10) when k approaches zero, is convex in the q−I1 plane for f=0, but not for f≠0, which provides and additional explanation for the numerical difficulties often reported in the literature.
Appendix B. Plane stress criterion: implementation for a backward Euler integration scheme
As most FE programmes adopt the classical Continuum Mechanics sign convention, whereby tensile stresses are considered positive, in this appendix signs are reversed as compared to the rest of the manuscript, and the classical sign convention is adopted.
A constitutive model for plane stress analyses based on the proposed failure criterion can be very efficiently integrated by means of a full backward Euler scheme following the procedure outlined in . However, it is shown here that after some algebraic manipulation, the internal Newton can be made to iterate on one equation in one unknown only, which results in an extremely efficient and fast integration, even when very large strain steps are applied.
As previously discussed, the explicit formulation in terms of invariants I1 and q of the plane stress failure criterion should be used for implementation purposes. Thus, adopting the classical Continuum Mechanics sign convention, equation (4.11) becomes B 1 with the rounding parameter k respecting the condition set by equation (A 2).
In plane stress conditions, the only non-vanishing out-of-plane strain component is ε3, which is not set by kinematic conditions and needs to be evaluated at the constitutive level. Therefore, the strain increment vector is B 2 where and are here adopted to indicate known and unknown quantities, respectively. Similarly, in the volumetric strain increment B 3 the contribution of the known strain components can also be highlighted, resulting in B 4 The deviatoric strain increment vector can then be expressed as B 5 where , i=1,2.
The plastic components of the volumetric and deviatoric strain increments are obtained via the flow rule B 6 and B 7 The isotropic component of the stress state at time t+Δt is given by B 8 where for the sake of simplicity, the subscript t+Δt is omitted, with the understanding that all quantities are evaluated at the end of the time step, and B 9 The deviatoric component at the end of the load increment is B 10 which is based on the result that the plastic strain tensor can be written as B 11 where B 12 and I and s are the identity and the deviatoric stress tensors, respectively. Moreover, we indicate with B 13 the analogous of the stress tensor predictor, which, however, cannot be directly evaluated as it is a function of which is unknown. As g is not a function of Lode's angle θ B 14 and equation (B 10) can be expressed as B 15 Substituting equations (B 13) and (B 7) into (B 15) leads to B 16 which allows the evaluation of the invariant q of both sides of the previous equation B 17 and solving for q gives B 18 The unknown out-of-plane strain increment can be evaluated from the plane stress condition σ3=s3+I1/3=0 by substituting the s3 component of the deviatoric stress tensor given by equation (B 16) and considering the q definition of equation (B 18) B 19 Finally, substituting equation (B 8) into the last expression, a relationship between Δε3 and Δε1,Δε2,Δλ is obtained B 20
Once Δε3 is evaluated by means of the previous equation, the stress invariants, I1 and q, at the end of the time-step t+Δt, are a function of Δε1, Δε2 and Δλ only. They can therefore be evaluated by means of equations (B 8) and (B 18), where B 21 B 22 B 23 Therefore, the only unknown Δλ can be obtained by imposing f=0.
- Received July 25, 2014.
- Accepted September 25, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.