Practical Soil-Shallow Foundation Model for Nonlinear Structural Analysis

Soil-shallow foundation interaction models that are incorporated into most structural analysis programs generally lack accuracy and efficiency or neglect some aspects of foundation behavior. For instance, soil-shallow foundation systems have been observed to show both small and large loops under increasing amplitude load reversals. This paper presents a practical macroelement model for soil-shallow foundation system and its stability under simultaneous horizontal and vertical loads. The model comprises three spring elements: nonlinear horizontal, nonlinear rotational, and linear vertical springs. The proposed macroelement model was verified using experimental test results from large-scale model foundations subjected to small and large cyclic loading cases.


Introduction
Several researchers ( [1][2][3][4][5], among others) have investigated extensively the subject of soil-structure interaction (SSI). Rayhani et al. [6] showed that soil-structure interaction might amplify or attenuate the base shear through inertial and kinematic interactions. Stewart et al. [7] concluded that the effects of SSI on rigid structures founded on soil are more significant, when compared to flexible structures. Bobet et al. [8] showed that the soil-structure system is a function of the relative rigidity of the structure compared to that of the ground. SSI is shown to be significant in the presence of soft soils or when structural mass is very large [9]. As a consequence, the SSI problem can be excluded from computations if the soil where the structure is founded is very rigid. When considering the effect of soil-structure interaction in base isolated multistoried structures on elastic layered soil, Spyrakos et al. [10] have found that SSI effects are significant for squat lightweight buildings on low stiffness soil-stratum.
Although their study dealt with harmonic excitations, it still gives an insight on the danger of neglecting SSI in the design of base isolated buildings.
In the seismic resistant design of structures, we are most interested in the strength reduction factors to account for the nonlinear behavior that might be experienced by a structure subjected to an earthquake ground motion. Few researchers [11,12] have recently attempted to assess the effect of SSI on the strength reduction factors. Eser et al. [11] have shown that the presence of soft soils reduces the strength reduction factors, which is primarily controlled by the changes in the structural period and displacement ductility.
Incorporation of SSI requires explicit modelling of soilfoundation system adequately. For instance, several models have been proposed depending on the foundation type, its embedment, and its rigidity ( [13][14][15], among others). The Federal Emergency Management Agency [16] required that the foundation stiffness should be determined with one of the following three methods: (i) uncoupled spring model comprising three spring elements, for shallow foundations that are stiffer than the supporting soil; (ii) a finite element formulation of linear (or nonlinear) foundation behavior using Winkler models, for shallow foundations that are less stiff than the supporting soil; (iii) decoupled Winkler model, for shallow foundations that are flexible with respect to the supporting soil.
However, El Ganainy and El Naggar [17] have demonstrated that the decoupled technique (Beam on a Nonlinear Winkler Foundation, BNWF) is not capable of predicting accurately the settlement seen in foundations on soft soils. Although it can be used to predict the overall deformation behavior of foundations, the BNWF requires a large number of nonlinear springs, which is considered as a major drawback [17]. To address some of the above-mentioned issues, macroelement formulations have been proposed. The first formulation has been developed by Nova and Montrasio [18] and later modified and/or extended by other researchers [19][20][21][22] and recently the formulations developed by Gajan et al. [23] and Shirato et al. [3]. The major advantages of such models are their simplicity and their ability to capture the global response of bearing foundations [22,24,25]. However, on one hand, calibrating macroelement parameters pose a constraint on their adoption for practical applications. And on the other hand, knowing that most of the available macroelement models are based on specified bounding surfaces poses another problem for their capability to cover a wide range of problems [17].
In a completely different modelling approach, El Shamy and Zamani [26] proposed a new 3D particle-based technique using the discrete element method (DEM) to analyze the seismic performance of soil-foundation-structure systems. In their model, the soil is idealized as a collection of spherical particles using DEM; the footing is considered as a rigid block, whereas the structure is modelled using a number of spherical particles in the form of a column, which can be clamped to simulate a rigid structure or bonded to simulate a flexible structure of predefined rigidity.
To overcome the difficulties in performing complete nonlinear simulations, Seylabi et al. [27] proposed an equivalent linearization of nonlinear soil-structure systems considering both the effect of SSI and the nonlinear behavior of the structure on equivalent linear parameters. In their model the structure is modelled as an elastoplastic single-degree-offreedom system (SDOF) and the soil beneath the structure is modelled by a discrete model combining different spring and dashpot elements.
In this paper, a new macroelement model is developed for the analysis of the nonlinear response of shallow foundations under cyclic loading. This model may easily be incorporated into available structural analysis programs such as OpenSees [28]. The soil-foundation system is simulated using three spring elements: horizontal and rotational nonlinear springs and a linear vertical spring. The nonlinear springs are assigned appropriate nonlinear model of plasticity with material degrading parameters.

Proposed Macroelement Model
The problem being studied here is that of a shallow foundation of any shape embedded in soil and subjected to simultaneous axial and lateral forces, as shown in Figure 1. The foundation is considered to be very stiff. The depth of embedment is . The proposed model incorporates three types of springs: (i) vertical translational elastic spring with stiffness V ; (ii) shear inelastic spring with preyield stiffness ; (iii) rotational inelastic spring with preyield stiffness .
These equivalent springs represent the foundation-soil system. The macroelement model replaces the system soilshallow foundation, thus decreasing considerably both the overall number of degrees of freedom and the computation effort required to run large models.

Constitutive Equations.
Two material models are considered in this study, namely, the Bouc-Wen model [29,30] for the shear and rotational springs and the linear model for the vertical spring.

Shear and Rotational Springs
Model Assumptions. The smooth Bouc-Wen model of hysteresis by Bouc [29] and Wen [30] has found many engineering applications. For instance, the use of original Bouc-Wen model and its extensions to soil-structure interaction include [31][32][33].
In this study, we considered the Baber and Noori [34] extension to the original Bouc-Wen model. This version includes the degrading behavior observed in many engineering materials.
The nonlinear behavior of the soil-shallow foundation system is modelled via the nonlinear shear and nonlinear rotational springs. A force ( ) is mobilized at the shear spring and a moment ( ) is mobilized at the rotational spring. In what follows, we use a general force ( ) to denote ( ) and ( ), depending on which spring we are modelling. The constitutive relationship for is expressed in the Bouc-Wen model as a linear part and a hysteretic part: where is the deformation; is the preyield stiffness; is the mobilized force at the beginning of yielding; is the post-to-preyield stiffness ratio; and is a hysteretic quantity controlling the nonlinear behavior. The latter is governed by the following differential equation with respect to time : Equation (2) represents the Baber and Noori [34] formulation of the rate of hysteretic deformation, which accommodates degradation. In (2) , , and are parameters that control the shape and transition from elastic to inelastic regions of the hysteretic loop, while , ], and are variables that control the stiffness degradation and material deterioration. Note that the yield deformation of the spring ( ) does not appear in (2), unlike the original Bouc-Wen model. In fact, the adoption of the Baber-Noori [34] formulation in calculating forces is easy by considering that and are expressed as in which and can take values as in the original Bouc-Wen model of hysteresis. In this case, is a dimensional quantity and is bounded between − and + . The evolution of material degradation is governed by the following equations [34]: where is defined by the following rate equation: And , ] , and are the parameters that control the degradation.
Incremental Response Equations. The tangent of ( ) (i.e., ( )/ ) is computed using (1): Note that this is the continuum tangent and not the algorithmically consistent tangent. It is clear from (6) that the derivative of with respect to is needed. To get it, (2) can be rewritten aṡ Hence, From this point, we will derive incremental response equations to obtain computer implementable equations. In what follows, the same procedure developed by Haukaas and Der Kiureghian [35] is used. The force at time +1 is obtained as The variable is next discretized by a Backward-Euler solution scheme: It is seen that Δ cancels from the equation, yielding a nonlinear equation in +1 . A Newton scheme of the form +1 = − ( )/ ( ) to solve the general nonlinear equation of the form ( ) = 0 is employed to solve for +1 in (10) [35].
The equations describing the degradation behavior are described as follows: where +1 is found by discretization of the rate equation in (5) using the Backward-Euler solution scheme: where Δ again cancels. Note that , , and are history variables that must be stored at each converged step.
For the complete implementable procedure to obtain +1 and the algorithmically consistent tangent +1 / +1 (referred to here as the Local Newton Routine), the reader is referred to [35].

Vertical
Spring. The vertical force V mobilized in the vertical translational spring is given by in which V is the elastic vertical stiffness of the soilfoundation system and is the vertical displacement of the spring. Note that an appropriate constitutive model could be assigned to the vertical spring; however, in this paper, the spring is considered linear elastic, but the effect of axial loads is present on the soil-shallow foundation global behavior.

Equilibrium Equations and Model
Implementation. The macroelement model described in the previous subsections is used to simulate the soil-shallow foundation interaction. Initially, we assume that the shear and rotational springs are linear; then we replace them by the general forces ( ) and ( ), respectively. The total potential energy Π of the model with reference to Figure 2 is expressed as follows: where and are the foundation's top horizontal and vertical displacements, respectively, and are dependent on the depth of embedment and the spring's displacements: = (1 − cos ( )) + sin ( ) + V cos ( ) .
In (15b), , positive in case of compression, is the sum of the vertical component of the displacement originating from the axial flexibility of the soil and the additional vertical displacement that happens in the laterally deformed configuration shown in Figure 2.
The above kinematic equations (see (15a), (15b)), which relate the total displacement and the vertical displacement to the internal displacements , V, and , assume large displacements and large rotations.
Considering the model subjected to the axial load with the resulting horizontal displacement , the vector of unknowns is then x = ⟨ , V, , ⟩ . Using Castigliano's second theory by imposing the stationary of Π with respect to , V, and yields the following governing equations: The following equation should also be considered to permit the displacement-controlled cyclic analysis by imposing the top lateral displacement : However, when the lateral force-controlled cyclic response is desired, is to be imposed. Hence, the system of governing equations to be solved is reduced to the three equilibrium equations (see (16a), (16b), and (16c)). The lateral and vertical displacements ( , ) may be calculated at the end of the analysis using the kinematic equations (see (15a), (15b)).
Considering the case of lateral displacement-controlled analysis, (16a), (16b), (16c), and (16d) are rewritten in the following form: Mathematical Problems in Engineering 5 Note that the resultant moment ( = ( − ) + ) found by the equilibrium of the macroelement in the deformed configuration needs also to be applied at the top of the foundation. Also note that in (17), the terms and were replaced by the general forces ( ) and ( ), respectively. These forces are assigned the constitutive model ( ) developed in Section 2.1.1.
The system g(x) can be solved using Newton's method following the pseudocode provided below (referred to here as the Global Newton Routine): Previous converged solution x = ⟨ , V, , ⟩ From the current estimation of the top lateral displacement and axial load , the code returns the vector x = ⟨ , V, , ⟩ that satisfies equilibrium and kinematics (see (17)). In each global Newton iteration the current values of , , / , and / based on the current values of and are updated using the Local Newton Routine (Section 2.1.1). The above pseudocode (Global Newton Routine) along with the Local Newton Routine (that implements the Bouc-Wen-Baber-Noori model of hysteresis) has been implemented numerically in MATLAB (Mathworks, Inc.).

Verification with Experimental Results
To verify the validity of the proposed macroelement model in predicting the cyclic behavior of the soil-shallow foundation system, its predictions are compared with experimental results. In the framework of the TRISEE Project (3D Site Effects and Soil-Foundation Interaction in Earthquake and Vibration Risk Evaluation) a program of large-scale 1 g model has been tested to investigate the response of soil-shallow foundation under cyclic and dynamic loads. The experiment was carried out at ELSA (European Laboratory for Structural Assessment) in Ispra, Italy; test results are reported in many references ( [36], among others).
The TRISEE experiment consists of three phases; only Phases I and III are considered in the comparison because the proposed macroelement model is restricted to 2D loading conditions only. However, the model may be extended to include the 3D case.
The experimental setup consists of a square steel shallow foundation (1 m × 1 m) mounted inside a stiff concrete caisson (4.6 m × 4.6 m × 3 m) filled with Ticino sand. The embedment of the foundation was about 1 m, with a steel framework placed around the foundation to retain the sand. Two different soil relative densities have been used, = 85% and = 45%, representing high (HD) and low density (LD) soil conditions, respectively.
The HD and LD specimens were loaded vertically by 300 kN and 100 kN, respectively. The static safety factor was about 5 in both tests. In the HD test, the load was applied at 0.9 m above the foundation and in the LD test at 0.935 m. In Phase I, a series of unidirectional force-controlled small amplitude cycles was applied to identify the significance of nonlinear soil behavior. In Phase III, displacementcontrolled, unidirectional, increasing amplitude cycles were imposed to the top of the foundation. For further information on experimental setup and results refer to [36].
The new macroelement model is used to simulate the foundation. The parameters of the numerical model are presented in Table 1. These parameters have been calibrated using experimental moment-rotation and horizontal forcehorizontal displacement curves. Parameters for the Bouc-Wen-Baber-Noori model of hysteresis are given in Tables 2  and 3, for HD and LD tests, respectively. The parameters were again calibrated with experimental results. Figures 3-6 compare the experimental and the numerical hysteretic horizontal force-horizontal displacement curves and moment-rotation curves for Phases I and III. The numerical model reproduces correctly the overall behavior of the foundation, verifying the ability of the proposed macroelement model to simulate the cyclic behavior of shallow foundations.   The uplift is important for the HD Phase III test, and this can be observed from the S shaped moment-rotation curve (see Figure 5). For the LD sand, only plasticity is developed and the uplift is not present.

Foundation Stiffness Matrix
In the previous section, large rotations have been considered in formulating the macroelement model. However, assuming that small rotations permits the construction of the foundation stiffness matrix. In this case, (17) can be rewritten as Ryan et al. [37,38] studied the stability of bearing isolators and constructed the lead-rubber bearing stiffness matrix relating a change in the isolator forces to the change in displacements. Following the same idea, the foundation matrix k , relating a change in foundation forces f = ⟨ , ⟩ to the change in displacements U = ⟨ , ⟩ , is derived in three steps: (1) Differentiate the equations of equilibrium (in (18): 1 , 2 and 3 ) resulting in in which V = ⟨ , , ⟩ and the matrices k eq and Γ are given by (2) Differentiate the equations of kinematics (see (15a) and (15b)): (3) Substitute V from (19) into (21): where f eq = k −1 eq . The resultant flexibility matrix, f = Γ f eq Γ, relates the displacement increment U to the force increment F and is given by ] .
Then foundation stiffness matrix k is the inverse of the flexibility matrix, f . The stiffness matrix is computed after calculation of displacements and computation of forces to satisfy the governing equations (see (17)). In a structural analysis program, the routine to compute the foundation stiffness matrix may be incorporated into existing structural analysis software and used for analysis of structures on shallow foundations.

Conclusions
A practical macroelement model is presented in this paper to simulate the response of soil-shallow foundation systems. The proposed model was verified against experimental test results of large-scale model foundations subjected to small and large loading cycles. A summary of the main points presented in this paper is given below: Ganainy and El Naggar [17] can be included by proper choice of degradation parameters of the Baber-Noori model. (4) The proposed model does not take into consideration full coupling between the different springs. Nevertheless, the model represents a first step for future improvements.