3D LBM NUMERICAL SIMULATION OF MIXED SAND SORTING UNDER OSCILLATORY FLOWS AND PROGRESSIVE WAVES

A 3D numerical model based on Lattice Boltzmann Method (LBM) has been developed to investigate the particle behavior of mixed-grain-size sediments under oscillatory flows and progressive waves. The model considers both particle-particle interaction and particle-fluid interaction with free surface flow and SGS turbulence model. By applying the present model, simulations of mixed-grain-size sediments with different influence factors (i.e., Shields parameters, periods, bottom layer thicknesses) under symmetric oscillatory flows as well as under progressive cnoidal wave are performed to understand the vertical sorting process of graded sands. The results of mixed sands under symmetric oscillatory flows show that the water particle semi-excursion is a key influence factor on the armoring effect. For the cnoidal wave, the vertical sorting proceeds until the wave crest passes by and it fully develops in a complete period. The concentration centroid of large particles becomes higher landward within 2 periods.


INTRODUCTION
The sorting process of mixed-grain-size sediments is essential for study of sand transport in sheetflow regime and is a fundamental element for predicting the beach profiles evolution, as naturally occurring sand is almost always graded in size.Numerous hydrodynamic experiments (Dibajnia and Watanabe 2000, Ahmed and Sato 2003, Hassan and Ribberink 2005) and numerical simulations (Gotoh and Sakai 1997, Dong and Zhang 1999, Thaxton and Calantoni 2006, Harada and Gotoh 2008) have been conducted to study the mechanisms of vertical sorting processes and the transport rates of mixed-grain-size sediments.Most of the above studies are carried out in the sheetflow regime, as a large amount of sediment is transported under sheetflow.
In the sheetflow regime, the concentration is so high that the existence of the sediment particles influences the flow characteristics.Hence, sub-mechanisms of both the particle-fluid interaction and particle-particle interaction should not be ignored in the numerical models.However, the interaction between fluid and particle still remains as a complex problem, and the fundamental physical phenomena are often described in an empirical fashion.The motion of the particles is driven collectively by gravity and the hydrodynamic forces exerted by the fluid, and may also be altered by the interaction between particles.On the other hand, the fluid flow pattern can be greatly affected by the presence of the particles, and is often of a turbulent nature (Feng et al. 2010).These complicated interactions and inter-granular stress make it difficult to establish a numerical model for sediment transport in oscillatory flows and progressive waves including both particle-fluid interaction and particle-particle interaction.
A great many efforts have been done to solve the above problem.For example, Dong and Zhang (1999) proposed a two-phase flow model for the simulation of sediment transport in oscillatory flow.Their model is able to calculate the time-varying concentrations and particle velocities and is superior to the single-phase flow model, which is only available for the computation of time-averaged concentrations and flow velocities.However, this kind of Euler-Euler coupling model (i.e., Hsu, Jenkins and Liu 2004) is good at dealing with fluid-particle interaction but weak at particle-particle interaction.The motion of particles is Langrangian in nature, while both flow model and particle model are Euler.Another numerical model, namely a granular model coupled with one-way or two-way flow models, is also proposed.This Langrangian-Euler coupling model (Gotoh andSakai 1997, Thaxton andCalantoni 2006) improves the deficiency of the Eulerian model by emphasizing the behaviors of particles, but simplifies the characteristics of fluid flow.Gotoh et al. (2002) and Harada and Gotoh (2008) proposed an improved 3D granular material model based on the Distinct Element Method (DEM), which considers particle-fluid and particle-particle interactions.
Recently, Lattice Blotzmann Method (LBM) has been used to solve the problems of particle-fluid interaction and particle-particle interaction.LBM has been recognized as an alternative solution scheme to the conventional Computational Fluid Dynamics methods employing Navier-Stokes equations with the advantages of stabile, ease of implementation and extensibility, and intrinsic parallel structure allowing high performance computing.Ladd (1994a, b) is the first to develop LBM for chemical particle-fluid suspension by considering both particle-fluid interaction and particle-particle interaction.Feng and Michaelides (2003) applied Ladd's particle-fluid suspension model with a modified lift force to study the behavior of particle motions in 2D space.Feng et al. (2010) combined LBM with DEM to investigate 3D contact of particles.In their model, DEM is applied for the motion of particles and discrete Boltzmann equation is used for the solution of fluid motion.He et al. (2013) expanded Ladd's LBM for particle-fluid suspension with free surface flow and SGS turbulence model to study the sediment transport of uniform particles under open channel flows.In this paper, the 3D numerical movable bed model based on LBM (He et al. 2013) is applied to investigate the processes of mixed-grain-size sediments under both symmetric sinusoidal oscillatory flows and asymmetric fluid particle motion based on cnoidal wave theory.The concentration centroids of mixed-grain-size sediments under symmetric oscillatory flow are compared to the experiment results of Harada and Gotoh (2008).

MODEL DESCRIPTION
Lattice Boltzmann Method (LBM) is originated from the Lattice Gas Automata (LGA), which can be considered as a simplified fictitious molecular dynamics model, where space, time and particle velocities are all discrete.With an appropriate choice of the lattice symmetry, lattice gas models (Frisch et al. 1987) in fact represent numerical solutions of the Navier-Stokes equations and are thus able to describe macroscopic hydrodynamic problems (Chen et al. 1992, Qian et al. 1992).
In LBM, the macroscopic fluid can be represented by an aggregation of microscopic virtual particles moving on lattices, and the flow of fluid can be obtained through the calculations of time evolution of streaming and collision processes of these virtual particles.According to Ladd's (1994a, b) LBM model for chemical particle-fluid suspension, the interaction force between fluid and movable solid particles can be calculated from the momentum exchange by applying the link-bounce-back scheme for the movable solid particles in the fluid.The basic LBM algorithm consists of two steps; the streaming step and the collision step.Before streaming, the fraction of virtual particles moving with a certain velocity of an independent lattice is represented by the DF t f i , x .During the streaming step, all DFs are advected to the adjacent fluid lattice with the corresponding velocity direction (see Fig. 1).During the collision step, the virtual particles are redistributed with a percentage of i Δ to a local equilibrium state.All the DFs after collision step can be obtained by the Boltzmann Equation (BE), (1) The virtual particle's DF, t f i , x shows the number of particles with a velocity vector of i direction in the lattice of x position at t time.With the particle distribution function, the density, velocity, momentum flux and local Equilibrium Distribution Function (EDF) of the fluid field can be expressed respectively with weighted coefficient of i e a as follows: 1 e e uu e j (3) The particle distribution function for the next step, namely the post-collision DF, where 1 is the unit tensor, and s c is the speed of sound, expressed with a function of virtual particle velocity e, . The pressure is determined by the ideal gas equation of state, The post-collision non-equilibrium momentum flux in the above equation is calculated as where the eigenvalues λ and v λ are related to shear and bulk viscosities, respectively, and the overbar indicates the traceless projection.The non-equilibrium momentum flux neq Π can be obtained directly by subtracting the equilibrium momentum flux eq Π from the total momentum flux Π .
The link-bounce-back collision rule is applied for the moving boundary condition of a solid particle.It is assumed that the surface of an individual particle is constructed by a series of boundary nodes which locate in the center of fluid nodes and solid nodes inside the rigid particle (Fig. 2).It is also assumed that there are no DFs in the solid nodes, and the momentum transfers between fluid and particle are through the link-bounce-back scheme.The velocities along links cutting the boundary surface are represented by arrows.Set the location of fluid node just outside the particle surface as x and the position of solid node adjacent to the particle surface as (the subscript b represents the solid boundary node), then each of the corresponding population densities is updated according to a simple rule which takes into account the motion of the particle surface.Figure 3 shows the momentum transfers between DFs of fluid nodes and boundary nodes.It should be emphasized that the velocity components transfer parallel to the 18 DF directions.Hence, the DFs of moving boundary node can be updated from x is the location of surface boundary node, . Originally, the bounce-back method (Ladd, 1994a) is assumed with fluid filled inside of the particle.However, it is impossible to treat a solid particle as a rigid body, as time lag of momentum exchange exists due to the inertia of fluid inside the particle.The method proposed by Aidun et al. (1998), in which no fluid is included inside particles and momentum exchange only occurs at the boundary nodes of particle surfaces, is applied in the present model.The force F and torque T exerted on particles are composed of external forces ext F (including drag forces, gravity force, buoyance force, etc.) depending on the incoming DFs, and forces related to the particle velocity U and angular velocity Ω.
where ζ FU , ζ FΩ , ζ TU and ζ TΩ are high-frequency friction coefficients (Nguyen and Ladd 2002).The velocity independent forces F 0 and torques T 0 are calculated at the half-time step, as Then the particle velocity U and angular velocity Ω at the next time step are updated implicitly: where m is the mass, and I is the inertia momentum.
With the updated velocity and angular velocity, the new position of an individual particle is finally updated with the following simple equations:

SGS Turbulence Term
Sub Grid Smagorinsky (SGS) turbulence model is rather easily applied by adding the modified total kinematic viscosity to the collision term of the general Lattice Boltzmann Equation with BGK approximation (He et al. 2013).Specifically, the eddy viscosity coefficient t ν of SGS turbulence term is added to the normal kinematic viscosity coefficient ν , resulting in a total kinematic viscosity . The relaxation time ( λ τ 1 ) related to kinematic viscosity is then revised with The equation of Smagorinsky eddy viscosity is given in the form where is the magnitude of the local stress tensor, C is Smagorinsky constant, and x Δ is the lattice length.Since the filter length is assumed here less than 2 x Δ , a Smagorinsky constant of 0.1 C is applied throughout numerical simulations.The local stress tensor in LBM model can be determined from

Lubrication Forces
When particles are getting close to each other or impinging on a wall, a strong repulsive force occurs, which is caused by the fluid squeezed between the solid bodies.These lubrication flows generate very high pressures in the gap and are generally difficult to be solved.However, Nguyen and Ladd (2002) successfully proposed a method to incorporate lubrication corrections into a lattice Boltzmann simulation as where 1 a and 2 a are the radii of the particles, c Δ is the cut off for the added lubrication force including the normal cN Δ , tangential cT Δ and rotational cR Δ components,

Free Surface Flows
Free surface condition is included in the model for the simulation of open channel flows, oscillatory flows, and progressive waves.The intrinsic characteristic of LBM, of which all the macroscopic variables (density, momentum, momentum flux of the fluid) can be simulated with a serious of microscopic DFs of virtual particles in each lattice, makes the expansion of free surface condition more conveniently without any huge computational cost.In the present paper, free-surface flows is simulated based on the method proposed by Thürey and Rüde (2009), in which a fluid fraction parameter ε is used to distinguish fluid ( 1 ε ), gas ( 0 ε ) and interface lattices ( 1 0 ε ).With this marking method, tracking of the free surface is feasible in the LBM simulation.This can be done by the following steps: (1) calculation of the mass exchange, (2) computation of the gas lattice next to an interface lattice, (3) the ordinary streaming and collision processes, (4) re-initialization of lattice types and redistribution of extra mass to the surrounding lattices.
From the fluid fraction computed with the lattice mass m and density ρ, ρ m ε , the movement of the fluid interface is traced by the calculation of the mass that is contained in each lattice.It is assumed that the mass exchange only happens between a fluid lattice and an interface lattice or between interface lattices.For an interface lattice and a fluid lattice at ( ), the exchanged mass is given as When the mass exchange occurs for two interface lattices, the area of the fluid interface between these two lattices is considered.This is approximated by averaging the fluid fraction values of the two lattices, and the mass exchanged between them becomes With the results from Eq. 22 and Eq. 23, the total mass for the next time step is computed from When a gas lattice is located next to an interface lattice, it is assumed that the atmosphere has a pressure of ρ ρ G and the fluid has a much lower kinematic viscosity.Hence, the gas at the interface simply flows along the direction which is pushed by the fluid.Then the DFs from the gas lattice are modified with With all DFs for the interface lattice, the density is calculated during the collision step and then the updated density is used to check whether the interface lattice filled or emptied by the following equations: where κ is an additional offset and set to 10 -3 in the present simulation.The extra mass caused by lattice type change is then redistributed to its neighbor lattices to ensure the mass conservation.
The extra mass can be calculated from  In order to study the influence factor to the vertical sorting process of mixed-grain-size sediments, cases with different Shields parameters, oscillatory flow periods and sediment layer thicknesses are conducted under symmetric sinusoidal oscillatory flows, as shown in Table 1.Case 1 is set as the standard case for comparisons.Details of the components of sediment mixture in Table 1 are listed in Table 2.The Shields parameter Ψ is calculated with the maximum flow velocity of the oscillatory flow,

Initialization of Cnoidal Wave
The second approximation of cnoidal wave theory is applied to generate a wave with a still water depth h of 20cm, wave height H of 14cm and period T of 1.4s.The initial wave level is the same as the still water depth and initial fluid velocity is zero.Periodic boundaries with time series of the water level and water velocities in x-and z-direction are given at left and right sides of a 211cm×16cm×50cm flume.The EDFs of each lattice at the left and right side of the computation domain is then updated with the given water velocities by using Eq. 3, and the fluid fraction parameter ε at the left and right boundary is controlled by the given water level.
Sediment mixture C in Table 2 with bottom layer thickness δ of 4cm and specific gravity s of 2.65 are placed in the numerical flume to obtain a view of the behavior of mixed-grain-size sediments.Using Eq. 32, the corresponding Shields parameter Ψ to the maximum flow-velocity amplitudes of 160cm/s is 0.19.

NUMERICAL RESULTS OF MIXED SANDS UNDER OSCILLATORY FLOWS
As listed in Table 1, the first four cases with different simulation characteristics are performed to study the influence of Shields parameters, periods of oscillatory flows and sediment layer thicknesses on the vertical sorting progress of mixed-grain-size sediments.Particles with large diameters, which is 10 times larger than usually used in experiments (i.e., Ahmed andSato 2003, Hassan andRibberink 2005) but the same as Harada and Gotoh's (2008), are used in the simulations.Oscillatory flows with relatively large maximum velocities and short periods of 1s and 2s, which also have been applied by Harada et al. (2010) and Gotoh et al. (2002)

Concentration Centroids of Mixed-Grain-Size Sediments under Oscillatory Flows
For the simulation of oscillatory flows, a function for the flow velocity of t ω U U b sin max is applied.The phase change of the movement of mixed-grain-size sediment particles follows that of the oscillatory flow.However, a relatively small phase lag of 16.5° is detected due to the turbulence effect.
Figure 5 shows the time series of both vertical (Fig. 5a, b, c) and transverse (Fig. 5d, e, f) concentration centroids of graded sediment particles for all four cases under symmetric oscillatory flows.The concentration centroid of particles is an important parameter to judge the vertical sorting process of mixed-grain-size sediments.It is defined as the deviation of the summary of centers of an individual group of one uniform-grain-size sediments to that of the initial condition.For instance, the vertical concentration centroids of particles are determined by 0 The subscribe i = 1, 2, 3 represents small, middle and large particles with diameter of d i .In Eq. 33, is the initial vertical concentration centroid of particles, and t i z is the vertical concentration centroid of particles calculated by  Sediment particles with all three diameters start to move in the acceleration phase, small particles move downward through the voids between the large particles in the deceleration phase, leave large particles on the top of the bed-load layer.Then the armoring coat is generated.It is observed from Fig. 5a, b, c, the vertical concentration centroids of particles in all cases except Case 3 (with double flow period) become nearly constant at after 2 periods, which indicates the ending of the developing stage of grading.There are predominant layers of large, middle and small particles from the surface layer to the bottom.For clarity, the process in the first period is denoted as the developing stage of grading and the process in the third period is marked as the fully developed stage of grading in the following sections.
The vertical concentration centroids of larger particles increase significantly and middle particles rise up with a small height, while the vertical concentration centroids of small particles decrease in Case 2 (Fig. 5a).The numerical result is quite similar to the experimental result in Harada and Gotoh's (2008) case with similar Shields parameter of 0.27.The vertical concentration centroids of large and middle particles of Case 2 with smaller Shields parameter are lower than that of Case 1 while the vertical concentration centroid of small particles is higher than that of Case 1.This indicates the phenomenon of vertical sorting is more distinct with larger Shields parameter.In Case 3, where the period of oscillatory flow increases to the double of Case 1, the concentration centroids of large particles fall down while the concentration centroid of middle particles rise up after 1.5T (Fig. 5b).Armor effect is not sufficiently working and the middle and large particles become suspended and mixed in the surface layer.The flow velocity max b

U
and Shields parameter for Case 1 and Case 3 are the same, while the angular velocity ω becomes half.Hence the water particle semi-excursion A (in the function ) increases, which induces the Keulegan-Carpenter number to increase, resulting the increase of drag force with the remained pressure force.As a result, sediment particles become more active, and the armoring effect breaks.For Case 3, the longer computational domain of 24cm in length with totally 5,265 particles is used to ensure that the maximum step length of particles during the half-period of oscillation is sufficiently shorter than the length of numerical flume.Particles in the bottom layer move more slowly than particles in the top layer in Case 4 with double sediment layer thickness, hence the centroids of concentrations of large and middle particles are lower than those of Case 1 while the centroid of concentrations of small particles is higher than that of Case 1 (Fig. 5c).
The transverse concentration centroids (Fig. 5d, e, f) of all particles show visible high-frequency fluctuations in comparison with the vertical concentration centroids.This is due to the repeated collisions between particles which are able to move freely in the transverse direction under periodic boundary conditions.However, the horizontal concentration centroids of particles in all cases show no sizable differences.The maximum deviation is less than the double diameter of the small particle.

Vertical Sorting Process of Mixed-Grain-Size Sediments under Oscillatory Flows
Figure 6 shows some computational snapshots with both side view and top view of Case 1 at different phases of the movement of mixed-grain-size sediments.In the developing stage of grading, a mixture of particles with all three sizes can be obviously detected at the maximum velocity phase (Fig. 6a) from the top view.At the decelerating phase in the developing stage of grading, more red large particles and white middle particles rise up to the surface layer (Fig. 6b).When the grading reaches to the fully developed stage, the surface layer is occupied with the large and middle particles at the maximum velocity phase (Fig. 6c) at this stage, a clear vertical sorting with large, middle and small particles locating from top to bottom can be seen from the side view (Fig. 6c).
Figure 7 shows both the side view and top views of computational snapshots for the other three cases at different phases in both developing stage and fully developed stage of grading.At the maximum velocity phase in the fully developed stage of grading, very few blue small particles can be seen from gaps between large and middle particles from the top view of Case 1(Fig.6c).From the side view of Case 2 with smaller Shields parameter, red large particles and blue small particles can be found between while middle particles in the middle layer (Fig. 7a), while large particles and small particles can hardly be seen in the middle layer in Case 1 (Fig. 6c).Hence, the phenomenon of vertical sorting is more pronounced with larger Shields parameter.
In the double period case, Case 3, vertical sorting is not fully attained.As shown in Fig. 7b, at the maximum velocity phase in the third period, armoring effect is not working.All particles, especially the middle and large particles are suspended, and as a result, more white middle particles rather than the red large particles are detected from the top view.This coincides with the result of vertical concentration centroids discussed in the above section.
From the top view of Case 4 (Fig. 7c), where the thickness δ of sediments increases to the double of Case 1, more large and middle particles can be found compared to Case 1 (Fig. 6b) at the decelerating phase in the developing stage of grading.Meanwhile, a vertical sorting with large, middle and a mixture of all particles from top to bottom can be seen clearly from the side view.This shows that armoring takes effect earlier when the sediment layer thickness increases.

Mean Velocity of Mixed-Grain-Size Sediments under Oscillatory Flows
Figure 8 shows the dimensionless mean velocity of sediment particles for Case 1, namely the average particle velocity with the same diameter (u pi ) over the maximum particle velocity in one period (u pmax ).At the accelerating phase (Fig. 8a), at the maximum velocity phase (Fig. 8b) and at the decelerating phase (Fig. 8c) in the developing stage of grading, the mean velocity of all particle shows almost the same distribution of downward convex velocity profile.At the maximum velocity phase (Fig. 8e) and the decelerating phase (Fig. 8f) in the fully developed stage of grading, the velocity profile of small particle does not reach the upper layer due to armoring effects.The distribution of all particles in the fully developed stage of grading shows a nearly linear profile at the accelerating phase (Fig. 8d) and the maximum velocity phase (Fig. 8e), and a tendency of upward convex profile at the decelerating phase (Fig. 8f).

NUMERICAL RESULTS OF MIXED SANDS UNDER PROGRESSIVE CNOIDAL WAVE
The second-order approximation of cnoidal wave theory is applied to generate a cnoidal wave with a still water depth of 20cm, wave height of 14cm and period of 1.4s.Figure 9 shows the free surface of water level elevation and flow velocity fields of the simulated progressive cnoidal wave in the x-z plane at different phases.A forward leaning cnoidal wave is generated (Fig. 9a) and the wave breaks at some phases (Fig. 9b) due to the water level difference induced during the simulation (He et al. 2014).

Vertical Sorting Process of Mixed-Grain-Size Sediments under Progressive Wave
The vertical sorting process of mixed sands in three areas of the whole computation domain, marked as seaward, center and landward areas (Fig. 9), with sub-computation domains of 16×16×16cm, are focused for comparisons.Figure 10 shows some snapshots of sorting process of mixed-grain-size sediment under the progressive cnoidal wave at the center area of the computational domain.At half wave period before the wave crest passing by (Fig. 9a) the flow velocity near the bed-load layer is so small that particles are unable to start moving (Fig. 10a).When the wave crest passes by, the large bottom flow velocity pushes the sediment particles with all three diameters to rising up.However, after the wave crest passed by, the flow velocity slows down causing the sediment particles to settle down.Small particles move down between the gaps of large particles during this process.Consequently, the vertical sorting is observed in a complete wave period (Fig. 10b).A clearer vertical sorting can be found near three periods in both side view and top view in Fig. 10c, although an armor coat is not formed as a result of vertical sorting.

Concentration Centroids of Mixed-Grain-Size Sediments under Progressive Wave
Figure 11 shows both the vertical and transverse concentration centroids of particles for landward and seaward, respectively.It can be seen that the vertical concentration centroid of particles for overall width is not uniform in wave progressing direction.The vertical concentration centroid of large particles is higher landward from around 0.7T, but becomes slightly higher seaward after 2.4T (Fig. 11a), due to the effects of spatially periodic boundary conditions.The transverse concentration centroids of particles show no significant difference for all particles located in the landward area and seaward area (Fig. 11b).

CONCLUDING REMARKS
A 3D numerical movable bed model based on LBM is established by including both particleparticle interaction and particle-fluid interaction with free surface flow and SGS turbulence model.The present model is applied to study the vertical sorting process of mixed-grain-size sediments under oscillatory flows and progressive waves.
In order to understand the influence factors on the vertical sorting process of graded sands, four cases with different Shields parameters, period, sediment layer thickness under symmetric oscillatory flow are performed.The simulation results show that the water particle semi-excursion A is one of the key influence factors in the sorting process, which increases with the increase of period.Sediment particles become more active with the increase of A, and as a result, the vertical sorting is not well established.The results also show that the phenomenon of vertical sorting can be more easily seen with larger Shields parameter.Armoring takes effect earlier when the thickness of sediments increases.However, the vertical concentration centroids of large and middle particles are lower and the vertical concentration centroid of small particles is higher with the increase of bottom layer thickness.
The results for mixed sand under asymmetric oscillatory flow generated from the second approximation of cnoidal wave theory show that the vertical sorting proceeds until the wave crest passes by and it fully develops in a complete period.The concentration centroids of particles for overall width are not uniform in the direction of the wave propagation.The concentration centroid of large particles becomes higher landward within 2 periods.

(
(Fig.1) is used in the present simulation, as it is necessary to determine the macroscopic variables, momentum flux in addition to density and velocity from the velocity Distribution Function (DF).The velocity vectors i e (i=0, 1, ..., 18) of each single lattice are given as 0 Δ are lattice interval in space and time, respectively).
* x is the post-collision DF at (x, t) in the direction b e , b e a is the coefficient of weights, 0 ρ is the mean density used to simplify the update procedure, and the subscript b ~ denotes the anti-parallel lattice vector, namely, b b e e ~.The local velocity of the particle surface b u is determined by the particle velocity U, angular velocity Ω and the center of the mass R, shown as

Figure 2 .
Figure 2. The surface boundary nodes of curved particles.

Figure 3 .
Figure 3. DFs before and after collision with stationary (a) and moving (b, c) boundary nodes.
is the shear viscosity.When gaps are smaller than c Δ , the lubrication forces work, otherwise 0 lub ij F .
Simulations of a mixture of grains with diameters of d 1 = 0.5cm, d 2 = 1.0cm and d 3 = 1.5cm are performed under sinusoidal oscillatory flows and progressive waves, respectively.Figure4shows the set-up of mixed-grain-size sediments under oscillatory flows with different sediment layer thicknesses.Small d 1 , middle d 2 and large d 3 sediment particles with specific gravity s = 2.65 are placed in a 16cm×16cm×16cm cubic flume for cases with the sediment layer thickness of 4cm and a 16cm×16cm×20cm flume for the case with the sediment layer thickness of 8cm.No-slip boundary conditions are used for the bottom of the computation domain and periodic boundary conditions for the others.The initial positions of particles are obtained from the settling down of all particles, which are initially suspended randomly without contact to the adjacent particles in the computation domain.

Figure 4 .
Figure 4. Initial set-up of mixed sands with (a) 4cm layer thickness and (b) 8cm layer thickness under oscillatory flows.
the present study.
, are used in the present simulations to save computation time.With a frictional coefficient 024 , 280cm/s give the Shields parameters of 0.33 and 0.58, respectively.
center of particles with the diameter of d i , and i n is the number of particles in the whole computation domain with the diameter of d i .

Figure 5 .
Figure 5. Comparisons of vertical concentration centroids (a, b, c) and transverse concentration centroids (d, e, f) of particles between Case 1 and other cases.

Figure 6 .Figure 7 .
Figure 6.Computational snapshots for Case 1 in developing stage (a), (b) and fully developed stage (c) of grading.

Figure 8 .
Figure 8. Mean velocity profiles of sediment particles in Case 1 in the developing stage of grading (a, b, c) and fully developed stage of grading (d, e, f).

Figure 9 .Figure 10 .
Figure 9. Water elevation and velocity field of wave generated from initial cnoidal wave.

Figure 11 .
Figure 11.Comparisons of (a) vertical and (b) transverse centroids of concentration of particles between the landward area and the seaward area.