Critical grid method: An extensible Smoothed Particle Hydrodynamics fluid general interpolation method for Fluid-Structure Interaction surface coupling based on preCICE

In the SPH method, for any continuous field function $f(\bm{r})$, the function value at a point on the definition domain $\Omega$ can be expressed by the following equation (1).

Where $\bm{r}$ is the vector in $n$-dimensional space ($n=1,2,3,\cdots$) and $\delta(\bm{r}-\bm{r}^{\prime})$ is Dirac $\delta$ Function, which satisfies the following properties.

However, Dirac $\delta$ function is discontinuous in reality, so the smooth function $W\left(\bm{r}-\bm{r}^{\prime},h\right)$ is used to approximately replace the $\delta$ function, which obtains the approximate kernel function expression of the field function $f(\bm{r})$.

$W(\bm{r}-{\bm{r}}^{\prime},h)$ represents the smooth kernel function in SPH method, and its value depends on the distance $\left|\bm{r}-{\bm{r}}^{\prime}\right|$ for two points and the smooth length $h$, and its influence region is determined by the smoothing factor $\kappa$ and the smoothing length $h$.

Based on the same principle, the kernel approximation of the derivative function of $f(\bm{r})$ can be evaluated according to the integration formula and the divergence theorem.

Equations (4) and (5) show that any field function and its derivatives in space can be expressed by the smooth function $W$. However, equation (5) can only show that the support region completely falls within the boundary range. When the support region is truncated by the boundary, boundary conditions need to be imposed. Therefore, imposing boundary conditions is also a hot and difficult topic in the theoretical research of the SPH method.

Equations (4) and (5) obtain the integral expressions of the field function and its derivatives based on the smooth kernel function $W$. But it can not discretize the governing equations of the fluid. Therefore, the SPH method uses the particle approximation method to interpolate and discretize the smooth kernel function. Its discretization principle is as follows.

In the solution domain $\Omega$, as shown in Fig.1, the integral expression of the kernel function can be further transformed into the form of weighted summation of a series of discrete particles in the support domain, which is called the particle approximation of the SPH method. If $V_{j}$ is used to represent the volume element $d{\bm{r}}^{\prime}$ of particle $j$ with arbitrary distribution, the mass of particle $j$ can be expressed by equation (6).

Where $\rho_{j}$ is the density of particle $j$. Based on equations (4) and (6), the particle estimation expression for the field function at particle $i$ can be solved as shown in equation (7).

Similarly, according to equations (5) and (6), the derivative expression of the field function at particle $i$ can be obtained as shown in equation (8 or 9).

Where

$r_{ij}$ represents the distance between particles $i$ and $j$, then, an excellent feature of equation (11) is that the right side of the equation appears in the form of paired particles. The advantage of this symmetric form is that it can improve the calculation accuracy of the SPH method 2006Restoring .

In the field of Computational Fluid Dynamics, the momentum conservation equation controlling fluid flow in a continuous field can be described by equation (12).

Where $\bm{\Gamma}$ is the dissipative term and $\bm{g}$ is the gravitational acceleration ($\bm{g}=9.8m/s^{2}$). Considering the effect of fluid viscosity and surface tension, artificial viscosity was proposed by Monaghan to simulate fluid flow because it is easy to realize. Equation (12) can be written as equation (13).

Where $P$ and $\rho$ are the pressure and density of the corresponding particle $i$ or $j$, respectively, and $\Pi_{ij}$ is the viscosity term, which is given by equation (14).

Where $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$, $\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{j}$, $\mu_{ij}=h\bm{v}_{ij}\cdot\bm{r}_{ij}/(r_{ij}^{2}+\eta^{2})$, $\bar{c_{ij}}=0.5(c_{i}+c_{j})$, $c$ is the speed of local sound, $h$ is the smooth length, $\eta^{2}=0.01h^{2}$, $\alpha$ is a coefficient associated with the problem.

In addition, in most application scenarios, the laminar viscous stresses also needs to be introduced into the momentum equation Edmond2002Simulation , which can be written as

Where $\upsilon_{0}$ is the kinetic viscosity, generally, its value is $10^{-6}m^{2}s$ for water. Therefore, according to the interpolation principle of the SPH method, equation (15) can be transformed into equation (16) with discrete form of the SPH method through equations (7) and (8).

Gotoh2001Sub first introduced the concept of Sub-Particle Scale (SPS) to describe the affected by turbulence for their Moving Particle Semi-implicit (MPS) model. Their momentum conservation equation is defined as follows.

Where equation (15) is used as the laminar term, and $\vec{\tau}$ is the SPS stress tensor. Finally, Rogers et al.dalrymple2006numerical introduce SPS into the weakly compressible SPH (WCSPH) method by using Favre averaging, and the SPH discrete form of equation (17) can be written as

In the process of solving weakly compressible SPH simulations, the mass of particles remains unchanged, the initial density solution is obtained by weighted summation of particle masses monaghan1992smoothed , which tends to cause violent fluctuations near the boundary and free surface. Therefore, the density can be solved currently by the conservation equation of mass, or continuity equation as shown in equation (19).

In the weakly compressible SPH method, the pressure is solved by the equation of state (EOS) monaghan1994simulating . Therefore, the relationship between pressure and density is described in equation (20).

Where $\gamma=7$ for water, $\rho_{0}=1000kg/m^{3}$ is the reference density, and $c_{0}=c(\rho_{0})=\sqrt{\partial P/\partial\rho}$ is the local speed of sound when $\rho=\rho_{0}$.

preCICE provides two coupling schemes to realize the coupling process: 1) explicit coupling schemes and 2) implicit coupling schemes. The difference lies in whether all physical quantities are obtained at the same time. The explicit coupling scheme needs to be called in a single time step, while the implicit coupling scheme needs to iterate the coupling equation until it converges, so it has iterative subcycles. Taking the solution of two solvers on FSI problem as an example, the fluid solver is annotated as $F_{1}$, and the solid solver is denoted as $S_{2}$. The data of vectorization type exchanged between the two solvers are $X_{1}$ and $X_{2}$, respectively. Therefore, they have the following corresponding relationship on the coupling surface.

Most of the partitioned coupling schemes belong to Dirichlet–Neumann type coupling denk2015structurally . In FSI problem, fluid solver $F_{1}$ needs to provide forces $X_{1}$ as the input of the wet interface and transmit it to solid solver $S_{2}$ through preCICE. Then, solid solver $S_{2}$ sends back displacements $X_{2}$ to fluid solver $F_{1}$. Therefore, boundary conditions and conservation equations are also satisfied.

preCICE supports the coupling of multiple solvers. Generally speaking, different solvers adopt different meshes on the coupling surface (many-to-many relationships). It requires data mapping of unmatched meshes on the coupling surface, and this mapping relationship can not destroy the basic conservation equations such as mass conservation and energy conservation, otherwise, the simulation results are difficult to converge. Therefore, preCICE provides consistency and conservatism principles to deal with data mapping between unmatched grids of Fig.5 (a).

Consistent mapping and conservative mapping are used in different situations. First of all, the consistency form just copies the data from one grid to another grid of the receiver, and its value does not change. For example, in heat transfer simulation, the change of temperature in a certain region is consistent. As shown in Fig.5 (b), the temperature between adjacent grids should use the consistent mapping principle, which ensures that the corresponding grid temperature within the same area is consistent (the grid is the basic unit of atomic operation). In addition, Newton’s third law states that the force and its reaction are equal in magnitude. For solving Fluid-Structure Interaction problem, the transfer of force should strictly follow the conservative mapping principle. As shown in Fig.5 (c), the resultant force on both sides of the fluid and solid should remain unchanged. Therefore, the forces between adjacent mesh interpolation points shall be evenly distributed, such as $F=4N=2N+2N$ and $F=6N=3N+3N$.

preCICE provides the following methods and their variants to implement the data mapping function, which can ensure the conservation conditions between the coupling surfaces.

• •

  Nearest neighbor: Its principle is to find the nearest mesh points for interpolation in space. It does not need any topology information. This method is simple and easy to implement, but its disadvantage is that it only has first-order accuracy.

• •

  Nearest projection: In this method, the target grid points are projected onto the source grid elements according to the grid topology information, and the source grid elements are linearly interpolated according to the coordinates of the projection points. Because the orthogonal distance between the coupling interface and the grids on both sides is very small, therefore, this method has second-order accuracy.

• •

  Radial Basis Function (RBF): This method does not need either grid topology information or projection information. It constructs an interpolation function (RBF) centered on the source grid points. Through this function, the target grid points can be interpolated and the interpolation results can be obtained. Although RBF is globally supported, considering the complexity and efficiency of calculation, its application in practice is only limited to a small range. At the same time, preCICE also provides a variety of interpolation functions that can be consulted in literaturebungartz2016precice , among which Gaussian and Thin Plate Spline are the most widely used RBF interpolation functions.

It is a great challenge to couple multiple solvers with parallel computing for efficient communication. Each participant may be executed by multiple processes and distributed on multiple nodes of the supercomputer. Therefore, preCICE provides MPI technology to solve the problem of large-scale parallel communication. It allows open-source programs to run efficiently across nodes. However, for closed-source commercial software, no corresponding interface is provided for calling. preCICE provides a lower-layer TCP/IP sockets communication protocol to support this function which is difficult to program. Each participant communicates point-to-point in parallel, and then it sets up a master process to steer other slave processes. These communication relationships are set at the time of initialization. They are configured through an XML file. preCICE will automatically recognize the mesh decomposition and establish contact with the corresponding processes. This communication method also has some disadvantages, for example, it is difficult to deal with adaptive mesh and FSI problems of fluid immersed boundaries.

When using preCICE for coupling development, each participant (solver) needs to make appropriate modifications. Among them, the solver interacts with preCICE through the adapter, as shown in Fig.2. Therefore, the corresponding adapter for fluid or solid must be developed to use preCICE for coupling. Although preCICE was written in C++. However, the interface it provides also supports C, Python, Fortran and other widely used programming languages. This highly modular software architecture enables it to be flexibly compatible with most of the current influential CFD/CSD software, including open-source and closed-source commercial software. At present, the compatible open source software has been released on the official website (www.precice.org/adapters-overview).

At present, preCICE has integrated many open/closed-source software including those enumerated in Table LABEL:tab01, and it has achieved many successful applications in FSI rodenberg2021fenics , conjugate heat transfer kim2020simulation , multi-phase flow bungartz2011 , and acoustic fields bungartz2016 . It has a strong function and rich development documentation, at the same time, software community users very active support the preCICE updating constantly, and it attracts more and more researchers worldwide to improve and strengthen the function of preCICE. In future work, the development team will provide the support to change the interface grid.

In section III of this paper, we mentioned that the preCICE coupling library and the solver are connected through adapters. This highly modular approach makes preCICE compatible with most solvers. Similarly, for the SPH method based on pure Lagrange, we use PMC method to develop an adapter connecting preCICE, which is called SPH-Adapter. It uses a critical grid as the medium to interact with preCICE and solid solver, and its process can be simplified as shown in Algorithm 1.

Lines 1-2 of Algorithm 1 are the configuration of preCICE, where precice-config.xml is a configuration file with XML format, which contains a lot of configuration information. Including the configuration of fluid and solid, the adopted data mapping method, data communication mode, setting time window and convergence conditions, etc. The configuration file needs to be loaded into a path recognized by the SPH solver. Then, set a critical mesh in line 3 as the medium of particle interpolation, which meets all the standards of the preCICE interface and is composed of only one layer of mesh as shown in Fig.7 (f). Its mesh vertices are loaded into the array in the form of coordinates according to the topological relationship, and preCICE provides the functions of mesh connection and data initialization to support these operations. Meanwhile, the adjacency relationship between grids has been determined at the beginning and will not be changed later. After the configuration is completed, line 4 calls the initialization function, and then preCICE will establish communication between the fluid and solid solver. SPH is solved in line 10, which is responsible for the calculation of fluid in the FSI problem. Line 7 is to obtain the displacement returned by the solid solver based on the critical grid, and line 12 is to transfer the force generated by fluid particles to the solid structure based on the critical grid. Data transmission and coupling communication only occur in the advance function in line 11. In addition, in line 9, it is necessary to ensure that the time step used between fluid and solid is consistent, which avoids the calculation error caused by the disordered time advance.

In summary, we established a link between the preCICE framework and SPH particles through a critical grid, and we will describe in detail the process by which the interpolation is achieved for the critical grid and particles, the interpolation results and solid structures are achieved for the exchange of force and displacement through preCICE.

In solving the FSI problem, fluid particles transfer the force to preCICE through interpolation with the critical mesh $\Upsilon$, which ensures that the force generated by fluid particles will act on the patch closest to the critical mesh. Therefore, the interpolation between fluid particles and the critical mesh has the following relationship.

Where $\phi_{i}$ represents the fluid particle $i$, and $\wp_{k}$ represents the pressure exerted by the fluid particle on the critical mesh patch $k$. $\Gamma$ represents the force generated by fluid particles in the $\alpha$ dimension ($\alpha=1,2,3$). $N$ represents the total number of fluid particles in contact with the patch as shown in Fig.10. Obviously, equation (32) shows that the force on the critical mesh can be obtained by summing all the contact particles, and then the resultant force obtained by particle interpolation is transferred to the solid structure through preCICE as the pressure generated by the fluid on the contact surface. In addition, the particle interpolation method can calculate the resultant force by using two methods. The first method is to calculate directly through the particle node using Newton’s second law. First, the acceleration of the particle is obtained through the momentum equation (13), and then the acceleration and particle mass are used to calculate the corresponding force, Finally, the forces of all particles are added and transferred to the critical grid as the force generated by the fluid on the coupling boundary. However, because this method relies on contact particles as the calculation carrier, there is an obvious disadvantage of this method, which is vulnerable to numerical errors caused by particle distribution disturbances. In order to make up for this defect, equation (33) describes the solution method based on particle pressure integration. It uses the basic principle of integration to solve the contact force generated by fluid particles on solid structures. This integral representation method based on continuous function has the advantages of robustness and stability, and it can better maintain numerical stability in the calculation process.

Where $P$ is the pressure calculated from equation of state (20), $\Omega$ is the range of contact particles, and $\mathrm{d}S$ is the area element. Therefore, equation (33) shows that the force on the patch $\wp$ can be obtained by solving the integral of the pressure on the contact surface by its contact particles, which is less affected by the particle distribution. Therefore, this calculation method is more accurate than the method based on a single particle to solve the resultant force for critical grid. In order to realize the interpolation of force, we first define the particle contact method (PCM) as shown in algorithm 2.

Algorithm 2 describes the contact algorithm between particles and patches, which is used in the interpolation process of particles and critical meshes, including force and displacement interpolation. Among them, the line 2 sets the contact threshold $\Theta$, which is generally set to 2$h$. The fourth line represents the process of searching the nearest critical mesh patch for particles close to the coupling surface. Due to the weak dependence, the module can adopt multi-threaded parallel technology to optimize the efficiency of calculation. In line 7, calculate the distance between particles and patches, which is compared with the contact threshold in line 8. When the distance is less than the set threshold, we judge that there is contact between fluid and solid. On the contrary, when the contact distance is greater than the threshold $\Theta$, it is considered that there is no contact relationship. Finally, in line 10 or 13, we save the particles in contact around the solid to array $\chi$, and it is called contact particles, which are provided for subsequent interpolation with the critical mesh. The particle mesh contact process is described in Fig.9 (d) and (e). It should be noted that the contact particles calculate the shortest distance, which ensures that each particle almost establishes a contact relationship with only one patch. This correspondence is a one-to-many mapping relationship, that is, a patch contains multiple contact particles, but a contact particle only corresponds to the target patch. Therefore, based on this contact relationship, the interpolation of particle force can be implemented by the following algorithm 3.

Algorithm 3 describes the simplified process of particle force interpolation, in which line 3 traverses all patches $\Re$ of the critical mesh. Then, according to the contact relation $\chi$ provided by algorithm 2, we calculate the acceleration of particles according to the momentum equation, then calculate the force generated by a single particle on the critical grid according to Newton’s second law (the particle mass is a constant), and finally calculate the resultant force $\wp$ generated by all particles on the critical grid as a whole. Then the resultant force $\wp$ is transferred to the solid structure through preCICE. At the same time, the force interpolation can also be implemented through equation (33). In theory, although it has high accuracy, the disadvantage is that it will be difficult to realize the pressure interpolation on the complex coupling surface, and this process will consume additional computing resources. To sum up, after realizing the critical mesh and particle contact algorithm, the force of particles on the solid mesh can be calculated, and the accuracy of interpolation will also have a certain impact on the accuracy of the final coupling results. Therefore, it is very important to choose an efficient particle mesh interpolation algorithm to obtain high-precision results.

In section V.2, we describe the force interpolation procedure of particles, and then the solid structure calculates the displacement of deformation according to the force exerted by particles on its surface, and the displacement is transmitted back to the critical grid through precCICE. Therefore, it is necessary to interpolate the displacement received by the critical mesh to the contact particles and update the position coordinate information of the particles. The displacement interpolation procedure can be described by algorithm 4.

Algorithm 4 describes the interpolation procedure of displacement between the critical grid and its corresponding contact particles. It includes two parts: 1) the calculation of critical grid displacement and 2) the correction and update of particle position information. It is necessary to calculate the last displacement $\xi_{last}$ of the critical grid and the current displacement $\xi_{current}$ in line 4, and correct the current contact particles according to the relative displacement $\xi_{relative}$ in line 7 to prevent the boundary particles from non-physical penetration of the coupling surface. However, the contact particles with high speed may still penetrate the coupling boundary in the system update stage of SPH (the information of physical quantities such as particle position and speed will be updated). Therefore, it is also necessary to correct the speed, it needs to satisfy the boundary conditions of equation (34).

Where $\bm{\xi}$ and $\bm{v}$ represent the displacement and velocity of the critical mesh respectively, and $\bm{n}$ represents the normal vector at the critical mesh patch. Equation (34) indicates that to prevent fluid particles from penetrating the boundary at the wet interface, the velocity of the contact particles in the normal vector direction of the critical mesh patch needs to be less than or equal to the moving velocity of the patch, the velocity of the contact particles is corrected using equation (34) in line 9. At the same time, in order to improve computing efficiency, the procedure can use multi-threaded acceleration technology in line 4. To sum up, the displacement interpolation needs to deal with the contact relationship between particles and the critical grid, avoid the non-physical penetration of the contact particles to the coupling boundary resulting in calculation errors.

This example gives the numerical simulation results of the classical FSI problem of the dam breaking flow impacting the elastic plate to verify the accuracy of the PMC method, the benchmark was completed by Liao et al.liao2015free . And it is an important example to verify the FSI problem. The setting of the model is shown in Fig.11, in which the bottom of the elastic plate is fixed and made of rubber, and its corresponding parameters are as follows, the density is $\rho=1161.54$ kg/m³, the Young’s modulus is $E=3.5\times 10^{6}$ Pa, and the Poisson’s ratio is $\nu=0.45$, the model is in a gravity field with a vertically downward direction, and the gravity acceleration is $\bm{g}=9.8$ m/s². The size of the model is shown in Table 1.

In this experiment, the width of the elastic plate is 0.2 m, the thickness is 0.004 m, and the height is 0.1 m. The spacing between the rubber plate and the right wall is 0.2 m, the length of the water is 0.2 m, and the height is 0.4 m, the rubber plate was marked with points, where $h_{3}$=0.087 m. The whole device is placed in a water tank with the length-width-height of $0.8\times 0.2\times 0.6$. This example is used to verify the accuracy of the numerical results of the PMC method. The dam breaking flow is solved by the SPH method, and the elastic baffle is solved by the FEM method. The simulation results are shown in Fig.12.

Fig.12 describes the impact process of dam break flow on the baffle, which can better verify the performance of the PMC method. It records the flow field and the deformation process of the baffle within 0.25 to 0.36 seconds. The color of the fluid indicates the flow velocity of the flow field, and the color of the elastic baffle indicates the deformation. First, within 0.25 s, the fluid contacted the elastic plate and began to impact it. When the time reaches 0.27 seconds, the elastic baffle begins to produce concave deformation due to the impact of water flow. Obviously, the fluid particles contacting the elastic baffle do not produce non-physical penetration to the baffle, which is the result of using equation (34) to constrain the boundary conditions. When the time is within 0.32 s, because the bottom of the baffle is fixed, its shape variable is smaller than that of the top. Therefore, the water velocity near the bottom of the baffle will drop sharply and the direction will also change. At the same time, the top of the baffle is movable, and it will produce greater deformation after being impacted by the water flow. And most of the fluid particles will continue to flow forward across the top of the baffle. In the time from 0.32 to 0.36 s, the flow velocity of the flow field has a large loss until it is finally reduced to 0. Finally, the elastic baffle will return to its original position due to the force balance. The contact radius we used in this simulation example is $\Theta=0.01$ m. Different contact radii will affect the simulation results. For example, if the contact radius is too small, even if the boundary conditions are imposed, the fluid particles can still produce non-physical penetration to the elastic baffle, which is caused by the updating of particles in the calculation process of the SPH method. In this way, the faster particles penetrate into the interior of the solid without modifying the coupling boundary conditions, which is illegal. Therefore, in the actual simulation, we need to comprehensively consider the impact of various factors on the simulation results.

In order to further verify the performance of the PMC method, we compared it with the actual experimental results of Liaoliao2015free . As shown in Fig.LABEL:fig131 and Fig.13, the experiment gives the results within 0.25 s to 0.3 s, and the flow procedure at different times is recorded by a high-speed camera. The fluid particles first touched the baffle within 0.25 s, and the structural module has been deformed. With the further increase of the fluid velocity, the structure part has depression, and the simulation results are consistent with the actual experimental results. The structure (elastic plate) is further deformed within 0.28 s after being impacted by water flow. Finally, the structure was deformed greatly within 0.3 s, and the fluid particles passed through its top. However, the experimental results show that some fluid particles are thrown out, and there are subtle differences in the location of local fluid. The reason is that the actual experiment also has an influence on airflow. In addition, as shown in Fig.14, the displacement in the Y-direction is marked (refer to the marker point in Fig.11), it can be seen that this procedure describes the deformation process of the marker point simulated by the PMC method within one second, and it is compared with the actual experimental results of Liao liao2015free . Firstly, the fluid particles come into contact with the rubber baffle within 0.25 s and force it to produce deformation displacement along the Y-direction. Secondly, when the simulation time reaches 0.4 s, the deformation displacement generated by the baffle reaches its maximum. Overall, between 0.3 s and 0.6 s, the deformation displacement amplitude simulated by the PMC method is slightly lower than the experimental results. Therefore, within the allowable error range, the simulation results are consistent with the actual physical results. The original experiment is to verify the multi-physical field coupling problem between water flow, air, and structure. This paper only considers the FSI problem between fluid and structure. The scale of simulation is relatively small, and the accuracy of the SPH solver itself is limited. Therefore, comprehensively excluding the influence of the above factors, it can be concluded that the PMC method proposed in this paper is effective in solving the FSI problem.

In order to further verify the performance of the PMC method, we give the numerical simulation results of water flow through the elastic gate. This example is also a representative benchmark for the FSI problem. Its specific experiment was completed by Antoci et al.antoci2007numerical . The experimental setup is shown in Fig.15. The fixed plate and the rubber gate fix the water on the right side of the water tank. The top of the rubber plate is fixed and its bottom is free to move. The whole device is placed in a water tank of size $4\times 12\times 8$ m. In this experiment, the FEM method is used to disperse the rubber gate, and the SPH method is used to disperse the water flow. The fluid adapter module interpolates data (forces and displacements) through the critical grid and interacts with preCICE.

The experimental dimensions used in this paper are shown in Table 2, where the young’s modulus of the elastic plate $E=1.2\times 10^{7}$ Pa. It determines the deformation of the elastic door under the pressure of water flow. The density of the elastic gate made of rubber is $\rho=1100$ kg/m³, and the Poisson’s ratio is $\nu=0.45$. The whole device has a vertical downward gravity field. The value of gravitational acceleration in this paper is $\bm{g}=9.8$ m/s², the actual results of this experiment are shown in Fig.16.

The experiment of Antoci et al.antoci2007numerical showed that the whole fluid is in the initial state at t=0 s, the rubber gate is fixed by the baffle, and the water generates pressure on the surrounding wall. With the removal of the baffle, the bottom of the rubber gate began to deform greatly under the pressure of water at t=0.04 s. As the water flows out of the bottom of the rubber gate, the rubber gate is further deformed by the pressure of the water flow, and the state is maintained until the water flows out from the right side. In this procedure, the elastic gate has been under the pressure of water from the right side, so the bottom of the elastic gate deforms to the left.

Fig.17 is the simulation result using the PMC method with a time span from 0.3 s to 0.8 s. Obviously, rubber gate begin to be subjected to the pressure of water flow and has a slight distortion. The bottom of rubber gate is movable, so it will be subjected to the pressure of water flow and show greater distortion at the bottom. As the pressure on the bottom of the right flow increases, the water flows out from its bottom and reaches a relatively stable state. As the fluid particles on the right decrease, the pressure on the bottom of the rubber gate decreases until a new force balance is reached and the elastic plate finally returns to its original position. The simulation results are in good agreement with the experiments of Antoci et al. Because the precision of the SPH solver is limited, the size used in this paper is slightly larger than that of the actual experiments. Therefore, the results of each phase are lagging behind those of the actual experiment. However, the simulation results show that the flow of water and the deformation procedure of the elastic plate is consistent with the experimental results.

Fluid-structure interaction (FSI) problem is a classical problem in multi-physical field coupling, which has been successfully applied in many fields liu2019smoothed ; zhang2021multi ; gong2016two . The fluid-structure interaction phenomenon is also widespread in practice. The study of fluid-structure interaction has very important applications in aerospace takizawa2020computational ; vijayanandh2019design , water conservancy engineering wang2017numerical ; guo2020fluid ; guo2020fluidpp , safety protection paik2018one , petrochemical lee2018fluid , oceanographic engineering wilkes2016numerical , and biological engineering samaee2017coupling . For example, the earthquake disaster area is prone to cause debris flow disasters, which can destroy houses, bridges, and other public facilities and pose a major threat to human life and property. In order to study the impact force of debris flow on the target buildings and establish stable safety protection measures, Dai et al. dai2017sph established the FSI model of simulated debris flow, in which the debris flow material is simulated as a viscous fluid and solved by SPH, while the check dam is simulated as elastic solid. They solved the governing equations of the two phases and calculated the interaction between them. Finally, they verified the model by simulating the sand flow model test and successfully applied its propagation evolution process to practical problem research. In addition, Zhang et al.zhang2017smoothed summarized the application of the FSI method in oceanographic engineering in detail, including the capture of phenomena such as water waves, impact, and splash jets. The formation of these phenomena is mainly related to free surface flow.

For most computational fluid dynamics (CFD) solvers based on the Euler framework principle, such as the finite element volume method and the finite difference method, it is a difficult and challenging task to simulate these complex free surface flows. As a Lagrangian and meshless method, smoothed particle hydrodynamics (SPH) has the ability to track different complex boundaries and provides direct support to meet different boundary conditions. Due to the robustness and high accuracy of the SPH method in simulating complex fluid dynamics problems such as free surface boundary, multiphase interface, or material discontinuity, it is very suitable to use the SPH method in fluid simulation of FSI problems. Therefore, another example of dam break impact is given to verify the free surface flow in FSI.

As shown in Fig.18, it describes the process of free surface flow and solid deformation in the FSI problem. The fluid module is simulated by SPH method and the solid structure is simulated by FEM method. The dam break water flows out from the left side, the top of the square baffle is fixed, and the bottom is free to move. The physical time corresponding to the evolution process from Fig.18 (a) to Fig.18 (f) is from 0.5s to 1s, with a time interval of 0.1s. First, from this evolution procedure, it can be seen that the fluid contacts the solid structure and compresses it, resulting in the deformation of the solid. The water flow at the bottom flows out from the bottom of the baffle. The water flow at the top is blocked by the baffle and generates broken waves. The flow field at the top splashes from above the baffle and flows forward across the baffle. This procedure can capture the actual physical phenomena such as broken waves, the impact of water flow on the baffle and the deformation of the baffle. This example shows that the PMC method proposed in this paper can also effectively simulate a series of physical procedures in FSI problems involving free surface flow.