Wavelet Element Modelling for Inviscid Fluid–Solid Coupling Problem based on Partitioned Approach

To provide a simple numerical formulation based on fixed grids, a wavelet element method for fluid–solid modelling is introduced in this work. Compared with the classical wavelet finite element method, the presented method can potentially handle more complex shapes. Considering the differences between the solid and fluid regions, a damping-like interface based on wavelet elements is designed, in order to ensure consistency between the two parts. The inner regions are constructed with the same wavelet function in space. In the time and spatial domains, a partitioned approach based on Jacobi iteration is combined with the pseudo-parallel calculation method. Numerical convergence analyses show that the method can serve as an alternative choice for fluid–solid coupling modelling.


Introduction
In the context of structural health monitoring, as well as seismic and acoustic exploration, the numerical modelling and analyses of fluid-solid interfaces and coupling have been considered crucial issues, due to the complexity of the related physical phenomena. For waves propagating in a single medium (i.e., fluid or solid), nearly all types of numerical methods can be used, such as the finite difference method [1,2] and its derivative formulations [3], boundary element methods [4], spectral element methods [5,6], pseudo-spectral element methods (i.e., time-domain spectral element methods) [7], wavelet spectral methods [8][9][10], and mixed formulations based on spectral methods [11]. Although the performances of these numerical methods typically differ, due to the corresponding superiorities or shortcomings, reasonable results can be usually obtained for the self-contained establishment of modelling and solving methodologies. Fluid-solid interfaces lead to significant numerical difficulties, as the correct implementation of physical matching conditions plays an important role in the convergence of numerical algorithms. Modelling such an interface involves many issues concerning the selection of grid techniques, simulation techniques, and time discretization.
In fluid-solid interaction problems, both the fluid domain and the solid domain are deformed. For elastic and incompressible problems, the solid region is usually modelled by a fixed grid technique, due to its small deformation, while the grid of the fluid part can be altered between moving and fixed grids [12]. The arbitrary Lagrangian-Eulerian formulation [13] is a typical moving grid formulation, in which the grid can move with arbitrary velocity; but not necessarily at the velocity of the fluid. The movement of the grid involves interpolation from the old to the new grid, inevitably causing errors. Moreover, the interaction between the fixed solid grid and the moving fluid grid is also complicated, although the arbitrary Lagrangian-Eulerian formulation has the advantage that the wall shear stress at the fluid-structure interface can be calculated accurately. To obtain a simpler formulation for the fluid-solid coupling problem, a fixed grid is used for the fluid. The immersed boundary method where vector u denotes the displacement of the solid. For the two-dimensional case, u is expressed by [u x , u y ] T . The symbols σ and ε denote the stress and strain tensors, respectively. The symbols λ and µ are the Lame parameters, and ρ s is the mass density of the solid waveguide, referring to the inertial force of the inspected region. The matrix I is the identity tensor. The vector f s defines the general external force in the solid region. The function tr(·) is employed to obtain the trace of the bracketed term. The fluid region, where only the acoustic wave propagates, is governed by the conservation and dynamics equations: .
where the vector v denotes the velocity of the fluid, the constant c = k/ρ f is the velocity of the acoustic wave, and k is the bulk modulus of fluid. The mass density of fluid is defined by ρ f . As the fluid is irrotational, the variable v is replaced by the gradient of the potential function ϕ. Substituting the relationship v = ∇ϕ into Equations (4) and (5), a second order system which only involves the potential function variable can be obtained: The above equations independently define the dynamic properties in the solid and fluid. To couple the two waveguides together, we need to find the connection between u and v. To address this problem, Komatitsch et al. [18] proposed a fluid-solid interface method. From Equation (4), it can be obtained as follows: Considering Equation (7) and the continuity of traction on the interface (the boundary between solid and fluid), τ = σ · n (where n denotes the unit normal to the interface), we obtain the following: The continuity of the normal component of speed is expressed as: Note that v = ∇ϕ. Equation (9) constructs the energy transport path between solid and fluid. In the present model, it should be noticed that the shear force is considered in the solid region but absent in the fluid region. The present interface fully ignores the viscosity on the acoustic wave propagation, which is not an accurate model as the equilibrium condition on the interface is not satisfied due to neglection of shear stress. The balance of tangential element of the stress on the interface is ignored. A more precise interface containing the shear effect should be considered to achieve more accurate results in practice.

The Weak Form (Numerical Models)
In finite element methods, problems are solved using their corresponding weak forms. By dotting the partial differential equation with the trial function (or interpolating function) ψ = [ψ x (x, y), ψ y (x, Materials 2020, 13, 3699 4 of 21 y)] T for the solid region, and ψ(x, y) for the fluid region, and integrating by parts over the field of interest, we obtain Equation (10) for the solid region: Ω s ψ f s dΩ, (10) and Equation (11) for the fluid region: where the inspected solid and fluid regions are denoted by Ω s and Ω f , respectively, and the interface between them is defined by the symbol Γ. In the above equations, we consider a wave source, f s , in the solid region, which can be implemented into the fluid region. Based on the interface components in Equations (10) and (11), the solid and the fluid system are coupled. In numerical formulation, the Dirichlet-Neumann decomposition method is used to connect the coupled regions: by which the velocities in the solid part can be converted to the velocity of the fluid. By . u = n x ∇ϕ, n y ∇ϕ the velocity vector in the solid can be calculated from the variables in the fluid part.

Wavelet Element Spatial Discretization
Numerical solution of the wave propagation in solid-fluid coupling problems involves the discretization of spatial and temporal domains. We chose the mth-order j scale B-spline wavelets on the interval (BSWI) interpolating function ψ j m,k (ξ) for spatial discretization. The region Ω, either for the solid or fluid, is divided into a set of non-overlapping sub-domains, Ω e , following which each sub-domain is mapped onto a unit interval, considering the dimension of the problem being analyzed. According to the mth-order 0 scale B-spline functions and the corresponding wavelets given by Goswami [20], the j scale mth-order BSWI scaling functions ψ j m,k (ξ), denoted by BSWIm j , can be derived as:  (14), the interpolating functions, in horizontal and vertical directions, are defined as: where ξ and η are restricted to the interval [0, 1], respectively, depicting the normalized x and y co-ordinates. The two-dimensional interpolating function is formulated based on the Kronecker product between the two vectors in Equation (15); namely, Ψ = ψ ξ ⊗ ψ η . In the classical finite element method, the unknown field function u is approximated by the interpolating function N and node vectors u e as: Materials 2020, 13, 3699

of 21
Note that the wavelet coefficients can be obtained when the interpolating function N is replaced by Ψ. The physical meaning of the wavelet coefficients is not evident, due to the overlap in support, and the boundary condition is not easy to implement. Therefore, an additional matrix T is required to transform the wavelet coefficients into the physical domain, in which the interpolating function N yields: where the transform matrix is Equations (10) and (11) can be rewritten as: where the superscript e depicts "elemental". Equilibrium is also satisfied in each element. Based on the Hamilton principle, we further obtain the wavelet element method in matrix form: Assembling the elementary matrices together, we further obtain: where M is the mass matrix and K is the stiffness matrix. The pseudo-damping C is the coupling matrix of the interface, which acts on the first order temporal derivatives of u and ϕ. It differs from the real damping, in that the energy is not dissipated in the layer but transported to the coupled region instead. The interface is integrated into the finite element formulation as a pseudo-damping layer, which couples the solid and fluid parts. The explicit forms of the corresponding global matrices for the two-dimensional case are: where J is the Jacobian matrix which transforms the normal region to a curved region, and the constants w i and w j are the weights referring to the Gaussian integral. As an isoparametric method, the shape function related to J is also defined by the same BSWI used for interpolation. Here, we give the explicit expression of the Jacobian matrix based on BSWI4 3 (the scale j = 3, and fourth order) element. The even order spline is usually used in structural analysis; thus, the fourth-order spline is selected. In order to ensure the existence of inner wavelet, scale j and order m should yield: Then scale j is fixed at 3. The BSWI4 0 is given by [20], then the BSWI4 3 function ψ ξ = { ψ 3 4,−3 (ξ) ψ 3 4,−2 (ξ) . . . ψ 3 4,7 (ξ)} can be derived from Equation (14) as: and Similarly, we can obtain The present element is constructed as the isoparametric element, thus Ψ is also used as the shape function. The direct interpolating by wavelet function induces the solution in the forms of wavelet coefficient a like: where F(x, y) depicts the unknow function to be interpolated. Thus, T matrix is employed to connect wavelet coefficient a with physical region u as: , and the ξ i , η i can be flexibly selected in the interval [0, 1]. As we mentioned, the interpolating function/shape function can be replaced by ΨT = N by aid of the T matrix. The Jacobian matrix: Here the location is denoted by d = (x, y) T to distinguish with deflection u. The main entries of J matrix are calculated by: The main components of ∂ψ ξ ∂ξ are: and ∂ ∂ξ ψ 3 4,5 ∂ψ η ∂η is obtained in the same manner, and the J matrix can be calculated according to Equation (37). The interpolating results of a quarter circle idealized by one BSWI4 3 element is presented in Figure 1. The element contains 11 × 11 inner nodes as shown in Figure 1a. Figure 1b,c present the distribution and values of error independently. The scattered points in Figure 1c illustrate the value of every arrow in Figure 1b. It is seen the interpolating error is limited to an acceptable level; with the use of more elements in modelling, the error can be further restrained. (37).
The interpolating results of a quarter circle idealized by one BSWI43 element is presented in Figure 1. The element contains 11 × 11 inner nodes as shown in Figure 1a. Figure 1b,c present the distribution and values of error independently. The scattered points in Figure 1c illustrate the value of every arrow in Figure 1b. It is seen the interpolating error is limited to an acceptable level; with the use of more elements in modelling, the error can be further restrained.

Temporal Discretization and Partitioned Approach
The central difference time integration scheme for temporal discretization is Equation (45): The subscripts related to the time t and time step Δt define the varying times in the numerical calculation. Two difficulties are implicated in solving the above equations directly: (1) the large dimension of the stiffness matrix and (2) the iteration between two mediums in different time steps.
To simulate wave propagation accurately, the grid should satisfy the rule that at least 10 nodes are required for each wavelength [21][22][23]. This requirement usually involves a large dimensional grid for high-frequency wave motion. On the other hand, a dense grid further compels the use of a tiny Δt (usually less than 1 × 10 −6 s for ultrasound wave simulation) in the central difference time integration scheme, in order to ensure convergence. In each time step, the matrix calculation and the inverse computation should be iteratively solved. Therefore, Equation (45) is difficult to solve on the level of the global matrix. We chose the pseudo-parallel method (as shown in Table 1) to address this problem.

Temporal Discretization and Partitioned Approach
The central difference time integration scheme for temporal discretization is Equation (45): The subscripts related to the time t and time step ∆t define the varying times in the numerical calculation. Two difficulties are implicated in solving the above equations directly: (1) the large dimension of the stiffness matrix and (2) the iteration between two mediums in different time steps.
To simulate wave propagation accurately, the grid should satisfy the rule that at least 10 nodes are required for each wavelength [21][22][23]. This requirement usually involves a large dimensional grid for high-frequency wave motion. On the other hand, a dense grid further compels the use of a tiny ∆t (usually less than 1 × 10 −6 s for ultrasound wave simulation) in the central difference time integration scheme, in order to ensure convergence. In each time step, the matrix calculation and the inverse computation should be iteratively solved. Therefore, Equation (45) is difficult to solve on the level of the global matrix. We chose the pseudo-parallel method (as shown in Table 1) to address this problem.
Secondly, the iteration between fluid and solid is also a key issue for partitioned approaches. Considering the last terms on the right-hand side of Equation (45): (a) in the solid region, in order to estimate the displacement field ut+∆t , the distribution of the potential function ϕt+∆t is required. However, ϕt+∆t is unknown for the current step; and (b) we can obtain ϕt+∆t from the fluid region to solve the issue. In the fluid region, the analogous obstacle is that u t+∆t is unknown. This problem can be addressed by using the partitioned approach based on Jacobi iteration [12], as shown in Table 2.
Based on the definition in Equation (45), the algorithms in Tables 1 and 2 are proposed. The presented method is a hybrid of the pseudo-parallel method and the partitioned approach, where these two methods are employed to address the memory requirement and iteration problems, respectively. It can be observed that the global stiffness matrix K, a semi-full matrix, is not used in the calculation; instead, it is established on an elementary level by the equationsf assembled in our process, the mass matrix M is a diagonal matrix, which can be treated as a vector, and its inverse can be easily calculated (using the reciprocal of each entry). The interface matrix C is similar to M but sparser. The presented calculations are all accomplished in only one grid system, where the absence of a staggered grid guarantees the simplicity of the structure of the developed algorithm.
The partitioned approach based on Jacobi iteration" End #2 End of the loop over instant t Output u = u u , ϕ = ϕ u The superscript "u" represents the "updating step". The superscript "p" represents the "predicting step". Table 2. The partitioned approach based on Jacobi iteration.

Calculate force vectors
The superscript "u" represents the "updating step". The superscript "p" represents the "predicting step". ω denotes the relaxation parameter, see Section 4.

Numerical Examples
Numerical validations were organized in the following manner: (1) Cases A-C were designed for different regions (i.e., from rectangular to solid circle) and different interfaces (i.e., from straight to curved). (2) Two sub-cases were set in each case (e.g., Case A-1 and Case A-2 for Case A). These two sub-cases were used to simulate wave propagation from solid to fluid and fluid to solid, respectively. (3) Cases A-C only present some qualitative comparisons with the theoretical wavefront in snapshots.
Quantitative analyses are given at the end of this section, in terms of convergence analysis.
(4) Considering the similarity of Cases A-C, convergence analysis was only conducted for Case C, the most complex case among them.

Case A: Rectangular with a Straight Interface
In order to demonstrate the performance of the presented method, a simple fluid-solid coupling example is conducted in this section. The inspected regions are treated as the plane strain problem and the material parameters are given in Table 3. In the following analysis, the primary wave (P-wave) and the secondary wave (S-wave) appear. The P-wave is a compressional wave which can move through both solid and fluid. The particles in the medium that the P-wave passes through are pushed and pulled by wave propagation. It travels faster than the S-wave. The S-wave can only move through the solid, not through any liquid medium (as it is a kind of shear wave). Particles are moved back and forth perpendicularly to the direction of wave motion. Due to the assumption of the present fluid medium, only the P-wave exists (1468.63 m/s) in the fluid part. The excitation is shown in Figure 2, which was defined by a modulated wave with the carrier frequency f 1 = 75 kHz and the modulation frequency f 2 = 75 kHz:  It should be noted that the BSWI elements used in this paper had 121 nodes per element, thus having 121 dofs and 242 dofs for fluid and solid, respectively. The number of dofs seemed to be not large enough, compared with some commercial software; however, it becomes further multiplied when time iteration (as shown in Tables 1 and 2) is considered. For each time step, enough memory was required to compute 60,903 dofs. Usually, more than 5000 steps are used in a typical simulation. The situation can be further compounded when a more complex shape and interface are investigated.
3.1.1. Case A-1: Wave Travels from Solid to Fluid Figure 3 presents four snapshots, from 10 to 50 μs, induced by a fixed source located at 0.191 m and 0.125 m in the solid region. The grey auxiliary curves present the theoretical wavefronts of the Pwave Cp, the S-wave Cs, and the secondary wavefronts. The auxiliary solid line located in the middle of the region depicts the interface. The symbol Cpp depicts the P-wave in fluid induced by the P-wave in the solid, and Csp depicts the P-wave in fluid induced by the S-wave in the solid. In the same manner, we can define the secondary wavefront Cps, which is presented in Figure 4. One can observe, from Figure 3, that the developed results qualitatively agreed with the wavefronts, either for the main modes Cp and Cs, or the secondary waves Cpp and Csp. The interface transformed the P-and S-waves The wave propagation is illustrated by contours in the following parts. In the solid region, the contours represent the displacement defined by u 2 x + u 2 y , while the contours in the fluid part are drawn based on ϕ . These contours are all normalized, for the sake of expression.
The solid part, with dimensions 0.125 m × 0.25 m, was divided by 10 × 20 BSWI elements-20,301 nodes and 40,602 dofs (degrees of freedom)-and the fluid part, with dimensions 0.125 m × 0.25 m, was also meshed by 10 × 20 BSWI elements-20,301 nodes and 20,301 dofs. It should be noted that the BSWI elements used in this paper had 121 nodes per element, thus having 121 dofs and 242 dofs for fluid and solid, respectively. The number of dofs seemed to be not large enough, compared with some commercial software; however, it becomes further multiplied when time iteration (as shown in Tables 1 and 2) is considered. For each time step, enough memory was required to compute 60,903 dofs. Usually, more than 5000 steps are used in a typical simulation. The situation can be further compounded when a more complex shape and interface are investigated.
3.1.1. Case A-1: Wave Travels from Solid to Fluid Figure 3 presents four snapshots, from 10 to 50 µs, induced by a fixed source located at 0.191 m and 0.125 m in the solid region. The grey auxiliary curves present the theoretical wavefronts of the P-wave C p , the S-wave C s , and the secondary wavefronts. The auxiliary solid line located in the middle of the region depicts the interface. The symbol C pp depicts the P-wave in fluid induced by the P-wave in the solid, and C sp depicts the P-wave in fluid induced by the S-wave in the solid. In the same manner, we can define the secondary wavefront C ps , which is presented in Figure 4. One can observe, from Figure 3, that the developed results qualitatively agreed with the wavefronts, either for the main modes C p and C s , or the secondary waves C pp and C sp . The interface transformed the P-and S-waves in the solid region to the corresponding secondary P-waves in fluid, where the wavefronts are also clear.

Case A-2: Wave Travels from Fluid to Solid
To bilaterally validate the effectiveness of the interface, the source was then moved to the fluid region for verification. The grid was kept the same. In Figure 4, the source was further implemented in the fluid region (at 0.083 m and 0.125 m), and the snapshots were recorded from 25 to 50 μs. Likewise, the numerical P-wave wavefront in fluid qualitatively agreed with the theoretical wavefront. When the wave passed through the interface, two wave modes were generated in the solid region (i.e., Cpp and Cps). The wavefronts of these two secondary waves also qualitatively agreed with the theoretical wavefronts. The above results validate the effectiveness of the present method in

Case A-2: Wave Travels from Fluid to Solid
To bilaterally validate the effectiveness of the interface, the source was then moved to the fluid region for verification. The grid was kept the same. In Figure 4, the source was further implemented in the fluid region (at 0.083 m and 0.125 m), and the snapshots were recorded from 25 to 50 μs. Likewise, the numerical P-wave wavefront in fluid qualitatively agreed with the theoretical wavefront. When the wave passed through the interface, two wave modes were generated in the solid region (i.e., Cpp and Cps). The wavefronts of these two secondary waves also qualitatively agreed with the theoretical wavefronts. The above results validate the effectiveness of the present method in

Case A-2: Wave Travels from Fluid to Solid
To bilaterally validate the effectiveness of the interface, the source was then moved to the fluid region for verification. The grid was kept the same. In Figure 4, the source was further implemented in the fluid region (at 0.083 m and 0.125 m), and the snapshots were recorded from 25 to 50 µs. Likewise, the numerical P-wave wavefront in fluid qualitatively agreed with the theoretical wavefront. When the wave passed through the interface, two wave modes were generated in the solid region (i.e., C pp and C ps ). The wavefronts of these two secondary waves also qualitatively agreed with the theoretical wavefronts. The above results validate the effectiveness of the present method in normalized regions.

Case B: Rectangular with a Curved Interface
Curved interfaces and boundaries are further investigated in this section. The material is presented in Table 4. The grid of Case B, a rectangular region divided by a curve, is shown in Figure 5. The left part was meshed by 560 wavelet elements (56,481 nodes) and the right part was meshed by 336 wavelet elements (34,001 nodes). Compared with the uniform grid for solid and fluid used in Figures 3  and 4, the mesh in Figure 5 was intentionally designed with different dimensions. As mentioned, the verification was hoped to be bidirectional. Thus, the left grid was used to model the solid region in Case B-1 (as shown in Figure 6), and then to model the fluid region in Case B-2 (as shown in Figure 7).       Figure 6 presents the wave propagation from solid to fluid. The left part was treated as solid and, thus, had 112,962 dofs; while the grid on the right of the interface was used to simulate the fluid, with 34,001 dofs. The source was fixed at 0.275 m and 0.125 m in the solid region. Due to the complexity of the theoretical secondary wavefronts, we only present the main modes in the following simulations. Figure 6a shows the overall view of the P-wave and S-wave, whose wavefronts qualitatively agreed with the theoretical ones. In Figure 6b, the P-wave passes through the interface and a clear secondary P-wave in the fluid region is generated after 37.5 µs. Figure 6c presents the transition from S-wave to the corresponding secondary P-wave in the fluid part after 50 µs. It can be seen that the snapshot is complex and blurred, due to reflections.

Case B-2: Wave Travels from Fluid to Solid
In the following simulation, shown in Figure 7, we exchange the solid and fluid grids, such that the left part was defined as fluid, with 56,481 dofs; while the right part was defined as solid, with 68,002 dofs. The source was still fixed at 0.275 m and 0.125 m, but was in the fluid region. Figure 7a illustrates the original P-wave. Figure 7b,c show the generation of the secondary waves when the original P-wave passes through the interface. As the source is near the focus of the interface, the reflection from the interface is nearly straight. The plane wave travels to the left and its shape remains as a line before 200 µs, as shown in Figure 7d-f.

Case C: Solid Circle with a Curved Interface
In Case C (Figure 8), a circular region with a curved interface was considered. The annulus was defined as the fluid region (960 elements, 160,800 dofs), while the inner circle was defined as the solid region (400 elements, 80,802 dofs). The grid in Figure 8a is not optimized, in view of the finite element method, as some elements nearly degenerate to triangles, as shown in Figure 8b. This means that the Jacobi matrix is nearly singular here. Furthermore, the dimension of the elements in the inner part of Figure 8a are not uniform, which involves the potential instability in temporal discretization and solution. Figure 8c presents an alternative grid mode for the inner region, with more uniform dimension and shape. However, in this section, we chose the worse grid for validation, based on the idea that, if the method performs well on a bad grid, it could perform well on better meshes. 3.3.1. Case C-1: Wave Travels from Solid to Fluid Figure 9 shows the wave propagating from solid to fluid in the inspected region. Figure 9a,b present the early development of the wave. The agreement between the theoretical wavefront and the numerical wavefront validated the effectiveness of the method in the solid region. In the snapshot at 25 μs (Figure 9b), the P-waves in fluid induced by the P-and S-waves in solid are concentrated near the interface. In Figure 9c-e, the secondary P-waves in fluid induced by the P-and S-waves in solid can be observed and distinguished. In addition, the presented result maintained the symmetry of the wavefield, as can be seen in Figure 9f. 3.3.1. Case C-1: Wave Travels from Solid to Fluid Figure 9 shows the wave propagating from solid to fluid in the inspected region. Figure 9a,b present the early development of the wave. The agreement between the theoretical wavefront and the numerical wavefront validated the effectiveness of the method in the solid region. In the snapshot at 25 µs (Figure 9b), the P-waves in fluid induced by the P-and S-waves in solid are concentrated near the interface. In Figure 9c-e, the secondary P-waves in fluid induced by the P-and S-waves in solid can be observed and distinguished. In addition, the presented result maintained the symmetry of the wavefield, as can be seen in Figure 9f. 3.3.1. Case C-1: Wave Travels from Solid to Fluid Figure 9 shows the wave propagating from solid to fluid in the inspected region. Figure 9a,b present the early development of the wave. The agreement between the theoretical wavefront and the numerical wavefront validated the effectiveness of the method in the solid region. In the snapshot at 25 μs (Figure 9b), the P-waves in fluid induced by the P-and S-waves in solid are concentrated near the interface. In Figure 9c-e, the secondary P-waves in fluid induced by the P-and S-waves in solid can be observed and distinguished. In addition, the presented result maintained the symmetry of the wavefield, as can be seen in Figure 9f.

Case C-2: Wave Travels from Fluid to Solid
Then, the source is further replaced at −0.0833 m and 0 m in the fluid region. The corresponding snapshots are given in Figure 10. We can see that creep waves were generated, due to the shape of the interface. The creep wave travelled along the interface, following it in a curved manner (as shown in Figure 10d-f). Induced by the creep wave, wave motions occurred in the area shadowed by the solid region.
Materials 2020, 13, x FOR PEER REVIEW 16 of 22 Then, the source is further replaced at −0.0833 m and 0 m in the fluid region. The corresponding snapshots are given in Figure 10. We can see that creep waves were generated, due to the shape of the interface. The creep wave travelled along the interface, following it in a curved manner (as shown in Figure 10d-f). Induced by the creep wave, wave motions occurred in the area shadowed by the solid region.

Convergence Analysis in Space
Convergence analyses were conducted on Case C for three types of grid: (1) Grid-1, the densest mesh composed of 1920 elements (321,600 dofs) for the fluid part and 800 solid elements (161,604 dofs) for the solid part.
(2) Grid-2, the grid used in Section 3.3, which was composed of 960 elements (160,800 dofs) for the fluid part and 400 solid elements (80,802 dofs) for the solid part. It should be mentioned that the grids were similar in structure, but had different densities. First, convergence analyses were conducted on Case C-1, in which the wave travelled from the solid part to the fluid part. The wave source was located at −0.0360 m and 0 m, and assigned vibrations along X-direction. The receiver was located at −0.0625 m and 0 m in the solid region. The fluctuation of the wave source determined that the X-direction responses at the receiver were dominated by the P-wave (as shown in Figure 11), and that the Y-direction here was dominated by the S-wave (as shown in Figure 12).
In order to illustrate the P-and S-waves clearly, the X-and Y-directional responses are shown in Figures 11 and 12, respectively. The highlighted part in Figure 11 defines the theoretical wavefront and duration of P-wave. It can be observed that the presented results, for all grids, agreed with the theoretical wavefront. A similar conclusion can be obtained from Figure 12 for the S-wave, in which the highlighted part defines the corresponding parameters for S-wave. Meanwhile, the fluctuation in

Convergence Analysis in Space
Convergence analyses were conducted on Case C for three types of grid: (1) Grid-1, the densest mesh composed of 1920 elements (321,600 dofs) for the fluid part and 800 solid elements (161,604 dofs) for the solid part. (2) Grid-2, the grid used in Section 3.3, which was composed of 960 elements (160,800 dofs) for the fluid part and 400 solid elements (80,802 dofs) for the solid part. It should be mentioned that the grids were similar in structure, but had different densities. First, convergence analyses were conducted on Case C-1, in which the wave travelled from the solid part to the fluid part. The wave source was located at −0.0360 m and 0 m, and assigned vibrations along X-direction. The receiver was located at −0.0625 m and 0 m in the solid region. The fluctuation of the wave source determined that the X-direction responses at the receiver were dominated by the P-wave (as shown in Figure 11), and that the Y-direction here was dominated by the S-wave (as shown in Figure 12).  Figure 13, in order to demonstrate the convergence of the method. It can be seen that the different grids achieved convergence in the fluid part.  Figure 13, in order to demonstrate the convergence of the method. It can be seen that the different grids achieved convergence in the fluid part. In order to illustrate the P-and S-waves clearly, the X-and Y-directional responses are shown in Figures 11 and 12, respectively. The highlighted part in Figure 11 defines the theoretical wavefront and duration of P-wave. It can be observed that the presented results, for all grids, agreed with the theoretical wavefront. A similar conclusion can be obtained from Figure 12 for the S-wave, in which the highlighted part defines the corresponding parameters for S-wave. Meanwhile, the fluctuation in the fluid at (−0.0700 m, 0 m) is presented in Figure 13, in order to demonstrate the convergence of the method. It can be seen that the different grids achieved convergence in the fluid part.
Further analysis was conducted on Case C-2, in which the wave travelled from fluid to solid. The source was located at −0.0833 m and 0 m, and the receiver was assigned at −0.125 m and 0 m in the fluid and at 0 m and 0.0625 m in the solid. Figures 14 and 15 show the convergence in the fluid part and solid part, respectively. For a better illustration of the creep wave, the response in the solid part is presented in the form of u 2 x + u 2 y , keeping consistency with the snapshots shown in Figure 10. It can be observed that the present method also converged for Case C-2.

Convergence Analysis in Time Domain
Besides the convergence analysis on space, the convergence analysis in time domain is also conducted in this section based on Gird-2. In the premise of avoiding the numerical diffusion in central difference time integration scheme, we vary the time step for presentation. Time steps are selected as 0.01 μs, 0.005 μs, and 0.0025 μs, denoted by "Step-1" to "Step-3", respectively. Figure 16 presents the results in solid and fluid. Compared with the convergence in space, the influence induced by the change of time step is tiny if it satisfies the convergency requirement of central difference convergence.

Convergence Analysis in Time Domain
Besides the convergence analysis on space, the convergence analysis in time domain is also conducted in this section based on Gird-2. In the premise of avoiding the numerical diffusion in central difference time integration scheme, we vary the time step for presentation. Time steps are selected as 0.01 μs, 0.005 μs, and 0.0025 μs, denoted by "Step-1" to "Step-3", respectively. Figure 16 presents the results in solid and fluid. Compared with the convergence in space, the influence induced by the change of time step is tiny if it satisfies the convergency requirement of central difference convergence.

Convergence Analysis in Time Domain
Besides the convergence analysis on space, the convergence analysis in time domain is also conducted in this section based on Gird-2. In the premise of avoiding the numerical diffusion in central difference time integration scheme, we vary the time step for presentation. Time steps are selected as 0.01 µs, 0.005 µs, and 0.0025 µs, denoted by "Step-1" to "Step-3", respectively. Figure

Discussions about Relaxation
As presented in Table 1, the relaxation method was not used; however, this does not mean that relaxation is not required in the present method. Thus, a discussion of the convergence of the partitioned approach is given in this section. The general formulation can be described by: The weighting constant ω is also called the relaxation parameter. The simplest method is to choose a fixed parameter ω for all time steps, in which the relaxation parameter has to be small enough to keep the iteration from diverging, but as large as possible to avoid unnecessary iterations. In addition, non-relaxation results can be obtained by setting ω = 1. Besides the fixed method, some self-adaptive relaxation methods are also available, such as Aitken relaxation [24]. The relaxation parameter is calculated by Equations (48) and (49) in the Aitken relaxation method: These relaxation methods can be integrated into the present method (see Table 2). To demonstrate the performance of these methods, the 1st-250th time steps of Case B-1 were selected. The convergence criteria es and ef (see Table 1) were used as inputs and the average number of iterations (denoted by "i" in Table 5) was defined as the output. Step-1 Step-2 Step - Step-1 Step-2 Step Step-1 Step-2 Step Step-1 Step-2 Step-3

Discussions about Relaxation
As presented in Table 1, the relaxation method was not used; however, this does not mean that relaxation is not required in the present method. Thus, a discussion of the convergence of the partitioned approach is given in this section. The general formulation can be described by: The weighting constant ω is also called the relaxation parameter. The simplest method is to choose a fixed parameter ω for all time steps, in which the relaxation parameter has to be small enough to keep the iteration from diverging, but as large as possible to avoid unnecessary iterations. In addition, non-relaxation results can be obtained by setting ω = 1. Besides the fixed method, some self-adaptive relaxation methods are also available, such as Aitken relaxation [24]. The relaxation parameter is calculated by Equations (48) and (49) in the Aitken relaxation method: These relaxation methods can be integrated into the present method (see Table 2). To demonstrate the performance of these methods, the 1st-250th time steps of Case B-1 were selected. The convergence criteria e s and e f (see Table 1) were used as inputs and the average number of iterations (denoted by "i" in Table 5) was defined as the output. It is clear that Aitken relaxation performed best among the inspected methods. From the convergence behavior of Aitken relaxation shown in Figure 17, one can observe that the iteration steps became stable swiftly after 100 steps, which is nearly the end instant of the wave source. These results demonstrate the adaptability of the Aitken method. By contrast, the average number of iterations for fixed relaxation methods were all linear with the value shown in Table 4 and, hence, are not presented here. For a simple problem such as that considered in this paper, even the non-relaxation method (ω = 1) can achieve a satisfying iteration speed, while under-relaxations converge slowly.  It is clear that Aitken relaxation performed best among the inspected methods. From the convergence behavior of Aitken relaxation shown in Figure 17, one can observe that the iteration steps became stable swiftly after 100 steps, which is nearly the end instant of the wave source. These results demonstrate the adaptability of the Aitken method. By contrast, the average number of iterations for fixed relaxation methods were all linear with the value shown in Table 4 and, hence, are not presented here. For a simple problem such as that considered in this paper, even the nonrelaxation method (ω = 1) can achieve a satisfying iteration speed, while under-relaxations converge slowly.

Conclusions
In this work, the inviscid fluid-solid coupling problem was solved by an iterative model using the wavelet element method and a partitioned approach. Based on the use of the Jacobian matrix in the wavelet element method, problems relating to curved boundaries can be solved. The interface was integrated into the model as a pseudo-damping layer, which can conduct transformations between force vectors in the solid region and pressure in the inviscid fluid region. As validated by our results, the presented method can serve as an alternative choice for the issue under consideration, both for normal regions and the curved region.
This work focused on inviscid fluid, but the problems relating to the influence of viscosity on the interface and numerical format (e.g., the implementation of boundary conditions and the iteration method used) can be further considered and investigated. To achieve efficient computation, fixed grid technology was used for the fluid region, by use of a Cartesian grid. As a result, the flow solver may be simple and fast; however, it sacrifices accuracy near the interface, due to the interpolations used. The use of a finite element method in both solid and fluid regions makes the formulation simple; however, the accuracy depends on the ratio between the fluid and solid grid size. Due to the effects of temporal iteration and numerical oscillation on convergence behavior, the sizes of elements are determined by the highest frequency domain considered, which may necessitate a memory-and time-related computational burden.

Conclusions
In this work, the inviscid fluid-solid coupling problem was solved by an iterative model using the wavelet element method and a partitioned approach. Based on the use of the Jacobian matrix in the wavelet element method, problems relating to curved boundaries can be solved. The interface was integrated into the model as a pseudo-damping layer, which can conduct transformations between force vectors in the solid region and pressure in the inviscid fluid region. As validated by our results, the presented method can serve as an alternative choice for the issue under consideration, both for normal regions and the curved region.
This work focused on inviscid fluid, but the problems relating to the influence of viscosity on the interface and numerical format (e.g., the implementation of boundary conditions and the iteration method used) can be further considered and investigated. To achieve efficient computation, fixed grid technology was used for the fluid region, by use of a Cartesian grid. As a result, the flow solver may be simple and fast; however, it sacrifices accuracy near the interface, due to the interpolations used. The use of a finite element method in both solid and fluid regions makes the formulation simple; however, the accuracy depends on the ratio between the fluid and solid grid size. Due to the effects of temporal iteration and numerical oscillation on convergence behavior, the sizes of elements are determined by the highest frequency domain considered, which may necessitate a memory-and time-related computational burden.