# A macro-element for a circular foundation to simulate 3D soil-structure interaction

S. Grange*, †, P. Kotronis and J. Mazars

Laboratoire 3S-R (Sols, Solides, Structures-Risques) and VOR Research Network, Grenoble Université, CNRS, Domaine Universitaire BP 53, 38041 Grenoble cedex 9, France

This paper presents a non-linear interface element to compute soil-structure interaction (SSI) based on the macro-element concept. The particularity of this approach lies in the fact that the foundation is supposed to be infinitely rigid and its movement is entirely described by a system of global variables (forces and displacements) defined in the foundation's centre. The non-linear behaviour of the soil is reproduced using the classical theory of plasticity. Failure is described by the interaction diagram of the ultimate bearing capacity of the foundation under combined loads.

The macro-element is appropriate for modelling the cyclic or dynamic response of structures subjected to seismic action. More specifically, the element is able to simulate the behaviour of a circular rigid shallow foundation considering the plasticity of the soil under monotonic static or cyclic loading applied in three directions. It is implemented into FedeasLab, a finite element Matlab toolbox. Comparisons with experimental monotonic static and cyclic results show the good performance of the approach.

KEY WORDS: soil-structure interaction; macro-element; foundation; plasticity; uplift; circular

## 1. INTRODUCTION

In structural engineering, soil-structure interaction (SSI) is an important phenomenon that has to be taken into account. Different experimental results [1] and observations after strong earthquakes [2, 3] showed that non-linearities at the soil level (plasticity) and between the soil and the foundation (uplift of the foundation) often result in an isolation of the structure and thus to a reduction in the

Contract/grant sponsor: European Contract LESSLOSS, Risk Mitigation for Earthquakes and Landslides, Sixth Framework Programme; contract/grant number: GOCE-CT-2003-505488

<footnote>*Correspondence to: S. Grange, Laboratoire 3S-R (Sols, Solides, Structures-Risques) and VOR Research Network, Grenoble Université, CNRS, Domaine Universitaire BP 53, 38041 Grenoble cedex 9, France.</footnote>

<footnote>†E-mail: stephane.grange@hmg.inpg.fr</footnote>
forces and the moments developed at its base during an earthquake. Maximum values of stresses are limited because of larger energy dissipation, but more important displacements are generated at the top.

In order to study SSI, three kinds of methods can be found in the literature [4]:

(i) Sub-structuring that aims at decomposing the problem into simpler problems (Kausel's superposition principle: kinematics interaction and inertial interaction [5]). This method is valid only for linear problems.

(ii) Direct methods that use a classical finite element approach (see [6, 7]). They provide very good results but require numerical expertise, good knowledge of the constitutive laws and are costly in terms of computational power.

(iii) Hybrid methods that are a combination of the two previous methods and therefore more attractive on the numerical level. The macro-element approach belongs to this last category.

The macro-element concept developed for geomechanics in [8] consists in condensing all non-linearities into a finite domain and works with generalized variables (forces and displacements) that allow simulating the behaviour of shallow foundations in a simplified way. Several 2D macro-elements exist in the literature [8–19]. The 2D macro-element developed in [10–12] can be used not only for static/cyclic loading but also for dynamic loading (i.e. earthquake) applied in the horizontal direction, considering the plasticity of the soil and the uplift of the foundation.

Inspired by this last work, a new macro-element is developed hereafter. The goal is to compute the behaviour of a circular shallow rigid foundation lying on an infinite space of soil submitted to static or dynamic loading in three directions. The macro-element takes into account the plasticity of the soil (resulting from the thermodynamics of irreversible processes), but not the influence of the uplift. It is implemented into FedeasLab, a finite element Matlab toolbox developed by Prof. F. Filippou and his co-workers in UC Berkeley [20].

After the mathematical description of the macro-element, numerical results compared with experimental tests under monotonic static [21] and cyclic [22] loading are provided to show the good performance of the approach.

## 2. MATHEMATICAL DESCRIPTION OF THE MACRO-ELEMENT

### 2.1. Shape of the foundation and associated dimensionless variables

In order to simplify the problem, the foundation studied hereafter is considered circular (Figure 1). Because of the axial symmetry, the horizontal loads in the directions x and y are computed in a similar way. Furthermore, it is also easier to reproduce the interaction between horizontal forces and moments. Being a macro-element, the foundation is supposed to be infinitely rigid and all non-linearities are condensed in a representative point: the centre of the foundation. Within that framework it is appropriate to work not only with the generalized (global) variables, the vertical force V, the horizontal forces $H_x$, $H_y$, and moments $M_x$, $M_y$, but also with the corresponding displacements, the vertical displacement $u_z$, the horizontal displacements $u_x$, $u_y$, and rotations $\theta_x$, $\theta_y$. The torque moment $M_z$ is not taken into account.

In order to simplify the equations and to be able to compare the results, an adimensional notation of the variables independent of the plastic properties of the soil is used as follows:

(i) reduced horizontal forces: $H'_x$, $H'_y$;
<figure><img src="image_1.png" /><figcaption>Figure 1. Shape of the foundation and generalized variables: (a) forces and (b) displacements.</figcaption></figure>

(ii) reduced vertical force: $V'$;

(iii) reduced moments: $M'_x, M'_y$.

$$ \underline{\underline{\mathbf{F}}} = \begin{pmatrix} V' \\ H'_x \\ M'_y \\ H'_y \\ M'_x \end{pmatrix} = \frac{1}{Sq_{\max}} \begin{pmatrix} V \\ H_x \\ \frac{M_y}{D} \\ H_y \\ \frac{M_x}{D} \end{pmatrix} \quad (1) $$

With $q_{\max}$ being the ultimate compression stress of the soil under a vertical centred load [23–25], $D$ the diameter of the foundation, and $S=\pi D^2/4$ the surface of the foundation.

The ultimate compression stress is given by the following relationship adapted for a circular footing:

$$ q_{\max} = \frac{0.6}{2} \gamma b N_{\gamma} + q_0 N_q + 1.3 c N_c \quad (2) $$

where $q_0$ is the vertical effective stress at the depth the foundation is lain, $N_\gamma$ is the surface term, $N_q$ is the term depending on the depth of the foundation, and $N_c$ is the cohesion term. The relationships allowing the calculation of $N_\gamma$, $N_q$, and $N_c$ are given in [26]. They depend only on the cohesion $c$ of the soil and the friction angle $\phi$.

In a similar way we obtain for the displacements as follows:

(i) reduced horizontal displacements: $u'_x, u'_y$;

(ii) reduced vertical displacement: $u'_z$;

(iii) reduced rotations: $\theta'_x, \theta'_y$.
$$ \underline{\mathbf{u}} = \begin{pmatrix} u'_z \\ u'_x \\ \theta'_y \\ u'_y \\ \theta'_x \end{pmatrix} = \frac{1}{D} \begin{pmatrix} u_z \\ u_x \\ D\theta_y \\ u_y \\ D\theta_x \end{pmatrix} \qquad (3) $$

## 2.2. Plasticity model

The macro-element developed in this study is able to reproduce the behaviour of a shallow circular rigid foundation under monotonic static or cyclic loading in three directions considering the plasticity of the soil.

The plasticity model is obtained within the framework of the classical thermodynamic laws of irreversible processes [27, 28]. The load surface and the failure criterion are inspired from the works of [10–12] and are expressed in generalized reduced variables $V'$, $H'_x$, $M'_x$, $H'_y$, and $M'_y$. The failure surface is the one developed in the works of [29–31].

They are described hereafter as

(i) the elastic behaviour;

(ii) the plastic behaviour (the failure criterion, the loading surface, the kinematic and isotropic hardening rules, the tangency rule, the flow rule).

### 2.2.1. Elastic behaviour.
The constitutive law can be expressed as: $\underline{\mathbf{F}}=\underline{\mathbf{K}}^{\text{el}}\underline{\mathbf{u}}^{\text{el}}$, where displacement and force vectors are dimensionless, with $\underline{\mathbf{K}}^{\text{el}}$ being the elastic stiffness matrix and $\underline{\mathbf{u}}^{\text{el}}$ the elastic displacement array.

The stiffnesses corresponding to both horizontal displacements are identical, thanks to the circular shape of the foundation. The same stands for the rotations. Using the dimensionless notation presented previously, the following dimensionless stiffness matrix is found:

$$ \underline{\mathbf{K}}^{\text{el}} = \begin{pmatrix} K_{zz}^{\text{el}} & 0 & 0 & 0 & 0 \\ 0 & K_{hh}^{\text{el}} & 0 & 0 & 0 \\ 0 & 0 & K_{00}^{\text{el}} & 0 & 0 \\ 0 & 0 & 0 & K_{hh}^{\text{el}} & 0 \\ 0 & 0 & 0 & 0 & K_{00}^{\text{el}} \end{pmatrix} \qquad (4) $$

with

$$ K_{zz}^{\text{el}} = \frac{K_{zz}^{\text{el}} D}{Sq_{\max}}, \quad K_{hh}^{\text{el}} = \frac{K_{hh}^{\text{el}} D}{Sq_{\max}} \quad \text{and} \quad K_{00}^{\text{el}} = \frac{K_{00}^{\text{el}}}{D Sq_{\max}} $$
Elastic stiffnesses are calculated using the real part of the static impedances defined by [32] and given in [4] for a circular shallow foundation:

$$ 
\begin{aligned}
K_{hh} &= \frac{4G_0D}{2-v} \left(1+0.5\frac{R}{2H}\right) \\
K_{zz} &= \frac{2G_0D}{1-v} \left(1+1.28\frac{R}{2H}\right) \\
K_{00} &= \frac{G_0D^3}{3(1-v)} \left(1+0.17\frac{R}{2H}\right)
\end{aligned} 
\qquad (5) $$

Equation (5) is a function of the geometrical properties of the foundation (diameter D), the elastic properties of the soil (shear modulus G₀, Poisson's ratio v), and the profile of the soil (H being the height of the layer of a soil with very good characteristics (bed rock) measured from the surface).

In agreement with the macro-element presented in [10], the cross-coupled elastic terms are neglected for a surface foundation (see also [4]). However, one has to keep in mind that the off-diagonal elastic terms are responsible for the volumetric behaviour; therefore, they are sometimes introduced in some numerical models [9].

### 2.2.2. Plastic behaviour.
The total displacement **u** is defined as $\mathbf{u}=\mathbf{u}^{\text{elpl}}=\mathbf{u}^{\text{el}}+\mathbf{u}^{\text{pl}}$, with $\mathbf{u}^{\text{pl}}$ the plastic displacement array.

### Failure criterion.
The failure criterion comes from the works in [29–31] and has already been used in the 2D macro-element of Crémer [10]. This criterion was initially developed for a shallow strip and rigid foundation in 2D lying on a half space with a homogeneous cohesion. However, it is shown in [21] that the shapes of the load and failure surfaces for a circular footing are very similar.

The adaptation of the previous criterion in 3D is feasible because the foundation studied presents an axial symmetry. The horizontal forces on the one hand and the moments on the other hand are therefore treated in a similar manner, i.e. the interactions between the two horizontal forces or the two moments are described by circles (the denominators for the horizontal forces (moments) are the same, Equation (6)). It results that the adaptation in 3D consists in adding two terms in relation with the horizontal force $H'_x$ and the moment $M'_y$ to obtain a 5D surface as follows:

$$ f_{\infty} = \left( \frac{H'_{x}}{aV'^{c}(1-V')^{d}} \right)^{2} + \left( \frac{M'_{y}}{bV'^{e}(1-V')^{f}} \right)^{2} + \left( \frac{H'_{y}}{aV'^{c}(1-V')^{d}} \right)^{2} + \left( \frac{M'_{x}}{bV'^{e}(1-V')^{f}} \right)^{2} - 1 = 0 \qquad (6) $$

with the coefficients:

(i) a, b defining the size of the surface in the plane ($H'-M'$);

(ii) c, d, e, and f defining the parabolic shape of the surface in the planes ($V'-H'$) and ($V'-M'$), respectively.
<figure><img src="image_2.png" /><figcaption>Figure 2. Evolution of the loading surface [11].</figcaption></figure>

These parameters can be fitted to different experimental results found in the literature (see, for example, the numerical simulations in Section 3).

*Loading surface:* The loading surface must be able to treat independently the vertical force on the one hand and the horizontal forces and moments on the other hand. It must evolve towards the failure surface when the forces increase. As in [11], its evolution has been guided by the need of reproducing the behaviour of the soil and the foundation initially submitted to the weight of the structure, and then loaded along any loading path. This is also understood by looking at the failure criterion (Equation (6)): the interaction between *H* and *M* is a perfect ellipse centred on the origin of the space 0 (Figure 2), whereas the interaction between *H* and *V* is a curve passing by the origin of the space 0, defined only for $V \ge 0$.

During the initialization phase of the gravity loads, the loading surface is reduced to a straight line segment along the $V'$ axis, starting at point 0 (Figure 2) and with a length equal to $\gamma$ (a hardening variable defined in the paragraph 'Evolution of $\gamma$').

$$f_c \equiv V' - \gamma = 0 \qquad (7)$$

In order to parametrize the size of the ellipse in the hyper-plane ($H'_x, M'_y, H'_y, M'_x$) during the subsequent steps of the loading, two kinds of hardening laws—one kinematic (τ) and another isotropic ($\rho$)—are used. This description is essential in order to reproduce the principal behaviour during a cyclic loading, i.e. the recovery of the plastic state and slope developed at the previous cycle [10], see also Figure 4 and the paragraph 'Kinematic hardening variable' (another possible approach is to use the hypoplasticity theory [13, 14, 18]). The analytical expression of the 5D loading surface $f_c$—represented in Figure 2 as a 3D projection in a space ($V', M', H'$)—takes the following form:

$$f_c(\mathbf{F}, \tau, \rho, \gamma) \equiv \left( \frac{H'_x}{\rho a V'^c (\gamma - V')^d} - \frac{\alpha}{\rho} \right)^2 + \left( \frac{M'_y}{\rho b V'^e (\gamma - V')^f} - \frac{\beta}{\rho} \right)^2 + \left( \frac{H'_y}{\rho a V'^c (\gamma - V')^d} - \frac{\delta}{\rho} \right)^2 + \left( \frac{M'_x}{\rho b V'^e (\gamma - V')^f} - \frac{\eta}{\rho} \right)^2 - 1 = 0 \qquad (8)$$

where $\tau = (\alpha, \beta, \delta, \eta)$ is the array containing the kinematic hardening variables, and $\rho$ the isotropic hardening variable, and $\mathbf{F} = (V', H'_x, M'_y, H'_y, M'_x)$ the force loading vector with $F$ the point at its
<figure><img src="image_3.png" /><figcaption>Figure 3. FEM simulations: relationships between $M'$ and $\theta'$, independent of $V'$ [10].</figcaption></figure>

extremity (Figure 2). The relationships driving the evolution of these hardening variables are given in the following paragraphs.

**Kinematic hardening rule:** The kinematic variables $\alpha$, $\beta$, $\delta$, and $\eta$ permit to determine the centre of the ellipse in the hyper-plane ($H'_x$, $M'_y$, $H'_y$, $M'_x$). The evolution of these variables has been obtained by studying the experimental and numerical behaviour of a foundation under a monotonic static loading. More specifically, Gottardi et al. [21] provide the relationships obtained from experimental tests for a circular footing and for different kinds of soils. Crémer [10] uses similar curves obtained with finite element simulations. The different relationships between forces/moments and associated displacements/rotations can be approached using exponential relationships, driven by a first-order differential equation independent of $V'$. For example, the relationship between $M'_y$ and $\theta'_y = \theta_y^{\text{elp}}$ can be approached as follows (Figure 3):

$$M'_y(\theta'_y) = M'_y^{\infty} \left( 1 - \exp \left( \frac{-K_{00}^{\text{el}}}{M'_y^{\infty}} \theta'_y \right) \right) \quad (9)$$

where $M'_y^{\infty}$ is the limit of the curve $M'_y(\theta'_y)$ when $\theta'_y$ tends to infinity. The corresponding first-order differential equation is as follows (where (.) is the symbol for the derivative with respect to time):

$$\dot{M}'_y = K_{00}^{\text{el}} \left( 1 - \frac{M'_y}{M'_y^{\infty}} \right) \dot{\theta}'_y \quad (10)$$

Assuming again the classical partition of the elastic-plastic total displacement $\mathbf{u}$ into an elastic part $\mathbf{u}^{\text{el}}$ and a plastic part $\mathbf{u}^{\text{pl}}$ ($\mathbf{u}=\mathbf{u}^{\text{el}}+\mathbf{u}^{\text{pl}}$) and considering that $\mathbf{F}=\mathbf{K}^{\text{el}}\mathbf{u}^{\text{el}}$, it is easy to link the increment of force with the increment of the associated plastic displacement:

$$\dot{M}'_y = K_{00}^{\text{el}} \left( \frac{M'_y^{\infty}}{M'_y} - 1 \right) \dot{\theta}'^{\text{pl}}_y \quad (11)$$
<figure><img src="image_4.png" /><figcaption>Figure 4. Influence of the rocking on the plasticity mechanism.</figcaption></figure>

Moreover, as $M'_y$ can be represented in the same space as $\beta b V'^e(\gamma - V')^f$ (in Equation (8) the two terms have the same size), the evolution of the kinematic hardening variable $\beta$ is driven by the following relationship:

$$ \dot{\beta} = \frac{1}{2bV'^e(\gamma - V')^f} K'_{00}^{\text{el}} \left( \frac{M_y^{\star\infty}}{\beta} - 1 \right) \dot{\theta}_y^{\text{pl}} \quad (12) $$

where $M_y^{\star\infty}$ is the limit of the curve $\beta(\theta_y^{\text{pl}})$ when $\theta_y^{\text{pl}}$ tends to infinity. It is clear at this point that the evolutions of the other kinematic hardening variables are driven by similar relationships.

For a real structure, when uplift is possible, loading of the foundation on one side almost does not influence its behaviour on the other side. This is clearly shown in Figure 4: when the foundation undergoes a loading in the right direction ($\dot{M}<0$), the soil under the right part of the foundation is plastified, and no stresses are developed under the left part. When the sign of the loading is reversed ($\dot{M}>0$), only the soil under the left part of the foundation is plastified. The soil under the right part remains in the same plastic state reached during the previous phase. If a third loading is applied at the right direction ($\dot{M}<0$), the soil under the right part of the foundation is again plastified, starting from the plastic state (and with the same slope) obtained during the first phase.

The behaviour being different depending on the sign of the loading, two families of kinematic hardening laws and variables are used to describe the evolution of each force. Eight relations and variables are therefore used in the model for the eight forces $\dot{H}_x'>0$, $\dot{H}_x'<0$, $\dot{M}_y'>0$, $\dot{M}_y'<0$, $\dot{H}_y'>0$, $\dot{H}_y'<0$, $\dot{M}_x'>0$, and $\dot{M}_x'<0$. For example, for a radial loading, each kinematic hardening variable has the following expression (for simplicity, only the case of $\beta$ is presented):

$$ \dot{\beta}^{\oplus} = \frac{1}{2bV'^e(\gamma - V')^f} K'_{00}^{\text{el}} \left( \frac{M_y^{\star\infty}}{\beta^{\oplus}} - 1 \right) \dot{\theta}_y^{\text{pl}} \quad (13) $$

$$ \dot{\beta}^{\ominus} = \frac{1}{2bV'^e(\gamma - V')^f} K'_{00}^{\text{el}} \left( -\frac{M_y^{\star\infty}}{\beta^{\ominus}} - 1 \right) \dot{\theta}_y^{\text{pl}} \quad (14) $$
The first equation of the system (Equation (13)) is activated when $\dot{\theta}_y'^{pl} \ge 0$ (corresponding to the sign $\oplus$), and the second equation is activated when $\dot{\theta}_y'^{pl} \le 0$ (corresponding to the sign $\ominus$). Finally, at each step $\beta$ is calculated as follows:

$$ \beta = \beta^{\oplus} + \beta^{\ominus} \qquad (14) $$

For a non-radial loading (i.e. when forces and moments are coupled with a non-linear relationship), the evolution of the hardening variables is more complicated. Indeed, in many cases $\beta_{lim}$ can decrease, whereas $\dot{M}'_y \ge 0$ ($\beta_{lim}$ being the limit value of $\beta$). In other words, even if $M'_y$ increases and $\dot{\theta}_y'^{pl} \ge 0$, $\beta^{\ominus}$ needs to be activated. This is the reason why for a non-radial loading the evolution of the kinematic hardening variables is driven by relationships of the form (Equation (15)):

$$ \begin{aligned} \dot{\beta}^{\oplus} &= \frac{1}{2bV'^e(y-V')f} K_{00}^{el} \left( \frac{M_y^{*\infty}}{\beta^{\oplus}} - 1 \right) |\dot{\theta}_y'^{pl}| \\ \dot{\beta}^{\ominus} &= \frac{1}{2bV'^e(y-V')f} K_{00}^{el} \left( -\frac{M_y^{*\infty}}{\beta^{\ominus}} - 1 \right) |\dot{\theta}_y'^{pl}| \end{aligned} \qquad (15) $$

The value of $M_y^{*\infty}$ is evaluated as follows:

$$ M_y^{*\infty} = \beta_{lim} - \beta^{\ominus} \quad \text{if } \dot{\beta} > 0 \qquad (16) $$

$$ M_y^{*\infty} = -(\beta_{lim} - \beta^{\oplus}) \quad \text{if } \dot{\beta} < 0 $$

The limit $M_y^{*\infty}$ is always positive, less than or equal to 1.

The first equation of the system (Equation (15)) is activated when $\dot{\beta} \ge 0$, whereas the second equation when $\dot{\beta} < 0$. The sign of $\dot{\beta}$ is identical to the sign of $\beta_{lim} - \beta$, where $\beta_{lim}$ is provided by the tangent rule explained in the paragraph 'Tangent rule', and $\beta$ is the one calculated during the previous step.

**Isotropic hardening rule:** The independence of the directions (for $\dot{M} > 0$ and $\dot{M} < 0$) is taken into account using the specific kinematic hardening laws described in the previous paragraphs. However, one can also link the isotropic with the kinematic hardening laws as in [10] in order to introduce the coupling between plasticity and uplift. Indeed, when a plastic state is reached during a new cycle, plastic behaviour is recovered at the same state and with the same slope as before (Figures 4 and 5). The evolution of the loading surfaces describing this property is given in Figure 5.

This property is translated into the mathematical relation:

$$ \dot{\rho} = \left| \frac{d \| \boldsymbol{\tau} \|}{dt} \right| = \left| \frac{\alpha \dot{\alpha} + \beta \dot{\beta} + \delta \dot{\delta} + \eta \dot{\eta}}{\sqrt{\alpha^2 + \beta^2 + \delta^2 + \eta^2}} \right| \qquad (17) $$

with $\|\boldsymbol{\tau}\|$ the norm of the vector $\boldsymbol{\tau}$ and $d./dt$ the symbol for the derivative with respect to time.

Note that the influence of the coupling between plasticity and uplift can be reduced or cancelled by choosing an arbitrary evolution of the isotropic hardening variable as follows:

$$ \dot{\rho} = k \left| \frac{d \| \boldsymbol{\tau} \|}{dt} \right| \quad \text{where } 0 \le k \le 1 \qquad (18) $$
<figure><img src="image_5.png" /><figcaption>Figure 5. Relationships between kinematic and isotropic variables in a plane (M'-H') [10].</figcaption></figure>

*Evolution of γ*: Its evolution is closely dependent on the evolution of the vertical force V'. During the initialization phase, where the foundation is submitted only to the weight of the structure, γ = V'. After this first phase, the evolution of γ is driven by the empirical relationship linking the vertical force and the vertical displacement given by [8] for a monotonic loading:

$$ \gamma = 1 - \exp(-K_{zz}^{\prime el}(u_z^{\prime})) \quad (19) $$

Equation (19) can be written as a first order differential equation:

$$ \dot{\gamma} = K_{zz}^{\prime el}(1-\gamma)\dot{u}_z^{\prime} \quad (20) $$

Nevertheless, the plastic displacements (horizontal displacements and rotations) can also increase the size of the loading surfaces in the direction of V'. Consequently, the evolution of γ also depends on them according to the following expression:

$$ \dot{\gamma} = K_{zz}^{\prime el}(a_1 \dot{u}_z^{\prime pl} + a_2 |\dot{u}_x^{\prime pl}| + a_3 |\dot{\theta}_y^{\prime pl}| + a_4 |\dot{u}_y^{\prime pl}| + a_5 |\dot{\theta}_x^{\prime pl}|) \frac{(1-\gamma)}{\gamma} \quad (21) $$

where a₁, a₂, a₃, a₄, and a₅ are parameters that permit to adjust the influence of each component of the plastic displacement array.

*Tangent rule*: As in [33], a tangent rule is used to ensure an evolution of the loading surface such as, at failure, the representative point of the force's state coincides with the meeting point of the loading surface with the failure surface. In other words, the loading surface never interpenetrates the failure criterion. It becomes tangent to the failure criterion at the final loading point. Failure occurs the first time the loading point intercepts it. The tangent rule has also another function: in the case of a non-radial loading, the final loading point is not known as it can be anywhere on the failure criterion. The tangent rule permits to determine the final loading point and the limits of the kinematic hardening variables.
<figure><img src="image_6.png" /><figcaption>Figure 6. Tangent rule in $E^*=(H_x^*, M_y^*, H_y^*, M_y^*)$.</figcaption></figure>

This rule is described in the space $E^*=(H_x^*, M_y^*, H_y^*, M_y^*)$, where

$$ 
\begin{align} 
H_x^* &= \frac{H_x'}{aV^{\prime c}(\gamma - V')} \nonumber \\ 
M_y^* &= \frac{M_y'}{bV^{\prime e}(\gamma - V')} \tag{22} \\ 
H_y^* &= \frac{H_y'}{aV^{\prime c}(\gamma - V')} \nonumber \\ 
M_x^* &= \frac{M_x'}{bV^{\prime e}(\gamma - V')} \nonumber 
\end{align} 
$$

From Equations (8) and (22) it is obvious that in this space the loading surfaces are circles with centre the extremity of the vector $\tau=(\alpha, \beta, \delta, \eta)$ and radius $\rho$. They evolve inside the failure criterion, which is also a circle when $\gamma=1$ (Figure 6).

The tangent rule adopted is based on the assumption that the final loading point $F^{*\infty}$ is given by the projection of the current loading point $F^*$ on the failure criterion in the direction of the current increment of force $dF^*$. Thus, the centre of the final circle is the extremity of the vector $\tau_{lim}=(\alpha_{lim}, \beta_{lim}, \delta_{lim}, \eta_{lim})$ and is situated on the radius of the failure criterion when $\gamma=1$. $\tau_{lim}$ is calculated using an iterative procedure: at the first iteration the size of the loading surface is kept constant and its centre is evaluated using the tangent rule. At the following iteration, a new position of the centre is calculated considering the increase of the size of the loading surface.

In the space $E^*=(H_x^*, M_y^*, H_y^*, M_y^*)$, the loading surfaces are circles; hence, the tangency is guaranteed. Nevertheless, for complex loadings (non-radial, variable vertical force) the algorithm presented in the previous paragraph must be improved.

**Flow rule:** The normality rule is defined as: $\dot{\mathbf{u}}^{\text{pl}}=\langle\dot{\lambda}\rangle(\partial g/\partial \mathbf{F})$, where $g$ is the flow rule (or plastic potential) and the plastic multiplier $\langle\dot{\lambda}\rangle=\dot{\lambda}$ if $\dot{\lambda}\ge0$ and $\langle\dot{\lambda}\rangle=0$ if $\dot{\lambda}<0$.
<figure><img src="image_7.png" /><figcaption>Figure 7. Associated flow rule: (a) in a plane ($H' - V'$) projection of $\mathbf{n}_f$ on $V'$ axis can be negative and (b) in a plane ($H' - M'$).</figcaption></figure>

In the 4D space defined by the hyper-plane ($H'_x, M'_y, H'_y, M'_x$), plasticity develops according to the normal at the loading surface $\mathbf{n}_f$. If the projection of the normal on a given axis i.e. ($H'_x$) is positive, then the projection of the plastic velocity vector ($u_x^{pl}$) is also positive (Figure 7(b)). The loading surface $f_c$ is sufficient to describe the direction of the plasticity velocity and the flow rule is thus in this hyper-plane associate.

In the other four planes, however, where the vertical force interferes (i.e. ($H'_x, V'$), ($M'_y, V'$), ($H'_y, V'$), ($M'_x, V'$)), the loading surface cannot be used to describe correctly the direction of the plastic velocity vector. Indeed, there are cases where the projection of the normal vector at this surface is negative, whereas the projection of the plastic velocity vector $u_z^{pl}$ must be always positive or equal to zero (Figure 7(a)). This is the reason why the following non-associated flow rule $g$ is used in these four planes:

$$ g(\mathbf{F}, \mathbf{\tau}, \rho, \gamma, V_c') \equiv \left( \frac{H_x'}{\rho \frac{aV_c^{\prime c}}{(\kappa\gamma + V_c')^d} (\gamma - V')^d (\kappa\gamma + V')^d} - \frac{\alpha}{\rho} \right)^2 + \left( \frac{M_y'}{\rho \frac{bV_c^{\prime e}}{(\xi\gamma + V_c')^f} (\gamma - V')^f (\xi\gamma + V')^f} - \frac{\beta}{\rho} \right)^2 + \left( \frac{H_y'}{\rho \frac{aV_c^{\prime c}}{(\kappa\gamma + V_c')^d} (\gamma - V')^d (\kappa\gamma + V')^d} - \frac{\delta}{\rho} \right)^2 + \left( \frac{M_x'}{\rho \frac{bV_c^{\prime e}}{(\xi\gamma + V_c')^f} (\gamma - V')^f (\xi\gamma + V')^f} - \frac{\eta}{\rho} \right)^2 - 1 = 0 \quad (23) $$
<figure><img src="image_8.png" /><figcaption>Figure 8. Representation of the computed flow rule *g*, the loading surface $f_c$, and the corresponding normal vectors for a given loading point *F* in a plane (*M*′−*V*′).</figcaption></figure>

where $V'_c$ is the current vertical force used in order to assure that the flow rule *g* and the loading surface *f*<sub>c</sub> have their intersection on the loading point *F* (*V*<sub>c</sub> represents the first compo- nent of the coordinates of the loading point *F*). *f*<sub>c</sub> and *g* are therefore identical in the planes (*H*'<sub>x</sub>, *M*'<sub>y</sub>), (*H*'<sub>y</sub>, *M*'<sub>x</sub>), (*H*'<sub>x</sub>, *H*'<sub>y</sub>), (*M*'<sub>y</sub>, *M*'<sub>x</sub>), (*M*'<sub>y</sub>, *H*'<sub>y</sub>), (*M*'<sub>x</sub>, *H*'<sub>y</sub>) (i.e. ellipse), whereas *f*<sub>c</sub> and *g* are non associated in (*H*'<sub>x</sub>, *V*'), (*M*'<sub>y</sub>, *V*'), (*H*'<sub>y</sub>, *V*'), (*M*'<sub>x</sub>, *V*'). *κ* and *ξ* are two parameters in order to modify the shape of the flow rule to present a horizontal tangent at a given vertical reduced force *V*′ chosen by the user. In the general case, the horizontal tangent in (*M*′−*V*′) planes is located on *V*'=γ(1−*ξ*)/2 and the horizontal tangent in (*H*′−*V*′) planes is located on *V*'=γ(1−*κ*)/2.

The representation of *g* in the plane (*M*', *V*') is given in Figure 8, where *κ* and *ξ* are taken by default equal to 1 to obtain always an horizontal tangent for *V*'=0. Similar figures are valid for the planes (*H*', *V*').

3. NUMERICAL SIMULATIONS

The macro-element is implemented into FedeasLab, a finite element Matlab toolbox [20]. The return mapping algorithm [34] is used for the plasticity mechanism. Two different simulations are provided hereafter:

(i) To see whether the macro-element is able to give good results under a monotonic static loading, numerical simulations are compared with experimental results obtained from [21].

(ii) The performance of the macro-element under a cyclic loading is tested using the experimental results obtained from the European program TRISEE [22].
<figure><img src="image_9.png" /><figcaption>Figure 9. Swipe test GG03. Step 1: a constant vertical displacement is applied. Step 2: an increasing horizontal displacement is applied.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10. Swipe test GG07. Step 1: once a vertical force $V_0=1600\text{N}$ is reached, it is reduced to $V_0=200\text{N}$. Step 2: an increasing horizontal displacement is applied.</figcaption></figure>

### 3.1. Monotonic static behaviour

Detailed presentation of the tests is presented in [21]. They concern a circular footing of diameter $2R=D=0.1\text{m}$ lying on a sand of a known density. The general procedure followed during the test is described hereafter.

At the beginning, a vertical displacement is applied at the foundation until a given vertical force is reached. Then, the vertical displacement is kept constant while another displacement (horizontal displacement or rotation or a combined displacement) starts increasing. The tests are thus completely displacement controlled. The experimental response of the foundation is represented in the space of forces. These tests allow a direct investigation of the failure surface for a specific permanent settlement, due to the fact that the load path is sweeping out the yield surface ('Swipe Test'). Of course, the failure surface implemented in the macro-element (Equation (6)) is more general, valid whatever be the vertical force, the horizontal forces, and the moments.

Among all the tests realized in [21], five of them are studied and reproduced hereafter numerically with the macro-element (GG03, GG04, GG06, GG07, and GG29).

During the GG03 test, an initial vertical displacement is applied until the vertical force equals $V_0=1600\text{N}$ (Figure 9). During the GG07 test, once the vertical force $V_0=1600\text{N}$ is reached, it
<table><caption>Table I. Parameters used for the Swipe tests.</caption><thead><tr><th>Parameters</th><th>Values</th><th>Remarks</th></tr></thead><tbody><tr><td><strong>Geometry</strong></td><td></td><td></td></tr><tr><td><em>D</em></td><td>0.1 m</td><td>Diameter of the foundation</td></tr><tr><td><strong>Elastic stiffness</strong></td><td></td><td></td></tr><tr><td><em>K</em><sub>zz</sub><sup>el</sup></td><td>15 MN/m</td><td></td></tr><tr><td><em>K</em><sub>hh</sub><sup>el</sup></td><td>0.4 MN/m</td><td></td></tr><tr><td><em>K</em><sub>θθ</sub><sup>el</sup></td><td>0.0035 MNm/rad</td><td></td></tr><tr><td><strong>Plasticity</strong></td><td></td><td></td></tr><tr><td><em>q</em><sub>max</sub></td><td>0.260 MPa</td><td>Ultimate bearing capacity</td></tr><tr><td><em>a</em></td><td>0.33</td><td>Failure criterion coefficients</td></tr><tr><td><em>b</em></td><td>0.28</td><td>Failure criterion coefficients</td></tr><tr><td><em>c</em></td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td><em>d</em></td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td><em>e</em></td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td><em>f</em></td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td><em>κ</em></td><td>0.608</td><td>Flow rule coefficients</td></tr><tr><td><em>ζ</em></td><td>0.412</td><td>Flow rule coefficients</td></tr><tr><td><em>a</em><sub>1</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td><em>a</em><sub>2</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td><em>a</em><sub>3</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td><em>a</em><sub>4</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td><em>a</em><sub>5</sub></td><td>1</td><td>Influence of each displacement component</td></tr></tbody></table>

is reduced to $V_0=200$ N (Figure 10). For both tests an increasing horizontal displacement is later applied to the foundation.

Swipe test GG04 is realized in a similar way as the GG03 test. The only difference is that once the initial vertical force of $V_0=1600$ N is reached, an increasing rotation is applied to the foundation.

During the swipe tests GG06 and GG29, a simultaneous horizontal displacement and rotation are applied at the centre of the foundation with a constant ratio of arctan($D d\theta/du$) = 126° for GG06 and arctan($D d\theta/du$) = 157° for GG29. This loading is again applied after the initial vertical load of $V_0=1600$ N is reached.

Numerically the foundation is discretized using the macro-element with two nodes. The node at the base of the macro-element is considered embedded (all the degrees of freedom are fixed), whereas the loading (displacements or forces) is applied at the upper node.

The parameters of the numerical model (see Table I) are calibrated using the tests GG03 (for the horizontal components) and GG04 (for the rotational components). Results concerning the other tests are calculated without any further calibration. Figure 11 shows the load path in the planes ($H_x-V$) and ($H_x-u_x$) for the experimental tests GG03 and GG07 ((a) and (b)) and the numerical tests ((c) and (d)). Figures 12 and 13 show the comparison between the experimental and numerical results for the GG04 (in ($M_y-V$) and ($M_y-\theta_y$) planes), GG06 and GG29 tests (in ($h=H'_x/V_0$; $m=M'_y/2RV_0$) plane).
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. (a), (b): Experimental results and (c), (d): numerical results for swipe tests GG03 and GG07 in planes $(H_x-u_x)$ and $(H_x-V)$.</figcaption></figure>

Figures 11–13 show that the macro-element is able to reproduce the trends of the behaviour of a circular footing under static loading under different combinations of displacements.

During the GG03 test, a decrease in the vertical force is observed. This is due to the fact that the horizontal loading creates an isotropic hardening which develops plasticity in all directions. The sand is glued under the foundation, implying a residual friction that can be computed through correct tuning of the flow rule parameters. The macro-element correctly reproduces this behaviour. It is interesting to note that the load path follows particularly well the failure criterion of the experiment.

During the GG07 test, the isotropic hardening developed during the application of the vertical load explains the quasi-linear behaviour observed during the first phase of the application of the horizontal displacement. Then, the numerically calculated horizontal force stays constant, whereas the experimental one decreases along the failure surface. This difference comes from the fact that during the experiments a coupling between the different directions exists; hence so the application of a horizontal displacement can generate moments. This behaviour is not reproduced by the macro-element because there is no link between the horizontal displacements and the moments (in the plastic domain the macro-element takes into account only two links: horizontal displacements with vertical forces and rotations with vertical forces).
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Figure 12. (a), (b): Experimental results and (c), (d): numerical results for swipe test GG04 in planes ($M_y-\theta_y$) and ($M_y-V$).</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 13. (a) Experimental results and (b) numerical results for swipe tests GG06 and GG29 in plane ($h-m$).</figcaption></figure>
During the GG04 test, the load path follows again particularly well the failure criterion of the experiment and reproduces correctly the evolution of the $M_y/2R - 2R\theta$ and $M_y/2R - V$ curves.

Finally, for the tests GG06 and GG29, the loop is well simulated. Some differences are obtained from the definition of the flow rule. Calibrating the parameters of the model could certainly lead to better numerical results.

In order to show the behaviour of the macro-element under a 3D loading, the following numerical 3D swipe test is performed. Figure 14 shows the load path in the ($M_x/2R - M_y/2R - V$) space. At the beginning, a vertical displacement is imposed till a constant value. After that, the foundation is driven with an increasing rotation $\theta'_y$ until the moment $M'_y$ reaches a given value. At this stage the rotation is halted and kept constant, and a new increasing rotation $\theta'_x$ is applied to the foundation. Moments in the two directions are clearly developed and at the end the load path is very close to the failure surface.

## 3.2. Cyclic behaviour

Within the European program TRISEE, experimental tests are performed on a shallow 1m x 1m rectangular foundation lying on two different sands: a ‘Low-density’ sand and a ‘High-density’ sand [22]. Sine-shaped horizontal displacement cycles of increasing amplitude are applied at the

<figure><img src="image_21.png" /><figcaption>Figure 14. Numerical simulations of a 3D swipe test: representation of the load path (a) in planes ($M_y/2R - M_x/2R$), (b) ($M_y/2R - V$), and (c) in space ($M_y/2R - M_x/2R - V$).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 15. TRISEE: scheme of the experimental mock-up (dimensions in m).</figcaption></figure>

top of a vertical beam embedded on the foundation (Figure 15). Numerical results concerning the low-density specimen using the macro-element are presented hereafter.

Numerically the structure is discretized using the macro-element and an elastic beam. The node at the base of the macro-element is, as before, considered embedded, whereas the loading is applied at the upper extremity of the beam. The diameter of the numerical footing is calculated to have the same area as that of the experimental one. This leads to a diameter *D* equal to 1.128 m. The parameters of the numerical model are presented in Table II. During the experiment, no uplift of the foundation was observed (only positive settlements—representative of a foundation lying on a low density soil or with low mechanical characteristics). By imposing $\dot{\rho}=0$ ($k=0$)—in other words by adjusting the evolution of isotropic hardening compared with kinematic hardening—it is possible to simulate this kind of behaviour. The other parameters of the macro-element have been calibrated using the experimental moment-rotation and horizontal force-horizontal displacement diagrams.

Figure 16 shows the relationship between the horizontal force and the horizontal displacement reproduced numerically and during the experiment. Results are satisfactory in terms of the maximum horizontal forces. Nevertheless, the experimental loops shift towards the negative direction of the horizontal displacements, although the loading and the geometry of the mock-up are symmetric. As mentioned in [22], during the experiment the foundation presented a permanent horizontal settlement and severe stiffness degradation. In our opinion, this is due to the development of high plastic strains in the soil leading to a local failure (abutment effect for cyclic loading, shear band formation). This is also the explanation of the small loops presented in the reloading phase at each cycle. This type of local behaviour cannot be reproduced by the macro-element.

Figure 17 shows the relationship between the moment and the rocking angle. Results are again very satisfactory. The main differences are that the moment reaches an asymptote ($M_{lim}=40$ kNm)
<table><caption>Table II. Parameters used for TRISEE modelling.</caption><thead><tr><th>Parameters</th><th>Values</th><th>Remarks</th></tr></thead><tbody><tr><td><strong>Geometry</strong></td><td></td><td></td></tr><tr><td>D</td><td>1.128 m</td><td>Diameter of the foundation</td></tr><tr><td><strong>Elastic stiffness</strong></td><td></td><td></td></tr><tr><td>K<sub>zz</sub><sup>el</sup></td><td>30 MN/m</td><td></td></tr><tr><td>K<sub>hh</sub><sup>el</sup></td><td>8.6 MN/m</td><td></td></tr><tr><td>K<sub>00</sub><sup>el</sup></td><td>6 MNm/rad</td><td></td></tr><tr><td><strong>Plasticity</strong></td><td></td><td></td></tr><tr><td>q<sub>max</sub></td><td>0.2 MPa</td><td>Ultimate bearing capacity</td></tr><tr><td>a</td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td>b</td><td>1</td><td>Failure criterion coefficients</td></tr><tr><td>c</td><td>0.9</td><td>Failure criterion coefficients</td></tr><tr><td>d</td><td>0.9</td><td>Failure criterion coefficients</td></tr><tr><td>e</td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td>f</td><td>0.8</td><td>Failure criterion coefficients</td></tr><tr><td>κ</td><td>1</td><td>Flow rule coefficients</td></tr><tr><td>ζ</td><td>1</td><td>Flow rule coefficients</td></tr><tr><td>a<sub>1</sub></td><td>0.8</td><td>Influence of each displacement component</td></tr><tr><td>a<sub>2</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td>a<sub>3</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td>a<sub>4</sub></td><td>1</td><td>Influence of each displacement component</td></tr><tr><td>a<sub>5</sub></td><td>1</td><td>Influence of each displacement component</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>Figure 16. TRISEE: comparison between (a) experimental results and (b) numerical results. Horizontal force vs horizontal displacement.</figcaption></figure>

more quickly with the macro-element than during the experiments. Once again the experimental curves are not symmetric but the phenomenon is less pronounced this time.

Finally, Figure 18 presents the vertical settlement of the structure measured during the test and reproduced by the macro-element. The different plateaus of the experimental curves are well simulated by the model, thanks to the good description of the plastic displacements by the flow rule. Nevertheless, experimental results show that the centre of the foundation rises from time to time slightly. Of course, the macro-element cannot reproduce this behaviour because during the numerical simulation uplift is prohibited.
<figure><img src="image_24.png" /><figcaption>Figure 17. TRISEE: comparison between (a) experimental results and (b) numerical results. Moment *vs* rocking angle.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 18. TRISEE: comparison between (a) experimental results and (b) numerical results. Vertical settlement *vs* time.</figcaption></figure>

During the 'High-density' sand test of the TRISEE program [22], the foundation developed a high negative vertical displacement component (uplift). The macro-element presented in this paper is not able to reproduce satisfactorily this test as it takes into account only material (plasticity of the soil) and not geometrical (uplift) non-linearities.

## 4. CONCLUSION

The macro-element developed within this work gives satisfactory results for simulating the non-linear behaviour of a circular rigid foundation lying on an infinite space submitted not only to a monotonic static loading but also to a cyclic loading. Using global variables it presents the advantage of inducing low computation costs. It is implemented into the Matlab toolbox FedeasLab providing the user with an easy-to-use interface and practical tools for the post treatment of the results.

Possible improvements of the macro-element are necessary, especially concerning the tangent rule. The extension of the macro-element for dynamic loadings is also straightforward. Some initial
results have already been presented in [35]. Nevertheless, the most challenging future directions of this work concern the implementation of the geometrically nonlinear mechanism of uplift by the introduction of a new loading surface and the extension of the macro-element approach in fully 3D conditions.

ACKNOWLEDGEMENTS

The authors are grateful for the financial support of the European Contract LESSLOSS, Risk Mitigation for Earthquakes and Landslides of the Sixth Framework Programme (Project No.: GOCE-CT-2003-505488).

REFERENCES

1. CAFEEL-ECOEST/ICONS. *Thematic Report No. 5. Shear Walls Structures*, Reynouard JM, Fardis MN (eds), Severn RT, Bairrão R (gen. eds). LNEC, September 2001; ISBN 972-49-1891-2.

2. Housner GW. The behaviour of inverted pendulum structures during earthquakes. *Bulletin of the Seismological Society of America* 1963; **53**(2):403-417.

3. Psycharis IN. Dynamic behavior of rocking structures allowed to uplift. *Journal of Earth Engineering and Structural Dynamics* 1983; **11**:57-76, 501-521.

4. Pecker A. *Dynamique des Sols*. Presse, ENPC: Paris, France, 1984.

5. Kausel E, Whitman A, Murray J, Elsabee F. The spring method for embedded foundations. *Nuclear Engineering and Design* 1978; **48**:377-392.

6. Lysmer J, Udaka T, Tsai CF, Seed HB. FLUSH—a computer program for approximate 3D analysis of soil structure interaction problems. *EERC Report. No. 75-30*, University of California, Berkeley, CA, 1975.

7. Prévost JH. DYNAFLOW. A finite element analysis program for static and transient response of linear and nonlinear two and three dimensional systems. *User Manual*, Department of Civil Engineering, Princeton University, 1999.

8. Nova R, Montrasio L. Settlements of shallow foundations on sand. *Géotechnique* 1991; **41**(2):243-256.

9. Cassidy MJ, Byrne BW, Houlsby GT. Modelling the behaviour of circular footings under combined loading on loose carbonate sand. *Géotechnique* 2002; **52**(10):705-712.

10. Crémer C. Modélisation du comportement non linéaire des fondations superficielles sous séismes. *Ph.D. Thesis*, LMT Cachan-ENS Cachan, France, 2001.

11. Crémer C, Pecker A, Davenne L. Cyclic macro-element for soil-structure interaction: material and geometrical non-linearities. *International Journal for Numerical and Analytical Methods in Geomechanics* 2001; **25**(13): 1257-1284.

12. Crémer C, Pecker A, Davenne L. Modelling of nonlinear dynamic behaviour of a shallow strip foundation with macro-element. *Journal of Earthquake Engineering* 2002; **6**(2):175-211.

13. di Prisco C, Nova R, Sibilia A. Analysis of soil-structure interaction of towers under cyclic loading. In *Proceedings of the NUMOG 8*, Rome, Pande GN, Pietruszczak S (eds). Swets & Zeitlinger: Lisse, 2002; 637-642 (invited lecture).

14. di Prisco C, Galli A. Mechanical behaviour of shallow foundations under cyclic loads. *ALERT Workshop 2006* (without proceedings), Technical University (Politecnico) of Milan, Italy, 2006.

15. Martin CM. Physical and numerical modelling of offshore foundations under combined loads. *D.Phil. Thesis*, University of Oxford, Oxford, 1994.

16. Montrasio L, Nova R. Settlements of shallow foundations on sand: geometrical effects. *Géotechnique* 1997; **47**(1):49-60.

17. Paolucci R, Pecker A. Seismic bearing capacity of shallow strip foundations on dry soils. *Soils and Foundations, Japanese Geotechnical Society* 1997; **37**(3):95-105.

18. Pedretti S. Nonlinear seismic soil-foundation interaction: analysis and modelling method. *Ph.D. Thesis*, Dpt Ing Strutturale, Politecnico di Milano, 1998.

19. Tan FSC. Centrifuge and theoretical modelling of conical footings on sand. *Ph.D. Thesis*, University of Cambridge, MA, 1990.

20. Filippou FC, Constandines M. *FedeasLab Getting Started Guide and Simulations Examples*. Department of civil and Environmental Engineering, University of California, Berkeley, CA, 2004.
21. Gottardi G, Houlsby GT, Butterfield R. Plastic response of circular footings under general planar loading. *Géotechnique* 1999; **49**(4):453-469.

22. TRISEE. 3D Site Effects and Soil-Foundation Interaction in Earthquake and Vibration Risk Evaluation, Part 4: Large-scale Geotechnical Experiments on Soil-Foundation Interaction. European Commission, Directorate General XII for Science, Research and Development, 1998.

23. Davis EH, Booker JR. The effect of increasing strength with depth on the bearing capacity of clays. *Géotechnique* 1973; **23**(4):551-563.

24. Matar M, Salençon J. Capacité portante des semelles filantes. *Revue française de Géotechnique* 1979; **9**:51-76.

25. Philipponnat G, Hubert B. *Fondations et Ouvrages en Terre*. Eyrolles: Paris, France, 2003.

26. Caquot A, Kérisel J. *Traité de Mécanique des Sols*. Gautiers-Villars: France, 1966.

27. Lemaitre J, Chaboche JL. *Mécanique des Matériaux Solides*. Dunod, Bordas: Paris, France, 1988.

28. Lemaitre J, Chaboche JL. *Mechanics of Solid Materials*. Cambridge University Press: Cambridge, 1990.

29. Pecker A. Analytical formulae for the seismic bearing capacity of shallow strip foundations. In *Seismic Behavior of Ground and Geotechnical Structures*, e Pinto S (ed.). Balkema: Rotterdam, 1997; 261-268.

30. Salençon J, Pecker A. Ultimate bearing capacity of shallow foundations under inclined and eccentric loads. Part I: purely cohesive soil. *European Journal of Mechanics - A/Solids* 1995; **14**(3):349-375.

31. Salençon J, Pecker A. Ultimate bearing capacity of shallow foundations under inclined and eccentric loads. Part II: purely cohesive soil without tensile strength. *European Journal of Mechanics - A/Solids* 1995; **14**(3):377-396.

32. Gazetas G. Foundations vibrations. In *Foundation Engineering Handbook*, Chapter 15, Fang H-Y (ed.). van Nostrand Reinhold: New York, 1991.

33. Prévost JH. Anisotropic undrained stress strain behaviour of clays. *Journal of the Geotechnical Engineering Division* 1978; **GT8**:1075-1090.

34. Simo JC, Hughes TJR. Mechanics and materials. *Computational Inelasticity*. Springer Interdisciplinary Applied Mathematics, vol. 7. Springer: New York, 1998.

35. Grange S, Kotronis P, Mazars J. Advancement of simplified modelling strategies for 3D phenomena and/or boundaries conditions for base-isolated buildings or specific soil-structure interactions. *Project LESSLOSS Risk Mitigation for Earthquakes and Landslides*, 2006. Available at: http://hal.ccsd.cnrs.fr/ccsd-00102513.