Use of interaction domains for a displacement-based design of caisson foundations

Caisson foundations, typically adopted for both onshore and offshore structures, are usually subject to combined loading acting during working conditions and exceptional events such as earthquakes. Assessment of their performance under general loadings is therefore fundamental, for both serviceability and ultimate limit states. In this study, a simplified displacement-based approach, aimed at preliminary designing caisson foundations subjected to combined loading, is presented. Such an approach requires the definition of both interaction domains (IDs) and generalised pushover curves, together with the assumption of an associative flow rule. The IDs and pushover curves are obtained by interpreting the results of a set of 3D finite element nonlinear static analyses, where the response of massive cylindrical onshore caisson foundations, embedded in a layered soil profile and subjected to both centred vertical (N) and combined loads (N–Q–M), is investigated. Following previous works, the influence of initial loading factor and caisson embedment ratio on both IDs shape and size is investigated. Additionally, the effect of soil drainage conditions on the IDs is discussed. Role of load reference point (LRP) is also assessed, since a suitable choice of LRP may strongly simplify the geometrical representation of the ID. Analytical expressions for dimensionless IDs and pushover curves are presented and used at a preliminary design stage to evaluate the maximum generalised load acting on the caisson for a given threshold generalised displacement, so as not to exceed either serviceability or ultimate limit states.


jFj
Dimensionless generalised force as defined in Eqs. (2) jFj lim Limit value of the dimensionless generalised force juj Dimensionless generalised displacement as defined in Eqs. (2) juj el ¼ jFj lim =K 0 Elastic value of juj corresponding to jFj ¼ jFj lim A c Caisson cross-section area a n , a l Ellipse semi-axes [Eqs. (9)] C Elastic compliance matrix of the soilcaisson system c' Soil effective cohesion c' red Soil reduced effective cohesion according to [10] c 11 , c 12 , c 13 Interpolating parameters [Eq. (8) Poisson's ratio n = Q/N lim,net Dimensionless horizontal force r v (z=H) Total vertical stress at z = H in lithostatic conditions u 0 Soil internal friction angle u 0 red Soil reduced internal friction angle according to [10] v = 1/F Sv Initial loading factor w Soil dilatancy angle x Ellipse rotation angle [Eq. (8)] 1 Introduction Caissons are embedded foundations characterised by large mass and stiffness, typically employed for both offshore and onshore structures, and usually subject to combined loading acting during working conditions (weights of superstructure and traffic, wind, etc.) as well as exceptional events (strong earthquakes, storms, etc.). According to an ultimate limit state approach, a foundation must be designed so that the load combinations expected during its service life are lower than the collapse loads, thus guaranteeing an adequate factor of safety. However, if the collapse load is attained for displacements incompatible with serviceability of the superstructure, a displacement-based approach is recommended when evaluating the soil-foundation system capacity, as in the case of large-diameter piles. Such an approach may be particularly useful for caisson foundations, as failure mechanisms involve large soil volumes due to their noticeable dimensions, leading to high axial and lateral loading capacity. This need is supported by these foundations being typically used for critical facilities characterised by tall superstructures, such as long-span bridges and wind turbines, whose limit conditions are typically defined in terms of generalised displacements and permanent rotations rather than generalised forces.
In view of the above, developing a simplified procedure aimed at computing the limit load acting on caissons for a given threshold displacement may prove useful, as an expeditious tool to be used at a preliminary design stage. Such a procedure would require the definition of: (1) the combined N-Q-M failure domain of the soil-foundation system, where N and Q are the vertical and horizontal forces and M is overturning moment; (2) the generalised load-displacement curve describing the response of the system from the onset of loading until its failure condition. For foundations subject to combined loads, interaction diagrams (IDs, i.e. three-dimensional failure envelopes in the N-Q-M space) are a useful tool to relate the different loading components at failure. The use of IDs would allow to define the factor of safety either as the minimum distance of the current N-Q-M combination from the envelope or the distance evaluated along the load path [16,27].
Interaction diagrams have been investigated for a variety of foundations by employing different approaches (i.e. experimental, analytical and numerical): shallow footings [7,9,18,30,33], solid and skirted shallow foundations for offshore structures typically characterised by an embedment ratio H/D \ 1 (H being the embedment depth and D the in-plane dimension) [5,37] and spudcan footings employed for jack-up units [23]. Less studies have been devoted to failure conditions of massive caisson foundations for onshore structures subjected to combined loads, concerning either cuboid-shaped [2,17,38] or cylindrical caissons [3,8], both with values of H/D C 1. Furthermore, most of the numerical studies [4,5,18,19,37] are based on 2D plane-strain numerical analyses, while just a few account for the three-dimensional stress/strain state properly [3,8,17,38]. Under undrained conditions, the analyses are always carried out in terms of total stresses [3,5,17,18,37], where the soil is described as an equivalent single-phase medium, therefore ignoring its two-phase nature.
In this study, the IDs of massive cylindrical caisson foundations are obtained for different embedment ratios H/ D = 0.5, 1, 2, installed in a layered alluvial deposit under both centred vertical and combined N-Q-M loads. The IDs are computed through a series of 3D finite element (FE) nonlinear static analyses carried out assuming both fully undrained and drained conditions. To account for the twophase nature of the saturated soil, the undrained analyses have been performed in terms of effective stresses, rather than the typically adopted total stress approach. For the sake of simplicity, a linear elastic-perfectly plastic behaviour is assumed for the foundation soils. This choice, in agreement with what already done by many other authors also in the recent past [3][4][5]19], is justified by the theoretical objectives of the paper, aimed at defining the mechanical response of the system, irrespective of the specific datum. In fact, when monotonic loads are applied and hydro-mechanical coupling is disregarded, the use of more sophisticated constitutive relationships could modify the values but not the nature of the system mechanical response.
In the first part of the paper, after giving insights on the problem layout (Sect. 2) and the numerical model (Sect. 3), the results of the numerical analyses are presented and discussed (Sect. 4). The FE results are validated against those coming from the application of limit analysis (LA) and those discussed by [5,17,38] for a homogeneous soil. Additionally, the choice of considering a two-layer foundation soil showed that, for a given profile of shear strength (either constant or linearly increasing with depth), the size of IDs is governed by the stratigraphic profile but not its shape: this allowed to generalise the obtained results to a large number of practical applications in which a colluvial stratum is positioned on a cohesive one. Following previous works by [17,38], influence of foundation geometry (H/D) and initial loading factor (v) on the IDs have been first investigated, showing a good agreement with their results. Furthermore, the comparison of the IDs obtained in the two limit drainage conditions is provided for a given geometry and initial loading factor, showing that drainage scales the ID size rather than influencing its shape. A short discussion on how the choice of the load reference point (LRP) affects the shape of the IDs is also given, with the caisson centroid being the most suitable choice as it strongly simplifies the geometrical representation of the ID.
In the second part of the paper, based on the results of the nonlinear static analyses, the authors propose an approach for the assessment of the performance of caisson foundations: (1) an analytical expression for ID in the nondimensional N-Q-M space for both drainage conditions (Sect. 5) and (2) the equation of dimensionless generalised pushover curves (Sect. 6). The expressions (1) and (2) are integrated in a displacement-based simplified approach; moreover, the hypothesis of an associative flow rule, confirmed by previous experimental and numerical works [35,38], allows to compute the displacement components, for an assigned load path and a given threshold for generalised displacements.

Problem layout
IDs have been obtained by performing nonlinear static (i.e. pushover) analyses on different cylindrical caisson foundations subject to combined N-Q-M loading until the collapse is attained. In Fig. 1a, a sketch of the adopted geometrical scheme is illustrated. A rigid cylindrical caisson of diameter D = 12 m and height H is embedded in an alluvial deposit consisting of a 5-m-thick layer of gravelly sand and a layer of silty clay. Water table is located at the bottom of the gravelly sand layer: an initial hydrostatic pore water pressure regime is imposed. While the diameter of the caisson is kept constant over the parametric study, three embedment ratios H/D = 0.5, 1, 2 are considered. In this study, the load reference point (LRP), defined as the point where loads and displacements are referred to, is chosen to be coincident with the caisson centroid unless specified differently; the sign convention adopted in the analyses is also illustrated in Fig. 1a.
The profiles of both overconsolidation ratio OCR and small-strain shear modulus G 0 are plotted in Fig. 1b. The silty clay is assumed to be slightly overconsolidated with an OCR profile consistent with a uniform erosion process (unloading vertical stress Dr v = -340 kPa), reproduced in the analyses by means of the stepwise profile of Fig. 1b. G 0 is assumed to increase with depth according to the empirical relationship proposed by [20] for gravelly sands (assuming a maximum void ratio e max = 0.8, a minimum e min = 0.4 and a relative density D r = 60%) and by [29] for the silty clay (assuming a plasticity index I P = 25%). The earth pressure coefficient at rest is evaluated by using the relationship proposed by [24]. A linear elastic-perfectly plastic behaviour is assumed for the foundation soils with a Mohr-Coulomb failure criterion. The influence of the choice of such a constitutive law on the numerical results is discussed in ''Appendix 1''. The physical and mechanical properties adopted in the analyses are listed in Table 1, where c is the unit weight, m the Poisson's ratio, G/G 0 the ratio between the current and the small-strain shear modulus, c' the cohesion, u' the internal friction and w the dilatancy angles.

Numerical modelling
The numerical analyses have been carried out by using the FE code Plaxis 3D AE [6]. In case of H/D = 0.5 and 1, the numerical model shown in Fig. 1c, d has been used: the 3D mesh consists of 95,500 10-node tetrahedral elements with 4 Gaussian points [12]. In case of H/D = 2, the depth z of the model has been doubled using a mesh of approximately 188,000 elements. Thanks to the problem symmetry with respect to the x-axis, only half domain has been modelled. Lateral boundaries are located at a distance x = y = 6.25 D from caisson axis, where the contours of both soil displacements and stresses have been checked not to be affected by the mesh boundaries (e.g.: [13,14]).  Table 1 Properties of the foundation soils and parameters adopted in the soil constitutive model Horizontal displacements on vertical boundaries, as well as horizontal and vertical displacements at the base, are not allowed.
In the analyses, the caisson is always ''wished in place'', neglecting the simulation of the construction process. A linear-elastic behaviour is assumed for the caisson with a Young's modulus E c = 30 GPa and a Poisson's ratio m = 0.15; the unit weight of reinforced concrete is assumed equal to c c = 25 kN/m 3 . At the soil-caisson interface, purely attritive interface elements with Mohr-Coulomb failure criterion are inserted to simulate relative sliding, with a friction angle d = tan -1 [2/3 tanu']. The choice of such a value of d, although consistent with common practice, does not affect the generality of the numerical results since, as shown by [17,38], the friction angle at the foundation-soil interface seems to have a slight influence on the IDs especially for high values of H/D.

Analyses and results
Two series of 3D FE analyses have been performed to investigate the bearing capacity of the caissons (1) under centred vertical loading and assuming drained conditions, (2) under a general combination of N, Q and M, assuming both drained and undrained conditions for the foundation soils. Specifically, in (2) drainage conditions have been varied only in the calculation phases during which horizontal forces and overturning moments are applied, bearing in mind that such components can be representative of seismic-induced inertial forces acting under undrained conditions. Conversely, drained conditions are always considered for the vertical load, assuming that the excess pore water pressures in the foundation soil, due to the construction of the superstructure, are fully dissipated at the end of the construction phase.

Centred vertical loading
To investigate the bearing capacity under a centred vertical load, three displacement-controlled analyses have been performed. They consist of the following calculation phases: (1) initialisation of the effective stress state; (2) wished-in-place caisson activation; (3) progressive application of a vertical displacement w. For each embedment ratio H/D, the resulting N-w curves are plotted in Fig. 2a. As was expected, the limit load increases with H/D and is attained for very large values of w. Ultimate loads are compared with the values of N lim provided by the code OPTUM G3 [26], where LA calculations are combined with the FE method [32]. The numerical model, characterised by the same dimensions and boundary conditions as the ones described in Sect. 3, comprises the soil and the caisson, the latter modelled as a rigid body. The soil domain is discretised by means of solid elements characterised by a rigid-perfectly plastic behaviour assuming a Mohr-Coulomb failure criterion and an associative flow rule. In each analysis, an adaptive mesh with automatic refinement in proximity of plastic strain localisations has been used. After a few iterations, the analysis stops when the range of the upper and lower bounds to the bearing capacity is small enough to give an accurate estimation of the exact solution.
The drained FE analyses, carried out by assuming a nonassociative flow rule (dilatancy is nil), cannot provide bearing capacities coincident with those obtained from LA calculations, where associativeness is imposed. For this reason, FE results have been compared with LA results computed using both the original strength parameters, c 0 and u 0 , and the reduced values, c 0 red and u 0 red , evaluated as proposed by [10]. Since the role of dilation is more crucial for deeper failure mechanisms, as it was expected, the agreement between FE and LA results, the latter obtained with the original strength parameters, is better for H/ D = 0.5, whereas the opposite is true for H/D = 2 ( Fig. 2a). H/D = 1 may be considered as an intermediate case.
In Fig. 2b, the N-w curves are plotted in the non-dimensional plane N/N lim -w/w el , where w el corresponds to the elastic vertical displacement calculated for N = N lim : therefore, w el = N lim /K 0 , where K 0 is the initial tangent stiffness of the N-w curves. By definition, all normalised curves exhibit an initial linear-elastic response with tangent stiffness equal to 1, followed by a nonlinear load-displacement curve until the ultimate condition (N/N lim = 1) is attained. The nonlinear and irreversible response of caisson foundations to vertical load is strongly influenced by the embedment ratio: as H/D increases, the ''structural hardening'' related to the progressive increase in the plastic zone developing in the soil surrounding the caisson's base and shaft becomes more and more evident. Similarly to what observed in [11] for the unloading occurring at the face of a deep tunnel, the progressive loading of the foundation causes an expansion of the plastic domain: yielding starts from the edges of the foundation's base to deepen and widen as loading increases, with plastic strains spreading over the soil surrounding the caisson. As soon as the mechanism reaches the upper boundary, the horizontal asymptote of the load-displacement curve is attained and failure occurs. In Fig. 3, the progressive evolution of the yielded volume during the unloading at the face of a deep tunnel (Fig. 3a) and that computed for a caisson with H/ D = 2 ( Fig. 3b) is shown. Under the assumption of a sufficiently large ratio of the cover to tunnel diameter, the mechanism ( Fig. 3a) involves the soil surrounding the tunnel face and lining without reaching the ground surface, this resulting in the absence of a plateau in the load-displacements curve [11].
The following analytical expression of the non-dimensional N-w curves is proposed to simulate the numerical results: This represents a modified hyperbole with the following properties: (1) the initial tangent stiffness is equal to 1; (2) the ultimate condition (N/N lim = 1) is reached for a finite value of w/w el = k, where k describes the ''structural hardening'' not governed by the upper boundary and uniquely determined once the constitutive relationship and the constitutive parameters are assigned; (3) the curve is discontinuous in the first derivative for w/w el = k; (4) parameter r = 1.9 is chosen to best-fit the curves obtained from the numerical analyses. The discontinuity of the curve for w/w el = k emphasises the sudden change in stiffness observed when the plastic zone reaches the upper boundary. In Table 2, the values of k evaluated for each embedment ratio are listed. The non-dimensional Nw curves computed using Eq. (1) are compared in Fig. 2b with those obtained numerically. The agreement is more than satisfactory. Bearing in mind that the limit vertical load is attained for high values of the vertical displacement, definitely incompatible with the superstructure operating conditions, Eq. (1) can be employed to introduce an alternative criterion to define the attainment of an ultimate Push-over

Combined N-Q-M loading
To investigate the bearing capacity under a general load combination, about 370 load-controlled pushover numerical analyses have been performed. After the initialisation of effective stresses and the caisson activation, the following calculation phases are completed: (1)  In contrast to what observed for shallow footings [7,25], in case of caisson foundations: (1) when N net = 0 bearing capacity is different from zero; (2) owing to the shaft resistance ID extends also for N net \ 0 (traction); and (3) similarly to what observed also by [17,38], the ID section in the Q-M/D plane shows a nonzero strength even for N net very close to N lim,net ; in any case, this portion of the load space has not been investigated further in this paper.
From the results of the numerical analyses, the collapse of the soil-foundation system under a general N-Q-M load combination (plateau of the pushover curves) may be inferred to be attained for values of displacements varying in a wide range: from centimetres to metres, depending on (1) the loading path (a G ), (2) the embedment ratio (H/D), (3) the vertical load (N net ) and (4) the drainage conditions. The highest values of displacements are attained for the deepest caissons (H/D = 1, 2) subject to high vertical loads and under drained conditions. Therefore, as was suggested by [17,38], a different criterion to define the ultimate condition is adopted here in the following: the attainment by the tangent stiffness of the pushover curve of the 1% of the initial stiffness (K tan /K 0 = 1%) (Fig. 5a). Such a condition leads to the evaluation of limit loads corresponding to computed displacements much smaller than those referred to the plateau of the pushover curves. For H/ D = 1, N net = 7 MN and drained conditions in Fig. 5b, IDs obtained by employing the two different criteria, plateau (solid line) and K tan /K 0 = 1% (solid line with crosses), respectively, are compared. The pushover curve for a G = 0 in Fig. 5a is represented in the non-dimensional generalised force-displacement plane |F|-|u|, defined as: In case the tangent stiffness criterion is adopted, a value |u|= 0.015 is obtained for the case under consideration, whereas |u|= 0.089 is computed at last step of convergence.
Owing to its symmetry with respect to the origin of the M/D-Q plane, half of the envelope is represented in Fig. 5b. IDs obtained by using the different criteria are characterised by the same shape and this is true for Fig. 4b-d as well, where the other sections are plotted. The same may be inferred in case of loci corresponding to assigned values of juj (juj= 0.005, 0.015, 0.025, 0.035, 0.045): an almost homothetic expansion of the envelopes correspond to an increase in juj, progressively less spaced as ultimate condition is approached. This observation, in view of modelling the soil-foundation response to general load combinations by means of a macro-element approach [22,25,28,31,35] based on the theory of elasto-plasticity, seems to suggest the adoption of an isotropic hardening rule in the Q-M plane. Finally, in Fig. 5b the results of FE pushover analyses are compared with those obtained by employing limit analysis and in particular code OPTUM G3, in which strength parameters are reduced as proposed by [10]: again, the agreement is generally good for all loading paths (for any values of a G ) although less satisfactory for high vertical loads. As was previously discussed (Sect. 4.1), the reduction to be applied to strength parameters is however problem dependent.
The model capability to predict the bearing capacity under a general loading combination has also been checked by comparing the IDs obtained for the H/D = 1 caisson subject to N = 0 assuming in turn (1) undrained and (2) drained conditions with published data (Fig. 6). The results are represented in the non-dimensional plane Q/Q u -M/M u , where Q u and M u denote the values of Q and M bringing the soil-caisson systems to collapse when the other load component is zero: Q u = Q (M=0) and M u = M (Q=0) . For case (1), the solutions obtained by [5,17] by means of 2D and 3D total stress analyses are shown, the former assuming a shear strength constant (S u = const.) or linearly increasing with depth (S u = k z), the latter assuming S u = const. only. For case (2), instead, the comparison is made with the ID obtained in [38] by means of 3D numerical analyses where the foundation soil consists of a layer of sand.
In case (1), the IDs are referred to the base of the caisson, chosen by [5] as LRP. The best agreement is obtained with the ID computed in [5] under the assumption of S u = kÁ z, as the same assumption of a soil shear strength increasing with depth is made. For both drainage conditions, the comparison is satisfactory showing that, for the case of a two-layer soil deposit hereby considered, the ID shape is very slightly affected by the assumed soil profile

Parametric study
In this section, the role in affecting the shape and size of IDs of (1) the initial loading factor v = 1/F Sv , where Prior to showing the results of the parametric study, a discussion on how the LRP location affects the shape of failure envelopes is presented. Given that LRP can be chosen arbitrarily, the most frequently adopted location is the top [17,38], the centroid [23] and the base of the foundation [5]. The effect of LRP location is presented in Fig. 7a Fig. 7a), the envelope is almost perfectly symmetric and crosses the x and y-axes perpendicularly: this is due to the decoupling between the rotational and horizontal degrees of freedom. In view of developing a macro-element model for caisson foundations, it is evident that both symmetry and decoupling are desirable for an easy analytical representation of ID. Decoupling is however lost for higher values of embedment (Fig. 7b). To highlight this item, vector plots of some failure mechanisms computed for the two caissons are shown in Figs. 8 and 9 together with the relevant contours of the deviatoric strain. IDs, relative to the caissons centroid, are represented in the nondimensional plane Q/Q u -M/M u to focus on their shape. In Fig. 8, the caisson with H/D = 0.5, subject to M only (a G = !, point A), undergoes an almost pure clockwise rotational mechanism around a point close to the centroid (''scoop'' mechanism, according to [38]). Similarly, when subject to Q only (a G = 0, point C), the caisson undergoes a pure translational mechanism with deviatoric strains developing in front and at the back of the caisson, where pseudo-active and pseudo-passive wedges are detected.
Between points A and B, a coupled (''scoop-slide'') mechanism is attained with the caisson rotation pole moving downwards from the centroid to infinite (in point C). Between points B and D, a translational mechanism is observed, while after point D the mechanism is coupled again, with an anticlockwise rotation of the caisson (''reverse-scoop'') and the rotation pole progressively approaching the caisson centroid again, moving from above as Q decreases (points E and F). Similar failure mechanisms are observed for the deeper caisson (H/D = 2, Fig. 9). However, a rotation about a point deeper than the centroid is observed in A and a coupled mechanism in C: both symmetry and orthogonality of the failure envelope with respect to both x and y-axes are lost. The way the coupling effect of sliding and rotation affects the depth of the rotation point has also been observed in [3,38]. This condition becomes more and more evident as H/D increases, as the lever arm of the resultant force transferred to the caisson centroid increases with H. Indeed, for the caisson with H/D = 2, the pure translational mechanism is attained when the horizontal force is applied to a point deeper than the centroid, at a depth z = H/2-a G = 16 m = 2/3ÁH (point D in Fig. 9).

Influence of the initial loading factor
For increasing vertical loads (v = 0.63, F sv = 1.6), symmetry and orthogonality of ID are lost even for H/D = 0.5 (Fig. 10). Higher values of vertical load induce a deepening of failure mechanisms, as shown by the comparison between the contour plots plotted in Figs. 8 and 10, with a more pronounced asymmetry in the Q-M plane, this resulting in a different shape of the envelope, especially To better appreciate the influence of the vertical load on both ID size and shape, in Fig. 11 these are plotted for caissons with H/D = 0.5 and 1, for undrained applications of Q and M, both in the dimensional Q-M/D and nondimensional Q/Q u -M/M u planes. Owing to the progressive change in the failure mechanisms geometry, which become deeper and more asymmetric (Figs. 8, 10), at increasing values of v failure envelopes in Fig. 11a, c, first harden and subsequently shrink, in accordance with [38]. Specifically, at increasing v values up to a certain threshold, the bearing capacity increases since larger volumes of soil are involved as mechanisms deepen. However, above a threshold value of v, bearing capacity reduces since the asymmetry of the mechanisms seems to prevail, with a contraction of the failure envelopes.
Comparison of Fig. 11b, d allows us to appreciate the influence of v on the shape of IDs: this is more pronounced for the shallower caisson (H/D = 0.5). A significant change

Influence of caisson embedment ratio
While the influence of embedment ratio H/D on IDs size has been already discussed (Sect. 4.3), here the focus is on the role of H/D in affecting the shape of IDs (Fig. 12) in case undrained conditions are accounted for. Consistently with the results shown by [17], a significant change in the shape of IDs is observed for small vertical loads, v B 0.21, as H/D increases (Fig. 12a, b): the envelope eccentricity progressively increases in the IV quadrant of the non-dimensional plane. By contrast, the influence of H/D vanishes for higher vertical loads (Fig. 12c, d).

Influence of drainage conditions
The influence of drainage conditions, on both size and shape of IDs, is illustrated in Fig. 13, where IDs obtained imposing both undrained and drained conditions are compared for v = 0.09 and different embedment ratios. As was expected for compression loading paths, undrained soil behaviour always results into smaller envelopes with respect to those computed under drained conditions (Fig. 13a, c, e). Indeed, due to the excess pore water pressure, foundation soils undergo an effective stress and in turn a shear strength reduction. Conversely, drainage conditions do not affect the shape of the failure envelopes (Fig. 13b, d, f). This observation may lead to assuming a failure envelope hardening homothetically as drainage occurs with time, until reaching the drained envelope.

Proposed relationships for interaction diagrams
An analytical expression for IDs is proposed to fit all the conditions investigated in this study. The general equation, valid for every initial loading factor v, embedment ratio H/ D and for both undrained and drained conditions, is that of a unit circle in the n 0 -l 0 plane: By scaling Eq. (3), the equation of an ellipse is obtained in the form n 0 a n where a n and a l , defined in the following, are a function of v, H/D and drainage conditions. Equation (4) is then mapped into the n-l plane by applying a rotation of an angle x (positive if counterclockwise) where n and l, following [7], are defined as Hence, the proposed equation in the n-l plane (Eq. 7) is that of a rotated ellipse having semi-axes, a n and a l , and rotation angle, x, depending on v, H/D and drainage conditions: The numerical IDs corresponding to the five values of v under consideration have been first interpolated by ellipses to obtain the values of x, a n and a l plotted by symbols in Fig. 14, for each H/D ratio and for both drainage conditions. In turn, these values have been best-fitted to define the following functions of v: and a n v ð Þ ¼ c 21 ðv À c min Þ c 22 ðc max À vÞ where the non-dimensional parameters c i are listed in Table 3 as a function of H/D and the drainage conditions. Figure 14 shows a satisfactory agreement between the proposed Eqs. (8)(9) and the values of x, a n and a l represented by symbols. The n-l plane locus is obtained from the analytical expression above according to the method illustrated in Fig. 15. Specifically, with reference to undrained conditions and H/D = 0.5, the unit circle (Fig. 15a) is first distorted to the ellipse (Fig. 15b), by applying Eq. (4), and then rotated by the values of x provided by Eq. (8) (Fig. 15c).
In Fig. 16a, the three-dimensional ID under drained conditions provided by Eq. analytical solution proposed by [7] for shallow footings resting on sand is also presented (dotted line). It is evident that: (1) IDs increase in size with H/D, and (2) in contrast to the case H/D = 0, the cross sections of the caissons IDs in planes n = 0 and l = 0 (Fig. 16c, d) do not close neither for v = 0 nor for v = 1 (Sect. 4.2). Conversely, in the solution proposed by [7], an almost perfectly symmetric to v = 0.5 parabola is suggested. Finally, in planes v = const., the solution proposed for shallow footings is an ellipse rotated by a constant angle x = -14°, while for embedded caisson foundations the rotation angle x is found to vary with v, especially for low values of H/D. The IDs provided by Eq. (7) for undrained conditions and each embedment ratio are plotted in Fig. 17. Similarly to what observed for drained conditions, the tendency of failure domains not to close for v = 1 becomes more evident as H/D increases: K 0 denotes the elastic generalised displacement corresponding to F lim and K 0 is the initial tangent stiffness of the F À u j j curves. As it is evident in Fig. 18, for deep caissons, in contrast to what is observed in case of low values of H/D (B 1), the dependency of LDCs on a G is less significant, being the scatter not marked. The same observation does apply for drained conditions, although the representation is here omitted for the sake of brevity. In view of defining a pre-design simplified approach, the following equations of LDCs are proposed for both drainage conditions: These simple expressions, neglecting the dependency of LDCs on parameters v, H/D, a G , are justified by the range of values of u j j u j j el considered for design purposes, far from F . F lim ¼ 1. In preliminary design computations, the user can evaluate F lim from the rotated ellipse discussed in Sect. 5 (Eq. 7) and K 0 from the elastic solutions proposed by [15,21,34] (see ''Appendix 2''), as a function of input parameters such as v, H/D, a G and drainage conditions. Then, the threshold displacement ratio u j j u j j el needs to be selected only in Eqs. (10).  Table 4, for three values of v and both drainage conditions, correspondingly to u j j lim ¼ 1&. As it is shown in Fig. 18, for undrained conditions, in the range of values of u j j u j j el here considered for design purposes, u j j u j j el 1, the analytical expressions [Eqs. (10)] fairly match the curves obtained from the numerical analyses. It is worth mentioning that one can also compute the

Conclusions
Rigid and massive caisson foundations are typically subject to combined loading conditions, thus withstanding vertical and horizontal forces at the same time, as well as moment. For this reason, interaction domains may be a useful tool to assess their safety against limit states when designing these foundations, as the distance of the image generalised stress point from the boundary of the interaction domain may be interpreted as a sort of safety factor. However, as full mobilisation of the caissons capacity under general loading conditions tends to be achieved for large displacements not compatible with working conditions, a displacement-based approach for assessing the safety factor at a preliminary design stage may be preferred. This approach has been followed in this paper, based on the non-dimensional expressions provided for both the IDs and generalised load-displacement curves obtained by performing FE analyses. Indeed, the response of massive cylindrical onshore caisson foundations in a two-layer soil, subject to general loading conditions, has been investigated by means of a series of 3D elastic-plastic FE numerical analyses.
Bearing capacity under a centred vertical load is first investigated by assuming a drained response for the foundation soil and three different values for the embedment ratio. Numerical results testify that load-displacement relationships are severely affected by the caisson embedment ratio: indeed, structural hardening related to spatial propagation of the plastic zone, from the caisson base to the upper horizontal boundary of the FE mesh, becomes more and more pronounced as H/D increases.
Ultimate response of the caissons under a general combination of vertical and horizontal loads (N, Q), as well as overturning moment (M), is also investigated. The 3D envelopes in the N-Q-M/D space present a rugby-ball shape, similar to that obtained for shallow footings in previous works. However, soil-caisson strength under Q-M combinations is different from zero when the vertical load is equal to zero and close to its limit value (N lim,net ), differently from shallow foundations.
From the comparison of the results obtained in this study (two-layer soil) with those from previous ones (homogeneous soil), it follows that IDs shape is very slightly affected by the assumed soil layering, whereas it is influenced by the profile of shear strength (constant or linearly increasing with depth).
A parametric study has been carried out to understand the factors mainly affecting both size and shape of failure envelopes in the N-Q-M space. The influence of initial loading factor (v) and caisson embedment ratio (H/D) has been evaluated for both fully undrained and drained conditions, showing that the assumed drainage condition scales up the ID size rather than modifying its shape. The influence of the load reference point location is also discussed, with the centroid of the caisson being usually the most suitable choice, as this strongly simplifies the shape of ID, leading to decoupling of rotational and horizontal degrees of freedom for low embedment ratios (H/D = 0.5).
Analytical expressions best-fitting the three-dimensional interaction diagrams computed in the non-dimensional N-Q-M/D load space, together with the generalised loaddisplacement curves are also proposed, for both undrained and drained conditions: these are integrated according to a simplified approach proposed as a useful tool at a preliminary design stage to ensure the foundation displacements to be compatible with the desired structural performance, for both serviceability and ultimate limit states. Future work will be oriented towards integrating the obtained results in an elastic-plastic isotropic hardening macroelement model, where the ID equation can be used to define the yield function.

Appendix 1
The influence of the soil constitutive law on the IDs is shown in Fig. 19. For the H/D = 1 caisson subject to v = 0.21, the IDs have been obtained from two sets of purely undrained numerical analyses where the following models have been used to describe the soil behaviour: (1) linear elastic-perfectly plastic and (2) nonlinear elasticplastic model with isotropic hardening (Hardening Soil with small-strain stiffness, [1]), both assuming a Mohr-Coulomb failure criterion. Table 5 summarises the parameters assumed for HS Small: details about the meaning and choice of such parameters can be found in [13]. The choice of different constitutive laws affects the evolution of the soil stiffness and, consequently, the excess pore pressure arising under undrained conditions. This affects the evolution of the generalised force-displacement curves attaining a different plateau, as shown in Fig. 19b for a G = 0. However, Fig. 19a shows that the choice of the constitutive law does not affect the shape of the IDs and that the overestimation of the Q-M strength does not exceed the 10%. Under drained conditions, instead, nor the shape neither the dimensions of the IDs are affected since the plateau attained by the curves is the same as a consequence of the same failure criterion assumed by the two models (Fig. 20). If a criterion different from the plateau is used to obtain the IDs, as those based on the K tan /K 0 ratio or on the generalised displacement |u| described in Sect. 4.2, its influence on the IDs shape and dimensions is negligible though the constitutive law affects the evolution of the generalised force-displacement curves. In view of the above and of the purpose of the paper, a less complex constitutive law has been adopted in the parametric numerical study, thus taking advantage of a lower computational effort.

Appendix 2
To calculate u j j el ¼ F lim . K 0 , stiffness K 0 is to be evaluated as follows (a) (b) Fig. 19 Influence of the constitutive law under undrained conditions: ID for the H/D = 1 caisson subject to v = 0.21 obtained assuming the elastic-perfectly plastic (black) and the HS Small (grey) model for the soil (a); pushover curve for the load path a G = 0 (b) where ||Á|| stands for the matrix norm. In Eq. (11) C represents the elastic compliance matrix of the soil-caisson system calculated as the inverse of the stiffness matrix The translational, rotational and coupled stiffness terms in Eq. (12) can be evaluated according to [15,21] for H/ D B 1 and to [34] for H/D = 2. The ± sign in Eq. (11) is introduced to properly exploit the previously cited elastic solutions since they assume the caisson base (?) as LRP for H/D B 1 and the caisson top (-) for H/D = 2, respectively.
Authors' contributions Not applicable.
Funding Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement.
Availability of data and materials Data will be made available on request.
Code availability Not applicable.

Declaration
Conflict of interest The authors have no conflict of interest do declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.