Macroscopic turbulent flow via hard sphere potential

In recent works, we proposed a hypothesis that the turbulence in gases could be produced by particles interacting via a potential, and examined the proposed mechanics of turbulence formation in a simple model of two particles for a variety of different potentials. In this work, we use the same hypothesis to develop new fluid mechanics equations which model turbulent gas flow on a macroscopic scale. The main difference between our approach and the conventional formalism is that we avoid replacing the potential interaction between particles with the Boltzmann collision integral. Due to this difference, the velocity moment closure, which we implement for the shear stress and heat flux, relies upon the high Reynolds number condition, rather than the Newton law of viscosity and the Fourier law of heat conduction. The resulting system of equations of fluid mechanics differs considerably from the standard Euler and Navier-Stokes equations. A numerical simulation of our system shows that a laminar Bernoulli jet of an argon-like hard sphere gas in a straight pipe rapidly becomes a turbulent flow. The time-averaged Fourier spectra of the kinetic energy of this flow exhibit Kolmogorov's negative five-thirds power decay rate.


Introduction
Reynolds [31] demonstrated that an initially laminar flow in a long straight pipe develops turbulent motions whenever the high Reynolds number condition is satisfied.More than half a century later, Kolmogorov [20,21,22] and Obukhov [27,28,29] observed that the time-averaged Fourier spectra of the kinetic energy of a turbulent flow possess a universal decay structure, corresponding to the negative five-thirds power of the Fourier wavenumber.To this day, the mechanism of turbulence formation in a laminar flow, as well as the origins of the power scaling of its kinetic energy spectrum, remain unknown.
In our recent work [3], we proposed a hypothesis that the turbulence in gases could be produced by their particles interacting via a potential.Later, in [4], we studied the ability of the potential interaction in a pair of particles to form turbulent structures with power decay spectra from an initially laminar shear flow.The key difference between the conventional fluid mechanics equations and our system was the closure for the velocity moment hierarchy.While the Newton law of viscosity and the Fourier law of heat conduction are typically used to truncate the velocity moment hierarchy emerging from the Boltzmann equation [8,12,13,18] (which subsequently leads to the Euler [6] and Navier-Stokes [15] equations), in our potential-driven system we used the high Reynolds number condition to obtain the closures for the shear stress and heat flux.
In [4], we examined the dynamics for the following common interaction potentials: electrostatic, gravitational, Thomas-Fermi [14,33], and Lennard-Jones [25].For each potential, we found via numerical simulations that the time-averaged kinetic energy of the flow decays as the negative five-thirds power of its Fourier wavenumber, thus corresponding to the Kolmogorov turbulence spectrum.Out of curiosity, in [4] we also simulated the dynamics with a generic large scale potential, which did not correspond to any known type of interaction, but rather "impersonated" a generic mean field forcing, which typically arises in a Vlasov-type equation [35].To our surprise, we found that, even for a large scale mean field potential, the Fourier transform of the kinetic energy of the flow had the Kolmogorov decay rate.
The present work complements [4], in the sense that here we focus entirely on the macroscopic dynamics of a single particle in a multiparticle system under the averaged forcing from all other particles.This is similar to the conventional kinetic theory approach, with the exception that we do not replace the potential interaction between particles with a Boltzmann collision integral, and instead leave it as is.Subsequently, we truncate the velocity moment hierarchy using the high Reynolds number condition, just as we did in [4].The result is the system of fluid dynamics equations for the density, velocity, temperature and pressure, which is markedly different from the conventional Euler equations.The first difference is that the momentum equation now includes the mean field forcing, which is absent from the Euler equations.The second difference is that our system preserves the pressure along the streamlines, while the Euler equations preserve the entropy.We subsequently demonstrate, via a numerical simulations that, in our system, a laminar Bernoulli jet develops into a turbulent flow which has the Kolmogorov energy spectrum.
The work is organized as follows.In Section 2, we consider a system of many particles which interact via a potential.For this system, we compute the Liouville equation, the Gibbs equilibrium state of the system, and show that a generic solution preserves the family of the Rényi metrics [30], including the Kullback-Leibler divergence [23], between itself and the equilibrium state.
In Section 3 we carry out the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) formalism [7,9,19] to obtain the transport equation for the distribution density of a single particle.We then close this transport equation by expressing the two-particle marginal density via the single-particle density in the same manner as in our recent work [2].The result is the Vlasov-type equation with a large scale mean field forcing.In the case of a short range interaction potential, we obtain a simplified expression for the mean field forcing, and, in particular, an explicitly computable formula for the mean field forcing derived from a hard sphere collision interaction.
In Section 4, we formulate a hierarchy of the transport equations for the velocity moments of the resulting Vlasov-type equation.Due to the fact that the potential forcing replaces the usual Boltzmann collision integral, a different closure must be used to truncate the moment hierarchy.Here, we use the same closure for the shear stress and heat flux, based on the high Reynolds number condition, as we introduced in our recent work [4].This closure leads to the system of transport equations for the density, momentum, and pressure (or, equivalently, the inverse temperature).
In Section 5 we show the results of a numerical simulation of the new transport equations in a realistic setting.Namely, we simulate the flow of argon in a straight pipe.The initial condition for the flow is a laminar Bernoulli jet, which is a steady state in the conventional Euler equations.However, in our model, this laminar Bernoulli jet quickly develops into a turbulent flow, which behaves visually similarly to turbulent flows observed in nature.To verify that the behavior of our model is indeed related to turbulent flows, we compute the time averages of the Fourier transform of the kinetic energy of the turbulent flow.We observe that these averages indeed decay with the rate of negative five-thirds power of the Fourier wavenumber, which is a famous Kolmogorov spectrum.Section 6 summarizes the results of this work.

Multiparticle dynamics
Here, we consider a dynamical system which consists of N identical particles, interacting via a potential φ(r).Denoting the coordinate and velocity of the i-th particle via x i and v i , respectively, we have the following system of equations of motion: Observe that adding an arbitrary constant to φ(r) leaves the equations in (2.1) unaffected.
In what follows, for convenience we assume that this additive constant is such that φ(r) → 0 as r → ∞.For the purposes of this work, we also assume that the potential φ(r) is repulsive, that is, φ(r) 1 as r → 0. We do not, however, assume that φ(r) necessarily becomes infinite at zero, to avoid technical issues associated with this singularity.
The total momentum and energy of all particles are preserved by the dynamics: (2.2) Observe that, for a given value of the momentum, it is always possible to choose the inertial reference frame in which the momentum becomes zero; thus, without much loss of generality, we will further assume that the total momentum of the system is zero.Following our work [3], we concatenate the coordinates as X = (x 1 , . . ., x N ), and velocities as V = (v 1 , . . ., v N ).In these notations, we can write In the variables X and V , the conservation of the energy can be expressed via Let F(t, X, V ) be the density of states of the dynamical system above.Then, the Liouville equation for F is given via Any suitable F 0 of the form is a steady state for (2.5).Among those, the canonical Gibbs state is where θ 0 is the kinetic temperature.As shown by Abramov [4], (2.5) preserves the family of general Rényi divergences [30]; indeed, for ψ 1 (F) and ψ 2 (F) we have and, if ψ 1 (x) = x α , ψ 2 (x) = x 1−α for some α > 0, the following Rényi divergence is preserved in time: (2.9) The Kullback-Leibler divergence [23] is a special case of the Rényi divergence for α = 1.As we remarked in [4], the conservation of the Rényi metrics plays a similar role to Boltzmann's H-theorem -namely, if the initial condition of (2.5) is close, in the sense of any Rényi metric, to (2.7), then the solution will also remain a nearby state.
2.1.Two-particle marginal density of the Gibbs state.The marginal density of the Gibbs state for the particles #1 and #2 is given via and Y N (r) is a multiple of the two-particle cavity distribution function [10]: (2.12) Y (2) The gas is dilute if the measure of the set where the integrand of (2.12) is not unity is negligible in comparison with the measure of the set over which the integration occurs.
In such a case, V 2 Y (2) N = 1, which leads to the simplified expression for the two-particle marginal density: (2.13) F (2) 3. The closure for a single particle (Vlasov equation) Here, we isolate a single particle (say, #1), and denote x = x 1 , v = v 1 .We then examine the transport of its marginal distribution f (t, x, v), given via Integrating the Liouville equation in (2.5) over all particles but the first one, in the absence of boundary effects we arrive at where 1,i is the two-particle marginal density for particles 1 and i.The equation above is a part of the BBGKY hierarchy [7,9,19], which needs to be closed, that is, F 1,i needs to be approximated via f .In order to do this, we assume that all pairs of particles are statistically identical, and thus their joint marginal distributions are equal.This leads to where F (2) is the two-particle marginal density for any pair of particles.Here, the simplest closure would, apparently, be the direct factorization of F (2) into the product of f 's, however, such closure would not be compatible with the form of the steady state marginal distribution in (2.13).Instead, we use the same type of closure as we previously did in [2].Namely, we take the exact expression for the two-particle equilibrium state of a dilute gas in (2.13), and replace, first, the single particle equilibrium marginals f G with f , and, second, the equilibrium kinetic temperature θ 0 with that of f , computed at the midpoint between x and y: . This leads to the closed transport equation for f : This is an equation of the Vlasov type [35], with the mean field forcing provided by g.Following the standard approach [2], we renormalize f so that it becomes a mass density, f → Nm f , where m is the mass of a particle.The transport equation for f remains the same, however, g becomes renormalized as follows: Next, we assume that N is large enough so that (N − 1)/N ≈ 1, which yields 3.1.Short range potential.If φ is a short range potential, the integral above can be simplified considerably.Let us change the dummy variable of integration y → x − y: Next, let us observe that where in the last identity we add 1 to the quantity under the differentiation, so that the expression in square brackets approaches zero as y → ∞.This leads, via integration by parts, to Observe that 1 − e −φ/θ ∼ φ/θ as φ → 0, yet is bounded above by 1 even if φ is arbitrarily large at zero.Thus, if the effective range of φ is short enough so that ρ, θ and their spatial derivatives can be treated as constants within it, we can ignore their dependence on the dummy variable of integration y: This leads to Denoting, for further convenience, the mean field potential φ via we obtain the following Vlasov-type transport equation for a short range interaction potential φ: Here, observe that, generally, the mean field potential φ is time-dependent, since ρ and θ are time-dependent.
3.2.Special case: the mean field potential for a hard sphere particle.If φ is the hardsphere potential (that is, it is zero outside the radius of the sphere, and arbitrarily large inside, with a rapid, yet smooth, transition), then the expression for the mean field potential in (3.13) can be simplified even further.Indeed, observe that, for a hard sphere potential, the integral over dy in (3.13) is simply the volume of the hard sphere V sp irrespectively of the value of θ(x): Substituting the above expression into (3.13), and denoting the mass density of such a spherical particle via ρ sp = m/V sp , we obtain the following explicit formula of the mean field potential for a hard sphere particle of the mass density ρ sp :

The transport equations for a turbulent flow
Just as done in the course of the conventional formalism in fluid mechanics, here we convert (3.14) into a hierarchy of the transport equations for the velocity moments, which is subsequently truncated at a suitable point.However, unlike the conventional fluid mechanics, here we will use a different closure to truncate the hierarchy, which is based on the high Reynolds number condition.
As usual, we denote the velocity average a of a function a(v) via (4.1) a (t, x) = a(v) f (t, x, v)dv.
The velocity moments are the averages of powers of v, that is v k , where v k is the k-fold outer product (which is a tensor of rank k).The first two velocity moments are the mass density and the momentum, respectively: The transport equation for any velocity moment is obtained by multiplying (3.14) by an appropriate power of v, integrating over dv, and using the integration by parts for the term with φ.In particular, for the density and momentum we have, respectively, where ∇ now denotes the x-differentiation operator (as the v-derivatives are no longer present).Observe that, due to the presence of v in the advection term of the Vlasov equation (3.14), each moment equation is coupled to a higher-order velocity moment (i.e., the equation for ρ is coupled to ρu, the equation for ρu is coupled to v 2 , and so forth), thus forming a hierarchy of the velocity moment transport equations.In order to truncate this hierarchy, a closure is needed.
To obtain a suitable closure, let us write the transport equation for the quadratic moment: where, as above, the terms with the potential forcing were integrated by parts.Next, we decompose the quadratic and cubic moments as follows: we introduce the temperature tensor T , and express where "⊗" denotes the outer product, while "T" and "TT" denote the two permutations of a symmetric 3-rank tensor.The kinetic temperature θ is the normalized trace of T , while the viscous shear stress S is the corresponding deviator: The product ρθ is known as the pressure.By construction, S is traceless.Likewise, we decompose the centered cubic moment by introducing the heat flux q, and the corresponding deviator Q: By construction, the contraction of Q along any pair of its indices is zero.In the new variables, the transport equation for the quadratic moment can be replaced by the transport equations for the pressure and shear stress: The equations above are obtained by subtracting appropriate combinations of the mass and momentum transport equations in (4.3) from the quadratic moment transport equation in (4.4) (for more details, see Grad [16], Abramov [1,2], or Abramov and Otto [5]).

4.1.
A turbulent closure for the moment transport equations.As we can see, the moment transport equations in (4.3), (4.4) and (4.8) are chain-linked together via the higherorder moment in the advection term, forming an infinite hierarchy of equations.For a practical computation, this hierarchy must be truncated at some point, so that the number of equations to be solved is finite.
For the conventional velocity moment hierarchy emerging from the Boltzmann equation, such truncation typically occurs at the equations for the shear stress S and heat flux q.The reason why it happens is that the time-irreversible collision integral of the Boltzmann equation creates linear damping in the shear stress and heat flux equations, which is inversely proportional to the viscosity.For a small viscosity, this damping becomes strong, and one can replace the shear stress and heat flux equations with their respective steady states, obtaining the Newton law of viscosity, and Fourier law of heat conduction, respectively.This results in the Navier-Stokes equations, or, if the viscosity is set to zero (thus making the shear stress and heat flux damping infinitely strong), in the Euler equations.For rarefied gas flows with a large Knudsen number, the truncation can be made at a higher order in the velocity moment hierarchy, resulting in the 13-moment Grad equations [16], or the regularized 13-moment Grad equations [32].
Here, however, observe that such collision damping is not present in the transport equation for the shear stress (4.8b), thus necessitating a different closure.Thus, we follow our recent work [4] and use the high Reynolds number condition to truncate the moment hierarchy.Let L and U denote, respectively, the characteristic length scale of the flow and its reference speed, such that the ratio U/L specifies the characteristic time scale.Additionally, T will specify the reference temperature of the gas, while V will be the reference viscosity constant.We rescale the time and space differentiation operators, as well as the velocity, temperature, the mean field potential, the shear stress, the heat flux, and its deviator, via In the rescaled variables, the mass, momentum, pressure and stress transport equations in (4.3) and (4.8) become where Ma and Re are the Mach and Reynolds numbers, respectively: Here, for convenience, our definition of the Mach number differs from the conventional one by the adiabatic constant factor.
As in [4], here we assume that Ma ∼ 1, Re → ∞.Subsequently, the terms with the shear stress S are discarded from the momentum and pressure equations (as they are divided by Re).The closure for the heat flux q follows from the stress equation, where we set the terms, multiplied by Re, to zero, and take into account the fact that the trace of Q along any pair of its indices is zero: Recall that the Newton law of viscosity and the Fourier law of heat conduction, which are used to close the conventional hierarchy of velocity moments of the Boltzmann equation, are obtained from the truncated Hilbert expansion of f in Hermite polynomials around the equilibrium [13,16,18].To be accurate, this truncation requires the distribution f to be sufficiently close to the Maxwell-Boltzmann equilibrium state.In contrast, our high Reynolds number closure above makes no assumption on the form of f in (3.14), and does not require f to be related to the equilibrium state in any particular way.Substituting (4.12b) into the pressure equation (4.10b), and reverting back to the dimensional variables by discarding the Mach number, we arrive at the following closed system of macroscopic turbulent equations for a short range interaction potential: In (4.13b), either of the two equations can be used to compute the solution; however, for a numerical simulation, the inverse temperature equation is more practical as it has the form of a conservation law.The notable difference between (4.13) and the conventional Euler equations is that the former have an additional forcing term ∇ φ in the momentum equation, and preserve the pressure ρθ along the streamlines.
In [4], we further imposed the hydrostatic balance assumption onto the respective analog of (4.13).This was necessitated by the fact that, in the context of [4], the transport equations were forced directly by the interaction potential φ, which usually peaks strongly near zero and thus can force the velocity to attain unrealistic values.In the present context, however, the momentum equation in (4.13) is forced by the mean field potential φ, which is given via (3.13), and behaves regularly throughout the domain.In addition to this, in practical scenarios, the effect of φ is expected to be weak, so one cannot realistically expect a hydrostatic balance to hold for (4.13).Thus, in the context of the present work, we abstain from further simplifications of (4.13).
We also have to note that the equations in (4.13) are "pathologically turbulent", in the sense that, by construction, they cannot model anything but a turbulent flow with a high Reynolds number.For example, an airfoil, being moved in an otherwise steady gas at a constant pressure, will be unable to create, in the context of (4.13), the aerodynamic lift (or the drag, for that matter), due to the fact that the pressure remains constant along the streamlines.The conservation of the pressure along the streamlines is the direct consequence of the high Reynolds number condition requirement, which is used to achieve the heat flux closure in (4.12b).4.2.Special case: turbulent flow at constant pressure.A notable special case of the turbulent flow equations in (4.13) is where the pressure is constant throughout the domain.This occurs for example, in the case of a simple channel flow at a constant rate (such as the Bernoulli flow), which also happens to be a steady state for the conventional Euler equations.For a constant pressure ρθ, the turbulent flow equations in (4.13) become Although the pressure is no longer a variable in the equations, it still affects the dynamics as a parameter.For example, in the case of a hard sphere mean field potential (3.16), the momentum equation is given via where p 0 is the constant pressure parameter.It is clear that the magnitude of the constant pressure parameter affects the potential forcing, and, therefore, the development of the resulting turbulent effects.

Numerical simulation
Here, we show the results of a numerical simulation of the turbulent flow of an argonlike hard sphere gas in a straight pipe of square crossection at normal atmospheric conditions, assuming that the interaction potential is that of a hard sphere (3.16).The constant sphere density parameter in (3.16) is set to ρ sp = 2716 kg/m 3 , which is estimated from the mass of an atom of argon (6.634 • 10 −26 kg) and its van der Waals diameter (∼ 3.6 • 10 −10 meters).The length of the pipe is 10 meters, while the width and the height are both set to 2.5 meters.
To carry out the numerical simulations, we use OpenFOAM [36].Since the turbulent flow equations in (4.13) comprise a system of conservation laws, we implement them with the help of an appropriately modified rhoCentralFoam solver [17], which uses the central scheme of Kurganov and Tadmor [24] for the numerical finite volume discretization, with the flux limiter due to van Leer [34].The time-stepping of the method is adaptive, based on the 20% of the maximal Courant number, which is evaluated over the whole domain at each time step.The spatial domain is discretized uniformly in all directions, comprising 400 points along the length of the pipe, and 100 points in each of the two transversal directions (4 million discretization cubes in total).
We choose the initial condition to be a steady (in the context of the conventional Euler equations), constant pressure, laminar Bernoulli jet, with the velocity profile given via Above, r is the distance to the central axis of the pipe, and the parameters are chosen to be u 0 = 10 m/s, and r 0 = 0.5 m.As we can see, the velocity of the jet equals 10 meters per second at the axis of the pipe, and smoothly decays to zero at the distance of half a meter away from the axis.Our set-up corresponds to the Reynolds number in the Re ∼ 10 6 range; indeed, the size of the domain is ∼ 1 meter, the reference velocity is ∼ 10 meters, while the kinematic viscosity of argon at normal conditions is ∼ 10 −5 m 2 /s.The constant pressure is set to p 0 = 10 5 kg/m s 2 , which corresponds to the normal atmospheric pressure at sea level.The temperature θ is calculated according to Bernoulli's law for monatomic gases: Above, R = 8.314 J/mol K is the universal gas constant, M = 3.99 • 10 −2 kg/mol is the molar mass of argon, and T 0 = 288.15K is the atmospheric temperature at normal conditions.The density of the initial flow is calculated according to .
The boundary conditions are specified as follows: • For the velocity u, the inlet and wall boundary conditions are fixed at the values which correspond to the initial flow, while the outlet boundary condition is set to zero normal gradient; • For the pressure p, the inlet and wall boundary conditions are set to zero normal gradient, while the outlet boundary condition is set to constant pressure corresponding to the initial flow; • For the inverse temperature θ −1 , the boundary condition is set to zero normal gradient everywhere.Effectively, this set-up corresponds to the constant-pressure equations in (4.14), with the density boundary condition set to zero normal gradient.
Upon conducting the numerical simulation, we found that, for the above specified conditions, the time scale of development of turbulent flow from the laminar initial condition is only a few seconds.In Figure 1 we show the two-dimensional snapshots of the magnitude of the flow velocity, measured within the plane z = 0 (that is, the plane that dissects the pipe into two symmetric halves along its axis) at elapsed times of 1, 2, 3 and 4 seconds.As we can see, at one second, the jet already has small turbulent fluctuations within it.At two seconds, the last third of the jet develops larger scale fluctuations, which eventually lead to a break-up.At three and four seconds, the jet shows three distinct stages: roughly, the first third of the jet is intact, the second third has small fluctuations, while in the last third the jet completely falls apart.Even though we have not observed this turbulent jet approach a steady state, it appears that its statistical steady state is practically reached at about three seconds of the elapsed time -that is, further simulation shows a similar principal structure of the jet, consisting of the same three stages.
In Figure 2 we show the show the two-dimensional snapshots of the density, measured within the plane z = 0 at elapsed times of 1, 2, 3 and 4 seconds.It is interesting that the most significant fluctuations of the density occur within the jet itself, and not in the turbulent flow around it (for example, compare the snapshots at t = 1, where the jet is mostly intact throughout the length of the pipe, and at t = 4, where the jet falls apart in the last third of the pipe).The fine "checkerboard" fluctuations on the spatial scale of the discretization above and below the jet appear to be an analog of the acoustic waves; indeed, combining the time-derivative of the density equation in (4.14) with the divergence of the momentum equation in (4.15) leads to (5.4) where the matrix which multiplies ρ is obviously positive definite.The "speed of sound" here is, however, quite low (in comparison with the true acoustic waves) -observe that, for the given pressure and the density of the atom of argon, we have p 0 /2ρ sp ∼ 20 m 2 /s 2 , plus u 2 contributes at most 100 m 2 /s 2 , which yields about 10-11 m/s for the speed of such a wave in the given conditions.While one can probably devise numerical methods to suppress these "acoustic waves" by filtering out the appropriate spatial or temporal scales, here we leave them as they are, since they do not seem to significantly affect the turbulent dynamics of the jet in Figure 1.
To verify that the behavior of the jet in Figure 1 indeed corresponds to the turbulent behavior observed in nature, in Figure 3 we show the time-averaged Fourier spectra of the kinetic energy of the flow.These spectra is computed as follows: first, the kinetic energy E(x) = u(x) 2 /2 is averaged over the cross-section of the pipe (thus becoming the function of the x-coordinate only).Then, the linear least-squares fit is subtracted from the result in the same manner as was done by Nastrom and Gage [26]; this is done to ensure that there is no sharp discontinuity between the energy values at the inlet and the outlet.Finally, the one-dimensional discrete Fourier transformation is applied to the result, so that the resulting Fourier transform of the kinetic energy, computed along the pipe, corresponds to the zero Fourier wavenumbers across the pipe.The subsequent the kinetic energy averages.Also, the plots are visually similar to the observations in the experimental work by Buchhave and Velte [11].

Summary
In the current work, we develop and test a fluid mechanical model for a macroscopic turbulent flow, which is triggered by a short range interaction potential.We start with a system of many identical particles, which interact via a short range potential.Applying the BBGKY formalism and a suitable closure, we obtain a Vlasov-type equation for the distribution density of a single particle with a mean field forcing.We show that the mean field forcing has an explicit and relatively simple form if the interaction potential corresponds to a hard sphere collision.
Next, for this Vlasov-type equation, we formulate the velocity moment hierarchy, similarly to how it is done in the conventional transition from kinetic theory to fluid mechanics.However, due to the fact that the mean field forcing replaces the conventional Boltzmann collision integral, the closure for the moment hierarchy differs from the one typically used in the conventional fluid mechanics; namely, we have to use the high Reynolds number condition for the closure instead of the Newton and Fourier laws.As a result, we arrive at a system of transport equations for the density, momentum and temperature, which is markedly different from the conventional Euler equations.First, the momentum equation in the new system has an additional forcing originating from the interaction potential, and, second, the new system preserves the pressure along the streamlines (while the conventional Euler equations preserve the entropy).
Finally, we numerically simulate the new equations in a realistic setting of argon flowing through a straight three-dimensional pipe.The laminar Bernoulli jet is chosen as an initial condition for the simulation, given that it is a steady state in the conventional Euler equations.Our simulation, however, shows that the laminar Bernoulli jet quickly develops into a fully turbulent flow, which visually resembles the observations in nature.We also confirm that the time averages of the kinetic energy of the simulated flow decay with the rate of negative five-third power of the Fourier wavenumber, which corresponds to the Kolmogorov spectrum.