Mean-field theory of turbulence from the variational principle and its application to the rotation of a thin fluid disk

... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A new mean-field theory of turbulence that treats the effective viscosity as a dynamical degree of freedom is presented on the basis of the stationary action principle, and is shown to reproduce some experiments for the mean flow profile of turbulence in the laboratory. Comparison with eddy viscosity models is also made. Then, with the help of the viscosity expansion method, the theory is applied to a rotating thin fluid disk for the purpose of evaluating the viscosity effect in the dynamics of a spiral galaxy or protoplanetary nebula. We find two types of physically interesting solution. In the first type, the rotation curve at long distances from the disk center is flat as a natural feature of the rotating viscous fluid with neither strong radial motion nor radial pressure gradient. The flow is gravitationally maintained only when a sufficient amount of matter other than viscous fluid is present. In the second type, the rotation is Keplerian with a centrally localized mass distribution. In both types of solution, the effective viscosity tends to act to stabilize perturbation in the region of shorter distances. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Subject Index E27, J14, J15, J18


Introduction
The structure and internal motion of spiral galaxies have been shown to be understood within the framework of fluid dynamics.This approach may be more pertinent in the early stages of galaxy evolution because the galaxy is a viscous gaseous disk constituted of gas, dusts, and asteroids which behave as molecules of viscous fluid at astronomical scales [1][2][3].As the temperature cools and the constituents aggregate into larger stellar bodies of smaller number density, the chance of collisions between them, and thereby the viscosity, is gradually lost at large scales and the "fluid" will become effectively inviscid.Such a view provides a natural basis for traditional studies of the astronomical fluid in terms of the Euler equations [4][5][6][7][8].
On the other hand, the role of viscosity in the evolution of differentially rotating gaseous bodies has been noted through studies of the mechanism that separates mass and angular momentum.The viscosity is effectively contributed from the shift of momentum by random molecular motion and the gravitational collapse of matter.Relevant major driving mechanisms may be the gravitational encounters of particles and clouds of mass, and the differential rotation of the disk [9,10].
Turbulence and viscosity are sources of complexity, and put obstacles in the way of analytical studies.On an astronomical problem such as the one mentioned above, direct numerical simulations have been intensively performed for quantitative studies [11,12], and the method is taken over to the current setting of the dark matter problem in the evolution of galaxies [13,14].
In fluid dynamics, it has been established that the Kolmogorov cascade from larger to smaller eddies takes place toward the viscous dissipation at the smallest scale.Because eddies are deviations from the mean flow, the concept of the eddy viscosity closely related to the Reynolds stress in turbulence has been put forward, and, using the conventional mean-field theory of the Reynolds stress, eddy viscosity models (EVMs) have been developed in fluid dynamics and hydraulic engineering.For details of the EVMs, see, e.g., reviews by Launder and Spalding [15], Rodi [16], Suga [17], Bredberg [18], and Launder and Sandham [19].For recent progress in understanding turbulence, see, e.g., Marusic et al. [20].
Eddies or turbulence with sizes ranging from very large [21] to protoplanetary [22] may be ubiquitous in nebulae and galaxies.Von Weizsäcker [23,24] has already applied the concept of turbulent viscosity to the evolution of astronomical rotating gas.Since the Kolmogorov-like cascade of eddies is also expected to inhere in such systems [25], we address the question of whether an EVM devised to understand the phenomena on earth is applicable to astronomical problems.However, there exists a problem to be acknowledged before this question is answered.Such being the case, let us take a look at the so-called K-ε model, a typical EVM.
The eddy viscosity is a representative quantity in dissipative eddy dynamics, and is expressed phenomenologically by ν t = C μ K 2 /ε, where K is the turbulent energy, ε the rate of energy dissipation, and C μ a constant.The K-ε model is a dynamical system of K, ε , and the mean velocity.Its equations of motion are derived on the basis of the Navier-Stokes equations by taking various moments of fluctuations.The number of equations is infinite.In order for such a mathematical system to work, the equations are truncated at the second-order moment and an assumption, the Boussinesq hypothesis (see Sect. 2.2.), is adopted that relates the second-order moment of fluctuations (Reynolds stress) to the mean velocity gradient via ν t .Arbitrariness in such a modeling is sifted by experiments in laboratory.
The method of truncation is, to some extent, arbitrary.Specifically, there exist situations in which the original Boussinesq hypothesis does not hold [16,19].The relations among higher-order moments and lower-order ones are derived sometimes by scaling arguments, which are effective but do not give unique solutions.This kind of unsatisfactory feature, called the closure problem, is inevitable in EVMs and is not yet handled in a unified way.Because of the difficulty inherent in EVMs, reproducing the observed pattern of a flow with boundaries needs somewhat ad hoc prescriptions like damping functions to be repeatedly employed [17][18][19]26].
In passing, we note that there have been endeavors to derive moment reduction formulae by eliminating the above arbitrariness in the EVMs or by seeking theories free from the closure problem, in principle.Such studies are usually implemented by taking advantage of the Fourier series expansion of the velocity field because the equations of motion assume relatively compact expressions for incompressible fluid.The direct interaction approximation proposed by Kraichnan [27] is a perturbative method in Fourier space for homogeneous turbulence that makes the special assumption of statistical independence of the interacting degrees of freedom, and has been improved or extended by many researchers.By the expansion method in a scale parameter that is supposed to conveniently separate the slow and fast modes, Yoshizawa [28] and Okamoto [29] systematically derived the moment reduction formula for the Reynolds stress, thereby improving the Boussinesq hypothesis.Their approach was extended by Yokoi and Yoshizawa [30] and Yokoi and Brandenburg [31] to incorporate the helicity, which is closely correlated with the inhomogeneous mean velocity, and to bring methodological consistency to the K-ε model.These works were aimed at obtaining formally exact expressions for the Reynolds stress for incompressible turbulent flows.For practical use, however, some assumptions on the statistical properties of fluctuating quantities together with calculation approximations are needed.A compact guide to this field is found in Ref. [32].
The structures of the astronomical system we are interested in are quite inhomogeneous.For example, the density of the visible mass in our galaxy exhibits a spiral structure and varies from ∼ 10 −16 kg m −3 in the central bulge to ∼ 10 −22 kg m −3 or less in the interstellar medium [33].Such structures are known as the grand design of spiral galaxies and seem to be common to numerous galaxies in the extragalactic domain of our universe.They must be the evolutional consequence of the instability of a rotating disk of compressible gas subjected to self-gravity.Once instability is triggered, the density wave and arms grow.The acoustic speed in our galaxy is of the order of mκ/k = O(10 km s −1 ), where m is the number of spiral arms, κ the epicycle frequency, and k the wavenumber [2,6].Since the rotational speed of matter is generally supersonic (of the order of 10 2 km s −1 ), shocks are generated and the rate of compression becomes non-negligible.Because of this feature of galaxies, the possibility of making use of existing EVMs for incompressible fluid seems dubious for the study of galactic spatiotemporal structures from the viewpoint of turbulence.If we aim at consistent understanding of the inhomogeneity in astronomical gaseous objects in terms of turbulent fluctuations, we instead need a dynamically closed turbulence model of compressible fluid that is applicable not only to static but also dynamical problems of a gas disk.For recent developments of the hydrodynamical simulation method in this area, see, e.g., the review by Sellwood [34] and the references cited therein.
EVMs that take account of compressibility have also been proposed by, e.g., Yoshizawa [35], Speziale [36], Gatski and Bonnet [37], and Germano et al. [38].The method, being effective with the large-eddy simulation [39], in which small-scale components are separated from large-scale ones by so-called grid filtering, consists of evaluating the fluctuations of such physical quantities as velocity, density, and Reynolds stress from the density-weighted ensemble averages.For correlation functions of fluctuations from filtered quantities, as in EVMs for incompressible fluid, some statistical assumptions are employed.Since the density is also treated as a dynamical variable, the method needs to model not only the second moment of fluctuations but also the third moment.Together with the need for calculation approximations, problems arise similar to the ones we encounter in EVMs for an incompressible fluid.Here, a question arises: "Is there a guiding principle for constructing a model of turbulence for compressible fluid that is free from these problems and is suited for analytical studies toward our aim?" This paper is divided into two parts.We first propose a method of the stationary action principle to overcome the above deficit in the prevalent EVMs.The stationary action principle requires the distinct entities in a system to mutually adjust their trajectories in phase space so as to keep the quantity called action stationary under infinitesimal deviations of their real paths.It is known that the Navier-Stokes equations are derived from the stationary action principle when a stochastic method for diffusive velocity [40][41][42][43] or a thermodynamic condition on dissipation for viscous fluid velocity [44,45] is introduced.In either method, the dynamical ingredient is the velocity field.
In the present problem, the entities are mean velocity and eddy viscosity.Given that they are represented by respective field functions, we do not have a priori unique knowledge of their mutual interactions.However, if there is an action for the velocity and the eddy viscosity, we will be able to derive their equations of motion in a closed form.In this paper, instead of asking how the "eddy" viscosity is generated and behaves in accordance with the fundamental dynamics of a fluid, we would PTEP 2017, 083J01 K. Takahashi like to primarily try to construct a mean-field model of velocity and viscosity on the basis of the stationary action principle.This is the precise meaning of a "dynamically closed" model.Since this mathematical model will be for the velocity and viscosity and does not involve the Reynolds stress, we do not a priori know if the "viscosity" of this model is equivalent to the eddy viscosity in EVMs constructed on the basis of the Navier-Stokes equations.Therefore, for the time being, we shall call our viscosity the effective viscosity.
Note that the concept of closure in our approach must be discriminated from the standard one in EVMs.There, closure usually implies that various model parameters in the transport equations for quantities characterizing fluctuations are determined by comparing the model with an experiment, and can be used to predict the results of different experiments.The complete resolution of the closure problem in this context does not exist.However, for a recent attempt, see, e.g., Ref. [46].
It will be desirable that this attempt is made by resorting to invariance principles that are as general as possible.For this purpose, we integrate the velocity and the effective viscosity to a single representation of the group GL(2,C) and, in terms of such a representation, write down an action that is invariant under rotation and Galilean transformation.At the same time, we require the resultant equations of motion to reduce to the Navier-Stokes equations when the effective viscosity is set equal to the molecular viscosity.We will see that the simplest version of such models, which we shall call the minimal dynamical effective viscosity model (MDEVM) for brevity, is possible and semi-quantitatively reproduces experimental results for the mean velocity in well-developed turbulent channel and pipe flows.
As the second step, having galaxies and protostar nebulae in mind, the rotation of a thin fluid disk is studied using the MDEVM.The viscosity-expansion (VE) method, which was useful in solving the Navier-Stokes equations for the problem of rotating fluid with small viscosity [47], will be employed to find the rotation curve and the density profile of the disk.
In the next section, the MDEVM is constructed and applied to turbulent flows in the laboratory.In Sect.3, using the VE method, the equations of motion in cylindrical coordinates are derived.In Sect.4, stationary solutions to the equations given in Sect. 3 are sought.In Sect.5, a brief comment on perturbations is presented.Section 6 is devoted to summary and conclusions.

MDEVM and its applications to turbulent channel and pipe flows 2.1. Construction of the model
We first present the model construction [48].The dynamical variables in the MDEVM are the complex mean velocity u, which is a vector, and the complex (pseudo-)scalar φ, which corresponds to the eddy viscosity in EVMs.Physical observables are given by the real parts of these variables.We integrate u and φ to form a complex 2 × 2 matrix as where σ i , i = x, y, z, are the Pauli matrices.The 's form GL(2,C). Its center is given by φ, which is required in order for the 's to close under products.In terms of , we construct the "pseudo-action" A for the MDEVM by ) 2 Tr , (2.2f ) The dot stands for a partial time-derivative.c V , c 3 , λ 0 , and ξ 0 are real constants.F is a matrix that includes the pressure gradient and the body force.L Ld + L f in L corresponds to the Euler equation, which is invariant under rotation and Galilean transformation.The remaining terms in L consist of Tr , ∇ , and their complex conjugates that are also invariant under the above transformations.Concerning the rotation, we can also say that L is obviously invariant under The "potential" V (ξ ) is chosen so as to be odd under ξ → −ξ in order to maintain the ν-inversion invariance of the equations which governs steady vortex solutions [47].(See also Sect.3.2.)In Eqs.(2.2), the kinetic and interaction terms in the Lagrangian density L are represented in the "minimal" way, i.e., by L φ kin , L φ V , and L (3) , although other higher-order terms may also be possible.For example, for the "potential" in Eq. (2.2e), a form like ξ 2 0 ξ 2m+1 /(2m + 1) − ξ 2m+3 /(2m + 3), which is odd in ξ and has a local minimum at ξ = ξ 0 , could be a candidate.However, it does not make much difference in the results given in what follows.
The equation of motion for is obtained by requiring A to be stationary under small variations of † .The equations of motion for real physical quantities (a concise exposition on this procedure is In the above, we defined a new variable and constants by ϕ ≡ φ/ξ 0 , ν 0 ≡ c 3 ξ 0 , λ 1 ≡ 2c V ξ 0 .In addition to Eqs. (2.3), we require the density to obey the continuity equation Notice that, by taking a limit ϕ → 1 and then letting ξ 0 → ∞, keeping other quantities finite, Eqs.(2.3) reduce to the Navier-Stokes equations for incompressible fluid.Thus, Eq. (2.3a) implies that ν 0 ϕ → ν 0 (with ξ 0 → ∞) corresponds to the molecular viscosity, ν.Comparisons of the outcomes of the model with experiments, however, are always to be done with finite ξ 0 , in which case ν 0 in Eq. (2.3a) does not necessarily coincide with ν.Later in this section, we will see that ν 0 ν is required to maintain consistency with the Reynolds stress equation.By requiring the total functional derivative of A given by Eqs.(2.2) to vanish, we obtained a set of equations that involves the Navier-Stokes equations as a special case.Note that ϕ behaves as

5/23
Downloaded from https://academic.oup.com/ptep/article-abstract/2017/8/083J01/4098172/Mean-field-theory-of-turbulence-from-the by CERN -European Organization for Nuclear Research user on 28 September 2017 a local viscosity, whose positive gradient acts as a resistance to the increase of the velocity.The consequential field-theoretical action-reaction relation is that u in Eq. (2.3b) acts on ϕ to advect it by the second term on the left-hand side, while the gradient term of ϕ 2 on the right-hand side of Eq. (2.3a) emerges as the reaction of ϕ to u.Another action-reaction relation is observed in the appearance of the viscous terms in Eqs.(2.3).There is no purely passive field, although ϕ can be treated as a passive scalar when the gradient term of ϕ 2 is sufficiently small and the reaction term in Eq. (2.3a) can be neglected.We called A a pseudo-action because its integrand does not have the canonical form of the potential energy density subtracted from the kinetic energy density.
The anticipated scale dependencies that are read from Eqs. (2.3) are as follows.The velocity shear, the second term on the right-hand side of Eq. (2.3b), will act destructively for ϕ.On the other hand, because of the positive "potential" term, i.e., the third term on the right-hand side of the same equation, the decrease of ϕ tends to bring about the growth of ϕ.The net effect will be determined by the relative sizes of λ 0 , λ 1 , and ν 0 together with the velocity shear and ϕ.Specifically, if ϕ and the velocity gradient are small (i.e., for large velocity scales at long distances), the destruction is small and ϕ will temporally grow.Contrariwise, if ϕ and the velocity gradient are large (i.e., for small velocities at short distances), the destruction is large and ϕ will decay.
In order for the system (2.3) to be a good model of fluid motions, it may have to inherit some fundamental properties of the Navier-Stokes equations.Apart from the symmetries associated with the coordinate transformations, scale invariance is the most important one.When length and velocity are measured in units of the system's characteristic length L and the characteristic velocity U (and, therefore, the characteristic time is L/U ), the explicit dependencies of the Navier-Stokes equations on L and U are eliminated and the pattern of dynamics is mainly specified by the Reynolds number Re = LU /ν 0 and Fr defined by (p,ρ, and f are not scaled) Fr is called the Froude number when f i is the acceleration of gravity and a pressure gradient is absent.We shall use the same name and symbol for any body force in this paper.
By the same scalings, Eqs.(2.3) are transformed as where Pr ≡ ν 0 /λ 0 is the Prandtl number.Equations (2.3a ) and (2.3b ) show that Re, Pr, and Fr are the three similarity numbers for the dominant component of the flow governed by the Navier-Stokes equations and the transport equation.This is realized if where the proportionality "constants" γ 0 and γ 1 may be functions of Re, Pr, and Fr.Equation (2.6) implies that ξ 0 and λ 1 are not universal constants but vary with U and L. Since we do not have any a priori knowledge of γ 0 and γ 1 , the model parameters ξ 0 and λ 1 remain undetermined.More detailed information on γ 0 and γ 1 will be obtained by comparing the outcome of Eqs.(2.3) with experimental data.
In (2.3), the field ϕ is transported by and interacts with u in an apparently analogous manner to the eddy viscosity in the EVMs.The flow equation (2.3a) is supplemented by an additional equation for ϕ.Such a system is an example of a so-called one-equation EVM.In the next subsection, therefore, we shall give a brief comment on the distinction between our model and the EVM.
It should be noticed that, since ϕ has space-time dependencies, the velocity field u in our equation system is not identical to the one appearing in the Navier-Stokes equations with a constant molecular viscosity.Specifically, the solutions to our system, even if they are stationary, generally involve the contributions from gradient-finite ϕ (i.e., the second term on the right-hand side of Eq. (2.3a)) and cannot be the laminar flows described by the Navier-Stokes equations.On the basis of this reasoning, together with the apparent similarity of our model to an EVM, we would like to interpret our u as the mean velocity in turbulent flow.This interpretation will be corroborated in what follows by applying the model to flows with boundaries.

Comparison with EVMs
Many EVMs have been proposed so far, in which, in addition to the mean flow equation, transport equations for Reynolds averaged fluctuating quantities together with their combinations are employed for the purpose of capturing the statistical property of turbulence for an incompressible fluid.In order to handle the closure problem, the higher-order moments of fluctuations in those models have to somehow be related to those of lower moments.The Boussinesq hypothesis for incompressible fluid is the one frequently referred to: where u i is the fluctuation of the velocity around the mean ūi .The eddy viscosity ν t is introduced in this way.Other relations between higher and lower moments are assumed, if need be.A large number of studies on turbulence have been reported in the literature.Unfortunately, the Boussinesq hypothesis itself is problematic in certain circumstances [16,19].Furthermore, given our goal of finding a model applicable also to astronomical compressible fluid, we cannot rely much on the accumulated knowledge of EVMs for incompressible fluid.Another limitation is inherent in our model: it involves only a scalar and a vector and not tensors, so that finding out any relation between the Reynolds stress and the mean velocity (and the effective viscosity) is beyond the scope of our model.Nevertheless, it may be instructive to compare (2.3) with an analogous transport equation in a known EVM.It is the one devised by Spalart and Allmaras [49] that employs the following Reynolds equations and local transport equation for ν t : (2.8b) We note the correspondence between the Reynolds stress term in Eq. (2.8a) and the gradient term of ϕ 2 in Eq. (2.3a).
In Eq. (2.8b), S is a measure of the rate of strain, and f w is a coordinate-dependent adjuster function for numerical tuning, which takes a finite value near the wall.The other quantities are positive constants.The second and third terms on the right-hand side are productive, while the fourth term is destructive.The role of f w is to suppress the destruction effect of the last term outside the boundary layer.Since f w is finite near the wall, Eq. (2.8b) admits a finite value for ν t near the wall, too.The Reynolds stress is assumed to be given by − The Spalart-Allmaras model has been applied to a turbulent channel flow [49].There, as the distance from the wall increases, ν t initially increases, and then the coupling term with shear starts to dominate over the destruction term; thereby, ν t decreases to a very small value at the midpoint of the flow.The method of tuning the model by introducing a more or less arbitrary function like f w was originally proposed by Van Driest [26].The model provided good fits to experimental data.
We readily notice an apparent similarity between ν t in Eq. (2.8b) and ν 0 ϕ in Eqs.(2.3).However, the role in transportation is different: In Eq. (2.3b), the destruction of ϕ takes place via the term i.e., the second-order term of the rate of strain and the rotation of the flow, while, in Eqs.(2.8), the rate of strain emerges in the first order in the production term.(The emergence of vorticity in the dynamics of fluctuation correlation is also understood from symmetry considerations [50,51], where its applications to anisotropic flows are also reported.Another application is found in Ref. [52].See also Sect.6.) Furthermore, Eq. (2.3b) involves a production term λ 1 (1 − ϕ 2 )/2, while the corresponding term c w1 f w (ν t /d) 2 in Eqs.(2.9) is destructive.In this sense, the correspondence between ν t and the Reynolds stress in the Spalart-Allmaras model will not be straightforwardly translated into a possible relation of ν 0 ϕ in our model to the Reynolds stress.

Calculations and comparisons with experiments
We are now ready to apply Eqs.(2.3) and (2.4) to well-developed turbulent flows.We first consider a channel flow of width 2d and averaged velocity ū in the x-direction under constant negative pressure gradient dp/dx and a constant body force f = (f , 0, 0), if any.Adopting the units c ≡ (λ 0 /λ 1 ) 1/2 for length and ξ 0 for velocity, the equations are written down in Cartesian coordinates as The prime denotes a derivative with respect to the dimensionless distance ŷ ≡ y/ c measured from one of the walls.Equation (2.9a) is the flow equation Eq. (2.3a) with i = x.α in Eq. (2.9a) is a constant defined using Eqs.(2.5) and (2.6) by (2.10) γ 0 and γ 1 , which were defined by Eq. (2.6), are as yet unknown functions of Re, Fr, and Pr.The flow equations for i = y and i = z give ξ 2 0 (ϕ 2 ) /2 = −p /ρ and 0 = 0, respectively.The former is used to determine the y-dependence of the pressure.
In Eq. (2.9b), the velocity shear always acts to increase ϕ .On the other hand, when ϕ is sufficiently small, the potential term in Eq. (2.9b) acts to decrease ϕ .At the same time, owing to symmetry, ϕ is negative at the midpoint as long as ϕ is smaller than unity, and ϕ continues to decrease to zero (symmetry value) until the midpoint is reached.Thus, in contrast to ν t in the Spalart-Allmaras The boundary conditions adopted are ū = ϕ = 0 at the midpoint.ϕ on the wall is very small but finite.Data are taken from Refs.[53] and [54] for channel flow with 14000 <Reynolds number (Re) <40000 and Ref. [55] for pipe flow with Re = 50000.model, ϕ will acquire a finite value at the midpoint.In this sense, our ϕ will behave like the eddy viscosity in the EVM near the wall, while its behavior far from the wall will be very different from those in the EVM.
For an axially symmetric pipe flow with diameter 2d, cylindrical coordinates are used to obtain Here, the prime denotes a derivative with respect to r = r/ c , the dimensionless distance from the central axis of the circular pipe.The body force f = (0, 0, f ) is assumed, if any.Concerning the role of ϕ, reasoning similar to that given above for the channel flow also applies to the pipe flow.
Pr and α are the essential parameters of the model.Note that Eqs.(2.9) and (2.11) are invariant under the scaling ū → sū, Pr → Pr /s 2 , α → sα.Therefore, in obtaining solutions to Eqs. (2.9) and (2.11), the boundary value of the scaled ū was fixed and the optimum values of Pr and α were sought by trial and error until global conformity with the experimental data was attained.Fitting the solutions to the experimental data was quite difficult when Pr appreciably differs from unity.The results with Pr = 1 are shown in Fig. 1. (As regards the effect of changing the parameters, decreases of Pr and/or α give rise to a tendency of ū that the increase in the transition region gets steeper and the one in the logarithmic region gentler, although this is not shown here.)Semi-quantitative agreement of the model calculations with the experiments over the viscous sublayer and the logarithmic region are obvious.Note that no additional assumption to connect both regions was employed.This result under the condition of Pr ≈ 1 may be interpreted as the equivalence of the dissipation of energy to the diffusion of ϕ, just like Kolmogorov's theory of energy cascade via formation of eddies envisages.One may call ϕ the viscosity function.
It has been known experimentally or by scaling arguments that the average velocity of turbulent channel and pipe flows has the limiting profile shown in Fig. 1  for large Re.Whether the relation (2.12) is universal or not is another problem.Even if it is universal, the functional forms of γ 0 and γ 1 together with the values of ξ 0 and λ 1 in Eq. (2.6) are still undetermined.In this sense, our model (2.3) does not have strong numerical predictability for other experimental circumstances.Nevertheless, some qualitative features of the mean flow can be understandable because of the definite model structure described by Eq. (2.3).For example, note that Eqs.(2.9a) and (2.11a) are easily integrated twice to yield (2.13) The wall is at ŷ = 0. ūc is the mean velocity at the center of the pipe.For pipe flow, the distance from the wall has been introduced by ŷ = d − r.At the very proximity of the wall, the viscosity function takes small but finite values (Fig. 1) and increases approximately linearly with ŷ.This behavior of ϕ is the cause of the distinction between the viscous sublayer and the logarithmic region in the mean velocity.
There has been a debate as to whether the mean velocity in the "logarithmic" region is logarithmic or power-like, and an experiment to resolve the problem was conducted by Zanoun and Nagib [56] (see also the references cited by them).From Eq. (2.13), the logarithmic behavior of the mean velocity is due to a power series expansion ϕ = ϕ 0 + ϕ 1 ŷδ + • • • with δ = 1.Values other than δ = 1 could also be chosen, depending on the forms of the remaining terms.In such a case, power-like behavior would result.
In Fig. 1, small discrepancies from experiments in the buffer layer are still present in the calculation for the channel flow.This unfavorable feature of the results is expected to be improved by continuing further the above trial and error method, or by adding new terms in our MDEVM, which may be a troublesome task.
As was already noted, our model does not involve tensors and cannot evaluate the Reynolds stress within the framework of the model.Here, we try to find how u and φ are related to the Reynolds stress with the help of the Reynolds stress equation.For the channel flow, the latter is given by where ν is the molecular viscosity.The meanings of the variables, constants, and symbols are the same as before.Dividing both sides by ν and subtracting (2.9a) from Eq. (2.14), we have Integrating Eq. (2.15) once, we have ξ 0 ū is the unnormalized mean velocity.For ν 0 ν, the second term in parentheses is neglected, and Eq. ( 2.18) has a form reminiscent to the Boussinesq hypothesis with ν 0 ϕ formally playing the role of the eddy viscosity.Note also that Eqs.(2.16) and (2.18) are compatible only when ū ∝ 1 − y/d near the middle point.

Mean field equations for a rotating thin fluid disk in cylindrical coordinates
In this section, we apply the MDEVM developed in the previous section to a rotating thin fluid disk in inviscid limit as a simplified model of a galaxy or protostar nebula.Since, as noted in the previous section, the model parameters in the MDEVM may generally be dependent on the system, some justification may be necessary for the tremendous leap in size from laboratory on the earth to universe before going into the detailed calculations.
We have seen, by comparing with the Reynolds stress equation, that ν 0 ϕ will correspond to the eddy viscosity ν t that connects the Reynolds stress and the mean velocity shear, provided that ν 0 ν.Thus, ν 0 will be a measure of dissipation rate due to the Kolmogorov cascade of eddy formations that is effectively viewed as the diffusion of turbulent viscosity.In this sense, Pr ∼ 1 was regarded as representing the equivalence of dissipation and diffusion.
The prominent feature of the eddy cascade is that it is mathematically expressed by Kolmogorov's scaling law in the spectrum of turbulent energy.Notice that the EVMs have effectively taken the process of eddy cascade into account by introducing the eddy viscosity ν t ≈ k 2 /ε deduced by scaling analysis in the k-ε model, or analogous scaling relations in other models.Here, k is the turbulent kinetic energy and ε the energy dissipation rate (see, e.g., Ref. [57]).Since scaling laws apply in principle to any system with arbitrary size, and since our model reproduced experiments that have been the subject of EVMs, it is natural and intriguing to apply our concept of model construction with Pr ∼ 1 to other systems with different sizes, e.g., astronomical systems, albeit with a lot of obvious differences, many of which are ignored.
There are two things we are going to take into consideration in the following.In fluid dynamical studies of astronomical disks, the effective kinematic viscosity is introduced as a measure of effectiveness of matter collisions and angular momentum transfer under differential rotation, and is treated as an increasing function of density [11,12,58].Within the MDEVM, if we pose a condition for the viscosity function ϕ → 1 at large distances from the center, then ϕ will be greater than unity in most parts of the area of the disk.In the previous section, we have seen that "steady" turbulence is realized as a balance of the productive contribution of the "potential" term proportional to λ 1 and the destructive contribution from the velocity shear, as was seen in Eqs.(2.3b) or (B.3d).In channel and pipe flows, the balance took place with λ 1 > 0 and ϕ < 1.For the astronomical disk with ϕ > 1, such a balance will be assured if λ 1 < 0 in the same equations.This is the principal change in the MDEVM dictated by the astronomical situations.Another one is the incorporation of compressibility.These will be done in the next subsection.

Equations of motion for a rotating thin disk of compressible flow
For the present purpose, it is convenient to adopt the cylindrical coordinate system (Appendix B).By "thin disk," we mean that the density and the pressure have a factor δ(z) as a z-dependence.Here, δ(z) is Dirac's delta function.Accordingly, we write ρ = σ δ(z) and p = πδ(z), with σ and π being functions of r, θ, and t for the density and pressure, respectively.All functions and equations are evaluated at z = 0. Furthermore, following Burgers [59], Sullivan [60], and Takahashi [61], u z at z = 0 is constrained as Since the density in a galaxy or nebula changes from place to place, we have to treat the system as compressible.This will be done by adding to the Lagrangian density (2.2b) a term or c 4 ∇(∇ • u) to the right-hand side of the flow equations, where c 4 is a constant, which is related to the standard notation, e.g., in Refs.[11,62], as We seek steady and axially symmetric flow solutions under the condition that the equations of motion do not explicitly depend on time or azimuthal angle.Namely, u r , u θ , ϕ, h, π/σ , and f r are non-vanishing and are functions of r only.The equations to be solved are written with the replacement Gravity appears solely in Eq. (3.3a).The equation for h is obtained by differentiating the equation for u z given in Appendix B with respect to z and by setting z = 0.Under the same condition, the continuity equation (2.4) takes the form The reason the time-derivative term is kept in Eq. (3.4) will be clarified in Sect. 4.
12/23 If we set ϕ = 1, Eqs. (3.3a)-(3.3c)reduce to the corresponding Eq. (3.1) in Ref. [47] for a constant viscosity.The above (3.3d) is the equation for ϕ, which is newly included in the present system.The correlations between the terms in these equations are generally intricate.The VE method [47,61] applicable to the case of small viscosity will give us a clearer perspective in this respect.

VE method
The VE method consists of noting first that the axially symmetric flow equations (3.3) are invariant under the following change of signs (ν 0 , λ 0 , λ 1 , c 4 , u r , u θ , h, ϕ, σ , π , f r ) → (−ν 0 , −λ 0 , (−ν 0 , −λ 0 , −λ 1 , −c 4 , −u r , u θ , −h, φ, σ , π , f r ), and that λ 0 , λ 1 ,c 4 ,u r , and h are odd functions of ν 0 , while u θ , ϕ, and π are even Refs.[47,61].Thus, the fields are subjected to the following power series expansions in ν 0 , or equivalently, in the inverse of the Reynolds number UL/ν 0 where U and L are the characteristic rotation velocity and the size of the disk, respectively: ) and so on.In Eq. 3.5, it has been implicitly assumed that the series are convergent under the condition that ν 0 is far greater than the astronomical molecular viscosity.Substituting these expressions into the equations of motion and comparing the terms with the same power of ν 0 , we have the equations for the expansion coefficients in Eq. 3.5.The equations of the expansion coefficients for the Navier-Stokes equations close up to the second order.Although this is not the case for the MDEVM, we keep the lowest-order terms only since we will eventually take the limit ν 0 → 0. Omitting hereafter the superscripts for u, π, and ϕ on the expansion coefficients (i.e., u θ → u θ , etc.), we have a set of ν 0 -independent equations: ) ) ) where the prime stands for a derivative in r, λ0,1 ≡ λ 0,1 /ν 0 , and c4 ≡ c 4 /ν 0 .λ−1 0 is the Prandtl number.By the reasoning presented at the beginning of this section, we will assume λ0 = 1 throughout the numerical calculations.Note that u θ has the dimension of velocity, while u r and h above have the dimensions of length −1 and length −2 , respectively.The last equation shows that the natural scale of length is given by ≡ λ−1/2 1 = (ν 0 /λ 1 ) 1/2 , so that the characteristic value of the radial (or vertical) velocity is given by ν 0 / = (ν 0 λ 1 ) 1/2 .On the other hand, Eq. (3.6e) implies that the natural scale for the azimuthal velocity is given by ξ 0 .
Equations Eq. (3.6e), which states that the advection of ϕ is determined by the diffusion term, the potential term, and the velocity shear term.If ϕ approaches unity from above at large distances, the signs of the potential and shear terms are opposite and may cancel each other.This canceling will take place if, for r → ∞, and other quantities decay sufficiently fast.The dynamically closed system Eqs.(3.6) bears the possibility of constant u θ at long distances.
In the following, we shall solve the equations of O(ν 1 0 ) and O(ν 2 0 ), i.e., Eqs.(3.6b)-(3.6e).Note that these equations do not explicitly depend on the pressure or the body force.Therefore, the resultant solutions reflect directly only the fluid dynamics of very small viscosity.Other dynamics like gravity enter in the balance equation of O(ν 0 0 = 1), i.e., Eq. (3.6a).

Solutions to Equations (3.6b)-(3.6e)
In this section, for convenience of the numerical calculations, we adopt the units = λ−1/2 1 for r, −1 for u r , and ξ 0 for u θ .

Viscosity function and velocity
We restrict ourselves to the case h = 0. From Eq. (3.6e), we have d dr If we reinterpret ϕ(r) as the "coordinate" at "time" rof a point mass in one-dimensional motion in the potential (ϕ 3 /3−ϕ)/2 λ0 , Eq. (4.1) expresses the rate of energy change due to the non-conservative force on the right-hand side.
Suppose that the point mass is initially at rest at a position right of the local minimum point ϕ = 1.The mass point will begin sliding down the slope of the potential toward the minimum point with the "velocity" ϕ < 0. If |u r | is sufficiently small, the first and the last terms on the right-hand side will dominate, thereby rendering the motion dissipative.Since this dissipation will fade as the time elapses, the point mass will either monotonically approach the local minimum point ϕ = 1 by losing the "energy" or oscillate near ϕ = 1.When the oscillation takes place with small or damping amplitude, its "period" will be given by 2π λ1/2 0 .ϕ will be approximated asymptotically by unity, so that (3.6a-d) are virtually equivalent to the Navier-Stokes equations with the effective viscosity ν 0 at long distances.
Once u r and h (= 0 for the present treatment of the equations) are determined, the axially symmetric density σ is also known by employing the continuity equation (3.4).We are interested in the adiabatic limit ν 0 → 0 in the VE method, so that σ in steady configuration can be of the order of ν 0 .Since Eq. (3.4) is linear in σ , it is generally allowed to have an exponential factor like exp(−ν 0 t/L 2 ) with a constant L. In this case, the time-derivative term is also of the order of ν 0 , and the continuity condition (3.4) reads The factor 2 in the last term on the left-hand side of Eq. (4.2) is due to measuring r and u r in units of and −1 , respectively.Equation (4.2) is used to calculate the density from the velocity field.Two types of physically interesting solution exist.Solutions of the first type are shown in Fig. 2. Quantities other than ϕ depend on λ1 very weakly.u θ vanishes at r = 0 and increases to some constant as r → ∞ as Such a behavior of u θ emerges owing to ϕ ≈ ϕ 0 − ϕ 1 r δ with positive ϕ 0 and ϕ 1 near r = 0 and ϕ ≈ 1 for r → ∞.Near r = 0, u r and σ behave as ≈ −ar −α with a > 0, α > 0, and ∝ r α−1 , respectively.The flow is inward.The total radial momentum flow is finite.The contribution to the kinetic energy within the radius r behaves as ∼ r 1−α if α < 1.In the example shown in Fig. 2, α ≈ 1/2 is realized.
u r vanishes as ∼ 1/r for r → ∞.In Fig. 2, /L is non-zero, so that, from Eq. (4.2), the decay of σ at long distances is Gaussian.The profile of σ in the intermediate region strongly depends on /L.A plateau is observed for ( /L) 2 0.05 in Fig. 2(d), which may be compared with the result of numerical simulation for compressible viscous fluid conducted by Laughlin and Bodenheimer [11].The fluid matter is localized within a well-defined radius, so that its self-gravity is insufficient to maintain the constant azimuthal velocity.Therefore, in order for this type of solution to be consistent with gravity theory, another kind of matter that does not interact with the fluid is required to exist, just like the dark matter hypothesis [62], or modification of the standard gravity theory will be necessary [64].
The second type of solution are those whose u θ exhibits Kepler rotation, as is shown in Fig. 3. u r decreases as ∼ 1/r for r → ∞, and therefore, as in the first solution type, the density shows the Gaussian localization around the disk center.Such a mass distribution is not gravitationally inconsistent with the Kepler motion, which is actually observed in the protoplanetary system HD 163296 [64].This second type of solution can be a primitive model of protoplanetary nebulae.The gross consistency of Figs.If we assume that masses in a protoplanetary disk on nearby orbits are collected by sweeping to eventually form a planet, as is supported by recent ALMA observation [64,65], then the mass of the planet will be proportional to r 1+ε σ with ε 0 (r 1 is due to the circumference of the orbit) and the angular momentum of the planet at radius r will roughly be given by l(r) ∝ r 2+ε σ (r)u θ (r).
(4.5) l(r) with ε = 0 (other values of ε > 0 are also possible), together with the orbital angular momenta of the planets in our solar system, is plotted in Fig. 4 for the solution drawn by the solid curves in Fig. 3.The model calculation grossly reproduces the tendency of the observed localization of angular momenta near the center of the solar system, if the scaling parameters are chosen appropriately.Discrepancies at large radii may indicate the significance of perturbation on the low-density region.Similar localization of angular momenta of fluid elements takes place for the first type of solution.

Perturbation
A uniformly rotating fluid disk generally reforms itself by generating a density wave [4].This mechanism has been regarded as responsible for the "grand design" of astronomical disks [6,8,10,12,[67][68][69][70].The next natural step toward understanding the role of this mechanism within our model is therefore to elicit the profile of small perturbations around the solutions found in the previous section.
The author previously studied the stability problem of the thin fluid disk within linear inviscid perturbations and found an exact algebraic equation obeyed by the coordinate-dependent perturbation frequency [47].Unfortunately, in the present problem, it is quite difficult to treat analytically the full set of equations of motion together with the continuity equation, even in inviscid linear approximation.In order to render the problem tractable, we assume that the radial and azimuthal perturbations are simply correlated.The perturbations are further assumed to have a common phase factor e i(ωt−nθ) , n = 0. Since the power behavior σ ∝ r , s > 0, holds for the unperturbed density, the perturbations turn out to exist when ω satisfies the relation where β with nβ > 0 is a constant of order unity, v s = dπ/dσ the sound velocity, and M (r) the total mass within a distance r from the disk center (see Appendix C for details).Equations 5.1 are compared to the result of the WKB approximation in the tight-winding limit for a spiral galaxy provided that n/r is identified with the radial gradient of the phase of the potential [6,70].By introducing the average density σ 0 by M (r) = πr 2 σ 0 , the right-hand side of Eq. (5.1a) is expressed as n(n + 2βs)v 2 s /r 2 − πn(n + 2βs)Gσ 0 /r + κ2 .The perturbation is neutral when this is positive for arbitrary r, or This is analogous to Toomre's stability condition [67,70] with κ corresponding to the epicycle frequency.What is peculiar to Eqs. (5.1) and (5.2) is that, recalling that φ in the central region significantly deviates from its value as r → ∞, and that φ < 0 there, giving κ2 larger and positive values, the contribution from the effective viscosity dominates near the central region of the disk and enhances the stability of the perturbed structure.

Summary and perspective
We first constructed a mean-field theory, MDEVM, of velocity and effective viscosity that is compatible with the stationary action principle.In MDEVM, dynamics is described as interactions among complex multiplets of velocity and effective viscosity.Incorporation of effective viscosity in the multiplet is required in order for MDEVM to be closed within a single kind of multiplet.It was shown that MDEVM can reproduce the experimentally established mean flow pattern of turbulent channel and pipe flows without employing damping functions.
We next applied MDEVM to a rotating thin fluid disk using the VE method.Two types of axially symmetric vortical flows of physical interest were found by requiring the azimuthal velocity component u θ to vanish at the center of the disk.In one type of solution, u θ approaches a nonvanishing constant value at large distances.This behavior is a natural consequence of interaction between velocity and effective viscosity, and is independent of the body force.This solution has observational support from data of galaxy rotation curves [71,72].Density perturbations generally alter the monotonic profile of the rotation curve to oscillatory ones.In this solution, since the fluid mass is quite localized, the existence of other, possibly non-viscous, matter whose density decays as 1/r is required to reconcile the solution with the Newtonian gravitational force, thereby being a possible model of galaxies with dark matter [62].Modification of the acceleration law has also been proposed [64,73].
Another type of solution exists that exhibits Kepler motion, which is also supported by the recent ALMA observation of a protoplanetary nebula [22,65,66].Such a rotation curve is compatible with the central gravity dominance because the fluid mass of this solution is quite localized around the center of the disk.
Appreciable deviations of the viscosity function from that in the small density region are observed near the center of the disk where the density is enhanced.This seems consistent with the view that the astronomical turbulent viscosity originates from cloud-cloud interactions of interstellar medium [2,9,24].
Evaluation of the Reynolds stress, one of the main concerns in EVMs, is beyond the scope of the present MDEVM, although a comparison of the model with the Reynolds stress equation provides an indication that φ plays the role of the eddy viscosity.In order to extract a relation between u, ϕ, and the Reynolds stress within our scheme of model construction, extension of the model is needed to incorporate a tensor as a representation of GL(2, C) together with interactions of the tensor with u and ϕ.A tensor R ij can be integrated into a GL(2, C) multiplet by R ij is of even parity if it is to represent the Reynolds stress.Just as for in the MDEVM, the set of R i closes under multiplication when they have the center.Thus, the above R i may be modified as If we further require R i to have a definite parity, the centre ω must be a pseudo-vector, whose candidate is the vorticity.Namely, the incorporation of the Reynolds stress into the theory necessarily means the incorporation of vorticity.In this context, it is intriguing that the significance of vorticity in the turbulent system with broken parity has been noticed [30,31].A model construction by this approach seems worthwhile.
We conclude that the MDEVM presented in this paper, with its natural extensibility, will facilitate studying the profile of mean turbulent flows with a wide variety of scales, including astrophysical phenomena.where βis a real constant of order unity (2β = 1 if the epicyclic orbit is circular at the rest frame of the fluid element in unperturbed motion).We choose the direction of the epicyclic motion of perturbation in such a way that nβ > 0. By this assumption, we consider the perturbations of u θ , φ, and σ by Eqs.(B.3b), (B.3d), and (2.The last term in the square root is the contribution from the effective viscosity.For solutions of the first type, u θ is ignored at large distances as compared with .Since σ ∝ r −s , the relation Eq. (5.1) is obtained.

Fig. 1 .
Fig. 1.Left panel: Mean velocity ū vs. distance y + from the wall.ū and y + are appropriately normalized.Right panel: Viscosity function ϕ vs. y/d with d being half the channel width or the pipe radius.Solutions are found in Cartesian coordinates (red solid curve) with α = 0.0124 and in cylindrical coordinates (blue dotted curve) with α = 0.008 for a common Prandtl number Pr = 1.The boundary conditions adopted are ū = ϕ = 0 at the midpoint.ϕ on the wall is very small but finite.Data are taken from Refs.[53] and[54] for channel flow with 14000 <Reynolds number (Re) <40000 and Ref.[55] for pipe flow with Re = 50000.
2(d) and3(d)  for the density profiles with the result of numerical simulation in Ref.[11] is because the time dependence of the density has been taken into consideration in the continuity equation (see Eqs. (3.4) and (4.2)).