A Novel Approach for the Numerical Simulation of Fluid-Structure Interaction
Problems in the Presence of Debris

A novel algorithm is proposed for the simulation of fluid-structure interaction problems. In particular, much attention is paid to natural phenomena such as debris flow. The fluid part (debris flow fluid) is simulated in the framework of the smoothed particle hydrodynamics (SPH) approach, while the solid part (downstream obstacles) is treated using the finite element method (FEM). Fluid-structure coupling is implemented through dynamic boundary conditions. In particular, the software “TensorFlow” and an algorithm based on Python are combined to conduct the required calculations. The simulation results show that the dynamics of viscous and non-viscous debris flows can be extremely different when there are obstacles in the downstream direction. The implemented SPH-FEM coupling method can simulate the fluid-structure coupling problem with a reasonable approximation.


Introduction
Fluid-solid coupling is universal, such as water dripping through rocks and surging against the bank. The dynamic characteristics of debris flow are sophisticated and easily affected by many factors, such as particle gradation, viscosity content, and bulk density; besides, these characteristics often involve complex and large deformations. In the process of movement, the strain rate is far more than 100%. The traditional finite element method (FEM) finite difference method, and other mesh-based calculation methods are easy to cause mesh distortion when solving the problem of large deformation. At present, some scholars have developed the technology of re-meshing in large deformation area and applied the finite element theory to the calculation of significant deformation problems. However, this not only increases the calculation time but also makes the calculation result easy to produce significant errors. With the development of technology, computer software and hardware have reached the needs of the simulation. For fluid simulation, there are two typical methods: The Eulerian grid method and the LaGrange particle method. The LaGrange particle method is widely used because it is more suitable for dynamic fluid simulation. Notably, the smoothed particle hydrodynamics (SPH) method is the most popular.
After its proposal, SPH was quickly applied to fluid simulation problems. Solid simulation is simple, and the FEM is widely used. The first application is to divide an entire stress surface into triangular elements for calculation, and finally, obtain the plane stress successfully. The FEM can divide an entire area arbitrarily, and boundary conditions do not limit it. It is convenient and efficient to calculate. Therefore, the advantage of the FEM is distinct when simulating the fluid-structure coupling model. These advantages make the FEM popular.
The debris flow fluid is a non-Newtonian fluid. Given the characteristics of the non-Newtonian fluid, the SPH method is used to simulate the debris flow. Starting from 2D non-viscous dam break fluid, the simulation is gradually extended to 3D viscous dam break fluid, and the governing equation of the fluid motion is derived. Based on the Python language, the 3D motion program of debris flow fluid is written. The calculation is carried out through TensorFlow, and the 3D motion control equation of debris flow fluid is obtained. FEM is used to simulate the deformable solid, and its finite element analysis is carried out. Finally, the fluid-structure coupling boundary is analyzed. The innovation of this simulation is to combine the problem of fluid-structure coupling with TensorFlow to simulate more quickly and accurately, which is of positive significance to the simulation of the process of debris flow.

Experiment and Method
Fluid-structure coupling, a common natural phenomenon, has been studied and explored by scholars for many years. All kinds of methods emerge in endlessly. Here, the reason that the problem of fluid-structure coupling has been studied is that this kind of natural disaster, i.e., debris flow, has threatened the lives of human beings for years. Thus, if humans can master the movement state of this kind of natural disaster through simulations, the prediction or early-warning of this disaster will be achieved, and humans may minimize the damages brought by it.
The initial ideas and methods of simulation experiments are as follows. (1) The fluid in the fluidstructure coupling problem is simulated by SPH, which is regarded as a single fluid particle. Each particle is analyzed for its motion state, and the results are extended to the entire fluid. (2) FEM simulates the solid in the fluid-solid coupling problem. The entire solid is divided into one small unit for simulation, which better adapts to the solid boundary with different shapes. (3) DBC simulates the fluid-solid boundary. The boundary particles are generated automatically from the solid wall boundary; therefore, it will be harder for the fluid particles to penetrate the solid boundary. Besides, in this case, these particles can better adapt to the solid boundary with different shapes, thereby improving the calculation efficiency. (4) The leapfrog scheme is used to solve the Navier-Stokes equation, update the motion state of particles in real-time, and control the time step of the coupling algorithm. The solutions of fluid and solid are limited to the same time step. (5) Python is used for programming. Through the interface, TensorFlow is used for calculation, which dramatically shortens the calculation time and better satisfies the real-time performance of the simulation process. (6) A complete fluid-structure coupling simulation program is formed. Given obstacles in the downstream, the non-viscous fluid dam break, viscous fluid dam break, and debris flow dam break are simulated, and the simulation results are consistent with the reality.

SPH
SPH regards continuous, flowing, and changing fluid as discrete and interacting particles. All physical quantities are attached to each particle separately, such as mass and velocity. By studying the motion state of each particle, the motion of the entire fluid is obtained. According to the obtained fluid situation, the entire fluid is simulated, as shown in Fig. 1.
The key of the algorithm is to solve the motion of each particle and analyze the particle by Newton's second law, F = ma. According to hydrodynamics, the force acting on particles consists of three parts [1], as shown in Eq. (1): In Eq. (1), F external represents gravity, F external ¼ qg. Since the fluid unit density determines the fluid mass, the mass is replaced by q. F pressure represents pressure, the force generated by the pressure difference inside the fluid; numerically, it is equal to the gradient of the pressure field, and the direction is from the area with high pressure to the area with low pressure, i.e., F pressure ¼ Àrp. F viscosity is the viscous force, which is caused by the velocity difference between particles [2]. The magnitude of this force is related to the viscosity coefficient μ of the fluid and the velocity difference, i.e., F viscosity ¼ lr 2ũ . The smooth length is expressed in h. Therefore, the force of particles is shown in Eq. (2): SPH algorithm also involves the concept of "smooth kernel" [3]. The particle's attributes "spread" around, and the effect diminishes as the distance increases. The function that changes with the distance between particles is "smooth kernel function," and the maximum influence radius is "smooth kernel radius." Here, the fluid is treated as a single particle. However, the fluid fills the whole space. The operation of each position in the fluid is related to every particle around [4]. Therefore, the parameters or values at each position are accumulated by the surrounding particles, as shown in Fig. 2.
An attribute A on point r is shown in Eq. (3): The density q ðr i Þ of point r is shown in Eq. (4): The pressure p of point r is shown in Eq. (5): The pressure acting on point r is shown in Eq. (6): The viscous force acting on point r is shown in Eq. (7): The equation of particle motion is shown in Eq. (8):

FEM
FEM is to divide a complete solution area into finite triangles, which are connected but not coincident [5]. The polynomial difference is used to solve each triangle region. By combining the solutions of all triangles, the approximate solution of the region can be obtained [6]. FEM can adapt to various shapes of the solution area, with fast speed and high accuracy.
Based on the dynamic equation of continuous medium, the equation of the FEM is similar to its definition. According to the weighted residual method, the velocity is taken as the weight function. The momentum conservation equation is shown in Eq. (9): Z V dvað @r ba @x b þ qb a À q€ u a ÞdV ¼ 0 In Eq. (9), @r ba is Cauchy stress, and b a is physical force. The following aspects should be noted for highly precise FEM solution: 1. Unit order. The first order, second-order, and third-order can summarize the current unit order [7]. A higher-order element is above the second-order. Because of the existence of curves and surfaces, the higher-order elements can better simulate the approximation. However, the higher the order is, the more grids need to be processed. In addition, the difficulty of calculation, the demand for hardware and software, and the calculation time will increase. 2. Grid density. The more polygons are divided, the closer the calculated area is to the area of the circle.
Here, the more grids are divided, the closer they are to the real situation [8]. However, the number of grids also affects the calculation cost. Therefore, to better adapt to various situations, different segmentation methods are adopted in different parts of the same model. Such processing methods can effectively avoid the situation of too high calculation cost [9]. 3. The quality of the grid. The quality of the grid is closely related to the size and shape of the mesh [10].
However, the mesh quality is usually challenging to control. Therefore, in the actual design of the finite element, it is assumed that there are constraints on the grid to ensure the grid quality [11].
Therefore, the accuracy of the finite element solution is closely related to the grid division method; thus, scholars have developed three high-quality division methods: Grid method, Delaunay triangulation method, and mapping method.
The quality of the generated grid is poor [12]. Delaunay triangulation is suitable for the simulation of 2D problems, as shown in Fig. 3. The mapping method is suitable for the simulation of 3D problems, as shown in Fig. 4.

Neighborhood Particle Search Method
At present, there are three commonly used methods for neighborhood particle search: Global pairing method, tree search method, and linked list search [13]. The basic principle of the global pairing method is to delimit a fixed radius area around a particle, and the particles in this area are the neighborhood particles of the chosen particle. Assuming that there are N particles in the defined area, the time complexity of the global pairing method is O (N2) [14]. The more particles are, the longer the calculation time of the global pairing method is, which is not suitable for real-time model simulation. Therefore, the global pairing method is generally not used for the model with many particles.
The tree search method needs to divide the whole model into several small blocks to guide each block to contain one particle; thus, it needs too much storage space and requires high software and hardware.
The method used in this study is a linked list search, which divides an entire area into several cells of the same size, and numbers each cell. Particles in each cell correspond to the cell number, and particles adjacent to the cell are searched. Currently, the search for neighborhood particles is much more convenient.

Boundary Treatment Method
Except for the simulations of solid and fluid, a severe problem in simulating the fluid-solid coupling problem is the boundary treatment. Since SPH is based on the LaGrange particle method, the motion of some particles in the fluid follows their LaGrange velocity. Hence, the fluid boundary can be captured automatically; however, on the interface where the fluid contacts with the solid, the fluid particles must not penetrate the solid. There are two solutions to these problems: the virtual particle method and the penalty function method. The virtual particle method is suitable for the situation with a simple boundary but high precision demand, and the penalty function method is suitable for the situation with complex boundaries [15]. Because it is difficult to simulate the process of debris flow movement, the existing boundary treatment methods are insufficient. Therefore, the original method is improved to achieve the expected results. The method used here is a dynamic boundary condition (DBC). More than two layers of solid wall particles are set, which are staggered with each other and utterly consistent with the motion state of the fluid, as shown in Fig. 5.
Boundary particles can be generated automatically by a solid wall boundary.

Overall Algorithm Flow
The procedure of SPH-FEM coupling algorithm is as follows: (1) The minimum element size and the minimum time step of the finite element are calculated.
(2) Combined with the time step of the SPH algorithm, the time step of the SPH-FEM coupling algorithm is calculated.
(3) The DBC method is used to process boundaries.  (4) The viscosity, density, and pressure of fluid particles are calculated.
(5) The leapfrog scheme is used to solve Navier-Stokes equations and update the motion state of particles in real-time.
At present, most of the incompressible viscous fluid simulation uses Navier-Stokes equations as the governing equations of fluid motion. Its vector form is as follows: ru ¼ 0 In Eqs. (10) and (11), u is the velocity of the fluid, q is the density of the fluid, p is the internal pressure of the fluid, f is the external force on the fluid, r is the gradient operator, v is the elastic viscosity coefficient, and r 2 is the Laplace operator. Eq. (10) shows the momentum conservation. The four terms on the right represent advection term, pressure term, viscous force, and external force in the process of fluid motion. Eq. (11) shows the fluid mass conservation, which is also known as the incompressible condition of the fluid.  It is shown in Fig. 6:

Hardware and Software Environment
The Python language is used to write programs. Because of the complexity of simulation and many fluid particles, the calculation process costs too much time. The calculation efficiency of standard software cannot meet the real-time research very well; so, the high-level interface provided by TensorFlow is used to realize the calculation.
TensorFlow, as an open-source, high-performance computing library, supports IOS, Windows, MacOS, Linux operating systems [16]. Also, it can realize large-scale local and distributed parallel computing on Fluid particle The boundary particle The solid wall boundary Figure 5: Dynamic boundary conditions handling solid wall boundaries FDMP, 2020, vol.16, no.5 GPU, CPU, TPU, FPGA, and other hardware architectures [17]. Hence, the same code can achieve highperformance computing on different hardware platforms without any modification, reducing the cost of software and time. The frame composition is shown in Fig. 7. The hardware environment required for the experiment is shown in Tab An obstacle with length, width, and height of 0.15 m, 0.5 m, and 0.4 m, respectively, and a water column with length, width, and height of 1 m, 1 m, and 0.5 m, respectively, are set. Obstacles and water columns are installed in a water tank, which is 3 m long, 1 m wide, and 1 m high. The initial distance between the obstacle and the water column is 1 m.

Start
In order to prevent fluid particles from penetrating, three-layer boundary particles are set. The specific parameters are shown in Tab. 2: The movement process of dam break fluid in the flume is as follows:   The whole process is shown in Fig. 8:

Simulation Results of the Viscous Fluid
The laminar viscosity model is transformed into the SPH form. The same model is simulated to compare the difference between viscous fluid and non-viscous fluid. The equivalent viscous fluid replaces only the non-viscous fluid, and the calculation parameters are shown in Tab. 3: For a viscous fluid, unlike a non-viscous fluid, particles are affected by viscous force. At the initial stage, the viscous fluid is contractive, the velocity of the viscous fluid is lower than that of the non-viscous fluid, and the time of the viscous fluid reaching the edge of the obstacle is shorter than that of the non-viscous fluid. After collisions, no "back splashing" phenomenon like a non-viscous fluid occurs, but a slow climbing  trend appears. Then, under the influence of gravity, it will fall and tend to be gentle, flowing forward along both sides of the obstacle, as shown in Fig. 9: It takes about 3.0 s from the beginning to the surface of the whole dam break fluid to become calm. It takes about 0.43 s for viscous fluid to contact the left edge of the obstacle. Moreover, due to the effect of viscous force, the dispersion between particles of viscous fluid is far less than that of non-viscous particles, which is also the reason that the viscous fluid does not collide with the obstacles to produce intense water bloom.

Simulation Results of Viscous Debris Flow
Debris flow is a kind of fluid formed by many poorly sorted solid particles and liquid water under the action of gravity. As a kind of natural disaster, the initial gravitational potential energy of debris flow is abundant, the movement speed is fast, and it is difficult to predict it. Therefore, the study of debris flow movement has practical significance, which can help to prevent or reasonably control debris flow.
According to the nature of the fluid, debris flow can be divided into three types: Strong viscosity (bulk density is more than 2 g/cm −3 ), sub-viscosity (bulk density is more than 1.85 g/cm −3 ) and rareness (bulk density is between 1.46-1.85 g/cm −3 ). In addition to debris flow, there are mudflow and water stone flow. Here, the viscous debris flow formed by the mixture of sediment and rainwater is studied.
Based on viscous fluid, shear stress is added into the fluid to simulate debris flow fluid. Again, the same model as before is used for comparison. The fluid is replaced by debris flow fluid. Besides, to simulate the situation of debris flow accompanied by rain, water is used to simulate rainwater during the simulation process.
According to the previous model, the debris flow model is analyzed. Compared with viscous fluid, the velocity of viscous debris flow is the same as that of viscous fluid. However, when the viscous debris flow passes over the obstacles, "overflow" will appear, that is, after the viscous debris flow collides with the obstacles, the viscous debris flow will not flow back into the fluid, but will flow forward through the obstacles, as shown in Fig. 10: Because the debris flow fluid is affected by multiple factors, such as sediment content and particle size, the debris flow simulation conducted in the laboratory may have some errors with the debris flow generated in different locations and different states. Here, only the continuous and uniform debris flow is discussed.
Under the same working condition and small model size, the three kinds of fluids present different motion patterns when they encounter obstacles. It is proved that the addition of viscous force and shear stress is hugely beneficial better to simulate the movement state of viscous debris flow.

Conclusion
The motion state of dam break fluid when encountering obstacles is mainly simulated. The cohesionless dam break fluid contacts the left boundary of the obstacle. Then, the dam break fluid hits the barrier, and its velocity decreases after hitting the barrier. Part of the fluid climbs up the obstacle. Part of the fluid continues to flow forward from both sides of the obstacle and forms an "open space" on the right side of the obstacle. The fluid that climbs the barrier falls back into the overall fluid. The entire fluid area is steady and covers the entire water tank. For a viscous fluid, unlike a non-viscous fluid, particles are affected by viscous force. At the initial stage, the viscous fluid is in the form of contraction, and the velocity of the viscous fluid is lower than that of the non-viscous fluid. The time for the viscous fluid to reach the edge of obstacles is significantly shorter than that for non-viscous fluid. After collisions, no "back splashing" phenomenon like a non-viscous fluid occurs, but a slow climbing trend appears. Then, it falls under the influence of gravity, tends to be gentle, and flows forward along both sides of the obstacle. Compared with viscous fluid, the velocity of viscous debris flow is the same as that of viscous fluid. However, when the viscous debris flow passes over the obstacles, it will appear "overflow," that is, after the viscous debris flow collides with the obstacles, the viscous debris flow will not flow back into the fluid but will flow forward through the obstacles. The indepth study carried out by scholars worldwide on the dynamics of natural disasters such as flood, dam break, landslide, and debris flow by using the SPH method is consistent at home and abroad.
SPH-FEM coupling method is still in the development stage. Thus, there are still some imperfections in this method. It is hoped that scholars can continue to study and apply the SPH-FEM coupling method to real life as soon as possible. The following aspects need to be further studied. (1) According to the movement model of different obstacles, such as smooth stones, trees, and grass, the impact of different obstacles on the movement speed of debris flow can be analyzed; (2) Because SPH is a meshless method, which is very similar to DEM, the two can be combined in the later study. SPH is used in the fluid part, and DEM Funding Statement: The author(s) received no specific funding for this study.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.