2017 Volume 93 Issue 8 Pages 525-546
For innovative designs of coastal structures, Numerical Wave Flumes (NWFs), which are solvers of Navier-Stokes equation for free-surface flows, are key tools. In this article, various methods and techniques for NWFs are overviewed. In the former half, key techniques of NWFs, namely the interface capturing (MAC, VOF, C-CUP) and significance of NWFs in comparison with the conventional wave models are described. In the latter part of this article, recent improvements of the particle method are shown as one of cores of NWFs. Methods for attenuating unphysical pressure fluctuation and improving accuracy, such as CMPS method for momentum conservation, Higher-order Source of Poisson Pressure Equation (PPE), Higher-order Laplacian, Error-Compensating Source in PPE, and Gradient Correction for ensuring Taylor-series consistency, are reviewed briefly. Finally, the latest new frontier of the accurate particle method, including Dynamic Stabilization for providing minimum-required artificial repulsive force to improve stability of computation, and Space Potential Particle for describing the exact free-surface boundary condition, is described.
To show a necessity of Numerical Wave Flume (NWF), current structurally-resistive design against wave action is overviewed. The current structurally-resistive design is based on empirical, experimental, or deterministic design concepts, and they cannot provide sufficient hydraulic information for a detailed design method such as a reliability design.
The structurally-resistive design against wave action requires two keys of information on (i) wave force and (ii) wave overtopping. For the wave force, design formulae exist for all of different types of structures, such as piles (columnar structure), vertical wall and armor materials. Morison formula (Morison et al., 1950)1) is famous for estimating wave force acting on piles. This formula displays wave force acting on a pile for the sum of drag in proportion to square of the approaching velocity and inertial force. It is an empirical formula in which stress distribution due to a complex flow around a pile is simplified into two coefficients, namely drag coefficient and inertia coefficient. These coefficients are experimentally estimated for a single pile in sufficiently large domain. In most of practical cases, group of the piles is often used as a foundation of structure, hence it is necessary to estimate the coefficients of piles group. As for a pile, modified formula for irregular wave and wave breaking are investigated experimentally. Considering infinite variation of arrangement patterns of piles, it is impossible to conduct hydraulic experiments to find universally appropriate coefficients of Morison formula for piles group.
A drag force can be estimated without a model constant by surface integral of stress on a pile surface, if the stress distribution on a pile surface is sufficiently resolved. In addition, a frictional resistance is approximately negligible in the most of practical cases, a drag force is computed by surface integral of pressure on a pile surface. For this computation, an NWF, which can directly calculate the time history of the pressure distribution around a pile, is useful.2)
An estimation of the wave force acting on a vertical wall is necessary for design of a caisson breakwater. Goda formula (Goda, 1973),3) which describes both of pressure due to standing wave and impact pressure by wave breaking, is widely used. The Goda formula is the superior one that can simply estimate maximum wave pressure. But it has been experimentally established for standard breakwater sections, without covering various cross-sections, such as a caisson with water chambers for wave dissipation. In addition, it describes impact or maximum pressure acting on a caisson breakwater, but it cannot treat time history of pressure. The NWF, which can calculate pressure distribution around the caisson and time history of boundary-neighboring velocity directly, is necessary to examine the stability of structure especially for cases under the existence of liquefaction and scouring around the caisson.2)
Hudson formula (Hudson, 1959)4) is used for assessing stability of armor material, or rubble, of the mound. This formula estimates a threshold of stability, namely an initiation of motion of rubble, from a static balance of lift force, which is proportional to square of the approaching velocity, and gravity force for rubble placed on a slope. Although a lift coefficient in Hudson formula is estimated by an experimentally-investigated coefficient for a single body, rubble often exists in swash zone with repeating dry-and-wet condition, and a sheltering effect must be significant for pick-up event of rubble.
As for a wave overtopping, many hydraulic experiments for estimating run-up height and wave overtopping quantities have been carried out and summarized as diagrams parameterized by wave conditions (wave height, and wavelength) and topographical conditions (water depth, and bottom slope). As for wave overtopping quantity onto a vertical wall, Goda diagram (Goda et al., 1975)5) is widely used for a design of coastal dikes. Because water surface profile is complex in run-up and wave overtopping processes, and even a topology of water surface may change in case of a plunging breaker, wave transformation models with a depth-integral cannot reproduce a time dependent change of water depth. In such a case, a water-surface tracking model, which can work flexibly and accurately under a topological change of water surface, must be introduced. Hence an NWF should be introduced.
As shown above, around a structure, there exists a violent flow with wave breaking, and empirical design formulae have limits in accuracy. In recent years, combination of wave-control measures and structures based on a concept of the integrated shore protection is taken into the coastal designs. In such situation, structures interact with each other complicatedly, and an evaluation only by a combination of empirical formula becomes insufficient. Then there is no choice to conduct a one-by-one hydraulic experiment using a scaled-down model of the structure. Each hydraulic experiment needs enormous cost and time. Furthermore, the similarity with a full-scale phenomenon is not guaranteed enough, and sufficiently accurate data may not be provided by a scaled-down model experiment. As a method to collect data effectively without taking costs, technologies based on Computational Fluid Dynamics (CFD), including the NWF, are hoped to be a core tool of the design of coastal structures.
Although the NWF is a computational tool for coastal wave, its target is not limited to a water-surface wave. The NWF can simulate a free-surface flow, the governing equation of which is Navier-Stokes equation. In case of an incompressible Newtonian viscous fluid, continuity equation and Navier-Stokes equation are as follows:
| \begin{equation} \frac{D\rho}{Dt} + \rho\nabla\cdot\boldsymbol{{u}} = 0 \end{equation} | [1] |
| \begin{equation} \frac{D\boldsymbol{{u}}}{Dt} = - \frac{1}{\rho}\nabla p + \nu\nabla^{2}\boldsymbol{{u}} + \boldsymbol{{g}} \end{equation} | [2] |
In hydraulic engineering, only water (liquid phase) is usually calculated, hence an equation of a single-phase flow is the governing equation. When motion of atmosphere has to be calculated in addition to a water flow, a governing equation of a gas-liquid multiphase flow is necessary. Furthermore, when smaller-scale fluid motion than the grid scale cannot be ignored to describe the grid-scale flow, a sub-model to describe turbulence is necessary.
Here, we presuppose that all information on flow behavior, as far as required in engineering viewpoint at least, is provided, when Navier-Stokes equation is solved in sufficiently-high resolution in both of spatial and time domains. Although multiphase flow or non-Newtonian fluid requires a modeling of the phase interaction or a constitutive law, we stand in the viewpoint that a computation of Navier-Stokes equation gives sufficiently accurate physical quantities at least for a single-phase flow.
Appropriate discretization is necessary to obtain a numerical solution of Navier-Stokes equation, which is a governing equation for continuum. The discretization applied to a fluid computation widely is Eulerian method with a calculation grid. In Eulerian method, the calculation grid is placed to cover a fluid-existing domain, and the physical quantities are defined discretely at grid points.
The operation to place the calculation grid is called pre-processing. The Navier-Stokes equation is transformed into simultaneous algebraic equations, which describe a relation of physical quantities at neighboring grid points. And velocity and pressure at the grid points are provided as a solution of the algebraic equation. The visualization processing, or post-processing, of the calculation results is carried out afterwards.
2.2 Poisson equation and pressure calculation.After the discretization of governing equations, velocity at the next time step is obtained from the explicit calculation of the discretized Navier-Stokes equation. However the resultant velocity field doesn’t satisfy the continuity equation because the pressure field is not updated simultaneously with the velocity field. In order to compensate this deficiency of the explicit scheme, Poisson pressure equation, which is derived by substituting the Navier-Stokes equation into the continuity equation, is introduced and solved in the implicit way.
There are several methods to solve the Poisson pressure equation. The common process is that it is numerically solved with an iteration process so that the velocity vectors are adjusted to satisfy the continuity equation with the minimum accumulation error under optimized distribution of scalar pressure. It should be noted here that in the process, calculated pressure values have a role to compensate numerical errors at the previous and current steps. This means that the obtained pressure may not be the physically exact pressure and it sometimes shows large fluctuation or a false value.
2.3 Detection of phase interface.Detection of the free surface is classified into an interface-tracking method and an interface-capturing method. Because the interface-tracking method defines the interface as a boundary of a deforming calculation domain, and the coordinate of the interfacial calculation point strictly coincides with the position of interface, the interface is tracked accurately. However, when interface-tracking method is applied to a grid method under the topological change of the water surface from simply connected to biconnected like a plunging breaker, neighboring grids of water surface are extremely distorted and a calculation breaks down due to a crush of the grid.
The application of the interface-tracking method is therefore limited to the condition with relatively small wave steepness without a wave breaking. When wave steepness increases, the ALE (Arbitrary Lagrangian-Eulerian) method (Hirt, et al., 1977)6) is applicable, but even ALE cannot calculate a plunging breaker with topological change of the water surface. In the particle method, which is a complete Lagrangian method free from a computational grid, interface-neighboring particles are given interfacial boundary condition. Thus, the particle method belongs to interface-tracking method.
The interface-capturing method places a grid to cover the domain where the interface is assumed to pass through during a calculation, and an interface is captured by the advection calculation of a function representing the interface. Interface-capturing method consists of three methods: a method with a height function, a method with Lagrangian markers, and a method solving an advection equation.
The height function is simple, but it is essential for water surface to be a single-valued function and is not applicable to the overhanging water-surface which is found at the wave peak at the moment of plunging breaker generation. The MAC (Marker and Cell) method (Harlow and Welch, 1965)7) is well known as a method with Lagrangian markers. Advection of a large number of markers placed near the interface is tracked to capture an interface. The MAC method is theoretically superior one that can deal with a topological change of the water surface, but it is not efficient because high-calculation load of marker tracking is required only for water-surface capturing and governing equations are solved by a grid method.
A method to solve an advection equation of a physical quantity (for example, volume fraction of fluid in a grid cell) for an identification of the water surface is effective to reduce computational load of the interface tracking. Both of VOF (Volume of Fluid) method (Hirt and Nichols, 1981)8) and Level-Set method (Sussman et al., 1994)9) are categorized as interface-capturing method.
Pioneering technique of the interface capturing is the VOF method, which computes the volume fraction of fluid in a grid cell. At first, this approach was devised as a method for improving the volume-tracking performance of MAC method with saving memory and CPU cost. For a gas-liquid two-phase flow, VOF method introduces volume fraction, F, in each cell as a variable in place of markers in MAC method. The interface exists between the cells filled with liquid and the cells filled with gas. Otherwise, the interface exists between the cells filled with liquid-gas mixture and the cells filled with liquid. In other words, the interface can be detected at the cells satisfying 0 < F < 1.
However, the interfacial direction, or the direction of uneven distribution of liquid in the cell, cannot be estimated only based on the F value. A Donor-acceptor method8),10) is necessary to solve this problem. Because the Donor-acceptor method calculates the advection with taking uneven distribution of liquid in the up-stream-side cell into account, an interface profile does not become discontinuous. In addition, the Donor-acceptor method satisfies volume conservation. However, in an actual computation, an interfacial cell (F < 1) arises inside of fluid due to computational errors, and volume conservation becomes incomplete.
As for the grid methods except for VOF method, there are techniques to solve gas-liquid two-phase flow directly: the density function method, the C-CUP (CIP-Combined and Unified Procedure) method (Yabe and Wang, 1991)11) and the Level Set method.8) In the density function method, the color function, whose value jumps from 0 to 1 across an interface, is defined, and the advection equation of the color function is solved to capture a gas-liquid interface. In the C-CUP method, the advection calculation of the color function is carried out in the same manner with the density function method. In the Level Set method, isosurface (closed surface) of the signed distance function (actual value scalar function) depending on the distance from the interface is introduced to identify an interface, and an advection equation of the distance function is solved to capture an interface.
In the methods in which an interface is captured by calculating advection of color function, an interface becomes unclear with calculation progress due to numerical diffusion. In the Level Set method, an interface is kept to be clear with attenuating numerical diffusion, but reinitialization, or reset, of distance function is necessary when an interfacial profile changes due to advection, and the algorithm becomes complicated. Since it is caused by the nonlinearity of the advection equation, the numerical diffusion is controlled effectively by a higher-order calculation scheme of advection. The CIP (Cubic Interpolation Pseudo-particle or Constrained Interpolation Profile) method (Yabe and Aoki, 1991;12) Yabe et al., 200113)) is a calculation scheme to overcome problems of a numerical diffusion and a phase error at the same time. The CIP method enables sharp interfacial capturing in a concise algorithm. However, the interface becomes unclear with calculation progress, because even the CIP method cannot completely control numerical diffusion. The unclear interface causes a domain classified as neither the liquid phase nor the gas phase and brings a problem of violation of volume conservation. Therefore, it is necessary to compensate error of volume of the gas-and-liquid phases during computation.
CADMAS-SURF (SUper Roller Flume for Computer Aided Design of MAritime Structure) (Isobe et al., 1999)14),15) is a solver of the Navier-Stokes equation of the single-phase flow based on VOF method. And in late years, CADMAS-SURF is used frequently in structurally-resistive design against wave action in coastal engineering. CADMAS-SURF can handle practical topographies such as a sloping sea bottom, a porous wave absorbing structure and so on. CADMAS-SURF furnishes a large number of functions specialized in coastal and harbor engineering: wave generation source, wave-making function for permanent progressive wave or the arbitrary wave pattern, non-reflecting boundary processing, setting of the arbitrary shape structure, k-ε turbulence model, and air bubbles and water drop processing.
2.4 Connecting wave models based on Boussinesq equation.Since numerical models based on the Navier-Stokes equation are very costly, wave models are frequently used to give necessary information at the offshore boundary of the calculation. Boussinesq (1872)16) derived a set of equations for weak nonlinear and dispersive waves in very shallow water by using Taylor expansion of velocity potential around the bottom. Boussinesq-type equation has many variations. Peregrine (1967)17) showed a derivation in terms of the depth-averaged velocity with taking into account the bottom slope. Modified Boussinesq equations (e.g., Nwogu, 1993;18) Beji and Nadaoka, 199619)), that are applicable to irregular wave at an arbitrary water depth, are widely used for computation of wave propagation. Boussinesq equation was derived for wave field without energy dissipation under the assumption of inviscid and irrotational flow, and is not theoretically applicable to wave fields with energy dissipation, such as breaking waves and turbulent flow field around coastal structures.
Being applied to wave breaking, a resistance term proportional to velocity, a dispersive term proportional to second-derivative of velocity in wave-propagating direction, or a forced-dissipation term depending on water-surface profile (e.g., Schäffer et al., 1992)20) is usually added for reproducing energy dissipation artificially. The energy dissipation by these additional terms does not have a physical basis, it should be therefore understood as an artificial damping. Treating a water-surface motion under energy dissipation due to wave breaking as a pure wave is questionable, hence, there exists a limit for calculation of breaking wave propagation and transformation by wave equations. The wave equation cannot describe three-dimensionality, unsteadiness, and intermittency, which are essences of a turbulent flow. As a result, it cannot sufficiently reproduce an energy dissipation due to turbulence. To deal with these problems, by getting rid of the limitation as the wave, Navier-Stokes equation must be solved directly. It is the NWF that enables this approach.
The NWFs are applicable to not only the structurally-resistive design against wave action but also various problems relating to the coastal hydraulics. In the field of sediment transport dynamics, particularly sediment movement in the surf zone, contribution from the computational science is of the great importance. Because methods to track intensely-fluctuated water surface were not well developed, the detailed calculation of the water surface and velocity field in the surf zone was difficult with the conventional wave models. However, rapid progress of the NWF allows detailed understanding of the time-and-space structure of the velocity field in the surf zone.
A Lagrangian model in which motion of sediment particles are directly calculated following the mechanism of sediment movement is effective for evaluation of sediment transport, calculation load, however, becomes too heavy in cases of simultaneous tracking of a large number of sediment particles. To decrease the load, a method to track the density of sediment based on an advection-diffusion equation is commonly applied instead of tracking each individual sediment particle. In the practical calculation of beach profile changes, it is even unrealistic to calculate advection-diffusion for density of sediment, because a target domain is too large and the tracking time for phenomenon is too long. For instance, CERC formula,21),22) which estimates longshore sediment transport rate by referring the longshore-directional component of energy flux at the wave breaking point, has been used for a long time, despite its ignorance of the detailed dynamics of sediment.
In the surf zone, which is defined as a domain on the onshore side of the breaker line, the wave energy dissipates through turbulence energy. The surf zone is a complex and unsteady turbulent field, in which sediment transport due to turbulent flow is quite active. The surf zone is a source of the suspended sediment, which is driven by organized unsteady turbulence due to a wave breaking. It is therefore important to understand the physics of the flow in the surf zone in order to predict the advection and diffusion of sediment accurately.
In the numerical models before the NWFs, net steady velocity due to wave breaking was of importance as well as the oscillatory motion of water. The onshore velocity was divided into four components: the steady-flow component, the wave component, the wave-breaking-induced near-surface vortices (surface roller) component and the turbulence component. The near-surface vortices, which are generated by wave breaking, are transported in onshore direction and bring strong turbulence near the water surface. Offshore-directed steady velocity (undertow) can be estimated with an empirically-given mass transport due to surface roller (e.g., Okayasu et al., 1988).23) The onshore sediment transport rate can be easily predicted from the near-bottom steady velocity, but an estimation of the velocity based on the Navier-Stokes equation, namely NWF, is required to evaluate suspended sediment transport in the turbulent flow field with strong unsteadiness in the surf zone.
For a calculation based on the Navier-Stokes equation, Direct Numerical Simulation (DNS) which calculates the velocity consistent with a micro-scale turbulence directly, is ideal. But the DNS needs an enormous calculation cost to ensure the resolution in Kolmogorov scale, hence its application is restrictive (e.g., Wijayaratna and Okayasu, 2000).24) Suspended sediment clouds on ripples in the offshore zone (Sunamura et al., 1978)25) characterize the onshore sediment transport driven by a periodic oscillatory flow. In order to calculate it, Reynolds-Averaged Navier-Stokes equations (RANS) has been applied before an adoption of Large Eddy Simulation (LES). Sato et al. (1986)26) simulated the advection diffusion of the suspended sediment on ripples by the k-ε turbulence model. To conduct a calculation in a scale of the actual surf zone, LES, which describes the motion of vortices smaller than the grid size by Sub-Grid-Scale (SGS) turbulence model is one of candidates (e.g., Christensen and Deigaard, 2001;27) Okayasu et al., 200528)).
Even if unsteady turbulent flow field due to a wave breaking was calculated by LES, a model describing motion of suspended sediment is necessary. Since the inertia of the suspended sediment is small, the suspended sediment should follow well the movement of surrounding fluid. Therefore, motion of suspended sediment can be estimated by solving an advection-diffusion equation in velocity field calculated by LES. Suzuki et al. (2007)29) simulated the density distribution of suspended sediment in the surf zone by using LES and the advection-diffusion equation. In the calculation of the advection-diffusion equation, estimation of the bottom boundary condition, that is source or sink, of the suspended sediment density is necessary. Nielsen (1992)30) proposed a conceptual model providing time-dependent density flux of suspended sediment at the bottom boundary. In this model, the density flux of suspended sediment from the bottom which is called “the pick-up function” is described as a function of the instantaneous Shields parameter, specific mass of sediment and sediment diameter.
However, in Nielsen’s model, the pick-up from a bed-load motion is estimated empirically, and the physics of transition mechanism from bed-load motion to suspension is still unclear. The Lagrangian model, namely the particle tracking method, which treats motion of sediment particle directly, is then necessary to discuss dynamics of suspended-sediment generation in detail. As for the particle tracking method, the sediment transport which was not affected by the wave breaking such as the formation of the suspended sediment cloud on vortex ripples has been targeted so far. However, now complex sediment transport in the surf zone can be directly computed by applying the sediment transport mechanics with the NWF which enables detailed understanding of the velocity-field structure in the surf zone.
The sediment particle is accelerated by hydrodynamic force, and the fluid loses a momentum. In other words, the solid-liquid multiphase flow model, which evaluates momentum transport between liquid-and-solid phases, is necessary. As for the bed load, the momentum transport between liquid-and-solid phases becomes active, since the relative velocity of sediment particle to the surrounding fluid cannot be negligible, due to a frequent collision of sediment particles with the bottom surface. Therefore, the interaction force between liquid-and-solid phases should be introduced (e.g., Gotoh et al., 1995).31) On the other hand, as for the mean velocity of the suspended sediment, the effect of the momentum transport between liquid-and-solid phases is quite small, because there is no significant difference in the velocities of sediment particle and the surrounding fluid. However, the change of the interaction force between liquid-and-solid phases gives a change to turbulence structure. Therefore, in transport equation of the turbulence, it is necessary to consider a correlation term of fluctuating components of fluid velocity and liquid-and-solid interaction force (e.g., Gotoh and Sakai, 2000).32)
A motion of gas-liquid interface is greatly affected by the turbulence structure of the surf zone. Particularly, in the plunging breaker, wave crest overturns and splashes down to form horizontal vortices (horizontal roller), and strong turbulence is produced at the plunging location. Until a jet lands on the water surface after wave breaking, a two-dimensional transformation of water surface is predominant, hence the water surface can be tracked by the two-dimensional model (e.g., Watanabe and Saeki, 2002).33) It changes afterwards to a bore while repeating splash-up. In the domain of the bore region, the organized structure of three-dimensional vortices called obliquely descending eddies (Nadaoka et al., 1989)34) are prominent. A measurement was difficult in this domain, because a large quantity of air bubbles are mixed into water. Watanabe et al. (2005)35) conducted the three-dimensional single-phase turbulent flow computation and showed that a reorientation of the initial two-dimensional spanwise vortices into the coherent streamwise-and-vertical vortices was produced by a strong stretching force and a shear instability. Furthermore, Watanabe et al. (2008)36) simulated the process when a finger-shaped jet characterizing splash-up was formed from the water surface with longitudinal counter-rotating vortices.
The solid-solid interaction would give a dominant effect to the high-density sediment-laden flow, thus, it is important to reproduce the physical solid-solid interaction directly with describing sediment particles as rigid bodies. For this problem, Distinct Element Method, abbreviated as DEM (Cundall and Strack, 1979),37) which is a soft sphere model, has been widely applied to express the physical solid-solid interaction, namely repulsive/contact force acting on sediment particles (e.g., Gotoh an Sakai, 1997).38) The easy approach to couple DEM with flow model is the one-way coupling by giving the constant drag force to solid particles (e.g., Reddy and Kumaran, 2010;39) Estep and Dufek, 201340)).
In late years, studies on a highly-concentrated solid-liquid multiphase flow have became active by coupling a solver of the Navier-Stokes equation with DEM (e.g., Tang et al., 2016;41) Mondal et al., 201642)). The NWF, or the solver of liquid phase, enables the development of the computational science of the sediment transport by being coupled with DEM, or the solver of solid phase (e.g., Harada et al., 2015,43) 201644)).
Grid methods suffer from numerical diffusion associated with the discretization of an advection term, while a particle method is free from numerical diffusion, because a substantial derivative is equivalent to a time derivative in Lagrangian description. In this context, a particle method possesses significant potentials to provide fairly accurate solutions with moderate computational costs for violent free-surface flows. In other words, a particle method can show robustness in reproduction of violent free-surface flows with fragmentation and coalescence of water.
Particle methods are categorized into two major methods, namely Smoothed Particle Hydrodynamics (SPH) method (Lucy, 1977;45) Gingold and Monaghan, 197746)) and Moving Particle Semi-implicit (MPS) method (Koshizuka and Oka, 1996).47)
In coastal and ocean engineering, particle methods had been applied to wave breaking (Gotoh and Sakai, 1999;48) Khayyer and Gotoh, 2008;49) Farahani and Dalrymple, 201450)), wave overtopping (Gotoh et al., 2005;51) Shao et al., 200652)), wave run-up (Shadloo et al., 2015),53) wave impact (Khayyer and Gotoh, 2009;54) Lee et al., 2011;55) Altomare et al., 201556)), wave-induced nearshore current (Farahani et al., 2014),57) violent sloshing (Delorme et al., 2009;58) Gotoh et al., 201459)), green water on ships (Shibata and Koshizuka, 2007;60) Le Touzé et al., 201061)), sediment transport (Gotoh and Sakai, 2006),62) landslide-generated waves (Panizzo and Dalrymple, 2004;63) Fu and Jin, 201564)) and fluid–structure interactions (Rafiee and Thiagarajan, 2009;65) Shibata et al., 2012;66) Hwang et al., 2014;67) Colagrossi et al., 2015;68) Wei et al., 201569)).
Either explicit algorithm or implicit algorithm can be applied to particle methods, but both of these algorithms require a stabilizer to ensure numerical stability. WCSPH (Weakly Compressible SPH) (Colagrossi and Landrini, 2003)70) method with an explicit algorithm must introduce artificial viscosity for stable computation. While MPS method and ISPH (Incompressible SPH) method (Shao and Lo, 2003),71) both of which are based on the projection method with semi-implicit algorithm, require artificial repulsive force. Hence, how to minimize stabilizer, namely artificial viscosity or artificial repulsive force, is one of the main issues in development of accurate particle method.
In semi-implicit algorithm, a projection onto the solenoidal field provides error clean-up in pressure field, while in explicit algorithm, re-initialization of density field to remove error accumulation is essential. In the context of error control, semi-implicit algorithm has superiority to explicit algorithm.72) Taking before-mentioned superiority of semi-implicit algorithm into account, in this paper, particle methods with semi-implicit algorithm are focused.
4.2 Standard MPS method.In MPS method, which is based on a projection method, the pressure is obtained by solving the Poisson Pressure Equation (PPE) as follows:
| \begin{equation} \langle\nabla^{2}p^{k + 1}\rangle_{i} = \frac{\rho}{n_{0}\varDelta t}\left(\frac{Dn}{Dt}\right)_{i}^{*} \end{equation} | [3] |
| \begin{equation} \nabla^{2}p^{k + 1} = \frac{\rho}{\varDelta t^{2}}\frac{n^{*} - n_{0}}{n_{0}}. \end{equation} | [4] |
| \begin{equation} \langle\nabla p\rangle_{i} = \frac{D_{s}}{n_{0}}\sum_{j \neq i}\frac{p_{j} - \hat{p}_{i}}{|\boldsymbol{{r}}_{ij}|^{2}}\boldsymbol{{r}}_{ij} w(|\boldsymbol{{r}}_{ij}|) \end{equation} | [5] |
| \begin{equation} \hat{p}_{i} = \min_{j \in J}(p_{i},p_{j}),\ J = \{j{:}\ w(|\boldsymbol{{r}}_{ij}|)\neq 0\} \end{equation} | [6] |
| \begin{equation} \boldsymbol{{r}}_{ij} = \boldsymbol{{r}}_{j} - \boldsymbol{{r}}_{i} \end{equation} | [7] |
| \begin{equation} w(r) = \left\{ \begin{array}{lr} \dfrac{r_{e}}{r} - 1 &(0 \leq r < r_{e})\\ 0 &(r_{e} \leq r) \end{array} \right. \end{equation} | [8] |
| \begin{equation} n_{i} = \sum_{i \neq j}w (|\boldsymbol{{r}}_{ij}|). \end{equation} | [9] |
| \begin{equation} \langle \nabla^{2}\phi\rangle_{i} = \frac{2D_{s}}{\lambda n_{0}}\sum_{j \neq i}\frac{\phi_{j} - \phi_{i}}{|\boldsymbol{{r}}_{ij}|^{2}}w(|\boldsymbol{{r}}_{ij}|) \end{equation} | [10] |
| \begin{equation} \lambda = \frac{1}{n_{0}}\sum_{j \neq i}|\boldsymbol{{r}}_{ij}|^{2} w(|\boldsymbol{{r}}_{ij}|). \end{equation} | [11] |
An existence of unphysical pressure fluctuations, which bring a numerical instability, is the most serious drawback of both SPH and MPS methods. An instantaneous distribution of hydrostatic pressure in a rectangular tank calculated by the original MPS method47) is shown in Fig. 1. A hydrostatic pressure includes very strong fluctuations. Although an interval averaging brings a physically sound distribution of pressure, physical pressure peaks are also flattened by an interval averaging, resultantly impact pressure is unable to be estimated. Pressure fluctuation is inevitable in particle methods, which employ moving calculation points. When a particle number density is too high, local pressure goes up to push neighboring particles away.73)
Hydrostatic pressure of water in a rectangular tank calculated by the original MPS method; the initial depth of the tank is 0.2 m filled with the particles with the diameter 4.0 × 10−3 m.
The more exactly the position of particle is updated, the less number density is disturbed. Hence, inter-particle repulsive force to adjust particle number density can be reduced. Theoretically speaking, there exist many of approximations in the local kernel-based interpolations of physical properties on moving calculation points. By improving methods for local kernel-based interpolations, pressure fluctuations can be controlled more effectively. In the next chapter, a digested review of the state-of-the-art of the improvements for attenuating pressure fluctuations in particle methods is presented.
All particle methods suffer from an existence of unphysical pressure fluctuation54),73) which is caused by the Lagrangian nature of moving calculation points, or particles. Hence, a reduction of unphysical pressure fluctuation has been a key issue to improve the performance of particle methods. Methodologies for reducing pressure fluctuation, which also enable us to improve stability in computation, are called accurate particle methods.
In the context of weakly compressible SPH with explicit schemes, Godunov SPH with Riemann solver (Inutsuka, 1994;74) Inutsuka, 200275)) as well as delta-SPH (Antuono et al., 2012)76) schemes have been proposed to enhance an accuracy. Several rigorous studies also investigate the accuracy of SPH with highlighting an importance of higher order interpolation schemes (Le Touzé, 2013).77) Corrective terms for SPH to restore the completeness or consistency of formulations (Randles and Libersky, 1996;78) Chen et al., 1999;79) Oger et al., 2007;80) Fatehi and Manzari, 201181)) as well as momentum conservation (Bonet and Lok, 1999)82) have also been developed and applied to coastal and ocean engineering problems (Sun et al., 2010;83) Xie et al., 201284)).
As for semi-implicit schemes, such as MPS47) and ISPH,71) based on a projection method, differential operator models have been refined for discretizing the source term and Laplacian of PPE as well as for restoring momentum conservation (Khayyer et al., 2008).85) A Higher-order Source term of PPE abbreviated as HS,54) a Higher-order Laplacian (HL) model (Khayyer and Gotoh, 2010;86) Khayyer and Gotoh, 201287)) were proposed for further enhancement of pressure-field accuracy. To attenuate errors of projection onto solenoidal field and resultant errors of volume conservation, the ECS (Error Compensating Source) of PPE was proposed for both MPS and ISPH methods (Khayyer and Gotoh, 2011).88)
Particle methods commonly encounter tensile instability, which occurs in presence of attractive inter-particle forces (Swegle et al., 1994).89) A key point in removal of tensile instability is an accurate estimation of derivatives (Dilts, 1999).90) Actually, differential operator models in the original particle method are not consistent with first-order Taylor series. The Taylor-series consistent pressure gradient model with Gradient Correction (GC)88) plays a significant role in controlling tensile instability.
5.1 CMPS method for momentum conservation.To conserve momentum, interparticle forces must be anti-symmetric, namely equal in magnitude and opposite in direction.82) Since the pressure-gradient force is not anti-symmetric, linear momentum is not preserved in the original MPS method. By considering an auxiliary calculation point on the midpoint of the position vector of particle i and its neighboring particle j, an anti-symmetric and co-linear pressure gradient model:
| \begin{equation} \langle \nabla p \rangle_{i} = \frac{D_{s}}{n_{0}}\sum_{j \neq i}\frac{(p_{i} + p_{j}) - (\hat{p}_{i} + \hat{p}_{j})}{|\boldsymbol{{r}}_{ij}|^{2}}\boldsymbol{{r}}_{ij} w(|\boldsymbol{{r}}_{ij}|) \end{equation} | [12] |
In ISPH method,71) which is the SPH for incompressible fluid, the algorithm is the same with MPS method. Linear momentum is exactly conserved in ISPH method, because both pressure-gradient force and viscous force are anti-symmetric. However, due to anisotropic nature of viscous stresses, viscous forces are not co-linear, thus, angular momentum is not conserved. Corrected ISPH (CISPH) method (Khayyer et al., 2008)85) was proposed based on the variational approach82) to derive a corrective matrix for the viscous stress to guarantee the invariance of potential energy with respect to rigid body rotation. CISPH method ensures the conservation of linear and angular momentums.
5.2 Higher-order Source (HS) in PPE.In the original MPS method, the time differentiation of the particle number density (Dn/Dt) in the PPE is discretized in the first-order forward-difference as shown in Eq. [4]. The ISPH method treats the time differentiation of the particle density (Dρ/Dt) in the PPE in the same manner. Khayyer and Gotoh (2009)54) derived the higher-order time-differentiation of the particle number density with revisiting the definition of the particle number density (Eq. [9]).
| \begin{equation} \left(\frac{Dn}{Dt}\right)_{i}^{*} {}={} \sum_{j \neq i}\left(\frac{\partial w_{ij}}{\partial r_{ij}}\right)\frac{\boldsymbol{{r}}_{ij}^{*} \cdot \boldsymbol{{u}}_{ij}^{*}}{|\boldsymbol{{r}}_{ij}|} \end{equation} | [13] |
| \begin{equation} \langle\nabla^{2}p^{k + 1}\rangle_{i} = -\frac{\rho}{n_{0}\varDelta t}\sum_{i \neq j}\frac{r_{e}}{|\boldsymbol{{r}}_{ij}|^{3}}\boldsymbol{{r}}_{ij}^{*} \cdot \boldsymbol{{u}}_{ij}^{*}. \end{equation} | [14] |
Laplacian must be calculated in both the viscous term of Navier-Stokes equation and the pressure term of PPE. As mentioned above, in the standard MPS method, Laplacian is modeled as exchanges of properties in neighboring particles, which indicates diffusion (Eq. [10]). Khayyer and Gotoh (2010)86) derived Higher order Laplacian (HL) for the MPS method by considering mathematical definition of Laplacian, as a divergence of gradient. Taking the divergence of SPH gradient model (Monaghan, 1992),91) Laplacian is written as:
| \begin{align} \langle\nabla^{2}\phi\rangle_{i} &= \frac{1}{n_{0}}\sum_{j \neq i}\biggl\{\frac{\partial\phi_{ij}}{\partial r_{ij}}\frac{\partial w_{ij}}{\partial r_{ij}} \\ & \quad+ \phi_{ij}\biggl(\frac{\partial^{2}w_{ij}}{\partial r_{ij}^{2}} + \frac{D_{s} - 1}{r_{ij}}\frac{\partial w_{ij}}{\partial r_{ij}}\biggr)\biggr\}. \end{align} | [15] |
| \begin{equation} \langle \nabla^{2}\phi\rangle_{i} = \frac{1}{n_{0}}\sum_{j \neq i}\frac{(5 - D_{s})r_{e}}{|\boldsymbol{{r}}_{ij}|^{3}}\phi_{ij}. \end{equation} | [16] |
In improved particle method, an accuracy of the source term in PPE plays a key role. A semi-implicit scheme can be applied on the mathematical ground of Helmholtz decomposition of a vector field, where an intermediate velocity field, u*, is composed of the divergence-free velocity field and the gradient of a scalar field.
A two-step prediction-correction solution process in the semi-implicit schemes is brought by this vector decomposition. In the prediction step, an intermediate velocity field u* is obtained explicitly by considering viscous and gravitational forces. At this step, the incompressibility of the fluid is not satisfied, i.e., the divergence of u* would not be equal to zero. Hence, a new corrective velocity field, u**, is computed to project u* onto a solenoidal, or divergence-free, field so that the final velocity field, uk+1, would be solenoidal. While, due to numerical errors, the updated velocity field uk+1 slightly deviates from the solenoidal field.
To minimize the deviation from the solenoidal field, two error-compensating terms should be introduced into the source term of the PPE. The first term is related to the instantaneous time variation of the particle number density at the time step k or the divergence of the velocity at this time step, both of which should be zero theoretically. The second term reflects the deviation of particle number density at the time step k from its reference n0, or the time rate of overall volumetric change at the time step k, i.e., it accounts for the accumulative error in the particle number density.
| \begin{equation} \left. \begin{array}{c} \langle \nabla^{2}p^{k + 1}\rangle_{i} = \dfrac{\rho}{\varDelta t}\left\{\dfrac{1}{n_{0}}\left(\dfrac{Dn}{Dt}\right)_{i}^{*} {}+{} S_{\textit{ECS}}\right\}\\ S_{\textit{ECS}} = \textit{func}\left[\dfrac{1}{n_{0}}\left(\dfrac{Dn}{Dt}\right)_{i}^{k}, \dfrac{1}{\varDelta t}\dfrac{n^{k} - n_{0}}{n_{0}}\right] \end{array} \right\} \end{equation} | [17] |
| \begin{equation} S_{\textit{ECS}} = \alpha\left[\frac{1}{n_{0}}\frac{n^{k} - n^{k - 1}}{n_{0}}\right] + \beta\left[\frac{1}{\varDelta t}\frac{n^{k} - n_{0}}{n_{0}}\right] \end{equation} | [18] |
To remove an inconvenience of calibration of these unknown constants, Khayyer and Gotoh (2011)88) modified the multi-term source for the PPE with describing constants α and β as functions of the instantaneous flow field.
| \begin{equation} \left. \begin{array}{c} S_{\textit{ECS}} = \alpha\left[\dfrac{1}{n_{0}}\left(\dfrac{Dn}{Dt}\right)_{i}^{k}\right] + \beta\left[\dfrac{1}{\varDelta t}\dfrac{n^{k} - n_{0}}{n_{0}}\right]\\ \alpha = \left|\dfrac{n^{k} - n_{0}}{n_{0}}\right|;\ \beta = \left|\dfrac{\varDelta t}{n_{0}}\left(\dfrac{Dn}{Dt}\right)_{i}^{k}\right| \end{array} \right\} \end{equation} | [19] |
The instantaneous distributions of hydrostatic pressure of water in a rectangular tank calculated by the original MPS method and CMPS-HS-HL-ECS method, are shown in Fig. 2. Although the original MPS method provides unphysical and unstable instantaneous pressure field, physically sound and stable profile of pressure is found in the result of CMPS-HS-HL-ECS method.
Hydrostatic pressure of water in the rectangular tank calculated by the CMPS-HS-HL-ECS method in comparison to the original MPS method.
In Fig. 3, wave breaking and run-up on a uniform slope calculated by the original MPS method and CMPS-HS-HL-ECS method, are shown. Only on the free surface, the dark-colored zone, which indicates a zero pressure, should exist. But in the result of the original MPS method, many of dark colored particles, which indicate zero pressure, exist in water. On the other hand, CMPS-HS-HL-ECS method provides physically sound profile of pressure. The time variation of pressure at a fixed measuring point is shown in Fig. 4. Less fluctuated and more physically sound time-variation curve is provided by CMPS-HS-HL-ECS method.
Numerical wave flume by the original MPS and the CMPS-HS-HL-ECS methods for simulating wave run-up on uniform slope.
Comparison of the time series of pressure calculated by the original MPS and the CMPS-HS-HL-ECS methods at the toe of the slope.
A tensile instability, which occurs in presence of attractive inter-particle forces, is a common knotty problem of particle methods (Swegle et al., 1994).89) As shown by Dilts (1999),90) accurate estimation of derivatives contributes to remove the tensile instability. Actually, differential operator models in the original MPS method do not satisfy consistency with first-order Taylor series.
The Taylor-series consistent pressure gradient model with Gradient Correction (GC)88) is expressed as:
| \begin{equation} \left\langle\frac{\nabla p}{\rho}\right\rangle_{i} {}={} \frac{D_{s}}{\rho n_{0}}\sum_{j \neq i}\frac{p_{j} - p_{i}}{|\boldsymbol{{r}}_{ij}|^{2}}\boldsymbol{{C}}_{i} \cdot \boldsymbol{{r}}_{ij}w_{ij} \end{equation} | [20] |
| \begin{equation} \boldsymbol{{C}}_{i} = \frac{1}{D_{s}}\left(V_{i}\sum_{j \neq i}\frac{\boldsymbol{{r}}_{ij} \otimes \boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|^{2}}w_{ij}\right)^{-1} \end{equation} | [21] |
| \begin{equation} V_{i} = \frac{1}{\displaystyle\sum_{j \neq i}w_{ij}} \cong \frac{1}{n_{0}}. \end{equation} | [22] |
Evolution process of a square patch of 2-dimensional inviscid fluid subjected to a rigid rotation, in which the initial square side is 1.0 m in length and angular velocity around the center of the square is 10.0 s−1.
Taylor series consistency also plays a key role in computation of density in the vicinity of the interface of a multi-phase flow. At the liquid-gas interface, there exist a mathematical discontinuity of the density, which is the source of numerical instability. In ordinary particle methods, the density is discretized in zero-order accuracy, which brings non-physical density diffusion. The description of density to satisfy the consistency with first-order Taylor series enables a violent sloshing simulation of liquid-gas two-phase flow, the density ratio of which is 1000 (Khayyer and Gotoh, 2013).94)
Consideration of Taylor series consistency also reveals an existence of artificial repulsive force in gradient model of MPS method (Tsuruta et al., 2013).95) The artificial repulsive force as well as the artificial viscosity in SPH method is essential for computational stability, but it should be minimized due to its unphysical nature. In SPH method with an explicit algorithm, Godunov SPH with Riemann solver74),75) and Balsara switch (Balsara, 1995)96) work effectively in minimizing artificial viscosity. On the other hand, in MPS method and other projection-based schemes, DS (Dynamic Stabilization) scheme95) provides a smart algorithm to find minimum artificial repulsive force.
Boundary condition is an unsettled issue of accurate particle method. Considering a single-phase flow, in particle method, free-surface particles are regarded as instantaneous fixed wall boundary to which dynamic, or Dirichlet, condition of pressure is given.45)–47) Although this boundary treatment is quite simple, it disregards volume conservation in the vicinity of free surface. Treating the free surface as the fixed wall in PPE prevents fluid-particles from flowing into a void space, because there is no calculation point to define a void. Due to this lack of information of existing voids, unphysical voids generated erroneously cannot be removed efficiently. In SPH method, so-called ghost (Colagrossi and Landrini, 2003)70) or mirror (Basa et al., 2009)97) particles as fictitious neighboring particles to complete the truncated kernel supports at boundaries, can be applied to improve a conservation of volume. But physical definitions of ghost and mirror particles are not clear. In MPS method, SPP (Space Potential Particle) (Tsuruta et al., 2015),98) which is placed at the centroid of void of each boundary particles to give an exact Dirichlet condition at the boundary, was proposed.
6.1 Dynamic stabilization.Since an overlapping of adjacent particles brings a failure of computation, in particle method, an inter-particle force should be repulsive. An artificial repulsive force is essential to prevent particles from overlapping each other.
Pressure gradient model consistent with 0th-order Taylor series is written as follows:
| \begin{equation} \langle \nabla p \rangle_{i} = \frac{D_{s}}{n_{0}}\sum_{j \neq i}\frac{p_{j} - p_{i}}{|\boldsymbol{{r}}_{ij}|}\frac{\boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|}w_{ij}. \end{equation} | [23] |
| \begin{equation} \left. \begin{array}{c} \langle \nabla p \rangle_{i} = \dfrac{D_{s}}{n_{0}}\displaystyle\sum_{j \neq i}\dfrac{p_{j} - p_{i}}{|\boldsymbol{{r}}_{ij}|}\dfrac{\boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|}w_{ij} + \boldsymbol{{F}}_{i}{}^{\textit{AR}}\\ \boldsymbol{{F}}_{i}{}^{\textit{AR}} = \dfrac{D_{s}}{n_{0}}\displaystyle\sum_{j \neq i}\dfrac{p_{i} - \hat{p}_{i}}{|\boldsymbol{{r}}_{ij}|}\dfrac{\boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|}w_{ij} \end{array} \right\}. \end{equation} | [24] |
| \begin{equation} \boldsymbol{{F}}_{i}{}^{\textit{AR}} = \frac{D_{s}}{n_{0}}\sum_{j \neq i}\dfrac{(p_{i} - \hat{p}_{i}) + (p_{i} - \hat{p}_{j})}{|\boldsymbol{{r}}_{ij}|}\frac{\boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|}w_{ij}. \end{equation} | [25] |
In order to find minimum-required artificial repulsive force, Tsuruta et al. (2013)95) proposed Dynamic Stabilization (DS) scheme, which defines the minimum-required artificial repulsive force as the force to push back the overlapped particles into the state of contacting each other.
The DS scheme provides the stabilizing force FDS in the followings equations:
| \begin{equation} \boldsymbol{{F}}_{i}{}^{\textit{AR}} = \frac{1}{n_{0}}\sum_{j \neq i}\boldsymbol{{F}}_{ij}{}^{\textit{DS}}w_{ij} \end{equation} | [26] |
| \begin{equation} \begin{cases} \boldsymbol{{F}}_{ji}{}^{\textit{DS}} = 0 & \text{($|\boldsymbol{{r}}_{ij}^{*}| \geq d_{ij}$)}\\ \boldsymbol{{F}}_{ij}{}^{\textit{DS}} = -\rho_{i}\varPi_{ij}\boldsymbol{{e}}_{ij\|} & \text{($|\boldsymbol{{r}}_{ij}^{*}| < d_{ij}$)} \end{cases} \end{equation} | [27] |
| \begin{equation} d_{ij} = \alpha_{\textit{DS}}\frac{d_{i} + d_{j}}{2};\quad\boldsymbol{{e}}_{ij\|} = \frac{\boldsymbol{{r}}_{ij}}{|\boldsymbol{{r}}_{ij}|};\quad |\boldsymbol{{e}}_{ij\|}| = 1 \end{equation} | [28] |
| \begin{equation} \alpha_{\textit{DS}} + \alpha_{\textit{dt}} = 1 \end{equation} | [29] |
| \begin{equation} \varPi_{ij} = \frac{\rho_{j}}{\varDelta t^{2}(\rho_{i} + \rho_{j})}(\sqrt{d_{ij}^{2} - |\boldsymbol{{r}}_{ij\bot}^{*}|^{2}} - |\boldsymbol{{r}}_{ij\|}^{*}|) \end{equation} | [30] |
Concept of Dynamic Stabilization (DS) scheme.
Although the GC scheme can calculate tensile force, or negative pressure, it tends to be unstable in the presence of violent motion of a free surface. To improve a stability of GC scheme, DS scheme, which provides minimum-required artificial repulsive force, is effective. GC-based pressure gradient stabilized by DS scheme is written as:
| \begin{align} \left\langle \frac{\nabla p}{\rho}\right\rangle_{i} &= \frac{D_{s}}{\rho n_{0}}\sum_{j \neq i}\frac{p_{j} - p_{i}}{|\boldsymbol{{r}}_{ij}|^{2}}\boldsymbol{{C}}_{i} \cdot \boldsymbol{{r}}_{ij}w_{ij} \\ &\quad+ \frac{1}{\rho n_{0}}\sum_{j \neq i}\boldsymbol{{F}}_{ij}{}^{\textit{DS}}w_{ij}. \end{align} | [31] |
Settlement process of heavy particles, the density of which is 2.65 × 103 kg/m3, in water; the tank is 0.25 m in length and 0.2 m in height, and is filled with 7000 fluid particles with the diameter 2.5 × 10−3 m.
On a free surface, a simple and robust but ad-hoc boundary condition is applied in particle method. Free-surface particles are treated as instantaneous fixed wall boundary, on which the dynamic, or Dirichlet, condition of pressure is given. Although this boundary treatment is quite simple and robust, but brings a problem about disregard of the volume conservation. Since the free-surface boundary is set as a fixed wall in the PPE, it is premised that particles (or fluids) do not flow into a void space. Such a calculation process without any information of existence of voids misses a potential of voids, consequently unphysical voids cannot be removed efficiently.
Space Potential Particle (SPP) was proposed by Tsuruta et al. (2015)98) to improve these problems in free-surface boundary condition. The SPP describes the existence of void space as ‘space potential’ by introducing virtual free-surface boundary particle additionally. This technique achieves that all the fluid particles can be handled as pure fluids with consistency with mass conservation.
Accurate estimation of differential operator requires kernel support of a targeted particle to be filled-up by neighboring particles. However, particles in the vicinity of a free surface do not have sufficient neighboring particles around them. Here, the kernel support of such a free-surface particle is interpreted to be filled-up with substantial particles, consequently non-substantial void space remains in their kernel support.
The SPP for the description of the void space is set for all the particles around free surface one by one. The number density of particle represented by SPP of the particle i, wispp is given as:
| \begin{equation} w_{\textit{ispp}} = \begin{cases} n_{0} - n_{i} & \text{($n_{0} > n_{i}$)}\\ 0 & \text{($n_{0} \leq n_{i}$)} \end{cases}. \end{equation} | [32] |
| \begin{equation} |\boldsymbol{{r}}_{\textit{ispp}} - \boldsymbol{{r}}_{i}| = \frac{r_{e}}{w_{\textit{ispp}} + 1} = \frac{r_{e}}{n_{0} - n_{i} + 1} \end{equation} | [33] |
| \begin{equation} \boldsymbol{{r}}_{\textit{ispp}} = \boldsymbol{{r}}_{i} - |\boldsymbol{{r}}_{\textit{ispp}} - \boldsymbol{{r}}_{i}|\frac{\boldsymbol{{r}}_{gi}}{|\boldsymbol{{r}}_{gi}|} \end{equation} | [34] |
| \begin{equation} \boldsymbol{{r}}_{\textit{gi}} = \frac{1}{n_{0}}\sum_{\substack{j \neq i\ \textit{and}\\ j \neq \textit{ispp}}}\boldsymbol{{r}}_{j} w_{ij} \end{equation} | [35] |
| \begin{equation} \boldsymbol{{u}}_{\textit{ispp}} = \boldsymbol{{u}}_{i} \end{equation} | [36] |
In high Reynolds number flows, all of the particle methods suffer from particle clumping and resultant void formation in negative pressure zone (Adami et al., 2013).99) At the downstream of a bluff body, streamlines separate from the body and a negative pressure zone is formed. Kármán vortex generated in the flow downstream of cylindrical column is well-known benchmark test of computational stability in existence of negative pressure. Kármán vortex simulated by GC scheme combined with SPP is shown in Fig. 8. Absence of SPP leads to a formation of unphysical void downstream of the cylinder, but the SPP contributes well to fill a void and adjust positions of fluid particles appropriately. Alternating vortices are well simulated by the SPP.
Performance of SPP in Kármán vortex simulation under Re = 1200, with time resolution Δt = 2.5 × 10−3 s.
In Fig. 9, relationship of the accurate particle methods is shown schematically. There are two major items for improving accuracy: the gradient model and the Poisson pressure equation (PPE). CMPS method ensures momentum conservation of a gradient model, while GC scheme provides Taylor-series consistent gradient model, whose computational stability is enhanced by DS scheme. HL scheme enables accurate estimation of Laplacian in PPE, HS and ECS schemes improve the accuracy of the source term of PPE. SPP scheme provides improved boundary condition of PPE.
Relationship of the accurate particle methods.
Figure 10 shows snapshots of the dam breaking simulation by the accurate particle methods. The standard MPS [M0] shows significantly irregular distribution of pressure with severe noises, on the other hand, CMPS [M1] moderately suppresses such noises. Then, HS scheme [M2] drastically improves the reproducibility of the pressure. HL and ECS schemes [M3, M4] further gain the reproducibility with enhancing the smoothness of free-surface line. At t = 0.75 s, GC scheme [M5] without stabilizer shows unphysical perturbation of the free-surface in the splash region and at the position x < 0.7 m. However, it leads to reproduction of negative pressure. Then, the free-surface line is further smoothed by introducing DS scheme [M6]. Comparing with the repulsion-based model corresponding to CMPS-HS-HL-ECS method [M4], the pressure noise is more effectively suppressed by the Taylor-series consistent model, namely MPS-GC series, which can be easily recognized from the snapshots at t = 0.28 s around the bottom left corner of the calculation domains. In the snapshots with SPP scheme [M7] at t = 0.75 s, it is shown that the particles around the free surface change their motions characterized by more continuity resulting in reduction of unphysical particle dispersiveness seen in other simulations.
Dam breaking simulation by the accurate particle methods.
The required computational time in these benchmarks is shown in Fig. 11. Each computational time is normalized by that of the standard MPS method [M0]. By introducing accurate schemes, the repulsion-based scheme, namely CMSP series, gradually increases the computational time. However, even CMPS-HS-HL-ECS scheme brings an increment of less than 20% of computational load. Taylor-series consistent schemes, or MPS-GC series, save the computational time by more than 20%. Although MPS-GC scheme has a problem of numerical instability in simulations of the violent flow, DS scheme enhances its numerical stability. Increments of the computational load with the introduction of DS and SPP schemes are quite small. The set of accurate schemes showing the best performance of dam breaking simulation in Fig. 10, namely MPS-HS-HL-ECS-GC-DS-SPP scheme [M7] saves the computational time by more than 20% in comparison to the standard MPS method.
Required computational time in dam breaking simulations.
Here, the Numerical Wave Flume, which is the key tool for an innovative design of coastal structures, has been introduced.
First, to show a necessity of NWF, current structurally-resistive design against wave action was overviewed. Around a structure, there exists a violent flow with wave breaking, and empirical design formulae have a limit in accuracy. CFD technologies including the NWF are hoped to be a core tool of the design of coastal structures. Secondly, techniques of the interface capturing (MAC, VOF, C-CUP) were shortly reviewed, as key items of the NWF. In all of the grid method, a numerical diffusion, which is caused by the non-linear advection term in Navier-Stokes equation, makes interface capturing difficult. Thirdly, the role of NWF in computational coastal hydraulics was discussed. The surf-zone dynamics to simulate mass-and-momentum transport in a surf zone, or multi-phase turbulent flow model for sediment transport in a surf zone, was overviewed.
The particle method, which is free form numerical diffusion due to the absence of the advection term in Navier-Stokes equation, shows robustness in computation of the violently-deforming liquid-gas interface including a topological change of the interface like a wave breaking. The latter part of this article described recent improvements of the particle method as core methods of the NWF. Accurate particle methods have been developed for attenuating unphysical pressure fluctuation, which is the common nature of Lagrangian methods. Key methods for improving accuracy, such as CMPS method for momentum conservation, Higher-order Source (HS) of Poisson Pressure Equation (PPE), Higher-order Laplacian (HL), Error-Compensating Source (ECS) in PPE, and Gradient Correction (GC) for ensuring Taylor-series consistency were reviewed briefly. Finally, Dynamic Stabilization (DS) for providing minimum-required artificial repulsive force to improve stability of computation and Space Potential Particle (SPP) for describing exact free-surface boundary condition were described as the latest new frontier of the particle method.
The NWF is expected to contribute to further more subjects. Bubble generation process by wave breaking is a conundrum from the scientific viewpoint. A plunging jet hits water surface to inject bulk air into water. But the process of fragmentation of bulk air into bubbles may contain a series of complex phenomena.2) A microscopic scale modeling of liquid-gas interface, including turbulence due to bubble motion, may bring some developments on these complex problems. In the design of the coastal structure against Tsunami in the Level-2, whose return period is 1000 years, a multi-physics modeling with computing interaction of flow, structure and soil-foundation draws attention. An overflow brings jets hitting sea bottom at the onshore-side of the structure, and generates scour hole. These complex processes include multi-physics problems. In both of grid method and particle method, the modeling of elastic and plastic bodies, which are essential for the multi-physics modeling, progresses rapidly.100)–102) A role of the CFD in the design of coastal structures will become more and more significant.
Hitoshi Gotoh was born in 1963 and graduated from the Faculty of Engineering at Kyoto University, specializing in civil engineering. He continued his study at the Graduate School of Engineering and obtained master’s and doctoral degrees under the supervision of Prof. Hiroji Nakagawa. He began his professional career as Research Associate in 1992 at the Prof. Tetsuo Sakai’s Laboratory at Kyoto University, where he started research in coastal engineering. Since then he has experienced Lecturer (1996) and Associate Professor (1997 to 2008) and Professor (2008 to date) at Kyoto University. In the meanwhile, from 1998 to 1999, he stayed at Technical University of Denmark as a visiting scientist. His major is coastal engineering, especially wave dynamics and sediment transport mechanics. He has done pioneering works on a multiphase-flow modeling of sediment transport. His recent research topic is concentrated on a development of computational methods for free-surface flow, e.g., particle method. He has developed many of methods to enhance an accuracy of Lagrangian computations of free-surface flow, some of which play key roles in the numerical wave flume. He won Coastal Engineering Journal Awards of 2005, 2009 and 2012.
Akio Okayasu was born in 1961 and graduated from the Faculty of Engineering at the University of Tokyo in 1984, specializing in civil engineering. He continued his study of coastal engineering at the Graduate School of Engineering and obtained master’s degree in 1986 under the supervision of Prof. Kiyoshi Horikawa. Then he received Doctor of Engineering from the University of Tokyo under the supervision of Prof. Akira Watanabe in 1989. He began his professional career as Research Associate in 1989 at Prof. Tomoya Shibayama’s Laboratory at Yokohama National University. He then became Associate Professor in 1996 and moved to Tokyo University of Fishery in 2001. In the meanwhile, he stayed at Delaware University of the United States of America in 1993 as Visiting Researcher and at Technical University of Denmark as Research Associate Professor in 2001. Due to the integration of universities in 2003, Tokyo University of Fisheries changed its name to Tokyo University of Marine Science and Technology where he was promoted to Professor in 2005. He was elected to be Dean of Graduate School of Marine Science and Technology at Tokyo University of Marine Science and Technology in 2012 and Dean of School of Marine Environment and Resources in 2017. His major is coastal engineering, especially surfzone dynamics and sediment transport. His recent research topics include tsunami disaster prevention and mitigation. He won JAMSTEC Nakanishi Award in 2013 and CEJ Best Citation Award in 2015.
We are particularly grateful to Prof. Kiyoshi Horikawa, M.J.A., who kindly invited and encouraged us to write this review article in the Proceedings of the Japan Academy, Series B.
