Reduced T-shaped soil domain for nonlinear dynamic soil-bridge interaction analysis

The one-directional three-component wave propagation in a T-shaped soil domain (1DT-3C) is a numerical modeling technique, in a finite element scheme, to investigate dynamic soil-structure interaction (SSI) coupled with seismic site effects, under the assumption of vertical propagation of three-component seismic motion along a horizontal multilayered soil. A three-dimensional elasto-plastic model is adopted for soils, characterized using their shear modulus reduction curve. In this research, the 1DT-3C wave propagation modeling technique is proposed as an efficient tool for bridge design to take into account directly the spatial variability of seismic loading. This approach, in the preliminary phase of bridge study and design, allows the reduction of the soil domain and the easier definition of boundary conditions, using geotechnical parameters obtained with only one borehole investigation for each pier. This leads to a gain in modeling and computational time.

macro-elements (Grange et al. 2009;Abboud 2017). These approaches reduce the computational time but increase the modeling time to calibrate the soil parameters and also the uncertainties due to the numerical calibration of springs and dashpots. On the other side, the three-dimensional (3D) modeling of both soil and structure and the direct solution of the dynamic equilibrium equation for the assembly (Elgamal et al. 2008;Mazzieri et al. 2013) has higher computational cost and uncertainties in soil properties, considering that generally only one borehole investigation per pier is undertaken in the case of a preliminary phase of bridge design.
The one-directional three component (1D-3C) wave propagation approach proposed by Santisi d' Avila and Lopez Caballero (2018) for soil-building interaction analysis strongly reduces the modeling and computational time, under the assumption of vertical propagation and horizontally layered soil. A one-dimensional soil model, of unit area, is assembled with a 3D building structure to consider inertial soil-structure interaction. Quadratic line elements are used for the soil profile, having zero horizontal strains. The input motion and the absorbing boundary condition are imposed in only one node, at the soil-bedrock interface, and the geotechnical parameters are easy to define, needing only one borehole investigation. However, the whole structure base is submitted to the same seismic loading transmitted by the soil column. Consequently, this model is suitable only for structures with rigid shallow foundation and negligible rocking effects. Fares et al. (2019) propose the 1DT-3C wave propagation modeling technique to take into account the foundation deformability and rocking effects in soil-building interaction analysis. The soil model is 3D until a fixed depth, where the SSI is not negligible, and one-dimensional for deeper soil layers (T-shaped soil domain). Fares et al. (2019) verify the 1DT-3C wave propagation approach by comparison with the results of a fully 3D soil domain, assuming in both models the same conditions of vertical propagation and horizontally layered soil. The T-shaped soil domain allows an important reduction in computational time, compared with a fully 3D soil domain, due to the reduced mesh.
The aim of this research is to propose the 1DT-3C wave propagation modeling technique that was first introduced by Fares et al. (2019) for buildings with shallow foundation, as an approach for bridge design, taking into account directly the spatial variability of seismic loading associated to the irregularity of soil layers and variation of soil mechanical properties. The dynamic equilibrium equation is solved directly for the bridge-soil system and a different seismic motion is expected at the base of each pier, depending on local seismic effect and SSI. Moreover, in this study, the proposed modeling technique is applied in the case of deep foundation.

1DT-3C wave propagation model for bridge design
In this research, the 1DT-3C wave propagation approach is extended to the case of bridge design and it is proposed as a valuable solution to overcome the difficulty of taking into account the spatial variability of the seismic loading.
The proposed modeling technique consists of assembling each pier with a T-shaped soil domain having different stratigraphy, for SSI analysis. Even if the same incident seismic motion is given at the soil-bedrock interface (at the same level for all soil profiles), as it propagates along soil profiles having different stratigraphy (Fig. 1), a different motion is estimated at the base of each pier. In this way, seismic site effects, coupled with SSI, are reproduced directly by the solution of the dynamic equilibrium equation for the bridge-soil system.
The 1DT-3C wave propagation approach is efficient for bridges having pier spacing much longer than the T-shaped soil domain length, so that the soil-pier interaction can be considered independent for each pier and the pier-soil-pier interaction is negligible.

Dynamic balance and time integration
The discrete dynamic equilibrium equation for the assembly of soil domain and structure ( Fig. 1), including compatibility conditions, three-dimensional nonlinear constitutive relation and the imposed boundary conditions, is solved directly. ΔD is the displacement increment vector, Ḋ and D are the increment of velocity and acceleration vectors, respectively, i.e. the first and second time derivatives of displacement vector. K and M are the stiffness and consistent mass matrix, respectively, for the assembled soil-structure system. C and F are the assembled damping matrix and load vector, respectively, derived from the imposed absorbing boundary condition at the soil-bedrock interface. The damping in the structure is also taken into account in the damping matrix C. According to Rayleigh's approach (Chopra 2001), the damping matrix for the structural part, representing the damping of structural and nonstructural elements, is assumed as mass and stiffness proportional, using coefficients depending on the first two natural frequencies.
The time integration of the equation of motion (1) is performed through the implicit Hilber-Hughes-Taylor algorithm (Hughes 1987). The three parameters α = − 0.1, β = 0.25 (1 − α) 2 = 0.3025 and γ = 0.5 − α = 0.6 guarantee an unconditionally numerical stability of the time integration scheme and numerical damping to reduce high frequency content, without having any significant effect on the meaningful, lower frequency response. The dynamic equilibrium equation is directly solved using an automatic time step lower than that used for the input signal sampling.

T-shaped soil domain
The soil in each layer is assumed to be a continuous and homogeneous medium, with nonlinear constitutive behavior. Shear and pressure waves propagate vertically in a horizontally layered soil domain, from the top of the underlying elastic bedrock to the surface. In this research, the finite element models are developed with Abaqus software by Dassault Systèmes, but the same modeling technique can be adopted in any commercial finite element code. The soil domain can be modeled using linear 8-node brick elements, quadratic 20-node solid finite elements (reduced integration) or quadratic 27-node solid finite elements (without reduced integration), having three translational degrees of freedom per node. The spatial discretization is done using one element every meter and verifying that this choice is more accurate than the minimum number of nodes per wavelength p = 10, to accurately represent the seismic signal. The maximum frequency, above which the spectral content of the input signal can be considered negligible, is fixed as f = 15 Hz. The number of finite elements per layer is checked according to the maximum reduction of the shear wave velocity occurred during the dynamic process, that affects the wavelength.
The soil extension along the horizontal directions implies zero axial strain in these directions. Tie constraints are imposed (for translational degrees only) at the nodes of soil domain lateral boundaries to obtain zero axial strains in horizontal directions. Imposing zero axial strains through the tie constraint between the lateral boundaries of soil domain limits the extension of the domain, maintaining the soil profile natural frequency. If it is decided to adopt absorbing lateral conditions, using dashpots at the lateral boundaries of soil domain for the three components of motion, the tie constraints (imposed to obtain zero horizontal strains) can be adopted right after. A constraint equation is used to condense out the degrees of freedom at the base of the 3D soil domain to those at the top of the unit area soil column (Fig. 2).
In a 1D wave propagation model, the area A of the soil column appears as a constant factor in each term of the equilibrium equation. Consequently, the free-field motion can be correctly obtained even if a unit area is adopted. This is not the case in SSI analyses where the area of each part (structure and soil domain) is different. In a commercial finite element software, the soil domain area A, adopted for the 3D soil domain, has to be taken into account in the balance by defining a corrected soil density ρA and elasticity modulus in compression EA in the unit area soil column.
The soil column is bounded at the bottom by a semi-infinite bedrock having elastic behavior. An absorbing boundary condition is imposed at the soil-bedrock interface (as adopted by Joyner and Chen 1975), to take into account the finite rigidity of the bedrock and allow energy to be radiated back into the underlying medium. The seismic loading is applied at the soil column base in terms of horizontal forces ρ b v sb A i ḋ x − 2 v bx and ρ b v sb A iḋy − 2 v by , and vertical force ρ b v pb A i ḋ z − 2 v bz . The first factor (ρ b v sb A i for horizontal directions x and y, ρ b v pb A i for the vertical direction z) is imposed as damping coefficient of dashpots at each node of the soil column base. The parameters ρ b , v sb and v pb are the bedrock density and shear and compressional wave velocities in the bedrock, respectively. A i = A/n is the influence area of each node, A is the soil domain area and n is the number of nodes at the soil-bedrock interface (Fares et al. 2019). The three terms ḋ x , ḋ y and ḋ z are the unknown velocities (incident and reflected motion), at the soil-bedrock interface, in x-, y-and z-direction, respectively, that are evaluated during the process. The seismic loading components are applied at the bottom of the soil column in terms of forces ( The three components of the incident seismic motion at the bedrock level in terms of velocity v bx , v by and v bz , in x-, y-and z-direction, respectively, can be obtained by halving the seismic motion at the outcropping bedrock. The soil model is assumed 3D until a fixed depth h (Fig. 2), where the SSI is not negligible, and one-dimensional for deeper soil layers. The depth h of the 3D soil domain is fixed after comparison of the maximum shear strain profile with depth obtained for both free-field condition and SSI analysis (using a unit area soil column assembled with the pier), in the linear elastic regime. The effect of SSI is considered not negligible where the maximum shear strain profiles with depth are remarkably different.
A guide to assist users in using the 1DT-3C wave propagation modeling technique in Abaqus is presented by Fares (2018).

Constitutive models
The proposed modeling technique is independent on the constitutive relationship selected for soils in a total or effective stress analysis. In this research, the Iwan's elastoplastic model (Iwan 1967;Joyner 1975; Santisi d' Avila and Lopez Caballero 2018) is adopted to represent the 3D hysteretic behavior of soil in total stress conditions.
The calibration of Iwan's constitutive model (Iwan 1967;Joyner 1975) is done using the elastic moduli in shear and compression, G 0 and E 0 , respectively, and the nonlinearity of soil is characterized using the shear modulus reduction curve, employed to deduce the size of the yield surface. The soil damping is purely hysteretic. The linear kinematic model approximates the hardening behavior with a constant rate of hardening (Prager hardening rule). The plasticity model assumes an associated plastic flow, which allows for isotropic yield. In the present study, the adopted shear modulus reduction curve is written in the form G(γ)/G 0 = 1/(1 + |γ/γ r0 |), where γ r0 is a reference shear strain corresponding to a tangent shear modulus G(γ) equivalent to 50% of the elastic shear modulus G 0 . This shear modulus reduction curve provides a hyperbolic stress-strain relationship (Hardin and Drnevich 1972), having asymptotic shear stress τ 0 = G 0 γ r0 in the case of simple shear. If no additional information is available, the normalized compressional modulus reduction curve E/E 0 is assumed equal to the shear modulus reduction curve G/G 0 , under the commonly agreed hypothesis of constant Poisson's ratio during the time history, in undrained conditions.
The 3D Iwan constitutive model is present in Abaqus software as an elasto-plastic model with isotropic and multilinear kinematic hardening (Fares 2018). The size of the yield surface is imposed by the half-cycle first loading curve in the uniaxial stress case, in terms of axial stresses and strains. The backbone curve is deduced by the compressive modulus reduction curve. It is discretized using almost 100 intervals and the nonlinear kinematic hardening with ratcheting is modeled using 10 backstresses (kinematic shift of the yield surface). If resonant column tests provide shear modulus reduction curves G/G 0 (γ), the first loading curve is evaluated as σ(ε) = E/E 0 (ε) E 0 ε, where the axial stress σ(ε) can be calculated from shear stress τ(γ) as σ (ε) = √ 3 τ (γ ) , E/E 0 (ε) is the reduction curve of normalized elastic modulus in compression versus axial strain ε that is assumed equal to G/G 0 (γ) and ε = √ 3 γ G 0 /E 0 . In this research, the nonlinear behavior of soil only is taken into account. According to Fares (2018), soil nonlinearity effect is predominant on the structural seismic response compared with the structure nonlinearity. In fact, in the case of a high level of attained soil nonlinearity, the effect of a nonlinear behavior material for the structure is negligible.
In this research, the foundation deformability is directly taken into account by the model and it is assumed that it is rigid enough to remain in a linear-elastic range.
If the nonlinear behavior of the structure is to be taken into account, steel bars have to be modeled. In this case, a bilinear elastic-plastic relationship could be adopted for the steel bars and a concrete plasticity model under cyclic loading, calibrated from a backbone curve for uniaxial compression and with linear behavior until the low tensile strength.

Soil-structure system
The bridge piers are rigidly connected to the foundation plate, using a tie constraint (including rotational degrees of freedom if applicable). As solid elements of soil domain do not have rotational degrees of freedom, generally, it is necessary to block the three rotations of nodes at the base of vertical structural elements. This is not the case for piers modeled using solid elements.
The bridge foundation plate is rigidly connected to the soil. It is a partition of the T-shaped domain to which the mechanical parameters of reinforced concrete (RC) are assigned as material properties (Fig. 3). RC piles are rigidly connected to the 3D soil domain. The pile depth represents a minimum threshold for the 3D soil domain depth h.

Soil-bridge interaction analysis
A short and regular beam bridge is modeled to illustrate the 1DT-3C wave propagation model, but the proposed approach can be adopted for different types of bridges and for piers with different height. The bridge model is assembled with different T-shaped soil domains (one per pier) in a finite element scheme (Fig. 3), using the stratigraphy of a horizontally layered soil. The incident motion propagates along the soil columns having different stratigraphy. The spatial variability of the seismic motion at the base of the piers is obtained by numerical simulation.

Bridge structure
A 200 m long RC beam bridge having only four spans of 50 m and three piers is selected as application of the proposed modeling technique. It is inspired by the B232 regular bridge tested by Pinto et al. (1996), but the same height for the three piers is assumed for simplification. Moreover, the horizontal and vertical stiffness of the four isolators (lead rubber bearings) between each pier and deck, modeled as springs, are selected arbitrarily equal to K h = 2.32 kN/mm and K v = 1424 kN/mm, respectively. According to Pinto et al. (1996), the deck is a 14 m wide hollow-core concrete girder. The deck extremities are supported at the abutments and they can rotate. The mechanical parameters and deal load are listed in Table 1. Each pier is 14 m high, with a constant 4 m × 2 m rectangular hollow section that is 0.4 m thick. The cross-section of deck and pier is shown in Fig. 4. A damping ratio of 1.6% is fixed for the bridge (Pinto et al. 1996). The first two fundamental frequencies obtained for the fixed base bridge are 0.58 Hz and 1.36 Hz. Each pier is rigidly connected to a RC foundation slab that is 6 m long, 4 m wide and 2 m deep.

Soil stratigraphy
The stratigraphy is freely inspired by the multilayered soil described by Gara et al. (2019). The soil depth is imposed by the position of the soil-bedrock interface, where the incident motion is imposed. In this analysis, the input motion is imposed at a depth of 50 m, for the three soil columns. A stiffer layer (AD5) is added between the limit of  investigations and the assumed soil column base (Table 2). This stiff layer above the bedrock is assumed to behave elastically. The geotechnical parameters are given in Table 3. The three soil profiles (Table 2) have first natural frequency f s1 = 3.4 Hz, f s2 = 3.2 Hz and f s3 = 2.2 Hz. The simulation of SSI effects requires the representation of an adequate soil volume. The soil domain area A = L 2 (Fig. 2) is selected by evaluating the pier base to bedrock transfer function for different soil areas, in the linear elastic regime. The selected soil area is the smallest that provides the soil column fundamental frequency equivalent to the free-field case (Santisi d' Avila and Lopez Caballero 2018). The soil area A = 20 m × 20 m is selected for this analysis. The depth h (Fig. 2) of the 3D soil domain is fixed after comparison of the maximum shear strain profile with depth obtained for both free-field condition and SSI analysis (using a unit area soil column assembled with the pier), in the linear elastic regime. The  Avila et al. Advances in Bridge Engineering (2022) 3:9 effect of SSI is considered not negligible where the maximum shear strain profiles with depth are remarkably different. SSI effects are observed in the first soil layer. Moreover, the depth of piles represents a minimum dimension for the 3D soil part in the T-shaped soil domain. In this example, 6 m deep piles are present under each slab in the first soil layer, having a circular cross-section of 0.4 m diameter. Consequently, the depth h = 8 m of the first soil layer is fixed for the 3D soil domain. As discussed in Section 3.4, displacement and acceleration time history at the structure base are not significantly different from them obtained with a slightly deeper 3D soil domain or adopting a 3D soil domain all along the soil profile.

Input motion
A synthetic wavelet is adopted as outcropping motion and applied, after deconvolution, as seismic loading at the base of each soil profile. The synthetic wavelet is used to highlight the amplitude and phase variation of the seismic motion when arriving at the soil surface and at the base of each pier. It has the following expression (Mavroeidis and Papageorgiou 2002) in terms of velocity: The interest of this specific wavelet is justified in order to control not only the amplitude but also the predominant frequency and the number of cycles. The motion duration is 2 t 0 , where t 0 = n/(2f q ) is the time of envelope peak, the predominant frequency is f q and n = 5 is the selected number of cycles. The peak outcropping acceleration in the (2) d g (t) =ḋ g max /2 1 + cos 2 π f q /n (t − t 0 ) cos 2 π f q (t − t 0 )   Avila et al. Advances in Bridge Engineering (2022) 3:9 direction of the longitudinal axis of the bridge is d g max = 1m/s 2 and 3 m/s 2 , for the analyses considering the linear and nonlinear behavior of soil, respectively. The peak acceleration is assumed equal to 0.8d g max in the direction of the transversal axis of the bridge and 0.5d g max in vertical direction. The sample rate of the synthetic signals is dt = 0.01 s.
The predominant frequency f q of the input signal, as represented in Fig. 5, is selected to excite the first natural frequency of the third soil profile (f q = 2.2 Hz) or the two other soil profiles (f q = 3.3 Hz).

Seismic response of soil and structure
The proposed numerical model allows the simulation of the seismic response of soil and structure taking into account site effects coupled with soil-bridge interaction. First, the seismic response of soil and structure is compared for the 1DT-3C wave propagation approach (Fig. 1) and the case of 3D soil domain (Fig. 6) having area A = 20 m × 20 m (as at the top of the T-shaped soil domain).
The amplification and frequency content for both cases are compared in Fig. 7 where the soil top to bottom transfer function is represented for the three soil profiles ( Table 2). The peak of each transfer function corresponds to the first natural frequency of the related soil profile (f s1 = 3.4 Hz, f s2 = 3.2 Hz and f s3 = 2.2 Hz). The difference between the obtained curves for the T-shaped and 3D soil domain are negligible. The differences in terms of acceleration and displacement time history at the soil-bedrock interface, pier bottom and top (represented in Fig. 8 for the T-shaped soil domain) are trivial.
A higher variation in the soil top to bottom transfer function (and therefore in the motion time histories) could be obtained using an entire 3D soil domain if the hypothesis of independent behavior of each soil column is not assured. This could be the case Soil-bridge interaction model using a 3D domain for the soil profiles d' Avila et al. Advances in Bridge Engineering (2022) 3:9 of bridges having low pier spacing and, at the same time, important soil heterogeneity in the horizontal direction. The computation time for an input signal duration of 4 s (Fig. 5), in linear elastic regime, is 9min for the 1DT-3C wave propagation model (Fig. 1) and 44 min for the 3D soil columns (Fig. 6), using 14 cores of an Intel Core i9 processor.
The variability of the horizontal motion at each pier base and top, in the direction of the longitudinal axis of the bridge, is represented in Fig. 8a-b in the case of linear elastic regime, for a predominant frequency of the input signal f q close to the fundamental frequencies of soil profiles (f q = 3.3 Hz and f q = 2.2 Hz). In both cases, the seismic motion at the base of each pier is amplified if compared with the same component of motion at the soil-bedrock interface. In the case of f q = 3.3 Hz, the peak acceleration at the pier base and peak displacement at the top are similar for the three piers but with a certain delay in soil profile 3. In the case of f q = 2.2 Hz, the ground and structure motion show both a difference in amplitude and a delay in soil profile 3 having a fundamental frequency f s3 = f q . The variability of the seismic motion at the base of each pier is affected by the site effects, resulting in different dynamic features of each soil profile, and soil-pier interaction. This variability could be increased in the case of piers having different height and longer bridges. The variability of seismic motion at the base of piers obtained when the predominant frequency of the input signal is close to the fundamental frequency of structure (f q = 0.5 Hz) is negligible.
In the case where the soil has a nonlinear behavior, for f q = 2.2 Hz, the horizontal acceleration time history at the base of each pier is shown in Fig. 9. The time delay is maintained in the soil column 3 even for nonlinear behavior of soil, nevertheless the resonance effect combined with the plastic soil behavior reduces the acceleration peak.

Conclusions
The aim of this research is to propose a one-directional three-component (1DT-3C) wave propagation model, in a T-shaped soil domain, for the numerical simulation of the response of bridges under seismic loading, taking into account the site effects and dynamic soil-structure interaction. Fig. 9 Time history of the acceleration at the soil-bedrock interface (top) and pier base (middle) and of the displacement at the pier top (bottom), in the direction of the longitudinal axis of the bridge, for a predominant frequency of the input signal f q = 2.2 Hz and nonlinear behavior of soil Whereas the structure is modeled in detail, the 1DT-3C wave propagation modeling technique strongly reduce the soil domain and simplifies the definition of boundary conditions at the soil-bedrock interface. This lead to a gain in modeling and computational time if compared with a fully 3D soil domain. On the other side, it also represents a gain in modeling time and uncertainty reduction if compared with approaches as spring-dashpot and macro-elements. The incident motion propagates along soil columns having different stratigraphy. Consequently, the seismic motion estimated at the base of each pier considers the site effects and soil-pier interaction, taking into account directly the spatial variability of the seismic motion. In this research, the 1DT-3C wave propagation model is applied in the case of deep foundation. The depth of the piles represents a minimum dimensions for the 3D soil part in the T-shaped soil domain.
The 1DT-3C wave propagation approach is efficient for bridges having pier spacing much longer than the T-shaped soil domain length, so that the soil-pier interaction can be considered independent for each pier and the pier-soil-pier interaction is negligible.
The proposed modeling technique can be adopted for different types of bridges and for piers with different height in any commercial finite element code. Moreover, it is independent on the constitutive relationship selected for soils in a total or effective stress analysis.