Time Domain Implementation of Transmitting Boundaries in ABAQUS for Discrete Soil-structure Interaction Systems

Soil-structure-interaction (SSI) analyses are essential to evaluate the seismic performance of important structures before finalizing their structural design. SSI under seismic condition involves much more complex interaction with soil compared to the dynamic loads having source on the structure. Seismic SSI analysis requires due consideration of sitespecific and structure-specific properties to estimate the actual ground motion (scattered motion) experienced at the base of the structure, and subsequently the effects of the scattered motion on the structure. Most challenging aspect of seismic SSI analysis is to implement transmitting boundaries that absorb the artificial reflections of stress waves at the truncated interface of the finite and infinite domains, while allowing the seismic waves to enter the finite domain. In this paper, the time domain implementation of seismic analysis of a soil-structure system is presented using classical discrete models of structure and interactive force boundary conditions for soil. These models represent typical SSI systemsa single Degree of Freedom (DOF) of a spherical cavity with mass attached to its wall, a two DOF system consisting of a mass attached by a nonlinear spring to a semi-infinite rod on elastic foundation, and a three DOF system with additional DOFs for modelling the structural stiffness and damping. The convolution integral representing the force boundary condition on the truncated interface, is evaluated interactively using UAMP user-subroutine in ABAQUS and applied as concentrated forces at the interface (truncated interface) nodes of the bounded domain or generalized-structure domain. The verification problems presented in the paper show the satisfactory performance of the developed MATLAB code and ABAQUS implementation with FORTRAN user-subroutines. The classical phenomena associated with the dynamic soil-structure systems are discussed through the present work. KeywordsSoil-structure interaction, Seismic analysis, Transmitting boundaries, ABAQUS, UAMP subroutine.


Introduction
The accuracy, and thus the usefulness, of a dynamic Soil-Structure Interaction (SSI) analysis depend on the ability to model the boundary at the truncated interface of the finite geometry (bounded domain) and the semi-infinite geometries (unbounded domain). The frequency dependence of the stiffness of the unbounded domain along with the stress waves and their multiple reflections makes the modelling of this domain the most challenging part. Although classical solutions of this SSI problem are obtained either in the frequency domain (Emani and Maheshwari, 2009;Spyrakos and Xu, 2003;Wolf, 1985), or hybrid domain (Emani and Maheshwari, 2010;Maheshwari and Emani, 2015), the dynamic analysis is intuitively done only in time domain using the numerical integration of the equations of motion. The displacement boundary conditions are highly inappropriate at the truncated interface since such boundaries cause considerable reflection of waves leading to chaotic results in the finite geometry (Kramer, 1996). So, a time domain implementation of the (interaction-) force boundary condition is needed motion (Sarkar and Maheshwari, 2012;Solberg et al., 2016;Wegner et al., 2005;Wolf and Song, 1996). That method which can integrate with the well-developed state-of-the-art finite element (FE) solvers like ABAQUS will be of special significance. Recently, researchers are using user subroutines to achieve customization of FE solvers for SSI analysis (Poul and Zerva, 2018;Wang and Yang, 2013). In the following section the results for three representative soil-structure systems-one modelled as a single DOF, another as 2DOF and one other system modelled as 3DOF are presented. The implementation in all three cases is done using UAMP user subroutine in ABAQUS.

Formulations of System and Method of Analysis
Classically a SSI system is modelled as a finite domain with the effect of the truncated soil accounted through an interactive force applied at the truncated interface called as "soil structure interface". Basic equation of motion for such a system in time domain is given by Eq. (1) (Wolf, 1988) [ where [ ]is the partitioned form of mass matrix at the structure nodes (with subscript ss) and at the interface nodes ( In time domain, the dynamic stiffness matrix of truncated unbounded domain has a general form given in Eq. (3) After time discretization as per Newmark-time integration scheme with implicit algorithm, the Eq. (4), at time step = Δ becomes The basic equations of Newmark time integration method for determining total displacements and velocities of the soil-structure interaction nodes 'b' at current time step are Using Eq. (6), (7), the velocities {̇} and accelerations {̈} at current step can be expressed in terms of displacements { } at = Δ . where Alternative to the process described from Eq. (11) to Eq. (14), the convolution integral can also be evaluated using displacement impulse response (or dynamic stiffness in time domain) as follows where [ , ] − is given by Eq. (22) On substituting back into the equation of motion Eq. (1) the formulation can be reduced to Eq.
(16) and can be solved for total displacements. Using Eq. (7) and Eq. (8), the total velocities and accelerations are also evaluated subsequently.

Flow Chart
The flow chart in Figure 1 depicts the implementation of the above equations in a programming language. In the present work, MATLAB is used for this purpose.

Implementation in ABAQUS
ABAQUS is general purpose Finite Element software by Dessault systems. The procedure described in the previous section can also be implemented in ABAQUS using an implicit dynamic analysis step. In this case, the time discretization (involving the Newmark's integration parameters, namely, γ and β) is internally performed by ABAQUS. Thus, only the terms involving convolution integral and the ground motion are to be considered in the interaction forces, which shall act on a fixed base system. This results in the solution in total displacements, velocities and accelerations.
The calculation of convolution integral in ABAQUS needs an interactive procedure at each increment of the dynamic step. This is because the convolution integral term of the applied interaction forces at any increment depends on the displacements, velocities and accelerations at all the previous steps. For this, sensors are created for displacements, velocities and accelerations at the SSI nodes. Using the sensor data at each step, the convolution terms and thus the total interaction forces on SSI nodes are calculated by a user subroutine UAMP. The UAMP subroutine is coded in FORTRAN 77 for calculating the interaction forces at each time increment. These forces when applied to the fixed base structure give the total displacements, velocities and accelerations at all nodes of the model. These steps are summarized in Table 1.

Analysis of SDOF Soil-Structure Interaction (SSI) System
In the simplest possible system representing general features of SSI, a mass is attached to a spherical cavity of radius and subjected to symmetric (radial) stress waves. The medium has a shear modulus of , mass density of and Poisson's ratio of . Such a system can be modelled as a SDOF system with mass supported by frequency dependent springs as shown in Here , represent the P-wave and S-wave velocities in the medium respectively, and is the static stiffness value of the spherical cavity. The regular part of dynamic stiffness in frequency domain is given by Eq. The verification of the developed code and its implementation is checked for an artificial earthquake used by Wolf (1988). The ground motion { ( )}, {̇( )} and the response in terms of the total acceleration of the mass are shown in the Figure 3. The ABAQUS implementation matches the result from MATLAB programme indicating that there are no hidden errors in implementation. The minor discrepancy in the result compared to the result reported by Wolf, 1988 is the sight numerical errors in digitizing the input ground motion used by Wolf, 1988.
The formulation and its implementation in MATLAB and ABAQUS are, thus, verified using the single DOF system of spherical cavity. The same procedure can be applied for any general soilstructure system. The soil modelling can be either continuum or discrete.

Analysis of 2DOF System
The selected system is a mass ( ) connected by a nonlinear spring (of stiffness and to a semiinfinite rod on elastic foundation (Figure 4 (a)). The mass and non-linear spring represents the structure and the semi-infinite rod represents the foundation-soil. The dynamic stiffness coefficient of semi-infinite rod and the expressions for the interaction force history ( ) in terms of convolution integral are given in Eq. (35) to Eq. (39). Here E is the elastic modulus of the soil, A is the cross-sectional area of the rod, is the elastic stiffness of the soil per unit length of the rod. The terms of ( ) containing current (unknown) displacements { } and unknown velocities {̇} will contribute to the stiffness [ ] and [ ] of the foundation-soil system (Eq. (18), Eq. (39). The remaining known terms of the convolution integral shall be calculated in ABAQUS using user-defined UAMP subroutine, and the corresponding amplitude is used to apply a concentrated force on mass 0 . Thus, the SSI system reduces to a fixed-base 2DOF dynamic system as shown in Figure 4 (b), where the ground motion is replaced by equivalent interaction force ( ). The two degrees of freedom are the (total) displacements of the two masses and 0 . The numerical properties of the system are as detailed in Table 2. In this case, the dynamics stiffness matrix of the foundation is a scalar value since there is only one DOF at the soil structure interface. The formulation is given by Eq. (35)-Eq. (39).
where = √ and 0 = , = √ (36) So that , = 0 , = and , ( ) = 1 ( ) The system is subjected to an artificial time history, whose displacement and velocity histories are shown in Figure 5 (a). The system on being solved for the seismic response gives an accurate response as can be seen in Figure 5 (b). From the figure, it can also be observed that there is a permanent deformation caused by the elastic-perfect plastic nonlinearity of the structure spring.  Table 2. Parameters for 2DOF system (Wolf, 1988) Parameters assumed Dimensionless parameters assumed Parameters calculated = 1.5 / 3 = 6.0 × 10 6 / 2 = 0.

Analysis of 3DOF SSI System
The system consists of a mass-spring-damper combination representing the dynamic properties of a structure, mounted on a rigid column of height 'h'. This column is supported by a roller which represents a rigid ground in vertical direction. The horizontal translational and rotational degrees of freedom at the bottom of the column structure are governed by the corresponding components of foundation elasticity and damping. For the analysis purpose a foundation consisting of a rigid circular disk on an elastic half-space is considered. From the Figure 6, it can be noted that there are three degrees of freedom for the whole system-horizontal degree of freedom of the structure mass 'm', the horizontal and rotational degrees of freedom of the foundation. The input values of the parameters and the expressions used to calculate the derived input values are shown in Table  3. Here, is the stiffness of the structure spring, is the structural damping ratio, is the structure mass, is the radius of the rigid circular disk foundation, , , , are the shear moduli, mass density, Poisson ratio, damping ratio of the soil respectively. Figure 6. Three-DOF dynamic system representing SSI The stiffness and damping values of foundation-soil system are obtained for a rigid circular disk of radius 'a', resting on a ground/ soil having a shear modulus of 'G', density of 'ρ', Poisson ratio of 'ν', material damping constant of 'ςg' This procedure is described in (Wolf, 1985) and the coefficients are shown in Figure 7 for horizontal and rocking components.
High frequency behaviour of rigid circular base mat resting on the surface of a layered half-space are (Wolf, 1988) : , This forms the singular part of the dynamic stiffness coefficient. The regular part of the dynamic stiffness coefficient of a rigid circular disk or base mat resting on an undamped half-space can be obtained by separating the singular part from the dynamic stiffness coefficients shown in Figure  7.
From the singular part of dynamic stiffness coefficients the initial value (i.e., at t = 0+) of the coefficient in time domain can be obtained. The regular part of dynamic stiffness coefficient is transformed to time domain using numerical Fourier transform to get , ( ). The system can then be modelled in ABAQUS into a three DOF system shown in Figure 8. It is subjected to artificial earthquake ground motion described in the previous section ( Figure 5). The results are shown in Figure 9 and Figure 10. Figure 9 gives the displacement history at base and structure mass, as well as the rotation at the base. It can be seen that the displacement response is considerably reduced with respect to the (free-field) ground motion in Figure 5. This is due to the high stiffness of soil. Further, the displacement of structural mass is quite magnified from that at the base. This is bound to happen for the selected ̅ = 0.1, which means that the structure is quite flexible compared to the foundation stiffness. It can also be seen that the rotation of the base is quite negligible.  The internal force in structural connector is shown in Figure 10 (b) and the applied force at the base in shown in Figure 10 (a). The applied force is the interaction force calculated using the convolution integral through the UAMP subroutine. The structure-connector being quite flexible in relation to the foundation experiences very low internal force as can be seen in Figure 10.

Conclusions
The paper successfully presents steps for carrying out seismic analysis of three different SSI systems in time domain using ABAQUS software, wherein the user defined subroutine UAMP is programmed to interactively evaluate the soil structure interaction forces that ensure realistic transmitting boundaries. The procedure presented in the paper can be applied to continuum SSI models also.