Simulation of Pressure–Velocity Correlations by Green’s Function Based on Reynolds Stress Model

Abstract

Cost-effective wind energy harvesting by wind turbines in urban areas needs to strengthen the required flow field properties, such as mean velocity, turbulence, and its distribution. This paper conducts a series of CFD simulations to investigate the characteristics and related mechanisms of flow within the cavity, considering the force–turbulence interactions at the RANS scales. The pressure–velocity correlation term is formulated and solved by the elliptic relaxation equation to compensate for the Reynolds stress overestimation. Numerical simulations of flow over an open cavity with the proposed model are compared with corresponding PIV data. The results show that the mean velocity and the fluctuation velocity along the streamwise direction exist a slightly favorable pressure gradient. While the fluctuation velocity and fluctuation pressure show different correlation characteristics along the streamwise direction. Moreover, the pressure–velocity fluctuation correlation becomes obvious near the upper corner of the cavity due to the favorable pressure gradient. Hence, the leading and trailing locations of the cavity are both obvious favorable regions and further emphasis should be put on both high-accurate simulation methods and practical applications.

1. Introduction

Due to the impacts on pressure, velocity, and density fluctuation, flow over an open cavity has been investigated extensively in a variety of areas, such as pantograph cavities, door gaps, open-window-induced air leakage in transportation, flow-induced structural vibrations in civil engineering, and pipe networks in other energy-related side branches [1]. Because of its extensive presence and effects on the above areas, cavity flow has received considerable investigations over the past few decades [2]. Especially as material constraints in the turbine and complex aerodynamic solutions, it is imperative to study the interaction between the secondary flow with the main circulation [3]. Thus, it is required to reduce turbine energy losses by understanding of the interaction between the primary and the secondary flow [4].

Most open-cavity-flow studies have carried out flow and pressure data collection by visualization facilities [5]. Corresponding results have shown that a feedback mechanism of the shear layer–cavity corner having obvious impacts on the unsteady turbulent flow in the cavity is the intrinsic characteristic over a considerable quantity of Reynolds (Re) and Mach numbers [6]. This feedback mechanism generates a stable recirculation zone in an open cavity. Meanwhile, it is found that the unsteady pressure is fundamentally correlated with the vorticity distribution within the cavity [7]. Therefore, accurate assessment of the pressure distribution is indispensable in flow-induced vibrations [8]. Moreover, a deep understanding of the fluctuations leads to a better understanding of the vorticity dynamics around the cavity region [9]. It also should be noted that cavity flow of non-Newtonian materials can extend from mineral oil [10,11], melting material [12,13], painting [14], polymer [15], or temperature effects [16].

Even the non-intrusive system, such as PIV, can serve high-precision flow field data; pressure distribution is required to be calculated by instantaneous material acceleration based on the Poisson equation after finishing the PIV tests. It is not only time-consuming but also needs a high amount of computer power [17]. Recently, numerical methods for science and engineering applications were also adopted to describe the complex flow such as shear layer effects in an open cavity [18]. Compared to laboratory test technology, there are several methods to model a flow field by numerical simulation. They are different in complexity, accuracy, and cost of use and are more or less specifically derived for certain types of flows. The large-eddy simulation (LES) method directly simulates turbulent flows by solving a time-dependent Navier–Stokes (N-S) equation. This technique is a powerful tool, but its application is limited to systems with moderate Re numbers and simple geometry shapes because of the formidable computation requirements. The Reynolds-averaged Navier–Stokes (RANS) offers a more practical option to simulate flows with a high Re number that are always encountered in engineering practices. The RANS method, which solves the Reynolds-averaged transport equations, can resolve the mean motion accurately [19]. A turbulence model is adopted to represent the small-scale eddy flows [20,21]. This method is simple to use and appealing for industrial problems. However, it is impossible to build a universal RANS model that is able to cope with all kinds of engineering turbulent flow [22]. RANS models should be modified according to various different flows. It has been, however, reported that there are still several unsolved key issues in the simulation of flow over an open cavity, which are large Re number, shear layer flow and recirculation zones existing along the streamwise direction, strong streamline curvatures, and favorable and adverse pressure gradients, etc., [23,24]. From a perspective of turbulence modeling, the pressure–velocity correlation term in N-S equations are of significant relevance, which is responsible for energy transfer between mean flow and fluctuations [25]. There still exists the closure problem in the Reynolds stress function, which plays an important role in cavity flow simulation [26].

It has been proved that the Green’s function can be reclaimed by cross-correlated fluctuations at two locations for fluid and waves [27]. The extraction of the Green’s function from fluctuations has recently received so much attention and been applied to acoustics [28], seismology [29], and structural engineering [30]. More importantly, the capability of the Green’s function in handling very thin boundary layers has been emphasized [31], such as channel flow, Couette flow, and pipe flow, which are characterized by shear flows at high Reynolds numbers. Those effects can be captured exactly by the analytic form of the Green’s function [29].

With the aim to solve the above-mentioned key issues, it is imperative to calculate the mean flow and turbulence quantities accurately. The pressure–velocity correlation that exists in the Reynolds stress transport equations is associated with the volume integrals of the two-point correlation between the velocity and pressure fluctuations [31]. Most of the turbulence models of RANS methods do not consider the pressure–velocity correlation term (denoted as $D_{ij}$ in the following sections), which are still based on a quasi-homogeneous and local-equilibrium assumption [32]. Several models have been proposed in previous works for modeling the $D_{ij}$ term by a realizability condition or Caley–Hamilton theorem, which fulfill the closure requirement. Most of the models separate the $D_{ij}$ term into three parts, which are simulated individually [33]. However, these models cannot be combined with the wall boundary conditions, and thus fail to calculate the near-wall turbulent quantities correctly [34]. With the purpose of improving this inadequacy, an ad hoc wall echo model was adopted along the normal direction to damp the velocity fluctuations. However, this method encounters great difficulties in numerical stability and applications. Therefore, it still remains a challenging problem to evaluate the pressure–velocity correlation by clear mathematical expression [35].

In this study, the Reynolds stress model (RSM) is utilized as a basic turbulence-modeling method and the pressure–velocity fluctuation correlation transport equations are reconstructed by the elliptic relaxation theory [36]. The corresponding discretization is carried out by Green’s function to make sure the whole sets of expressions are in accordance with N-S functions [29]. Cavity flow calculation is performed to confirm the proposed model by computation. The prominent advantage of this new method is that pressure–velocity correlation terms are solved without any discarding. Corresponding results show good agreement with PIV tests [37].

2. Methodology

2.1. Basic Model and Discussions on Pressure–Velocity Correlation Term

Before presenting the transport equation of the pressure–velocity correlation term, a brief derivation of the pressure–velocity correlation term will be extracted based on the TKE and turbulence energy dissipation equations. The TKE ($k=\frac{u^2+v^2+w^2}{2}$) equation can be deduced straightforwardly by Reynolds stresses transport equation. The steady TKE equation is expressed as:

$$\rho U_i\frac{\partial k}{\partial x_i}=\underset{D_k}{\underbrace{\frac{\partial }{\partial x_i}\left(\mu \frac{\partial k}{\partial x_i}-p{u}_i-\rho k^{\prime }{u}_i\right)}}-\underset{G_k}{\underbrace{\rho u_i{u}_j\frac{\partial U_i}{\partial x_j}}}-\underset{\epsilon }{\underbrace{\mu {\left(\frac{\partial u_i}{\partial x_j}\right)}^2}}$$

(1)

where $\rho$ is the fluid density (kg/m³), $U_i$ is the mean wind velocity components (m/s), $\mu$ is the fluid viscosity (kg/m·s⁻¹), and $p$ and $u_j$ are fluctuation pressure and fluctuation velocity; $k^{\prime }$ denotes the instantaneous TKE. In other words, Equation (1) expresses the fact that a small fluid element carries along its TKE changes among the processes. It is the result of a net imbalance of diffusivity transport due to viscosity and turbulence (associated with both velocity and pressure–velocity interactions),$D_{k^{\prime }}$ production of TKE $G_{k^{\prime }}$ and dissipation of TKE by fluid viscosity $\epsilon$.

Among all the terms mentioned in Equation (1), the term $p{u}_i$ is denoted as $D_{ij}^p$ in the following. Especially in the buffer and outer regions of boundary layer, the $D_{ij}^p$ term makes the main contribution here. Additionally, in the central core of the flows, the $D_{ij}^p$ term should be paid attention to according to the direct numerical simulation (DNS) results. However, the instantaneous turbulence energy term $\rho k^{\prime }{u}_i$, contributing not that much, can be combined into the turbulence kinetic viscosity $u_t$ in the TKE transport equation without further modeling and calculation. Therefore, the reorganized TKE equation can be expressed as:

$$\rho U_i\frac{\partial k}{\partial x_i}=G_k-\epsilon +\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]-\rho \frac{\partial }{\partial x_i}\left(p{u}_i\right)$$

(2)

where ${\mu }_t$ is the modified turbulence dynamic viscosity.

Previous methods usually use a third-order momentum function for the $D_{ij}^p$ term, which can be applied in homogeneous turbulent flows with the standard turbulence models. However, it is challenging to apply the above method in inhomogeneous turbulent flows. In this paper, the $D_{ij}^p$ term will be evaluated by solving a second-moment closure model expression of the pressure correlations in TKE equation. The pressure field of incompressible turbulent flows satisfies the following Poisson equation

$${\nabla }^{2}p\left(x\right)=-\rho \left(\frac{\partial U_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}+\frac{\partial U_j}{\partial x_i}\frac{\partial u_i}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}-\overline{\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}}\right)$$

(3)

where $U_i$ and $u_i$ ($i=x$,$y$,$z$) denote the three mean and fluctuation velocity components along the streamwise, spanwise, and vertical directions, respectively. $p\left(x\right)$ denotes the kinematic pressure. The last term in Equation (3) is the average value of velocity gradients.

Equation (3) can be solved by the Green function $G_{\Omega }$, which expresses $p\left(x\right)$ in terms of an integral of the entire volume of the velocity field, together with integrals over the boundaries. The Green function for Laplacian operator is

$$G\left(x,x^{\prime }\right)=\frac{-1}{4\pi \left|x-x^{\prime }\right|}$$

(4)

The general solution of Equation (3) by the Green’s function in Equation (3) takes the following expression

$$p\left(x\right)={\displaystyle {\int }_{\Omega }{\nabla }^{2}p\left(x^{\prime }\right)G_{\Omega }\left(x,x^{\prime }\right)dV\left(x^{\prime }\right)}$$

(5)

where $x$ and $x^{\prime }$ denote position vectors and $dV\left(x^{\prime }\right)$ the elementary volume. Based on Equations (4) and (5), the fluctuating pressure is usually decomposed into three components [9], corresponding to rapid pressure, slow pressure, and wall influence terms, while in this paper, $D_{ij}^p$ will be treated as a whole term, which can be derived from the divergence $\partial /\partial x_i-\partial /\partial r_i$ of Reynolds stress transport equation, which leads to a Poisson equation for $D_{ij}^p$. The integration of Equation (5) is:

$$\rho p\left(x\right)u_j=\rho D_{ij}^p={\displaystyle {\int }_{\Omega }{D}_{ij}^p{G}_{\Omega }\left(x,x^{\prime }\right)dV\left(x^{\prime }\right)}$$

(6)

where $D_{ij}^p$ suggests the correlation between the velocity and the Laplacian of the pressure:

$$D_{ij}^p=\overline{u_i\left(x\right)p\left(x^{\prime }\right)}$$

(7)

The right-hand part of the Poisson Equation (3) consists of linear and nonlinear velocity fluctuation terms; hence, the pressure–velocity correlation terms in the evolution equation are the n-point velocity moments.

2.2. The Elliptic Relaxation Equation

In Equation (7), the correlations between the velocity and the Laplacian of the pressure need to be modeled by mathematical expression. Durbin defined a correlation function [32]:

$$D_{ij}^p=D_{ij}{}^{\prime }f\left(x,x^{\prime }\right)$$

(8)

$f\left(x,x^{\prime }\right)$ is modeled by:

$$f\left(x,x^{\prime }\right)=\exp \left(\frac{-\left|x-x^{\prime }\right|}{L}\right)$$

(9)

where $L$ is the turbulence integral length scale, which will be defined in the following section. The transport equation of $D_{ij}^p$ can be constructed on the basis of the above expressions. In a free space $V$ using Equation (9), the pressure–velocity correlation term can be written as:

$$\rho D_{ij}^p={\displaystyle {\int }_{\Omega }{D}_{ij}{}^{\prime }\frac{\exp \left(-\left|x^{\prime }-x\right|/L\right)}{4\pi r}dV\left(x^{\prime }\right)}$$

(10)

In the above equation, $D_{ij}^p$ appears as a convolution product between $D_{ij}{}^{\prime }$ and $\frac{\exp \left(-\left|x^{\prime }-x\right|/L\right)}{4\pi r}$, which is the free-space Green’s function associated with the operator $-{\nabla }^2+1/L^2$.

Based on the Durbin’s work [32], it is concluded that the existence of open cavity, wall boundary, and the corresponding roughness effects not only extend the turbulent structures, but also change Reynolds stresses normal to the wall. It appears that the correlation function cannot be expressed by an elementary exponential function. Actually, the two-point fluctuation quantities’ correlations are much stronger toward the wall boundary than apart from the wall boundary. It is rationally derived that the fluctuation velocity and the fluctuation pressure share an identical behavior. In order to modify the elliptic relaxation equation, a new model is derived based on DNS analysis results, which yields the elliptic relaxation expression as:

$$D_{ij}^p-L^2{\nabla }^2{D}_{ij}^p=-\frac{L^2}{\rho }\left(\overline{u_i\left(x\right){\nabla }^{2}p\left(x^{\prime }\right)}+\overline{u_j\left(x\right){\nabla }^{2}p\left(x^{\prime }\right)}\right)$$

(11)

The first purpose is to find ways to solve the above pressure–velocity correlation Equation (11). For instance, the IP model assumed isotropization of the TKE production. Rotta split the correlation term into rapid and slow parts. The corresponding simulation results show that both of the IP and Rotta’s models overestimate the energy redistribution in the logarithmic layer. Following Gibson and Launder’s method, the wall echo terms are introduced into the TKE transport equation, which reimburse the TKE overestimation in the logarithmic layer. These terms are expressed by the elliptic relaxation equation. For the sake of simplification and practical application, Durbin proposed a quasi-homogeneous model $D_{ij}^h$ instead of the right-hand side of Equation (11), which leads to the original model in [32]:

$$D_{ij}^p-L^2{\nabla }^2{D}_{ij}^p=D_{ij}^h$$

(12)

In the logarithmic layer, the Reynolds stresses are consistent with $y^{-1}$, while the production and the dissipation behave with that value. Thus, Durbin’s proposal can be treated as a simple way of extending the quasi-homogeneous models from the core of the flow domain to the wall boundary with appropriate boundary conditions for $D_{ij}^h$. Thus, in order to model the near-wall characteristic flows correctly, it is desirable to choose a proper expression on the right-hand side in Equation (12), and this will be introduced in Section 2.3.

2.3. Length Scales

The turbulence characteristic scales, which traditionally define the Re number, describe the eddy influence on the mean flows. In isotropic turbulence, the characteristic scale has only a single value. However, almost all the wind-engineering-related practical turbulent flows are anisotropic, especially at the large scales. To model such flows, it is required to find its predictable features.

In the near-wall region, from the DNS data, it is noted that the correlation length scale $L$ and the turbulent integral length scale behave differently. However, the correlation length scale $L$ is very close to the integral length scale in the main part of the flows. The correlation function has great influence on the flow separations. To clarify the asymptotic log-layer behavior, new formulations will be established to explain strong anisotropy and inhomogeneity in the cavity and the wall boundary. Additionally, an asymmetric transport equation will be proposed using spatially variable length scales $L_{ij}$, which can be expressed as:

$$L_{ij}=\frac{3}{2}\frac{u_i{u}_j}{k}{L}_K$$

(13)

where $u_i{u}_j$ is defined in the equation as Reynolds stress transport equations. $L_K$ is the Kolmogorov length scale, which is adopted in the near-wall flows:

$$L_K=C_L\max \left(C_{\eta }\frac{{\nu }^{3/4}}{{\epsilon }^{1/4}};\frac{k^{3/2}}{\epsilon }\right)$$

(14)

The above analysis addressed some important issues in the pressure–velocity correlation calculations. It confirmed the necessity of restricting the length scale around the wall by Kolmogorov length scale. It is also noted that modeling the correlation function $f\left(x,x\right)$ by a simple isotropic exponential function is too crude. Reformulations of the original model will be introduced in the following part.

2.4. Correction of the $D_{ij}^p$ Transport Equation

The original correlation function $f\left(x,x\right)$ does not consider the anisotropic effects. This shortcoming will overestimate the turbulence quantities, especially for the TKE generation and pressure terms. Taking the asymmetric normal to the wall into account, due to t, from the DNS data, it can be noted that $f\left(x^{\prime },x^{\prime }\right)$ behaves with ${y^{\prime }}^{-2}$. With flow driving away from the wall, $f\left(x^{\prime },x^{\prime }\right)$ becomes increasingly isotropic, and is only slightly elongated in the streamwise direction. It is specified that the original pressure–velocity correlation term does not account for the elongation in the streamwise direction, close to the near-wall region. Meanwhile, the correlation expression (a function of distance $r$) cannot distinguish different directions. Especially for complex practical flows, this shortcoming becomes significant. Thus, a modified elliptic relaxation equation will be proposed in this section, which combines different length scales in each direction to account for the anisotropy of the turbulent flow structures.

When the integral is a convolution product, Equation (10) can be inverted to give an appropriate expression of the above elliptic relaxation equation, which can be used in CFD simulation. Moreover, in order to take length scale depending on the location into account, the length scales of Equation (13) are introduced into the elliptic relaxation equation.

The $D_{ij}^p$ equation can be expressed as:

$$D_{ij}^p-L_{\alpha }^2{A}_{kl}{A}_{ml}\frac{{\partial }^2{D}_{ij}^p}{\partial x_k\partial x_m}-L_{\alpha }^2{A}_{ml}\frac{\partial A_{kl}}{\partial x_m}\frac{\partial D_{ij}^p}{\partial x_k}=D_{ij}^h$$

(15)

The only difference between Equation (15) and Equation (11) is the third term, which accounts for the length scales’ variations and near-wall flow redistribution.

To introduce the asymmetric characteristics to the correlation function and also to satisfy the original transport equation, gradient of the length scale is adopted in this paper. A matrix of length scales can be defined as $A_{ij}=L_{ij}/L_{\alpha }$, which yields a new sort of the elliptic relaxation expression:

$$D_{ij}^p-L_{kl}{L}_{ml}\frac{{\partial }^2{D}_{ij}^p}{\partial x_k\partial x_m}-L_{ml}\frac{\partial L_{kl}}{\partial x_m}\frac{\partial D_{ij}^p}{\partial x_k}=D_{ij}^h$$

(16)

The new formulation can consider length scale anisotropy, which has been found to be very significant in the cavity flows. It also does overcome the spurious errors in the logarithmic layer.

In fact, no matter what expression of $D_{ij}^h$ is adopted, the same analysis will lead to an expression of the form:

$$D_{ij}^p=\Gamma D_{ij}^h$$

(17)

The amplification factor $\Gamma$, which characterizes the effect of near-wall re-distribution in the elliptic relaxation equation, can be expressed as:

$$\Gamma =\frac{1}{1+\gamma C_L^2{C}_{\mu }^{-3/2}{\kappa }^2}$$

(18)

where $C_L$, $C_{\mu }$, and $\kappa$ are the turbulence model coefficients. Manceau and Hanjalic derived three different models ($M_1$, $M_2$, $M_3$) of coefficient $\gamma$ involving the gradient of the length scales, aiming to account for the anisotropy. Among the three models, the parameter $\beta$ of model $M_3$ can be chosen by the amplification factor from 0 to 1.51. The coefficient $\beta$ enables adjustment of weights of the diffusion and source terms. It can be easily concluded that the model becomes neutral ($\Gamma =1$) when $\beta =1/12$.

The detailed results will show the merits using the proposed form of the elliptic relaxation equation in cavity flow.

3. Solution Strategy

3.1. Turbulence-Modeling Method and Mesh Schemes

As derived in the above sections, the pressure–velocity correlation model is based on the RSM by introducing the pressure–velocity correlation term $p{u}_i$. The new model is implemented into FLUENT, and the corresponding numerical results are compared with the PIV data. The empirical constants of the transport equations and the $k-\omega$ SST model adopted here are the same as the original ones. The modified functions are summarized here:

$$\left\{\begin{array}{c}\frac{D\overline{u_i{u}_j}}{Dt}=-\overline{u_i{u}_k}\frac{\partial U_j}{\partial x_k}-\overline{u_j{u}_k}\frac{\partial U_i}{\partial x_k}+k{D}_{*ij}-\frac{\overline{u_i{u}_j}}{k}\epsilon \hfill \\ +\frac{\partial }{\partial x_l}\left(\frac{C_{\mu }}{{\sigma }_k}\overline{u_l{u}_m}T\frac{\partial \overline{u_i{u}_j}}{\partial x_m}\right)+\nu {\nabla }^2\overline{u_i{u}_j}\hfill \\ \frac{D\epsilon }{Dt}=\frac{C_{\epsilon 1}P-C_{\epsilon 2}\epsilon }{T}+\frac{\partial }{\partial x_l}\left(\frac{C_{\mu }}{{\sigma }_{\epsilon }}\overline{u_l{u}_m}T\frac{\partial \epsilon }{\partial x_m}\right)+\nu {\nabla }^2\epsilon \hfill \\ D_{ij}^p-L_{kl}{L}_{ml}\frac{{\partial }^2{D}_{ij}^p}{\partial x_k\partial x_m}-L_{ml}\frac{\partial L_{kl}}{\partial x_m}\frac{\partial D_{ij}^p}{\partial x_k}=D_{ij}^h\hfill \end{array}\right.$$

(19)

The cavity model considered in this study for comparison between the numerical predictions and water tunnel test results is a scaled rigid model described by Liu and Katz [34]. The numerical model for the pressure–velocity correlation simulation is shown in Figure 1. The L (444.5 mm), W (50.8 mm), and H (120.0 mm) are length, width, and height of the overall model, respectively. A 38.1 mm-long, 50.8 mm-wide, and 30.0 mm-deep cavity is constructed, which is about 190.5 mm from the inlet boundary and 177.8 mm from the outflow boundary. A contraction and a diffusion ramp are arranged along the leading wall and downstream wall of the cavity, which is sketched in Figure 1. The mean velocity above the cavity is about 10 m/s, which corresponds to Reynolds numbers of about 3.4 × 10⁵.

Although the geometrical shape of the model is simple, the mesh scheme should take into account not only the near-wall flow characteristics, but also the computation efficiency. The mesh adjacent to the ground wall and aligned with cavity surfaces must be refined and located within the viscous boundary layer. A hexahedral mesh arrangement is applied to the overall domain, including the cavity area. Since the local Re number around the cavity is relatively large, it is required to ensure that the local Re number $y^{+}$ satisfies $5<y^{+}=\rho uy/\mu <12$. Thus, the height of the first-layer cell centroid is set to be 0.1 mm and the growth ratio is 1.15, as shown in Figure 2.

3.2. Numerical Discretization and Boundary Conditions

All the transport equations are solved by the pressure-based method with the semi-implicit method for pressure-linked equations (SIMPLE). The third-order scheme for convective kinematics (QUICK) difference scheme is used for spatial discretization. Moreover, second-order discretization is adopted for the turbulence transport quantities and pressure–velocity correlation term.

It is obvious that different inlet profiles can cause resulting discrepancies, whether in experiment measurements or numerical simulations. In order to obtain better agreement between the experiment and numerical results, the boundary conditions adopted in the numerical simulations should be the same as those in the experiments, especially for the inlet profiles of mean and turbulent flow quantities. In this simulation, a uniform velocity inlet boundary condition is adopted to simulate to water inflow. The pressure outlet boundary condition is used for water outflow. Turbulence has great influences on fluctuation pressure and fluctuation velocity, mainly on the dynamics of the shear layer. In this study, turbulence is generated by the roughness robs in front of the flow domain. It could be predicted that a turbulence intensity spike value would exist in the near-wall region, which is larger than that generated in the water tunnel test. All the boundary conditions are listed in

Table 1. Boundary conditions summary.

Position	Boundary Condition	Expressions
Inflow	Velocity–inlet	$U_{in}=7.5$ m/s
Outflow	Pressure–outlet	\
Cavity walls	No-slip boundary wall	$U_x=U_y=U_z$
Top, side and ground walls	Slip boundary wall	\

.

3.3. Grid Independence

Two structured mesh plans with different resolutions were tested to examine the grid independence, namely mesh plan G1 and G2, as listed in

Table 2. Grid independence tests summary.

Mesh Plan	Total Grid Numbers	Total Mass (m3/s)
G1	1,936,990	0.0463
G2	3,228,316	0.0462

. An O-type hexahedron-structured grid was applied in the above two mesh schemes as shown in Figure 2. Attention was focused on the grid sizes along the open cavity length and depth directions ($\Delta x$ and $\Delta y$). Therefore, the domain is divided uniformly into 15 grid points for all the two meshes in the z direction. The open cavity between the front and rear walls is discretized by Nx × Ny = 30 × 40, 40 × 50 grids, respectively, for G1 and G2. The first-layer grid distance $\Delta z$ equals 5 × 10⁻⁴ for both cases, which obeys the empirical function $0.1H/\sqrt{\mathrm{Re}}$.

The results of the grid sensitivity analysis are shown with regard to the mass flow $m_a$ across the central surface in Table 2. It is apparent that $m_a$ values are almost insensitive to the grid plan G1 and G2. Therefore, mesh scheme G1 is considered fine enough, which is chosen to resolve the present Reynolds number and the flow recirculation in the cavity. The G1 mesh scheme is shown in Figure 2.

4. Results and Discussion

For examining the performance of the new model, the numerical results of mean velocity distribution, Reynolds stresses, and pressure–velocity correlations are compared with the corresponding water tunnel measurement data in this section.

4.1. Mean Flow Distributions

To illustrate the whole flow field of the cavity flow, Figure 3 gives the mean velocity distributions along the streamwise direction. Due to the contraction of the section area around the cavity, the maximum velocity is about 11 m/s, which is defined as reference velocity $U_e$. The approaching flow from upstream of the cavity separates at the cavity lip and begins to develop and reattaches downstream over the cavity space. The boundary layer grows during this process. Compared with the PIV data, the recirculation center of CFD is higher than that of the test as shown in Figure 3. It is more far away from the ground of the cavity.

Figure 4 shows the mean streamwise velocity distribution overlapped with the streamlines in the open cavity. There exists an obvious large recirculation region within the cavity. The maximum mean flow velocity is about 3.7 m/s around the recirculation core. Above the recirculation region, there exists a small shear layer with large velocity gradients, the thickness of which is about 3 mm. Above the shear layer, the mean horizontal velocity increases slightly in the horizontal direction due to the contraction of the section area, the growth of the shear layer, and the expansion of the boundary layer on the upper wall.

4.2. Reynolds Stress Distributions

The fluctuation part of the flow is represented by considering the normalized Reynolds stress components $\overline{u^2}/U_e^2$, $\overline{v^2}/U_e^2$, and $\overline{-uv}/U_e^2$ distributions, which are shown in Figure 5a,b,c, respectively. It can be summarized that intense turbulent transports occur along the small shear layer. The maximum Reynolds stress values exist at the upper corner of the trailing wall (right wall) of the cavity, and remain high until close to the trailing edge. Along the downstream direction, the turbulent transports extend into the cavity to the leading wall (left wall) of the flow recirculation center.

Generally, the maximum $\overline{u^2}/U_e^2$, $\overline{v^2}/U_e^2$, and $\overline{-uv}/U_e^2$ values are almost the same as the test measurement data, compared with the PIV results, while it could be noted that the peak value area is larger than that of the test results, which is caused by the thicker shear layer predicted by RSM. The peak values and areas of the $\overline{u^2}/U_e^2$ component are similar to the test data, while for the $\overline{v^2}/U_e^2$ component, the maximum value is about 27% larger than that of the test data and it remains large along the trailing wall to the cavity ground. The CFD data give smaller values, especially on the leading wall of the cavity. Similar to the distribution of $\overline{u^2}/U_e^2$, $\overline{-uv}/U_e^2$ reaches its maximum value along the shear layer to the trailing edge wall, while the peak value area is smaller than that of the test data.

4.3. Pressure–Velocity Correlation Distributions

Distributions of $pu/\rho U_e^{-3}$ and $pv/\rho U_e^{-3}$, the pressure–velocity correlations, are presented in Figure 6a,b, respectively. As is obvious, fluctuation velocity $u$ and fluctuation pressure $p$ are negatively correlated in most parts of the shear layer. This negative correlation explains why the pressure increases with the fluctuation flow velocity decelerating in the streamwise direction and vice versa. As shown in Figure 6a, the $pu/\rho U_e^{-3}$ component changes its sign from the shear layer to the trailing edge of the cavity. Meanwhile, the $pu/\rho U_e^{-3}$ value decreases its magnitude to zero gradually, and then increases to a positive value of about 0.005, while for the $pv/\rho U_e^{-3}$ correlation shown in Figure 6b, fluctuation velocity $v$ and fluctuation pressure $p$ are positively correlated in most parts of the shear layer. When approaching the trailing edge wall of the cavity, the $pv$ correlation also changes its sign and magnitude.

Generally, the pressure–velocity correlations $pu/\rho U_e^{-3}$ and $pv/\rho U_e^{-3}$ give similar results to those of the test measurement data, compared with the PIV results. These trends are associated with the adverse mean pressure gradients near the upper corner of the trailing edge of the cavity wall. When the mean velocity increases due to the contraction of the section area, the adverse pressure gradients also increase simultaneously. Then, there exists a positive pressure–velocity correlation. Additionally, this procedure can be found in the Reynolds stress component $\overline{-uv}/U_e^2$, which is presented in Figure 5c.

5. Conclusions

By solving the pressure–velocity correlations individually and introducing these terms into the RSM, the elliptic relaxation method is a promising way to model the pressure–velocity correlation terms within the whole flow field correctly. The first modification is based on the length scales in the elliptic relaxation equation that cannot be considered locally as constant. A new formulation, considering this feature, includes the effects of anisotropy of turbulence through a rather complicated tensor expression. In the second modification, an asymmetric correlation function was introduced, by using the gradient of the length scales to identify the directions of inhomogeneity. These modifications result in a new formulation of the elliptic relaxation equation to model the pressure–velocity correlations. Based on the physical insights gained through the water tunnel test PIV data, CFD data are compared with the corresponding results.

Numerical simulations of water flows over an open cavity model were conducted as a validation study. The performance of the modified turbulence model with a high Reynolds number was investigated in detail. The corresponding numerical results were compared with the water tunnel test PIV data and the conclusions are summarized below:

The mean velocity distributions along the streamwise direction suggest that there exists a slightly favorable pressure gradient in the free stream region due to the contraction of the section, with the growth of the shear layer. The fluctuation velocities and fluctuation pressure show different correlation characteristics along the streamwise direction, which also can be approved by Reynolds stress distributions.

For most flow around an open cavity, the pressure is negatively correlated with the streamwise velocity. Despite the favorable pressure gradient along the shear layer to the trailing edge of the cavity, the pressure–velocity correlation becomes obvious near the upper corner of the cavity. Thus, the pressure–velocity correlation should be emphasized in RANS simulations of flow with obvious impact and separation flows.

Further studies, including the simulation of flows other models using the new formulations, will be carried out using RANS and LES. It is predicted that the pressure–velocity correlation could improve the flow and pressure distributions around the wind turbine applications.