The role of wave kinematics in turbulent flow over waves

TL;DR

This study investigates turbulent flow over waves using LES, analyzing wave kinematics' impact on form drag across different wave ages and Reynolds numbers, and introduces a novel split system approach for modeling.

Contribution

A new split system method is developed to separate laminar and turbulent responses, improving understanding of wave-induced stresses in turbulent flows.

Findings

Form drag sensitivity to Reynolds number decreases at high wave ages.

Wave kinematics significantly influence turbulent stresses and form drag.

Fundamental Reynolds stress mode becomes more critical at higher Reynolds numbers.

Abstract

The turbulent flow over monochromatic waves of moderate steepness is studied by means of wall resolved LES, and cover a range of wave ages and Reynolds numbers. We compute the Fourier modes of the flow and analyse the momentum balance for the mean and fundamental mode. At low wave ages, the form drag displays a large sensitivity to changes in Reynolds number, and the interaction between turbulent and wave-induced stresses increases with Reynolds number. At higher wave ages, a quasi-laminar regime is entered, where the wave-induced stress is primarily balanced by viscous stresses. To exploit the increasing importance of the wave kinematics observed in the intermediate to high wave age regime, a novel split system approach is introduced, where a laminar response to the wave forces a turbulent RANS type of shear flow . To account for the effects of turbulence, we force the system using Reynolds stresses from the corresponding LES. We give an analytic functional dependence for the form drag associated with the laminar solution. For intermediate to high wave ages, the form drag of the shear flow exhibits relatively simple behaviour, and we derive approximate functional dependencies for the quasi-laminar regime. The high sensitivity of the form drag to variation in Reynolds number at low wave ages is more challenging to evaluate using this approach. This is primarily due to the increased importance of nonlinearity in the shear flow, which are inherently coupled with higher harmonics in the turbulent stresses. Nevertheless, the split system approach can be utilized to quantify the importance of different harmonics in the turbulent stresses by explicitly choosing which modes to include in the split system forcing. We demonstrate that the fundamental mode of the Reynolds stresses becomes increasingly important for accurate prediction of the form drag as the Reynolds number increases.

Tables

Table 1. Table 1: Computation cases. All cases use a geometric stretching in the vertical direction of r = 1.03 𝑟 1.03 r=1.03 and the spacing from the surface to the first point in the domain is Δ ​ ξ min ​ u τ / ν = 0.7 Δ subscript 𝜉 min subscript 𝑢 𝜏 𝜈 0.7 \Delta\xi_{\mathrm{min}}u_{\tau}/\nu=0.7 .

$R{e}_{\tau }$	$\Delta x{u}_{\tau }/\nu$	$\Delta z{u}_{\tau }/\nu$	$\Delta {\xi }_{\max }{u}_{\tau }/\nu$	$N_x$	$N_z$	$N_y$	$c/u_{\tau }$
200	8.0	8.0	7.1	150	75	70	0,2,4,8,12,16,24,36
260	8.7	7.8	8.1	200	100	85	4
395	9.5	7.9	12.4	250	150	94	0,2,4,8,12,16,24,36
950	9.5	9.5	27.4	600	300	130	2,4,8,12,16,24

Equations

$$(u,v,w{)}_{y=\eta (x)}=(akc\sin (kx)-c,-akc\cos (kx),0).$$

$${\partial }_t{u}_i+u_j{\partial }_j{u}_i$$

$${\partial }_i{u}_i$$

$$u_i(x_j,t)=⟨u_i(x_j,t)⟩+u_i^{\prime }(x_j,t)=U_i(x,y)+u_i^{\prime }(x_j,t),p(x_j,t)=P(x,y)+p^{\prime }(x_j,t).$$

$$(u-c){\partial }_{x}u+v{\partial }_{y}u$$

$$(u-c){\partial }_{x}v+v{\partial }_{y}v$$

$${\partial }_{x}u+{\partial }_{y}v=0,$$

$$\chi =x,\xi =y-g(\chi ,\xi )=y-\frac{ai}{2}(\exp (-ik\chi )-\exp (ik\chi ))\exp (-k\xi ).$$

$$(u-c){\partial }_{\chi }u+v{\partial }_{\xi }u$$

$$(u-c){\partial }_{\chi }v+v{\partial }_{\xi }v$$

$${\partial }_{\chi }u+{\partial }_{\xi }v$$

$${\nabla }_{\chi }^2=\frac{{\partial }^2}{\partial {\chi }^2}+\frac{{\partial }^2}{\partial {\xi }^2},{\mathcal{L}}_N=g_{\chi }^2\frac{{\partial }^2}{\partial {\xi }^2}-2g_{\chi }\frac{{\partial }^2}{\partial \chi \partial \xi }+(g_{\chi }{g}_{\chi \xi }-g_{\chi \chi })\frac{\partial }{\partial \xi },$$

$${\mathcal{L}}_u=(u-c)g_{\chi }{\partial }_{\xi }u+g_{\chi }{\partial }_{\xi }p+g_{\chi }{\partial }_{\xi }⟨u^{\prime }{u}^{\prime }⟩,{\mathcal{L}}_v=(u-c)g_{\chi }{\partial }_{\xi }v+g_{\chi }{\partial }_{\xi }⟨u^{\prime }{v}^{\prime }⟩.$$

$$f(\chi ,\xi )=\sum _{m=-M}^M{\hat{f}}_m(\xi )\exp (ikm\chi ),{\hat{f}}_{-m}={\hat{f}}_m^{\ast },$$

$$\mathcal{F}(f{\partial }_{x}h)\simeq \mathcal{F}(f{\partial }_{\chi }h-f{g}_{\chi }{\partial }_{\xi }h)=\mathcal{F}(f)\ast (\mathcal{F}({\partial }_{\chi }h)-\mathcal{F}(g_{\chi }{\partial }_{\xi }h)).$$

$$\mathcal{F}(g_{\chi }{\partial }_{\xi }h{)}_m=\frac{ak}{2}\exp (-k\xi )\mathcal{D}({\hat{h}}_{m-1}+{\hat{h}}_{m+1})={\mathcal{D}}_g({\hat{h}}_{m-1}+{\hat{h}}_{m+1})={\mathcal{D}}_g[{\hat{h}}_m].$$

$$ikm{\hat{u}}_m-{\mathcal{D}}_g[{\hat{u}}_m]+\mathcal{D}{\hat{v}}_m=0,$$

$$\begin{gathered} \mathcal{F}(f{\partial }_{\chi }h-f{g}_{\chi }{\partial }_{\xi }h{)}_m={\hat{f}}_0\left(imk{\hat{h}}_m-{\mathcal{D}}_g[{\hat{h}}_m]\right)+ \\ {\hat{f}}_{-1}\left(i(m+1)k{\hat{h}}_{m+1}-{\mathcal{D}}_g[{\hat{h}}_{m+1}]\right)+{\hat{f}}_1\left(i(m-1)k{\hat{h}}_{m-1}-{\mathcal{D}}_g[{\hat{h}}_{m-1}]\right)+\dots \end{gathered}$$

$$ikm({\widehat{uu}}_m-c{\hat{u}}_m)+\mathcal{D}{\widehat{uv}}_m$$

$${\hat{u}}_{\tau }^2/h{\delta }_{0m}+{\mathcal{D}}_g[{\widehat{⟨u^{\prime }{u}^{\prime }⟩}}_m+{\hat{p}}_m+{\widehat{uu}}_m-c{\hat{u}}_m]$$

$$ikm({\widehat{uv}}_m-c{\hat{v}}_m)+\mathcal{D}{\widehat{vv}}_m$$

$${\mathcal{D}}_g\left[{\widehat{⟨u^{\prime }{v}^{\prime }⟩}}_m+{\widehat{uv}}_m-c{\hat{v}}_m\right]$$

$${\widehat{uv}}_m={\hat{u}}_0{\hat{v}}_m+{\hat{u}}_{-1}{\hat{v}}_{m+1}+{\hat{u}}_1{\hat{v}}_{m-1}+{\hat{u}}_{-2}{\hat{v}}_{m+2}+{\hat{u}}_2{\hat{v}}_{m-2}+\dots$$

$$x=\chi ,y=\eta (\chi )+\xi (1-\eta (\chi )/h),z=\zeta ,$$

$$\mathcal{D}{\widehat{uv}}_0-2{\mathcal{D}}_g({\widehat{uu}}_1^r-c{\hat{u}}_1^r)=\mathcal{D}(-{\widehat{⟨u^{\prime }{v}^{\prime }⟩}}_0+\nu \mathcal{D}{\hat{u}}_0)+{\hat{u}}_{\tau }^2/h+2{\mathcal{D}}_g({\widehat{⟨u^{\prime }{u}^{\prime }⟩}}_1^r+{\hat{p}}_1^r),$$

$$u_{\tau }^2=\frac{1}{\lambda }{\int }_0^{\lambda }\left(\nu ({\partial }_{y}u+{\partial }_{x}v-2{\partial }_{x}u{\partial }_x\eta )+p{\partial }_x\eta \right)\mathrm{d}x=D_v+D_p.$$

$${\eta }_x=0.5ak(\exp (ikx)+\exp (-ikx)),$$

$$D_p={\lambda }^{-1}{\int }_0^{\lambda }{\eta }_{x}p\mathrm{d}x=\frac{1}{2}ak({\hat{p}}_1+{\hat{p}}_{-1})=ak{\hat{p}}_1^r.$$

$${\hat{p}}_1^r\approx {\left\{{\int }_0^{1}ik({\hat{u}}_0-c){\hat{v}}_1\mathrm{d}\xi \right\}}_{re}={\int }_0^1-k({\hat{u}}_0-c){\hat{v}}_1^{im}\mathrm{d}\xi ,$$

$$\begin{gathered} \underset{\text{Wave induced}}{\underbrace{\mathcal{D}(c{\hat{v}}_1+{\widehat{uv}}_1)+ik{\widehat{uu}}_1-{\mathcal{D}}_g[{\widehat{uu}}_1]}}=\underset{\text{Viscous}}{\underbrace{\nu {\mathcal{D}}^2{\hat{u}}_1}}-\underset{\text{Turbulent}}{\underbrace{(\mathcal{D}{\widehat{⟨u^{\prime }{v}^{\prime }⟩}}_1+ik{\widehat{⟨u^{\prime }{u}^{\prime }⟩}}_1-{\mathcal{D}}_g[{\widehat{⟨u^{\prime }{u}^{\prime }⟩}}_1])}}- \\ \underset{\text{Pressure}}{\underbrace{(ik{\hat{p}}_1+{\mathcal{D}}_g[{\hat{p}}_1])}}. \end{gathered}$$

$$\frac{c}{u_{\tau }}=\frac{c}{k\nu }\frac{\nu }{u_{\tau }\lambda }(2\pi )=\frac{R{e}_w}{R{e}_{\tau }}(2\pi ),$$

$$u_j{\partial }_j{u}_i=\underset{\text{particular}}{\underbrace{u_j^p{\partial }_j{u}_i^p}}+\underset{\text{homogeneous}}{\underbrace{u^p{\partial }_j{u}_i^h+u_j^h{\partial }_j{u}_i^p+u_j^h{\partial }_j{u}_i^h}}.$$

$$u_i^p=u_i^{\varphi }+u_i^{\psi },$$

$$-c{\partial }_x{u}_i^{\psi }+u_j^{\varphi }{\partial }_j{u}_i^{\psi }+u_j^{\psi }{\partial }_j{u}_i^{\varphi }+\underset{\text{nonlinear}}{\underbrace{u_j^{\psi }{\partial }_j{u}_i^{\psi }}}=-{\partial }_i{p}_{\psi }+\nu {\partial }_{jj}{u}_i^{\psi },{\partial }_i{u}_i^{\psi }=0,$$

Full text

To exploit the increasing importance of the wave kinematics observed in the intermediate to high wave age regime, a novel split system approach is introduced. The key ingredient in this approach is that in the absence of background shear, the flow response to the wave kinematics is laminar. This laminar solution acts as a forcing in a RANS type model for the shear flow that is subject to homogeneous boundary conditions. To account for the effects of turbulence, we force the system using Reynolds stresses from the corresponding large eddy simulation. We give an analytic functional dependence for the form drag associated with the laminar solution. For intermediate to high wave ages, the form drag of the shear flow exhibits relatively simple behaviour, and we derive approximate functional dependencies for the quasi-laminar regime.

The high sensitivity of the form drag to variation in Reynolds number at low wave ages is more challenging to evaluate using the split system approach. This is primarily due to the increased importance of nonlinearity in the shear flow. These nonlinearites are inherently coupled with higher harmonics in the turbulent stresses. Nevertheless, the split system approach can be utilized to quantify the importance of different harmonics in the turbulent stresses by explicitly choosing which modes to include in the split system forcing. We demonstrate that the fundamental mode of the Reynolds stresses becomes increasingly important for accurate prediction of the form drag as the Reynolds number increases.

1 Introduction

Stresses from turbulent wind over water generate waves, which in turn alter the wind. The exchange of momentum, energy and mass occurring at air-water interfaces is to a large extent determined by this wind-wave interaction. Despite receiving constant attention from the research community in the form of field measurements (Hasselmann et al., 1973; Hristov et al., 1998; Edson et al., 2007), laboratory experiments (Plant, 1982; Mastenbroek et al., 1996; Grare et al., 2013; Buckley & Veron, 2016), theoretical studies (Phillips, 1957; Miles, 1957; Janssen, 1991; Belcher & Hunt, 1993) and numerical simulations (Gent & Taylor, 1976; Van Duin & Janssen, 1992; Li et al., 2000; Meirink & Makin, 2000; Sullivan et al., 2000, 2008; Yang et al., 2013), there are still some very simple questions that remain only partly answered. One of these questions is how the form drag, exerted by the waves on the air, depends on the properties of the air flow and the waves.

Naturally, early attempts at numerical modelling of this flow were devoted to the use of Reynolds averaged Navier-Stokes (RANS) models for monochromatic waves (Gent & Taylor, 1976; Van Duin & Janssen, 1992). The results turned out to be sensitive to the choice of turbulence closure scheme (Van Duin & Janssen, 1992), and suitable closures were adopted (Mastenbroek et al., 1996; Li et al., 2000; Meirink & Makin, 2000). Increased computing capabilities led to the advent of direct numerical simulation (DNS) studies (Sullivan et al., 2000; Kihara et al., 2007; Yang & Shen, 2009, 2010). DNS studies have provided valuable insight into the shape of the wave correlated motion (Sullivan et al., 2000) and the structure of the turbulence (Yang & Shen, 2009, 2010). Due to the strict resolution requirements inherent to DNS, only low Reynolds numbers have been considered so far. Meirink & Makin (2000) used RANS models suitable for low Reynolds number flows to discuss the Reynolds number effects on the drag. They found reasonable agreement with Sullivan et al. (2000) at low Reynolds number, but further showed that there are considerable differences between the form drag obtained at low and intermediate Reynolds numbers. Large eddy simulation (LES) studies have primarily been aimed at closing the gap to operational models, with special attention given to subgrid near-surface modelling and the inclusion of a realistic wave spectrum (Sullivan et al., 2008; Yang et al., 2013; Hara & Sullivan, 2015). However, Yang et al. (2013) used wall resolved LES at an intermediate to high Reynolds number as a baseline for their discussions on subgrid modelling.

Waves that are slower than the bulk wind extract momentum from the wind through the action of the wave-correlated pressure against the slope of the surface. These waves are responsible for most of the momentum transfer from the atmosphere to the sea. In this part of the spectrum, breaking is frequent, which increases both drag and aerosol transport (Donelan, 1998). Operational wave forecasting relies on the celebrated critical layer mechanism of Miles (1957, 1959), where the discontinuity in wave-correlated Reynolds stresses at the critical layer is responsible for resonance. The existence of the critical layer, and the corresponding jump in wave-correlated Reynolds stresses, has been demonstrated both numerically (Sullivan et al., 2000) and in field observations (Hristov et al., 2003). However, Belcher & Hunt (1993, 1998) argued on the basis of a rapid distortion framework, that the Miles mechanism could only be active for relatively fast waves and reintroduced a non-separated version of the Jeffreys (1925) sheltering mechanism for slow waves.

Waves that propagate faster than the bulk wind typically emerge as the result of nonlinear interaction of slower waves, or as the result of distantly storm-generated waves that enter domains of lighter winds. The momentum flux between the air and the sea is dramatically altered in the presence of fast moving waves (Sullivan et al., 2008; Kahma et al., 2016). In severe situations, fast waves even have the ability to generate wind (see for instance Semedo et al., 2009). Examples of simplified models for the air-side response in this regime include Cohen & Belcher (1999); Kudryavtsev & Makin (2004); Semedo et al. (2009).

In the present work, we perform wall resolved LES over monochromatic waves for different wave ages and Reynolds numbers. With wall resolved LES the requirements for resolving the turbulence in an adequate manner are less severe than in DNS, and as such, it may be a valuable tool for studying the flow dynamics in the intermediate Reynolds number range. The layout of this paper is as follows: Section 2 contains the problem formulation and a description of the mathematical tools used. In section 3 we analyse the LES results in terms of Reynolds number and wave age dependence. Both the mean and the fundamental mode momentum balances are investigated to quantify the relevance of the different stresses involved. The LES results point towards the existence of a quasi-laminar regime at intermediate to high wave ages, characterised by a balance between viscous and wave induced stresses. Therefore, in section 4, we tentatively perform a splitting of the flow field into a laminar wave-generated response (Lamb, 1932), and a turbulent shear flow. The former is driven by the non-homogeneous boundary conditions, whereas the latter is driven by a constant volume force, turbulent stresses and the laminar solution. We approach the problem in the opposite direction compared to for instance Miles (1957), where the wave-correlated motion is found by linearisation about an estimated mean flow. The aim of the present formulation is to emphasise the role of the simple wave induced motion for intermediate to high wave ages and explore the possibility of linear solutions in the resulting wave correlated flow fields. With linear solutions, we mean linear in the sense of a flat boundary layer; while turbulent fluctuations are active at all scales, we investigate whether or not the wave-correlated mean flow associated with the turbulent shear flow interacts with itself.

2 Problem formulation

A fully developed pressure-driven open channel flow with a moving lower wavy surface is considered. Figure 1 gives an illustration of the set up. The domain size is $(L_{x},L_{y},L_{z})=(6,1,3)$ wave lengths in the streamwise, vertical and crosstream directions, respectively. The lower surface is described by a propagating wave $y=\eta(x,t)=a\sin(k(x-ct))$ with amplitude $a$, wave number $k=2\pi/\lambda$, and phase velocity $c$. At this boundary, the flow is subject to kinematic boundary conditions from the wave orbital movement imposed by the propagating wave. In the frame of reference that moves with the propagating wave, the wave shape and the wave orbital velocities are stationary, *i.e. *,

[TABLE]

At the top boundary, slip conditions are prescribed, whereas in the streamwise and crosstream directions periodicity is imposed. The fluid motion is described by the incompressible Navier–Stokes equations

[TABLE]

where $x_{i}=(x_{1},x_{2},x_{3})=(x,y,z)$ denotes the streamwise, vertical and crosstream coordinates, respectively. The velocities in the corresponding directions are $u_{i}=(u_{1},u_{2},u_{3})=(u,v,w)$. In the above equations, the summation convention applies to repeated indices. The constant density is absorbed into the pressure, $p$, and the volume force, $\hat{u}_{\tau}^{2}/h$. This volume force is defined as in Kihara et al. (2007), using a virtual friction velocity $\hat{u}_{\tau}$ at the surface. As in that work, we keep $\hat{u}_{\tau}=1$ for all Reynolds numbers and wave ages, and the actual computed friction velocity, $u_{\tau}$, as given later in equation (19), is almost identical to this virtually defined friction velocity for all cases considered. Furthermore, the traceless strain-rate tensor is $s_{ij}=1/2(\partial_{j}u_{i}+\partial_{i}u_{j})$ and we have allowed for a varying kinematic viscosity, $\nu_{t}(x_{i},t)$, since the equations will be solved using LES to account for unresolved turbulence.

2.1 Averaging, RANS and Fourier transform

For analysis, the simulation results are averaged in the crosstream direction and time, thus enabling the Reynolds decomposition

[TABLE]

For simplicity, throughout the rest of the paper we will let $(u,v,p)$ refer to the mean velocity field $(U(x,y)+c,V(x,y),P(x,y))$ and $(u^{\prime}_{i},p^{\prime})$ refer to the turbulent fluctuations. In the frame of reference moving with the wave, the Reynolds Averaged Navier-Stokes formulation of (2) is

[TABLE]

where we have used the same averaging as above, and $\langle u^{\prime}_{i}u^{\prime}_{j}\rangle$ are the Reynolds stresses. In the above equations, the subgrid eddy viscosity terms have not been included, because their contribution to the the momentum balances considered in the rest of the paper is small.

To aid the analysis, the RANS equations are Fourier transformed in the streamwise direction. This may be achieved by first mapping the physical curved domain into a rectangular domain by the transformation

[TABLE]

In this coordinate system the RANS equations (4) can be approximated by

[TABLE]

where

[TABLE]

and

[TABLE]

Here, we have used the approximation $\partial_{y}=(1-g_{\xi})^{-1}\partial_{\xi}\approx\partial_{\xi}$. Strictly, this requires a lower wave steepness than what is considered here, but for analysis purposes it is more convenient to use the simplified expression. In section 4.3, where a numerical solution to the system is sought, the full expression is used.

The variables are expanded in multiples of the fundamental mode,

[TABLE]

where $\hat{f}$ denotes the complex Fourier amplitude function and $M$ is the number of modes. The fundamental mode is $m=\pm 1$, the purely real zero component, or mean, is $m=0$, and higher harmonics are $m=\pm 2,\pm 3,\ldots$.

Nonlinear advection terms and terms involving $g_{\chi}$, henceforth referred to as geometric terms due to their link to the surface shape, both give rise to convolution terms. Note specifically that advection terms involving $\partial_{x}$ transforms to

[TABLE]

Since the function $g(\chi,\xi)$ is periodic with wave number k, the Fourier transform of the last term leads to a coupling with higher and lower harmonics,

[TABLE]

Throughout the rest of the paper the $\xi-$derivative will be denoted $\partial_{\xi}=\mathcal{D}$ and the shorthand notation for the geometric coupling contribution $\mathcal{D}_{g}[\hat{h}_{m}]$ will be used. The continuity equation (4c) for mode $m$ is then simply

[TABLE]

and the $m$-th mode of an advection term can be expressed as

[TABLE]

By application of the above two relations, the momentum equations (4a) and (4b) can be transformed to

[TABLE]

These equations have been casted to conservative form in order to simplify the identification of the various stresses in the flow. Since the wave-correlated stresses are nonlinear terms, they couple higher and lower harmonics. This can be seen from the expansion of the wave correlated shear stress

[TABLE]

The geometrical correction, $\mathcal{D}_{g}$, to the streamwise derivative also introduces a coupling between higher and lower harmonics, even for linear terms such as pressure, as seen from equation (11). In equation (14) and (15), the geometry corrections to the Laplacian have been omitted. We have verified numerically that their contributions to the analysis performed in subsequent sections are negligible.

There is a difference in the way the wave-correlated stresses and the Reynolds stresses are treated. The wave-correlated stresses are inherently mean flow quantities. Therefore, these stresses are obtained as products of the individual Fourier components of the mean flow. On the other hand, the Reynolds stresses are turbulent quantities and their individual modes, $\widehat{\langle u^{\prime}_{i}u^{\prime}_{j}\rangle}_{m}$, can only be accessed from a Fourier transform of $\langle u^{\prime}_{i}u^{\prime}_{j}\rangle$.

3 Large Eddy Simulation

3.1 Numerical method

The Navier–Stokes equations (2) are solved by means of LES, using the unstructured finite-volume node-based solver VIDA from Cascade Technologies (see for instance Ham et al., 2006, 2007). LES subgrid terms are modelled by means of a dynamic Smagorinsky procedure suitable for unstructured grids, as described in Mahesh et al. (2004). The equations are advanced in time using the second order BDF-2 scheme, and a fractional step predictor-corrector procedure is employed to ensure conservation of mass. The computational grid is generated by first considering a rectangular domain with coordinates $(\chi,\xi,\zeta)$. Uniform spacing is imposed in the two horizontal directions ($\chi$ and $\zeta$) and a grid stretching is applied in the vertical direction ($\xi$). For all cases, the average grid spacing in viscous units along the surface is $\Delta\xi^{+}=\Delta\xi_{\mathrm{min}}u_{\tau}/\nu=0.7\nu/u_{\tau}$. A transformation to physical space is obtained by means of

[TABLE]

where $\xi\in[0,h]$, and $\eta(\chi)$ is the surface deformation. This results in an $\mathcal{O}(a/h)$ variation of the grid spacing over a wavelength close to the surface. Furthermore, since there is surface stress variation over a wavelength, the effective viscous grid spacing close to the surface varies from $\Delta y^{+}\approx 0.2$ to $\Delta y^{+}\approx 1.5$. Note that this transformation is different from the one used in to transform the equations for analysis. The grid spacing is comparable to previous DNS studies Sullivan et al. (2000); Kihara et al. (2007); Yang & Shen (2010). However, we use a second order finite volume framework, which for the same grid spacing is not expected to resolve the physics as accurately as the pseudo-spectral methods used in previous DNS studies (see for instance Boyd (2001) for differences in global versus local approximation methods). Therefore, we have used LES to account for unresolved turbulent fluctuations.

The flow over the waves depend on both the wave friction Reynolds number, $Re_{\tau}=u_{\tau}\lambda/\nu$, henceforth referred to as the Reynolds number, and the wave age, $c/u_{\tau}$. The simulations cover a range of wave ages for the three Reynolds numbers $Re_{\tau}=\{200,395,950\}$. In addition, a simulation at $c/u_{\tau}=4$ with $Re_{\tau}=260$ was performed in order to compare with the shear-driven flow of Sullivan et al. (2000). Table 1 shows the details of the different simulations.

3.2 Mean flow description

An illustration of how the flow field varies with Reynolds number and wave age is found in figure 2. It contains pressure contours and mean flow streamlines for combinations of two wave ages ($c/u_{\tau}=4$ and $16$) and two Reynolds numbers ($Re_{\tau}=200$ and $950$). The critical layer, *i.e. *where the streamwise velocity $u$ matches the wave speed $c$, is shown as (   ). Above (below) this line, the mean flow is faster (slower) than the wave. Around the critical layer, the so called Kelvin (1880) cat’s-eyes appear, as seen by the closed streamlines. The cat’s-eye plays an integral role in the theory of Miles (1957), and its role in turbulent flows has been discussed in Sullivan et al. (2000); Kihara et al. (2007).

An increase in wave age shifts the location of the cat’s-eye outwards. On the other hand, an increase in Reynolds number compresses the mean flow profile towards the surface as the boundary layer becomes steeper. In viscous units, however, the critical layer position is almost constant for low wave ages. Specifically, for the wave age $c/u_{\tau}=4$, the critical layer location is $\xi u_{\tau}/\nu\approx 5$ for all Reynolds numbers considered. For higher wave ages, the critical layer location increases with Reynolds number, and for $c/u_{\tau}=16$ its location ranges from $\xi u_{\tau}/\nu\approx 40$ for $Re_{\tau}=200$ to $\xi u_{\tau}/\nu\approx 60$ for $Re_{\tau}=950$. Furthermore, the vertical extent of the cat’s-eye decreases with Reynolds number, thus approaching the inviscid form assumed in Miles (1957).

The form drag is a result of asymmetry in the pressure about the wave crest. Specifically, a positive (negative) form drag is present when the low pressure is behind (in front of) the crest. This is the case for both Reynolds numbers at the low wave age, with a more optimal distribution (for generating form drag) and higher amplitude in the high Reynolds number case. At this wave age, the critical layers are close to the surface and appear highly correlated with the pressure minima. This suggests that the interaction between the background mean flow and the wave motion is dynamically important. For the intermediate wave age cases this correlation is no longer present, and the surface pressure variation appears to be contained in a region well beneath the critical layer. The pressure distribution is also more or less symmetric, resulting in low form drag. Consequently, the variation in form drag with Reynolds number is small at intermediate wave ages.

To see how the wave-correlated motion affects the mean flow balance, we consider the streamwise momentum equation (14) for the zero mode

[TABLE]

where, as before, $\mathcal{D}_{g}=\frac{1}{2}ak\exp(-k\xi)\mathcal{D}$ is the geometric correction associated with the streamwise derivative. The various contributions to the wave-induced stresses can easily be obtained from the decomposition in equation (16). In standard channel flow, there is a balance between the first three terms on the right-hand side of equation (18) (see for instance Pope, 2000). The presence of the wave modifies this balance, resulting in additional wave-induced stresses on the left-hand side of the equation, as well as geometry induced turbulence and pressure terms on the right-hand side. The first term on the left-hand side contains the wave correlated shear stress as defined in Hussain & Reynolds (1970, 1972a, 1972b). However, Hussain & Reynolds (1970) explicitly subtracted the mean before forming their wave correlations. In the present work, the zero mode contribution is retained in the stress. We believe that it should be considered as part of the wave-induced stress, since it is a result of the wave motion and surface geometry.

To illustrate the dependence of the various stresses on Reynolds number and wave age, we plot their respective values, for the same cases as in figure 2, in figure 3. Since the streamwise wave-induced stress in equation (18) is not on a form that allows direct comparison with the shear stresses, we have approximated its influence $ak\exp(-k\xi)\mathcal{D}(\widehat{uu}_{1}^{r}-c\hat{u}_{1}^{r})$ as $\mathcal{D}(ak\exp(-k\xi)(\widehat{uu}_{1}^{r}-c\hat{u}_{1}^{r}))$. Consequently, the streamwise wave-induced stress is presented in premultiplied form. Note that the wave-induced streamwise and shear stresses have opposite signs in the equation. Hence, in the figure, it is the difference of these stresses that yields their dynamical relevance. For all cases, the turbulent shear stress is dominant in the outermost part of the flow, and the momentum balance is similar to that of a channel flow. For the low wave age cases, there is a remarkable dependence in the amplitude of the wave-induced stresses with Reynolds number. For $Re_{\tau}=200$, they have moderate amplitudes and peak in a region where viscous stresses are dominant. On the other hand, for $Re_{\tau}=950$ the wave-induced stresses are large in most of the flow domain. They peak well outside the viscous sublayer, and interact with the turbulent shear stress, as seen from the peaks in both terms at $\xi u_{\tau}/\nu\approx 50$. For the intermediate wave age case, the situation is quite different. The wave-induced stresses have similar amplitudes for the two Reynolds numbers, and their imbalance primarily occurs in the viscous sublayer. There thus appears to be less interaction between turbulent and wave-induced stresses as the wave age increases.

3.3 Form drag

Figure 4(a) shows a comparison of the form drag, $D_{p}/u_{\tau}^{2}$, obtained from the LES at different Reynolds numbers as a function of wave age. In addition, the low Reynolds number DNS results of Sullivan et al. (2000) and the intermediate Reynolds number RANS results of Meirink & Makin (2000) are also plotted. The form drag, $D_{p}$, is defined as the pressure contribution to the surface stress

[TABLE]

For wave ages above $c/u_{\tau}\approx 10$ there is a small upward shift in form drag with increasing Reynolds number for the current simulations. The same behaviour is seen in Meirink & Makin (2000). This implies that the wave age at which the form drag switches from positive (wave growth regime) to negative (wave decay regime) shifts towards higher wave ages. All results point towards a low Reynolds number sensitivity in this regime.

Below $c/u_{\tau}\approx 10$ both the present results and other published work show considerable variation in form drag as a function of Reynolds number, which is consistent with the large Reynolds number sensitivity of the wave-induced stresses discussed in section 3.2. Before we describe the differences, we first note that the Couette flow results of Sullivan et al. (2000) and Yang & Shen (2010) are in good agreement, despite the large differences in friction Reynolds number based on the channel half height, $Re_{\tau}^{h}=u_{\tau}h/\nu$. On the other hand, their wave friction Reynolds numbers are almost identical. This suggests that the wave length, as opposed to the half channel height, is the appropriate outer length scale in this flow. This observation is further supported by comparing our $Re_{\tau}^{h}=Re_{\tau}=200$ results with the results of Kihara et al. (2007), where $Re_{\tau}^{h}=150$ and $Re_{\tau}=200$. Indeed, our results are in excellent agreement with Kihara et al. (2007).

At $Re_{\tau}=200$ the maximum form drag is found at $c/u_{\tau}=8$, but in contrast to the higher Reynolds numbers, there is relatively small wave age sensitivity below $c/u_{\tau}\approx 10$. Although the Reynolds number in the shear-driven flow of Sullivan et al. (2000) is only $30\%$ larger, there is a distinct difference between the two at $c/u_{\tau}=4$. In fact, the results of Sullivan et al. (2000) resemble our $Re_{\tau}=395$ results, although with a slightly lower form drag. At first glance, this seems to support the well known observation that shear-driven flows inherently contain more high Reynolds number dynamics than a pressure-driven flow at the same Reynolds number. To check this, we performed a simulation at $Re_{\tau}=260$ for this wave age. As seen from the figure, our form drag is almost identical to the one in Sullivan et al. (2000). Therefore, it appears that the large variation across data sets at low wave ages is not primarily a flow configuration issue. Instead, the results point towards a high Reynolds number sensitivity at low wave ages.

Both at $Re_{\tau}=395$ and $950$ there is a distinct peak in form drag at $c/u_{\tau}=4$. The $Re_{\tau}=1000$ ( $\vartriangle$ ) results of Meirink & Makin (2000) has a peak at $c/u_{\tau}=8$, but their variation with wave age in this regime is small. They also found that maximum form drag was obtained at $Re_{\tau}=800$, and linked this to an optimal cooperation of viscous and turbulent stresses. For $c/u_{\tau}=5$ they also report a higher form drag at $Re_{\tau}=260$ (figure 9 in that paper) than at $Re_{\tau}=1000$. This is clearly contrasted by our $c/u_{\tau}=4$ results, where $Re_{\tau}=950$ has a substantially higher form drag than $Re_{\tau}=260$. Our results do not contradict the observed $Re_{\tau}=800$ maximum in Meirink & Makin (2000), but the agreement between our $Re_{\tau}=950$ results and their $Re_{\tau}=1000$ results is rather poor. Note also that their $Re_{\tau}=260$ results are only partly in agreement with Sullivan et al. (2000). Although Meirink & Makin (2000), with support from Gent & Taylor (1976), claim that the growth rate is independent of wave steepness below $ak<0.1$, we are not convinced that this is the case.

From equation (19) we see that the form drag is a result of the action of the pressure against the slope of the surface. Since the slope is given as

[TABLE]

we end up with the simple relation,

[TABLE]

From this we observe that the form drag is determined by the out-of-phase pressure at the surface. We follow Mastenbroek et al. (1996); Meirink & Makin (2000), where the out-of-phase surface pressure component is found by integrating the vertical momentum equation (15) from the surface to the freestream. They neglected nonlinear and geometric terms, which could be justified by their low wave steepness. Even for the present steepness of $ak=0.1$, we find that the surface pressure to a large extent may be computed from the wave-induced shear stress component

[TABLE]

with the remaining contribution mainly stemming from the turbulent stress component $\langle v^{\prime}v^{\prime}\rangle$. A comparison of the pressure drag computed using equation (22) and the LES results can be found in figure 4(b).

The vertical dependence of the wave-induced shear stress, $-k(\hat{u}_{0}-c)\hat{v}_{1}^{im}$, for different Reynolds numbers and wave ages, is seen in figure 5. For wave ages below $c/u_{\tau}=12$, all Reynolds numbers have predominantly positive wave-induced shear-stress profiles, resulting in positive form drag. For $c/u_{\tau}=16$, there is a negative region close to the surface which is approximately balanced by a positive outer region, resulting in low form drag. For the highest wave age, $c/u_{\tau}=24$, all Reynolds numbers have strictly negative profiles. As the wave age increases, the support of the wave-induced stresses are increasingly confined close to the surface. The Reynolds number dependence is also simplified in the sense that the shape becomes similar. The amplitude decreases with increasing Reynolds number. On the other hand, for low wave ages the wave-induced stress has support throughout the boundary layer and the dependence of peak locations and amplitudes are non-monotonic in the parameters.

3.4 Reynolds number dependence at a low wave age

To gain insight into the high Reynolds number sensitivity at low wave ages, we consider the wave age $c/u_{\tau}=4$, where simulations have been performed at all four Reynolds numbers. Figure 6(a) and (b) show the mean flow $(\hat{u}-c)/u_{\tau}$ and the out-of-phase vertical velocity $-kh\hat{v}_{1}^{im}/u_{\tau}$, respectively. Close to the surface, the mean flow profiles display the expected behaviour of increasing gradients with increasing Reynolds number. In the outer part of the flow, however, there is no evidence of a direct link between the form drag and the mean flow magnitude. In other words, a lower freestream velocity does not follow from an increased form drag, as would be the case for standard channel flow. This may seem counter intuitive, but as seen from the presence of wave-induced stresses in the zero mode streamwise momentum balance (18), there is no a priori reason for one to follow from the other. Naturally, given the increasing form drag with Reynolds number, the out-of-phase vertical velocity also exhibits a non-monotonic behaviour with increasing Reynolds number. For all Reynolds numbers, the out-of-phase vertical velocity has support up to $\xi k\approx 3$. It therefore has the ability to interact with the mean flow in a large part of the flow domain.

Overall, the Reynolds number dependence can not be attributed to the individual behaviour of the mean flow or the out-of-phase vertical velocity, since both of these vary in a non-monotonic manner. In order to assess which physical processes that are involved in determining the wave-induced stresses, we consider the streamwise momentum balance for the fundamental mode,

[TABLE]

Here, we have used the approximation $(\mathcal{D}^{2}-k^{2})\approx\mathcal{D}^{2}$. Specifically, we consider the imaginary part of this equation, since it contains $\hat{u}_{0}\hat{v}_{1}^{im}$.

Figure 7 shows the imaginary part of (a) the viscous, (b) the turbulent and (c) the wave-induced stress gradients for all four Reynolds numbers. Notice that all stress terms contribute significantly to the momentum balance even in the interior of the domain. The contribution of the pressure is almost constant in most of the domain, and it is thus not shown in the figure.

For all Reynolds numbers, the viscous forces are active from the surface well beyond the critical layer. While the $Re_{\tau}=200$ case exhibits an almost constant behaviour below the critical layer, an increasing variation in the viscous force is seen with increasing Reynolds number. For $Re_{\tau}=950$ there is a negative peak of similar amplitude as the surface value. This indicates that for the fundamental mode, the importance of viscosity around the critical layer increases with Reynolds number.

The turbulent stress gradients display a large variation with Reynolds number. It is interesting that the turbulent and viscous stress gradients have opposite signs and similar amplitudes around the critical layer for all cases. Given the small variation in the wave-induced stresses in this part of the domain, the results hint at a balance between viscous and turbulent forces close to the surface. In the figure, the turbulent stress gradient is the sum of the streamwise and shear-stress contributions. However, at the critical layer, the shear stress is much larger than the streamwise stress for $Re_{\tau}=950$. With decreasing Reynolds number, the streamwise stress becomes more important. This type of behaviour can be explained by comparing the streamlines in figure 2(a) and (b). For $Re_{\tau}=200$ the vertical extent of the cat’s eyes region is large. Therefore, a fluid particle located in its centre will experience a smaller shear anisotropy than in the $Re_{\tau}=950$ case.

In the outer part of the flow, the amplitudes of the turbulent and wave-induced stress gradients increase with Reynolds number in a correlated manner. Notice in particular the extremal points for $Re_{\tau}=950$ around $\xi u_{\tau}/\nu\approx 20$ and $\xi u_{\tau}/\nu=60$. This behaviour supports the rapid-distortion theory of Belcher & Hunt (1993).

3.5 Reynolds number dependence at an intermediate wave age

We next consider the case $c/u_{\tau}=16$. At this wave age, all simulations have close to zero form drag. Previously it was conjectured that the flow exhibits a simpler behaviour in this regime. From figure 5 we also observed that the out-of-phase wave-induced stress term $(\hat{u}-c)\hat{v}_{1}^{im}$ behaved similarly for all Reynolds numbers.

To investigate why this is the case, we consider the stresses in the out-of-phase streamwise momentum equation as a function of Reynolds number. The resulting stress gradients are found in figure 8. At this wave age, all the dynamics in the fundamental mode occur well within the critical layer, and the main contributing terms are the viscous and wave-induced stresses. The turbulent stress gradient becomes more pronounced with increasing Reynolds number, and for $Re_{\tau}=950$ it has a significant contribution above $\xi u_{\tau}/\nu=10$. The dominance of viscous and wave-induced stresses is even more pronounced at $c/u_{\tau}=24$ and $36$. It thus seems that a quasi-laminar regime is entered in the intermediate to high wave age regime.

Note that the wave age can be written as a ratio of Reynolds numbers

[TABLE]

where $Re_{w}$ is the wave Reynolds number (Lamb, 1932). A large wave age implies that the wave Reynolds number is large compared to the friction Reynolds number, or that the smallest scales in the flow are governed by the imposed wave. It thus seems that the simple wave age dependence at higher wave ages is the result of a separation of scales. Turbulence mainly has its dynamical relevance in maintaining the mean flow, while the flow response to the propagating wave is maintained by a near-surface balance between viscous and wave-induced forces. These observations, coupled with the seemingly simple dependence of $D_{p}$ on $Re_{\tau}$ at intermediate to high wave ages, suggest that splitting the flow fields into a wave response and a remainder may be fruitful.

4 A split system approach

Consider a splitting of the flow fields into a “shear” contribution $(u_{i}^{h},p^{h})$ and a wave contribution $(u_{i}^{p},p^{p})$. The latter, hereafter referred to as the particular solution, satisfies the non-homogeneous Dirichlet boundary conditions imposed by the travelling wave, and is the solution to Navier-Stokes equations in the absence of a background shear flow. The former, hereafter referred to as the homogeneous solution, is subject to homogeneous boundary conditions, and is driven by the particular solution and the constant body force, as well as turbulent stresses. The system is hence built in a bottom-up approach (from the known boundary conditions), instead of the common top-down approach (Miles, 1957; Belcher & Hunt, 1993; Kudryavtsev et al., 2001), where linear solutions about a prescribed estimate of the mean flow are sought.

To illustrate the splitting of the flow fields, we consider the decomposed nonlinear advection terms

[TABLE]

The first advection term is solved separately as the laminar flow response to the non-homogeneous boundary conditions, whereas the three latter terms are parts of the homogeneous system that is forced by the particular solution. It is expected that the velocities of the particular solution scale as $ack$, and that the homogeneous solution has a mixed scaling, involving both $ack$ and $u_{\tau}$. We may not a priori state the importance of the different terms, but linearity in the homogeneous system equations would result from the terms $u_{j}^{h}\partial_{j}u_{i}^{h}$ being small compared to the others.

4.1 The particular solution

For the particular solution we use a Helmholtz decomposition

[TABLE]

where $u_{i}^{\phi}=\partial_{i}\phi$ is irrotational and $u_{i}^{\psi}$ is divergence free. The velocity potential can be obtained by solving the Laplace equation $\partial_{ii}\phi=0$ subject to $n_{i}\partial_{i}\phi=n_{i}u_{i}$ at the surface and $\partial_{y}\phi=0$ in the freestream. Here, the surface velocity $u_{i}$ is given by linear Airy wave theory. The equation for the steady streamfunction velocity is

[TABLE]

where $u_{i}^{\psi}=(u^{\psi},v^{\psi})$ is subject to the boundary conditions $u_{i}^{\psi}=u_{i}-u_{i}^{\phi}$ at the lower boundary $y=\eta$, and $u_{i}^{\psi}=0$ in the freestream. When the nonlinear terms are skipped, the above equations represent a linear system forced by a periodic base flow with wave number $k$. The Helmholtz decomposition enables interaction between the velocity potential and the streamfunction velocity through the advection terms. An interesting observation is that the introduction of the geometry transformation (5) leads to the following cancellation

[TABLE]

where we have used $cg_{\chi}=v^{\phi}$. This implies that the single mode behaviour to leading order is well described by a linear system with $u_{i}^{\psi}=\epsilon_{ij3}\partial_{j}\psi=(\mathcal{D}{\psi},-\partial_{\chi}\psi)$,

[TABLE]

subject to $\mathcal{D}{\psi}=-2ack\sin(k\chi)$ and $\partial_{\chi}\psi=0$ at the lower surface $\xi=0$, and $\mathcal{D}{\psi}=\partial_{\chi}\psi=0$ at the top boundary $\xi=1$. This system is the one-phase version of the system in Lamb (1932); Harrison (1908), where we match only the kinematic properties of the lower phase. Due to the viscous scaling of the streamfunction velocities, the standard linear wave theory approximation of evaluating at $y=0$ is prone to increasing errors with increasing wave Reynolds number. Therefore, the boundary conditions must be evaluated at the interface ($\xi=0$).

The main motivation for resorting to a coupled formulation is that we are concerned with the multi-mode behaviour of the particular solution. Specifically, we are dealing with the transfer of momentum from the fundamental mode to lower and higher harmonics, which is enabled by the action of the velocity potential on the viscous solution. The particular solution will in turn act as a multi-modal base flow for the homogeneous solution, and the effect of small mean or higher harmonic contributions can not be neglected a priori.

4.2 The homogeneous solution

Once the particular solution is obtained, one may solve for the homogeneous solution, which has the following RANS formulation

[TABLE]

The velocity field is subject to a no-slip condition at the surface, and a slip condition in the freestream. Once nonlinear terms ($u_{j}^{h}\partial_{j}u_{i}^{h}$) are skipped, these equations represent a linear system forced by a multi-modal base flow. The body force term, $\hat{u}_{\tau}^{2}/h$, is chosen in accordance with the LES set up, and the turbulent Reynolds stresses are taken from the LES simulations.

The choice of using turbulent stresses from the LES renders this framework of little practical use as a RANS model, and indeed this is not the aim of the present work. Instead, we aim at using the split system approach to analyse the LES results and better understand the underlying dynamics of the flow. It should be pointed out that we also could have obtained the homogeneous solution by subtracting the particular solution from the LES results directly, but this approach is not as convenient for exploring the dependencies of the equation system on various assumptions.

There has been some debate regarding how to properly model turbulent stresses in flow over waves. Belcher & Hunt (1993) argued, on the basis of a rapid-distortion framework, that eddy viscosity type closures are suitable only in the near surface region, since turbulent eddies far from the surface are transported too quickly to be strained by the local shear. This has been confirmed numerically by Mastenbroek et al. (1996) who compared different RANS closures to experimental results. In the present work, we have tested the simple eddy-viscosity closure of Reynolds & Tiederman (1967), based on the Cess (1958) framework, which has proven successful in studying linear energy amplification in turbulent channel flow Pujals et al. (2009). Indeed, we found that no simple relation between the mean shear and the turbulent stresses is present.

4.3 Solution procedure

The equation systems (27) and (30) have a similar form that can be written

[TABLE]

where $\mathcal{L}=\nu\nabla^{2}+(c-U)\partial_{x}-V\partial_{y}$, with $\nabla^{2}=\partial_{xx}+\partial_{yy}$ being the Laplacian operator. The above system represents the equations for the particular solution when the base flow is $(U,V)=(u^{\phi},v^{\phi})$ and the forcing is zero ($f=0$). It represents the homogeneous solution when $(U,V)=(u^{p},v^{p})$ and the forcing terms are

[TABLE]

Here, the additional forcing term $f_{b}$ is added to correct for incompatibility of the divergence free criterion used in the LES and the spectral collocation method, which leads to a momentum imbalance of the LES results on the spectral collocation grid (see Ham & Iaccarino, 2004, for a discussion on the conservation properties of the LES). This term is generally small, but the streamwise mean velocity in the outer part of the domain is highly sensitive to any pointwise momentum imbalance.

In the linear case, it is possible to solve system (31) explicitly by introducing a coordinate transformation to map the physical domain to a rectangular domain and then Fourier transforming along the streamwise direction, as was done in previous sections. However, both the periodic base flow and the geometrical terms lead to coupling of the different harmonics, which makes this a tedious and time consuming effort. Instead, we have used a more attractive approach based on two-dimensional Fourier-Chebyshev collocation Weideman & Reddy (2000); Hoepffner et al. (2019). The physical domain $0\leq x\leq 2\pi/k$ and $\eta(x)\leq y\leq 1$ is mapped to a rectangular computational domain $0\leq\chi\leq 2\pi/k$ and $0\leq\xi\leq 1$ using the transformation (17), and the corresponding derivative operators are redefined numerically to account for this mapping. Since the base flow is periodic, the resolution in the streamwise direction can be kept rather low. For the particular solution, we found $(N_{x},N_{y})=(8,100)$ to be sufficient at all Reynolds numbers. The homogeneous solution has an additional forcing from turbulent stresses, which results in stricter resolution requirements. We found $(N_{x},N_{y})=(12,120)$ to be sufficient for the cases considered.

In the nonlinear configuration, the system is solved using simple Newton iterations. The linear solution is used as an initial guess. For all but the highest Reynolds number, we found that convergence could be reached in less than ten iterations. At the highest Reynolds number, convergence is sensitive to both initial guess and the choice of iterative solver. We found that convergence could be ensured by successively adding Fourier components of the turbulent forcing.

4.4 Discussion on the particular solution

Figure 9 shows vertical profiles of the fundamental mode velocities of the particular solution for the four wave Reynolds numbers $Re_{w}=\{10,100,1000,10000\}$. To put these wave Reynolds number in context, one may use relation (24). For instance, at $c/u_{\tau}=2$ the corresponding wave Reynolds numbers of $Re_{\tau}=200$ and $950$ are $Re_{w}\approx 60$ and $300$, respectively. The in-phase parts of the streamwise and vertical velocities are shown in figure 9(a), and the out-of-phase velocities are shown in figure 9(b). The vertical velocity approaches the velocity potential with increasing wave Reynolds number, which means that the out-of-phase part tends to zero, and the in-phase part approaches the velocity potential $-akc/2\exp(-k\xi)$. In the outer part of the domain, the streamwise velocity is well described by the velocity potential, but the streamfunction velocity is needed close to the surface to match the boundary conditions.

As previously stated, the fundamental mode of the particular solution may be approximated from the analytical solution of equation (29). For high wave Reynolds numbers ($Re_{w}>100$) we may write the solution as

[TABLE]

The velocity components are then approximately

[TABLE]

Here, the streamwise velocity amplitude is independent of Reynolds number, while the vertical velocity amplitude decays with Reynolds number as $\sqrt{2/Re_{w}}$. We also observe that vertical profiles collapse when scaled with $k\sqrt{Re_{w}}$. The same behaviour is seen for the full particular solution, with minor corrections to the amplitude of the streamwise velocity. The out-of-phase streamwise velocity reflects the adjustment from the imposed boundary condition to the velocity potential in the outer flow. The phase shift occurs in the region $k\xi\sqrt{Re_{w}}=[0,10]$. The only clear difference for the fundamental mode behaviour between the current formulation and the uncoupled analytical approach given in (29), which is similar the the one in Lamb (1932); Harrison (1908), is that the peak amplitude is slightly higher for the coupled system at high Reynolds numbers.

We find that the coupling between the velocity potential and the streamfunction velocity first and foremost gives rise to a positive valued zero mode with a peak amplitude of approximately half of the fundamental mode, as seen in figure 9(c). The amplitude of the second harmonic is approximately $ak$ times that of the fundamental mode for all wave Reynolds numbers considered. This is indicative of linear behaviour. The present case, with a wave steepness of $ak=0.1$, is transitional; for less steep waves ($ak<0.05$), coupling terms are negligible, and for steeper waves ($ak>0.15$) both coupling terms and nonlinear interactions of the streamfunction velocity with itself become important.

4.5 Discussion on the total solution

Figure 10 shows a comparison of LES results with both linear and nonlinear split-system solutions for $Re_{\tau}=395$ at the two wave ages $c/u_{\tau}=4$ and $24$. The split-system solution is the sum of the homogeneous and particular solution. For the low wave age case, seen in the top frames of the figure, neither the streamwise mean velocity nor the out-of-phase fundamental mode velocities are captured by the linear system. On the other hand, the nonlinear solution is in good agreement with the LES results. In the outermost part of the domain, the streamwise mean velocity deviates significantly from the LES. This discrepancy, which is present for all cases considered, is most likely due to slight differences in implementation of the top boundary condition. If we explicitly enforce the freestream velocity from the LES as the top boundary condition, the profiles are in close agreement. Despite the deviation of the streamwise mean profiles in the outer part of the domain, the agreement in form drag between the two is excellent for all wave ages and Reynolds numbers. The reason for this, is that the deviation occurs in a region where the out-of-phase vertical velocity has negligible amplitude.

For the high wave age case, the streamwise mean velocities of the linear and nonlinear solutions are almost identical. For the out-of-phase velocity components, the linear solution captures the main features of the flow, but the amplitudes are somewhat underpredicted. For the even higher wave age of $c/u_{\tau}=36$, the linear and nonlinear solutions are in excellent agreement. In the present formulation, the linearity requirement is very strict, since it implies that the nonlinear term $v_{h}\mathcal{D}u_{h}$, which includes the mean shear, is negligible in the momentum balance. The wave age at which the linear regime is entered increases with Reynolds number.

To explore the functional dependence of the form drag with Reynolds number and wave age, we consider the particular and homogeneous surface pressure separately in figure 11(a). The surface pressure of the particular solution is strictly negative, and for low wave ages its contribution is almost negligible compared to its homogeneous counterpart. With increasing wave age this situation changes, and the particular solution eventually outweighs the positive contribution of the homogeneous solution, resulting in a net negative form drag. The behaviour of the homogeneous pressure is non-monotonic, with a minimum around $c/u_{\tau}=16$. Beyond this, we observe an approximate linear growth with wave age.

As previously stated, the main contributor to the surface pressure is the wave-induced shear stress term (see equation (22)), which decomposes to

[TABLE]

The first term on the right hand side is the contribution to the particular solution pressure. Since the mean flow of the particular solution is at least $ak/2$ smaller than the phase velocity (see section 4.4), $\hat{u}_{0}^{p}\hat{v}_{1}^{p,\mathrm{im}}$ is small and can be neglected. With this simplification, a scaling for the particular solution pressure can be obtained using the analytical solution (34). The analytical solution gives the following estimate for the magnitude of the out-of-phase vertical velocity,

[TABLE]

This scaling is consistent with the profiles in figure 12(b). It implies that the particular solution pressure scales as

[TABLE]

The particular solution pressure decreases in magnitude with increasing friction Reynolds number and grows rapidly with increasing wave age. This functional form is clearly seen from the blue lines in figure 11(a).

The rest of the terms on the right hand side of (35) contribute to the homogeneous solution pressure. The term containing the zero mode of the particular solution is small and can be neglected, whereas the three remaining terms are shown for $Re_{\tau}=395$ in figure 11(b). Of the velocity components involved in these expressions, only the out-of-phase vertical velocity of the homogeneous solution displays a complicated wave age dependence (see fig. 12(c)). For low wave ages, $\hat{v}_{1}^{h,\mathrm{im}}$ is negative and large, with peak amplitude encountered at $c/u_{\tau}=4$. Consequently, the surface pressure is dominated by the term $(c-\hat{u}_{0}^{h})\hat{v}_{1}^{h,\mathrm{im}}$. At approximately $c/u_{\tau}=12$ there is a switch in behaviour of the out-of-phase vertical velocity, and a change of sign occurs close to the surface. As the wave age increases, the profiles become strictly positive, until finally, the change with wave age is slow (as seen by the overlapping lines of $c/u_{\tau}=24$ and $c/u_{\tau}=36$). This dictates the high wave age behaviour seen in figure 11(b). Firstly, the integral of $c\hat{v}_{1}^{h,\mathrm{im}}$ ends up being positive and scales linearly with wave age. Secondly, the integral of $-\hat{u}_{0}^{h}\hat{v}_{1}^{h,\mathrm{im}}$ reaches a constant state, since $\hat{u}_{0}^{h}$ (see fig. 12a) is unaltered in the region where $\hat{v}_{1}^{h,\mathrm{im}}$ is non-negligible. Lastly, the contribution from, $\hat{u}_{0}^{h}\hat{v}_{1}^{p,\mathrm{im}}$ is positive and scales approximately as $(c/u_{\tau})^{1/2}$. The end result is an almost linear dependence of the homogeneous surface pressure for high wave ages. Based on the above scaling, we can estimate the pressure drag in the quasi-laminar regime given its inception value. The existence of the quasi-laminar regime is good news, as it reduces the parameter space of interest for future investigations to a smaller region of the wave age Reynolds number plane. As previously noted, the onset of this quasi-laminar regime is Reynolds number dependent. Exploring this Reynolds number dependence is a topic for future work.

The sensitivity of the form drag to variation in Reynolds number at low wave ages is more challenging to evaluate using the present decomposition. From the results in figure 10, we know that the full nonlinear system is required to describe the flow dynamics. The analysis in section 3.4 established that the interaction of the wave-induced stress with turbulence increases with Reynolds number. To quantify the importance of the wave correlated turbulent stresses, we solve the nonlinear split-system using only the zero mode of the Reynolds stresses as a forcing. The resulting wave-induced stresses, for $Re_{\tau}=200$ and $950$ at the wave age $c/u_{\tau}=4$ are found in figure 13. For $Re_{\tau}=200$, using only the zero mode yields a fairly good estimation of the wave induced stress and the form drag is slightly overpredicted. On the other hand, for $Re_{\tau}=950$ the stress is severely underpredicted, which results in a fifty percent underestimation of the form drag. The fundamental mode turbulence is thus an essential part of the dynamics responsible for the form drag.

Another interesting observation is that at the lowest wave age, $c/u_{\tau}=2$, and $Re_{\tau}=200$, using only the forcing of the zero component yields almost identical results to the full system. In addition, the coupling with the particular solution can be excluded altogether, which implies that the kinematics of the underlying wave is not important. The wave kinematics and turbulence is therefore fully decoupled. However, as the Reynolds number is increased, the exclusion of coupling with the particular solution leads to an increasing inability to describe the flow. The main difference in the particular solution with increasing Reynolds number is that the shear becomes stronger (due to the scaling $\sqrt{Re_{w}}$). This means that at low wave ages, the effect of both the laminar forcing from the wave kinematics and turbulence increase with Reynolds number. Since the laminar forcing also contributes to generating turbulence in multiple harmonics, isolating the effect of one from the other is a challenge.

5 Conclusions

In this paper we have presented results from wall resolved large eddy simulations of turbulent flow over a simple propagating wave at a moderate wave steepness of $ak=0.1$. We have varied both Reynolds number and wave age over a range previously not considered in such detail.

A key quantity of interest is the dependence of form drag on the governing parameters. In the intermediate to high wave age regime, previous studies have found a somewhat simplified dependence of form drag on both Reynolds number and wave age, but to the authors’ knowledge, conclusions have been limited to an observed collapse when scaling the problem using an outer velocity scale instead of the friction velocity. In the low wave age regime, the form drag depends in a complex manner on Reynolds number, wave age and wave steepness.

We analysed the results by Fourier-decomposing the computed flow variables, and evaluate the contributions to the momentum balance in both the mean and fundamental mode. Using this decomposition, the form drag can be found by integrating the vertical momentum equation from the surface to the free stream. The dominant contribution to this integral was found to be the wave-induced shear stress which is the product of the streamwise mean velocity and the out-of-phase vertical velocity. We observed that for low wave ages, this wave-induced stress has a significant contribution in a large part of the domain. On the other hand, at high wave ages the stress is increasingly confined to a region close to the surface.

By examining the results in terms of Reynolds number dependence at a low wave age, we found that the interaction between the wave induced and turbulent stresses plays an increasingly important role in the streamwise momentum balance of the fundamental mode as the Reynolds number increased. At an intermediate wave age, we found that the fundamental mode balance was primarily dominated by viscous and wave induced stresses. This indicates that a quasi-laminar regime, where wave kinematics become increasingly important, is entered as the wave age increases. Therefore, we introduce a novel split system approach, where we first construct a laminar flow response to the wave kinematics. This in turn acts as a forcing on a shear flow subjected to homogeneous boundary conditions. To account for the effects of turbulence, we force the system using Reynolds stresses from the large eddy simulation. By splitting the flow in this manner, we were able to formulate an analytic functional dependence for the form drag associated with the laminar response. The results also show that the form drag of the shear flow simplifies at high wave ages, and an approximate wave age dependence of the full solution could be constructed. The reason for this simplicity is that, at high wave ages, the forcing from the laminar wave-induced flow overwhelms the geometry induced interaction of the shear flow with itself.

The high sensitivity of the form drag to variation in Reynolds number at low wave ages is more challenging to evaluate using the split system approach. This is primarily due to the increased importance of nonlinearity in the shear flow. These nonlinearities are inherently coupled with higher harmonics in the turbulent stresses. Nevertheless, the split system approach can be utilized to quantify the importance of different harmonics in the turbulent stresses by explicitly choosing which modes to include in the split system forcing. In support of Belcher & Hunt (1993), we demonstrated that the fundamental mode of the Reynolds stresses becomes increasingly important for accurate prediction of the form drag as the Reynolds number increases. However, the importance of nonlinearity in the low wave age regime suggest that any linear low wave age model for the form drag can only work for a very low wave steepness and/or Reynolds number. The dependence of form drag on wave steepness and Reynolds number in the low wave age regime is an interesting topic for future work.