A systematic framework for formulating convex failure envelopes in multiple loading dimensions

The failure envelope approach is widely used to assess the ultimate capacity of shallow foundations for combined loading, and to develop foundation macro-element models. Failure envelopes are typically determined by fitting appropriate functions to a set of discrete failure load data, determined either experimentally or numerically. However, current procedures to formulate failure envelopes tend to be ad hoc, and the resulting failure envelopes may not have the desirable features of being convex and well-behaved for the entire domain of interest. This paper describes a new systematic framework to determine failure envelopes – based on the use of sum of squares convex polynomials – that are guaranteed to be convex and well-behaved. The framework is demonstrated by applying it to three data sets for failure load combinations (vertical load, horizontal load and moment) for shallow foundations on clay. An example foundation macro-element model based on the proposed framework is also described.

A failure envelope for a specific shallow foundation concept can be determined by measuring, or computing, a set of discrete failure load combinations for different loading conditions (e.g. the experimental data reported in Martin & Houlsby (2000) and the numerical data in Taiebat & Carter (2000)). These failure envelope data are then approximated using an appropriate mathematical formulation (referred to in this paper as a 'failure envelope formulation') together with an optimisation process to achieve a best fit.
In the current paper, V, H and M refer to the applied vertical and horizontal loads, and moment, respectively. Furthermore,H ¼ H=H 0 ,M ¼ M=M 0 andṼ ¼ V =V 0 refer to the normalised loads, where H 0 , M 0 , V 0 are the respective uniaxial capacities. (Uniaxial capacity refers to the failure load for one load component when all other load components are maintained at zero throughout the loading.) In the notation employed in this paper, vectors are denoted in bold face, italic, lower case (e.g. x) and matrices are denoted in bold face, upper case, upright font (e.g. A).
Current failure envelope formulations can generally be categorised into one of two groups (although the wideranging literature in this area also includes a number of formulations that do not fit into either group, e.g. Taiebat & Carter (2000)). One group represents the failure envelope using a sum of power functions of the applied loads. This approach has been used to represent the 'cigar-shaped' failure envelopes determined from experimental data for shallow foundations on sand (Nova & Montrasio, 1991;Gottardi & Butterfield, 1993;Cassidy, 1999;Gottardi et al., 1999;Bienen et al., 2006) and clay (Martin & Houlsby, 2000). Published formulations of this type are mainly defined for planar VHM loading, although extensions to six degrees-of-freedom loading have been described by Martin (1994), Byrne & Houlsby (2005) and Bienen et al. (2006). An example of a formulation in this group is 'Model A' (Martin, 1994), which was developed to represent the failure envelope for a spudcan foundation in clay. Model A is defined as although for numerical implementation, Martin & Houlsby (2001) recommend the alternative form f ðH;M;Ṽ Þ ¼ ðH 2 þM 2 Þ 1=2 À 4Ṽ ð1 ÀṼ Þ ¼ 0 The second group of failure envelope formulations, termed 'HM-based', has primarily been used in connection with numerically determined failure envelope data (Feng et al., 2014;Vulpe et al., 2014;Vulpe, 2015;Shen et al., 2016Shen et al., , 2017. Formulations in this group are defined in terms of composite horizontal and moment loads, where the composition provides a means of incorporating the influence of additional load components within a basic HM framework. An example of a formulation in this group (from Vulpe et al. (2014) whereH Ã andM Ã are composite horizontal and moment loads defined byH Ã ¼H=ð1 ÀṼ 4Á69 Þ andM Ã ¼M=ð1 ÀṼ 2Á12 Þ; and a and b are parameters (a = 2·13 and b = À0·09 forṼ 0Á5, or a = 1·63 and b = À0·05 forṼ . 0Á5).
Formulation of foundation macro-element models Failure envelope formulations are a convenient basis for the development of plasticity-based macro-element models for foundations. For example, based on experimental data, Martin & Houlsby (2000) develop a modified version of Model A (termed 'Model B', described further in the later subsection within 'Example applications' entitled 'Model B') as the yield function for a macro-element model of a spudcan foundation in clay.
A macro-element model typically incorporates a yield function and a plastic potential function; these functions govern the admissible load states and the plastic displacements, respectively. The yield function may also contain a hardening parameter that governs the evolution of the yield surface (see Appendix for discussion on the distinction between a yield surface and a yield function). However, the current paper is concerned mainly with developing failure envelope formulations to provide an approximate fit to discrete failure load data. For simplicity, for the macroelement formulation described later in this paper, associated flow (i.e. where the yield function is also the plastic potential function) and perfect plasticity (i.e. no hardening) are adopted.
The yield function adopted in macro-element models should ideally be convex throughout the domain and it is insufficient to have convexity at just the admissible loads boundary (i.e. a convex yield surface). Franchi et al. (1990) has discussed the key differences between the yield function and the yield surface and the theoretical relevance of the convexity of the yield function. A summary of the characteristics of convex functions and the difference between a yield surface and a yield function is given in the Appendix. The importance of convexity of the yield function has been largely overlooked, especially in geotechnics, over the past few decades. This particular issue was raised by Panteghini & Lagioia (2014, 2018a, 2018b, who reported numerical problems (e.g. lack of convergence in boundary value problems) when using implicit integration algorithms with non-convex yield functions. To avoid such problems, they propose an approach for providing convexity to the yield functions. For rate-independent, non-hardening/softening materials, a convex yield function has the added benefit of thermodynamic consistency, provided that an associated flow rule is adopted and the yield surface contains the origin of the load space (Ottosen & Ristinmaa, 2005).
Yield (and plastic potential) functions employed in macroelement models must be continuous and real-valued. Ideally, they should also be differentiable with a continuous gradient and Hessian throughout the load domain. Furthermore, they should not have a restricted domain, as it is not known a priori the load states at which implicit integration algorithms will require the evaluation of the yield function. Also, functions should be avoided that have singularities (either in the function itself or its derivatives) at certain load states. For the purposes of the current paper, a function with these ideal attributes is referred to as 'well-behaved'. It is noted that successful plasticity models can be developed on the basis of yield functions that are not well-behaved; the conventional Mohr-Coulomb yield function, for example, provides a common basis for soil constitutive models despite not being differentiable at the 'edges' of the yield surface (although rounding approximations (e.g. Sloan & Booker, 1986) can be employed to address the numerical difficulties associated with these edges). In general, employing non-well-behaved functions is inconvenient from a numerical implementation perspective.
Failure envelope formulations can be used in two separate types of practical application: (a) assessment of ultimate foundation capacity and (b) representing the yield surface (and the plastic potential for associated flow) in macroelement models. Failure envelope formulations that are suitable for both purposes are therefore particularly convenient. However, not all of the current failure envelope formulations are applicable to both applications, either due to lack of convexity or because they are not wellbehaved. For example, the extensive use of fractional exponents for the HM-based formulations means that the yield function may not be real-valued in parts of the domain (e.g. equation (3) is not real-valued forM , 0 or V , 0, since raising negative numbers by fractional powers results in complex values). Furthermore, the composition approach within this group of formulations results in singularities in the function and its gradient atṼ ¼ 1, which makes this form of failure envelope inconvenient for macro-element model implementation. It is noted that the HM-based formulations were not developed with macroelement modelling in mind, and thus suitability for this use should not necessarily be expected. However, failure envelopes developed using the formulation framework described in the current paper are guaranteed to be convex and wellbehaved; they are therefore applicable to a range of practical applications.
A new framework for convex, well-behaved, failure envelopes For the purposes of the current paper, a 'globally convex' function is defined as a function that is convex throughout the domain R n (real coordinate space of n dimensions, where n is interpreted in the context of the current work as the number of individual load and moment components applied to the foundation). A new framework is proposed for the formulation of failure envelopes that are guaranteed to be globally convex and well-behaved for all possible failure load combinations. Furthermore, the framework is systematic and general. This is significant, as previous researchers (e.g. Gourvenec, 2007) have suggested that difficulties in deriving closed-form expressions to approximate failure envelopes present a practical obstacle to the adoption of the failure envelope approach in design. The proposed framework provides a systematic process that can be applied to failure envelope data sets consisting of multiple load (and moment) components. The framework is also able to represent various failure envelope shapes that correspond to different types of shallow foundation. It therefore facilitates the further development of the failure envelope approach for practical applications.
In developing this framework, it is recognised that proving convexity for a general function is typically a difficult task. For example, a twice differentiable function with a positive semi-definite Hessian throughout its domain is guaranteed to be convex. However, a heuristic approach that demonstrates that a Hessian is positive semi-definite at a finite number of points is insufficient to prove convexity for the entire domain. Conversely, demonstrating non-convexity is often relatively straightforward; a single case where the Hessian is not positive semi-definite is sufficient.
The approach employed in the current framework is therefore to restrict the search for failure envelope formulations to well-behaved polynomial functions that are known a priori to be globally convex. In developing this procedure, use is made of a mathematical technique known as 'sum of squares' (SOS) programming, which provides a computationally tractable method for determining globally convex polynomials.

POLYNOMIAL FAILURE ENVELOPE FORMULATIONS General homogeneous polynomials
Ellipsoids are a class of well-behaved functions that are guaranteed to be convex, as outlined below. A general ellipsoid can be expressed in the quadratic form where x, c [ R n are vectors and A is a positive definite square constant matrix. The Hessian in this case is r 2 f x ð Þ ¼ 2A and therefore guaranteed to be positive definite. Since positive definite matrices are also positive semi-definite, it follows that ellipsoids are convex. As ellipsoids are also well-behaved, they are candidates for representing failure envelopes, for which the components of x are the foundation loads (in an appropriate normalised form) and the components of A and c are selected to provide a best fit with the available failure load data. However, ellipsoids have the significant disadvantage in this application that they are not expressive enough to represent the varied shapes of failure envelopes that correspond to the range of plausible shallow foundation configurations.
Ellipsoids are second degree polynomials; higher degree polynomials provide options for developing failure envelope functions with increased expressiveness. Since polynomials of odd degree ! 3 are guaranteed to be non-convex (Ahmadi et al., 2013), odd polynomials are excluded as candidates for the current application. Polynomials of even degree ! 4 can either be convex or non-convex, but it is typically computationally difficult to prove their convexity in the general case (Ahmadi et al., 2013). Thus, this paper considers a subset of polynomials of even degree ! 2, called sum of squares convex (SOS-convex) polynomials, for which convexity can be demonstrated in a computationally tractable manner.

Sum of squares polynomials
where g j (x) are polynomials of degree 4 d and npoly is the number of individual polynomials, g j (x), employed in the summation. Equivalently, s(x) is SOS if and only if it can be represented as s(x) ¼ z T Bz where B is a positive semi-definite matrix and z is a vector of all monomials of degree up to and including d. This can be shown by noting that (5) and Q T Q is the Cholesky decomposition of B). It is evident from equation (5) that s(x) ! 0 for x [ R n . A necessary and sufficient condition for a twice differentiable polynomial f (x) to be convex is for its Hessian r 2 f x ð Þ to be positive semi-definite (see Appendix). This requires that A subset of polynomials that satisfy equation (6) is defined by the restricted but more tractable condition Polynomials f (x) that satisfy equation (7) are convex (although convex polynomials exist that do not satisfy equation (7)). The subset of convex polynomials that satisfy equation (7) is referred to as 'SOS-convex' polynomials (Ahmadi & Parrilo, 2012).
A key benefit of employing SOS-convex polynomials for failure envelope applications is that the requirement in equation (7) can readily be incorporated within the search process for failure envelope formulations using semi-definite programming (Parrilo, 2003). Additionally, SOS-convex polynomials are well-behaved and support general n-dimensional loading. Thus, SOS-convex polynomials form the basis of the failure envelope formulations in the proposed framework.

DETERMINATION OF SOS-CONVEX POLYNOMIAL FAILURE ENVELOPES
Procedures to determine SOS-convex polynomials to provide a fit with failure load data sets are outlined below.

Standardise the data
Failure envelope data typically consist of sets of discrete failure load combinations. These failure load combinations are standardised bȳ where x i is a load component, x i,ref is a reference load and x i,c is a shift parameter. The purpose of equation (8) is to shift and scale the data such that the values ofx i for uniaxial loading in the positive and negative directions (according to the selected coordinate system) are 1 and À1, respectively. Geometrically, this ensures that the failure envelope intersects eachx i axis at ± 1. Importantly, if fx ð Þ is convex, then convexity in f (x) is preserved as the mapping in equation (8) is affine (Boyd & Vandenberghe, 2004).
Define the failure envelope functional form SOS-convex polynomial functions of the standardised loading variablesx, of even degree 2d where d is a positive integer, are now sought to represent the failure envelope fx ð Þ. Although SOS-convex polynomials are in general nonhomogeneous, homogeneous forms are adopted in the current framework as they have the important advantage that the coefficients for uniaxial loading can be identified straightforwardly, as indicated in step 3 below.
As an example, a homogeneous polynomial fx ð Þ with degree 2d = 4 and number of loading dimensions n ¼ 2, with coefficients a i (that are subsequently determined to ensure that the polynomial is SOS-convex) is The failure envelope is defined by fx ð Þ À 1 ¼ 0.

Apply the uniaxial conditions
The standardisation process in step 1 requires that the coefficients of all monomials comprising a single loading variable are 1. (This is to represent failure correctly for uniaxial loading conditions.) For the example in equation (9), this condition gives a 1 ¼ 1 and a 5 ¼ 1.

Apply symmetry principles
Other coefficients can be determined by exploiting any available symmetry conditions. For the example in equation (9), if the problem is symmetric such that fx 1 ;x 2 ð Þ¼f Àx 1 ;x 2 ð Þ , then the coefficients of all monomials that are odd inx 1 must be zero; in the current example this implies a 2 ¼ 0 and a 4 ¼ 0.

Set up a convex optimisation problem
To identify the remaining unknown coefficients (there is only one unknown a 3 in the current example but a greater number of coefficients will typically need to be determined at this stage), a convex optimisation problem is set up. This is based on the conditions: (a) fx ð Þ is SOS-convex and (b) fx ð Þ provides a best fit with the failure envelope data by minimising the objective function C wherex data i is a set of failure load combinations and ndata is the size of the data set. The optimisation problem is therefore To solve the SOS program defined in equation (11), it must first be converted into an equivalent semi-definite program. In the current work, the Matlab toolbox 'YALMIP' (Löfberg, 2004) was used to convert SOS problems such as equation (11) automatically to their equivalent semi-definite programs (Löfberg, 2009). Thereafter, the unknown coefficients in f ðxÞ can be readily determined using this toolbox.

EXAMPLE APPLICATIONS
Three examples are provided below to demonstrate the application of the proposed formulation framework. The first example concerns the failure envelope for a surface circular footing on homogeneous clay; the second example is the failure envelope for a suction caisson foundation embedded in homogeneous clay; the last example is a new failure envelope formulation to represent Model B, which demonstrates that an existing failure envelope formulation can be approximated within the proposed framework to achieve the advantages of global convexity and well-behaved function attributes.
In all the examples, YALMIP was used, in conjunction with the SeDuMi semi-definite solver (Sturm, 1999), to solve the SOS programs. All the failure envelopes for the examples are plotted in the normalised load space (e.g.HM). This enables convenient comparison with figures in previous publications (for the case of Model B).

Surface and caisson foundations
The first two examples use the VHM failure load combinations for surface and caisson foundations in homogeneous clay, determined as described in Suryasentana et al. (2019) using the 'sequential swipe test' in three-dimensional (3D) finite-element analysis. The soil in these cases was modelled as a linear elastic-perfectly plastic von Mises material. Fig. 1 shows the foundation configurations being considered, including the assumed location of the loading reference point (LRP) (which is needed to provide a datum for the definition of the applied moment).
Based on the numerical failure load data, standardised forms of the failure load data ½H;M;V are obtained, from equation (8) Two separate SOS-convex polynomial failure envelopes have been determined: f 4 (with degree 2d ¼ 4) and f 6 (with degree 2d ¼ 6). The general forms of these polynomials are Application of the uniaxial condition gives a 1 , a 5 , a 15 ¼ 1 in f 4 and a 1 , a 7 , a 28 ¼ 1 in f 6 . Additionally, symmetry inH andM requires that a 6 , a 7 , a 8 , a 9 , a 13 , a 14 ¼ 0 in f 4 and a 8 , a 9 ,  a 10 , a 11 , a 12 , a 13 , a 19 , a 20 , a 21 , a 22 , a 26 , a 27 ¼ 0 in f 6 . The remaining coefficients are determined from the optimisation in equation (11). Results. The coefficients for f 4 and f 6 determined from this process are listed in Tables 1 and 2, respectively. Figs 2 and 3 compare the normalised VH, VM and HM failure envelopes determined by f 4 À 1 ¼ 0 and f 6 À 1 ¼ 0 for the surface and caisson foundation, respectively; also shown are the corresponding 3D finite-element results from Suryasentana et al. (2019). Differences between f 4 and f 6 can be seen more clearly in Fig. 4 (which omits the finite-element results).
Both f 4 and f 6 are seen to provide reasonably good approximations of the computed data, although there is a small under-prediction of the VM failure envelope for the suction caisson (Fig. 3(b)) and f 4 is unconservative compared to f 6 for the surface foundation at high vertical loads (Fig. 2(c)). It is seen that f 6 provides a better fit with the numerical data than f 4 for the surface foundation (Fig. 2(c)), but f 4 is marginally better for the caisson (Fig. 3(c)). This is confirmed by the data in Table 3, which shows values of the objective function (and their root-mean) at the end of the optimisation process (the smaller the better).
Failure envelopes have been determined on the basis of f 4 and f 6 (i.e. equations (12) and (13)) incorporating the same uniaxial and symmetry conditions as adopted in the earlier section 'Surface and caisson foundations' examples. The remaining coefficients were determined using the optimisation in equation (11).
The computed coefficients are listed in Tables 4 and 5. Fig. 5 shows the normalised VH, VM and HM failure envelopes for f 4 À 1 ¼ 0 and f 6 À 1 ¼ 0, together with comparisons with the original Model B equation. Evidently, both f 4 and f 6 provide reasonably good approximations of the Model B equation (although neither of them predicts the correct value ofṼ at which the peaks in theṼH andṼM envelopes occur). Table 6 shows the minimised objective function values (and their root-mean) for f 4 and f 6 at the end of the optimisation; these data suggest that f 4 provides a better fit than f 6 . Figs 5(a) and 5(b) show that, atṼ ¼ 0 andṼ ¼ 1, the failure envelopes implied by f 4 and f 6 are more 'rounded' than Model B, and the surface normal is directed along theV axis. This means that, when employed as a plastic potential, the model correctly predicts vertical plastic displacement for pure vertical loading; this is significant, as special numerical procedures are required to achieve such behaviour with Model B (Martin, 1994). In general, the results show that the SOS-convex Model B formulation provides a failure envelope that is reasonably close to the formulation in equation (15), although it is evident from Fig. 5(c) that the SOS-convex formulation is unconservative for relatively high values of normalised vertical load; the SOS-convex formulation, however, is guaranteed to be globally convex and well-behaved.
Increasing Ṽ Macro-element model for a surface foundation To demonstrate the application of the proposed SOSconvex failure envelope formulations in a macro-element model, the 4th degree version of the failure envelope (equation (12)) was employed, using the coefficients for a surface foundation in Table 1. Failure envelopes of the form fH;M;V À Á such as equation (12) are implicit functions of the loading variables H, M, V by way of the affine transformation in equation (8) The failure envelope f 4 (H, M, V ) À 1 ¼ 0 was implemented as the yield surface and plastic potential of a macro-element model, with the elasto-plastic integration performed using the implicit closest point projection method (Simo & Hughes, 2006). The uniaxial capacities (V 0 , H 0 , M 0 ) were determined from the finiteelement data in Suryasentana et al. (2019), while the elastic stiffness matrix for the macro element was calibrated using pre-failure data obtained from the finite-element results.
The macro-element model is used, as an example, to simulate a sequential swipe test in the HM load space for planar HM loading withṼ = 0·25, to replicate the corresponding finite-element simulations in Suryasentana et al. (2019). Fig. 6 provides a comparison of the loaddisplacement behaviour and the normalised HM failure envelope, computed using the macro-element model and the 3D finite-element model (Suryasentana et al., 2019). It is evident that the macro-element predictions agree well with the finite-element results.

DISCUSSION
The procedures presented in the paper provide a convenient means of generating failure envelopes that are convex and well-behaved. This is achieved by selecting functions from a restricted set of SOS-convex polynomials. The process is demonstrated by applying it to previously published numerical data on failure load combinations for a surface footing and a caisson on soil that is modelled as an elasticperfectly plastic von Mises material. Additionally, it is used to reformulate a previously determined failure envelope, Model B.
The results indicate that an SOS-convex polynomial of degree 6, f 6 , does not necessarily provide a better fit to a prescribed data set than an SOS-convex polynomial of degree 4, f 4 . The results indicate that f 6 is a better fit for the surface foundation failure envelope than f 4 , but the opposite is true for the caisson and Model B failure envelopes. This may seem counterintuitive, but it is a consequence of the fact that  Fig. 6. Comparison of load-displacement and the normalised HM failure envelopes predicted by the macro-element and 3D finite-element results (Suryasentana et al., 2019), for a sequential swipe test in the HM load space whenṼ = 0·25. S h and θ m refer to the horizontal and rotational displacements of the surface foundation, respectively. D is the diameter of the surface foundation the current approach is restricted to the use of homogeneous polynomials. Although a general polynomial can be used to represent any polynomial of lower degree, this is not the case for the homogeneous polynomials employed here. For practical applications of this method, it is recommended that a relatively low degree homogeneous polynomial (e.g. degree = 4) is initially adopted. Higher degree polynomials would then be considered if the low degree polynomial is found to provide an unsatisfactory fit to the data.
The YALMIP toolbox provides a useful SOS decomposition feature to double check that the failure envelopes are, indeed, convex throughout the entire domain. For example, for f 4 employed for Model B, the Hessian can be decomposed into SOS form using YALMIP as This equation demonstrates thatȳ T r 2 f 4x ð Þȳ is SOS, hencē y T r 2 f 4x ð Þȳ ! 0 for its entire domain; f 4 is therefore shown to be SOS-convex.
The values of the coefficients in the failure envelope polynomials determined using the proposed framework do not necessarily have a physical interpretation. This disadvantage is outweighed by the new approach being able to represent a wide range of failure envelope shapes, while maintaining desirable theoretical attributes. A simple way to interpret the polynomial coefficients would be to view their (absolute) magnitudes as indicative of the strength of coupling between the respective loading components at failure.
The framework presented in the paper has certain limitations, as follows. The application of the SOS-convex failure envelope formulation to macro-element modelling has only been demonstrated for perfect plasticity with associated flow. Further work is needed to develop the modelling approach for the case of non-associated flow and hardening plasticity. Furthermore, the current study only explores the application of homogeneous polynomials. Non-homogeneous SOSconvex polynomials are not explored, but they potentially provide greater expressiveness than homogeneous SOS-convex polynomials (since homogeneous polynomials are subsets of non-homogeneous polynomials). However, based on initial analyses using non-homogeneous polynomials, the SeDuMi solver employed in the current work encountered convergence problems for some of the above examples; this contrasts with the homogeneous polynomials where convergence issues did not occur. Further work is needed to explore the implementation of non-homogeneous SOS-convex polynomials and their potential benefits.

CONCLUSIONS
This paper introduces a novel SOS-convex polynomialbased framework to systematically formulate failure envelopes that are suitable for use in both ultimate capacity calculations and macro-element modelling. The approach is applicable to failure envelope data sets of multiple loading dimensions, and it is demonstrated in the current paper for VHM loading. The principal advantage of the proposed framework is that it leads to formulations for failure envelope functions (and yield functions) that are guaranteed to be globally convex and well-behaved. As discussed in this paper, these attributes lead to significant advantages in terms of theoretical considerations as well as the numerical stability of macro-element models.
Examples are presented to demonstrate the application of the framework to three different shallow foundation cases, with significantly different failure envelope shapes. These examples demonstrate the expressiveness and general capabilities of the proposed framework. The framework rationalises and expedites the formulation of failure envelopes and therefore supports the development of the failure envelope approach for design applications.
Although not shown here, this framework may have potential uses in other fields of solid mechanics, such as determining convex yield functions from experimentally determined yield data for different materials.