# Reduced T-shaped soil domain for nonlinear dynamic soil-bridge interaction analysis – Advances in Bridge Engineering

#### ByMaria Paola Santisi d’Avila, Luca Lenti, Stefania Gobbi and Reine Fares

Aug 20, 2022

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

$$mathbf{M}kern0.5em Delta ddot{mathbf{D}}+mathbf{C};Delta dot{mathbf{D}}+mathbf{K};Delta mathbf{D}=Delta mathbf{F}$$

(1)

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, (Delta dot{mathbf{D}}) and (Delta ddot{mathbf{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 non-structural 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 ({rho}_bkern0.5em {v}_{sb};{A}_i;left({dot{d}}_x-2;{v}_{bx}right)) and ({rho}_bkern0.5em {v}_{sb};{A}_i;left({dot{d}}_y-2;{v}_{by}right)), and vertical force ({rho}_bkern0.5em {v}_{pb};{A}_i;left({dot{d}}_z-2;{v}_{bz}right)). The first factor (ρbvsbAi for horizontal directions x and y, ρbvpbAi for the vertical direction z) is imposed as damping coefficient of dashpots at each node of the soil column base. The parameters ρb, vsb and vpb are the bedrock density and shear and compressional wave velocities in the bedrock, respectively. Ai = 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 ({dot{d}}_x), ({dot{d}}_y) and ({dot{d}}_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 (ρbvsbAi) 2vbx, (ρbvsbAi) 2vby and (ρbvpbAi) 2vbz. The three components of the incident seismic motion at the bedrock level in terms of velocity vbx, vby and vbz, 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 elasto-plastic 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, G0 and E0, 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(γ)/G0 = 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 G0. This shear modulus reduction curve provides a hyperbolic stress-strain relationship (Hardin and Drnevich 1972), having asymptotic shear stress τ0 = G0γr0 in the case of simple shear. If no additional information is available, the normalized compressional modulus reduction curve E/E0 is assumed equal to the shear modulus reduction curve G/G0, 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/G0(γ), the first loading curve is evaluated as σ(ε) = E/E0(ε) E0ε, where the axial stress σ(ε) can be calculated from shear stress τ(γ) as (sigma left(varepsilon right)=sqrt{3};tau left(gamma right)), E/E0(ε) is the reduction curve of normalized elastic modulus in compression versus axial strain ε that is assumed equal to G/G0(γ) and (varepsilon =sqrt{3};gamma;{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.