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

Abstract

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.

Keywords:
Dynamic soil structure interaction; viscoelastic soils; time-domain iterative coupling procedure; convolution integral

Graphical Abstract

1 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.

2 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 F₁(t), u₁(t), F₂(t) and u₂(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:

$u_1\left(t_i\right)=u_2\left(t_i\right)$ (1)

and

$F_1\left(t_i\right)+F_2\left(t_i\right)=0$ (2)

In the present article, the soil response u₂(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).

Figure 1:
Statement of problem.

3 COUPLING PROCEDURE

3.1 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 (t_i) of the response. It is assumed that, for every time step t = t_i , 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 F₂(t_i) (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:
Input-Output Modes: Newman-Dirichlet (Mode 1), Dirichlet-Newman (Mode 2).

Figure 2 should be read in the following way. For every time step (t_i), the forces acting upon the structure are the external excitation F(t_i) and the interface force F₁(t_i). In this case the Newmark algorithm can render the displacement u₁(t_i) that fulfills the dynamic equilibrium equations of motion (Damasceno et al., 2013). For the soil, the convolution integral may be so formulated that given a force applied at the soil interface F₂(t_i), the convolution algorithm delivers the corresponding displacement u₂(t_i) (Damasceno et al., 2013; 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.

To properly formulate the iterative algorithm, the inverse operation, in which, for every time step the interface displacements are the (Dirichlet) input quantities, u₁(t_i) and u₂(t_i), and the output are the forces F₁(t_i) and F₂(t_i) (Newman quantities) that satisfies the equations of motion of both subdomains. This operation is called Mode 2.

3.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 (t_i), 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, F₁(t_i) = F₂(t_i) = 0. The outputs are the displacements for each subsystem u₁(t_i , j) and u₂(t_i , j). The internal iterative loop within every time step (t_i) is designated by the index ‘j’.

Next, the resulting displacements are compared within a pre-determined error err_u level to check if the continuity at the soil-foundation interface is satisfied, u₁(t_i , j) - u₂(t_i , j) < err_u. If this is the case, the coupling has been attained. Both interface displacements are equal, within an error bound err_u, satisfying Equation (1) and both interface forces satisfy Equation (2), in this case with F₁(t_i) = F₂(t_i) = 0.

In general Equations (1) and (2) are not fulfilled at the first iteration step. In this case, new displacement for the iterative inner counter j=j+1, u₁(t_i , j+1) and u₂(t_i , j+1) are determined from weighting the previous displacement results, u₁(t_i , j) and u₂(t_i , j):

$u_1\left(t_i,j+1\right)=u_2\left(t_i,j+1\right)=\left[{\alpha }_u{u}_1\left(t_i,j\right)+\left(1-{\alpha }_u\right)u_2\left(t_i,j\right)\right]$ (3)

These new displacements u₁(t_i , j+1) and u₂(t_i , j+1) are now fed into the Mode 2 solvers, having as an output the new interface forces F₁(t_i , j+1) and F₂(t_i , j+1). To start a new iterative cycle (j+2), the forces F₁(t_i , j+1) and F₂(t_i , j+1) are weighted to obtain the input forces F₁(t_i , j+2) and F₂(t_i , j+2):

$F_1\left(t_i,j+2\right)=F_2\left(t_i,j+2\right)=\left[{\alpha }_f{F}_1\left(t_i,j+1\right)+\left(1-{\alpha }_f\right)F_2\left(t_i,j+1\right)\right]$ (4)

These new forces, F₁(t_i , j+2) and F₂(t_i , j+2), are now fed into the Mode 1 algorithms, delivering the displacements u₁(t_i , j+2) and u₂(t_i , j+2), which are compared within the limits of an established error bound u₁(t_i , j+2) - u₂(t_i , j+2) < err_u. Again, if this condition is fulfilled the subsystems are considered coupled within this numerical error bound. Otherwise, the cycle continues until at step ‘n’ the displacement compatibility condition u₁(t_i , n) - u₂(t_i , n) < err_u is satisfied, or the inner counter ‘j ‘ reaches a predefined maximum value, ‘j=n_max’, and the procedure stops without having reached the required convergence.

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).

4 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.

4.1 Formulation of the Input-Output Algorithm for the Structure Transient Response

The equation of motion for the subsystem 1, the structure, is given by:

$M_1\cdot {\ddot{u}}_1+C_1\cdot {\dot{u}}_1+K_1\cdot u_1=F+F_1$ (5)

The Newmark algorithm may be written in the input-output Mode 1 and Mode 2, as indicated in Figure 2.

Figure 3:
The iterative coupling algorithm.

4.1.1 Newmark Mode 1

For the Mode 1, considering the time instant (t_i) and the inner iteration ‘j’, having as input the external force F(t_i) and the interface iterative force F₁(t_i , j), the output, i.e. the displacement u₁(t_i , j+1) at the inner iteration (j+1), may be expressed using the Newmark algorithm as (Damasceno, 2013):

$u_1\left(t_i,j+1\right)={\left[\frac{1}{\alpha \Delta t^2}M+\frac{\beta }{\alpha \Delta t}C+K\right]}^{-1}\cdot \left\{F\left(t_i\right)+F_1\left(t_i,j\right)+A+B\right\}$ (6)

with

$A=M\left[\frac{1}{\alpha \Delta t^2}{u}_1\left(t_i,j\right)+\frac{1}{\alpha \Delta t}{\dot{u}}_1\left(t_i,j\right)+\left(\frac{1}{2\alpha }-1\right){\ddot{u}}_1\left(t_i,j\right)\right]$ (7)

$B=C\left[\frac{\beta }{\alpha \Delta t}{u}_1\left(t_i,j\right)+\left(\frac{\beta }{\alpha }-1\right){\dot{u}}_1\left(t_i,j\right)+\left(\frac{\beta }{\alpha }-2\right)\frac{\Delta t}{2}{\ddot{u}}_1\left(t_i,j\right)\right]$ (8)

The Newmark integration parameters are α = 0.25 and β = 0.5.

4.1.2 Newmark Mode 2

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

$F_1\left(t_i,j+1\right)=M\cdot {\ddot{u}}_1\left(t_i,j+1\right)+C\cdot {\dot{u}}_1\left(t_i,j+1\right)+K\cdot u_1\left(t_i,j+1\right)-F\left(t_i\right)$ (9)

4.2 Formulation of the Input-Output Algorithm for the Soil Transient Response

The soil transient displacement response u₂(t_i), is obtained by the convolution integral relating the actual force applied at the soil foundation interface F₂(t_i) with the soil impulse response function h₂ ^δ(t_i) (Wolf, 1988; Mesquita et al., 2012):

$u_2\left(t_i\right)={\displaystyle \underset{\tau =0}{\overset{\tau =t_i}{\int }}{F}_2\left(\tau \right)\cdot h_2^{\delta }\left(t_i-\tau \right)}d\tau$ (10)

4.2.1 Convolution Integral Mode 1

For the case of Mode 1, considering the time instant (t_i) and the inner iteration ‘j’, having as input the interface iterative force F₂(t_i , j), the output, i.e. the displacement u₂(t_i , j+1) at the inner iteration (j+1) may be expressed in discretized as (Tovo, 2018):

$u_2\left(t_i,j+1\right)={\displaystyle \sum _{k=0}^{k=i-1}{F}_2\left(t_k\right)h_2^{\delta }}\left(t_{i-k}\right)\Delta t+F_2\left(t_i,j\right)h_2^{\delta }\left(0\right)\Delta t$ (11)

4.2.2 Convolution Integral Mode 2

Equation (11) can easily be reformulated to have, as an input, the displacement u₂(t_i , j+1) and, as output, the soil interface F₂(t_i , j+1) at the next inner iteration cycle.

$F_2\left(t_i,j+1\right)=\frac{1}{h_2^{\delta }\left(0\right)\Delta t}\left[u_2\left(t_i,j+1\right)-{\displaystyle \sum _{k=0}^{k=i-1}{F}_2\left(t_k\right)h_2^{\delta }\left(t_{i-k}\right)\Delta t}\right]$ (12)

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).

4.3 The Transient Soil Impulse Response

4.3.1 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 H₂(ω) due to a unitary load F₂(ω) = 1(ω) can be seen in Figure 4a. The parameters used in these calculations are: G = 1.0 [Pa], ρ = 1.0 [kg/m³], ν = 0.25, η = 0.01. It is costumary to use unit constitutive values to test numerical methodologies (Soares, 2008).

4.3.2 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 H₂(ω) to render the transient soil impulse response h₂ ^δ(t) with very small time steps shown in Figure 4b.

Two remarks should be mentioned about this topic. The transient impulse response h₂ ^δ(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.

Figure 4:
(a) Soil frequency domain displacement solution H₂(ω) for a unitary load F₂(ω) = 1 (ω). (b) Soil transient impulse response h₂ ^δ(t) due to a unit impulse excitation f(t)=δ(t).

4.4 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 H₁(ω) and the soil frequency response H₂(ω) 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 H_T(ω). The structure frequency response H₁(ω) 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 H₂(ω) is given, exemplarily, in Figure 4a. The expression for coupled solution H_T(ω) can be shown to be:

$H_T\left(\omega \right)=\left[\frac{H_1\left(\omega \right)\cdot H_2\left(\omega \right)}{H_1\left(\omega \right)+H_2\left(\omega \right)}\right]$ (13)

Figure 5a shows the frequency domain response of the direct coupled soil-structure system for the parameters M = 0.4 [kg], C = 0.01 [Ns/m], K = 0.01 [N/m], G = 1.0 [Pa], ρ = 1.0 [kg/m³], ν = 0.25, η = 0.01. Using the FFT algorithm on Equation (13) will furnish the transient unit impulse response of the coupled system h_T ^δ(t), which is given in Figure 5b. This response h_T ^δ(t) will be used in conjunction with the convolution integral to render the transient response of the direct coupled soil-structure system:

$u_T\left(t_i\right)={\displaystyle \underset{\tau =0}{\overset{\tau =t_i}{\int }}F\left(\tau \right)\cdot h_T^{\delta }\left(t_i-\tau \right)}d\tau$ (14)

Figure 5:
(a) Frequency domain response of the coupled soil-structure system H_T(ω). (b) Transient unit impulse response of the coupled soil-structure system h_T ^δ(t).

5 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 F₁(t) and F₂(t) are also shown. The number of iterative steps ‘j’ within each time step Δt is also presented. Throughout these examples the system parameters are: M = 0.4 [kg], C = 0.01 [Ns/m], K = 0.01 [N/m], G = 1.0 [Pa], ρ = 1.0 [kg/m³], ν = 0.25, η = 0.01. The value of the convergence criteria for the displacements of both subsystems u₁(t_i , j) - u₂(t_i , j) < err_u in the iterative process is established to be err_u = 10^-12. The number of maximum iterative steps within a time step is set to n_max = 200. The displacement u units are given in meters [m].

Two other error definitions are used in this article. They are, respectively, the absolute, err_abs, and relative, err_rel, errors between the subsystems u₁(t) = u₂(t) = u_1,2(t) and the response of the coupled system u_T(t):

$er{r}_{abs}=\left|u_T-u_{\mathrm{1,2}}\right|$ (15)

$er{r}_{rel}=\left|\frac{u_T-u_{\mathrm{1,2}}}{u_T}\right|$ (16)

5.1 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 u_T(t) and the iteratively coupled systems u₁(t) = u₂(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 u₁(t_i , j) - u₂(t_i , j) < err_u is shown in Figure 8b. This Figure 8b indicates that at some time steps the maximum number of iterations per time step, n_max = 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 u₁(t_i , j) - u₂(t_i , j) < err_u fulfill the convergence criteria, this iterated solution still presents larger absolute and relative errors when compared to the originally coupled system displacement u_T(t).

Figure 6:
(a) Trapezoidal excitation - F(t) [N]. (b) Displacements: u₁ , u₂ , u_T [m].

Figure 7:
Displacement errors u_1,2(t)- u_T(t) : (a) Absolute error. (b) Relative error.

Figure 8:
(a) Coupling interface forces F₁(t), F₂(t). (b) Iterations to convergence.

5.2 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.

Figure 9:
(a) Triangular excitation - F(t) [N]. (b) Displacements: u₁ , u₂ , u_T [m].

Figure 10:
Displacement errors u_1,2(t)- u_T(t) : (a) Absolute error. (b) Relative error.

Figure 11:
(a) Coupling interface forces F₁(t), F₂(t). (b) Iterations to convergence.

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.

6 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.

ACKNOWLEDGEMENTS

The research leading to this work was funded by the São Paulo Research Foudation (Fapesp) through grant FAPESP CEPID Process 2013/08293-7 and Maranhão Research Foundation (Fapema). The support of Capes, CNPq and Faepex/Unicamp is also gratefully acknowledged. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

References

• Barros, P.L.A., Mesquita, E., (1999). Elastodynamic Green’s Functions for Orthotropic Plane Strain Continua with inclined Axis of Symmetry. International Journal of Solids and Structures 36:4767-4788.
• Barros, P.L.A., Labaki, J., Mesquita, E., (2019). IBEM-FEM model of a piled plate within a transversely isotropic half-space. Engineering Analysis with Boundary Elements 101:281-296.
• Bathe, K.J. (1996). Finite Element Procedures. Englewood Cliffs, NJ: Prentice-Hall.
• Carrion, R., Amilcar, D.O.S., Mesquita, E., (2007). The influence of soil damping mechanisms and geo-profiles on the stationary response of rigid block foundations. In Proceedings of the 19th Cobem 01-10.
• Coda, H.B., Venturini, W.S., (1995). Non-singular time-stepping BEM for transient elastodynamic analysis. Engineering Analysis with Boundary Elements 15:11-18.
• Coda, H.B. (2001). Dynamic and static non-linear analysis of reinforced media: a BEM/FEM coupling approach. Computers & Structures 79:2751-2765.
• Crisfield, M.A. (1991). Non-linear Finite Element Analysis of Solids and Structures. Chichester: John Wiley & Sons. Vols. 1 & 2.
• Damasceno, D.A. (2013). Transient analysis of systems with soil-structure interaction through iterative coupling techniques. Master Thesis (in Portuguese), University of Campinas, Brazil.
• Damasceno, D.A., Labaki, J., Mesquita, E., (2013). Dynamic transient analysis of partitioned systems through an iterative coupling technique. In Proceedings of the XXXIV Iberian Latin-American Congress on Computational Methods in Engineering 01-20.
• Dangla, P., Semblat, J.-F., Xiao, H., Delépine, N., (2005). A Simple and Efficient Regularization Method for 3D BEM: Application to Frequency-Domain Elastodynamics. Bulletin of the Seismological Society of America 95:1916-1927.
• Dominguez, J. (1978). Dynamic Stiffness of Rectangular Foundations. MIT Research Report R78-20. NSF-RANN Grant No. ENV 77-18339.
• François, S., Coulier, P., Degrande, G., (2015). Finite element-boundary element coupling algorithms for transient elastodynamics. Engineering Analysis with Boundary Elements 55:104-121.
• Gaul, L., Schanz, M., Fiedler, C., (1992). Viscoelastic formulations of BEM in time and frequency domain. Engineering Analysis with Boundary Elements 10:137-141.
• González, J.A., Park, K.C., Felippa, C.A., (2007). FEM and BEM coupling in elastostatics using localized Lagrange multipliers. International Journal for Numerical Methods in Engineering 69:2058-2074.
• Hughes, T.J.R. (1987). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs, NJ: Prentice-Hall.
• Labaki, J., Mesquita, E., Rajapakse, R.K.N.D., (2014). Vertical vibrations of an elastic foundation with arbitrary embedment within a transversely isotropic, layered soil. Computer Modeling in Engineering & Sciences 103:281-313.
• Lima, L.F.V., Labaki, J., Mesquita, E., (2019). Studies on the Stationary Dynamic Response of a Foundation Supported by Flexible Pile and Soil to a Vertical Incident Wave Field or External Force. In Proceedings of the XL Iberian Latin-American Congress on Computational Methods in Engineering 01-11.
• Mansur, W.J., Brebbia, C.A., (1982a). Numerical implementation of the boundary element method for two dimensional transient scalar wave propagation problems. Applied Mathematical Modelling 6:299-306.
• Mansur, W.J., Brebbia, C.A., (1982b). Formulation of the boundary element method for transient problems governed by the scalar wave equation. Applied Mathematical Modelling 6:307-311.
• Mazzotti, M., Bartoli, I., Marzani, A., Viola, E., (2013). A coupled SAFE‐2.5D BEM approach for the dispersion analysis of damped leaky guided waves in embedded waveguides of arbitrary cross‐section. Ultrasonics 53:1227-1241.
• Mesquita, A.D., Coda, H.B., Venturini, W.S., (2001). Alternative time marching process for BEM and FEM viscoelastic analysis. International Journal for Numerical Methods in Engineering 51:1157-1173.
• Mesquita, E., Antes, H., Thomazo, L.H., Adolph, M., (2012). Transient wave propagation phenomena at visco-elastic half-spaces under distributed surface loadings. Latin American Journal of Solids and Structures 9:453-474.
• Rizos D.C., Karabalis D.L., (1998). A time domain BEM for 3-D elastodynamic analysis using the B-Spline fundamental solutions. Computational Mechanics 22:108-115.
• Rizos D.C. (2000). Rigid surface boundary element for soil-structure interaction analysis in the direct time domain. Computational Mechanics 26:582-591.
• Rizos, D.C., Wang, Z., (2002). Coupled BEM-FEM solutions for direct time domain soil-structure interaction analysis. Engineering Analysis with Boundary Elements 26:877-888.
• Schanz, M., Antes, H., (1997). Application of ‘Operational Quadrature Methods’ in Time Domain Boundary Element Methods. Meccanica 32:179-186.
• Soares, D., von Estorff, O., Mansur, W.J., (2004). Iterative coupling of BEM and FEM for nonlinear dynamic analysis. Computational Mechanics 34:67-73.
• Soares, D. (2008). An optimised FEM-BEM time-domain iterative coupling algorithm for dynamic analysis. Computers & Structures 86:1839-1844.
• Soares, D., Araújo, F.C., (2019). An explicit direct FEM-BEM coupling procedure for nonlinear dynamics. Engineering Analysis with Boundary Elements 103:94-100.
• Tovo, O.A. (2018). A contribution to obtain the transient response of soil-structure systems by an iterative coupling. Master Thesis (in Portuguese), University of Campinas, Brazil.
• Tovo, O., Lima, L.F.V., Mesquita, E., (2019). Study of dynamic response of multidegree of freedom systems through an iterative procedure. In Proceedings of the XL Iberian Latin-American Congress on Computational Methods in Engineering 01-15.
• Von Estorff, O., Prabucki, M.J., (1990). Dynamic response in the time domain by coupled boundary and finite elements. Computational Mechanics 6:35-46.
• Von Estorff, O., Antes, H., (1991). On FEM-BEM coupling for fluid-structure interaction analysis in the time domain. International Journal for Numerical Methods in Engineering 31:1151-1168.
• Wang, Y., Rajapakse, R.K.N.D., (1993). Transient Fundamental Solutions for a Transversely Isotropic Elastic Half Space. Proceedings of the Royal Society of London, 442:505-531.
• Wolf, J.P. (1988). Soil-Structure Interaction Analysis in Time Domain. Englewood Cliffs, NJ: Prentice-Hall 21:363-363.

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  09 Nov 2020
• Date of issue
  2020

History

• Received
  30 Mar 2020
• Reviewed
  31 July 2020
• Accepted
  28 Aug 2020
• Published
  02 Sept 2020