Large–Eddy Simulation of the Flow Inside a Scour Hole Around a Circular Cylinder Using a Cut Cell Immersed Boundary Method

Abstract

We present a novel symmetry-preserving cut cell finite volume method which is a three-dimensional generalisation of the method by Dröge and Verstappen (Int J Numer Method Fluids 47:979–985, 2005). A colour-coding scheme for the three-dimensional cut momentum cell faces reduces the number of possible cut cell configurations. A cell merging strategy is employed to alleviate time step constraints. We demonstrate the energy conservation property of the convective and pressure gradient terms, and the second-order spatial convergence with suitable benchmark cases. We used the scheme to perform highly resolved large–eddy simulations of the flow inside a scour hole around a circular cylinder mounted vertically in a flume. The simulation results are extensively compared to a stereoscopic particle image velocimetry experiment of the same configuration performed by Jenssen and Manhart (Exp Fluids 61:217, 2020). We demonstrate that for the investigated Reynolds numbers (20,000 and 40,000) nearly converged solutions are obtained; however at large computational efforts (up to 2.35 billion cells for the higher Reynolds number). It turns out that the flow topology of the horseshoe vortex system is strongly dependent on the grid resolution. For simulation results obtained on the finest grid, the mean flow and turbulence quantities agree well with the experiment. We investigate the shape and turbulence structure of the horseshoe vortex based on three-dimensional fields, and discuss the distribution of the mean and standard deviation of the wall shear stress in the scour hole and the implications for the physics of the scouring process over a sand bed.

1 Introduction

Local disturbances of the flow around hydraulic structures such as bridge piers can result in increased wall shear stress levels and the development of local scour holes. Such scour holes can be a severe threat to the safety of the hydraulic structure. The estimation of the depth of a scour hole is a difficult task and empirical formulae for the design of pier foundations go in hand with large uncertainties (Link 2006). Numerical models coupling a flow simulation with a sediment transport model promise a better accuracy for predicting the evolution of scour holes (Roulund et al. 2005). However, there is a long error chain starting with a prediction of the flow field and ending with a sediment transport model. In every part of the modelling chain a number of uncertainties arises. It is the purpose of this research to shed light on the uncertainties in the flow simulation and to provide highly resolved and accurate data of the turbulent flow field inside the scour hole and thus alleviate one of the numerous uncertainties in the modelling chain. A particular focus will be placed on the reliability of the wall shear stress as it can be used as the basis of a sediment transport model (Roulund et al. 2005; Kraft et al. 2011).

The so-called horseshoe vortex was found to be the main agent in the development of the scour hole around a circular cylinder (Melville and Raudkivi 1977). It has been suggested to give rise to large local wall shear stresses on the mobile bed. A large number of investigations has been devoted to the dynamics of the horseshoe vortex (Dargahi 1990; Dey and Raikar 2007; Oliveto and Hager 2002; Escauriaza 2008; Kirkil et al. 2005, 2008, 2009; Link et al. 2012; Khosronejad et al. 2012; Baykal et al. 2015). Unfortunately, very few experimental data have been published which spatially resolved the turbulence structure and dynamics of this horseshoe vortex. Unger & Hager 2006 applied monoscopic particle image velocimetry (PIV) to measure the flow inside a half model of a scour hole. They used a half model to gain optical access to the flow whereas in a full model, the scour hole is hidden to an observer from the side. This optical access-problem has been solved by Jenssen and Manhart (2020) who used a stereoscopic PIV. The stereoscopic setup has the additional advantage that it provides three-dimensional velocity data in a two-dimensional plane. Fully three-dimensional data have been reported so far only from numerical simulations.

Previous simulations of the flow inside a scour hole have employed Reynolds-averaged simulations (Roulund et al. 2005), detached-eddy simulations (DES) and large–eddy simulations (LES) (Kirkil et al. 2008) with block-structured or overset boundary-fitted curvilinear grids (Escauriaza 2008; Kirkil et al. 2009) or on boundary-fitted unstructured meshes (Baykal et al. 2015; Song et al. 2022). In all these simulations, however, no direct comparison to experimental data – time-averaged flow fields, turbulence statistics and wall shear stress distribution – has been presented. This is partly due to the unavailability of reliable experimental data at the time the simulation was performed or due to incompatibility of the scour geometry used. In this paper, we use the recently published measurements of Jenssen and Manhart (2020) to validate a highly-resolved LES of the same flow configuration.

We performed LES using a three-dimensional cut cell finite volume method on a block-structured Cartesian grid. The Cartesian grid method coupled with a variant of the immersed boundary method has certain advantages and certain disadvantages compared to methods which employ body-fitted grids. The main advantage of Cartesian methods is the numerical efficiency which allows for a considerably larger number of grid cells than unstructured mesh methods (Strandenes et al. 2015), leading to a higher numerical accuracy due to the absence of grid distortion effects. The Achilles heel of Cartesian immersed boundary methods is the near-wall resolution. In general, a proper wall-normal grid resolution requires a refinement in all three coordinate directions which is much more expensive in terms of grid cell numbers than an anisotropic wall-normal refinement in a body-fitted grid approach. On the other hand, an efficient refinement strategy in a Cartesian method using zonal grids (Manhart 2004) helps to keep the numerical effort moderate and simultaneously provides grid resolution in off-wall critical flow zones. The choice of a cut cell approach allows to determine the wall shear stress consistently with the simulation. In contrast, in a ghost cell method the wall shear stress can only be obtained in a post-processing step, introducing additional uncertainty.

In the first part of this work (Sect. 2), we describe a novel cut cell method that generalises the method of Dröge and Verstappen (2005) to three dimensions. The large variety of cut momentum cells is categorised using a colour-coding scheme on which the selection of the appropriate discretisations of the fluxes is based. A cell merging scheme is used to alleviate the time step constraints that would be introduced by small cut cells. The energy conservation property of the discrete convective and pressure gradient terms, which is a distinguishing feature of our method, is demonstrated for an inviscid flow. Moreover, second order convergence to an analytical solution is established and the results for a benchmark case (flow around a sphere at $Re =100$) are compared to those reported in the literature.

In the second part of this work, we intend to clarify the question which computational effort is required to achieve reliable and accurate results from a large–eddy simulation of the flow inside a scour hole around a circular cylinder using our cut cell method. In particular, we aim to address the following research questions. What grid resolution is necessary to obtain (i) The correct flow topology and (ii) Distributions of flow quantities that agree with experimental results? How reliable are LES results at a given grid resolution? What is the three-dimensional shape of the turbulent horseshoe vortex in the time-averaged sense? What is the structure, shape and coherence of the horseshoe vortex in the instantaneous representations? How are mean values and fluctuations of the wall shear stress distributed and linked to the horseshoe vortex?

The paper is structured as follows. In the next section, the numerical method and its verification are documented. The flow configuration and the simulation setup are presented in Sect. 3. The results Sect. 4 Contains a detailed assessment of grid effects on the computed flow topology 4.2 and velocity profiles 4.2). Some flow and turbulence quantities are compared with measured results in Sects. 4.3 and 4.4 presents three-dimensional representations of the time-averaged and instantaneous horseshoe vortex configuration. The accuracy and distribution of the wall shear stress is discussed in Sect. 4.5. Finally in Sect. 5 we conclude and discuss our results. In the appendix some details of the cut cell method are documented.

2 The Numerical Method

2.1 Basic Numerical Method

We solved the conservation laws for incompressible Newtonian flow by an in-house Finite Volume solver (MGLET). This solver uses a Cartesian grid with a staggered arrangement of the variables. The fluxes are approximated by second order central schemes which lead to a skew symmetric convective term that is energy conserving (Verstappen and Veldman 2003). For time integration, a third order Runge-Kutta scheme is embedded in a pressure–velocity coupling according to Chorin’s projection method (Chorin 1968). The resulting Poisson equation for the pressure update is solved by the strongly-implicit procedure (SIP) of Stone (1968). For local mesh adaption, we use the method of zonal grid refinement introduced by Manhart (2004). The code is parallelised using the Message Passing Interface (MPI).

The sub-grid scale stresses (SGS stresses) are parameterised by the WALE model (Nicoud and Ducros 1999) and the wall shear stress is computed by the Werner-Wengle wall function (Werner et al. 1993) which switches between the linear law of the wall and a $\frac{1}{7}$-th power law depending on the wall distance.

The non-Cartesian boundaries are represented by a cut cell immersed boundary method that will be described in the next section. A verification simulation of the implementation will be presented in Sect. 2.3.2.

2.2 The Cut Cell Method

For the simulations presented here, we developed and implemented a variant of the cut cell method originally proposed by Dröge and Verstappen (2005), Dröge (2006) and extended upon by Cheny and Botella (2010), Nikfarjam et al. (2018), Portelenelle et al. (2019). The advantage of the cut cell method compared to the ghost cell method (Peller et al. 2006; Peller 2010) that we have used for previous simulations (Schanderl and Manhart 2016; Schanderl et al. 2017) is that the whole fluid volume is covered by control volumes. This gives improved conservation of mass and momentum as no leakage flux can occur near the immersed surface. The wall shear stress – which is one of the quantities in focus – is an intrinsic flux quantity used in the momentum balance of the cut cells whereas in a ghost cell method, the wall shear stress needs to be evaluated in a post-processing step.

The underlying design principle of our implementation is the cut cell method proposed by Dröge and Verstappen (2005) which derives the approximations based on the skew-symmetry of the convective and the symmetry of the viscous operator. This approach has the advantage of strict energy conservation on the cut cells. However, the original proposal of Dröge and Verstappen (2005), Dröge (2006) was only for two-dimensional geometries. While there have been proposals for three-dimensional cut cell methods for compressible flows on collocated grids (Hartmann et al. 2009, 2011), for incompressible flow by Meyer et al. (2010) or by Kirkpatrick et al. (2003) or Xie et al. (2020), the symmetry-preserving method of Dröge and Verstappen (2005) has been considered as too complicated for general 3D cases. For example, Nikfarjam et al. (2018) presented a partially three-dimensional extension of the two-dimensional cut cell method of Cheny and Botella (2010), but state that “[o]ne of the major difficulty in cut cell discretization for 3D geometries is the large number of cases to consider: we have enumerated 16 different types of cut cells in 3D whereas there are only 3 in 2D.”, and also Beltman et al. (2020) considers the “[...] extension of the method [of Dröge and Verstappen (2005); Cheny and Botella (2010)] to 3D [as] problematic due to the many possible cut cell configurations.” We address these problems by introducing a colour-coding scheme to classify the complicated cell faces of the three-dimensional momentum cells. Based on the assigned face colour, the appropriate discretisation of the momentum fluxes are selected. The coding scheme essentially reduces the 3-dimensional cut cell configurations to two-dimensional ones. Another novelty of our implementation is the three-dimensional cell merging scheme to alleviate the time step constraints in explicit time integration schemes. The method was first presented at a conference (Kurz et al. 2014) and further extended and validated for aero-acoustical computations by Kreuzinger et al. (2015), Ewert and Kreuzinger (2021); however, it has not been described in detail yet. In what follows, we formulate our scheme for uniform and cubic grids which we have used in the simulations presented here.

2.2.1 Pressure Cells and Continuity Equation

The computational grid is composed of pressure cells which are labelled with three Cartesian indices (i, j, k). The pressure $p_{i,j,k}$ is located at the cell centre and the velocities $u_{i,j,k}$, $v_{i,j,k}$ and $w_{i,j,k}$ are staggered in the positive x-, y- and z-coordinate direction, respectively.

The body geometry is intersected with the edges of the pressure cells. Assuming linear cuts, the open face areas $S^{u}_{i-1,j,k}$, $S^{u}_{i,j,k}$, $S^{v}_{i,j-1,k}$, $S^{v}_{i,j,k}$, $S^{w}_{i,j,k-1}$ and $S^{w}_{i,j,k}$, the corresponding face centre points and the open cell volume $V^{p}_{i,j,k}$ are obtained. Depending on the curved polygonal wall face, 128 different forms of three-dimensional pressure cells are possible if rotated and mirrored cases are taken into account.

For a cut pressure cell, the velocity is defined at the midpoint of the open face area of the pressure cell. The discrete continuity equation is formulated according to Cheny and Botella (2010), in which the midpoint rule is used to approximate the volume flux across the face areas of the pressure cell by multiplying the midpoint velocity of the neighbouring momentum cell with the cell-face area. For example, $\overline{u}_{i,j,k}$ is the volume flux between the cells (i, j, k) and $(i+1,j,k)$ and is computed as $\overline{u}_{i,j,k}=S^{u}_{i,j,k}u_{i,j,k}$.

Fig. 1

Geometric definition of cut momentum cells. The blue lines delimit the control volume boundaries of the staggered u-component. Left:non-balanced momentum cells (red); Right: cell-to-boundary formulation

2.2.2 Momentum Cells and Momentum Balance

As described by Dröge (2006) there are two ways to formulate momentum cells which are cut by an immersed boundary. They can be distinguished by the way how small cut cells are treated. In the first option, small cut cells are dropped, which leads to cells which are not part of the overall momentum or mass balance. We chose the cell-to-boundary formulation in which small momentum cells are merged with neighbouring cells, thus extending the momentum balance to the wall, see Fig. 1. This comes with the drawback of introducing additional errors due to larger aberrations between volume and face midpoints. However, the advantage of this formulation is that the boundary created from staggered momentum cells conforms to the three-dimensional geometry. This simplifies the calculation of boundary forces and facilitates the application of wall models.

The momentum cell is determined according to Dröge and Verstappen (2005) by halving the two adjacent pressure cells in the direction of the velocity component to which the cell belongs. When the neighbouring momentum cell is blocked, the cell is extended to the boundary. The volume of a momentum cell can then be easily computed as

$$\begin{aligned} V^u_{i,j,k}= & {} \frac{1}{2}\left[\underbrace{V^p_{i,j,k}+V^p_{i+1,j,k}}_{\text {default formulation by Droge and Verstappen (2005)}}\right.\\ &\left.+\underbrace{(1-B^{u}_{i-1,j,k})V^p_{i,j,k}+ (1-B^{u}_{i+1,j,k}) V^p_{i+1,j,k}}_{\text {cell-to-boundary}}\right] \end{aligned}$$

(1)

where the indicator field $B^{u}_{i,j,k}$ is 1 if the face belonging to the velocity $u_{i,j,k}$ is open and 0 if it is closed. Using this definition of the momentum cell ensures that the balanced volume is the same for the mass and all momentum components.

The semi-discretised Navier–Stokes equations read

$$\begin{aligned} \mathbf {M}\mathbf {u}_h&= 0 \end{aligned}$$

(2a)

$$\begin{aligned} \mathbf {V}\frac{\mathrm {d}\mathbf {u}_h}{\mathrm {d}t} + \mathbf {C}(\mathbf {u}_h)\mathbf {u}_h&= -\mathbf {G}\mathbf {p}_h - \mathbf {D}\mathbf {u}_h \end{aligned}$$

(2b)

where $\mathbf {M}$ is the discrete divergence operator, $\mathbf {V}$ is a diagonal matrix containing momentum cell volumes according to Eq. (1), $\mathbf {C}$ is the convective operator, $\mathbf {G}$ is the discrete gradient operator and $\mathbf {D}$ is the discrete diffusive operator. The symmetry-preserving method of Verstappen and Veldman (2003) requires that $\mathbf {G} = \mathbf {M}^\mathrm {T}$, $\mathbf {C}(\mathbf {u}_h) = - \mathbf {C}(\mathbf {u}_h)^\mathrm {T}$ and $\mathbf {D}+\mathbf {D}^\mathrm {T}$ is positive definite.

The convective and diffusive approximations are selected with the help of a colour-coding scheme for the face areas of the momentum cell. We define 7 different cell-face categories based on their orientation and the way they are cut by the boundary. The use of the colour-coding effectively reduces the discussion of the cell shapes from three to two dimensions. The formulation of the fluxes is described in the appendix. The convective approximations are implemented according to the method of Dröge and Verstappen (2005) and the viscous terms are implemented using the non-orthogonality correction of Kirkpatrick et al. (2003). This correction is needed because for cut cells the straightforward two-point finite difference between the velocities located on the faces of the pressure cell is not a second order accurate representation of the gradient projected on the normal vector at the centre of the momentum cell face.

In the following, we describe the colour-coding of the cell faces. The faces of a momentum cell in negative and positive x-direction are abbreviated with w (west) and e (east); the faces in negative and positive y-direction with s (south) and n (north) and the faces in negative and positive z-direction with b (bottom) and t (top). Every cell-face is subdivided at the boundary of the pressure cell and, in the cell-to-boundary case, at the midpoint of the pressure cell. Each segment is then marked with one or more colours. All possible colour-coded face categories are shown by using the example of different u-momentum cells in Fig. 2. The coding scheme uses a total of 7 colours (orange, yellow, red, pink, black, grey and green). The different cases can be discriminated based on the logical information whether the neighbouring u-velocity cells are open or closed. This information is represented by the $B^u$ indicator field.

The colour codes for an open u-momentum cell ($B^u_{i,j,k} = 1$) are assigned based on the following rules:

• (1) Each wall face area is assigned the colour red.

• (2) If the neighbouring momentum cell at the e-face (w-face) is open, i.e. $B^u_{i+1,j,k} = 1$, the part of the e-face area which is parallel to the grid is represented by the colour orange. For the skewed parts of the corresponding face the colour yellow is used.

• (3) If the neighbouring momentum cell at the t-face (n-, s- or b-face) is open, i.e. $B^u_{i,j,k+1} = 1$, the face part parallel to the grid is marked with the colour pink. Moreover, if the neighbouring momentum cell at the w-face is closed ($B^u_{i-1,j,k}=0$), the t-face area of a momentum cell is expanded towards the wall on the w side. For this segment the additional colour grey is used if the t-w neighbouring momentum cell is open ($B^u_{i-1,j,k+1}=1$) and the additional colour black if the t-w neighbouring momentum cell is closed ($B^u_{i-1,j,k+1}=0$).

• (4) Else the neighbouring momentum t-face (n-, s- or b-face) is closed, i.e. $B^u_{i,j,k+1} = 0$. If both neighbours at the w- and the t-w-side are open ($B^u_{i-1,j,k}=1$ and $B^u_{i-1,j,k+1}=1$), the face is marked with the colour green.

The configurations for the third and fourth rule are depicted in Fig. 3 for a t-face. Figure 4 displays two possible three-dimensional configurations of the momentum cell with colour coding.

Fig. 2

Cell-face types for a u-momentum cell. Each type is represented by a different colour

Fig. 3

Assignment of cell-face types (colour codes) at a T-face of a u-momentum cell (i, j, k) based on the indicator field $B^u$. The value of $B^u_{i,j,k}$ is 1 if the face area belonging to the velocity $u_{i,j,k}$ is open and 0 if it is closed. The shaded area schematically represents the original momentum cell (i, j, k)

Fig. 4

Cell-face types for a u-momentum cell in two three-dimensional configurations. Each type is represented by a different colour

2.2.3 Treatment of Small Cells

When an arbitrary surface is intersecting a Cartesian grid, one cannot avoid that arbitrarily small cells emerge. In the two-dimensional case, the stability of the skew-symmetric discretisation of Dröge and Verstappen (2005); Dröge (2006) is not affected by the presence of small cells. Contrarily, in the three-dimensional case small cells lead to large imaginary eigenvalues of the discrete convective operator matrix $\mathbf {V}^{-1}\mathbf {C}(\mathbf {u}_h)$ and thus result in a severe time step constraint (Dröge 2006). In our implementation, we use cell merging to avoid those instabilities. The cell merging approach described in Hartmann et al. (2011) is extended for staggered momentum cells.

A cell is considered to be small when the cell volume is smaller than $0.25 \varDelta x\varDelta y\varDelta z$, which turned out to be a reasonable compromise between accuracy and efficiency. Each small cell needs an appropriate master cell to combine with for which the largest direct neighbour is selected. The merging starts with the pressure cells and from the merged pressure cells the staggered momentum cells are formed. Figure 5 shows a simple example of a merged pressure cell and the associated staggered momentum cells. The velocity associated with a merged staggered momentum cell is located in the midpoint of the cell face of the merged pressure cell as shown in Fig. 5b and c.

The fluxes across the surface of a merged cell are computed by evaluating the fluxes on both the master and the small cells individually and then adding the results. For each sub-cell the same formulation is used as for a non-merged cut cell; however, the modified position of the merged velocity is taken into account when gradients are computed. After the global pressure correction, divergence-free volume fluxes are computed for all sub-cells by solving a local Poisson problem.

Fig. 5

Merging of small cells

2.3 Verification of the Cut Cell Implementation

As the scheme we used in this project has not been verified before, we present results for three test cases. The first test case establishes the energy conservation of the skew-symmetric formulation of the cut cell method. The second test case demonstrates that the implementation including cell merging is second order in all terms of the Navier–Stokes equation and the third test case places the method in relation to other schemes published.

2.3.1 Conservation of Energy for the Inviscid Taylor Vortex Flow Around a Body

In this section, we test the energy conservation properties of the cut cell implementation. Following (Morinishi 2010), we simulate an inviscid flow in a periodic domain. Due to the skew symmetry of the discrete convective operator and the adjoint relation between the discrete gradient and divergence operators, the kinetic energy must be conserved if the initial velocity field is divergence-free and satisfies the no-penetration boundary condition on the immersed boundary.

We construct the body geometry and the initial condition from the two-dimensional Taylor vortex (Taylor 1923) which remains stationary in the limit of zero viscosity. It is given by the stream function

$$\begin{aligned} \psi (x,y) = -\frac{U_0}{\kappa } \cos (\kappa x) \cos (\kappa y) \end{aligned}$$

(3)

with the wavenumber $\kappa$. We set the body surface at the isocontour of the stream function $\psi _\mathrm {body} = - U_0/(2\kappa )$, which extends from $\pi /3 \le \kappa x \le 2\pi /3$, $\pi /3 \le \kappa y \le 2\pi /3$. Since this is a closed streamline, the flow is always tangent to the body and thus the no-penetration wall boundary condition for inviscid flow is satisfied. The discrete velocity fields were initialised with the specific volume fluxes across the cell faces that can be easily computed from the stream function and the open face areas. The velocity field and the body are shown in Fig. 6a.

The flow was simulated for a total time $\kappa U_0 T=2\pi$. Four grids were considered that cover one wavelength of the flow with a cell height $\kappa \varDelta z= 2\pi /10$ and cell widths of $\kappa \varDelta x = \kappa \varDelta y = 2\pi /20$, $2\pi /40$, $2\pi /80$ and $2\pi /160$, respectively. We investigated multiple time step sizes between $U_0\,\varDelta t/\varDelta x=5\times 10^{-4}$ and $5\times 10^{-1}$. The simulations were performed in double precision floating point arithmetic.

Figure 6b displays the relative error in the conservation of kinetic energy over the time interval [0, T] for the different grid resolutions. It can be seen that the error depends on the time step size as $O (\varDelta t^3)$ and can be brought down to floating point accuracy for very small time steps. This error is caused by the numerical dissipation of the explicit third-order Runge–Kutta scheme. Moreover, a non-systematic dependency of the error on the grid spacing is apparent. Please note that we have also considered rotated variants of the case with identical results.

As the error in the conservation of the kinetic energy for any given grid can be reduced to machine precision by choosing a small enough time step, we conclude that the implementation of the convective and pressure gradient terms in the cut cell method indeed conserves the kinetic energy in an inviscid flow.

Fig. 6

Inviscid Taylor vortex flow around a body defined by $\psi = -U_0/(2\kappa )$. a Streamlines and u-component of the initial velocity field. b Relative error in the total kinetic energy $E_\mathrm {kin}=\frac{1}{2}\mathbf {u}_h^\mathrm {T} \mathbf {V}\mathbf {u}_h$ over the time interval $[0,T]=[0,2\pi /(\kappa U_0)]$. The dashed line indicates third order convergence with the time step $\varDelta t$

2.3.2 Convergence to an Analytical Solution for Transient Flow through an Oblique Pipe

In this section, we demonstrate the spatial convergence of the cut cell implementation for transient laminar flow in a circular pipe. This test case is a three-dimensional and unsteady generalisation of the oblique channel flow considered by Dröge and Verstappen (2005), Meyer et al. (2010). The pipe of radius a is oriented the in the (1, 1, 1)-direction and is embedded into a cubic domain of length 4a with triple periodic boundary conditions (cf. Fig. 7a). A pressure gradient is applied along the x-direction. For this flow, an analytical solution exists in the form of a series expansion (Pozrikidis 2017):

$$\begin{aligned} u_s(r, t) = -\frac{1}{4\mu }\frac{1}{\sqrt{3}}\frac{\mathrm {d}p}{\mathrm {d}x}\left[ a^2 - r^2 - 8 a^2 \sum _{n=1}^{\infty } \frac{1}{\alpha _{n}^3}\frac{\mathrm {J}_0(\alpha _n r/a)}{\mathrm {J}_1(\alpha _n)} \exp \left( -\alpha _n^2 \frac{\nu t}{a^2}\right) \right] \end{aligned}$$

(4)

where $\mathrm {J}_0$ and $\mathrm {J}_1$ are the Bessel functions of the first kind and $\alpha _n$ is the n-th zero of the Bessel function $\mathrm {J}_0$. The dimensionless pressure gradient magnitude was chosen as $-\mathrm {d}p/\mathrm {d}x\,a^3/( \rho \nu ^2 )=20$. This would eventually result in a steady state Reynolds number $Re =u_b\,2 a/\nu =2.89$ with the bulk velocity $u_b$.

The oblique orientation of the pipe ensures that the convective terms in the Cartesian frame are nonzero. The pressure gradient (body force) is applied in the (1, 0, 0)-direction and forms a $54.7^\circ$ angle with the pipe axis. This has the consequence that the pressure iterations need to redirect this pressure gradient and that the pressure gradient terms in the momentum equation are nonzero. Consequently, in this test case all terms of the incompressible Navier–Stokes equations are active. Moreover, the oblique placement of the pipe with respect to the Cartesian grid leads to a large variety of cut cell geometries and sizes.

We simulated the flow with 4000 time steps until $\nu T/ a^2 = 0.04$ was reached. At this time the viscous layer thickness is $20\%$ of the pipe’s radius. We used cubic cells with a size of a/8, a/16, a/32, and a/64, respectively. The time step $\nu \,\varDelta t/a^2=10^{-5}$ was identical for all cases and was chosen small enough such that the temporal error is negligible compared to the spatial error. We sampled the velocity field at probe points with $r/a = 0, 0.1, \dots , 0.9$ and 36 points in azimuthal direction. For each radial position, we computed the maximum error with respect to the analytical solution. Figure 7b shows the maximum error at each radius over the grid spacing; we observe that the error decreases towards the centre of the pipe and we see a second order convergence with the grid spacing. Please note that the computations were performed in single precision arithmetic.

Fig. 7

Transient flow through an oblique pipe. a Configuration of the pipe in a periodic domain and the applied pressure gradient. b Grid convergence of the $L^\infty$-norm of the streamwise velocity $u_s(r, T)$ at probe points with $r/a = 0, 0.1, \dots , 0.9$. The dashed line indicates second order convergence with the grid spacing $\varDelta x$

2.3.3 Steady Flow Around a Sphere at Re = 100

Following (Kim et al. 2001; Marella et al. 2005; Hartmann et al. 2011), we simulated steady laminar flow around a sphere of diameter D at $Re = U_0\,D/\nu =100$, where $U_0$ denotes the velocity of the uniform approach flow. This flow separates from the sphere and forms a closed recirculation region.

The extent of the simulation domain was chosen as $[-15 D; 49 D]$ in the x-direction and $[-16 D; 16 D]$ for the y- and the z-direction. Like (Hartmann et al. 2011), we used a zonal grid with a total of 5 levels with a cell size $\varDelta x = 0.5 D$ on the coarsest and $\varDelta x = 0.03125 D = D/32$ on the finest level (grid #1). Moreover, we performed a full refinement by a factor of 2 (grid #2). Due to the block-structured grid in our code compared to the octree grid of Hartmann et al. (2011), we had to modify the grid layout. The grid definition is given in Table 1. A constant velocity $\mathbf {u}=U_0\,\mathbf {e}_x$ was set at the inlet and $p=0$ and homogeneous Neumann boundary conditions for the velocity were set at the outlet. At the sides of the domain, slip boundary conditions were applied. The velocity field was initialised with $\mathbf {u}=U_0\,\mathbf {e}_x$ and the flow was integrated over a time $U_0\,T/D = 100$ with 16000 and 64000 time steps for grids #1 and #2, respectively.

In Table 2 the friction drag coefficient $C_\mathrm {f}$, the pressure drag coefficient $C_p$, the total drag coefficient $C_\mathrm {d}=C_\mathrm {f}+C_p$ and the length of the recirculation zone from our simulations are compared to results from the literature. The friction and pressure drag coefficients were computed from the wall shear stress and the pressure in the cut cells as

$$\begin{aligned} C_\mathrm {f}&= \frac{1}{\frac{1}{2}\rho U_0^2 \pi \frac{D^2}{4}} \int _{A_\mathrm {sphere}} \tau _{\mathrm {w}x}\,\mathrm {d}A \end{aligned}$$

(5)

$$\begin{aligned} C_p&= \frac{1}{\frac{1}{2}\rho U_0^2 \pi \frac{D^2}{4}} \int _{A_\mathrm {sphere}} -p n_x\,\mathrm {d}A \end{aligned}$$

(6)

The length of the recirculation zone $L_\mathrm {r}$ was defined as the distance from the rear end of the sphere at $x=0.5 D$ to the zero crossing of the u-component on the x-axis. It can be seen from Table 2 that the results of our simulation lie within the range of values in the literature.

Moreover, we investigated the dependency of the time step on the cell merging threshold that represents the open volume fraction of a cut cell below which a cut cell is merged to one of its neighbours. Figure 8 displays the maximum time step for which the numerical solution remains stable over an integration time of $\varDelta x/U_0$ and the minimum time step for which the numerical solution blows up as a function of the cell merging threshold. For a cell merging threshold larger or equal than 0.25, the time step is limited by the diffusion number of an uncut cell. For smaller cell merging thresholds (leading to fewer merged cells and smaller cut cells), the maximum stable time step is reduced severely. As mentioned earlier we therefore decided for a threshold value of 0.25.

Table 1 Grid configuration for the flow around a sphere at $Re =100$. Grid #1 contained 3.3m cells and grid #2 contained 26m cells

Table 2 Comparison of results for the drag coefficent $C_\mathrm {d}$, the friction drag coefficient $C_\mathrm {f}$, the pressure drag coefficient $C_p$, and the length of the recirculation zone $L_\mathrm {r}/D$ for the flow around a sphere at $Re =100$. The abbreviations FD and FV stand for finite difference and finite volume method, respectively

Fig. 8

Dependency of the maximum stable time step on the cell merging criterion for flow around a sphere at $Re =100$ (grid #1). The time step is normalised with the diffusion number limit of an uncut cell; the numerical value 2.5127 represents the stability limit of the Runge–Kutta scheme along the negative real axis

Fig. 9

Geometry of the scour hole and definition of the global coordinate system

3 Flow Configuration, Experimental and Simulation Setup

3.1 Flow Configuration

The flow under investigation is a fully developed turbulent flow in a water channel (width $b=1.17 m$) with a free surface and a sandy bed (water depth $h=0.15 m$). We performed an experiment at the hydraulic laboratory of the Professorship of Hydromechanics at the Technical University of Munich. We placed a cylinder (diameter $D=0.1 m$) vertically in an open channel. This serves as a model for a bridge pier in a river with a mobile bed. In a preliminary study, we let a scour hole develop in a sediment bed with a grain diameter $d=2 mm$ and measured its geometry after certain development times (Pfleger et al. 2011; Pfleger 2011). The measurements were performed using a laser distance sensor on a $5\,mm$ grid. For the measurements used herein, the scour geometry after one hour of development was taken and milled into aluminium to provide a stable and well controlled surface which could also be taken as a boundary for the LES. The geometry of the scour hole is depicted in Fig. 9 and is provided in the supplementary data. The surface of the scour cannot be regarded as smooth. Instead, there are undulations at a length scale between $2\,mm$ and 5mm resulting from the grain size and the laser scanning resolution. The main flow is along the x-direction. The y-direction is horizontal and the z-direction is vertical and constitutes the axis of the cylinder.

The study targeted two nominal Reynolds numbers, 20, 000 and 40, 000 which are close to the ones used by Dargahi (1990) and Jenssen et al. (2021) for the case of the flow of a cylinder mounted on a flat plate. These target Reynolds numbers were not fully met in either experiment and simulation. In the experiment, the precision of the discharge measurement was limited. The simulation was performed at exactly the same friction Reynolds number as in Schanderl and Manhart (2016), but the flow rate changed slightly due to a change in the grid spacing compared to the latter simulation. When using the depth-averaged velocity in the symmetry plane as reference, the Reynolds numbers are $Re _{D}=u_\mathrm {sym}D/\nu =21,256$ and 39, 662 in the experiment and $Re _{D}=22,618$ and 45, 958 in the LES, respectively. The Reynolds numbers based on the cross-sectionally averaged velocity $u_\mathrm {cs}$ were $Re _{D,\mathrm {cs}}=21,251$ and 41, 342 for the LES and 20, 331 and 38, 806 for the experiment. As the simulations and the experiments were conducted in parallel, the discrepancy in Reynolds numbers could not be compensated for. As the flow changes with Reynolds number are relatively small in this Reynolds number range (Jenssen et al. 2021), we will refer to the nominal Reynolds numbers in the remainder of this paper.

3.2 Experimental Setup

In this section, we describe selected aspects of the experimental setup. A detailed description of the experiment can be found in Jenssen and Manhart (2020). The measurements were done by a stereoscopic PIV. A streamwise/vertical plane was illuminated by a laser light sheet entering the water body from above. An acrylic glass plate was applied at the water surface in front of the cylinder in order to prevent refractions of the light sheet by surface waves. Thus, the bow wave in front of the cylinder was suppressed. Consequently, despite of the relatively high Froude numbers of $Fr =u_\mathrm {sym}/\sqrt{g h}=0.17$ and 0.33 for the low and high Reynolds number, respectively, the water surface can be considered as flat. The effect of this glass plate onto the flow was investigated numerically for the flat bed case at $Re _{D,\mathrm {cs}}=20,000$ in Alfaya (2016). It was found that near the horseshoe vortex the velocity and the turbulence statistics are only weakly affected by the presence of the plate.

The PIV images were taken by two cameras looking from above at a $45^\circ$ angle. The optical access was realised through water-filled prisms at the water surface. After the usual stereoscopic transformations the $16 \times 16$ pixel interrogation window size was $5.88 \times 10^{-3}\,D=0.588 mm$ with an overlap of $50 \%$. In order to reduce the light reflections the scour surface was coated by a fluorescent Rhodamine-B varnish which returned light in a shifted wavelength which was filtered out by a band-pass filter in front of the cameras. However, the efficiency of the Rhodamine-B varnish was not perfect and reflections were still present in a layer along the scour surface affecting the measurements at the wall. Considerable difficulties were associated with locating the wall position in the acquired images with an accuracy related to the pixel size ($\approx 39 \upmu m$). Therefore, the position of the wall was placed in the centre of the reflection layer which resulted in an uncertainty of about 25 pixels distance. For the computation of the wall shear stress, the wall position was refined by assuming a linear variation through the five wall-nearest points (Jenssen 2019). The wall shear stress was then determined from the velocity vector closest to the wall. Please note that in the flat bed case the wall shear stress obtained in this way was underpredicted compared to the LES and a single-pixel PIV evaluation (Schanderl et al. 2017; Strobl 2017; Jenssen et al. 2021).

3.3 Simulation Setup

As in the experiment, the cylinder and the scour hole were placed inside an open channel of width $b=11.7\,D$ and height $h=1.512\,D$. The domain consisted of a precursor domain which was used to generate a fully turbulent inflow condition and the main domain. The precursor grid had a length $22.4\,D=14.8h$ with periodic boundary conditions applied in streamwise direction. The main simulation domain had a length of $22.4\,D$ as well with the cylinder placed in the centre ($x/D=0$). The time-dependent velocity profiles from the precursor simulation were imposed at the inflow plane located at $x/D=-11.2$. At the outflow plane, a zero-gradient condition was used for the velocities. The water surface was modelled by a free-slip condition, which gives a non-deformable free surface. This corresponds to the limit of zero Froude number, see also discussion in Sect. 3.2. At rigid walls (bottom wall, scour geometry, cylinder and side walls) the flow satisfies a no-slip boundary condition and the Werner-Wengle wall function was employed to compute the wall shear stress.

The flow in the precursor domain was driven by a constant pressure gradient such that $Re _\tau = u_\tau \,h/\nu =1522$ and 2762, respectively, with the wall shear velocity

$$\begin{aligned} u_\tau =\sqrt{\frac{\tau _\mathrm {w}}{\rho }} = \sqrt{-\frac{1}{\rho }\frac{\mathrm {d}p}{\mathrm {d}x}\frac{b\,h}{b+2 h}}. \end{aligned}$$

(7)

Fig. 10

Grid configuration around the cylinder and inside the scour hole. Every grid box sketched corresponds to $20^3$ cells at $Re _{D}=20,000$ and to $34^3$ cells at $Re _{D}=40,000$

Table 3 Details of the computational grid. The spatial extent of the refinement levels was the same for both Reynolds numbers. “+” indicates grid spacing based on the wall friction velocity $u_\tau$ in the precursor domain according to Eq. (7)

The flow domain was represented by a block-structured Cartesian grid with local refinement on five levels each reducing the cubic grid cell spacing by a factor of two. For both Reynolds numbers we performed a grid study by successively adding more refinement levels. The base grid (grid #1) covered the whole computational domain, grid #2 added refinement level #2 and so on. The grid was refined at the bottom wall, the cylinder and inside the scour geometry (see Fig. 10). Please note that the refinement also included the rim of the scour hole where a flow separation takes place. At the finest resolution, a total number of 478 million and 2.35 billion cubic grid cells was used for $Re _{D}= 20,000$ and 40, 000, respectively. In Table 3 the main parameters of the computational grids are given.

In a first step, the appropriateness of the grid resolution can be judged based on the grid spacing ${\varDelta x}^+= u_\tau {\varDelta x}/\nu$ in wall units. On the one hand, the channel flow in grid #1 is resolved with a cell size ${\varDelta x}^+ \approx 30$ wall units which constitutes a wall-modelled LES. As the SGS model and wall function have been calibrated and validated with channel flow, we expect little uncertainty in this region of the flow. On the other hand, the grid spacing in the scour hole in grid #5 is at most $\varDelta x^+ =5$ wall units with the friction velocity based on the local wall shear stress (cf. the peak friction coefficient in Fig. 19). As the Werner-Wengle wall function reduces to the linear law of the wall for ${\varDelta x}^+ \le 11.81$, the flow in the scour hole can be considered a wall-resolved LES. Hence, the refinement from grid #2 to #5 gradually moves the simulation from a wall-modelled to a wall-resolved LES. The adequacy of the grid resolution will be discussed in more detail in the Sects. 4.1 and 4.2.

The simulations were run until a statistically steady state was achieved before sampling was started. The samples cover a time of $450D/u_\mathrm {sym}$ at $Re _{D}=20,000$ and $118D/u_\mathrm {sym}$ at 40, 000, respectively. At the finest grid resolution, this led to a computational cost of 1.0 and 2.3 million CPU hours for the two Reynolds numbers.

4 Results

In this section, we first evaluate the reliability of the LES results by investigating the development of representative quantities with grid refinement and by comparing the finest grid results with the ones measured by PIV. Based on this analysis we discuss certain features of the flow. Furthermore, we put a focus on the three-dimensional shape of the horseshoe vortex inside the scour hole as well as the wall shear stress field, since these quantities are challenging to obtain experimentally.

4.1 Flow Topology

Before evaluating the ability of the LES to represent the flow topology in an adequate sense, we introduce the main characteristics of the flow in the scour hole in front of the cylinder. The flow inside the scour hole is dominated by the horseshoe vortex wrapping around the cylinder. This vortex is the only vortex in the symmetry plane upstream of the cylinder which can be made visible in the streamlines (Fig. 11) for the time-averaged velocity field from the experiment of Jenssen and Manhart (2020). The downflow in front of the cylinder is a result of the vertical gradient in the velocity profile of the approaching open channel flow which leads to a higher stagnation pressure at the upper part of the cylinder than at the lower part. This drives the flow downwards. At the bottom of the scour hole, we can find a stagnation point near $x/D = -0.55$. The flow upstream of the stagnation point is directed into the horseshoe vortex; the flow on the other side moves around the cylinder. The measurements could identify a secondary flow detachment between $x/D=-1.1$ and $-1.0$ at $Re _{D}=20,000$ and a local near-wall velocity peak directly underneath the horseshoe vortex. This velocity peak could give rise to a large wall shear stress; this will be discussed in Sect. 4.5.

Fig. 11

Time-averaged streamlines of the flow field inside the scour hole in front of the cylinder from the PIV experiment (Jenssen and Manhart 2020) super-imposed with the velocity magnitude normalised by $u_\mathrm {sym}$ for a $Re _{D}= 20,000$ and b 40, 000

Fig. 12

Time-averaged streamlines of the flow field inside the scour hole in front of the cylinder taken from LES super-imposed with the velocity magnitude. The “+” indicates the position of the horseshoe vortex observed in the experiment. Left column: grids #1 – #5 for $Re _{D}= 20,000$. Right column: grids #1 – #5 for $Re _{D}= 40,000$

In the following, we discuss the flow topology resulting from our simulations and its variation with grid resolution. Figure 12 depicts the time-averaged streamlines in the symmetry plane as well as the magnitude of the time-averaged velocity for the different grid resolutions. We can observe a strong variation of the flow topology with the grid spacing. In the coarsest simulations (LES #1) the flow remains attached downstream of the rim and only separates near $x/D= -1.15$. At all higher resolutions, it separates from the rim of the scour hole. At both Reynolds numbers a large single horseshoe vortex can be observed in the grids #2 and #3 at a distance of approximately 0.1D from the position observed in the experiment, which is marked by a “$+$” in the plots. In grid #4 the horseshoe vortex moves down into the scour hole closer to the position observed in the experiment. This gives space for a small secondary vortex in the detached shear layer upstream of the horseshoe vortex. A further grid refinement along the wall (grid #5) does not result in a significant change of the position of the horseshoe vortex. However the secondary vortex becomes weaker at the lower or vanishes at the higher Reynolds number. The shift of the vortex between grids #3 and #4 gives room for a small recirculation zone at the wall at $x/D=-1.05$ which also is visible in the measured flow field at $Re _{D}= 20,000$. While position and shape of the horseshoe vortex do not change significantly between grids #4 and #5, there remains some uncertainty about the secondary vortex that is visible in grid #5 at $Re _{D}= 20,000$. It has not been measured in the experiment. However, such a second vortex has been observed by Kirkil et al. (2008), Kirkil et al. (2009) in a deeper scour hole at equilibrium depth which is obtained after a long time of scour development. It might be conjectured that the change in flow topology that is observed between the LES #3 and #4 is due to a better resolution of the flow separation at the rim of the scour hole and the detached shear layer developing upstream of the main vortex.

Table 4 Main parameters of the horseshoe vortex. $x_\mathrm {HV}/D$ is the location of the horseshoe vortex

In Table 4, some global parameters of the horseshoe vortex obtained on different grids are compared to the experimental values. The horseshoe vortex positions were determined by the focus of the streamlines and confirmed with the swirl criterion of Kida and Miura (1998). The maximum vorticity and turbulent kinetic energy (TKE) were determined in the symmetry plane in the neighbourhood of the horseshoe vortex centre. From grid #2 on, there is a roughly monotonic trend in the horseshoe vortex positions with differences between the finest and the second finest grid below $2.5\%$ for both Reynolds numbers. The measured positions are close to the simulated ones at the finest grids. The maxima of the time averaged spanwise vorticity $\overline{\omega }_{y,max}$ and TKE in the horseshoe vortex centre show a similar trend although no monotonous convergence.

Generally, we can conclude that the simulation at the coarsest grid $\#1$ is far from a grid-converged solution. It is evident that there is no monotonic convergence of the global quantities, which show larger variations between resolutions differing in flow topology. An application of the Grid Convergence Index (GCI) (Roache 1997) would result in inconsistent values among the various grid refinements. However, the simulation at the second finest grid $\#4$ is close to the finest one #5 from which we can conclude that those simulations are close to a grid-converged solution. We would also like to emphasise here, that the statistics we could gather with our limited computing time are probably not fully converged.

Fig. 13

Time-averaged velocity profiles in the symmetry plane ($y=0$) at $Re _{D}= 40,000$. a u-velocity profile through the horseshoe vortex centre from the experiment, and b w-velocity above the scour hole ($z=0$). The grey shaded area represents the $2.8 \%$ uncertainty associated with the experimental results

4.2 Velocity Profiles

In this section we discuss the grid dependency of selected profiles of the mean velocity at the Reynolds number $Re _{D}= 40,000$. This will provide further evidence to judge the uncertainty in the simulations. We do not show this comparison for the lower Reynolds number as there were some uncertainties in the normalisation of the measured values.

First, we investigate a vertical profile of the time-averaged streamwise velocity component $\overline{u}$ through the measured horseshoe vortex centre. In Fig. 13a the velocity profiles obtained with the different grid resolutions are compared to the corresponding profile from the experiment (PIV). The backflow under the horseshoe vortex is well resolved in all grids. The agreement between grids #4 and #5 is fully satisfactory. Note that grid #5 added further refinement only along the wall compared to grid #4 (cf. Fig. 10). The simulations reach the measured streamwise velocity magnitude above the scour hole and there is a slight undulation in the profiles of grids #4 and #5 around $z/D=-0.1$, which we cannot explain at this moment. However, the vertical position of and the velocity gradient around the horseshoe vortex centre fully agree between the PIV and the two simulations with the finest grids.

Second, we present in Fig. 13b various simulated and measured horizontal profiles of the time-averaged vertical velocity $\overline{w}$ at the upper end of the scour hole ($z=0$). This quantity can be understood as the convective velocity transporting horizontal and vertical momentum into the scour hole. The curves indicate that (i) The results on the coarsest grid are far from the others, (ii) The results on the two finest grids are very close and (iii) The experimental results are not far from the simulated ones on the finest grid. We can also observe a qualitative deviation of the profiles grids $\#1$ – $\#3$ from the experiment and the simulations with the two finest grids. While the latter have a fairly monotonic decrease of $\overline{w}$ towards the maximum downflow close to the cylinder, the solutions on the intermediate grids show a local minimum of the downward velocity near $x/D=-1.1$. This seems to be the footprint of the larger and upward shifted horseshoe vortex predicted by grids #2 and #3. While the peak amplitude of the downward velocity agrees between the simulations and the experiment, the experiment shows a stronger downflow than the simulations around $x/D=-1.0$. Nevertheless, the shape of the profiles predicted by grids #4 and #5 agrees quite well with the measured one.

In conclusion, we can observe good quantitative agreement between results obtained on grids #4 and #5 for both velocity profiles and with the results presented in Sect. 4.1, we consider those simulations sufficiently converged in terms of grid resolution. In the following, we only discuss the results for the simulations on grid #5.

4.3 Vorticity and Turbulence Quantities in the Symmetry Plane

In this section, we compare the spatial distributions of the mean vorticity, TKE, Reynolds stresses and production of TKE in the symmetry plane between LES and experiment at $Re _{D}= 40,000$.

Fig. 14

Time-averaged vorticity $\overline{\omega _y} = \frac{\partial \overline{w}}{\partial x} - \frac{\partial \overline{u}}{\partial z}$ normalised by $u_\mathrm {sym}/D$ inside the scour hole in front of the cylinder taken from a LES (grid $\#5$), and b PIV at $Re _{D}= 40,000$. The “+” indicates the position of the horseshoe vortex core observed in the LES and PIV, respectively

The distributions of $\overline{\omega _y}$ in the symmetry plane show a high level of agreement between LES and PIV (Fig. 14). Large values of negative vorticity (corresponding to a clockwise rotation) can be found in the shear layer detaching from the scour rim and around the horseshoe vortex centre. Large positive vorticity is found upstream of the horseshoe vortex at the wall. This can be attributed to the small recirculation region visible in the streamlines of the LES near $x/D=-1.05$ (Fig. 12).

Fig. 15

Turbulent kinetic energy (a,b), Reynolds normal stresses $\overline{u'u'}$ (c,d), Reynolds normal stresses $\overline{w'w'}$ (e,f) and production of TKE (g,h) inside the scour hole in front of the cylinder taken from LES (grid $\#5$) (a,c,e,g) and PIV (b,d,f,h) at $Re _{D}= 40,000$. Data are normalised by $u_\mathrm {sym}/D$. The “+” indicates the position of the horseshoe vortex observed in the LES and PIV, respectively

Figures 15a and b compare the TKE predicted by our LES with the measured one. In both data sets, a TKE maximum can be found at the vortex centre with values of 0.08 to $0.09 u_\mathrm {sym}^2$. The distributions are very similar to each other. However, the near-wall peak at $x/D=-1.0$ is better resolved by the LES, which can be attributed to the better near-wall resolution of grid #5 compared to the interrogation window size in the experiment (cf. Sect. 3.2). It is important to note that in the deepest region near the cylinder, the TKE is almost zero.

The detailed turbulence structure in front of the cylinder has been discussed by Jenssen and Manhart (2020). They showed that when the coordinate system is rotated to align with the surface of the scour hole, the turbulence structure is pretty similar to the one in front of a cylinder on a flat bed. This includes a ”C” -shaped distribution of the TKE in which the region around the horseshoe vortex is dominated by wall-normal fluctuations and the foot of the “C” is dominated by wall-parallel fluctuations. As shown in Fig. 15c–f, the distributions and levels of the streamwise and vertical Reynolds normal stresses are well predicted by grid $\#5$ in comparison to the PIV. The same can be concluded with some reservations from Fig. 15g,h for the production of TKE, $P=-\overline{u'_i u'_j}\, \partial \overline{u_i}/\partial x_j$. The PIV measured the largest levels of TKE production upstream of the horseshoe vortex centre and could not resolve the near-wall peak which is well visible in the simulated results. The same could be observed in the flat-plate case (Schanderl et al. 2017), so we are convinced that this is not an artifact of the simulation.

In summary, a highly satisfactory agreement the mean vorticity as well as second order velocity statistics can be observed. This suggests that the LES represents a flow very similar to the one measured in the experiment.

4.4 Mean and Instantaneous Appearance of the Horseshoe Vortex

In this section, we discuss the shape of the time-average horseshoe vortex based on streamlines, contours of the mean pressure and the TKE. Furthermore, we investigate the appearance of the horseshoe vortex at $Re _{D}= 40,000$ visible in instantaneous shapshots of the pressure distributions.

Fig. 16

Side view (a,b) and top view (c,d) of the mean pressure isosurface with $\overline{p}-\overline{p}_\mathrm {ref}= -0.72 \rho u_\mathrm {sym}^2$ and streamlines near the cylinder with scour hole for (a,c) $Re _{D}= 20,000$ and (b,d) $Re _{D}= 40,000$. The reference pressures were computed at the location $x/D=-0.5$, $y/D=0$ and $z/D=0.5$

Fig. 17

Top view of TKE isosurface with $0.05 u_\mathrm {sym}^2$ and streamlines near the cylinder with scour hole for a $Re _{D}= 20,000$ and b $Re _{D}= 40,000$

Fig. 18

Top view of instantaneous pressure isosurfaces with $p-\overline{p}_\mathrm {ref}= -0.72 \rho u_\mathrm {sym}^2$ and streamlines near the cylinder with scour hole for different times at $Re _{D}= 40,000$. The reference pressure was computed at the location $x/D=-0.5$, $y/D=0$ and $z/D=0.5$

Figure 16 shows isosurfaces of the mean pressure and streamlines of the mean velocity field in side and top view, respectively. The streamlines were seeded on the line $x/D=-1.6$, $z/D=0.05$. A low pressure region in the core of the horseshoe vortex consistent with the streamlines is clearly visible. In front of the cylinder, the horseshoe vortex core remains at an approximately constant depth of $z/D=-0.2$. The elevation of the core increases as the horseshoe vortex passes around of the cylinder. Finally, the horseshoe vortex leaves the scour hole around $x/D = 0.3$ for the lower and $x/D=0.2$ for the higher Reynolds number. Another interesting feature of the isosurfaces is that the length of the low pressure region behind the cylinder decreases with the depth. This suggests that the flow separation from the cylinder is less pronounced near the scour hole.

We also determined the core line of the horseshoe vortex from the mean velocity field using the parallel vectors method of Roth and Peikert (1998), Haimes and Kenwright (1999) implemented in ParaView. The core lines lie in the centre of the mean pressure isosurfaces (cf. Fig. 16). We provide the isolated and sorted core lines for both Reynolds numbers as supplementary data.

The horseshoe vortex is known to be a region of high turbulence intensity (Devenport and Simpson 1990; Jenssen and Manhart 2020). Therefore, also the spatial distribution of the TKE can provide further evidence on the shape of the horseshoe vortex. Figure 17 shows an isosurface of TKE and streamlines of the mean velocity field. It is clearly visible that the TKE is concentrated in the horseshoe vortex and in the wake of the cylinder. A peculiar feature of the TKE surface associated with the horseshoe vortex is the indentation at both ends of the vortex. These indentations appear approximately where the surface reaches the elevation of the channel bed ($z/D \approx -0.1$). It can be seen from both the shape of the TKE isosurface and the streamlines that the width of the cylinder wake and the distance between the branches of the horseshoe vortex is larger for $Re _{D}= 40,000$ than for 20, 000. Moreover, the comparable volume enclosed by the isosurfaces indicates that the values of the TKE for both Reynolds numbers are quite similar when normalised with $u_\mathrm {sym}^2$.

In the following, we analyse four snapshots of the flow at $Re _{D}= 40,000$ that were sampled with a temporal separation of at least $12 D/u_\mathrm {sym}$. Figure 18 shows snapshots of an isosurface of the instantaneous pressure with the same contour value as in Fig. 16 together with instantaneous streamlines. The plots indicate that the horseshoe vortex is also present in the instantaneous flow fields. The isosurfaces at the location of the horseshoe vortex also show smaller turbulent vortices superimposed on the horseshoe vortex. Low pressure regions can be observed on the sides and in the wake of the cylinder; in some snapshots the vortices of the von Kármán vortex street are clearly discernible. There seems to be little variation in the location of the horseshoe vortex as it remains close to the mean position in all the snapshots. It can be seen that sometimes the pressure isosurface near the horseshoe vortex is broken into smaller elements. A further breakup of the vortex can be observed near the sides of the cylinder. A possible interpretation of this is that the vortex loses its coherence and disintegrates into smaller turbulent structures.

Please note that the Q-criterion (Hunt et al. 1988) that has been employed in comparable studies to visualise the coherent structures in the vicinity of the horseshoe vortex (Kirkil et al. 2008); Link et al. (2012) turned out not to be appropriate to isolate the horseshoe vortex in our simulations. The reason for this is that as the Q-criterion is proportional to the Laplacian of the pressure which highlights much finer scales than the pressure field itself. Due to the high Reynolds numbers and the comparably much finer grid, there appears to be a scale separation between the horseshoe vortex and the sea of turbulent vortices identified by the Q-criterion (with a strong clustering inside the scour hole.

4.5 Wall Shear Stress

The wall shear stress (defined as the stress $\tau _{\mathrm {w}}$ exerted on the wall) is the key quantity for modelling the sediment transport in the scour hole (Roulund et al. 2005) and large efforts have been made to determine it in experimental studies (Melville 1975; Graf and Istiarto 2002; Dey and Raikar 2007) and by numerical simulations (Dey and Bose 1994; Kirkil et al. 2008). There is a large scatter in reported experimental values which have been given in terms of [Pa] and when transformed into a friction coefficient $c_\mathrm {f}=\overline{\tau }_{\mathrm {w}}/\frac{1}{2}\rho u_{\mathrm {sym}}^2$ are around $c_f \lesssim 3\times 10^{-3}$ (Graf and Istiarto 2002) and $c_f\lesssim 0.15$ (Dey and Raikar 2007). Maximum values found in the LES by Kirkil et al. (2009) are $c_f\sim 0.013$. In our previous experimental investigation (Jenssen and Manhart 2020) we have determined the wall shear stress from the velocity gradient between the wall-nearest point and the wall. The values obtained by this procedure, $c_f\lesssim 0.01$, can be regarded as lower limit only. We have seen that our LES provide us with a much better near-wall resolution in the scour hole than the PIV. The experiment uses a spatial resolution of $16\times 16$ pixels in an interrogation window which compiles into a size of 0.608mm, i.e. D/165. This is considerably larger than the wall resolution achieved in grid #5 (D/1000) and lies between the resolutions of grids #2 and #3. Thus, we can expect that the values provided by our LES establish a better estimation of the wall shear stress in the scour hole.

Fig. 19

Friction coefficient $c_{\mathrm {f}}=\overline{\tau }_{\mathrm {w}x}/\frac{1}{2}\rho u_{\mathrm {sym}}^2$ in the symmetry plane from LES and PIV for a $Re _{D}= 20,000$ and b $Re _{D}= 40,000$

In Fig. 19 we compare the time-averaged friction coefficients $c_\mathrm {f}$ from LES on grid #5 with the experimental estimations. For the lower Reynolds number the computed magnitudes are about two times larger than the measured ones. This factor increases at the higher Reynolds number which is not surprising as the limited near-wall resolution of the experiment is more serious. There is a high level of scatter in the computed wall shear stress. The scatter with a wavelength of $\varDelta x$ can be attributed to the cut cell method whereas the scatter with a wavelength of $0.05\,D$ is a consequence of the piecewise triangular representation of the scour geometry obtained from the experiment by Pfleger (2011), see Sect. 3.1.

There is a very short region with a positive wall shear stress (related to near-wall flow in the positive x-direction) in the corner between cylinder and scour surface (at $x/D=-0.5$). This could be the trace of a small anticlockwise corner vortex. Moving away from the cylinder (at $x/D\approx -0.8$, see Fig. 12). The edge can explain the formation of the sharp peak in the magnitude of the wall shear stress which has not been observed in such an intensity in the flat plate case (Schanderl et al. 2017). It is likely that the sediment material is eroded at the edge and plateau and lifted up to be entrained by the horseshoe vortex. The wall shear stress between the peak and the cylinder opposes gravity and can therefore stabilise the steep slope of the scour hole. Comparing the wall shear stress distributions with the streamlines (Fig.  12), we also can identify the footprint of the small recirculation region near $x/D=-1.1$, which leads to a short region where the time-averaged wall shear stress is positive. The erosive potential in this region is higher than the magnitude of the wall shear stress indicates, as the gravity forces together with the positive wall shear stress can lead to local avalanches of bed material in this region. This process has been described by e.g. Dargahi (1990).

Furthermore, we document the grid-resolution dependency of the friction coefficient for $Re _{D}= 40,000$ in Fig. 19b. As could be expected from the discussion of the flow topology (cf. Sect. 4.1), there is a qualitative difference between the distribution of $c_\mathrm {f}$ on grid #3 on the one hand and on grids #4 and #5 on the other hand. On grid #3 the sharp wall shear stress peak at $x/D\approx -0.8$ and its magnitude around the plateau ($-0.9< x/D < -0.8$) are considerably smaller. On the other hand the LES on grid #3 produces large negative values at the region where the finer grids have small positive values at $x/D\approx -1.05$. This can be attributed to the wrong prediction of the horseshoe vortex position and topology on the coarser grids. Overall there is a strong dependence of the predicted wall shear stress on the grid resolution at the wall. The grid spacing of grid #4 does not yet seem sufficient to give grid-converged levels of the wall shear stress inside the scour hole. For $Re _{D}= 40,000$ this means that a grid spacing of D/1000 or even finer is required in wall-normal direction. This poses an extreme challenge to Cartesian grid methods as grid refinement has to be done in all three Cartesian directions. Body-fitted grids would have an advantage in this aspect.

Fig. 20

Distribution of the friction coefficient $c_{\mathrm {f}}=|\overline{\mathbf {\tau }}_{\mathrm {w}}|/\frac{1}{2}\rho u_{\mathrm {sym}}^2$ over whole scour hole (top view) from LES a $Re _{D}= 20,000$ and b $Re _{D}= 40,000$. The yellow line represents the core line of the time-average horseshoe vortex (cf. Sect. 4.4)

Fig. 21

Distribution of the magnitude of the component-wise root mean square value of the wall shear stress fluctuations $|{\tau }_{\mathrm {w},\mathrm {rms}}^\prime |/\frac{1}{2}\rho u_{\mathrm {sym}}^2$ over whole scour hole (top view) from LES a $Re _{D}= 20,000$ and b $Re _{D}= 40,000$. The yellow line represents the core line of the time-average horseshoe vortex (cf. Sect. 4.4)

In Fig. 20 we document the spatial distribution of the friction coefficient obtained on grid #5. The distribution is similar for both Reynolds numbers and the magnitude of $c_\mathrm {f}$ is slightly smaller for the larger Reynolds number than for the lower Reynolds number. Please note that the triangular patches that can be made out in Fig. 20 originate from the geometric representation of the scour hole which is based on a grid with size $0.05\,D$. The cell size of the LES in these regions is 28 and 50 times smaller for the lower and higher Reynolds number, respectively. As the scour geometry is composed of flat triangular pieces, the wall shear stress varies little within a triangle and strongly across the boundaries between triangles.

The wall shear stress is enhanced in a ring-like region following the core of the horseshoe vortex which is marked as a thick yellow line in the plots. At the sides of the cylinder a region of increasingly high wall shear stress emerges and extends up to an angle of approximately $145^\circ$ from the negative x-axis. On the sides of the cylinder the wall shear stress acts mostly in the positive x-direction. This suggests the possibility that sediment erosion can take place on the side of the cylinder and that the eroded material would be transported downstream. Apart from the spatial variations due to the triangulation of the scour surface, our distributions are similar to time-averaged wall stress distributions reported by Kirkil et al. (2008).

The spatial distribution of the standard deviation of the wall shear stress is documented in Fig. 21. The most intense fluctuations can be observed directly behind the cylinder, where the mean wall shear stress is close to zero. This is due to the alternating vortex shedding from the cylinder in this region. Moreover, fluctuations with a standard deviation of approximately $30-50 \%$ of the mean wall shear stress can be observed below the horseshoe vortex (especially below the outer branches). The location of these fluctuations is consistent with the near-wall maximum of the turbulent kinetic energy (cf. Fig. 15a), the Reynolds normal stress $\overline{u'u'}$ (cf. Fig. 15c) and of the production of turbulent kinetic energy (cf. Fig. 15g). The turbulent dynamics under the horseshoe vortex seem to give rise to increased wall shear stress fluctuations which act in the plateau region of the scour geometry and amplify the already large mean magnitudes of the wall shear stress in this area. On the other hand, the wall shear stress fluctuations are small between the horseshoe vortex and the cylinder as well as on the sides of the cylinder. It seems that these deepest parts of the scour hole are excavated by the action of the time-averaged wall shear stress caused by the strong velocity overshoot at the sides of the cylinder.

5 Conclusions

We performed highly resolved LES to investigate the flow inside a scour hole around a circular cylinder using a cut cell immersed boundary method. The scour hole which was originally determined experimentally was fixed at a depth corresponding to about a half of the equilibrium scour depth which establishes after a very long time.

A detailed description of the cut cell method is provided. It is a three-dimensional generalisation of the method by Dröge and Verstappen (2005); Dröge (2006) and conserves mass and momentum on consistent balance volumes. The cut cell approach allows to directly determine the wall shear stress and wall pressure consistently with the momentum balance. The convective term preserves the kinetic energy which makes the scheme well suited for scale-resolved simulations. A cell merging scheme is used to circumvent the time step constraint imposed by small cut cells. The energy conservation property of the convective and pressure gradient terms was demonstrated for an inviscid flow case. Moreover, a second-order convergence to an analytical solution was shown for a transient oblique pipe flow. Finally, we simulated flow around a sphere at $Re =100$ and the drag coefficient and length of the recirculation zone lie close to the values from the literature.

In this contribution, we examined a series of simulations with increasingly fine grids and at two Reynolds numbers $Re _{D}=20,000$ and 40, 000. We first compare the LES results to a stereoscopic PIV experiment (Jenssen and Manhart 2020) of the identical configuration. We clarified how close to the experimental values we can get with the LES and which numerical resolution is required to obtain results that are accurate enough to draw conclusions on flow physics. Furthermore, we investigated the three-dimensional shape of the horseshoe vortex in the instantaneous and time-averaged fields. Finally, we discussed the wall shear stress field and compared it to results from the literature.

In the simulations presented herein, we used a refinement by zonal grids within a Cartesian solver. Therefore, refinement was not limited to the immediate wall layer. The grid spacings in the most refined regions varied from 0.028D (36 cells per diameter) for the coarsest grid at $Re _{D}= 20,000$ to 0.001D (1000 cells per diameter) for the finest grid at $Re _{D}= 40,000$. We found that the topology of the time-averaged flow and the position of the horseshoe vortex are very sensitive to the grid resolution – a finding that should be valid also for other numerical methods. For a more detailed representation of the flow physics and a quantitative prediction, our results indicate that at least the second finest grid – using 285 cells per diameter at the lower and 500 cells at the higher Reynolds number – was necessary. This results in a relatively large computational effort of around $10^8$ cells for the lower and $5\times 10^8$ cells for the higher Reynolds number, respectively. For an accurate prediction of the wall shear stress inside the scour hole the finest resolutions near the wall are definitely required. This can be explained by the extremely thin viscous wall layer forming in such flows, which has been demonstrated for the flow around a cylinder mounted on a flat plate (Schanderl and Manhart 2016; Schanderl et al. 2017; Jenssen et al. 2021). It is also unlikely that traditional wall models could resolve this problems (Schanderl et al. 2017). It should be kept in mind that the statements on necessary grid resolution are highly dependent on the numerical order of accuracy of the code and the grid structure, topology and refinement strategy. One has to bear in mind, however, that the grid resolution should be fine enough to resolve the shear layer separating from the scour rim, the horseshoe vortex and its dynamics, the velocity field between the horseshoe vortex and the wall and the wall shear stress. A refinement strategy that concentrates solely on the wall would fail in our point of view as the development of the shear layer and the interaction of the horseshoe vortex with entrained small scale vorticity would not be captured.

Based on the LES results, the three-dimensional shape of the horseshoe vortex was investigated. We provide the core line of the time-average horseshoe vortex as supplementary data. Moreover, the horseshoe vortex has a clear signature in the instantaneous pressure field. The low pressure region of the horseshoe vortex is visible in instantaneous pressure fields although it doesn’t always seem to render a fully coherent vortex structure.

We have analysed the wall-shear stress distribution inside the scour hole. This quantity is strongly grid dependent and needs a proper representation of the flow topology and a proper wall-normal grid resolution. The distribution is qualitatively similar to distributions found in the literature. We found that reported results for LES at low Reynolds numbers are of a similar order of magnitude as ours, although a direct comparison is not possible as our scour hole geometries and depths are not the same as the ones used in other references. Our distributions can be used to explain the character of the erosion process in the initial phase of the scour development in which near-cylinder material is transported around the cylinder, and in the later phases in which material is picked up under the horseshoe vortex and entrained into the same. The sediment avalanches observed in the upper parts of the scour hole are supported by downward pointing wall shear stresses found in the secondary recirculation zone upstream of the horseshoe vortex.

On a final note, we would like to point out that at the grid resolution $\varDelta x=0.0021\,D$ which we found to be the minimum requirement for reliable results at a Reynolds number $Re _{D}=40,000$, each sand grain in the original experiment of Pfleger et al. (2011), Pfleger (2011) would be resolved by 10 cells. This resolution has for example been employed in Mazzuoli et al. (2019), Scherer et al. (2020) to perform grain-resolved direct numerical simulation of sediment erosion and transport in turbulent channel flow. Therefore, a direct simulation of the scouring process in the present configuration seems to be within the reach of currently available high performance computing technology.

6 Supplementary Data

We provide the core lines of the horseshoe vortex obtained from the mean velocity field of our LES at $Re _{D}=20,000$ and 40, 000. The dataset includes the arclength s/D and the coordinate values x/D, y/D and z/D of the points on the core line. Moreover, we provide the scour geometry that has been used in the simulations and experiments.

Appendix

1.1 Convective Fluxes

In this section, we describe the computation of the convective fluxes in the u-momentum equation. The equations for v and w follow by interchanging u with the respective velocity component and i with either j or k.

For the e- and w-cell faces (orange and yellow), the convective fluxes are written as the product of the volume flux across the momentum cell face and a transported velocity interpolated with the weight $\frac{1}{2}$

$$\begin{aligned} \left. \left( \mathbf {C}(\mathbf {u}_h)\mathbf {u}_h\right) _{i,j,k}\right| _{\textsc {w}}= -\underbrace{\frac{1}{2} \left( \overline{u}_{i-1,j,k}+\overline{u}_{i,j,k}\right) }_{\text {volume flux}} \underbrace{\frac{1}{2} \left( u_{i-1,j,k}+u_{i,j,k}\right) }_{\text {transported velocity}} \end{aligned}$$

(8)

The interpolation of the volume flux is chosen such that the transporting velocity field is divergence-free on the momentum cell. The weight $\frac{1}{2}$ in the interpolation of the transported velocity follows from the skew-symmetry requirement of Verstappen and Veldman (2003).

The cell faces parallel to the velocity vector of the momentum cell are denoted as n- s-, t- and b-boundaries with the possible colours pink, green, grey and black. We compute the convective flux by summing over all colours:

$$\begin{aligned} \left. \left( \mathbf {C}(\mathbf {u}_h)\mathbf {u}_h\right) _{i,j,k}\right| _{\textsc {t}} = \sum _{\mathrm {colour}} \left( \overline{w}^\mathrm {colour}_{i,j,k} u^\mathrm {colour}_{i,j,k} + \overline{w}^\mathrm {colour}_{i+1,j,k} u^\mathrm {colour}_{i+1,j,k}\right) \end{aligned}$$

(9)

The definition of the volume fluxes and transported velocities is detailed in Table 5. For every half-face, the volume flux is computed from the adjacent velocity. The transported velocities are obtained using the rule of skew-symmetric approximation (weighting by $\frac{1}{2}$) where we use the velocity of the directly neighbouring cell for a pink or black face and the velocity from the diagonally neighbouring cell for a green or grey face.

The convective fluxes over the wall faces (red) are zero due to the impermeability condition at the wall.

Table 5 Convective fluxes at a T-face area for the momentum cell $u_{i,j,k}$

1.2 Diffusive Fluxes

The viscous fluxes at immersed boundaries (red faces) are non-zero due to the no-slip condition. The wall shear stress at the red cell faces is computed according to Meyer et al. (2010). First the velocity vector is interpolated to an interpolation point at a specified wall distance $\varDelta n = \varDelta x$ from the red face. For this a weighted linear least squares method is applied to neighbouring velocities. Then the wall shear stress is computed from an assumed velocity profile for the wall-parallel velocity component (here we use the Werner-Wengle wall function).

The geometry of the yellow face is generally quite complicated. Therefore, the viscous flux across a yellow face is approximated as (i) An increased viscous flux across the orange face and (ii) Half of the wall shear flux. The enlarged orange face has the area $\frac{1}{2}(S^{u}_{i-1,j,k} + S^{u}_{i,j,k})$ instead of $\min (S^{u}_{i-1,j,k}, S^{u}_{i,j,k})$ for the original orange face. Due to these approximations, the viscous flux across a yellow face-area is not evaluated separately, but in combination with the fluxes across the orange and red face areas, respectively.

The combined viscous flux due to the orange and yellow parts of a w- or e-face is formulated as

$$\begin{aligned} \left. (\mathbf {D}\mathbf {u}_h)_{i,j,k}\right| _{\mathrm {orange}+\mathrm {yellow}} = 2 \nu _{i,j,k} \left[ \frac{u_{i,j,k}-u_{i-1,j,k}}{{\varDelta x}_i} + (C^{u}_x )_{i,j,k} \right] \frac{1}{2}\left( S^{u}_{i-1,j,k} + S^{u}_{i,j,k}\right) \end{aligned}$$

(10)

where $(C^{u}_x )_{i,j,k}$ is a correction term for the non-orthogonal alignment between the cell face and the line connecting the velocities. The effective viscosity $\nu$, which is the sum of molecular and SGS viscosity, is interpolated for the whole face area using a harmonic average (Patankar 1980; Werner 1991; Peller 2010).

Fig. 22

Formulation of the diffusive fluxes for the red, yellow and orange cell-faces

The combined viscous flux due to the red and yellow faces is formulated as

$$\begin{aligned} \begin{aligned} \left. (\mathbf {D}\mathbf {u}_h)_{i,j,k}\right| _{\mathrm {red}+\mathrm {yellow}}&= -\frac{1}{2} \big [1 + \underbrace{(1-B^{u}_{i-1,j,k})}_{\text {cell-to-boundary}}\big ] B^{u}_{i,j,k} S^{\mathrm {IB}}_{i,j,k} (\tau _{\mathrm {w}x})_{i,j,k} \\&- \frac{1}{2} \big [1 + \underbrace{(1-B^{u}_{i+1,j,k})}_{\text {cell-to-boundary}}\big ] B^{u}_{i,j,k} S^{\mathrm {IB}}_{i+1,j,k} (\tau _{\mathrm {w}x})_{i+1,j,k} \end{aligned} \end{aligned}$$

(11)

with the wall shear stress $\tau _{\mathrm {w}x}$ and the cut wall area inside the cell (Hartmann et al. 2011)

$$\begin{aligned} S^{\mathrm {IB}}_{i,j,k} = \sqrt{(S^u_{i,j,k}-S^u_{i-1,j,k})^2 + (S^v_{i,j,k}-S^v_{i,j-1,k})^2 + (S^w_{i,j,k}-S^w_{i,j,k-1})^2} \end{aligned}$$

(12)

which is zero if a cell does not contain an immersed boundary. Summing Eq. (11) for all cells with $B^u_{i,j,k}=1$, the entire friction force $-\sum _{i,j,k} S^{\mathrm {IB}}_{i,j,k} (\tau _{\mathrm {w}x})_{i,j,k}$ is recovered.

Figure 22 displays two configurations of a u-momentum cell. For the first case (Fig. 22a), equation (11) assigns half the wall shear flux of each adjacent pressure cell to the depicted momentum cell whereas in the second case (Fig. 22b) the entire wall shear flux of the left and half the wall shear flux of the right cell is assigned to the depicted momentum cell.

The viscous fluxes across the tangential (n, s, t and b) cell faces are calculated according to

$$\begin{aligned} \begin{aligned} \left. \left( \mathbf {D}\mathbf {u}_h\right) _{i,j,k}\right| _{\textsc {t}}&= \sum _{\mathrm {colour}} \frac{1}{2} S^w_{i,j,k} \nu ^\mathrm {colour}_{i,j,k} \left[ \left( \frac{\varDelta u}{\varDelta z}+ C^{u}_z\right) ^\mathrm {colour}_{i,j,k} + \left( \frac{\varDelta w}{\varDelta x}+ C^{w}_x\right) ^\mathrm {colour}_{i,j,k}\right] \\&+ \sum _{\mathrm {colour}} \frac{1}{2} S^w_{i+1,j,k} \nu ^\mathrm {colour}_{i+1,j,k} \left[ \left( \frac{\varDelta u}{\varDelta z}+ C^{u}_z\right) ^\mathrm {colour}_{i+1,j,k} + \left( \frac{\varDelta w}{\varDelta x}+ C^{w}_x\right) ^\mathrm {colour}_{i+1,j,k}\right] \end{aligned} \end{aligned}$$

(13)

The formulation of the effective viscosity and the velocity gradient for the half-faces can be found in Tables 6 and 7. For pink, black and grey faces, the viscosity and the velocity gradients are constant over the entire face of the momentum cell.

The non-orthogonality correction terms in Eqs. (10) and (13) are formulated according to Kirkpatrick et al. 2003. The velocity gradient values for the correction terms are computed from the assumed velocity profile that is used in the calculation of the wall shear stress at a wall distance $\varDelta n = \varDelta x$.

Table 6 Effective viscosity at a T-face area for the momentum cell $u_{i,j,k}$

Table 7 Velocity gradients at a T-face area for the momentum cell $u_{i,j,k}$