LATTICE BOLTZMANN SIMULATION TO CHARACTERIZE ROUGHNESS EFFECTS OF OSCILLATORY BOUNDARY LAYER FLOW OVER A ROUGH BED

The 3-D lattice Boltzmann method was applied to characterize roughness effects of oscillatory boundary layer flow over a rough bed. The direct numerical simulation was carried out and the flow resistance of the flat and fixed bed was investigated. The position of the theoretical bed, equivalent roughness height and the behavior of friction factor at small values of relative roughness were obtained using the log-fit method.


INTRODUCTION
In coastal area, the seabed is composed of different roughness elements in different places, such as sand, stone, coral reef and so forth.The appearance of roughness enhances the interaction between the wave boundary layer and the seabed.It increases the impact on the sediment transport, flow resistance, energy dissipation, velocity distribution and so on.Plenty of theoretical, experimental and numerical investigations have been conducted to understand oscillatory boundary layer flow phenomena associated with different types of roughness elements and parameterize their effects (Jonsson 1966;Kajiura 1968;Sleath 1987;Jensen et al. 1989;Blondeaux et al. 2004;Dixen et al. 2008;Sana et al. 2009).However, the understanding of this problem is not sufficient because of its difficulties.
There are two commonly used parameters to characterize the resistance, one is the equivalent roughness height k s and the other is friction factor f w .For the flat and fixed beds, k s is expected to be on the order of the grain diameter d.The proportional coefficient is in the range of 1 ~ 5 according to the literature (Jonsson 1966;Kamphuis 1975;Nielsen 1992;You and Yin 2006;Camenen et al. 2009).For the fixed rippled beds, it is often suggested k s is in proportion to ripple height.Some investigation indicates it also depends on the ripple steepness with the coefficient varying from 7.4 to 27.7 (Jonsson and Carlsen 1976;Grant and Madsen 1982;Nielsen 1992;Van Rijn 1993).It can be seen that the equivalent roughness height remains quite difficult to determine and the existing results differ considerably from each other under the same conditions, even for the simple case.Take the flat and fixed beds for instance, although the coefficient 1 or 2.5 is often recommended, it still remains difficult to choose which one should be applied in practice.For f w , it is often considered as the function of Reynolds number Re a and the relative roughness a / k s (Jonsson 1966), where a is the amplitude of freestream motion.When the oscillatory boundary layer flow becomes rough turbulent which often takes place in field, it only varies with respect to the a / k s .Some previous investigators (Jonsson 1966;Kajiura 1968;Swart 1974;Kamphuis 1975;Fredsøe and Deigaard 1992;Soulsby et al. 1993;Simons et al. 2000) had provided a few implicit and explicit relationship between f w and a / k s in the rough turbulent regime.The existing expressions indicate that there are mainly two different point of view in the behavior of f w at small values of a / k s .On the one hand, f w approaches a constant value, for example, 0.24 suggested by Bagnold (1946), 0.25 by Kajiura (1968) and 0.30 by Jonsson (1966).On the other hand, it has such a relationship with a / k s which means that an increase in f w even at very small values of a / k s with decreasing a / k s .As the value of the f w in this range has a close relationship with the research on the stability of stones, rock and armour blocks in coastal engineering, it should be paid more attention to.So the purpose of the present paper is to shed more light on the understanding of two problems about the oscillatory boundary layer over a rough bed by numerical simulation.One is the equivalent roughness height for the fixed beds, here the flat beds condition will be taken into consideration.The other is the variation of the f w at small values of a / k s .They can be considered to be complementary to the earlier work.
To the best of authors' knowledge, the numerical simulation of oscillatory boundary layer over rough beds before were mainly based on Reynolds-averaged Navier-Stokes (RANS) equations, e.g., Puleo et al. (2004), Sana et al. (2009).Only a few works were carried out with large eddy simulation (LES) (Lohmann et al. 2006) or direct numerical simulation (DNS) (Blondeaux et al. 2004).Despite the studies quoted above, the understanding of the flow resistance is not sufficient because the grid resolution is not fine enough to represent the roughness elements and their effects are only included empirically.Recently, Fornarelli and Vittori (2009) (referred to as FV hereinafter) performed DNS of oscillatory boundary layer close to a rough bed which was composed of a layer of semi-spheres placed in a hexagonal pattern.It extended the research with considering the shape and distribution of the roughness elements.In this paper, the 3-D fully resolved simulation was carried out by a new promising method-the lattice Boltzmann (LB) method.It has the ability of dealing with an arbitrary complex curved boundary easily and straightforwardly using a fixed Eulerian mesh.What's more, it can be extended to simulate moving boundary problems without any extra effort.Meanwhile, the intrinsic parallelism feature also makes it easier to fully utilize the high performance computers.
The rest of the paper is organized as follows.In section 2, a brief introduction of LB method is given.The numerical validation of the model is presented in section 3. Section 4 and 5 are devoted to the simulation setup and discussion of the computational results.It ends with conclusions in section 6.

LATTICE BOLTZMANN METHOD
Different from conventional numerical methods, which solve the discrete macroscopic Navier-Stokes equations, LB method aims at modeling fluids in terms of the density distribution function of fictitious particles at the so-called mesoscopic level (Chen and Doolen 1998;Succi 2001).The fundamental concept is to construct simplified kinetic models incorporating mass and momentum conservation principles so that the macroscopic averaged properties obey the desired macroscopic equations.The time evolution of the density distribution function is described by the LB equation (Ladd and Verberg 2001).
where f i (r,t) is the density distribution function at location r at time t, ∆ i is the linearized collision operator (Ladd 1993), F i is the external force density term, e i is the discrete velocity of the simplified kinetic model.In the simulation, the three dimensional 19-velocity lattice model D3Q19 which was shown to be both stable and efficient (Mei et al. 2000) was used.The schematic diagram of it is shown in Fig. 1.The discrete velocity set for the D3Q19 model is defined as 11,11,11, 0,0,00, ,0,00,,00,0,16, ,,0,0,0,,718 The hydrodynamic properties of mass density ρ, momentum density ρu and momentum flux Π are moments of f i .
In order to take the particles into consideration, the link-bounce-back rules (Ladd and Verberg 2001), which make the boundary locate at the middle of the links between lattice nodes, have been used to match the velocity at the solid boundary and to account for the momentum transfer along the link direction. where a is the weight coefficient of the velocity direction, u b is the velocity of each boundary node and c s 2 = 1 / 3 in D3Q19 model.It can be seen from the work of Ladd (1994) as well as Feng and Michaelides (2002) that choosing the diameter of the particle to be long as 20 lattice units could represent the curved boundary of particle well and yield the reasonable results.

NUMERICAL VALIDATION
At first, the computation was performed in order to validate the reliability and accuracy of the LB model.The values of the parameters were chosen from the number 41 experiment conducted by Keiller and Sleath (1976) (referred to as KS hereinafter).The rough bed was composed of a layer of spheres placed in a hexagonal pattern (See Fig. 2).The diameter of the sphere d was fixed to 6.95δ, where δ was the Stokes layer thickness.It should be noted that the maximum free-stream velocity U 0 and d were taken as the characteristic velocity and length separately.The other parameters were made dimensionless with them and they were represented by a superscript asterisk.There was a conversion between the physical units and LB units in terms of similarity law.They were all shown in Table 1, where L x , L y , and L z were the length, height and width of the computational domain separately, T was the oscillatory period, Re a = U 0 a / ν.The boundary layer flow is often generated by an oscillatory pressure gradient in the streamwise direction 0 sin() where P was the pressure, ω = 2π / T was the angular frequency.The pressure gradient was converted to an external harmonic mass force in LB.The boundary conditions in streamwise and spanwise directions were periodic that was equivalent to considering an infinite rough bed with the roughness elements placed in a hexagonal pattern.The free-slip condition was imposed at the upper boundary.The non-slip boundary condition was implemented on the bottom boundary as well as the particles.
The ensemble averaging procedure was introduced to compute the mean quantities.Take the streamwise velocity u for instance, it was defined as follows where N was the number of averaged period.It was worth noting that in the present simulation the rough bed was kept fixed and the fluid moved to and fro, while in the KS experiment the bed oscillated with respect to the fluid otherwise at rest.When comparing the experimental and numerical results, the quantity U 0 cos(ωt) was added to the LB results in order to partially account for the difference.
The time development of the dimensionless magnitude of the projection of the velocity vector on a vertical plane parallel to the oscillatory direction |<V * >| at different distances from a sphere crest was shown in Fig. 3. Phase was taken to be zero as the moment when the rough bed velocity was maximum in the experiment.It indicated that there were two maxima of |<V * >| near the rough bed during each half period.One was in phase with that of maximum free-stream velocity and the other took place close to flow reversal.It was in accordance with the observation of KS.What's more, it showed the maximum of the secondary peak was equal to 0.49U 0 , which was attained at a distance 0.12d (0.82δ) from the sphere crest and at a phase of 93 degrees.While the KS results showed a corresponding maximum equal to 0.49U 0 , at a distance 0.70δ from the crest and characterized by a phase of 88 degrees.Comparing the LB and the measured values of the maximum of |<V * >| and of the phase of the secondary peak, a good agreement was obtained.The variation of maxima of |<V * >| and their phases with distance from a sphere crest was illustrated by Fig. 4. The LB results were in overall agreement with the KS experimental data.It could also be seen from Fig. 4 that the corresponding phases showed a marked change in the vicinity of 0.084d (0.6δ).Below this point, the phase value of the maximum |<V * >| was relative small.Above it, the maximum |<V * >| happened at around the phase of 90 degrees.The phase of this peak remained almost constant with increasing distance which was different from that close to the bed.KS also observed the similar phenomenon.The FV results were presented in Fig. 4 as well.It should be mentioned that the overall geometry of the bed in FV simulation was similar to the experiment but the details were different.The magnitude of the maximum |<V * >| was affected more by the difference.

SIMULATION SETUP
Once the reliability of the LB model had been tested, the numerical simulation of the oscillatory boundary layer over a rough bed would be further carried out.The rough bed was still composed of a layer of spherical particles placed regularly in two patterns.One was the hexagonal packing which was shown in Fig. 2. The other was the cubic packing which could be seen in Fig. 5.The choice of such regular roughness elements resulted from a compromise between computational efficiency and reproduction of a realistic geometry.The computational parameters were listed in Table 2.It should be mentioned that dimensional quantities are nondimensionalized and reported using lattice units.In terms of Jonsson's work (1966), the flows were all in the rough turbulent regime.The same boundary conditions mentioned above were imposed and DNS was carried out in the computation.

RESULTS AND DISCUSSION
The equivalent roughness height and the friction factor could be obtained from the log-fit method, which was based on Eq. 9.The term on the left-hand side was the ensemble-and space-averaged horizontal velocity.u * was the friction velocity, κ = 0.4, y t was the distance from the theoretical bed.The least squares method was used to determined the unknown parameters in the expression.

Equivalent Roughness Height
Fig. 6 showed the ensemble-and space-averaged velocity profiles of case c1 at different phases in semi-logarithmic graph.Phase was taken to be zero when the flow reversed.In order to ascertain the boundary layer developed quite substantially, the results with the phase between 50 and 100 degrees in the first half of period, and from 230 to 280 in the second half were chosen.The velocity profiles for other cases would not be given for the page limit.
Table 3 showed the time-averaged results of dimensionless distance y '* from the theoretical bed to the bottom boundary and the equivalent roughness height k s * , where Sd y'* and Sd ks* were the standard deviations of y '* and k s * separately.For both patterns, it can tell that the y '* as well as the k s * varied little around its average value in the oscillatory period according to the standard deviation, especially for the former one.The values of y '* in all cases were less than 1.0.It meant that the theoretical bed located at 0.19 ~ 0.25 times diameter below the crests of spherical particles.It was in good agreement with the experimental data reported by Dixen et al. (2008).The equivalent roughness heights were in the range of 2.51 ~ 3.41 times diameter but in most cases they were nearly 2.8 times diameter.They were close to the recommended value of 2.5 times diameter.

Friction Factor
The relationship between friction factor and friction velocity could be expressed by Eq. 10 (Lundgren and Jonsson 1961).The two maxima of dimensionless friction velocity u *mc / U 0 and u *mt / U 0 in both first and second half of period were given separately with the corresponding phase φ c and φ t in Table 4.The magnitude of the two maxima was close to each other and the larger one would be used to calculate the friction factor for each case.It was also shown that there existed a phase difference between the maximum friction velocity u *m and U 0 as the free-stream velocity reached its maximum U 0 at a phase of 90 or 270 degrees.The maximum friction velocity u *m led over U 0 , similar to those reported in earlier studies.Fig. 8 illustrated the phase lead ∆φ of u *m over U 0 .Although there was a considerable scatter in the data, it seemed that the ∆φ did not change markedly with a / k s .The LB results as well as experimental data were in the range of 10 ~ 30 degrees.When the a / k s decreased, especially when it was less than 1.0, the phase lead appeared a decline trend from its maximum 30 degrees to 20 degrees.The friction factor f w obtained from Eq. 10 was shown in Table 5.It could be seen that the f w mainly depended on the a / k s which was in good agreement with Jonsson's results (1966) in the rough turbulent regime.It was also not sensitive to the packing pattern.It should be paid more attention to the results that the f w did not seem to approach a constant value as suggested by Jonsson et al. but showed an increase when the a / k s reached small values, such as the case c4, c5, c10.There existed some expressions to represent such a behavior of f w for small values of a / k s .For example, Eq. 11 by Kamphuis (1975), Eq. 12 by Simons et al. (2000), Eq. 13 by Dixen et al. (2008).Eq. 14 was the one we got on the basis of LB results.It had a similar form with other expressions.
The LB results as well as experimental data were plotted in Fig. 9.The LB results were generally consistent with the existing data.Bagnold's data for small values of a / k s called for special attention.They deviated from the data of others and tended to a constant value.That may be the reason why Jonsson (1966) and Kajiura (1968) made such a suggestion about the tendency of the f w which was based on Bagnold's data.Eq. 11 to Eq. 14 were also shown in Fig. 9.It could be seen that Eq. 11 overpredicted the f w for the range indicated in the figure.Eq. 12 and Eq. 13 were close to each other.The Eq. 14 we got agreed well with Eq. 12 and Eq. 13 when the a / k s was less than 1.0 as well as most of the experimental data, especially for the one reported in Simons et al. (2000).

CONCLUSIONS
The 3-D lattice Boltzmann method was employed to characterize roughness effects of oscillatory boundary layer flow over a sphere-covered bed.The flow resistance including the equivalent roughness height and friction factor was investigated in the paper.The main conclusions were summarized as follows.
1.The benchmark problem of oscillatory boundary layer over a rough bed at low Reynolds number was solved.It indicated the LB model was feasible to this kind of problem from the mesoscale view.2. Direct numerical simulation of rough turbulent oscillatory boundary layer over a rough bed was carried out.The log-fit method based on the least squares was employed to analyze the position of the theoretical bed, equivalent roughness height and friction velocity.It could give reasonable results.3. The rough bed in computational domain was composed of a layer of spherical particles placed regularly in two patterns.One was cubic packing and the other was hexagonal packing.For both of them, it showed that the theoretical bed located at 0.19 ~ 0.25 times diameter below the crests of spherical particles.The dimensionless equivalent roughness height was nearly 2.8 in most cases, which was in good agreement with the recommended value of 2.5.4. The log-fit results indicated that the dimensionless friction velocity for the fixed bed appeared a sinusoidal-like behavior in the oscillatory period.The maximum friction velocity led over the maximum free-stream velocity.The phase lead was found to be in the range 10 ~ 30 degrees for the computational values of relative roughness a / k s = 0.47 ~ 6.58.The friction factor for small values of a / k s did not seem to tend to a constant value suggested by some previous investigators, but constantly increased with decreasing a / k s .

Figure 2 .
Figure 2. Left: Sketch of the computational domain Right: Hexagonal packing (top view).

Figure 3 .Figure 4 .
Figure 3. |<V * >| at different phases and distances from a sphere crest on the vertical plane.

Figure 5 .
Figure 5. Left: Sketch of the computational domain Right: Cubic packing (top view).

Figure 9 .
Figure 9. fw for small values of a / ks.