Turbulence revisited

Turbulence revisited

Chaos, Solitons and Fractals 41 (2009) 1136–1149 Contents lists available at ScienceDirect Chaos, Solitons and Fractals journal homepage: www.elsevi...

460KB Sizes 0 Downloads 136 Views

Chaos, Solitons and Fractals 41 (2009) 1136–1149

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals journal homepage: www.elsevier.com/locate/chaos

Turbulence revisited Michail Zak Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, United States

a r t i c l e

i n f o

Article history: Accepted 28 April 2008

a b s t r a c t A new approach to turbulence based upon the Stabilization Principle is introduced. Onset of turbulence is interpreted as loss of stability of solutions to the Navier–Stokes equations in the class of differentiable functions. Developed turbulence is considered as postinstability motion in the class of non-differentiable functions. A non-linear version of the Liouville equation is proposed for describing postinstability motions of dynamical systems with exponential divergence of trajectories such as those leading to chaos and turbulence. As a result, velocities are represented by a set of stochastic invariants found from the Stabilization Principle. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Unsolved problems in Newtonian mechanics Turbulence is one of the most fundamental problems in theoretical physics that is still unsolved. Although applicability of the Navier–Stokes equations as a model for fluid mechanics is not in question, the instability of their solution for flows with supercritical Reynolds numbers raises a more general question: Is Newtonian mechanics complete? During several centuries, Newtonian mechanics enjoyed an unprecedented success being considered as the most general, the most accurate, and the most elegant branch of science. It became consensual that its structure is perfect and complete as far as the non-relativistic macroworld is concerned. However the problem of turbulence (stressed later by the discovery of chaos) demonstrated that the Newton’s world is far more complex than those represented by classical models. It appears that the Lagrangian or Hamiltonian formulations do not suggest any tools for treating postinstability motions, and this is a major flaw of the classical approach to Newtonian mechanics. The explanation of that limitation was proposed by Zak [13,14]: the classical formalism based upon the Newton’s laws exploits additional mathematical restrictions (such as space–time differentiability, and the Lipschitz conditions) that are not required by the Newton’s laws. The only purpose for these restrictions is to apply a powerful technique of classical mathematical analysis. However, in many cases such restrictions are incompatible with physical reality. This statement can be illustrated by a trivial example presented below, Fig. 1. Consider a surface of a tangential jump of velocity V2  V1 in a horizontal unidirectional flow of an inviscid incompressible fluid. Applying the principle of virtual work to a small volume V containing both flows of the fluid and the surface of the tangential jump of velocities separating the flows, one obtains

Z

ðqa1  dU 1 þ qa2  dU 2 ÞdV ¼ 0;

ð1Þ

V

where q, U1, U2, a1, a2 are density, displacements, and accelerations of the fluid. The displacements U1 and U2 are mutually independent in the region that does not contain the separating surface, but they are dependent at the surface due to its impenetrability

E-mail address: [email protected] 0960-0779/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2008.04.059

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1137

Fig. 1. Loss of differentiability in fluid mechanics.

U 1  n ¼ U 2  n ¼ U:

ð2Þ

Hence, as follows from Eq. (1), at the surface the following equality holds:

"

2  2 # o o o o U þ a ¼ 0; þ þ V1 þ V2 ot oS ot oS

ða1 þ a2 Þ  n ¼ 0; i:e:

ð3Þ

where a is the term that does not contain the second order derivatives of U, and therefore, it does not effect the characteristic equation

ðk  V 1 Þ2 þ ðk  V 2 Þ2 ¼ 0

ð4Þ

and its characteristic roots



1 ½ðV 2 þ V 1 Þ  iðV 2  V 1 Þ: 2

ð5Þ

Let us start with studying the propagation of high frequency oscillations of the transverse displacements U. Recall that Eq. (3) being linear with respect to the second order time–space derivatives, strictly speaking, is non-linear with respect to U and its first time–space derivatives that are contained in the term a. For small amplitudes and their first derivatives, this term can be linearized:

a ¼ a1

oU oU þ a2 þ a3 U þ a4 : ot oS

ð6Þ

For further simplifications, all the coefficients in Eq. (3) can be linearized with respect to an arbitrarily selected point S0 and instant of time t0 = 0. Then Eq. (3) takes the form of a linear elliptic PDE with constant coefficients

"

# 2  2 o o 0 o 0 o 0 o 0 o 0 þ þ a1 þ a2 þ a3 U þ a04 ¼ 0: þ V1 þ V2 ot oS ot oS ot oS

ð7Þ

Obviously this equation is valid only for small amplitudes with small first derivatives, the small area around the above selected point S0, and within a small period of time Dt. Let us derive the solution to Eq. (7) subject to the following initial conditions:

U 0 ¼

1 k0 Si e ; k0

at t ¼ 0

ð8Þ

assuming that k0 can be made as large as desired, i.e. k0 > N ? 1. Consequently, the initial disturbances can be made as small as desired, i.e. U 0 < N 1 ! 0. The corresponding solution is written in the form

U  ¼ C 1 ek0 ðk1 tSÞi þ C 2 ek0 tSÞi ;

ð9Þ

where k1 and k2 are the roots of the characteristic Eq. (4) (see Eq. (5)). Since the characteristic roots are complex, the solution (9) will contain the term

1 jIm k1;2 jDt e sin k0 S; k0

k0 ! 1

ð10Þ

that leads to infinity within an arbitrarily short period of time Dt and within an infinitesimal area around the point S0. Hence one arrived at the following situation: jU*j ? 1 in spite of the fact that jU0j ? 0. To obtain a geometrical interpretation of the above described instability, let us note that if the second derivatives in Eq. (7) are of order k0, then the first derivatives are of order of 1, and U is of order of 1/k0. Hence, the period of time Dt can be selected in such a way that the second derivatives will be as large as desired, but U and its first derivatives are still sufficiently small. Taking into account

1138

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

that the original governing equation (3) is quasi-linear with respect to the second derivatives, and therefore, the linearization does not impose any restrictions on their values, one concludes that the linearized equation (7) is valid for the solution during the above-mentioned period of time Dt. Turning to the term (10) of the solution (9), one can now interpret it as being represented by a function having an infinitesimal amplitude and changing its sign with an infinite frequency (k0 ? 1). The first derivatives of this function can be small and change their signs by finite jumps (with the same infinite frequency), so that the second derivatives at the points of such jumps are infinite. From the mathematical viewpoint, this kind of function is considered as continuous, but non-differentiable. The result formulated above was obtained under specially selected initial conditions (8), but it can be generalized to include any initial conditions. Indeed, let the initial conditions be defined as

Ujt¼0 ¼ U  0 ðXÞ

ð11Þ

and the corresponding solution to Eq. (3) is

U ¼ U  ðX; tÞ:

ð12Þ

Then, by altering the initial conditions to

Ujt¼0 ¼ U 0 ðXÞ þ U  0 ðXÞ;

ð13Þ

U 0

where is defined by Eq. (8), one observes the preceding argument by superposition that vanishingly small change in the initial condition (13) leads to unboundedly large change in the solution (12) that occurs during an infinitesimal period of time. That makes the solution (12) non-differentiable, but still continuous. However, the Euler’s equation that governs the flow discussed above was derived under condition of space–time differentiability of the velocity field. This discrepancy manifests inapplicability of the Euler’s equation for description of postinstability motion of an inviscid fluid, and therefore, for a developed turbulence. But this does not imply the incompleteness of Newtonian mechanics: it only means that the mathematical formalism that expresses the Newton’s laws should include non-differentiable components of the velocity field. Incorporation of thermodynamics into Newtonian mechanics does not eliminate the loss of differentiability: the Navier–Stokes equations appear to have even more sophisticated patterns of instability than the Euler equations if the Reynolds number is supercritical, and from the point of view of mathematical formalism, it leads to the phenomenon described above. As a matter of fact, the Navier–Stokes equations impose even stronger mathematical restriction on the velocity field requiring it’s twice differentiability with respect to space coordinates. Therefore, the Navier–Stokes equations require a modification that would allow one to include non-differentiable velocity field similar to those in the Euler’s equations. It should be noticed that the same kind of phenomena occurs in other branches of continuum mechanics, and in particular, in theory of elasticity. Indeed, consider an ideal filament in the gravity field stretched in the vertical direction with the additional tension T0 as shown in Fig. 2. Let us cut it at the middle point and observe a behavior of its lower part. As shown by Zak [13,17,18] the lower part has imaginary characteristic speeds

sffiffiffi S ; k ¼ i g

ð14Þ

where S is the current length (see Fig. 2), and the solution to the corresponding governing equation contains the same destabilizing term (10) as those in Eq. (9). As shown by Zak [13,17,18], the same loss of differentiability occurs in two- and three-dimensional elastic bodies (wrinkles in films, buckling in soft shells, and fractures in composite materials). In mathematics, all these phenomena are associated with the Hadamard’s instability. Remark 1. Returning to Fig. 1, it should be noticed that the motion of the upper part of the stretched and cut-at-the-middle string is also non-trivial. As shown by Zak [9,13,17,18], the configuration of the string remains differentiable within the open

Fig. 2. Loss of differentiability in elasticity.

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1139

interval that excludes the free end; however, the ‘‘price” for such a mathematical convenience is a loss of the snap-of-a-whip effect well known from experiments. Only inclusion of the free end allows one to observe this effect, but that leads to violation of another mathematical restriction: the Lipschitz condition. Therefore, two mathematical restrictions: space–time differentiability and the Lipschitz condition that do not follow from the Newton’s laws are imposed on the formalism of Newtonian mechanics, and that leads to its apparent incompleteness. A removal of both of these restrictions will be discussed below. 1.2. Reynolds approach Over 100 years ago, Reynolds proposed a modification of the Navier–Stokes equations by decomposing the fluid velocity  and fluctuating component v0 [6]. The new form of the Reynolds-averaged Nainto the mean (time-averaged) component v vier–Stokes (RANS) equations for an incompressible fluid is the following:

i ov ¼ 0; oxi      v  ov ov o ovi ovj dij þ l q i þ q i j ¼ qfi þ p þ  qvi vj Þ : oxj ot oxj oxj oxi

ð15Þ ð16Þ

; andl are density, mean pressure, and coefficient of viscosity, while Eqs. (5) and (6) represent conservation of mass Here q; p and momentum, respectively. The last term in Eq. (16) (generally referred as Reynolds stresses) represents the contribution of the inertial forces of fluctuations into momentum an energy of a turbulent flow. However, along with the qualitative advantages, this term brings a fundamental quantitative disadvantage: the system Eqs. (15) and (16) becomes open since the number of unknowns now is larger than the number of equation; that creates the so called ‘‘closure” problem that is to find additional relationships between fluctuations and mean velocities. A typical closure was proposed by Prandtl 2

vi vj ¼ l j

i ov i ov j ; oxj oxj

ð17Þ

where a non-local parameter l (called mixing length) is supposed to be found from experiments. However, the problem with this and other closures is that they are effective only for a narrow class of boundary conditions that are close to those for which the corresponding experiment was performed, while they fail otherwise. This failure could be expected. Indeed, Reynolds did not change the physics of the viscous fluid: he changed only a mathematical representation of this physics. Namely, he enlarged the class of functions in which turbulence should be described by moving from differentiable to multi-valued functions. Therefore, a closure should be searched not in new physical laws, but rather in relationships between fluctuations and the rate of instability that causes the failure of differentiability. These relationships follow from the Stabilization Principle introduced by Zak [10,13]. It will be discussed in the next section. Remark 2. More consistent and elegant way to enlarge the class of differentiable functions in physics has been discussed by Nottale [5] in his review paper. The approach is based upon the theory of scale relativity. The alternative to smoothness is fractal dimensionality, while non-differentiability is represented by stochastic variables describing fluctuations. However, this approach requires the same kind of reference to the rate of instability that the Reynolds equations do. Indeed, according to the governing equations written in scale relativity form, the non-smooth space–time components of motion can be generated only by non-smooth force or non-smooth boundary–initial conditions. The only source of such conditions in a laminar flow is inertial force of unstable fluctuations. Therefore, prior to application of the relative scale equations, one has to perform the stability analysis of the laminar flow and find the fluctuations. 1.3. Stabilization Principle A physical quantity is invariant if it has the same magnitude in all inertial reference frames. This definition disqualifies the concept of stability from being a physical invariant. Indeed, it can be shown that stability depends upon the frame of reference, upon the class of functions in which the motion is described, and upon the metrics of configuration space, and in particular, upon the way in which the distance between the basic and perturbed solution is defined. The dependence upon the frame of reference can be illustrated by the following example [1]: consider an inviscid stationary flow with a smooth velocity field

vx ¼ A sin z þ C cos y;

vy ¼ B sin x þ A cos z;

vz ¼ C sin y þ b cos x:

ð18Þ

Surprisingly, the trajectories of this flow are unstable (Lagrangian turbulence). It means that this flow is stable in the Eulerian coordinates, but is unstable in the Lagrangian coordinates. Several examples of dependence of the stability of a motion upon the metric of the configuration space were demonstrated by Synge [8]. Finally, the dependence of the stability of a motion upon the class of functions was illustrated by Zak [9,13].

1140

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

Thus, the concept of stability is an attribute of a mathematical formalism rather than physics, and consequently, unsolved problems in Newtonian mechanics do not imply an incompleteness of the Newton’s laws: they manifest inadequacy of the corresponding mathematical formalism, i.e. inconsistency between idealized models and reality. 1.4. Stabilization Principle for closure problem Based upon comments made above, we will now formulate the stabilization principle for closure of the Reynolds-averaged equations. Consider a dynamical model that in some domain of its parameters loses its space–time differentiability, i.e. it becomes unstable in the class of differentiable functions. As noticed earlier, this means that the corresponding physical phenomenon cannot be adequately described without an appropriate modification of the mathematical formalism representing the model. The modification should be based upon an enlarging the original class of functions in such a way that the instability is eliminated. The mathematical formulation of this principle can be expressed in the following symbolic form:

AR  X ¼ A0  X þ rR :

ð19Þ

If the original dynamical model

A0  X ¼ 0

ð20Þ

is unstable, but the Reynolds-averaged model

A0  X þ rR ¼ 0

ð21Þ

is stable, obviously the stabilization is performed by the Reynolds stresses rR that represent contribution of all the nonsmooth components into the modified model. Indeed, driven by the mechanism of instability of the original model, they grow until the instability is suppressed down to a neutral stability. From that viewpoint, the Prandil’s closure equation (17) can be considered as a feedback that stabilizes an originally unstable laminar flow. Turning, for instance, to a plane Poiseuille flow with a parabolic velocity profile, one arrives at its instability if the Reynolds number is larger than R ffi 5772. Experiments show that the new steady turbulent profile is no longer parabolic: it is very flat near the center and is very steep near the walls. The same profile follows from the Prandtl solution based upon the closure (17). But since this profile can be experimentally observed, it must be stable, and this stabilization is carried out by the feedback (17). Mathematical justification of the neutral stability is presented by Zak [13]. Experimental verification of neutral stability of free turbulent jets was reported in Lessen and Paillet [3]. The meaningfulness of the Stabilization Principle formulated above and its application to the closure problem had been illustrated for inviscid and viscous turbulent flows by Zak [13,14]. However, a major limitation of this approach is a necessity to find a rate of instability of the original laminar flow prior to application of the Stabilization Principle, and this pre-condition is very complex and laborious. 2. Extension of the Liouville formalism to postinstability dynamics 2.1. Background Newtonian dynamics describes processes in which future can be derived from past, and past can be traced from future by time inversion. Because of such determinism and reversibility, classical dynamics becomes fully predictable, and therefore, it cannot explain the emergence of new dynamical patterns in nature. Dissipative version of Newtonian dynamics that includes elements of thermodynamics is not reversible, but is still fully deterministic and predictable. The only phenomenon which breaks down determinism and creates ‘‘intrinsic stochasticity” is instability of the mathematical models that may lead to chaos and turbulence. The major limitation of Newtonian dynamics is that its mathematical formalism does not discriminate between stable and unstable motions, and therefore, an additional stability analysis is required for that. However, such an analysis is not constructive: in case of instability, it does not suggest any model modifications to efficiently describe postinstability motions. This flaw in classical dynamics has attracted the attention of many outstanding scientists (Gibbs, Planck, and Prigogine, among others). From the mathematical viewpoint, the instability of solutions manifests incompleteness of the corresponding model. For instance, the derivation of the Navier–Stokes equation (as well as other models of continua) is based upon the assumption that the velocity field is twice-differentiable. Obviously, this assumption is violated when, driven by instability, a laminar flow is turning into a turbulent one. However, instability is the attribute of a mathematical model rather than a physical phenomenon, i.e. it is associated with a certain class of function in which the solution is sought. Indeed, a stochastic description of the same turbulent flow does not suffer any instability, and this gives a hint: enlarge the class of function to include stochasticity into the original mathematical model. In order to do that, we will use another hint: turning to quantum mechanics and representing the Schrödinger equation in the Madelung form, we observe the feedback from the Louville equation to the Hamilton–Jacobi equation in the form of the quantum potential [2]. Preserving the same topology, we will replace the quantum potential by other types of feedbacks to introduce intrinsic stochasticity in the governing equation of Newtonian dynamics, and thereby, to enlarge the class of functions in which the Newton’s laws are to be implemented.

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1141

2.2. Destabilizing effect of Liouville feedback We will start with derivation of an auxiliary result that illuminates departure from classical models. For mathematical clarity, we will consider here a one-dimensional motion of a unit mass under action of a force f depending upon the velocity v and time t

v_ ¼ f ðv; tÞ:

ð22Þ

If initial conditions are not deterministic, and their probability density is given in the form

q0 ¼ q0 ðVÞ; where q P 0; and

Z

1

q dV ¼ 1;

ð23Þ

1

while q is a single-valued function, then the evolution of this density is expressed by the corresponding Liouville equation

oq o ðqf Þ ¼ 0: þ ot oV

ð24Þ

The solution of this equation subject to initial conditions and normalization constraints (23) determines probability density as a function of V and t: q = q(V, t). In order to deal with the constraint (23), let us integrate Eq. (24) over the whole space assuming that q ? 0 at jVj ? 1 and jfj < 1. Then

o ot

Z

1

q dV ¼ 0;

1

Z

1

q dV ¼ const:

ð25Þ

1

Hence, the constraint (25) is satisfied for t > 0 if it is satisfied for t = 0. Let us now specify the force f as a feedback from the Liouville equation

f ðv; tÞ ¼ u½qðv; tÞ

ð26Þ

and analyze the motion after substituting the force (26) into Eq. (22)

v_ ¼ u½qðv; tÞ:

ð27Þ

This is a fundamental step in our approach. Although the theory of ODE does not impose any restrictions upon the force as a function of space coordinates, the Newtonian formalism does: equations of motion are never coupled with the corresponding Liouville equation. Moreover, it can be shown that such a coupling leads to new properties of the underlying model. Indeed, substituting the force f from Eq. (26) into Eq. (25), one arrives at the non-linear equation for evolution of the probability density

oq o fqu½qðV; tÞg ¼ 0: þ ot oV

ð28Þ

Let us now demonstrate the destabilizing effect of the feedback (27). For that purpose, it should be noted that the derivative oq/ov must change its sign, at least once, within the interval 1 < v < 1, in order to satisfy the normalization constraint (23). But since

Sign

ov_ du oq Sign ¼ Sign ov dq ov

ð29Þ

there will be regions of v where the motion is unstable, and this instability generates randomness with the probability distribution guided by the Liouville equation (28). It should be noticed that the condition (29) may lead to exponential or polynomial growth of v (in the last case the motion is called neutrally stable; however, as will be shown below, it causes the emergence of randomness as well if prior to the polynomial growth, the Lipschitz condition is violated). 2.3. Emergence of randomness In order to illustrate mathematical aspects of the concepts of Liouville feedback, as well as associated with it instability and randomness, let us take the feedback (26) in the form

f ¼ r2

o ln q ov

ð30Þ

to obtain the following equation of motion:

v_ ¼ r2

o ln q: ov

ð31Þ

This equation should be complemented by the corresponding Liouville equation (in this particular case, the Liouville equation takes the form of the Fokker–Planck equation)

1142

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

oq o2 q ¼ r2 2 : ot oV

ð32Þ

Here v stands for the particle velocity, and r2 is the constant diffusion coefficient. The solution of Eq. (32) subject to the sharp initial condition is

! 1 V2 q ¼ pffiffiffiffiffiffi exp  2 : 4r t 2r p t

ð33Þ

Substituting this solution into Eq. (22) at V = v, one arrives at the differential equation with respect to v(t)

v_ ¼

v 2t

ð34Þ

and, therefore,

pffiffi v¼C t

ð35Þ

Here C is an arbitrary constant. Since v = 0 at t = 0 for any value of C, the solution (35) is consistent with the sharp initial condition for the solution (33) of the corresponding Liouville equation (22). The solution (35) describes the simplest irreversible motion: it is characterized by the ‘‘beginning of time” where all the trajectories intersect (that results from the violation of Lipschitz condition at t = 0, Fig. 3), while the backward motion obtained by replacement of t with (t) leads to imaginary values of velocities. One can notice that the probability density (33) possesses the same properties (see Fig. 4). For a fixed C, the solution (35) is unstable since

dv_ 1 ¼ >0 dv 2t

ð36Þ

and, therefore, an initial error always grows generating randomness. Initially, at t = 0, this growth is of infinite rate since the Lipschitz condition at this point is violated

dv_ ! 1 at t ! 0: dv

ð37Þ

This type of instability had been introduced and analyzed by Zak [11,12]. Considering first Eq. (35) at fixed C as a sample of the underlying stochastic process (33), and then varying C, one arrives at the whole ensemble characterizing that process (see Fig. 3). One can verify that, from Eq. (33) [7] the expectation and the variance of this process are, respectively,

MV ¼ 0;

DV ¼ 2r2 t:

ð38Þ

The same results follow from the ensemble (35) at 1 6 C 6 1. Indeed, the first equality in (38) results from symmetry of the ensemble with respect to v = 0; the second one follows from the fact that:

DV / v2 / t:

ð39Þ

It is interesting to notice that the stochastic process (35) is an alternative to the following Langevin equation [7]

v_ ¼ CðtÞ;

M C ¼ 0;

DC ¼ r

ð40Þ

Fig. 3. Stochastic process and probability density.

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1143

Fig. 4. Effect of terminal attractor.

that corresponds to the same Fokker–Planck equation (22). Here C(t) is the Langevin (random) force with zero mean and constant variance r. 2.4. Negative diffusion Let us consider the Langeven equation (40) and couple it with the corresponding Liouville equation in the same fashion as in Eq. (31), and present this equation in a dimensionless form assuming that the parameter a absorbs the velocity and the time scale factors

v_ ¼ CðtÞ þ a

o ln q; ov

CðtÞCðt 0 Þ ¼ qdðt  t 0 Þ:

ð41Þ

If one chooses a = q2, then the corresponding Liouville equation (that takes the form of the Fokker–Planck equation) will change from Eq. (32) to the following:

 2  oq o2 q o q oq ¼ 0; q ¼ q2 2  oV ot q oV oV

q ¼ q0 ðVÞ ¼ const:

ð42Þ

Thus, the Liouville feedback stops the diffusion. However, the feedback force can be even more effective: it can reverse the diffusion process and push the probability density back to the sharp value in finite time. Indeed, suppose that in the Liouville feedback

pffiffiffiffi

a ¼ q2 exp D; where DðtÞ ¼

Z

1

qV 2 dV:

ð43Þ

1

Then the Fokker–Planck equation takes the form

pffiffiffiffi o2 q oq ¼ ½q2 1  exp DÞ 2 : ot oV

ð44Þ

Multiplying Eq. (44) by V2, then integrating it with respect to V over the whole space, one arrives at ODE for the variance D(t)

pffiffiffiffi D_ ¼ 2q2 ð1  exp DÞ; i:e: D_ 6 0 if D P 0:

ð45Þ

Thus, as a result of negative diffusion, the variance D monotonously vanishes regardless of the initial value D(0). It is interesting to note that the time T of approaching D = 0 is finite



1 q2

Z

0

Dð0Þ

dD

1 pffiffiffiffi 6 1  exp D q2

Z

1 0

dD p pffiffiffiffi : ¼ exp D  1 3q2

ð46Þ

This terminal effect is due to violation of the Lipschitz condition, [12] at D = 0. Let us turn to a linear version of Eq. (44)

oq o2 q ¼ q2 2 ot oV

ð47Þ

and discuss a negative diffusion in more details. From the linear equivalent of Eq. (45)

dD_ ¼ q2 ; i:e: D ¼ D0  q2 t < 0 at t > D0 =q2 : dD

ð48Þ

Thus, eventually the variance becomes negative, and that disqualifies Eq. (47) from being meaningful. It had been shown [16] that the initial value problem for this equation is ill-posed: its solution is not differentiable at any point. (Such an ill-posedness expresses the Hadamard instability similar to those studied by Zak [13].) Therefore, a negative diffusion must be non-lin-

1144

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

ear in order to protect the variance from becoming negative (see Fig. 6). One of possible the realization of this condition is placing a terminal attractor [12] at D = 0 as it was done in Eq. (36). It should be emphasized that negative diffusion represents a major departure from Newtonian formalism. 3. Modified Newtonian formalism 3.1. General model In this section, based upon the stabilizing effect of negative diffusion considered above, we will introduce a general approach to modeling postinstability dynamics. The modified Newtonian formalism is based upon coupling the classical governing equations with the corresponding Liouville equation and suppression exponential divergence of trajectories by the effect of negative diffusion introduced above. The idea of the proposed approach is in introduction such a Liouville feedback that as a fictitious force acts only upon the erratic component of a trajectory without affecting its ‘‘expected” value. For that purpose, introduce a system of n first order ODE with n unknowns

v_ i_ ¼ fi ½fvðtÞg; t;

fvg ¼ v1 ; . . . vn ; i ¼ 1; 2; . . . ; n

ð49Þ

subject to initial conditions

vi ð0Þ ¼ v0i

ð50Þ

Due to finite precision, the values (50) are not known exactly, and we assume that the error possesses some joint distribution

ErrðV 0i Þ ¼ qðV 01 ; . . . ; V 0n Þ ¼ q0 :

ð51Þ

It is reasonable to assume that the initial conditions (50) coincide with the initial expectations i.e. that P0 has a maximum at V 0i ¼ v0i , i = 1, 2, . . ., n. This means that

oq0 ¼ 0; oV 0

o2 q0 < 0; oV i oV j

i ¼ 1; 2; . . . ; n:

ð52Þ

This is true for any symmetric initial density (for instance, the normal distribution) when the expected values have the highest probability to occur. The Liouville equation describing the evolution of the joint density q is

oq þ r  ðqf Þ ¼ 0; ot

f ¼ f1 ; . . . ; fn ;

f i ¼ fi ðfVg; tÞ;

q ¼ qðfVg; tÞ:

ð53Þ

Its formal solution

 Z t  P ¼ P0 exp  r  f ds

ð54Þ

0

suggests that the flattening of the error distribution is caused by the divergence of the trajectories of the governing equation (49) from the target trajectory that starts with the prescribed initial conditions (50), Fig. 5. Let us introduce the following Liouville feedback:

F i ¼ ai

o ln q; ovi

ai > 0:

ð55Þ

Then the system (49) in modified to the following one:

Fig. 5. Divergence of trajectories.

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

v_ i ¼ fi þ ai

o ln q; ovi

1145

ð56Þ

that should be complemented by the corresponding Liouville equation

  oq X o oq ai qfi þ þ ¼ 0: oV i ot oV i i

ð57Þ

The coupled ODE–PDE equations of this type had been studied by Zak [15–18], in connection with model of Livings. Here we will summarize only mathematical aspects of these systems. Firstly, the force Fi makes the Liouville equation non-linear, while ODE becomes dependent upon PDE. Secondly, this force introduces to PDE a negative diffusion that changes the type of the PDE from the hyperbolic to the parabolic one. At the same time, the behavior of the solution to Eq. (57) is fundamentally different from its Fokker–Planck analog. Thirdly, from Eq. (55), the force Fi does not affect the motion along that trajectory vi ¼ vi which has the maximum probability of occurrence since

oq ðvi ¼ vi Þ ¼ 0 ovi

ð58Þ

and that property makes this force fictitious. Before formulating the proposed model in the final form, we will consider a trivial, but instructive example. 3.2. Example Let us consider an unstable linear ODE

v_ ¼ v:

ð59Þ

In this particular case, the expected trajectory is known in advance:

 ¼ 0: v

ð60Þ

However, any small error in initial conditions leads to a different trajectory that diverge exponentially from those in Eq. (59):

v ¼ v0 exp t:

ð61Þ

Similar result follows from the corresponding Liouville equation:

oq o ðqVÞ; ¼ oV ot V ¼ V 0 exp t:

ð62Þ ð63Þ

Let us now introduce the fictitious force as



pffiffiffiffi o D ln q; ov

ð64Þ

where D is the variance

DðtÞ ¼

Z

1

V 2 qðV; tÞdV

ð65Þ

1

and obtain the following modified version of Eq. (59):

v_ ¼ v þ

pffiffiffiffi o ln q: D ov

ð66Þ

Due to this Liouville feedback, Eq. (60) is modified to the following Fokker–Planck equation:

pffiffiffiffi o2 q oq o ðqVÞ  D 2 : ¼ oV ot oV

ð67Þ

Multiplying Eq. (67) by V, then using partial integration, one obtains for expectations the same Eq. (60) and its solution (63). Similarly one obtains for variances

pffiffiffiffi D_ ¼  D:

ð68Þ

For the initial condition

D ¼ D0

at t ¼ 0

the solution to Eq. (68) is

ð69Þ

1146

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149



pffiffiffiffiffiffi 2 pffiffiffiffiffiffi pffiffiffiffiffiffi t D0  for t < 2 D0 and D 0 for t P 2 D0 : 2

ð70Þ

It is easily verifiable that the Lipschitz condition at D = 0 is violated since

oD_ 1 ¼  pffiffiffiffi ! 1 at D ! 0: oD 2 D

ð71Þ

As will be shown later, this property of the solution is of critical importance for multi-dimensional case. Now the solution to the non-linear version of the Fokker–Planck equation in (67) can be approximated by the first term in the Gram–Charlier series represented by the normal distribution with the variance D. For the case close to a sharp initial value at V = 0

1

V2

!

q ¼ pffiffiffiffiffiffiffi exp  2 : D 2p 2D

ð72Þ

Substituting Eq. (72) (with reference to the solution (70)) into Eq. (66) one obtains

" # 1 v_ ¼ v 1  pffiffiffiffiffiffi t D0  2

ð73Þ

whence for v = v0 at t = 0 the solution is



 2 v0 t pffiffiffiffiffiffi t D0  e 0 6 t 6 2D0 ; 2 D0

v 0 t > 2D0 :

ð74Þ

For sufficiently small variance of initial error distribution D0 1, an exponential growth of initial error v0 is totally elimipffiffiffiffiffiffi nated after t > 2 D0 , Fig. 6. It should be noticed that a finite time of approaching equilibrium is special properties of the terminal attractors discussed by Zak [9,11,12]. One has to recall again that although the example we just considered is trivial, the stabilization mechanism performed by the negative diffusion-based Liouville-feedback forces is the same. It is also important to learn from this example that the true expected solution is given by Eq. (74) rather than by Eq. (63) despite the fact that Eq. (63) directly follows from the Liouville equation (67). Indeed, the solution (63) is identical to the original solution (61), and any initial error will grow exponentially. This means that both these solutions are unstable in the class of differentiable functions. But the same physical phenomenon described by Eq. (74) is stable in the enlarged class of functions that includes stochastic components. Obviously the stochastic components are found from the Stabilization Principle discussed in the previous section. 3.3. Final form of modified Newtonian formalism Based upon the example considered above, we can now specify the coefficients ai in Eqs. (56) and (57)

pffiffiffiffiffiffi o Dii ln q; ovi ! X pffiffiffiffiffiffi oq oq o Dii ¼ 0; qfi þ þ ot oV i oV i i

v_ i ¼ fi þ

ð75Þ ð76Þ

Fig. 6. Suppression of instability.

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1147

where Dii are principal variances

Dii ¼

Z

Z

1

1

... 1

ðV i  V i Þ2 qðdVÞn :

ð77Þ

1

In order to verify the stabilizing effect of negative diffusion for the n-dimensional case, let us linearize equation (75) with respect to the initial state vi = 0. Then the linearized versions of Eqs. (75) and (76) will be, respectively,

v_ i ¼ aij vj þ

pffiffiffiffiffiffi o Dii ln q; ovi

 aij ¼

X pffiffiffiffiffiffi oq oq o Dii qaij vj þ þ ot oV i oV i i

!

ofi ovj

 ;

ð78Þ

vj ¼0

¼ 0:

ð79Þ

An n-dimensional analog of Eq. (68) can be obtained by multiplying Eq. (79) by Vi and then using partial integration

D_ ij ¼ ail Dlj  ajl Dli 

qffiffiffiffiffiffi Dij :

ð80Þ

_ ij =oDlk . Its diagonal eleLet us first analyze the effect of terminal attractor and, turning to Eq. (80), start with the matrix ½oD ments become infinitely negative when the variances vanish

oD_ ij ¼ oDij

! 1 2aij  pffiffiffiffiffiffi ! 1 at Dii ! 0; 2 Dij

ð81Þ

while the off-diagonal elements are bounded. Therefore, due to the effect of terminal attractor (81), the system equation (80) has infinitely negative characteristic roots, i.e. it is infinitely stable with respect to small errors regardless of the parameters aij of the original dynamical system (78). In addition to that, the terminal attractor (as well as any attractor) guarantees ‘‘impenetrability” of the state Dii = 0, i.e. if the principle variances initially were non-negative, they will never become negative, and that prevent ill-posedness of the problem for the PDE (76). Thus, all the properties of the modified model discovered in one-dimensional case are preserved in the n-dimensional case, namely: a simultaneous solution of the coupled ODE–PDE system (75) and (76) describes a stable ‘‘expected” motion regardless of the original instability. 3.4. Representation of higher moments Although expected values of the state variables play an important role in description of postinstability motions, they do not expose the full dynamical picture. Indeed, measured velocities of turbulent flows look truly random, and therefore, the behavior of the higher moments is required within the framework of stochastic formalism. For that purpose, let us turn to Eq. (49) and introduce new variables

vij ¼ vi vj :

ð82Þ

After trivial transformations, the system (49) can be rewritten in an equivalent form being expressed via new variables

v_ ij ¼ fij

ðv11 ; . . . ; vnn Þ;

ð83Þ

in which

fij ¼

pffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffi pffiffiffiffiffiffiffi vjj fi ð v11 ; . . . ; vnn g þ vii fj ð v11 ; . . . ; vnn Þ:

ð84Þ

Let us now augment Eq. (84) with the Liouville feedback similar to those in Eq. (75)

v_ ij ¼ fij þ

qffiffiffiffiffiffiffiffi o Dijij ln q : ovij

ð85Þ

Then the corresponding Liouville equation will be similar to (76)

ffi oq X qffiffiffiffiffiffiffi oq o qfij þ Dijij þ oV ij ot oV i ij

! ¼ 0;

ð86Þ

where q* and Dijij are probability density and principal variances for the new variables. Solving Eqs. (85) and (86) simultaneously, one obtains the evolution of the expectations of the new state variables that are equivalent to the second moments of the old variables (see Eq. (82))

ij ¼ vi vj : v

ð87Þ

It should be noticed that q and q* are different: for instance, if initially q is normally distributed, q* must be recalculated by applying the rules for the change of variables (82); that is why the expectations and the second moments must be found from different equations.

1148

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

The higher moments can be found in a similar way by introducing new variables vijk, vijkl, etc. Based upon the expectation and higher moments, one can reconstruct the joint probability distribution of state variables, and therefore, to obtain a complete information about dynamics of the underlying physical process in a stable form. 3.5. Application to postinstability models The proposed approach can be stated as follows: in order to find the solution to a dynamical system (49)

v_ i_ ¼ fi ½fvðtÞg; t;

fvg ¼ v1 ; . . . ; vn ;

i ¼ 1; 2; . . . ; n

ð49Þ

subject to the initial conditions (50)

vi ð0Þ ¼ v0i ;

ð50Þ

and avoid possible computational errors due to exponential divergence of the neighboring trajectories, it is sufficient to modify Eq. (49) by applying a fictitious stabilizing force in the form of the Liouville feedback (55) and to solve the following system

v_ i ¼ fi þ

pffiffiffiffiffiffi o Dii ln q; ovi

ð75Þ

subject to the same initial conditions (50), simultaneously with the modified Liouville equation

X pffiffiffiffiffiffi oq oq o Dii qfi þ þ ot oV i oV i i

! ¼ 0:

ð76Þ

Actually, since Eq. (76) does not depend upon Eq. (75), the latter can be solved in advance, so that the Liouville feedback becomes known. Eq. (76) is to be solved subject to the normalization constraint and an initial density that can be established based upon the accuracy of the computations. Obviously if stability of the system (49) in known in advance, the modification (75) and (76) is not necessary, although application of Eqs. (75) and (76) will only reaffirm this stability. However, in general, there is no analytical criterion that would predict chaos based only upon the system (49) without actual numerical runs. In view of that, the proposed approach seems quite universal. 3.6. Application to turbulence In the case of turbulence, prior to application of the proposed methodology, the Navier–Stokes equations must be approximated by a system of ODE. Such an approximation can be performed using finite differences, finite elements, or the Galerkin method. The relevance of the finite-dimensional approximation to solutions of the fluid dynamics has been successfully demonstrated by Lorenz [4], who applied the Galerkin method to the Rayleigh–Benard convection model keeping only three Fourier components; as a result, he arrived at a strange attractor that now bears his name. Despite a sharp truncation of the Galerkin expansion, there are obvious qualitative similarities between the original PDE model and its tree-dimensional approximation. However, it would be unreasonable to expect a quantitative similarity since both the models are unstable. At the same time, a finite-dimensional approximation for laminar flows exhibit a good quantitative approximation, and that justifies application of the proposed approach to modeling turbulence. 4. Conclusion Thus, a non-linear version of the Liouville equation is proposed for describing postinstability motions of dynamical systems with exponential divergence trajectories such as those leading to chaos and turbulence. The approach is based upon introduction of stabilizing forces that couple equations of motion and the evolution of the probability density of errors in initial conditions. These stabilizing forces create a powerful terminal attractor in the probability space that corresponds to occurrence of the target trajectory with the probability one. In configuration space, this effect suppresses exponential divergence of the close neighboring trajectories without affecting the target trajectory. As a result, the postinstability motion is represented by a set of functions describing the evolution of the statistical invariants such as expectations and higher moments, while this representation is stable. General analytical proof has been introduced. Since the proposed approach is not restricted by any special assumptions about the original dynamical system, it can be applied to both conservative and dissipative systems. The main applications for conservative systems are in celestial mechanics (for instance, many-body problems). The broad class of dissipative systems to which the proposed approach can be applied includes chaotic attractors and turbulence. It should be noticed that the proposed approach combines several departures from the classical methods. Firstly, it introduces a non-linear version of the Liouville equation that is coupled with the equation of motion (in Newtonian dynamics they are uncoupled). Secondly, it introduces terminal attractors characterized by violation of the Lipschitz conditions (in Newtonian dynamics as well as in the theory of differential equations these conditions are preserved). Finally, the idea of

M. Zak / Chaos, Solitons and Fractals 41 (2009) 1136–1149

1149

a forced stabilization of unstable equations follows from the Stabilization Principle formulated by Zak [13]. This is the most fundamental conceptual departure from the classical approach to mathematical modeling. Acknowledgements Copyright 2008 California Institute of Technology. Government sponsorship acknowledged. The research described in this paper was performed at Jet Propulsion Laboratory California Institute of Technology under contract with National Aeronautics and Space Administration. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

Arnold V. Mathematical methods of classical mechanics. New York: Springer; 1988. Landau LD, Lifshitz EM. Quantum mechanics. Oxford: Heineman; 1989. Lessen M, Paillet F. Phys Fluids 1976;19(7). Lorenz EN. J Atmos Sci 1963;20:130. Nottale L. Scale relativity: a fractal matrix for organization in nature. EJTP 2007;4(16(II)):187–274. Reynolds O. On the dynamical theory of incompressible viscous fluid and the determination of the criterion. Phil Trans R Soc Lond, Ser A 1895;186:123–64. Risken H. The Fokker–Planck equation. NY: Springer; 1989. Synge J. On the geometry of dynamics. Phil Trans R Soc Lond, Ser A 1926;226:31–106. Zak M. Uniqueness and stability of the solution of the small perturbation problem of a flexible filament with a free end. J Appl Math Mech Moscow 1970;34(6):1048–52 [Transl. Engl. pp. 988–992]. Zak M. Closure in turbulence theory using stabilization principle. Phys Lett A 1986;118:139–43. Zak M. Terminal attractors for associative memory in neural networks. Phys Lett A 1989;133(1–2):18–22. Zak M. Terminal model of Newtonian dynamics. Int J Theor Phys 1992;32:159–90. Zak M. Postinstability models in dynamics. Int J Theor Phys 1994;33(11):2215–80. Zak M, Zbilut J, Meyers R. From instability to intelligence. NY: Springer; 1997. Zak M. Self-supervised dynamical systems. Chaos, Solitons, & Fractals 2004;19:645–66. Zak M. From reversible thermodynamics to life. Chaos, Solitons, & Fractals 2005:1019–33. Zak M. From quantum entanglement to mirror neuron. Chaos, Solitons, & Fractals 2007;34:344–59. Zak M. Physics of life from first principles. EJTP 2007;4(16(II)):11–96.