Abstract
This work presents a three-dimensional two-way coupled method to simulate moving solids in viscous free-surface flows. The fluid flows are solved by weakly compressible smoothed particle hydrodynamics (SPH) and the displacement and rotation of the solids are calculated using the multisphere discrete element method (DEM) allowing for the contact mechanics theories to be used in arbitrarily shaped solids. The fluid and the solid phases are coupled through Newton’s third law of motion. The proposed method does not require a computational mesh, nor does it rely on empirical models to couple the fluid and solid phases. To verify the numerical model, the floating and sinking processes of a rectangular block in a water tank are simulated, and the numerical results are compared with experimental results reported in published literatures. The results indicate that the method presented in this paper is accurate and is capable of modelling fluid-solid interactions with a free-surface.
1. Introduction
Fluid-solid interaction in free-surface flows has gained a lot of interests in recent years due to its involvement in a wide variety of industrial processes and engineering disciplines such as hydraulic, ocean, and environmental engineering. Accurate numerical simulations are helpful in obtaining detailed information of the interaction process and in improving the engineering design while avoiding tedious and time-consuming physical experiments. There are generally two types of coupling model [1] applied to calculate the interaction between fluid and solids: one-way coupling and two-way coupling. Substantial efforts have been made in the past to develop numerical methods for two-way coupling of fluid and solid. For example, Singh et al. [2] used the distributed Lagrange multiplier method for particulate flows with collisions. Takizawa et al. [3] modelled fluid-solid interactions with free-surface flows using CIP method. Latham et al. [4] simulated discrete solids in motion within fluids by coupling of DEM and CFD based on the arbitrary Lagrangian-Eulerian (ALE) frame. Wu et al. [5] presented a two-way coupled Moving Solid Algorithm (MSA) to calculate the motion of the solids in free-surface flows. In contrast to the above-mentioned mesh-based techniques, meshless methods do not require an interface capturing scheme or the moving mesh technology, which is a clear advantage when the problem involves breaking waves and moving boundary. Of the meshless techniques now available, smoothed particle hydrodynamics (SPH) is proving to be popular and robust. It was originally developed for astrophysical simulations [6, 7] and has been extended to simulate free-surface flows by Monaghan [8]. Instead of using a mesh, the SPH method uses a set of interpolation nodes placed arbitrarily within the fluid. This gives several advantages in comparison to mesh-based methods when simulating nonlinear flow phenomena. More complete reviews on SPH can be found at [9, 10]. Attempts have been made in most recently to develop SPH based method for two-way coupling of solid and fluid in a free-surface flow. Hashemi et al. [11] modelled free floating bodies using a ghost particle method which is usually used for fairly regular boundaries. Liu et al. [12] simulated fluid-structure coupling problems, especially for moving structures, using incompressible smoothed particle hydrodynamics (ISPH) model. Ren et al. [13] developed a two-dimensional SPH-DEM model to simulate the wave-structure interaction by describing the movement of the solids based on multisphere DEM introduced by Favier et al. [14]. Canelas et al. [15] coupled SPH and DEM to describe the motion of arbitrarily shaped solids in viscous fluids based on the concept of DCDEM introduced by Cummins and Cleary [16]. Yang et al. [17] developed an improved SPH-EBG coupling method for modelling the interaction of fluid with free-surface and flexible structure with free and fixed ends.
This work presents a 3D coupled DEM-WCSPH method to simulate moving solids in viscous free-surface flows. The fluid flow was solved by a weakly compressible smoothed particle hydrodynamics (WCSPH) which has an advantage of avoiding solution of the pressure Poisson equation. The motion of the solids was determined by the discrete element method and does not need to be prescribed beforehand. The DEM is also a Lagrangian method which was developed by Cundall and Strack [18] to describe granular material and nowadays is widely used in particulate flows. The DEM model in our study has been implemented using a multisphere approach introduced by Favier et al. [14] to enable solids with different complex shapes to be modelled. The accuracy of the simulation results was verified by comparison with published laboratory experiments of falling and floating solids in water.
2. Methodology
2.1. SPH for Fluid Phase
For the fluid phase, the following continuity and momentum equations are employed and given bywhere is the time, is the velocity, is the pressure, is gravity acceleration, is the kinetic viscosity, and is the density.
The key idea of SPH is that the concerning fluid domain is discretized by use of the moving particles which carry all the properties, such as mass, density, velocity, and pressure. And any quantity of a particle can be approximated by summing over the relevant quantities of its neighboring particles within its support domain [9, 10]:where , , and are, respectively, the mass, density, and position of a particle, is the smoothing length, being the distance between particles and , and is the smoothing kernel. In general, the higher the order of the kernels, the greater the accuracy of the SPH scheme. The Wendland kernel has some advantages over the Gaussian and the Cubic Spline [19]. So, it is used in this study and written by [20]where is in 2D and in 3D, .
Following Morris et al. [21] and Lo and Shao [22], (1) can be written in SPH discretization aswhere , and denotes the gradient taken with respect to the coordinates of particle .
As mentioned before, the fluid phase is treated as weakly compressible where the density in the fluid is allowed to vary slightly. The closure problem is solved by using an equation of state to relate changes in pressure to changes in density. Following Monaghan [8] and Batchelor [23], the relationship between pressure and density for a particle is given bywhere constant , is the reference density, is the sound speed at the reference density, and for a fluid like water. The speed of sound is generally chosen as 10 times the maximum velocity in the fluid to ensure the fluctuations in density are less than 1%.
When solving any boundary value problem, the proper implementation of boundary conditions is crucial. Due to its Lagrangian description and to kernel truncation, it is challenging to implement solid boundaries in SPH. In this work, the so-called dynamic boundary conditions [24] are used for both the moving and fixed boundary. It consists of creating boundary particles that satisfy the same equations of continuity (4) and state (6) as the fluid particles, but their positions are not updated using the momentum equation (5) and they remain fixed in position (fixed boundaries) or move according to some externally imposed function (moving objects).
Finally, particle positions updated every time step using the XSPH variant [10]:where and is a constant, whose values range between zero and unity, and is often used. This method is a correction for the velocity of particle . This correction lets particles to be more organized and, for high fluid velocities, helps to avoid particle penetration.
2.2. DEM for Solid Phase
DEM coupled with SPH has been used to numerically model the moving solids in viscous free-surface flows. The motion of each solid is tracked based on Newton’s laws of motion as follows:where solid possesses mass , velocity , and angular velocity . is gravity acceleration. Force is due to fluid-solid interaction and force represents any solid contact that might occur. Integrating (8) into time advances the linear motion of the solid, whereas (9) accounts for the rotational motion.
In order to calculate the fluid-solid interaction force , the boundary of a moving solid is represented by groups of SPH particles and these moving boundary particles have an interparticle spacing that is equal to half of the fluid particle spacing to prevent the fluid particles from penetrating the moving solid boundary. The fluid-solid interaction is ensured through the interactive force balance condition satisfying Newton’s third law of motion. Following Ren et al. [13], the fluid-solid interaction force can be determined aswhere the inner summation means the total force on moving boundary particle belonging to solid due to the neighborhood fluid particle .
Following Latham et al. [4], a multisphere approach originally proposed by Favier et al. [14] was used for modelling the complex-shaped multibody dynamics in which the surface of each solid is represented by a cluster of small spheres whose diameter is equivalent to the spacing of the moving boundary particles. The solids are allowed to interact via contact forces when the surface spheres of different solids overlap. No relative movement between spheres of the same body is allowed. Based on the multisphere approach, the resultant contact forces and torque acting upon solid are evaluated aswhere is the radius of surface sphere belonging to and represents the normal unit vector pointing from the center of sphere belonging to to the center of sphere belonging to . and are the normal and tangential forces, respectively, on surface sphere of solid due to surface sphere of solid . Based on the model by Cundall and Strack [18], the normal contact forces are calculated bywhere is the normal spring stiffness, the overlap , where and are the surface sphere radii, and and are the surface sphere centers. is the normal relative velocity, and is the normal damping coefficient given bywhere is the restitution coefficient. The effective mass , where and are the mass of surface sphere.
The tangential component of the contact force is calculated by a Coulomb friction law using a coefficient of friction . It can be expressed aswhere , , and are the tangential spring stiffness, tangential overlap, and tangential damping coefficient, respectively. The tangential relative velocity and the tangential unit vector . The tangential damping coefficient given bywhere and . The value of causes the normal and tangential springs to react on similar time scales [24].
3. Numerical Examples and Discussion
In this section, numerical examples are presented to validate both solid-solid and fluid-solid interactions. For the cases involving water, the initial velocity of the fluid particles is considered zero. The fluid particles are assigned an initial density based on hydrostatic pressure. Therefore, the density of particle (located at depth ) is calculated taking in account of the water column height as where is the initial water depth at position and is the vertical distance from the bottom. The pressure is calculated using equation of state (6) following the density value. In all the simulations presented herein, the modelled fluid is water with reference density kg/m3 and kinetic viscosity m2/s. For all simulations, the smoothing length is defined as , where represent the initial particle spacing.
3.1. Solid Dropping onto a Flat Floor
In this simulation, a sphere is dropped under the influence of gravity onto a flat floor. The purpose of this simulation is to verify that the momentum exchange between a solid and a flat boundary agrees with theory based on the principle of impulse and also to see if the motion of the solid subject to a gravitational force is accurate. The problem contains one sphere whose bottom is 0.5 m above a flat floor. The solid sphere has a density of 2.6 × 103 kg/m3, a diameter of 0.1 m, and a normal spring stiffness of 1 × 105 N/m. The collision is assumed to be frictionless and purely elastic. So, coefficient of friction and the restitution coefficient . A time step of 1.0 × 10−4 s and an initial fluid particle spacing m are used in the simulation. The sphere is dropped and falls under the force of gravity until it contacts the level surface. The translational motion of the sphere can be described in three stages: free fall, contact, and rebound. Figure 1 shows snapshots of the falling sphere at different time. Figure 2 compares the vertical velocity and the height of the solid as a function of time. The results agree very well with the analytical solution given in the work of Chen et al. [25].
(a)
(b)
(c)
(d)
3.2. Falling Solid in Water
The aim of this simulation is to evaluate the fluid-solid coupling of a solid in a free-surface flow. The simulation involves a cube falling into calm water in a tank. The trajectory of the cube is compared to laboratory experiments proposed by Wu et al. [5]. In their experiments, a cube was released by a clamp at the water surface and the trajectory of the cube was recorded by a high speed camera. In order to reduce the water splash caused by the cube falling into water, the bottom of the cube was immersed in the water before the cube was released. The clamp was controlled by an electromagnetic device so the cube could be released instantly (see Figure 3).
The problem domain contains a water tank of 150 mm × 140 mm × 140 mm, and the water depth in the tank was 131 mm. The cube with a side length of 20 mm was assumed be a rigid as in the experiment. The density of the cube was 2120 kg/m3, the normal spring stiffness was 1 × 105 N/m, the coefficient of friction , and the restitution coefficient . A time step of 1.0 × 10−4 s and an initial fluid particle spacing mm are used in the simulation. Figure 4 compares the vertical displacements of the cube to the experimental results in Wu et al. [5]. The simulation results are in good agreement compared with the measured trajectory in the vertical direction. Figure 5 presents some snapshots at different instants of falling cube evolution with the velocity magnitude () on the central cutting section. The color bar is common to all snapshots and velocities are in meters per second.
(a) sec
(b) sec
(c) sec
(d) sec
3.3. Floating Solid in Water
In this simulation, a wooden block rising from the bottom of a water tank was considered. Figure 6 shows the initial configuration of the problem. The numerical results are also compared with the experimental data provided by Wu et al. [5]. The size of the water tank was 150 mm × 140 mm × 140 mm, and the water depth in the tank was 52 mm. A rectangular wooden block (width = 48 mm, length = 49 mm, height = 24 mm, and density 800.52 kg/m3) was released from the bottom of the water tank. In our simulation, the normal spring stiffness was 1 × 105 N/m, the coefficient of friction , and the restitution coefficient . A time step of 1.0 × 10−4 s and an initial fluid particle spacing mm are used in the simulation. Figure 7 shows the comparison of measured and predicted vertical displacements of the block. It shows a good agreement between the simulated results and experimental data by Wu et al. [5]. Figure 8 shows the snapshots of the configuration of the floating solid at different time and the calculated free-surface profiles are very similar to that of the laboratory results.
(a) sec
(b) sec
(c) sec
4. Conclusions
A three-dimensional numerical model has been developed to simulate moving solids interacting with a free-surface flow by coupling the smoothed particle hydrodynamics and the discrete element method. The motions of the solids do not need to be prescribed beforehand and the mesh-free method is used for the computational fluid dynamics. The DEM model has been implemented using a multisphere approach to enable solids with arbitrarily complex shapes to be modelled. The comparison of results from the coupled approach with known analytical solutions and experimental data indicates that the proposed methods could be an efficient way of simulating complex moving solids in viscous free-surface flows. Our numerical examples focus mainly on fluid-solid interaction; however, multibody systems with solid-solid contacts are also within the capabilities of the model. Further investigations might be needed to consider the solid deformable and turbulent effects.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
The authors are grateful for funding from the National Natural Science Foundation of China (Grant nos. 11172321 and 51509248), the Open Fund of the State Key Laboratory of Hydroscience and Engineering of China, Tsinghua University (Grant no. sklhse-2015-C-03), the Scientific Research and Experiment of Regulation Engineering for the Songhua River Mainstream in Heilongjiang Province, China (Grant no. SGZL/KY-12), and the Open Fund of the State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research (Grant no. IWHR-SKL-201505).