Dynamic interaction between two 3D rigid surface foundations subjected to oblique seismic waves
Artikel-Kategorie: Original Study
Online veröffentlicht: 28. März 2025
Seitenbereich: 103 - 120
Eingereicht: 04. Okt. 2023
Akzeptiert: 03. Jan. 2025
DOI: https://doi.org/10.2478/sgem-2025-0004
Schlüsselwörter
© 2025 Messioud Salah, published by Sciendo
This work is licensed under the Creative Commons Attribution 4.0 International License.
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.
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.

Geometry of two rigid foundations subjected to harmonic seismic waves.
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:
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:
{Δ} = {Δ
The load vector {
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:
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).
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
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.
where
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
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

Vertical discretization.

Calculation model.
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
The presence of a foundation on the surface of a half-space results in a diffraction of the total displacement field {
Substituting equation (5) into equation (12) written in terms of the scattered field {
The coefficients of the input motion matrix [
In the case of an angle of incidence
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:
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.
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°.
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 (
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).
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).
In Figure 3, we can observe the influence of the vertical angle of incidence on the displacements Δx of both foundations (1) and (2).

Influence of the vertical angle of incidence on the horizontal displacement
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.

Influence of the vertical angle of incidence on the vertical displacement
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:

Influence of the vertical angle of incidence on the rotation
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.
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.

Influence of the vertical angle of incidence on the horizontal displacement

Influence of the vertical angle of incidence on the vertical displacement
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
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 Δ
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.

Influence of the vertical angle of incidence on the rotation
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.
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.

Influence of the vertical angle of incidence on the horizontal displacement

Influence of the vertical angle of incidence on the torsion
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.
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.