A quadratic constitutive equation for the turbulent Kolmogorov flow
Wenwei Wu, Francois G. Schmitt, Enrico Calzavarini, Lipo Wang
AA quadratic constitutive equation for the turbulent Kolmogorov flow
Wenwei Wu,
1, 2, ∗ Fran¸cois G. Schmitt, † Enrico Calzavarini, and Lipo Wang Univ. Lille, CNRS, ULCO, Laboratory of Oceanology and Geosciences, UMR LOG 8187, Wimereux, France UM-SJTU Joint Institute, Shanghai JiaoTong University 200240, Shanghai, People’s Republic of China Univ. Lille, Unit´e de M´ecanique de Lille - J. Boussinesq (UML) ULR 7512, F59000 Lille, France (Dated: January 14, 2021)We study the three-dimensional turbulent Kolmogorov flow, i.e. the Navier-Stokes equationsforced by a low-single-wave-number sinusoidal force in a periodic domain, by means of direct nu-merical simulations. This classical model system is a realization of anisotropic and non-homogeneoushydrodynamic turbulence. Boussinesq’s eddy viscosity linear relation is checked and found to beapproximately valid over half of the system volume. A more general nonlinear quadratic constitu-tive equation is proposed and its parameters estimated at varying the Taylor scale-based Reynoldsnumber in the flow up to the value 200. This provides a Reynolds number-dependent quadraticclosure for the Kolmogorov flow. The case of a forcing with a different shape, here chosen Gaussian,is considered and the differences with the sinusoidal forcing are emphasized.
I. INTRODUCTION
In the late 1950s, A.N. Kolmogorov proposed to study thestability properties of an incompressible flow described bythe Navier-Stokes equation forced by a sinusoidal shearforce in a periodic domain. An answer was put for-ward soon after [1], indicating the existence of a criticalReynolds number R = √
2, confirmed also later in [2](we will provide below the definition of Reynolds num-ber which was used in these works). Such a model sys-tem, since then dubbed Kolmogorov flow (KF), is notstraightforward to be realized in experiments. However,an important work published in the Russian literature[3] and discussed by A. N. Obukhov [4] describes an ex-periment using a thin layer of electrolyte in an exter-nal force field capable to generate an analogous flow,for which the stability properties as well as the transi-tion to the turbulent state were studied. The resultsof this experiment, supported by a previous theoreticalwork [5], have been taken as an evidence that there ex-ists in the KF a succession of instabilities with increas-ing Reynolds numbers, until reaching a fully turbulentstate for Reynolds numbers of the order of 1000 [3, 4].One the other hand, since the advent of computers andin particular since the introduction of the Fast FourierTransform (FFT) algorithm the Kolmogorov flow has be-come amenable to be explored via numerical simulations,even in high-Reynolds number conditions [6]. This makesthe turbulent Kolmogorov flow (TKF) possibly the sim-plest and most accessible prototype of open flow, i.e. aflow without boundaries, which is at the same time sta-tistically stationary, anisotropic, and non-homogeneous(along one direction) [6–11]. As discussed by Musacchioand Boffetta [11], the TKF can be considered, to somerespect, as a turbulent channel flow without boundaries.In recent years, the TKF has been mainly studied the- ∗ wenwei˙[email protected] † [email protected] oretically and numerically. The large-scale forcing wasoriginally proposed as a sinusoidal force, which, however,is not a requirement, and rather constitutes a convenientsimplification for the theoretical analysis and for numer-ical implementations (which are mostly based on FFT).Other shapes for the large-scale forcing could be imag-ined as well [12]. Many numerical studies devoted on KFand TKF have adopted a two-dimensional configurations[13–16]. In this work we focus on the three-dimensionalcase, i.e. fully resolved Navier-Stokes incompressible tur-bulence forced along the x direction by a large-scale sinu-soidal force depending solely on the z coordinate. As ob-served by Borue and Orszag [6], such flow is a convenienttest ground for transport models. Such a considerationmotivates the present study.The structure of this article is as follows: After a sectionintroducing the present notations and numerical details,the theoretical framework of linear and nonlinear closureequation for the Reynolds stress tensor is presented andadapted to the TKF model system. First, we considerthe classical turbulence closure based on eddy-viscosityBoussinesq’s approach, where the traceless stress tensoris assumed to be proportional to the mean strain-ratetensor. Such closure is at the basis of many turbulencemodels including k - (cid:15) , k - ω and all eddy viscosity trans-port models [17]. We show the range of applicability andlimitations of this assumption in the context of TKF.Second, a nonlinear quadratic closure that make use oftensor invariants is directly tested on TKF. This is per-formed along the lines of previous direct test done forchannel flows [18, 19], and for various Reynolds num-bers. In a following section a model flow system with adifferent forcing, non-sinusoidal shape, is considered andits differences with the original KF forcing are consid-ered. The last section is devoted to a discussion of themain findings and conclusions. a r X i v : . [ phy s i c s . f l u - dyn ] J a n II. THE KOLMOGOROV FLOW MODELSYSTEMA. Equations of motion and numericalimplementations
The governing equations for velocity field u ( x , t ) are theincompressible Navier-Stokes equations, ∂ u ∂t + ( u · ∇ ) u = − ρ ∇ p + ν ∆ u + f , (1)where p is the hydrodynamic pressure, ρ the fluid densityand ν is the kinematic viscosity. This flow is sustainedby a constant in time and spatially dependent force f ofthe form: f = A sin (cid:16) π zH (cid:17) e x , (2)where A is a constant, H is the length of the side of thecubic domain, chosen here as the characteristic lengthscale. Such a force, directed along the x direction anddepending only on the z coordinate, makes the turbulentflow statistically anisotropic and non-homogeneous alongthe z direction (but statistically homogeneous in the x - y planes). It is convenient to introduce the followingreference scales for velocity and time: U = ( AH ) / ; (3) T = HU = (cid:18) HA (cid:19) / . (4)From this, one can construct the Reynolds number as Re = U Hν , (5)which thus becomes the only dimensionless control pa-rameter in the system. Let us mention that in sta-bility analyses [2, 5] another Reynolds number is used,in the present notation reads R = AH / ( ν (2 π ) ). Itthus yields R = Re / ((2 π ) ). The stability criterion R > √ Re > / (2 π ) / (cid:39) .
7. In the fol-lowing we will also use the Reynolds number based onTaylor microscale, Re λ = λu (cid:48) /ν , where λ = u (cid:48) (cid:112) ν/(cid:15) , ν is the kinematic viscosity, u (cid:48) = (cid:112) u (cid:48) is the globalroot-mean square of single component velocity, (cid:15) = ν (cid:80) i (cid:80) j ( ∂ i u j + ∂ j u i ) is the global energy dissipationrate, and the overbar · · · denotes the global average (intime and all over the spatial domain). The latter numberis more convenient to quantify the degree of turbulencerealized in the system. In the rest of this article, all thereported quantities are dimensionless with reference tothe units defined in (3) and (4).The Kolmogorov flow model system is numerically sim-ulated in a cubic tri-periodic domain. The dynami-cal equations (1) are solved numerically by means ofa pseudo-spectral code using a smooth dealiasing tech- nique for the treatment of non-linear terms in the equa-tions [20]. The spatial resolution is chosen such that thecondition | (cid:126)k | max · η >
1, where | (cid:126)k | max ≈ . N is themaximum wave number amplitude kept by the dealias-ing procedure, is always satisfied. The time-marchingscheme adopts a third order Runge-Kutta method. Theglobal non-dimensional values of the key parameters forthe simulations are reported in table I. Two criteria havebeen proposed for the convergence of Kolmogorov flowsimulations [21]: first the mean energy injection shouldbe equal (within numerically accuracy) to the total dissi-pation, and second the left-hand-side and right-hand-sideof equation (8) must be equal. It has been checked herethat these two criteria are satisfied in our simulations(the right-hand-side of equation (8) is shown in figure2(b)). The total integration time is chosen in such a wayto have comparable datasets for each run and to ensurethe statistical convergence of the measurements (see Ta-ble I). B. Reynolds decomposition and velocity moments
Let us consider a Reynolds decomposition of the velocityinto mean and fluctuating quantities u = (cid:104) u (cid:105) + u (cid:48) ( (cid:104)·(cid:105) denotes the average over time and spatially along x and y directions) and note the three cartesian components ofthe velocity u = ( u, v, w ). Because of the periodicity in x and y directions, the derivatives with respect to x and y of mean quantities are 0. By taking the average (cid:104)·(cid:105) of the Navier-Stokes equations, one obtains the followingrelations: − ∂ z τ = 1 Re ∂ z U ( z ) + sin(2 πz ) ,∂ z (cid:104) w (cid:48) (cid:105) = − ∂ z (cid:104) p (cid:105) , (6)where U = (cid:104) u (cid:105) = (cid:104) u (cid:105) and the shear stress is τ = −(cid:104) u (cid:48) w (cid:48) (cid:105) .Using the first line of this relation, one can verify thatfor laminar flows when τ = 0, the mean velocity profileis sinusoidal while the pressure field is constant [1].For turbulent flows, it is well-known that the mean ve-locity profile is also sinusoidal [6]. However, this is anumerical result which has so far, to our knowledge, nodirect analytical explanation. We obtain the following z -dependence for U : U ( z ) = κ sin(2 πz ) . (7)where κ is a coefficient whose numerical estimation, atvarying the Reynolds number, is plotted in figure 1(b).The maximum value of the mean turbulent velocity isof the order of the characteristic velocity built using theforcing values, since we obtain values of κ between 1.01and 1.12, increasing with the Reynolds number (figure1(a) and Table II). The values of κ found here are com-patible with the value of κ = 1 . Re Re λ Re ∗ ν (cid:15) η N | (cid:126)k | max · η T total /T l ∆ t Re = HU ν , Re λ = λu (cid:48) ν the Taylor-scale based number, Re ∗ = UH πν = HκU πν (same as in Ref [11]); the kinematic viscosity ν ; (cid:15) = ν (cid:80) i (cid:80) j ( ∂ i u j + ∂ j u i ) is the global energy dissipation rate, · denotes the global average (in time and all overthe spatial domain). N is the grid size; η = ( ν /(cid:15) ) / is the global Kolmogorov scale; | (cid:126)k | max · η is the spatialresolution condition, where | (cid:126)k | max ≈ . N is the maximum wave number amplitude kept by the dealiasingprocedure; T total is the total simulation time and T l is the large eddy turnover time [17], which thus implies that T total /T l denotes the number of large eddy turnover time spanned by the simulation in statistically steadyconditions; ∆ t is the numerical time step.cosity of the 8th order and the value of the Reynoldsnumber is not provided). In the work by Musacchio andBoffetta [11], the dependence of the friction coefficient f (written as f = AH πκ in the present notation), on theReynolds number based on the forcing scale and meanvelocity (denoted here Re ∗ , which writes Re ∗ = κ π Re in our notation) was investigated. Such dependence isalso plotted in the inset of figure 1(a): the friction coeffi-cient values obtained in our simulations are comparablewith those reported by Musacchio and Boffetta [11] inthe same range of Reynolds numbers.For large Reynolds numbers, using equations (6) and (7)we find that ∂ z τ is proportional to sin(2 πz ), obtainingfinally: τ = 12 π (cid:16) − (2 π ) κRe (cid:17) cos(2 πz ) . (8)The first and second moments of the velocity obtainedafter averaging Navier-Stokes equations are shown in fig-ures 2 a) and b). We observe that only one componentof the mean velocity is non-zero; concerning second mo-ments, only the shear stress term (cid:104) u (cid:48) w (cid:48) (cid:105) is non-zero. Theturbulence is globally anisotropic since all normal stresscomponents of the stress tensor are different. Specifically, (cid:104) u (cid:48) (cid:105) > (cid:104) w (cid:48) (cid:105) > (cid:104) v (cid:48) (cid:105) (see figure 3). The diagonal termshave twice the spatial frequency of the forcing. Since cos(2 θ ) = 2 cos θ −
1, they can be written as: (cid:104) u (cid:48) (cid:105) = α + β cos (2 πz ) , (cid:104) v (cid:48) (cid:105) = α + β cos (2 πz ) , (cid:104) w (cid:48) (cid:105) = α + β cos (2 πz ) , (9)where ( α i , β i ) are numerical coefficients. The estimatedvalues of such coefficients for the different runs are listedin table II, which will be important for the quadraticclosure done in the next section. Consequently, we canwrite also the evolution of the mean kinetic energy: K ( z ) = α + β cos (2 πz ) . (10)The mean kinetic energy profiles are presented in figure4, from which we can see that the increase of Reynoldsnumber leads to the global growth of the kinetic energy.The coefficients α and β as functions of Re λ are plottedin figure 1(b). The values found for the highest Reynoldsnumber values are in good agreement with the valuesreported in [6] ( α = 0 .
391 and β = 0 . III. CLOSURES FOR THE TURBULENTKOLMOGOROV FLOW
In this section we aim at deriving a relation express-ing the Reynolds stress in terms of the gradients of themean velocity flow. Such a relation, also known as tur-No. α β α β α β κ . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . T = −(cid:104) u (cid:48) ⊗ u (cid:48) (cid:105) (with ⊗ denoting the dyadic product).The anisotropic stress tensor is R = − T + K I , where K is the kinematic energy and I is the identity tensor.The mean velocity gradient tensor A = ∂ (cid:104) u i (cid:105) /∂x j , andthe mean strain-rate S and rotation-rate W tensors arealso introduced as: S = 12 (cid:18) ∂ (cid:104) u i (cid:105) ∂x j + ∂ (cid:104) u j (cid:105) ∂x i (cid:19) , (11) W = A − S . (12)A closure for the turbulence equations corresponds toexpressing the Reynolds stress tensor using mean quan-tities, e.g. when the closure is local, using the tensors S and W . Below we first consider the simplest linearclosure and estimate the eddy-viscosity, and later on weaddress a nonlinear closure using a quadratic constitutiveequation. A. Boussinesq’s eddy-viscosity hypothesis and itsassessment
It is seen from equations (7) and (8) that the only non-zero non-diagonal term in the stress tensor has the same z -dependence as the mean gradient term. This leads toan eddy-viscosity of the form: ν T = τU (cid:48) ( z ) = (cid:18) π ) κ − Re (cid:19) . (13)The eddy-viscosity does not depend on z , but depend onthe Reynolds number and the coefficient κ . For Run 7 wefind a value of ν T = 0 . z . This implies that ν T /ν = O (10 ), in agree-ment with [11] results at comparable Reynolds number. However, the estimation of an eddy-viscosity does notvalidate the linear closure. The Boussinesq’s hypothe-sis, which is at the basis of all eddy-viscosity turbulencemodels, corresponds to a linear proportionality betweentensors [22] : R = 2 ν T S . (14)For the flow considered here, there are some symmetriesso that the strain as well as the stress have a simplifiedform: S = a (15)and R = K − σ u τ K − σ v τ K − σ w , (16)where a = U (cid:48) ( z ), σ u = (cid:104) u (cid:48) (cid:105) and the same for σ v and σ w .It is then clear, as also the case for turbulent channelflows [17, 23, 24], that such linear relation between ten-sors can be realized only when diagonal terms are zero,i.e. in an isotropic situation. However, the Kolmogorovflow is anisotropic and as seen in figure 3, the three nor-mal stresses are all different, which means that a preciseproportionality does not exist.In order to quantify this alignment, we consider the innerproduct between tensors: A : B = { A t B } = A ij B ij ,where { X } is a notation for the trace of X . The normis then || A || = A : A . As a direct test of Boussinesq’shypothesis, we first represent here the normalized innerproduct of R and S tensors (which is similar to the cosine Re Re * I U L F W L R Q F R H I I L F L H Q W AH 0 % Re F R H I I L F L H Q W V Figure 1: (a) Amplitude of mean velocity profile, κ inEq. 7, as function of Re λ . The inset panel shows thedependence of the friction coefficient ( f ) on Re ∗ (blue),in comparison with the result obtained by Musacchioand Boffetta [11] (red). The dashed black line in theinset panel shows the curve of f = 0 .
124 + 5 . /Re ∗ ,which fits the red dots in the large Re ∗ range Musacchioand Boffetta [11]. (b) The global dissipation rate ( (cid:15) )and the coefficients obtained by fitting the profiles ofkinetic energy, α and β in Eq. 10, as function of Re λ .The horizontal dotted red and green lines represent thevalues α = 0 .
391 and β = 0 . Re ∗ = 2000.of an “angle” between vectors, see [18, 25]): ρ RS = R : S || R || || S || . (17)The ratio ρ RS is thus a number between − z P H D Q Y H O R F L W \ uvw sin(2 z ) z V K H D U V W U H V V H V u v u w v w Td udz
Figure 2: (a) The adimensional mean quantities of eachcomponent of the velocity of Run 8. The only non-zeroterm is (cid:104) u (cid:105) , having a maximum value of κ , where κ = 1 .
12. The black dotted line shows the curve of κ sin(2 πz ). (b) The different adimensional shear stressterms of Run 8. The only non-zero term is (cid:104) u (cid:48) w (cid:48) (cid:105) , whose z dependence is given by relation (8). The black dottedline shows the function − ν T d (cid:104) u (cid:105) dz , where ν T = 0 .
023 isthe turbulent viscosity (Eq. (13)) for the Run 8.is 1 when this hypothesis is valid, and when close to 0it corresponds to the case of “orthogonal” tensors. Thebehaviour of this quantity is shown in figure 5. It isseen that a plateau close to the value one appears incertain regions; in particular the Boussinesq’s hypothesisis approximately valid when the mean velocity gradientis large, whereas it fails dramatically for some range ofvalues around the positions where the mean velocity gra- z Q R U P D O V W U H V V H V u v w (a) Re λ = 38 . z Q R U P D O V W U H V V H V u v w (b) Re λ = 49 . z Q R U P D O V W U H V V H V u v w (c) Re λ = 66 . z Q R U P D O V W U H V V H V u v w (d) Re λ = 122 . Figure 3: The different normal stresses with (cid:104) u (cid:48) (cid:105) > (cid:104) w (cid:48) (cid:105) > (cid:104) v (cid:48) (cid:105) . The z -dependence is given by thefits of equation (9). z N L Q H W L F H Q H U J \ ( K ) Re = Re = Re = Re = Figure 4: The mean kinetic energy profile K ( z ) = (cid:104) u i u i (cid:105) for different Reynolds numbers.dient vanishes.More quantitatively, from Run 7, we find ρ RS = 0 .
93 for z = 1 / ρ RS = 0 . . ≤ ρ RS ≤ z ∈ [0 , . ∪ [0 . , . ∪ [0 . , ρ RS larger than 0.9, while for the restof the flow such linear relation fails to a large extent. z RS RS u (a) Re λ = 38 . z RS RS u (b) Re λ = 49 . z RS RS u (c) Re λ = 66 . z RS RS u (d) Re λ = 122 . Figure 5: Simulation results for the test of the validityof Boussinesq’s hypothesis, representing the alignment ρ RS between R and S . The mean velocity profile issuperposed in dotted line for reference. B. Test of a quadratic constitutive equation
We have seen above that the linear closure model can-not produce an anisotropic Reynolds stress tensor foranisotropic flows such as the Kolmogorov flow. Pope [26]has proposed to use the invariant theory in turbulencemodeling, to represent the stress tensor as a develop-ment into a tensor basis. Originally it was on the form R = (cid:80) i =1 a i T i with 10 basis tensors. By consideringa quadratic development, only three tensors are used,which is complete for two dimensional flows [26], and isalso a good approximation for fully 3-dimensional flows[27]. As it was used for channel flows [19, 28], we proposehere to use it also for the TKF.In this framework, the anisotropic stress tensor writes asa three-terms development R = a T + a T + a T ,where the three tensors of the basis are all symmetricand traceless [26]: T = S , T = SW − WS , T = S − η I . (18)The coefficients a i can be written using scalar invariantsof the flow, which correspond to scalar fields whose val-ues are independent of the system of reference. Invari-ants can be defined as the traces of different tensor prod-ucts [29]. Some of the first invariants are the following: η = { S } , η = { W } , η = { S } , η = { SW } , η = { S W } , µ = { R } , µ = { RS } , µ = { RSW } , and µ = { RS } . All these invariants can be here estimatednumerically. The coefficients a , a and a can be ex-pressed using the above invariants by projecting the con-stitutive equation onto the tensor basis: successive in-ner products of this equation with tensors T i provides asystem of scalar equations involving the invariants [27].For two-dimensional mean flows such as the KF, we have η = 0 and η = η η /
2, and the system of scalar equa-tions is inverted to provide finally the quadratic consti-tutive equation using invariants: R = µ η S − µ η η T + 6 µ η T (19)where the invariants write for the TKF: η = a η = − a µ = aτ (20) µ = a (cid:0) σ u − σ w (cid:1) , (21) µ = a (cid:18) σ v − K (cid:19) . (22)The two remaining tensors of the tensor basis are: T = a − (23)and T = a − . (24)The quadratic constitutive equation finally writes, re-placing invariants in equation (19): R = 2 τa S + (cid:0) σ u − σ w (cid:1) a T + (cid:0) σ v − K (cid:1) a T . (25)Equation (25) is a quadratic constitutive equation whichexpresses a nonlinear closure of the turbulent Kol-mogorov flow; the first constant coefficient is twice theeddy-viscosity ( τ /a = ν T = const. ), whereas the othercoefficients are space-dependent because they involve ra-tio of const. + cos divided by cos terms.When a = U (cid:48) ( z ) (cid:39)
0, for z (cid:39) / z (cid:39) / πz ) = 0 and all S , T and T vanish, but in thethree-terms development of R , the second term and thethird are non-zero constants, since the coefficients diverge(the a terms cancel). In those positions, we see that R is a diagonal tensor which is not vanishing: figure 6)shows that the second term is also very small and thatthe third term is dominant. This means that in thosepositions, the Boussinesq’s linear eddy-viscosity approxi-mation is no longer appropriate and the anisotropic stress z D P S O L W X G H V R I W K H W H U P V a || S || u w a || T || v Ka || T || u (a) Re λ = 38 . z D P S O L W X G H V R I W K H W H U P V a || S || u w a || T || v Ka || T || u (b) Re λ = 49 . z D P S O L W X G H V R I W K H W H U P V a || S || u w a || T || v Ka || T || u (c) Re λ = 66 . z D P S O L W X G H V R I W K H W H U P V a || S || u w a || T || v Ka || T || u (d) Re λ = 122 . Figure 6: The amplitudes of the terms at the right handside of Eq. (25) as function of z . The mean velocityprofile is also represented as a dotted line, for reference.The horizontal red dotted lines mark the 0 value for theamplitudes. The simulation results of Run 1, 3, 5, 7 areshown here.tensor is a constant perpendicular to the linear term andapproximately proportional to T = S − η I . IV. PERIODIC FLOW WITHNON-SINUSOIDAL FORCING
As we have seen, the Kolmogorov flow, in its originaldefinition, is sustained by a monochromatic sinusoidalforcing. The resulting mean flow profile is also sinusoidalwith the same shape of the force term, both in laminarand in turbulent flow conditions. This peculiar propertyis of great help in the analysis of the turbulent flow and,as we have seen, it simplifies the formulation of a closurerelation. It is therefore of interest to ask what happenswhen the shape of the force is changed to other periodicor quasi-periodic shapes. Here we investigate a forcinghaving a Gaussian shape. Although this function is non-periodic and has an unbounded support, one can adjustits width in such a way that its value and its derivativebecomes sufficiently small at the borders. Furthermore,the Gaussian has the advantage to be easily implementedin spectral space. Here we test only one given Reynoldsnumber, as this is sufficient to contrast the qualitativedifferences with the sinusoidal forcing case. The dimen-sionless viscosity for this case are 0 . . f = ( A exp (cid:18) − ( z − . (cid:96) (cid:19) + C ) e x , (26)where A is the forcing parameter, C is a constant to make f = 0 in equation (1) and (cid:96) = 0 . H units), whichis equivalent to the standard deviation, and controls thewidth of the Gaussian shape.Figure 7 shows the Reynolds stress components: theshear stress and normal stresses. The normal stresseshave a shape that can not be fitted with known func-tions and only the numerical result is shown here. It isvisible that normal stresses are again anisotropic, with σ u > σ w > σ v at all positions. The shear stress τ = −(cid:104) u (cid:48) w (cid:48) (cid:105) is as in the TKF the only non-zero shearstress, and is proportional to U (cid:48) ( z ), with a coefficient ν T = 0 .
016 (it was found of the same order, 0 . τ = ν T U (cid:48) ( z ) weobtain U (cid:48)(cid:48) ( z ) = − Reν T /ν + 1 f . (27)This shows that the mean velocity profile is proportionalto twice the integral of the forcing. The primitive of theGaussian is non-analytical and involves the error func-tion; it can be estimated numerically, as shown in figure8. An excellent superposition is found. The shape of themean velocity is still Gaussian-like, but its precise ana-lytical expression is given by the function whose secondderivative is a Gaussian.The alignment between the tensors of R and S (as de-fined in Eq. (17)) is also examined for the case withGaussian forcing, as shown in Fig. 9. We observe also inthis case the plateaus obtained at the position where themean velocity gradient is large. Qualitatively, we findthat 0 . ≤ ρ RS ≤ z ∈ [0 . , . ∪ [0 . , . V. CONCLUSION
We have considered here the closure relation for theReynolds stress in a numerically simulated turbulent Kol-mogorov flow. As the simplest realization of turbulencewith a spatially dependent mean flow, such model sys-tem is a convenient test ground for turbulent transportmodels. With a forcing of the form sin(2 πz ), it was found z Q R U P D O V W U H V V H V u v w (a) z V K H D U V W U H V V H V D Q G P H D Q G H U L Y D W L Y H u v u w v w Td udz (b)
Figure 7: (a) The normal stresses; (b) The shearstresses and the function of − ν T d (cid:104) u (cid:105) dz (black dotted line),where ν T is the turbulent viscosity (Eq. (13)) andnumerically found as 0 .
016 here, for the case withGaussian forcing (same parameters with Run 4). z P H D Q Y H O R F L W \ uvw Re T / + 1 f x dz Figure 8: Comparison of the mean profile ( (cid:104) u (cid:105) , bluesolid line) and the double integral of the forcing (blackdotted line) for the case with Gaussian forcing (sameparameters with Run 4). The coefficient of Reν T /ν +1 isnumerically found as 61.27. The black dashed linerepresents the Gaussian profile forcing (26).that the mean velocity profile has the same form, with adamping of a factor κ with respect to the mean velocityvalue calculated from forcing terms. The value of κ wasfound to increase with the Reynolds number, and of theorder of 1.01 to 1.12 for the range of Reynolds numbersconsidered here. The only non-zero shear stress termis proportional to cos(2 πz ) as expected, and the normalstress components all involve a square cosine expressionof the form α + β cos (2 πz/H ), where the parameters α and β are numerically estimated and found to saturatefor the largest Reynolds numbers considered here. The z RS RS u Figure 9: Estimation of the alignment ρ RS between R and S , for the case with Gaussian forcing (sameparameters with Run 4). The mean velocity profile issuperposed in dotted line for reference. z D P S O L W X G H V R I W K H W H U P V a || S || u w a || T || v Ka || T || u Figure 10: The amplitudes of the terms at the righthand side of Eq. (25) as function of z for the case withGaussian forcing (same parameters with Run 4). Themean velocity profile is also represented as a dottedline, for reference. The horizontal red dotted lines markthe 0 value for the amplitudes.normal stresses, i.e. (cid:104) u (cid:48) (cid:105) , (cid:104) w (cid:48) (cid:105) , (cid:104) v (cid:48) (cid:105) , are all different inamplitude, showing that the turbulence is anisotropic.It was also shown that a quadratic nonlinear constitutiveequation can be proposed for this flow. Specifically a lin-ear term and two nonlinear terms in the form of tracelessand symmetric tensors SW − WS and S − { S } I areinvolved and their coefficients are here numerically esti-mated. For about half of the flow domain, the linear term is dominating, whereas for the vanishing mean velocityregions a constant term is the only one remaining. Hencean effective viscosity coefficient can indeed be estimatedfor the Kolmogorov flow, but contrary to what has beenindicated previously [12] this type of turbulence withoutboundaries does not generate an effective diffusion of mo-mentum, since nonlinear terms are needed: globally alllinear and nonlinear terms are needed for the completeReynolds stress closure.The values obtained here are in agreement with previ-ous works [6, 11]. Using 8 different runs with differentgrid sizes from 128 to 512 , and with Reynolds numbersfrom Re λ = 39 to 198, the Reynolds number dependenceof involved parameters has been checked with expectedconvergence toward the largest Reynolds number consid-ered.Finally a periodic flow with non-sinusoidal forcing hasbeen considered, with the choice of a Gaussian shape. Itwas found that the shear stress term τ is proportionalto the mean velocity derivative, indicating that for suchforcing also the eddy-viscosity does not depend on z .With such numerical result, we obtain that the meanvelocity profile is twice the integral of the forcing. Theshape of the normal stresses in this case is non-trivial andcan not be precisely fitted. A quadratic development ofthe constitutive equation can also be applied to this flow.We may qualitatively compare here the turbulent flowssustained by the sinusoidal forcing (TKF) and the Gaus-sian forcing. In both cases the eddy-viscosity is found tobe constant. We observe that the relation stating that U (cid:48)(cid:48) is proportional to the forcing (27), is in fact also valid forthe TKF case. There are however, two main differencesbetween the two forcing cases. The first lies in the shapeof the mean velocity profile. For the TKF, since the sec-ond integral of the forcing is proportional to the forcing,equation (27) directly gives the mean velocity profile, asbeing proportional to the forcing. In the Gaussian forc-ing case, the mean velocity is a non-analytical function,obtained as the second integral of the Gaussian. The sec-ond difference is in the shape of normal stresses, whoseexpression could be fitted using cos terms for the TKF,while no known analytical fit for the Gaussian forcing isavailable.As a perspective let us mention the recent work [30] de-scribing a flow behind a grid in a wind tunnel as hav-ing locally, close to the grid, sinusoidal variations. Thiscould provide ideas to perform measurements of a TKFand to check experimentally the closure proposed in thepresent study. A more systematic numerical study ofnon-sinusoidal forcing in future work may also help toprovide a general expression for normal stresses, whichwould be valid for all kind of forcing. It remains also tobe understood from analytical arguments why the eddy-viscosity does not depend on z for such periodic flow,contrary to what is found in similar bounded flow suchas the channel flow [19] or the boundary-layer flow.0 Acknowledgements
This work is under the joint supportof Shanghai Jiao Tong University and the French Region“Hauts-de-France” in the framework of a cotutella PhDprogramme. We thank Dr. Michael Gauding (CORIA(CNRS UMR 6614), Rouen, France) for providing a nu- merical code which was adapted for our specific presenttopic. We acknowledge the computing resources includ-ing the High Performance Computing Center (HPCC) atUniversit´e de Lille, CALCULCO of Universit´e du LittoralCˆote d’Opale and the National Supercomputer Center inGuangzhou, China. [1] L. D. Meshalkin and I. G. Sinai, Investigation of the sta-bility of a stationary solution of a system of equations forthe plane movement of an incompressible viscous liquid,Journal of Applied Mathematics and Mechanics , 1700(1961).[2] J. S. A. Green, Two-dimensional turbulence near the vis-cous limit, Journal of Fluid Mechanics , 273 (1974).[3] N. F. Bondarenko, M. Z. Gak, and F. V. Dolzhan-skii, Laboratory and theoretical models of plane periodicflow, Akademiia Nauk SSSR, Izvestiia, Fizika Atmosferyi Okeana , 1017 (1979).[4] A. M. Obukhov, Kolmogorov flow and laboratory simula-tion of it, Russian Mathematical Surveys , 113 (1983).[5] V. I. Kliatskin, On the nonlinear theory of stability ofperiodic flows, Journal of Applied Mathematics and Me-chanics , 243 (1972).[6] V. Borue and S. A. Orszag, Numerical study of three-dimensional Kolmogorov flow at high Reynolds numbers,Journal of Fluid Mechanics , 293 (1996).[7] J. V. Shebalin and S. L. Woodruff, Kolmogorov flow inthree dimensions, Physics of Fluids , 164 (1997).[8] L. Biferale and F. Toschi, Anisotropic homogeneous tur-bulence: hierarchy and intermittency of scaling expo-nents in the anisotropic sectors, Physical Review Letters , 4831 (2001).[9] N. J. Balmforth and Y.-N. Young, Stratified Kolmogorovflow, Journal of Fluid Mechanics , 131 (2002).[10] G. Boffetta, A. Celani, A. Mazzino, A. Puliafito, andM. Vergassola, The viscoelastic Kolmogorov flow: eddyviscosity and linear stability, Journal of Fluid Mechanics , 161 (2005).[11] S. Musacchio and G. Boffetta, Turbulent channel with-out boundaries: The periodic Kolmogorov flow, PhysicalReview E , 023004 (2014).[12] B. Rollin, Y. Dubief, and C. R. Doering, Variations onKolmogorov flow: turbulent energy dissipation and meanflow profiles, Journal of Fluid Mechanics , 204 (2011).[13] Z. S. She, Metastability and vortex pairing in the Kol-mogorov flow, Physics Letters A , 161 (1987).[14] S. Berti and G. Boffetta, Elastic waves and transition toelastic turbulence in a two-dimensional viscoelastic Kol-mogorov flow, Physical Review E , 036314 (2010).[15] D. Lucas and R. R. Kerswell, Spatiotemporal dynamicsin two-dimensional Kolmogorov flow over large domains,Journal of Fluid Mechanics , 518 (2014). [16] D. Lucas and R. R. Kerswell, Recurrent flow analysisin spatiotemporally chaotic 2-dimensional Kolmogorovflow, Physics of Fluids , 045106 (2015).[17] S. B. Pope, Turbulent flows (Cambridge UniversityPress, 2000).[18] F. G. Schmitt, About boussinesq’s turbulent viscosity hy-pothesis: historical remarks and a direct evaluation of itsvalidity, Comptes Rendus M´ecanique , 617 (2007).[19] F. G. Schmitt, Direct test of a nonlinear constitutiveequation for simple turbulent shear flows using dns data,Communications in Nonlinear Science and NumericalSimulation , 1251 (2007).[20] T. Y. Hou and R. Li, Computing nearly singular solu-tions using pseudo-spectral methods, Journal of Compu-tational Physics , 379 (2007).[21] I. E. Sarris, H. Jeanmart, D. Carati, and G. Winckel-mans, Box-size dependence and breaking of translationalinvariance in the velocity statistics computed from three-dimensional turbulent Kolmogorov flows, Physics of Flu-ids , 095101 (2007).[22] J. Boussinesq, Essai sur la th´eorie des eaux courantes(M´emoires pr´esent´es par divers savants `a l’Acad´emie desSciences vol. XXIII, Imprimerie Nationale, 1877).[23] C. G. Speziale, On nonlinear kl and k- ε models of turbu-lence, Journal of Fluid Mechanics , 459 (1987).[24] S. Nisizima and A. Yoshizawa, Turbulent channel andcouette flows using an anisotropic k-epsilon model, AIAAjournal , 414 (1987).[25] F. G. Schmitt and C. Hirsch, Experimental study of theconstitutive equation for an axisymmetric complex tur-bulent flow, ZAMM/Zeitschrift f¨ur Angewandte Mathe-matik und Mechanik , 815 (2000).[26] S. B. Pope, A more general effective-viscosity hypothesis,Journal of Fluid Mechanics , 331 (1975).[27] T. Jongen and T. B. Gatski, General explicit alge-braic stress relations and best approximation for three-dimensional flows, International Journal of EngineeringScience , 739 (1998).[28] D. Modesti, A priori tests of eddy viscosity models insquare duct flow, Theoretical and Computational FluidDynamics , 713 (2020).[29] A. J. M. Spencer, Part iii. theory of invariants, Contin-uum Physics , 239 (1971).[30] W. T. Bos, Production and dissipation of kinetic en-ergy in grid turbulence, Physical Review Fluids5