Second Moment Closure Modeling and Direct Numerical Simulation of Stratified Shear Layers

Wake flows differ markedly in stratified compared to nonstratified environments. Specifically, numerous complex physical phenomena arise, including strong Reynolds stress anisotropy, wake flattening/collapse, countergradient fluxes, the coupling between kinetic and potential energies, internal gravity waves, and regime transition.

Many researchers have studied stratification in homogeneous flows, and these studies have led to a community focus on salient physics and flow parameters that characterize these complex turbulent systems. Numerical analysis focusing on 2D stratified mixing layers has been performed in Refs. [1–4], including flows with initially turbulent perturbations [5]. However, these turbulence resolving direct numerical simulations (DNSs) and large eddy simulations (LESs) require significant computational resources at the higher Reynolds and Richardson numbers observed in many practical applications. Accordingly, Reynolds-averaged Navier– Stokes (RANS) methods are attractive in these flows, provided they can return sufficient accuracy.

Reynolds-averaged Navier–Stokes-based eddy viscosity models (EVMs) are inherently incapable of accommodating the physics of stratified flow where the stratification is of importance due to underlying stress anisotropy and attendant wake stabilization and regime transition, particularly at high Richardson numbers [6]. An alternate approach is to use second moment closure (SMC) methods. SMCs abandon the Boussinesq approximation that underpins eddy viscosity models, and solve individual transport equations for the Reynolds stresses, fluxes, and the density variance. SMC-based methods have demonstrated superior predictive capability to EVMs for many flows, see, e.g., Refs. [7–10]. Historically, these methods have been less widely applied than EVMs due to increased stiffness, complexity, and computational cost associated with solving additional transport equations, and for being numerically nonrobust. However, significant increases in available computational resources have allowed for adequate mesh density for 3D applications, leading to stable, practical, and robust applications of SMC methods [11–14].

To date, experimental measurements of important modeling terms (pressure-redistribution, for example) have been very challenging and unreliable [15], and detailed DNS/LES studies that yield such terms “exactly” have been restricted to simple flow configurations. Consequently, calibration of the numerous SMC model constants has been limited to flows using reduced forms of the RANS equations [16], rather than on the accuracy with which the individual submodels approximate the exact terms. Although these key simplifications enable calibration, they do not generalize well to more complex flows, e.g., separated flows and simple shear flows [17,18].

Research on Reynolds stress RANS modeling has focused primarily on improving model performance in reproducing important first and second moment statistics, turbulence energetics [19–21]. Because of the recent availability of DNS and LES predictions, it is now feasible to conduct detailed comparisons of the individual termwise performance of SMC submodels with the exact terms in the governing equations, even in complex dynamical systems like stratified flows. The aim of this work is to perform such detailed analysis.

Here, the flow dynamics of initially turbulent, temporally evolving stratified and nonstratified shear layers are studied. Full Reynolds stress RANS and DNS are applied, with the goal of quantifying model shortcomings and inconsistencies. The paper is organized as follows: A summary of the governing equations, numerical schemes, simulation parameters, and the initial and boundary conditions is provided in Sec. 2. The RANS SMC submodels investigated are delineated in Sec. 3. The results from the two different flow configurations are presented in Sec. 4. Finally, conclusions of the current work are provided in Sec. 5.

The flow configuration considered is a temporally evolving shear layer that develops when two miscible fluids with velocities equal in magnitude but opposite signs are brought together. The mixing layer generated is developed in a stably stratified fluid with a linear density gradient as shown in Fig. 1.

where in Eqs. (9), (10), and (11), production terms, $Pij,PiT,PTT$⁠, and buoyancy terms, $Bij,BiT$⁠, are retained in their exact form at the level of SMC. The remaining terms require modeling as summarized in Sec. 3.

The NPhase solver employs a segregated pressure-based methodology with a collocated variable arrangement and lagged coefficient linearization method (see Ref. [24], for example). A diagonal dominance preserving, finite volume spatial discretization scheme is selected for the momentum and turbulence transport equations. Continuity is introduced through a pressure correction equation, based on the SIMPLE-C algorithm [25]. The cell face fluxes are generated through a momentum interpolation scheme [26], which introduces damping in continuity equation. At each iteration, the discrete momentum equations are solved approximately, followed by a more exact solution of the pressure correction equation. The turbulence scalar and enthalpy equations are then solved in succession. For further details, refer to Ref. [23].

The AFiD-DNS solver is an open-source code that uses a second-order finite difference scheme and a staggered grid to solve for fluid velocities, and a second-order Adams-Bashforth method for time discretization. Further details can be found in Ref. [22]. For a better convergence of higher order statistics, an ensemble average of 100 AFiD-DNS realizations was found sufficient and used here for results.

The UCSD-DNS solver employs a staggered grid with normal velocities stored at cell faces, and density and pressure at cell centers. A second-order central difference scheme for spatial discretization and low storage third-order Runge–Kutta method for temporal integration is employed. A sponge region is employed at the top and bottom boundaries to control spurious reflections. A detailed description can be found in Refs. [3] and [5].

Details of the different simulation parameters including the computational domain size $(L1,L2,L3)$⁠, node count $(N1,N2,N3)$⁠, Prandtl number (Pr), initial gradient Richardson number $(Rig,0)$ and Reynolds number $(Re0)$ are listed in the Table 1.

For NPhase, the periodic cyclic and symmetric boundary conditions are specified along the streamwise (x₁) and spanwise (x₃) directions, respectively. In the vertical (x₂) direction, in addition to Eq. (18), a symmetry boundary is specified for Reynolds stresses, temperature fluxes, variance, and specific dissipation rate.

A short description of the models considered to close the system of Eqs. (9), (10), and (11) is presented below, and all model coefficients are listed in Tables 2 and 3.

where for R, a constant value of 1.5 recommended by Ref. [33] was used here.

The RANS and DNS predictions for two different test cases are considered here: a nonstratified shear layer (S0) and a stratified shear layer (S1). In the S0 case, the mixing layer is evolved with a linear density gradient (with Boussinesq approximation). For this case, the buoyancy terms are removed from the equations to make the flow dynamics identical to a nonstratified system. The resulting system develops with the density/temperature fields behaving as a passive scalar. In the stratified shear layer (S1) case, the buoyant effects are fully included resulting in a complex dynamic system. The simulation details are listed in Table 1.

where $Cδ=0.16$ and $Dω=5$⁠. In Fig. 2, the nondimensional momentum thickness is plotted versus nondimensional time for the RANS model (only SSG with DH model results are presented for the sake of brevity), Eq. (28), and the two DNS simulations.

The SMC model successfully captures the linear momentum thickness growth rate and thus the self-similar state (for SD diffusion case as well). In Fig. 3, profiles of the three normalized RMS velocity fluctuations and square root of shear stress are compared between the RANS (SSG with DH model) and DNS calculations, and two different experimental results, at a nondimensional time-step, $tΔU/4δθ,0=50$⁠, which is in the self-similar region. In these plots, the velocity difference, $ΔU$⁠, and length scale, $δω=Dωδθ$ are used for nondimensionalization.

The shear stress profile, plotted in Fig. 3(d), matches the DNS and experimental predictions very well, attendant to the linear growth rate prediction. The $Rij$ profiles are well-predicted by the SMC model with a slight overprediction of $R11$ and underprediction of $R22$⁠. These results also serve to validate the model implementation and the ability of SMC model to capture anisotropy and predict the self-similar state. The centerline evolution of the nondimensional temperature fluxes and RMS of temperature fluctuations for the three simulations (SSG with DH model) are shown in Figs. 4(a) and 4(b). Good agreement between the RANS and DNS results can be seen, and the flux anisotropy is well captured by the SMC model.

Next, the performance of the SMC submodels is studied by comparing, in Fig. 5, the nondimensional budget terms for the $u′1u′1¯, u′2u′2¯$ and $u′1u′2¯$ Reynolds stresses with the exact terms obtained by AFiD-DNS, at $tΔU/4δθ,0=30$ (again, SSG with DH for SMC).

The model performs well except, importantly, for the terms in the vertical and streamwise budgets, where the peak production and pressure-redistribution terms are seen to be overpredicted by RANS. To further understand SSG submodel performance, the exact pressure-redistribution terms from DNS are compared in Fig. 6 (labeled DNS), to what DNS returns for the pressure-redistribution model, Eq. (19) (labeled DNS, SSG). Qualitative agreement is observed with an under-prediction of peak values for the nonzero stress budget components. This gives an assessment of the accuracy of the SSG pressure-redistribution model independent of the full RANS implementation.

The diffusion term in the current flow configuration does not play an important role close to the centerline but is dominant at the shear layer edges. The widely used diffusion models applied in this work (SD and DH) have historically been formulated to model the exact triple-velocity correlation terms (turbulent-diffusion) in the context of assuming that the pressure-diffusion term is negligible [30,32]. In order to assess this approximation for the present flow, in Fig. 7, the exact turbulent-diffusion and pressure-diffusion terms returned by AFiD-DNS are presented, with time extruded onto the horizontal axis.

It is observed that the pressure-diffusion term is of the same order of magnitude as the turbulent-diffusion and hence non-negligible. Despite this, the SMC models (SSG with DH models) predicted the shear layer growth rate and Reynolds stresses accurately. Thus, it appears that neglecting pressure-diffusion and the over-prediction of pressure-redistribution and diffusion terms offset one another in the calibrated SMC model leading to an overall good prediction of bulk properties.

The introduction of buoyancy in the transport equations has a significant impact on the flow dynamics. Compared to the S0 case, the shear layer growth rate is reduced significantly in the stably stratified S1 case due to the damping of vertical velocity fluctuations. The SMC nondimensional momentum thickness growth rate and centerline evolution of gradient Richardson number are compared with the DNS predictions in Fig. 8. The nondimensional momentum thickness from both DNS models exhibits a decreasing growth rate over time and eventually reaches a near-asymptotic state. The SMC models (SD and DH results included here) also return a decreasing growth rate compared to the S0 case, indicating that the effect of stratification is being accounted for. However, the growth rate returned by the SMC is much higher than DNS. Also, RANS underpredicts Ri_g and its growth rate, and fails to capture the plateauing behavior. These S1 shear layer thickness and gradient Richardson number results are quantitatively far less accurate than the shear layer evolution results presented for the S0 case. These incorrect predictions by RANS SMC models imply the need to examine the accuracy of second-order statistics.

The centerline evolution of the nondimensional temperature fluxes and RMS of velocity and temperature fluctuations for RANS (SD and DH models) and AFiD-DNS are compared in Fig. 9. The growth rates predicted by SMCs (both models) match well with the DNS results. For the peak values of RMS velocity fluctuations and streamwise temperature flux, the error in SD model predictions is over 10%, while the DH model is under 5%. The error in vertical temperature flux and temperature variance predictions by both models is over 20%. The decay rate following this peak is significantly underpredicted by SMCs. The DNS results exhibit two different decay rates due to the laminarization of the flow around, $tΔU/4δθ,0≈80$⁠. In the SMC predictions, an initial rapid decay transitions to a slower decay rate around, $tΔU/4δθ,0≈30$⁠. Additionally, the decay rates predicted by the two diffusion models are almost identical, indicating RANS predictions during the decay period is independent of diffusion model choice. The incorrect decay rate predictions signify the need to improve the SMC models to capture the flow dynamics accurately.

The SMC models considered here employ the isotropic dissipation rate, ε, as a key parameter in formulation of all submodels considered. The validity of this assumption is checked by comparing the exact components of the nondimensional dissipation rate tensor, $εij$⁠, with the corresponding isotropic value, $(23ε)$⁠, estimated by AFiD-DNS, at $tΔU/4δθ,0=60$⁠, in Fig. 10.

At early time instants, the flow evolves like a nonstratified system, and the departure from isotropy was found negligible (not presented for the sake of brevity). The buoyancy term increases with time, and a time instant when its impact on the flow is significant is chosen for comparison in Fig. 10. A significant departure from isotropy can be observed, highlighting an important inconsistency in all submodel formulations where the isotropic value is used as a key parameter. This conclusion can be extended to future times where the buoyancy effects are significant. This invalid assumption is hypothesized to be a reason behind inaccurate SMC predictions. Departure from isotropy in the S0 case was found negligible, consistent with the computational and numerical studies at high Reynolds number [38,39]. The anisotropy and complex dynamics introduced by stratification render the assumption invalid, resulting in incorrect SMC decay rate predictions. Thus, this result implies the need to account for anisotropy in the dissipation rate tensor to improve RANS predictions.

The performance of RANS-SMC models was compared against DNS results of the same case. Despite the reasonable prediction of first and second-order statistics for the nonstratified case, individual submodels were shown to depart from DNS in the Reynolds stress budget termwise comparison. The omission of pressure-diffusion term offset by calibration of model coefficients is hypothesized to be a reason for the inconsistencies. For the stratified case, SMC models were shown to capture the complex flow dynamics less accurately. The local isotropy assumption in dissipation modeling was shown to be invalid for the stratified case and hypothesized to be a reason for incorrect RANS predictions.

• U.S. Office of Naval Research (Contract No. N00014-19-1-2057; Funder ID: 10.13039/100000006).

a_ij =

anisotropy tensor, $aij=Rijk−23δij$

$Bij$ =

Reynolds stress buoyancy production tensor

$BiT$ =

temperature flux buoyancy production vector

$Cμ$ =

Prandtl–Kolmogorov constant, $Cμ=0.09$

$Dij$ =

Reynolds stress diffusion tensor

$DiT$ =

temperature flux diffusion vector

$DTT$ =

temperature variance diffusion scalar

g_i =

gravitational vector

k =

turbulent kinetic energy, $k=12(u′iu′i¯)$

L_i =

domain length along x_i direction

m =

wavenumber

N =

Brunt–Vaisala frequency, $N=−g2ρ0∂ρ∂x2$

N_i =

grid points in the ith Cartesian direction

P =

instantaneous pressure

$p′$ =

fluctuating pressure

$P¯$ =

mean pressure

$Pij$ =

Reynolds stress production tensor

$PiT$ =

temperature flux production vector

$PTT$ =

temperature variance production scalar

Pr =

Prandtl number, $Pr=να=1$

q² =

turbulence intensity, $q2=u′iu′i¯$

R_ij =

Reynolds stress tensor

Re₀ =

initial Reynolds number, $Re0=ΔU×4δθ,0ν=640$

$Rig,0$ =

initial gradient Richardson number, $Rig=N2S2$

S =

shear, $S=∂u1¯∂x2$

t =

time

T =

instantaneous fluid temperature

$T′$ =

fluctuating fluid temperature

T₀ =

reference fluid temperature, 1K

$T¯$ =

mean fluid temperature

$T′T′¯$ =

temperature variance

u_i =

instantaneous velocity vector

U₁ =

higher streamwise velocity magnitude, $U1=0.5ms$

U₂ =

lower streamwise velocity magnitude, $U2=−0.5ms$

$u′i$ =

fluctuating velocity vector

$ui¯$ =

mean velocity vector

$u′iT′¯$ =

temperature fluxes

$u′iu′j¯$ =

Reynolds stresses

ν =

fluid kinematic viscosity

x_i =

Cartesian coordinate

$x2,0$ =

centerline vertical coordinate

α =

fluid thermal diffusivity

β =

expansion coefficient at constant pressure

$δθ$ =

momentum thickness at time t, $δθ=∫−∞∞(14−u1¯2ΔU2)dx2$

$δθ,0$ =

momentum thickness at time t = 0, 0.25 m

δ_ij =

Kronecker delta

$ΔT$ =

temperature difference across $4δθ,0$ at t = 0

$ΔU$ =

velocity difference between bottom and top layer, $ΔU=U1−U2=1ms$

ε =

isotropic turbulence dissipation rate, $ε=ν∂u′i∂xj∂u′i∂xj¯$

$εij$ =

turbulence dissipation rate tensor, $εij=2ν∂u′i∂xk∂u′j∂xk¯$

$εiT$ =

temperature flux dissipation rate, $εiT=(ν+α)∂u′i∂xk∂T′∂xk¯$

$εTT$ =

temperature variance dissipation rate, $εTT=2α∂T′∂xk∂T′∂xk¯$

Π_ij =

pressure-redistribution tensor

Π_iT =

pressure-scrambling vector

ρ =

instantaneous fluid density

$ρ′$ =

fluctuating fluid density

$ρ¯$ =

mean fluid density

ρ₀ =

reference fluid density, $1kgm3$

ω =

specific turbulence dissipation rate, $ω=εCμk$