Printer Friendly

Computational wave dynamics for innovative design of coastal structures.

1. Introduction

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.

2. Numerical wave flume by grid methods

2.1 Solver of Navier-Stokes equation.

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:

[D[rho]/Dt] + [rho][nabla] x u = 0 [1]

Du/Dt = -[1/[rho]] -[nabla]p + v[[nabla].sup.2]u + g [2]

where u is velocity vector, t is time, [rho] is fluid density, p is pressure, g is gravitational acceleration vector and v is laminar kinematic viscosity.

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 method (8,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., 2001 (13)) 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-[epsilon] 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, 1996 (19)), 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., Schaffer 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.

3. Calculation of sediment transport with numerical wave flume

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-[epsilon] 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., 20 05 (28)).

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, 20 00). (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, 2013 (40)).

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., 2016 (42)). 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) 2016 (44)).

4. Numerical wave flume by particle method

4.1 particle method. 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, 1977 (46)) 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, 2014 (50)), wave overtopping (Gotoh et al., 2005; (51) Shao et al., 2006 (52)), wave run-up (Shadloo et al., 2015), (53) wave impact (Khayyer and Gotoh, 2009; (54) Lee et al., 2011; (55) Altomare et al., 2015 (56)), wave-induced nearshore current (Farahani et al., 2014), (57) violent sloshing (Delorme et al., 2009; (58) Gotoh et al., 2014 (59)), green water on ships (Shibata and Koshizuka, 2007; (60) Le Touze et al., 2010 (61)), 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., 2015 (69)).

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:

[mathematical expression not reproducible] [3]

where n is particle number density, [n.sub.0] is the reference particle number density as a criteria of the fluid volume, which is defined by the initial particle number density in a regular distribution of particles, and [DELTA]t is the time increment of computation. Subscripts k, i denote calculation time step and particle number, respectively. Superscript * refers to pseudo time step k + [1/2]. In the original MPS method, (47) PPE is given as:

[mathematical expression not reproducible]. [4]

In the original MPS method, to ensure the stability of computation, the pressure gradient is defined by replacing [p.sub.i] with the minimum pressure value in the kernel support domain:

[mathematical expression not reproducible] [5]

[mathematical expression not reproducible] [6]

[r.sub.ij] = [r.sub.j] - [r.sub.i] [7]

where [D.sub.s] is the number of the spatial dimension, [r.sub.i] the position vector of particle i and w the kernel function. In the original MPS method, to avoid particle clustering, kernel function which is infinite at r = 0 is used.

[mathematical expression not reproducible] [8]

where [r.sub.e] is the radius of kernel support. The particle number density is defined as follows:

[n.sub.i] = [summation over (i[not equal to]j)]w([absolute value of [r.sub.ij]]). [9]

In the original MPS method, Laplacian is modeled to provide exchanges of properties in neighboring particles, which indicates diffusion.

[mathematical expression not reproducible] [10]

[mathematical expression not reproducible]. [11]

In SPH method, differential operations are defined as a differentiation of the kernel function, hence vector differential operators are mathematically consistent. But to ensure the accuracy of differential operations, a higher-order kernel function, which needs rather heavier computational load, must be introduced. While, the MPS method has a specific model for each differential operator, to provide stability with less computational load.

4.3 Pressure fluctuation. 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 method (47) 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)

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.

5. Accurate particle method

All particle methods suffer from an existence of unphysical pressure fluctuation (54,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, 2002 (75)) 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 Touze, 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, 2011 (81)) 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., 2012 (84)).

As for semi-implicit schemes, such as MPS (47) 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, 20 08). (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, 2012 (87)) 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:

[mathematical expression not reproducible] [12]

can be derived. (49) This formulation is called Corrected MPS (CMPS) method.

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 approach (82) 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[rho]/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]).

[mathematical expression not reproducible] [13]

where [w.sub.ij] is the kernel function, [r.sub.ij], [u.sub.ij] are the relative position and velocity vectors of the particle j to the particle i, respectively. The superscript * denotes that the physical quantities are calculated at the prediction step. By applying standard kernel of MPS method (Eq. [8]), the PPE with higher order source terms for the MPS method is given as follows:

[mathematical expression not reproducible] [14]

The CMPS methods with Higher order Source terms was named CMPS-HS method.

5.3 Higher-order Laplacian (HL). 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:

[mathematical expression not reproducible] [15]

By applying the standard MPS kernel (Eq. [8]), above equation is rewritten as:

[mathematical expression not reproducible] [16]

Replacing the SPH gradient model by Taylor-series consistent gradient model, further-accurate Laplacian model, namely Corrected HL (CHL) is derived. (92)

5.4 Error-Compensating Source terms (ECS) for PPE. 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.sup.*], 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.sup.*] 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.sup.*] would not be equal to zero. Hence, a new corrective velocity field, [u.sup.**], is computed to project [u.sup.*] onto a solenoidal, or divergence-free, field so that the final velocity field, [u.sup.k+1], would be solenoidal. While, due to numerical errors, the updated velocity field [u.sup.k+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 [n.sub.0], 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.

[mathematical expression not reproducible] [17]

Kondo and Koshizuka (2011) (93) suggested a multi-term source for the PPE to obtain smooth time variation of the particle number density while keeping it invariant.

[mathematical expression not reproducible] [18]

where, [alpha] and [beta] are unknown constants, that require appropriate calibration.

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 [alpha] and [beta] as functions of the instantaneous flow field.

[mathematical expression not reproducible] [19]

Certainly, for particles at-and-close to the free surface, the initial particle number density would be much smaller than its reference [n.sub.0], which is calculated for homogeneously spaced particles. In violent flows computation, [n.sub.0] should be updated frequently, e.g., every 20 time steps, in the vicinity of free surface. (88)

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.

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.

5.5 Taylor series consistent scheme. 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:

[mathematical expression not reproducible] [20]

where the Gradient Correction matrix, [C.sub.i], which is given as:

[mathematical expression not reproducible] [21]

[V.sub.i] = 1/[summation over (j[not equal to]i)][w.sub.ij] [congruent to] 1/[n.sub.0]. [22]

Figure 5 shows the evolution process of a square patch of fluid subjected to a rigid rotation, which is one of the typical benchmark problems for tensile instability. In a short duration after the beginning of the rotation, MPS method shows a spurious pressure field with dispersed particles at the patch corners. The result of MPS-HS shows that the higher order source term improves stability. A combined use of the ECS term further improves the stability in some degree, however, the MPS-HS-ECS simulation failed at t = 0.20 s. On the other hand, an introduction of gradient correction, namely MPS-HS-GC, provides a more accurate and stable reproduction of the fluid patch at t = 0.20 s. The best snapshot is provided by MPS-HS-HL-ECS-GC method. These comparisons clarify that the gradient correction plays a significant role in controlling tensile instability. Under the absence of tensile stress, or negative pressure, CMPS-HS-HL-ECS method ensures accurate and stable computation. The MPS-GC, instead of the CMPS, enables accurate computation in existence of tensile stress.

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 nonphysical 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)

6. The latest new frontier of the accurate particle method

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 solver (74,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) scheme (95) 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:

[mathematical expression not reproducible]. [23]

Pressure gradient model of original MPS method, namely Eq. [5], can be re-written with separating the 0th-order-Taylor-series consistent term as follows:

[mathematical expression not reproducible]. [24]

The first term corresponds to Eq. [23], and the second term corresponds to the artificial repulsive force. As for the CMPS method, the artificial repulsive force in pressure gradient model is given as follows:

[mathematical expression not reproducible]. [25]

All of the pressure gradient models in semi-implicit schemes must include artificial repulsive force in addition to the Taylor-series consistent term. Needless to say, the artificial repulsive force should be minimized, similar to an artificial viscosity, which is inevitable in stabilization of the WCSPH method.

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 [F.sup.DS] in the followings equations:

[F.sub.i.sup.AR] = [1/[n.sub.0]][summation over (j[not equal to]i)][F.sub.ij.sup.DS][w.sub.ij] [26]

[mathematical expression not reproducible] [27]

[mathematical expression not reproducible] [28]

[[alpha].sub.DS] + [[alpha].sub.dt] = 1 [29]

where [e.sub.ij[parallel]] is unit vector of [r.sub.ij], [d.sub.i] is diameter of particle i. The constant [[alpha].sub.DS] is for adjusting active range of [F.sub.ij.sup.DS] according to the Courant stability condition for the time resolution. The constant [[alpha].sub.dt] is the ratio of time increment to Courant number. The superscript * refers to a state after the advection by the Taylor series consistent gradient term, Eq. [20]. Through a geometrical consideration, the parameter [[PI].sub.ij] for describing a magnitude of [F.sub.ij.sup.DS] is written in the following form.

[mathematical expression not reproducible] [30]

where [DELTA]t is time increment, [r.sup.*.sub.ij[parallel]] is the parallel vector of [r.sup.*], and [r.sup.*.sub.ij[perpendicular to]] is the normal vector of [r.sup.*]. The concept of DS is shown in Fig. 6.

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:

[mathematical expression not reproducible]. [31]

Figure 7 shows the settlement process of heavy particles, the density of which is 2.65 x [10.sup.3] kg/[m.sup.3], in water. The CMPS--HS method reproduces a stable free surface, however, velocity vectors are unphysical with perturbations. In addition, some heavy particles, which are fully bounded within the light fluid particles, do not settle down. Excessive repulsive forces may cause these unphysical results. On the other hand, MPS-HS-GC-DS method portrays a physically sound reproduction of the phenomenon. A stable free surface and a circulating flow driven by the settlement of heavy particle are reproduced. The settlement of heavy particles is well simulated. The DS scheme also contributes to attenuate unphysical void generation, which is one of the common crucial drawbacks of particle method. The purely repulsive interaction works so as to get rid of unphysical voids. Being pushed by the opposite-side particle of the void due to repulsive forces, adjacent particles to the void likely move towards the void.

6.2 Space Potential Particle. 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, [w.sub.ispp] is given as:

[mathematical expression not reproducible]. [32]

Considering the maldistribution of particles due to the void, the SPP is set at the opposite side of centroid of neighboring fluid particles in the support domain. Referring standard kernel function of MPS method (Eq. [8]), the distance between target particle i and its SPP is formulated as follows:

[mathematical expression not reproducible] [33]

where [r.sub.ispp] is position vector of SPP. As SPP exists in the opposite side of centroid of neighboring fluid particles, and the target particle should be at the centroid of particles in the kernel support, the position vector of SPP is written as follows:

[r.sub.ispp] = [r.sub.i] - [absolute value of [r.sub.ispp] - [r.sub.i]] []/[absolute value of []] [34]

[mathematical expression not reproducible] [35]

where [] is position vector of the centroid of neighboring particles of particle i. In accordance with the free-surface boundary condition, SPP velocity is defined not to give any shear stress to the targeted particle i as:

[u.sub.ispp] = [u.sub.i] [36]

where [u.sub.i] is the velocity of particle i.

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. Karman vortex generated in the flow downstream of cylindrical column is well-known benchmark test of computational stability in existence of negative pressure. Karman 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.

6.3 Benchmarks of accurate particle methods. 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.

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.

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.

7. Concluding remarks

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.

doi: 10.2183/pjab.93.034


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.


(1) Morison, J.R., O'Brien, M.P., Johnson, J.W. and Schaaf, S.A. (1950) The force exerted by surface wave on piles. Petroleum Transactions, American Institute of Mining Engineers 189, 149-154.

(2) Gotoh, H., Okayasu, A. and Watanabe, Y. (2013) Computational Wave Dynamics. World Scientific, Singapore; Hackensack, N.J.

(3) Goda, Y. (1973) A new method of wave pressure calculation for the design of composite breakwater. Report of the Port and Harbour Research Institute 12, 31-70 (in Japanese).

(4) Hudson, R.Y. (1959) Laboratory investigations of rubble-mound breakwaters. Proc. ASCE 85, 93-121.

(5) Goda, Y., Kishira, Y. and Kamiyama, Y. (1975) Laboratory investigation on the overtopping rates of seawalls by irregular waves. Report of the Port and Harbour Research Institute 14, 3-44 (in Japanese).

(6) Hirt, C.W., Amsden, A.A. and Cook, J.L. (1977) An Arbitrary Lagrangian-Eulerian computing method for all flow speeds. J. Comput. Phys. 135, 203-216.

(7) Harlow, F.H. and Welch, J.E. (1965) Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys. Fluids 8, 2182-2189.

(8) Hirt, C.W. and Nichols, B.D. (1981) Volume of fluid VOF method for the dynamics of free boundaries. J. Comput. Phys. 39, 201-225.

(9) Sussman, M., Smereka, P. and Osher, S. (1994) A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146-159.

(10) Nichols, B.D., Hirt, C.W. and Hotchkiss, R.S. (1980) SOLA-VOF--A Solution Algorithm for Transient Fluid with Multiple Free Boundaries. Report LA-8355, Los Alamos Scientific Laboratory, University of California, Los Alamos, N.M.

(11) Yabe, T. and Wang, P.Y. (1991) Unified numerical procedure for compressible and incompressible fluid. J. Phys. Soc. Jpn. 60, 2105-2108.

(12) Yabe, T. and Aoki, T. (1991) A universal solver for hyperbolic equation by cubic-polynomial interpolation, I. One-dimensional solver. Comput. Phys. Commun. 66, 219-232.

(13) Yabe, T., Xiao, F. and Utsumi, T. (2001) The constrained interpolation profile method for multiphase analysis. J. Comput. Phys. 169, 556-593.

(14) Isobe, M., Takahashi, S., Yu, S.P., Sakakiyama, T., Fujima, K., Kawasaki, K., Jiang, Q., Akiyama, M. and Ohyama, H. (1999) Interim development of a numerical wave flume for maritime structure design. Proc. Civil Eng. in the Ocean. 15, 321-326 (in Japanese).

(15) Coastal Development Institute of Technology (2001) Research and Development of Numerical Wave Flume: CADMAS-SURF. Tokyo, Japan (in Japanese).

(16) Boussinesq, J. (1872) Theorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. Journal de Mathematique Pures et Appliquees. Deuxieme Serie 17, 55-108.

(17) Peregrine, D.H. (1967) Long waves on a beach. J. Fluid Mech. 27, 815-827.

(18) Nwogu, O. (1993) Alternative form of Boussinesq equations for nearshore wave propagation. J. Waterw. Port Coast. Ocean Eng. ASCE 119, 618-638.

(19) Beji, S. and Nadaoka, K. (1996) A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth. Ocean Eng. 23, 691-704.

(20) Schaffer, H.A., Deigaard, R. and Madsen, P.A. (1992) A two-dimensional surfzone model based on the Boussinesq equations. Proc. 23rd Int. Conf. Coastal Eng., ASCE, 576-589.

(21) Komar, P.D. and Inman, D.L. (1970) Longshore sand transport on beaches. J. Geophys. Res. 75, 5514-5527.

(22) U.S. Army Corps of Engineers, Coastal Engineering Research Center. (1984) Shore Protection Manual. Dept. of the Army, Waterways Experiment Station, Corps of Engineers, Coastal Engineering Research Center, Vicksburg, Miss.

(23) Okayasu, A., Shibayama, T. and Horikawa, K. (1988) Vertical variation of undertow in the surf zone. Proc. 21st Int. Conf. Coastal Eng., ASCE, 478-491.

(24) Wijayaratna, N. and Okayasu, A. (2000) DNS of wave transformation, breaking and run-up on sloping beds. Proc. 4th Int. Conf. Hydrodynamics, 527-532.

(25) Sunamura, T., Bando, T. and Horikawa, K. (1978) An experimental study of sand transport mechanism and rate over asymmetrical sand ripples. Proc. 25th Japanese Conf. Coastal Eng., 250-254 (in Japanese).

(26) Sato, S., Uehara, H. and Watanabe, A. (1986) Numerical simulation of the oscillatory boundary layer flow over ripples by a k-[epsilon] turbulence model. Coast. Eng. Jpn., JSCE 29, 65-78.

(27) Christensen, E.D. and Deigaard, R. (2001) Large eddy simulation of breaking waves. Coast. Eng. 42, 53-86.

(28) Okayasu, A., Suzuki, T. and Matsubayashi, Y. (2005) Laboratory experiment and three-dimensional Large Eddy Simulation of wave overtopping on gentle slope seawalls. Coast. Eng. J. 47, 71-90.

(29) Suzuki, T., Okayasu, A. and Shibayama, T. (2007) A numerical study of intermittent sediment concentration under breaking waves in the surf zone. Coast. Eng. 54, 433-444.

(30) Nielsen, P. (1992) Coastal Bottom Boundary Layers and Sediment Transport. Advanced Series on Ocean Eng. Vol. 4, World Scientific, Singapore; River Edge, N.J.

(31) Gotoh, H., Tsujimoto, T. and Nakagawa, H. (1995) Refined PSI-cell model for interphase and interparticle momentum transfer in bed-load layer. J. Hydroscience Hydraulic Eng., JSCE 13, 13-24.

(32) Gotoh, H. and Sakai, T. (2000) Behavior of bed-material particles as a granular material in a bed-load transport process. J. Hydroscience Hydraulic Eng., JSCE 18, 141-151.

(33) Watanabe, Y. and Saeki, H. (2002) Velocity field after wave breaking. Int. J. Numer. Methods Fluids 39, 607-637.

(34) Nadaoka, K., Hino, M. and Koyano, Y. (1989) Structure of the turbulent flow field under breaking waves in the surf zone. J. Fluid Mech. 204, 359-387.

(35) Watanabe, Y., Saeki, H. and Hosking, R.J. (2005) Three-dimensional vortex structures under breaking waves. J. Fluid Mech. 545, 291-328.

(36) Watanabe, Y., Saruwatari, A. and Ingram, D.M. (2008) Free-surface flows under impacting droplets. J. Comput. Phys. 227, 2344-2365.

(37) Cundall, P.A. and Strack, O.D.L. (1979) A discrete numerical model for granular assembles. Geotechnique 29, 47-65.

(38) Gotoh, H. and Sakai, T. (1997) Numerical simulation of sheetflow as granular material. J. Waterw. Port Coast. Ocean Eng. 123, 329-336.

(39) Reddy, K.A. and Kumaran, V. (2010) Dense granular flow down an inclined plane--A comparison between the hard particle model and soft particle simulations. Phys. Fluids 22, 113302.

(40) Estep, J. and Dufek, J. (2013) Discrete element simulations of bed force anomalies due to force chains in dense granular flows. J. Volcanol. Geotherm. Res. 254, 108-117.

(41) Tnag, T., He, Y., Ren, A. and Zhao, Y. (2016) Investigation on wet particle flow behavior in a riser using LES-DEM coupling approach. Powder Technol. 304, 164-176.

(42) Mondal, S., Wu, C.-H. and Sharma, M. (2016) Coupled CFD-DEM simulation of hydrodynamic bridging at constrictions. Int. J. Multiph. Flow 84, 245-263.

(43) Harada, E., Gotoh, H. and Tsuruta, N. (2015) Vertical sorting process under oscillatory sheet flow condition by resolved discrete particle model. J. Hydraul. Res. 53, 332-350.

(44) Harada, E., Tsuruta, N. and Gotoh, H. (2016) Two-phase flow LES of the sedimentation process of a particle cloud. J. Hydraul. Res. 51, 186-194.

(45) Lucy, L.B. (1977) A numerical approach to the testing of the fission hypothesis. Astron. J. 82, 1013-1024.

(46) Gingold, R.A. and Monaghan, J.J. (1977) Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181, 375-389.

(47) Koshizuka, S. and Oka, Y. (1996) Moving particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123, 421-434.

(48) Gotoh, H. and Sakai, T. (1999) Lagrangian simulation of breaking waves using particle method. Coast. Eng. J. 41, 303-326.

(49) Khayyer, A. and Gotoh, H. (2008) Development of CMPS method for accurate water-surface tracking in breaking waves. Coast. Eng. J. 50, 179-207.

(50) Farahani, R.J. and Dalrymple, R.A. (2014) Three dimensional reversed horse-shoe vortex structures under broken solitary waves. Coast. Eng. 91, 261-279.

(51) Gotoh, H., Ikari, H., Memita, T. and Sakai, T. (2005) Lagrangian particle method for simulation of wave overtopping on a vertical seawall. Coast. Eng. J. 47, 157-181.

(52) Shao, S., Ji, C., Graham, D.I., Reeve, D.E., James, P.W. and Chadwick, A.J. (2006) Simulation of wave overtopping by an incompressible SPH model. Coast. Eng. 53, 723-735.

(53) Shadloo, M.S., Weiss, R., Yildiz, M. and Dalrymple, R.A. (2015) Numerical simulation of long wave runup for breaking and nonbreaking waves. Int. J. Offshore Polar Eng. 25, 1-7.

(54) Khayyer, A. and Gotoh, H. (2009) Modified moving particle semi-implicit methods for the prediction of 2D wave impact pressure. Coast. Eng. 56, 419-440.

(55) Lee, B.H., Park, J.C., Kim, M.H. and Hwang, S.C. (2011) Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads. Comput. Methods Appl. Mech. Eng. 200, 1113-1125.

(56) Altomare, C., Crespo, A.J.C., Dominguez, J.M., Gomez-Gesteira, M., Suzuki, T. and Verwaest, T. (2015) Applicability of smoothed particle hydrodynamics for estimation of sea wave impact on coastal structures. Coast. Eng. 96, 1-12.

(57) Farahani, R.J., Dalrymple, R.A., Herault, A. and Bilotta, G. (2014) Three-dimensional SPH modeling of a bar/rip channel system. J. Waterw. Port Coast. Ocean Eng. 140, 82-99.

(58) Delorme, L., Colagrossi, A., Souto-Iglesias, A., Zamora-Rodriguez, R. and Botia-Vera, E. (2009) A set of canonical problems in sloshing, part I: pressure field in forced roll-comparison between experimental results and SPH. Ocean Eng. 36, 168-178.

(59) Gotoh, H., Khayyer, A., Ikari, H., Arikawa, T. and Shimosako, K. (2014) On enhancement of incompressible SPH method for simulation of violent sloshing flows. Appl. Ocean Res. 46, 104-105.

(60) Shibata, K. and Koshizuka, S. (2007) Numerical analysis of shipping water impact on a deck using a particle method. Ocean Eng. 34, 585-593.

(61) Le Touze, D., Marsh, A., Oger, G., Guilcher, P.M., Khaddaj-Mallat, C., Alessandrini, B. and Ferrant, P. (2010) SPH simulation of green water and ship flooding scenarios. J. Hydrodyn. Ser. B 22, 231-236.

(62) Gotoh, H. and Sakai, T. (2006) Key issues in the particle method for computation of wave breaking. Coast. Eng. 53, 171-179.

(63) Panizzo, A. and Dalrymple, R.A. (2004) SPH modelling of underwater landslide generated waves. Proc. 29th Int. Conf. Coastal Eng., Lisbon, ASCE, 1147-1159.

(64) Fu, L. and Jin, Y.C. (2015) Investigation of non deformable and deformable landslides using mesh-free method. Ocean Eng. 109, 192-206.

(65) Rafiee, A. and Thiagarajan, K.P. (2009) An SPH projection method for simulating fluid-hypoelastic structure interaction. Comput. Methods Appl. Mech. Eng. 198, 2785-2795.

(66) Shibata, K., Koshizuka, S., Sakai, M. and Tanizawa, K. (2012) Lagrangian simulations of ship-wave interactions in rough seas. Ocean Eng. 42, 13-25.

(67) Hwang, S.C., Khayyer, A., Gotoh, H. and Park, J.C. (2014) Development of a fully Lagrangian MPS-based coupled method for simulation of fluid-structure interaction problems. J. Fluids Structures 50, 497-511.

(68) Colagrossi, A., Marrone, S., Bouscasse, B. and Broglia, R. (2015) Numerical simulations of the flow past surface-piercing objects. Int. J. Offshore Polar Eng. 25, 13-18.

(69) Wei, Z., Dalrymple, R.A., Herault, A., Bilotta, G., Rustico, E. and Yeh, H. (2015) SPH modeling of dynamic impact of tsunami bore on bridge piers. Coast. Eng. 104, 26-42.

(70) Colagrossi, A. and Landrini, M. (2003) Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J. Comput. Phys. 191, 448-475.

(71) Shao, S.D. and Lo, E.Y.M. (2003) Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv. Water Resour. 26, 787-800.

(72) Gotoh, H. and Khayyer, A. (2016) Current achieve ments and future perspectives for projection-based particle methods with applications in ocean engineering. J. Ocean Eng. Mar. Energy. 2, 251-278.

(73) Gotoh, H. (2009) Lagrangian particle method as advanced technology for numerical wave flume. Int. J. Offshore Polar Eng. 19, 161-167.

(74) Inutsuka, S. (1994) Godunov-type SPH. J. Ital. Astron. Soc. 65, 1027-1031.

(75) Inutsuka, S. (2002) Reformulation of smoothed particle hydrodynamics with Riemann solver. J. Comput. Phys. 179, 238-267.

(76) Antuono, M., Colagrossi, A. and Marrone, S. (2012) Numerical diffusive terms in weakly-compressible SPH schemes. Comput. Phys. Commun. 183, 2570-2580.

(77) Le Touze, D., Colagrossi, A., Colicchio, G. and Greco, M. (2013) A critical investigation of smoothed particle hydrodynamics applied to problems with free surfaces. Int. J. Numer. Methods Fluids 73, 660-691.

(78) Randles, P.W. and Libersky, L.D. (1996) Smoothed particle hydrodynamics: some recent improvements and applications. Comput. Methods Appl. Mech. Eng. 139, 375-408.

(79) Chen, J.K., Beraun, J.E. and Jih, C.J. (1999) An improvement for tensile instability in smoothed particle hydrodynamics. Comput. Mech. 23, 279-287.

(80) Oger, G., Doring, M., Alessandrini, B. and Ferrant, P. (2007) An improved SPH method: towards higher order convergence. J. Comput. Phys. 225, 1472-1492.

(81) Fatehi, R. and Manzari, M.T. (2011) Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives. Comput. Math. Appl. 61, 482-498.

(82) Bonet, J. and Lok, T.S. (1999) Variational and momentum preservation aspects of smooth particle hydrodynamic formulation. Comput. Methods Appl. Mech. Eng. 180, 97-115.

(83) Sun, J.W., Liang, S.X., Sun, Z.C. and Zhao, X.Z. (2010) Simulation of wave impact on a horizontal deck based on SPH method. J. Mar. Sci. Appl. 9, 372-378.

(84) Xie, J., Nistor, L. and Murty, T. (2012) A corrected 3-D SPH method for breaking tsunami wave modelling. Nat. Hazards 60, 81-100.

(85) Khayyer, A., Gotoh, H. and Shao, S.D. (2008) Corrected incompressible SPH method for accurate water-surface tracking in breaking waves. Coast. Eng. 55, 236-250.

(86) Khayyer, A. and Gotoh, H. (2010) A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method. Appl. Ocean Res. 32, 124-131.

(87) Khayyer, A. and Gotoh, H. (2012) A 3D higher order Laplacian model for enhancement and stabilization of pressure calculation in 3D MPS-based simulations. Appl. Ocean Res. 37, 120-126.

(88) Khayyer, A. and Gotoh, H. (2011) Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230, 3093-3118.

(89) Swegle, J.W., Attaway, S.W., Heinstein, M.W., Mello, F.J. and Hicks, D.L. (1994) An analysis of smooth particle hydrodynamics. Sandia Report SAND93-2513.

(90) Dilts, G.A. (1999) Moving least squares hydrodynamics: consistency and stability. Int. J. Numer. Methods Eng. 44, 1115-1155.

(91) Monaghan, J.J. (1992) Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys. 30, 543-574.

(92) Ikari, H., Khayyer, A. and Gotoh, H. (2015) Corrected higher order Laplacian for enhancement of pressure calculation by projection-based particle methods with applications in ocean engineering. J. Ocean Eng. Mar. Energy 1, 361-376.

(93) Kondo, M. and Koshizuka, S. (2011) Improvement of stability in moving particle semi-implicit method. Int. J. Numer. Methods Fluids 65, 638-654.

(94) Khayyer, A. and Gotoh, H. (2013) Enhancement of performance and stability of MPS mesh-free particle method for multi-phase flows characterized by high density ratios. J. Comput. Phys. 242, 211-233.

(95) Tsuruta, N., Khayyer, A. and Gotoh, H. (2013) A short note on dynamic stabilization of moving particle semi-implicit method. Comput. Fluids 82, 158-164.

(96) Balsara, D.S. (1995) Von Neumann stability analysis of smoothed particle hydrodynamics--Suggestions for optimal algorithms. J. Comput. Phys. 121, 357-372.

(97) Basa, M., Quinlan, J.N. and Lastiwka, M. (2009) Robustness and accuracy of SPH formulations for viscous flow. Int. J. Numer. Methods Fluids 60, 1127-1148.

(98) Tsuruta, N., Khayyer, A. and Gotoh, H. (2015) Space potential particles to enhance the stability of projection-based particle methods. Int. J. Comput. Fluid Dyn. 29, 100-119.

(99) Adami, S., Hu, X.Y. and Adams, N.A. (2013) A transport-velocity formulation for smoothed particle hydrodynamics. J. Comput. Phys. 241, 292-307.

(100) Li, S. and Liu, W.K. (2004) Meshfree Particle Methods. Springer, Heidelberg.

(101) Liu, P.L.-F., Yeh, H. and Synolakis, C. (2008) Advanced Numerical Models for Simulating Tsunami Waves and Runup. World Scientific, Singapore; Hackensack, N.J.

(102) Liu, M.B. and Liu, G.R. (2016) Particle Methods for Multi-Scale and Multi-Physics. World Scientific, Singapore; Hackensack, N.J.

(Received Mar. 6, 2017; accepted May 30, 2017)


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.

By Hitoshi GOTOH * [1], [dagger]] and Akio OKAYASU * [2]

(Communicated by Kiyoshi HORIKAWA, M.J.A.)

* [1] Department of Civil and Earth Resources Engineering, Kyoto University, Kyoto, Japan.

* [2] Department of Marine Resources and Energy, Tokyo University of Marine Science and Technology, Tokyo, Japan.

([dagger]) Correspondence should be addressed: H. Gotoh, Department of Civil and Earth Resources Engineering, Kyoto University, Kyoto daigaku-Katsura, Nishikyo-ku, Kyoto 615-8540, Japan (e-mail:

Abbreviations: ALE: Arbitrary Lagrangian-Eulerian; CADMAS-SURF: SUper Roller Flume for Computer Aided Design of MAritime Structure; C-CUP: CIP-Combined and Unified Procedure; CERC: Coastal Engineering Research Center; CFD: Computational Fluid Dynamics; CHL: Corrected HL; CIP: Cubic Interpolation Pseudo-particle or Constrained Interpolation Profile; CISPH: Corrected ISPH; CMPS: Corrected MPS; DEM: Distinct Element Method; DNS: Direct Numerical Simulation; DS: Dynamic Stabilization; ECS: Error Compensating Source; GC: Gradient Correction; HL: Higher-order Laplacian; HS: Higher-order Source term; ISPH: Incompressible SPH; LES: Large Eddy Simulation; MAC: Marker and Cell; MPS: Moving Particle Semi-implicit; NWF: Numerical Wave Flumes; PPE: Poisson Pressure Equation; RANS: Reynolds-Averaged Navier-Stokes; SGS: Sub-Grid-Scale; SPH: Smoothed Particle Hydrodynamics; SPP: Space Potential Particle; VOF: Volume of Fluid; WCSPH: Weakly Compressible SPH.

Caption: Fig. 1. 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 x [10.sup.-3]m.

Caption: Fig. 2. Hydrostatic pressure of water in the rectangular tank calculated by the CMPS-HS-HL-ECS method in comparison to the original MPS method.

Caption: Fig. 3. Numerical wave flume by the original MPS and the CMPS-HS-HL-ECS methods for simulating wave run-up on uniform slope.

Caption: Fig. 4. 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.

Caption: Fig. 5. Evolution process of a square patch of 2-dimensional inviscid fluid subjected to a rigid rotation, in which the initial square side is 1.0m in length and angular velocity around the center of the square is 10.0[s.sup.-1].

Caption: Fig. 6. Concept of Dynamic Stabilization (DS) scheme.

Caption: Fig. 7. Settlement process of heavy particles, the density of which is 2.65 x [10.sup.3] kg/[m.sup.3], 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 x [10.sup.-3]m.

Caption: Fig. 8. Performance of SPP in Karman vortex simulation under Re = 1200, with time resolution [DELTA]t = 2.5 x [10.sup.-3]s.

Caption: Fig. 9. Relationship of the accurate particle methods.

Caption: Fig. 10. Dam breaking simulation by the accurate particle methods.

Caption: Fig. 11. Required computational time in dam breaking simulations.
COPYRIGHT 2017 The Japan Academy
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2017 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Gotoh, Hitoshi; Okayasu, Akio
Publication:Japan Academy Proceedings Series B: Physical and Biological Sciences
Article Type:Report
Geographic Code:9JAPA
Date:Aug 1, 2017
Previous Article:Proceedings at the 1110th general meeting.
Next Article:Biological control for grapevine crown gall using nonpathogenic Rhizobium vitis strain ARK-1.

Terms of use | Privacy policy | Copyright © 2020 Farlex, Inc. | Feedback | For webmasters