Calibrating passive scalar transport in shear-flow turbulence
aa r X i v : . [ a s t r o - ph . S R ] J u l Calibrating passive scalar transport in shear-flow turbulence
Enik˝o J. M. Madarassy and Axel Brandenburg NORDITA, AlbaNova University Center, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, SE 10691 Stockholm, Sweden (Dated: November 19, 2018, Revision: 1.32 )The turbulent diffusivity tensor is determined for linear shear flow turbulence using numerical simulations.For moderately strong shear, the diagonal components are found to increase quadratically with Peclet andReynolds numbers below about 10 and then become constant. The diffusivity tensor is found to have compo-nents proportional to the symmetric and antisymmetric parts of the velocity gradient matrix, as well as productsof these. All components decrease with the wave number of the mean field in a Lorentzian fashion. The compo-nents of the diffusivity tensor are found not to depend significantly on the presence of helicity in the turbulence.The signs of the leading terms in the expression for the diffusion tensor are found to be in good agreement withestimates based on a simple closure assumption.
PACS numbers: PACS Numbers : 47.27.tb, 47.27.ek, 95.30.Lz
I. INTRODUCTION
In a turbulent flow, chemicals tend to be mixed more ef-fectively than in the absence of turbulence. Indeed, turbu-lence disperses chemicals by advecting particles along chaotictrajectories. This rapidly causes large concentration gradi-ents that speed up their mixing down toward the smallestscales. Turbulent mixing is a complicated and rich process;see Ref. [1] for a comprehensive review on this subject. Themathematical treatment of the description of turbulent mixingis closely related to that of turbulence itself, but it is in manyways much simpler and provides therefore an ideal tool formaking conceptual progress in that field [2].Here we are mainly interested in cases where it is mean-ingful to define a mean concentration whose scale of variationis large compared with the scale of the energy-carrying ed-dies. In such cases it can be useful to describe the change inthe mean concentration by an effective turbulent diffusion ten-sor. On smaller scales the change in the mean concentrationcan still be described in such a way, but in that case the mul-tiplication with a turbulent diffusivity must be replaced witha convolution. The turbulent diffusion tensor quantifies theeffective exchange of chemicals or other passive scalar quan-tities advected by the flow. If there is a gradient in the meanconcentration C of chemicals, there will be a net mean flux F = u c of chemicals resulting from a systematic correlationof fluctuations in the concentration c and the turbulent velocity u . Here, overbars denote averaging. Under isotropic condi-tions with sufficient scale separation, this mean flux will bedown the gradient of concentration, with F = − κ t ∇ C, (1)where κ t is the turbulent diffusivity. However, modificationsare expected when the turbulence is anisotropic. In that casethis relation takes the form F i = − κ ij ∇ j C, (2)where κ ij is now the turbulent diffusion tensor. In this paperwe are interested in the anisotropy caused by the presence ofshear. One of the results one expects to see is a suppression of turbulent transport in the cross-stream direction. This effect isdiscussed in various physical circumstances such as geophys-ical flows [3], turbulent plasmas [4], and solar physics [5, 6].Much of this research is done using analytical techniquessuch as the first-order smoothing approximation and the renor-malization group analysis. However, in recent years it hasbecome possible to calculate turbulent transport coefficientsusing numerical realizations of turbulence from direct simula-tions. Turbulent transport coefficients can then be determinedby imposing a gradient in the passive scalar concentration andmeasuring the resulting concentration fluxes [7]. By imposinggradients in three different directions it is possible to assembleall components of the turbulent diffusion tensor.In recent years such a technique has been applied to the caseof magnetic fields whose evolution is controlled not just byturbulent magnetic diffusion, but also by non-diffusive contri-butions known as the α effect [8, 9]. In this way it has beenpossible to investigate numerically the effects of shear and ro-tation in regimes that cannot be treated analytically. The tech-nique is known under the name test-field method, which refersto the fact that this approach involves the analysis of correla-tions for a set of different pre-determined test fields. In theanalogous case of passive scalars, this method is now oftenreferred to as test-scalar method [10].Using this method, it has recently been possible to deter-mine the turbulent diffusion tensor in cases where the turbu-lence is anisotropic owing to the presence of either rotationor an imposed magnetic field [10]. In the case of rotation theangular velocity vector Ω provides a new element for con-structing an anisotropic rank-2 tensor of the form [11] κ ij = κ δ ij + κ Ω ǫ ijk ˆΩ k + κ ΩΩ ˆΩ i ˆΩ j , (3)where ˆ Ω = Ω / | Ω | is the unit vector along the rotation axisand κ , κ Ω , and κ ΩΩ are functions of the flow parameters.Note that Ω is a pseudo vector while κ ij is a proper tensor,so all three coefficients in Eq. (3) are proper scalars. In thecase of a shear flow, an obvious possible ansatz is obtainedby replacing Ω with the vorticity W = ∇ × U , which isalso a pseudovector (or axial vector), and U is the mean shearflow. However, such an ansatz would be incomplete, becauseit only captures the antisymmetric part of the velocity gradientmatrix U i,j , where a comma denotes partial differentiation.A more natural approach would therefore be to invoke bothsymmetric and antisymmetric parts of the velocity gradientmatrix by writing it as U i,j = S ij + A ij , where S ij = ( U i,j + U j,i ) , (4) A ij = ( U i,j − U j,i ) . (5)The latter can also be written as A ij = − ǫ ijk W k . A properrank-2 tensor can then be expressed as κ ij = κ t δ ij + κ S S ij + κ A A ij + κ SS ( S S ) ij + κ AS ( A S ) ij , (6)where κ t , κ S , κ A , κ SS , and κ AS are proper scalars that areagain functions of the flow parameters. In the absence of he-licity, no further rank-2 tensors can be constructed from a lin-ear shear flow. We return to the case with helicity in Sec. III D.An important goal of this work is to determine the coeffi-cients in Eq. (6) for a linear shear flow of the form U = (0 , Sx, , (7)where S = const is the shear rate, which is not to be confusedwith the tensor S . For a linear shear flow given by Eq. (7), thetensors S and A are constants, and their only non-vanishingcomponents are S xy = S yx = − A xy = A yx = S/ . (8)Note also that S = − A = ( S/ diag (1 , , , (9) A S = − S A = ( S/ diag ( − , , . (10)With these preparations we can now express all nine compo-nents of κ ij in terms of the five coefficients in Eq. (6) as fol-lows: κ = κ t + S ( κ SS − κ AS ) , (11) κ = κ t + S ( κ SS + κ AS ) , (12) κ = κ t , (13) κ = S ( κ S − κ A ) , (14) κ = S ( κ S + κ A ) , (15) κ = κ = κ = κ = 0 . (16)Given that all nine components of κ ij can be determined fromsimulation data using the test-scalar method, we can use the relations above to compute the five unknown coefficients inEq. (6) via κ S S = κ + κ , κ A S = κ − κ , (17) κ SS S / κ + κ − κ , (18) κ AS S / κ − κ , (19) κ t = κ . (20)Note that combinations such as κ S S and κ SS S / have stillthe same dimension as κ ij , so in the following we shall quotethese combinations in that form.In principle it is possible to construct κ ij using also the ve-locity vector U itself. However, U varies in x and vanishes at x = 0 . On the other hand, we expect the components of κ ij not to depend explicitly on position, making a constructionin terms of U less favorable. Furthermore, the tensor U i U j ,which has only one component in the yy position, can alreadybe constructed from S − A S = diag (0 , , , so no new in-formation would be added. However, this changes when wealso admit helical turbulent flows, because then there couldbe tensors of the form W i U j and W j U i which have compo-nents in the yz and zy directions. For this reason we shall alsoinvestigate helical turbulence in some cases.A comment regarding the case of rotation without shear ishere in order. In hindsight it might have been more naturalto write Eq. (3) in terms of the antisymmetric matrix A ij = − ǫ ijk ˆΩ k , i.e. κ ij = κ Ωt δ ij + κ Ω A A ij + κ Ω AA ( A ) ij , (21)with coefficients that are related to those in Eq. (3) via κ Ωt = κ + κ ΩΩ , κ Ω A = − κ Ω , κ Ω AA = 4 κ ΩΩ . (22)Evidently, this representation is equivalent to that of Eq. (3).In the rest of this paper we continue with the case of a pureshear flow. The aim is to determine the coefficients in Eq. (6)as functions of flow parameters such as the Peclet number andthe shear parameter. II. SIMULATIONS
We simulate turbulence by solving the compressible hydro-dynamic equations with an imposed random forcing term andan isothermal equation of state, so that the pressure p is relatedto ρ via p = ρc , where c s is the isothermal sound speed.We consider a periodic Cartesian domain of size L . In thepresence of shear the hydrodynamic equations for ρ and thedeparture U from the imposed shear flow U take the form, D ln ρ D t = − ∇ · U , (23) D U D t = − SU x ˆ y − c ∇ ln ρ + f + F visc , (24)where D / D t = ∂/∂t + ( U + U ) · ∇ is the advective deriva-tive with respect to the full velocity, F visc = ρ − ∇ · ρν S is the viscous force, ν is the kinematic viscosity, S ij = ( U i,j + U j,i ) − δ ij ∇ · U is the traceless rate of strain ten-sor of the departure from the shear flow, and f is a randomforcing function consisting of plane transversal waves withrandom wave vectors k such that | k | lies in a band arounda given forcing wave number k f . The vector k changes ran-domly from one timestep to the next, so f is δ correlated intime. We have carried out simulations with helical and non-helical forcings using the modified forcing function f k = R · f (nohel) k with R ij = δ ij − i σǫ ijk ˆ k k √ σ , (25)where f (nohel) k is the non-helical forcing function. In the fullyhelical case ( σ = ± ) we recover the forcing function usedin Ref. [12], and in the non-helical case ( σ = 0 ) this forc-ing function becomes equivalent to that used in Ref. [13].The forcing amplitude is chosen such that the Mach number,Ma = u rms /c s , is about 0.1. We use triply-periodic boundaryconditions, except that the x direction is shearing–periodic,i.e. U ( − L, y, z, t ) = U ( L, y + LSt, z, t ) , (26)where L is the side length of the cubic domain. This condi-tion is routinely used in numerical studies of shear flows inCartesian geometry [14, 15].In this paper we are interested in the turbulent mixing of apassive scalar concentration C . Its evolution is governed bythe equation ∂C∂t = − ∇ · ( U C ) + κ ∇ C, (27)where κ is the microscopic (molecular) passive scalar diffu-sivity. In the absence of any sources, the dynamics of C depends essentially on initial conditions. For example, if C is initially concentrated in a plane with its normal pointingin one of the three coordinate directions, turbulence tends tospread this initial distribution away from the plane – regard-less of its orientation. Only the speed of spreading will bedifferent in the different directions. The spreading is then bestdescribed by introducing planar averages over the same direc-tions as the initial distribution. These averages are denotedby overbars and they depend only on time and the directionnormal to the plane of averaging, i.e. C = C ( x j , t ) , where x j denotes x , y , or z for j = 1 , ..., 3, just depending on the ini-tial distribution. This allows us then to quantity the speed ofspreading by the different components of the diffusion tensor κ ij in Eq. (2). We do this by introducing different ‘test scalars’and calculating the evolution for each case separately..In the following we are interested in the fluxes of the pas-sive scalar concentration, F = u c , where c = C − C is thefluctuation around the mean concentration and u = U − U is the velocity fluctuation around the mean flow U . The test-scalar equation is obtained by subtracting the averaged pas-sive scalar equation from the original one and applying it to apredetermined set of six different mean fields, C ic = C cos kx i , C is = C sin kx i , (28)where C is a normalization factor. Again, the overbars de-note planar averaging over the directions that are perpendic-ular to the direction in which the mean field varies. For eachtest field C pq we obtain a separate evolution equation for thecorresponding fluctuating component c pq , ∂c pq ∂t = − ∇ · ( U c pq + u C pq + u c pq − u c pq )+ κ ∇ c pq , (29)where p = 1 , ..., 3, and q = c or s . In this way, we calcu-late six different fluxes, F pq = u c pq , and compute the ninerelevant components of κ ij , κ ij = −h cos kx j F jsi − sin kx j F jci i /k, (30)for i, j = 1 , ..., . Here, angular brackets denote volume aver-ages. A visualization of c s , c s , and c s on the periphery ofthe computational domain is shown in Fig. 1 after about oneturnover time for a run with k/k f = 0 . , which is smaller thanin most of the runs analyzed in this paper. This ratio is cho-sen here for visualization purposes only, because this way thelarge-scale modulation compared with the scale of the turbu-lence becomes evident.We emphasize that Eq. (29) is an inhomogeneous equationin c pq . The term u C pq can be regarded as a forcing termthat guarantees that the direction of the turbulent concentra-tion flux will not change with time.In this paper we present the values of κ ij in non-dimensional form by normalizing with κ t0 = u rms / k f , (31)which is the expected value for large values of Pe. Here wehave defined the root-mean-square value of the velocity fluc-tuation as u rms = h u i / .Our simulations are characterized by two important non-dimensional control parameters, the shear parameter Sh andthe Peclet number Pe, defined asSh = S/ ( u rms k f ) , Pe = u rms / ( κk f ) . (32)In addition, there is the Schmidt number Sc = ν/κ , but wekeep it equal to unity in all cases reported below. Note alsothat in most cases we use negative values of S , so we haveSh < . The smallest wave number that fits into the compu-tational domain is k = 2 π/L . In most of the cases reportedbelow we choose the forcing wave number to be 3 times larger,i.e. k f /k = 3 .The simulations have been carried out using the P ENCIL C ODE [34] which is a high-order finite-difference code (sixthorder in space and third order in time) for solving the com-pressible hydrodynamic equations. The test-scalar equationswhere already implemented into the public-domain code, but
FIG. 1: (Color online) Visualization of c s , c s , and c s on the periphery of the computational domain after about one turnover time for arun with k/k f = 0 . . In the middle panel, arrows indicate the direction of the shear flow with negative S , i.e. d U y / d x < . Note the clearsinusoidal modulation in the x , y , and z directions for the three panels, respectively. In the middle panel this modulation is already smearedout by the shear. have now been generalized to determining all nine compo-nents of κ ij . The numerical resolution used in the simulationsdepends on the Peclet number and reaches meshpointsfor runs with Pe ≈ . In this paper we restrict ourselves totime spans short enough so that the so-called vorticity dynamohas no time do develop; see Refs. [16, 17] for details on thiseffect. III. RESULTSA. Dependence on the shear parameter
We begin by discussing the dependence of the coefficientsin Eq. (6) on the shear parameter Sh. The result is shown inFig. 2 for Pe = 25 . It turns out that all five coefficients arepositive. We find that κ t /κ t0 = const = 2 for non-helicalturbulence and 3 for helical turbulence, independent of thevalue of shear, provided | Sh | < . . The other coefficientsshow the following approximate scaling behavior: κ S S/κ t0 ≈ | Sh | , κ SS S / κ t0 ≈ Sh , (33) κ A S/κ t0 ≈ | Sh | , κ AS S / κ t0 ≈ | Sh | . (34)The fact that κ A and κ AS scale with the third power of Shsuggests that these are higher order effects that are not easilycaptured by perturbative approaches.A comment regarding the values of Sh is here in order. Al-though values of Sh larger than unity have not yet been ex-plored, it is unlikely that the uprise of κ t continues. Further-more, one might speculate that all coefficients in Eqs. (33)and (34) should eventually decrease as | Sh | → ∞ .In Fig. 2 we have also shown results for cases where theforcing function has maximum helicity. No significant depen-dence can be seen, except for κ t which is slightly enhancedin the helical case with weak shear. This suggests that thisdependence is not connected with the presence of shear. B. Dependence on Peclet number
We have performed simulations for different values of thePeclet number and have determined the coefficients in Eq. (6)for each simulation. The results are shown in Fig. 3 for fixedSh = 0 . . It turns out that the first four coefficients can wellbe approximated by simple algebraic functions, κ t κ t0 = 2 κ Sh Pe Pe + Pe , κ S Sκ t0 = κ Sh Pe ( Pe + Pe ) / , (35) κ SS S κ t0 = κ Sh Pe ( Pe + Pe ) , κ A Sκ t0 = κ Sh Pe ( Pe + Pe ) . , (36)where κ Sh = 0 . κ t0 and Pe = 3 . are fit parameters. Inthe case of κ AS the error bars are so large that no conclusivestatements can be made. Likewise, the error bar on the firstdata point is quite large too. This is caused by the numericaltime step becoming rather short at large diffusivities, so therun is short and the statistics poor.In Fig. 4 we show the dependence of the diagonal com-ponents of κ ij on Pe. Over the range of parameters shownhere, the difference between the three components is small,although there is a tendency for κ yy to be somewhat enhancedaround Pe = 20 , while κ zz is slightly smaller than κ xx . C. Wavenumber dependence
We consider now the dependence of the diagonal com-ponents of κ ij on the wave number k of the test scalar inEq. (28). A dependence of κ ij on k reflects the fact that thereis poor scale separation, i.e. k/k f is no longer small. In such acase, the multiplication with a turbulent diffusivity in Eqs. (1)and (2) must be replaced by a convolution with an integralkernel [10]. In Fourier space the convolution corresponds toa multiplication. The full integral kernel can be assembledby determining the full k dependence and then Fourier trans-forming back into real space. FIG. 2: Dependence of the coefficients in Eq. (6) on Sh for Pe = 25 .The dashed line in the first panel is for a run with maximum helicity.All runs with helicity are marked with open symbols. Filled symbolsindicate runs without helicity. Solid lines represent the fits given byEqs. (33) and (34).
The resulting dependence on k is shown in Fig. 5 for twovalues of the shear parameter and Pe around 50. In agreementwith earlier findings, the components of κ ij show a Lorentziandependence on k , i.e. κ ij = κ (0) ij ak/k f ) , (37)where a ≈ . for the κ and κ components, and a ≈ . for the κ component. Here, κ (0) ij is the value for k = 0 ,which is approximately equal to κ t0 , defined in Eq. (31).Given that the Schmidt number is always kept equal tounity, there will be a fully developed cascade in the passivescalar concentration when the Peclet number is large. The va-lidity of Eq. (37) has only been tested for values of Pe up to60. It is unclear whether this equation holds also for large FIG. 3: Dependence of the coefficients in Eq. (6) on Pe for Sh = − . . The symbols give the numerical results and the solid linesrepresent fits given by Eqs. (35) and (36).FIG. 4: Dependence of the diagonal components of κ ij on Pe. values of Pe when contributions from the high wave numberdynamics may become important in the mixing of the mean FIG. 5: Dependence of the diagonal components of κ ij on k forSh = − . at Pe = 40 (upper panel) and Sh = − . at Pe = 60 (lower panel). concentration.The case of high wave numbers is interesting in view ofpossible applications of our results to subgrid scale modelingin large-eddy simulations of turbulence. The highest possi-ble wave number is the Nyquist wave number, k Ny = π/δx ,where δx is the mesh scale. In the Smagorinsky model [18]the subgrid scale viscosity is proportional to the modulus ofthe rate of strain tensor times δx . For a turbulent flow wherethe local velocity difference δu ℓ over a distance ℓ is propor-tional to ℓ / we expect the subgrid scale viscosity to be ef-fectively proportional to ℓ / , suggesting an asymptotic k − / scaling for k ≫ k f . Here we have identified ℓ with δx andthus k with k Ny . Only for a smooth velocity field, where δu ℓ scales linearly with the separation ℓ , the subgrid scale viscos-ity would be proportional to ℓ , justifying an asymptotic k − scaling. This uncertainty warrants further studies of the valid-ity of Eq. (37) for k ≫ k f . D. Effects of helicity
As discussed in the Introduction, the presence of helicity al-lows one in principle to construct proper tensors proportionalto W i U j and W j U i , because we have now access to a pseu-doscalar given by the kinetic helicity of the turbulence. If thisdoes indeed have an effect, one would expect finite yz and zy components. In Fig. 6 we present results for κ yz and κ zy us-ing Pe = 25 . We see that κ yz = κ zy = 0 within error bars, so FIG. 6: Plot of κ yz (dotted line) and κ zy (dashed line) versus Sh formaximally helical turbulence and Pe = 25 . No significant depen-dence can be seen. there is no evidence for the presence of additional terms whenthe turbulence is helical. IV. EXPECTATIONS FROM THE τ APPROXIMATION
Passive scalar transport is closely related to the transportof a mean magnetic field. Commonly applied techniquesfor computing turbulent transport coefficients in mean-fieldelectrodynamics are the first order smoothing approximation[19, 20] and the τ approximation [21–23]. The τ approxi-mation consists in writing down an evolution equation for thequadratic correlations which, in the case of mean-field elec-trodynamics, is the mean electromotive force E . Its solutiongives then an expression for E in terms of the mean magneticfield and its derivatives. For a recent review see Ref. [24].This technique has also been used to compute the Reynoldsand Maxwell stress in rotating shear flows [25–27]. In thepresent case of passive scalar transport one starts with theevolution equation for the mean flux F = u c , as is done inRefs. [7, 28]. Thus, we write ∂ F i ∂t = ˙ u i c + u i ˙ c, (38)where dots denote time derivatives that are given essentiallyby Eqs. (24) and (29). This results in quadratic and triple cor-relations. The sum of all triple correlations is substituted bya damping term of the form − F /τ on the right-hand side ofthe evolution equation for F . Here, τ = St /u rms k f is theturnover time and St is a positive dimensionless parameter oforder unity (referred to as Strouhal number). This is a closureassumption that cannot be motivated rigorously [29], but ithas been found numerically that the triple-correlations are in-deed locally and temporally proportional to the negative fluxterm divided by τ ; see Ref. [7] for passive scalar diffusion andRef. [30] for the case of mean-field electrodynamics.As a first orientation, and in order to gain some understand-ing of our numerical results, we make the additional assump-tion that we can subsume the effects of the pressure term inour closure assumption. Since our forcing function f is δ cor-related in time we have f c = and thus obtain ˙ u i c = − Sδ i δ k u k c + triple correlations , (39) u i ˙ c = − u i u j ∇ j C + triple correlations . (40)The triple correlation terms result from the nonlinearities inthe evolution equations, Eqs. (24) and (29). In the τ ap-proximation one substitutes the sum of the triple correlationsby quadratic correlations, i.e. in the present case by − u i c/τ [21, 31]. We write the resulting equation in matrix form, τ ∂ F i ∂t = − L ik F k − τ u i u j ∇ j C, (41)where L ik = δ ik + Sτ δ i δ k . We solve this equation for F and obtain F i = − ( L − ) ij (cid:18) τ u j u k ∇ k C + τ ∂ F i ∂t (cid:19) , (42)where ( L − ) ik = δ ik − Sh δ i δ k with Sh = Sτ . In the pres-ence of shear, the Reynolds stress tensor u j u k is no longerdiagonal, but it has finite xy and yx components. Also thethree diagonal components are no longer the same. In the fol-lowing we represent u j u k in the form uu = u x − δ − δ ǫ
00 0 1 + ǫ z ! , (43)where δ = − u x u y /u x characterizes the value of the off-diagonal components, while ǫ = u y /u x − and ǫ z = u z /u x − characterize the change in the two lower diago-nal components. The dependence of δ and ǫ on Sh is shown inFig. 7, while ǫ z is found to be small. Inserting this expressioninto Eq. (42), we obtain κ κ t0 = − δ − δ − Sh ǫ + δ Sh
00 0 1 + ǫ z ! . (44)In the stationary state we may ignore the time derivative andrecover Eq. (3) with κ S Sκ t0 = − δ − Sh , κ A Sκ t0 = − Sh , (45) κ SS S / κ t0 + 2 ǫ z = κ AS S / κ t0 = ǫ − δ Sh . (46)We recall that Sh is negative, and that δ changes sign with Sh.Therefore we expect κ S and κ A to be positive, which agreeswith the simulations. Furthermore, we expect κ yy to be en-hanced, which also agrees with the simulations. However, theslight suppression of κ zz cannot be explained by the simpletheory, because ǫ z is small and perhaps even positive, sug-gesting at best an opposite trend. FIG. 7: Dependence of δ and ǫ on Sh for nonhelical turbulence (solidsymbols) and helical turbulence (open symbols). The solid line in thesecond panel has a slope of 5. V. CONCLUSIONS
The present work has shown that shear introducesanisotropies in the diffusivity tensor for passive scalar dif-fusion. These additional components are proportional to theeven and odd parts of the velocity gradient tensor, as well asproducts of these tensors. Those components that are con-nected with the antisymmetric part of the velocity gradienttensor scale with the third power of the shear parameter, sug-gesting that these effects cannot be captured perturbatively.Given that Sc = 1 in all our runs, we always have Re = Pe,which is at most about 100, so the inertial range of the turbu-lence is not very big yet. It is therefore important to investigatethe dependence of the various transport coefficients on the val-ues of Re and Pe, as was done in Fig. 3. The results availableso far suggest that the first three coefficients ( κ t , κ S , and κ SS )do not change with Re for Re > . If there were indica-tions that the resulting coefficients change beyond Re = 100 ,it would be important to make an effort to increase the valuesof Re even further. This would require more resolution and isobviously expensive. In view of the constancy of the first threecoefficients, this may not be well justified. The fourth coeffi-cient ( κ A ) seems to tend to zero, and the fifth one ( κ AS ) showslarge error bars. The situation regarding these last two coeffi-cients may not improve significantly towards larger Reynoldsnumbers, unless the simulations are run for long enough time.In general, turbulent transport tends to be enhanced in thedirection of the shear, i.e. κ yy tends to be larger than κ xx and κ zz . Furthermore, κ zz tends to be suppressed relative to κ xx .This is a result that is not reproduced by a simple analyticalclosure in which triple correlations are being replaced withquadratic ones. In particular, there is no evidence for a sup-pression of turbulent transport in the cross-stream or x direc-tion. Instead, there is a suppression in the spanwise directionout of the plane of the shear flow.We recall that the moduli of the diagonal components of theturbulent diffusivity tensor are found to decrease with increas-ing wave number of the mean concentration in a Lorentzianfashion. This is in agreement with earlier findings both inthe contexts of mean-field electrodynamics with and withoutshear [32, 33], as well as passive scalar transport in the ab-sence of shear [10]. The limit of high wave numbers maybe of interest for subgrid scale modeling in large-eddy sim-ulations of turbulence. However, it still needs to be clarifiedwhether the effective diffusivity is proportional to the inverseNyquist wave number to the second power, as suggested byour current results, or to some smaller power, ∼ k − / , asexpected for Kolmogorov turbulence. In order to address thisquestion, simulations at larger Peclet and Reynolds numbers are required. Such simulations do not require the presence ofshear. This is however beyond the scope of the present paper.Finally, we note that, in shear flows, the passive scalartransport properties are not affected by the presence of he-licity. In other words, there is no evidence for the existenceof components to the turbulent diffusivity tensor κ ij that areproportional to W i U j and W j U i . Acknowledgments
We thank Alexander Hubbard and Karl-Heinz R¨adler forsuggestions and stimulating discussions. We acknowledge theuse of computing time at the Center for Parallel Computers atthe Royal Institute of Technology in Sweden. This work wassupported in part by the European Research Council underthe AstroDyn Research Project No. 227952 and the SwedishResearch Council Grant No. 621-2007-4064. [1] G. Falkovich, K. Gawedzki, and M. Vergassola, Rev. Mod.Phys. , 913 (2001).[2] B. Shraiman and E. D. Siggia, Nature , 639 (2000).[3] P. Terry, Rev. Mod. Phys. , 109 (2000).[4] K. Burrell, Phys. Plasmas , 1499 (1997).[5] E. Kim, Astron. Astrophys. , 763 (2005).[6] N. Leprovost and E. Kim, Astron. Astrophys. , 617 (2006).[7] A. Brandenburg, P. K¨apyl¨a, A. Mohammed, Phys. Fluids ,1020 (2004).[8] M. Schrinner, K.-H. R¨adler, D. Schmitt, M. Rheinhardt, U.Christensen, Astron. Nachr. , 245 (2005).[9] M. Schrinner, K.-H. R¨adler, D. Schmitt, M. Rheinhardt,U. Christensen, Geophys. Astrophys. Fluid Dynam. , 81(2007).[10] A. Brandenburg, A. Svedin, G. M. Vasil, Mon. Not. R. Astron.Soc. , 1599 (2009).[11] L. L. Kitchatinov, G. R¨udiger, V. V. Pipin, Astron. Nachr. ,157 (1994).[12] A. Brandenburg, Astrophys. J. , 824 (2001).[13] N. E. L. Haugen, A. Brandenburg, and W. Dobler, Astrophys.J. , L141 (2003).[14] J. Wisdom, S. Tremaine, Astronom. J. , 925 (1988).[15] J. F. Hawley, C. F. Gammie, S. A. Balbus, Astrophys. J. ,742 (1995).[16] T. Elperin, N. Kleeorin, and I. Rogachevskii, Phys. Rev. E ,016311 (2003).[17] P. J. K¨apyl¨a, D. Mitra, and A. Brandenburg, Phys. Rev. E ,016302 (2009).[18] J. Smagorinsky, Monthl. Weather Rev. , 94 (1963).[19] H. K. Moffatt, Magnetic field generation in electrically con- ducting fluids . Cambridge University Press, Cambridge (1978).[20] F. Krause and K.-H. R¨adler,
Mean-field magnetohydrodynamicsand dynamo theory . Pergamon Press, Oxford (1980).[21] S. I. Vainshtein, L. L. Kitchatinov, Geophys. Astrophys. FluidDynam. , 273 (1983).[22] N. Kleeorin, M. Mond, I. Rogachevskii, Astron. Astrophys. , 293 (1996).[23] K.-H. R¨adler, N. Kleeorin, I. Rogachevskii, Geophys. Astro-phys. Fluid Dynam. , 249 (2003).[24] A. Brandenburg and K. Subramanian, Phys. Rep. , 1 (2005).[25] Ogilvie, G. I., Mon. Not. R. Astron. Soc. , 969 (2003).[26] P. Garaud, G. I. Ogilvie, J. Fluid Mech. , 145 (2005).[27] Liljestr¨om, A. J., Korpi, M. J., K¨apyl¨a, P. J., Brandenburg, A.,& Lyra, W., Astron. Nachr. , 92 (2009).[28] E. G. Blackman and G. B. Field, Phys. Fluids , L73 (2003).[29] K.-H. R¨adler, M. Rheinhardt, Geophys. Astrophys. Fluid Dy-nam. , 11 (2007).[30] A. Brandenburg, K. Subramanian, K., Astron. Astrophys. ,835 (2005).[31] N. I. Kleeorin, I. V. Rogachevskii, and A. A. Ruzmaikin, Sov.Phys. JETP , 878 (1990).[32] A. Brandenburg, K.-H. R¨adler, and M. Schrinner, Astron. As-trophys. , 739 (2008).[33] D. Mitra, P. J. K¨apyl¨a, R. Tavakol, and A. Brandenburg, Astron.Astrophys. , 1 (2009).[34] http://pencil-code.googlecode.com/, 1 (2009).[34] http://pencil-code.googlecode.com/