A UNIFIED HIGH-ORDER APPROACH TO WAVE PROPAGATION IN BOUNDED AND UNBOUNDED DOMAINS USING THE SCALED BOUNDARY FINITE ELEMENT METHOD

Abstract. A uniform high-order time-domain approach for wave propagation in bounded and unbounded domains is proposed. It is based on improved continued-fraction expansions of the dynamic stiffness. The coefficient matrices of the continued-fraction expansion are determined recursively from the scaled boundary finite element equations in dynamic stiffness. The resulting solution is suitable for systems with many degrees of freedom as it converges over the whole frequency range, even for high orders of expansion. In the time-domain, the continued-fraction solutions correspond to equations of motion with symmetric, banded and frequency-independent coefficient matrices.


INTRODUCTION
The numerical modelling of wave propagation in unbounded domains is required in a number of engineering applications.Consider for example the design of machine foundations, earthquake analysis of large structures and dam-reservoir systems or the vibro-acoustic optimization of passenger vehicles.In these examples, displacement or pressure waves propagate in semi-infinite media such as soil, air or water.Here, it is of crucial importance to accurately model the effect of radiation damping in a numerical simulation.The well-established finite element method cannot be used straightforwardly in these types of problems, since outgoing waves are reflected at the artificial boundaries of the finite element mesh.The consistent modelling of wave propagation in unbounded domains has been a major research topic for more than 30 years, and is still challenging [7,14].The numerical approaches which have been developed in this context include absorbing boundaries [8,10,15], the boundary element method [3,4], infinite elements [5], perfectly matched layers [2] and the scaled boundary finite element method [16].
Wave propagation in bounded domains is of interest for example in non-destructive testing using wave based methods or in dynamic fracture analysis, when simulating propagating cracks.Here, the finite element method is not competitive if singularities occur inside the bounded domain due to cracks or material interfaces.Typically, very fine finite element meshes are required around crack tips.In dynamic crack propagation analyses or in the context of inverse analyses the numerical effort associated with remeshing can be prohibitive.
In this paper, a high-order mass or damping matrix approach for wave propagation analyses in bounded or unbounded domains is presented, which overcomes these drawbacks of the finite element method.It is based on a unified time-domain formulation of the scaled boundary finite element method.This is a semi-analytical technique which excels in modelling time-dependent problems in unbounded domains and in modelling bounded domains with singularities.
The scaled boundary finite element method is based on a coordinate transformation which allows the governing equations to be discretized in the circumferential directions, while the solution in the direction of wave propagation is obtained analytically.The method was originally formulated in the frequency domain, leading to a nonlinear differential equation in dynamic stiffness.The scaled boundary finite element method in the time-domain is obtained by expanding the dynamic stiffness into a series of continued fractions [1,12].An improved continued-fraction solution is presented in this paper, which converges over the whole frequency range with increasing order of expansion.Compared to an existing approach, it leads to numerically more robust solutions for large-scale systems and arbitrarily high orders of expansion.By using the continued-fraction solution and introducing auxiliary variables, the equation of motion of a bounded domain is expressed in high-order static stiffness and mass matrices.For an unbounded domain, a formulation in terms of high-order static stiffness and damping matrices is obtained.Both formulations correspond to equations of motion with symmetric, banded and frequency-independent coefficient matrices, which can be coupled seamlessly with finite elements.Standard procedures in structural dynamics are then directly applicable.
The further outline of this paper is as follows.In Section 2, the concept of the scaled boundary finite element method is briefly summarized.In Section 3, the improved continued fraction solutions of the scaled boundary finite element equations for the dynamic stiffness matrix of a bounded and unbounded domain, respectively are derived.In Section 4, the corresponding equations of motion are constructed by using the continued fraction solutions and introducing auxiliary variables.Numerical examples are analysed in Section 5.

SUMMARY OF THE SCALED BOUNDARY FINITE ELEMENT METHOD
The scaled boundary finite element method is described in detail in the book [16] and in Reference [11].For completeness, the equations necessary for the development of the high-order time-domain formulations are summarized briefly in the following.
In the scaled boundary finite element method, a so-called scaling centre O is chosen in a zone from which the total boundary, other than the straight surfaces passing through the scaling centre, must be visible (Figure 1 where ξ, η and ζ are called the scaled boundary coordinates. The nodal unknown displacements {u(ξ)} are introduced along the radial lines passing through the scaling centre O and a node on the boundary.(The dependency on the excitation frequency ω in a frequency-domain analysis or on time t is omitted from the argument for simplicity when it is not explicitly required.)The unknown displacements at a point (ξ, η, ζ) are interpolated from the nodal displacements {u(ξ)} as In Equation ( 2), the size of the identity matrix [I] is n × n, where n is the number of degrees of freedom per node.In a next step, Galerkins weighted residual technique or the virtual work method is applied in the circumferential directions η, ζ to the governing differential equations.
In the frequency domain, the scaled boundary finite element equation in displacements {u(ξ)} results, where s (=2 or 3) denotes the spatial dimension of the domain.
] and [M 0 ] are coefficient matrices obtained by assembling the element coefficient matrices as in the finite element method.For three-dimensional elastodynamic problems the element coefficient matrices are expressed as where [B In an elastodynamic problem, the internal nodal forces {q(ξ)} on a surface with a constant ξ are obtained by integrating the surface traction over elements.This yields The internal nodal forces are related to the nodal forces {R} on the boundary by {R} = −{q(ξ = 1)} for an unbounded domain.The dynamic-stiffness matrix [S ∞ (ω)] of an unbounded domain is defined by In Equation ( 6), the notations {R(ω)} = {R(ξ = 1, ω)} and {u(ω)} = {u(ξ = 1, ω)} are introduced to denote the force and displacement amplitudes at the boundary.Using Equations ( 5) and (6), the relationship between the nodal displacements and the radial derivatives of the displacements on the boundary is expressed as Using Equation (7), the scaled boundary finite element equation (3) in displacement {u(ξ)} can be transformed into the so-called scaled boundary finite element equation in dynamic stiffness (8) (see [11,16]).
Equation ( 8) is valid for an unbounded domain with dynamic stiffness [S ∞ (ω)].For a bounded domain, the internal nodal forces are related to the nodal forces {R} on the boundary by {R} = +{q(ξ = 1)}.The corresponding scaled boundary finite element equation in dynamic stiffness (9) is derived analogously to the unbounded case, where [S b (ω)] is the dynamic stiffness matrix of a bounded domain.Equation ( 8)/(9) is a system of non-linear differential equations in the independent variable ω.For ω → ∞, it can be solved using an asymptotic power expansion [16].In the rigorous scaled boundary finite element method, the dynamic stiffness matrix at intermediate and low frequency is obtained by numerical integration of Equation ( 8)/(9).This computationally expensive task is avoided by constructing a continued-fraction solution of the scaled boundary finite element equation in dynamic stiffness.This is described in detail in the following section.

CONTINUED FRACTION SOLUTIONS OF DYNAMIC STIFFNESS MATRIX
In this section, continued-fraction solutions for the dynamic stiffness matrix are determined from the SBFE equation in dynamic stiffness.Similar approaches have been originally derived in References [1] and [12] for unbounded and bounded domains, respectively.These methods, however, have only been used for the analysis of small problems.For systems with many DOFs and high-orders of expansion, the numerical steps involved in these approaches may become ill-conditioned.In Reference [6], the reason of these numerical problems has been identified studying a simple analytical example and an improved algorithm for unbounded domains has been proposed.This algorithm is presented in Section 3.1 and extended to the bounded case in Section 3.2.

Unbounded domain
In a first step, the continued-fraction solution (10) is assumed, The first two terms are the constant dashpot and spring matrix, respectively.The term R (1) (ω) denotes the yet unknown residual of the two-term expansion at high frequency.Substitution of Equation (10) in Equation ( 8) leads to The terms in Equation ( 11) can be sorted in descending order of powers of (iω).Equation ( 11) is satisfied when the two terms corresponding to (iω) 2 and (iω) and the remaining lower-order term are equal to zero.Setting the terms corresponding to (iω) 2 and (iω) 1 equal to zero yields, (iω) In the solution process, the eigenvalue problem ( 14) is used.
The eigenvalues ⌈Λ 2 ⌋ are positive, since both [M 0 ] and [E 0 ] are positive definite.Normalizing the eigenvectors [Φ] with respect to the matrix Using Equation ( 17), pre-and post-multiplying Equation ( 12) by [Φ] T and [Φ], respectively, and introducing yields The matrix [k ∞ ] results from: Equation ( 20) can be solved directly by back substitution as the coefficient matrix at the lefthand side is diagonal.The remaining part of Equation ( 11) is an equation for The unknown residual [R (1) (ω)] is expressed as In Equation ( 23), the terms [Y 1 ] are constants corresponding to the constant and linear term of the i-th continued fraction and [R (i+1) ] is the residual of the order i expansion.
[X (i) ] is a yet undetermined factor.If the factor [X (i) ] is chosen as [X (i) ] = [I], then Equations ( 22) and ( 23) are identical to the decomposition used in Reference [1].In this paper, [X (i) ] will be selected such that the robustness of the numerical algorithm is improved.
Here and in the following, the superscript −T denotes the transpose of the inverse of a matrix.Using Equations ( 18) and (19), Equation ( 25) is written as the case i = 1 of the following equation: with [b [b [c (1) ] = [X (1) Using Equation (23), Equation ( 26) is again expanded to (iω) 2 , (iω) and remaining lowerorder terms, As for Equation (11), Equation ( 28) is satisfied when all the three terms in the power series are equal to zero.Setting the (iω) 2 term to zero leads to an equation for [Y Pre-and post-multiplying the above equation with [Y Equating the terms corresponding to (iω) 1 to zero yields . The remaining lower-order term is reformulated using Equation (22).Pre-and post-multiplying the resulting equation by Equation ( 32) is the (i + 1)-case of Equation ( 26), [b Equation ( 33) can be solved by following the same steps as for solving Equation (26).For given coefficients [X (i) ], the coefficient matrices [a (i) ], [b 1 ] and [c (i) ] are evaluated recursively using Equation (34), starting from those at i = 1.An order M continued fraction terminates with the approximation [R (i+1) (ω)] = 0. Increasing the order of continued fraction does not require the recalculation of the coefficient matrices determined previously for a lower order.
The coefficients [X (i) ] are yet undetermined.It is worth noting that the algorithm presented above reduces to the method presented in Reference [1] if these coefficients are equal to [X (i) ] = [I].In the following, these coefficients are chosen such that the robustness of the approach is improved.The analytical study of the scalar wave equation in a fullspace bounded by a spherical cavity in Reference [6] has shown that the original continuedfraction approach breaks down if the corresponding scalar coefficient c (i) becomes zero.The same study has revealed that this problem can be overcome by choosing the corresponding scalar coefficient X (i) such that |c (i) | = 1.In the multidimensional case, the coefficient [c (i) ] approaches a singular matrix if the original continued-fraction procedure of Reference [1] is used.The idea of choosing the coefficient X (i) such that |c (i) | is equal to one is extended to the matrix case in the following equations.In each step of the recursive procedure, the coefficient [c (i) ] can be expressed as follows: with The matrix [c (i) ] is symmetric.It can be expressed as the product of a lower triangular matrix [L (i) ], a diagonal matrix [D (i) ] and an upper triangular matrix [L using the so-called LDL T -decomposition, see [9], Sec.5.1, page 82.Here, [L (i) ] is normalized such that the entries of the diagonal matrix [D (i) ] are ±1.0.The LDL T -decomposition of a matrix [A] is a generalization of the Cholesky decomposition, which is applicable even if the matrix Thus, the coefficient [c (i) ] is diagonal with entries c (i) kk = ±1.0 and thus perfectly wellconditioned.

Bounded domain
In order to facilitate the derivation of the continued-fraction solution of the dynamic stiffness of a bounded domain, Equation ( 9) is rewritten as where the independent variable has been changed to The derivation is started by assuming In Equation (43), the matrices [K] and [M ] are the static stiffness and mass matrix of the bounded domain, respectively.The term [R (1) ] denotes the yet unknown residual of the continued-fraction expansion, representing the high-frequency response.Substituting Equation (43) in Equation (41) yields The terms of Equation ( 44) can be written in ascending orders of powers of x.Setting the constant term equal to zero yields: Equation ( 45) is the scaled boundary finite element equation in static stiffness [K] of a bounded domain.The solution of this algebraic Riccati equation is described in detail in Reference [13] and not repeated here.Equating the linear terms in x to zero leads to: Equation ( 46) is a Lyapunov equation for the mass matrix [M ] .It solution is described in detail in Reference [11].Setting the remaining terms equal to zero yields an Equation for the unknown residual [R (1) (x)]: The unknown residual term [R (1) (x)] is decomposed as The derivative [R (i) (x)] ,x is determined as Substituting Equations ( 48) and (50) in Equation ( 47) and pre-and post-multiplying the resulting expression by [S (1) (x)][X (1) ] −1 and [X (1) ] −T [S (1) (x)], respectively, leads to an equation for [S (1) ], [S (1) Equation ( 51) is written as the case i = 1 of the following equation [b [b [c
Equation ( 65) is a standard equation of motion of a linear system in structural dynamics written in the frequency-domain.It is expressed in the time-domain as a system of first-order differential equations with high-order stiffness and damping matrices.
For a bounded domain, Equations (60b) and (64b) can be combined into the following matrix form: with 0 ], [S and {Z(ω)}, {F (ω)} given in Equation (66c).Equation (68) can be expressed in the timedomain as a standard equation of motion with high-order stiffness and mass matrices.
The matrices h are symmetric and sparse.The natural frequencies and vibration modes of a bounded domain can be determined from the eigenproblem corresponding to Equation (70).

NUMERICAL EXAMPLES
In this section, the accuracy of the proposed improved high-order formulations in both frequency and time domains is evaluated by numerical examples.Its superiority with respect to the original continued-fraction approach [1,12] is demonstrated.

Three-dimensional elastc foundation embedded in homogeneous isotropic halfspace
As a 3D vector-valued problem, vertical motion of a square foundation 2b × 2b embedded with depth e = 2/3b in a homogeneous isotropic halfspace is analysed.The system is One quarter of a square foundation embedded in halfspace; geometry and mesh shown in Figure 2.Only one quarter of the symmetric system is modelled using the SBFEM.The foundation-soil interface is meshed with 12 8-node SBFEs, leading to a total of 129 DOFs.The continued-fraction solution for the dynamic stiffness matrix [S ∞ ] is constructed using the proposed method or the original approach presented in Reference [1].The accuracy in the frequency domain is evaluated in Figures 3-5.As an example, the diagonal term S 1,1 and the off-diagonal term S 1,2 of the 129 × 129 dynamic stiffness matrix [S ∞ ] are shown in a dimensionless form.The continued-fraction solutions of order M = 3, M = 7 and M = 10 and 17 are compared to the dynamic stiffness obtained by numerical integration of Equation ( 8) in Figures 3, 4 and 5, respectively.For M = 3, the two continued-fraction solutions obtained using the original approach [1] or the proposed method are identical, as is shown in Figure 3.This low-order continued-fraction solution, however, differs strongly from the reference solution in the low-frequency range.The agreement between the dynamic stiffness obtained by numerical integration and the continued-fraction solution is improved significantly, if the proposed method is used to calculate the coefficients [Y 1 ] of an order M = 7 approximation, as can be seen in Figure 4. On the contrary, the continued-fraction solution of order M = 7 calculated using the original approach [1] is clearly erroneous.Since the continued-fraction solution of [1] diverges for M ≥ 5, it is not shown for higher orders of M .Figure 5 confirms that the proposed continued-fraction expansion converges to the exact solution for increasing order M .
The transient response of the three-dimensional soil halfspace with excavation (no foundation) initially at rest is evaluated.In the time-domain, the unbounded domain is described by the system of first-order differential equations (67).It is assumed that a vertical force P (t) acts at the bottom centre of the foundation (x = 0, y = 0, z = e).The timedependence of the excitation is prescribed as a Ricker wavelet.The time history of the Ricker wavelet is given as where t s is the time when the wavelet reaches its maximum, 2/t 0 is the dominant angular frequency of the wavelet and P 0 is the amplitude.A Ricker wavelet with the parameters  Figure 6 shows the computed vertical displacement at the bottom centre of the foundation.The results are non-dimensionalized with P 0 /Gb.The displacement response obtained using the proposed continued-fraction solution of order M = 10 and of order M = 17 is compared to the numerical result calculated using the rigorous SBFEM based on convolution.In the rigorous analysis, the time-step ∆t = 0.005b/c s is selected.The agreement between the result of the rigorous analysis based on convolution and that of the present transmitting boundary is excellent for t < 1.5b/c s .After that, very small differences occur.These deviations are due to the fact that the proposed singly-asymptotic transmitting boundary approximates the high-frequency behaviour with higher accuracy than the static stiffness.The accuracy in the time-domain can be further improved by further increasing the order of continued-fraction expansion M .

Regular polygon
As a bounded domain example, a regular octagon is considered.The scaling center is chosen at the geometry center.The elastic medium has the following parameters: modulus of elasticity E = 1 N m 2 , Poisson's ratio ν = 0.25, mass density ρ = 1 kg m 3 .Plane stress state is considered.Two nine-node elements on each edge are used.The system and scaled boundary finite element mesh are shown in Figure 7.The eigenfrequencies of the octagon are calculated solving the eigenvalue problem corresponding to Equation (70).A converged finite element solution is used as a reference solution.Figures 8 and 9 show the first 150 natural frequencies obtained with M = 6 and M = 19 using the original approach [12] and the proposed improved method, respectively.A parameter study in Reference [12] has shown that a continued-fraction expansion of order M = 6 leads to accurate results.Figure 8(b) shows that the error corresponding to M = 6 for higher modes is lower than 3%.Nevertheless, it is expected that this error decreases further for increasing order of approximation. Figure 8(b), however, shows that this is not the case if M is increased to M = 19.In fact, the error corresponding to M = 19 is considerable higher than that associated with M = 6, if the original approach [12] is used.For high degrees of approximation ill-conditioning is observed.The improved continued-fraction expansion proposed in this paper leads to more robust results, as can be seen in Figure 9. Here, the higher-order expansion clearly leads to a smaller error and strict convergence.

Transient analysis of two-dimensional bounded domain under strip loading
The in-plane motion of the two-dimensional domain shown in Figure 10 with shear modulus G, mass density ρ and Poisson's ratio ν = 0.25 is considered.A uniformly distributed strip loading P (t) is applied on the free surface.Its time-dependence is considered as a triangular impulse with an amplitude P and a duration 3b/cs, see Figure 11  Each edge of the subdomains is discretized with one ten-node high-order element.The total calculation time is 20b/c s .The fixed time step t = 0.02b/c s is selected.
The vertical displacement responses at points A, B are plotted in Figure 12.The results are non-dimensionalised with P/Gb.The agreement between the results of the finite el- ement method and those of the improved method is excellent when the order of the continuedfraction expansion M is increased to 9. On the other hand, the results of the original method calculated using M = 9 at points A, B are NaNs (Not a Number), due to ill-conditioning, and thus not plotted in Figure 12.

CONCLUSIONS
High-order time-domain formulations for the modelling of wave propagation in unbounded and bounded domains of arbitrary geometry are developed.These formulations are based on continued-fraction expansions of the dynamic stiffness.Improved continuedfraction solutions are proposed, which are computationally more robust than previous procedures [1,12].These are characterized by an additional, matrix-valued factor, which is chosen such that singularities are removed.Numerical examples demonstrate that the improved continued-fraction expansion converges to the exact dynamic stiffness for increasing orders of expansion M .No ill-conditioning is observed, even for very high orders of approximation.

Figure 1 .
Figure 1.Concept of the scaled boundary finite element method.(a) unbounded domain, (b) three-node line element on boundary of 2D problem, (c) eight-node surface element on boundary of 3D problem

Figure 3 .
Figure 3. Continued-fraction solution of order M = 3 for dynamic stiffness matrix of embedded square foundation (diagonal term S 1,1 and off-diagonal term S 1,2 )

Dimensionless
(a).The corresponding Fourier transform is plotted in Figure 11(b).In the scaled boundary finite Frequency ωb/c s Force Amplitude P(ω)/P (b) frequency domain

Figure 12 .
Figure 12.Vertical displacement due to strip loading 1 (η, ζ)] and [B 2 (η, ζ)] represent the strain-displacement relationship, [D] and |J(η, ζ)| are the elasticity matrix and the determinant of the Jacobian matrix on the boundary, respectively.The matrices [B 1 ] and [B 2 ] depend on the geometry of the boundary only.The coefficient matrices [E 0 ] and [M 0 ] are positive definite.[E 2 ] is symmetric.