## Abstract

The uniform flow over a nominally two-dimensional normal thin flat plate with blockage ratio 0.214 was numerically investigated in three dimensions by three methods: unsteady Reynolds-averaged Navier–Stokes (URANS) based on the realizable k–epsilon (RKE) turbulence model, URANS based on the k–omega shear stress transport (SST) turbulence model, and detached eddy simulation (DES). The Reynolds number based on the inlet flow velocity and the chord width of the plate was 117,000. A comprehensive comparison against earlier experimental results showed that URANS-SST method only could give a correct Strouhal number but overestimated the mean base pressure distribution and mean drag coefficient, while URANS-RKE and DES methods succeeded in giving accurate predictions of all. Moreover, by comparing the instantaneous vorticity contours and three-dimensional (3D) turbulent flow structures, it is found that DES is better suited for the present case because it can capture irregular small-scale structures and reproduce the three-dimensionality and low-frequency unsteadiness of the vortex shedding. Finally, through the volume-of-fluid (VOF)-based simulation of the free surface, it is demonstrated that the free surface has no significant effect on the mean drag coefficient and Strouhal number.

## 1 Introduction

In hydrodynamics, flow past a normal thin flat plate is a classic and fundamental topic. Besides, structures and components similar to flat plates are widely used in the field of offshore and marine engineering, such as heave damping plates of truss spar platforms, interceptors in high-speed crafts, and ship stabilizers. Despite its simple geometry, the wake is very complex since it is very sensitive to boundary conditions such as blockage ratios and end-plate conditions, especially at high Reynolds numbers (Re) due to the turbulence phenomena. Therefore, the detailed numerical study of this classic subject not only will enhance my understanding of this phenomenon including blockage effect and free surface effect but also can provide ground to verify which is the more precise and robust numerical approach considering contemporary computational methods.

In the present paper, the flow past a normal thin flat plate with a blockage ratio (*S/C*, where *S* is the reference area of a model and *C* is the cross-sectional area of a test-section) of 0.214 is studied numerically by three methods: unsteady Reynolds-averaged Navier–Stokes (URANS) in conjunction with the realizable *k*–epsilon model (URANS-RKE), URANS with Menter's *k*–*ω* shear stress transport (SST) model (URANS-SST) and detached eddy simulation based on Menter's *k*–*ω* SST model (DES-SST). Grid verification based on the grid convergence index (GCI) method is conducted. The results from URANS-RKE, URANS-SST, and DES-SST are compared and validated with data available in the literature. Whether it is possible to predict flow characteristics such as mean drag coefficient $(Cd\xaf)$ and Strouhal number (St) with reasonable accuracy by these methods is identified. Computational settings such as grid resolution, boundary conditions, and the size of the time-step needed for accurately capturing the correct flow field are explored. Besides, the free surface effect is studied numerically by the volume-of-fluid (VOF) method.

In Sec. 2, a brief introduction to the background of this problem is provided. In Sec. 3, a brief description of the mathematical formulation is given. In Sec. 4, computational details including the grid and numerical schemes used for the simulations are outlined. Section 5 contains the results and discussion. The free surface effect is explained in detail in Sec. 6. Finally, Sec. 7 summarizes the findings of the study.

## 2 Background

In the 20th century, the investigation of flow characteristics over bluff bodies and the effect of blockage was mainly concentrated upon experiments. One of the earliest comprehensive experimental studies was carried out by Fage and Johansen [1] who obtained detailed data such as pressure distribution, normal force coefficient, Strouhal number (St), and longitudinal spacing of vortices in the wake through systematic tests of the flow of air behind an inclined flat, sharp-edged plate of infinite span at 18 angles of incidence from 0.15 deg to 90 deg with a blockage ratio of 0.0714 at Re = 1.5 × 10^{5}. Additional experiments to determine the magnitude of the blockage effect showed that drag coefficients $(Cd\xaf)$ increase as *S/C*s increase, being 1.928, 2.05, and 2.144 corresponding to *S/C* = 0.0238, 0.0476, and 0.0952, respectively. Then, Maskell [2], El-sherbiny and Modi [3,4], Raju and Garde [5], Raju and Singh [6], Awbi [7,8], and Takeuchi and Okamoto [9] conducted a series of experiments for different bluff body shapes including two-dimensional (2D) circular cylinders, flat plates, and rectangular prisms of various aspect ratios with and without splitter plates. The Re in these experimental studies far exceeded 10^{3} and, namely, was limited to the subcritical Reynolds number regime wherein the pressure distribution and the drag coefficient are well known to be independent of Re for a normal plate in uniform flow. The *S/C* in all tests varied from 0.075 to 0.355. Takeuchi and Okamoto [9] concluded that the blockage effect of sidewalls can be ignored for *S/C* < 0.05 for a 2D flat plate.

At the same time, numerous semi-empirical methods for blockage correction on flow quantities such as $Cd\xaf$ and St were introduced based on the analysis of these experimental data. One of the well-known classical theories is the Maskell equation [2] based on an approximate relation describing the momentum balance and two empirical auxiliary relations describing invariance of the form of pressure distribution under constraint and governing wake distortion, respectively. Due to the success of Maskell's theory at small blockage ratios, a variety of modified Maskell equations, which extended its applicability to a wider range of bluff body shapes and flow conditions, were developed by other investigators such as El-sherbiny and Modi [3], Awbi [7], and Ramamurthy and Balachandar [10]. In addition to the aforementioned semi-empirical methods, other analytical methods, e.g., the potential flow model [3] and the purely empirical formula [5], were proposed.

With the rapid advances in computing technology over the past three decades, computational fluid dynamics (CFD) has become an important research tool in fluid mechanics. Regarding numerical simulations of flow over a normal flat plate, the most used methods include URANS, large eddy simulation (LES), detached eddy simulation (DES), and direct numerical simulation (DNS). Among them, DNS is currently only used in the range of about 80 < Re < 1200 due to the limitation of computer capacity. Najjar et al. [11–13] performed three-dimensional (3D) DNSs for two types of basic geometrical configurations, which were the separated-reattaching flow past a normal flat plate with an attached downstream splitter plate and the separated flow past a flat plate held normal to a uniform stream, on massively parallel processing computers at low Re, i.e., 250 and 1000. A low-frequency unsteadiness in the wake of a normal flat plate at Re = 250, in addition to the primary shedding frequency corresponding to Karman vortices, was captured in the time histories of instantaneous span-averaged drag and lift coefficients. The wake switched alternately between high and low drag regimes and reasonably better-organized streamwise vortices were seen in the high-drag regime compared with the low-drag regime. It was further shown that the flow at Re > 1000 is inherently 3D even if it is 2D extrinsically, so 2D flow simulations are untrustworthy. Then, Narasimhamurthy and Andersson [14], and Hemmati et al. [15] preformed 3D DNSs for the wake behind a 2D thin flat plate placed normally to the free stream at Re = 750 and 1200. The flow features such as the mean recirculation length, drag, pressure distributions and Reynolds-averaged statistics, were available in both cases. Furthermore, the hypothesis that there is a strong Re-effect only at Re < 1000 was confirmed.

Some LES and URANS numerical investigations in addition to Hemmati et al. [15] are listed below. Lasher [16] computed a 2D blocked flow (*S/C* = 0–0.33) normal to a flat plate at Re = 32,200 by both quasi-steady and transient methods, in combination with several versions of the *k*–*ɛ* turbulence model. They found that the quasi-steady computations provide an excellent prediction of the change in drag under blockage but underpredict the exact value of drag, while the transient simulations can predict the actual value of drag but are sensitive to the choice of the turbulence model. Tian et al. [17–20] also carried out a series of numerical studies on the flow past a normal flat plate at some high Reynolds numbers. First [17,19], the 2D flow normal to rectangular cylinders with different aspect ratios (*R* = 0.05–1) was numerically calculated using a *k–ω* SST turbulence model at Re = 21,400. The main conclusion was that the predicted hydrodynamic characteristics were in good agreement with published experimental data and numerical results for the rectangular cylinder with *R* > 0.6 while the 2D simulations were not able to give a satisfactory result for *R* < 0.6. The 3D LES/DES/DNS, which are known to capture 3D flow features well, were recommended for *R* < 0.6. Then [18,20], the 3D effect of the flow normal to a flat plate with the thickness ratio (thickness/height) of 0.02 at Re = 1.5 × 10^{5} was investigated using 2D URANS and 3D LES methods. The result showed that the 3D effect is significant at a high Reynolds number even though the plate is nominally 2D, and a low-frequency unsteadiness was found in the flow pattern as well, which both agree well with Najjar's study. All previous studies mentioned above are summarized by the authors in Table 1 along with the crucial characteristic parameters and methods. Among them, most numerical simulations did not consider blockage effect, and the study on the free surface effect is also lacking, so further numerical investigation of the flow over a flat plate with blockage and free surface effects is essential.

Study | Methods | Re | Blockage ratio (S/C) | Geometry | Thickness ratio | Reported results |
---|---|---|---|---|---|---|

Fage and Johansen (1927) [1] | Exp. | 1.5 × 10^{5} | 0.024–0.095 | Nominally 2D | 0.03 | $Cd\xaf$, St |

El-sherbiny (1972) [3] | Exp. | 10^{4}–1.5 × 10^{5} | 0.09–0.355 | Nominally 2D | 0.042–0.083 | $Cd\xaf$, St |

Raju and Singh (1975) [6] | Exp. | 1.64 × 10^{4} | 0.075–0.24 | Nominally 2D | 0–2.0 | $Cd\xaf$ |

Awbi (1978) [7] | Exp. | 9.2 × 10^{4} | 0.05–0.25 | Nominally 2D | 0.2–5 | $Cd\xaf$, St |

Takeuchi and Okamoto (1984) [9] | Exp. | 3.22 × 10^{4} | 0–0.3 | Nominally 2D | — | $Cd\xaf$, St |

Ringuette (2004) [21] | Exp. | 3000 | 0.052–0.093 | Nominally 2D/3D | 0.036–0.068 | $Cd\xaf$ |

Ramamurthy et al. (1989) [10] | Analysis | (0.2–3) × 10^{5} | 0.025–0.7 | Nominally 2D | – | $Cd\xaf$ |

Lasher (2001) [16] | URANS | 3.22 × 10^{4} | 0–0.33 | 2D | – | $Cd\xaf$, St |

Mannini (2015) [22] | URANS | 2.64 × 10^{4} | 0 | 2D/3D | 0.2 | $Cd\xaf$, St |

Tian et al. (2011) [17] | URANS | 2.14 × 10^{4} | 0.05 | 2D | 0.05–1.0 | $Cd\xaf$, St |

Tian et al. (2013) [19] | URANS | 2.14 × 10^{4} | 0.05 | 2D | 0.05–1.0 | $Cd\xaf$, St |

Najjar (1994) [11] | DNS | 1000 | 0.066 | 3D | 0 | $Cd\xaf$, St |

Najjar and Balachandar (1998) [13] | DNS | 250 | 0.0625 | 2D | 0 | $Cd\xaf$, St |

Narasimhamurthy et al. (2009) [14] | DNS | 750 | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Tian et al. (2012) [18] | LES | 1.5 × 10^{5} | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Tian et al. (2014) [20] | LES | 1.5 × 10^{5} | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Hemmati et al. (2018) [15] | DNS/LES | 1200 | 0.0625 | 3D | 0 | $Cd\xaf$, St |

Study | Methods | Re | Blockage ratio (S/C) | Geometry | Thickness ratio | Reported results |
---|---|---|---|---|---|---|

Fage and Johansen (1927) [1] | Exp. | 1.5 × 10^{5} | 0.024–0.095 | Nominally 2D | 0.03 | $Cd\xaf$, St |

El-sherbiny (1972) [3] | Exp. | 10^{4}–1.5 × 10^{5} | 0.09–0.355 | Nominally 2D | 0.042–0.083 | $Cd\xaf$, St |

Raju and Singh (1975) [6] | Exp. | 1.64 × 10^{4} | 0.075–0.24 | Nominally 2D | 0–2.0 | $Cd\xaf$ |

Awbi (1978) [7] | Exp. | 9.2 × 10^{4} | 0.05–0.25 | Nominally 2D | 0.2–5 | $Cd\xaf$, St |

Takeuchi and Okamoto (1984) [9] | Exp. | 3.22 × 10^{4} | 0–0.3 | Nominally 2D | — | $Cd\xaf$, St |

Ringuette (2004) [21] | Exp. | 3000 | 0.052–0.093 | Nominally 2D/3D | 0.036–0.068 | $Cd\xaf$ |

Ramamurthy et al. (1989) [10] | Analysis | (0.2–3) × 10^{5} | 0.025–0.7 | Nominally 2D | – | $Cd\xaf$ |

Lasher (2001) [16] | URANS | 3.22 × 10^{4} | 0–0.33 | 2D | – | $Cd\xaf$, St |

Mannini (2015) [22] | URANS | 2.64 × 10^{4} | 0 | 2D/3D | 0.2 | $Cd\xaf$, St |

Tian et al. (2011) [17] | URANS | 2.14 × 10^{4} | 0.05 | 2D | 0.05–1.0 | $Cd\xaf$, St |

Tian et al. (2013) [19] | URANS | 2.14 × 10^{4} | 0.05 | 2D | 0.05–1.0 | $Cd\xaf$, St |

Najjar (1994) [11] | DNS | 1000 | 0.066 | 3D | 0 | $Cd\xaf$, St |

Najjar and Balachandar (1998) [13] | DNS | 250 | 0.0625 | 2D | 0 | $Cd\xaf$, St |

Narasimhamurthy et al. (2009) [14] | DNS | 750 | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Tian et al. (2012) [18] | LES | 1.5 × 10^{5} | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Tian et al. (2014) [20] | LES | 1.5 × 10^{5} | 0.0625 | 3D | 0.02 | $Cd\xaf$, St |

Hemmati et al. (2018) [15] | DNS/LES | 1200 | 0.0625 | 3D | 0 | $Cd\xaf$, St |

## 3 Mathematical Formulation

### 3.1 Unsteady Reynolds-Averaged Navier–Stokes.

*x*

_{i}is the spatial coordinate,

*t*is the time, $u\xafi$ is the mean velocity, $ui\u2032$ is the fluctuating velocity, $P\xaf$ is the mean pressure,

*ρ*is the fluid density,

*μ*is the fluid viscosity, and $\tau ij=\u2212ui\u2032uj\u2032\xaf$ is Reynolds stress tensor. The Reynolds stresses are additional unknowns introduced by the averaging procedure, which are modeled using a turbulent viscosity

*ν*

_{t}based on the Boussinesq hypothesis.

*k*is the turbulent kinetic energy.

#### 3.1.1 Realizable k–Epsilon Model.

Compared with the standard *k*–epsilon model, the realizable *k*–epsilon model uses a new formulation for the dissipation rate *ɛ*. In the realizable *k*–epsilon model, transport equations of turbulent kinetic energy *k* and its dissipation rate *ɛ* are as follows:

*ν*is the molecular kinematic viscosity rate, and $Pk=\u2212ui\u2032uj\u2032\xaf\u2202u\xafi/\u2202xj$ is the rate of production of turbulent kinetic energy.

*ω*

_{k}.

The model constants specified in Openfoam are: *C*_{2} = 1.9, *σ*_{k} = 1.0, and *σ*_{ɛ} = 1.2. For further details of the realizable *k*–epsilon model, please see the study by Shih et al. [23]

#### 3.1.2 Menter's k–ω Shear Stress Transport Model.

*k–ω*model and the

*k–*epsilon model such that the

*k–ω*formulation is used in the near-wall region and switches to the

*k–*epsilon formulation in the free shear flow. The equations for the turbulence kinetic energy

*k*and the turbulence specific dissipation rate

*ω*are

*γ*,

*β*,

*σ*

_{k}, and $\sigma \omega $ are blended between the corresponding constants

*ϕ*

_{1},

*ϕ*

_{2}of the

*k–*epsilon and the

*k–ω*model via

*d*is the distance to the nearest wall,

*CD*

_{kω}= max(2

*ρσ*

_{ω2}1/

*ω*∂

*k*/∂

*x*

_{j}∂

*ω*/∂

*x*

_{j}, 10

^{−10}).

The constants of the SST model are *γ*_{1} = 5/9, *γ*_{2} = 0.44, *σ*_{k1} = 0.85, *σ*_{ω1} = 0.5, *β*_{1} = 0.075, *σ*_{k2} = 1.0, *σ*_{ω2} = 0.856, *a*_{1} = 0.31, *β*_{2} = 0.0828, and $\beta *=0.09$. For further details of the SST model, please see the study by Menter et al. [24]

### 3.2 Detached Eddy Simulation.

*F*

_{DES}is included in the dissipation term of the

*k*equation as follows:

*x*

_{1}, Δ

*x*

_{2}, Δ

*x*

_{3}) is the largest dimension of the local grid cell, and

*C*

_{DES}= 0.6.

## 4 Numerical Methods

The flow was assumed as 3D, unsteady, incompressible, viscous, and turbulent. The open-source CFD code Openfoam based on the finite-volume method (FVM) was used here for solving all corresponding Navier–Stokes equations in parallel. The PIMPLE (PISO-SIMPLE merged) scheme was used in the present study. Pressure Implicit Splitting of Operator (PISO) is an algorithm for time-dependent flows while SIMPLE stands for semi-implicit method for pressure linked equations and is used for steady-state problems. All spatial schemes and time schemes are in the second order. In all simulations, an adjustable time-step was used so that the maximum Courant–Friedrichs–Lewy (CFL) number did not exceed a pre-set value (1 for URANS, 0.5 for DES, respectively). The residual of 10^{−8} was used as the convergence criterion for all parameters such as pressure and velocity at each time-step.

The computational domain was defined in a 3D Cartesian coordinate system as shown in Fig. 1. A fixed thin flat plate with a thickness of *t* = 0.005 m and a height of *h* = 0.3 m is centered at the origin, normal to a uniform incoming flow with velocity *U*_{0} = 0.39 m/s. The distances from the plate center to the inlet and outlet boundaries are 4 m and 12 m, which are large enough to eliminate the far-field effects from the boundaries. The top and bottom boundaries are positioned at *y* = +0.7 m and −0.7 m and the spanwise length is 0.5 m.

The boundary conditions used for the simulations were set as follows:

- At the inlet boundary, a uniform flow profile was prescribed:
*u*=*U*_{0}= 0.39 m/s and*v*=*w*= 0; and the pressure was set as a Neumann boundary condition. The turbulence kinetic energy at the inlet was estimated bywhere(18)$k=32(I|uref|)2$*I*is the turbulence intensity and equal to 10% here, and the reference velocity*u*= 0.39 m/s. The turbulence dissipation rate_{ref}*ɛ*and specific dissipation rate*ω*were calculated aswhere $C\mu $ is a constant equal to 0.09, and the reference length scale used was(19)$\epsilon =C\mu 0.75k1.5L\omega =k0.5C\mu 0.25L$*L*= 0.07*H*. At the outlet boundary, Neumann boundary conditions were used for

*u*,*v*,*w*,*k*,*ω*,*ɛ*and the dynamic pressure was set to zero.On the flat plate, no-slip boundary condition was prescribed, i.e.,

*u*=*v*=*w*= 0;*p*was set as zero normal gradient;*k*was fixed at 10^{−10}.Free-slip boundary conditions were applied on the sidewalls, and symmetry boundary conditions were used for all flow variables on chord-wise planes.

The unstructured polyhedral mesh was constructed using starccm+ meshing software and then converted to Openfoam format by running the ccm26ToFoam utility. In the *XY*-plane, the computational domain was divided into three parts so that the grid resolutions could be controlled separately. All the cases in this study use the same topology. Furthermore, a very fine-grid resolution layer near the plate surface was employed to ensure that the maximum dimensionless height of the first layer $y+=yu\tau /\nu =y\tau w/\rho /\nu $, where *y* is the absolute distance from the wall, *u*_{τ} is the wall friction velocity, *τ*_{w} is the wall shear stress, *ρ* is the density of the fluid, and *ν* is the kinematic viscosity, was less than 1. In the *z*-direction, a uniform grid distribution was employed. Figure 2 shows the mesh of case A3.

All simulations were carried out on the 252 nodes SGI ICE-X computer cluster at the Advanced High-Performance Computing Nucleus (NACAD). Each node comprised two twelve-core 2.3-GHz Intel Xeon E5-2670v3 processors with 64-GB memory. The statistical average calculation of hydrodynamic quantities was performed after the periodic vortex shedding began. In all cases, the sampling time 100 *h/U _{0}* (approximately 20 shedding cycles) was used according to previous studies [14,18–20], which was found to be long enough for the results to converge to a statistically stationary state.

## 5 Results and Discussion

### 5.1 Verification and Validation Studies.

For the verification of URANS simulation, the GCI method derived by Roache [25], which has been evaluated over numerous CFD cases, is widely accepted and recommended. However, for LES and DES, the verification procedure is still an area of active research and there is not a uniform and widely used method until now due to their inherent feature that the resolved length scale is smaller as the mesh is refined. In this study, for the convenience of research, the GCI method was uniformly adopted for all cases using the mean drag coefficient (i.e., $Cd\xaf=Fd/0.5\rho U2S$) as a verification variable.

The hydrodynamic results of the cases with different grid resolutions and simulation methods are listed in Table 2, where CFL represents Courant number. $Cd\xaf$ was calculated from a non-dimensional time of 77*U _{0}/h*, as shown in Figs. 4(a), 5(a), and 6(a). St was estimated from the Fast Fourier Transform (FFT) of the instantaneous cross-stream velocity at points A and B in Fig. 1, which are located at

*h*and 2

*h*downstream along the

*x*-axis, respectively. Taking case A4 for example, the time series of

*y*-component velocity at point A is shown in Fig. 4(c), and the related spectrum is shown in Fig. 4(d). The dominant peak in the spectrum corresponds to St and was the same at two monitoring points. Cases A1, B1, and C1 were used to verify the influence of the number of grids in the spanwise direction on the results, i.e., $Cd\xaf$ and St. It was found that the grid resolution in the spanwise direction has a small influence on the mean drag coefficient for all simulation methods but almost has no effect on St. Cases C4 and C6 are to explore the influence of different time schemes on DES result.

Cases | Methods | Time scheme | Grids in XY-plane | Grids in z-direction | Elements | CFL | $Cd\xaf$ | St |
---|---|---|---|---|---|---|---|---|

A1 | URANS-RKE | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <1 | 3.3 | – |

A2 | CrankNicolson 0.9 | 109,546 | 32 | 3,505,472 | <1 | 3.1759 | 0.2 | |

A3 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <1 | 3.163 | 0.2 | |

A4 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <1 | 3.17 | 0.2 | |

B1 | URANS-SST | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <1 | 4.72 | 0.2 |

B2 | CrankNicolson 0.9 | 109,546 | 32 | 3,505,472 | <1 | 6.97 | 0.2 | |

B3 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <1 | 6.86 | 0.2 | |

B4 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <1 | 6.72 | 0.2 | |

C1 | DES-SST | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <0.5 | 2.95 | 0.19 |

C2 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <0.5 | 2.874 | 0.2 | |

C3 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <0.5 | 3.116 | 0.19 | |

C4 | Backward | 244,927 | 50 | 12,246,350 | <0.5 | 3.152 | 0.2 | |

C5 | CrankNicolson 0.9 | 388,396 | 63 | 24,468,948 | <0.5 | 3.118 | 0.2 | |

C6 | Backward | 388,396 | 63 | 24,468,948 | <0.5 | 3.169 | 0.2 |

Cases | Methods | Time scheme | Grids in XY-plane | Grids in z-direction | Elements | CFL | $Cd\xaf$ | St |
---|---|---|---|---|---|---|---|---|

A1 | URANS-RKE | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <1 | 3.3 | – |

A2 | CrankNicolson 0.9 | 109,546 | 32 | 3,505,472 | <1 | 3.1759 | 0.2 | |

A3 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <1 | 3.163 | 0.2 | |

A4 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <1 | 3.17 | 0.2 | |

B1 | URANS-SST | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <1 | 4.72 | 0.2 |

B2 | CrankNicolson 0.9 | 109,546 | 32 | 3,505,472 | <1 | 6.97 | 0.2 | |

B3 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <1 | 6.86 | 0.2 | |

B4 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <1 | 6.72 | 0.2 | |

C1 | DES-SST | CrankNicolson 0.9 | 158,323 | 20 | 3,166,460 | <0.5 | 2.95 | 0.19 |

C2 | CrankNicolson 0.9 | 158,323 | 40 | 6,332,920 | <0.5 | 2.874 | 0.2 | |

C3 | CrankNicolson 0.9 | 244,927 | 50 | 12,246,350 | <0.5 | 3.116 | 0.19 | |

C4 | Backward | 244,927 | 50 | 12,246,350 | <0.5 | 3.152 | 0.2 | |

C5 | CrankNicolson 0.9 | 388,396 | 63 | 24,468,948 | <0.5 | 3.118 | 0.2 | |

C6 | Backward | 388,396 | 63 | 24,468,948 | <0.5 | 3.169 | 0.2 |

Table 3 presents the parameters of verification based on the data gathered in Table 2. St is not included in this verification study because the differences of St among all cases are very small and almost all of them equal to 0.2. The verification of $Cd\xaf$ for URANS-SST cases was also not conducted as their values were seriously incorrect and the simulation could not achieve a stable periodic state as shown in Fig. 5(a). Refinement ratio is *h*_{coarse}/*h*_{fine}, where *h* is the grid size. The convergence ratio is (*S*_{2} − *S*_{1})/(*S*_{3} − *S*_{2}) where the solutions for the fine, medium, and coarse grids are *S*_{1}, *S*_{2}, and *S*_{3}, respectively. *P* is the apparent order of convergence. Whether the solutions are in the asymptotic range of convergence (i.e., a converged solution is being approached asymptotically) can be checked with the relationship *GCI*_{23}/*r*^{p}*GCI*_{12} ≅ 1. Hence, according to Table 3, the numerical uncertainties in the fine-grid solution of the mean drag coefficient for all three cases were small and the grid independence was achieved.

Variable | Grids | Refinement ratio | Convergence ratio | P | Asymptotic solution | Grid uncertainty (GCI_{21}) % | Grid uncertainty (GCI_{32}) % | GCI_{32}/r^{p}GCI_{21} |
---|---|---|---|---|---|---|---|---|

$Cd\xaf$ | A234^{a} | 2^{1}/^{3} | −0.512 | 2.9 | 3.177 | 0.27 | 0.53 | 1.002 |

C235^{b} | 2^{1}/^{3} | 0.0083 | 20.8 | 3.118 | 0.0007 | 0.8 | 1.0006 | |

C246^{c} | 2^{1}/^{3} | 0.06 | 12.1 | 3.17 | 0.04 | 0.7 | 1.005 |

Variable | Grids | Refinement ratio | Convergence ratio | P | Asymptotic solution | Grid uncertainty (GCI_{21}) % | Grid uncertainty (GCI_{32}) % | GCI_{32}/r^{p}GCI_{21} |
---|---|---|---|---|---|---|---|---|

$Cd\xaf$ | A234^{a} | 2^{1}/^{3} | −0.512 | 2.9 | 3.177 | 0.27 | 0.53 | 1.002 |

C235^{b} | 2^{1}/^{3} | 0.0083 | 20.8 | 3.118 | 0.0007 | 0.8 | 1.0006 | |

C246^{c} | 2^{1}/^{3} | 0.06 | 12.1 | 3.17 | 0.04 | 0.7 | 1.005 |

A234 represents cases A2, A3, and A4.

C235 represents cases C2, C3, and C5.

C246 represents cases C2, C4, and C6.

In order to evaluate the performance of different numerical methods, the results in the fine-grid solution for all cases were compared with published experimental data, as shown in Table 4. The experimental data from El-sherbiny (1972) [3] was used as the benchmark data *D*. Validation was performed using both the simulation value *S* and the extrapolated value *Se*, as summarized in Table 5, where *E* is the comparison error and *Ee* is the corrected comparison error. The table shows that all the simulation results were in good agreement with the experimental data except $Cd\xaf$ of URANS-SST cases in which the relative error reached 116.8%.

Case/method | $Cd\xaf$ | St | Blockage ratio | Re | Thickness ratio |
---|---|---|---|---|---|

Present A4/RKE | 3.17 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present B4/SST | 6.72 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C5/DES | 3.118 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C6/DES | 3.169 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

El-sherbiny (1972)/Exp.^{a} | 3.1 | 0.193 | 0.205 | 1.2 × 10^{5} | 0.054 |

RANGA (1975)/Exp. | 2.9 | – | 0.185 | >10^{3} | 0 |

RANGA (1975)/corrected Exp. | 3.1 | – | 0.214 | >10^{3} | 0 |

AWBI (1978) /Exp. | 3.27 | 0.195 | 0.2 | 9.2 × 10^{4} | 0.2 |

Takeuchi (1984)/Exp. | 3.1 | – | 0.214 | 0.3 × 10^{5} | – |

Case/method | $Cd\xaf$ | St | Blockage ratio | Re | Thickness ratio |
---|---|---|---|---|---|

Present A4/RKE | 3.17 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present B4/SST | 6.72 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C5/DES | 3.118 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C6/DES | 3.169 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

El-sherbiny (1972)/Exp.^{a} | 3.1 | 0.193 | 0.205 | 1.2 × 10^{5} | 0.054 |

RANGA (1975)/Exp. | 2.9 | – | 0.185 | >10^{3} | 0 |

RANGA (1975)/corrected Exp. | 3.1 | – | 0.214 | >10^{3} | 0 |

AWBI (1978) /Exp. | 3.27 | 0.195 | 0.2 | 9.2 × 10^{4} | 0.2 |

Takeuchi (1984)/Exp. | 3.1 | – | 0.214 | 0.3 × 10^{5} | – |

This experimental result will be used as benchmark data in the validation study.

### 5.2 Hydrodynamic Forces and Frequency Analyses.

The distribution of the mean pressure coefficient along the front and back sides of the plate from the present simulations compared against the previous experimental studies is shown in Fig. 3. The mean pressure coefficient is defined as $Cp\xaf=2(P\xaf\u2212P0)/(\rho U02)$, where the reference pressure *P*_{0} is the inlet pressure, the reference velocity *U*_{0} is the free stream velocity, and the over-bars signify averaging in time and spanwise direction. The pressure distribution on the upstream side is almost the same regardless of the different numerical methods and matches well with the experimental results of El-sherbiny [3] and Awbi [7]. However, on the backside, the pressure predicted by URANS-SST is largely lower than those in the URANS-RKE and DES-SST simulations and the experiments, which are consistent with each other and nearly constant. As shown in Table 4, the values of $Cd\xaf$ calculated by 3D URANS-RKE and SST-DES are in good agreement with the experimental results, whereas the present URANS-SST overpredicts $Cd\xaf$ almost twice, compared with experimental results. Due to the above discrepancy, it can be concluded that the URANS-SST method failed to simulate the present case. URANS-SST maybe not applicable for simulating large separation flows behind bluff bodies especially at the high blockage. The possible cause will be analyzed and discussed later in Sec. 5.4.

The time histories of *C*_{d} of cases A4, B4, and C3 are plotted in Figs. 4(a), 5(a), and 6(a). It appears that *C*_{d} in 3D URANS simulations exhibits a regularly and periodically repeating feature while *C*_{d} in 3D DES displays an irregular fluctuation. A low-frequency unsteadiness is observed in 3D DES according to the spectrum of *C*_{d} in Fig. 6(b), which was also reported by many previous numerical studies [11,13,20]. The time series of the cross-stream velocity at point A in cases A4, B4, and C3 are plotted in Figs. 4(c), 5(c) and 6(c). Similarly, it appears that the instantaneous velocities in 3D URANS simulations exhibit a regularly and periodically repeating feature while that in 3D DES displays more random fluctuations. As shown in Figs. 4(d), 5(d), and 6(d), St, which corresponds to the dominant peak in the spectra, was nearly the same for cases A4, B4, and C3. In fact, according to Table 4, both 3D URANS and 3D DES can give reasonable predictions of St.

As presented in Fig. 6(d), the frequency power spectrum of the cross-stream velocity fluctuation conforms well to the −5/3 slop power law [26], which seems to indicate that the resolved scales reach an inertial subrange. In order to justify that the present DES simulations are accurate enough to capture Taylor microscale [27] flow variation, the grid size was compared with Taylor microscale $\lambda =10\nu k/\epsilon $, where *ɛ* is the energy dissipation rate defined as $\epsilon =2\nu SijSij\xaf$, *k* is the turbulent kinetic energy. Figure 7 shows the contour of the spanwise averaged ratio of the grid size $\Delta cell=\Delta x\Delta y\Delta z3$ to Taylor microscale, and Table 6 gives the average grid resolution relative to Taylor microscale along six different downstream positions, demonstrating that Δ_{cell}/*λ* is typically less than 2 in most of the region behind the plate except for a small part of the high shear area where the ratio is relatively large and close to 5.9, and the average ratios are less than 0.35 from h to 6 h downstream the plate for cases C3 and C6.

### 5.3 Time and Spanwise Averaged Streamlines and Velocity Profiles.

The streamlines of the time and spanwise averaged flow field in the near wake for cases A4, B4, and C6 are shown in Figs. 8(a)–8(c). The nearly symmetrical recirculation bubbles behind the plate are demonstrated for those three cases, confirming the statistical convergence of the present cases. The mean recirculation lengths, which are defined as the x-axis distance from the plate to the position in the centerline where the mean streamwise velocity changes from negative to positive, correspond to 1.92 *h*, 1.65 *h*, and 2.2 *h* for cases A4, B4, and C6, respectively. In addition, the positions of the cores of the bubbles can be determined, approximately located at X ≈ 0.83 *h*, 0.304 *h* and 1.04 *h* for cases A4, B4, and C6, respectively. The closer distance between the core and the plate in the URANS-SST case is believed to be directly related to the inconsistent low mean pressure distribution shown in Fig. 3 and high drag force.

The profiles of the mean streamwise velocity in the near wake at six different positions downstream of the plate for cases A4 and C6 are shown in Figs. 9(a) and 9(c). When the sections are within the recirculation zone, the streamwise velocity is negative near the centerline as shown in the velocity profiles at the location X = *h* in cases A4 and C6. In Figs. 9(b) and 9(d), the mean cross-stream velocity profiles perform anti-symmetric variations, which also coincides with earlier studies. The mean velocity distributions of cases A4 and C6 are basically same, except for the regions very close to the plate such as X = *h* and X = 2 *h*.

### 5.4 Spanwise Averaged Vorticity and Time and Spanwise Averaged Turbulent Viscosity Contours.

Figures 10–12 show spanwise averaged vorticity contours in one period for cases A4, B4, and C6, respectively. The observation of the figures indicates that URANS predicts regular von Karman vortex-shedding while DES reveals more small-scale disordered structures. By comparing Figs. 10 and 11, we can find that the vorticity behind the plate in case B4 is stronger than that in case A4.

Figure 13 compares the mean turbulent viscosity from the 3D URANS-RKE with that from 3D URANS-SST. It is noteworthy that the eddy viscosity predicted by URANS-SST is significantly lower than that predicted by URANS-RKE. This underestimation of the turbulent viscosity may directly result in reduced dissipation and enhanced convection of the flow velocity field behind the plat as shown in Fig. 11. Finally, a lower pressure field in the wake than the real value is produced, and an incorrect higher drag is estimated in consequence.

The SST *k–ω* turbulence model is a hybrid model combining the Wilcox *k*–*ω* and the *k*–epsilon models because the *k*–*ω* model is more suited for simulating flow in the viscous sublayer and the *k*–epsilon model is ideal for predicting flow behavior in regions away from the wall. A blending function F1, as shown in Eq. (11), activates the *k–ω* model near the wall and the *k–*epsilon model in the free-stream. Therefore, in theory, the calculation result of URANS-SST should be better than URANS-RKE or at least nearly the same. In order to figure out the reason why URANS-SST failed in this study, it is very necessary to study its formula specifically. In the SST model, the kinematic eddy viscosity limiter (*SF*_{2}), which aims to limit the increase of eddy viscosity, was introduced to help the *k–ω* model better predict flow separation in aerodynamics, as shown in Eq. (13). This item is very likely to be the direct cause of the underestimation of the turbulent viscosity for a URANS-SST simulation of a bluff body at a large Reynolds number like the present case.

### 5.5 Three-Dimensional Flow Structures.

*Q*criterion proposed by Hunt et al. [28] was used here to identify 3D vortex structures in the wake. The

*Q*formulation computes the second invariant of the velocity gradient tensor as

_{ij}is the vorticity tensor and

*S*

_{ij}is the strain rate tensor.

The flow structures resolved by DES-SST and URANS-RKE are shown in Fig. 14, in which Figs. 14(a)–14(d) present an example of the 3D and spanwise views of the instantaneous iso-surfaces of *Q* = 0.1 for cases C6 and A4, respectively. The primary coherent shedding vortices can be resolved by the two methods. However, only a series of relatively large parallel 2D vortices are reproduced in the result of URANS-RKE simulation while DES-SST captures highly irregular small-scale structures. Moreover, in DES, the 3D vortex tubes are observed to meander widely and their evolution processes including vortex stretching, twisting, and breaking are also clearly displayed. The reproduced irregular dynamical features of the flow structures by DES correspond to the low-frequency unsteadiness of force coefficients as explained in Sec. 5.2. Finally, the vortex around the plate looks like less complex, compared to pure LES [20]. This is because unsteady vortex motions in near-wall regions are not resolved in DES but calculated with the RANS model. If an accurate investigation of near-wall regions behavior is needed, pure LES will be more appropriate.

## 6 The Effect of Free Surface

In this paper, interFoam solver is used to model the free surface, which solves the Navier Stokes equations for two incompressible, isothermal immiscible fluids using the VOF method. Figure 15 shows the grids of cases with the free surface. Cases E1 and E2 are used to verify the influence of the number of grids in the spanwise direction on the results. Their grid distributions in the *XY*-plane are the same as case A3, but the total number of grids in the *z*-direction of case E2 is 60 which is double that of case E1.

The results including $Cd\xaf$ and St for cases E1 and E2 which simulated free surface were compared with those in cases C5 and C6 which did not consider free surface and published experimental data, as shown in Table 7. The time histories of *C*_{d} of cases E1 and E2 are plotted in Figs. 16(a) and 16(d). It appears that *C*_{d} in case E1 displays a bigger fluctuation than that in case E2. Besides, according to Table 7, the value of $Cd\xaf$ in case E1 is much larger than that in case E2, which are 3.32 and 2.89, respectively; however, both of them are very close to the values 3.118 and 3.169 in cases C5 and C6 and the experimental data 3.1. It is concluded that the free surface has a relatively small effect on $Cd\xaf$ because the relative differences of $Cd\xaf$ in cases E1 and E2 compared with the cases without free surface and the experiment are less than 8%, but the value of $Cd\xaf$ is kind of sensitive to the number of grids in the *z*-direction in the present study since the relative difference between case E1 and case E2 is approximately 15%. The time series of the cross-stream velocity at point B in cases E1 and E2 are plotted in Figs. 16(b) and 16(e). Similarly, it appears that the instantaneous velocity in case E1 exhibits a stronger fluctuation than that in case E2. As shown in Figs. 16(e) and 16(f), St, which corresponds to the dominant peak in the spectra, are 0.19 and 0.2 for cases E2 and E2, respectively. The values of St in cases E1 and E2 are nearly the same and in concert with those in cases C5 and C6 and the experimental data, indicating that the free surface and the number of grids in the *z*-direction has almost no effect on St.

Case/method | $Cd\xaf$ | St | Blockage ratio | Re | Thickness ratio |
---|---|---|---|---|---|

Present C5/DES | 3.118 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C6/DES | 3.169 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

El-sherbiny (1972)/Exp. | 3.1 | 0.193 | 0.205 | 1.2 × 10^{5} | 0.054 |

Present E1/DES and VOF | 3.32 | 0.19 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present E2/DES and VOF | 2.89 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Case/method | $Cd\xaf$ | St | Blockage ratio | Re | Thickness ratio |
---|---|---|---|---|---|

Present C5/DES | 3.118 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present C6/DES | 3.169 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

El-sherbiny (1972)/Exp. | 3.1 | 0.193 | 0.205 | 1.2 × 10^{5} | 0.054 |

Present E1/DES and VOF | 3.32 | 0.19 | 0.214 | 1.17 × 10^{5} | 0.017 |

Present E2/DES and VOF | 2.89 | 0.2 | 0.214 | 1.17 × 10^{5} | 0.017 |

Figure 17 shows some surface results in cases E1 and E2. It is observed that the stronger vortex just behind the plate leads to the larger height difference of free surface between the front and rear of the plate in Fig. 17(a), while there are no large vortexes just behind the plate and the height difference of free surface between the front and rear of the plate is relatively small in Fig. 17(b). However, overall, the fluctuation of the free surface in Fig. 17 is small compared to the height of the plate immersed into the water surface. Thereby, it is reasonably justified that the free surface effect on the statistical flow field characteristics including $Cd\xaf$ and St is small when the fluctuation of the free surface is small compared to the height of the plate immersed into the water surface. Figure 18 shows some details of the calculation results including the velocity contour and the pressure contour for cases E1 and E2.

## 7 Conclusions

This paper did a comprehensive CFD analysis including verification and validation on the problem of flow over a flat plate with three-dimensional, blockage, and free surface effects, using URANS and DES methods. Detailed results such as $Cd\xaf$, St, streamlines, and turbulent viscosity contour are presented and discussed. To sum up, the following concluding remarks can be made:

The URANS-RKE and DES-SST methods can predict the mean drag coefficient $(Cd\xaf)$ and St accurately, but URANS-SST overestimates $Cd\xaf$. The cores of the bubble of the streamlines in the URANS-SST simulation are closer to the plate than URANS-RKE and DES-SST simulations as shown in Fig. 8. In addition, the vorticity behind the plate in the URANS-SST simulation is significantly stronger than that in URANS-RKE and DES-SST simulations, as shown in Figs. 10–12. It is probably because the introduction of the kinematic eddy viscosity limiter in the SST model causes eddy viscosity to be underestimated. Then, this underestimation of the turbulent viscosity results in reduced dissipation and enhanced convection of the flow velocity field behind the plate. Further confirmation is still awaited.

The low-frequency unsteadiness of the drag force coefficient in the wake behind the plate was found in DES, which is in line with earlier numerical studies. Besides, the DES method can capture highly irregular small-scale structures and 3D transient flow behavior, but URANS-RKE only can reproduce the quasi-periodic regular shedding corresponding to St. Therefore, for the unsteady turbulence simulation of the separated flow behind a blunt body at large Reynolds numbers, the DES method is a more reliable and better-suited method.

In general, when the fluctuation of the free surface is small compared to the height of the plate immersed into the water, the free surface almost has no effect on the statistical flow field characteristics including $Cd\xaf$ and St.

## Acknowledgment

This research has been developed with the support of the Advanced High-Performance Computing Nucleus (NACAD) at COPPE, Federal University of Rio de Janeiro (UFRJ). This support is gratefully acknowledged. The authors would also like to express their thanks to CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), CNPq (Brazilian Research Agency) and LOC/COPPE/UFRJ (Laboratory of Waves and Current of COPPE, Federal University of Rio de Janeiro).

## Conflict of Interest

There are no conflicts of interest.

## References

*ɛ*Eddy Viscosity Model for High Reynolds Number Turbulent Flows