A fast-multipole accelerated method of fundamental solutions for 2-D broadband scattering of SH waves in an infinite half space

The traditional method of fundamental solution (T-MFS) is known as an effective method for solving the scattering of elastic waves, but the T-MFS is inefficient in solving large-scale or broadband frequency problems. Therefore, in order to improve the performance in efficiency and memory requirement for treating practical complex 2-D broadband scattering problems, a new algorithm of fast multi-pole accelerated method of fundamental solution (FM-MFS) is proposed. Taking the 2-D scattering of SH waves around irregular scatterers in an elastic half-space as an example, the implementation steps are presented in detail. Based on the accuracy and efficiency verification, the FM-MFS is applied to solve the broadband frequency scattering of plane SH waves around group cavities, inclusions, a V-shaped canyon and a semi-elliptical hill. It shows that, compared with T-MFS, the FM-MFS has great advantages in reducing the consumed CPU time and memory for 2-D broadband scattering. Besides, the FM-MFS has excellent adaptability both for broad-frequency and complex-shaped scattering problems.

Among different boundary-type methods, the method of fundamental solution (MFS) not only reduces the dimension of problem and automatically satisfy the boundary conditions at infinite, but also is a meshless method.In the MFS, the approximate solution of the partial differential equation is constructed by the linear combination of a series of Green's functions.In order to avoid singularities, the fictitious wave sources are placed at some distance to the physical boundaries of scatters.The densities of the fictitious sources are firstly obtained by the boundary integral equations discretized through boundary collocation points.So, the MFS can also be regarded as a special indirect boundary integral equation method which is very similar to the Trefftz method [23].
Due to the excellent numerical accuracy and simplicity of usage, the MFS (also named indirect boundary integral equation method) has been widely applied in the field of wave scattering in an elastic medium [24][25][26][27][28][29][30][31][32].However, just as other boundary-type numerical methods, the MFS suffers a large difficulty in treating large-scale models.A significant weakness of MFS is that the calculation efficiency of solving high-DOF problem (such as large scale, broadband frequency or multi-body scattering) is seriously affected by its dense characteristic of the equation matrix.Chen et al. [33] greatly improved the computational efficiency of MFS using a sparse storage method, but the method still has some limitations to overcome this weakness.Some researchers tried to combine the FMM with the MFS in recent years.In essence, The FMM separates source points from collocation points of the kernel function expressions to change the Green's function calculation from one-to-one interaction to group-to-group interaction.Then the FMM can substantially reduce the memory requirement and improve the computational efficiency [34][35][36].
The FM-MFS has been developed for solving many large-scale practical engineering problems [37][38][39].The memory of MFS can be reduced to  () by the latest development of FMM, and the computation amount is reduced to  (log2) order.It should be noted that the references mentioned above are all about the potential problems, the static analysis or the acoustic simulation.
According to the author's knowledge, the FM-MFS applied to elastic waves scattering problem in solid and the corresponding research are still rare at present.It is worth mentioned that there have been some reports related to the FMM-BEM [40,41], but the associated BEM is largely different from the meshless MFS.
In this paper, in order to overcome the disadvantage of traditional MFS and solve large-scale elastic wave-motion problems in a half-space, a new algorithm of the FM-MFS is proposed and validated.The method is based on the single layer potential theory, and introduces the fast multi-pole expansion to construct the Green's functions matrix.The virtual source density is obtained through GMRES iteration, of which the solving process is direct, and it is convenient for numerical implementation.Taking the 2-D scattering of SH waves around an irregular scatter in an elastic half-space as an example, the implementation steps are presented in detail, and then several typical examples indicate the excellent accuracy and efficiency of the method.Finally, the scattering of plane SH waves by randomly distributed cavity or inclusion group, by a canyon and hill topography in a half-space are solved and discussed.

The traditional MFS
Fig. 1 shows the calculation model.Assuming the incident waves are plane SH waves ,  are the two parameters of Lame constant, and  is the medium density  is the shear wave number, defined as  = /, where  is the circular frequency of seismic waves,  is the velocity of shear waves.The steady-state motion equation in isotropic homogeneous medium  is as follows: where  is the anti-plane displacement.Note that the time factor exp(i) is omitted in Eq. ( 1) and thereafter.According to numerical experience, the wave sources can be placed within 0.4-0.6 times of the relevant radius of the scatterer for low frequency waves (non-dimensional frequency  ≤ 2 ,  is the scatterer radius); while for high frequency wave ( > 2), the wave sources should better be placed within 0.7-0.9times of the relevant radius of the scatterer.
Firstly, the wave field should be separated.And the total displacement field and the stress field in the elastic half-space can be expressed as: where  and  are anti-plane displacement and shear stress of the free-field under incident plane SH waves. and  are displacement and stress of the scattered wave field.The scattered field in  is constructed by the virtual sources .Consequently, the boundary condition can be written as: According to the single layer potential theory, virtual sources on  produce the scattered field in the half-space.Related displacement and stress are expressed as follows: () =  () (, ) .
In which  ∈ ,  ∈ . () is the strength of shear-wave line source corresponding to  position on virtual sources . (, ) and  (, ) are the displacement Green's functions and stress Green's functions in the elastic half-space which can be introduced as:  Set  observed points on boundary , and  points on the virtual sources .We choose the unit center physical quantity instead of unit physical quantity, then the linear system of equations become: Eq. ( 9) can be reformatted as: where  is the dense non-symmetric coefficients matrix,  is the densities of fictitious sources on virtual boundary. denotes the response related to the free field on the boundary.
To speed up the calculation, a new algorithm of fast multi-pole fundamental solution method (FM-MFS) is presented.This method can greatly shorten the calculation time and reduce the needed memory.Combining General Minimize Residual Methods (GMRES) [42] with the fast multi-pole expansion technique [43], large-scale or high-frequency scattering of elastic waves can be solved effectively.

The new fast multi-pole for MFS
In traditional method of fundamental solution, it is difficult to calculate the approximate solutions for large-scale or high-frequency problems using the GMRES iterative algorithm.As shown in Fig. 2, consider a large number of cracks and inclusions existing in the actual complex terrain.For such problems, integration equations with high DOFs need to be calculated.

Fig. 2. The calculation model of large-scale or multi-scatter problem
To solve the Eq. ( 10), the large coefficients matrix  is needed to storage.But the memory requirement is so high that this algorithm loses its practical value.The tree structure is used as the main storage and computing object in the FM-MFS, which can expand and transfer the kernel function.By taking the GMRES iterative algorithm, we can multiply  and the tree structure instead of the large coefficients matrix .Then results are obtained by controlling iterative precision.Different kernel functions apply to different expansion methods.In this paper, the plane wave expansion theorem is taken to expand the kernel function  ( ) (•) in Eq. ( 7) as: where, | ⃗| < |  ⃗| and |  ⃗| < |  ⃗|,  is the center of multi-pole expansion,  = (cos , sin ),  = 2/(2 + 1),  is the angle between   ⃗ and  axis,  = − − ,  = |  ⃗|, and the number  of truncated.Therefore, it can be shown that Eq. ( 6) can be rewritten as: In which  ( ⃗) and    ⃗ are the multi-pole moments and can be expressed as: ( ⃗) and    ⃗ in Eq. ( 15) are only related to  and  .Furthermore, it needs to calculate only once for arbitrary collocation points  for satisfying conditions.

Fig. 3. Related points for FM-MFS expansion
As shown in Fig. 3, new multi-pole moments can be obtained when points for expansion moved from  , ′ to  , ′ .The M2M transfer can be written as: From equations above, the multi-pole moments do not need to be calculated again, but the calculated coefficient should be translated into new multi-pole moments.
Assuming a point  ,which is close to point  , satisfying | ⃗| < |  ⃗|.And the local expansion to Eq. ( 15) can be expressed as: where,  ( ⃗) and    ⃗ can be obtained by M2L transformation of multi-pole moments  ( ⃗) and    ⃗ , as follows: The Finally, the local expansion transmission of L2L can be expressed as: Take the scattering of SH waves by a cavity in a half space as an example, the implementation procedures is listed as follows.
Step 1. Discretization.Discrete the boundary  and , the same as in the T-MFS, and then automatically generate the adaptive tree structure (See Fig. 4).Step 2. Upward pass.Translate information from the source points, ′ to the leaf cells.Eq. ( 16) and Eq. ( 17) are expanded in the center of leaves to generate multi-pole expansion moments.And then the multi-pole expansion point can be moved to the corresponding parent cell (M2M).Eq. ( 18) and Eq. ( 19) constitute the transfer relationship.In this step, transfer can be repeated several times to the parent cells.If satisfying the accuracy conditions, upward pass can be carried out ( ≥ 2) (see Fig. 4(a), (b)).
Step 3. Downward pass.The multi-pole expansion moments of each level cell on the tree

Non-leaf cell
Non-leaf cell structure can be transferred into corresponding local expansion moments (M2L), and Eq. ( 21) and Eq. ( 22) can be applied to build this relationship.Transfer the local expansion moments of each level cell in the tree structure to the leaf cell by using Eq. ( 23).So far, all local expansion moments at leaf cells can be translated to the boundary points  through local expansion.The total expansion and transfer processes are completely given.
Step 4. Near-source calculation.The contributions of the near field virtual source points can be calculated by the T-MFS.
Step 5. Iterations.Update the unknown vector in the system  =  until the solution converges within a given tolerance using the GMRES solver.

Accuracy and computational efficiency
Taking the wave scattering around a cavity in half-space as an example, the accuracy and numerical stability of FM-MFS are tested and compared with analytical solution [44].As is shown in Fig. 1, assuming the radius of the circular hole is , the number of collocation points on boundary  is  = 2000 ( = ), the embedment depth is  = 2, shear modular is , surface width is  = 8, the incident angle of SH waves is  = 90°,and the dimensionless frequency is  = 0-20, the expansion terms  = 10, the iterative error is  =10 -3 .Fig. 5 illustrates the displacement amplitude on the surface by FM-MFS compared with that by analytical solution.Table 1 shows the displacement amplitude on the surface center for the vertical incident SH waves ( = 90°) under different degrees of freedom.The comparison illustrates that the results by FM-MFS agree well with those by analytical solution, and have good numerical stability.To test the efficiency of FM-MFS, the scattering of plane SH wave is considered by a cylindrical cavity in an elastic half space.Model parameters are as above and the DOFs increase from 100 to 50000.As shown in Fig. 6, complexity of the T-MFS is ( ) about total CPU time, and the calculation time has reached 9147s when the DOFs are only 10000.However, the complexity of the FM-MFS is (), and the calculation time is only 30 s when the DOFs reach 50000.Near field calculation becomes unnecessary due to low frequency, sources with small radius and the tree structure with fine mesh.Fig. 7 illustrates comparison of memory between T-MFS and the present FM-MFS.The storage requirement of T-MFS increases rapidly with the increase of DOFs.When the DOFs reaches 10000, the memory is more than 1.5 GB.However, the required memory will be greatly reduced when using the FM-MFS.For 50000 DOFs, the occupied memory is only 66 MB for the FM-MFS.Therefore, compared with the T-MFS, both the CPU time and memory performance are greatly improved by FM-MFS.These advantages speed up fast solution of 2D large-scale or broadband scattering problem of elastic waves.

The scattering of SH wave by cavities group
The FM-MFS can be used to calculate the scattering of seismic waves for shallow underground The storage capacity(Mb) cavities distributed randomly in a half space.30 m under the surface of the ground, 100 random cavities are established to simulate fracture within the scope of 60 m×120 m.Radiuses of circular cavities are ranging from 1.5 m to 3 m and 10,00 boundary points are set on the surface of each cavity.The total number of DOFs reaches 10 5 .The dimensionless frequency is set as  = 150-300.The iterative error is  = 10 -3 , and  = 0°, 30°, 60°,90°.Fig. 8 shows the surface displacement amplitude of cavities.The consumed CPU time is 803 s and the storage of main matrix is 2011 MB.Numerical results show that random cavities have significant effect on the scattering of SH waves.As shown in Fig. 8, for vertical incidence of the SH waves, the displacement amplitudes of the surface is greatly decreased, for example, the displacement amplitude near  = 0 is only about 0.1.For the oblique incidence of SH waves, there is obvious amplification effect on the incident side and the amplitude close to two times of free-field displacement, but the displacement on the side behind the wave is significantly decreased.Therefore, random cavities group has a strong "barrier effect" on incident wave.Fig. 9 shows displacement amplitudes around the cavities under SH wave incidence.It can be seen that the coherence effect of scattering waves among cavities is very strong.The displacement amplitude reached more than 4.0 when the incident wave is close to the horizontal incidence.For the vertical incidence, the displacement amplitude around cavities is more than 1.5 times of the free field displacement amplitude.Therefore, for the safety of tunnel and underground pipeline engineering, crossing the cavities zone should be avoided, especially in the side of the incident wave.

The scattering of SH wave by inclusions group
The scattering of plane SH wave is considered by inclusions group in a 30×40 domain, where  is radius of the inclusion.Fig. 10 illustrates the results of displacement amplitude with variation of incident wave frequency.Calculation parameters are set to be: density ratio between the inclusion and the half space  / = 2/3, Poisson's ratio  = 0.25,  = 0.30.DOFs = 61000.The tolerance for convergence is set to  = 10 -3 .In this calculation, the CPU time is about 205 seconds and the memory usage is about 1.3 GB,  = 0.5, 1.0 , 2.0 , 5.0 .
Assume incident plane waves propagating along -direction.Fig. 10 plots the contours of displacement amplitude around inclusions group under SH wave incidence for  = 0.5, 1.0 , 2.0 , 5.0 , and shear-wave speed ratio  / = 1/2.It demonstrates that the scattering effect of elastic wave by inclusions seems more obvious as the frequency increases.For  = 1.0 , the normalized displacement amplitude reaches up to 5.3, more than five times that of incident wave, which is mainly due to coherent effects of multiple scattering wave among inclusions group.

The scattering of SH wave by a V-shaped canyon
Displacement amplitude on the surface of a V-shaped canyon under incident high-frequency SH waves are shown in Fig. 12 and Fig. 13.The discrete number of canyon surface is  = 200359, the calculation time is about 732 s, and the consumed memory is 2.4 GB.The parameters are set to be: the wideness of the V-shaped canyon 2 km, the incident frequency 20 Hz, the shear wave velocity 400 m/s, the dimensionless frequency  = 100, / = 1.0, 2.0, and the iterative error  = 10 -6 , incident angle  = 0°, 45°, 90° (Fig. 11).Numerical results show that the left of the canyon appears significant displacement amplification effect with the incidence of high-frequency SH waves (left incidence), about twice the size of free-field displacement amplitude ratio.And the right side shows a very long "shadow zone", indicating canyon has a shielding effect under the incident waves.Under oblique incidence, the amplification effect of canyon is particularly remarkable.With vertical incidence, the

The scattering and focusing of SH wave by a semi-elliptical hill
Fig. 14 illustrates the contour of displacement amplitude around a semi-elliptical hill under incident high-frequency SH waves.The parameters are set to be: the wideness of the hill 1km, the incident frequency 5 Hz, the shear wave velocity 500 m/s, the dimensionless frequency  = 10, / = 1.0, 2.0,5.0 and the iterative error  = 10 -6 , incident angle  = 60°, 90°.
The results show that the seismic response characteristics of the hill are significantly dependent on the incident wave angle and the aspect ratio of the hill.For the high steep hill, the focal effect of the seismic wave near the top is the most significant, such as for / = 5.0, the displacement peak reaches 8.3, and the corresponding semi-circular hill displacement peak is 5.5.As for the spatial distribution feature, displacement amplification area is widely distributed on the semi-circular hill, but is mainly concentrated in the hill surface for the high steep hill.Thus, in actual violent earthquakes, the damage of high and steep slopes is usually serious.In addition, due to the coherence effect of the surface creeping wave, reflected wave and direct wave, the wave crest and wave node of the displacement response appear alternately on the surface of the hill, which is more significant for a high-steep hill.Considering the change of incident angle, the focusing position is also greatly shifted for oblique incidence, and the range of displacement amplification of a high and steep hill is enlarged.Correspondingly, the peak displacement is smaller than that of vertical incidence.The similar focusing phenomenon has also been observed by Chen et al. [45].Note that the accurate solution to the high-frequency wave scattering by a high-steep hill also illustrate excellent adaptability of the FM-MFS for treating practical problems.Mr. Zhongxian Liu formulated the overall research program and guided the preparation of code.Details as follows: Zhikun Wang was responsible for writing essays and programming code, data sorting and analysis was finished by Linping Guo and the code was prepared by Fengjiao Wu, additionally, Dong Wang is in charge of formula derivation and data management.

Conclusions
Combined with the new fast multi-pole expansion technique, the large-scale or high-frequency scattering problem of SH waves by 2-D superficial irregularities or heterogeneity in a solid half-space is solved by the proposed FM-MFS.By using the plane wave expansion of SH wave potential function, requirement of the CPU time and memory of multi-pole and local expansion is substantially reduced.The accuracy of this method is verified by comparing with exact solution.It is demonstrated that the FM-MFS can effectively solve practical complex scattering problem.Numerical results also show several interesting scattering characteristics of the high-frequency elastic waves generated by a V-shaped canyon, a semi-elliptical hill and inclusions (cavities) group.All these numerical examples indicate that the proposed FM-MFS has favorable numerical accuracy and efficiency, providing a novel effective method for the simulation of large-scale and where  ( ) (•) is the Hankel function of the second kind with order 0,  and  are the unit normal vectors,  = | − |, ′ = | − ′|.As is shown in Fig. 2, ′ is the mirror point of .

4 .
a) Construction of quad-tree structure b) Virtual wave source point c) Boundary point Fig.The adaptive tree structure

Fig. 5 .
Fig. 5. Comparisons between the displacement amplitude on the surface by FM-MFS and that of analytical solution

Fig. 11 .
Fig. 11.Diffraction of plane SH waves by a V-shaped canyon

Table 1 .
Numerical stability of displacement amplitudes