© ESO, 2013

1. Introduction

As many astrophysical measures are calibrated on stellar evolution theory, realistic stellar evolution models are crucial for astrophysics. One of the major issues is the influence of the macroscopic motions induced by the rotation of a star on its internal structure (see the recent review of Maeder & Meynet 2012). In particular, chemical mixing can significantly extend the lifetime of the stars by continuously providing fresh combustible matter to the core, thus feeding nuclear reactions for a longer time.

A key point of these stellar evolution models is the radial turbulent mixing of chemical species induced by radial differential rotation. According to Zahn’s model (Zahn 1974), this turbulent mixing is described by a diffusion coefficient that reads as (1)where κ is the thermal diffusivity, Ri_cr the critical value of the Richardson number, N_T the Brunt-Väisälä frequency, and rsinθdΩ/dr the shear associated with radial differential rotation. This expression is taken from Eq. (2.14) of Zahn (1992). While other versions of this diffusion coefficient exist in the literature, notably those accounting for additional effects like the μ-gradient (Maeder & Meynet 1996; Talon & Zahn 1997; Mathis et al. 2004), Zahn’s prescription is the basis of most implementations of the turbulent mixing driven by differential rotation in stellar evolution codes (Endal & Sofia 1978; Pinsonneault et al. 1989; Heger et al. 2000, see a detailed discussion in Meynet & Maeder 2000). Moreover, this mixing plays a significant role as compared to other transport processes, such as the meridional circulation combined with horizontal turbulence, especially for massive stars (Meynet & Maeder 1997, 2000; Denissenkov et al. 1999) and low-mass stars with strong μ-gradients (Palacios et al. 2003, 2006). It thus has a direct impact on the confrontations of stellar evolution models with observations. While chemical abundances at the star surface have been widely used for such confrontations (Michaud & Zahn 1998; Brott et al. 2011; Potter et al. 2012), asteroseismology is now providing direct constraints on these processes, especially for subgiant and giant stars (Deheuvels et al. 2012; Eggenberger et al. 2012).

The derivation of Zahn’s prescription is based on the hypothesis that turbulent flows are statistically stationary and that this is achieved when the mean shear flow is close to the linear marginal stability criterion. Following Townsend (1958), Zahn modified Richardson’s instability criterion to take the destabilizing effect of thermal diffusivity into account in such a way that marginal stability is controlled by RiPe ~ 1, where the Péclet number Pe compares the diffusive time scale with the dynamical one (as later confirmed through numerical linear stability analyses, e.g. Lignières et al. 1999). Then, he assumed that, in a turbulent state, the relevant Péclet number in the expression RiPe ~ 1 is the turbulent one , defined with the turbulent scales of velocity and length (u and ℓ, respectively). The diffusion coefficient is then estimated through the relation D = uℓ/3, which finally leads to Eq. (1).

The purpose of our present work is to use local 3D direct numerical simulations (DNS) to test Zahn’s modeling of turbulent mixing. Previous numerical simulations performed in a geophysical context have shown that a steady stably stratified homogeneous sheared turbulence can be obtained by imposing both a uniform mean shear and a uniform stratification and by tuning the Richardson number to reach statistical stationarity (Jacobitz et al. 1997). Such a flow is well adapted to studying chemical turbulent transport and enables one to relate the transport to the local gradients of velocity and entropy as in Zahn’s criterion. Imposing the shear and the stratification amounts to an assumption of scale separation between the variations in the mean and turbulent flows. The main difficulty of performing such numerical simulations in a stellar context is the very low Prandtl number of the stellar fluid. Indeed, this introduces a huge gap between the dynamical time scale and the diffusive one, which forces use of a large number of numerical time steps to accurately compute the effect of thermal diffusion on a dynamical time. This would require a prohibitive amount of computation time for DNS, where the whole range of scales of turbulence down to the dissipative scales are resolved. An attempt to test Zahn’s prescription was made by Brüggen & Hillebrandt (2001), but in this case the simulations relied on the numerical diffusivity of the ZEUS code that is strongly dependent on the grid size (Fromang & Papaloizou 2007). Fortunately, an asymptotic form of the Boussinesq equations in the limit of the small Péclet numbers exists, where the numerical constraint associated with the dynamical effect of the high thermal diffusivity disappears (Lignières 1999).

In this letter, we perform DNS of 3D steady, stably stratified, homogeneous sheared turbulence at decreasing Péclet numbers and use both the Boussinesq equations and the so-called small-Péclet-number approximation (SPNA) to explore the asymptotic Pe ≪ 1 regime. The mathematical model of the considered turbulent flow is presented in Sect. 2, followed by the numerical simulations (Sect. 3) and the methods used to study the turbulent transport (Sect. 4). Then, in Sect. 5, constraints on the diffusion coefficient obtained from our numerical simulations are described and compared to Zahn’s prescription. The results are discussed in Sect. 6.

	Fig. 1 Sketch of the flow configuration.
Open with DEXTER

2. Mathematical model of the flow

To obtain results that are as generic as possible, we consider the simplest flow likely to produce stably stratified, homogeneous, sheared turbulence. Its characteristics are a uniform vertical shear of horizontal velocity and a uniform vertical temperature gradient, as shown in Fig. 1. The governing equations of the flow are presented in Sect. 2.1, and the boundary conditions along with the forcing and the initial conditions in Sect. 2.2.

2.1. Governing equations

We use the Boussinesq equations in which density fluctuations are neglected, except in the buoyancy term. For the usually low Mach-number flows inside stars, this approximation is justified as long as the motion vertical length scale is smaller than the pressure scale height. Using L, ΔU, and ΔT = (T_max − T_min)/2 as length, velocity, and temperature units (see Fig. 1), their non-dimensional form reads as where is the velocity, the pressure deviation from hydrostatic equilibrium, the temperature deviation from the mean profile, and and forcing terms (see Sect. 2.2).

These equations show three non-dimensional parameters: (i) the Richardson number Ri = (N_T/S)², noting S = dU/dz = ΔU/L the shear rate, N_T² = αgΔT/L the Brunt-Väisälä frequency, and α the thermal expansion coefficient, (ii) the “macroscopic” Reynolds number Re = LΔU/ν = SL²/ν, where ν is the molecular viscosity, and (iii)  the “macroscopic” Péclet number Pe = LΔU/κ = SL²/κ.

To explore the regime of very low Péclet numbers, we use the SPNA (Lignières 1999). The basic principle of this approximation is to Taylor-expand all variables up to first order in Pe. Thus, Eqs. (3)and (4)become with . The background thermal stratification remains fixed because temperature deviations become infinitely small as Pe approaches zero. There are two non-dimensional parameters left: the Reynolds number and the “Richardson-Péclet” number RiPe, which compares the effect of the stratification modified by the thermal diffusion (with a time scale of τ_B²/τ_κ, where τ_B = 1/N_T and τ_κ = L²/κ are the buoyancy and diffusive time scales, respectively) to that of the shear. This number is known to control the linear stability of shear flows at low Péclet number (Lignières et al. 1999, and references herein).

An expansion similar to the SPNA has been proposed for an unstable thermal stratification (Thual 1992). However, it does not allow the feedback of the convective turbulent motions on the stratification, as observed in the solar convective zone.

The numerical code uses a Fourier colocation method in the horizontal directions, compact finite differences in the vertical, and a projection method plus Runge-Kutta for time advancing.

2.2. Boundary conditions, forcing, and initial conditions

Periodic boundary conditions apply in the horizontal directions, whereas in the vertical direction the fluid is bounded by two horizontal surfaces that are impenetrable (w = 0), stress-free in the spanwise direction (∂_zu = 0), uniformly sheared in the streamwise direction (∂_zv = S), and thermally controlled (θ = 0). The mean shear profile U(z) is imposed thanks to a body force applied at every time step Δt, where is the horizontal average of the streamwise velocity. The mean vertical temperature profile is imposed in the same way. With such a set-up, turbulence is not rigorously homogeneous in the vertical direction. Nevertheless, as already found in Schumacher & Eckhardt (2000) using similar forcing and boundary conditions, non-homegeneity is limited to thin layers (~10% of the domain vertical extent) close to the upper and lower boundaries.

As initial conditions, we construct a statistically homogeneous and isotropic Gaussian random velocity fluctuation field (Orszag 1969) with an energy spectrum proportional to k²e^−k2/k₀2, where k is the norm of the wave vector and k₀ corresponds to the maximum of the spectrum (Jacobitz et al. 1997). The initial temperature deviation θ is set to zero.

3. Numerical simulations

For each value of the Péclet number, the Richardson number can be tuned to obtain statistical steady-state flows. In the regime where thermal diffusivity has no dynamical effect, various studies (Holt et al. 1992; Jacobitz et al. 1997) had already shown there is a critical Richardson number such that the turbulent kinetic energy increases if Ri < Ri_cr and decreases if Ri > Ri_cr, as illustrated in Fig. 2. We could also find such critical Richardson numbers in the low-Péclet-number regime.

	Fig. 2 Evolution of the turbulent kinetic energy K for three different Richardson numbers at turbulent Péclet number Peℓ = 52. After a transient regime, the turbulent kinetic energy is stationary for Ri = 0.124, increases for Ri = 0.05 and decreases for Ri = 0.2.
Open with DEXTER

Moreover, to obtain reliable statistics of the turbulent transport, the simulation domain must contain several correlation length scales of the turbulence and the duration of the simulation must be long compared to its correlation time scale. Denoting ℓ as the horizontal integral length scale defined by (7)where E(k_h) is the horizontal spectrum of the turbulent kinetic energy averaged over the vertical direction, we find that the ratio ℓ/L_h never exceeds 0.16 in our simulations. This indicates that there are several large structures in the simulated domain. For illustration, Fig. 3 displays 2D snapshots of the temperature fluctuations and of the vertical velocity.

	Fig. 3 Snapshots of temperature fluctuations scaled by ΔT (left) and vertical velocity scaled by ΔU (right) in the yOz plane at Peℓ = 0.34. Vertical velocity presents smaller scales than temperature fluctuations; both are anti-correlated.
Open with DEXTER

To properly resolve all scales of turbulence, we use a resolution of 128² × 257 and an aspect ratio L_h/L_v = π/4. We have verified that in our simulation k_maxη > 1, where k_max is the largest wavevector present in the flow, η = (ν³/ε)^1/4 the Kolmogorov dissipation scale, and ε the dissipation rate. This criterion, which states that all scales down to the dissipation scale are resolved, has been found to be adapted for the study of turbulent transport in isotropic turbulence simulations (Gotoh et al. 2002). In addition, simulations performed at a higher resolution have confirmed that our resolution is sufficient.

4. Turbulent transport

The passive scalar is introduced once a statistical steady state is reached. We then consider and compare two complementary approaches to determine its vertical turbulent transport: one Lagrangian, by following fluid particles; the other Eulerian, by solving an advection/diffusion equation for a concentration field.

Owing to the stationarity and the spatial homogeneity of the turbulence, a fluid particle encounters statistically homegeneous conditions as it moves with the flow (as long as it keeps away from the upper and lower boundaries). The turbulence is thus homogeneous from a Lagrangian point of view. This is an important property because it enables us to apply Taylor’s turbulent transport theory (Taylor 1921). Accordingly, the transport is diffusive, and the vertical turbulent diffusion coefficient D_t reads as (8)where T_L is the Lagrangian correlation time and  ⟨  ⟩  the ensemble average over the particles. If z(t) denotes the vertical position of a particle, the vertical dispersion is then given by (9)We computed the turbulent diffusivity D_t using either Eq. (8)or a linear regression of the mean square displacement . In both cases the main source of error is the temporal fluctuations of the averaged turbulent quantities. Then, depending on the starting point of the time average, this creates a significant dispersion (up to 20%) in the values of D_t. We favor the turbulent diffusivity D_t obtained with Eq. (8)because in this case the dispersion is generally lower than using the linear regression of the mean square displacement.

In the Eulerian approach, the concentration field c is governed by the equation (10)where the molecular diffusivity D_m is such that Pe_c = LΔU/D_m = 10⁴. To determine D_t, we impose a mean concentration gradient dC/dz and compute (11)where c′ refers to the concentration deviation from the mean profile C(z).

We found that the relative difference between D_t computed from the Lagrangian approach (Eq. (8)) and the Eulerian one does not exceed 15%. This difference is partly due to molecular diffusion, which is present in the Eulerian approach but not in the Lagrangian one. Nevertheless, this effect is limited because D_m/D_t remains always less than 5%. Again, the difference is mostly due to temporal fluctuations of the turbulent quantities that generate dispersion in the time averages. A way to reduce this dispersion is to average over a longer time. While it is not possible in the Lagrangian approach since the displacements of particles are limited by the upper and lower boundaries of the domain, a much longer integration time can be used in the Eulerian one. We have thus used the Eulerian determination of the turbulent diffusion coefficient, keeping the error linked to the time average to a value lower than 5%.

5. Results

We present the results of five simulations corresponding to five different values of the turbulent Péclet number, Pe_ℓ = 52,0.90,0.72,0.34, and  ≪ 1 (the last using the SPNA). Table 1 displays the turbulent diffusion coefficients, together with other relevant physical parameters.

The most striking feature is that the full computation at Pe_ℓ = 0.34 and the SPNA computation give similar results. Indeed, at Pe_ℓ = 0.34 the critical Richardson number multiplied by the turbulent Péclet number is very close to the critical “Richardson-Péclet” number of the SPNA computation. This is strong evidence that at small Péclet number, the steady, stably stratified, sheared turbulence is characterized by a critical “Richardson-Péclet” number (RiPe_ℓ)_cr independent of the Péclet number. This regime has not been reached in the simulations from Pe_ℓ = 52 to Pe_ℓ = 0.7 where Ri_crPe_ℓ > (RiPe_ℓ)_cr.

Table 1 further shows the parameter β = D_t/(uℓ), which in the regime where RiPe_ℓ = (RiPe_ℓ)_cr, verifies (12)This parameter clearly has the same value for Pe = 0.34 and Pe ≪ 1. Thus, in the small-Péclet-number regime, the turbulent diffusion coefficient can be written in the form of (12)with a constant β. As a consequence, D_t is proportional to κRi^-1.

It is obvious that the final expression found for the vertical turbulent diffusion coefficient D_t (see Eq. (12)) has the same form as Zahn’s prescription given in Eq. (1), thus validating in the Pe ≤ 0.34 range the hypotheses made to derive it. The present numerical simulations also provide a quantitative determination of the critical turbulent “Richardson-Péclet” number (RiPe_ℓ)_cr = 0.426 and of the coefficient β = D_t/(uℓ) ≃ 0.131. This has to be compared with the order of magnitude estimates proposed by Zahn (1992), namely Ri_c ≲ 1/4 and β = 1/3 ≃ 0.333. Nevertheless, it is the proportionality constant between D_t and κRi^-1, i.e. the product of β and (RiPe_ℓ)_cr, which is relevant in practice in stellar evolution codes. We find a value of 5.58    ×    10^-2 for this quantity, whereas Zahn’s prescription leads to 8.33    ×    10^-2, which is a relative difference of 33%.

6. Conclusion and discussion

We thus conclude that the results of our DNS agree with the form of Zahn’s prescription and provide a first quantitative determination of D_t beyond Zahn’s order of magnitude estimate. In general, the main drawback of numerical simulations of turbulent flows is the unrealistically low value of the Reynolds number. While it may be interesting to test the dependence of our results on the Reynolds number, we note that estimates of effective turbulent Reynolds numbers associated with the radial diffusion of chemicals in stars lead to surprisingly low values Re_ℓ ~ 100 (Michaud & Zahn 1998), which are comparable to the values obtained in the present work.

Other relevant issues can be addressed numerically. In particular, the concentration field is considered here as a passive scalar so that the mean gradient has no effect on the flow, whereas in reality such a stratification would enhance the effect of temperature stratification. This effect is considered differently by Maeder & Meynet (1996) and Talon & Zahn (1997), and simulations might help in deciding between them. Another important aspect of the models of rotationally-induced mixing is the magnitude of the horizontal turbulent transport driven by the horizontal differential rotation. It controls the efficiency of the meridional circulation transport and the departure from shellular rotation (Zahn 1992). Again, the turbulent transport can be investigated with DNS similar to the present ones by varying the angle between the velocity and temperature gradients.

To study the turbulent transport in less idealized flow configurations would require a better understanding of the shear flows induced by the large-scale motions driven by rotation.

Acknowledgments

This work was granted access to the HPC resources of CALMIP under the allocation 2012–P0507.

References

All Tables

All Figures

	Fig. 1 Sketch of the flow configuration.
Open with DEXTER
In the text

	Fig. 2 Evolution of the turbulent kinetic energy K for three different Richardson numbers at turbulent Péclet number Peℓ = 52. After a transient regime, the turbulent kinetic energy is stationary for Ri = 0.124, increases for Ri = 0.05 and decreases for Ri = 0.2.
Open with DEXTER
In the text

	Fig. 3 Snapshots of temperature fluctuations scaled by ΔT (left) and vertical velocity scaled by ΔU (right) in the yOz plane at Peℓ = 0.34. Vertical velocity presents smaller scales than temperature fluctuations; both are anti-correlated.
Open with DEXTER
In the text