Direct Numerical Simulations of the Pulsating Flow over a Plane Wall

Abstract

The results of direct numerical simulations of the flow generated in a plane duct by a pressure gradient which is the sum of two terms are described. The first term of the pressure gradient is constant in space but it oscillates in time whereas the second term is constant both in space and in time. Therefore, a pulsating flow is generated, similar to that generated at the bottom of a monochromatic propagating surface wave when nonlinear effects are taken into account. The simulations are carried out for values of the parameters similar to those considered in previous investigations. It is shown that even a small constant pressure gradient influences the flow regime in the bottom boundary layer. In particular, turbulence strength is damped when the steady velocity component has the direction opposite to the oscillating velocity component whereas turbulence strength increases when the steady and oscillating components point in the same direction. Even though the flow is not exactly equal to that generated at the bottom of sea waves, where second order effects in the wave steepness induce a steady streaming in the direction of wave propagation, our results provide information on the interaction of the steady streaming with the oscillatory flow and are also relevant for investigating the dynamics of sediment close to the sea bottom. Indeed, since the turbulent eddies tend to pick-up the sediment from the bottom, it can be inferred that the triggering of turbulence enhances sediment transport towards the shore.

1. Introduction

A fair approximation of the flow field generated by a propagating surface wave is obtained by assuming the flow to be irrotational and the fluid to be inviscid. However, this approach cannot describe the flow close to the bottom where vorticity is generated because of the no-slip condition and viscous effects turn out to be significant. Hence, the flow provided by the irrotational and inviscid approach should be corrected by introducing a bottom boundary layer. When surface waves of small amplitude are considered, the flow field can be expanded as a power series of the wave steepness and, at the leading order of approximation, the flow close to the bottom is described by the well-known Stokes’ solution that considers the laminar regime [1].

However, turbulence in an oscillatory boundary layer (Stokes boundary layer) appears when the Reynolds number is larger than a critical value falling approximately between 500 and 600. Hereinafter, the Reynolds number $R_{\delta }$ is defined using the amplitude $U_0^{*}$ and the angular frequency ${\omega }^{*}$ of the oscillations of the irrotational velocity evaluated at the bottom along with the kinematic viscosity ${\nu }^{*}$ of the seawater:

$$R_{\delta }=\frac{U_0^{*}{\delta }^{*}}{{\nu }^{*}}.$$

(1)

In (1), ${\delta }^{*}=\sqrt{\frac{2{\nu }^{*}}{{\omega }^{*}}}$ denotes the order of magnitude of the thickness of the viscous bottom boundary layer.

Indeed, the measurements of Hino et al. [2] show that the velocity field is characterized by the appearance of turbulent oscillations towards the end of the accelerating phases of the cycle when the Reynolds number is larger than a value approximately equal to 500. However, small perturbations of the Stokes flow are superimposed to the laminar flow already for values of $R_{\delta }$ falling around 100, but the average velocity profile exhibits only small deviations from the laminar solution. Only when the Reynolds number becomes larger than about 500, turbulence becomes significant and it affects the flow field during larger parts of the oscillatory cycle as the Reynolds number is increased. The fully turbulent regime is observed, characterized by turbulence presence throughout the whole cycle, when the Reynolds number becomes larger than a value that is difficult to quantify exactly but the results of Jensen et al. [3] suggest to range between 2500 and 3500. Therefore, the experimental measurements show the existence of four flow regimes, i.e., the laminar regime, the disturbed laminar regime, the intermittently turbulent regime and the fully developed turbulent regime. The disturbed laminar regime is characterized by the presence of small amplitude perturbations superimposed on the laminar flow. The intermittently turbulent regime is characterized by turbulence bursts which appear during the decelerating phases of the cycle but the flow recovers a laminar like behavior during the remaining parts of the cycle. The fully turbulent regime is characterized by turbulence presence throughout the whole cycle.

The theoretical and numerical investigations of both the transition process and the dynamics of vorticity in the oscillatory flow at the bottom of the sea waves are difficult because of the unsteadiness of the basic flow [4,5,6,7,8,9,10,11,12,13,14] and because of the presence of bedforms [15,16,17]. When the bottom is flat, recent theoretical investigations and numerical simulations [18,19,20,21,22,23,24,25], show that turbulence appearance is triggered by either wall imperfections [18,20] or by the bottom roughness, because of the vortex structures shed by the roughness elements [23,25] or by external perturbations [12,26].

Even though the oscillatory boundary layer generated close to a wall by the harmonic oscillations of a fluid has been deeply investigated, much less is known on the boundary layer generated at the bottom of sea waves. Indeed, no accurate measurements have been made because of the difficulties to generate, in a laboratory, propagating surface waves characterized by large values of the Reynolds number of the bottom boundary layer.

The numerical simulations are difficult to be made because of the different length scales involved in the phenomenon. Indeed, the thickness ${\delta }^{*}$ of the bottom boundary layer is of the order of ${10}^{-3}$ m while the length of the surface waves is of order ${10}^2$ m. Taking into account that at least 10 grid points or more are needed every ${\delta }^{*}$, it turns out that the numerical simulations would require a number of grid points of order ${10}^{15}$.

The phenomenon becomes more complex when the surface waves interact with steady currents, since the average flow at time $t^{*}$ is no longer the mirror image of that at time $t^{*}+T^{*}/2$ and turbulence production during half a cycle becomes different from that during the previous or following half cycle.

Herein, to provide some information on the effect that a steady velocity component has on the transition process and on turbulence dynamics, direct numerical simulations are made of the flow generated in a plane duct by a pressure gradient which is the sum of two terms. The first term of the pressure gradient is constant in space but it oscillates periodically in time whereas the second term is constant both in space and in time. Hence, a pulsating flow is generated as that considered by [10,27,28]. The numerical simulations are carried out for values of the parameters similar to those considered by [18,19].

Since both theoretical analyses and previous numerical simulations show that wall imperfections play a key role in triggering turbulence appearance, in the present numerical simulations the wall is not perfectly plane but characterized by the presence of a small waviness, the amplitude of which is so small that the wall is flat from a macroscopic point of view.

The structure of the rest of the paper is the following: in the next section, we formulate the problem and briefly describe the numerical approach. The third section is devoted to describing the results and in particular the velocity field, the spatial and temporal distribution of the turbulent kinetic energy and the time development of the bottom shear stress. The last section is devoted to the conclusions.

2. Problem

An incompressible viscous fluid, bounded by a horizontal wall, is set in motion by a pressure gradient described by

$$\frac{\partial p^{*}}{\partial x_1^{*}}=-{\rho }^{*}{U}_0^{*}{\omega }^{*}\left(sin\left({\omega }^{*}{t}^{*}\right)+C\right),\frac{\partial p^{*}}{\partial x_2^{*}}=0,\frac{\partial p^{*}}{\partial x_3^{*}}=0,$$

(2)

where $x_1^{*},x_2^{*}$ and $x_3^{*}$ are the streamwise, cross-stream and spanwise coordinates of a Cartesian coordinate system, ${\rho }^{*}$ denotes the fluid density, $U_0$ and $T^{*}=2\pi /{\omega }^{*}$ are related to the amplitude and the period of the fluid velocity oscillations induced by the oscillating pressure far from the wall and C is an arbitrary constant. The rigid wall, located at $x_3^{*}=0$, is not perfectly flat but is characterized by a waviness of small amplitude as in [18,19]. More in detail, the wall profile is described by

$$x_2^{*}={ϵ}^{*}\eta (x_1^{*},x_3^{*})={ϵ}^{*}\sum _{n=1}^N{a}_n^{*}cos\left({\alpha }_n^{*}{x}_1^{*}+{\gamma }_n^{*}{x}_3^{*}+{\varphi }_n\right),$$

(3)

where $a_n$ is the amplitude of the nth component of the wall waviness characterized by wavenumbers ${\alpha }_n^{*}$ and ${\gamma }_n^{*}$ in the $x_1^{*}$- and $x_3^{*}$-directions respectively and by a random phase. Moreover, the flow is assumed to be symmetric with respect to the plane $x_2^{*}=L_{x2}^{*}$.

Hereinafter, a star denotes dimensional quantities and the same symbols without the star denote their dimensionless counterparts which are defined by

$$t=t^{*}{\omega }^{*},(x_1,x_2,x_3)=\frac{(x_1^{*},x_2^{*},x_3^{*})}{{\delta }^{*}},(u_1,u_2,u_3)=\frac{(u_1^{*},u_2^{*},u_3^{*})}{U_0^{*}},p=\frac{p^{*}}{{\rho }^{*}{U}_0^{*2}}.$$

(4)

In (4), $t^{*}$ is time and $u_1^{*},u_2^{*},u_3^{*}$ are the fluid velocity components. A value of $L_{x2}^{*}$ much larger than ${\delta }^{*}$ is chosen so that the ghost wall located at $x_2^{*}$ equal to $2L_{x2}^{*}$ does not affect the flow close to the wall described by (3).

The hydrodynamic problem is posed by Navier–Stokes and continuity equations

$$\frac{\partial \mathbf{u}}{\partial t}+\frac{R_{\delta }}{2}\nabla ·\left(\mathbf{u}\mathbf{u}\right)=-\frac{R_{\delta }}{2}\nabla p+i_{x1}\left(sin\left(t\right)+C\right)+\frac{1}{2}{\nabla }^2\mathbf{u},$$

(5)

where $i_{x1}$ denotes the unit vector in the streamwise direction and $R_{\delta }=\frac{U_0^{*}{\delta }^{*}}{{\nu }^{*}}$ is the Reynolds number of the flow.

The equations are solved by a numerical approach in a computational domain of size $L_{x1},L_{x2}$ and $L_{x3}$ in the streamwise, vertical and spanwise directions, respectively. The numerical approach uses finite difference approximations of the spatial and temporal derivatives. An exhaustive description of the numerical approach can be found in [19] where the tests carried out to validate the numerical code are also described.

The no-slip condition is enforced at the wall

$$\left(u_1,u_2,u_3\right)=0\mathrm{at}{x}_2=ϵ\eta \left(x_1,x_3\right).$$

(6)

Since the value of ${ϵ}^{*}$ is assumed to be much smaller than ${\delta }^{*}$, the parameter $ϵ$ turns out to be much smaller than 1 and the boundary condition at the wall can be approximated by

$$\left(u_1,u_2,u_3\right)=-ϵ\eta \frac{\partial \left(u_1,u_2,u_3\right)}{\partial x_2},\mathrm{at}{x}_2=0$$

(7)

since the accuracy of (7) is within the accuracy of the numerical approach itself. Indeed the numerical approach is second order accurate in space and in the simulations the value of $ϵ$ is smaller than the size of the first computational step in the $x_2^{*}$-direction. Moreover, at $x_2^{*}=h_0^{*}$, a symmetry condition is enforced, which is equivalent to enforce the vanishing of tangential stresses far from the wall.

The computational grid is uniform in the streamwise and spanwise directions, with $n_{x1}$ and $n_{x3}$ grid points, respectively. In the cross-stream direction a non-uniform mesh with $n_{x2}$ points is used to cluster the gridpoints close the wall where velocity gradients are expected to be larger.

Costamagna et al. [29], for a flow similar to that investigated here and for similar values of the Reynold number, discussed the adequacy of a computational domain with size $L_{x1}=25.13{\delta }^{*}$, $L_{x2}=25.13{\delta }^{*}$, $L_{x2}=12.57{\delta }^{*}$ and number of gridpoints equal to $n_{x1}=64$, $n_{x2}=64$, $n_{x3}=32$ (small box). Considering the spatial autocorrelation functions and comparing with results obtained in control runs made with a larger computational domain (big box, $L_{x1}=50.27{\delta }^{*}$, $L_{x2}=25.13{\delta }^{*}$, $L_{x2}=25.14{\delta }^{*}$, $n_{x1}=192$, $n_{x2}=64$, $n_{x3}=96$), Costamagna et al. [29] concluded that as in [30], the small computational domain does not allow accurate predictions of the high-order statistical quantities, but it is adequate to isolate the basic process generating turbulence. On the other hand, Costamagna et al. [29] showed that the big box allows us to investigate turbulence structure both near and far from the wall. Moreover Costamagna et al. [29], by considering power spectra of the velocity, showed that the number of grid points is adequate to resolve the smallest scales of the flow, even at the phases of the cycle when turbulence appears.

Considering the size of the computational box and the number of grid points used in the present simulations (see

Table 1. Values of the hydrodynamics parameters of the numerical simulations.

${\mathit{R}}_{\mathit{\delta }}$	C	$\frac{{\mathit{L}}_{\mathit{x}1}^{*}}{{\mathit{\delta }}^{*}}$	$\frac{{\mathit{L}}_{\mathit{x}2}^{*}}{\mathit{\delta }}$	$\frac{{\mathit{L}}_{\mathit{x}3}^{*}}{\mathit{\delta }}$	${\mathit{n}}_{\mathit{x}1}$	${\mathit{n}}_{\mathit{x}2}$	${\mathit{n}}_{\mathit{x}3}$
100	${10}^{-3}$	50.27	25.13	25.13	129	65	129
500	0	50.27	25.13	25.13	129	65	129
500	${10}^{-3}$	50.27	25.13	25.13	129	65	129
600	0	50.27	25.13	25.13	129	65	129
600	${10}^{-3}$	50.27	25.13	25.13	129	65	129
600	$5\times {10}^{-3}$	50.27	25.13	25.13	129	65	129
600	${10}^{-2}$	50.27	25.13	25.13	129	65	129
800	0	50.27	25.13	25.13	129	65	129
800	${10}^{-3}$	50.27	25.13	25.13	129	130	129
800	$5\times {10}^{-3}$	50.27	25.13	25.13	129	130	129
800	${10}^{-2}$	50.27	25.13	25.13	129	130	129
1000	0	50.27	25.13	25.13	129	65	129
1000	${10}^{-3}$	50.27	25.13	25.13	129	130	129
1000	$5\times {10}^{-3}$	50.27	25.13	25.13	129	130	129
1000	${10}^{-2}$	50.27	25.13	25.13	129	130	129
1200	0	50.27	25.13	25.13	129	65	129
1200	${10}^{-3}$	50.27	25.13	25.13	129	130	129
1200	$5\times {10}^{-3}$	50.27	25.13	25.13	129	130	129
1200	${10}^{-2}$	50.27	25.13	25.13	129	130	129

), which are comparable to those of the big box by [29], it appears that the numerical parameters used in the present simulations are adequate to describe the average flow. The initial condition is either the steady velocity profile, induced by the steady pressure gradient or the final flow field of a run with a similar Reynolds number. In any case, the simulations are carried out till a regime configuration is attained and the results do not depend on the initial condition. The numerical code was run for different values of the Reynolds number $R_{\delta }$ ranging from 100 to 1200 and for values of the steady pressure gradient of order ${10}^{-2}$ (see Table 1).

The parameters describing the wall imperfections are the same as those used by [19]. In particular, a wall profile with two harmonic components ($N=2$) is used and the first component of the wall waviness is two-dimensional ($a_1=1.0,{\alpha }_1=0.5,{\gamma }_1=0,{\varphi }_1=0$) whereas the second component is a three-dimensional undulation ($a_2=0.1,{\alpha }_2=0,{\gamma }_2=1.0,{\varphi }_2=0$). Moreover, $ϵ$ is fixed equal to $0.005$ which is the order of magnitude of the imperfections of a mirror-shine smooth wall. The numerical code was extensively validated by [19,29] by comparing the results it provides with analytical solutions, available for low Reynolds numbers, and with experimental data.

3. Results

In a purely oscillatory boundary layer, for values of the Reynolds number larger than a first critical value ranging between 500 and 600, the transition to turbulence takes place every half cycle and it is followed by a phase during which turbulence decays and the flow recovers a laminar-like behaviour [19]. This scenario is partly modified by the presence of a constant positive pressure gradient, even though turbulence dynamics is similar to that observed in [19].

In Figure 1, the instantaneous spanwise velocity component in the horizontal plane $x_2^{*}=0.65{\delta }^{*}$ is plotted at two different phases of the cycle for $R_{\delta }=800$ and $C=1.0\times {10}^{-3}$. The results clearly show that turbulence intensity is rather small at $t^{*}\simeq 19.26T^{*}$ whereas it is significant at $t^{*}\simeq 19.51T^{*}$ (the reader should notice the different scales on the right of the panels of Figure 1).

Moreover, the time development of the turbulent kinetic energy ($K^{*}$), averaged over horizontal planes and integrated into the vertical direction, shows that the pressure gradient has different effects on the flow, depending on the direction of the oscillatory velocity component with respect to the direction of the current induced by the steady pressure gradient. Indeed, as shown in Figure 2 for $R_{\delta }=1000$, whereas in purely oscillatory flow $K^{*}$ has two similar peaks per cycle, the presence of a moderate steady current decreases the intensity of one of the maxima. In particular, the maximum located in the time interval during which the oscillatory velocity component opposes the steady current is damped. On the other hand, when the oscillatory velocity component has the same direction as the steady current, the turbulence strength increases with respect to the case characterized by $C=0$. If the steady current increases in intensity ($C={10}^{-2}$), only one maximum per cycle is observed. In Figure 2, K is equal to $K^{*}$ divided by ${\rho }^{*}{U}_0^{*2}{\delta }^{*}$ where $K^{*}$ is the turbulent kinetic energy per unit bottom area, averaged over the horizontal plane and integrated over the vertical direction.

It is worthwhile to mention that K does not appear to reach a periodic behavior because of the size of the computational box which is not large enough to provide an accurate description of the dynamics of the largest vortex structures. Indeed, the size of the computational domain and the number of grid points are the results of a compromise between the accuracy of the numerical results and the computational costs. The evaluation of the spectra of the velocity field at different phases during the oscillatory cycle shows that the grid size is small enough to capture the dynamics of the smallest turbulent eddies but the size of the computational domain does not allow to describe accurately the largest vortex structures. However, the energy of these large vortices is rather small and we are confident that the results described in the paper provide a reliable description of turbulence dynamics.

For values of the Reynolds number close to the first critical value for the transition from the disturbed laminar to the intermittently turbulent regime in the purely oscillatory boundary layer [19], the steady current appears to have a different effect of K.

Figure 3 shows that for $R_{\delta }=600$ and $C=0$ two maxima of K appear during one oscillating cycle but for large values of C, only one maximum of the turbulent kinetic energy can be clearly observed. The phase within the cycle of the maximum of K depends significantly on the value of C and it is difficult to draw definitive conclusions. However, it appears that for moderate values of the Reynolds number, the opposite current tends to stabilize the flow for all the values of C considered and the flow tends to recover the laminar regime. On the other hand, the turbulent kinetic energy increases when the oscillatory velocity component points in the same direction as the steady velocity.

The vertical distribution of the plane-averaged turbulent kinetic energy $k^{*}$, plotted in Figure 4 and Figure 5 for $R_{\delta }=600$ and $R_{\delta }=1000$ respectively, shows that turbulence reaches larger distances when the flow has the same direction as the current. The intensity of the steady current does not appear to alter the vertical profile of the turbulent kinetic energy from a qualitative point of view. The turbulent kinetic energy profiles for $R_{\delta }=800$ and $R_{\delta }=1200$ are qualitatively similar to those of $R_{\delta }=1000$.

Figure 6 shows the time development of the bottom shear stress. Vittori and Verzicco [19] showed that, for the smallest value of $R_{\delta }$ and $C=0$, the bottom shear stress significantly deviates from the values typical of the laminar regime just after its maximum and minimum values, when turbulence is generated. Of course when the Reynolds number is increased, turbulence becomes stronger and appears earlier within the cycle and the peaks of the bottom shear stress become larger. Consistently with the time evolution of the turbulent kinetic energy, when values of C different from zero are considered, the peaks of the bottom shear stress in the direction of the steady current become larger whereas those in the opposite direction decrease.

So far, we have described the influence of the steady pressure gradient on the oscillatory components of the flow. However, the presence of a steady pressure gradient generates a steady velocity component too. If no oscillatory velocity component is present ($U_0^{*}=0$), in the laminar regime the velocity profile can be easily evaluated. Such a velocity profile was used as the initial condition for the numerical simulations. After a large number of cycles, the velocity profile, averaged over horizontal planes and over one period should not change from one cycle to the other. Presently, due to the limited statistical sample, the velocity profile shows little difference form one cycle to the following ones. Figure 7 shows the time-averaged velocity profiles obtained by considering many cycles. As expected, it can be seen that the steady velocity component $u_1$ is always negative and the values of the velocity increase for increasing values of C. It is also worthwhile to notice that, as the Reynolds number is increased, the intensity of the steady velocity component tends to decrease.

4. Conclusions

The direct numerical simulations of the flow generated in a plane duct by a pressure gradient which is the sum of a constant pressure gradient plus a temporal oscillating contribution show that, for moderate values of the Reynolds number, the turbulence strength, evaluated when the oscillating velocity component points in the same direction as the steady velocity component, is stronger than that observed when the two velocity components have opposite directions. At the bottom of sea waves a steady velocity component, pointing onshore, is generated by nonlinear effects [31]. Even though the flow provided by the present direct numerical simulations is not exactly equal to that generated at the bottom of sea waves, the obtained results provide information on the effect that the steady streaming [31], that is present in the boundary layer at the bottom of a sea wave, has on the transition to turbulence and on the sediment transport. Indeed, it can be inferred that, the turbulence strength observed when the oscillatory velocity points onshore is stronger than when it points offshore. Because the sediment transport is significantly increased by turbulence presence, it appears that the onshore sediment transport rate is significantly increased by turbulence appearance.