Abstract

In this paper, we present a grid-free modelling based on the finite particle method for the numerical simulation of incompressible viscous flows. Numerical methods of interest are meshless Lagrangian finite point scheme by the application of the projection method for the incompressibility of the Navier–Stokes flow equations. The moving least squares method is introduced for approximating spatial derivatives in a meshless context. The pressure Poisson equation with Neumann boundary condition is solved by the finite particle method in which the fluid domain is discretized by a finite number of particles. Also, a continuous particle management has to be done to prevent particles from moving into configurations problematic for a numerical approximation. With the proposed finite particle technique, problems associated with the viscous free surface flow which contains the study on the liquid sloshing in tanks with low volumetric fluid type, solitary waves movement, and interaction with a vertical wall in numerical flume as well as the vortex patterns of the ship rolling damping are circumvented. These numerical models are investigated to validate the presented grid-free methodology. The results have revealed the efficiency and stability of the finite particle method which could be well handled with the incompressible viscous flow problems.

1. Introduction

In the last several decades, various meshless techniques [16] have been developed as alternatives to traditional grid-based methods such as finite volume method [7, 8] and finite element method [9, 10] as well as boundary element method [11, 12], which have attracted plenty of research fellow for the potential interest and requirement. Compared to the traditional mesh generation method which is of consumption and difficulty and even hardly solving the thorny problem with the large deformation viscous flow or within complex geometry boundary, the grid-free modelling reflects superior ability to complete simulation and adapt it locally by constructing space arbitrary point distribution without using any mesh to connect them. This approach could well put an end to distortion occurrence in the numerical grid solution. And, the other significant consideration factor of computation effort, especially for the situation of sophisticated time-dependent geometries, the meshless method shall balance better accuracy and less time. Therefore, the satisfactory performance of the meshless or particle method has been exploited and demonstrated to be conveniently implemented and computationally efficient.

Amongst the numerous grid-free and particle methods, the longest established meshless method is smoothed particle hydrodynamics that was originally developed in the late 1970s [13, 14], which has been extensively studied, explored, and conducted in many fields. With its wide attention in the area of computational fluid dynamics, large numbers of applications have been implemented which range over free surface flows [15], viscous flows [16, 17], multiphase flows [18, 19], geophysical flows [20, 21], the coastal engineering [22], turbulence modelling [23], viscoelastic flows [24, 25], and free-surface viscoelastic flows [26]. Though there are several advantages of SPH over grid-based methods, such as its Lagrangian and adaptive nature for the issue on large deformation and complex inconstant flows or easy programming, this meshless method is still confronted with numerical shortcomings which involve inconsistency and instability as well as incompatible boundary conditions treatment. The approach for enforcing the incompressibility condition in governing flow equations, namely, the incompressible SPH (ISPH), has been explored [27, 28]. In addition, an enhanced ISPH-SPH coupled method has also been proposed and studied for the simulation of incompressible fluid-elastic structure interactions [29]. Also, a great many studies have been devoted to improve the performance of this particle method, and the related research progress can refer to the literature cites [30, 31].

As another Lagrangian meshless method, moving-particle semi-implicit method (MPS) has been adopted in solving the hydrodynamics phenomenon, where the computational domain is similarly discrete and reconstructed by particles without the use of mesh in traditional methods. Hereby, there exist a number of studies on the flow problems resolved by MPS [32, 33], especially for the violent free-surface flows accompanied by large deformations and fragmentations, and it appears to be an appropriate candidate due to the particle characteristics [34, 35]. However, the unphysical pressure fluctuations in numerical computation procedure have also bothered the development of this meshless method.

Keeping this in mind and with the purpose to gain some insight into the real capabilities of meshless approach, a technique known as the finite particle method [36] has been extensively brought into effect for the numerical solution of a wide range of problems in computational mechanics. And, similar meshless discretization technique such as the finite-volume particle method (FVPM) has been employed for convergence analysis for compressible as well as incompressible flows [37, 38]. In addition, the validations of the finite particle method have been conducted through classical benchmark tests with exact analytical solutions, in terms of accuracy as well as conservation properties [39, 40]. The truly meshless method presented in references [4143] has been developed and implemented for the flow problems with complicated and rapidly changing geometry [44], free-surface flows [45], and multiphase flows [46]. In the finite particle method, particles which schlep the density and pressure as well as necessary physical characteristics move with fluid velocities. Similarly, boundaries are also approximated by a finite number of boundary particles, where the boundary conditions are also applied to these particles. Meanwhile, these points need not be regularly distributed and can even be quite arbitrary. Hence, the numerical treatment technique makes the finite particle method valid and efficient in the flow problems resolution.

In this thesis, we consider numerical simulation of incompressible viscous flow problems by the finite particle method. Numerical methods of interest are meshless Lagrangian particle methods. This approach is in some sense natural for fluid flows, since the governing equations themselves can be derived by considering small fluid particles which move with the flow. The incompressible Navier–Stokes equations are resolved by the projection method. This requires the solution of Poisson problems in each time step in which the spatial derivatives generated can be approximated by the moving least squares (MLS) method. As for the Neumann boundary condition, it is settled by placing stabilized finite points. In detail, the particle data management of the remeshing technique is adopted to circumvent the problems associated with irregular particle distribution. The proposed finite particle method is tested for three benchmark problems associated with incompressible viscous flows containing liquid sloshing of tanks with low volumetric fluid type, solitary waves movement, and interaction with a vertical wall in numerical flume and vortex patterns for the cross sections of ship rolling damping, in which the capability of finite particle method is demonstrated as a meshless method to solve the problems on incompressible viscous flow accurately and robustly.

The paper is organized as follows. In the next section, the theoretical principle of the finite particle method is expounded firstly. The application of the moving least squares (MLS) method for derivatives approximation is introduced in Section 2. We describe the explanation and discretization of the incompressible Navier–Stokes equations in Section 3. And, the particle data management as well as the resolution of time step has been explored in Section 4. Then, the numerical tests are presented in Section 5. The paper ends up with concluding remarks eventually.

2. Derivatives Approximation by the Moving Least Squares (MLS) Method

2.1. Moving Least Squares (MLS) Method

The moving least squares (MLS) method was developed as an approach to reconstruct the function for approximation in the combination form of polynomial basis, weighted with respect to the values of the scattered interpolation particles. Since the weights move through over the domain but with staid polynomial basis, this method is also designated by the name of moving weighted least squares (MWLS) method.

Firstly, let a function be given on a domain . Consider a point cloud . The function values are given at the data points by . The task is to approximate the function value at an arbitrary point , i.e., in the first place, we are confronted with a meshless interpolation problem.

u is assigned to be the function to be sought at control point x and its approximation function ; then, there iswhere the operator acts on a function and yields a new polynomial function which is constructed only from the values of at the neighbouring data points . denotes the coefficient vector, is the basis function vector, and is the number of basis functions. For the condition of 2D, the canonical basis is

While the canonical basis in 3D can be as follows:

The next step of the MLS method is to assign weights to the points via a distance weight function . The distance weight function is chosen smooth and decaying with increasing distance. The weight function used in this calculation takes the distance weight function and the specific expression asor splineswhere is a positive distance between any two particles in Euclidian norm, since it is rotationally invariant. The radius of the support domain h0 determines , the number of neighbouring particles around to be used for MLS approximation.

In the moving least squares (MLS) method, for the aim of minimum error square sum, a specific point is considered and the functionalshould be minimized, and the polynomial is defined by the coefficient vector , which minimizes the moving least squares distance at the data points . Since the weights are taken with respect to the point , for every point , one obtains a different polynomial and thus a different coefficient vector .

The coefficient vector is obtained by the normal equationwhere

Hence, for each point , a small linear system could be solved with matrix .

2.2. Approximation of the Derivatives Scheme by MLS

In the process of the finite point method, the derivative term of velocity and pressure arises in the numerical scheme. The MLS method with distance weight function is implemented to approximate derivatives for the pressure Poisson equation. While taking the weight function, such as equations (4) and (5), the MLS approximation can be differentiated at data points. The resulting expression will be a linear combination of function values at the data points .

Approximating function equation (1) by MLS, we can find that the coefficient vector is the solution to the systemwhere the matrix isand the right hand side is

Note that both and depend on . Consequently, as well as the basis vector depends on . In the following, we omit the -dependence in the expressions. For equation (9), the coefficient vector is calculated as

First of all, we are going to take the first-order derivative of both sides of equation (9); that is,where the prime denotes a differential operator (such as ), and thus,Then, the second-order derivative of both sides of equation (12) is taken; it yieldswhere the double prime denotes the operator applied twice, and then,

According to equation (14), the first derivative approximation of the MLS function can be obtained as

Similarly, using equations (14) and (16), we obtain the MLS function’s second derivative approximation as

Therefore, the spatial derivatives could be approximated by the above approach described in this subsection. The final sparse linear algebraic equations for the unknown pressure values are solved by employing an iterative scheme known as the stabilized biconjugate gradient (BiCgStab) method [47].

3. Discretization of the Incompressible Flow Governing Equations

3.1. Incompressible Flow Governing Equations

In the Lagrangian form, the incompressible Navier–Stokes equations can be written aswhere is the Cartesian components of the velocity field, is the Cartesian components of the position vector, is the time, is the fluid density, is the pressure, is the kinematic viscosity and , is the dynamic viscosity, and is the source term (normally the gravity ). denotes the total or material time derivative of a function .

3.2. Numerical Discretization by the Projection Method

The projection method which can be referred in reference [48] is employed in the finite particle method for the solution of equations (19) and (20) in time, which is completed in two fractional steps. The first step is to explicitly calculate the particle update position and the intermediate velocity from the momentum conservation equation (20):

The second step is to correct the intermediate speed and get the speed at the next time step by solving the equation

Equation (22) is observed, where the pressure gradient term is added and the new arising velocity can be handled by the incompressibility of the fluid:

Meanwhile, the pressure is calculated implicitly with the following Poisson equation for pressure deduced from equations (22) and (23); that is

The boundary condition for is obtained by projecting equation (22) on the outward unit normal vector to the boundary . Thus, the pressure Neumann boundary condition is written as

In the projection method, the particle positions change only in the first step. The intermediate velocity, the pressure, and the final divergence free velocity are all computed on the new particle positions. The spatial derivatives appearing in the above equations are approximated by the moving least squares (MLS) method proposed above which will directly be used to solve the pressure Poisson equation (24). Eventually, we use this mathematical model to numerically simulate the incompressible viscous flow problems.

4. Treatment of Particle Data and Time Step Size

4.1. Particle Data Management

In the finite particle method, points move with the fluid velocity. Since the velocity field is the solution itself, one cannot predict in advance where the particles will move. Hence, even if the particles are well distributed initially, they will in general cluster in some places and become scarce in others. The former can in best case result in an unnecessarily high local resolution and thus computational effort and in the worst case, yield numerical instabilities. The latter can spoil the accuracy and cause problems in the point neighbourhood search.

In order to remedy these problems, a constant particle data management has to be done. Particles too close to each other have to be removed, or better merged, while new particles have to be inserted into large holes and gaps. This requires on the one hand efficient detection procedures for nearby particles and holes and on the other hand, correct interpolation of the field data, when particles are merged and inserted.

When finding some points’ shortest distances closer than a given distance dmin by looping over all points, removing close points by the method of interpolating using all neighbouring points is worth effort. As for the condition with too many largest holes generating during flow, searching for all the Voronoi cell of point xi in 2D or in 3D and constructing the full Voronoi diagram [49] are the core missions. Once those largest holes are identified by Voronoi diagrams, new points will be inserted and participant calculation is done.

4.2. Determination of Time Step

For numerical stability, several time step constraints must be satisfied, including a Courant–Friedrichs–Lewy (CFL) condition,where is the maximum value of velocity for a given problem. And, additional constraints due to the hydrodynamical force acting on the particle ,

5. Numerical Examples

In this section, three numerical examples concerning free deformation of a viscous fluid patch with free surface are presented to test the performance of the finite particle method proposed in this paper to simulate incompressible viscous fluid flow. The numerical results are compared with the available experimental data for validating the presented grid-free methodology.

5.1. Liquid Sloshing of Tanks with Low Volumetric Fluid Type

The liquid sloshing of tanks is a typical strongly nonlinear flow problem, accompanied by traveling wave, breaking wave, splashing and high-speed lashing, and some other complicated flow phenomena. Especially for the high-pressure impacting force on the tank structure, this condition may result in local structural deformation or damage; therefore, it is well worth intensively studying. This numerical example is used for the simulation of liquid sloshing of tanks with low volumetric fluid by the proposed finite particle method, analysis of sharp motion characteristics and large deformation of the fluid flow, periodic slamming pressure on tank walls, and comparison with test results.

For consistency with the experiment condition in reference [50], liquid sloshing of tanks with low volumetric fluid type of the same liquid capacity (40%) but with two different sway excitation periods is selected and numerically simulated. Model I has the period of excitation 1.5 s and Model II has excitation period 1.3 s, while the resonant period of tank is found to be 1.3 s. At the same time, the pressure monitoring point O is set in the calculation to compare with the measured pressure acting on the tank bulkhead in the test, quantitatively analysing the accuracy of the numerical calculation results. The calculation model of the tank which is shown in Figure 1 takes a two-dimensional rectangular chamber with a width of 0.6 m and a height of 0.3 m.

The sloshing of the liquid is generated by the regular sway motion of the tank, and its expression iswhere is the excitation amplitude and is the corresponding maximum amplitude of 0.05 m.

5.1.1. Liquid Sloshing of Tank Model I

In Model I, the excitation frequency of the tank with a filling rate of 40% is less than the resonance frequency. Figure 2 shows the calculation results of the turbulent Model I at different times in a cycle (t = 0.1 T, 0.2 T, 0.3 T, and 0.4 T), given in the form of the snapshots of water particles together with the pressure field.

It can be seen from Figure 2 that, at t = 0.1 T, the fluid in the tank moves to the right wall under the action of the excitation force, accompanied by a traveling wave. When it reaches 0.2 T, the wall-piercing fluid starts to rise along the wall under the action of inertia. As t = 0.3 T, the fluid moving along the right side wall hits the tank ceiling again and the liquid at the bottom accumulates more and more. It brings about the opposite movement “water bag” under the excitation movement of the tank. While approaching 0.4 T, the fluid slamming the right side wall begins to move to the left side of the bulkhead in a raised “water bag” shape. Under the action of the bulkhead excitation force, the fluid has sloshed in different forms with strong nonlinear phenomena such as traveling wave, water jump, and rolling. The free surface flow of these complex deformations is accurately simulated by the finite particle method. Compared with the traditional numerical simulation method based on mesh, the meshless method not only eliminates the numerical dissipation problem of convection terms but also has strong ability of capturing viscous flow, which has guiding significance in a certain degree to some other similar solution methods for the resolution of large-scale motion of free surface flow.

In a further step, to quantitatively verify the accuracy in liquid sloshing problems by the finite particle method, the pressure measuring point O of Model I is numerically calculated and obtained. To be consistent with the experimental measurements [50], the calculation pressure monitoring point O is in the same height of 0.1 m, whose time-history pressure is presented in Figure 3.

From the presented Figure 3, it is observed that the numerically calculated pressure peaks produce a certain range of nonphysical oscillations. For the first four cycles, the pressure peaks are larger than those of the experimental result. From the comparative analysis of the two curves in the first cycle, it has been shown that the peak value of the test measurement is about 1/2 of the peak value of the numerical calculation. But, in the subsequent cycle, the numerical calculation results are basically consistent with the experimental measurements. Comparing the numerical result and test data, we can see that there are basically two pressure peaks in each cycle, and the second pressure peak remains the same with the experimental values in general. The overall trend of the pressure fluctuation is consistent on the whole. Hence, through quantitative comparison of the time histories of pressure, it has been demonstrated that the finite particle method has a good reliability to solve liquid sloshing of tanks problems.

5.1.2. Liquid Sloshing of Tank Model II

The sloshing Model II is with the same low volumetric rate of 40% and an excitation period of 1.3 s. It is noted that this model excitation period remains in accordance with the resonant period of tank. Figure 4 presents a qualitative comparison corresponding to Model II of the free surface profiles from the experiment [50] and numerical methods in which the simulation results are depicted by the snapshots of water particles together with the pressure field at two instants (t = 0.2 T and 0.4 T).

As can be seen from Figure 4, the sway motion of the tank excites the liquid to violently slosh, in accordance with noticeably convex free surface because of the same excitation frequency with the resonance frequency. In particular, the time reaches 0.2 T, the fluid pressure near the right side wall of the tank increases sharply, forming a large pressure gradient. While t = 0.4 T, due to the strong slap of the top and left bulkheads with the liquid, the fluid is caused to roll down and accompanied by splashing water. By qualitative comparison between the numerical solution and the experiment, these complex deformations are accurately captured by the finite point method and are in good agreement with the experimental results.

5.2. Solitary Waves Movement and Interaction with a Vertical Wall in Numerical Flume

Solitary wave motion and interaction with solid as a strong nonlinear physical process, accompanying with wave climbing and fragmentation, is endowed with significant research value in the field of coastal engineering. Due to the complexity of its physical process, it is difficult to fully explain it mathematically. With the development of experimental techniques and numerical calculation techniques, relevant studies have achieved certain progress. The subsection will establish the numerical flume based on the finite particle method to perform numerical simulation of the problem on solitary waves.

In numerical wave flume, the solitary wave is generated by pushing a wavemaking plate on the wave source position, whose surface equation can be defined by the first-order solution of Boussinesq equation:where presents the vertical coordinates in wave section, the wave height, the water depth, the horizontal coordinates, the time, the wave velocity and , and is the acceleration of gravity. By deducing from equation (29), the displacement expression of wavemaking plate is

The displacement values are passed to the limit, respectively, at and , and they are subtracted each other to get the maximum displacement of wavemaking plate (stroke) as follows:

It is transparent that as , . Let T be the truncation period of wavemaking plate motion and there is at . This relational expression can be brought into equation (30), so the period can be expressed as

The zero displacement point is set at , whereupon the displacement is redescribed by the following equation under different values of relative wave heights .

The dimension of the numerical wave tank model is of , and the computational domain is instead by particles with the water depth . For the numerical conditions of and , Figure 5 presents the analytical free surface at a different time through solitary wave propagation and interaction with a vertical wall.

As can be seen in Figure 5, either for or condition, when the solitary wave begins to climb the wall, the shape of the solitary wave has changed significantly. The front peak shape becomes steeper and the back remains flat. As the solitary wave climbs to the highest on the vertical wall, it has been broken and the water particles begin to fall along the straight wall. Because of the gross movement, the fluid possesses a large potential energy. It results in a reflective tail wave after the particles fall back. It could plainly be found that the wave height of the reflective solitary wave is obviously smaller than that before climbing and the wavelength becomes larger, indicating that the climbing motion makes it lose a large amount of energy. This numerical case indicates that the meshless method has a strong ability of simulation of the entire physical process of the solitary wave propagation and interaction of climbing along the vertical wall and falling back, where the waveform is kept well and visible. In a further step, the finite particle method has good applicability and accuracy especially for the large deformation of the viscous free surface flow.

5.3. Vortex Patterns for the Cross Sections of Ship Rolling Damping

The numerical simulation of vortex patterns for the cross sections of ship rolling damping has been conducted by the finite particle method proposed in this subsection. The numerical simulation has effectively depicted the nonlinear characteristics, such as flow separation phenomenon and vortices outflowing. Also, the results of comparison between the numerical simulation and experiment are investigated.

5.3.1. Rolling Computation Model Representation

In the calculation of the viscous flow field, the two-dimensional rolling object is selected by the cross sections of the S60 ship. The corresponding dimension is that the length between perpendiculars of the ship model Lpp = 3.569 m, the breadth moulded B = 0.477 m, and the draught d = 0.190 m. Ikeda has conducted a series of experiments on two-dimensional rolling of S60 ship, and the specific test data and detailed discussion can be found in reference [51]. In this paper, three cross sections in different positions of the ship model based on large differences in shape are determined for numerical calculation of two-dimensional forced roll motion. The ship model and each cross section selected are shown in Figure 6.

The rolling centre of each two-dimensional cross section is located at the intersection of the static fluid surface and the longitudinal section. The rolling angular displacement equation iswhere the maximum angular displacement and the circular frequency . The specific geometric dimensions and related calculation parameters of each cross section are shown in Table 1.

In Table 1, S denotes sectional area, H0 and σ denote Lewis profile parameters, H0 = B/2d, and σ = S/Bd.

5.3.2. Results Analysis of Vortex Patterns for Ship Rolling Damping

Through the rolling computation by the finite particle method, the velocity vector distribution of the flow field around the three cross sections (S.S.1.1/2, S.S.5, and S.S.8.1/2) of the S60 ship in each period is depicted in Figures 79. Observaing the velocity vector distribution of each cross section in Figures 68, it is known that the vortex is periodically formed and exterminated near the ship bilge with the roll motion of the hull, which is basically consistent with the test conducted by Ikeda. Through the analysis of rolling motion, it can be seen that there are two vortexes on sides of the bilge for the cross section S.S.5 located in the middle of ship. However, only one vortex in relatively high lever is generated at the bottom either for section S.S.8.1/2 close to the stem or section S.S.1.1/2 near the stern. The numerical result is very similar to that of the literature [51].

Ikeda had concluded that, through a series rolling test of cross sections for S60, there are two types of vortexes usually formed in ship sections undergoing rolling motion which are the double bilge vortexes near the midship section and the single bottom vortex at the stem or stern section. The fundamental principle of the vortex patterns classification can be determined according to the Lewis profile parameters (H0, σ) (Figure 10). The rule diagram of Ikeda’s S60 series rolling test is also drawn in Figure 10 [51].

Associating Figure 10 and the data in Table 1, as the Lewis profile parameter H0 deviates far from 1, the cross section S.S.8.1/2 will only form a strong vortex at the bottom, which is the same with that numerical condition in Figure 9. When H0 = 1, it is necessary to combine σ to further determine the number of vortexes. If the value is near 1, a vortex will be formed on both sides of the cross section bilge. As the Lewis profile parameter σ decreases, the vortex discharge position would move along the downstream and the vortex strength shall also reduce. The experiment conclusion has verified the calculation results of S.S.5 cross section in Figure 8. When the value of σ continues to decrease below 0.6, only one vortex will exist near the bottom of the cross section, such as S.S.1.1/2 (Figure 7). By exploring the vortex patterns for the cross sections in Figures 710, it can be concluded that the numerical results are in good agreement with Ikeda’s experimental phenomenon, which validates the effectiveness and accuracy of the finite point method for solving the ship rolling motion problems.

6. Conclusions

A grid-free modelling based on the finite particle method is proposed in this paper for accurate simulation of incompressible viscous flows. The moving least squares is introduced for approximating the derivatives. And, the governing equations of fluid flow are solved using a projection method. The pressure Poisson equation with Neumann boundary condition is handled by the biconjugate gradient method. This concept is shown to improve the accuracy of the finite particle method for simulating complex-free surface flows. The meshless method has also been used to simulate three different numerical examples including liquid sloshing of tanks with low volumetric fluid type, solitary waves movement, and interaction with a vertical wall in numerical flume and vortex patterns for the cross sections of ship rolling damping whose results are presented and compared with corresponding experimental data. The results keep a good agreement between numerical calculation and the test. In addition, the method is also shown to be able to produce accurate and stable results for simulating strongly nonlinear flow problems, which indicates the finite particle method is an accurate and efficient numerical method for the solution of incompressible viscous flows.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (grant no. 51809029) and the Fundamental Research Funds for the Central Universities under grant nos. 3132019113 and 3132018207.