Revision as of 19:10, 8 December 2009

2D Periodic Hill

Underlying Flow Regime 3-30

Test Case

Brief Description of the Test Case Studied

The measures of the geometry introduced by Mellen et al. (2000) relate to the hill height h. The hill constricts the channel height of 3.036h by about one third, whereas the inter hill distance is 9h. The contour of the 3.857h long two-dimensional hill is described by the following six polynomials.

$x\;\epsilon \;\left[0;\;0.3214h\right];\quad y(x)=\min(1;1+0\cdot x+2.420\cdot 10^{-4}\cdot x^{2}-7.588\cdot 10^{-5}\cdot x^{3});$ $x\;\epsilon \;\left[0.3215h;\;0.5h\right];\quad y(x)=0.8955+3.484\cdot 10^{-2}\cdot x-3.629\cdot 10^{-3}\cdot x^{2}+6.749\cdot 10^{-5}\cdot x^{3};$ $x\;\epsilon \;\left[0.5h;\;0.7143h\right];\quad y(x)=0.9213+2.931\cdot 10^{-2}\cdot x-3.234\cdot 10^{-3}\cdot x^{2}+5.809\cdot 10^{-5}\cdot x^{3};$ $x\;\epsilon \;\left[0.7144h;\;1.071h\right];\quad y(x)=1.445-4.927\cdot 10^{-2}\cdot x+6.950\cdot 10^{-4}\cdot x^{2}-7.394\cdot 10^{-6}\cdot x^{3};$ $x\;\epsilon \;\left[1.071h;\;1.429h\right];\quad y(x)=0.6401+3.123\cdot 10^{-2}\cdot x-1.988\cdot 10^{-3}\cdot x^{2}+2.242\cdot 10^{-5}\cdot x^{3};$ $x\;\epsilon \;\left[1.429h;\;1.929h\right];\quad y(x)=\max(0;2.0139-7.180\cdot 10^{-2}\cdot x+5.875\cdot 10^{-4}\cdot x^{2}+9.553\cdot 10^{-7}\cdot x^{3});$

At x/h=0 the hill height is maximal, whereas the boundary is flat in the range between x/h=1.929 and x/h=7.071. Between x/h=7.071 and x/h=9 the contour follows the above equations but the hill geometry is mirrored at x/h = 4.5. Besides the geometry, the following figure shows streamlines at a Reynolds number, that is based on the hill height h and the bulk velocity above the crest, of 5,600 [Rapp 2009)].

$Re={\frac {u_{b}h}{\nu }}$

$u_{b}={\frac {1}{2.035}}\int _{h}^{3.035h}u\left(y\right)dy$

The mean flow separates at the curved hill crown. In the wake of the hill the fluid recirculates before it attaches naturally at about x/h=4.5.

Experimental Setup

A water channel has been set up in the Laboratory for Hydromechanics of the Technische Universität München to investigate the flow experimentally. In total 10 hills with a height of 50 mm were built into the rectangular channel to accomplish periodicity whilst the measurement range lies between hills seven and eight. To achieve homogeneity in the spanwise direction an extent of 18 hill heights was appointed. The following figure sketches the experimental setup.

The 2D PIV measurements were undertaken between hills seven and eight - and to investigate the periodicity of the flow - between the hill pair six and seven through vertical laser light sheets. The homogeneity in the spanwise direction was controlled by 2D PIV measurements in horizontal planes. The PIV field data was thoroughly validated through 1D LDA measurements. Experiments were done at four Reynolds numbers: Re=5,600; Re=10,600; Re=19,000 and Re=37,000.

CFD Methods

The numerical part of the present study relies on two completely independent codes based on either curvilinear body-fitted grids with a colocated variable arrangement or Cartesian non-uniform grids using a staggered configuration. The objective is to present highly reliable results obtained by carefully cross-checking between the outcome of both numerical schemes and additional experimental data. Afterwards the investigations concentrate on the physical aspects of the flow considered.

In the following, both codes are described briefly.

Finite-volume code LESOCC

LESOCC solves the (filtered) Navier-Stokes equations based on a three-dimensional finite-volume method for arbitrary non-orthogonal and non-staggered block-structured grids (see, e.g. Fig. ??). The spatial discretization of all fluxes is based on central differences of second-order accuracy. Time advancement is performed by a predictor--corrector scheme. A low-storage multi-stage Runge-Kutta method (three sub-steps, second-order accuracy) is applied for integrating the momentum equations in the predictor step. Within the corrector step the Poisson equation for the pressure correction is solved implicitly by the incomplete LU decomposition method. Explicit time marching works well for DNS and LES with small time steps which are necessary to resolve turbulence motion in time. In order to ensure the coupling of pressure and velocity fields on non-staggered grids, the momentum interpolation technique is used. For modeling the non-resolvable subgrid scales, two different models are implemented, namely the well-known Smagorinsky model (1963) with Van Driest damping near solid walls and the dynamic approach with a Smagorinsky base model proposed by Germano et al. (1991) and modified by Lilly (1992). In order to stabilize the dynamic model, averaging of the numerator and the denominator in the relation for the determination of the Smagorinsky value was carried out in the spanwise homogeneous direction and also in time using a recursive digital low-pass filter (Breuer and Rodi 1996, Breuer 2002). The code and the implemented SGS models were validated on a variety of different test cases. For more information on this issue, please refer to Breuer and Rodi (1996), Breuer (1998, 2000, 2002).

Finite-volume code MGLET

MGLET is based on a finite-volume formulation for non-uniform Cartesian grids with a staggered arrangement of the spatially filtered variables (see, e.g.\ Fig. ??). The spatial discretization of the convective and diffusive fluxes is based on second-order central differences. The momentum equations are advanced in time by a fractional time stepping using either an explicit second-order central leapfrog scheme or a third-order Runge-Kutta scheme. For the solution of the Poisson equation for the pressure the ``Strongly Implicit Procedure (SIP) is implemented. For the representation of the hill geometry in the Cartesian grid an immersed boundary technique is used. All Cartesian cells lying inside the body are excluded from the computation. The excluded cells are determined by the intersection of the hill geometry with the Cartesian cells. The geometry of the hills is represented by a triangle mesh. The immersed boundary technique provides a smooth representation of the body surface in the Cartesian mesh by using third-order least squares interpolation for the interface cells (Peller 2006). This method prevents instabilities which are present in high-order Lagrange interpolation schemes. The code is used for DNS and LES simulations. It has been shown by several authors that second-order accuracy can be sufficient for DNS of flows provided the grid resolution is sufficient (Manhart and Friedrich 2002, Peller et al. 2006).

Boundary conditions and simulation parameters

Since the grid resolution in the vicinity of the wall is sufficient to resolve the viscous sublayer, the no-slip and impermeability boundary condition is used at the wall in both codes. The flow is assumed to be periodic in the streamwise direction and thus periodic boundary conditions are applied. Similar to the turbulent plane channel flow case the non-periodic behavior of the pressure distribution can be accounted for by adding the mean pressure gradient as a source term to the momentum equation in streamwise direction. Two alternatives exist. Either the pressure gradient is fixed which might lead to an unintentional mass flux in the configuration or the mass flux is kept constant which requires to adjust the mean pressure gradient in time. Since a fixed Reynolds number can only be guaranteed by a fixed mass flux, the second option is chosen here.

Furthermore, the flow is assumed to be homogeneous in spanwise direction and periodic boundary conditions are applied, too. For that purpose the use of an adequate domain size in the spanwise direction is of major importance in order to obtain reliable and physically reasonable results. To assure this criterion the two-point correlations in the spanwise direction have to vanish in the half-width of the domain size chosen. Based on the investigations by Mellen et al. (2000) a spanwise extension of the computational domain of $L_{z}=4.5h$ is used in all computations presented. It represents a well-balanced compromise between spanwise extension and spanwise resolution.

The Table summarizes the most important parameters of the simulations available. $N_{tot}$ denotes the total number of grid points used; the corresponding number of control volumes is slightly lower. Although a direct comparison of the number of grid points used by LESOCC and MGLET in one x-y plane is not reasonable, at least the number of points equidistantly distributed in the spanwise direction, $N_{span}$ , can be compared.

The dimensionless time-step size $\Delta t$ is also tabulated in Table~\ref{tab:simtab}, where the time is normalized by the ratio of the hill height $h$ and the bulk velocity ${U_{B}}$ taken at the crest of the hill. To reduce statistical errors due to insufficient sampling to a reasonable minimum, the flow field was averaged in spanwise direction and in time over a long period of $\Delta T_{avg}$ which is also given in Table~\ref{tab:simtab}. Partially $\Delta T_{avg}$ covers a time interval of about ${140}$ $ flow-through times.

Resolution issues

A detailed discussion about this concern is provided in Breuer et al. (2009). Here only a few issues should be mentioned. We start with the curvilinear grid design for the wall-resolved LES prediction at the highest Reynolds number chosen, i.e. Re = 10,595.

Curvilinear grids

For the simulation at Re = 10,595 LESOCC applies a curvilinear block-structured grid consisting of $N_{tot}=13.1$ million grid points corresponding to a total of about 12.4 million control volumes. The grid points are clustered in the vicinity of the lower wall, the upper wall, and the region where the free shear layer appears. Besides classical quality criteria such as orthogonality and smoothness, two main issues motivated the distribution of the grid points in space. These are the resolution of the near-wall region and of the inner domain.

To evaluate the first concern, the most important quality criterion is the distribution of non-dimensional $y^+$ values defined by $y^+ = \Delta y_{cc} \, u_{\tau} / \nu$ where $ \Delta y_{cc}$ denotes the distance of the cell center from the wall and $u_{\tau} = \sqrt{\tau_w

 / \rho}$ describes the shear stress velocity. Note that due to the

cell-centered variable arrangement $\Delta y_{cc}$ is half of the corresponding cell height $\Delta y$. Figure~\ref{fig:yplus} depicts the $y^+$ distribution along the lower wall at all Reynolds numbers considered. At $Re = 10,595$ the values are below $0.45$ with a mean value of about $0.2$ except at the windward side of the hill. Here the largest values of the wall shear stress are observed and the $y^+$ value reaches its maximum of about $1.2$. Hence the lower wall is well resolved. Regarding the wall-normal resolution the grid satisfies the requirements of a wall-resolved LES prediction. Compared to \cite{froehlich_05} who employed in their highly resolved simulations a curvilinear grid with about $4.6$ million CVs ($196 \times 128 \times 186$) especially the number of grid points in the wall-normal direction was increased to $220$ in the present investigation. Furthermore, the simulations resolve not only the lower wall (the hills) in more detail but also resolve the upper wall by a DNS-like representation ($y^+ \leq 0.95$ at $Re = 10,595$). Thus in contrast to \cite{froehlich_05} the application of wall functions is avoided. That allows to establish the influence of the resolution of the upper wall on the results. To prove the enhanced resolution, some numbers are provided. For instance, the cell sizes at the hill crest, which is a key region for the periodic hill flow are in the current case $\Delta x_{crest}/h = 0.026$ and $\Delta y_{crest}/h = 2.0 \times 10^{-3}$ whereas the corresponding values in \cite{froehlich_05} are $\Delta x_{crest}/h = 0.032$ and $\Delta y_{crest}/h = 3.3 \times 10^{-3}$, respectively. Owing to the increased resolution in streamwise and spanwise direction the cell sizes expressed in wall units are below $\Delta x^+ = 20$ and $\Delta z^+ = 9$ and thus lower than in \cite{froehlich_05} and substantially lower than the recommendations for wall-resolved LES given by Piomelli and Chasnov~\cite{piomelli96}. That also holds at the windward slope of the hill where the largest shear stresses are found.

%----------------------------------------------------------------------------- \begin{figure}

 \psfrag{x}[c][][2.5][0]{$x / h$}
 \psfrag{yplus}[c][][2.5][-90]{$y^{+}$}
 \begin{center}

%% \epsfig{file=./eps/yplus_all_re.eps,angle=-90,width=0.65\textwidth,clip=}

      \epsfig{file=./eps_col/yplus_all_re.eps,angle=-90,width=0.65\textwidth,clip=}
 \end{center}
 \vspace*{-0.5cm}
 \caption{Distribution of $y^+$ along the lower wall at different
          $Re$ using \lesocco.\label{fig:yplus}}

\end{figure} %-----------------------------------------------------------------------------

To evaluate the second concern, the resolution of the inner region, it is reasonable to estimate the size of the smallest scales given by the Kolmogorov length $\eta = (\nu^3 / \epsilon)^{1/4}$. In order to determine this quantity within the wall-resolved LES prediction, the dissipation tensor $\epsilon_{ij}$ was predicted and averaged in time and in spanwise direction. It has to be mentioned that this procedure represents only a rough estimate of $\epsilon_{ij}$ since on the one hand the diminutive SGS contribution is not included and on the other hand the present grid might slightly underestimate the dissipation tensor. Based on the dissipation rate $\epsilon = 1/2 \, \epsilon_{ii}$ and the kinematic viscosity $\nu$ of the fluid, $\eta$ can be determined and compared to the filter width $\Delta = (\Delta x \times \Delta y \times \Delta z)^{1/3}$ applied. In Fig.~\ref{fig:eta} profiles of $\Delta / \eta $ are shown at three different locations, $x/h = 0.5$, $1.0$ and $6.0$, respectively. With respect to the estimation given by Pope~\cite{pope2000} that the maximum dissipation takes place at a length scale of about $24 \, \eta$, these structures are resolved by at least 4-5 grid points at $x/h = 0.5$. At $x/h = 1.0$ and $6.0$ the maximum of $\Delta / \eta$ is about $6.3$ and thus at least about 4 times smaller than the decisive scales found by Pope~\cite{pope2000}. Consequently, the grid allows to resolve a substantial part of the dissipation. Overall the values $\Delta / \eta$ are smaller than in \cite{froehlich_05}. Especially in the vicinity of the upper wall the situation is strongly improved.

%----------------------------------------------------------------------------- \begin{figure}

 \psfrag{y/h}[r][][2.5][-90]{$y / h$}
 \psfrag{delta / eta}[c][][3.][0]{$\Delta / \eta$}
 \begin{center}
     \epsfig{file=./eps/prof_deltaeta_10595.eps,angle=-90,
                   width=0.6\textwidth,clip=}
 \end{center}
 \vspace*{-0.5cm}
 \caption{Profiles of $\Delta / \eta$ at three different vertical
          positions: $x / h = 0.5$, $1.0$, and $6.0$; curvilinear 
          grid used by \lesocc for the 
          wall-resolved LES prediction at $Re = 10,595$.
 \label{fig:eta}}

\end{figure} %-----------------------------------------------------------------------------

As explained above the generation of the grid for \lesocc was designed for the wall-resolved LES predictions at $Re = 10,595$. Regarding the lower Reynolds numbers considered, the grid was not modified when $Re$ was reduced. The reason is twofold. On the one hand a grid which is sufficiently fine for a certain Reynolds number should also be adequate for a lower $Re$. That is visible for example in Fig.~\ref{fig:yplus} which depicts the distributions of $y^+$ along the lower wall for all $Re$. The shear stress increases with decreasing Reynolds numbers, but for fixed $U_B$ and $h$ the viscosity also increases with decreasing $Re$. As a result the average and maximum $y^+$ values are strongly reduced. Consequently, with decreasing $Re$ the resolution becomes better and better. Thus by applying the same grid for the comparison of two $Re$, the effect of the grid is completely excluded from the consideration. On the other hand the intention was to perform DNS predictions for the lowest Reynolds numbers considered, i.e.\ at $Re = 700$, $1400$, and $2800$. Consequently, the grid should be sufficiently fine for a DNS at $Re = 2800$. That is exactly the Reynolds number in the classical plane channel flow predictions \cite{kim87,moser99} who applied a grid of about 2 million points. In the present case a six times finer grid is applied which accounts for the more complex flow field and the lower accuracy of the numerical method. Besides the criteria discussed in detail above, a further evidence of the adequacy of the grid for DNS at $Re \leq 2800$ is provided by two simulations carried out at $Re = 5600$ (see Table~\ref{tab:simtab}). One simulation was done as a wall-resolved LES using the dynamic SGS model (case \LDL) and the other was carried out without a subgrid-scale model (case \LDD). As will be discussed below, only marginal deviations were found between both cases. That is a clear hint that for the further reduced Reynolds numbers $Re \leq 2800$ the grid used delivers results which can be regarded as DNS.

In accordance with \cite{froehlich_05} the simulation at $Re = 10,595$ (case \LEL) was performed with the dynamic Smagorinsky model. Owing to the increased resolution the ratio of $\nu_t / \nu$ found in the present prediction is smaller than in the previous study. Moreover, by applying two different SGS models which delivered strongly differing eddy-viscosity values, Fr\"ohlich et al.~\cite{froehlich_05} have shown that the influence of this deviation on the LES prediction is low if a very fine grid is used as in the present case. That confirms that the present simulation is not materially inferior to a DNS near the walls.

Besides, for some additional simulations at very low Reynolds number between $100$ and $700$ to be discussed in Sect.~\ref{sect:instat}, a coarser grid with $164 \times 100 \times 64$ CVs was used.

Contributed by: (*) Christoph Rapp, (**) Michael Breuer — (*) Technische Universität München, (**) Helmut-Schmidt Universität Hamburg