Buoyancy-driven flow through a bed of solid particles produces a new form of Rayleigh-Taylor turbulence

TL;DR

This study uses high-resolution simulations to explore how buoyancy-driven flow through a bed of solid particles alters Rayleigh-Taylor turbulence, revealing distinct mixing physics due to porosity.

Contribution

It introduces a new understanding of RT turbulence in porous media and proposes an effective continuum model accounting for porosity effects.

Findings

Mixing layer growth becomes linear in time for low porosity

Velocity fluctuations saturate instead of growing linearly

Porosity significantly alters classical RT turbulence behavior

Abstract

Rayleigh--Taylor fluid turbulence through a bed of rigid, finite-size, spheres is investigated by means of high-resolution Direct Numerical Simulations (DNS), fully coupling the fluid and the solid phase via a state-of-the art Immersed Boundary Method (IBM). The porous character of the medium reveals a totally different physics for the mixing process when compared to the well-known phenomenology of classical RT mixing. For sufficiently small porosity, the growth-rate of the mixing layer is linear in time (instead of quadratical) and the velocity fluctuations tend to saturate to a constant value (instead of linearly growing). We propose an effective continuum model to fully explain these results where porosity originated by the finite-size spheres is parameterized by a friction coefficient.

Equations

$$\frac{\partial {\mathbf{u}}_{\mathbf{f}}}{\partial t}+{\mathbf{u}}_{\mathbf{f}}\cdot \nabla {\mathbf{u}}_{\mathbf{f}}=-\nabla p+\nu {\nabla }^2{\mathbf{u}}_{\mathbf{f}}-\beta \mathbf{g}T+\mathbf{f}$$

$$\frac{\partial T}{\partial t}+{\mathbf{u}}_{\mathbf{cp}}\cdot \nabla T=\nabla \cdot ({\kappa }_{cp}\nabla T)$$

$$\mathcal{M}=\int 4\rho Y_1{Y}_2{d}^{3}x$$

$$\alpha =\nu \frac{45(1-\varphi {)}^2}{r_p^2{\varphi }^2}.$$

$$h$$

$$U$$

Full text

Introduction. – Rayleigh-Taylor (RT) turbulence is strongly influenced by physical phenomena such as rotation Carnevale et al. (2002); Baldwin et al. (2015); Boffetta et al. (2016) surface tension Chertkov et al. (2005), viscosity variations and/or viscoelastic effects Boffetta et al. (2010a); Zhou (2017a, b); Boffetta and Mazzino (2017). Little is know about RT turbulence, and buoyancy-driven turbulence in general, in porous media. This is in spite of the importance of the problem in a variety of environmental applications. Among the many possible examples, we mention here the geological storage of $CO_{2}$ in saline aquifers Michael et al. (2010) in order to mitigate the effects of emissions on climate changes, and all processes involving the injection of a hot fluid into a cooler, fluid-saturated, subsurface rock including the so-called thermal Enhanced Oil Recovery Obembe et al. (2016).

All these applications have renewed the interest to understand RT-induced mixing in porous media Huppert and Neufeld (2014); Slim (2014); Nakanishi et al. (2016); Kalisch et al. (2016); Binda et al. (2017). Accurate prediction of the performance of these processes requires a model describing the fully-coupled dynamics of the rock–fluid system. Our aim here is to propose a simple model of porous buoyancy-driven turbulent flow whose universal properties are extracted via a phenomenological theory. The simple model we study here shares with all real complex systems two key properties: i) the fluid motion is triggered by the RT instability; ii) the fluid motion evolves in a porous medium. We will show that the fluid-structure interaction problem radically changes the classical RT scenario giving rise to new physics which can be captured by simple theoretical arguments.

Here, we address the problem of RT turbulence in porous media by extensive numerical simulations of a fully resolved two-phase flow, representing a disordered distribution of solid particles (spheres) in the computational domain for different values of porosity. A state-of-the-art immersed-boundary method is employed to simulate the presence of the particles Breugem (2012); Lashgari et al. (2014); Picano et al. (2015). We find that the growth of the mixing layer is strongly affected by the presence of particles and, for sufficiently large concentrations, the mixing layer grows linearly in time. Velocity fluctuations are reduced and saturate to a constant value in the limit of large concentrations, with increasing anisotropy. The presence of particles also suppresses the turbulent heat transfer. The resulting phenomenology is in sharp contrast with the well-known quadratic growth rate of the mixing layer, accompanied by the linear growth in time of velocity fluctuations occurring in classical RT turbulence Boffetta and Mazzino (2017). Moreover, we compare the results of the fully resolved model with an effective continuous model in which the porosity is parameterized by a friction coefficient, a model for which simple theoretical predictions are possible, and find a good agreement with the results of the fully resolved model.

Model for porous RT turbulence. – We consider the Boussinesq model for the buoyancy-driven incompressible flow with velocity ${\bf u_{f}}({\bf x},t)$ and temperature $T({\bf x},t)$ in the presence of gravity ${\bf g}=(0,0,-g)$

[TABLE]

where $\nu$ is the kinematic viscosity of the fluid, $p$ the pressure, $\beta$ the thermal expansion coefficient and ${\bf f}({\bf x},t)$ is the immersed-boundary forcing that accounts for the presence of the particles. The temperature equation is solved in all the computational domain for both fluid and solid phases

[TABLE]

where $\bf u_{cp}$ and $\kappa_{cp}$ are the velocity and thermal diffusivity of the combined phase. These last two quantities can be expressed in a Volume of Fluids (VoF) formulation Hirt and Nichols (1981), based on the combined single phase values and on the local volume fraction. Since in the case considered here the particles do not move, i.e. the forcing term ${\bf f}$ does not depend on time, ${\bf u_{cp}}=(1-\xi){\bf u_{f}}$, where $\xi({\bf x})$ is a phase indicator field equals to [math] in the fluid phase and to $1$ in the solid phase. Similarly the combined thermal diffusivity is written as $\kappa_{cp}=(1-\xi)\kappa_{f}+\xi\kappa_{s}$ where $\kappa_{f}$ and $\kappa_{s}$ are the thermal diffusivity of the fluid and solid phases respectively Ardekani et al. (2018a). The computational domain contains a random distribution of $N$ solid spherical particles (obstacles) of macroscopic radius $r_{p}$. Particles are fixed in space and the no-slip and no-penetration boundary conditions on their surface are imposed indirectly via the forcing term ${\bf f}({\bf x})$ in (1). Further details on the numerical method can be found in the Supplemental Material Sup and in Ardekani et al. (2018a, b); Roma et al. (1999); Kempe and Fröhlich (2012).

The velocity and temperature fields in (1-2) are defined in a domain of volume $V=L_{x}\times L_{y}\times L_{z}$, with periodic conditions on the domain boundaries. The porosity of the domain, the ratio of the void volume over the total volume, is $\phi=1-NV_{p}/V=1-\langle\xi({\bf x})\rangle$ where $V_{p}=(4/3)\pi r_{p}^{3}$ is the volume of a single particle and $\langle\cdot\rangle$ represent the volume average.

We perform direct numerical simulations of Eqs.  (1-2) at different values of porosity. For simplicity and numerical convenience, the simulations described in this Letter assume $\kappa_{f}=\kappa_{s}=\kappa_{cp}=\nu$. The domain size has horizontal dimensions $L_{x}=L_{y}=32r_{p}$ and vertical height $L_{z}=128r_{p}$. A resolution of 16 points per particle diameter is used, giving a total of $N_{x}=N_{y}=256$ and $N_{z}=1024$ grid points on a regular grid. Numerical results are averaged over $4$ independent runs starting with different initial perturbations and are presented as dimensionless quantities using $L_{z}$ and $\tau=(L_{z}/Ag)^{1/2}$ as space and time units respectively.

The initial condition for Rayleigh-Taylor instability and turbulence is a layer of cooler (heavier) fluid over a warmer (lighter) layer at rest, i.e. $T({\bf x},0)=-(\theta_{0}/2)sgn(z)$ ($T=0$ is the reference temperature) and ${\bf u_{f}}({\bf x},0)=0$ where $\theta_{0}$ is the initial temperature jump which defines the Atwood number $A=(1/2)\beta\theta_{0}$. This initial condition is unstable and after the linear instability phase, the system develops a turbulent mixing zone that grows in time starting from the plane $z=0$ Boffetta and Mazzino (2017).

The phenomenology of the pure fluid case ($\phi=1$) is well known Boffetta and Mazzino (2017); Zhou (2017a). After the initial linear instability, the flow enters into a nonlinear phase where a turbulent mixing layer is produced and evolves in the vertical direction. The mixing layer amplitude can be defined in terms of the mean vertical temperature profile $\overline{T}(z,t)\equiv{1\over L_{x}L_{y}}\int T({\bf x},t)dxdy$ as the region of width $h$ in which $|\bar{T}(z)|\leq{\theta_{0}\over 2}r$ where $r<1$ is a threshold (typically $r=0.9$). In the turbulent phase the width of the mixing layer grows asymptotically as $h(t)=cAgt^{2}$, while vertical and horizontal velocity fluctuations grow linearly in time, with vertical fluctuations about two times larger than horizontal fluctuations and isotropic velocity gradients Boffetta et al. (2010b). The determination of the dimensionless coefficient $c$ has been the object of many numerical and experimental studies both in $3D$, where it is in the range $0.02-0.04$ Cabot and Cook (2006); Vladimirova and Chertkov (2009); Boffetta et al. (2010c); Zhou (2017a), and in $2D$ Young et al. (2001); Clark (2003); Celani et al. (2006).

Figure 1 shows a section of the temperature field for classic RT turbulence and a case with porosity $\phi=0.6$. Qualitative differences between the two cases are evident; in particular, the presence of strongly anisotropic, vertically elongated, plumes in the porous case. These differences are quantified in Fig. 2, where we plot the time evolution of the mixing layer $h(t)$ for different values of the porosity, starting from the standard case $\phi=1$. We observe that the presence of solid particles strongly reduces the growth of the mixing layer. While in the pure fluid case the mixing layer at late times follows the classical $t^{2}$ law Boffetta and Mazzino (2017), already for $\phi=0.8$ it shows a different scaling law and, for the smallest values of porosity $\phi=0.7$ and $\phi=0.6$, the growth becomes linear (see inset of Fig. 2). Moreover also the coefficient of the linear growth depends on the porosity. We notice that at short times, $t/\tau<0.5$, the presence of the particles has no effects on the evolution of $h(t)$ since the width of the mixing layer is here comparable with interparticle scale.

The reduced growth of the mixing layer is associated to the suppression of the turbulent velocity fluctuations. In Fig. 3, we show the horizontal $u_{x}(t)$ and vertical $u_{z}(t)$ rms velocities in the mixing layer. These are computed in a phase averaged sense as ${u_{x}}(t)=\langle({\bf u_{f}}\cdot{\hat{\bf x}})^{2}\rangle^{1/2}$, $\hat{\bf x}$ being the unit vector along the $x$-axis, and similarly for $u_{z}$, where brackets indicate average over the mixing layer. Figure 3 shows that both components are reduced in the presence of particles. For the smallest value of porosity, the velocity fluctuations become almost constant at large times. This is in agreement with the linear growth of the mixing layer observed in Fig. 2. We observe also a small increment of the anisotropy of the velocity components $u_{z}/u_{x}$ with respect to the case of pure fluid $\phi=1$, which is not surprising given the elongated structures observed in Fig. 1.

The growth of the mixing layer is a basic measurement of the amount of mass mixed by the turbulent flow. Recently, a more direct indicator of the mixed mass, ${\cal M}$, has been introduced which has the advantage of being a conserved inviscid quantity Zhou et al. (2016). It is defined by the integral

[TABLE]

where $\rho$ is the mixture density and the mass fractions, in the present case of a symmetric temperature jump, are $Y_{1}({\bf x})=(\theta_{0}/2-T)/\theta_{0}$ and $Y_{2}({\bf x})=(\theta_{0}/2+T)/\theta_{0}$.

While for the higher values of porosity ${\cal M}$ follows the $t^{2}$ behavior observed in the standard RT turbulence Zhou et al. (2016), for the lower values $\phi=0.7$ and $\phi=0.6$ it displays a clear linear behavior (see inset of Fig. 3).

Figure 4 shows the dimensionless turbulent heat transfer $Nu=1+\langle{\bf u_{f}}\cdot{\hat{\bf z}}\,T\rangle h/(\kappa_{f}\theta_{0})$ as a function of the Rayleigh number, defined for RT turbulence as $Ra=Agh^{3}/(\nu\kappa_{f})$. We observe large fluctuations for all the values of porosity, even after averaging over realizations. Nonetheless it is possible to observe a reduction of $Nu$, for given $Ra$, by decreasing the value of porosity. A similar behavior has been observed in the case of rotating Rayleigh-Taylor turbulence where the reduction of $Nu$ is produced to the decoupling of velocity and temperature fluctuations due to the bi-dimensionalization of the flow Boffetta et al. (2016). For large porosity, $\phi=1$ and $\phi=0.8$, the scaling is in agreement with the so-called ultimate state regime $Nu\simeq Ra^{1/2}$ already observed in the pure fluid case Boffetta and Mazzino (2017). For smaller values of $\phi$ there is a clear indication of a transition to a different regime compatible with a $Ra^{1/3}$ scaling. Indeed, assuming that $h\sim t$ and $u_{rms}\sim t^{0}$ we obtain $Nu\sim t$ and $Ra\sim t^{3}$ which imply $Nu\sim Ra^{1/3}$.

Interpretation in terms of an effective model – Let us now show that the features of the porous RT turbulence can be obtained by an effective continuous model (without particles), in which the porous medium is parameterized by a friction coefficient. The model is obtained by averaging the microscopic equation over a volume which includes many particles and therefore filters the discrete nature of the porous medium. In the limit of small particles, the porous medium is considered as an homogeneous fluid with an additional effective friction term $-\alpha{\bf u}$ added to the momentum equation (1) Breugem et al. (2006). The friction coefficient $\alpha$ is:

[TABLE]

We remark that the use of a continuous model for the problem discussed in this Letter is not justified a priori, since there is no large-scale separation between particles size and box size. Moreover, particles are not very small to guarantee the presence of a Stokes flow in the pores. This is why the fluid inertia contribution is retained. The continuous model can be corrected taking into account finite particle Reynolds number $Re_{p}=r_{p}u/\nu$ (where $u$ represents the magnitude of the flow velocity around the particle) by the factor $\left(1+{\phi\over 50(1-\phi)}Re_{p}\right)$ Breugem et al. (2006). For simplicity, in the following we consider the extension of the continuous model to the Boussinesq equations in the limit of small particles with linear friction (4) only and we find that it is able to reproduces many of the results of the full microscopic model and sheds light upon the mechanism at the basis of the results discussed in the previous Section.

In the limit of large porosity, $\phi\simeq 1$, the friction coefficient (4) vanishes and therefore we expect that the standard RT turbulence phenomenology holds. Therefore in this limit we can assume that $h\simeq\beta g\theta_{0}t^{2}$ and $U\simeq\beta g\theta_{0}t$. On dimensional grounds, by using these scaling laws, one sees that $\alpha{\bf u}$ becomes dominant over ${\bf u}\cdot{\bf\nabla u}$ in (1) after a time $t_{\alpha}\simeq 1/\alpha$. Therefore, for $t>t_{\alpha}$ we expect a different phenomenology given by the balance of the buoyancy term which injects energy and the friction term which removes the energy in the system. This balance gives the new scaling laws

[TABLE]

Therefore, already at the level of dimensional analysis, the effective model is able to reproduce the behavior observed in the fully microscopic model, i.e. the saturation of velocity fluctuations and the linear growth of the mixing layer.

In Fig. 5 we plot the time evolution of the rms of the horizontal and vertical velocities for three different simulations of the effective model: one for the standard RT with $\alpha=0$ and two with larger values of the friction coefficient. The case $\alpha\tau=3$ corresponds to the case $\phi=0.6$ according to (4) and will be used to make a quantitative comparison of the homogeneous model with the full microscopic model. As in the microscopic model, we observe that while for $\alpha=0$ the large-scale velocity grows linearly in time (after an initial transient), in the simulations with friction the velocity saturates to a constant value. Moreover, anisotropy increases with $\alpha$, as the horizontal velocity is suppressed more than the vertical one, a feature also observed in the microscopic model (see Fig. 3). In the inset, we report the vertical rms velocity multiplied by the friction coefficient $\alpha$, which, according to (6), gives a constant value independent on $\alpha$.

Figure 6 shows the evolution of the mixing layer width, $h(t)$, from the three simulations of the effective model at increasing values of the friction coefficient. The presence of friction slows down the growth of the mixing layer but also changes its slope. For the largest values of $\alpha$, the growth becomes linear (see inset of Fig. 6) in agreement with the prediction (5).

The case $\alpha\tau=3$ is also plotted in the inset of Fig. 2 for a direct comparison with the full model. It is evident that the simple homogeneous model is able to reproduce quantitatively the behavior and the transition observed in the microscopic model.

Conclusions. – We have numerically studied Rayleigh-Taylor turbulence in the presence of fixed macroscopic solid particles, for different values of porosity coefficient. We have shown that the presence of particles reduces the growth of the mixing layer, which for small porosity follows asymptotically a linear behavior. In this regime turbulent velocity fluctuations saturate to a constant value. We have interpreted these results in terms of a continuous homogeneous model with an additional linear friction term representing the effective porosity of the medium. Dimensional analysis predicts that the friction term at late times modifies the asymptotic growth of the mixing layer. This is confirmed by extensive simulations of the effective model which is shown to reproduce the main features observed in the full model.

Acknowledgments. G.B. acknowledges financial support by the project CSTO162330 Extreme Events in Turbulent Convection and from the Departments of Excellence grant (MIUR). HPC center CINECA is gratefully acknowledged for computing resources.