Turbulence via Intermolecular Potential: Viscosity and Transition Range of the Reynolds Number

Abstract

Turbulence in fluids is an ubiquitous phenomenon, characterized by spontaneous transition of a smooth, laminar flow to rapidly changing, chaotic dynamics. In 1883, Reynolds experimentally demonstrated that, in an initially laminar flow of water, turbulent motions emerge without any measurable external disturbance. To this day, turbulence remains a major unresolved phenomenon in fluid mechanics; in particular, there is a lack of a mathematical model where turbulent dynamics emerge naturally from a laminar flow. Recently, we proposed a new theory of turbulence in gases, according to which turbulent motions are created in an inertial gas flow by the mean field effect of the intermolecular potential. In the current work, we investigate the effect of viscosity in our turbulence model by numerically simulating the air flow at normal conditions in a straight pipe for different values of the Reynolds number. We find that the transition between laminar and turbulent flow in our model occurs, without any deliberate perturbations, as the Reynolds number increases from 2000 to 4000. As the simulated flow becomes turbulent, the decay rate of the time averaged Fourier spectrum of the kinetic energy in our model approaches Kolmogorov’s inverse five-thirds law. Both results are consistent with experiments and observations.

1. Introduction

Historically, the phenomenon of turbulence in fluids was noted by Leonardo da Vinci. However, its first scientifically documented account was due to Boussinesq [1], while the word “turbulence” itself was suggested by Thomson, Lord Kelvin [2]. In his famous experiment, Reynolds [3] demonstrated that a laminar flow of water spontaneously developed turbulent motions without any measurable external disturbances and that the onset of turbulence was reliably associated with a sufficiently high value of the Reynolds number. Later, Kolmogorov [4,5,6] found that the time-averaged Fourier spectrum of the kinetic energy of turbulent air flow decayed as the inverse five-thirds power of its wavenumber.

Despite overwhelming research efforts spanning multiple decades [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23], the phenomenon of turbulence in gases and liquids remains unexplained, namely, an adequate fluid-mechanical model does not exist, for either a liquid or a gas, where, at appropriate values of the Reynolds number, turbulent flow naturally emerges from laminar initial and boundary conditions in the absence of artificial external disturbances. As examples, one can refer to relatively recent works [24,25,26,27], where turbulent-like motions in a numerically simulated flow had to be created artificially by deliberate perturbations. In reality, turbulence emerges spontaneously by itself, even if all reasonable measures are taken to preserve the laminarity of the flow (e.g., the experiment of Reynolds [3])—moreover, it was precisely the spontaneity of the onset of turbulence that attracted world-wide scientific interest to this intriguing phenomenon.

A different approach to model turbulence is based on the concept of “eddy viscosity”, first proposed by Boussinesq [1]. This approach has been actively explored in numerical turbulence modeling, with methods such as Reynolds-Averaged Navier–Stokes equations (RANS) or large eddy simulation (LES) gaining popularity (refer, e.g., to Wilcox [28] or Ferziger and Perić [29] for a detailed exposition). However, evidence has emerged lately [30,31] that this concept fails to describe a freely developing turbulence.

In our recent works [32,33,34,35], we proposed a theory of turbulence in a gas, where turbulent motions in an initially laminar, inertial (that is, constant pressure) flow were created via the average effect of gas molecules interacting by means of their intermolecular potential $\varphi \left(r\right)$. In our theory, this average effect is expressed via the mean field potential $\overline{\varphi }$ (which, of course, depends on $\varphi \left(r\right)$), whose gradient enters the equation for the transport of momentum. According to our theory, in the absence of the pressure gradient, the effect of $\overline{\varphi }$ becomes the key driving force of turbulent dynamics. This novel effect is absent from the conventional Euler and Navier–Stokes equations of fluid mechanics because; in the Boltzmann–Grad limit [36], it is tacitly assumed that the effect of $\overline{\varphi }$ is negligible. However, we found that, in our model of inertial gas flow, $\overline{\varphi }$ produces turbulent flow with Kolmogorov decay of the Fourier spectra of the kinetic energy. Evidently, its effect is non-negligible and quantifiable.

Remarkably, Tsugé [37] attempted to explain the creation of turbulence via long-range correlations between molecules, but Tsugé’s result was restricted to an incompressible flow. On the other hand, from what we discovered thus far, it appears that density fluctuations are instrumental in the creation of turbulent dynamics. In our past work [38], we also considered long-range interactions as a possible reason for the manifestation of turbulence. However, we later found that even the hard sphere potential creates turbulence in our model, which means that a typical intermolecular potential is also capable of the same effect.

Thus far, our model of turbulent gas flow did not include viscous effects. In the current work, we equip the momentum transport equation of our model with the standard viscous term, which counteracts the effect of the mean field potential (i.e., $\overline{\varphi }$ creates turbulent motions, while the viscosity dissipates them). We numerically simulate the air flow at normal conditions within a straight pipe at different values of the Reynolds number and find that the transition between the laminar and turbulent flow in our model occurs within the same range of values of the Reynolds number as observed in practice [39,40]. Additionally, we find that, for turbulent values of the Reynolds number, the rate of decay of the time-averaged Fourier spectrum of the kinetic energy of the flow approaches the famous Kolmogorov decay rate of the inverse five-thirds power of the wavenumber.

The work is organized as follows. In Section 2, we present the inviscid model of inertial flow with the mean field potential and demonstrate that, at a low Mach number, the effect of the mean field potential in nondimensional variables is of the same order as the rest of the terms. In Section 3, we add viscosity into the momentum equation in a standard fashion. In Section 4, we present the results of a numerical simulation of the air flow at normal conditions in a straight pipe and show that the laminar-to-turbulent transition occurs as the Reynolds number increases from 2000 to 4000. Section 5 summarizes the results of this work.

2. Our Model of Inertial Turbulent Gas Flow

In the context of our theory [32,33,34], the turbulent flow of a gas with density $\rho$ and velocity $\mathit{u}$, at a constant pressure $p_0$ (inertial flow) and in the absence of viscous effects, is described by the following mass and momentum transport equations

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0,$$

(1)

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla ·\left(\rho {\mathit{u}}^2\right)+\nabla \overline{\varphi }=\mathbf{0}.$$

(2)

Above, observe that Equation (2) for the transport of momentum possesses a novel term $\nabla \overline{\varphi }$, which replaces the pressure gradient in the conventional Euler or Navier–Stokes equations. This term quantifies the average (or “mean field”) effect of the motion of molecules of mass m, interacting via a potential $\varphi \left(r\right)$. Therefore, in our preceding works we referred to $\overline{\varphi }$ as the mean field potential. For a general short-range intermolecular potential $\varphi \left(r\right)$, in our work [34] we computed $\overline{\varphi }$ via

$$\overline{\varphi }=\frac{2\pi p_0\rho }{3m}{\int }_0^{\infty }\left(1-e^{-\varphi \left(r\right)/\theta }\right)\frac{\partial }{\partial r}\left(r^{3}Y\left(r\right)\right)\mathrm{d}r,$$

(3)

where $\theta =p_0/\rho$ is the kinetic temperature, and $Y\left(r\right)$ is the pair cavity distribution function [41]. Under the assumption that the gas is sufficiently dilute (i.e., $Y\approx 1$) and that the intermolecular potential $\varphi \left(r\right)$ can be approximated via the hard sphere potential with the effective range $\sigma$,

$${\varphi }_{HS}\left(r\right)=\left\{\begin{array}{cc}0\hfill & r\ge \sigma ,\hfill \\ \infty \hfill & r<\sigma ,\hfill \end{array}\right.$$

(4)

the general formula in Equation (3) simplifies to

$${\overline{\varphi }}_{HS}=\frac{4p_0\rho }{{\rho }_{HS}},$$

(5)

where ${\rho }_{HS}=6$ m/$\pi {\sigma }^3$ is the density of the equivalent hard sphere of mass m and diameter $\sigma$. Substituting ${\overline{\varphi }}_{HS}$ from Equation (5) into Equation (2) yields

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla ·\left(\rho {\mathit{u}}^2\right)+\frac{4p_0}{{\rho }_{HS}}\nabla \rho =\mathbf{0}.$$

(6)

This novel effect is absent from the conventional Euler and Navier–Stokes equations of fluid mechanics because, in the Boltzmann–Grad hydrodynamic limit (see Grad [36], pp. 352–353), it is assumed that ${\sigma }^2/$m∼constant as m and $\sigma$ are taken to zero, and therefore ${\rho }_{HS}\to \infty$. However, it can be shown that, in the nondimensional variables, the effect of the mean field potential in the inertial flow is comparable to the rest of the terms at sufficiently low Mach numbers. Indeed, let us rescale the variables in Equations (1) and (6) in a standard fashion by introducing the reference values of spatial scale L, flow speed U, and density ${\rho }_0$:

$$\tilde{t}=\frac{U}{L}t,\tilde{\mathit{x}}=\frac{\mathit{x}}{L},\tilde{\rho }=\frac{\rho }{{\rho }_0},\tilde{\mathit{u}}=\frac{\mathit{u}}{U}.$$

(7)

In addition, we denote the packing fraction $\eta$ and the Mach number $\mathit{Ma}$ via

$$\eta =\frac{{\rho }_0}{{\rho }_{HS}},\mathit{Ma}=U\sqrt{\frac{{\rho }_0}{\gamma p_0}},$$

(8)

where $\gamma$ is the adiabatic exponent of the gas. In the nondimensional variables, Equation (1) for density remains the same, while Equation (6) for momentum becomes

$$\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial \tilde{t}}+\tilde{\nabla }·\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)+\alpha \tilde{\nabla }\tilde{\rho }=\mathbf{0},\alpha =\frac{4\eta }{\gamma {\mathit{Ma}}^2}.$$

(9)

For air, $\gamma =1.4$, and ${\rho }_{HS}=1850$ kg/m${}^3$ (see Equation (64) in our work [33] for details). Under normal conditions (sea level, 20 °C), we have ${\rho }_0=1.204$ kg/m${}^3$, and $p_0=101.3$ kPa. As a result, the packing fraction $\eta \approx 6.5·{10}^{-4}$. Additionally, taking $U=30$ m/s, as in the numerical simulations below, we obtain $\mathit{Ma}\approx 8.74·{10}^{-2}$. Combining the estimates, we arrive at $\alpha \approx 0.24$, for the settings of our numerical simulations below. Moreover, reducing the speed of the flow to 15 m/s leads to $\alpha \approx 1$. This confirms that the mean field effect of the intermolecular potential at normal conditions and relatively low Mach numbers is not negligible and certainly does not vanish as presumed in the Boltzmann–Grad limit. Additionally, we found via numerical simulations in [32,33,34] that the system of Equations (1) and (6) produces distinctly turbulent flow with Kolmogorov decay of the Fourier spectra of the kinetic energy with the mean field potential present and develops purely laminar flow with no signs of turbulence when the mean field potential is removed (compare Figures 2 and 3 in our work [34]). Evidently, the effect of mean field potential is non-negligible and quantifiable.

3. A Model of Inertial Turbulent Gas Flow with Viscosity

Due to the absence of viscosity, numerical solutions of Equations (1) and (6) always develop turbulent dynamics from an initially laminar flow [32,33,34]. In the current work, we add the standard viscous term with the dynamic viscosity $\mu$ into Equation (6):

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla ·\left(\rho {\mathit{u}}^2\right)+\frac{4p_0}{{\rho }_{HS}}\nabla \rho =\nabla ·(\mu \nabla \mathit{u}).$$

(10)

According to kinetic theory of gases [42,43], $\mu$ is proportional to the square root of the temperature. In a gas, the product of temperature and density is proportional to the pressure, which, in turn, is constant in an inertial gas flow. Therefore, in our setting,

$$\mu ={\mu }_0\sqrt{\frac{{\rho }_0}{\rho }},$$

(11)

where ${\mu }_0$ is the reference value of viscosity. In the nondimensional variables, Equation (10) becomes

$$\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial \tilde{t}}+\tilde{\nabla }·\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)+\alpha \tilde{\nabla }\tilde{\rho }=\frac{1}{\mathit{Re}}\tilde{\nabla }·\left({\tilde{\rho }}^{-1/2}\tilde{\nabla }\tilde{\mathit{u}}\right),$$

(12)

where the Reynolds number $\mathit{Re}$ is given via

$$\mathit{Re}=\frac{{\rho }_{0}UL}{{\mu }_0}.$$

(13)

It is clear that, in this model, turbulent flow cannot emerge if the coefficients of the forcing and dissipative terms are balanced, that is, $\mathit{Re}={\alpha }^{-1}\approx 4.1$ in our setting. However, for $\mathit{Re}\gg {\alpha }^{-1}$, the manifestation of turbulent flow depends on the geometry of the domain. In particular, it is known from observations and practical engineering knowledge [39] that the transition between laminar and turbulent gas flow in straight pipes occurs within the range of values of the Reynolds number between 2000 and 4000, with the parameters L and U in Equation (13) being the width of the pipe and the maximum speed of the flow, respectively.

We have to note that, in general, the presence of viscosity in a flow induces pressure variations due to viscous friction; yet, above we introduced viscosity into Equation (10) for the transport of momentum rather formally, without accounting for such an effect. Thus, the system of Equations (1) and (10) should not be viewed as a practical method for accurate prediction of a real-world gas flow. Instead, it should be treated as a “proof-of-concept” model, whose purpose is to investigate how the relation between the mean field potential forcing and viscous dissipation alone results in the manifestation of turbulence in a gas flow, and at which values of the Reynolds number such a transition occurs.

To confirm that our model is suitable for the above stated purpose, first observe that, in an established laminar solution of the conventional Navier–Stokes system, expressed in the nondimensional variables, the term with the pressure gradient, induced by the viscous friction, must be comparable to the viscous term itself (that is, ∼${\mathit{Re}}^{-1}$). In the numerical simulations below, we find that a laminar flow still develops for $\mathit{Re}=1000$, which is also the smallest value of the Reynolds number used. This, in turn, means that, for $\alpha \approx 0.24$ in Equation (12), the effect of the induced pressure gradient would be ∼240 times smaller than that of the mean field potential term and could, therefore, be disregarded.

4. Numerical Simulations

To investigate whether or not our model predicts the transition to a turbulent flow within a realistic range of values of the Reynolds number, here we present numerical simulations of an inertial air flow at normal conditions in a straight pipe. As in our recent works, we use the appropriately modified rhoCentralFoam solver [44], which uses the central discretization scheme of Kurganov and Tadmor [45], with the flux limiter of van Leer [46]. The rhoCentralFoam solver is a standard component of the OpenFOAM suite [47].

Here, we simulate the inertial air flow in a straight pipe of a square cross-section, using Equations (1) and (10), with $p_0=101.3$ kPa, and ${\rho }_{HS}=1850$ kg/m${}^3$. The size of the pipe is $36\times 5.2\times 5.2$ cm${}^3$. The domain is uniformly discretized in all directions with the step of $0.8$ mm, which comprises $450\times 65\times 65=$ 1,901,250 finite volume cells in total. The pipe is open-ended at the outlet side and has a wall at the inlet side, with the circular inlet of 1 cm in diameter located in the middle of this wall. The longitudinal section of the domain is shown in Figure 1.

The boundary conditions are as follows. The density is set to $\rho =1.204$ kg/m${}^3$ at the outlet and has zero normal derivative at the inlet and the walls. The velocity has zero normal derivative at the outlet, no-slip condition at the walls, and a radially symmetric parabolic profile at the inlet, with a maximum of 30 m/s in the middle of the inlet, directed along the axis of the pipe. The diameter of the inlet, and the speed of the entering flow, correspond to the experiment by Buchhave and Velte [48]. Initially, the gas inside the pipe is at rest, with zero velocity and uniform density set to 1.204 kg/m${}^3$.

We conducted numerical simulations of Equations (1) and (10) for the values of the Reynolds number $\mathit{Re}=$ 1000, 2000, 3000, and 4000, by setting the reference viscosity ${\mu }_0$ to $1.872·{10}^{-3}$, $9.36·{10}^{-4}$, $6.24·{10}^{-4}$, and $4.68·{10}^{-4}$ kg/m s, respectively. We integrated Equations (1) and (10) forward in time using the explicit (forward Euler) scheme for the advection and mean field potential forcing terms, and the implicit (backward Euler) scheme for the viscous term. The forward time stepping of the scheme was adaptive, with the time step set to 20% of the maximal allowed by the Courant number.

Here, observe that the goal of the simulation is not the local accuracy of the solution but rather the accurate capture of the statistical regime of the dynamics (and, in particular, the transition to nonlinear chaos). This means that the numerical integration scheme should be chosen to avoid the introduction of artificial damping into the advection part of the system. Due to this reason, the simple forward Euler method appears to be better for such a specific purpose than more advanced numerical integration schemes, such as the 4th order Runge–Kutta method, since the latter tend to introduce artificial damping into numerical solutions as a result of their superior stability properties [49].

Results

We found that the flow fully developed by the elapsed time $t=0.07$ s in all simulations. In Figure 2, we show four snapshots of the speed of the flow, taken in the longitudinal symmetry plane of the pipe at $t=0.15$ s, which illustrates the transition between the laminar and turbulent flow. For $\mathit{Re}=1000$, the flow is laminar, as evidenced by the smoothness of the level curves, and symmetric relative to the axis of the pipe. For $\mathit{Re}=2000$, the symmetry of the flow is broken, and small intermittent fluctuations appear in the otherwise laminar flow. These fluctuations become larger and more numerous for $\mathit{Re}=3000$; for $\mathit{Re}=4000$, the flow is fully turbulent. Clearly, the range of values of the Reynolds number, at which the turbulent transition occurs in our model, agrees with observations [39]. The breaking of the flow symmetry during the transition to turbulence is likely associated with the onset of chaos in the dynamics, and, in numerical simulations, happens due to exponentially growing machine round-off errors. The hypothesis that turbulence is a manifestation of nonlinear chaos has also been discussed in the literature (see Letellier [40] and references therein).

In addition to the snapshots of the speed of the flow, we computed the time averages of its kinetic energy spectrum. The computation was done within the central core of the pipe of $3.6\times 3.6$ cm${}^2$ in cross-section, extending between 0 and 24 cm of the length of the pipe (shown in red in Figure 1) and thus largely containing the jet stream. The time averaging was carried out in the interval between $0.1$ and $0.2$ s of the elapsed time. For the detailed description of the energy spectrum computation, see our recent works [32,33,34].

In Figure 3, we show time averages of the computed Fourier spectra of the streamwise component of the kinetic energy for the same simulated flows, which are displayed as functions of their Fourier wavenumber $k_x$ in the longitudinal direction of the pipe on a logarithmic scale. In addition, we show the kinetic energy spectrum for the simulation with ${\mu }_0=1.825·{10}^{-5}$ kg/m s ($\mathit{Re}$∼${10}^5$), which corresponds to the viscosity of air at normal conditions. For the reference, we show two power slope lines, given via $E_0{k}_x^{-8/3}$ and $E_0{k}_x^{-5/3}$, which share the same empirically chosen scaling constant $E_0=20$ m${}^2$/s${}^2$.

At the large scale Fourier wavenumbers, the structures of the kinetic energy spectra of all computed flows are similar, while the major differences are observed at moderate and small scales. For $\mathit{Re}=1000$, 2000, and 3000, the kinetic energy spectrum at small scales generally appears to match the $k^{-8/3}$-decay slope, with the following variations. In the $\mathit{Re}=2000$ regime, the spectrum decays along the $k^{-8/3}$-slope rather monotonously, whereas in the $\mathit{Re}=1000$ regime (which is fully laminar) the spectrum also exhibits oscillations around this slope. As we hypothesized in [34], the latter could be a manifestation of the quasi-periodic dynamics at the unstable Fourier wavenumbers, with the periodicity of orbits destroyed by chaos as $\mathit{Re}$ increases to 2000. In the regime with $\mathit{Re}=3000$, an unusual growth of the energy spectrum is observed at small scales. Remarkably, the rate of decay of the kinetic energy spectrum for the turbulent regime $\mathit{Re}=4000$, as well as that of the air, approaches Kolmogorov’s $k^{-5/3}$-slope.

5. Summary

In the current work, we formally introduce viscosity into our model of turbulence via an intermolecular potential [32,33,34], to investigate the transition between laminar and turbulent air flows with varying Reynolds number. We numerically simulate the air flow at normal conditions in a straight pipe at different values of the Reynolds number and find that the transition into turbulent flow occurs when the Reynolds number increases from 2000 to 4000. This appears to be consistent with observations, experiments, and practical knowledge. Additionally, we find that, in our model, the corresponding rate of decay of the time-averaged Fourier transform of the streamwise kinetic energy at small scales changes from the $k^{-8/3}$-slope towards the $k^{-5/3}$-slope (Kolmogorov’s law), as the flow transitions from laminar to turbulent.

The results of our work seem to be encouraging. For the first time in history, we created a model of compressible gas flow in the form of fluid mechanics equations, where, first, turbulent dynamics emerge from an initially laminar flow at appropriate values of the Reynolds number naturally and without the help of artificial disturbances, and, second, the rate of decay of the Fourier spectrum of the kinetic energy of turbulent flow in our model matches that observed in nature.

At the same time, our model has its limitations as it describes the inertial gas flow; in a realistic flow, the pressure generally fluctuates, even if slightly. However, other models of fluid mechanics have their own limitations, for example, both the incompressible and compressible Euler equations are incompatible with the process of convection (the density is constant in the former and increases when the air warms up in the latter [34]). Yet, such models are widely used because they describe specific features of the flow that are needed for relevant practical applications. Similarly, our model may find its own use, perhaps as a “stepping stone”, improving our general understanding of the fluid mechanics of turbulence.