# Practical Soil-Shallow Foundation Model for Nonlinear Structural Analysis

Authors: Moussa Leblouba; Salah Al Toubat; Muhammad Ekhlasur Rahman; Omer Mugheida
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

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.


---

## Body

## 1. Introduction

Several researchers ([1–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 soil-foundation system adequately. For instance, several models have been proposed depending on the foundation type, its embedment, and its rigidity ([13–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–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-of-freedom 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.


## 2. 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 Figure1. The foundation is considered to be very stiff. The depth of embedment is H. The proposed model incorporates three types of springs:(i)
vertical translational elastic spring with stiffnessk
v;
(ii)
shear inelastic spring with preyield stiffnessk
s;
(iii)
rotational inelastic spring with preyield stiffnessk
θ.These equivalent springs represent the foundation-soil system. The macroelement model replaces the system soil-shallow foundation, thus decreasing considerably both the overall number of degrees of freedom and the computation effort required to run large models.Figure 1
Proposed macroelement model for soil-shallow foundation system.


### 2.1. 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.


#### 2.1.1. 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–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 forcef
s
(
s
) is mobilized at the shear spring and a moment f
θ
(
θ
) is mobilized at the rotational spring. In what follows, we use a general force f
(
x
) to denote f
s
(
s
) and f
θ
(
θ
), depending on which spring we are modelling.The constitutive relationship forf is expressed in the Bouc-Wen model as a linear part and a hysteretic part:(1)
f
=
α
k
x
+
1
-
α
f
y
z
,where x is the deformation; k is the preyield stiffness; f
y is the mobilized force at the beginning of yielding; α is the post-to-preyield stiffness ratio; and z is a hysteretic quantity controlling the nonlinear behavior. The latter is governed by the following differential equation with respect to time t:(2)
z
˙
=
A
x
˙
-
β
x
˙
z
z
n
-
1
+
γ
x
˙
z
η
ν
η
.Equation (2) represents the Baber and Noori [34] formulation of the rate of hysteretic deformation, which accommodates degradation. In (2) β, γ, and n are parameters that control the shape and transition from elastic to inelastic regions of the hysteretic loop, while A, ν, and η are variables that control the stiffness degradation and material deterioration. Note that the yield deformation of the spring (
x
y
) 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(3)
β
=
β
′
x
y
η
,
γ
=
γ
′
x
y
ηin which β
′ and γ
′ can take values as in the original Bouc-Wen model of hysteresis. In this case, z is a dimensional quantity and is bounded between -
x
y and +
x
y. The evolution of material degradation is governed by the following equations [34]:(4a)
A
=
A
0
-
δ
A
e
,
(4b)
ν
=
1
-
δ
ν
e
,
(4c)
η
=
1
+
δ
η
e
,where e is defined by the following rate equation:(5)
e
˙
=
1
-
α
f
y
x
˙
z
.And δ
A, δ
ν, and δ
η are the parameters that control the degradation.Incremental Response Equations. The tangent of f
(
x
) (i.e., ∂
f
(
x
)
/
∂
x) is computed using (1):(6)
∂
f
x
∂
x
=
α
k
+
1
-
α
f
y
∂
z
∂
x
.Note that this is the continuum tangent and not the algorithmically consistent tangent. It is clear from (6) that the derivative of z with respect to x is needed. To get it, (2) can be rewritten as(7)
z
˙
=
A
x
˙
-
x
˙
z
n
β
sgn
⁡
x
˙
z
+
γ
ν
η
=
∂
z
∂
x
∂
x
∂
t
.Hence,(8)
∂
z
∂
x
=
A
-
z
η
β
sgn
⁡
x
˙
z
+
γ
ν
η
.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 t
i
+
1 is obtained as(9)
f
i
+
1
=
α
k
x
i
+
1
+
1
-
α
f
y
z
i
+
1
.The variable z is next discretized by a Backward-Euler solution scheme:(10)
z
i
+
1
=
z
i
+
1
+
Δ
t
A
i
+
1
-
z
i
+
1
n
β
sgn
⁡
x
i
+
1
-
x
i
/
Δ
t
z
i
+
1
+
γ
ν
i
+
1
η
i
+
1
x
i
+
1
-
x
i
Δ
t
.It is seen that Δ
t cancels from the equation, yielding a nonlinear equation in z
i
+
1. A Newton scheme of the form x
i
+
1
=
x
i
-
f
x
m
/
f
′
x
m to solve the general nonlinear equation of the form f
x
=
0 is employed to solve for z
i
+
1 in (10) [35].The equations describing the degradation behavior are described as follows:(11a)
A
i
+
1
=
A
0
-
δ
A
e
i
+
1
,
(11b)
ν
i
+
1
=
1
-
δ
ν
e
i
+
1
,
(11c)
η
i
+
1
=
1
+
δ
η
e
i
+
1
,where e
i
+
1 is found by discretization of the rate equation in (5) using the Backward-Euler solution scheme:(12)
e
i
+
1
=
e
i
+
1
+
Δ
t
1
-
α
f
y
x
i
+
1
-
x
i
Δ
t
z
i
+
1
,where Δ
t again cancels. Note that z
i, x
i, and e
i are history variables that must be stored at each converged step.For the complete implementable procedure to obtainf
i
+
1 and the algorithmically consistent tangent ∂
f
i
+
1
/
∂
x
i
+
1 (referred to here as the Local Newton Routine), the reader is referred to [35].


#### 2.1.2. Vertical Spring

The vertical forcef
v mobilized in the vertical translational spring is given by(13)
f
v
=
k
v
Vin which k
v is the elastic vertical stiffness of the soil-foundation system and V 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.


### 2.2. 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 forcesf
s
(
s
) and f
θ
(
θ
), respectively. The total potential energy Π of the model with reference to Figure 2 is expressed as follows:(14)
Π
=
1
2
k
θ
θ
2
+
1
2
k
s
s
2
+
1
2
k
v
v
2
-
P
V
-
F
U
,where U and V are the foundation’s top horizontal and vertical displacements, respectively, and are dependent on the depth of embedment H and the spring’s displacements:(15a)
U
=
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
,
(15b)
V
=
H
1
-
cos
⁡
θ
+
s
sin
⁡
θ
+
v
cos
⁡
θ
.In (15b), V, 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.Figure 2
Undeformed and deformed configurations of the proposed macroelement model.The above kinematic equations (see (15a), (15b)), which relate the total displacement U and the vertical displacement V to the internal displacements s,  v, and θ, assume large displacements and large rotations.Considering the model subjected to the axial loadP with the resulting horizontal displacement U, the vector of unknowns is then x
=
s
,
v
,
θ
,
F
T. Using Castigliano’s second theory by imposing the stationary of Π with respect to s, v, and θ yields the following governing equations:(16a)
g
1
=
k
s
s
-
F
cos
⁡
θ
-
P
cos
⁡
θ
=
0
,
(16b)
g
2
=
k
v
v
+
F
sin
⁡
θ
-
P
cos
⁡
θ
=
0
,
(16c)
g
3
=
k
θ
θ
-
P
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
-
F
H
-
v
cos
⁡
θ
-
s
sin
⁡
θ
=
0
.The following equation should also be considered to permit the displacement-controlled cyclic analysis by imposing the top lateral displacement U:(16d)
g
4
=
U
-
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
=
0
.However, when the lateral force-controlled cyclic response is desired, F 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 (
U
,
V
) 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:(17)
g
=
g
1
g
2
g
3
g
4
=
f
s
s
-
F
cos
⁡
θ
-
P
sin
⁡
θ
k
v
v
+
F
sin
⁡
θ
-
P
cos
⁡
θ
f
θ
θ
-
P
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
-
F
H
-
v
cos
⁡
θ
-
s
sin
⁡
θ
U
-
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
.Note that the resultant moment M  (
M
=
F
H
-
V
+
P
U
) 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 k
s
s and k
θ
θ were replaced by the general forces f
s
(
s
) and f
θ
(
θ
), respectively. These forces are assigned the constitutive model f
(
x
) developed in Section 2.1.1.The systemg
(
x
) can be solved using Newton’s method following the pseudocode provided below (referred to here as the Global Newton Routine):Previous converged solutionx
=
s
,
v
,
θ
,
F
T
(1)
While (g
(
x
)
<
t
o
l
e
r
a
n
c
e)
(2)
Updatef
s, f
θ, ∂
f
s
/
∂
s, and ∂
f
θ
/
∂
θ using the Local Newton Routine
(3)
Computeg
(
x
)
(4)
Compute the Jacobian matrixJ
(
x
)
(5)
d
x
=
J
x
-
1
g
(
x
)
(6)
x
=
x
+
d
x
endwhere the Jacobian matrixJ
(
x
) is defined as J
i
j
=
∂
g
i
/
∂
x
j. The tangents ∂
f
s
/
∂
s and ∂
f
θ
/
∂
θ required in the Jacobian matrix should be algorithmically consistent tangents (Section 2.1.1).From the current estimation of the top lateral displacementU and axial load P, the code returns the vector x
=
s
,
v
,
θ
,
F
T that satisfies equilibrium and kinematics (see (17)). In each global Newton iteration the current values of f
s, f
θ, ∂
f
s
/
∂
s, and ∂
f
θ
/
∂
θ based on the current values of s 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.).


## 2.1. 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.


### 2.1.1. 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–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 forcef
s
(
s
) is mobilized at the shear spring and a moment f
θ
(
θ
) is mobilized at the rotational spring. In what follows, we use a general force f
(
x
) to denote f
s
(
s
) and f
θ
(
θ
), depending on which spring we are modelling.The constitutive relationship forf is expressed in the Bouc-Wen model as a linear part and a hysteretic part:(1)
f
=
α
k
x
+
1
-
α
f
y
z
,where x is the deformation; k is the preyield stiffness; f
y is the mobilized force at the beginning of yielding; α is the post-to-preyield stiffness ratio; and z is a hysteretic quantity controlling the nonlinear behavior. The latter is governed by the following differential equation with respect to time t:(2)
z
˙
=
A
x
˙
-
β
x
˙
z
z
n
-
1
+
γ
x
˙
z
η
ν
η
.Equation (2) represents the Baber and Noori [34] formulation of the rate of hysteretic deformation, which accommodates degradation. In (2) β, γ, and n are parameters that control the shape and transition from elastic to inelastic regions of the hysteretic loop, while A, ν, and η are variables that control the stiffness degradation and material deterioration. Note that the yield deformation of the spring (
x
y
) 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(3)
β
=
β
′
x
y
η
,
γ
=
γ
′
x
y
ηin which β
′ and γ
′ can take values as in the original Bouc-Wen model of hysteresis. In this case, z is a dimensional quantity and is bounded between -
x
y and +
x
y. The evolution of material degradation is governed by the following equations [34]:(4a)
A
=
A
0
-
δ
A
e
,
(4b)
ν
=
1
-
δ
ν
e
,
(4c)
η
=
1
+
δ
η
e
,where e is defined by the following rate equation:(5)
e
˙
=
1
-
α
f
y
x
˙
z
.And δ
A, δ
ν, and δ
η are the parameters that control the degradation.Incremental Response Equations. The tangent of f
(
x
) (i.e., ∂
f
(
x
)
/
∂
x) is computed using (1):(6)
∂
f
x
∂
x
=
α
k
+
1
-
α
f
y
∂
z
∂
x
.Note that this is the continuum tangent and not the algorithmically consistent tangent. It is clear from (6) that the derivative of z with respect to x is needed. To get it, (2) can be rewritten as(7)
z
˙
=
A
x
˙
-
x
˙
z
n
β
sgn
⁡
x
˙
z
+
γ
ν
η
=
∂
z
∂
x
∂
x
∂
t
.Hence,(8)
∂
z
∂
x
=
A
-
z
η
β
sgn
⁡
x
˙
z
+
γ
ν
η
.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 t
i
+
1 is obtained as(9)
f
i
+
1
=
α
k
x
i
+
1
+
1
-
α
f
y
z
i
+
1
.The variable z is next discretized by a Backward-Euler solution scheme:(10)
z
i
+
1
=
z
i
+
1
+
Δ
t
A
i
+
1
-
z
i
+
1
n
β
sgn
⁡
x
i
+
1
-
x
i
/
Δ
t
z
i
+
1
+
γ
ν
i
+
1
η
i
+
1
x
i
+
1
-
x
i
Δ
t
.It is seen that Δ
t cancels from the equation, yielding a nonlinear equation in z
i
+
1. A Newton scheme of the form x
i
+
1
=
x
i
-
f
x
m
/
f
′
x
m to solve the general nonlinear equation of the form f
x
=
0 is employed to solve for z
i
+
1 in (10) [35].The equations describing the degradation behavior are described as follows:(11a)
A
i
+
1
=
A
0
-
δ
A
e
i
+
1
,
(11b)
ν
i
+
1
=
1
-
δ
ν
e
i
+
1
,
(11c)
η
i
+
1
=
1
+
δ
η
e
i
+
1
,where e
i
+
1 is found by discretization of the rate equation in (5) using the Backward-Euler solution scheme:(12)
e
i
+
1
=
e
i
+
1
+
Δ
t
1
-
α
f
y
x
i
+
1
-
x
i
Δ
t
z
i
+
1
,where Δ
t again cancels. Note that z
i, x
i, and e
i are history variables that must be stored at each converged step.For the complete implementable procedure to obtainf
i
+
1 and the algorithmically consistent tangent ∂
f
i
+
1
/
∂
x
i
+
1 (referred to here as the Local Newton Routine), the reader is referred to [35].


### 2.1.2. Vertical Spring

The vertical forcef
v mobilized in the vertical translational spring is given by(13)
f
v
=
k
v
Vin which k
v is the elastic vertical stiffness of the soil-foundation system and V 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.


## 2.1.1. 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–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 forcef
s
(
s
) is mobilized at the shear spring and a moment f
θ
(
θ
) is mobilized at the rotational spring. In what follows, we use a general force f
(
x
) to denote f
s
(
s
) and f
θ
(
θ
), depending on which spring we are modelling.The constitutive relationship forf is expressed in the Bouc-Wen model as a linear part and a hysteretic part:(1)
f
=
α
k
x
+
1
-
α
f
y
z
,where x is the deformation; k is the preyield stiffness; f
y is the mobilized force at the beginning of yielding; α is the post-to-preyield stiffness ratio; and z is a hysteretic quantity controlling the nonlinear behavior. The latter is governed by the following differential equation with respect to time t:(2)
z
˙
=
A
x
˙
-
β
x
˙
z
z
n
-
1
+
γ
x
˙
z
η
ν
η
.Equation (2) represents the Baber and Noori [34] formulation of the rate of hysteretic deformation, which accommodates degradation. In (2) β, γ, and n are parameters that control the shape and transition from elastic to inelastic regions of the hysteretic loop, while A, ν, and η are variables that control the stiffness degradation and material deterioration. Note that the yield deformation of the spring (
x
y
) 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(3)
β
=
β
′
x
y
η
,
γ
=
γ
′
x
y
ηin which β
′ and γ
′ can take values as in the original Bouc-Wen model of hysteresis. In this case, z is a dimensional quantity and is bounded between -
x
y and +
x
y. The evolution of material degradation is governed by the following equations [34]:(4a)
A
=
A
0
-
δ
A
e
,
(4b)
ν
=
1
-
δ
ν
e
,
(4c)
η
=
1
+
δ
η
e
,where e is defined by the following rate equation:(5)
e
˙
=
1
-
α
f
y
x
˙
z
.And δ
A, δ
ν, and δ
η are the parameters that control the degradation.Incremental Response Equations. The tangent of f
(
x
) (i.e., ∂
f
(
x
)
/
∂
x) is computed using (1):(6)
∂
f
x
∂
x
=
α
k
+
1
-
α
f
y
∂
z
∂
x
.Note that this is the continuum tangent and not the algorithmically consistent tangent. It is clear from (6) that the derivative of z with respect to x is needed. To get it, (2) can be rewritten as(7)
z
˙
=
A
x
˙
-
x
˙
z
n
β
sgn
⁡
x
˙
z
+
γ
ν
η
=
∂
z
∂
x
∂
x
∂
t
.Hence,(8)
∂
z
∂
x
=
A
-
z
η
β
sgn
⁡
x
˙
z
+
γ
ν
η
.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 t
i
+
1 is obtained as(9)
f
i
+
1
=
α
k
x
i
+
1
+
1
-
α
f
y
z
i
+
1
.The variable z is next discretized by a Backward-Euler solution scheme:(10)
z
i
+
1
=
z
i
+
1
+
Δ
t
A
i
+
1
-
z
i
+
1
n
β
sgn
⁡
x
i
+
1
-
x
i
/
Δ
t
z
i
+
1
+
γ
ν
i
+
1
η
i
+
1
x
i
+
1
-
x
i
Δ
t
.It is seen that Δ
t cancels from the equation, yielding a nonlinear equation in z
i
+
1. A Newton scheme of the form x
i
+
1
=
x
i
-
f
x
m
/
f
′
x
m to solve the general nonlinear equation of the form f
x
=
0 is employed to solve for z
i
+
1 in (10) [35].The equations describing the degradation behavior are described as follows:(11a)
A
i
+
1
=
A
0
-
δ
A
e
i
+
1
,
(11b)
ν
i
+
1
=
1
-
δ
ν
e
i
+
1
,
(11c)
η
i
+
1
=
1
+
δ
η
e
i
+
1
,where e
i
+
1 is found by discretization of the rate equation in (5) using the Backward-Euler solution scheme:(12)
e
i
+
1
=
e
i
+
1
+
Δ
t
1
-
α
f
y
x
i
+
1
-
x
i
Δ
t
z
i
+
1
,where Δ
t again cancels. Note that z
i, x
i, and e
i are history variables that must be stored at each converged step.For the complete implementable procedure to obtainf
i
+
1 and the algorithmically consistent tangent ∂
f
i
+
1
/
∂
x
i
+
1 (referred to here as the Local Newton Routine), the reader is referred to [35].


## 2.1.2. Vertical Spring

The vertical forcef
v mobilized in the vertical translational spring is given by(13)
f
v
=
k
v
Vin which k
v is the elastic vertical stiffness of the soil-foundation system and V 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.


## 2.2. 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 forcesf
s
(
s
) and f
θ
(
θ
), respectively. The total potential energy Π of the model with reference to Figure 2 is expressed as follows:(14)
Π
=
1
2
k
θ
θ
2
+
1
2
k
s
s
2
+
1
2
k
v
v
2
-
P
V
-
F
U
,where U and V are the foundation’s top horizontal and vertical displacements, respectively, and are dependent on the depth of embedment H and the spring’s displacements:(15a)
U
=
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
,
(15b)
V
=
H
1
-
cos
⁡
θ
+
s
sin
⁡
θ
+
v
cos
⁡
θ
.In (15b), V, 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.Figure 2
Undeformed and deformed configurations of the proposed macroelement model.The above kinematic equations (see (15a), (15b)), which relate the total displacement U and the vertical displacement V to the internal displacements s,  v, and θ, assume large displacements and large rotations.Considering the model subjected to the axial loadP with the resulting horizontal displacement U, the vector of unknowns is then x
=
s
,
v
,
θ
,
F
T. Using Castigliano’s second theory by imposing the stationary of Π with respect to s, v, and θ yields the following governing equations:(16a)
g
1
=
k
s
s
-
F
cos
⁡
θ
-
P
cos
⁡
θ
=
0
,
(16b)
g
2
=
k
v
v
+
F
sin
⁡
θ
-
P
cos
⁡
θ
=
0
,
(16c)
g
3
=
k
θ
θ
-
P
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
-
F
H
-
v
cos
⁡
θ
-
s
sin
⁡
θ
=
0
.The following equation should also be considered to permit the displacement-controlled cyclic analysis by imposing the top lateral displacement U:(16d)
g
4
=
U
-
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
=
0
.However, when the lateral force-controlled cyclic response is desired, F 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 (
U
,
V
) 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:(17)
g
=
g
1
g
2
g
3
g
4
=
f
s
s
-
F
cos
⁡
θ
-
P
sin
⁡
θ
k
v
v
+
F
sin
⁡
θ
-
P
cos
⁡
θ
f
θ
θ
-
P
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
-
F
H
-
v
cos
⁡
θ
-
s
sin
⁡
θ
U
-
H
-
v
sin
⁡
θ
+
s
cos
⁡
θ
.Note that the resultant moment M  (
M
=
F
H
-
V
+
P
U
) 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 k
s
s and k
θ
θ were replaced by the general forces f
s
(
s
) and f
θ
(
θ
), respectively. These forces are assigned the constitutive model f
(
x
) developed in Section 2.1.1.The systemg
(
x
) can be solved using Newton’s method following the pseudocode provided below (referred to here as the Global Newton Routine):Previous converged solutionx
=
s
,
v
,
θ
,
F
T
(1)
While (g
(
x
)
<
t
o
l
e
r
a
n
c
e)
(2)
Updatef
s, f
θ, ∂
f
s
/
∂
s, and ∂
f
θ
/
∂
θ using the Local Newton Routine
(3)
Computeg
(
x
)
(4)
Compute the Jacobian matrixJ
(
x
)
(5)
d
x
=
J
x
-
1
g
(
x
)
(6)
x
=
x
+
d
x
endwhere the Jacobian matrixJ
(
x
) is defined as J
i
j
=
∂
g
i
/
∂
x
j. The tangents ∂
f
s
/
∂
s and ∂
f
θ
/
∂
θ required in the Jacobian matrix should be algorithmically consistent tangents (Section 2.1.1).From the current estimation of the top lateral displacementU and axial load P, the code returns the vector x
=
s
,
v
,
θ
,
F
T that satisfies equilibrium and kinematics (see (17)). In each global Newton iteration the current values of f
s, f
θ, ∂
f
s
/
∂
s, and ∂
f
θ
/
∂
θ based on the current values of s 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.).


## 3. 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, D
r
=
85
% and D
r
=
45
%, representing high (HD) and low density (LD) soil conditions, respectively.The HD and LD specimens were loaded vertically by300
k
N and 100
k
N, 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, displacement-controlled, 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 Table1. These parameters have been calibrated using experimental moment-rotation and horizontal force-horizontal 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.Table 1
TRISEE: parameters of the numerical model (units: kN, m).
Phase
Shear spring
Rotational spring
Vertical spring
k
/
10
3
f
y
α (%)
k
/
10
3
f
y
α (%)
k
v
/
10
3
HD
I
132.2
12.5
60
58.6
10
60
120
III
70
99
9
35
111
2
80
LD
I
54
3.8
47
25.4
3.8
4
65
III
35
40.4
1
8
33.3
1
27Table 2
Parameters of the Bouc-Wen-Baber-Noori model, HD test.
Phase
Shear spring
Rotational spring
A
β
′
γ
′
n
δ
A
δ
v
δ
η
A
β
′
γ
′
n
δ
A
δ
v
δ
η
I
1
0.5
0.5
1
0
−0.01
0
1
0.5
0.5
1
0
−0.01
0
III
1
0.5
0.5
0.7
0
0
0.1
1
0.1
0.9
0.7
0
0
0.1Table 3
Parameters of the Bouc-Wen-Baber-Noori model, LD test.
Phase
Shear spring
Rotational spring
A
β
′
γ
′
n
δ
A
δ
v
δ
η
A
β
′
γ
′
n
δ
A
δ
v
δ
η
I
1
0.5
0.5
1
0
−0.01
0
1
0.5
0.5
1
0
−0.01
0
III
1
0.5
0.5
0.3
0
−0.01
0.1
1
0.33
0.67
0.2
0
−0.01
0.1Figures3–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.Figure 3
Comparison of experimental and numerical results: HD test, Phase I.Figure 4
Comparison of experimental and numerical results: LD test, Phase I.Figure 5
Comparison of experimental and numerical results: HD test, Phase III.Figure 6
Comparison of experimental and numerical results: LD test, Phase III.The uplift is important for the HD Phase III test, and this can be observed from the S shaped moment-rotation curve (see Figure5). For the LD sand, only plasticity is developed and the uplift is not present.


## 4. 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(18)
g
=
g
1
g
2
g
3
g
4
=
f
s
s
-
F
-
P
θ
k
v
v
-
P
f
θ
θ
-
P
H
θ
+
s
-
F
H
U
-
H
θ
-
s
.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
f, relating a change in foundation forces d
f
f
=
d
F
,
d
P
T to the change in displacements d
U
f
=
d
U
,
d
V
T, is derived in three steps:(1)
Differentiate the equations of equilibrium (in (18): g
1, g
2 and g
3) resulting in(19)
k
e
q
d
V
=
Γ
d
F
f
in whichd
V
=
d
s
,
d
V
,
d
θ
T and the matrices k
e
q and Γ are given by(20)
k
e
q
=
∂
f
s
∂
s
0
-
P
0
k
v
0
-
P
∂
f
θ
∂
θ
-
P
H
,
Γ
=
1
θ
0
1
H
H
θ
+
s
.
(2)
Differentiate the equations of kinematics (see (15a) and (15b)):(21)
d
U
f
=
Γ
T
d
V
.
(3)
Substituted
V from (19) into (21):(22)
d
U
f
=
Γ
T
f
e
q
Γ
d
F
f
,
wheref
e
q
=
k
e
q
-
1. The resultant flexibility matrix, f
f
=
Γ
T
f
e
q
Γ, relates the displacement increment d
U
f to the force increment d
F
f and is given by(23)
f
f
=
∂
f
θ
/
∂
θ
+
P
H
+
∂
f
s
/
∂
s
H
2
∂
f
s
/
∂
s
∂
f
θ
/
∂
θ
-
P
H
-
P
2
∂
f
s
/
∂
s
H
+
P
H
θ
+
s
∂
f
s
/
∂
s
∂
f
θ
/
∂
θ
-
P
H
-
P
2
∂
f
s
/
∂
s
H
+
P
H
θ
+
s
∂
f
s
/
∂
s
∂
f
θ
/
∂
θ
-
P
H
-
P
2
∂
f
s
/
∂
s
H
θ
+
s
2
+
P
H
θ
2
+
2
P
s
θ
∂
f
s
/
∂
s
∂
f
θ
/
∂
θ
-
P
H
-
P
2
+
1
k
v
.Then foundation stiffness matrix k
f is the inverse of the flexibility matrix, f
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.


## 5. 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:(1)
The proposed macroelement model can simulate with a good accuracy the lateral response and rocking of shallow foundations under quasistatic cyclic loadings.
(2)
The uplift could be simulated adequately using well-chosen parameters of the Bouc-Wen-Baber-Noori model of hysteresis.
(3)
The soil squeeze-out phenomenon observed by El 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.



---

*Source: 4514152-2016-08-01.xml*