A new global nonreflecting boundary condition with diagonal coefficient matrices for analysis of unbounded media

Accepted Manuscript

A new global nonreflecting boundary condition with diagonal coefficient matrices for analysis of unbounded media M. Mirzajani , N. Khaji , M.I. Khodakarami PII: DOI: Reference:

S0307-904X(15)00617-4 10.1016/j.apm.2015.09.083 APM 10782

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

5 June 2014 27 June 2015 22 September 2015

Please cite this article as: M. Mirzajani , N. Khaji , M.I. Khodakarami , A new global nonreflecting boundary condition with diagonal coefficient matrices for analysis of unbounded media, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.09.083

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights A new semi-analytical global nonreflecting boundary condition (BC).



The proposed BC is employed at medium-structure interface.



The coefficient matrices of BC and dynamic-stiffness matrix are diagonal.



Only the boundary of the problem’s domain is discretized.



The proposed method leads to diagonal differential equations in the frequency domain.

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT

A new global nonreflecting boundary condition with diagonal coefficient matrices for analysis of unbounded media

CR IP T

M. Mirzajani1, N. Khaji*,1, M.I. Khodakarami2

1- Faculty of Civil and Environmental Engineering, Tarbiat Modares University, P.O. Box 14115–397, Tehran, Iran

AN US

2- Faculty of Civil Engineering, Semnan University, P. O. Box 35131-19111, Semnan, Iran

Abstract

In this paper, a new semi-analytical method is developed with introducing a new global nonreflecting boundary condition at medium-structure interface, in which the coefficient

M

matrices, as well as dynamic-stiffness matrix are diagonal. In this method, only the boundary of

ED

the problem’s domain is discretized with higher-order sub-parametric elements, where special shape functions and higher-order Chebyshev mapping functions are employed. Implementing the

PT

weighted residual method and using Clenshaw-Curtis quadrature leads to diagonal Bessel’s differential equations in the frequency domain. This method is then developed to calculate the

CE

dynamic-stiffness matrix throughout the unbounded medium. This method is a semi-analytical method which is based on substructure approach. Solving two first-order ordinary differential

AC

equations (i.e., interaction force-displacement relationship and governing differential equation in dynamic stiffness) allows the boundary condition of the medium-structure interface and radiation condition at infinity to be satisfied, respectively. These two differential equations are then diagonalized by implementing the proposed semi-analytical method. The interaction force* Corresponding author. Tel.: +98-21-82883319, Fax: +98-21-82884914. E-mail address: [email protected] (N. Khaji).

2

ACCEPTED MANUSCRIPT

displacement relationship may be regarded as a nonreflecting boundary condition for the substructure of bounded domain. Afterwards, this method is extended to calculate the asymptotic expansion of dynamic-stiffness matrix for high frequency and the unit-impulse response coefficient of the unbounded media. Finally, six benchmark problems are solved to illustrate

other numerical methods available in the literature.

CR IP T

excellent agreements between the results of the present method and analytical solutions and/or

Key words: semi-analytical method; global nonreflecting boundary condition; Chebyshev

AN US

polynomials; sub-parametric elements; soil-structure interaction; unbounded media

1. Introduction

M

Modeling of unbounded media may be regarded as one of the difficult challenges in the seismic analysis of structures. In the dynamic analysis of unbounded media, the radiation

ED

condition at infinity should be satisfied. Various methods such as finite element method (FEM), boundary element method (BEM), and scaled boundary finite element method (SBFEM) are the

PT

most popular methods for modeling the unbounded media. Each of these methods has their own

CE

advantages and disadvantages. In the FEM, the bounded medium is discretized with finite elements, while the remaining part of unbounded medium is elaborated by particular boundary

AC

conditions (BCs). These BCs are introduced to avoid the artificial reflection of incident waves [1–6]. The dynamic analysis of soil-structure interaction (SSI) problems, based on substructure method, requires investigating the dynamic response of the unbounded medium. Hence, the dynamic-stiffness matrix of the unbounded medium is required to be determined. The unbounded domain far from the soil-structure interface is usually assumed to be isotropic and homogeneous.

3

ACCEPTED MANUSCRIPT

In the substructure method, two substructures are modeled independently. The dynamic property of unbounded domain could be described by force-displacement relationship formulated at the medium-structure interface. By implementing of this relationship in the bounded domain formulation, the total dynamic system of equations may be obtained. In the substructure method,

CR IP T

this relationship could be regarded as a BC for bounded substructure. In the frequency domain, this BC is demonstrated as dynamic-stiffness matrix on the medium-structure interface. Global open BCs have been extensively constructed using the BEM [7], thin-layer method [8], exact non-reflecting BCs [9], and the SBFEM [10]. In the frequency domain SBFEM, the dynamic-

AN US

stiffness matrix is obtained from asymptotic solution and numerical integration of the SBFEM equation. In the time domain SBFEM, the unit-impulse response matrix is calculated from the convolution integral, in which the computational cost for large scale problems with respect to

M

time and space are dramatically expensive, while the accuracy of these BCs are very high [11]. Consequently, proposing suitable and accurate BCs to decrease the computational costs of

ED

practical large scale (e.g., SSI and dam-reservoir interaction) problems may play a crucial role in this regard.

PT

Moreover, various local BCs have been proposed, which are simple and approximate

CE

[12–13], and should be generally applied to a boundary sufficiently far from the region of interest, in order to reach acceptable results. In order to increase the accuracy of these simple

AC

BCs, high-order local boundaries have been proposed [11, 14–16], among which a nonreflecting BC has been proposed for an infinite reservoir [17], dam-reservoir interaction [18, 19] and SSI problems [20]. In the frequency domain analysis, the dynamic-stiffness matrix of the unbounded medium is required to be determined, which is applicable for SSI and dam-reservoir system problems. The soil stiffness matrix may be derived by implementing the FEM formulation by

4

ACCEPTED MANUSCRIPT

idealizing the foundation as an elastic half space [21]. The BEM is an alternative method for the modeling of unbounded media. In this method, no domain discretization is required. The radiation condition at infinity and the governing differential equations in the unbounded domain are exactly satisfied by using fundamental solutions. The main disadvantage of the BEM is that

CR IP T

this method needs a fundamental solution which depends on underlying problems [22–23]. The SBFEM combines the advantages of the FEM and the BEM. This method is a semi-analytical algorithm, for which no fundamental solution is required, and the radiation condition at infinity is satisfied exactly. This method uses the substructure method for modeling the SSI system. A

AN US

disadvantage of this method is that the coefficient matrices are coupled [24–30].

The Spectral Element Method (SEM) is another method which combines the flexibility of the FEM with the accuracy of spectral techniques. This method implements the higher-degree

diagonal mass matrix [31, 32].

M

Lagrange polynomials on Gauss-Lobbatto-Leggendre interpolation points which yield to a

ED

The boundary value problem of an embedded rigid foundation in an orthotropic elastic soil was examined by implementing the indirect BEM [33]. A multi-cell cloning algorithm based

PT

on the FEM formulation was developed to compute the dynamic-stiffness matrix of the

CE

unbounded soil [34]. Wolf and Song [35] developed the consistent infinitesimal finite-element cell method to compute the dynamic-stiffness matrix at the structure-medium interface of an

AC

unbounded medium, in the frequency domain. The infinitesimal finite-element cell method was also developed to estimate the unit-impulse response matrix of the unbounded medium, in the time domain [36]. De Barros and Luco [37] investigated an indirect BEM to obtain the dynamic response of two-dimensional (2D) rigid strip foundation of semi-circular cross-section embedded in a uniform or layered viscoelastic half-space.

5

ACCEPTED MANUSCRIPT

The goal of this paper is to introduce a semi-analytical method for diagonalization of dynamic-stiffness matrix, in the frequency domain, and proposing a new global nonreflecting BC for unbounded domains. The present method has been recently developed for solving potential [38], elastostatic [39, 40] and elastodynamic problems in the time [41] and the frequency [42]

CR IP T

domains. In this method, only the boundaries of the problem’s domain are discretized with

higher-order sub-parametric elements. Using the weighted residual method and implementing Clenshaw-Curtis quadrature leads to a diagonal system of Bessel’s differential equation, in the frequency domain. In other words, the governing differential equation for each DOF is

AN US

independent from other DOFs. In this method, the far field radiation BC at infinity is satisfied, exactly. The basic concepts of the new semi-analytical method are summarized in section 2. The diagonalization procedure of dynamic-stiffness matrix and the high frequency asymptotic

M

expansion are then described in section 3. Afterwards, the unit impulse response coefficient is

ED

presented in section 4. Six benchmark problems are examined and discussed in section 5.

2. Summary of the semi-analytical method

PT

To explain the principals of the semi-analytical method, a general configuration of an

CE

unbounded medium is addressed in Figure 1. As shown in this figure, the side faces that are not passing from the Local Coordinate Origin (LCO) are not discretized. In the semi-analytical

AC

method [38–42], the LCO is selected at a point, from which the whole problem’s domain is visible. For bounded domains, the LCO is chosen inside the domain or on the boundary; while for unbounded domain, the LCO is located outside of the domain. After choosing the LCO, the origin of global 2D Cartesian coordinate is selected at the LCO. The structure-medium interface which is visible from the LCO should be discretized using new sub-parametric elements. In the

6

ACCEPTED MANUSCRIPT

present method, Chebyshev polynomials are used to transform global Cartesian coordinates into local dimensionless coordinates of (  , ). In the local coordinates, the elements’ geometry is interpolated using n  1 mapping functions of φi (η) as follows nη 1

i 1

i 1

x(η)   φi (η) xi , y (η) 

 φ (η) y i

i

(1)

.

CR IP T

nη 1

Besides, any given point in the domain may be related to the corresponding point on the elements of boundaries by using the radial coordinate  as given by

ˆx( , )   x(η) , ˆy( , )   y(η) .

AN US

(2)

It should be noted that the radial coordinate  varies between 1 and infinity for an unbounded medium (see Figure 1).

M

In the present method, the proposed mapping functions φi (η) are higher-order Chebyshev polynomials which are equal to either zero or one at any given control point

ED

(Chebyshev-Lobatto-Legendre or CLL points). Using Jacobian matrix, global coordinates may be transformed into local coordinates as

PT

ˆx (ξ, η) ˆy ,ξ (ξ, η) ˆJ (ξ, η)   ,ξ  ˆx (ξ, η) ˆy (ξ, η) . ,η  ,η 

CE

(3)

Using Eq. (3), the relation between a differential element of area in the global coordinates

AC

( d ˆx d ˆy ) and a differential element of area in the local coordinates ( d ξ d η ) may be written as

d Ω  d ˆx d ˆy  ˆJ(ξ, η) d ξ d η  ξ J(η) d ξ d η

(4)

2.1. In-plane motion

7

ACCEPTED MANUSCRIPT

In the case of in-plane motions such as P and/or SV waves, the displacement field may be depicted as vectors. So, the spatial derivatives in the two coordinate systems are related by the following equation  1 2   b I (η) ξ ξ η

(5)

CR IP T

L  b1I (η)

where 0   ˆy    ˆx , L   ˆy   ˆx   0 T

(6)

AN US

 y ,η (η) 0  1  b ( η)   0  x ,η (η) , J  x ,η (η) y ,η (η)   

0    y ( η) 1 b (η)   0 x(η)  . J  x(η)  y (η)

(7)

(8)

ED

2 I

M

1 I

In the present method, new shape functions N(η) are used for the approximation of

PT

displacement function and its derivatives, which have two specific characteristics: (a) the shape functions have Kronecker Delta property, and (b) their first derivatives are equal to zero at any

CE

given node. In the frequency domain analysis [42], the strain vector may be obtained from the

AC

displacement field at any point (ξ, η) and exciting frequency ω as





1 ˆε(ξ, η,ω)  εˆ x (ξ, η,ω) εˆ y (ξ, η,ω) ˆγ xy (ξ, η,ω) T  b1I (η)Nu ˆ ,ξ  b 2I (η)N ,η u ˆ ξ

(9)

in which





ˆ (ξ, η,ω)  Nu ˆ  N uˆ x (ξ,ω) uˆ y (ξ,ω) T , u

(10)

8

ACCEPTED MANUSCRIPT

B1I (η)  b1I N ,

(11)

B 2I (η)  b 2I N .

(12)

ˆ (ξ, η,ω) is the displacement field at any point (ξ, η) , which may be obtained by where u

Hook’s law, the stress-strain relation may be expressed as

ˆ (ξ,η,ω)  Dˆε(ξ,η,ω) σ

CR IP T

interpolation of the displacement function using special shape functions. Using Eq. (9) and

(13)

AN US

where D is the elasticity matrix of plane strain condition.

The governing equation of elastodynamic problems in the frequency domain is formulated as [42]

σˆ I ij, j  ˆf I i  ρω2uˆ I i  0

(14)

M

where σ I ij shows the stress tensor components, ˆf I i refers to the external source of exciting

ED

forces generated per unit volume for in-plane motion, and ρ is the mass density. After some algebraic elaborations [42], weighted residual method is used in the circumferential direction η

PT

which results in the following system of Bessel’s differential equations

ˆ b (ξ, ω)  0 , ˆ ,ξξ  D1Iu ˆ ,ξ  ω 2 ξ M I u ˆ ξF ξ D 0I u I

CE

(15)

AC

where the coefficient matrices and vectors are given by D 0I 



1

T

B1I D B1I J d η ,

(16)

D1I   B1I,η D B 2I J d η ,

(17)

1

1

T

1

9

ACCEPTED MANUSCRIPT

1

MI   NT ρ N J d η ,

(18)

ˆ b  1N T ˆf b J d η . F I I 

(19)

1

1

CR IP T

2.2. Anti-plane motion

In the case of anti-plane motion such as SH wave, the displacement field is scalar. The

     ˆx   1  1  b 2A (η) ,     b A (η)  ξ ξ  η     ˆy  where

 y ,η (η)  1 1   , B A (η)  b A N ,  x ( η )  ,η 

1 J

b 2A (η) 

1  y (η) 2 2   , B A (η)  b A N . J  x(η) 

PT

ED

M

b1A (η) 

AN US

spatial derivatives for two coordinate systems are considered as the follow relation

(20)

(21)

(22)

In this case, the strain vector may be obtained by the following equation





CE

1 ˆε(ξ, η,ω)  εˆ x (ξ, η,ω) εˆ y (ξ, η,ω) T  b1A (η)Nu ˆ z ,ξ  b 2A (η)N ,η u ˆ z. ξ

(23)

AC

Moreover, the stress-strain relation may be written as

ˆ (ξ,η,ω)  Gˆε(ξ,η,ω) , σ

(24)

where G denotes the shear modulus of the problem’s medium. For the case of anti-plane motion, the governing equation in the frequency domain is

10

ACCEPTED MANUSCRIPT

σˆ zij, j  ˆf zi  ρω2uˆ zi  0 .

(25)

where σˆ zij denotes the stress tensor components, ˆf zi refers to the external source of exciting forces generated per unit volume for anti-plane motion. Similar to in-plane motion, the

ˆ b (ξ, ω)  0 , ˆ z,ξξ  D1A u ˆ z,ξ  ω 2 ξ M A u ˆ z ξF ξ D 0A u A in which 1

D 0A  G  B1A B1A J d η , T

AN US

D1A  G  B1A ,η B 2A J d η , T

1

1

MA  ρ  NT N J d η ,

M

1

ˆ b  1 N T ˆf b J d η . F A A 

(28)

(29)

(30)

ED

1

(26)

(27)

1

1

CR IP T

governing equation may be derived as [42]

PT

2.3. Diagonal coefficient matrices

As already discussed, the present method employs special higher-order sub-parametric

CE

elements, based on Clenshaw-Curtis quadrature and implementing weighted residual method

AC

[38–42]. Consequently, the present method consists of diagonal coefficient matrices, which results in Bessel’s differential equations in the frequency domain [41, 42]. For in-plane motion, the coefficient matrices are transformed into the following forms DI0ij  2 ijWi B1I ( i )D B1I ( i ) J( i ) , T

(31)

11

ACCEPTED MANUSCRIPT

DI1ij  2 ijWi B1I, ( i ) D B 2I ( i ) J( i ) ,

(32)

M Iij  2 ijWi N T ( i )  N( i ) J( i ) .

(33)

T

For the anti-plane motion, the corresponding coefficient matrices may be computed as DA0 ij  2 ijWi B1A (i ) G B1A (i ) J(i ) , T

CR IP T

(34)

DA1 ij  2 ijWi B1A, ( i )G B 2A ( i ) J( i ) , T

(35)

M Aij  2 ijWi N T ( i )  N( i ) J( i ) ,

AN US

(36)

where δij denotes Kronecker delta, and Wi implies the weighting values of Clenshaw–Curtis quadrature. Lastly, Eqs. (15) and (26) become diagonal as a set of single Bessel’s differential equations for ith DOF, as given by the following relations

(37)

 DI0ii uˆ yi ,  DI1ii uˆ yi ,   2 M Iii uˆ yi   Fˆ Iyib  0 ,

(38)

ED

M

 DI0ii uˆ xi ,  DI1ii uˆ xi ,   2 M Iii uˆ xi   Fˆ Ixib  0 ,

PT

for in-plane motion. For anti-plane motion, one may write (39)

CE

 DA0ii uˆ zi ,  DA1 ii uˆ zi ,   2 M Aii uˆ zi   Fˆ Abzi  0 .

AC

Eqs. (37)–(39) are decoupled Bessel’s differential equations which depend only on the ith DOF.

3. Dynamic-stiffness matrix The present semi-analytical method is based on substructure method. In the substructure method, the bounded structure and unbounded soil are modeled independently. The bounded domain may be modeled with the present method [41, 42] or other numerical methods (e.g., the FEM). The unbounded soil is currently being modeled with various numerical methods such as the FEM, the 12

ACCEPTED MANUSCRIPT

BEM, and the SBFEM. In this paper, the new semi-analytical method is used to diagonalize the dynamic-stiffness matrix of the unbounded medium. To this end, the medium-structure interface is discretized with sub-parametric elements. The dynamic behavior of unbounded medium is described by the interaction force-displacement relationship on the medium-structure interface

CR IP T

(see Figure 2). The interaction force-displacement relationship could be regarded as a global nonreflecting BC. In the present method, it is assumed that the nodes are connected to each other with springs of certain stiffness. As schematically shown in Figure 2, all physical and numeral features of the problem may be gathered in the LCO. This point will be discussed in the coming

AN US

sections. In the present method, the LCO is identical for all nodes. In other words, the LCO has the same displacement components for all nodes. Therefore, the physical concept of this phenomenon may be considered as a few parallel springs with distinct stiffnesses, which adjoin

M

to each other at the LCO (see Figure 2).

ED

3.1. Dynamic-stiffness matrix derivation

ˆ ( ) is related to interaction forces vector In the frequency domain, the displacement vector u

PT

R( ) by the dynamic-stiffness matrix S  ( ) as follows

ˆ ( ) . R( )  S  ( )u

CE

(40a)

AC

In fact, the force-displacement relationship describes the contribution of unbounded medium in the governing equation of motion of dynamic SSI and dam-reservoir system problems. In other words, computation of dynamic-stiffness matrix is a crucial step in these problems. As shown in Figure 2, only the soil-structure interface is required to be discretized, in the present method. Therefore, the free surface of soil does not require any discretization. Eq. (40a) could be regarded as a global nonreflecting BC at the present method. As linear behavior of the domain is 13

ACCEPTED MANUSCRIPT

considered in the present paper, the Fourier and Inverse Fourier transformation may be applied. Applying the Inverse Fourier transformation to Eq. (40a) leads to t

R(t )   S  (t -  )u( )d ,

(40b)

0

CR IP T

which is a convolution integral that is diagonalized using the present method. In other word, each DOF is independent from other DOFs.

The governing differential equations of in-plane (I) and anti-plane (A) motions have been obtained in the previous section. Now for simplicity of formulations, the subscripts of I and A

AN US

would not be mentioned in the following sections. To apply the BCs at infinity and at the

medium-structure interface, the governing displacement differential equation is not solved, directly. It is replaced by two first-ordered differential equations: (1) interaction force-

method in dynamic stiffness [14].

M

displacement relationship, and (2) the governing differential equation of the new semi-analytical

ED

In order to determine the interaction force-displacement relationship as in the FEM, the internal nodal forces Q  Q( ) , which are equal to the resultants of the stresses  , on a line with

PT

a constant  plus the redistribution effect, which is similar to the moment distribution at classic

CE

structural analysis [45], for internal nodal forces appeared in the present method are taken into account. Implementing the principle of virtual work results in the following expression w T Q   ξ wn ξ σdS ξ , T

(41)

AC

S

in which d S ξ and w are determined as (see Figure 3)

d S ξ  ˆx,2 (η)  ˆy ,2 (η)d   x,2 (η)  y ,2 (η)d ,

(42)

w  w ( , ,)  N(η)w( ,) .

(43)

14

ACCEPTED MANUSCRIPT

Substituting Eqs. (42), (43), (21) into Eq. (41) for an arbitrary w( , ) leads to 1

ˆ J d . Q    N T b1 σ T

(44)

1

The first derivative of shape functions used in the present method is equal to zero at given nodes

CR IP T

( N ,  0 ). As a result, using Eq. (11) and by substituting Eqs. (13) and (24) into Eq. (44), and considering the redistribution effects of the present method (as shown in Figure 2), the following relation may be obtained for 1D anti-plane problems 1

1

ˆ A ,ξ J d  ς  B1 GB1u ˆ A J d , Q    B1 GB1u T

T

1

AN US

1

which may be rewritten as

ˆ A ,ξ  ςD0A u ˆ A. Q A  ξ D0A u

(45)

(46)

M

In the present method, a coefficient matrix ς is introduced whose components denote the

ς Aii 

DA0 ii   A i DA1 ii m

D

(47)

PT

j 1

,

0 Ajj

ED

redistribution coefficients, and may be obtained for 1D problems as given by

CE

in which m is the number of DOFs of the problem and

Aii 

DA1 ii

m

D

(48)

AC

j 1

.

1 Ajj

Using Eqs. (16) and (27), Eq. (46) may be rewritten as a set of single differential

equations. On the other hand, each DOF may be independent from other DOFs as given by

15

ACCEPTED MANUSCRIPT

 D110 z Q1z          Q   0  nz  A 

 D110 z  0  uˆ 1z             ς Aii   0     0  Dnnz  A uˆ nz  A , 

 0  uˆ 1z           , 0     Dnnz  A uˆ nz  A

(49)

or

QA i  ξ DA0 ii uˆ Ai,ξ  ς Aii DA0 ii uˆ A i .

CR IP T

(50)

The relation between nodal forces’ vector R(ξ ) and internal nodal forces’ vector Q(ξ ) may be introduced as follows (see Figure 4)

R  Q ,

AN US

(51)

in which positive and negative signs represent the bounded and unbounded media, respectively.

ˆ ( , ) and nodal forces vector R( , ) may be The relation between displacement vector u

ˆ ( ,) R( ,)  S( ,)u

M

rewritten in the frequency domain as

(52)

ED

where S( , ) denotes the dynamic-stiffness matrix on a line with constant ξ . Using Eq. (50),

PT

Eq. (52) may be rewritten as

ˆ ( ) R F ( ) , R( )  S( , )u

(53)

CE

in which R F ( ) implies the amplitude of nodal force due to the body load and surface traction.

AC

Substituting Eqs. (46), (53), and (45) into Eq. (51) results in the following equation

 Suˆ A  R F  ξ D0A uˆ A ,ξ  ςD0A uˆ A

(54)

Differentiating Eq. (54) with respect to ξ leads to

ˆ A  Su ˆ A,  R F,  ςD0A u ˆ A,   D0A u ˆ A,  D0A u ˆ A,  0 .  S , u

16

(55)

ACCEPTED MANUSCRIPT

Adding Eq. (26), and Eq. (46) multiplied by ξ , and considering the mode effects [44], results in the following relation





ˆ b  0, ˆ A ,ξξ  ξ χ A D1A  D 0A u ˆ A ,ξ  n 2 D 0A u ˆ A  ω2 ξ 2 M Au ˆ A ξ2 F ξ 2 D 0A u A

(56)





CR IP T

or 1 1 1 b ˆ  0, ˆ A    I  χ A D0A D1A u ˆ A   n 2u ˆ A   2 2 D0A M A u ˆ A  ξ 2 D0A F  2u A ,

,

(57)

where n is the number of modes. Multiplying Eq. (55) by ξ and adding Eq. (56) gives





AN US

ˆ b  0 .(58) ˆ A  ξ  S   A D1A  ς D0A u ˆ A ,ξ  R F,  n 2 D0A u ˆ A  ω2 ξ 2 M Au ˆ A ξ2 F  ξ S , u A ˆ A ,ξ and substituting into Eq. (58) yields Solving Eq. (54) for u 1







 ξ S , uˆ A  D 0A  S  χ A D1A  ςD 0A  S  ς D 0A uˆ A 1





.

(59)

M

ˆ b 0  R F,  D 0A  S  χ A D1A  ςD 0A R F  n 2 D 0A uˆ A  ω 2 ξ 2 M A uˆ A  ξ 2 F A

1

ED

ˆ A , the remaining coefficient matrix of Eq. (59) should vanish as For arbitrary u







(60)

PT

  S ,  D0A  S  χ A D1A  ςD0A  S  ς D0A  n 2 D0A  ω 2 ξ 2 M A  0 .

The matrix S may be determined using Eq. (60). The remaining part of Eq. (60) yields to

ˆ b . The boundary condition is R F (  0)  0 for R F,  D0A  S  χ A D1A  ςD 0A R F   ξ 2 F A

CE

1

AC

bounded domain, and R F (  )  0 for unbounded domain. As  S ,   S ,  applies, and

D0A and M A are independent of ξ and ω , it may be concluded from Eq. (60) that S is a function

of ωξ only. In other words, S is a function of dimensionless frequency introduced for the line with constant ξ as given below

17

ACCEPTED MANUSCRIPT

a

ωr0 , cs

(61)

where c s  G ρ is shear wave velocity, and r0 denotes the characteristic length of the boundary ( ξ  1 ). The derivative with respect to a may be performed either with varying ω

CR IP T

(with fixed ξ ) or with varying ξ (with fixed ω ). For 2D cases, S  S(a) may be written as

a S ,a (a)   S ,   S , .

(62)

As a result in Eq. (60), the derivation with respect to ω may be replaced by the derivation with

1



AN US

respect to  as given below





  S ,  D0A  S  χ A D1A  ςD0A  S  ς D0A  n 2 D0A  ω 2 ξ 2 M A  0 .

(63)

The dynamic-stiffness matrix for the bounded medium can be demonstrated by the following

1

M

relation 1

 Sb,  D0A Sb  (χ A D0A D1A  2ς )Sb  ς(χ A D1A  ςD0A )  n2D0A   2 2MA  0 ,

(64)

ED

2

where superscript b stands for the bounded medium. It should be noted that Eq. (64) is written on

1

PT

the boundary ( ξ  1 ). Similarly, for the unbounded medium one may write 1

  S,  D0A S  (χ A D0A D1A  2ς )S  ς(χ A D1A  ςD0A )  n2D0A   2 2MA  0 . (65)

CE

2

AC

Eqs. (64) and (65) are the governing differential equations in dynamic stiffness in the frequency domain for 1D bounded and unbounded media. These equations are the system of non-linear first-order ordinary differential equations with the frequency ω , which are diagonalized with the new semi-analytical method.

3.2. Diagonal system of governing differential equations in dynamic stiffness 18

ACCEPTED MANUSCRIPT

As mentioned in section 2, the new semi-analytical method uses Clenshaw-Curtis quadrature and special shape functions for diagonalization of the coefficient matrices. As shown in Eqs. (31)–(36), the coefficient matrices of Eqs. (64) and (65) are diagonal. This means that the coupled system of differential equations in dynamic stiffness has been transformed into

CR IP T

decoupled differential equations. On the other hand, the dynamic stiffness at a given point may be evaluated by solving the governing equation corresponding to that point only. In this regard, the system of differential Eqs. (64) and (65) may be converted to a single differential equation

S

b ii ,

AN US

related to the ith DOF, as given by the following forms

DA1 ii 1 b2  0 Sii  (  A ii 0  2 A ii ) Siib   A ii (  A ii DA1 ii   A ii DA0 ii )  n 2 DA0 ii   2 2 M A ii  0 (66) DA ii DA ii DA1 ii 1 2 S  (   2 Aii ) Sii   Aii (  Aii DA1 ii   Aii DA0 ii )  n 2 DA0 ii   2 2 M Aii  0 . (67) ii Aii 0 0 DAii DAii

M

  Sii, 

ED

Eq. (66) is usually not computed. Instead, by ignoring damping effects, the low frequency expansion of S b is evaluated using the static-stiffness matrix K b , and the mass matrix M b as

PT

given below

S b  K b   2 M b  O( 4 ) .

CE

(68)

AC

In the following sections, the detailed calculations of K b and M b are discussed.

3.3. Analytical solution for diagonal system of differential equations 3.3.1. Statics

As Eqs. (56), (64) and (65) are diagonal, decoupled form of equations for the ith DOF may be taken into account. Therefore, in static status (   0 ) of the present method, for

19

ACCEPTED MANUSCRIPT

vanishing body forces, the governing differential equations in displacement and stiffness are given by

 2 DA0 uˆ A    A DA1  DA0 uˆ A  n 2 DA0 uˆ A  0 , ii

i ,

ii

ii

ii

i ,

ii

(69)

i

(70)

DA1 ii 1 2 K  (   2 A ii ) K Aii   A ii (  A ii DA1 ii   A ii DA0 ii )  n2 DA0 ii  0 , A ii A ii 0 0 DA ii DA ii

(71)

CR IP T

DA1 ii 1 b2 K  (   2 Aii ) K Ab ii   Aii (  Aii DA1 ii   Aii DA0 ii )  n 2 DA0 ii  0 , A ii A ii 0 0 DAii DAii

AN US

where K iib and K ii are diagonal components of static-stiffness matrices. The general solution of Eq. (69) may be given as the following form for 1D problems

uˆ i ( )    i .

(72)

M

Substituting Eq. (72) into Eq. (69) yields a quadratic algebraic equation, whose two solutions for

 i1,2

ED

 i may be derived as  1   A DA1 1   A DA 2    4n    DA0 2 DA0   ii

ii

PT

ii

ii

ii

ii

   

2

  .  

(73)

CE

As a result, the general solution of Eq. (69) may be written as

uˆ i ( )  c1 i1  c2 i 2 ,

(74)

AC

where c1 and c 2 imply the integration constants. As for the bounded case, uˆ i ( ) involves finite values at the LCO (   0 ), the constant c 2 should be zero. In other words,

uˆ i ( )  c  i1 .

(75)

20

ACCEPTED MANUSCRIPT

Considering the positive root of Eq. (70), the static-stiffness coefficients at the bounded case may be easily determined by the following form 1 K iib     A DA1  2 A DA0  4n 2 DA0 2 ii

ii

ii

ii

2   A2 DA1  

2 ii

ii

(76)

ii

CR IP T

For the unbounded case when ξ   , the displacement leads to zero and therefore c1  0 . As a result,

uˆ i ( )  c i 2 .

(77)

AN US

Again, the positive root of Eq. (71) is equal to the diagonal coefficient of static-stiffness matrix of unbounded domain as 1 1 0 2 0   A DA  2 A DA  4n DA 2

K ii 

ii

ii

ii

2 ii

2   A2 DA1  .  ii

ii

(78)

M

ii

3.3.2. Dynamics

ED

Since Eq. (57) is diagonal, decoupled form of equations for the ith DOF may be

PT

considered. Consequently, in dynamic status for vanishing body forces, the governing differential equations in displacement in the present method is written as

 2 uˆ A, I    i  1uˆ A, I  n 2 uˆ A, I  aA,2 I uˆ A, I  0 , i ,

CE

i ,

AC

where  A, I  i

 A, I DA,1 I ii

DA,0 I

ii

i

(79)

i

.

(80)

ii

3.3.2.1.1D problems The solution of Eq. (79) for the bounded case is as the following form 21

ACCEPTED MANUSCRIPT

uˆ i ( )  c1



i 2

J  i (a)  c 2



i 2

(81)

Yi (a) ,

where J i (a) and Yi (a) are Bessel’s functions of the first and second kind, respectively. 1 2

 i2  4n 2 indicates the order of Bessel’s functions for each DOF. In

CR IP T

Furthermore,  i 

addition, c1 and c 2 are the integration constants. Assuming a finite displacement in the domain, Eq. (81) is rewritten as given by 

i 2

(82)

J  i (a) .

AN US

uˆ i ( )  c

Using Eq. (82), the ith component of nodal force vector R may be derived as Ri   Dii0 uˆ i ,   ii Dii0 uˆ i .

(83)

The diagonal components of dynamic-stiffness matrix may be determined for the bounded case

Ri . uˆ i

(84)

ED

S iib 

M

as

PT

The mass matrix M b of the bounded domain is defined using the low frequency expansion of the dynamic-stiffness matrix (see Eq. (68)). Substituting Eq. (68) into Eq. (64), and

CE

considering Eq. (76), the diagonal coefficient of mass matrix may be found after some algebraic

AC

manipulations as M iib 

M ii . 2(2 ii  1)   i   i 

(85)

For the unbounded medium, the solution of Eq. (79) is given by

uˆ i ( )  c1



i 2

H  i (a)  c 2 (1)



i 2

(86)

H (2i ) (a) ,

22

ACCEPTED MANUSCRIPT

where H (1i) (a) and H (2i ) (a) denote the Hankel functions of the first and second kinds of order  i for each DOF, respectively. The radiation condition is represented by implementing the asymptotic behavior of Hankel functions for ξ   as 

 



2 exp  i a  a

H (2) (a) 

2 exp  i a  a

2

2



 4 ,

(87)

CR IP T

 

H (1) (a) 



 4 ,

(88)

AN US

where i  (1)1/2 .

Eqs. (87) and (88) are corresponding to the wave propagating from infinity and toward infinity, respectively. As only wave propagation toward infinity are permitted, c1  0 and then 

i 2

H (2i ) (a) .

ED

The nodal force equals Ri   Dii0 uˆ i ,   ii Dii0 uˆ i .

(89)

M

uˆ i ( )  c

(90)

PT

Therefore, using Eq. (90), the diagonal coefficients of dynamic-stiffness matrix for the

CE

unbounded medium may be derived as Ri uˆ i

(91)

AC

S ii  

As Eqs. (84) and (91) obviously indicate, the dynamic-stiffness components are functions of dimensionless frequency a . On the boundary of the problem’s domain ( ξ  1 ), Eq. (91) is rewritten as

23

ACCEPTED MANUSCRIPT

S ii (a0 )  

Ri , uˆ i

(92)

where the dimensionless frequency on the boundary equals

ω . cs

(93)

CR IP T

a0 

3.3.2.2. 2D problems

In order to solve 2D problems, the lame potentials  ( ) and  ( ) are introduced [1, 46].

AN US

Therefore, the displacement fields of Eq. (79) may be written as (see Figure 1b) 1 uˆ I ( ,  )  , ( ,  )   ( ,  ) ,



1

1

 ( ,  )   , ( ,  ) .  

(94b)

M

uˆ I ( ,  ) 

(94a)

ED

Substituting Eqs. (94) into Eq. (79) yields to the potential governing differential equations as given by

(95a)

  DI1(i 1) (i 1)  ωξ 2 2     1 )  (i 1)  0 , ( i 1) ,  n  ( i 1)  ( 0  DI ( i 1) (i 1)  cs  

(95b)

CE

  (i 1) ,

( i 1 )

AC

2

i

PT

   DI1ii  ωξ  1 i ,  n 2 i  ( ) 2 i  0 , 0  DIii  cp  

 2 i ,   

where c p 

2 (1 ) 1 2 s

c is dilatational wave velocity. As shown in Figure 1b, the BCs on the

structure-medium interface are assumed as uˆ I (  1)  u 0 uˆ I (  1)   u 0

(96)

24

ACCEPTED MANUSCRIPT

The general solutions of Eq. (95) for the bounded medium may be given by the following form 

i 2

i

J ( i

 (i 1) ( )  d1



 ωξ ωξ )  c2 2 Y ( ) , cp cp

 ( i  1 ) 2

(97a)

i

J

( i 1 )

 ωξ ( )  d 2 cs

 ( i  1 ) 2

Y

( i 1 )

(

ωξ ) , cs

while for unbounded case are 

i 2

i

H (1) ( i

 (i 1) ( )  d1



 ωξ ωξ )  c2 2 H (2 ) ( ) , cp cp

  ( i 1 ) (1)

2

(98a)

i

H

( i 1 )

 ωξ ( )  d 2 cs

 ( i 1 ) 2

AN US

i ( , )  c1

(97b)

CR IP T

i ( , )  c1

H (2) (  ( i 1 )

ωξ ) . cs

(98b)

Similar to the previous section, the selected solutions for bounded and unbounded cases are

2

i

 (i 1) ( )  d1



2

i 2

J

( i 1 )

(

ωξ ) , cs

H (2) (

AC

 (i 1) ( )  d 2

i



 ( i  1 ) 2

H (2) (  ( i 1 )

(99a)

(99b)

ωξ ωξ )  c2  B H (A2 ) ( ) , cp cp

CE



ωξ ), cp

 ( i 1 )

and

i ( , )  c2

J (

M

i

ED



PT

i ( , )  c1

(100a)

ωξ ωξ )  d 2  Z H (2) ), E ( cs cs

(100b)

in which c2 and d 2 are the integration constants, which are obtained using the BCs introduced in Eq. (96). Implementing these BCs, the unknown coefficients for unbounded domain are obtained as the following relations 25

ACCEPTED MANUSCRIPT

c2 

 2     

d2 

u0  2 , ( 2)      H E (a0 )

u0 , a c H (A2 ) ( 0 s ) cp

(101a)

CR IP T

(101b)

where   a0

H (E2 )1 (a0 ) , H (E2 ) (a0 )

(102a)

cs a0 ) cp cs   a0 . c cp H (A2 ) ( s a0 ) cp

AN US

H (A2)1 (

(102b)

After determining these coefficients, the ith and (i+1)th components of nodal force vector

M

R may be derived as

1 1 Ri   DI0ii (i , ) ,  PDI0ii (  (i 1) ) ,  DI0ii (  (i 1) ) ,

ED





1 1 Ri 1   DI0(i 1)(i 1) (  (i 1) , ) ,  TDI0(i 1)(i 1) ( i ) , ,

PT



cs 2 ( A  B  1)( A  B  1) ) (Z  1)  , cp A B3

AC

P  2(

CE

where

T



  (

(103b)

(104)

cp 1 ( ) 2 A 2  A(2 B  1)  B( B  5)  2 , A  B  3 cs



(103a)



cs 2 2 ( A  B  1)( A  B  1)(E  Z  3) ) ( E  1  ( Z  1) 2 )  . cp A B3

26

(105)

(106)

ACCEPTED MANUSCRIPT

The diagonal coefficients of dynamic-stiffness matrix for the unbounded medium may be obtained as      1 1  S ii  (  ) R  ( ) R  i i 1  0 c u0  D p I ( i 1)( i 1)  DI0ii ( ) 2   c s  

CR IP T

(107a)

in which [1] 

   cos2 ( )d .

(107b)

AN US

0

3.4. Asymptotic expansion for high frequency

In the first step, the governing differential equation in dynamic stiffness for 1D unbounded medium (Eq. (65)), which introduces the radiation condition at infinity, should be

M

solved. The radiation condition is calculated at a high but finite frequency ( ωh ) from high

ED

frequency asymptotic expansion. In fact, the high frequency asymptotic expansion of the dynamic-stiffness matrix permits the radiation condition to be satisfied exactly. The governing

PT

differential equation is then solved starting from ωh . For this purpose, the dynamic-stiffness

CE

matrix is extended in a power series of order m as given by [24] m

1 Aj . j j 1 (i )

S  ( )  i C   K   

AC

(108)

The first two terms on the right hand side of Eq. (108) shows the singular part S s ( ) with the constant dashpot matrix C  and the constant spring matrix K  . The third term indicates the regular part S r ( ) , which is expanded implementing m terms of the power series with the unknown coefficient matrices A j . 27

ACCEPTED MANUSCRIPT

The radiation condition may be satisfied by formulating a positive definite C  from an eigenvalue problem. The starting point of the eigenvalue problem is presented as

MΦ  D0 ΦΛ 2

(109)

CR IP T

in which the coefficient matrices D 0 and M have been diagonalized using the new semianalytical method. In addition, the eigenvectors Φ are normalized as given by

Φ T D0 Φ  I .

(110)

In addition

AN US

Φ T MΦ  Λ 2 , 1

D0  ΦΦT .

(111) (112)

By pre-multiplying Eq. (65) by Φ T , and post-multiplying by Φ , the following results

M

may be obtained

  s ,  s   (χe1  2ς ) s   ς(-χe1  ς)  n 2 I   2  2  0

(114) (115)

CE

e1  Φ T D1Φ .

PT

where

s   ΦT S Φ ,

(113)

ED

2

AC

Substituting Eq. (108) into Eq. (114) leads to m

1 aj , j j 1 (i )

s  ( )  i c   k   

(116)

where

c   Φ T C Φ ,

(117)

28

ACCEPTED MANUSCRIPT

k   ΦT K Φ ,

(118)

a j  ΦT A j Φ .

(119)

Considering only two terms (i.e. m  2 ), and substituting Eq. (116) into Eq. (113) results

2

CR IP T

in the following relation

or







AN US

 a 2a 2   a a     i c   k   1  2 2    i c   1  2  i (i )   i (i )    a a    2ς  χe1  i c   k   1  2 2   ς  χe1  ς   n 2 I   2 Λ 2  0 i  (i )  

  k  2c a   2ς  χe k  ς χe  ς   n I  1  a  2c a  2k a   2ς  χe a   0 i

(120)

(i ) 2 c 2  Λ 2  (i ) 2c  k    2ς  χe1 c   c  2 

1



1

2



1

1





2

1

1

M

1

(121)

The coefficient matrix of each term of the power series in i (Eq. (121)) should be set equal to

ED

zero in a descending order. The first term gives (122)

PT

c 2  Λ 2 .

Choosing the positive roots of the diagonal Λ 2 yields

CE

c  Λ .

(123)

AC

The dashpot matrix C  is determined using Eq. (117)

C   Φ T ΛΦ1 .

(124)

As can be inferred from Eq. (124), the radiation condition is exactly satisfied, when a positive definite matrix is constructed. So, the unbounded medium acts as a sink of energy (radiation condition). 29

ACCEPTED MANUSCRIPT

Implementing Eq. (123) in the second term of Eq. (121), and setting it equal to zero results in the following relation

k 



 

1  2ς  χe1  I . 2

(125)

a1 

1 1 2 Λ n I  k 2   2ς  χe1 k   ς χe1  ς  , 2

a2 

1 1 Λ  a1  2k a1   2ς  χe1 a1 , 2



 

(126)

(127)

AN US



CR IP T

Following the same procedure for the third and fourth terms of Eq. (121) results in

Afterwards, the spring matrix K  and coefficient matrix A j may be calculated using Eqs. (118) and (119) as

M

K   Φ T k  Φ 1 ,

(129)

ED

A j  Φ T a j Φ 1 .

(128)

As shown in the previous sections, the system of governing differential equations in

PT

dynamic stiffness becomes diagonal, using the new semi-analytical method. This means that, each DOF is independent of other DOFs and can be solved separately. Furthermore, all

CE

coefficient matrices of the present paper are diagonal. For solving the problems using the

AC

proposed method, after calculating the matrices C , K  , and A j , the asymptotic behavior of dynamic-stiffness matrix is firstly determined using Eq. (108). At the second step, Eq. (67) is  solved with known S ( h ) from the first step, for each DOF. However, in the present method

the dynamic-stiffness matrix may be calculated analytically after redistribution. This may be achieved by implementing Eq. (67) for unbounded domain with asymptotic expansion at infinity

30

ACCEPTED MANUSCRIPT

as BC. For bounded domain, the mentioned procedure is performed using Eq. (66) (or Eq. (68) for low frequency expansion, instead).

4.1. Inverse Fourier transform

CR IP T

4. Unit-impulse response

Using Eq. (52), the force-displacement relationship on the boundary in the frequency domain may be written by the following equation

ˆ 0. R 0  Su

AN US

(130)

Moreover, Eq. (108) with m  1 is rewritten as S  ()  i C   K   S r ( ) .

(131)

The regular part of dynamic-stiffness matrix is square integrable. Using the inverse Fourier

(132)

ED

S  (t )  C (t )  K  (t )  S r (t ) ,

M

transform of Eq. (131), the unit-impulse response may be obtained from

PT

where  (t ) indicates the Dirac delta function, and S r ( ) and S r (t ) are Fourier transform pairs. Inverse Fourier transform of Eq. (130) with substituting Eq. (131) results in [6] t

CE

R(t )  C  u (t )  K  u(t )   S r (t   )u( ) d . 0

(133)

AC

The first two terms on the right-hind side of Eq. (133) depicts the instantaneous response, and the convolution integral indicates the delaying part of the response. At the first step and using Eq. (65), the corresponding equation in the regular part of

dynamic stiffness is formulated. To this end, by implementing Eqs. (65) and (131), and after some algebraic procedures, the following relations may be obtained

31

ACCEPTED MANUSCRIPT

1



1

1

1



  S r ,  D 0 S r  2iD 0 C  S r  2D 0 K   χD 0 D1  2ς S r



2



 D 0 K 2  χD 0 D1  2ς K   ς χD1  ς D 0   n 2 D 0  0 1

1

(134)

The numerical discretization is simplified using the following abbreviation

1  S r ( ) , i

(135)

CR IP T

Vr ( ) 

which may be rewritten in the time domain as t

Vr (t )   S r ( )d

(136)

0

AN US

or

  (t ) . S r (t )  V r

(137)

2 Substituting Eq. (135) into Eq. (134) and dividing by (i ) yields

1



1

1

1



2





0 1





(138)

D1  2ς K   ς χD1  ς D 0   n 2 D 0  0

ED

 Vr,  (i )  2 D K 2  χD

0 1

M

D 0 Vr  2D 0 C  Vr  (i ) 1 2D 0 K   χD 0 D1  2ς Vr  (i ) 1 Vr

The inverse Fourier transform is applied to Eq. (138), which results in 1

V t

0

 r



1

1

1

(t   )Vr ( )d  2D 0 C  Vr (t )  2D 0 K   χD 0 D1  2ς

PT

D0







 V t

0

 r

( )d



  Vr ( )d  tVr (t )  t D 0 K 2  χD 0 D1  2ς K   ς χD1  ς D 0   n 2 D 0 H (t )  0 t

1

(139)

CE

0

1

where H (t ) denotes the Heaviside step function. Eq. (139) corresponds to the equation of new

AC

semi-analytical method in the time domain for Vr (t ) . On the other hand, the high frequency asymptotic expansion of the dynamic-stiffness matrix, Eq. (108), is related to the early-time asymptotic expansion of the unit-impulse response matrix. Finally, the inverse Fourier transform is applied to Eq. (108) with m  1 , which leads to the following relation

32

ACCEPTED MANUSCRIPT

S (t )  C(t )  K  (t )  A1H (t )

(140)

4.2. Time discretization

CR IP T

By comparing Eqs. (132) and (140) at t  0 , the regular part of the unit-impulse response matrix is evaluated as given by

S r (0)  A1 ,

(141)

AN US

 which indicates that S r (0) is finite. Consequently, using Eq. (139) for t  0 results in

Vr (0)  0 .

(142)

In order to solve the integral equation (139), the equation should be discretized with respect to   time to get an equation for Vr m at each time step m ( m  1 ). In addition, it is assumed that Vr m is

0

j 1

PT



mt

0

m 1

Vr ( )d  t  Vrj  12 t Vrm ,

(143)

m 1

Vr (mt   )Vr ( )d  t  Vrj Vrm  j .

CE



mt

ED

terms in Eq. (139) are as

M

piecewise constant over the time of (m  12 )t  t  (m  12 )t . It may be shown that the integral

(144)

j 1

AC

 and finally, a linear equation for Vr m may be obtained by substituting Eqs. (143) and (144) into

Eq. (139) formulated for t  mt , as given below

33

ACCEPTED MANUSCRIPT



2 t



D 0 C   D 0 K   12 χD 0 D1  ς  m  12 I Vr 1

 D 0



1

1

m 1

V

 rj

j 1

1



1

1

m

Vr  2D 0 K   χD 0 D1  2ς  I m j





 V m 1 j 1

 rj

 m D 0 K 2  χD 0 D1  2ς K   ς  χD1  ς D 0   n 2 D 0 1

1

(145)



CR IP T

Eq. (137) is then employed to calculate the regular part of the unit-impulse response matrix S r m as

m

1  V  Vr t r m

m 1

.

(146)

AN US

S r 

5. Numerical Examples

In order to assess the accuracy of the present method for analyzing the problems involving dynamic-stiffness with diagonal coefficient matrices and examine the efficiency of

M

proposed global nonreflecting BC, six benchmark examples are selected and solved by the present method. The results of the present method are compared with analytical solutions and/or

ED

other numerical methods available in the literature. No physical damping is considered in all

PT

examples. The measurements are in the SI units.

CE

5.1. Semi-infinite layer with constant depth The first example is a semi-infinite layer with constant depth which has a cut-off

AC

frequency. The geometry of this example is shown in Figure 5. In order to discretize the semiinfinite layer, the LCO is located at infinity and only vertical side of semi-infinite layer is discretized using three-node line elements [47]. As the LCO is selected at infinity, the geometry of the elements in local coordinate may be expressed as

34

ACCEPTED MANUSCRIPT

x   xb   , y   y b ,

(147)

in which b refers the coordinates of elements on the boundary. Using Eq. (147), one may write 0  1 J     , 0 y  , 

CR IP T

(148)

and

AN US

1 b 1      0 1 1 b 2      J 0

(149)

Following the procedure described in this paper and after a couple of algebraic manipulations,

ˆ ,  n 2 D 0 u ˆ   2 Mu ˆ  0. D0 u

M

the governing differential equation in the frequency domain may be derived as (150)

ED

Implementing the present method for diagonalization of coefficient matrices leads to a

PT

diagonal dynamic-stiffness matrix, whose diagonal components may be normalized as given by 2

S i a   1  0  . n n

CE

(151)

This dynamic stiffness coefficient is shown in Figure 6 and compared with analytical solution

AC

[11], which shows excellent agreement. Furthermore, in order to examine the global nonreflecting BC (i.e., force-displacement

relationship) proposed in this paper, the unit-impulse response coefficient is calculated using the present method and compared with exact solutions [11] in Figure 7. In addition, Figure 9 illustrates the response of semi-infinite layer to a surface traction prescribed as the Ricker 35

ACCEPTED MANUSCRIPT

wavelet (Figure 8). The response is computed for n  5 , 10, 15 and compared with exact solutions [11], with excellent agreement as may be observed in Figure 9.

5.2. Truncated semi-infinite wedge

CR IP T

The second example examines the anti-plane motion of a truncated semi-infinite wedge with a free and a fixed boundary extending to infinity (see Figure 10). The side 1 2 may be considered as a structure-medium interface. The boundary of the problem’s domain is discretized

AN US

using one two-node element. The node positions are as: x1  x2  r0 , y1   r0

5 , y2  0 ,

where r0 is the characteristic length of the boundary. In this example, the redistribution coefficients are calculated using Eqs. (47)-(48) as

 A  0 , DA1  0 , 22

ςA 

1 . 2

(152a)

(152b)

ED

22

M

22

Then, the governing differential equation in displacement is written as

2 ,

2 ,



PT

 2 uˆ A   uˆ A  n 2 uˆ A  ( 2

cs

) 2 uˆ A  0 .

(153)

2

CE

The static stiffness coefficient is calculated using Eq. (78) as

AC

K 22  12 (2 2  1) D220 .

(154)

The dynamic-stiffness coefficient for 2nd DOF may be obtained using Eq. (92) as

S 22 (a0 )   12 D220  12 a0 D220

1 H (2)1 (a0 )  H (2)1 (a0 ) . H  (a0 ) ( 2)

2

2

2

36

(155)

ACCEPTED MANUSCRIPT

 In this example n=4 is selected. The dynamic-stiffness coefficient S 22 (a0 ) may be written in a

distinct form as S 22 (a0 )  K 22 k (a0 )  i a0 c (a0 ) .

(156)

Eqs. (155) and (156) as given below

H 3( 2 ) (a 0 )  H 5( 2 ) (a 0 )  1   . k (a 0 )  ia 0 c(a 0 )   1  a 0 7 H (42 ) (a 0 ) 

CR IP T

The dimensionless spring and damping coefficients are determined by setting  2  4 and using

(157)

AN US

At the second step, the governing differential equation in dynamic-stiffness in the frequency domain, Eq. (65), is solved. For this purpose, the dynamic-stiffness coefficient S 22 (h ) at a high 

but finite frequency  h  8 c s r0 is calculated from Eq. (108), and the differential equation is

M

then solved with this starting value. Finally, the regular part of unit-impulse response coefficient

ED

S r (t ) is numerically computed step-by-step, by using Eqs. (145) and (146) (  t  0.05 r0 c s ).

The results of the present method are compared with the results of the SBFEM [27] in Figures 11

PT

and 12. In Figure 12, the unit-impulse response coefficient and time are normalized as

S r (t ) (Gcs r0 ) where t  tc s r0 denotes the dimensionless time variable. Obviously, there are

AC

CE

very good agreements between the results of the present method and the results of the SBFEM.

5.3. Anti-plane motion of circular cavity embedded in full-plane In the third example, the anti-plane motion of a circular cavity of radius r0 embedded in a

full-plane is considered as shown in Figure 13. In this example, the propagation of each mode is one-dimensional and can be explained analytically. The Poisson’s ratio of the problem’s material

37

ACCEPTED MANUSCRIPT

is   13 . The anti-plane displacement uˆ 0 ( ) with exciting frequency  is applied to the structure-medium interface as shown in Figure 13. In this example, due to symmetry, only one-quarter of the structure-medium interface is

coefficients are calculated using Eq. (47) and (48) as

 A  0.0406161867 ,

CR IP T

discretized using 4 three-node elements of equal length. For the 1st DOF the redistribution

(158a)

11

ς A  0.0202812947 .

(158b)

11

 2 uˆ A  1.03873424  uˆ A  n 2 uˆ A  ( 3 ,

3 ,

3

AN US

Then, the governing differential equation for displacement is written as

 cs

) 2 uˆ A  0 . 3

(159)

M

The analytical solution of this example may be obtained using Eq. (92). The solution of this example is again computed using Eq. (67). The starting point of

ED

solving the governing differential equation for dynamic-stiffness is achieved by using the asymptotic expansion of dynamic stiffness at  h  3c s r0 (see Eq. (108)). The governing differential equation for dynamic stiffness and S11 (h ) may be obtained as

PT



  S ii,  15.0875693 S ii  0.0018283 S ii

CE

2

 0.0000248  0.0662797n 2  0.06628 2  0

AC

and

(160)

S11 (3)  0.032160099  0.20159636i ; for n  0 S11 (3)  0.035842319 + 0.19054972i ; for n  1

38

(161)

ACCEPTED MANUSCRIPT

The results of the present method are compared with the exact solution [14, 25] for the mode numbers of n  0 and n  1 , in Figures 14, 15, respectively. Again, excellent agreement is

5.4. In-plane motion of circular cavity embedded in full-plane

CR IP T

obtained for the results of the present method and those of the exact solution.

The fourth example demonstrates a circular cavity embedded in a homogeneous fullplane as plotted in Figure 16. Again, the Poisson’s ratio of the problem’s material is   13 . A

AN US

translational displacement uˆ 0 ( ) is applied to the structure-medium interface (see Figure 16). Because of symmetry, only one-quarter of the structure-medium interface is discretized using 4 three-node elements with equal length (see Figure 17). This example is a 2D problem which is solved using Eqs. (94)-(107). To this end, the coefficients are calculated for the 1st and 2nd DOFs

M

as shown in Table 1. The dynamic-stiffness coefficient may be normalized using the shear modulus G as

2

(162a) (162b)

PT

  0 cos 2 ( )d   ,

ED

S  (a0 )  Gk (a0 )  ia0 c(a0 ) ,

CE

with dimensionless spring coefficient k (a0 ) and damping coefficient c(a 0 ) . The results of the present method are compared with analytical solution [1, 14] in Figure

AC

18. As may be observed from Figure 18, there are excellent agreements between the results of the present method and analytical solutions.

5.5. In-plane motion of a semi-circular rigid foundation

39

ACCEPTED MANUSCRIPT

In the fifth example, the vertical in-plane motion of a semi-circular rigid foundation in a half-plane is examined (see Figure 19). As no analytical solution of this problem is available, comparison of the results of the present method is made with the BEM [33] and the finite element multi-cell cloning [34] results. In this problem, the Poisson’s ratio of the material is

CR IP T

  14 .

In this example, only the boundary of semi-circular foundation is discretized using 8 three-node elements. This example is solved by implementing Eqs. (94)-(107), by which the

AN US

coefficients are obtained for 17th and 18th DOFs as illustrated in Table 1. The dynamic-stiffness coefficient is normalized with the shear modulus G as given by

S  (a0 )  Gk (a0 )  ia0 c(a0 ) , 

 2

.

(163b)

M

  0 cos 2 ( )d 

(163a)

The results of the present method are plotted in Figure 20 with good agreements with the

PT

ED

results of other numerical methods.

5.6. Concrete gravity dam

CE

The last example is devoted to dynamic behavior of a concrete rigid dam with vertical upstream face as shown in Figure 21. The geometry and mechanical properties of the problem

AC

are as follow: the depth of the reservoir is h  180 m, the pressure wave velocity is c  1438.7 m/s, and the fluid density is   1000 kg/m3. The interface between the concrete

gravity dam and the reservoir is discretized using 2 three-node elements (see Figure 21). The dam is excited by the horizontal acceleration of Ricker wavelet as shown in Figure 22. The results of the present study are shown in Figure 23. As may be seen from Figure 23, excellent 40

ACCEPTED MANUSCRIPT

agreement is obtained compared to the results of exact solution [19, 48]. To present the performance of the proposed method for practical earthquake loads, the El Centro earthquake ground motion is applied to the dam’s support (see Figure 24). Again, the obtained results of the present method show excellent agreement with those of analytical solution [19, 48], as depicted

CR IP T

in Figure 25.

6. Conclusions

In this paper, a new semi-analytical method has been summarized and developed for

AN US

diagonalization of dynamic-stiffness matrix in the frequency domain, and proposing a new

global nonreflecting BC. In this method, only the problem’s boundary at the free surface of a half-space should be discretized. This method has leaded to a set of decoupled Bessel’s

M

differential equations in the frequency domain. The Bessel’s differential equations have then been used to derive the dynamic-stiffness matrix in the present method. This method is based on

ED

substructure method, in which the dynamic-stiffness matrices have been derived for both statics and dynamics. Using the present method has caused two principal first-order differential

PT

equations of the problem become diagonal. These principal equations include interaction force-

CE

displacement relationship and governing differential equation in dynamic stiffness. In other words, these two differential equations for each DOF have become independent from other

AC

DOFs of the problem’s domain. The force-displacement relationship may be regarded as global nonreflecting BC in the medium-structure interface. The high frequency asymptotic expansion for the dynamic-stiffness matrix and unit-impulse response coefficient of an unbounded media has been presented using the proposed method. Six benchmark examples have been successfully carried out by implementing the proposed method. The results of the present method may be

41

ACCEPTED MANUSCRIPT

directly used in the SSI and dam-reservoir system problems. The coupled SSI problems may be decoupled using the proposed method, which leads to solving the structure and soil media directly, for each DOF. These objects are currently being followed by the authors, and the results may appear soon.

CR IP T

It is worthwhile remarking that the main limitation of the proposed method is its

application for layered media, where it is required to discretize more boundaries. In other words, the number of DOFs may be increased for layered problems, which is a common case for

boundary-wise methods. However, it should be noted that the increasing in computation cost due

AN US

to more DOFs is somehow compensated by low computational cost thanks to the proposed

diagonal matrices. Finally, as already indicated and discussed, the unbounded domain far from the soil-structure interface is usually assumed to be isotropic and homogeneous. For more

M

practical cases of anisotropic and non-homogeneous media, the proposed method should be

ED

further developed.

References

PT

[1] Wolf JP, Song C. Finite-element modelling of unbounded media. John Wiley & Sons:

CE

Chichester, 1996.

[2] Sommerfeld A. Partial differential equations in physics. Academic Press: New York, 1949.

AC

[3] Kausel E. Local transmitting boundaries. Journal of Engineering Mechanics 1988; 114, 1011–1027.

[4] Givoli D. Non-reflecting boundary conditions. Journal of Computational Physics 1991; 94, 1–29.

42

ACCEPTED MANUSCRIPT

[5] Givoli D. High-order local non-reflecting boundary conditions: a review. Wave motion 2004; 39, 319–326. [6] Tsynkov SV. Numerical solution of problems on unbounded domains. A review. Applied Numerical Mathematics 1998; 27, 465–532.

Heidelberg, 1984.

CR IP T

[7] Brebbia CA, Telles JCF, Wrobel LC. Boundary Element Techniques. Springer Berlin

[8] Kausel E. Thin‐layer method: Formulation in the time domain. International Journal for Numerical Methods in Engineering 1994; 37, 927–941.

AN US

[9] Keller JB, Givoli D. Exact non-reflecting boundary conditions. Journal of Computational Physics 1989; 82, 172–192.

[10] Song C, Wolf JP. The scaled boundary finite-element method—alias consistent infinitesimal

M

finite-element cell method—for elastodynamics. Computer Methods in Applied Mechanics and Engineering 1997; 147, 329–355.

ED

[11] Prempramote S, Song C, Tin-Loi F, Lin G. High-order doubly asymptotic open boundaries

79, 340–374.

PT

for scalar wave equation. International Journal for Numerical Methods in Engineering 2009;

CE

[12] Smith WD. A nonreflecting plane boundary for wave propagation problems. Journal of Computational Physics 1974; 15, 492–503.

AC

[13] Liao ZP, Wong HL. A transmitting boundary for the numerical simulation of elastic wave propagation. International Journal of Soil Dynamics and Earthquake Engineering 1984; 3, 174–183.

43

ACCEPTED MANUSCRIPT

[14] Song C, Bazyar MH. A boundary condition in Padé series for frequency-domain solution of wave propagation in unbounded domains. International Journal for Numerical Methods in Engineering 2007; 69, 2330–2358. [15] Bazyar MH, Song C. A continued-fraction-based high-order transmitting boundary for wave

CR IP T

propagation in unbounded domains of arbitrary geometry. International Journal for Numerical Methods in Engineering 2008; 74, 209–237.

[16] Birk C, Song C. A local high-order doubly asymptotic open boundary for diffusion in a semi-infinite layer. Journal of Computational Physics 2010; 229, 6156–6179.

AN US

[17] Wang Y, Lin G, Hu Z. Novel nonreflecting boundary condition for an infinite reservoir based on the scaled boundary finite element Method. Journal of Engineering Mechanics 2015; 141, 04014150.

M

[18] Lin G, Wang Y, Hu Z. An efficient approach for frequency‐domain and time‐domain hydrodynamic analysis of dam–reservoir systems. Earthquake engineering and structural

ED

dynamics 2012; 41, 1725–1749.

[19] Wang X, Jin F, Prempramote S, Song C. Time-domain analysis of gravity dam–reservoir

PT

interaction using high-order doubly asymptotic open boundary. Computers and Structures

CE

2011; 89, 668–680.

[20] Chen D, Birk C, Song C, Du C. A high‐order approach for modelling transient wave

AC

propagation problems using the scaled boundary finite element method. International Journal for Numerical Methods in Engineering 2014; 97, 937–959.

[21] Sharma KG. Condensation of elastic half space soil stiffness matrix considering symmetry. International Journal for Numerical and Analytical Methods in Geomechanics 1985; 9, 391– 395.

44

ACCEPTED MANUSCRIPT

[22] Dominguez J. Boundary Elements in Dynamics. Computational Mechanics Publications: Southampton, 1993. [23] E. Kausel. Fundamental solutions in elastodynamics. A compendium. Cambridge University Press, 2006.

CR IP T

[24] Wolf JP. The Scaled Boundary Finite Element Method. Wiley: New York, 2004.

[25] Song C, Wolf JP. The scaled boundary finite-element method: analytical solution in

frequency domain. Computer Methods in Applied Mechanics and Engineering 2000; 164, 249–264.

AN US

[26] Wolf JP, Song C. The scaled boundary finite-element method – a primer: derivations. Computers and Structures 2000; 78, 191–210.

[27] Song C, Wolf JP. The scaled boundary finite-element method – a primer: solution

M

procedures. Computers and Structures 2000; 78, 211–225.

[28] Wolf JP. Response of unbounded soil in scaled boundary finite-element method. Earthquake

ED

Engineering and Structural Dynamics 2002; 31, 15–32. [29] Song C. Dynamic analysis of unbounded domains by a reduced set of base functions.

PT

Computer Methods in Applied Mechanics and Engineering 2006; 195, 4075–4094.

CE

[30] Deeks AJ, Wolf JP. Semi-analytical elastostatic analysis of unbounded two-dimensional domains. International Journal for Numerical and Analytical Methods in Geomechanics 2002;

AC

26, 1031–1057.

[31] Komatitsch D, Ritsema J, Tromp J. The spectral-element method, Beowulf computing and global seismology. Science 2002; 298(5599):1737–1742.

45

ACCEPTED MANUSCRIPT

[32] Y. Maday and A. T. Patera. Spectral element methods for the incompressible Navier-Stokes equations. In A. K. Noor and J. T. Oden, editors, State-of-the-art surveys on computational mechanics, chapter 3, pages 71–143. ASME, 1989. [33] Wang Y, Rajapakse RKND. Dynamics of rigid strip foundations embedded in orthotropic

CR IP T

elastic soils. Earthquake Engineering and Structural Dynamics 1991; 20, 927–947.

[34] Wolf JP, Song C. Dynamic-stiffness matrix of unbounded soil by finite-element multi-cell cloning. Earthquake Engineering and Structural Dynamics 1994; 23, 233–250.

[35] Wolf JP, Song C. Consistent infinitesimal finite-element cell method in frequency domain.

AN US

Earthquake Engineering and Structural Dynamics 1996; 25, 1307–1327.

[36] Wolf JP, Song C. Unit-impulse response matrix of unbounded medium by infinitesimal finite-element cell method. Computer Methods in Applied Mechanics and Engineering 1995;

M

122, 251–272.

[37] de Barros FCP, Luco JE. Dynamic response of a two-dimensional semi-circular foundation

1995; 14, 45–57.

ED

embedded in a layered viscoelastic half-space. Soil Dynamics and Earthquake Engineering

PT

[38] Khaji N, Khodakarami MI. A new semi-analytical method with diagonal coefficient

CE

matrices for potential problems. Engineering Analysis with Boundary Elements 2011; 35, 845–854.

AC

[39] Khodakarami MI., Khaji N. Analysis of elastostatic problems using a semi-analytical method with diagonal coefficient matrices. Engineering Analysis with Boundary Elements 2011; 35, 1288–1296.

46

ACCEPTED MANUSCRIPT

[40] Khaji N, Khodakarami MI. A semi-analytical method with a system of decoupled ordinary differential equations for three-dimensional elastostatic problems. International Journal of Solids and Structures 2012; 49, 2528–2546. [41] Khodakarami MI, Khaji N, Ahmadi MT. Modeling transient elastodynamic problems using

CR IP T

a novel semi-analytical method yielding decoupled partial differential equations. Computer Methods in Applied Mechanics and Engineering 2012; 213, 183–195.

[42] Khaji N, Mirzajani M. Frequency domain analysis of elastic bounded domains using a new semi-analytical method. Acta Mechanica 2013; 224, 1555–1570.

AN US

[43] Wolf JP. Soil–Structure-Interaction Analysis in Time Domain. Prentice-Hall: Englewood Cliffs, 1988.

[44] R. Courant, D. Hilbert, Methods of Mathematical Physics, Vol. I, Interscience, New York,

M

1953.

[45] Hibbeler RC. Structural Analysis, 8th edition, Pearson Prentice Hall, 2012.

ED

[46] Graff KF. Wave motion in elastic solids. Courier Dover Publications, 1975. [47] Li B, Cheng L, Deeks AJ, Teng B. A modified scaled boundary finite-element method for

PT

problems with parallel side-faces. Part I. Theoretical developments. Applied Ocean Research

CE

2005; 27, 216–223.

[48] Chopra AK. Hydrodynamic pressures on dams during earthquakes. Journal of the

AC

Engineering Mechanics Division 1967; 93(EM6), 205–223.

47

ACCEPTED MANUSCRIPT

A le ich

ace e-f t BC

Dir

Sid

A

CR IP T

e fac BC e Sid ann

LCO

AN US

S

M

ξ →∞

AC

CE

PT

ED

(a)

48

ξ →∞

ξ →∞

um Ne

ACCEPTED MANUSCRIPT

∞

→ ∞

ξ→

ξ

ξ

>1 >1

LCO

M

∞ → ξ

→ ∞

AN US

SS

CR IP T

ϴ

ξ

ξ =1

u0ξ u0 -u0η

ED

(b)

Figure 1. Soil-structure interface of unbounded medium. (a) Unbounded medium with Dirichlet

AC

CE

PT

and Neumann BCs without spatial discretization; (b) Cavity embedded in full space.

49

ACCEPTED MANUSCRIPT

ŷ

LCO Neumann BC

Neumann BC

û R

Diric hlet B C

Unbounded soil

C hlet B Diric

AN US

∞ Ω

ξ→

∞

ξ→

CR IP T

ξ =1

ξ =1

Soil-structure interface

AC

CE

PT

ED

M

Figure 2. Interaction force-displacement relationship on the soil-structure interface.

50

xˆ

ACCEPTED MANUSCRIPT

ξ

η 2 (ξ=1,η=+1)

nξ nη

Sξ

CR IP T

ξ

ŷ

1 (ξ=1,η=-1) Sη xˆ

LCO (ξ=0)

0

M

1 -1

AN US

(a)

2 +1

η

ED

(b)

Figure 3. (a) Schematic view of a sample 2D domain, in which the LCO is shown in local

AC

CE

PT

coordinate systems; (b) A sample two-node element.

51

η= +1

ACCEPTED MANUSCRIPT

η

ξ→

∞

CR IP T

ξ

2

η=-1

1

AN US

LCO (ξ=0)

(a)

η

η= +1

ED

ξ

1 ξ=

M

2

1

η=-1

PT

LCO (ξ=0)

the LCO

AC

CE

Figure 4. (a) Unbounded and (b) bounded media with two-node line element and the location of

52

ACCEPTED MANUSCRIPT

z η Free surface η

CR IP T

h

ξ

LCO at infinity

AN US

Fixed boundary

(a)

+


+1 0 -1

x

(b)

AC

CE

PT

ED

M

Figure 5. Example 1: (a) Geometry and Discretization of semi-infinite layer with constant depth, (b) A sample three-node element.

53

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

Figure 6. dynamic stiffness coefficient of semi-infinite layer. 54

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

Figure 7. Unit-impulse response of semi-infinite layer for (a) the first 100, and (b) the total 300 steps of the analysis.

55

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 8. Ricker wavelet traction of the first numerical example.

56

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

AC

CE

PT

ED

M

(a)

57

AN US

CR IP T

ACCEPTED MANUSCRIPT

(c)

Figure 9. The response of semi-infinite layer to Ricker wavelet traction: (a) n  5 , (b) n  10 , and

AC

CE

PT

ED

M

(c) n  15 .

58

ACCEPTED MANUSCRIPT

ŷ 2

∞

1

∞

CR IP T

xˆ LCO

Figure 10. Example 2: Truncated semi-infinite wedge and the proposed mesh including one two-

AC

CE

PT

ED

M

AN US

node element.

59

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

Figure 11. Dynamic stiffness coefficients of truncated semi-infinite wedge of the second example. (a) spring coefficient, (b) damping coefficient.

60

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 12. Regular part of unit-impulse response coefficient of the second example.

61

ACCEPTED MANUSCRIPT




r0


(ω) u0

ŷ

ξ xˆ

u

AN US

LCO

τ zξ

CR IP T




Structure-Medium Interface

AC

CE

PT

ED

M

Figure 13. Example 3: Anti-plane motion of circular cavity embedded in a full-plane.

62

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

Figure 14. Dynamic stiffness coefficient of the third example, for the mode number of n  0 . (a) real part, and (b) imaginary part.

63

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b) Figure 15. Dynamic stiffness coefficient of the third example, for the mode number of n  1. (a) real part, and (b) imaginary part. 64

ACCEPTED MANUSCRIPT




ση




uξ




-u0η

u0ξ
u0(ω)

CR IP T

r0 ŷ




u

-uη










σξ

ξ LCO

xˆ

AN US

Structure-Medium Interface

AC

CE

PT

ED

M

Figure 16. Example 4: In-plane motion of circular cavity embedded in a full-plane.

65

ACCEPTED MANUSCRIPT

LCO

AN US

CR IP T

r0

Figure 17. The LCO and the proposed mesh employing three-node elements for the fourth

AC

CE

PT

ED

M

example.

66

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

Figure 18. Dynamic-stiffness coefficient of the fourth example. (a) real part, and (b) imaginary part.

67

ACCEPTED MANUSCRIPT

ˆy

xˆ

LCO (0,0)

 1

CR IP T



AN US

 

Figure 19. Vertical motion of a semi-circular foundation embedded in a half-plane. In this example, only the semi-circular boundary is discretized by 8 three-node elements. In this figure,

AC

CE

PT

ED

M

the solid circles indicate the location of beginning and ending nodes of each element.

68

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

69

ACCEPTED MANUSCRIPT

Figure 20. Dynamic-stiffness coefficient of the fifth example: (a) spring coefficient, and (b)

AC

CE

PT

ED

M

AN US

CR IP T

damping coefficient.

70

ACCEPTED MANUSCRIPT

z

Rigid dam

Dam-Reservoir interface

x

CR IP T




Reservoir

AC

CE

PT

ED

M

AN US

Figure 21. Example 6: Rigid dam with a semi-infinite reservoir of constant depth.

71

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 22. Ricker wavelet load of the last numerical example.

72

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 23. Hydrodynamic pressure on the dam’s heel under Ricker wavelet acceleration.

73

AC

CE

PT

ED

M

Figure 24. El Centro earthquake time history.

AN US

CR IP T

ACCEPTED MANUSCRIPT

74

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 25. Hydrodynamic pressure on the dam’s heel due to El Centro earthquake.

75

ACCEPTED MANUSCRIPT

Table 1. Coefficients of the present method for the fourth and fifth examples

th

4 example

th



A

0.035365885 0.0421817261 7

0.020308049 0.0609241495 8

B

1.00323117

1.00379045

E

M ED PT CE AC

76

Z

0.0804536475

0.0871508770

P

T

1.00000888

0.004215912 50

0.57957067 3

4

1.00000521

0.003227840 00

0.75248231 3

3

AN US

5 example



CR IP T

Example