A flexible and robust yield function for geomaterials

This paper presents a new Modified Cam Clay (MCC) type yield function, that is designed for robust and efficient use with mplicit stress integration algorithms. The proposed yield function attains non-elliptical (e.g. tear and bullet) shapes, as well s the typical elliptical shape of the MCC model. Like that of MCC, and unlike most other yield functions with non-elliptical hapes available in literature, it is non-singular and unique throughout stress space. The experimental yielding stresses of a wide ange of geomaterials have been accurately simulated using the yield surface. The yield function can be used in constitutive odels based on classical elasto-plasticity theory. c 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license http://creativecommons.org/licenses/by/4.0/). eywords: Flexible; Robust; Non-singular; Uniqueness; Yield surface; Implicit stress integration


Introduction
To numerically simulate the mechanical performance of large-scale engineering infrastructures, constitutive equations that represent the mechanical behaviour of materials are used. These can be implemented in boundaryvalue solvers, such as the finite element method or the material point method, using several approaches (see [1][2][3] for detailed information), of which the most efficient belong to the implicit stress integration algorithm category [1,3]. In this category, after the boundary-value solver has calculated a strain increment, a corresponding (trial) stress state is calculated and, when needed, this is mapped back onto the yield surface. More specifically, the trial stress state may temporarily go beyond the admissible stress states that are bounded by the yield surface. When this happens, the stress (located outside the yield surface) is returned back onto the yield surface using the plasticity component of the constitutive equations. This return mapping algorithm is generally an iterative procedure, with the stress being updated in each iteration. Usually, the return mapping procedure requires the first derivatives of the yield function and plastic potential with respect to the updated stress state. Therefore, it is essential for robust computation that computational methods employ a yield function (and plastic potential) that is defined uniquely in the entire stress space to ensure the updated stress is correctly mapped onto the yield surface [4]. Otherwise, this may result in Stress tensor at time step n + 1 φ c Friction angle limitations on the time/load steps used, or worse, to non-convergence or incorrect converged stresses which impairs the predictive capability of the constitutive equations [4].
Owing to a wide range of behaviours observed for geomaterials, researchers have paid particular attention to developing yield surfaces and plastic potentials with a variety of shapes in meridian and deviatoric stress space, in order to capture their behaviour as accurately as possible. Examples can be found in [5][6][7][8][9][10][11][12][13][14]. These constitutive equations have been derived from different theories, with classical plasticity and thermodynamics principles being commonly accepted examples.
As will be seen later, the Modified Cam Clay (MCC) [15] yield function and plastic potential are inherently robust and therefore highly suitable for implicit stress integration algorithms. However, although the MCC yield function therefore provides a suitable framework for developing constitutive equations, its elliptical surface does not represent the observed behaviour of many geomaterials. In fact, MCC has no flexibility to cope with experimentally determined loci of yield stress points with non-elliptical shapes (e.g. so-called tear, bullet and egg shapes). Many alternative surfaces can be found in the literature that have addressed this limitation (see for example [5,7,[11][12][13]16]). However, most of them have non-unique elastic domains, i.e., locations in stress space outside the intended elastic domain where the yield function value implies elasticity, a so-called undesired (false) elastic domain or nucleus. This was investigated by Coombs & Crouch [17] for the yield surface proposed by Collins & Hilder [11]. One may use the convexification technique to overcome this problem by developing an alternative formulation of the yield function which resembles and preserves the non-elliptical shape of the yield surface at all values; for example, Stupkiewicz et al. [18] used such an approach to extend the formulation of a previously proposed yield surface [13]. Moreover, Panteghini & Lagioia [19] corrected the undesired elastic domain of a distorted MCC-type yield surface with a single shape-parameter via a convexification technique [20], to develop a fully convex and unique yield function expression. However, in both of these works, the analytical expression for the yield function is complex and the approach they pursue may not be practical or convenient for all yield surfaces, as imaginary solutions may arise.
In this study, a new yield function is proposed which has the flexibility of attaining non-elliptical shapes. As for MCC, it is robust for implicit stress integration algorithms, without having the shortcomings of false elastic domains, singularities or erratic gradients in meridian stress space. Note that the formulations presented are in accordance with geotechnical conventions, where compressive stresses and contractive strains are considered to be positive.

Requirements of the yield function
Yield surface functions are generally defined by the effective stress tensor, σ, and hardening variables, q, which may be tensors (such as kinematic hardening variables) or scalars (e.g. isotropic hardening variables). Here, as is common in geotechnical engineering, yield surfaces defined in meridian stress space, i.e. by stress invariants p and q (respectively, the mean effective stress, p = tr(σ)/3 (kPa) and deviatoric stress, q = (3/2 s:s) 1/2 (kPa), where s = σ − tr(σ)/3:1 is the deviatoric stress tensor and 1 is the second order identity tensor), and specifically those with isotropic hardening, defined by a pre-consolidation pressure p c , are considered. This is particularly relevant for history dependent materials such as soils. A function y is defined which is used to control the constitutive behaviour. Hereafter, the terms yield surface and yield function are distinguished by level sets of y [4]. The yield surface is defined by the level set "zero", i.e., y(σ, q) = 0 and the yield function is defined for "non-zero" level sets, i.e., y k (σ, q) = k where k > 0.
In elasto-plasticity, a yield surface is employed to divide the stress space into three regions [4]: • stresses inside the yield surface, that represent the elastic domain, y < 0; • stress states on the yield surface, y = 0, that undergo elasto-plastic deformation; • inadmissible stresses outside the yield surface, which satisfy y > 0.
Yield surfaces should uniquely define the admissible/elastic stress space (y ≤ 0), i.e., there should only be one yield surface. Moreover, the yield surface is required to be convex over its admissible stress space, which means that stresses cannot go outside the yield surface from one point inside the yield surface to another point inside the yield surface. When the yield surface is implemented in a boundary-value solver with an implicit stress updating procedure, additional restrictions for the yield function are needed. When used in boundary-value solvers, usually a so-called return mapping procedure is used in order to correctly calculate the stress and hardening parameters associated with a certain strain (increment) (see Fig. 1). Here, a closest point projection method (CPPM) is employed as the return mapping approach, where it is assumed that for a time interval [t n , t n+1 ], at time step n, the stress tensor, σ n , strain tensor, ε n , hardening variables, q n , and yield surface (y n ), are converged and defined, and that the stress state is located on the yield surface. To calculate the corresponding variables at time step n+1, the return mapping technique applies two steps: (i) by the use of the total strain increment (obtained from the boundary-value solver), ∆ε, an incremental stress predictor, ∆σ predictor , is calculated, assuming only elastic behaviour. Consequently, the trial stress is calculated (σ trial = σ n + ∆σ predictor ); (ii) when σ trial goes beyond the yield surface (i.e. y > 0), the stress state must be corrected (mapped) back onto the yield surface. This is performed by the application of a stress increment corrector, ∆σ corrector . In this step, the plastic behaviour is estimated, giving a new prediction of the updated stress, the yield surface and the hardening variables. This may be done iteratively, in order to return the new stress (σ n+1 ) onto the updated yield surface (y n+1 ), updated via its hardening variables (such as p c ). As a result of step (i), σ trial sits on a non-zero level set of the yield function (i.e. y(σ trial ) > 0), and σ corrector in step (ii) is typically calculated iteratively by using the first and second derivatives of the yield function (y j n+1 , j = 1,2,3,. . . ) (in the case of an associated flow rule) or plastic potential function (g j n+1 , j = 1,2,3,. . . ) at every state that the return path experiences. To have a correct return mapping path, the yield function (y k (σ, q) = k) should meet the criteria considered for the yield surface (y = 0), including uniqueness and convexity. It should be noted that these restrictions are due to numerical implementation considerations and not material behaviour/energy considerations. Moreover, even if the yield surface satisfies these conditions, it is not automatic that the yield function also satisfies them. It is also possible to have mathematical expressions that provide the same yield surface but very different yield functions in terms of their level set variation.
The majority of this paper focuses on the uniqueness of yield functions and how to develop a robust and efficient form for implicit stress updating algorithms. Non-uniqueness may appear from several sources. Three of these are identified and their detrimental consequences summarised below: • Non-unique elastic domains: in the stress space outside the desired yield surface (y(σ, q) = 0) it may be possible to identify domains where the condition y > 0 is not satisfied, i.e., the yield surface defines two (or more) enclosed areas where the condition y < 0 is satisfied; this is a so-called false (undesired) elastic domain or nucleus ( Fig. 2(a)). Potentially, this means that the stress predictor could move the trial stress into the false elastic domain; then the behaviour is falsely predicted to be elastic, and consequently would not trigger the return mapping procedures and, therefore, not calculate plasticity, material hardening or the return to the yield surface.
• Singularity: in the stress space outside the yield surface, the yield function or its derivatives may not be mathematically defined, which is hereafter called a "singularity" (Fig. 2(b)). As a result, there are sudden  [13]); (b) singularity and discontinuity in the gradient for same level sets (yield function of [11]); (c) high curvature of the yield function (yield function of [11]). changes in the behaviour of the yield function (i.e., discontinuities), for example in its gradients, resulting in divergence of the solution. • Erratic gradients: in the stress space outside the yield surface, high curvatures of the yield function may be observed ( Fig. 2(c)), even though it is mathematically defined at those stress states (in contrast to singularity points). In these regions, convergence can be difficult.
Other issues, including local minima (the stress state getting trapped and being unable to return to the yield surface), also result in an inefficient stress integration algorithm. Although these issues are not covered in this paper, they may be linked to the sources of non-uniqueness mentioned above and these issues may be solved by following the approach proposed here.

Analysis of flexible elliptical yield functions
Many constitutive modelling approaches and yield surfaces have been proposed for geomaterials. Collins & Hilder [11], using the principles of thermodynamics, derived a family of isotropic yield surfaces for geomaterials (readers may refer to [11,21] for detailed descriptions of the derivation procedure). These yield surfaces have a quadratic-type formulation, similar to that of an ellipse: or, alternatively, in the form Both these forms will be used later to identify the geometrical constraints required to formulate a unique and non-singular yield surface. Eq. (1) is the zero level set of the general yield function (y k ): or, alternatively, where k ≥ 0 represents the level set of interest.
A and B are, respectively, the semi-major and semi-minor axes of the yield surface ( Fig. 3) and C controls the horizontal distance (or shift) between the semi-minor axis of the surface and the origin (denoted by point O in Fig. 3). All the surfaces in Fig. 3 are bounded on the hydrostatic axis at p c (the pre-consolidation pressure) and the centre of the surface (point c) is obtained as the intersection of the semi-major and semi-minor axes. A 1 is the horizontal distance of the centre of the yield surface from the tensile apex (which is here coincident with the origin, so that A 1 = C), whereas A 2 is the horizontal distance between the centre of the yield surface and the compression apex of the yield surface (p c ). Hence the sum of A 1 and A 2 is the size of the major axis of the yield surface (i.e., A 1 + A 2 = p c ). For an ellipse ( Fig. 3(a)), it holds that C/(A 1 + A 2 ) = 0.5 (i.e., A 1 /A 2 = 1) and the gradients (demonstrated by red arrows) at the co-vertexes (shown by cv 1 and cv 2 ) are zero with respect to the p-axis. Flexible non-elliptical shapes can be formed by several approaches. Here, two of those possible methods are discussed.
The first method fixes the semi-minor axis on the p-axis in the same place as an ellipse (i.e. A 1 = A 2 ) ( Fig. 3(a)) and at this point the yield surface is inclined to the p-axis. Two such surfaces are depicted in Fig. 3(b) and (c). Note that, for these surfaces, the gradients with respect to the p-axis at points cv 1 and cv 2 are non-zero (in contrast to the ellipse). In the second approach, the location of the minor-axis on the p-axis is changed compared to an ellipse (and the gradient at this point is maintained at zero). This implies that the ratio of A 1 /A 2 changes with respect to elliptical shapes (i.e. A 1 /A 2 ̸ =1). Fig. 3(d) and (e), respectively, illustrate yield surfaces with A 1 /A 2 < 1 and A 1 /A 2 > 1, that have been developed via this approach. For these surfaces, similar to an ellipse, the gradients at points cv 1 and cv 2 are zero with respect to the p-axis. A combination of these two methods may also be employed to develop non-elliptical surfaces of an appropriate form to present the material under consideration.
To mathematically incorporate these surfaces in Eq. (1), A, B and C are defined in terms of stress invariants (for example, the hydrostatic pressure p) and hardening variables (for example, the pre-consolidation pressure p c ), i.e., A = A(p, p c ), B = B(p, p c ) and C = C(p c ). Therefore, hereafter they are called 'stress-like functions'. Note that the value of A = A(p, p c ) varies between A 1 and A 2 depending on p. By modifying the definition of the stress-like functions (A, B and C) and using the general quadratic-form of the yield surface presented in Eq. (1), non-elliptical (flexible) MCC-type yield surfaces compatible with the observed yielding stresses in geomaterials have been proposed [5,11,12]. These models employ two-parameter (α and γ ) functions to define A, B and C, in order to form 'tear' and 'bullet' shape MCC-type surfaces. The functions for four yield surfaces of the MCC-type family, applicable to soils with zero tensile strength, are summarised in Table 1. Regardless of the different formulations proposed for these functions, they all comply to certain geometrical restrictions. Assuming an isotropic yield surface, the yield surface at zero deviatoric stress, q, should be defined at two points, (p, q) = (−p t , 0) and (p, q) = (p c , 0), where p c and p t are, respectively, the compression and tensile apex of the yield surface along the hydrostatic axis (p-axis), and p t = 0 for soils. Substitution of these two points in Eq. (1) results in Fig. 3 and shows that the functions A and C are inter-related in such a way that function C is the same as function A when the dependency on p is omitted. Thus, they have similar signs. In addition, Eq. (4) represents A 2 in Fig. 3.
Assuming that C represents the hydrostatic pressure at Critical State (CS) conditions, i.e.
where subscript 'cs' indicates a CS-related variable, the third geometrical constraint is then derived by considering that the yield surface is required to be defined at (p, q) = (p cs , q cs ) = (p cs , Mp cs ): where M is the CS stress ratio (the gradient of the CS Line (CSL) in p-q stress space). By substituting Eq. (5) in Eq. (6) and combining Eq. (3) with Eq. (4), the following relations are extracted: Note that the dependency of function B on function C (and consequently on function A) in Eq. (7c) demonstrates how the stress-like functions are linked. The above criteria provide a general method to define stress-like functions in order to develop non-elliptical (flexible) MCC-type yield surfaces; i.e., any proposed stress-like functions that Table 1 Different yield surface formulations.

MCC
Collins & Hilder [11] Chen & Yang [12] Zhang et al. [5] Function Roots of the yield function of the form Eq. (2b) Note: N.A. stands for not applicable. satisfy these conditions result in a flexible MCC-type yield surface in meridian stress space (see the stress-like functions for several MCC-type yield surfaces presented in Table 1). Note that Eq. (1) can easily be extended to account for the influence of the third invariant, the Lode angle θ , by multiplying the function B with an appropriate Lode angle dependent (LAD) function ζ (θ ): or, alternatively, in the form A number of these LAD functions are reported in [13,[22][23][24]. Since these functions are independent of p and p c , the generality of the above approach is not compromised. Since the focus of the current work specifically concerns numerical deficiencies in the meridian plane, the incorporation of a LAD function into the yield surface is hereafter ignored.

Analysis of the Modified Cam-Clay (MCC) model
To use the same yield surface as in the Modified Cam Clay (MCC) model, the functions A, C and B are defined as A = C = p c /2 and B = Mp c /2, respectively. Thus, for the MCC model the yield surface (Eq. (1)) becomes Level sets of the MCC yield function in normalised meridian stress space are shown in Fig. 4(a). All the level sets are defined in meridian (p, q) stress space and are concentric with respect to the yield surface, i.e., transformations of the yield surface by an increase in the yield function have the same shape. In addition, the yield function at any level set has identical properties, such as identical gradients. Since the MCC yield surface is unique and convex, its level sets (k > 0) also adhere to these properties and are ensured to be unique and convex, as is obvious from Fig. 4(a). It can also be inferred from Fig. 4(a) that, during the return mapping procedure of an implicit stress integration approach, any stress state outside the MCC yield surface (resulting from the stress predictor) can be correctly returned back onto the yield surface. This is because, during the return mapping, the stress will be located on surfaces that all have their gradients towards the yield surface (y = 0).

Analysis of the general formulation of the yield function (Eq. (2)) and non-elliptical MCC-type models
The stress-like functions of the MCC and three flexible MCC-type yield surfaces [5,11,12] are presented in Table 1. The functions A, B and C defined for MCC are independent of the stress invariants p and q (see Table 1). On the other hand, for the other three yield functions, the shape parameters α and γ are related to p to form non-elliptical MCC-type surfaces. It should be noted that by assigning specific values for the shape parameters, the dependency of these functions on p can be eliminated to resemble the MCC yield function. The respective shape parameters for the yield functions from [11,12] and [5] are (α, γ ) = (1, 1), (α, γ ) = (1, 0) and (α, γ ) = (1/2, 0) to recover the MCC yield surface.
The level set contours of these four yield surfaces are illustrated in Fig. 4. It is observed that all the yield functions (except for MCC) have either singularities or are non-unique for some level sets. Although the defined stress-like functions result in a unique yield surface (for a specifically defined range of α and γ ), the yield functions lose their uniqueness or may be undefined for a specific value of p (due to singularity) at non-zero level sets (see Fig. 4). These drawbacks can be investigated in the roots and singularity points of the general yield function equation (Eq. (2a)) and will be explored in the following subsections. The general formulation of the yield function, Eq. (2a), may have several roots in meridian stress space (p, q), i.e., on the p-axis (when q = 0). The roots are values of p that satisfy the following three equations: The roots corresponding to R 1 , R 2 and R 3 (obtained from Eq. (10)) for each yield function [5,11,12] are presented in Table 1 and are shown, respectively, by solid circle, open circle and red star in Fig. 4.

Uniqueness
By substituting A, B and C for the MCC yield function in Eq. (10), the yield surface (y = 0 and k = 0) intersects the p-axis at two points; (p, q) = (0, 0) and (p, q) = (p c , 0). These points can easily be determined by considering q = 0 in the yield function. In other words, the yield surface in meridian stress space has only two roots (the roots of R 1 and R 2 in Eqs. (10a) and (10b)). Since the function B is independent of p, R 3 has no root, which keeps the total number of roots of the yield surface function as two. When the roots of the yield function at non-zero level sets (for a given k) in meridian stress space are analysed (by solving for p at a given k > 0), again only two roots, , are obtained (shown by the solid and open circles in Fig. 4).
For the other three yield functions, by substituting the stress-like functions in Eq. (10), the roots of the yield function for a given k are extracted (see Table 1). Since these functions are defined by p and p c , there is a possibility that the number of roots (t) for Eq. (2b) can exceed two. In fact, t roots (in Eq. (2b)) enclose the meridian stress space with t − 1 domain(s) and each domain represents an elastic domain (bounded by a yield surface). For the MCC yield function, t = 2 and, therefore, the number of elastic domains is one. However, for the other yield functions, t ≥ 2 (due to the root of R 3 ) and, as a result, there are one (t = 2) or more (t > 2) elastic domains, which means that it may be possible to identify the existence of at least one false (undesired) elastic nucleus, indicating a loss of uniqueness. It is therefore concluded that, to avoid false elastic nuclei, it is required to keep the number of roots of the yield function Eq. (2b) (for a given level set k) equal to two (t = 2).

Singularity (discontinuity)
By substituting q = 0 in Eq. (2a), the yield function reduces to y k = (p−C) 2 /A 2 -(1 + k) = 0. This function is not defined when A = 0, i.e., at p values that fulfil At these values of p, the yield function is undefined (as indicated by the solid red squares in Fig. 4), thereby causing a sudden change in behaviour such as gradients, i.e., resulting in a discontinuity. This type of behaviour is due to the lack of a fully defined yield function in meridian stress space. Of the four yield functions investigated, only the yield functions of MCC and [12] have no singularities (see Table 1).

Observations and conclusions
From the above analysis of both forms of the general yield function (Eqs. (2a) and (2b)), the following observations are made: • To form a non-elliptical MCC-type yield surface, the functions A or B (or both) should be a function of p, and controlled via shape parameters such as α and/or γ .
• The yield function is required to be defined for all values of its shape parameters.
• To preserve the uniqueness of the admissible elastic region in stress space, the number of roots (t) of the yield function Eq. (2b), obtained from Eq. (10), is required to be equal to two.
• To avoid singularities (discontinuities), the yield function Eq. (2a) should be defined for all stresses in meridian stress space.

Proposed yield surface
A countless number of equations can be defined for stress-like functions to form non-elliptical MCC-type yield functions. By accounting for the observations reported in Section 3.5, general geometric constraints are designed which result in the development of a non-elliptical MCC-type yield function that is unique and non-singular, and is robust for a return-mapping technique using the implicit stress integration approach. Two shape parameters, α and γ , are used to define the A, B and C functions, following existing methods. Moreover, an additional parameter, β, which represents the back-stress, is introduced, allowing the yield surface (and its non-zero level sets) to be sheared/distorted with respect to the p-axis. The yield surface is described by a quadratic-type expression which is derived from a thermodynamic basis by defining an appropriate dissipation function (see [11,21,25] for the derivation of a quadratic-type yield surface with β). The yield surface can therefore be incorporated into constitutive models which have either a thermodynamic or a phenomenological basis. The modified yield surface and yield function are respectively expressed as The roots of Eq. (13b), when q = 0 and β̸ =0, are

Geometrical constraints
As observed in Section 3, in order to avoid singularities and to preserve uniqueness, a yield function in the form of Eq. (13b) is required to have only two roots (for a given level set k) on the hydrostatic axis (p-axis) and a yield function in the form of Eq. (13a) is required to be defined over the entire meridian stress space. However, from Eqs. (14c) and (14d), it is shown that there is the possibility to have more than two roots in Eq. (13b), and/or to have discontinuities in Eq. (13a). Note that, amongst these roots, level set values (k) are only incorporated within the roots R 1 and R 2 . To overcome these discrepancies, the general strategy (followed here) is to define functions A and B individually so as to have no roots on the p-axis. In this way, the roots R 3 and R 4 in Eqs. (14c) and (14d), respectively, are directly eliminated and the other two roots corresponding to R 1 and R 2 (for k = 0) satisfy the Eqs. (7a) and (7b) criteria.
Functions A and B are defined using the shape parameters γ and α, respectively, and are proposed to have the following general forms: These functions are specifically selected in order to have no roots (no intersection with the p-axis) and to be always positive, i.e., A, B ∈ (0, ∞). The parameters a 1 and b 1 are determined by satisfying the criteria in Eq. (7), which leads to the stress-like functions being defined as The MCC yield function is retrieved by setting α = γ = 0. The variation of functions A and B for different values of α and γ are presented in Fig. 5(a) and (b), respectively. As can be seen, these functions are always positive, regardless of the (positive or negative) values assigned to the shape parameters. In addition, the stress-like functions asymptote towards zero when p, depending on the value of α and γ , approaches ±∞, but never intersect with the p-axis. This ensures that the functions do not have any roots, thereby avoiding the possibility of non-uniqueness or singularities. Furthermore, R 1 and R 2 (for k = 0) in Eqs. (14a) and (14b), for three different γ , are depicted in Fig. 5(c), showing that, independently of γ , R 1 and R 2 have only one root on the hydrostatic axis (p-axis), intersecting at (p, q) = (p c , 0) and (p, q) = (0, 0) (the origin), respectively. This is similar to the MCC yield surface, but, for different values of α and γ , giving non-elliptical (flexible) MCC-type surfaces.

Convexity and range of shape parameters α and γ
The quadratic yield surface (yield function with level set of zero), Eq. (12a), with the addition of a LAD function, can be re-written as Eq. (20) demonstrates the formulation of the yield surface that is divided into two separate functions defined in the meridian, h(p), and deviatoric, z((q−βp), θ ), planes. Note that q−βp is the radius of the yield surface, at a shear level of β, in the deviatoric plane. For such yield surfaces that are distinguished by meridian and deviatoric functions, through the use of the theorem of convex analysis [26], Bigoni and Piccolroaz [13] demonstrated that the convexity of the yield surface can be separately investigated in meridian and deviatoric stress spaces; i.e., convexity of the functions, h(p) and z((q−βp), θ), implies that the yield surface is also convex. Note that the convexity of the deviatoric function, z((q−βp), θ), is directly related to the form of the LAD function, ζ (θ ). When convex LAD functions such as those proposed in [13,[22][23][24] are used, the convexity of the yield surface in deviatoric space is ensured and, consequently, the convexity of the yield surface reduces to the convexity of the meridian function.
Here, the convexity in meridian stress space is considered, as it relates to the main topic of the paper. The convexity of the meridian function, h(p) (Eq. (20)), is satisfied when h ′′ ( p) ≥ 0. h(p) can be rewritten as Consequently, the convexity of the meridian function results in a condition that satisfies where h ′ 1 ( p) and h ′′ 1 ( p) are, respectively, the first and second derivatives of the function h 1 (p) with respect to p. Although the yield surface is defined for the entire range of values of α and γ , it is only convex if a certain combination of these parameters is used, and it is possible to determine these ranges by applying the convexity criterion (Eq. (21b), which, for the current yield surface, is shown in Fig. 6. In practice, the non-elliptical yielding stresses in geomaterials can be captured by using −2≤α, γ ≤2. The convexity of the yield surface over this range of shape parameters is shown in Fig. 6(b). Note that the above analysis is limited to the convexity of the yield surface (yield function at the level set of zero). The convexity of the yield function for a non-zero level set is investigated in Section 7.

Extension of formulation to account for tensile strength
Geomaterials, such as sandstones and concrete, in general may have tensile strength, i.e., at zero mean effective stress, p = 0, the material has a shear strength of q t . Alternatively, the tensile strength can be characterised by a corresponding tensile pressure, p t ; i.e., the decompression pressure (initial isotropic pressure) is shifted from (p, q) = (0, 0) to (p, q) = (−p t , 0). The stress-like functions can be modified to extend the formulation to account for tensile strength. Hence, the yield surface should intersect the p-axis at p = −p t and p = p c . By applying the aforementioned conditions in Eq. (12a), Eq. (7) becomes To account for p t in functions A and B, the approach of Bigoni & Piccolroaz [13] is followed, where the normalised pressure term (p/p c ) is replaced by (p + p t )/(p c + p t ). Then, by satisfying the criteria in Eq. (23), functions A, B and C are modified to The tensile strength, q t , can be determined by substituting p = 0 in Eq. (12a).

Parametric study of the proposed yield function
The results of a sensitivity analysis of the parameters on the yield surface, Eq. (12), and level set contours of the yield function, Eq. (13), defined by Eqs. (24)- (26), are shown in Figs. 7 to 9.
In geomechanics, the ratio of the mean effective stress at the CS condition, p cs , to the over-consolidation pressure, p c , is often used to calibrate the yield surface. This ratio, r, is called the "spacing ratio" [11,16] and, by using Eqs. (5) and (19), it is calculated as As the spacing ratio is determined from experimental data, Eq. (27) can be used to calibrate the shape parameter γ as The effect of the shape parameter γ on the yield surface with α = 0 is demonstrated in normalised stress space (p/p c , q/p c ) in Fig. 7(a). In the figure β = p t = 0, and the intersection of the CS stress ratio (M) with the yield surfaces are indicated by red dots and their corresponding p/p c represents the normalised mean effective stress at the CS (p cs / p c ,), i.e., indicating the spacing ratio, r. As the parameter γ increases, r increases, as shown by Eq. (27). In addition, the yield surface is stretched with respect to the q-axis as γ increases, which means the material can withstand more deviatoric stress for a given p. Fig. 7(b) shows how the yield surface changes as the shape parameter α varies. At constant spacing ratio, r (or rather, constant γ ; see Eq. (28)), the yield surface attains a tear shape for α̸ =0. This is shown for r = 0.5 (γ = 0) in Fig. 7(b), where the yield surface curvature changes with respect to the CS stress point, (p, q) = (p cs , q cs ). In Fig. 8 the level sets of the proposed yield function for a range of shape parameters α and γ are plotted. Fig. 8(e) with α = γ = 0 resembles the MCC model, in which all the level sets intersect the hydrostatic axis (p-axis) at two points. With the new formulations defined for the stress-like functions (Eqs. (24)-(26)), which have no singularities or roots, the level sets of the proposed yield surface with different shape parameters α and γ , similar to MCC, intersect the p-axis at two points, i.e., the yield function for a given k has only two roots. Therefore, the yield function is defined over the entire domain of p (i.e., there is no singularity) and does not produce any undesired (false) elastic nuclei.
The CS stress ratio, M, indicates the ultimate strength of the material and controls the size of the yield surface along the q-axis. In terms of the CS friction angle of the material, φ c , this parameter can be determined on the compression and extension sides, respectively, as The change in the shape of the yield surface, as parameter M varies, is depicted in Fig. 9(a) for yield surfaces with (α, γ , β) = (0.5, 0.5, 0). It is observed that, as M increases, the size of the yield surface increases, i.e., the yield surface stretches along the q-axis while the p value of the intersection point of the yield surface with the CSL (shown by red dots) remains unchanged. Moreover, the level sets of the yield function ( Fig. 9(b)) indicate that the uniqueness of the yield function is preserved as M changes. Geomaterials, such as sands and clays, may be consolidated under K 0 conditions (i.e. anisotropically consolidated), where the yield surface (and plastic potential) have been observed to be inclined along a non-zero stress ratio (η̸ =0). It is important to capture this feature of geomaterials, as their mechanical behaviour may be significantly influenced by it. The yield surface (Eq. (12)) can be sheared/distorted with respect to the p-axis by the parameter β and a suitable hardening rule can be assigned for it during plastic deformation. β may be related to K 0 by [8] where η K 0 is the stress ratio at K 0 conditions. Shearing/distorting of yield surfaces with shape parameters (α, γ ) = (0.5, 0.5), for β = 0, 0.3, 0.6 and 0.9, as well as the original CSL, are plotted in Fig. 9(c). As the shearing/distortion level (β) increases, the intersection point of the yield surface with the CSL changes. This implies that the CS may not remain unique. However, one may follow the approach of Coombs [10] to relate the shape parameters, α and γ , to β to develop a unique CS yield surface (see Appendix A). Note that the uniqueness of the CSL does not place any requirements on the yield surface formulation, although some formulations of the CSL utilise features of the yield surface. The influence of β on the level sets of the proposed yield function is investigated in Fig. 9(d).
The level sets of the yield function with shape parameters (α, γ ) = (0.5, 0.5) and β = 0.3 (in Fig. 9(d)) indicate   that the proposed yield function maintains its uniqueness as β is varied and does not produce false elastic nuclei or singularities. The influence of parameter p t on yield surfaces with shape parameters (α, γ , β) = (0.5, 0.5, 0) is presented in Fig. 9(e). The yield surface stretches along the p-axis by increasing p t , whereas it shrinks along the q-axis; i.e., the yield surface flattens as p t increases. The level sets of the yield surface, for the same shape parameters and |p t /p c | = 0.3 (in Fig. 9(f)), indicate that the uniqueness is not lost as p t varies.

Efficiency and robustness
This section aims to demonstrate the advantages that the proposed yield surface brings in implicit stress integration algorithms. Therefore, some comparative numerical investigations are presented using MCC, the models of Collins & Hilder [11], Chen & Yang [12] and Zhang et al. [5], and the proposed model. These analyses consider stress changes as might be experienced, for example, due to loading of a soil layer near the ground surface, for a soil close to normally consolidated conditions as might be modelled with a relatively small yield surface. The applied loads may cause large strain increments and, when these are passed to a constitutive model, trial stresses are computed that may be a few times bigger than the hydrostatic extent of the yield surface. All of the yield equations were integrated via an implicit CPP algorithm and common constitutive components, including (hyper)elasticity, isotropic hardening rule and the assumption of associated plastic flow, were used. The complete set of constitutive equations and the numerical implementation are provided in Appendix B. For two types of tear-shape yield surfaces, the stress update algorithm using isotropic non-linear and linear elasticity was subjected to various trial stress states and the number of iterations to find the updated stress state were recorded. The yield surface includes hardening, whilst remaining pinned at the origin of p-q stress space. The model's hardening law is based on the MCC volumetric hardening relationship, such thaṫ where λ and κ are the bi-logarithmic elasto-plastic and elastic compressibility indices, respectively, andε p v is the rate of plastic volumetric straining. Common material parameters based on Lower Cromer Till are provided in Table 2.
The shape parameters to form two types of tear-shape surfaces, for each of the yield equations of Collins & Hilder [11], Chen & Yang [12], Zhang et al. [5], and the proposed form are presented in Table 3. The shape parameters for the models of Chen & Yang [12], Zhang et al. [5] and the proposed form were determined via regression analysis to provide a best fit to Collins & Hilder [11]. For Sections 5.1 and 5.2, the major axis of each yield surface is aligned with the p-axis (that is, β = 0). Analyses involving return mapping when β̸ =0 are investigated in Section 5.3. Table 4 Stress-space iterations for non-linear elasticity: maximum and total number of iterations, number of points that failed to converge (hit the maximum number of iterations), number of elasto-plastic points (trial states outside the yield surface with y > 0) and average number of iterations for elasto-plastic points. The maximum number of iterations was set to 25.

Iteration-stress maps for non-linear elasticity
For non-linear (hyper)elasticity, the Cauchy stress (σ) is linked to the elastic strains, ε e , via where e e is the elastic deviatoric strain, e e = ε e − 1 3 ε e v 1, ε e v is the elastic volumetric strain, G is the shear modulus and p r is the reference pressure.
The trial states were related to the hydrostatic extent of the yield surface, p c , such that the following region of stress space was explored: where they were applied from the same starting stress state, (p, q) = (p c /2,0), for all the models. The strain increment to cause these elastic trial states was determined and applied in a single step. The number of iterations for each of the trial states is shown in Fig. 10, for the two series of tear-shape yield surfaces and the MCC model. A total of 120,400 trial states were explored for each model. Table 4 gives the maximum and total number of iterations, and the number of points that failed to converge (i.e., hit the maximum number of iterations, which was set to 25 for this analysis) for the MCC model and for the four other yield equations forming tear-shape surfaces. The table also provides the number of elasto-plastic trial states (i.e., those trial states that were outside the yield surface and required iterations to update the stress state), and the average number of iterations for the elasto-plastic stress updates.

Tear-shape type I yield surface analysis
For trial stresses with p/p c > 1.5, the yield function proposed by Collins & Hilder [11] encountered difficulties to return back to the yield surface and 10,081 points failed to converge, as shown by the red points in Fig. 10(b), and it is this region that caused the average number of iterations to be double that of the MCC model. The model of Zhang et al. [5] predicts the lowest number of elasto-plastic states (47,423) compared to other models, and is nearly half of those predicted by Collins & Hilder [11], Chen & Yang [12] and the proposed yield equation. This is due to the false elastic region when p/p c > 1.15 (Fig. 10(d)). The model utilising the Chen & Yang [12] yield equation, although on average requiring a lower number of iterations for elasto-plastic states to converge, struggles to converge for trial stress states around p/p c = 0 and q/p c = 1.5 (Fig. 10(c)), or failed to converge (one failure point). The proposed form on average requires only one more iteration than the MCC model, despite having far more control over the shape of the yield envelope, and all the trial states converged.  Table 5 Stress-space iterations for linear elasticity: maximum and total number of iterations, number of points that failed to converge (hit the maximum number of iterations), number of elasto-plastic points (trial states outside the yield surface with y > 0) and average number of iterations for elasto-plastic points. The maximum number of iterations was set to 25. Tear-shape type II yield surface analysis Fig. 10(f) shows the contours of number of iterations required to return the trial stress back onto the yield surface of Collins & Hilder [11]. For trial stresses with p/p c > 1.5, the model faced difficulties to converge and, 157 points failed to return to the yield surface. Similar behaviour is observed for the model of Chen & Yang [12] (Fig. 10(g)) where 319 trial stresses did not converge. The model of Zhang et al. [5] performed the worst, with some trial stresses falsely predicted to be in an elastic domain (false elastic domain indicated in Fig. 10(h)), and 557 points failed to converge. The proposed yield surface performed the best, as all the trial stresses were converged within a maximum of 9 iterations. The total number of iterations required for all trial stresses to return to the yield surface using this yield function (667,343) was nearly 20% less compared to those of Collins & Hilder [11] and Chen & Yang [12], which required 816,518 and 802,520 iterations, respectively. Compared to the MCC model, on average one extra iteration is required to return all of the trial stresses.

Iteration-stress maps for linear elasticity
Here, the previous analysis is repeated but with a linear elastic formulation (that allows negative p values for trial states to be explored). The Cauchy stress is obtained from and I is the fourth-order identity tensor. The following region of p−q stress space was explored: with a total of 180,900 equally distributed trial states. The strain increment to cause these elastic trial states (starting from (p, q) = (p c /2,0)) was determined and applied in a single step. The number of iterations for each of the trial states is shown in Fig. 11, for the MCC model and for the four other yield surfaces for the two different tear-shape surfaces (with the same shape parameters as analysed for the return mapping using non-linear elasticity). Table 5 gives the corresponding results of Table 4, including the maximum and total number of iterations, the number of points that failed to converge, the number of elasto-plastic trial states, and the average number of iterations for elasto-plastic stress updates, for the five yield equations and linear elasticity. Only those models using MCC and the yield function proposed in this paper converged for all the trial stresses. Note that the model utilising the proposed yield equation only required, on average, 1 additional iteration compared to the MCC model.

Tear-shape type I yield surface analysis
It can been seen from Fig. 11(b) that the yield equation of Collins & Hilder [11] failed to converge for some trial states at p/p c < −0.5, and for all trial stresses at p/p c > 1.5 (a total number of 33,792 points shown by the red region). Although the model of Chen & Yang [12] performed relatively well using non-linear elasticity (Fig. 10(c)), the model using linear elasticity struggled to converge for trial states at p/p c < −0.5 and failed to converge for 330 trial stresses. The model of Zhang et al. [5] contains a spurious y < 0 region when p/p c > 1.15, as indicated in Fig. 11(d). Although the yield equation of Zhang et al. [5] has the lowest number of total iterations, this is a spurious result because the yield equation incorrectly predicts that a number of trial states are in a (false) elastic region of stress space (this point is explained in more detail below). In addition, for some trial stresses at −1.00 ≤ p/p c ≤ −0.75 (542 points) the model was not able to converge. Fig. 12(a) shows the distribution of the number of iterations within the stress update algorithm for the MCC model and for the four other models of tear-shape type I, in terms of the proportion of the total number of trial stress states. For example, for the MCC model around 37% of the trial states require 6 iterations to converge. Fig. 12(b) shows the normalised cumulative distribution, which indicates the total proportion of trial states equal to or below the given number of iterations for each of the models. It is clear from the cumulative distribution ( Fig. 12(b)) that the Zhang et al. [5] model spuriously predicts around 40% of the trial states as being within the yield surface, even though the actual percentage is around 12%, due to an additional region where y < 0 for normalised pressures above p/p c ≈ 1.1, as shown in Fig. 11(d). It is also clear from Fig. 12 [12] model with that proposed in this paper ( Fig. 12(b)), the proposed formulation has a shorter tail, with 99% of states converging within 9 iterations, whereas the Chen & Yang [12] model requires 18 iterations before 99% of the states have converged and only 88% have converged within 9 iterations.
Tear-shape type II yield surface analysis 14,249 trial stresses at −1.00 ≤ p/p c ≤ −0.75 failed to converge for the yield equation of Collins & Hilder [11] (shown in red in Fig. 11(f)) and, consequently, this increased the average number of iterations to 8.28. The model of Chen & Yang [12] faced difficulties in converging with trial stresses at p/p c > 1.5 and p/p c < −1.0 ( Fig. 11(g)), resulting in 1004 non-converged stresses. For the studied stress space, two undesired elastic regions were identified for the model of Zhang et al. [5], where trial stresses were not returned back on to the yield surface ( Fig. 12(h)). This is also demonstrated in Fig. 12(d) where 25% of trial stresses were identified to remain in the elastic domain, while the MCC model showed approximately 8.5% of trial stresses were located inside the yield surface. Moreover, 2,244 trial stresses were not returned back on to the yield surface. The proposed yield surface only required 8 iterations to converge for more than 96% of all the trial stresses ( Fig. 12(d)) and 9 iterations for all of the trial states to converge.

Iteration-stress maps for the proposed yield surface at sheared/distorted states
In Sections 5.1 and 5.2, the analysis of returned mapping stresses were limited to isotropic yield surfaces where β = 0. Here, to further demonstrate the robustness of the proposed yield surface for implicit stress integration algorithms, the convergence of trial stresses when the yield surface is sheared at β = 0.0, β = 0.1, β = 0.3, β = 0.5 and β = 0.7 for two types of tear-shape yield surface are investigated. The same material properties as reported in Table 2 and non-linear elasticity are used. The iteration-stress maps are shown in Fig. 13. The shape parameters for the tear-shape type I yield surface analysed in Fig. 13(a) to (e) are α = γ = 0.5, and for the tear-shape type II yield surface in Fig. 13(f) to (j) they are α = γ = −0.5. Detailed analysis, including the distribution and normalised cumulative distribution of the number of iterations are presented in Fig. 14 for both types of tear-shape yield surface.
For both types of tear-shape yield surface sheared at various β values, all the trial stresses were converged with a maximum of 8 iterations; the only exception was for the tear-shape type II yield surface with β = 0, where 9 iterations were needed to return less than 3% of trial stresses back on to the yield surface ( Fig. 14(d)). Moreover, most of the trial stresses were mapped on to the yield surface in 6 iterations for both types of tear-shape surface. Note that the number of elasto-plastic points reduced as β increased. This is because the yield surface bounds a larger portion of the stress space as it shears-off the p-axis.
In conclusion, while other flexible yield surfaces studied in this work have either singularities which result in difficulties in returning trial stresses back on to the yield surface, or have false elastic nuclei, the proposed yield surface is unique, without any singularities, and is demonstrated to be robust for implicit stress integration algorithms for different types of tear-shapes and levels of shearing/distortion.

Comparison with experimental data
The applicability of the proposed yield surface to resemble the yield stress points of a wide range of geomaterials, including clays, sands, carbonate sand, sandstones and gravel, are investigated here. The calibrated parameters for examples of each geomaterial are presented in Table 6.
The yield stresses of two clays, namely Boom Clay [27] and Bothkenner Clay [28], along with the calibrated yield surfaces, are presented in Fig. 15(a) and (b), respectively. The natural Boom Clay samples were initially compressed isotropically to 9 MPa, followed by isotropic unloading to form samples with a wide range of over-consolidation ratios. Then, the samples were sheared under drained triaxial stress paths. The yield stresses were calculated as the intersection of bilinear tangents to the ε v -lnp curves. For Bothkenner Clay [28], undisturbed samples were consolidated under K 0 conditions to resemble in-situ stress states, and then unloaded along the same stress path to form over-consolidated states. Triaxial drained probing tests were then performed and the yield stresses were calculated as the intersection of bilinear tangents of the p−ε v and q−ε s curves. It can be seen that the proposed non-elliptical yield surfaces successfully capture the yield stresses of both Boom Clay and Bothkenner Clay, which have respectively been isotropically and anisotropically loaded.   In Fig. 15(c), fully-saturated sand samples collected from a decomposed granite deposit were isotropically consolidated to 600 kPa and then unloaded isotropically to attain different over-consolidated stress states [29]. Then, the samples were subjected to drained triaxial shear stress paths and the yield stresses were measured by an acoustic emission technique. In this approach, the onset of plastic deformation, i.e., the yield stress, was detected from the sounds produced by the sliding of grains during shearing. The calibrated yield surface fits the experimental data with a good accuracy. Fig. 15. Comparison of the yield surface with experimental yielding stresses for: (a) Boom Clay [27]; (b) Bothkenner Clay [28]; (c) decomposed granite sand [29]; (d) Aio Sand [30]; (e) carbonate sand [31]; data taken from [32]; (f) Rothbach Sandstone [33,34]; (g) Darley Dale Sandstone [33,34]; (h) Pancrudo Slate gravel [35]. Fully-saturated dense Aio Sand [30] samples ( Fig. 15(d)) were loaded and unloaded anisotropically, and then subjected to several drained triaxial stress paths to determine the yield stress points. The yield stresses form a non-symmetric and non-elliptical shape which the proposed yield surface captures accurately with the calibrated parameters.
Many offshore-related structures are constructed on carbonate sands. These geomaterials exhibit distinctive behaviour, as they have a similar volumetric behaviour to fine-grained soils (e.g., clay) while their shear behaviour follows that of sands [32]. This feature is mainly due to carbonate sands being formed from the skeletal remains of marine organisms. The yield stress points and the calibrated yield surface are presented in Fig. 15(e), which shows the capability of the proposed yield surface to bound the non-elliptical shape of the experimental data.
The effectiveness of the proposed yield surface to represent the onset of shear compaction of sandstones is demonstrated for Rothbach Sandstone and Darley Dale Sandstone [33,34] in Fig. 15(f) and (g), respectively. The experimental data form non-elliptical envelopes that indicate the initiation of compaction during shearing. These envelopes are well represented by the calibrated yield surface. It should be noted that tensile pressures have been accounted for in simulating the yield surfaces.
Alonso et al. [35] determined the yield loci of Panncrudo Slate gravel by performing a single multistage triaxial test. The yield locus in normalised stress space is presented in Fig. 15(h). The yield stress points form a non-elliptical yield surface, for which the proposed yield function captures its curvature via its calibrated parameters.
7. Discussion on the potential loss of simple convexity at higher level sets The uniqueness of the proposed yield function has been proven and demonstrated in Sections 4 and 5. By referring to level set contours, for example in Fig. 8(c), one may argue that for non-zero value shape parameters, the convexity of the yield function may be lost as the level sets increase; i.e., although the yield surface is (simply) convex, an increase of the curvature of the yield surface at a higher level set may result in non-convexity (so that the yield function cannot be considered to be fully convex). For example, the range of shape parameters resulting in convex surfaces at level set k = 30 is shown in Fig. 16, where, compared to Fig. 6 for k = 0, a narrower range is indicated.
The loss of convexity does not affect the ability of the return mapping algorithm for the proposed yield function. In fact, it is still robust and is generally efficient, although the level of efficiency is likely to be reduced. Here, to demonstrate the robustness of the yield function for an implicit stress integration approach, a return mapping analysis for large strain increments is investigated.
For the simulations, the material parameters G = 55 (MPa), λ = 0.05, κ = 0.01, and shape parameters α = γ = 0.5 were used, and p c and β were respectively set to 200 kPa and 0.5. A linear elastic and associated plastic flow constitutive behaviour, similar to Section 5.2, were considered. Starting from an isotropic stress state at the reference pressure, (p, q) = (p r ,0), the model was subjected to trial stresses (in a single step) which were related to the hydrostatic extent of the yield surface, p c , such that a large region p trial / p c ∈ [0, 5] and q trial / p c ∈ [0, 5] in the stress space (p − q) was explored (with p trial and q trial being as large as 1000 kPa).
The level sets of the yield surface up to 1000, in the studied stress space, are shown in Fig. 17(a) and the map of the number of iterations to return the trial stresses back on to the yield surface is presented in Fig. 17(b). It is obvious that for the considered shape parameters (α = γ = 0.5), as the level sets of the yield function increase, they become more non-convex. However, as can be seen in Fig. 17(b), even for large trial stresses that sit on such non-convex level sets, they are returned back on to the yield surface, albeit with an increased number of iterations (although all trial stresses are converged within a maximum of 13 iterations).
The implicit stress integration is an iterative procedure in which, in each iteration, the updated stress converges towards the solution that satisfies the CPP energy equation by using the gradients of the yield surface at that stress value. This means that, in each iteration, the updated stress sits on a lower level set with a lower curvature intensity, and eventually returns back onto the yield function with a level set of zero. Because the intensity level of the curvature of the yield surface reduces iteratively, convexity is retrieved as the level sets reduce. Therefore, the iteration solution of the implicit stress integration algorithm ensures the correct returning path of the updated stress onto the yield surface. Thus, the proposed yield function, while not being fully convex (with respect to non-zero level sets), facilitates reliable return mapping paths for the trial stress.

Conclusion
A new yield function is proposed for geomaterials which is non-singular and uniquely defined in meridian stress space, and is therefore suitable and robust for implicit stress integration algorithms. The components and possible roots of the general yield function in meridian stress space were analysed in depth, and geometrical constraints are designed in such a way that the yield function preserves its uniqueness and non-singularity at any level set (value). Compared with other yield surfaces studied here, the proposed form of the yield function presented in this paper is robust and efficient for return mapping for large strain increments. It does not have undesired elastic domains, and therefore avoids the possible prediction of spurious false elastic behaviour or domains with erratic and divergent gradients, which could have resulted in a high number of iterations required to return stress states back onto the yield surface. Moreover, the yield function can attain elliptical and a wide range of non-elliptical shapes compatible with the experimental data of geomaterials available in literature.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. and are added to the previous increment values. The iteration procedure continues until the residual values are within an acceptable tolerance range. The components of the Jacobian matrix are  The derivative of the yield surface with respect to p c is For an associated flow rule: