Reactive Rayleigh-Taylor systems: flame propagation and non-stationarity
Reactive Rayleigh–Taylor systems are characterized by the competition between the growth of the instability and the rate of reaction between cold (heavy) and hot (light) phases. We present results from state-of-the-art numerical simulations performed at high resolution in 2d by means of a self consistent lattice Boltzmann method which evolves the coupled momentum and thermal equations and includes a reactive term. We tune the parameters affecting flame properties, in order to address the competition between turbulent mixing and reaction, ranging from slow to fast-reaction rates. We also study the mutual feedback between turbulence evolution driven by the Rayleigh-Taylor instability and front propagation against gravitational acceleration. We quantify both the enhancement of flame propagation due to turbulent mixing for the case of slow reaction-rate as well as the slowing down of turbulence growth for the fast reaction case, when the flame quickly burns the gravitationally unstable phase. An increase of intermittency at small scales for temperature characterizes the case of fast reaction, associated to the formation of sharp wrinkled fronts separating pure hot/cold fluids regions.
Many natural and industrial processes involve fluid transport and mixing of passive or active scalar fields; examples include concentration fields of chemicals or biological species as well as the temperature field in natural convection. While many of these phenomena have been the subject of in-depth studies, the cases where chemical reactions are involved, presenting an even richer phenomenology, have received considerably less attention. We address the problem of the interplay of reaction and turbulent mixing in Rayleigh-Taylor (RT) systems (a situation occurring, for example, in thermonuclear burning of type Ia supernovae [1, 2, 3] or in the inertially confined nuclear fusion ) focusing on the different regimes which develop as we vary the ratio between the characteric time scales of underlying turbulence, , and the reaction time, . We limit this study to the case of single-step reaction, i.e. the two reactant scalar fields are distinguished by a reaction progress variable, proportional to the temperature (see Fig. 1). The two different temperatures in the hot and cold blobs of fluid of our numerical setup, mimick the combustion of a cold mixture of actual reactants into a hot mixture of burnt products [5, 6, 7]. The interesting point in this set-up is given by the natural competition between gravitational forces, which tends to mix the fluid and to produce a larger and larger mixing layer with uniform temperature, and combustion, which works against this mixing, trying to burn the whole volume and producing a propagating flame of given tickness and velocity. Moreover, the global phenomenology is complicated by the natural unsteadiness of the underlying RT problem. The Damkhler number, , is the natural control parameter and is identified by the ratio between the turbulent time scale, and the reaction time scale . Notice that because of RT unsteadyness the number depends on time:
where , as of standard RT phenomenology .
We perform highly resolved numerical simulations in 2d, with a resolution up to grid points (see table 1).
The 2d set up allowed us to reach a wide scale separation and a time-span large enough to address problems at both small and large Damkhler numbers, something still unfeseable in 3d. Our study has also a methodological motivation. We adopted a fully consistent thermal lattice Boltzmann method to evolve simultaneously the momentum equations and the advection-diffusion-reaction equations for temperature. The novelty here is to show that the method works well also in a non trivial case where the thermal modes are directly forced by the combustion terms.
The main result of the paper concerns the quantification of the front propagation due to turbulence enhancement for the slow reaction case, , and the clear signature of a strong feedback on the fluid evolution induced by the flame propagation when . In the latter case, we also measure an important increase of the temperature intermittency at small scales.
Equations of motion and numerical setup
We adopt a numerical scheme based on a recently proposed thermal lattice Boltzmann algorithm [8, 9], which is able to reproduce the correct thermohydrodynamics of an ideal gas with good numerical accuracy . To do that, the probability densities for a particle with velocity (belonging to a discrete set, with the index running over 37 values ) at space location and time evolve according to the lattice Boltzmann BGK equation [11, 12, 13]
the lhs stands for the free streaming of particles and the rhs represents the relaxation process towards Maxwell equilibrium with a characteristic time ( is the simulation time step). Once the density (), velocity () and temperature () fields are defined in terms of the lattice Boltzmann populations as
where is the thermal conductivity, the reaction time, the gravity and the specific heat at constant volume. Eqs ( ‣ Equations of motion and numerical setup) are obtained provided that the thermohydrodynamic fields appearing in the equilibrium density functions are properly shifted:
The novelty here is in the extra reaction term introduced for the temperature field 111This shift represents a kind of implicit equation, since, in principle, should be a function of the “real” thermodynamic temperature, which must be shifted itself ; however, we can overcome this problem, observing that for (always true in actual situations), and since the other shift is (hence ) we can safely assume that . .
In the third equation of ( ‣ Equations of motion and numerical setup) we have already subtracted the compression term to avoid effects due to a varying heat capacity or global heating of the system coming from steady increase of the underlying mean pressure.
The reaction rate must be zero in the pure phases, which we set at temperatures , for the hot fluid at the bottom, and , for the cold fluid on top (see figure 1), so that ; it must also transforms, irreversibly, the pure cold phase (unstable) into the hot one (stable).
A simple model for with these properties is given by a logistic-type expression
Results and discussion
Any RT system, even in the case of , will eventually reach the fast reaction limit, i.e. a situation where . This is due to the fact that the underlying turbulence slows down adiabatically, . As a consequence, sooner or later the flame tends to become active, burning at a rate faster than the turbulence stirring/mixing. Here we study the two regimes and and the transition between them.
Mean temperature profiles evolution
For large , the mixing is effective only at very small scales (where the characteristic times of the fluid motion are shorter), while the reaction tends to make uniform the mixed regions: as a result we get a topology of the temperature field which is made of “patches” separated by rather thin interfaces, which are smoother than the non-reacting RT case ; in addition, the front of the hot phase moves, on average, with a non zero mean drift velocity towards the top. These preliminary features can be better understood, at a pictorial level, looking at figure 2, where we show the magnitude of the temperature gradient at three different times in the evolution for the fastest reaction rate that we have studied (top panel), and compare it with the non-reacting case (bottom panel).
On the other hand, the larger the reaction time the closer is the phenomenology to the standard RT case: to see this we compare in figure 3 the evolution of the mean temperature profile
for the two extreme cases in our database, runs A and C: while for the evolution is basically undistinguishable from the usual RT dynamics[9, 10], in the fast rection case () the center of mass of the system clearly moves upwards, due to the burning processes, causing a shift –and an asymmetry– of the mixing region.
in figure 4 we plot the function for the three different reaction rates: the growth of in time is greatly enhanced when going from small to large .
Front propagation speed
If we integrate eqs. ( ‣ Equations of motion and numerical setup) over the whole volume, and divide by , we get an exact equation for the propagating front speed:
(where ) since the boundary terms vanish, owing to the periodic conditions on the lateral walls and to the adiabatic condition at top and bottom plates (). For the laminar flame (that is without gravity producing turbulence), the integral can be evaluated exactly (using, for instance, the usual hyperbolic tangent profile) to give an explicit expression for the speed, that is , where is the flame thickness: as the latter can be estimated to be , we end up with the well known result:
that is the flame propagates at constant speed.
We now ask what changes when turbulence sets is. In the small limit, when turbulence has the time to mix the fluids before reaction becomes active, we are in the so-called pre-mixed combustion. In this case, it has been conjectured [18, 20] that the simplest way to extend the result of the laminar case is to replace in expression ( ‣ Front propagation speed) the molecular diffusivity with an effective (turbulent) eddy diffusivity . If we use for the latter the dimensional estimate , where and are large scale characteristic velocity and length (in our case, e.g. the root mean square velocity and the mixing region length), and plug it into ( ‣ Front propagation speed), we get:
This prediction, probably valid to describe the evolution of slow flames in stationary turbulent flows is unlikely to be relevant for RT turbulence. The reason is that in order to observe a “eddy-diffusivity” driven flames one needs also scale separation between the turbulent eddies and the flame tickness, something that is not realized by the evolving RT system. On the other hand, we can rewrite ( ‣ Front propagation speed) exactly as:
where with we denote the fluctuations with respect to the mean vertical profile. It is clear now that for , the flame cannot have any strong influence on the underlying RT evolution and we can identify the first term on the rhs as the mixing layer length . Moreover, we know that in RT temperature fluctuations are almost constant in time and homogeneous inside the mixing layer , so also the second term on the rhs is proportional to the mixing layer extension. A natural prediction for is therefore:
In standard (stationary) turbulent reacting systems, one can check this prediction against experiments/simulations at various , obtained changing the reaction rate or the underlying turbulent intensity, while in our reacting RT setup we can exploit the fact that varies in time. In figure 5 we plot the front speed (normalized with the root mean square velocity) as a function of (which is itself a function of the simulation time) for the three runs.
As one can see, our prediction ( ‣ Front propagation speed) works satisfactorily in a wide range of , showing deviations only for very small times, where turbulence is not yet developed and the flame evolution is strongly influenced by the initial configuration, and for where it cannot be expected to be valid. In the latter case our data point flatten, as we clearly observe the feedback of the flame on the turbulent evolution, with a sort of synchronization between flame propagation and evolution of the turbulent kinetic energy toward a value where . Such a behaviour turns out to be in agreement with recent theoretical results obtained through a mean-field approach .
Small scale intermittency
When the reaction rate is fast (), there are no extended regions which are well mixed, since the cold material, as soon as it is slightly entrained through the hot one is suddenly burnt. As a result, the temperature field organizes in patches of pure reactants/products separated by sharp interfaces (being in the so called “segregated regime”), and, consequently, it has been conjectured that an increased intermittency develops at the small scales . The authors in  also derived a phenomenological prediction for the scaling laws of fluid temperature (and velocity) structure functions, according to which, in the asymptotics of , they should follow the relation
(where is the mixing length), irrespective of the order . From Eqn. ( ‣ Small scale intermittency) the expression for the flatness reads:
and so it increases with decreasing for all orders, a clear indication of strong small scales intermittency.
In figure 6 we show the growth of as function of the mixing length , for the three runs: the flatness for run A, corresponding to the smallest reaction time, is in good agreement, within error bars, with the prediction of equation ( ‣ Small scale intermittency), ; instead, at increasing , intermittency is depleted and the flatness grows more slowly, at a rate comparable (within error bars) with the non-reacting RT case, whose data are also reported for comparison.
We used a self consistent thermal lattice Boltzmann algorithm to perform numerical simulations of 2D Rayleigh–Taylor turbulence, in presence of chemical reactions between hot and cold fluids. The reaction was modelled by means of a Fisher-Kolmogorov-Petrovsky-Piskunov source term in the temperature equation; this term has been introduced by application of a suitable shift of the temperature field appearing in the equilibria of the lattice Boltzmann equation.
We analyzed the crossover among the various regimes emerging from the competition of turbulent mixing and reaction, going from the segregated () to the well mixed one (). We showed that, in the latter case, the effect of turbulence is to enhance the reaction front speed leading to an homogeneous burning in the whole mixing layer region. On the other hand, for moderate and large Damkhler, there is a feedback of the reaction on the statistical properties of the temperature field, resulting in increased intermittency at small scales in reasonable accordance with the prediction of .
We thank A. Celani and M. Cencini for useful suggestions and the QPACE development team for support during the implementation of our code and execution of the runs. We acknowledge access to the QPACE and eQPACE systems.
-  Zingale M., Woosley S.E., Rendleman C.A., Day M.S. & Bell J.B., Astrophys. J., 1021 (2005) 632
-  Khokhlov A.M., Astrophys. J., 449 (1995) 695
-  Gamezo V.N., Khokhlov A.M., Oran E.S., Chtchelkanova A.Y. & Rosenber R.O.,, Science, 299 (2003) 77
-  Freeman J.R., Clauser M.J. & Thompson S.L., Nucl. Fus, 223 (1997) 17
-  Cetegen B.M. & Kasper K.D., Phys. Fluids, 8 (1996) 2974
-  Tieszen S.R., Annu. Rev. Fluid Mech., 33 (2001) 67
-  Chertkov M., Lebedev V. & Vladimirova N., J. Fluid Mech., 633 (2009) 1
-  Sbragaglia M., Benzi R., Biferale L., Chen H., Shan X. & Succi S., J. Fluid Mech., 628 (2009) 299
-  Scagliarini A., Biferale L., Sbragaglia M., Sugiyama K. & Toschi F., Phys. Fluids, 22 (2010) 055101
-  Biferale L., Mantovani F., Sbragaglia M., Scagliarini A., Toschi F. & Tripiccione R., Phys. Fluids, 22 (2010) 115112
-  Bathnagar P.-L., Gross E. & Krook M., Phys. Rev., 94 (1954) 511
-  Wolf-Gladrow D., Lattice-gas cellular automata and lattice Boltzmann models (Springer, New York) 2000
-  Succi S., The lattice Boltzmann equation and its applications (Oxford Science publications) 2001
-  Fisher R., Ann. Eugen., 7 (1937) 355
-  Kolmogorov A.N., Petrovsky I.G. & Piskunov N.G., Bull. Moskovskogo Gosudartsvennogo Univ. Mat. Mekh., 1 (1937) 1
-  G. Goldrian et al., Computing in Science & Engineering, 10 (2008) 46
-  Biferale L., Mantovani F., Pivanti M., Sbragaglia M., Scagliarini A., Sschifano S.F., Toschi F. & Tripiccione R., Procedia Comp. Science, 1 (2010) 1069
-  Koudella C.R. & Neufeld Z., Phys. Rev. E, 70 (2004) 026307
-  Constantin P., Kiselev A., Oberman A. & Rhyzik L., Arch. Ration. Mech. Anal., 154 (2000) 53
-  Damkhler G., Z. Elektrochem. Angew. Phys. Chem., 46 (1940) 601
-  Chertkov M., Phys. Rev. Lett., 91 (2003) 115001
-  Brandenbur A., Erland L. Haugen N. & Babkovskaia N., Phys. Rev. E, 83 (2011) 016304