MPS-FEM Coupled Method for Study of Wave-Structure Interaction

Nowadays, an increasing number of ships and marine structures are manufactured and inevitably operated in rough sea. As a result, some phenomena related to the violent fluid-elastic structure interactions (e.g., hydrodynamic slamming on marine vessels, tsunami impact on onshore structures, and sloshing in liquid containers) have aroused huge challenges to ocean engineering fields. In this paper, the moving particle semi-implicit (MPS) method and finite element method (FEM) coupled method is proposed for use in numerical investigations of the interaction between a regular wave and a horizontal suspended structure. The fluid domain calculated by the MPS method is dispersed into fluid particles, and the structure domain solved by the FEM method is dispersed into beam elements. The generation of the 2D regular wave is firstly conducted, and convergence verification is performed to determine appropriate particle spacing for the simulation. Next, the regular wave interacting with a rigid structure is initially performed and verified through the comparison with the laboratory experiments. By verification, the MPS-FEM coupled method can be applied to fluid-structure interaction (FSI) problems with waves. On this basis, taking the flexibility of structure into consideration, the elastic dynamic response of the structure subjected to the wave slamming is investigated, including the evolutions of the free surface, the variation of the wave impact pressures, the velocity distribution, and the structural deformation response. By comparison with the rigid case, the effects of the structural flexibility on wave-elastic structure interaction can be obtained.


Introduction
Nowadays, an increasing number of ships and marine structures are being manufactured and inevitably operated in rough sea, the safe operation of which are strongly related to the wave motion. In the process of wave slamming, the impact pressures with short duration and high peaks will occur, causing the destabilization of the superstructure and local fatigue failure. Especially the very large floating structure (VLFS), whose stiffness is relatively small, may present behavior similar to that of an elastic plate suspended over the water surface. When encountering extreme condition, such as tsunami, the plate structures may produce elastic vibration and considerable deformation, leading to new challenges of structural safety. In the above-described situation, some non-negligible coupling effects exist in many fluid-structure interaction problems. Thus, the studies of fluid-structure interaction problems become of paramount importance in the field of ocean engineering, requiring precise and comprehensive modeling of these important fluid-structure interaction (FSI) problems.
Until very recently, only a few existing laboratory experiments focused on the interaction of periodic waves with a horizontal plate, such as those conducted by Cuomo et al.

Article Highlights
• A coupled method of MPS-FEM is proposed for the interaction between a regular wave and a horizontal suspended structure.
• A more flexible weak coupling strategy is employed between the fluid and structure solver.
• At the interface between the fluid and structure domain, the displacement and force equivalence condition are satisfied.
• By comparison with the rigid structure, the effect of structure flexibility is investigated. (2007) and Schumacher et al. (2008), or on the solitary wave forces for a flat plate, such as those conducted by Seiffert et al. (2014) and Lo and Liu (2014). These experiments are of great value for further research; however, few of them considered structural flexibility. Nelli et al. (2017) conducted an experiment to study wave impact on a thin flexible plate to analyze the reflection and transmission of regular incident water waves. A conclusion can be drawn that the wave height is attenuated by the green water and wave breaking and that the attenuation is low when the plate is more flexible. As is well known, laboratory experiments and theoretical research studies on the field of FSI problems have some limitations because of the characteristics of obvious instantaneity and non-linearity; thus, CFD has been proven to be an effective tool to solve FSI problems. For the past decades, mesh-based methods have predominated studies on this topic, some of which have been successfully applied to model the intricate phenomena involved in FSI problems. For example, Idelsohn et al. (2008) simulated the interaction of sloshing flow with a bottom clamped elastic baffle using the particle finite element method (PFEM). The simulation results were found to be consistent with the experimental data. He and Kashiwagi (2010) studied the hydroelastic behavior of the interaction between a wave and a linear elastic plate based on boundary element method (BEM) for fluid analysis coupled with finite element method (FEM) for structure analysis. A finite differential method (FDM) and FEM coupled approach was applied by Liao and Hu (2013) to investigate the interaction between surface flow and a thin elastic plate. The results showed that the coupled method has great accuracy. Despite the effectiveness in mesh methods, taking violent fluid flows and deformable structures into account, these methods may suffer from many difficulties, such as the adjustment and regeneration of mesh while coordinating the interface between fluid and solid domain.
In view of the characteristics of FSI problems, Lagrangian mesh-free methods, e.g., smoothed particle hydrodynamics (SPH) developed by Gingold and Monaghan (1977) and Lucy (1977) or moving particle semi-implicit (MPS) proposed by Koshizuka and Oka (1996), have shown great capabilities for addressing such problems. In particular, the fluid flow is discretized into moving particles in mesh-free methods, thereby allowing some difficulties, such as the calculation of the convection term, the treatment of mesh, the capture of free surface, and large structural deformation, to be handled easily. In recent years, several studies have been conducted on implementation of these Lagrangian particle methods into FSI simulations. Antoci et al. (2007) performed simulations of dam-break flow interacting with an elastic gate using the SPH method. Khayyer et al. (2018aKhayyer et al. ( , 2018b presented enhanced fully Lagrangian mesh-free MPS-and SPH-based FSI analysis methods. Because of the multiple physics fields involved in FSI problems, the use of both a flow field and a structure domain is necessary to couple the dynamics of the two domains. Thus, many researchers have taken advantage of meshfree methods along with other robust computational methods to develop coupled FSI solvers. Since the 1990s, various SPH-FEM coupling methods have been proposed (Johnson 1994;Attaway et al. 1994) and subsequently used in solving FSI problems (Fourey et al. 2017;Canelas et al. 2016). Compared with the SPH, the pressure of the particle is obtained by solving the pressure Poisson equation (PPE) in the MPS method. Thus, the obtained pressure field using the MPS method is relatively smooth. The MPS method integrated with the FEM exhibits good performance in the FSI problem according to the numerical benchmark test of dam-break flow interacting with an elastic structure (Mitsume et al. 2014;Hashimoto and Le Touzé 2014;Sun et al. 2016).
In the present paper, we aim to apply the MPS-FEM coupled method to FSI problems of a regular wave slamming on a flexible structure. A fully Lagrangian solver MPS-FEM-STJU is developed for fluid-structure interaction problems based on the abovementioned method. The solver has been proven to be valid in 2D FSI problems, such as liquid sloshing with an elastic baffle and dam-break through an elastic gate (Zhang et al. 2016a;Zhang and Wan 2017), as well as 3D FSI problems, including solitary wave-plate interaction (Rao et al. 2017a;Rao et al. 2017b) and sloshing in an elastic tank (Zhang and Wan 2018). Theories of the MPS and FEM method together with the  coupling strategy are first presented. Next, the numerical wave tank is developed to calculate the interaction between waves and structure. The time history of wave propagation is measured and compared with the analytical solution to validate the accuracy of wave making. The interaction between the wave and the rigid plate is conducted firstly, and then the numerical impact pressure results are compared with the experimental results of Wang and Ren (2002) and Gao et al. (2012). Next, we simulate the regular wave impinges on a horizontal flexible plate. The case is compared with that of the rigid case to investigate the effects of flexibility on the slamming process.

Numerical Method
In this study, the MPS-FEM coupled method is adopted to address the wave-structure interaction problem. The MPS method for fluid domain analysis and the FEM for structural domain analysis are firstly introduced. The theories for the MPS, which have been presented in detail in our previous papers (Zhang and Wan 2012;Zhang et al. 2014;Tang and Wan 2015;Tang et al. 2016;Zhang et al. 2016b;Wan 2017, 2018), will be introduced briefly in this section. Next, a partitioned coupling strategy is adopted for the treatments of data transformation on the fluid-structure interface.

MPS Model for the Fluid
The governing equations for viscous incompressible fluid including the continuity equations and the Navier-Stokes equations are expressed in Lagrangian form as follows: where V, ρ, P, ν, and g denote the velocity vector, the fluid density, the pressure, the kinematic viscosity, and the gravitational acceleration, respectively. In the particle method, the governing equations should be expressed by the particle interaction models based on the kernel function. The kernel function in this paper can be formulated as where r = | r j − r i | is the distance between particles i and j and r e denotes the influence radius of the target particle.
In the MPS, the particle interaction models including the gradient, divergence, and Laplacian models are defined as Eqs. (4)-(6).
where ϕ is an arbitrary scalar function, Ф is an arbitrary vector, D is the number of spatial dimensions, r is the position vector, λ is a parameter that is expressed as Eq. (7), and n 0 is the initial density of the particle number defined as In this paper, the incompressible condition of the MPS method is represented by two parts: the particle number density and the divergence of velocity (Tanaka and Masunaga 2010;Lee et al. 2011). Pressure is implicitly calculated by solving the pressure Poisson equation (PPE), and the velocity and position of the particles are updated according to the obtained pressure. The PPE in the present MPS solver is defined as where γ is the weight of the particle number density term, which has a value in the range between 0 and 1. The range of 0.01 ≤ γ ≤ 0.05 is better according to numerical experiments conducted by Lee et al. (2011). In this paper, γ = 0.01 is adopted for all numerical experiments. The pressure of the fluid domain is closely affected by the accuracy of free surface detection. In this paper, we employ a free surface detection method by Zhang et al. (2014) that is defined as where F represents the asymmetry of the arrangements of neighboring particles. Particles satisfying < | F | > i > α|F| 0 will be considered a free surface particle, where |F| 0 is the initial value of |F| for surface particles. MPS is a fully Lagrangian method, where the boundary condition of the free surface is naturally satisfied. For the solid boundary, the treatment of multilayer particles is shown in Fig. 1. This treatment can ensure a smooth and accurate pressure field around the solid surface and prevent fluid particles from penetrating into the impermeable boundary. The solid boundary is represented by one layer of wall particles. The pressures of these particles and fluid particles are solved by PPE. The function of two layers of ghost particles is to fulfill the particle number density so that the particle interaction can be properly calculated near the solid boundary. The pressures of these ghost particles are obtained by interpolation.

FEM Model for the Structure
Based on FEM theory, the spatially discretized structural dynamic equations that govern the motion of structural nodes can be formulated as where M is the mass matrix, C represents the Rayleigh damping matrix, and K denotes the stiffness matrix. F represents the external force acting on the structure and varies with computational time. y denotes the nodal displacement of the structure. Coefficients α 1 and α 2 correspond to the natural frequency and the damping ratios of the structure, respectively. To solve Eq. (10), another two groups of functions should be added to form the equation system. Herein, Taylor's expansions of velocity and displacement developed by Newmark (1959) where β and γ are paramount parameters in the Newmark-β method and are set as β = 0.25 and γ = 0.5 for all simulations. Next, the displacement at t = t + Δt can be solved by the following formula proposed by Hsiao et al. (1999): where K À and F À denote the effective stiffness matrix and the effective force vector, respectively. Subsequently, the accelerations and velocities related to the next time step are updated as follows:

Fluid-Structure Coupling Scheme
In general, the coupling strategies for the FSI simulations can be classified into the strong coupling approaches and the weak coupling approaches. In the present study, the more flexible weak coupling strategy is employed for the FSI analysis because the separated fluid and the structure codes for each computational domain are admitted. Khayyer et al. (2018a) highlighted the importance of precise implementation of fluid-structure coupling schemes for stable and accurate simulations, especially in longterm simulations. The concept of the coupling system is shown in Fig. 1; the time step sizes for structure and fluid analyses are Δt s and Δt f , respectively. Because of the stability of the Newmark-β scheme, Δt s can be k multiples of Δt f , where k is an integer. The interaction procedure can be summarized as follows: 1) The positions of both fluid and structure particles (Fig. 2) will be updated based on the displacement response y n calculated by the structure module as well as velocity ẏ n and acceleration ÿ n .
2) At each fluid time step, the pressure of the fluid wall boundary particle is calculated. To obtain more accurate external loads for the structural response analysis, the mean values of particle pressures on the fluid-structure interface are calculated by where p n + i is the pressure of the fluid particle on the structural boundary at the instant t n + i .
3) The force vector F tþΔt s acting on the structural boundary particle is calculated by multiplying the average pressure p nþ1 and the influential area. Next, the structural displacement response can be determined based on the Newmark-β scheme.
4) The positions of particles should be updated by Eq. (21) at every fluid time step Δt f , not only at each structural time step Δt s , to avoid the instability of fluid field produced by a larger displacement within a larger time interval Δt s .
where ẏ nþi−1 and ÿ nþi−1 are the structural nodal velocity and acceleration at the time t n+i−1 , respectively.

Data Interpolation on the Fluid-Structure Interface
In this paper, the fluid domain is dispersed as particles and the structure is dispersed as beam elements. Thus, special schemes are required for data transformation on the fluidstructure interface, including applying the external force carried by the fluid particles onto the beam nodes and updating the position of boundary particles corresponding to the displacements of the beam elements. More specifically, the displacement equivalence condition and the force equivalence condition must be satisfied. A particle group scheme developed by Hwang et al. (2016), in which the structural particles located within the same section are grouped, is considered. For the force, the external load on the structure boundary is equal to the load on the fluid boundary at the interface. The concepts of the numerical considerations are shown in Fig. 3, where the vectors F Gi, upper and F Gi, lower represent the force acting on the upper and lower boundary particle of the structural group i, respectively. As stated above, the pressure of boundary particle is calculated using the MPS method initially. Next, the force acting on the structural boundary particle is calculated by the integration of the average pressure acting on the interface. Afterwards, the resultants within the same group are applied onto the structural FEM node as the external force for the structural physics analysis. For the deformation, the deformation of the fluid boundary is equal to the deformation of the structure boundary at the interface. Particles within a group move as one body according to the nodal linear velocities u i and v i . Next, the final position of structure particles can be updated according to the angular velocity ω i . The concepts of the deformation of structural particle model are shown in Fig. 4.

Validation of Wave Making
In this subsection, the convergent verification of a regular wave is conducted to determine the appropriate particle spacing. The model for simulation is depicted in Fig. 5. A pistontype wave generator is incorporated in the left side of the 2D numerical wave tank. A sponge layer 1.5 m long is installed at the end of the wave tank to absorb waves and avoid reflection. The NWT is of 7.5 m width and 1.5 m height, with an initial water depth of 1.05 m. The rigid fixing plate 3.5 m away from the piston paddle will be placed in the next section. The wave conditions used in the present numerical test are shown in Table 1. The regular wave adopted in the simulations is based on linear wave theory. The wave period (T) is 1.2 s, and the wave height (H) is 0.15 m. The particle spacing values used in the regular wave simulation are 0.0075 m, 0.01 m, and 0.015 m, with the corresponding total particle numbers of 144 336, 82 720, and 39 339, respectively. The simulation time is 30 s, and calculation time step is 0.0002 s. The waves are measured at the location of the leading edge of the plate structure without the plate to measure the incoming wave height and period and determine the repeatability of the generated waves. Figure 6 shows the results of numerical wave elevation with different particle spacings versus

Wave Impacting a Rigid Plate
In this subsection, the problem of a wave interacting with a rigid horizontal plate is investigated using the in-house MPS-FEM-SJTU solver. The structural model is 1.0 m long and 0.03 m thick, on which eleven probe points are set, marked as #1 to #11 starting from the leading side of the structure, as shown in Fig. 7. Some structural parameters are shown in Table 2. The clearance (S) signifies the distance between the still water level and the subface of the structure. In this paper, the clearance is set to 0.015 m, 0.03 m, and 0.06 m, such that the relative clearance (S/H) is 0.1, 0.2, and 0.4, respectively. In the simulation, the water depth is set as 1.05 m, which is different from the experiment. The difference is presented in Table 3. The relative error is approximately 5%. Nevertheless, the error may not be ignored. Therefore, in this paper, we focus on the characteristics of the impact process, including the features of the free surface and the impact pressure, rather than the value of the impact pressure. Figures 8, 9, and 10 show the evolution of the free surface, including the instantaneous pressure and instantaneous velocity distributions of the fluid particles within an impacting period, when the relative clearance (S/H) is 0.1, 0.2, and 0.4, respectively. In Fig. 8, note that, when the fluid flow impinges onto the subface of the structure, the wave flow has a significant reflection, the velocity vector of which near the leading side of the structure is almost perpendicular to the structure. Thus, evident green water can be observed. When the wave is entirely under the structure, the fluid particles move forward along the structure model according to the horizontal velocities, and the velocity vector of the anterior flow is inclined to the structure due to the slamming, whereas the rear flow appears to be away from the structure. Finally, when the wave crest separates from the structure, an evident spray can be observed at the trailing edge of the structure. With the increase of the relative clearance, evident differences can be observed during the slamming process. In Fig. 10, when wave impinges on the structure, there is no evident reflection; thus, the velocity vector is only inclined to the structure. In addition, the  Fig. 13 The comparison of impact pressure by simulation with experiment data (S/H = 0.1) green water is slight and intermittent. With increasing relative clearance, the larger the distance, the earlier the fluid detaches away from the structure (as Figs. 8(d),9(c), and 10(c)). Finally, in Fig. 10(d), the wave crest cannot slam against the trailing edge of the structure; this behavior may be attributed to the energy dissipation during the slamming process.
The time variations of the impact pressure on the subface of the structure with the relative clearances of 0.1, 0.2, and 0.4 are shown in Figs. 11 and 12 and the left column in Fig. 13, respectively. The detailed drawings of the impact pressure with S/H of 0.2 and 0.4 are also given. When the wave interacts with the structure, an impulsive impact pressure with high peak is generated. Next, the pressure drops suddenly, followed by a period of small oscillations of the pressure. The closer the pressure probe is to the trailing edge of the structure, the shorter the duration of the oscillations. In the end, the pressure returns to the initial zero pressure. An evident impact pressure period of approximately 1.2 s can be observed that is consistent with the incident wave. The peak value of wave impact pressure in different periods is randomized.
In the experiments of Wang and Ren (2002), a plate structure (L = 1 m, B = 0.69 m, and t = 0.02 m) is located at different clearances, and the impact pressures of a series of periodic waves of different heights and periods are considered. Here the case with T = 1.2 s, H = 0.15 m, and S/H = 0.1 is chosen for comparison. The comparison of the results in the present paper and the laboratory measurements by Wang and Ren (2002) and Gao et al. (2012) are shown in Fig. 13. In the simulation, the calculation time is 30 s, and the time starts from 10 s after the start of the wave paddle process to generate the regular waves shown in Figs. 10, 11, and 12. However, in the experiment of Wang and Ren (2002), the time started from the moment that the wave reached at the plate. The time spans are consistent with each other. In the process of the wave slamming mentioned above, periodicity, randomness, and impulsiveness are the obvious characteristics of impacting pressure, which are consistent with previous works (Wang and Ren 2002). A good agreement is observed overall. The trend of wave-induced pressure history agrees with the experimental results. This agreement indicates that the solver adopted in the present paper can be effectively applied to the wavestructure interaction problems.
The average peak impact pressure on the subface of the structure can be obtained by averaging the peak value in each period. The values versus different clearances are depicted in Fig. 14 to investigate the impact pressure distribution along the subface of the structure and its variation with the relative clearance (S/H). The variation of the peak value with relative clearance appears to be intricate. When S/H is 0.1 and 0.2, the peak value of impact pressure is distributed unevenly and varies greatly from point to point. Most of the maximum impact pressure is found to appear in #3 and #4, approaching the leading side. When S/H is 0.4, the peak value presents a uniform distribution. When S/H is 0.1, the peak values are significantly higher than those of the other two cases. However, when S/H is 0.2 and 0.4, the peak value exhibits no obvious relationship. In the future, we will investigate more cases to better understand the variation of the peak value for different relative clearances.

Wave Impact on a Flexible Plate
To study the effects of structural flexibility on the slamming process, the simulation of the regular wave interacting with a flexible plate is conducted in this subsection. The elastic modulus (E) of the structure is 10 MPa. The arrangements of the numerical flume and the plate are the same as those of the previous section. The relative clearance (S/H) is 0.4; thus, the primary impact is bottom slamming. The simulation time is 20 s. Figure 15 shows the evolution of the free surface while taking the structural flexibility into consideration, including the instantaneous pressure and instantaneous velocity distributions within an impact period.
When the wave propagates to the structure, the leading side of the structure is bottom slammed and the velocity vector of water particles is inclined to the structure, resulting in an upward deformation of the leading side of the structure. With the wave transmitting forward along the subface and impinging onto the rear structure, the upward deformation moves forward along the wave transmission. The wave can slam against the trailing edge of the structure slightly, despite the energy dissipation. Before the next incoming wave arrives at the structure, a vibration can be observed in Fig. 15(e). The structural flexibility has a great influence on the evolution process. When the structure is elastic, the peak pressure is too small to be observed. Moreover, the energy dissipation in the flexible case appears to be smaller than that in the rigid case; this observation is similar to the conclusion of Nelli et al. (2017).
The displacement history in the middle of the structure and the related impact pressure history are shown in Fig. 16. When the wave arrives at the structure, the structure possesses an upward deformation immediately because the water impinges on the subface. With the wave transmitting forward and separating from the structure, the structure vibrates to the equilibrium position under the action of the elastic restoring force of the structure. Because of the lack of storage of the green water, the equilibrium position of the structure is always the initial position. The vibration with evident attenuation will remain a fixed number of times until the next incoming wave impacts the structure. It can be observed that the peak value of wave impact pressure is more randomized than the rigid case; the displacement of the structure might be the reason for the difference. The moment that the structural deformation response occurs is consistent with the time of force action and presents an obvious periodicity. Furthermore, the peak value of the impact pressure is much lower than that of the rigid case for the #6 point in Fig. 12, and the impact duration is less than half of that of the rigid case.

Conclusions
In this paper, a fully Lagrangian FSI solver was primarily implemented based on the MPS-FEM coupled method. Numerical simulations of the interaction between the regular wave and the horizontal flexible structure were conducted. The interaction between the wave and rigid structure was conducted firstly. The feasibility of the present solver for FSI problems was validated using the test of a regular wave impacting onto a rigid structure. Good agreements between the present numerical results and the laboratory measurements were achieved. By analyzing the evolution of free surface, including the instantaneous pressure and instantaneous velocity distributions of the fluid particles and the time variations of the impact pressure on the subface of the structure, the characteristics of the wave slamming structure were clearly observed. Moreover, the effects of the relative clearance on the slamming process and the impact pressure were initially investigated. The results indicate that the solver adopted in this paper can be effectively applied to the wave-structure interaction problem.
Subsequently, the FSI problem of the regular wave interacting with a flexible structure was studied. We chose the case with the relative clearance (S/H) of 0.4 that has the primary impact of bottom slamming. The deformation of the structure was found to present an obvious periodicity, with the period related to the wave slamming. By comparison, the impact duration time when slamming on a flexible structure was found to be less than half of that regarding the rigid case. By analyzing the quantitative results, the impact pressure induced by the wave impacting onto the flexible structure was found to be much smaller than that regarding the rigid case. Furthermore, the peak value of wave impact pressure appeared to be more randomized, which could be caused by the displacement of the structure. By comparison of the velocity field, the energy dissipation in the flexible case appeared to be smaller.
In summary, the MPS-FEM coupled method can be applied to FSI problems with waves. Moreover, the phenomena and characteristics of wave-structure interaction can be observed more clearly by means of numerical simulation.