A Constitutive Equation of Turbulence

Even though applications of direct numerical simulations are on the rise, today the most usual method to solve turbulence problems is still to apply a closure scheme of a defined order. It is not the case that a rising order of a turbulence model is always related to a quality improvement. Even more, a conceptual advantage of applying a lowest order turbulence model is that it represents the analogous method to the procedure of introducing a constitutive equation which has brought success to many other areas of physics. First order turbulence models were developed in the 1920s and today seem to be outdated by newer and more sophisticated mathematical-physical closure schemes. However, with the new knowledge of fractal geometry and fractional dynamics, it is worthwhile to step back and reinvestigate these lowest order models. As a result of this and simultaneously introducing generalizations by multiscale analysis, the first order, nonlinear, nonlocal, and fractional Difference-Quotient Turbulence Model (DQTM) was developed. In this partial review article of work performed by the authors, by theoretical considerations and its applications to turbulent flow problems, evidence is given that the DQTM is the missing (apparent) constitutive equation of turbulent shear flows.


Introduction
Today there is still no consensus on the correct analytical solutions of elementary turbulent flow problems as the turbulent flows in the wake behind a cylinder, the axisymmetric turbulent jets, plane turbulent Couette and Poiseuille flows, the turbulent flows parallel to a wall, etc. The absence of finally established correct solutions finds its reason at the basis of the mathematical-physical problem setting; it is the result of a lack of knowledge of the correct constitutive equation of turbulence. In this article, the authors demonstrate that a first order turbulence model, the Difference-Quotient Turbulence Model (DQTM) [1], fulfils all the requirements of an apparent constitutive equation of turbulence. It is the analogous law of turbulent to constitutive equations of other areas of physics, e.g., micro-mechanical materials, nematic liquid crystals, non-Newtonian fluids, etc. Furthermore, it is the analogue law of turbulent flows to the constitutive equation of laminar flows which is Newton's law of viscosity [2]. However, this law is linear and local, whereas the new proposed law for turbulent flows is nonlinear and nonlocal and its Reynolds shear stress mathematically presents itself as a fractional derivative [3]. Newton's law of viscosity can be derived by a microscopic theory, the kinetic theory of gases [4] which as basic constituents considers flying single-size molecules with a constant number of entities with infinite lifetimes. The analogous theory of turbulence generalizes this concept by a cascade of multi-size fractal eddies with birth rates given by Lévy statistics [5] and lifetimes determined by the fractal β-model [6], bringing also into play the aspects of self-similarity [7] and intermittency [8] which are experimentally well observable phenomena of finite Reynolds number turbulent flows. The analogy

The General Form of Constitutive Equations
Usually, a physical problem is described by a set of equations or differential equations that have been, e.g., derived by integral balance equations. In ideal cases these problems, with initial and boundary values, may be uniquely solved. Such examples are electromagnetic problems, described by the Maxwell equations for electromagnetic fields in vacuum, or the Euler equations, describing inviscid fluid flows [10]. However, when electromagnetic fields enter material domains or the fluid is viscous, the response of the material on the field's forcing must additionally be known in advance. In a physical problem this important additional information is contained in a constitutive equation that completes the problem setting. In mechanics and fluid dynamics these equations are usually stress-strain relations. The most general description of a constitutive equation is nonlocal and shows memory effect and is also called to be dispersive, viz., where → x is the spatial variable (coordinate) and → x the spatial integration variable, mostly distant from → x , µ ijkl is containing information on the material under consideration, χ kl denotes the applied physical field (external forcing of the system), Ω is the normalization constant (for its importance, e.g., in turbulent flows, see below) and summations over repeated indices k and l are assumed. Determinism only allows a time connection of the physical variables from the past to the present (−∞, t], whereas nonlocality obeys no such restriction and, in the most general case, shows long-range interactions between spatial points in the entire domain of physical relevance. Such complex nonlocal descriptions are found in micro-mechanical materials, elastic multi-phase composites, elasto-viscoplastic samples, nematic liquid crystals, non-Newtonian fluids, etc. Moreover, the physical problem may be non-isotropic and then the mathematical formulation of the constitutive equation is based on tensor calculus. If the dependence on the applied field of the physical system is weak, a Taylor approximation leads to a linear and local gradient law which in this special case satisfies the demand to possess an accurate constitutive equation. The most famous examples of this type are Hooke's law of elastic materials, in which the stress is proportional to the material's dilatation, and Newton's law of viscosity or friction, respectively, in which the fluid's shear stress is proportional to the shear velocity or the local velocity field gradient. This law is only valid for the description of laminar but not transitional or turbulent shear flows. In some cases, these laws have been derived from microscopic theories that yield physical relations for the material dependent quantities µ ijkl . In other cases, from first physical principles such results could not be achieved, and relations had to be extracted by experimental observations and mathematical curve fitting procedures which delivered approximate equations that are called phenomenological or empirical relations.

Newton's Law of Viscosity
In 1687 in his Philosophiae Naturalis Principia, Isaak Newton [11] published the constitutive equation for steady two-dimensional laminar shear flows that is called Newton's law of viscosity. To introduce it, we assume a flow with velocity u 1 in the x 1 -direction that is sheared in the perpendicular x 2 -direction and shows a (linear) velocity profile u 1 (x 2 ). Furthermore, the flow does not depend on the x 3 -direction (u 3 ≡ 0, ∂/∂x 3 ≡ 0) (see Figure 1). Then, Newton's law of viscosity states that the shear stress between fluid layers at coordinates x 2 and x 2 + dx 2 is proportional to the spatial velocity derivative du 1 /dx 2 at location x 2 . For a given pressure p and absolute temperature T, the ratio of the shear stress τ 21 to shear rate du 1 /dx 2 is a constant defining the dynamic viscosity µ of the fluid (see Equation (2) below). By substituting the relation for the generalized dynamic viscosity µ 2112 (x 2 , x 2 ) = µδ(x 2 − x 2 ), where δ denotes the Dirac distribution, and the forcing field gradient χ kl = ∂u k /∂x l into Equation (1), furthermore, applying the two simplifying conditions for shear flows: |u 1 | >> |u 2 |, |∂/∂x 1 | << |∂/∂x 2 | and setting Ω = 1, the complex integral relation reduces to Newton's single-term law of viscosity, Note that in Equation (1) the Dirac distribution constricts the large-scale spatial dependence to a single point and makes the law local. Furthermore, for constant the relation is proportional to and, therefore, this law is also linear: ( + ) = ( ) + ( ).

Some Results of the Kinetic Theory of Gases
In a microscopic theory this important law can be derived by first physical principles of flying atoms or molecules, respectively obeying Brownian motion, that exchange momentum from the high-fluid-velocity regions to such of low velocity (see Figure 1). This theory is well established and is called the kinetic theory of gases. Because of space limits the derivation cannot be presented here (see, e.g., Reference [4]). However, one of its main results is the formula for the dynamic viscosity, being = 1/3 which contains the product of a characteristic length λ * (* following conventions in the physical literature, we leave away average and root mean square (rms) signs), being the mean distance between two collisions of molecules, also called mean free path length, and the characteristic (rms) flight velocity of the molecules *. The reason that the theory is valid in its linear form is that a single molecule's diameter d is several orders of magnitude smaller than the mean free path length λ, that, on the other hand, is usually smaller than the size of the gas containing vessel L: d << λ << L (see Figure 2a). Note that in Equation (1) the Dirac distribution constricts the large-scale spatial dependence to a single point x 2 and makes the law local. Furthermore, for constant µ the relation is proportional to du 1 and, therefore, this law is also linear: τ 21 (αu 1 + u 2 ) = ατ 21 (u 1 ) + τ 21 (u 2 ).

Some Results of the Kinetic Theory of Gases
In a microscopic theory this important law can be derived by first physical principles of flying atoms or molecules, respectively obeying Brownian motion, that exchange momentum from the high-fluid-velocity regions to such of low velocity (see Figure 1). This theory is well established and is called the kinetic theory of gases. Because of space limits the derivation cannot be presented here (see, e.g., Reference [4]). However, one of its main results is the formula for the dynamic viscosity, being µ = 1/3ρλu mol which contains the product of a characteristic length λ * (* following conventions in the physical literature, we leave away average and root mean square (rms) signs), being the mean distance between two collisions of molecules, also called mean free path length, and the characteristic (rms) flight velocity of the molecules u mol *. The reason that the theory is valid in its linear form is that a single molecule's diameter d is several orders of magnitude smaller than the mean free path length λ, that, on the other hand, is usually smaller than the size of the gas containing vessel L: d << λ << L (see Figure 2a).

The Closure Problem
In fluid mechanics there are strong analogies between the constitutive equations of laminar and that of turbulent flows. However, there are also essential differences in these equations and their microscopic theories that must be explained in detail. We start with the momentum equation of fluid dynamics, also called Navier-Stokes Equation (NSE) [12]. After a proposal of Reynolds, the main three variables in this equation, namely the three velocity components, are split into two parts each: = + ′, ∈ 1,2,3 , where denotes an ensemble average (or in a quasi-steady flow the time average) and ′ is its fluctuation quantity for which ′ = 0. Then, these two-fold velocities are substituted into the NSE and the average operator is applied to the entire equation. Formally, this produces the same equation as the NSE( ), but now it depends on the averaged velocities, NSE( ) (why it is called Reynold's Averaged Navier Stokes equation (RANS) [13]). As a result of this procedure, it also contains an additional term. Because of the quadratic products / in the inertia terms of the NSE, with help of the mass conservation equation (continuity equation), the just described recipe produces nine additional second order correlations ′ ′ . These are regrouped to define a symmetric 3 × 3 tensor, being similar as the NSE's dissipation tensor . Therefore, in analogy to the dissipation stress tensor, the additional turbulent stress tensor, also called Reynolds stress tensor, = − ′ ′ , is obtained. For this tensor, the general constitutive Equation (1) with an averaged driving term = / ′ leads to the nonlocal and fractional derivative closure [14,15] where, with the relation = / , the generalized turbulent kinematic viscosity , in turbulence research also called eddy viscosity (index m = momentum), was introduced. Note the formal equivalence of Equation (3) with the general definition of a fractional derivative, where the kernel or weighting function/distribution (in our case )

The Closure Problem
In fluid mechanics there are strong analogies between the constitutive equations of laminar and that of turbulent flows. However, there are also essential differences in these equations and their microscopic theories that must be explained in detail. We start with the momentum equation of fluid dynamics, also called Navier-Stokes Equation (NSE) [12]. After a proposal of Reynolds, the main three variables in this equation, namely the three velocity components, are split into two parts each: u i = u i + u i , i ∈ {1, 2, 3}, where u i denotes an ensemble average (or in a quasi-steady flow the time average) and u i is its fluctuation quantity for which u i = 0. Then, these two-fold velocities are substituted into the NSE and the average operator is applied to the entire equation. Formally, this produces the same equation as the NSE(u i ), but now it depends on the averaged velocities, NSE(u i ) (why it is called Reynold's Averaged Navier Stokes equation (RANS) [13]). As a result of this procedure, it also contains an additional term. Because of the quadratic products u j ∂u i /∂x j in the inertia terms of the NSE, with help of the mass conservation equation (continuity equation), the just described recipe produces nine additional second order correlations u i u j . These are regrouped to define a symmetric 3 × 3 tensor, being similar as the NSE's dissipation tensor τ ij . Therefore, in analogy to the dissipation stress tensor, the additional turbulent stress tensor, also called Reynolds stress tensor, τ Tij = −ρ u i u j , is obtained. For this tensor, the general constitutive Equation (1) with an averaged driving term χ kl = ∂u k /∂x l leads to the nonlocal and fractional derivative closure [14,15] where, with the relation ν = µ/ρ, the generalized turbulent kinematic viscosity ν T , in turbulence research also called eddy viscosity ε m (index m = momentum), was introduced. Note the formal equivalence of Equation (3) with the general definition of a fractional derivative, where the kernel or weighting function/distribution (in our case ν Tijkl ) characterizes the type of the derivative; if it is, for example, a Riemann-Liouville function, the resulting fractional derivative is of Riemann-Liouville type, etc.

Kraichnan's Direct Interaction Approximation (DIA)
It is highly satisfying that Kraichnan * (* Around 1950 Robert Kraichnan was one of the last assistants of Albert Einstein at Princeton University. He performed research on quantum field theories including the transfer of some of its ideas to the field of turbulence.) in his Direct Interaction Approximation (DIA)-by completely different means-also derived Equation (3) [16,17]. The output we got from the application of the Reynolds method is the new occurrence of six additional second order correlations u i u j for which no equations exist. Note that the three non-trivial symmetries, u i u j = u j u i , reduce the number of independent correlations from nine to six. There are two applied solutions to this problem. The first method is today most frequently applied in turbulence modeling: by manipulating the NSEs it is possible to introduce a whole cascade of higher order correlations and equations relating them up to the order infinity. Over the years a multitude of different methods for handling high order systems of equations have been developed. However, all these models are irrelevant to our purposes; we believe that they are not wrong, but possibly superfluous, why we here do not discuss them any further. −The second method was introduced by early turbulence researchers: in their lowest order approach, the second order correlations u i u j are directly described by functions and/or operators of the first order mean values u i and u j . If these (or in Equation (3) all the generalized kinematic viscosities ν Tijkl ) are known, they are substituted into the RANS equations and then the system of equations contains only mean first order variables u i or/and u j . By this the system of equations becomes closed and solvable. This is the reason why in turbulence research the equation u j u i = F u i , u j is called a first order closure scheme. Now, it is recognized that the DQTM-which is a first order closure equation of turbulence and a zero-equation turbulence model * (* The name 'zero-equation turbulence model' stems from the absence of transport equations of the turbulent kinetic energy k, the dissipation rate ε or the vorticity ω, as they occur in high order turbulence models, for example, in the k-ε model or in the k-ω model.) (see Equation (5) below)-is analogous to constitutive equations in other areas of physics and, therefore, is the constitutive equation of turbulence. A small variation in the analogy is that the turbulent stress tensor is not dependent on the dynamic viscosity, determined by the molecules of the fluid; it depends on the analogous quantities which in turbulence are virtual particles called eddies. Therefore, this dynamic viscosity is called eddy viscosity, ε m , and it is an apparent viscosity. The kernel in Kraichnan's nonlocal stress tensor must be generalized by ν Tijkl where, as we will experience later, the dependence on the spatial variable with an apostrophe, → x , in this term, but not in the derivatives ∂u k /∂x l → x , t , may be removed. This replacement can be physically explained and motivated. In laminar flows the eddy viscosity is zero: It is an increasing agitation (e.g., stirring or splashing) in the fluid that makes the eddy viscosity increase. So, the eddy viscosity is not a (constant) physical property, it is the result of the dynamics of the fluid. Therefore, it depends on the difference of the Reynolds numbers Re and Re c : Re − Re c . With Re = (ud/υ), where the characteristic length d and the fluid's molecular viscosity υ are constant, in the difference of the Reynolds number Re − Re c only the velocity difference u − u c , in our case, is assumed to be a variable. In the eddy viscosity of a Reynolds stress the velocity must be Galilean invariant and, therefore, it only appears in velocity derivatives or differences, in which the constant critical velocity → u c drops out. Therefore, the dependence ν Tijkl → x , → u holds. Finally, we hope that these arguments sufficiently convince also critical readers of the applicability of the above introduced replacement ν Tijkl u which was already applied in ancient times of turbulence research by numerous scientists, as, for example, Prandtl, Görtler, von Kàrmàn, Taylor, Wieghardt, etc. The dependence of the turbulent kinematic viscosity on the velocity makes the turbulent shear stress τ Tijkl nonlinear (this is seen by a substitution of ν Tijkl → x , → u into Equation (3)), but not nonlocal. It is the integral representation of the Reynolds stresses in Equations (1) and (3) that is responsible for their nonlocal properties.

The Boussinesq Closure Approximation
Joseph Boussinesq was the first who modeled turbulent momentum transfer in analogy to molecular momentum transfer. This method led to a turbulent shear stress that we now derive by making use of our general approach (3). By applying an idea of Prandtl . Then this viscosity term, the averaged forcing field gradient χ kl = ∂u k /∂x l and the same simplifying conditions as in the derivation of Newton's viscosity law (see above), but now for the averaged velocities, are substituted into Equation (3) which, in analogy to Equation (2), yields: With a constant eddy viscosity ε m this is Boussinesq's closure hypothesis. It is a linear and local "gradient" constitutive equation of turbulence which he published in 1877 [18]. The nonlinear version (4) also rests local and is still a main core part of present high order turbulence modeling work (see first method above). Therefore, one is not astonished that with this local approach so many problems in the past had occurred. To derive the correct constitutive equation of turbulence all conventional mathematical methods of physics * (* We exclude fractional calculus from the conventional methods of mathematical physics applied in the past. However, today in numerous fields of physics this calculus is becoming increasingly important and wider recognized. It slowly becomes exposed and evident that fractional calculus is also the right mathematical tool to describe Reynolds stresses of turbulent flows.) had failed which is one reason that Richard Feynman stated that turbulence is one of the most important remaining unsolved problems of classical physics [19], a judgment that is still true today!

Prandtl's Mixing-Length Turbulence Model
The microscopic theory leading to closure (4) goes back to Ludwig Prandtl, who in 1925 published his mixing-length turbulence model [20]. As a result of the similarity of Equations (2) and (4), he applied the ideas of the kinetic theory of gases to derive a kinetic theory of turbulence. We do not present the development of the model by Prandtl, but will demonstrate an analogy of the two approaches in a very brief manner. Prandtl's derivation led him to the following result for the eddy viscosity ε m (x 2 , u 1 ) = l 2 |du 1 /dx 2 |. Here l denotes the mixing length of the turbulent momentum transfer by the eddies. Rewritten this equation yields: ε m (x 2 , u 1 ) = l |du 1 /dx 2 |l = l u eddy (x 2 ). This is the analogous result to ν = 1/3λu mol of the kinetic theory of gases (the constant 1/3 may be thought to be absorbed by l). The mixing length of the eddies l is analogous to the mean free path length of the molecules λ and as the molecules have their characteristic velocity u mol also the eddies show such a velocity, u eddy . From the second last equation, we see that, if the eddy viscosity ε m and the mixing length l are constant, u eddy also has a single value and is also constant. In research of turbulence, it is a well-known fact that there is a strong monotonic power law dependence between the eddy diameter and the eddy velocity. Therefore, a single eddy velocity corresponds to a single eddy size. As a result, one has a strong analogy between this turbulence model and the kinetic theory of gases, presented in Figure 2a. In this figure, one may replace the molecules by small single-size eddies. In this simplest approach, a single class of eddies of equal size transfers all the turbulent momentum which is surely not a realistic image of real facts. In a slightly generalized turbulence model, the velocity u eddy (x 2 ) depends on the spanwise coordinate x 2 of the flow. As a result of the strong monotonic relation the eddy size is now also x 2 -dependent. Figure 2b shows a picture of the eddies of this generalized turbulence model, where they are of multiple sizes and show a stratification. This is a good step toward a more successful turbulence model. However, studying drawings by Leonardo Da Vinci of turbulent flows [13] shows that he already had recognized that the situation in these flows is even more complex. As shown in Figure 2c, he had drawn pictures of interpenetrating eddies of multiple sizes, that act simultaneously at a given point. In the past the first author had taken this picture as a guideline to develop a 'microscopic theory' of turbulence leading to the DQTM (see [13,21,22]). The main idea was to realize a modeling of turbulent momentum transfer based on multi-size, self-similar and fractal eddies [23].

The Difference-Quotient Turbulence Model
To discover the constitutive equation of turbulence, one may make a hypothesis and then prove it by a microscopic theory. In the development of Newton's law of viscosity, a Dirac distribution occurs which leads to the (most) local description. Thinking innovatively, we may now ask the question: "In the opposite case, which kernel leads to the least possible locality?" In Kraichnan's nonlocal model ν Tijkl x , t is usually a decaying function of the distance between two points in space, x . Fact is that the less the decrease is, the more nonlocal the problem will be. Therefore, a kernel, described by a Heaviside distribution θ, that in a fluid carrier domain is 1 (constant) and outside 0 (zero), shows the least decay and leads to the strongest nonlocal approach and this shall be further investigated. Now, we assume a generally large one-dimensional domain x 2min ≤ x 2 ≤ x 2max (see Figure 1). The position where the averaged velocity u 1 (x 2 ) is at its minimum is x 2min , u 1 (x 2min ) =: u 1min , and where it is at its maximum it is x 2max , u 1 (x 2max ) =: u 1max . With an eddy viscosity dependent on the down-slope diffusivity regime [x 2min , x 2 ) and a driving force dependent on the up-slope driving region [x 2 , x 2max ], we postulate a generalized eddy viscosity to be: . Furthermore, we keep to the driving field χ kl = ∂u k /∂x l and assume simplifying conditions similar as those in the derivation of Equation (4). Now, the integration interval can be reduced from the entire flow domain [x 2min , x 2max ] to [x 2 , x 2max ], because in the interval [x 2min , x 2 ) the Heaviside integrand is zero. On the other hand, in the reduced region [x 2 , x 2max ] it is equal to unity. Therefore, in the restricted integral of Equation (3) the spatial integration directly acts on the derivative of the average velocity leading to du 1 /dx 2 dx 2 = du 1 . By this the integration simply produces the velocity difference of the mean velocities evaluated at the two bounds, x 2 and x 2max , of the integral: u 1max − u 1 (x 2 ). In the denominator with a similar integral (with the same Heaviside carrier domain, but without a velocity derivative), the normalization condition leads to: Ω = x 2max − x 2 . Substituting these results into Equation (3), with analogies to the derivation of Equation (4), the constitutive equation of turbulence occurs with σ being a constant. The modified turbulent kinematic viscosity, in analogy to the kinetic theory of gases, was set to be υ T2112 (x 2 , u 1 ) = σλ T (u 1 − u 1min ). Also in this expression, we can recognize the product of a typical length (scale) λ T = x 2 − x 2min , that now can be a large-scale nonlocal quantity, because it varies between 0 and x 2max − x 2min which may span the entire flow domain, and the averaged velocity difference u 1max − u 1 (x 2 ) which at x 2min takes its maximum, u 1max − u 1min . Note that, just as gradients, space and velocity differences are Galilean invariant. It is remarkable that the nonlocal Difference-Quotient Turbulence Model [1,24]-which was discovered by completely different ideas-coincides with the constitutive equation of turbulence (5). This equation has a geometrical interpretation [13,25]. The Reynolds stress is proportional to the shear velocity: τ T2112 = ρε m dγ 21 /dt (see Figure 1). Therefore, with and by combining two of the above equations with ε m = σ(x 2 − x 2min )(u 1 − u 1min ) Equation (5) directly follows. Note that instead of the limes ∆x 2 → 0 the limes toward the value ∆x 2 → x 2max − x 2 applies which is a quantity that tends toward its maximum to a large-scale spatial difference. Note that in the difference quotient the velocity u 1 (x 2 ) may be nonlinear in x 2 . Nonlinearity of the Reynolds stress in u 1 appears only in the combination of the generalized eddy viscosity, ε m [x 2 , u 1 (x 2 ) ], with the difference quotient.

A 'Microscopic Theory' of Turbulence
It is tempting to derive the constitutive law of turbulence (5) also by a 'microscopic theory'. Because of space reasons, we can only sketch this theory by discussing the main ideas. A reader interested in the probabilistic/statistical details is guided to references [13,21]. In a turbulent flow, the higher its forcing and the lower the viscosity of the fluid are (high Reynolds number), the smaller will be the size of the smallest eddies. Kolmogorov's dissipation length defines the low-wave-length cut-off value of the eddy probability distribution [26]. Eddies with diameters smaller than this length dissipate their turbulent kinetic energy immediately into heat. Therefore, the smallest eddies are not microscopic like atoms and molecules, but they still may be exceedingly small. On the other hand, the largest eddies may have a size comparable to the overall flow domain (see Figure 2c) which is macroscopic. Therefore, we classify this theory to be more a meso-macroscopic model *(* To keep the analogy also of the name to the molecular case, we denote this theory 'microscopic theory,' where the apostrophes shall remind us of the highly stressed use of this terminology.) than a microscopic one as, for example, the analogous kinetic theory of gases. To build this theory, at first for a quasi-steady turbulent flow the steady state of the spatial eddy distribution is required which we call the space occupation number of eddies. Other than for the molecules, for which we assume a birth rate zero (a certain number is already present and rests constant) and lifetimes which are infinite, in a turbulent flow field eddies are frequently born and die again. Advantageous is that in the field of mathematics and turbulence a theory of the birth rates of eddies, namely Lévy statistics [5], and a theory of their lifetimes, the fractal β-model [6], have been developed, so that they only must be combined. It is appropriate to identify the lifetimes of eddies with their turnover times. Now, we split up the eddies into discrete classes n with eddies of multiple sizes d n, n ∈ {0, 1, 2 . . . , N}, where d 0 denotes the diameter of the largest and d N of the smallest eddies of a cascade. The combination of these two theories leads to the occupation number of eddies of class n, o n , and O is the total sum of all eddies. In this modeling the main relations are self-similar power laws which in their exponents also contain the Euclidean dimension D and the fractal or Hausdorff-Besicovitch dimension d, respectively. Then, the momentum transfer of eddy class n to the neighboring class n+1 is determined, leading to a result of impressive simplicity: p n = (a/b) n p 0 [21]. The momentum transfer depends only on the quantity a which is the mean number of Lévy flights of length b n that occurs before a rarer flight of length b n+1 is observed, the dimensionless flight length of the smallest flights, b > 1, the number of the class of Lévy flights (see Figure 3a), n, and the turbulent momentum transfer of the largest eddies, p 0 . To obtain the total momentum transfer at a certain space location this simple scaling relation adds up for a part of the n s ∈ {0, . . . , N} as the development of model (5) [22] requires. The summations of the self-similar expressions lead to results that contain amplification factors times the momentum transfers only of the eddies being members of the largest classes which at maximum are p 0 . These amplification factors mathematically represent the action of all medium and smaller scale eddies (see Figure 3b) with the result that in the constitutive Equation (5) only two large scales are observable. However, this does not mean that only large eddies are transporting momentum; all eddies of smaller sizes are also actively and simultaneously at work by transferring momentum in the x 2 -direction. Studying this equation, we see that the largest eddies' size depends on the coordinate x 2 which is very meaningful as, for example, the active eddies at a location x 2 above a wall, which is positioned at x 2 = x 2min (by the no-slip boundary condition the mean velocity at the wall is zero: u 1min = 0), can only have the maximum size being identical to the distance from the wall, x 2 − x 2min . smaller scale eddies (See Figure 3b) with the result that in the constitutive Equation (5) only two large scales are observable. However, this does not mean that only large eddies are transporting momentum; all eddies of smaller sizes are also actively and simultaneously at work by transferring momentum in the -direction. Studying this equation, we see that the largest eddies' size depends on the coordinate which is very meaningful as, for example, the active eddies at a location above a wall, which is positioned at = (by the no-slip boundary condition the mean velocity at the wall is zero: = 0), can only have the maximum size being identical to the distance from the wall, − .

Introduction
With the new constitutive equation of turbulence (5), the following elementary turbulent shear flows have analytically been solved (where each example listed is accompanied by the corresponding reference where the flow problem was solved): the turbulent flow in the wake behind a cylinder [24], the axisymmetric turbulent jet into a quiescent surrounding [24,27] and into a high-velocity constant co-flow [24], plane turbulent Couette flow [28] and plane turbulent Poiseuille flow [29] and finally the turbulent flow parallel to a wall [21]. To convince readers of the power of the new constitutive equation, in Sections 4.7.2-4.7.4 and 4.8 we present three elementary flow problems with their essentials and solutions.

Turbulent Wake Flow
The turbulent wake flow behind a cylinder is in a certain sense special, because it is possible to derive its Reynolds shear stress by rigorous physical methods. To calculate the drag on a cylinder in a parallel flow in the -direction, the mass (continuity) and the momentum conservation equation (NSE) are applied in their integral forms. This leads to the following rigorous physical relation [24,30] where is the density of the incompressible fluid, p is a pole distance upstream of the cylinder's center and the constant upstream velocity of the parallel horizontal main flow in the -direction. It is highly satisfying that this result reveals the same large-scale

Introduction
With the new constitutive equation of turbulence (5), the following elementary turbulent shear flows have analytically been solved (where each example listed is accompanied by the corresponding reference where the flow problem was solved): the turbulent flow in the wake behind a cylinder [24], the axisymmetric turbulent jet into a quiescent surrounding [24,27] and into a high-velocity constant co-flow [24], plane turbulent Couette flow [28] and plane turbulent Poiseuille flow [29] and finally the turbulent flow parallel to a wall [21]. To convince readers of the power of the new constitutive equation, in Sections 4.7.2-4.7.4 and 4.8 we present three elementary flow problems with their essentials and solutions.

Turbulent Wake Flow
The turbulent wake flow behind a cylinder is in a certain sense special, because it is possible to derive its Reynolds shear stress by rigorous physical methods. To calculate the drag on a cylinder in a parallel flow in the x 1 -direction, the mass (continuity) and the momentum conservation equation (NSE) are applied in their integral forms. This leads to the following rigorous physical relation [24,30] where is the density of the incompressible fluid, p is a pole distance upstream of the cylinder's center and U G the constant upstream velocity of the parallel horizontal main flow in the x 1 -direction. It is highly satisfying that this result reveals the same large-scale difference in space, x 2 − x 2min , and mean velocity difference, u 1max − u 1 , which also appear in Equation (5); but more than this, if Equation (5) is applied to the wake flow problem and an order-of-magnitude approximation is taken into consideration, the DQTM rigorously produces this Reynolds shear stress [24]. To the authors best knowledge, until today the DQTM is the only turbulence model that produces this first-physical-principle Reynolds stress and this even with the right multiplicative constant! The lack of any gradient in this Reynolds shear stress and the occurrence of nonlocal large-scale spatial lengths and averaged velocities, as occurring in Equation (5), show that nonlocal and fractional calculus is the right choice to close turbulent flow problems with a corresponding constitutive equation of turbulence.

The Turbulent Axisymmetric Jet
In this section, we study an axisymmetric turbulent jet flow into a quiescent surrounding, where the x 1 -direction is identical to the round jet's axis (the detailed calculations are in Refs. [24,27]). The flow width of the jet from a virtual origin behind the nozzle linearly increases with b(x 1 ) = βx 1 . Now, the radial spatial coordinate x 2 is normalized by the expression η = x 2 /b(x 1 ). Furthermore, the average downstream velocity is made dimensionless by f 1 (η) = u 1 (x 1 , x 2 )/u 1 (x 1 , 0) with its value on the jet axis, x 2 = 0, where the mean velocity shows its maximum u 1 (x 1 , 0) = u 1max (x 1 ) (see in Figure 4a the right half of the symmetric mean velocity profile with its maximum at η = 0). By combining the two partial differential equations (the continuity equation and the NSE) with the new closure of turbulence (5) which is a difference equation, substituting λ T = x 2max − x 2min and the self-similar coordinate η and function f 1 , the quasi-two-dimensional flow problem becomes one-dimensional which is essential to determine analytical solutions of this turbulent flow problem. By this procedure, the following highly nonlinear ordinary differential equation results: where an apostrophe, , denotes the derivative d/dη. Its solution is the Gaussian function which is self-similar with respect to the dimensionless coordinate in the transverse direction η and the mean downstream velocity and represented by the equation (see Figure 4a). By a substitution of Equation (9) into differential Equation (8), a reader may prove this by her-/himself. Note the extraordinary complexity of the highly nonlinear differential Equation (8) and, on the other hand, the high simplicity of its solution (9). difference in space, − , and mean velocity difference, − , which also appear in Equation (5); but more than this, if Equation (5) is applied to the wake flow problem and an order-of-magnitude approximation is taken into consideration, the DQTM rigorously produces this Reynolds shear stress [24]. To the authors best knowledge, until today the DQTM is the only turbulence model that produces this first-physical-principle Reynolds stress and this even with the right multiplicative constant! The lack of any gradient in this Reynolds shear stress and the occurrence of nonlocal large-scale spatial lengths and averaged velocities, as occurring in Equation (5), show that nonlocal and fractional calculus is the right choice to close turbulent flow problems with a corresponding constitutive equation of turbulence.

The Turbulent Axisymmetric Jet
In this section, we study an axisymmetric turbulent jet flow into a quiescent surrounding, where the -direction is identical to the round jet's axis (the detailed calculations are in Refs. [24,27]). The flow width of the jet from a virtual origin behind the nozzle linearly increases with ( ) = . Now, the radial spatial coordinate is normalized by the expression = / ( ). Furthermore, the average downstream velocity is made dimensionless by ( ) = ( , )/ ( , 0) with its value on the jet axis, = 0, where the mean velocity shows its maximum ( , 0) = ( ) (see in Figure 4a the right half of the symmetric mean velocity profile with its maximum at η = 0). By combining the two partial differential equations (the continuity equation and the NSE) with the new closure of turbulence (5) which is a difference equation, substituting = − and the self-similar coordinate and function , the quasi-two-dimensional flow problem becomes one-dimensional which is essential to determine analytical solutions of this turbulent flow problem. By this procedure, the following highly nonlinear ordinary differential equation results: where an apostrophe, ′, denotes the derivative d/dη. Its solution is the Gaussian function which is self-similar with respect to the dimensionless coordinate in the transverse direction and the mean downstream velocity and represented by the equation (see Figure 4a). By a substitution of Equation (9) into differential Equation (8), a reader may prove this by her-/himself. Note the extraordinary complexity of the highly nonlinear differential Equation (8) and, on the other hand, the high simplicity of its solution (9).  The dimensionless Reynolds shear stress is directly derivable from Equation (5) and is which is shown as Figure 4b. In the figure inset "D" denotes direct experimentally measured values of − f 21 and "I" are measurements of f 1 , shown in Figure 4a, which, by applying Equation (10), describing the Reynolds shear stress, were transformed to give indirectly (I) further values of − f 21 . The mean velocity f 1 and the Reynolds shear stress f 21 show an excellent agreement between the theoretical data [24] and the experimental observations by Wygnanski and Fiedler [31] (see Figure 4a,b). Also, excellent quality is obtained for the three normal Reynolds stresses, the turbulent kinetic energy, the production and dissipation of turbulent kinetic energy and the turbulent convection (see References [12,13,27]).

Plane Turbulent Couette Flow
This type of flow is created by shearing two parallel plane plates of distance 2a with a fluid in between by keeping the lower plate at rest and moving the higher positioned one with velocity U > 0 in the positive x 1 -direction. In the laminar case this is the prototype flow to demonstrate Newton's law of viscosity. It is well-known that the continuity and a reduced and linearized Navier-Stokes equation, supplemented by the constitutive equation of Newton (Newton's law of viscosity), lead to a linear differential equation and a linear velocity profile between the two different velocities of the plates (see Table 1). Note that (because of the no-slip boundary condition) the plate velocities are identical to the fluid velocities at the boundaries of these plates, where the dimensionless mean velocity g 1 fulfils the boundary conditions given as Equation (11). In these equations η = x 2 /a, where a is the half distance between the two plates. If instead of Newton's law of viscosity the DQTM is implemented-in this more general case-into the nonlinear, self-similar NSE, one obtains the following quadratic differential equation with the two boundary conditions (these and further following results and their calculations are found in [28]) Differential equation Char. velocity u mol u 1 − u 1min Boundary conditions Mean stream-wise velocity The solutions of the time-averaged turbulent velocity profiles are S-shaped curves (an example is shown in Figure 5a by the full line). With the dimensionless averaged downstream velocity g 1 (η) and the relation χ = β/4, it is described by the following equation:  Figure 6, respectively) with experimental data from Reference [32]. The mean velocity profile is nonlinear and S-shaped (here in a mirrored picture). In the infinite Reynolds number limit the distribution is a constant ≡ 1/2, fulfilling the boundary conditions (11) by the occurrence of two discontinuities at the boundaries. (b) In the second figure a comparison of the theoretical Reynolds stress (see Equation (14) below with β = 3.96 and α from Equation (13) or Figure  6, respectively) with numerical simulation results from [33] is presented. In the infinite Reynolds number limit the distribution is a constant ≡ 1, fulfilling the boundary conditions (±1) = 0 by the occurrence of two discontinuities at = ±1. Reprinted with permission from Egolf and Weiss (1995a) © 2021, American Physical Society.
By decreasing the Reynolds number, Re, from high values down toward the critical Reynolds number, , these mean velocity profiles loose curvature and converge to the linear laminar profile (see Figure 5a, the linear function is represented by a dashed line). One may easily show that for the Reynolds shear stress the DQTM directly leads to the formula = 4 ( ) [1 − ( ) ]. By a substitution of the right-hand side of Equation (12) into this equation, it follows that which is shown together with numerical data from Reference [33] as Figure 5b.

Turbulence: A Critical Phenomenon
In turbulent flow fields universality of rare fluctuations were experimentally observed and a lack of theory of critical phenomena of turbulence [34,35] was occasionally complained (see e.g., Refs. [36,37]). Analytical solutions, derived by an application of the DQTM, in a natural manner reveal theories of critical phenomena with criticalities and in which the two phases of the fluid are the low-entropy vorticity-rich regions and the high-entropy laminar streaks. Transition from laminar to transitional and turbulent flows at critical Reynolds numbers are known since the time of Osborne Reynolds (1842-1912). Theoretical vortisation curves * (* The name vortisation was proposed by the authors [25] from the word vorticity, in analogy to the name magnetisation which both are order parameters as function of related stress parameters of the respective physical systems.), in analogy to magnetisation  Figure 6, respectively) with experimental data from Reference [32]. The mean velocity profile is nonlinear and S-shaped (here in a mirrored picture). In the infinite Reynolds number limit the distribution is a constant g 1 ≡ 1/2, fulfilling the boundary conditions (11) by the occurrence of two discontinuities at the boundaries. (b) In the second figure a comparison of the theoretical Reynolds stress (see Equation (14) below with β = 3.96 and α from Equation (13) or Figure 6, respectively) with numerical simulation results from [33] is presented. In the infinite Reynolds number limit the distribution is a constant g 21 ≡ 1, fulfilling the boundary conditions g 21 (±1) = 0 by the occurrence of two discontinuities at η = ±1. Reprinted with permission from Egolf and Weiss (1995a) © 2021, American Physical Society. curves in magnetism, have been determined for numerous turbulent flow types. Figure 6 presents the vortisation curve of plane turbulent Couette flow. Figure 6. Order parameter χ, determined by Equation (13) and compared with experimental results from Reference [32]. There is a lack of good experimental results at excitations just slightly above the criticality . This is explained by the inverse relation between α and Re and the difficulty for experimentalists to perform good experimental set-ups and experiments to measure the turbulent fields at values just slightly above the critical Reynolds number (for more details see in Reference [28]). Reprinted with permission by Egolf and Weiss (1995a) © 2021, American Physical Society.
In References [13,38], in analogy to the mean field theory of magnetic materials, a mean field theory of turbulence was developed. The results are meaningful and accurate. This is explained by the fact that a quasi-steady turbulent flow, with a constant turbulent kinetic energy flow into and out of the system's Kolmogorov energy cascade, leads to behavior that is very similar to that of physical systems in thermodynamic equilibrium. In analogy, a Curie law of turbulence and a Curie-Weiss law of turbulence were discovered. From the Curie law of turbulence, in analogy to the Curie law of magnetism, the turbulence in- Figure 6. Order parameter χ, determined by Equation (13) and compared with experimental results from Reference [32]. There is a lack of good experimental results at excitations just slightly above the criticality Re c . This is explained by the inverse relation between α and Re and the difficulty for experimentalists to perform good experimental set-ups and experiments to measure the turbulent fields at values just slightly above the critical Reynolds number (for more details see in Reference [28]). Reprinted with permission by Egolf and Weiss (1995a) © 2021, American Physical Society.
In the mathematical process of determining solution (12), the relation between α and β is revealed. The parameter α is the stress parameter (control parameter) of the flow system, which is tightly related to the Reynolds number, and the second parameter β, respectively χ (normalized β), is the systems order parameter (see Figure 6). This mathematically derived order parameter relation is By decreasing the Reynolds number, Re, from high values down toward the critical Reynolds number, Re c , these mean velocity profiles loose curvature and converge to the linear laminar profile (see Figure 5a, the linear function is represented by a dashed line). One may easily show that for the Reynolds shear stress the DQTM directly leads to the formula g 21 = 4χ g 1 (η) [1 − g 1 (η)]. By a substitution of the right-hand side of Equation (12) into this equation, it follows that which is shown together with numerical data from Reference [33] as Figure 5b.

Turbulence: A Critical Phenomenon
In turbulent flow fields universality of rare fluctuations were experimentally observed and a lack of theory of critical phenomena of turbulence [34,35] was occasionally complained (see e.g., Refs. [36,37]). Analytical solutions, derived by an application of the DQTM, in a natural manner reveal theories of critical phenomena with criticalities and in which the two phases of the fluid are the low-entropy vorticity-rich regions and the high-entropy laminar streaks. Transition from laminar to transitional and turbulent flows at critical Reynolds numbers are known since the time of Osborne Reynolds (1842-1912). Theoretical vortisation curves * (* The name vortisation was proposed by the authors [25] from the word vorticity, in analogy to the name magnetisation which both are order parameters as function of related stress parameters of the respective physical systems.), in analogy to magnetisation curves in magnetism, have been determined for numerous turbulent flow types. Figure 6 presents the vortisation curve of plane turbulent Couette flow.
In References [13,38], in analogy to the mean field theory of magnetic materials, a mean field theory of turbulence was developed. The results are meaningful and accurate. This is explained by the fact that a quasi-steady turbulent flow, with a constant turbulent kinetic energy flow into and out of the system's Kolmogorov energy cascade, leads to behavior that is very similar to that of physical systems in thermodynamic equilibrium. In analogy, a Curie law of turbulence and a Curie-Weiss law of turbulence were discovered. From the Curie law of turbulence, in analogy to the Curie law of magnetism, the turbulence intensity has been determined to be the response function of turbulence, as, for example, compressibility is of a gas and susceptibility of a para-ferromagnetic phase-change material. Note that results of fluid dynamic considerations-involving the Difference-Quotient Turbulence Model (DQTM)-leading to critical phenomenon, are of higher physical relevance than results extracted only from an ad hoc model like the mean field theory of turbulence.

The Analogy between Laminar and Turbulent Flows
First order turbulence modeling has the great advantage that it follows the widely accepted method of introducing a generally valid constitutive equation to solve physical problems, as it is convention in other areas of physics. This is also the method of treating laminar flows, namely by applying Newton's law of viscosity. In turbulence, we propose the introduction of a new nonlinear, nonlocal and fractional constitutive equation, which was discovered in April 1985 and is called the Difference-Quotient Turbulence Model (DQTM). Some results stemming from the application of this model to the elementary plane turbulent Couette flows are presented in Section 4.7.4. They show strong analogies to the laminar case which in this section both are outlined and compared in Table 1.
As may be seen in Figure 6, at criticality Re = Re c ⇒ Re/Re c = 1, the order parameter is χ = 0. If these values are introduced into the differential equation describing turbulent plane Couette flows (see in Table 1, third line, sixth row), it follows that dg 1 /dη = 1/2. This is the differential equation for laminar plane Couette flows (third line, fifth row) and leads to the correct linear velocity profile g 1 = 1/2 (1 + η) (fifth line, fifth row), fulfilling the boundary conditions (11) (see also in Table 1, fourth line, fifth and sixth row).

Conclusions
This article promotes nonlocal and fractional calculus to model Reynolds shear stresses with a first order turbulence closure scheme and makes a new proposal for a constitutive equation of turbulence. This equation, in the limit of the Reynolds number converging to its criticality from above, conserves the full information of Newton's law of viscosity of laminar flows. It also extends Boussinesq's and Prandtl's models of turbulent momentum transfer by a generalization from eddies of a single size to such of a cascade with numerous classes of fractal eddies of multiple sizes. Six applications of the Difference-Quotient Turbulence Model, DQTM (constitutive equation of turbulence) to elementary turbulent flow problems lead to analytical results which, compared with experimental and (direct) numerical simulation results, are of convincing precision. All the derived non-empirical analytical solutions are of high simplicity and accuracy.