Turbulent Energy Redistribution in Spanwise Rotating Channel Flows

ARTICLE pubs.acs.org/IECR

Turbulent Energy Redistribution in Spanwise Rotating Channel Flows Karuna S. Koppula,† Andre Benard,‡ and Charles A. Petty*,† †

Department of Chemical Engineering and Materials Science ‡Department of Mechanical Engineering Michigan State University, East Lansing, Michigan 48824, United States ABSTRACT: For fully developed turbulent flows of a Newtonian fluid in a nonrotating channel with a large aspect ratio (i.e., spanwise length scale . transverse length scale), the longitudinal component of the mean velocity and the mean pressure are symmetrical functions of the transverse coordinate. The three unequal normal components of the Reynolds stress are also symmetrical functions about the midplane whereas the shear component of the Reynolds stress is spatially antisymmetric. However, if the channel rotates around the spanwise axis, the Coriolis acceleration breaks the spatial symmetry of the flow field. This phenomenon occurs because the angular velocity of the frame directly couples with the linear momentum of the flow field to redistribute the turbulent kinetic energy. Consequently, rotating channel flows are ideal for testing the efficacy of closure models for the Reynolds averaged Navier-Stokes (RANS) equation. It is noteworthy that none of the classical “eddy” viscosity models are able to predict the Coriolis symmetry breaking phenomenon in spanwise rotating fully developed channel flows. This fundamental deficiency occurs because the “eddy” viscosity hypothesis assumes that the Reynolds stress is a frame-indifferent dyadic-valued operator (i.e., objective). Unlike the viscous transport of momentum, the transport of momentum by turbulent fluctuations is directly influenced by the angular velocity of the reference frame. In this paper, a recently developed algebraic closure for the Reynolds stress [see Koppula, K.; Benard, A.; Petty, C. Realizable Algebraic Reynolds Stress Closure. Chem. Eng. Sci. 2009, 64, 4611], referred to hereinafter as the universal realizable anisotropic prestress (URAPS) closure, is used to predict the spatial distribution of the normalized Reynolds (NR) stress (i.e., R__ /tr() for rotating and for nonrotating fully developed channel flows. The new closure is formulated as an algebraic mapping of the NR stress into itself and is, thereby, realizable for all turbulent flows regardless of the specific benchmark flows used to estimate model parameters. The fixed points of the non-negative mapping depend explicitly on three specific time scales associated with the local statistical state of turbulence. The a priori predictions of the NR stress are consistent with previously published direct numerical simulations of the Navier-Stokes equation.

1. INTRODUCTION For constant density fluids, the Reynolds averaged NavierStokes (RANS) equation requires a closure model for the Reynolds stress,-F or, equivalently, the kinematic Reynolds momentum flux, = 2kR__ , where k( tr/2) is the turbulent kinetic energy. In this paper, a new algebraic closure for the NR stress, referred to hereinafter as the universal realizable anisotropic prestress (URAPS) closure, is used to predict the influence of rotation on the redistribution of the turbulent kinetic energy among the three components of the fluctuating velocity field. The new closure stems directly from the Navier-Stokes equation and is formulated as a nonlinear mapping of R__ into itself (i.e., R__ ðR__ , τR < F__ > Þ f R__ ). The mapping operator depends on the mean field operator < F__ > ð r þ 2Ω__ Þ and a scalar-valued turbulent transport time τR. As demonstrated in this paper, the URAPS closure is consistent with the anisotropic redistribution of the turbulent kinetic energy in rotating fully developed channel flows. By contrast, an “eddy” viscosity type closure, which relates the NR stress to the local mean strain rate (i.e., R__ ¼ __I =3 - 2τE < S__ > ), predicts that the turbulent kinetic energy is equally distributed among the three fluctuating components of the velocity for fully developed channel flows in inertial frames as well as in noninertial frames. This unphysical prediction is independent of the “eddy” time scale τE. Clearly, an “eddy” viscosity model is qualitatively inconsistent with the underlying dynamic equation governing turbulent fluctuations. r 2011 American Chemical Society

The objective of this paper is to calculate R__ for a given value of K__ τR < F__ > in the outer region of nonrotating and rotating channel flows. Previously published direct numerical simulation (DNS) results are used to estimate K__ at different positions in the channel. The components of R__ are calculated by solving the nonlinear algebraic equation R__ ¼ R__ ðR__ , K__ Þ by successive substitution. With K__ =K__ DNS, the predicted redistribution of the normal components and the shear component of the NR stress are qualitatively consistent with DNS predictions of the components of R__ . This prevalidation result confirms that the nonlinear algebraic structure of the URAPS equation has captured the important statistical physics that rely on the underlying coupling between the local mean field and the instantaneous velocity. A brief synopsis of the theoretical ideas associated with the URAPS closure is given in Section 2 below. The URAPS-closure depends on a scalar-valued turbulent transport time τR that depends on the local turbulent kinetic energy, k, and the local turbulent dissipation, ε v. The closure also depends on two coefficients, R and β, that depend on the eigenvalues of R__ . Two additional phenomenological coefficients Special Issue: Churchill Issue Received: October 7, 2010 Accepted: January 19, 2011 Revised: January 16, 2011 Published: April 04, 2011 8905

dx.doi.org/10.1021/ie1020409 | Ind. Eng. Chem. Res. 2011, 50, 8905–8916

Industrial & Engineering Chemistry Research

ARTICLE

2. THEORY RANS Equation. The unclosed, albeit exact, RANS equation

for a constant property Newtonian fluid in a rotating rigid-body frame-of-reference with a constant angular velocity Ω can be written as follows1-4 F

D ¼ - F 3 < F__ > Dt þ r 3 ð - < pH > __I þ 2μ < S__ > - F < u0 u0 > Þ

ð1Þ

In the above equation, < F__ > r þ 2Ω__

ð2Þ

< pH > - F½x 3 ðΩ__ 3 Ω__ T Þ 3 x=2 - Fg 3 x

ð3Þ

< S__ > ½r þ ðr ÞT =2

ð4Þ

The right-hand-side of eq 1 shows that the rate-of-change of mean momentum per unit volume of fluid in a rotating frame-of-reference depends on five factors, viz., a coupling between the mean velocity field and the mean velocity gradient, -F 3 r; a coupling between the mean velocity field and the Coriolis operator Ω__ , -2F 3 Ω__ ; the divergence of the mean hydrodynamic pressure stress, -I__ ; the divergence of the mean viscous stress, 2 μ; and, the divergence of the Reynolds stress, -F. In the above equations, the position vector x is relative to the rotating frame-ofreference. The rotational operator Ω__ is defined as follows, Ω__ ε__ 3 Ω. The operator _ε__ is the permutation triadic and the operator __I _ is the unit dyadic. The angular velocity of the frame Ω is constant and Ω__ is an antisymmetric operator. The hydrostatic forces per unit volume due to gravitational acceleration and Coriolis acceleration are included in the definition of . For constant density fluids, the continuity equation implies that r3 ¼ 0

ð5Þ

The mean kinematic operator is traceless. If Ω__ = 0__ , then eq 1 reduces to the unclosed RANS-equation in an inertial frame-of6 0__ , the instantaneous velocity u(x,t)is defined reference. For Ω__ ¼ relative to the rotating frame. The ensemble averaged velocity in the rotating frame is (x,t) and the fluctuating velocity u0 (x,t) is defined by the Reynolds decomposition, u0 u - . Clearly, the RANS-equation requires a closure model for the Reynolds stress, -F or, equivalently, the Reynolds momentum flux,

þF, which can be written in terms of the turbulent kinetic energy, k, and the NR-stress as follows: F < u0 u0 > ¼ 2FkR__

ð6Þ

Thus, a closure model for eq 1 requires a model for k and a model for R__ . In this paper, the focus is on the closure model for R__ and its ability to predict the distribution of turbulent energy among the three components of the fluctuating velocity. NR Stress. The NR stress, which is an ensemble average of a dyad, must satisfy the following fundamental inequality Q R R__ : z z < ðu0 3 zÞ2 > =ðtr < u0 u0 > Þ g 0 for " z ∈ E3 , z ¼ 1 ) )

are introduced by closure models for the “production” and dissipation of turbulent dissipation. The three coefficients introduced by the URAPS-mapping (i.e., τR, R, and β) have been previously estimated by using benchmark flows related to asymptotic nonrotating homogeneous shear and the decay of nonrotating and rotating homogeneous turbulence. Thus, the results reported herein for the NR stress for different values of K__ DNS are a priori predictions and provide a critical prevalidation test of the algebraic URAPS mapping. The component equations associated with the RANS equation and the URAPS equation for spanwise rotating fully developed channel flows are presented in Section 3. This is followed by Section 4, which compares the URAPS predictions of the NR stress with previously published direct numerical simulation (DNS) results. Section 5 summarizes the salient conclusions supported by the prevalidation results.

ð7Þ

Inequality 7 is a necessary and sufficient condition for the eigenvalues of R__ to be non-negative.5,6 All of the eigenvalues are in the first sextet of the non-negative orthant of a twodimensional hyperplane inasmuch as 3

trðR__ Þ ¼

3

∑ R jj ¼ i∑¼ 1 λRi ¼ 1 j¼1

0 e λR1 e λR2 e λR3 e 1

ð8Þ

The normal components of R__ are also in the non-negative orthant of a two-dimensional hyperplane (i.e., energy simplex) and, depending on the flow, may be distributed over all six sextets. For fully developed channel flows in an inertial frame-ofreference, the energy states are in the second sextet of the energy simplex (i.e., 0 e Ryy e Rxx e Rzz < 1). The primary normal stress difference, Rzz - Ryy, is non-negative at all spatial positions; and, the secondary normal stress difference, Ryy - Rxx, is nonpositive at all spatial positions. However, in a rotating frame, the energy states cover the entire energy simplex and, most significantly, the primary and secondary normal stress differences flip signs in the outer region of the flow. This complex behavior in a simple rotating mean shear flow is due to the explicit influence of the Coriolis acceleration on the NR stress. By definition, a turbulence closure model for the NR stress is realizable if all solutions to the model satisfy ineq 7. Furthermore, if ineq 7 is satisfied in all rigid-body frames (i.e., nonrotating and rotating), then the turbulence model is “universally” realizable. Although these two conditions are fundamentally important for any turbulence model, they do not guarantee the accuracy of any solutions other than the benchmark flows used to estimate model parameters. Recently, Koppula4 (also see ref 8) developed an algebraic closure model for the NR stress as a realizable mapping of R__ into itself (i.e., R__ ðR__ , K__ Þ f R__ , where K__ τR < F__ > ). The parameter τR is a scalar-valued turbulent transport time associated with the relaxation of the turbulent space-time correlation, and the relaxation of a Green’s function associated with a convective/diffusive parabolic differential operator.7,8 This closure for the RANS-equation is referred to hereinafter as the universal, realizable, anisotropic prestress (URAPS-) closure. Unlike other algebraic closure models for the NR-stress (see refs 9-12), the URAPS closure satisfies ineq 7 for all turbulent flows regardless of the set of benchmark flows selected to calibrate closure parameters. The components of the RANS equation and the URAPS equation for spanwise rotating fully developed channel flows are presented in Section 3 below. Previously published DNS results (see refs 13-15) are used to estimate the dimensionless groups that compare the relative importance of the turbulent time scale, 8906

dx.doi.org/10.1021/ie1020409 |Ind. Eng. Chem. Res. 2011, 50, 8905–8916


ARTICLE

)

)

)

)

k/ε; the viscous time scale, ν/k; and, the mean field time scale, 1/ , where (:T)1/2. The components of R__ based on the URAPS closure are determined by using the DNS operator DNS, the DNS turbulent kinetic energy kDNS, and the DNS turbulent dissipation εDNS. Preclosure Mapping and Normalized Prestress. The URAPS closure (see refs 2, 4, 7, 8,16-18) is based on the following preclosure mapping of the normalized prestress B__ into the normalized Reynolds stress R__ R__ ¼ A__ T 3 B__ 3 A__ =trðA__ T 3 B__ 3 A__ Þ ð9Þ In the above equation B__ < f 0 f 0 > =trð < f 0 f 0 > Þ where Ff 0 r 3 ½p0__I þ Fðu0 u0 - < u0 u0 > Þ

ð10Þ

The preclosure operator A__ in eq 9 is defined as follows: A__ ½I__ þ K__ -1 where K__ τR < F__ > and < F__ > r þ 2Ω__ 4

ð11Þ 2

)

)

Koppula developed eq 9 for rotating frames. Earlier, Parks and Weispfennig18 used the preclosure equation (i.e., eq 9) to study homogeneous shear and fully developed channel flows in inertial frames (i.e., Ω = 0). Equation 9, which stems directly from an analysis of the Navier-Stokes (NS) equation and the continuity equation, shifts the phenomenological closure problem from the NR stress R__ to the normalized prestress operator B__ . Like the NR stress, B__ is symmetric, non-negative, and normalized. The eigenvalues of B__ are all located in the first sextet of the positive orthant of a hyperplane defined by tr(B__ ) = ∑3i=1λBi = 1. As indicated by eq 10, the prestress is related to a fluctuating force Ff 0 caused by two factors: the divergence of the fluctuating pressure, p0__I ; and, the divergence of the fluctuating Reynolds momentum flux, F(u0 u0 - ). The operator K__ (see eq 11 above) depends on the mean velocity gradientr, the Coriolis operator 2Ω__ , and, the turbulent transport time τR. Turbulent Transport Time. The turbulent transport time τR in eq 11 depends on three local time scales: a viscous time scale ν/k; a turbulent time scale k/ε; and, a mean field time scale 1/ . Dimensional reasoning implies that

)

)

k τR CR1 ~τ R ðNF , Ret Þ ε where NF ðk=εÞ and Ret ðk=εÞ=ðν=kÞ ð12Þ The URAPS closure uses the following semiempirical equation for ~τR (see Koppula4,8) ! 3=2 ð1 þ CR3 NF Þ ~τ R ¼ Υwall 3=2 ð1 þ CR2 NF Þ where Υwall ½1 þ CR4 ð30=Ret Þ3=4 expð - Ret =30Þ ð13Þ Equation 13 stems from an analysis of rotating homogeneous decay by Park and Chung.19 The four parameters in eqs 12 and 13 are CR1 = 0.036, CR2 = 0.076, CR3 = 0.053, and CR4 = 7. For Ret . 30 and 0 e NF < ¥, eq 13 implies that Υwall f 1 and CR3/CR2 < ~τR e 1. For Ret f 0 and 0 e NF < ¥, Υwall f ¥, and consequently, ~τR f ¥. Near a solid/fluid interface (located at y = 0), the no-slip boundary condition and the continuity equation imply that :ezey y3, k y2, ε f εW > 0, R__ :ezey = Ryz y and

Ret y4(see Monin and Yaglom9 and Weispfennig et al.16). These results are consistent with the URAPS closure provided τR y-1 as y f 0. Equation 13 captures this behavior for NF f 0 and Ret f 0 inasmuch as ~τR(NF,Ret) f CR4(30/Ret)3/4 k-3/2 and τR = CR1kτ~R/ε k-1/2. The empirical exponential function in eq 13 was selected to transition the shear component of the NR stress across the wall sublayer wherein viscous transport of momentum is competitive with turbulent transport of momentum. With CR4 = 7, the DNS and URAPS shear components of the NR stress are about the = RDNS yþ for yþ yu*/ν ,10, where u* same, i.e., RURAPS yz yz 1/2 (τw/F) . The near wall transition function, defined by eq 13, implies that Υwall = 7(30/Ret)3/4 k-3/2 f ¥ as yþ f 0. In this paper, the URAPS closure is used to determine the behavior of the NR stress for Ret . 30; therefore, Υwall = 1. An assessment of the URAPS closure in the near wall region will be presented in a subsequent paper. Normalized Prestress Hypothesis. For nonrotating frames, Parks et al.7 and Weispfennig et al.16 developed a class of algebraic closure models for R__ based on eq 9 by assuming that the anisotropic component of B__ is a dyadic-valued function of the dimensionless strain rate k /ε. This hypothesis implies that the eigenvectors of the prestress operator B__ are collinear with the eigenvectors of the mean strain rate for all turbulent flows. Similarly, the “eddy” viscosity closure20 and its nonlinear generalizations (see refs 3, 11, 21, and 22) have the feature that the eigenvectors of R__ are collinear with the eigenvectors of the mean strain rate for all turbulent flows. For rotating (and nonrotating) frames, Koppula4 developed a closure model for the prestress operator B__ by assuming that B__ =B__ ðR__ Þ. This hypothesis is similar to the 1951 Rotta conjecture that the statistical correlation between the local fluctuating pressure and the local fluctuating strain rate can be represented in terms of the anisotropic component of the NR stress (see Pope,11 page 393). The primary mathematical motivation for the foregoing conjecture for the prestress operator is the observation that both the prestress and the NR stress are nonnegative operators. Thus, with B__ =B__ ðR__ Þ, the preclosure equation for R__ (i.e., eq 9) yields an algebraic equation for R__ that is parametrized by the kinematic/dynamic operator K__ , defined by eq 11. Furthermore, the foregoing hypothesis implies that the eigenvectors of the normalized prestress are collinear with the eigenvectors of the NR stress. A Cayley-Hamilton (CH) representation for the prestress operator in terms of the NR stress can be written as B__ ¼ B__ ðR__ Þ ¼ þR__ þ C1 ½R__ - __I =3 þ C2 ½R__ 3 R__ - ðR__ : R__ ÞR__ T

ð14Þ

T

Note that if R__ =R__ and tr(R__ ) = 1, then B__ = B__ and tr(B__ ) = 1. Moreover, if R__ is isotropic, then eq 14 implies that B__ is also isotropic. The “extra” anisotropic CH coefficients, C1 and C2, in eq 14 depend on the invariants of R__ . Without loss of generality, C1 and C2 can be expressed as C1 ¼ þ ½27detðR__ ÞβðIIR , IIIR Þ

ð15Þ

C2 ¼ - ½IIR - 1=3RðIIR , IIIR Þ

ð16Þ

8907



ARTICLE

In the above equations, IIR tr(R__ 3 R__ ) and IIIR tr(R__ 3 R__ 3 R__ ). If R__ =I__ /3 and if R and β are bounded, then C1 = þβ(1/3,1/9) and C2 = 0. Most significantly, it follows from eq 9 that if B__ : _z_z g0 for all z ∈ E3, then Q R ¼ R__ : z z ¼ ðz 3 A__ T Þ 3 B__ 3 ðA__ 3 zÞ=trðA__ T 3 B__ 3 A__ Þ ¼ ðA__ 3 zÞ 3 B__ 3 ðA__ 3 zÞ=trðA__ T 3 B__ 3 A__ Þ g 0

ð17Þ

Thus, the URAPS-closure (i.e., eqs 9 and 14) is realizable for all rotating and nonrotating turbulent flows provided the “extra” anisotropy coefficients R(IIR,IIIR) and β(IIR,IIIR) satisfy the following universal inequalities4,8 - 3=2 ¼ Rmin < RðIIR , IIIR Þ < Rmax ¼ 9

ð18Þ

- 1 ¼ βmin < βðIIR , IIIR Þ < βmax R=27 þ 4=9

ð19Þ

The above two inequalities guarantee that the CH mapping of R__ into B__ defined by eq 14 is non-negative. In other words, if the eigenvalues of R__ are in the first sextet of the positive orthant within the two-dimensional hyperplane defined by ΣiλRi = 1, then the eigenvalues of B__ are also within the first sextet of the positive orthant defined by ΣiλBi = 1. The upper and lower bounds on R and β noted above are universal bounds and apply for all rotating and nonrotating frames-ofreference. URAPS Closure for NR Stress. A combination of eqs 9 and 14 defines the URAPS closure for the NR stress R__ ðR__ , K__ Þ ¼ R__

ð20Þ

Figure 1. Spanwise rotation of fully developed channel flows.

the turbulent kinetic energy, and the turbulent dissipation. Therefore, R__ nþ1 = R__ (R__ n,K__ m), where K__ m denotes the hydrodynamic/kinematic operator determined by a previous iteration. In Section 4 below, the results of a prevalidation study of eq 20 for spanwise rotating fully developed channel flows are presented. The hydrodynamic/kinematic operator K__ was determined from previously published DNS statistics. The fixed points of eq 20 were determined by using the following successive substitution algorithm R__ n þ 1 ¼ R__ ðR__ n , K__ DNS Þ where ½ðR__ n þ 1 - R__ n Þ : ei ej e ½R__ n : ei ej 10-5

ð22Þ

3. RANS EQUATION FOR SIMPLE SHEAR FLOWS WITH TRANSVERSE ROTATION For fully developed, spanwise rotating channel flows, the mean velocity field has one nontrivial component (i.e., = (y)ez) . The mean hydrodynamic pressure varies in the transverse, longitudinal, and spanwise directions, = (x,y,z). As illustrated by Figure 1, the symmetry plane of the channel is located at y = δ; and, the solid/fluid interfaces are located at y = 0 and y = 2δ. The angular velocity of the channel is constant and has one component: Ω = Ωx ex ; the rotation operator and the transverse component of the angular velocity are related by Ω__ = Ωx[eyez-ezey]. The angular velocity is collinear with the mean vorticity, r∧ = ex, and the acceleration due to gravity g = gxex . The mean velocity gradient has only one dyad, r = Γyz(y)eyez. For this application, the transverse and longitudinal components of the RANS equation (see eq 1 above) are

~__ > and k/ε. It is most where K__ τR g CR1~τR d d < uz > þ μ - 2FkR yz 0¼ Dz dy dy

R__ n þ 1 ¼ R__ ðR__ n , K__ Þ where R__ n : _z_z g 0

The hydrodynamic pressure (x,y,z) is defined as follows

ð21Þ

Thus, not only are all the fixed points of eq 20 non-negative operators, but the sequence of operators R__ 0 , R__ 1 , R__ 2 , 3 3 3 are also non-negative for K__ ( τR). In practice, eq 21 must be solved simultaneously with the RANS equation (i.e., eq 1) and transport equations for the turbulent kinetic energy and the turbulent dissipation (see Section 5 below). Therefore, a solution to the RANS equation based on the algebraic URAPS closure would also involve iterations of K__ inasmuch as the operator K__ depends implicitly on the NR stress through the mean velocity gradient,

0 ¼ þ 2F < uz > Ωx -

D < pH > d þ ½ - 2FkR yy Dy dy

ð23Þ

ð24Þ

ðy2 þ z2 ÞΩ2x < pH > ðx, y, zÞ ðx, y, zÞ - F - Fgx x ð25Þ 2 The mean velocity and the mean axial pressure gradient satisfy the following conditions: < uz > ð0Þ ¼ < uz > ð2δÞ ¼ 0 D < pH > ¼ constant > 0 Dz 8908

ð26Þ



ARTICLE

The no-slip boundary conditions at the channel walls imply that the components of the Reynolds stress, -F, are zero for y = 0 and y = 2δ. Transverse Force Balance. For spanwise rotation, the Coriolis force in eq 1 acts in the transverse (cross-flow) direction and is given by - 2F 3 Ω__ ¼ - 2FΩ_ ∧ ¼ þ 2FΩx < uz > ðyÞey

ð27Þ

Equation 23 implies that the mean hydrodynamic pressure varies in the cross-flow direction due to two factors: changes in the local Coriolis force, and changes in the local transverse normal component of the Reynolds stress. For no-slip boundary conditions, transverse changes in the mean hydrodynamic pressure (i.e., ΔP(y) (y,z) - (0,z)) can be expressed as Z y ð^yÞd^y ð28Þ ΔP ðyÞ ¼ - 2FkðyÞR yy ðyÞ þ 2FΩx 0

For no rotation (Ωx = 0), Δp(y) is symmetric about the symmetry plane with Δp(0) = 0, min[Δp(y)] = Δp(δ) 0 and (y) g 0. Equation 28 shows that the mean hydrodynamic pressure (x,y,z) is not transversally uniform for either nonrotating or rotating fully developed turbulent channel flows. Longitudinal Force Balance. Equations 24 and 26 imply that the viscous stress plus the Reynolds stress is a linear function of the transverse (or cross-flow) coordinate d < uz > D < pH > - 2FkðyÞR yz ðyÞ ¼ - μ y þ τHPW dy Dz d < uz > where τHPW μ >0 ð30Þ dy

The friction velocity together with the kinematic viscosity of the fluid are often used to scale the mean axial velocity: uþ /u* = uþ(yþ;δþ;Ωþ). Thus, solutions to eq 24 depend on three dimensionless groups (yþ yu*/ν; δþ δu*/ν; and, Ωþ = Ωxν/(u*)2). A macroscopic longitudinal force balance follows directly from eq 30 by integrating over the flow domain. With ξ y/δ, the integration implies that ξmax ¼ þ 2=ð1 þ φÞ þ 2k þ R yz ðξmax Þ φ τLPW =τHPW duþ ξmax ’ þ Γþ ðξmax Þ ¼ 0 dy

For no rotation (Ωþ = 0), the mean velocity field is symmetric about y = δ; therefore, φ=1 and ξmax = 1. For this case, eq 33 anticipates that Ryz(1) = 0. DNS results14 for δþ = 300 and Ωþ = 0 support this conclusion. For Ωþ < 0, eq 33 anticipates that 1 < ξmax < 2 if φ < 1 and Ryz(ξmax) = 0. This means that transverse mixing of momentum on the HPW side of the channel is more effective than transverse mixing of momentum on the LPW side of the channel. DNS results15 for δþ = 300 and Ωþ = -0.0042 also support this 6 0 seem to show that inference. The DNS results for Ωþ ¼ 6 0; however, it is difficult to conclude that this occurs Ryz(ξmax) ¼ for Ωþ = -0.0042 in the absence of an analysis of numerical inaccuracies associated with the DNS results. In addition to Table 1 and Figure 3b, other complementary DNS results23,24 at higher values of |Ωþ| support the conclusion that there exist an ξ0 > ξmax ~ ~ for which Ryz(ξ0) = 0, Γ (ξ0) < 0, Ryz(ξmax) < 0 and Γ (ξmax) = 0, ~ where Γ kΓyz/ε. This would imply that the “production” of turbulent kinetic energy, Pk -:r, is negative within the region ξmax < ξ < ξ0 for |Ωþ| > 0. See Churchill25 for a discussion of negative “production” of turbulent kinetic energy in nonrotating annular flows. URAPS Equation for Spanwise Rotating Channel Flows. For fully developed channel flows, the NR stress and the normalized prestress can be represented as follows R__ ¼ R xx ex ex þ R yy ey ey þ R zz ez ez þ R yz ey ez þ R zy ez ey

y ¼ 2δ

As with the no-rotation case, the friction velocity u* for spanwise rotation is defined as follows: τHPW þ τLPW 1=2 δ D < pH > 1=2 u ¼ ð32Þ F Dz 2F

ð34Þ

B__ ¼ Bxx ex ex þ Byy ey ey þ Bzz ez ez þ Byz ey ez þ Bzy ez ey

y¼0

Unlike eq 23, there is no explicit influence of the Coriolis force on the longitudinal force balance; however, an implicit influence of rotation does occur through the shear component of the NRstress and the turbulent kinetic energy. Because k(0)Ryz(0) = 0 and k(2δ)Ryz(2δ) = 0, eq 30 implies that D < pH > þ 2δ ¼ τHPW þ τLPW 2FðuÞ2 Dz d < uz > where τLPW - μ >0 ð31Þ dy

ð33Þ

ξmax

ð35Þ

The component equations for eq 14 are Bxx ¼ R xx þ C1 ðR xx - 1=3Þ þ C2 ðR 2xx - IIR R xx Þ

ð36Þ

Byy ¼ R yy þ C1 ðR yy - 1=3Þ þ C2 ðR 2yy þ R yz R zy - IIR R yy Þ

ð37Þ

Bzz ¼ R zz þ C1 ðR zz - 1=3Þ þ C2 ðR 2zz þ R zy R yz - IIR R zz Þ

ð38Þ

Byz ¼ ½1 þ C1 þ C2 R yy þ R zz - IIR R yz

ð39Þ

The CH coefficients C1 and C2 are defined by eqs 15 and 16; and, the second invariant is IIR R__ :R__ = R2xx þ R2yy þ R2zz þ RyzRzy þ RzyRyz. 8909



ARTICLE

Table 1. DNS Groups in the Outer Region for No Rotation14 (δþ = 300 and Ωþ = 0) and for Rotation15 (δþ = 300 and Ωþ = -0.0042) Ωþ = 0 ξ

Ret

0 0.10a

Ωþ = -0.0042

~ Γ

Ret

Γ~

~ Ω

~ þ 2Ω ~ Γ

0

0

-0.4

þ7.2

0

0

0

0

177

þ5

454

þ8

0.15

193

þ4

492

þ5

-0.5

þ4.0

0.20

212

þ3

532

þ4

-0.7

þ2.6

0.25 0.30

228 241

þ3 þ3

567 600

þ4 þ4

-0.8 -0.9

þ2.4 þ2.2

0.40

256

þ3

651

þ4

-1.2

þ1.6

0.50

255

þ3

682

þ3

-1.4

þ0.2

0.60b

240

þ3

690

þ3

-1.5

0

0.70

214

þ3

686

þ3

-1.6

-0.2

0.80

186

þ2

669

þ3

-1.8

-0.6

0.90

158

þ1

642

þ3

-1.8

-0.6

1.00 1.10

147 158

0c -1

604 563

þ3 þ3

-1.9 -1.9

-0.8 -0.8

1.20

186

-2

525

þ2

-1.9

-1.8

1.30

214

-3

487

þ0.6 c

-1.9

-3.2

1.40

240

-3

460

-2

-1.8

-5.6

1.50

255

-3

429

-5

-1.6

-8.2.

1.60

256

-3

389

-8

-1.2

-10.4

1.70

241

-3

329

-10

-0.9

-11.8

1.75 1.80

228 212

-3 -3

291 255

-11 -12

-0.7 -0.6

-12.4 -13.2

1.85

193

-4

225

-15

-0.5

-16.0

1.90a

177

-5

223

-21

-0.4

-21.8

0

0

0

0

0

0

2.00

Note that for 0.1 e ξ e 1.9, Ret . 30; thus, Υwall = 1 in the outer region for no rotation and for rotation. b For rotation, the intrinsic vorticity is zero over a small region around ξ = 0.6. Also note that the turbulent Reynolds number has a local maximum at ξ = 0.6 and that the shear group has a plateau on the HPW side of the channel. c The shear group is zero on the symmetry plane for no rotation, and for rotation is zero on the LPW side of the channel. a

The preclosure operator A__ (see eq 11 above) for spanwise rotating channel flows can be written as

1 þ NΓ NΩ þ N2Ω A__ ¼

þ 1 þ NΓ NΩ þ N2Ω e x e x þ eyey þ ezez - ðNΓ þ NΩ Þe y e z þ NΩ e z e y

ð40Þ

~ and NΩ 2τRΩx = 2CR1~τRΩ ~. The where NΓ τRΓyz = CR1~τRΓ dimensionless transport time ~τR is defined by eqs 12 and 13. Equation 40 implies that the component equations for the preclosure mapping of B__ into R__ (see eq 9 above) are kR xx ¼ ð1 þ NΓ NΩ þ N2Ω Þ2 Bxx

ð41Þ

kR yy ¼ Byy þ N2Ω Bzz þ 2NΩ Byz

ð42Þ

kR zz ¼ þ Bzz þ ðNΓ þ NΩ Þ2 Byy - 2NΩ Byz - 2NΓ Byz

ð43Þ

kR yz ¼ þ ð1 - ðNΓ þ NΩ ÞNΩ ÞByz - NΓ Byy þ NΩ ðBzz - Byy Þ

ð44Þ

The normalization factor κ is defined as follows k þ ð1 þ NΓ NΩ þ N2Ω Þ2 Bxx þ ð1 þ ðNΓ þ NΩ Þ2 ÞByy þ ð1 þ N2Ω ÞBzz - 2NΓ Byz

ð45Þ

Solutions to eqs 36-39 and eqs 41-44 at different positions in the channel are determined by successive substitution for two cases: δþ = 300 and Ωþ = 0; and δþ = 300 and Ωþ = -0.0042. All fixed point solutions satisfy the convergence criterion ðR__ n þ 1 - R__ n Þ : ei ej e ½R__ n : ei ej 10-5 for i ¼ x, y, or z and j ¼ x, y, or z

ð46Þ

The dimensionless groups ΝΓ(ξ) and NΩ(ξ) are specified by ~ (ξ). Each solution was using DNS results for Γ~(ξ) and for Ω developed by starting with an isotropic operator, i.e., R__ 0 = __I /3. For Ret(ξ) g 30, the local operator K__ DNS is defined as follows ~__ >DNS DNS Γe where < F y z y z z y 3=2

¼ and ~τ DNS R

ð1 þ CR3 NF Þ 3=2

ð1 þ CR2 NF Þ

ð47Þ

DNS ~ In the above equation, Γ~ kDNSΓDNS ; Ω kDNEΩDNE / yz /ε x DNS 2 2 2 ~ þ 4Γ ~Ω ~ þ 8Ω ~. ε ; and NF = Γ

4. RESULTS DNS Temporal Groups. Table 1 shows how the turbulent Reynolds number Ret(ξ) (k/ε)/(ν/k), the shear group ~(ξ) kΩx/ε change Γ~(ξ) kΓyz/ε, and the rotation group Ω in the outer region of the channel for δþ = 300. The three groups were calculated for discrete points ξ( y/δ) in the interval 0.10 e ξ e 1.9 by using previously published DNS results.14,15 The results, which are rounded to one or two significant figures in the table to highlight trends, are used to predict the components of the NR stress in the outer region of the flow field based on the algebraic URAPS equation. The behavior of the URAPS equation in the near wall region will be presented in a subsequent paper. For no rotation (Ωþ= 0), the turbulent Reynolds number ~(ξ) is antisymRet(ξ) is symmetric and the local shear group Γ ~(1) = metric about the symmetry plane. On the symmetry plane, Γ ~(ξ) 0 because Γyz = 0 and k/ε > 0. The shear group Γ ~(0) = 0, changes rapidly in the near wall region with Γ ~(0.03) = 17 (local maximum, not listed in the table), and Γ ~(0.1) = 5. In the “outer” region, the shear group is approximately Γ ~(ξ) = 3 for 0.2 e ξ e 0.7. Near the wall, the constant with Γ turbulent Reynolds number Ret(ξ) increases monotonically from zero at ξ = 0 to Ret = 177 at ξ = 0.1. Ret(ξ) has a local maximum (max Ret = 256) at ξ = 0.4, and a local minimum (min Ret = 147) on the symmetry plane ξ = 1. For spanwise rotation (Ωþ = -0.0042), the shear group Γ~(ξ) is no longer antisymmetric about the symmetry plane. The DNS results in Table 1 show that the peak mean velocity shifts from the symmetry plane (for no rotation) to the LPW side of the channel at ξmax = 1.38 where Γ~(1.38) = 0. For ξ = 0 and ξ = 2, 8910



ARTICLE

Figure 2. Shear component of the NR stress predicted by DNS14,15 (δþ = 300): (a) Ωþ = 0; (b) Ωþ = -0.0042.

the turbulent Reynolds number and the dimensionless shear ~ = 0. On the HPWgroup are zero due to no slip: Ret = 0 and Γ ~(0.03) = 35 (local maximum, not listed in side of the channel, Γ the table) and Γ~(0.1) = 8. On the LPW side of the channel, Γ~(1.95) = -36 (local minimum, not listed in the table) and Γ~(1.9) = -21. For 0.5 e ξ e 1.1, the shear group is ~(ξ) = 3. The turbulent Reynolds approximately constant, Γ number Ret(ξ) increases monotonically from zero on the high pressure wall at ξ = 0 to Ret = 454 at ξ = 0.1; and, from zero on the low pressure wall at ξ = 2 to Ret = 223 at ξ = 1.9. In the outer region, the turbulent Reynolds number has a local maximum at ξ = 0.6 (max Ret = 690). The DNS results show that the turbulent transport time, defined by eq 13 above, does not depend on the turbulent Reynolds number for 0.1 e ξ e 1.9 (i.e., Υwall = 1). ~ þ 2Ω ~) is Table 1 also shows that the intrinsic vorticity (i.e., Γ approximately zero in a small region around ξ = 0.6. The ~(ξ) is zero on the boundaries and has local rotation number Ω extremes (not listed in the table) in the near wall regions (i.e., 0 < ξ < 0.10 and 1.90 < ξ < 2). In the outer region of the flow field, ~ = -0.4 at ξ = 0.1 and ξ = 1.9. It is noteworthy that Table 1 Ω ~ = -1.9 for 1.0 e ξ e 1.3. shows that Ω DNS NR Stress. Figures 2 and 3 show the influence of rotation on the shear component and the normal components of the NRstress in the outer region of the channel (0.1 e ξ e 1.9) based on the DNS results reported by Iwamoto14 and by Wu and Kasagi.15 Table 2 complements Figures 2 and 3 by giving specific numerical results for (δþ, Ωþ) = (300, 0) and for (δþ, Ωþ) = (300, -0.0042).

Figure 3. Normal components of the NR stress predicted by DNS14,15 (δþ = 300): (a) Ωþ = 0; (b) Ωþ = -0.0042; Rxx, Ο; Ryy, •; Rzz, 0).

For Ωþ = 0, the shear component of the NR-stress is zero on the walls and on the symmetry plane: Ryz(0) = Ryz(1) = Ryz(2) = 0 (see Figure 2a and Table 2). Local extremes in the magnitude of the shear stress occur at ξ = 0.40 and at ξ = 1.60 (i.e., |Ryz| = 0.15). In the near wall region, |Ryz |= 0.11 at ξ = 0.10 and at ξ = 1.90. For Ωþ = -0.0042, the shear component of the DNS NR stress is no longer antisymmetric about the symmetry plane (see Figure 2b and Table 2). The magnitude of Ryz near the HPW side of the channel increases a little (i.e., |Ryz| = 0.13); however, the magnitude near the LPW side of the channel decreases by a factor of 2 (i.e., |Ryz| = 0.06). With rotation, the shear stress is zero at ξ = 1.4þ and has local extremes at ξ=0.30 (i.e., |Ryz| = 0.20) and at ξ = 1.80 (i.e., |Ryz| = 0.09). Figure 3 (also see Table 2) shows the influence of rotation on the DNS distribution of turbulent kinetic energy. The upper graph (Figure 3a) shows the results for Ωþ = 0 and the lower graph (Figure 3b) shows the results for Ωþ = -0.0042. In the absence of rotation, all the energy states are in the second sextet of the energy simplex (i.e., 0 e Ryy < Rxx < Rzz < 1). The primary normal stress difference, Rzz - Ryy, is positive for 0 e ξ e 2 . The secondary normal stress difference, Ryy - Rxx, is negative for 0 e ξ e 2. With rotation, the lower graph (see Figure 3b) shows that Rzz - Ryy and Ryy - Rxx are positive and negative, respectively, near the high-pressure and the low-pressure walls. 8911



ARTICLE

Table 2. DNS Components of the NR Stress in the Outer Region for No Rotation14 (δþ = 300 and Ωþ = 0) and for Rotation15 (δþ = 300 and Ωþ = -0.0042) Ωþ = 0 ξ

Rxx

0.10

0.20

0

Ryy

Rzz

Ryz

0.71

-0.11

0 0.09

Ωþ = -0.0042 Rxx

Ryy

0.25

0.10

0

Rzz

Ryz

0.65

-0.13

0

0

0.15

0.24

0.14

0.62

-0.13

0.29

0.15

0.56

-0.15

0.20

0.27

0.16

0.57

-0.14

0.31

0.20

0.49

-0.17

0.25 0.30

0.27 0.28

0.18 0.19

0.55 0.53

-0.14 -0.15

0.31 0.31

0.23 0.27

0.46 0.42

-0.18 -0.19

0.40

0.28

0.12

0.60

-0.15

0.30b

0.33

0.37

-0.20

0.50

0.28

0.21

0.51

-0.15

0.29b

0.38

0.33a

-0.20

0.60

0.28

0.22

0.50

-0.14

0.28b

0.42

0.31a

-0.20

b

0.70

0.27

0.23

0.50

-0.13

0.27

0.44

0.29a

-0.19

0.80

0.27

0.25

0.48

-0.11

0.27b

0.45

0.28a

-0.17

0.90

0.27

0.27

0.46

-0.06

0.28b

0.46

0.26a

-0.16

1.00 1.10

(0.27) 0.27

(0.28) 0.27

(0.45) 0.46

0 þ0.06

0.29b 0.32b

0.45 0.43

0.26a 0.25a

-0.14 -0.11

1.20

0.27

0.25

0.48

þ0.11

0.34b

0.40

0.26a

-0.08

b

-0.03

1.30

0.27

0.23

0.50

þ0.13

0.37

0.34

0.29a

1.40

0.28

0.22

0.50

þ0.14

0.37b

0.27

0.36

þ0.009

1.50

0.28

0.21

0.51

þ0.15

0.35

0.21

0.44

þ0.05

1.60

0.28

0.12

0.60

þ0.15

0.32

0.16

0.53

þ0.07

1.70

0.28

0.19

0.53

þ0.15

0.28

0.13

0.59

þ0.08

1.75 1.80

0.27 0.27

0.18 0.16

0.55 0.57

þ0.14 þ0.14

0.26 0.24

0.11 0.10

0.63 0.66

þ0.09 þ0.09

1.85

0.24

0.14

0.62

þ0.13

0.21

0.08

0.71

þ0.08

1.90

0.20

0.09

0.71

þ0.11

0.15

0.05

0.80

þ0.06

2

0

0

0

0

With rotation, the primary normal stress difference, Rzz-Ryy, is spatially asymmetric; positive near both walls; and, negative for 0.50 e ξ e 1.30. With no rotation, the primary normal stress difference is spatially symmetric; and, non-negative for all spatial positions. b With rotation, the secondary normal stress difference, Ryy-Rxx, is spatially asymmetric; negative near both walls; and, positive for 0.40 e ξ e 1.40. With no rotation, the secondary normal stress differences is spatially symmetric, and nonpositive for all spatial positions. a

However, both normal stress differences flip sign in the outer region on the LPW side of the channel. Clearly, a significant redistribution of energy occurs between the longitudinal component of the fluctuating velocity and the transverse component of the fluctuating velocity. Table 2 also shows that the intrinsic vorticity is approximately zero on the HPW side of the channel. URAPS NR Stress. Figures 4 and 5 show the influence of rotation on the components of the URAPS NR-stress in the outer region of the channel based on the DNS time scales given in Table 1. Table 3 complements Figures 4 and 5 by giving specific numerical results for the components of the URAPS ~(ξ) were used to NR-stress. The DNS groups Γ~(ξ) and Ω determine the preclosure operator A__ for spanwise fully developed channel flows (see eq 40 above). The components of the NR stress were calculated by solving the nonlinear URAPS equation (see component equations given by eqs 36-39 and eqs 41-44 above) at discrete points by using successive substitution (see eq 22 above).

Figure 4. Shear component of the NR stress predicted by URAPS (δþ = 300): (a) Ωþ = 0; (b) Ωþ = -0.0042; see Table 1 for DNS values of ~(ξ). Γ~(ξ)and Ω

~(ξ) and the dimenIn practice, the dimensionless shear group Γ ~(ξ) would be predicted by simultaneously sionless rotation group Ω solving (subject to an appropriate set of boundary conditions) the RANS equation, the nonlinear algebraic URAPS equation, a closed transport equation for the turbulent kinetic energy k, and a closed transport equation for the turbulent dissipation ε. A favorable quantitative comparison of the results with a DNS simulation would validate the full RANS/URAPS model for spanwise rotating channel flows. The objective herein is to present the results of a prevalidation of the URAPS-equation by using DNS results for the shear group and the rotation group. This strategy is independent of the specific closures for the k-ε transport equations. The specific results presented in Figures 4 and 5 (also see Table 3) depend on three phenomenological model parameters, viz., the turbulent transport time τR (defined by eq 13 above) and the two dimensionless “extra” anisotropic prestress coefficients R and β. Because the turbulent Reynolds number is large in the outer region of the flow field (see Table 1) for rotating and for nonrotating turbulent flows, the turbulent transport time de~,Ω ~). The following pends only on the dimensionless group NF(Γ calibration parameters were used in the calculations: CR1 = 0.0036; CR2 = 0.076; and, CR3 = 0.053 (see eq 13 above). For the results presented herein, the “extra” anisotropic coefficients were assumed to be spatially uniform. Therefore, R = 0.1 and β = -0.01, which were estimated based on asymptotic homogeneous shear (see Koppula4). 8912



ARTICLE

Table 3. URAPS Components of the NR Stress in the Outer Region Based on DNS Time Scales with and without Rotation (δþ = 300; Ωþ = 0; and Ωþ = -0.0042; see Table 1 for DNS Groups) Ωþ = 0 ξ

Rxx

Ryy

—

0

—

0.10

0.24

0.18

0.58

0.15

0.26

0.22

0.52

0.20

0.26

0.23

0.25

0.26

0.30

Rxx

Ryy

0

—

0

—

0

-0.19

0.22

0.15

0.63

-0.25

-0.18

0.23b

0.24

0.53

-0.29

0.51

-0.18

0.23b

0.29

0.48

-0.30

0.23

0.51

-0.18

0.23b

0.34

0.43

-0.31

0.26

0.23

0.51

-0.18

0.23b

0.38

0.39

-0.30

0.40 0.50

0.26 0.26

0.23 0.23

0.51 0.51

-0.18 -0.18

0.24b 0.25b

0.45 0.49

0.31a 0.26a

-0.28 -0.23

0.60

0.27

0.23

0.50

-0.18

0.26b

0.51

0.23a

-0.20

b

0

Figure 5. Normal components of the NR stress predicted by URAPS (δþ = 300): (a) Ωþ = 0; () Ωþ = -0.0042 ; Rxx, Ο; Ryy, •; Rzz, 0 ; see ~(ξ). Table 1 for DNS values of Γ~(ξ)and Ω

For Ωþ = 0, the shear component of the NR-stress is zero on the walls and on the symmetry plane: Ryz(0) = Ryz(1) = Ryz(2) = 0(see Figure 4a and Table 3). In the near wall region, |Ryz| = 0.19 at ξ = 0.10. This prediction is almost a factor of 2 larger than the DNS result of |Ryz|DNS = 0.11 at ξ = 0.10. This occurs because the Cayley-Hamilton (CH) coefficients R and β were estimated by using asymptotic homogeneous shear data. Clearly, an improved quantitative comparison with nonrotating DNS results would occur by selecting fully developed channel flows as the benchmark flow. Nevertheless, it is most significant that the URAPS equation based on a calibration using unbounded flow statistics produces a realizable NR stress for all values of ~(ξ). Γ~(ξ) and Ω þ For Ω = -0.0042, the shear component of the NR-stress is no longer antisymmetric about the symmetry plane (see Figure 4b and Table 3). Near the HPW side of the channel, Ryz = -0.25 . Near the LPW side of the channel, Ryz = þ0.10. The URAPS-closure predicts that the shear stress is zero at ξ = 1.38, which is consistent with the DNS results portrayed in Figure 2b. Figure 4b (and Table 3) also shows that the URAPS closure predicts that rotation causes |Ryz|max to increase from |Ryz|max = 0.19 at ξ = 0.1 (for no rotation) to |Ryz|max = 0.31 at ξ = 0.25 (for rotation), which represents an increase of 63% in

Rzz

Ωþ = -0.0042 Ryz

Rzz

Ryz

0.70

0.27

0.24

0.49

-0.17

0.26

0.52

0.22a

-0.18

0.80

0.28

0.27

0.45

-0.16

0.27b

0.52

0.21a

-0.17

0.90

0.31

0.31

0.38

-0.11

0.27b

0.52

0.21a

-0.14

1.00

1/3

1/3

1/3

0

0.28b

0.51

0.21a

-0.13

1.10

0.31

0.31

0.38

-0.11

0.29b

0.48

0.23a

-0.09

1.20 1.30

0.28 0.27

0.27 0.24

0.45 0.49

þ0.16 þ0.17

0.32b 0.33b

0.43 0.36

0.25a 0.31a

-0.05 -0.01

1.40

0.27

0.23

0.50

þ0.18

0.33

0.28

0.39

þ0.02

1.50

0.26

0.23

0.51

þ0.18

0.30

0.21

0.49

þ0.05

1.60

0.26

0.23

0.51

þ0.18

0.27

0.16

0.57

þ0.07

1.70

0.26

0.23

0.51

þ0.18

0.25

0.13

0.62

þ0.09

1.75

0.26

0.23

0.51

þ0.18

0.24

0.12

0.64

þ0.10

1.80

0.26

0.23

0.51

þ0.18

0.24

0.10

0.66

þ0.11

1.85 1.90

0.26 0.24

0.22 0.18

0.52 0.58

þ0.18 þ0.19

0.23 0.22

0.08 0.06

0.69 0.72

þ0.11 þ0.10

0

—

0

—

—

0

2

—

0

With rotation, the primary normal stress difference, Rzz-Ryy, is spatially asymmetric; positive near both walls; and, negative for 0.40 e ξ e 1.30. With no rotation, the primary normal stress difference is spatially symmetric; and, non-negative for all spatial positions. b With rotation, the secondary normal stress difference, Ryy-Rxx, is spatially asymmetric; negative near both walls; and, positive for 0.15 e ξ e 1.30. With no rotation, the secondary normal stress difference is spatially symmetric, and, nonpositive for all spatial positions. a

the shear component of the NR stress. The DNS results show a 33% increase in the maximum magnitude of the shear stress (see Table 2). On the low pressure side of the channel, Figure 4a (see Table 3) shows that the URAPS closure predicts that rotation significantly suppresses the shear component of the NR stress. (i.e., ξ = 1.80, |Ryz| = 0.11 with rotation, and |Ryz| = 0.18 without rotation, a 39% decrease in the shear stress). The DNS results predict a 45% decrease in the shear stress at ξ = 1.80. Finally, it is notable that the shear component of the NR stress is zero near ξ = 1.4 for the DNS results (see Figure 2b) and for the URAPS results (see Figure 4b). Figure 5a (also see Table 3) shows the distribution of turbulent kinetic energy predicted by the URAPS-closure in the absence of rotation (i.e., Ωþ = 0). All of the energy states are in the second sextet within the positive orthant of the energy simplex (i.e., 0 < Ryy < Rxx < Rzz < 1). The primary normal stress difference, Rzz - Ryy, is positive for 0 e ξ e 2 and the secondary 8913



ARTICLE

normal stress difference, Ryy - Rxx, is negative for 0 e ξ e 2. These predictions are qualitatively consistent with the DNS results summarized by Figure 3a (and Table 2). For no rotation, the DNS results show that the energy state on the symmetry plane is anisotropic (i.e., Ryy = 0.27 < Rxx = 0.46), whereas the URAPS results show that Ryy = Rxx = Rzz = 1/3 at ξ = 0. Note, however, that the URAPS energy states just off the symmetry plane at ξ = 0.9 are qualitatively similar to the DNS results inasmuch as Ryy = Rxx = 0.31 < Rzz = 0.38 for the URAPS closure (see Tables 2 and 3). Near the symmetry plane, the turbulent energy is transferred equally from the longitudinal component of the fluctuating velocity to the other two components. With rotation (i.e., Ωþ = -0.0042), the URAPS closure predicts that Rzz - Ryy > 0 and Ryy - Rxx < 0 on the HPWside of the channel (i.e., 0 < ξ < 0.5) and on the LPW-side of the channel (i.e., 1.3 < ξ < 2; see Figure 5b). However, for 0.5 e ξ e 1.3, significant energy is transferred from the longitudinal component of the fluctuating velocity to the transverse component of the fluctuating velocity. In this region, the signs of the primary and secondary normal stress differences flip. It is notable that the NR stress is nearly isotropic for ξ = 1.4 (see Figure 3b and Table 2 for the DNS results, and Figure 5b and Table 3 for the URAPS results).

5. CONCLUSIONS URAPS Closure. The URAPS equation, defined by eqs 9 and 14 above, is a nonlinear dyadic-valued mapping of R__ into itself that depends on the mean velocity gradientr, the frame rotation operator Ω__ , and a turbulent transport time τR

R__ ¼ R__ ðR__ , K__ ; R, βÞ where K__ τR < F__ > and < F__ > r < u_ > þ 2Ω__

ð48Þ

)

)

The coefficients, R and β, depend on the eigenvalues of R__ . The turbulent time scale τR depends on the Reynolds number Ret( k2/(νε)) and the shear group NF( k /ε). The turbulent kinetic energy k and the turbulent dissipation ε are governed by the following transport equations Dk þ 3 rk ¼ Dt þ r 3 ½ðνI__ þ τR < u0 u0 > Þ 3 rkÞ - < u0 u0 > : < S__ > - ε

ð49Þ

Dε þ 3 rε ¼ Dt h

i þ r 3 νI__ þ τR Þ 3 rε þ CP

-: ε - CD τR τR

ð50Þ

The “production” coefficient CP and the dissipation coefficient CD also depend on Ret and NF. Equations 48, 49, and 50 define the theoretical framework of the universal realizable anisotropic prestress (URAPS) closure for the NR stress. Benchmark flows are used to determine the functional dependence of R and β on IIR and IIIR. Solutions to eq 48 are non-negative operators, provided R and β satisfy the

universal ineqs 18 and 19. The phenomenological coefficients CP and CD are selected so that solutions to eqs 49 and 50 are positive within the flow domain. Benchmark flows are also used to determine the functional dependence of τR, CP, and CD on Ret and NF. Continuity and the no slip boundary condition imply that the turbulent kinetic energy and its normal gradient at a solid/fluid interface are zero. Equations 48, 49, and 50 provide a closure for the RANS equation for all rigid body frames-of-reference (see eq 1). The NR-stress occurs explicitly in all four equations. Solutions to eq 48 are realizable for all K__ ; however, the accuracy of predictions for a wide class of flows depends on selecting benchmark flows that cover a wide range of scales (i.e., 0 e Ret < ¥, 0 e NF < ¥, 0 < IIR < 1, and 0 < IIIR < 1). Redistribution of Turbulent Kinectic Energy. A comprehensive validation of any closure model for the RANS equation is an incredibly important step in turbulence modeling. Section 4 above gives a summary of the results of a prevalidation study of eq 48 for R = þ0.1 and β = -0.01. Previously published DNS results for spanwise rotating, fully developed channel flows were used to support the prevalidation study. The spatial distributions of the mean velocity gradient, ΓDNS yz (ξ), the turbulent kinetic energy, kDNS(ξ), and the turbulent dissipation, ~DNS(ξ), and (ξ), Γ εDNS(ξ), were used to calculate ReDNS t DNS þ þ ~ Ω (ξ) for δ = 300 in an inertial frame (Ω = 0) and in a noninertial frame (Ωþ = -0.0042). Table 1 summarizes the DNS DNS . Table 3 gives a statistics used to calculate K__ DNS τDNS __ > R : __ 3

ð51Þ

For fully developed channel flows (either rotating or nonrotating), eq 51 implies that the turbulent kinetic energy is equally distributed among the three components of the fluctuating velocity for all positions in the flow domain: Rxx = Ryy = Rzz = 1/3 for 0 < ξ < 1. This prediction is clearly inconsistent with the underlying physical principles associated with turbulence. By contrast, the URAPS equation, which stems directly from the Navier-Stokes equation, provides a clear and physically unambiguous closure for the NR-stress in inertial and noninertial frames-of-reference. Although the URAPS equation is a nonlinear algebraic equation for the NR-stress, the solution to this equation can be determined by the method of successive substitution. Thus, the computational burden required to determine the mean velocity field and the mean pressure field based on the RANS equation and the continuity equation (i.e., eqs 1 and 4,

respectively) with the URAPS closure (i.e., eqs 48, 49, and 50) is comparable to the computational burden required for an “eddy” viscosity model. Most significantly, however, is the assurance that the URAPS closure will always produce realizable solutions for all inertial and noninertial frames-of-reference.

’ AUTHOR INFORMATION Corresponding Author

*Phone: (517) 353-5486. Fax: (517) 432-1105. E-mail: petty@ msu.edu.

’ ACKNOWLEDGMENT This research was partly funded by the National Science Foundation Industry/University Cooperative Research Center for Multiphase Transport Phenomena (NSF/ECC 0331977; NSF/IIP 093474). The authors gratefully acknowledge the members of the Center for their support. The authors sincerely thank Mr. Satish Muthu for preparing Figures 2-5.

8915

) )

’ NOMENCLATURE A__ preclosure operator (see eq 11) normalized prestress operator (see eq 14) B__ extra anisotropy prestress coefficients, i = 1, 2, 3 Ci (see eq 14) closure parameters, i = 1, 2, 3, 4 (see eqs 12 and 13) CRi e x, e y , e z Euclidean base vectors E3 Euclidean vector space mean field operator, 1/s f0 fluctuating acceleration (see eq 10), m/s2 g acceleration due to gravity, m/s2 unit dyadic-valued operator __I IR first invariant of the NR-stress, tr(R__ ) IIR second invariant of the NR-stress, tr(R__ 3 R__ ) IIIR third invariant of the NR-stress, tr(R__ 3 R__ 3 R__ ) k turbulent kinetic energy, tr/2, m2/s2 K__ kinematic operator (see eq 11) K__ DNS kinematic operator based on DNS statistics NR dimensionless time, τR NF dimensionless time, k /ε NΓ dimensionless time scale, τRΓyz NΩ dimensionless time scale, 2τRΩx p,p0 instantaneous and fluctuating pressure fields, N/m2
Reynolds average pressure, N/m2 mean hydrodynamic pressure (see eq 25), N/m2 Pk “production” of turbulent kinetic energy, m2/s3 QR quadratic form associated with R__ (see eq 7) R__ normalized Reynolds (NR) stress, /tr() Ret turbulent Reynolds number, k2/(νε) mean strain rate, [r þ (r)T]/2, 1/s t time, s mean velocity, m/s u0 fluctuating velocity, m/s u* friction velocity (see eq 32), m/s ub cross sectional average velocity, m/s uþ mean velocity scaled with u* ) )

) )

(4) It is noteworthy that in the presence of rotation (see Figure 3b for DNS and Figure 5b for URAPS) Rxx = Ryy = Rzz = 1/3 at ξ = 1.38. Also, note that Figure 2b (DNS) and Figure 4b (URAPS) show that Ryz = 0 at ξ = 1.38 for both DNS and URAPS. The above observations support the conclusion that the current calibration of the URAPS equation captures the underlying physics associated with the convective coupling between the Coriolis acceleration and the instantaneous velocity field. This phenomenon is responsible for breaking the symmetry of the mean velocity field and for the redistribution of the turbulent kinetic energy among the three components of the fluctuating velocity. This is a significant conclusion based on the results of the prevalidation analysis presented in Section 4 above. An improved quantitative agreement between the DNS and URAPS statistics for rotating and nonrotating channel flows is expected by using “extra” anisotropic coefficients (i.e., R and β) based on benchmark flows that cover a broader set of anisotropic turbulent states (i.e., 0 < IIR < 1, and 0 < IIIR < 1). “Production” of Turbulent Kinetic Energy. As noted in Section 3 above, the “production” of turbulent kinetic energy in eq 49 (i.e., Pk -:r) may be either positive or negative.11,25 The no-slip boundary condition implies that Pk is zero at a solid/fluid interface. For Ωþ = -0.0042, the DNS results for the shear component (see Figure 2b and Table 2) together with the results listed in Table 1 support the idea that the shear component of the NR stress and the mean velocity gradient are zero at different locations on the low pressure side ~DNS(ξmax) = 0; ξo 6¼ ξmax. of the channel: RDNS yz (ξo) = 0; Γ Grundestam et al.23,24 have published DNS results that also show this behavior at higher values of |Ωþ|. Thus, DNS statistics for transverse rotating channel flows seem to support the idea that the “production” of turbulent kinetic energy is positive near the solid/fluid interface and negative in the outer region on the low pressure side of the channel. As noted by Churchill,25 there is no a priori theoretical requirement that the turbulent “production” of turbulent kinetic energy be non-negative. As noted above, the current formulation of the URAPS equation does not predict a negative “production” of the turbulent kinetic energy. Comparison with Other Algebraic Closure Models. If K__ NR , 1 and if the prestress operator B__ is isotropic (i.e., B__ = __I /3) for all positions in the flow domain, then eq 48 reduces to an “eddy” viscosity type model7


Industrial & Engineering Chemistry Research - þ x (x, y, z) yþ z

kinematic Reynolds stress, m2/s2 kinematic Reynolds momentum flux, m2/s2 space-time velocity correlation, m2/s2 mean vorticity, 1/s position vector, m spanwise, transverse, and longitudinal directions, m transverse coordinate scaled with ν/u* arbitrary unit vector in E3

Greek symbols

R(IIR,IIIR) β(IIR,IIIR) δ δþ ε _ε_ φ κ Δp Γyz ~ Γ Γþ λRi λBi μ ν F ξ ξmax ξ0 τE τR ~τR τHPW τLPW Υwall Ω ~ Ω Ωþ Ω_

“extra” anisotropic prestress coefficient (see eq 16) “extra” anisotropic prestress coefficient (see eq 15) half-width of channel (see Figure 1), m half-width scaled with ν/u* turbulent dissipation, ν, m2/s3 permutation triadic ratio of wall shear stresses (see eq 33) normalization factor (see eq 45) relative pressure (see eq 28), N/m2 shear component of the mean velocity gradient, 1/s Γyz scaled with k/ε Γyz scaled with ν/(u*)2 eigenvalue of R_ eigenvalue of B_ fluid viscosity, Pa-s fluid kinematic viscosity, m2/s fluid density, kg/m3 transverse coordinate scaled with δ position where Γþ(ξmax) = 0 position where Ryz(ξ0) = 0 “eddy” transport time (see eq 51), s URAPS transport time (see eq 12), s URAPS transport time scaled with k/ε wall shear stress on the high pressure wall (see Figure 1), N/m2 wall shear stress on the low pressure wall (see Figure 1), N/m2 empirical wall function (see eq 13) angular velocity of reference frame, Ωxex, 1/s Ωx scaled with k/ε Ωx scaled with ν/(u*)2 rotational (or Coriolis) operator, _ε_ 3 Ω_ , [=] 1/s

Operators

a3b A__ 3 B__ A__ T 3 B__ T tr(A__ ) A__ : B__

Reynolds average operator scalar product between two vectors (scalar) inner scalar product between two dyadics (dyadic) outer scalar product between two dyadics (dyadic) trace of the operator A__ (scalar) trace of the operator A__ 3 B__ (scalar)

ARTICLE

(3) Piquet, J. Turbulent Flows: Models and Physics; Springer: New York, 1999. (4) Koppula, K. S. Universal Realizable Anisotropic Prestress Closure for the Normalized Reynolds Stress. Ph.D. Dissertation; Michigan State University, East Lansing, MI, 2009. (5) Schumann, U. Realizability of Reynolds Stress Turbulence Models. Phys. Fluids 1977, 20 (5), 721. (6) Lumley, J. Computational Modeling of Turbulent Flows. Adv. Appl. Mech. 1978, 18, 123. (7) Parks, S. M.; Weispfennig, K.; Petty, C. A. An Algebraic Preclosure Theory for the Reynolds Stress. Phys. Fluids 1998, 10 (3), 645. (8) Koppula, K. S.; Benard, A.; Petty, C. A. Realizable Algebraic Reynolds Stress Closure. Chem. Eng. Sci. 2009, 64, 4611. (9) Monin, A. S.; Yaglom, A. M. Statistical Fluid Mechanics: Mechanics of Turbulence; MIT Press: Boston, MA, 1965; Vol. 1. (10) Hanjalic, K. Advanced Turbulence Closure Models: A View of Current Status and Future Prospects. Int. J. Heat Fluid Flow 1994, 15 (3), 178. (11) Pope, S. B. Turbulent Flows; Cambridge University Press: Cambridge, U.K., 2000. (12) Gatski, T. B. Constitutive Equations for Turbulent Flows. Theoret. Comput. Fluid Dyn. 2004, 18, 345. (13) El-Samni, O. A.; Kasagi, N. The Effect of System Rotation with Three Orthogonal Rotating Axes on Turbulent Channel Flow. Proceedings of the 7th International Conference on Fluid Dynamics and Propulsion; Cairo, Egypt, Dec 19-21, 2001; ASME International: New York, 2001. (14) Iwamoto, K. Database of Fully Developed Channel Flow. THTLAB Internal Report No. ILR-0201. THTLAB, The University of Tokyo: Tokyo, 2002. (15) Wu, H.; Kasagi, N. Effects of Combined Streamwise and Spanwise Rotations on Turbulent Channel Flow. The Fifth JSME-KSME Fluids Engineering Conference; Nagoya, Japan, Nov 17-21, 2002; Japan Society of Mechanical EngineersJ: Tokyo, 2002. (16) Weispfennig, K.; Parks, S. M.; Petty, C. A. Isotropic Prestress for Fully Developed Channel Flows. Phys. Fluids 1999, 11 (5), 1262. (17) Parks, S. M.; Petty, C. A. Anisotropic Prestress Theory for Homogeneously Sheared Turbulent Flows. In Turbulent Shear Flow Phenomena; Banerjee, S., Eaton, E. K., Eds.; Begell House: New York, 1999; p 227. (18) Weispfennig, K. Relaxation/Retardation Model for Fully Developed Turbulent Channel Flow. Ph.D. Dissertation; Michigan State University, East Lansing, MI, 1997. (19) Park, J. Y.; Chung, M. K. A Model for the Decay of Rotating Homogeneous Turbulence. Phys. Fluids 1999, 11 (6), 1544. (20) Boussinesq, J. Theorie de l’Ecoulement Tourbillant. Memoires Presentes Par Divers Savants a l’Academie des Sciences, 1877, 23, 46. (21) Speziale, C. G. Some Important Remarks Concerning Recent Work on Rotating Turbulence. Phys. Fluids 1998, 10 (12), 3242. (22) Speziale, C. G.; Younis, B. A.; Rubinstein, R; Zhou, Y. On Consistency Conditions for Rotating Turbulent Flows. Phys. Fluids 1998, 10 (8), 2108. (23) Grundestam, O.; Wallin, S.; Johansson, A. V. Direct numerical simulations of rotating turbulent channel flows. J. Fluid Mech. 2008, 598, 177. (24) Grundestam, O.; Wallin, S.; Johansson, A. V. A priori evaluations and least-squares optimizations of turbulence models for fully developed rotating channel flow. Eur. J. Mech., B/Fluids 2008, 27, 75. (25) Churchill, S. W. A Critique of Predictive and Correlative Models for Turbulent Flow and Convection. Ind. Eng. Chem. Res. 1996, 35, 3122.

’ REFERENCES (1) Greenspan, H. P. Theory of Rotating Fluids; Cambridge University Press: Cambridge, U.K., 1968. (2) Parks, S. M. Relaxation Model for Homogeneous Turbulent Flows. Ph.D. Dissertation; Michigan State University, East Lansing, MI, 1997. 8916

Turbulent Energy Redistribution in Spanwise Rotating Channel Flows

Recommend Documents