Mean-field closure parameters for passive scalar turbulence
Direct numerical simulations of isotropically forced homogeneous stationary turbulence with an imposed passive scalar concentration gradient are compared with an analytical closure model which provides evolution equations for the mean passive scalar flux and variance. Triple correlations of fluctuations appearing in these equations are described in terms of relaxation terms proportional to the quadratic correlations. Three methods are used to extract the relaxation timescales from direct numerical simulations. Firstly, we insert the closure ansatz into our equations, assume stationarity, and solve for . Secondly, we use only the closure ansatz itself and obtain from the ratio of quadratic and triple correlations. Thirdly we remove the imposed passive scalar gradient and fit an exponential decay law to the solution.
We vary the Reynolds (Re) and Péclet (Pe) numbers while keeping their ratio at unity and the degree of scale separation and find for large Re fair correspondence between the different methods. The ratio of the turbulent relaxation time of passive scalar flux to the turnover time of turbulent eddies is of the order of three, which is in remarkable agreement with earlier work. Finally we make an effort to extract the relaxation timescales relevant for the viscous and diffusive effects. We find two regimes which are valid for small and large Re, respectively, but the dependence of the parameters on scale separation suggests that they are not universal.
pacs:47.27.E-, 47.27.tb, 47.40.-x
Fluid flows in astrophysical bodies are most often highly turbulent. Direct modeling of such high-Reynolds-number flows is currently impossible. Consequently, greatly enhanced diffusivities or modified diffusion operators are often applied in simulations . Such models are challenging in terms of the required computational resources, so wide-ranging parameter studies cannot be performed.
An alternative approach is to separate the large and small scales, and derive equations for the former in which correlations of small-scale quantities are parameterized. This is usually referred to as mean-field theory, see e.g. [17, 13, 21, 22]. Various schemes have been introduced to close the equations for the correlations of small-scale quantities. In astrophysics, the second-order correlation approximation (SOCA), and the minimal approximation (MTA, see e.g. [1, 2]) are widely used. The relevant relaxation time in MTA has been determined numerically for passive scalar transport  as well as for the effect in mean-field electrodynamics [7, 8]. Another approach, related to the MTA, where a relaxation term is invoked to describe the higher order correlations has been introduced in . In this ‘Ogilvie approach’ several nondimensional coefficients are invoked to describe physically motivated parameterizations of the higher order correlations in terms of relaxation and isotropization terms. This model has been applied to different physical setups in order to calibrate the coefficients [9, 16, 14, 10].
The validity of the various approximations can in principle be tested by comparing with direct simulations in the same parameter regime. In practise this is often not easy due to the limited parameter range accessible by the simulations. The starting point for such studies has been isotropically forced homogeneous turbulence under the influence of rotation  and/or shear . In a recent work, the timescales related to diffusion and isotropization that appear in the Ogilvie approach have been studied . In the present study we extend the work of  to the passive scalar case.
2 Mean-field modelling
2.1 Ideal case
Let us consider the transport of a passive scalar under the influence of a turbulent fluid motion. For simplicity we assume a homogeneous, incompressible fluid and neglect at first diffusion and viscous dissipation. Then the governing equations for the concentration of the passive scalar, , and the fluid velocity read
where is the pressure, is the constant density and is a forcing function with (with the unit ‘force per mass’). Upon introduction of a Reynolds averaging procedure, indicated by an overbar, and are decomposed into mean and fluctuating parts, , . The fluctuating fields, represented by lowercase letters, are then governed by
where the prime indicates extraction of the fluctuating part, e.g., . Simplifying further, we stipulate the absence of a mean velocity and assume that the forcing has no mean part, i.e., . In the present case, the goal of mean-field modeling consists in deriving a closed equation for the mean concentration . From (1) and (3), together with we obtain directly
with the mean density of the passive scalar flux, . So the task of closing (5) reduces to representing by the mean concentration . In the standard mean-field approach, (3) is solved for a prescribed fluctuating velocity , usually under some simplifying assumptions which inevitably limit the applicability of the obtained results. The solution is employed to express in terms of . Alternatively, one can abstain from deriving such an explicit solution for the fluctuating concentration and instead strive for establishing an evolution equation for which of course again has to be closed in the sense that the only variables occurring are the mean quantities and themselves. Such an equation is obtained by multiplying (3) with and (4) with , summing up and averaging, arriving at
By virtue of the incompressibility of the fluid the fluctuating pressure can be expressed by the velocity fluctuations:
which for an infinitely extended medium and vanishing at infinity is readily solved by
Now we can conclude that the second, third, and fourth terms on the r.h.s. of (6) are quadratic in and linear in , hence represent third order correlations. Following  we introduce here the closure assumption
with a relaxation time . Upon further neglect of the correlation , we arrive at
which governs the evolution of . Here, stands for the Reynolds stress tensor. Since a passive scalar would not act back onto the velocity, can here be considered given and the closure is completed. Nevertheless, in the presence of rotation, shear, or gravity, (4) contributes quadratic correlations to (6) even in the kinematic case; see, e.g., [12, 24, 15]. Equation (10) is similar to the penultimate row of Eq. (53) in  when replacing temperature perturbation by and neglecting the buoyancy term. Note that for small , that is, fast relaxation, will follow the inhomogeneity in (10) almost instantaneously, hence
and we may interpret as a turbulent diffusivity tensor. For a discussion of its relationship to traditional mean-field results, see Sect. 2.5.
To facilitate further comparisons to [10, 4], where an additional evolution equation for the mean temperature perturbation variance is derived, we give here an analogous equation for , although it is not necessary for completing the closure:
We note in passing that the term becomes important in reacting flows . Setting
with another relaxation time , the closed equation reads
and we have full analogy to the last row of Eq. (53) in .
Until now we have not constrained the properties of the turbulence, in particular we have not required isotropy or homogeneity. For example, inhomogeneous turbulence could be thought of giving rise to position-dependent relaxation times. However, from a strict point of view, the ansatz (9) is only consistent with an isotropic or uni-axial turbulence where the preferred direction of the latter coincides with the direction of . Consequently, turbulent properties of must not change along any other direction. The same restrictions should of course hold for the concentration fluctuations , but this is in conflict with the presence of the second preferred direction in this turbulence. Hence, (9) can be strictly justified only under very specific circumstances.
For that reason and for the sake of simplicity, we specify now the mean as horizontal average, i.e., as average over all and . Consequently, all mean quantities depend on only and only the component of is relevant. If we further restrict to have at best a anisotropy, then there is only a single preferred direction, namely that of and the ansatz (9) is legitimate. The system of mean field equations then simplifies to
where , and could still depend on and . Assuming now further homogeneous and statistically stationary fluctuations and and a uniform gradient of , , a stationary regime of (15) should be given by
Let us now assume that in a direct numerical simulation (DNS) Eqs. (1) and (2) with an appropriately defined forcing and an imposed uniform are integrated in time until a statistically stationary regime is established. Extracting now all mean quantities occurring in (16) from the numerical solution and assuming validity of the model (10) and (14), it is obviously possible to determine the crucial relaxation times from such runs (method M1). On the other hand, should of course also obey their defining relations (9) and (13). The third-order correlations , , and are again accessible in the DNS results and open up an independent path for determining the relaxation times (method M2). At the same time, it can also be checked to what extent the neglect of is justified.
Another approach to extract is available from decay experiments, for which, after having reached a stationary state in the DNS, the imposed gradient of is switched off. Then, according to (10) and (14), and should decay uniformly in space and exponentially in time with the increment and , respectively, and can be identified with the decay rates measured in the DNS (method M3). An exemplification of this method is shown in Figure 1.
The goal of this paper consists in systematically testing the validity of the presented closure assumptions for a range of Reynolds and Péclet numbers as well as different levels of scale separation. From this we expect hints with respect to the validity of the Garaud model of turbulent convection .
2.2 Non-ideal effects
Admitting now diffusion and viscous dissipation, we have to add the term with the diffusivity on the r.h.s. of (1) and with the kinematic viscosity on the r.h.s. of (2). Consequently, in the evolution equation (6) for , the additional terms and occur on the r.h.s.. Rewriting their sum as
or more symmetric, as done in , as
does nevertheless not allow a representation entirely by the mean flux. Even in the (very particular) case the second terms of (17) and (18) remain. As a skyhook, the second and third terms are replaced by the -ansatz-like expression, although they contain second order rather than third order correlations. Analogously, on the r.h.s. of (12), diffusion requires a term
and the second term is replaced by . Note that diffusion of and modelled in this way is obviously determined by the molecular (or microscopic) diffusivities.
In astrophysical applications the deviation from ideal conditions is usually small, and quantities expressing this smallness are given by the Reynolds and Péclet numbers, Re and Pe, which reflect the strength of advection relative to diffusion:
where is a characteristic scale of the turbulence. We will further make use of the Schmidt number .
2.3 Summary of method M2
2.4 Scaling of the relaxation times
For the relaxation times some reasonable scaling assumptions are in order and we follow essentially the choices of : , belonging to third-order correlation terms, are expressed as , and , belonging to diffusive second-order correlation terms are written as
respectively. The first of these expressions seems to be appropriate only for , hence in general the scaling ansatz should read instead
Here is the smallest wavenumber consistent with the box size, . The crucial question is now: are the constants universal, at least for a given type of turbulence, and in particular, are they independent of the dimensionless numbers of the problem, i.e., Re and Pe, and of the degree of scale separation? A preliminary answer to this question was given in , where the timescales were found to show a slightly increasing trend with increasing scale separation (see their Fig. 4).
Methods M1 and M3 for determining the relaxation times described in Sect. 2.1 have now to be modified in the following way: In (16) we have to replace by and by . Both methods then deliver only these aggregates and we have to employ the different scalings of the relaxation times to figure out the individual constants . In contrast, method M2 has merely to be extended to include also the additional second-order correlations showing up in (18) and (19), that is, to use expressions (21) and (22).
2.5 Comparison with traditional results
from which, under the assumption of good temporal scale separation,
can be concluded. can readily be identified as turbulent diffusivity tensor. Here the correlation time is just defined by the last identity. This clearly resembles the result (10) with being identified with the correlation time , the more so as for the validity of both (26) and (10) the relevant time parameter has, in a sense, to be small.
Relaxing the assumption of good temporal scale separation, that is retaining (25), we observe the presence of the so-called memory effect , that is, the influence of at earlier times on the mean flux at time by virtue of a convolution. Performing a Fourier transform with respect to time, this convolution turns into a simple multiplication and we can write
with a frequency-dependent turbulent diffusivity tensor . This quantity is directly accessible to the test-field method as described, e.g., in . Based on numerical simulations, and without resorting to SOCA, it has been found that for homogeneous isotropic turbulence a satisfactory approximation is accomplished already by
where is the turbulent diffusivity for stationary fields and is independent of . A slightly better fit is provided by
with a constant . Turning back to the physical space the first approximation (28) is equivalent to
For the general case of imperfect scale separation both with respect to space and time, we refer here to , albeit this work deals with the mean electromotive force of MHD rather than with the mean flux of passive scalar transport. In that work, non-locality due to imperfect spatial scale separation shows up in the form of a diffusion term in the evolution equation for with a diffusivity occasionally even larger than, but of the order of the SOCA estimate of the turbulent diffusivity in the high-conductivity limit, . Clearly, this value can be very different from the molecular diffusivity. When comparing with the diffusion term for identified in (17) or (18) where only the microscopic diffusivities occur we have to state that the Ogilvie approach deviates in this respect significantly from what we expect from the traditional approach. To reconcile them, possible diffusion terms and had to be taken into account in the parameterizing ansatzes.
2.6 Significance of method comparisons
Let us finally discuss what is really “tested” by comparisons of the results of the different methods M1–M3. For an incompressible fluid with the specific conditions of our model and under the assumption that no higher than the first order temporal derivative occurs, an ansatz for analogous to (16),
is exhaustive, as
and hence is quite generally a linear and homogeneous functional of (even when the correlation is not neglected), and
and hence are spatially constant, both in the statistically steady state and during the decay of . (That is why the diffusive term is absent.)
Note, however, that as in (16) is an assumption because a contribution proportional to or can also be provided by the triple correlations, the diffusive terms or by .
Since the passive scalar does not influence the turbulent velocity, the coefficients are completely determined by and hence true constants. Consequently, any comparison of the methods M1, M2, and M3 tests the influence of
the weak compressibility of the fluid in our simulations, and
deviations of the simulated velocity turbulence from homogeneity, isotropy and statistical stationarity.
Note, that the first influence can in principle be made arbitrarily small by increasing the sound speed in the numerical model, likewise the second by increasing scale separation and extending time ranges for averaging.
A comparison of M1 and M3 tests in addition the justification of
the neglect of higher temporal derivatives of ,
the assumption which was employed in calculating from the steady state.
On the other hand, a comparison of M1 and M2 tests again the assumption and specifically to what extent the neglect of the correlation is legitimate. This affects only, so we expect that here the discrepancies between M1 and M2 are more pronounced than in .
With respect to the evolution equation for in (15), the same conclusions hold true as far as only the steady state and the free decay with are taken into account. However, considering a general transient as a consequence of switching between two non-zero constants, a richer behavior should appear which perhaps allows further reaching conclusions.
3 Numerical setup
In order to take benefit of the capabilities of the Pencil Code 111freely available at http://pencil-code.googlecode.com/ we solve
instead of the incompressible system (1) and (2) the corresponding equations for a compressible, but isothermal fluid
where we employ the pseudo enthalpy instead of the density; is the constant speed of sound and the trace-less rate-of-strain tensor . Interpretation of the results of such simulations in terms of the incompressible model of course requires to keep the Mach number small, typically . Then it is particularly justified to replace the correlation , included in the ansatz, by .
The equations are solved by equidistant sixth order finite differences in space and an explicit third order time stepping scheme with step size control for stability. The computational domain is a cube with dimension and grid resolutions ranging from to according to the requirements raised by the values of the Reynolds and Péclet numbers and the forcing wavenumber. Boundary conditions are periodic throughout. The fluctuating force is specified such that it generates an approximately homogeneous, isotropic and statistically stationary fluctuating velocity . During any integration timestep, is a frozen-in linearly polarized (i.e., non-helical) plane wave with a wave-vector which is consistent with the periodic boundary conditions and whose modulus is close to a chosen average value . The wave amplitude is kept fixed whereas the wavevector is randomly changing between time steps and hence is approximately -correlated in time. For further details, see .
In general, the parameter space is spanned by the dimensionless numbers (degree of scale separation), Re and Pe, in whose definitions (20) we specify the characteristic length for simplicity by . In the following, we will however restrict ourselves to , that is, and leave the more general cases to future work. By this, we avoid in particular the complication with the scaling ansatz (23) which occurs for . For methods M1 and M2, all averaged quantities were derived from the simulations by performing in addition to the averaging and a temporal averaging over an interval in which in particular the correlation was found to be statistically stationary. Statistical errors were estimated by dividing this interval into three equally long parts and calculating averages over each of them. These individual averages were compared to the average over the whole interval, and the largest deviation was taken for the error estimate.
We performed a number of simulations, and extracted and using methods M1, M2 and M3. The results are summarized in Table 1, where the several simulations are grouped together into different sets by the values used for the scale separation . Within each set the Reynolds numbers were varied by changing (). The timescales and listed in Table 1 are also illustrated in Figures 2 and 3, respectively. We see that all three methods give quite similar results, with
There are some exceptions, however. When , methods M2 and M3 yield , unlike method M1 which stays in range given above. Moreover, although method M2 always yields positive values for , those for are sometimes negative. The absolute values of the M2 results can also be very high. This is because the sum of the correlations used to calculate may become very small. This problem manifests itself mainly in the high Reynolds number runs.
4.1 Universality of the closure ansatz
Multiplying by the viscous time yields
where a tilde indicates normalization by . Hence, when considering and as functions of Re (or Pe), the wanted parameters should be obtainable by a linear regression analysis. Figure 4 shows both functions (39) for different values of . Obviously, a linear relation is clearly present both for large and small values of Re, but with very different fit parameters for the two ranges.
Guided by these functional dependencies we hence propose as an alternative for (39)
and analogously for . This ansatz allows to model linear dependencies on Re both for small and large arguments, but with different coefficients:
where the slopes may well be different in sign. The constants in (40) were determined by a standard fitting procedure and are given in Table 2. For the fit is surprisingly good throughout, whereas for this holds true merely for . For larger scale separation the ansatz fails to model the pronounced minima at visible in Figure 4 (right panel).
The slopes for high Re, that is, and , show a possible saturation with growing scale separation , but the other coefficients do not. Hence, universality of the ansatz (40) with respect to scale separation is questionable.
4.2 Comparison of methods
Let us next compare the relaxation timescales and obtained with the different methods in more detail.
4.2.1 Methods M1 and M3.
Figure 5 shows the ratio of the values of and determined by either of the two methods in dependence on Re and . The M3 values differ from the M1 ones by up to %. For the deviations generally diminish with decreasing Re and with growing scale separation, falling beneath 10% for . However, neither of these tendencies can be confirmed for .
In view of the discussion in Sect. 2.6 it is satisfactory to find improving agreement between M1 and M3 with increasing scale separation for which our forced turbulence is more and more approaching the desired target of isotropic stationary turbulence.
4.2.2 Methods M1 and M2.
As can be seen in Figure 6, for small scale separations the values of from both methods coincide fairly. Deviations lie within errors. For large , , 10, however, we find the deviations grow with falling Re, being still within errors around . We have to conclude that the neglect of has its strongest effects for low Re and high . In contrast, the differences between the values from methods M1 and M2 are much smaller, reaching a significant magnitude (exceeding errors ) only for the higher and . A possible reason for this is insufficient numerical resolution.
5 Comments and extensions
5.1 Alternative scaling
As an alternative to (38) one might consider
Then by multiplying with the dynamic time we arrive at
where the normalization is now with respect to . Figure 7 shows the same results as Figure 4, but now with the altered scaling. Again, a linear fit is viable on each curve, but only separately for low and high Re. The corresponding fit parameters can be found in Table 3. Unfortunately, an overall fit analogous to (40) does here not work satisfactorily.
|s||small Re||large Re|
5.2 An Ogilvie approach for compressible hydrodynamics?
One could be tempted to treat the system (36) and (35) in the spirit of the Ogilvie approach quite analogously to what was shown in Sect. 2.1. In the ideal case and with , we have for the fluctuating fields
where the second order correlation can only partly be expressed by . The remaining parts could be modelled by a ansatz as used for the diffusion terms in Sect. 2.2, but note that here the “diffusivity” is and we have no argument to consider the not properly modelled terms as small.
The main conclusion to be drawn from the present work is that the time scales used to model closure terms in the equations for the mean flux and the mean square concentration are nearly independent of Re for and also nearly independent of the scale separation ratio for . Expressed in terms of dynamical times, the resulting non-dimensional time scales can be referred to as Strouhal numbers whose values are around 3 for the closure term and around 7 for the closure term. The former value is in good agreement with earlier work using the approximation .
Equipped with this knowledge, we may now be better justified in using the closure hypotheses discussed here for the quantities and . On the other hand, as explained in the present paper, it is quite clear that these closure hypotheses lack thorough justification . One should therefore in future strive to find systematic discrepancies from the anticipated scalings. One example that we alluded to in the present paper is the inhomogeneous case in which the approach my break down. Future work in that direction seems now to be highly desirable, because in virtually all astrophysical applications the turbulence is inhomogeneous or at least anisotropic.
The authors acknowledge the hospitality of NORDITA during the program ‘Dynamo, Dynamical systems and Topology.’ We acknowledge the allocation of computing resources by CSC - IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education. This work was supported in part by the Finnish Cultural Foundation (JES), the Finnish Academy grants 121431, 136189 (PJK, MR) and 218159, 141017 (MJK), the European Research Council under the AstroDyn Research Project 227952, and the Swedish Research Council under Project 621-2011-5076.
-  E. G. Blackman and G. B. Field. New Dynamical Mean-Field Dynamo Theory and Closure Approach. Physical Review Letters, 89(26):265007, 2002.
-  E. G. Blackman and G. B. Field. A new approach to turbulent transport of a mean scalar. Physics of Fluids, 15:L73–L76, 2003.
-  A. Brandenburg. The Inverse Cascade and Nonlinear Alpha-Effect in Simulations of Isotropic Helical Hydromagnetic Turbulence. Astrophys. J., 550:824–840, 2001.
-  A. Brandenburg, N. E. L. Haugen, and N. Babkovskaia. Turbulent front speed in the Fisher equation: Dependence on Damköhler number. Phys. Rev. E, 83(1):016304, 2011.
-  A. Brandenburg, P. J. Käpylä, and A. Mohammed. Non-Fickian diffusion and tau approximation from numerical turbulence. Physics of Fluids, 16:1020–1027, 2004.
-  A. Brandenburg and Å. Nordlund. Astrophysical turbulence modeling. Rep. Progr. Phys., 74(4):046901, 2011.
-  A. Brandenburg and K. Subramanian. Minimal tau approximation and simulations of the alpha effect. A&A, 439:835–843, 2005.
-  A. Brandenburg and K. Subramanian. Simulations of the anisotropic kinetic and magnetic alpha effects. Astronomische Nachrichten, 328:507, 2007.
-  P. Garaud and G. I. Ogilvie. A model for the nonlinear dynamics of turbulent shear flows. J. Fluid Mech., 530:145–176, 2005.
-  P. Garaud, G. I. Ogilvie, N. Miller, and S. Stellmach. A model of the entropy flux and Reynolds stress in turbulent convection. Mon. Not. R. Soc., 407:2451–2467, 2010.
-  A. Hubbard and A. Brandenburg. Memory Effects in Turbulent Transport. Astrophys. J., 706:712–726, 2009.
-  P. J. Käpylä and A. Brandenburg. Lambda effect from forced turbulence simulations. A&A, 488:9–23, 2008.
-  F. Krause and K.-H. Rädler. Mean-field magnetohydrodynamics and dynamo theory. Pergamon Press, Ltd., Oxford, 1980.
-  A. J. Liljeström, M. J. Korpi, P. J. Käpylä, A. Brandenburg, and W. Lyra. Turbulent stresses as a function of shear rate in a local disk model. Astron. Nachr., 330:92–99, 2009.
-  E. J. M. Madarassy and A. Brandenburg. Calibrating passive scalar transport in shear-flow turbulence. Phys. Rev. E, 82(1):016304, 2010.
-  N. Miller and P. Garaud. A Practical Model of Convective Dynamics for Stellar Evolution Calculations. In R. J. Stancliffe, G. Houdek, R. G. Martin, & C. A. Tout, editor, Unsolved Problems in Stellar Physics: A Conference in Honor of Douglas Gough, volume 948 of Am. Inst. Phys. Conf. Ser., pages 165–169, 2007.
-  H. K. Moffatt. Magnetic field generation in electrically conducting fluids. Cambridge University Press, Cambridge, England, 1978.
-  G. I. Ogilvie. On the dynamics of magnetorotational turbulent stresses. Mon. Not. R. Soc., 340:969–982, 2003.
-  K.-H. Rädler and M. Rheinhardt. Mean-field electrodynamics: critical analysis of various analytical approaches to the mean electromotive force. Geophys. Astrophys. Fluid Dyn., 101:117–154, 2007.
-  M. Rheinhardt and A. Brandenburg. Modeling spatio-temporal nonlocality in mean-field dynamos. Astron. Nachr., 333:80–86, 2012.
-  G. Rüdiger. Differential rotation and stellar convection. Sun and the solar stars. Akademie-Verlag, Berlin, 1989.
-  G. Rüdiger and R. Hollerbach. The magnetic universe: Geophysical and Astrophysical Dynamo Theory. Wiley-VCH, Berlin, 2004.
-  J. E. Snellman, A. Brandenburg, P. J. Käpylä, and M. J. Mantere. Verification of Reynolds stress parameterizations from simulations. Astron. Nachr., 332:883–888, 2011.
-  J. E. Snellman, P. J. Käpylä, M. J. Korpi, and A. J. Liljeström. Reynolds stresses from hydrodynamic turbulence with shear and rotation. A&A, 505:955–968, 2009.