Acceso abierto

Dynamic interaction between two 3D rigid surface foundations subjected to oblique seismic waves

  
28 mar 2025

Cite
Descargar portada

Introduction

The dynamic interaction between two or more foundations is typically analyzed independently for different types of loads, assuming harmonic variation over time. For external loads, the dynamic response of the foundation is evaluated by determining the dynamic impedance matrix of the soil–foundation system. In the case of seismic wave loads, the dynamic response is determined using the input motion matrix of the soil–foundation system (referred to as kinematic interaction).

To investigate the impact of soil–structure interaction on the dynamic response of the soil–foundation system, various methods have been proposed. One method to simplify this problem involves linear analysis. An effective and alternative approach is the substructuring method. In this approach, the dynamic response of the superstructure elements and the substructures is analyzed separately, as described by Kausel et al. (1978), Aubry and Clouteau (1992), and Messioud et al. (2017).

The analysis of the foundation system can be simplified by focusing on two key factors: the dynamic stiffness at the soil–foundation interface, often referred to as impedance functions, and the driving forces generated by the incident waves, as discussed by Messioud et al. (2016). The kinematic interaction of the foundation under incident waves is represented using a motion matrix, a concept introduced by Wong (1978) and further elaborated by Qian and Beskos (1996). Determining these dynamic impedance functions and the driving forces resulting from incident waves can be particularly intricate, especially when dealing with complex geometric configurations of the foundation and the mechanical properties of the soil, as highlighted in the works of Qian and Beskos (1996) and Messioud et al. (2016).

The substructuring method has been widely employed by numerous researchers to investigate the dynamic interaction within foundation–soil–foundation systems. Qian and Beskos (1996) developed a boundary element method (BEM) to explore the dynamic response of a system comprising two rigid, massless or massive surface foundations with arbitrary plan forms, subjected to various harmonic waves in three-dimensional (3D) conditions. Tham et al. (1998) presented a dynamic response analysis of a finite number of flexible surface foundations exposed to harmonic incident Rayleigh or shear horizontal waves. Karabalis and Mohammadi (1998) examined the interaction between the adjacent square foundations over multi layered isotropic half-space.

Chen (2016) introduced a numerical approach to calculating the dynamic response of a group of rigid surface foundations. This formulation is unconditionally stable and computationally efficient, involving only algebraic calculations. Sbartai (2016) utilized the BEM in conjunction with the thin layer method (BEM-TLM) to determine the dynamic compliance of two rigid surface foundations. More recently, Keawsawasvong and Senjuntichai (2017, 2020), as well as Han et al. (2017), investigated the dynamic response of two rigid foundations resting on a multilayered medium under time-harmonic loading. Zhenning et al. (2018) developed an indirect boundary element method to study the dynamic impedance functions (stiffness coefficients) of a rigid strip foundation.

In these studies, the foundations are assumed to be rigid, either placed on the surface or embedded in a viscoelastic multilayer soil, and subjected to harmonic dynamic loads in three directions. Within the discretized medium encompassing the foundation–soil–foundation system, the relationship between the distributed loads on the elements' borders and the displacements is established through integral equations. The presence of the two foundations imposes compatibility conditions on the border of elements, ensuring that the displacements are consistent with the movement of the rigid bodies. The dynamic response can be determined by enforcing compatibility between the two foundations and the surrounding soil.

This study introduces a 3D modeling approach to examine the interaction between foundation–soil–foundation systems under seismic loading conditions. The analysis employs the BEM-TLM to investigate the dynamic interaction within the foundation–soil–foundation system. The dynamic interaction between the two foundations is considered using the interaction coefficients proposed by Dobry and Gazetas (1988).

A comprehensive analysis of the 3D seismic response of two surface foundations, situated in either a heterogeneous or a homogeneous soil layer and underlain by a rigid substratum, was conducted. The results are presented in terms of displacements, rotations, and torsion for both rigid surface foundations. This study highlights the influence of the vertical angle of incidence of different waves on displacements, rotations, and torsion.

It is important to note that this work builds upon the research conducted by Messioud et al. in 2012, 2016, and 2019. The method proposed in this study can be readily applied by engineers familiar with the use of interaction factors in foundation design.

Dynamic response of foundations to oblique harmonic seismic waves
Physical model and basic equations

The critical phase of this study involves calculating the dynamic response of a foundation exposed to oblique seismic waves. This response is determined by calculating both the dynamic impedance matrix and the driving forces associated with the seismic waves applied at the center of the foundation.

In a continuous medium, calculating the impedance matrix is particularly intricate due to the transcendental nature of wave propagation and the presence of mixed boundary conditions. However, when the medium is discretized both vertically and horizontally, the problem can be transformed into an algebraic one. This transformation relies on the assumption that the displacement variation at the soil–foundation interfaces follows a linear pattern, as discussed in the works of Messioud et al. (2016, 2019).

In this study, we consider two rigid foundations that are square in shape and positioned on the surface of a homogeneous soil layer, bounded by a substratum (as shown in Figure 1). The soil layer has a height denoted by Ht and is assumed to exhibit linear viscoelastic behavior, characterized by its density (ρ), shear modulus (μ), damping coefficient (β), and Poisson's ratio (ν). The first foundation is subjected to incident harmonic oblique seismic waves, which vary over time and include P waves, shear vertical (SV) waves, and shear horizontal (SH) waves.

Figure 1:

Geometry of two rigid foundations subjected to harmonic seismic waves.

Mathematical model of calculation
Dynamic impedance

In the following sections, we will determine the dynamic response within the first foundation, and the interaction between the two foundations will be accounted for using the interaction coefficients proposed by Gazetas (1988). To calculate the total soil displacement matrix, we employ a method that involves applying unit loads successively across the individual elements of the discretized soil. Subsequently, the displacements within the soil can be expressed as follows: {u}=[G]{t} \left\{u \right\} = \left[G \right] \cdot \left\{t \right\} The vectors {u} and {t} are the nodal values of the amplitudes of displacements and traction, respectively, at the soil–foundation interface. [G] is the soil flexibility matrix.

The foundations impose on the different elements of the soil a displacement compatible with rigid body movements. The soil and the foundation are connected by the transfer matrix: {u}=[R].{Δ} \left\{u \right\} = \left[R \right].\left\{\Delta \right\} [R] is the transfer matrix of dimension (6N×3). [R]=[1000zy000z0x001yx0] \left[R \right] = \left[{\matrix{1 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & z \hfill & y \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & {- z} \hfill & 0 \hfill & x \hfill \cr 0 \hfill & 0 \hfill & 1 \hfill & y \hfill & {- x} \hfill & 0 \hfill \cr}} \right] N is the number of elements at the soil–foundation interface.

{Δ} = {Δx, Δy, Δz, Φx, Φy, Φz}T is the displacement vector of dimension (6×1); Δi (i = x, y, z) represents the translations and Φi (i = x, y, z) the rotations.

The load vector {P} applied to the center of the foundation, the equilibrium between the latter, and the forces (traction) distributed over the discretizing elements of the foundation are expressed by the following relation: {P}=[RT].{t} \left\{P \right\} = \left[{{R^T}} \right].\left\{t \right\} Combining equations (2), (3), and (4), the dynamic impedance matrix is obtained by the compatibility condition. {P}=([RT][G]1[R]){Δ}=[K(ω)]{Δ} \left\{P \right\} = \left( {[{R^T}] \cdot {{[G]}^{- 1}} \cdot \left[R \right]} \right) \cdot \left\{\Delta \right\} = \left[{K\left( \omega \right)} \right]\left\{\Delta \right\} K(ω)=[Kx(ω)000Kx,y(ω)00Ky(ω)0Ky,x(ω)0000Kz(ω)0000Kx,y(ω)0Kx(ω)00Ky,x(ω)000Ky(ω)000000Kz(ω)] K(\omega ) = \left[{\matrix{{{K_x}(\omega )} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{K_{x,y}}(\omega )} \hfill & 0 \hfill \cr 0 \hfill & {{K_y}(\omega )} \hfill & 0 \hfill & {{K_{y,x}}(\omega )} \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & 0 \hfill & {{K_z}(\omega )} \hfill & 0 \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & {{K_{x,y}}(\omega )} \hfill & 0 \hfill & {{K_x}(\omega )} \hfill & 0 \hfill & 0 \hfill \cr {{K_{y,x}}(\omega )} \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{K_y}(\omega )} \hfill & 0 \hfill \cr 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & {{K_z}(\omega )} \hfill \cr}} \right] where ω is the circular frequency and [K (ω)] the impedance function matrix of dimension (6×6).

Determination of Green's functions by the TLM

Obtaining Green's functions by the proposed method will subsequently allow us to form the displacement matrix (flexibility) of the discretized soil. The viscoelastic behavior of the soil is easily introduced in this study. According to the correspondence principle, by replacing the elastic constants λ and μ by their complex values: λ=λ(1+2) {{\rm{\lambda}}^ \star} = {\rm{\lambda}}\left( {1 + 2{\rm{i\beta}}} \right) μ=μ(1+2) {{\rm{\mu}}^ \star} = {\rm{\mu}}\left( {1 + 2{\rm{i\beta}}} \right) where β is the hysteresis damping coefficient.

Calculation of the flexibility matrix

The crux of the study concerning dynamic (seismic) response lies in determining the soil flexibility matrix. Once the flexibility matrix for the specific soil volume in question is known, the impedance matrix can be derived based on the geometric characteristics of the discretized model and the physical attributes of the foundation (degrees of freedom). However, calculating the flexibility matrix requires prior knowledge of Green's functions at the center of each element. The Green's function for a layered stratum is obtained by inverting the thin-layer stiffness matrix using a spectral decomposition procedure. The advantage of the thin-layer stiffness matrix technique over the classic transfer matrix technique for finite layers and the finite-layer stiffness matrix technique is that the transcendental functions in the layered stiffness matrix are linearized (Kausel and Peek, 1982; Boumekik, 1985; Messioud et al., 2012, 2016, 2019).

Discretization of the model

The vertical discretization of the soil mass involves dividing each soil layer into a specific number of sublayers, each with a height denoted by hj, and these sublayers are assumed to have identical physical characteristics, as shown in Figure 2b. These sublayers are considered to be horizontally oriented, viscoelastic, and isotropic, characterized by complex constants of Lamé (λj), complex shear moduli (μj), and densities (ρj). It is important to note that this discretization does not encompass the substratum located at a depth of Ht, which is considered to be infinitely rigid. It is assumed that waves reflecting off the substratum experience total reflection, resulting in zero displacements in that region.

Within each sublayer, the displacement is assumed to follow a linear function of the interface displacements above and below it. This assumption holds true when the height of the sublayer is much smaller than the wavelength being considered (typically in the order of λ/10). This approach bears similarities to the finite elements method in that it derives the displacements within each sublayer from the displacements at the middle of the interfaces. The interaction between the sublayers occurs exclusively through the nodes. The overall equilibrium of the soil mass is understood in terms of finite elements, and the degrees of freedom for the soil mass are reduced to the degrees of freedom of the nodes. The stiffness matrix of the soil mass is determined in a manner similar to the approach used in the structural analysis by finite elements. This technique, known as the thin layer method, was developed by Waas and Lysmer (1972) and is primarily employed when analyzing horizontal soil layers.

The displacements U(n), V(n) and W(n) in each sublayer vary linearly between planes and continue in the x, y, and z directions. Thus, the displacements in each sublayer are obtained by linear interpolation of the nodal displacements at the interface of the sublayer (n) as follows: U(n)(z)=(1η)Un+ηUn+1 {U^{\left( n \right)}}\left( z \right) = \left( {1 - \eta} \right){U^n} + \eta {U^{n + 1}} V(n)(z)=(1η)Vn+ηVn+1 {V^{\left( n \right)}}\left( z \right) = \left( {1 - \eta} \right){V^n} + \eta {V^{n + 1}} W(n)(z)=(1η)Wn+ηWn+1 {W^{\left( n \right)}}\left( z \right) = \left( {1 - \eta} \right){W^n} + \eta {W^{n + 1}} where η=(zZn)hn \eta = {{\left( {z - Zn} \right)} \over {{h_n}}} and 0 ≤ η ≤ 1, and U(n), V(n), and W(n) are the displacements along the x-axis, the y-axis and the z-axis as functions of z in the layer j, and with Un, Vn, and Wn, which are their nodal values at the interface layer z = Zn.

This approach has enabled the determination of 3D Green's functions and the flexibility matrix for the half-space. This method relies on an exact solution of the wave equations in the frequency-wavenumber domain, where only axisymmetric Green's functions are required to facilitate a transformation into a cylindrical coordinate system. The connection between the double Fourier transform in a Cartesian coordinate system and the zero-order or first-order Hankel transform in a cylindrical coordinate system allows for derivation of the transformed solution in a Cartesian coordinate system. The axisymmetric Green's functions used in this method are provided by Kausel and Peek in their work in 1982. Gijmn=l=12Naαβϕimlϕjmlk2kl2 G_{ij}^{mn} = \sum\limits_{l = 1}^{2N} {{{{a_{\alpha \beta}} \cdot \phi_i^{ml} \cdot \phi_j^{ml}} \over {{k^2} - k_l^2}}} with: aαβ = 1 si α = β and aαβ=kkl {a_{\alpha \beta}} = {k \over {{k_l}}} si αβ i ,j = x, y, z.

where k and kl are wave numbers; m represents the interface where the load is applied; and n represents the interface where Green's functions are calculated.

Green's functions obtained in this process are complex and serve as the initial data for calculating the flexibility matrix of any given soil mass. However, when considering the geometry of the foundation, a system of Cartesian coordinates is adopted. The expressions for U, V, and W and the Green's functions provided earlier essentially constitute the elements of the flexibility matrix of the soil. Determining this flexibility matrix allows for computation of the impedance functions for one or multiple foundations and provides information about the vibration amplitudes in the vicinity of a foundation. This method has the advantage of transforming the problem into an algebraic form, enabling the determination of Green's functions when combined with the BEM applied to the soil–foundation interface. Consequently, a horizontal discretization of the soil–foundation interface is necessary for this approach to work effectively.

The horizontal discretization involves subdividing any horizontal interface of the ground mass into square elements, denoted by Sk. These elements serve as constant boundary elements where the average displacement of the element is approximated by the displacement at its center, and it is assumed that the stress distribution within the element is uniform. For ease of integration and computational efficiency, the square elements are replaced by disks.

When distributed unit loads (in the x, y, and z directions) are applied to disk j, the Green's functions at the center of the disk I can be determined. By successively applying these loads to all the disks, we can construct the complex flexibility matrix of the soil at a given frequency ω. The discretized model, which allows us to calculate the impedance functions of the foundation, is presented in Figures 2a and 2b. Within this discretized model, equation (2) is expressed in algebraic form as follows: uj=i=1NSGjtids {u_j} = \sum\limits_{i = 1}^N {\int\limits_S {{G_j}{t_i}ds}} Where N is the total number of elements at the soil–foundation interface.

Figure 2a:

Vertical discretization.

Figure 2b:

Calculation model.

Dynamic response of the foundation subjected to seismic excitation

In this study, the first foundation is subjected to seismic excitation. Consider the plane-harmonic waves SH, SV, P, and R of vertical and horizontal incidence θV and θH (Figure 1). The free field motion of the half-space under seismic waves can be expressed at the soil surface (z = 0) by the following relation: {uf}={Uf}e[iω(x. cos θH+y. sin θH)/c] \left\{{{u^f}} \right\} = \left\{{{U^f}} \right\}{e^{\left[{- i\omega \left( {x.\cos {\theta_H} + y.\sin {\theta_H}} \right)/c} \right]}} where {Uf}={Uxf,Uyf,Uzf} \left\{{{U^f}} \right\} = \left\{{U_x^f,U_y^f,U_z^f} \right\} is the vector of the free field motion; and c is the apparent velocity of the incident waves ( c=cP cos θV c = {{{c_P}} \over {\cos {\theta_V}}} or c=cS cos θV c = {{{c_S}} \over {\cos {\theta_V}}} ) for the P and S waves, respectively and being equal to CR for R waves (CR/CS = 0.9325 for a Poisson's ratio v = 1/3).

The presence of a foundation on the surface of a half-space results in a diffraction of the total displacement field {u}, which can be expressed by the following relation: {u}={uf}+{us} \left\{u \right\} = \left\{{{u^f}} \right\} + \left\{{{u^s}} \right\} where {us} represents the scattered wave field, which satisfies the equation of motion (2).

Substituting equation (5) into equation (12) written in terms of the scattered field {us}, the applied force can be obtained by the following: {P}=[K]{Δ}[K*]{Uf} \left\{P \right\} = \left[K \right]\left\{\Delta \right\} - [{K^*}] \cdot \left\{{{U^f}} \right\} where [K*] is the driving force matrix (6×3N) given by the following relation: [K*]=[F][G]1{e[iω(xcos θH+ysinθH)/c]} [{K^*}] = \left[F \right]{[G]^{- 1}}\left\{{{e^{\left[{i\omega \left( {x \cdot \cos {\theta_H} + y \cdot \sin {\theta_H}} \right)/c} \right]}}} \right\} Equation (9) can be replaced with the following alternative form: {Δ}=[C]{P}+[S*]{Uf} \left\{\Delta \right\} = [C] \cdot \left\{P \right\} + [{S^*}]\left\{{{U^f}} \right\} where [C]=[K]−1 the dynamic compliance matrix (6×6) and [S*] the motion matrix (6×3) given by the following formula: [S*]=[C][K*] [{S^*}] = [C] \cdot [{K^*}] When the foundation is loaded only by the seismic waves, the external forces are zero (P = 0), and the seismic response of the foundation is obtained using equation (16) or by the following expression: {Δ}=[S*]{Uf} \left\{\Delta \right\} = \left[{{S^*}} \right]\left\{{{U^f}} \right\} With: [S*]={SxxoooSyySyzoSzySzzoRxySxzRyxooSzxoo} \left[{{S^*}} \right] = \left\{{\matrix{{Sxx} \hfill & o \hfill & o \hfill \cr o \hfill & {Syy} \hfill & {Syz} \hfill \cr o \hfill & {Szy} \hfill & {Szz} \hfill \cr o \hfill & {Rxy} \hfill & {Sxz} \hfill \cr {Ryx} \hfill & o \hfill & o \hfill \cr {Szx} \hfill & o \hfill & o \hfill \cr}} \right\} The displacements and rotations at the center of the foundation are given by the following relations: Δx=SxxUxfΔy=SyyUyf+SyzUzfΔz=SzyUyf+SzzUzfΦx=RxyUyf+SxzUzfΦy=RyxUxfΦz=SzxUxf} \left. {\matrix{{{\Delta_x} = SxxU_x^f} \hfill \cr {{\Delta_y} = SyyU_y^f + SyzU_z^f} \hfill \cr {{\Delta_z} = SzyU_y^f + SzzU_z^f} \hfill \cr {{\Phi_x} = RxyU_y^f + SxzU_z^f} \hfill \cr {{\Phi_y} = RyxU_x^f} \hfill \cr {{\Phi_z} = SzxU_x^f} \hfill \cr}} \right\} An important part of the equation of motion of the soil–structure interaction is the force vector due to incident waves. Physically, it is defined as the vector of driving forces when the foundation is subjected to an incident wave. The vector of the driving forces is replaced by the input motion matrix; it is linked by the dynamic impedance matrix. The input motion matrix represents, physically, the response of a foundation subjected to incident waves (kinematic interaction). This response is considered an important parameter in seismic analysis because it reduces the horizontal motion of the foundation to that of the free field motion.

The coefficients of the input motion matrix [S*] (Equation 18) are valid only for plane waves with a horizontal angle of incidence θH = 90° (x direction). The procedure presented in this work is applicable when determining the input motion matrix [S*] with other angles of horizontal incidence and for other types of excitations, but with caution because the elements constituting the input motion matrix change their positions when the seismic wave changes its direction.

In the case of an angle of incidence θH = 0°, the motion matrix takes the following form: [S*]={SxxoSxzoSyyoSzxoSzzoRxyoRyxoSyzoSzyo} \left[{{S^*}} \right] = \left\{{\matrix{{Sxx} \hfill & o \hfill & {Sxz} \hfill \cr o \hfill & {Syy} \hfill & o \hfill \cr {Szx} \hfill & o \hfill & {Szz} \hfill \cr o \hfill & {Rxy} \hfill & o \hfill \cr {Ryx} \hfill & o \hfill & {Syz} \hfill \cr o \hfill & {Szy} \hfill & o \hfill \cr}} \right\} Depending on the position of the elements in the input motion matrix, the influence of the horizontal angle of incidence on the seismic response is clearly evident through the effect of the coupling of the different coefficients. The terms of displacements, torsion, and rotations are given by the following relations: Δx=SxxUxf+SxzUzfΔv=SyyUyfΔZ=SzxUyf+SzzUzfΦx=RxyUxf+SyzUzfΦy=RyxUxf+SyzUzfΦz=RzyUyf} \left. {\matrix{{{\Delta_x} = SxxU_x^f + SxzU_z^f} \hfill \cr {{\Delta_v} = SyyU_y^f} \hfill \cr {{\Delta_Z} = SzxU_y^f + SzzU_z^f} \hfill \cr {{\Phi_x} = RxyU_x^f + SyzU_z^f} \hfill \cr {{\Phi_y} = RyxU_x^f + SyzU_z^f} \hfill \cr {{\Phi_Z} = RzyU_y^f} \hfill \cr}} \right\} The proposed approach is based on the concept of substructures, which involves dividing the problem into multiple stages, each of which is easier to solve than the overall problem. This method typically separates kinematic and inertial analyses. In this study, we focused on determining the foundation movement by solving the kinematic interaction problem (where the mass of the superstructure is considered to be zero). In this context, the interaction requires considering kinematic interactions in cases where the foundation soil profile has weak characteristics or is of lower quality, with successive layers having very different stiffness, especially in areas of moderate to high seismicity. In conclusion, this method can be successfully applied to soft soils and soils with weak characteristics. It is also suitable for heterogeneous soils or stratified environments limited by a rigid substratum. It is only necessary to introduce the appropriate parameters for each specific case.

Approximate study of the foundation–soil–foundation interaction

The problem under investigation involves two weightless surface foundations positioned on a uniform stratum or half-space, experiencing harmonic seismic waves originating from the midpoint of the first foundation. However, analyzing this intricate system requires dealing with a substantial number of degrees of freedom. A simplification can be made by assuming that the motions of the foundations remain unaffected by the excitation frequency.

In the case of a system comprising foundations and soil, the displacements of both foundations, denoted by Δ1 and Δ2, are identical. These displacements can be derived using the principle of superposition. Consequently, the total displacement of a pair of foundations is given by: Δ1=Δ11+Δ12=Δ11(1+αv(d)) {\Delta _1} = {\Delta _{11}} + {\Delta _{12}} = {\Delta _{11}}\left( {1 + \alpha v\left( d \right)} \right) with Δ12 = αv(d)Δ11.

where αv(d) is the interaction factor for the vertical vibration.

The dynamic response of the two foundations can be readily computed once the dynamic response of a single foundation has been determined. The impact of the vibrations from foundation 1 on the response of foundation 2 can be described using the dynamic interaction coefficient αv, which varies with frequency. αV=αV(ω)=Δ12/Δ11 {\alpha_V} = {\alpha_V}\left( \omega \right) = {\Delta_{12}}/{\Delta_{11}} The dynamic displacement around a vibrating foundation 1 is characterized by the following asymptotic wave expression presented in Dobry and Gastetas (1988), Gazetas and Makris (1991), and Makris and Gazetas (1992): Δ(a)A1a exp (βωa/Cs)exp[iω(tacs)] \Delta \left( a \right) \approx A{1 \over {\surd a}}\exp \left( {- \beta \omega a/{C_s}} \right)exp\left[{i\omega (t - {a \over {{c_s}}})} \right] where a represents the horizontal distance measured from the center of foundation 1, and A = A(z) signifies the amplitude of displacement. The remaining three terms in equation (24) indicate that the amplitude of motion decreases inversely proportional to a(1/2), which is influenced by the hysteretic damping factor, exp(−βωa/Cs). The phase shift in the motion of a specific point within the soil is solely reliant on distance a. On the axis of the adjacent foundation 2, located at a distance a = d, equation (24) yields: Δ12=Δ(r=d)=A1d exp (βωd)exp[iω(td/Cs)] \matrix{{{\Delta_{12}} = \Delta \left( {r = d} \right) =} \cr {A{1 \over {\surd d}}\exp \left( {- \beta \omega d} \right)exp\left[{i\omega \left( {t - d/{C_s}} \right)} \right]} \cr} The displacement Δ2 remains nearly identical at the periphery of foundation 1. To determine the displacement Δ1 of foundation 1 when subjected to its own dynamic loading, equation (24) continues to be employed to yield an approximate expression, as described in the works of Dobry and Gastetas (1988), Gazetas and Makris (1991), and Makris and Gazetas (1992): Δ11A1γ exp (iωt) {\Delta_{11}} \approx A{1 \over {\surd \gamma}}\exp \left( {i\omega t} \right) By dividing the results of equation (25) by those of equation (26), we can derive an approximate expression for the vertical interaction function. Here, a is defined as a = Bx/2, representing half the width of the foundation. αV(da)1/2 exp (βωd/CS) exp (iωd/CS) {\alpha_V} \approx {({d \over a})^{- 1/2}}\exp \left( {- \beta \omega d/{C_S}} \right)\exp \left( {- i\omega d/{C_S}} \right) This simple expression is essential for computing the vertical dynamic response of any group of foundations when the dynamic response of a single foundation is known. It is assumed based on both static and dynamic solutions that no interaction occurs due to the rotational motion of each foundation (Poulos and Davis, 1980; Dobry and Gazetas, 1988). Movements (without lateral displacement of the foundation) are virtually imperceptible a few foundation lengths away and result in a rapidly diminishing stress field around it. Therefore, although both axial and rotational movements occur in foundation 2 due to its own load, only the influence of displacements from foundation 1 on the axial displacements of foundation 2 is considered.

In the case of laterally loaded foundations, the interaction factor αv is contingent on the frequency ω, the distance d, and the angle θ, which denotes the angle between the foundation line and the direction of the applied force. However, it is sufficient to compute ah for just two specific angles: 0° and 90°. αh(θo)αh(0o)cos2θ+αh(90o)sin2θ {\alpha_h}\left( {{\theta^o}} \right) \approx {\alpha_h}\left( {{0^o}} \right)co{s^2}\theta + {\alpha_h}\left( {{{90}^o}} \right)si{n^2}\theta To improve estimations for an arbitrary angle θ, the approximate values of αh(0°) and αh(90°) can be used, which are derived from the formulas presented by Gazetas and Dobry (1984). Specifically, when considering foundation 2 at an angle of 90°, it is influenced by the Swaves originating from foundation 1, which have a phase velocity Cs. On the other hand, when foundation 2 is at an angle of 0°, it is affected by the extension of the compression wave emanating from foundation 1, which propagates with an apparent phase velocity CP. αh(90)αV {\alpha_h}\left( {{{90}^ \circ}} \right) \approx {\alpha_V} αh(0o)αh0(da)1/2 exp (βωd/CP)× exp (iωd/CP) \matrix{{{\alpha_h}\left( {{0^o}} \right) \equiv {\alpha_{h0}} \approx} \cr {{{({d \over a})}^{- 1/2}}{\rm{exp}}\left( {- \beta \omega d/{C_P}} \right) \times {\rm{exp}}\left( {- i\omega d/{C_P}} \right)} \cr} These coefficients play a crucial role in calculating the dynamic response of foundation 2 under two distinct types of loadings. The coefficients for rotational and torsional interactions can be determined based on the vertical and horizontal interaction coefficients, following the methods outlined in Dobry and Gastetas (1988), Gazetas and Makris (1991), and Makris and Gazetas (1992).

When dealing with oblique loading modes, where the horizontal incidence angle θH is set to 0° for the P and SV waves and 90° for the SH shear wave, the interaction coefficients can be approximated through interpolation, considering the apparent velocity for each type of wave ( c=CPcosθV c = {{{C_P}} \over {\cos {\theta_V}}} and c=Cs cos θV c = {{{C_S}} \over {\cos {\theta_V}}} ).

In this work, a parametric study was conducted to define the parameters of the calculation model. The influence of the discretization of the soil–foundation interface was studied. The thickness of a sublayer h must be small enough that the discrete model can transmit waves in an appropriate manner and without numerical distortion. This size depends on the frequencies involved and the velocity of wave propagation. The frequency of loading and the velocity of wave propagation affect the precision of the numerical solution. Kausel and Peek (1982) showed that the thickness of the sublayer must be smaller than a quarter of the wave length λ. Consequently, the maximum dimensionless frequency must not exceed the number of sublayer N divided by 4.

The validity of the BEM-TLM method utilized to analyze the 3D response of foundations under plane-harmonic waves with varying angles of incidence and vibration frequencies, denoted by ao, is confirmed through comparisons with findings from Luco and Wong (1977) and from Qian and Beskos (1996), conducted on semi-infinite ground (Messioud et al., 2016). Furthermore, the BEM-TLM method proposed in our study has been validated against finite element method numerical computations in the research conducted by Messioud and Dias (2023).

Results and discussion
Surface foundations

In this section, a parametric analysis is conducted for a pair of square foundations, each with a side length of Bx = 2a. These foundations are subjected to plane-harmonic waves with varying angles of incidence and vibration frequencies, as illustrated in Figure 1. The soil is characterized by a bedrock depth of Ht = 10 meters (simulating a semi-infinite soil medium), a Poisson's ratio of ν = 1/3, a coefficient of hysteretic damping of β = 0.05, a shear modulus of μ = 1, and a density of ρ = 1.

The critical step in this study involves determining the dynamic response of the first foundation using the BEM-TLM method. Specifically, only the first foundation is exposed to harmonic oblique seismic waves, while the second foundation is linked to the first one through interaction coefficients. The analysis varies the vertical angle of incidence (θV) at values of 30°, 45°, and 60°, while keeping the horizontal angle of incidence (θH) fixed according to the characteristics of the free field motion (θH = 90° for P waves and θH = 0° for SV waves). The displacements, rotations, and torsion are computed by multiplying the input motion matrix with the vector of free field motion amplitudes. The following figures illustrate the changes in displacement vectors as a function of the dimensionless frequency (ao).

Compression P waves

In Figure 3, we can observe the influence of the vertical angle of incidence on the displacements Δx of both foundations (1) and (2).

Figure 3:

Influence of the vertical angle of incidence on the horizontal displacement ΔxH = 0°; wave P). F1, foundation 1; F2, foundation 2.

The figure clearly demonstrates that the displacements Δx are significantly affected by the vertical angle of incidence, particularly when it comes to P waves. This effect is evident across all values of the vertical angle of incidence (θv) and various distances (d). At vertical angles of incidence of 45° and 60°, the horizontal displacement is quite substantial, especially at lower frequencies, while at higher frequencies, the displacement experiences strong attenuation. On the other hand, at a 30° angle of incidence, the displacement variation is relatively uniform, and the displacement amplitude is lower compared to the other angle values.

It's worth noting that the results indicate the second foundation (2) experiences significantly less displacement compared to the first foundation (1). This suggests the emergence of resonance peaks in the imaginary part, with their amplitudes increasing in relation to the displacement of the first foundation. This increase in the imaginary part is particularly pronounced at lower frequencies and is attributed to the rise in damping forces with frequency.

Figure 4 presents the impact of the vertical angle of incidence on the vertical displacement Δz of both foundations. Similar to the horizontal displacement Δx, the vertical displacement Δz is highly influenced by the angle of incidence (θV). The findings in Figure 4 indicate that the displacement Δz increases as the vertical angle of incidence increases. However, it is important to note that the amplitude of displacement for the second foundation is significantly smaller than that of the first foundation. Furthermore, the amplitudes of the displacements for the second foundation are amplified at lower frequencies, and their behavior becomes unpredictable due to wave reflections occurring at the level of the second foundation. This reflection phenomenon complicates the evolution of these amplitudes as the frequency changes.

Figure 4:

Influence of the vertical angle of incidence on the vertical displacement ΔzH = 0°, P wave). F1, foundation 1; F2, foundation 2.

In Figure 5, we observe the impact of the vertical angle of incidence on rotation and the influence of the interaction between the two foundations on the seismic response of the second foundation. The results obtained in Figure 5 highlight several key points:

Figure 5:

Influence of the vertical angle of incidence on the rotation ΦyH= 0°; P wave). F1, foundation 1; F2, foundation 2.

The imaginary part of the response is strongly affected by the increase in the vertical angle of incidence. As the angle of incidence significantly increases, the system exhibits higher damping at high frequencies when compared with an incidence angle of 30°. However, at lower frequencies, an opposite effect is noticeable: the imaginary part of the rotation weakens.

The interaction between the two foundations results in the emergence of resonance peaks, which are influenced by the distance between the foundations. Consequently, the imaginary part of the response is reduced, and the real part increases in comparison to that of the first foundation. This suggests that the presence of the second foundation alters the dynamic behavior of the system, causing a change in the seismic response characteristics.

SV shear wave

Figures 6 and 7 show the influence of the vertical angle of incidence on the displacement vector when the SV waves are applied in the x direction, with soil particles moving in the same direction as the SV wave. These figures likely depict how the angle of incidence affects the displacement behavior under SV wave loading.

Figure 6:

Influence of the vertical angle of incidence on the horizontal displacement ΔxH = 0°; SV wave). F1, foundation 1; F2, foundation 2.

Figure 7:

Influence of the vertical angle of incidence on the vertical displacement ΔzH = 0°; SV wave). F1, foundation 1; F2, foundation 2.

Figure 6 underscores the significant role played by the vertical angle of incidence, demonstrating how the displacements Δx are influenced by the vertical angle of incidence θV. The maximum displacement value is observed when θV is 60°, while the minimum value is reached when θV is 45°. It is noteworthy to observe the influence of distance (d) on the dynamic response of foundation 2: the displacement Δx2 decreases proportionally with the increase in distance.

The variation in the movement of foundation 2 with respect to the angle of incidence is quite significant. Specifically, at an angle of incidence equal to 60°, the displacement Δx2 reaches its maximum value. However, for the other values of the vertical angle of incidence, such as 30° and 45°, the real part of the displacement tends toward zero. This indicates that the amplitudes of the free field motion for the SV waves are nearly zero for these angles, implying that the movement of foundation 2 closely follows the movement of the free field. Consequently, the displacement amplitude of the second foundation differs significantly from that of the first foundation under these conditions.

Figure 7 highlights the significant influence of the vertical angle of incidence on the vertical displacement along the z-axis. The vertical displacements of both foundations, denoted by Δz1 and Δz2, are notably impacted by the obliquity of the seismic wave. There is a substantial increase in displacement as the vertical angle of incidence increases, especially at the angle of 30°. However, for these angles of incidence, the vertical displacement experiences strong attenuation at high frequencies, with a change in sign for frequencies ao ≥2.5. This results in the variation of vertical displacement Δz at an angle of incidence θV = 45° becoming uncontrolled, featuring some relatively weak peaks.

Figure 7 also demonstrates a considerable increase in the imaginary part compared with the displacement Δz due to the compression waves (P waves). This phenomenon occurs because the free field motion of the vertical shear wave is complex, and in this case the amplitude of displacement is influenced by the presence of the imaginary part of the free field motion.

Figure 8 demonstrates the impact of the vertical angle of incidence on rotation. In contrast to rotation caused by compression waves, there is a substantial increase in rotation as the vertical angle of incidence increases, particularly at 60°. This effect is not observed to the same extent for foundation 1. However, for a vertical angle of incidence of 30°, both foundations tend toward zero rotation. This behavior may be attributed to the disparity in the free field motion amplitudes between shear and compression waves.

Figure 8:

Influence of the vertical angle of incidence on the rotation ΦyH=0°; SV wave). F1, foundation 1; F2, foundation 2.

Furthermore, it is notable that the rotation of the second foundation becomes less predictable due to wave reflections within the foundation. The amplitudes of rotation are significantly influenced by the distance between the two foundations. The results indicate that the variation in rotations is inversely proportional to the increase in distance, suggesting that a larger distance between the foundations leads to less pronounced rotations.

SH shear wave

Figures 9 and 10 present the seismic response of the two square foundations when subjected to seismic oblique wave SH loading. In this case, the seismic wave is applied in the y direction, while the particles move along the x-axis.

Figure 9:

Influence of the vertical angle of incidence on the horizontal displacement ΔxH=90°; SH wave). F1, foundation 1; F2, foundation 2.

Figure 10:

Influence of the vertical angle of incidence on the torsion ΦzH = 90°; SH wave). F1, foundation 1; F2, foundation 2.

Figure 9 underscores the significant impact of the non-verticality of the seismic wave on the dynamic displacements of the two foundations. The displacements along the x-axis experience strong attenuation for low vertical angles of incidence, particularly at high frequencies. Moreover, the results indicate that increasing the distance between the foundations has an influence on the variation of displacement amplitudes for foundation 2. This suggests that the separation distance between the foundations affects how their dynamic displacements behave under the influence of the seismic wave.

Figure 10 illustrates the influence of the vertical angle of incidence on the torsion of the foundation under the influence of the SH wave. This figure reveals the impact of the SH wave on the dynamic response of the foundation. At low angles of incidence (30°), the real part of the torsion becomes negative, especially at high frequencies (ao ≥4). Conversely, for larger values of the vertical angle of incidence (θV ≥60°), the imaginary part of torsion approaches zero for all frequencies. This phenomenon highlights how the non-verticality of the SH shear wave affects the dynamic response of the two foundations.

Figure 10 also demonstrates that the torsion of the second foundation is strongly influenced by the presence of the first foundation. The displacement amplitudes for all angles of incidence are significantly reduced due to the radiation effect. The increase in torsion amplitudes is inversely proportional to the increase in the distance between the foundations. The variation of the imaginary part is not controlled, likely due to the appearance of resonance peaks, which can lead to unpredictable behavior in certain frequency ranges.

The study of the distance between the two foundations reveals a strong influence on the dynamic response of foundation 2. The amplitudes of displacements, rotations, and torsions experienced by the second foundation are consistently reduced for various types of harmonic seismic waves. The results indicate that as the distance between the foundations increases, the amplitudes of the displacement vector for the second foundation decrease proportionally. This suggests a clear inverse relationship between the separation distance and the displacement amplitudes of foundation 2.

Furthermore, the results demonstrate that the incidence of seismic waves affects the displacement vector amplitudes of the second foundation in a manner consistent with the behavior observed for the first foundation. In other words, both foundations respond similarly to the influence of seismic wave incidence and distance separation, although the specific amplitudes may differ. This underscores the significance of considering these factors in the analysis of dynamic responses in such foundation systems.

Conclusion

In the context of designing earthquake-resistant structures, the substructural method emerges as a highly appealing approach to addressing soil–structure interaction challenges. An essential component of this approach involves assessing the kinematic response of foundations. This study leverages the substructural method to precisely evaluate the kinematic interaction of 3D foundations situated on a viscoelastic half-space, constrained by a rigid substratum. This research developed numerical tools that combine the BEM and the TLM thin-layer theory. These tools were applied to calculate the input motion experienced by foundations under various traveling seismic waves. The solution was derived using the BEM in the frequency domain, employing the formalism of Green's functions. A constant quadrilateral element was utilized to investigate the seismic response of the foundations.

This study provides insights into the intricate phenomenon of kinematic interaction between two foundations. The dynamic interaction between these foundations is accounted for using the interaction coefficients proposed by Gazetas (1988). The research also implements a seismic analysis of two rigid square foundations situated on the surface of a homogeneous viscoelastic soil, subjected to different types of oblique harmonic plane waves.

This study underscores the significance of considering the angle of incidence of seismic waves in the dynamic behavior analysis of both square and rectangular foundations. Furthermore, it delves into the interactions between a group of foundations situated on the surface of a homogeneous half-space.

The research calculates displacements, rotations, and torsion by taking the product of the motion matrix and the vector of free field amplitudes of motion, with a focus on the vertical angle of incidence. The findings reveal that seismic activity, influenced by kinematic interaction, results in rotations of the foundation and reduction in translational movements during wave passage. Particularly noteworthy is the more pronounced reduction observed in transfer functions when subjected to S waves (SH and SV) compared with P waves.

The harmonic response of foundation groups is significantly impacted by the dynamic interaction between individual foundations. The amplitudes of displacements, rotations, and torsions for the second foundation are notably reduced across various types of harmonic seismic waves. Additionally, the study explores the effect of the distance between two foundations, demonstrating an inverse relationship between distance and displacement vector amplitudes for the second foundation. Furthermore, the incidence of seismic waves affects the displacement vector amplitudes of the second foundation in a manner similar to that of the first foundation.

This research highlights the importance of accounting for foundation–soil–foundation interaction in the design of foundation systems. It reveals that vibrations in the first foundation induce significant displacements, rotations, and torsion in the second foundation, with the amplitudes of these displacements being notably lower for the second foundation compared with the first. This straightforward approach can be applied during preliminary and extended design phases to account for inertial interaction, thus contributing to more robust foundation system designs.

In this work, the Green's function for a layered stratum is obtained by an inversion of the thin-layer stiffness matrix using a spectral decomposition procedure. The advantage of the thin-layer stiffness matrix technique over the classic transfer matrix technique for finite layers and the finite-layer stiffness matrix technique is that the transcendental functions in the layered stiffness matrix are linearized. The linear approximations inherent in the TLM method and the assumptions of linear elasticity in the BEM might restrict the method's accuracy in capturing complex non-linear behaviors. Theoretical constraints should be transparently acknowledged to guide future research directions and improvements. On the practical side, limitations related to computational efficiency and applicability to real project scenarios should be discussed. The method's scalability and performance on large-scale models, as well as its adaptability to irregular geometries, should be addressed, such as incorporating advanced constitutive models, enhancing computational algorithms, or addressing challenges in handling discontinuities or irregularities.