Iterative coupling procedure applied to the transient response of dynamic soil-structure systems

This paper proposes a methodology to obtain the transient response of structural system interacting with soil-foundation schemes supported by viscoelastic soils. The structure and soil are divided into sub-systems. The time domain solution for each subsystem is formulated by an appropriated methodology. The equations of motion of structure are solved by Newmark integration algorithm. The transient response of the soil is obtained by a convolution integral. The convolution integral uses transient impulse response of viscoelastic soils. Newmark and convolution algorithms are formulated as input and output schemes, which, in turn, are plugged to the time stepping iterative algorithm. The scheme is applied to vertical response of a dynamical system interacting with a massless foundation laying on a soil modelled as a three-dimensional homogeneous viscoelastic half-space. For two distinct external forces, the resulting coupled displacements, interface forces, errors and number of iterations within each time step are provided.


INTRODUCTION
The dynamic analysis of complex systems with parts presenting distinct properties constitutes an immense challenge to the scientific and professional community. Dynamic soil-structure interaction problems (DSSI) represent one such class of complex systems. The structures are bounded, heterogeneous, anisotropic and frequently present local non-linear behavior. The soil is usually less stiff, more homogeneous and usually is unbounded in two or three directions. The soil unboundedness introduces a geometric damping due to wave propagating into the infinite region and carrying energy with them.
Due to the complexity of constitutive equations, geometries, boundary conditions and applied loads it is not possible to tackle these problems with analytical tools. It is also very difficult to solve DSSI problems with a single numerical methodology. Presently, a very efficient methodology to describe DSSI problems is the coupling of the Finite Element Method (FEM) with the Boundary Element Method (BEM). The FEM can model, in a standard way, the heterogeneity and the non-linear behavior in localized parts of the structure. The FEM can also model the soil at the vicinity of the structure, where non-linear deformation or non-linear contact behavior may be expected (Hughes, 1987). On the other hand, the BEM can describe accurately the linear wave propagation process in the soil and the radiation condition. Thus, a coupling FEM-BEM has been an efficient and frequently used technique in the DSSI analysis (Von Estorff and Prabucki, 1990).
Dynamic analysis can be performed in frequency or in time domain. The FEM has become the standard numerical tool in linear and non-linear dynamic analysis (Bathe, 1996;Crisfield, 1991). A large set of very efficient numerical solvers are available to the FEM researcher or practitioner, that make full use of the sparsity and symmetry of the resulting FE system matrices. The BEM analysis for DSSI in the frequency domain was formulated by Dominguez (1978). An important feature of the BEM is that, when formulated based on a proper fundamental solution, the resulting Integral Equation only requires a boundary discretization, reducing thus the dimensionality of the problem by one. On the other hand, the resulting system matrices are fully populated and non-symmetric. A direct time domain formulation and implementation was first given by Mansur and Brebbia (1982a,b). Many contributions followed to improve the BEM in time domain (Von Estorff and Antes, 1991;Coda and Venturini, 1995;Rizos and Karabalis, 1998;Rizos and Wang, 2002). All these formulations represented distinct implementations of the so-called Stoke's time fundamental solution for linear elastic domains.
Transient response of unbounded soil problems and the interaction with supported structures have also been obtained with the use of Laplace transform domains (Wang and Rajapakse, 1993). A hybrid time-frequency domain procedure for transient BEM has been proposed by Gaul et al. (1992). A transient BEM formulation based on the Convolution Quadrature Method -CQM was proposed by Schanz and Antes (1997). A general formulation and implementation to transient linear viscoelatic models for FEs and BEs was presented by Mesquita et al. (2001). It requires a domain integration and cannot be directly applied to unbounded domains. A non-singular transient influence functions for semi-unbounded domains based on viscoelastic frequency domain solutions and the application of the FFT algorithm was synthesized by Mesquita et al. (2012). These solutions, which do not present singularities at the loading point (Barros and Mesquita, 1999) and require no integration regularization procedures (Dangla et al., 2005), were able to render accurate wave propagation phenomena for general linear viscoelastic models.
There are distinct strategies to couple FE and BE formulation in direct time domain. The first one is the direct coupling formulated by Von Estorff and Prabucki (1990). It is efficient but the integration time steps for both BE and FE domains are the same. This may result in instabilities due to the different stiffness and, consequently, distinct wave propagation speeds in both domains. Recent studies, nevertheless, stressed the efficiency of the direct coupling method face to face to the other methodologies (François et al., 2015). One possibility of the direct coupling methodology is to transform the BE into FE-like elements at the interfaces. This has been applied to determine dispersion relations in wave guides embedded in semi-unbounded media (Mazzotti et al., 2013). Coupled FE-BE systems with localized non-linearities have been presented by Coda (2001).
The second possibility is the so-called staggered coupling (Rizos and Wang, 2002). The technique has the advantage that each subdomain may be solved by the best available technique. In this procedure, for every time step, the interface solutions of one domain is fed as input to the second domain. The procedure implies in the same time step for the integration of both domains and require small integration steps to maintain stability of the solution (François et al., 2015).
The third FE-BE coupling mechanism is the iterative coupling. The first time domain FE-BE iterative coupling strategy was presented by Soares et al. (2004) and has been subjected to systematic improvements. An iterative coupling strategy shows many advantages over the other techniques. As in the staggered approach, every subdomain may be solved by the most suitable technique. The time steps for each domain may not be the same, potentializing the stability of the method (Soares, 2008;François et al., 2015). The iterative coupling strategy may incorporate a relaxation parameter which may accelerate the convergence of the iterative method (Soares, 2008). The iterative approach is, in principle, also prone to handle non-linear contact problems, such as foundation partial uplift (Wolf, 1988). A fourth possibility to coupled the FE and BE subdomains is through the Lagrange multipliers strategy (González et al., 2007). Recently, Soares and Araújo (2019) proposed an advanced coupling procedure that does not require the iterative convergence process.
All the previous works on iterative coupling used direct BEM time formulations for elastic domains. There is up to the present time, to the best of the authors' knowledge, no general time domain fundamental solution for linear viscoelastic unbounded domains. In this article an iterative coupling strategy is presented that allows linear or non-linear structures to be coupled directly in time domain with time response of rigid structures interacting with viscoelastic soils.
The time domain solution for both FE and BE methodologies are formulated in the so-called Newman-Dirichlet and Dirichlet-Newman modes (François et al., 2015). In this paper these solution schemes are called mode 1 and mode 2, respectively. For the FE domain, the Newmark integration algorithm is formulated to render mode 1 and mode 2 schemes for the iterative process.
An indirect version of the BEM is applied to synthesize the frequency response of a circular rigid and massless foundation interacting with a viscoelastic half-space (soil model) (Barros and Mesquita, 1999;Labaki et al., 2014). An accurate integration scheme allows to determine the frequency response for the soil-foundation at very high frequencies (Mesquita et al., 2012). This high frequency solution together with the FFT algorithm yield very accurate time response for the soil system, in which linear viscoelastic effects can be incorporated.
The time response of the BE domain to a general excitation is obtained by a superposition convolution integral, relating the time impulse response of the soil-foundation system to the time loads at the interface in every time step. The convolution integral has also been casted into mode 1 and mode 2 schemes. At every time step, the iterative coupling algorithm updates the interface nodal values, alternating in mode 1 and mode 2, until a prescribed convergence is reached. A relaxation parameter α is introduced for the interface variables update in mode 2, but in this study its value is kept constant α = 0.5. In the present article only the vertical response of the soil and of the structure were considered. So only the vertical structure and foundation degree freedom at the soil-structure interface is being coupled. The extension for the coupling FE interface variables to the multiple degrees of freedom of a rigid foundation can be obtained by standard procedures (Rizos, 2000).
The main accomplishment of the present article is related to the formulation of the iterative scheme for the viscoelastic unbounded domain, the BE discretized domain. A previously synthesized time domain solution for rigid foundation interacting with a linear viscoelastic unbounded domain is incorporated into a discretized convolution integral, which in turn, can, for every time step of the iterative procedure, deliver Newman or Dirichlet values of the variables at the domains interface. To the best of the authors' knowledge, iterative time domain coupling of unbounded viscoelastic domains with structures has not been reported in the literature.

STATEMENT OF PROBLEM
The problem being analyzed in this article is shown in Figure 1. This figure describes structure modelled as a linear system with mass M, stiffness K and damping coefficient C, interacting with a soil, described in this case as a homogeneous half-space, characterized by a density ρ, shear módulus G, Poisson ratio ν and internal damping ratio η. The structure is excited by an external load F(t). In the present case, the soil response is characterized by the vertical displacement at the geometric center of a massless and rigid foundation laying on a homogeneous half-space (Labaki et al., 2014).
The system described is subdivided into two sub-systems, the structure and the soil. The forces and the displacements at the soil-structure interface are designated by F1(t), u1(t), F2(t) and u2(t), respectively, as shown in Figure 1.
Fully bonded coupling conditions of the sub-systems require that at every time step, displacement kinematic compatibility and force equilibrium at the soil-foundation interface expressed, respectively, as: Iterative coupling procedure applied to the transient response of dynamic soil-structure systems In the present article, the soil response u2(t) is the vertical displacement at the geometric center of a massless and rigid foundation interacting with a homogeneous half space (Labaki et al., 2014). But the soil profile may be a layered one (Labaki et al., 2014), the foundation may be embedded in the soil (Carrion et al., 2007) or the response of a pile embedded in the soil may also be considered (Barros et al., 2019;Lima et al., 2019).

Input-output responses
To obtain the transient response of the coupled system, the equilibrium and kinematic compatibility equation must be fulfilled at every time step (ti) of the response. It is assumed that, for every time step t = ti , the response of each subsystem may be formulated as an input-output 'box'. In this article the equations of motion of the structure will be integrated by the Newmark time stepping algorithm (Damasceno, 2013). As detailed in the next session, the Newmark algorithm can be formulated in the form of a Newman-Dirichlet (force-displacement), or Dirichlet-Newman (displacement-force) scheme to be applied at the iterative coupling methodology (François et al., 2015). These two solution schemes are termed mode 1 and mode 2 in the present article and are illustrated in Figure 2.
The transient soil response will be obtained by means of a time convolution integral, relating the soil response due to unit (time) concentrated impulse to the actual loading function, in this case, the interface force F2(ti) (Mesquita et al., 2012). Both integration procedures will be detailed in the next session. In this session, for the purpose of explaining the iterative coupling procedure, it suffices to say that for both subsystems it is possible to formulate the response in the modes shown in Figure 2. Figure 2 should be read in the following way. For every time step (ti), the forces acting upon the structure are the external excitation F(ti) and the interface force F1(ti). In this case the Newmark algorithm can render the displacement u1(ti) that fulfills the dynamic equilibrium equations of motion . For the soil, the convolution integral may be so formulated that given a force applied at the soil interface F2(ti), the convolution algorithm delivers the corresponding displacement u2(ti) Tovo, 2018;Tovo et al., 2019) that corresponds to the solution of the soil differential equations of motion (Mesquita et al., 2012). This operation is called Mode 1.
Latin American Journal of Solids and Structures, 2020, 17(Thematic Section 8 MecSol), e309 5/13 To properly formulate the iterative algorithm, the inverse operation, in which, for every time step the interface displacements are the (Dirichlet) input quantities, u1(ti) and u2(ti), and the output are the forces F1(ti) and F2(ti) (Newman quantities) that satisfies the equations of motion of both subdomains. This operation is called Mode 2.

The coupling algorithm
With the numerical input-output algorithms established, the coupling algorithm, shown in Figure 3, can be explained. The first action is, for a given time step (ti), to solve both iterative algorithms in Mode 1, for the case that all interface forces are given an initial value. In the present case this value is zero, F1(ti) = F2(ti) = 0. The outputs are the displacements for each subsystem u1(ti , j) and u2(ti , j). The internal iterative loop within every time step (ti) is designated by the index 'j'.
Next, the resulting displacements are compared within a pre-determined error erru level to check if the continuity at the soil-foundation interface is satisfied, u1(ti , j) -u2(ti , j) < erru. If this is the case, the coupling has been attained. Both interface displacements are equal, within an error bound erru, satisfying Equation (1) and both interface forces satisfy Equation (2), in this case with F1(ti) = F2(ti) = 0.
Two remarks should be made about this algorithm. First, throughout this article the relaxation or weighting factors are kept constant with the value α = αu = αf = 0.5. This means that we are averaging displacements and forces within the j-cycles. Second, convergence is only tested for interface displacements and not for forces. The literature reports techniques to optimize the relaxation parameter, but this has not been incorporated in the present analysis (Soares, 2008).

FORMULATION
In this session the numerical algorithms used to render the transient response of the subsystems shown in the Problem Statement session will be given a formulation which, in turn, will render them amenable to be incorporated in the iterative scheme presented in Figure 3.

Formulation of the Input-Output Algorithm for the Structure Transient Response
The equation of motion for the subsystem 1, the structure, is given by: The Newmark algorithm may be written in the input-output Mode 1 and Mode 2, as indicated in Figure 2.

Newmark Mode 1
For the Mode 1, considering the time instant (ti) and the inner iteration 'j', having as input the external force F(ti) and the interface iterative force F1(ti , j), the output, i.e. the displacement u1(ti , j+1) at the inner iteration (j+1), may be expressed using the Newmark algorithm as (Damasceno, 2013): The Newmark integration parameters are α = 0.25 and β = 0.5.

Newmark Mode 2
The Newmark integrator may also be written in terms of Mode 2, in which the input is the displacement u1(ti , j+1) and the output is the interface force F1(ti , j+1) satisfying the equilibrium equation of motion (Damasceno, 2013):

Formulation of the Input-Output Algorithm for the Soil Transient Response
The soil transient displacement response u2(ti), is obtained by the convolution integral relating the actual force applied at the soil foundation interface F2(ti) with the soil impulse response function h2 δ (ti) (Wolf, 1988;Mesquita et al., 2012):

Convolution Integral Mode 1
For the case of Mode 1, considering the time instant (ti) and the inner iteration 'j', having as input the interface iterative force F2(ti , j), the output, i.e. the displacement u2(ti , j+1) at the inner iteration (j+1) may be expressed in discretized as (Tovo, 2018):

Convolution Integral Mode 2
Equation (11) can easily be reformulated to have, as an input, the displacement u2(ti , j+1) and, as output, the soil interface F2(ti , j+1) at the next inner iteration cycle.
In this article, both sub-systems are discretized with the same time step. Nevertheless, a linear and quadratic interpolation of the variables in the convolution integral have been implemented by Tovo (2018) and can be used to expand the iterative algorithm to use distinct time steps in each sub-system (François et al., 2015).

Frequency Domain Response
In this session, the strategy used to obtain the transient soil impulse response is described. The soil model is given by a rigid massless rigid foundation resting upon a viscoelastic homogeneous half-space. The frequency domain solution for this problem has been given by Labaki et al. (2014). The frequency domain vertical displacement of the foundation geometric center H2(ω) due to a unitary load F2(ω) = 1(ω) can be seen in Figure 4a. The parameters used in these calculations are: G = 1.0 [Pa], ρ = 1.0 [kg/m 3 ], ν = 0.25, η = 0.01. It is costumary to use unit constitutive values to test numerical methodologies (Soares, 2008).

Transient Impulse Response
The procedure developed to determine the frequency domain solution shown in Figure 4a can be used to calculate this response up to very high frequencies (Mesquita et al., 2012). This means that the FFT algorithm can be applied to the frequency domain response H2(ω) to render the transient soil impulse response h2 δ (t) with very small time steps shown in Figure 4b.
Two remarks should be mentioned about this topic. The transient impulse response h2 δ (t) shown in Figure 4b will be used in the iterative coupling algorithm, both in Mode 1, Equation (11), and Mode 2, Equation (12). In this article, the chosen soil model is a homogeneous half-space, but many other soils or foundation supporting systems like layered soils (Labaki et al., 2014), embedded foundations (Carrion et al., 2007) and pile supporting schemes (Barros et al., 2019;Lima et al., 2019) may be used.

Formulation of the Frequency Domain Solution of the Coupled Soil-structure System
In order to have a solution that allows a validation of the proposed iterative coupling scheme, a transient solution of the complete, direct coupled, soil-structure system is obtained. The idea is to use the frequency response solution of the structure H1(ω) and the soil frequency response H2(ω) and to couple both solutions through force equilibrium and kinematic continuity condition, given in Equations (1) and (2), to obtain the frequency domain solution of the coupled soil-foundation system HT(ω). The structure frequency response H1(ω) corresponds to a one degree of freedom oscillator with mass M, damping C and stiffness K and can be obtained in a rather straightforward manner. The more intricate soil solution H2(ω) is given, exemplarily, in Figure 4a. The expression for coupled solution HT(ω) can be shown to be:  (13) will furnish the transient unit impulse response of the coupled system hT δ (t), which is given in Figure 5b. This response hT δ (t) will be used in conjunction with the convolution integral to render the transient response of the direct coupled soil-structure system: Figure 5: (a) Frequency domain response of the coupled soil-structure system HT(ω). (b) Transient unit impulse response of the coupled soil-structure system hT δ (t).

NUMERICAL EXAMPLES
In this session, the responses obtained by the iterative coupling procedure are compared to those obtained from the direct coupled system, given in Equations (13) and (14). Two external excitation forces F(t) are used as examples, namely, trapezoidal and triangular excitation. The values of the absolute and relative displacement errors between both approaches are given. The resulting interface forces F1(t) and F2(t) are also shown. The number of iterative steps 'j' within each time step Δt is also presented. Two other error definitions are used in this article. They are, respectively, the absolute, errabs, and relative, errrel, errors between the subsystems u1(t) = u2(t) = u1,2(t) and the response of the coupled system uT(t):

Trapezoidal excitation
Figures 6 to 8 present the performance of the iterative coupling procedure for the case of a trapezoidal external excitation F(t), shown in Figure 6a. The displacement response for the both solutions, the initially coupled system uT(t) and the iteratively coupled systems u1(t) = u2(t), is given in Figure 6b. The errors as defined in Equations (15) and (16) are shown in Figures 7a and 7b, respectively. The resulting interface forces obtained in the iterative process for every time step is given in Figure 8a. The resulting interface forces are consistent with the excitation profile.
The number of inner iterations 'j' required to obtain convergence between the displacements of the two subsystems u1(ti , j) -u2(ti , j) < erru is shown in Figure 8b. This Figure 8b indicates that at some time steps the maximum number of iterations per time step, nmax = 200, is reached without displacement convergence within the prescribed error bound of 10 -12 . It should be mentioned that this is a very high convergence requirement. Even when the iterated displacement solution u1(ti , j) -u2(ti , j) < erru fulfill the convergence criteria, this iterated solution still presents larger absolute and relative errors when compared to the originally coupled system displacement uT(t).

Triangular excitation
The iterative coupling procedure is tested now for a triangular external excitation, shown in Figure 9a. The results for the displacements, errors, interface forces and number of iterations within each time step are given in Figures 9b through 11b. They follow basically the same pattern of the results obtained in the previous example. The idea is to show that, at least what concerns the external excitation source, the proposed methodology is consistent.   The results presented in Figures 6 to 11 are the main outcomes of the methodology proposed and implemented in this article. These initial results tend to indicate that the proposed transient iterative coupling procedure has the potential to obtain transient responses of more complex structural and foundation systems interacting with many soil profiles or soil-supported foundation schemes. It may be applied to obtain transient solutions of coupled structures with viscoelastic soils.

CONCLUSION
This article presented an iterative methodology to obtain the coupled transient response of structural systems interacting with viscoelatic soil-foundation arrangements. Structure and soil are considered two distinct subsystems. For the structure, an input-output scheme, relating Newman and Dirichlet interface quantities, is developed based on the Newmark time integration algorithm. The transient response of the viscoelastic soil-foundation sub-system is obtained by a discretized convolution integral relating the soil transient impulse response to the actual interface forces and displacements. The discretized convolution integral has also been formulated as input-output procedures relating Newman and Dirichlet quantities at the soil-structure interface. The kernel of the convolution integral for the soil is the transient impulse response of the viscoelastic soil model. It is obtained applying the FFT algorithm to an accurately synthesized frequency domain solution of the soil problem. The proposed iterative strategy couples the Newmark algorithm with a discretized version of the convolution integral.
The method is applied to obtain the vertical response of a mechanical structure interacting with a rigid and massless foundation resting on a viscoelastic half-space. The results obtained by the iterative coupling procedure presented a good agreement with the solution obtained for the complete, directly coupled, soil-structure system. The interface displacement, forces and the number of iterations within each step to achieve displacement convergence have been determined. Relative and absolute error measures were provided. The methodology may be further developed to consider structures with localized non-linear behavior, non-linear contact conditions at the soil-structure interface, including partial uplift, and also distinct integration time steps for each subdomain. It may also contribute to advance the challenging problem of obtaining transient responses for viscoelastic unbounded domains.