- Email: [email protected]
A Coupling Procedure of Finite Element and Scaled Boundary Finite Element Methods for Soil–Structure Interaction in the Time Domain1
C H A P T E R
6 A Coupling Procedure of Finite Element and Scaled Boundary Finite Element Methods for SoileStructure Interaction in the Time Domain1 Junyi Yan, Chuhan Zhang, Feng Jin Department of Hydraulic Engineering, Tsinghua University, Beijing
O U T L I N E 6.1 Introduction
118
6.2 Motion Equations of Coupling System
120
6.3 Realization and Model Approximation 6.3.1 Diagonalization 6.3.2 Decomposition 6.3.3 Truncation 6.3.4 Realization and Model Reduction
121 121 122 123 125
6.4 Evaluation of Interaction Forces
127
6.5 Numerical Verification 6.5.1 Accuracy 6.5.1.1 Spherical Cavity Embedded in Full Space 6.5.1.2 Dynamic Boussinesq Problem
127 127 127 129
International Journal for Numerical Methods in Engineering 2004; 59: 1453e71. Ó 2004 John Wiley & Sons Ltd. 1
Seismic Safety Evaluation of Concrete Dams http://dx.doi.org/10.1016/B978-0-12-408083-6.00006-4
117
Copyright Ó 2013 Tsinghua University Press. Published by Elsevier Inc. All rights reserved.
118
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
6.5.1.3 Massless Square Plate Resting on Half-Space 6.5.2 Efficiency
131 133
6.6 Conclusions
137
Acknowledgments
137
References
138
6.1 INTRODUCTION Dynamic analysis of soilestructure systems is an important topic in earthquake engineering. During recent decades many numerical procedures have been developed to deal with this kind of problem. Usually, the structure including a finite region of adjacent soil is modeled by the finite element (FE) method while the remaining unbounded soil is represented by boundary conditions on its interface with the structure, which leads to a model with a finite number of degrees of freedom. The boundary conditions describe the forceemotion relationship on the soilestructure interface, whose rigorous form is global in space and time. The boundary element (BE) method has become a popular procedure for constructing such a representation of boundary conditions. Another kind of procedure is categorized as various transmitting boundaries connected with finite elements to prevent a spurious wave reflection from the artificial boundaries. This procedure is usually constructed to be local in space and time to reduce the computational effort. Each of these two representations has its own merits and disadvantages. The BE method is well known as a powerful procedure for modeling the unbounded medium since the radiation condition is satisfied automatically as a part of the fundamental solution. However, the fundamental solution is not always available or it can be very complicated even if it exists. The scaled boundary finite element (SBFE) method [1,2], however, combines the advantages of the FE method and the BE method. Only the boundary is discretized and no fundamental solution is required. As a semi-analytical algorithm, it is exact in the radial direction and converges to an exact solution in the circumferential directions when the FE size approaches infinity. The SBFE method can be applied to static, diffusion and wave propagation problems in both frequency and time domain [3]. Therefore, developing a coupling procedure of FE and SBFE in the time domain to perform linear and nonlinear seismic analysis of the unbounded soilestructure system is the main objective of this chapter.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.1. INTRODUCTION
119
For numerical algorithms in the time domain such as a BE or SBFE procedure, the dynamic property of the unbounded medium is normally represented by a unit-impulse response matrix and the interaction-force vector is evaluated by the convolution integrals of the unit-impulse response matrix and the motion vector. A large computational effort (proportional to the square of the number of time stations) and storage requirement result in an expensive cost in the analysis of the practical problems with large degrees of freedom on the soilestructure interface. To reduce the computational effort, Wolf and Motosaka [4] presented two approximate techniques to evaluate the interaction forces recursively. Later, the concept of linear system theory was introduced by Paronesso and Wolf [5] for constructing rational functions in the frequency domain. The unit-impulse response matrices using expansion of a series of Legendre polynomials are obtained approximately, by which the equivalent stiffness, damping and mass matrices of the unbounded medium can be produced. Furthermore, a system realization algorithm implemented in the time domain is employed by Paronesso and Wolf [6] to calculate the interaction forces with a system of linear equations instead of convolution integrals. From physical aspects of the unbounded soil, a piecewise linear function is used by Zhang et al. [7] to describe the acceleration unitimpulse response matrix by which the implementation of the convolution integrals is simplified. However, the accuracy of this scheme is uncontrollable, which reduces the primary advantage of the SBFE method. In this chapter, a new coupling procedure of FEeSBFE is developed for three-dimensional (3D) dynamic analysis of unbounded soilestructure interaction in the time domain. The SBFE is used to represent the unbounded soil, while the bounded structure including a portion of the adjacent soil is modeled by FE. The concept of linear system theory [8] is applied to calculate the soilestructure interaction forces instead of using the time-consuming convolution integrals. To calculate the unit-impulse response matrix using the SBFE method, the convolution integral and Lyapunov function must be treated at each time station, resulting in enormous computational effort. Since a partial sequence of the unitimpulse response matrix is sufficient to represent the main character of the unbounded medium, only the terms before a cut-off time station are required and numerical cost can thus be reduced. In addition, a powerful system realization algorithm called FastERA [9] is employed to further reduce the computational effort. The accuracy of the presented procedure can be controlled by prescribed parameters. Thus, the advantage of high accuracy of the SBFE method is retained. Numerical examples of foundation response show good agreement compared with previous results and high efficiency compared with direct convolution integral implementation. The new scheme presented herein can be applied to the
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
120
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
dynamic analysis of unbounded soilestructure interaction of large-scale systems.
6.2 MOTION EQUATIONS OF COUPLING SYSTEM To demonstrate the accuracy and efficiency of the presented numerical scheme, only the unbounded soilestructure coupling system subjected to external time-varying loading is considered here. As shown in Fig. 6.1, the generalized structure including the real structure and a finite region of adjacent soil, either of which may exhibit nonlinear behavior, is modeled by FE, while the unbounded soil, which is assumed to be homogeneous and linear without losing its generality, is represented by SBFE. The equation of motion of the structure is written in the time domain as €s u Css Csb u_ s Kss Ksb us ps Mss Msb þ þ ¼ €b u Mbs Mbb Cbs Cbb Kbs Kbb ub rb u_ b (6.1) and the contribution of the unbounded soil can be represented by the interaction forces as Zt rb ðtÞ ¼
Gbb ðt sÞ€ ub ðsÞds
(6.2)
0
where M, C and K are the mass, damping and stiffness matrices of the structure, respectively; Gbb (t) is the acceleration unit-impulse response matrix of the unbounded soil which is determined via the SBFE method; u, u_ and u¨ are displacement, velocity and acceleration vectors, respectively, p is the external load vector and r is the interaction force vector.
FIGURE 6.1 Problem definition of dynamic unbounded soilestructure-interaction analysis: (a) coupled system; (b) unbounded soil.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
121
The subscripts b and s denote, respectively, the degrees of freedom associated with the soilestructure interface and the remaining ones associated with the structure. In general, the interaction forces can be expressed as a function of either the displacements or the accelerations. To avoid singularity, acceleration formulation is used herein. Using the time-integral scheme, the seismic response of the coupling system can theoretically be obtained by solving Eqs (6.1) and (6.2).
6.3 REALIZATION AND MODEL APPROXIMATION As a rigorous representation, interaction forces are evaluated by convolution integrals as shown in Eq. (6.2), which is exact but very expensive. The cost of the direct convolution integral implementation is proportional to the square of the number of time stations, and requires enormous storage capacities. Herein, concepts of linear system theory such as realization and model reduction are applied and the computational effort can be reduced enormously. To ensure that the present scheme is efficient and stable, diagonalization, decomposition and truncation are performed before the system realization. For the sake of conciseness, the subscript b of Gbb(t) is dropped in the following sections.
6.3.1 Diagonalization As the number of the degrees of freedom associated with the soilestructure interface is usually large in practical analysis, it is very inefficient and even unrealistic to apply the realization algorithm directly to the total matrix G(t). Therefore, the so-called diagonalization procedure is performed before the realization is conducted. Details of the diagonalization procedure are presented by Paronesso and Wolf [10,11]. With a time-independent transformation matrix T, the symmetric matrix G(t) of order n n is converted into the corresponding diagonal matrix Gd(t). The relationship between G(t) and Gd(t) can be expressed as GðtÞ ¼ TGd ðtÞTT
(6.3)
Premultiplying the transformation T and T , respectively, on rdb ðtÞ and u¨b(t), i.e., T
€ db ðtÞ ¼ TT u € b ðtÞ u
(6.4a)
rb ðbÞ ¼ Trdb ðtÞ
(6.4b)
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
122
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
and substituting Eqs (6.3) and (6.4) into Eq. (6.2) yields Zt rdb ðtÞ
Gd ðt sÞ€ udb ðsÞds
¼
(6.5)
0
€ db ðtÞ and rdb ðtÞ are, where u rb(t), which are defined
respectively, the transformed vector of u¨b(t) and by Eqs (6.4a) and (6.4b). The transformation process from Eq. (6.2) to Eq. (6.5) is called diagonalization. In Eq. (6.5), each term of Gd(t) can be addressed separately, which means that a multi-input multi-output (MIMO) system can be performed as a number of single-input single-output (SISO) systems. Therefore, approximation and realization are discussed in terms of the individual scalar term Gd(t) in the following sections.
6.3.2 Decomposition In view of the fact that the Hankel matrix realization algorithm is stable and robust when the corresponding Markov parameter sequence has the convergence-to-zero property [12], Gd(t) is decomposed as shown in Fig. 6.2, yielding Gd ðtÞ ¼ Cd HðtÞ þ Kd tHðtÞ þ Gm ðtÞ
(6.6)
with coefficients Kd and Cd chosen in such a way to achieve Gm(t) converging to zero as time increases to infinity, where H(t) is the Heaviside-step function. Then, Eq. (6.5) can be rewritten as Zt rdb
¼
Kd udb ðtÞ þ
Gm ðt sÞ€ udb ðsÞds
Cd u_ db ðtÞ þ
(6.7)
0
FIGURE 6.2
Decomposition of a typical term of the diagonalized unit-impulse response
matrix.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
123
where Kd, Cd and Gm (t) are diagonal matrices with the terms as Kd, Cd and Gm(t), respectively. The computation of the first two terms of the right-hand side of Eq. (6.7) is straightforward and is not discussed further. The evaluation of the third term expressed as a convolution integral can be accomplished with ease when each of scalar convolution integrals of Gm(t) and u€db ðtÞ is properly calculated. The existence of the decomposition is guaranteed by the fact that the second derivative of the acceleration unit-impulse response matrix is equal to the regular part of the displacement unit-impulse response matrix which converges to zero in a physical sense when time increases to infinity. The coefficient matrix Kd in Eq. (6.7) corresponds to the static-stiffness matrix of the unbounded medium K0. The relationship between these two matrices can be expressed with the transformation defined by Eq. (6.3) as K0 ¼ TKd TT
(6.8)
The static-stiffness matrix K0 can be obtained by solving an algebraic Riccati equation as mentioned by Wolf and Song [3]. The coefficient matrix Cd cannot be obtained directly, and is achieved approximately in Section 6.3.3.
6.3.3 Truncation Considering the convergence-to-zero property, Gm(t) is approximated as d G ðtÞ Cdc Kd t; 0 t < tc : (6.9) Gm ðtÞz 0; t tc where tc is the cut-off time and the truncated function Gm(t) is used. Coefficient Cdc can then be calculated from Eq. (6.9) with the known values of Kd and Gd(t). Whether the truncated function of the unit-impulse response matrix satisfies the prescribed tolerance depends on the cut-off time tc. For determination of the cut-off time tc the following criterion is proposed, which has a proven applicability. In the process of calculating the unit-impulse matrix via SBFE with a prescribed tolerance parameter "tol", the normalized Frobenius norm F(t)/F(0) is computed at each time station, where 11=2 0 n X n X 2 € ðtÞ A G (6.10) FðtÞ ¼ @ ij i¼1 j¼1
Once the inequality FðtÞ=Fð0Þ tol
(6.11)
is satisfied, the time station is chosen to be the cut-off time tc.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
124
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.3 Finite element meshes of soilestructure interfaces of typical excavations in half-space: (a) prismatic excavation; (b) hemispherical excavation.
FIGURE 6.4 Variation of the normalized Frobenius norms against dimensionless time.
To demonstrate the decaying behavior of the normalized norm F(t)/F(0), two typical cases of excavated foundation are analyzed. One is the half-space with a square prismatic excavation with length 2r0 and depth r0, the other is the half-space with a hemispherical excavation with radius r0. Figure 6.3 shows the FE meshes of one quarter of the soilestructure interfaces of these two cases. The normalized Frobenius norms F(t)/F(0) are plotted in Fig. 6.4 as functions of the dimensionless time t ¼ tcp =r0 ; where cp is the velocity of the dilatational wave. It is evident that the values of F(t)/F(0) converge to zero monotonously.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
125
Complete description of a system can be expressed by an infinite Hankel matrix, which is the starting point of the Eigensystem realization algorithm (ERA) [13]. When truncation is performed, the infinite Hankel matrix can be written as 2 3 ~m 2 / G ~ m kc 1 0 / ~m 1 G G 6 7 6 7 ~m 3 / ~m 2 6 G 7 G 0 0 . 6 7 6 7 6 « « 1 « « « 7 7 (6.12) HNN ð0Þ ¼ 6 6 m 7 ~ kc 1 6G 7 0 / 0 0 / 6 7 6 7 6 7 0 0 / 0 0 / 4 5 « « / « « 1 ~ m ðkÞ is the Markov parameter sequence sampled from in which G m G (t) with Dt and the total number of time stations kc equals tc/Dt. As the order of the system is equal to the rank of the Hankel matrix theoretically, the Hankel matrix can be reduced dramatically as 2 6 6 6 Hðkc 1Þðkc 1Þ ð0Þ ¼ 6 6 6 4
~m 1 G ~m 2 G
~m 2 G ~m 3 G
«
«
1
0
/
~ m kc 1 G
/ /
~ m kc 1 3 G 7 7 0 7 7 7 « 5
(6.13)
0
Since the order of the Hankel matrix is reduced, the following realization process involving eigensystem decomposition of the inner product of the Hankel matrix is more feasible and economical.
6.3.4 Realization and Model Reduction Considering the third term of the right-hand side of Eq. (6.7) in which Gm(t) has been diagonalized, each SISO system can be expressed as Zt rdb ðtÞ
¼
€db ðsÞds Gm ðt sÞu
(6.14)
0
In linear system theory, such convolution integral expression, called inputeoutput description, can be transformed into an equivalent system of linear equations, named state-variable description. In a discrete-time
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
126
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
system, this state-variable description is represented by finite-difference equations of first order at distinct time station k as ~ e xe ðkÞ þ B ~ eu €db ðkÞ xe ðk þ 1Þ ¼ A
(6.15a)
~ e xe ðkÞ þ D ~ eu €db ðkÞ rbd ðkÞ ¼ C
(6.15b)
where xe (t) is called the state-variable vector, the coefficient matrix ~ e; D ~ e; C ~ e g is named the realization of the unit-impulse ~ e; B quaternion fA m ~ response matrix G ðkÞ which is sampled from Gm(t) with time step Dt; and the subscript e corresponds to the specified SISO system. The ~ m ðkÞ using Hankel ~ e; D ~ e; C ~ e g can be calculated from G ~ e; B quaternion fA matrix realization algorithms. To reduce computational effort and eliminate the effects of unavoidable errors such as numerical round-off, it is desirable to construct a reduced model of the original system. This is generally achieved by singular value decomposition because the rank of the Hankel matrix, equivalent to the order of the system, can be evaluated numerically by non-zero singular values. Contributions associated with the smallest singular values below a prescribed tolerance threshold are neglected, thus the approximate model with reduced order is realized according to the remaining singular values. Herein, a powerful realization algorithm named FastERA is employed, in which only partial eigensystem decomposition of the inner product of the Hankel matrix is computed, rather than the full singular value decomposition of the Hankel matrix as processed in ERA. The reduced order of the system is determined with a prescribed allowance parameter "tole". Define a threshold of eigenvalue lc, which satisfies , !1=2 N X li (6.16) tole ¼ lc i¼1
where N and li denote, respectively, the order and the eigenvalue of the inner product of the Hankel matrix. The sum of all the eigenvalues in Eq. (6.16) can be calculated through the trace of the above inner product of the Hankel matrix. When the threshold lc is determined by Eq. (6.16) and the upper bound of the largest eigenvalue lmax is estimated with Gerschgorin disk theorem, the partial eigensystem decomposition corresponding to the eigenvalues between lc and lmax can be accomplished by the Sturm sequence method. The reduced ~ e and D ~ e; C ~ e can be ~ e; B order is then determined and the matrices A obtained. ~ e is Since the maximum magnitude of the eigenvalues of matrix A smaller than 1, the procedure is stable.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
127
6.5. NUMERICAL VERIFICATION
6.4 EVALUATION OF INTERACTION FORCES When each element of matrix Gm bb ðtÞ is approximated and realized, Eq. (6.15) is assembled and the interaction forces shown in Eq. (6.7) are evaluated in the discrete-time system as ~ ~ ud ðkÞ xðk þ 1Þ ¼ AxðtÞ þ B€ b
(6.17a)
~ ~ ud ðkÞ € db ðkÞ þ CxðtÞ rdb ðkÞ ¼ Kd udb ðkÞ þ Cdb u þ D€ b
(6.17b)
where the state-variable vector x is assembled from xe; and the coefficient ~ and D ~ B; ~ C ~ are block diagonal matrices with respect to matrices A; ~ e and D ~ e; C ~ e , respectively. Taking Eq. (6.4) into account ~ e; B the terms A gives g
g
g
€ b ðkÞ þ yðkÞ rb ðkÞ ¼ Kbb ub ðkÞ þ Cbb u_ b ðkÞ þ Mbb u
(6.18)
in which ~ ~ T ; YðkÞ ¼ TCxðkÞ Kbb ¼ TKd TT ; Cbb ¼ TCdc TT ; Mbb ¼ TDT (6.19) g
g
g
Discretizing Eq. (6.1) with time step Dt and substituting Eq. (6.18) into Eq. (6.1), the seismic response of the coupling system can be solved using a time-integral scheme such as the Newmark method.
6.5 NUMERICAL VERIFICATION The coupling procedure proposed in this chapter can be used in the analysis of the soilestructure interaction system subjected to direct external loading as well as seismic free-field input. Only external loading is taken into account in the following examples. The seismic free-field formulation will be reported later.
6.5.1 Accuracy This section verifies the numerical accuracy of the presented procedure compared with previous results by other researchers. To demonstrate the efficiency of the presented procedure, solutions calculated by a conventional convolution method are also presented and compared. 6.5.1.1 Spherical Cavity Embedded in Full Space A spherical cavity embedded in the full space subjected to a uniform normal pressure is analyzed. As a one-dimensional problem, its analytical
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
128
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.5 Finite element mesh of one octant of a spherical shell.
solution is available [3]. The rounded triangular pressure pulse acting on the wall of cavity is represented by
8 < p0 1 cos 2pt ; 0 t t 0 2 t0 pðtÞ ¼ (6.20) : 0; t > t0 where t0 ¼ 3.46r0/cp, in which r0 is the radius of the cavity and cp denotes the velocity of the dilatational wave. Owing to the symmetry of the problem, only one octant of the cavity is considered. As shown in Fig. 6.5, the spherical shell with thickness r0, a portion of the soil adjacent to the wall, is treated as the structure and discretized by three 20-node isoparametric hexahedral elements and the soilestructure interface is discretized by three eight-node isoparametric quadrangular elements. The total number of degrees of freedom for the whole system is 78, of which 33 are on the soilestructure interface. In Fig. 6.5, the nodes on the soilestructure interface are denoted by black circles and the remaining nodes of the structure by blank circles. For comparison with previous results, the dimensionless time is defined as t ¼ tcp =r0 ; and the dimensionless radial displacement is uðtÞ ¼ u0 ðtÞKN =ð4pr20 p0 Þ; where u0 is the radial displacement of the wall and KN ¼ 16pGr0 is the static-stiffness coefficient with shear modulus G. Poisson’s ratio is selected as 0.25 here. The results of the present method with two pairs of tolerance parameters are plotted in Fig. 6.6 (in which the first and the second parameters denote tol and tole, respectively) together with those of the convolution method and the analytical solution. Excellent agreement is observed, especially between the solutions of the convolution method and the present procedure, with smaller tolerance parameters. It is worth noting that the present procedure converges to the convolution method theoretically when tolerance parameters approach infinitesimal values.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
FIGURE 6.6
129
Displacement response of a spherical cavity embedded in full space.
6.5.1.2 Dynamic Boussinesq Problem A 3D, homogeneous, isotropic, elastic half-space subjected to a concentrated normal harmonic loading acting on its surface, named the dynamic Boussinesq problem, is presented. Owing to symmetry, only one-quarter of the half-space is modeled, as shown in Fig. 6.7. Similarly, one octant of a spherical portion of the soil is modeled as the structure and discretized by seven 20-node hexahedral elements, and the structureemedium interface is discretized by three eight-node quadrangular elements similar to the one shown in Fig. 6.5.
FIGURE 6.7
Finite element mesh of one quarter of half-space.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
130
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.8 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.2): (a) real part; (b) imaginary part.
For comparison with previous solutions, the steady-state response of the vertical displacement at point A (Fig. 6.7) is calculated in the time domain and transformed into the form of frequency domain. The dimensionless vertical displacement is defined as uz ¼ uz Gr=P and the dimensionless frequency a0 ¼ ur/cs, where G is the shear modulus and cs is the velocity of the shear wave. Comparisons between the results are drawn in Figs 6.8 and 6.9 for values of Poisson’s ratio of 0.2 and 0.4, respectively. Both the parameters tol and tole are selected as 0.005. Again, good agreement is seen. As a whole, the results obtained by the present procedure and the convolution method are close to those of Kobayashi and Nishimura [15]. Once
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
131
FIGURE 6.9 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.4): (a) real part; (b) imaginary part.
again, excellent agreement between the solutions of the convolution method and the present procedure is evident. As a typical example, Fig. 6.10 demonstrates that the results calculated from the present procedure are close to those of the convolution method when smaller tolerance parameters tole are selected. Apparently, in contrast to that of low frequency, the response of high frequency is more sensitive to the parameter tole. 6.5.1.3 Massless Square Plate Resting on Half-Space As a 3D problem, the response of a massless square plate resting on half-space subjected to uniform normal harmonic pressure is analyzed. Herein, the stiffness of the plate is assumed to be zero, which
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
132
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.10 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.2): (a) real part; (b) imaginary part.
corresponds to a uniform pressure acting on the surface of the half-space directly. Similar to the previous example, one-quarter of the half-space is modeled as shown in Fig. 6.11, where the shadowed square region denotes the area of the imposed loading. A square prism portion of the soil is modeled as the structure and discretized by 32 20-node hexahedral elements with 553 degrees of freedom, and the soilestructure interface is discretized by 12 eight-node quadrangular elements with 129 degrees of freedom as shown in Fig. 6.3(a). Poisson’s ratio is selected as {1/3} here. Similarly, the steady-state response calculated in the time domain is transformed into the form of frequency domain. The vertical
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
133
FIGURE 6.11 Finite element mesh of one quarter of half-space.
FIGURE 6.12 Vertical compliances of a square foundation under a uniform stress distribution (n ¼ 1/3).
dimensionless displacement amplitude is defined as uz ¼ uz Gb=P and the dimensionless frequency a0 ¼ ub/cs with the length of the square plate 2b and the amplitude of total load P. Solutions of the presented method with tole ¼ 0.005 and tole ¼ 0.01 are illustrated in Fig. 6.12 compared with those given by Wong and Luco [16]. Good agreement is also achieved.
6.5.2 Efficiency This section compares the efficiency of the presented procedure with the direct convolution integral method. Contributions of each
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
134
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.13 Comparison of numerical cost of different method: (a) total central processing unit (CPU) time cost; (b) numerical cost used in evaluating interaction forces.
approximation strategy in numerical cost savings are also demonstrated. As a typical example, the dynamic Boussinesq problem with a Poisson’s ratio of 0.2, as mentioned in Section 6.5.1.2, is addressed herein. The relationship of the numerical cost to the total computational time step is analyzed. It is obvious that the total computational effort spent in the convolution method is proportional to the square of the total time step, while that of the present procedure increases linearly, as demonstrated in Fig. 6.13(a). The computational effort expended in evaluating the interaction forces varies in the same way, as shown in Fig. 6.13(b). When the number of total time step is large, as encountered in practical seismic analysis, the presented procedure is far more economical.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
135
FIGURE 6.14 Numerical cost structure: (a) convolution method; (b) present procedure. SBFE: scaled boundary finite element.
Compared with the convolution method, an additional computation is needed to accomplish the realization implementation in the present procedure. As shown in Fig. 6.14(b), this part of the cost (denoted by Realization) is trivial, because the partial sequence of the Markov parameters is usually quite short, which implies that the order of the Hankel matrix is low. The benefit of the model reduction is demonstrated as follows. For the analyzed problem with 38 degrees of freedom on the soile structure interface, the unbounded medium can be represented by 741 SISO systems in a diagonalization sense (see Section 6.3.1). With the prescribed parameter tol ¼ 0.005, the order of each system is
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
136
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.15
Distribution of the single-input single-output (SISO) systems versus the
reduced order.
determined to be 41. When the model reduction is performed, the order of the system can be reduced dramatically. Figure 6.15 shows the distribution of the SISO systems versus the reduced order with different values of tole. Effects of the parameter tole on the computational cost are illustrated in Fig. 6.16, where the numerical cost used in realization implementation is denoted by Realization and that spent in solving the final coupling equations is termed Solving. Note that when tole is selected as 0.005, the numerical cost for these two aspects can be
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.6. CONCLUSIONS
137
FIGURE 6.16 Effects of the parameter "tole" on computational effort. CPU: central processing unit.
cut down by more than 80% compared with that without model reduction.
6.6 CONCLUSIONS A coupling procedure of FE and SBFE methods is developed in the time domain for 3D dynamic analysis of unbounded soilestructure interaction. The dynamic response of the structure and adjacent soil can be evaluated considering the effects of the surrounding unbounded soil. To enable the application of SBFE to practical engineering analysis, several techniques are used to achieve computational efficiency, such as diagonalization, sequence decomposition and truncation, system realization and model reduction. Furthermore, the balance between accuracy and efficiency can be controlled by prescribed parameters, and thus the main advantage of SBFE is retained. The accuracy and efficiency of the new scheme are demonstrated by numerical examples. The coupling procedure presented in this chapter, which is directly implemented in the time domain, can be extended to analyze the dynamic response of practical structures with nonlinearity.
Acknowledgments The authors would like to express their cordial thanks to Dr Chongmin Song of the University of New South Wales, Sydney, for providing the computer program SIMILAR and
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
138
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
for his beneficial discussions. Financial support by the National Key Basic Research and Development Program (grant no. 2002CB412709) is appreciated.
References [1] Wolf JP, Song C. Dynamic-stiffness matrix in time domain of unbounded medium by infinitesimal finite element cell method. Earthquake Eng Struct Dyn 1994;23(11): 1181e98. [2] Song C, Wolf JP. The scaled boundary finite-element method e alias consistent infinitesimal finite-element cell method e for elastodynamics. Comput Methods Appl Mech Eng 1997;147(3e4):329e55. [3] Wolf JP, Song C. Finite-element modeling of unbounded media. New York: Wiley; 1996. [4] Wolf JP, Motosaka M. Recursive evaluation of interaction force of unbounded soil in the time domain. Earthquake Eng Struct Dyn 1989;18(3):345e63. [5] Paronesso A, Wolf JP. Property matrices identification of unbounded medium from unit-impulse response functions using Legendre polynomials: formulation. Earthquake Eng Struct Dyn 1996;25(11):1231e45. [6] Paronesso A, Wolf JP. Recursive evaluation of interaction forces and property matrices from unit-impulse response functions of unbounded medium based on balancing approximation. Earthquake Eng Struct Dyn 1998;27(6):609e18. [7] Zhang X, Wegner JL, Haddow JB. Three-dimensional dynamic soilestructure interaction analysis in the time domain. Earthquake Eng Struct Dyn 1999;28(12):1501e24. [8] Chen CT. Linear system theory and design. New York: Holt, Rinehart and Winston; 1984. [9] Peterson LD. Efficient computation of the eigensystem realization algorithm. J Guid Control Dyn 1995;18(3):395e403. [10] Paronesso A, Wolf JP. Global lumped-parameter model with physical representation for unbounded medium. Earthquake Eng Struct Dyn 1995;24(5):637e54. [11] Paronesso A, Wolf JP. Property matrices identification of unbounded medium from unit-impulse response functions using Legendre polynomials: formulation. Earthquake Eng Struct Dyn 1996;25(11):1231e45. [12] Kung S. A new identification and modal reduction algorithm via singular value decompositions. In: Proceedings of the 12th Asilomar conference on circuits, system and computers. New York: IEEE; 1978. p. 705e14. [13] Juang J-N, Pappa RS. An eigensystem realization algorithm for modal parameter identification and model reduction. J Guid Control Dyn 1985;8(5):620e7. [14] Banerjee PK, Mamoon SM. A fundamental solution due to a periodic point force in the interior of an elastic half-space. Earthquake Eng Struct Dyn 1990;19(1):91e105. [15] Kobayashi S, Nishimura N. Green’s tensors for elastic half-spaces e an application of boundary integral equation method. Mem Fac Eng Kyoto Univ 1980;42(Pt 2):228e41. [16] Wong HL, Luco JE. Dynamic response of rigid foundations of arbitrary shape. Earthquake Eng Struct Dyn 1976;4(6):579e87.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6 A Coupling Procedure of Finite Element and Scaled Boundary Finite Element Methods for SoileStructure Interaction in the Time Domain1 Junyi Yan, Chuhan Zhang, Feng Jin Department of Hydraulic Engineering, Tsinghua University, Beijing
O U T L I N E 6.1 Introduction
118
6.2 Motion Equations of Coupling System
120
6.3 Realization and Model Approximation 6.3.1 Diagonalization 6.3.2 Decomposition 6.3.3 Truncation 6.3.4 Realization and Model Reduction
121 121 122 123 125
6.4 Evaluation of Interaction Forces
127
6.5 Numerical Verification 6.5.1 Accuracy 6.5.1.1 Spherical Cavity Embedded in Full Space 6.5.1.2 Dynamic Boussinesq Problem
127 127 127 129
International Journal for Numerical Methods in Engineering 2004; 59: 1453e71. Ó 2004 John Wiley & Sons Ltd. 1
Seismic Safety Evaluation of Concrete Dams http://dx.doi.org/10.1016/B978-0-12-408083-6.00006-4
117
Copyright Ó 2013 Tsinghua University Press. Published by Elsevier Inc. All rights reserved.
118
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
6.5.1.3 Massless Square Plate Resting on Half-Space 6.5.2 Efficiency
131 133
6.6 Conclusions
137
Acknowledgments
137
References
138
6.1 INTRODUCTION Dynamic analysis of soilestructure systems is an important topic in earthquake engineering. During recent decades many numerical procedures have been developed to deal with this kind of problem. Usually, the structure including a finite region of adjacent soil is modeled by the finite element (FE) method while the remaining unbounded soil is represented by boundary conditions on its interface with the structure, which leads to a model with a finite number of degrees of freedom. The boundary conditions describe the forceemotion relationship on the soilestructure interface, whose rigorous form is global in space and time. The boundary element (BE) method has become a popular procedure for constructing such a representation of boundary conditions. Another kind of procedure is categorized as various transmitting boundaries connected with finite elements to prevent a spurious wave reflection from the artificial boundaries. This procedure is usually constructed to be local in space and time to reduce the computational effort. Each of these two representations has its own merits and disadvantages. The BE method is well known as a powerful procedure for modeling the unbounded medium since the radiation condition is satisfied automatically as a part of the fundamental solution. However, the fundamental solution is not always available or it can be very complicated even if it exists. The scaled boundary finite element (SBFE) method [1,2], however, combines the advantages of the FE method and the BE method. Only the boundary is discretized and no fundamental solution is required. As a semi-analytical algorithm, it is exact in the radial direction and converges to an exact solution in the circumferential directions when the FE size approaches infinity. The SBFE method can be applied to static, diffusion and wave propagation problems in both frequency and time domain [3]. Therefore, developing a coupling procedure of FE and SBFE in the time domain to perform linear and nonlinear seismic analysis of the unbounded soilestructure system is the main objective of this chapter.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.1. INTRODUCTION
119
For numerical algorithms in the time domain such as a BE or SBFE procedure, the dynamic property of the unbounded medium is normally represented by a unit-impulse response matrix and the interaction-force vector is evaluated by the convolution integrals of the unit-impulse response matrix and the motion vector. A large computational effort (proportional to the square of the number of time stations) and storage requirement result in an expensive cost in the analysis of the practical problems with large degrees of freedom on the soilestructure interface. To reduce the computational effort, Wolf and Motosaka [4] presented two approximate techniques to evaluate the interaction forces recursively. Later, the concept of linear system theory was introduced by Paronesso and Wolf [5] for constructing rational functions in the frequency domain. The unit-impulse response matrices using expansion of a series of Legendre polynomials are obtained approximately, by which the equivalent stiffness, damping and mass matrices of the unbounded medium can be produced. Furthermore, a system realization algorithm implemented in the time domain is employed by Paronesso and Wolf [6] to calculate the interaction forces with a system of linear equations instead of convolution integrals. From physical aspects of the unbounded soil, a piecewise linear function is used by Zhang et al. [7] to describe the acceleration unitimpulse response matrix by which the implementation of the convolution integrals is simplified. However, the accuracy of this scheme is uncontrollable, which reduces the primary advantage of the SBFE method. In this chapter, a new coupling procedure of FEeSBFE is developed for three-dimensional (3D) dynamic analysis of unbounded soilestructure interaction in the time domain. The SBFE is used to represent the unbounded soil, while the bounded structure including a portion of the adjacent soil is modeled by FE. The concept of linear system theory [8] is applied to calculate the soilestructure interaction forces instead of using the time-consuming convolution integrals. To calculate the unit-impulse response matrix using the SBFE method, the convolution integral and Lyapunov function must be treated at each time station, resulting in enormous computational effort. Since a partial sequence of the unitimpulse response matrix is sufficient to represent the main character of the unbounded medium, only the terms before a cut-off time station are required and numerical cost can thus be reduced. In addition, a powerful system realization algorithm called FastERA [9] is employed to further reduce the computational effort. The accuracy of the presented procedure can be controlled by prescribed parameters. Thus, the advantage of high accuracy of the SBFE method is retained. Numerical examples of foundation response show good agreement compared with previous results and high efficiency compared with direct convolution integral implementation. The new scheme presented herein can be applied to the
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
120
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
dynamic analysis of unbounded soilestructure interaction of large-scale systems.
6.2 MOTION EQUATIONS OF COUPLING SYSTEM To demonstrate the accuracy and efficiency of the presented numerical scheme, only the unbounded soilestructure coupling system subjected to external time-varying loading is considered here. As shown in Fig. 6.1, the generalized structure including the real structure and a finite region of adjacent soil, either of which may exhibit nonlinear behavior, is modeled by FE, while the unbounded soil, which is assumed to be homogeneous and linear without losing its generality, is represented by SBFE. The equation of motion of the structure is written in the time domain as €s u Css Csb u_ s Kss Ksb us ps Mss Msb þ þ ¼ €b u Mbs Mbb Cbs Cbb Kbs Kbb ub rb u_ b (6.1) and the contribution of the unbounded soil can be represented by the interaction forces as Zt rb ðtÞ ¼
Gbb ðt sÞ€ ub ðsÞds
(6.2)
0
where M, C and K are the mass, damping and stiffness matrices of the structure, respectively; Gbb (t) is the acceleration unit-impulse response matrix of the unbounded soil which is determined via the SBFE method; u, u_ and u¨ are displacement, velocity and acceleration vectors, respectively, p is the external load vector and r is the interaction force vector.
FIGURE 6.1 Problem definition of dynamic unbounded soilestructure-interaction analysis: (a) coupled system; (b) unbounded soil.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
121
The subscripts b and s denote, respectively, the degrees of freedom associated with the soilestructure interface and the remaining ones associated with the structure. In general, the interaction forces can be expressed as a function of either the displacements or the accelerations. To avoid singularity, acceleration formulation is used herein. Using the time-integral scheme, the seismic response of the coupling system can theoretically be obtained by solving Eqs (6.1) and (6.2).
6.3 REALIZATION AND MODEL APPROXIMATION As a rigorous representation, interaction forces are evaluated by convolution integrals as shown in Eq. (6.2), which is exact but very expensive. The cost of the direct convolution integral implementation is proportional to the square of the number of time stations, and requires enormous storage capacities. Herein, concepts of linear system theory such as realization and model reduction are applied and the computational effort can be reduced enormously. To ensure that the present scheme is efficient and stable, diagonalization, decomposition and truncation are performed before the system realization. For the sake of conciseness, the subscript b of Gbb(t) is dropped in the following sections.
6.3.1 Diagonalization As the number of the degrees of freedom associated with the soilestructure interface is usually large in practical analysis, it is very inefficient and even unrealistic to apply the realization algorithm directly to the total matrix G(t). Therefore, the so-called diagonalization procedure is performed before the realization is conducted. Details of the diagonalization procedure are presented by Paronesso and Wolf [10,11]. With a time-independent transformation matrix T, the symmetric matrix G(t) of order n n is converted into the corresponding diagonal matrix Gd(t). The relationship between G(t) and Gd(t) can be expressed as GðtÞ ¼ TGd ðtÞTT
(6.3)
Premultiplying the transformation T and T , respectively, on rdb ðtÞ and u¨b(t), i.e., T
€ db ðtÞ ¼ TT u € b ðtÞ u
(6.4a)
rb ðbÞ ¼ Trdb ðtÞ
(6.4b)
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
122
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
and substituting Eqs (6.3) and (6.4) into Eq. (6.2) yields Zt rdb ðtÞ
Gd ðt sÞ€ udb ðsÞds
¼
(6.5)
0
€ db ðtÞ and rdb ðtÞ are, where u rb(t), which are defined
respectively, the transformed vector of u¨b(t) and by Eqs (6.4a) and (6.4b). The transformation process from Eq. (6.2) to Eq. (6.5) is called diagonalization. In Eq. (6.5), each term of Gd(t) can be addressed separately, which means that a multi-input multi-output (MIMO) system can be performed as a number of single-input single-output (SISO) systems. Therefore, approximation and realization are discussed in terms of the individual scalar term Gd(t) in the following sections.
6.3.2 Decomposition In view of the fact that the Hankel matrix realization algorithm is stable and robust when the corresponding Markov parameter sequence has the convergence-to-zero property [12], Gd(t) is decomposed as shown in Fig. 6.2, yielding Gd ðtÞ ¼ Cd HðtÞ þ Kd tHðtÞ þ Gm ðtÞ
(6.6)
with coefficients Kd and Cd chosen in such a way to achieve Gm(t) converging to zero as time increases to infinity, where H(t) is the Heaviside-step function. Then, Eq. (6.5) can be rewritten as Zt rdb
¼
Kd udb ðtÞ þ
Gm ðt sÞ€ udb ðsÞds
Cd u_ db ðtÞ þ
(6.7)
0
FIGURE 6.2
Decomposition of a typical term of the diagonalized unit-impulse response
matrix.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
123
where Kd, Cd and Gm (t) are diagonal matrices with the terms as Kd, Cd and Gm(t), respectively. The computation of the first two terms of the right-hand side of Eq. (6.7) is straightforward and is not discussed further. The evaluation of the third term expressed as a convolution integral can be accomplished with ease when each of scalar convolution integrals of Gm(t) and u€db ðtÞ is properly calculated. The existence of the decomposition is guaranteed by the fact that the second derivative of the acceleration unit-impulse response matrix is equal to the regular part of the displacement unit-impulse response matrix which converges to zero in a physical sense when time increases to infinity. The coefficient matrix Kd in Eq. (6.7) corresponds to the static-stiffness matrix of the unbounded medium K0. The relationship between these two matrices can be expressed with the transformation defined by Eq. (6.3) as K0 ¼ TKd TT
(6.8)
The static-stiffness matrix K0 can be obtained by solving an algebraic Riccati equation as mentioned by Wolf and Song [3]. The coefficient matrix Cd cannot be obtained directly, and is achieved approximately in Section 6.3.3.
6.3.3 Truncation Considering the convergence-to-zero property, Gm(t) is approximated as d G ðtÞ Cdc Kd t; 0 t < tc : (6.9) Gm ðtÞz 0; t tc where tc is the cut-off time and the truncated function Gm(t) is used. Coefficient Cdc can then be calculated from Eq. (6.9) with the known values of Kd and Gd(t). Whether the truncated function of the unit-impulse response matrix satisfies the prescribed tolerance depends on the cut-off time tc. For determination of the cut-off time tc the following criterion is proposed, which has a proven applicability. In the process of calculating the unit-impulse matrix via SBFE with a prescribed tolerance parameter "tol", the normalized Frobenius norm F(t)/F(0) is computed at each time station, where 11=2 0 n X n X 2 € ðtÞ A G (6.10) FðtÞ ¼ @ ij i¼1 j¼1
Once the inequality FðtÞ=Fð0Þ tol
(6.11)
is satisfied, the time station is chosen to be the cut-off time tc.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
124
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.3 Finite element meshes of soilestructure interfaces of typical excavations in half-space: (a) prismatic excavation; (b) hemispherical excavation.
FIGURE 6.4 Variation of the normalized Frobenius norms against dimensionless time.
To demonstrate the decaying behavior of the normalized norm F(t)/F(0), two typical cases of excavated foundation are analyzed. One is the half-space with a square prismatic excavation with length 2r0 and depth r0, the other is the half-space with a hemispherical excavation with radius r0. Figure 6.3 shows the FE meshes of one quarter of the soilestructure interfaces of these two cases. The normalized Frobenius norms F(t)/F(0) are plotted in Fig. 6.4 as functions of the dimensionless time t ¼ tcp =r0 ; where cp is the velocity of the dilatational wave. It is evident that the values of F(t)/F(0) converge to zero monotonously.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.3. REALIZATION AND MODEL APPROXIMATION
125
Complete description of a system can be expressed by an infinite Hankel matrix, which is the starting point of the Eigensystem realization algorithm (ERA) [13]. When truncation is performed, the infinite Hankel matrix can be written as 2 3 ~m 2 / G ~ m kc 1 0 / ~m 1 G G 6 7 6 7 ~m 3 / ~m 2 6 G 7 G 0 0 . 6 7 6 7 6 « « 1 « « « 7 7 (6.12) HNN ð0Þ ¼ 6 6 m 7 ~ kc 1 6G 7 0 / 0 0 / 6 7 6 7 6 7 0 0 / 0 0 / 4 5 « « / « « 1 ~ m ðkÞ is the Markov parameter sequence sampled from in which G m G (t) with Dt and the total number of time stations kc equals tc/Dt. As the order of the system is equal to the rank of the Hankel matrix theoretically, the Hankel matrix can be reduced dramatically as 2 6 6 6 Hðkc 1Þðkc 1Þ ð0Þ ¼ 6 6 6 4
~m 1 G ~m 2 G
~m 2 G ~m 3 G
«
«
1
0
/
~ m kc 1 G
/ /
~ m kc 1 3 G 7 7 0 7 7 7 « 5
(6.13)
0
Since the order of the Hankel matrix is reduced, the following realization process involving eigensystem decomposition of the inner product of the Hankel matrix is more feasible and economical.
6.3.4 Realization and Model Reduction Considering the third term of the right-hand side of Eq. (6.7) in which Gm(t) has been diagonalized, each SISO system can be expressed as Zt rdb ðtÞ
¼
€db ðsÞds Gm ðt sÞu
(6.14)
0
In linear system theory, such convolution integral expression, called inputeoutput description, can be transformed into an equivalent system of linear equations, named state-variable description. In a discrete-time
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
126
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
system, this state-variable description is represented by finite-difference equations of first order at distinct time station k as ~ e xe ðkÞ þ B ~ eu €db ðkÞ xe ðk þ 1Þ ¼ A
(6.15a)
~ e xe ðkÞ þ D ~ eu €db ðkÞ rbd ðkÞ ¼ C
(6.15b)
where xe (t) is called the state-variable vector, the coefficient matrix ~ e; D ~ e; C ~ e g is named the realization of the unit-impulse ~ e; B quaternion fA m ~ response matrix G ðkÞ which is sampled from Gm(t) with time step Dt; and the subscript e corresponds to the specified SISO system. The ~ m ðkÞ using Hankel ~ e; D ~ e; C ~ e g can be calculated from G ~ e; B quaternion fA matrix realization algorithms. To reduce computational effort and eliminate the effects of unavoidable errors such as numerical round-off, it is desirable to construct a reduced model of the original system. This is generally achieved by singular value decomposition because the rank of the Hankel matrix, equivalent to the order of the system, can be evaluated numerically by non-zero singular values. Contributions associated with the smallest singular values below a prescribed tolerance threshold are neglected, thus the approximate model with reduced order is realized according to the remaining singular values. Herein, a powerful realization algorithm named FastERA is employed, in which only partial eigensystem decomposition of the inner product of the Hankel matrix is computed, rather than the full singular value decomposition of the Hankel matrix as processed in ERA. The reduced order of the system is determined with a prescribed allowance parameter "tole". Define a threshold of eigenvalue lc, which satisfies , !1=2 N X li (6.16) tole ¼ lc i¼1
where N and li denote, respectively, the order and the eigenvalue of the inner product of the Hankel matrix. The sum of all the eigenvalues in Eq. (6.16) can be calculated through the trace of the above inner product of the Hankel matrix. When the threshold lc is determined by Eq. (6.16) and the upper bound of the largest eigenvalue lmax is estimated with Gerschgorin disk theorem, the partial eigensystem decomposition corresponding to the eigenvalues between lc and lmax can be accomplished by the Sturm sequence method. The reduced ~ e and D ~ e; C ~ e can be ~ e; B order is then determined and the matrices A obtained. ~ e is Since the maximum magnitude of the eigenvalues of matrix A smaller than 1, the procedure is stable.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
127
6.5. NUMERICAL VERIFICATION
6.4 EVALUATION OF INTERACTION FORCES When each element of matrix Gm bb ðtÞ is approximated and realized, Eq. (6.15) is assembled and the interaction forces shown in Eq. (6.7) are evaluated in the discrete-time system as ~ ~ ud ðkÞ xðk þ 1Þ ¼ AxðtÞ þ B€ b
(6.17a)
~ ~ ud ðkÞ € db ðkÞ þ CxðtÞ rdb ðkÞ ¼ Kd udb ðkÞ þ Cdb u þ D€ b
(6.17b)
where the state-variable vector x is assembled from xe; and the coefficient ~ and D ~ B; ~ C ~ are block diagonal matrices with respect to matrices A; ~ e and D ~ e; C ~ e , respectively. Taking Eq. (6.4) into account ~ e; B the terms A gives g
g
g
€ b ðkÞ þ yðkÞ rb ðkÞ ¼ Kbb ub ðkÞ þ Cbb u_ b ðkÞ þ Mbb u
(6.18)
in which ~ ~ T ; YðkÞ ¼ TCxðkÞ Kbb ¼ TKd TT ; Cbb ¼ TCdc TT ; Mbb ¼ TDT (6.19) g
g
g
Discretizing Eq. (6.1) with time step Dt and substituting Eq. (6.18) into Eq. (6.1), the seismic response of the coupling system can be solved using a time-integral scheme such as the Newmark method.
6.5 NUMERICAL VERIFICATION The coupling procedure proposed in this chapter can be used in the analysis of the soilestructure interaction system subjected to direct external loading as well as seismic free-field input. Only external loading is taken into account in the following examples. The seismic free-field formulation will be reported later.
6.5.1 Accuracy This section verifies the numerical accuracy of the presented procedure compared with previous results by other researchers. To demonstrate the efficiency of the presented procedure, solutions calculated by a conventional convolution method are also presented and compared. 6.5.1.1 Spherical Cavity Embedded in Full Space A spherical cavity embedded in the full space subjected to a uniform normal pressure is analyzed. As a one-dimensional problem, its analytical
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
128
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.5 Finite element mesh of one octant of a spherical shell.
solution is available [3]. The rounded triangular pressure pulse acting on the wall of cavity is represented by
8 < p0 1 cos 2pt ; 0 t t 0 2 t0 pðtÞ ¼ (6.20) : 0; t > t0 where t0 ¼ 3.46r0/cp, in which r0 is the radius of the cavity and cp denotes the velocity of the dilatational wave. Owing to the symmetry of the problem, only one octant of the cavity is considered. As shown in Fig. 6.5, the spherical shell with thickness r0, a portion of the soil adjacent to the wall, is treated as the structure and discretized by three 20-node isoparametric hexahedral elements and the soilestructure interface is discretized by three eight-node isoparametric quadrangular elements. The total number of degrees of freedom for the whole system is 78, of which 33 are on the soilestructure interface. In Fig. 6.5, the nodes on the soilestructure interface are denoted by black circles and the remaining nodes of the structure by blank circles. For comparison with previous results, the dimensionless time is defined as t ¼ tcp =r0 ; and the dimensionless radial displacement is uðtÞ ¼ u0 ðtÞKN =ð4pr20 p0 Þ; where u0 is the radial displacement of the wall and KN ¼ 16pGr0 is the static-stiffness coefficient with shear modulus G. Poisson’s ratio is selected as 0.25 here. The results of the present method with two pairs of tolerance parameters are plotted in Fig. 6.6 (in which the first and the second parameters denote tol and tole, respectively) together with those of the convolution method and the analytical solution. Excellent agreement is observed, especially between the solutions of the convolution method and the present procedure, with smaller tolerance parameters. It is worth noting that the present procedure converges to the convolution method theoretically when tolerance parameters approach infinitesimal values.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
FIGURE 6.6
129
Displacement response of a spherical cavity embedded in full space.
6.5.1.2 Dynamic Boussinesq Problem A 3D, homogeneous, isotropic, elastic half-space subjected to a concentrated normal harmonic loading acting on its surface, named the dynamic Boussinesq problem, is presented. Owing to symmetry, only one-quarter of the half-space is modeled, as shown in Fig. 6.7. Similarly, one octant of a spherical portion of the soil is modeled as the structure and discretized by seven 20-node hexahedral elements, and the structureemedium interface is discretized by three eight-node quadrangular elements similar to the one shown in Fig. 6.5.
FIGURE 6.7
Finite element mesh of one quarter of half-space.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
130
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.8 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.2): (a) real part; (b) imaginary part.
For comparison with previous solutions, the steady-state response of the vertical displacement at point A (Fig. 6.7) is calculated in the time domain and transformed into the form of frequency domain. The dimensionless vertical displacement is defined as uz ¼ uz Gr=P and the dimensionless frequency a0 ¼ ur/cs, where G is the shear modulus and cs is the velocity of the shear wave. Comparisons between the results are drawn in Figs 6.8 and 6.9 for values of Poisson’s ratio of 0.2 and 0.4, respectively. Both the parameters tol and tole are selected as 0.005. Again, good agreement is seen. As a whole, the results obtained by the present procedure and the convolution method are close to those of Kobayashi and Nishimura [15]. Once
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
131
FIGURE 6.9 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.4): (a) real part; (b) imaginary part.
again, excellent agreement between the solutions of the convolution method and the present procedure is evident. As a typical example, Fig. 6.10 demonstrates that the results calculated from the present procedure are close to those of the convolution method when smaller tolerance parameters tole are selected. Apparently, in contrast to that of low frequency, the response of high frequency is more sensitive to the parameter tole. 6.5.1.3 Massless Square Plate Resting on Half-Space As a 3D problem, the response of a massless square plate resting on half-space subjected to uniform normal harmonic pressure is analyzed. Herein, the stiffness of the plate is assumed to be zero, which
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
132
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.10 Dynamic Boussinesq problem (Poisson’s ratio ¼ 0.2): (a) real part; (b) imaginary part.
corresponds to a uniform pressure acting on the surface of the half-space directly. Similar to the previous example, one-quarter of the half-space is modeled as shown in Fig. 6.11, where the shadowed square region denotes the area of the imposed loading. A square prism portion of the soil is modeled as the structure and discretized by 32 20-node hexahedral elements with 553 degrees of freedom, and the soilestructure interface is discretized by 12 eight-node quadrangular elements with 129 degrees of freedom as shown in Fig. 6.3(a). Poisson’s ratio is selected as {1/3} here. Similarly, the steady-state response calculated in the time domain is transformed into the form of frequency domain. The vertical
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
133
FIGURE 6.11 Finite element mesh of one quarter of half-space.
FIGURE 6.12 Vertical compliances of a square foundation under a uniform stress distribution (n ¼ 1/3).
dimensionless displacement amplitude is defined as uz ¼ uz Gb=P and the dimensionless frequency a0 ¼ ub/cs with the length of the square plate 2b and the amplitude of total load P. Solutions of the presented method with tole ¼ 0.005 and tole ¼ 0.01 are illustrated in Fig. 6.12 compared with those given by Wong and Luco [16]. Good agreement is also achieved.
6.5.2 Efficiency This section compares the efficiency of the presented procedure with the direct convolution integral method. Contributions of each
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
134
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.13 Comparison of numerical cost of different method: (a) total central processing unit (CPU) time cost; (b) numerical cost used in evaluating interaction forces.
approximation strategy in numerical cost savings are also demonstrated. As a typical example, the dynamic Boussinesq problem with a Poisson’s ratio of 0.2, as mentioned in Section 6.5.1.2, is addressed herein. The relationship of the numerical cost to the total computational time step is analyzed. It is obvious that the total computational effort spent in the convolution method is proportional to the square of the total time step, while that of the present procedure increases linearly, as demonstrated in Fig. 6.13(a). The computational effort expended in evaluating the interaction forces varies in the same way, as shown in Fig. 6.13(b). When the number of total time step is large, as encountered in practical seismic analysis, the presented procedure is far more economical.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.5. NUMERICAL VERIFICATION
135
FIGURE 6.14 Numerical cost structure: (a) convolution method; (b) present procedure. SBFE: scaled boundary finite element.
Compared with the convolution method, an additional computation is needed to accomplish the realization implementation in the present procedure. As shown in Fig. 6.14(b), this part of the cost (denoted by Realization) is trivial, because the partial sequence of the Markov parameters is usually quite short, which implies that the order of the Hankel matrix is low. The benefit of the model reduction is demonstrated as follows. For the analyzed problem with 38 degrees of freedom on the soile structure interface, the unbounded medium can be represented by 741 SISO systems in a diagonalization sense (see Section 6.3.1). With the prescribed parameter tol ¼ 0.005, the order of each system is
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
136
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
FIGURE 6.15
Distribution of the single-input single-output (SISO) systems versus the
reduced order.
determined to be 41. When the model reduction is performed, the order of the system can be reduced dramatically. Figure 6.15 shows the distribution of the SISO systems versus the reduced order with different values of tole. Effects of the parameter tole on the computational cost are illustrated in Fig. 6.16, where the numerical cost used in realization implementation is denoted by Realization and that spent in solving the final coupling equations is termed Solving. Note that when tole is selected as 0.005, the numerical cost for these two aspects can be
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
6.6. CONCLUSIONS
137
FIGURE 6.16 Effects of the parameter "tole" on computational effort. CPU: central processing unit.
cut down by more than 80% compared with that without model reduction.
6.6 CONCLUSIONS A coupling procedure of FE and SBFE methods is developed in the time domain for 3D dynamic analysis of unbounded soilestructure interaction. The dynamic response of the structure and adjacent soil can be evaluated considering the effects of the surrounding unbounded soil. To enable the application of SBFE to practical engineering analysis, several techniques are used to achieve computational efficiency, such as diagonalization, sequence decomposition and truncation, system realization and model reduction. Furthermore, the balance between accuracy and efficiency can be controlled by prescribed parameters, and thus the main advantage of SBFE is retained. The accuracy and efficiency of the new scheme are demonstrated by numerical examples. The coupling procedure presented in this chapter, which is directly implemented in the time domain, can be extended to analyze the dynamic response of practical structures with nonlinearity.
Acknowledgments The authors would like to express their cordial thanks to Dr Chongmin Song of the University of New South Wales, Sydney, for providing the computer program SIMILAR and
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS
138
6. A COUPLING PROCEDURE OF FINITE ELEMENT AND SCALED BOUNDARY
for his beneficial discussions. Financial support by the National Key Basic Research and Development Program (grant no. 2002CB412709) is appreciated.
References [1] Wolf JP, Song C. Dynamic-stiffness matrix in time domain of unbounded medium by infinitesimal finite element cell method. Earthquake Eng Struct Dyn 1994;23(11): 1181e98. [2] Song C, Wolf JP. The scaled boundary finite-element method e alias consistent infinitesimal finite-element cell method e for elastodynamics. Comput Methods Appl Mech Eng 1997;147(3e4):329e55. [3] Wolf JP, Song C. Finite-element modeling of unbounded media. New York: Wiley; 1996. [4] Wolf JP, Motosaka M. Recursive evaluation of interaction force of unbounded soil in the time domain. Earthquake Eng Struct Dyn 1989;18(3):345e63. [5] Paronesso A, Wolf JP. Property matrices identification of unbounded medium from unit-impulse response functions using Legendre polynomials: formulation. Earthquake Eng Struct Dyn 1996;25(11):1231e45. [6] Paronesso A, Wolf JP. Recursive evaluation of interaction forces and property matrices from unit-impulse response functions of unbounded medium based on balancing approximation. Earthquake Eng Struct Dyn 1998;27(6):609e18. [7] Zhang X, Wegner JL, Haddow JB. Three-dimensional dynamic soilestructure interaction analysis in the time domain. Earthquake Eng Struct Dyn 1999;28(12):1501e24. [8] Chen CT. Linear system theory and design. New York: Holt, Rinehart and Winston; 1984. [9] Peterson LD. Efficient computation of the eigensystem realization algorithm. J Guid Control Dyn 1995;18(3):395e403. [10] Paronesso A, Wolf JP. Global lumped-parameter model with physical representation for unbounded medium. Earthquake Eng Struct Dyn 1995;24(5):637e54. [11] Paronesso A, Wolf JP. Property matrices identification of unbounded medium from unit-impulse response functions using Legendre polynomials: formulation. Earthquake Eng Struct Dyn 1996;25(11):1231e45. [12] Kung S. A new identification and modal reduction algorithm via singular value decompositions. In: Proceedings of the 12th Asilomar conference on circuits, system and computers. New York: IEEE; 1978. p. 705e14. [13] Juang J-N, Pappa RS. An eigensystem realization algorithm for modal parameter identification and model reduction. J Guid Control Dyn 1985;8(5):620e7. [14] Banerjee PK, Mamoon SM. A fundamental solution due to a periodic point force in the interior of an elastic half-space. Earthquake Eng Struct Dyn 1990;19(1):91e105. [15] Kobayashi S, Nishimura N. Green’s tensors for elastic half-spaces e an application of boundary integral equation method. Mem Fac Eng Kyoto Univ 1980;42(Pt 2):228e41. [16] Wong HL, Luco JE. Dynamic response of rigid foundations of arbitrary shape. Earthquake Eng Struct Dyn 1976;4(6):579e87.
II. DYNAMIC SOILeSTRUCTURE AND FLUIDeSTRUCTURE INTERACTIONS