A novel approach to model turbulence Part I: theoretical formulation.

ABSTRACT

This paper describes a novel unifying formulation of turbulence modeling via the optimal control theory. In this formulation, Reynolds stresses are considered control variables while the averaged velocities are considered state variables. The Reynolds stresses are selected to optimize a performance index with the averaged Navier-Stokes (N-S) equations as constraints. The main problem of the new approach is the selection of the performance index to model diverse turbulence problems. The key feature in this approach is the selection of a performance index to represent the information about flow field and geometry implicitly. This stands in contrast to the classical turbulence modeling that proposes explicit models, which differ according to the flow problem and the geometry. The entropy concept is utilized to formulate the performance index. Entropy is related to the turbulence data, as described by probability density function of the velocity field. Entropy is taken as a measure of the information losses that result from averaging. Turbulence model is calculated to optimize the entropy and to satisfy the averaged N-S equations. The paper describes the mathematical formulation of turbulence modeling and the idea of information content of differential equations. Also, it outlines the new research problems related to this novel approach.

BACKGROUND

It is postulated that N-S equations, the geometry of the flow field, boundary and initial conditions describe the turbulent flow field completely. Turbulent flow field is computed numerically either by integrating the full N-S equations or by integrating an averaged N-S equations. The averaged N-S equations contain unknown terms called Reynolds stresses. Reynolds stresses are described in advance, in specific mathematical forms called turbulence models, to solve the closure problem of turbulence. Models range from an algebraic and differential equation models to the sub-grid models of Large Eddy Simulation (LES) technique (Brachet et al., 1986; Durbin and Reif, 2001; Dwoyer et al., 1985; Ferziger, 1987; George and Arndt, 1989; Godeferd et al., 2001; Hunt et al. 2001; Leonard and Hill, 1988). Selecting a turbulence model depends on the problem and differs from application to application. It is a heuristic procedure that combines physical insight, engineering sense and experience. There is no rigorous procedure to select a model for a given fluid flow problem.

The closure problem of turbulence is the result of averaging the convective nonlinear terms of N-S equations. Reynolds stresses appear as new source terms in the averaged equations. To integrate these equations, models of Reynolds stresses are selected empirically. The assumed terms add information to the equations from sources other than N-S equations. Successful models must be consistent with the information contents of N-S equations. This means that, to solve the averaged N-S equations properly the turbulence models must be consistent with information content of the full N-S equations. In an attempt to establish a systematic and unified modeling approach for turbulent model selection, the optimal control methodology is proposed. In this approach, a performance index must be specified in advance to complete the problem formulation. It is suggested that a performance index, which measures the information loss caused by averaging, is a good candidate.

When averaging the flow field equations, certain information disappears and other terms appear. The additional terms (Reynolds stresses) must somehow include part or all of the information that disappeared by averaging. To have a consistent model of turbulence, an information measure of N-S equations should be defined. This information measure may explain where information goes when N-S equations are averaged and how to retrieve it through proper selection of Reynolds stresses. The Probability Density Function (PDF) of the velocity field is suggested as the information measure as it is used in the entropy definition of information theory. The maximum entropy method is employed by many investigators of spectrum analysis (Childers, 1978) to generate more data points from a limited record of data in such a way that keeps the information content of the limited record of data unchanged. The novel approach presented in this paper outlines the possibility of modeling turbulence by utilizing the maximum entropy concept to justify the turbulence models.

The turbulence-modeling problem is formulated using the optimal control theory. The averaged N-S equations are considered the dynamical constraints, the Reynolds stresses are considered control variables, and the information measure is considered as the performance index. Since the averaging of N-S equations does not conserve information, a more general formulation of a measure of the information content of N-S equations is required. In this case, the information measure of N-S equations may represent the losses of information due to averaging. If the information measure is expressed as a functional of Reynolds stresses, the averaged velocity, and the gradient of the averaged velocity, then the turbulence model is computed to minimize information losses subject to the averaged N-S equation. Formulating turbulence modeling using optimal control techniques will help to unify the turbulence-modeling problem and suggests a scheme that can be used to understand the turbulent models used in literature. But the main issue here is the quantitative formulation of information measure.

TURBULENCE AND CLOSURE PROBLEM

It is postulated that for homogeneous and isotropic fluids (Newtonian fluid), N-S equations describe completely the turbulence. This means that given the geometry, the flow conditions, and the proper initial and boundary condition, the solution of the N-S equations gives the full picture of turbulence. Solving full N-S equations for 3D turbulence is a very difficult problem. This difficulty is the result of the lack of a general theory of nonlinear phenomena and the difficulty of resolving the different scales of turbulence computationally. For Eulerian formulation, let the spatial coordinates be [x.sub.i], i=1, 2, 3. Also, let the velocity components, the density and the pressure be [V.sub.i], [rho] and P respectively. For incompressible, unsteady 3D flow, the equations of fluid motion are given as follows:

Mass conservation:

[[partial derivative][V.sub.i]]/[[partial derivative][x.sub.i]] = 0 (1)

Momentum conservation:

[[[partial derivative][V.sub.i]]/[[partial derivative]t]] + [V.sub.j][[[partial derivative][V.sub.i]]/[[partial derivative][x.sub.J]]] + [1/[rho]][[[partial derivative]P]/[[partial derivative][x.sub.i]]] - v[[[[partial derivative].sup.2][V.sub.i]]/[[partial derivative][x.sub.j][partial derivative][x.sub.i]]] = 0,i = 1,2,3 (2)

Let the velocity components and the pressure be separated into mean ([U.sub.i] and PP) and fluctuating quantities ([u.sub.i] & p):

[V.sub.i] = [U.sub.i] + [u.sub.i]

P = PP + p

Where the mean quantity is defined as follows:

[U.sub.i] = [1/[[t.sub.2] - [t.sub.1]]][[t.sub.2].[integral].[t.sub.1]]w[V.sub.i]dt (3)

In the above equation, w is a weighting function. From the definition of averaging, the average of any fluctuating quantities is zero. The averaging time is long compared with the time scale of the turbulent motion. In transition problems, [t.sub.2]-[t.sub.1] has to be small compared with the time scale of the mean flow. Using the above decomposition, with w = 1, the averaged N-S equations can be written as follows:

[[partial derivative][U.sub.i]]/[[partial derivative][x.sub.1]] = 0 (4)

[[[partial derivative][U.sub.i]]/[[partial derivative]t]] + [U.sub.j][[[partial derivative][U.sub.i]]/[[partial derivative][x.sub.J]]] + [1/[rho]][[[partial derivative]P]/[[partial derivative][x.sub.i]]] - v[[[[partial derivative].sup.2][U.sub.i]]/[[partial derivative][x.sub.j][partial derivative][x.sub.i]]] + [[partial derivative]/[[partial derivative][x.sub.j]]]([u.sub.i][u.sub.j]) = 0,i = 1,2,3 (5)

Averaging nonlinear terms leads to the appearance of additional Reynolds stresses terms [u.sub.i][u.sub.j]. Since equations (1) & (2) contain more variables, more equations are required to describe the additional terms and to close the system of equations. Closing the problem by selecting Reynolds stresses is called turbulence modeling. Researchers use empirical, semi-empirical, or analytical formulae for different models. These models include algebraic models, one-equation models, one and half-models, k-[epsilon] or two-equation models, and LES sub-grid models (Mathiew, and Scott, J. 2000; Braza and Dussauge, 1998; Monin, and Yaglom, 1971; Orszag and Kells, 1980; Piquet, 1999; Pope, 2000; Rodi, 1980; Wilcox, 1994). Selection of the model depends on the experience, the flow conditions, and the geometry. Turbulence models add information from outside sources other than N-S equations. This addition, consistent or inconsistent with the full N-S equations, may explain the difficulties of selecting the proper model.

TURBULENCE MODELING VIA OPTIMAL CONTROLS THEORY

The variational approach will be proposed here as a tool to unify turbulence modeling problems. Define the mean values [U.sub.i] & PP as the state variables, and x, y, z & t as the independent variables. Let the control variables be [[gamma].sub.1] = [u.sub.1.sup.2], [[gamma].sub.2] = [u.sub.1][u.sub.2], [[gamma].sub.3] = [u.sub.1][u.sub.3], [[gamma].sub.4]=[u.sub.2.sup.2], [[gamma].sub.5]=[u.sub.2][u.sub.3] & [[gamma].sub.6]=[u.sub.3.sup.2]. Let [GAMMA] be a control vector with six elements, and U is the velocity vector with three elements:

[GAMMA] = [[[gamma].sub.1][[gamma].sub.2][[gamma].sub.3][[gamma].sub.4][[gamma].sub.5][[gamma].sub.6]][.sup.T]

U = [[U.sub.1][U.sub.2][U.sub.3]][.sup.T]

The turbulence modeling problem is formulated as follows. Find the six control variables ([[gamma].sub.1]-[[gamma].sub.6]) to optimize the functional:

I = [integral][integral][integral][integral][PHI] dx dy dz dt (6)

Subject to the constraints:

[[[partial derivative][U.sub.1]]/[[partial derivative]x]] + [[[partial derivative][U.sub.2]]/[[partial derivative]y]] + [[[partial derivative][U.sub.3]]/[[partial derivative]z]] = 0 (7)

-[[[partial derivative][U.sub.1]]/[[partial derivative]t]] = [U.sub.1][[[partial derivative][U.sub.1]]/[[partial derivative]x]] + [U.sub.2][[[partial derivative][U.sub.1]]/[[partial derivative]y]] + [U.sub.3][[[partial derivative][U.sub.1]]/[[partial derivative]z]] + [1/[rho]][[[partial derivative]PP]/[[partial derivative]x]] - v[[nabla].sup.2][U.sub.1] + [[partial derivative]/[[partial derivative]x]]([[gamma].sub.1]) + [[partial derivative]/[[partial derivative]y]]([[gamma].sub.2]) + [[partial derivative]/[[partial derivative]z]]([[gamma].sub.3]) (8)

-[[[partial derivative][U.sub.2]]/[[partial derivative]t]] = [U.sub.1][[[partial derivative][U.sub.2]]/[[partial derivative]x]] + [U.sub.2][[[partial derivative][U.sub.2]]/[[partial derivative]y]] + [U.sub.3][[[partial derivative][U.sub.2]]/[[partial derivative]z]] + [1/[rho]][[[partial derivative]PP]/[[partial derivative]y]] - v[[nabla].sup.2][U.sub.2] + [[partial derivative]/[[partial derivative]x]]([[gamma].sub.2]) + [[partial derivative]/[[partial derivative]y]]([[gamma].sub.4]) + [[partial derivative]/[[partial derivative]z]]([[gamma].sub.5]) (9)

-[[[partial derivative][U.sub.3]]/[[partial derivative]t]] = [U.sub.1][[[partial derivative][U.sub.3]]/[[partial derivative]x]] + [U.sub.2][[[partial derivative][U.sub.3]]/[[partial derivative]y]] + [U.sub.3][[[partial derivative][U.sub.3]]/[[partial derivative]z]] + [1/[rho]][[[partial derivative]PP]/[[partial derivative]z]] - v[[nabla].sup.2][U.sub.3] + [[partial derivative]/[[partial derivative]x]]([[gamma].sub.3]) + [[partial derivative]/[[partial derivative]y]]([[gamma].sub.5]) + [[partial derivative]/[[partial derivative]z]]([[gamma].sub.6]) (10)

The LaGrange multiplier method is used to formulate the augmented functional, and the maximum principle (cf. appendix) is employed to derive the equations for the Reynolds stresses. The resulting equations for the Reynolds stresses are generally nonlinear. To complete the modeling problem, a well-defined [PHI] of the performance index (6) must be specified.

The closure problem is mainly due to nonlinearity of the original equations. Hence, the nonlinear terms, convective terms of the N-S equations, will affect the turbulence model. Formulating the turbulence-modeling problem using the optimal control methodology reveals that the general mathematical structure of the model should be nonlinear. The general nonlinear structure of the turbulence model can be seen from the mathematical formulation of equations and constraints. This can be demonstrated by the following observations. The optimal control problem, which is described by a linear system of differential equations and quadratic performance index, leads to the nonlinear Riccati equation as solution for the control. The Reynolds stresses appear as a result of averaging the nonlinear convective terms of the N-S equations; therefore, the most suitable and general form of the model may be nonlinear to be consistent with the N-S equations.

A COMPUTATIONAL ALGORITHM OF THE TURBULENCE MODEL

In the direct methods, the solution of the averaged N-S equations is advanced in steps, from time [t.sub.k] to time [t.sub.k+1], with integration time step [DELTA]t(k) = [t.sub.k+1] - [t.sub.k]. At every time step, initial and boundary conditions (IC/BC) are supplied and a turbulence model [u.sub.i][u.sub.j] is required to integrate the system of equations (7)-(10). At each time step [DELTA]t(k), the following procedure is repeated from r=0 to a fixed value R:

1. Assume any form for [u.sub.i][u.sub.j] ([t.sub.k])[.sup.r] to advance the solution from [U.sub.i] ([t.sub.k+1])[.sup.r].

2. Calculate [I.sup.r] = I ([u.sub.i][u.sub.j] ([t.sub.k])[.sup.r], [U.sub.i] ([t.sub.k+1])[.sup.r])

3. Select another form for [u.sub.i][u.sub.j] ([t.sub.k])[.sup.r+1] to advance the solution from [U.sub.i] ([t.sub.k]) to [U.sub.i] ([t.sub.k+1])[.sup.r+1]

4. Calculate [I.sup.r+1] = I (([u.sub.i][u.sub.j] ([t.sub.k])[.sup.r+1], [U.sub.i] ([t.sub.k+1])[.sup.r+1])

5. Calculate [delta][I.sup.r] = [I.sup.r+1] - [I.sup.r]

6. Update [u.sub.i][u.sub.j] ([t.sub.k])[.sup.r+2] = (1 - [theta]) [u.sub.i][u.sub.j]([t.sub.k])[.sup.r] + [theta] [u.sub.i][u.sub.j] ([t.sub.k])[.sup.r+1] + [gamma] [S.sup.r], 0 [less than or equal to] [theta] [less than or equal to] 1

[S.sup.r], in step 6, is the direction of search for the optimum value of I, and [gamma] is the value of travel in S- direction to reduce [delta][I.sup.r] to zero.

Let

[DELTA][U.sub.i]([t.sub.k+1])[.sup.r] = [U.sub.i] ([t.sub.k+1])[.sup.r+1] - [U.sub.i] ([t.sub.k+1])[.sup.r]

[DELTA][u.sub.i][u.sub.j] ([t.sub.k])[.sup.r] = [u.sub.i][u.sub.j] ([t.sub.k+1])[.sup.r+1] - [u.sub.i][u.sub.j]([t.sub.k+1])[.sup.r]

Repeat steps 1 -6 for many values of r, the best solution is obtained when

Limit [DELTA][U.sub.i] ([t.sub.k+1])[.sup.r] [right arrow] o([epsilon]), r[right arrow] R

Limit [DELTA][u.sub.i][u.sub.j] ([t.sub.k])[.sup.r] [right arrow] o([epsilon]), r[right arrow] R

and

[delta][I.sup.r] [right arrow] o([epsilon])

The behavior of [DELTA][u.sub.i][u.sub.j]([t.sub.k])[.sup.r] can be investigated during the iterations using the following formula:

[DELTA][u.sub.i][u.sub.j] ([t.sub.k])[.sup.r] = [[[partial derivative][u.sub.i][u.sub.j]]/[[partial derivative][U.sub.p]]][DELTA][U.sub.p] + [[[partial derivative][u.sub.i][u.sub.j]]/[[partial derivative][nabla][U.sub.p]]][DELTA][nabla][U.sub.p] + o([epsilon]), p = 1,2,3

The initial selection of [u.sub.i][u.sub.j] ([t.sub.k])[.sup.0] can use a simple algebraic model. If a converged solution occurs, the same procedure is repeated for [t.sub.k+2], [t.sub.k+3],..., [t.sub.k+N]. A good starting guess for [u.sub.i][u.sub.j] ([t.sub.k+1])[.sup.0] for the next time step is [u.sub.i][u.sub.j] ([t.sub.k])[.sup.R]. In this dynamic modeling, a heuristic relation between [u.sub.i][u.sub.j] ([t.sub.k+1]) and [u.sub.i][u.sub.j] ([t.sub.k]) may speed up the convergence.

SEARCHING FOR A POSSIBLE [PHI]

Since there are many choices of [PHI], the main problem in the optimal control formulation is the best selection of [PHI]. Three selections of [PHI] are considered to demonstrate the application of the variational approach.

An Unconstrained Case Model

For unsteady, 3D incompressible flow the Reynolds stresses (uiuj) is determined in such a way to minimize the functional equation (6). To minimize the above integral expression, [PHI] must be specified in advance. The argument of [PHI] may contain the mean velocity, the mean velocity gradient, mean convected terms in addition to Reynolds stresses. From the literature it appears that there is no unique model for turbulence but there are different models depending on the application. The existence of a unique model depends of the definition of averaging and the existence and uniqueness solution of the N-S equations. Since the N-S equations are nonlinear there is no guarantee that there exists only one solution. On the other side it can be shown that for a large category of flow problems governed by incompressible assumption there is a unique solution, Ladyzhenskaya, 1969, Temam, 1981. In this case let us assume that there exists a unique solution of the N-S equations and hence an optimum turbulence model. Let the best selection of Reynolds stress be as follows:

[u.sub.i][u.sub.j] = [u.sub.i][u.sub.j.sup.opt]

If we substitute [u.sub.i][u.sub.j.sup.opt] into equation (5) the right hand side will be zero. For other selections the right hand side will not be zero. And the solution of (4) and (5) will not give the correct solution of the averaged values. Assume that for any selection other than [u.sub.i][u.sub.j.sup.opt], the residuals of the momentum and continuity equations will be [[epsilon].sub.x], [[epsilon].sub.y], [[epsilon].sub.z] and [epsilon]. Here, [[epsilon].sub.x], [[epsilon].sub.y], [[epsilon].sub.z] and [epsilon] are the right hand side of equations (4) and (5). The Reynolds stresses can be computed to minimize the total squared residual errors as expressed by [PHI] = [[epsilon].sub.x.sup.2] + [[epsilon].sub.y.sup.2] + [[epsilon].sub.z.sup.2] + [[epsilon].sup.2], i.e.

I = [integral][integral][integral][integral] ([i=3.summation over (i=1)] [[[partial derivative][U.sub.i]]/[[partial derivative][x.sub.i]]])[.sup.2] + [i=3.summation over (i=1)] ([[[partial derivative][U.sub.i]]/[[partial derivative]t]] + [U.sub.j][[[partial derivative][U.sub.i]]/[[partial derivative][x.sub.j]]] + [1/[rho]][[[partial derivative]PP]/[[partial derivative][x.sub.i]]] - v[[[[partial derivative].sup.2][U.sub.i]]/[[partial derivative][x.sub.j][partial derivative][x.sub.i]]] + [[partial derivative]/[[partial derivative][x.sub.j]]]([u.sub.i][u.sub.j]))[.sup.2] dx dy dz dt

The only implied assumption here is that the Reynolds stresses are expressed as functions of the mean values of the velocities, the gradient of the mean values, and the diffusion terms. This is in contrast to Boussinesq assumption of using the linear gradient of velocity. Letting [U.sub.i] and [u.sub.i][u.sub.j] be the unknown functions, the Euler-LaGrange equations will produce the following equations for the Reynolds stresses [[gamma].sub.p]:

[[[partial derivative][PHI]]/[[partial derivative][[gamma].sub.p]]] - [[partial derivative]/[[partial derivative]x]]([[partial derivative][PHI]]/[[partial derivative][[gamma].sub.xp]]) - [[partial derivative]/[[partial derivative]y]]([[partial derivative][PHI]]/[[partial derivative][[gamma].sub.yp]]) - [[partial derivative]/[[partial derivative]z]]([[partial derivative][PHI]]/[[partial derivative][[gamma].sub.zp]]) = 0, p = 1,2,3,4,5,6.

Also, the Euler-LaGrange equations for Ui are:

[[[partial derivative][PHI]]/[[partial derivative][U.sub.i]]] - [[partial derivative]/[[partial derivative]t]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.t]]) - [[partial derivative]/[[partial derivative]x]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.x]]) - [[partial derivative]/[[partial derivative]y]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.y]]) - [[partial derivative]/[[partial derivative]z]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.z]]) + [[[partial derivative].sup.2]/[[partial derivative][x.sup.2]]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.xx]]) + [[[partial derivative].sup.2]/[[partial derivative][y.sup.2]]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.yy]]) + [[[partial derivative].sup.2]/[[partial derivative][z.sup.2]]]([[partial derivative][PHI]]/[[partial derivative][Ui.sub.zz]]) = 0,i = 1,2,3.

Turbulence Generation and Dissipation Model

One of the goals of the present study is to find a quantitative formulation for [PHI]. Since no specific performance index exists at present, the turbulent energy production (P) and dissipation ([epsilon]) terms are considered for investigation. The turbulence budget is a good start for this search. Some results from the literature (Pope, 2000; Durbin and Reif, 2001; Mathiew and Scott, 2000) suggest that the dissipation and production distribution are similar in shape and have symmetrical distributions. Therefore, a suitable integrand of equation (6) should have the form:

[PHI] = f ([epsilon], P)

To account for the published results of turbulence budgets, a possible function can be written as a difference between two functions [chi] & [psi]:

[PHI] = f ([epsilon], P) = [chi] ([epsilon]) - [psi] (P)

A suggested simple form is

[PHI] = f ([epsilon], P) = [alpha] [[epsilon].sup.2] - [beta] [P.sup.2]

where [alpha] and [beta] are suitable positive constants. More studies are needed to find the suitable function f. At present, the author is conducting pilot studies to check this form of [PHI]. Results of the pilot study focus on investigating the behavior of the dissipation term ([epsilon]) and production term (P) for different turbulent models of different geometries.

Information Losses Model

In the optimal control formulation of turbulence modeling, a reasonable selection of a performance index 'I' is required. Hence, the focus is changed from selecting a proper mathematical representation of Reynolds stresses to selecting a proper [PHI]. In information losses model, 'I' is selected to represent information losses produced by averaging. Part of the losses appears as new terms 'Reynolds stresses.' Minimizing 'I' will produce appropriate Reynolds stresses that minimize information losses. Here [PHI] represents information losses of the averaged N-S equations. Then the problem of turbulence closure becomes the problem of the choice of [PHI] as information measure, which leads to the new concept of information content of differential equations. To use this model, a definition of the information measure of the N-S equations is necessary. In this case, the following questions must be answered: What is the information measure of the N-S equations? Does it exist? Can we find it?

Information Measure: Does It Exist?

A mathematical model represented by the equations, boundary conditions, initial conditions, and geometry describe a physical phenomenon. It is postulated that an information measure of the system of differential equations exists that conveys basic features of the physics. Since the information about the field conveys physical characteristics, the information measure must be invariant under mathematical transformation of equations describing the physics. Averaging is a special type of transformation of equations. If the averaging is an invariant transformation with respect to the information, then we should be able to select Reynolds stresses to satisfy the invariance properties. If the average is not an invariant operation, we can select the Reynolds stress to minimize the amount of information lost by averaging. In the first case, finding a measure of information content requires finding the invariant quantities of differential equations and proving that information does not change by transformation. From the mathematical definition of averaging, it appears that, for non-linear equations, part of the information is lost by averaging. Since averaging does not conserve information, another measure of information losses needs to be defined. Finding [PHI] depends on answering the following: What are the invariants of nonlinear differential equations? And if the invariants do not exist, then how to define the information losses?

Information Measure and the Invariants of Differential Equation

In turbulence modeling, the nonlinearity of the N-S equations necessitates looking for invariant properties of nonlinear differential equations. Since invariants of a nonlinear equation is an unsolved mathematical problem, extending concepts from linear analysis may shed some light on finding the information measure of nonlinear differential equations. Then, as a first trial to find the invariants, the concepts of the eigen value and the trace of a matrix, as invariants of linear systems, can be extended to nonlinear systems. This issue requires more theoretical research.

Averaging and Information Losses

To study the information losses mechanism, an investigation of the process of averaging is required. In addition to that, investigating the different ways of representing a function using different mean values is essential. Averaging operation includes selection of time interval [t.sub.1]-[t.sub.2] and weighting function (w in equation 3). It can be assumed that there is an optimum selection of the time interval and weighting function that lead to minimum information losses. On the other side, higher order averages or moments will represent more closely the real function and reduce the information losses. These averages may be included in the performance index that represents the information content of differential equations. Here again, since the deterministic set of the N-S equations becomes indeterminate by averaging, the additional Reynolds stresses terms can be seen as a source of information losses. In this respect, an invariant property may be expressed as a function of different mean values of the flow field variables. This expresses information as an additive property by using more moments.

Elements of Information Measure

The assumed information measure is composed mainly of three terms. One term represents the geometry; the second term represents the initial and boundary conditions; the third term represents some invariant properties of the N-S equations. Since Reynolds stresses appear as a result of the nonlinear terms, some mathematical form derived from the convective term must be included in this term.

Maximum Entropy Concept

There is no single correct technique to calculate [u.sub.i][u.sub.j] in absence of knowledge or information about the type of process that generates the turbulence. Information about a flow field and the probability density function of velocities are connected through the concept of entropy, Childers, 1978; Katz, 1967. This is because the probability of occurrence of an event is related to information about the event. Therefore, entropy, as a measure of the uncertainty described by a set of probabilities, will be taken as a guide to formulate the performance index in turbulence modeling. If [u.sub.i][u.sub.j] is perfectly determined, i.e. no uncertainty exists, the entropy is zero and in all other cases the entropy is positive. Hence, entropy is taken as a measure of uncertainty in selecting the Reynolds stresses as it appears in the averaged N-S equations. By this selection, minimization of the entropy increases the possibility of selecting the right [u.sub.i][u.sub.j]. A mathematical form of entropy that describes the turbulent field is required. Taking the definition of the entropy from information theory, Childers, 1978, we can demonstrate the similarities between the data representing the turbulence in fluids and the data representing the signal in information theory. The entropy as an information measure of a random signal is defined by the equation,

S = -[SIGMA][p.sub.i] log [p.sub.i] (11)

Entropy, as a measure of the uncertainty described by a set of probabilities pi, is a nonlinear function. Since it is hypothesized that the full N-S equations, geometry, initial and boundary conditions determine turbulence, then the entropy should be based on these equations as a source of generating flow field data. A heuristic definition of the S may be formulated if we consider the discrete form of the N-S equations as a source of generating velocity and pressure signals Vi and Pi with sampling intervals [DELTA]x, [DELTA]y, [DELTA]z for each [DELTA]t. Given the PDF of the flow field, S in the above equation represents an information measure of the N-S equations. Some information will be lost as a result of averaging the nonlinear convective terms. Minimization means that we select the [u.sub.i][u.sub.j] terms in such a way as to reduce the uncertainties or information not related to the N-S equations. That is, we minimize the losses in information due to averaging or due to inappropriate selection of Reynolds stresses. The use of the logarithm in equation (11) indicates that the information is an additive quantity. The entropy is zero for a deterministic system; i.e., all the probabilities [p.sub.i] are zero except, one which is unity. In that case, the system is perfectly determined and no uncertainty exists. This corresponds to the case of direct solution of the N-S equations. In all other cases, the entropy is positive. Entropy is then a measure of disorder in the averaged N-S system, and hence it suggests our ignorance about the mathematical structure of the turbulence stresses. In this sense, optimizing entropy is a procedure to select the turbulence model without adding any other information other than that given by the N-S equations.

TURBULENCE MODELING ALGORITHM USING INFORMATION THEORY

Maximum entropy concept (Childers, 1978) will be utilized to formulate the proper performance index. It is known that the probability of occurrence of an event ([p.sub.i]) is related to information. If P(V) is the probability density function of velocity distribution (Monin and Yaglom, 1971; Pope, 2000) then

1 = [integral] P(V) dV

[U.sub.i] = [integral][V.sub.i] P(V) dV

[U.sub.i][U.sub.j] + [u.sub.i][u.sub.j] = [integral][V.sub.i][V.sub.j] P(V) dV

The last equation can be written as

[u.sub.i][u.sub.j] = [integral][V.sub.i][V.sub.j] P(V) dV - [integral][V.sub.i] P(V) dV [integral][V.sub.j] P(V) dV

It can be shown that a formula of P(V) that satisfies the above equations is given by

P(V) = exp([i=3.summation over (i=1)][[alpha].sub.i][V.sub.i] + [i=3.summation over (i=1)][j=3.summation over (j=1)][[beta].sub.ij][V.sub.i][V.sub.j])/ZZ

ZZ = [integral][integral][integral] exp([i=3.summation over (i=1)][[alpha].sub.i][V.sub.i] + [i=3.summation over (i=1)][j=3.summation over (j=i)][[beta].sub.ij][V.sub.i][V.sub.j])d[V.sub.1]d[V.sub.2]d[V.sub.3]

where [[alpha].sub.i] and [[beta].sub.ij] are arbitrary functions. From the above equations, the averaged velocities and the Reynolds stresses are given by:

[U.sub.i] = [[partial derivative]/[[partial derivative][[alpha].sub.i]]]log(ZZ)

[u.sub.i][u.sub.j] = [[partial derivative]/[[partial derivative][[beta].sub.ij]]]log(ZZ) - [U.sub.i][U.sub.j], i [greater than or equal to] j

The turbulence modeling problem can be stated as finding the nine coefficients [[alpha].sub.i] and [[beta].sub.ij] at each point (x, y, z; t), to optimize the entropy S locally, where

S(P) = -[integral] P(V) log P(V) dV

subject to equations (7)-(10) as constraints. The function S is a function of nine coefficients.

S = S ([[alpha].sub.i], [[beta].sub.ij])

The function S will be used in an iterative procedure to compute the Reynolds stress at each point.

NEW AREAS OF RESEARCH

In addition to the unifying mathematical formulation, the new approach has demonstrative capabilities that can be used to explain multitudes of turbulence models existing in the literatures. Despite the fact that the optimal control formulation of turbulence modeling is rational, logical, and consistent with basics of mathematics and physics, it needs proof. Additional measures must be proposed to prove the variational approach of turbulence modeling. To investigate the validity of the above idea, the following studies are new research topics:

1) Since no specific entropy function exists at present to represent the information content of the N-S equations, the following plan is taken. A kernel of the integral is formulated to measure the loss of the value of the velocity and the loss of the value of the direction cosines of the velocity vector due to averaging at specific space-time point (x, y, z, t). Assuming the validity of Boussinesq assumption, a selected number of problems from the literature using this assumption should be solved by a CFD code. For each specific problem, different models are used and the assumed entropy functional is computed. The pilot studies will be focused on investigating the behavior of the entropy functional for different turbulent models of the same problem. The same study will be repeated for other problems.

2) Another study will use analogy to form [PHI] similar to the entropy and define the probability using the discrete velocities. In this formulation, turbulence is considered a random phenomenon and the entropy is taken as a measure of turbulence randomness.

3) A third study will focus on the probabilistic formulation as a formal mathematical approach and solve for the PDF using the maximum entropy formulation.

4) A fourth study will focus on discovering the ways by which the information disappears by averaging. Simple turbulence problems will be selected and numerical integration of the full and the averaged N-S equations will be employed to solve the problem. Repeating the procedure using different models will help answering the following question: where does the information go when we average the N-S equations, and how can it be recovered through the controlled selection of the model to minimize the losses of information during averaging process?

5) The ultimate mathematical proof focuses on the selection of [PHI] to produce the familiar algebraic and differential equation models. Another problem is the derivation of Boussinesq hypothesis using quadratic form of the Reynolds stresses.

6) A computational algorithm and CFD code to solve the turbulence modeling based on equations (6)-(10) will be developed.

CONCLUSION

This paper presents a general mathematical formulation of the problem of turbulence modeling using optimal control theory and entropy as a measure of information. It is postulated that turbulence can be described fully by the N-S equations with the proper initial and boundary conditions. Assuming that there is information measure of the N-S equations, and assuming that averaging is a transformation, then the information measure should not be changed by averaging N-S equations. Through averaging, certain amounts of information disappear and other sources of information, Reynolds stresses, appear. The Reynolds stresses must contain part or all of the information that disappeared by averaging. The key concept is the calculation of turbulent stresses to minimize the information losses resulted from averaging N-S equations. Most of the existing turbulence models could be viewed from the above formulation. The new concept of information measure may be used to select the turbulence model in such a way as to make the information measure of the averaged N-S equations as close as possible to the information measure of the original N-S equations. The main problem in this formulation is the selection of information measure of the N-S equations. In the next paper, a heuristic approach will be used to study different measures computationally.

APPENDIX

Maximum principle for distributed parameter systems 

Suppose a dynamical system is described by the partial vector differential equation, Schwarz 1971;

[[partial derivative]U(X,t)]/[[partial derivative]t] = f(X,U(X,t), [[[[partial derivative].sup.k]U(X,t)]/[[partial derivative][X.sup.k]]],[GAMMA](X,t)),

where X is the vector of independent spatial coordinates (x, y, z), t is the time, U(X,t)[member of][OMEGA] is the vector of n-dimensional state vector, [GAMMA](X, t) is r-dimensional control vector, and f is a vector-valued function of its arguments. The notation [U.sub.X]k = [[[partial derivative].sup.k]U(X,t)]/[[partial derivative][X.sup.k]] takes care of all possible partial derivatives where k is the highest partial derivatives of U w. r. t. X. The solution of the above equation needs specification of initial and boundary conditions:

U([t.sub.0],X)=[U.sub.0](X), and [[[partial derivative].sup.k-1]U(X,t)]/[[partial derivative][X.sup.k-1]], i=1,2,3,...,k-1

at [t.sub.0] and on the boundary [partial derivative][OMEGA].

We desire to find a control vector [GAMMA](X, t) from among the admissible controls that minimizes the performance index I for fixed initial and final times [t.sub.0] and [t.sub.f]:

I = [[t.sub.f].[integral].[t.sub.0]][[integram].[OMEGA]]F(U, X, [U.sub.X]k, [GAMMA](X,t))d[OMEGA]dt

The solution of this problem, if it exists, can be obtained from the following equations:

[[partial derivative]U]/[[partial derivative]t] = [[partial derivative]H]/[[partial derivative][lambda]] = f

[[partial derivative][lambda]]/[[partial derivative]t] = -[[[partial derivative]H]/[[partial derivative]V]] - [summation over k] (-1)[.sup.k][[[partial derivative].sup.k]/[[partial derivative][X.sup.k]]]([[partial derivative]H]/[[partial derivative][[[[partial derivative].sup.k]U]/[[partial derivative][X.sup.k]]]])

where the Hamiltonian H is defined as

H = F + [[lambda].sup.T](X, t)f

The optimal controls [GAMMA]* are determined by minimizing H w.r.t. choice of [GAMMA] (X, t) such that

H* [greater than or equal to] H,

where the asterisk denotes that the optimal admissible control vector [GAMMA](X, t) is used. Under certain condition this reduces to:

[[partial derivative]H]/[[partial derivative][GAMMA]] = 0

ACKNOWLEDGMENTS

This research was supported by NASA Stennis Space Center (SSC) grant NAG 13 - 03009.

LITERATURE CITED

Brachet, M. E., M. Meneguzzi, Politano, H. and Sulem, P. L. 1986. Computer Simulation of Decaying Two-Dimensional Turbulence, In G. Comte-Bellot and J. Mathieu [eds.], Advances in Turbulence, pp. 245-254. Springer-Verlag, New York, USA.

Braza, M. and Dussauge, J.-P. 1998. Computation and Comparison of Efficient Turbulent Models for Aeronautics, European Research Project ETMA. Verlag-Vieweg, Wiesbaden, Germany.

Childers, D. G. 1978, Modern Spectrum Analysis. IEEE Press, The Institute of Electrical and Electronics Engineers, Inc., New York, USA.

Durbin, P. A. and Reif, B. A. P. 2001. Statistical Theory and Modeling for Turbulent Flows. John Wiley & Sons, Ltd., USA.

Dwoyer, D. L., Hussaini, M. Y., and Voigt, R. G. 1985. Theoretical Approaches to Turbulence. Springer-Verlag, New York, USA.

Ferziger, J. H. 1987. Simulation of Incompressible Turbulent Flows: A Review. Journal of Computational Physics 69, 1-48, USA.

George, W. K., and Arndt, R. 1989. Advances in Turbulence, Hemisphere Publishing Corporation, New York, USA.

Godeferd, F. S., Cambon, C., and Scott, J. F. 2001. Two-Point Closures and their Applications: report on a workshop. JFM, Vol. 436, pp. 393-407, UK.

Hunt, J. C. R. et al. 2001. Development in Turbulence Research: a review based on the 1999 Programme of the Isaac Newton Institute, Cambridge. JFM, Vol. 436, pp.353-391. UK

Katz, A. 1967. Principles of Statistical Mechanics: The Information Theory Approach, W. H. Freeman and company, San Francisco and London, USA.

Ladyzhenskaya, O. A. 1969. Incompressible Navier-Stokes Equations. Pergamon Press, New York, USA.

Leonard, A. D. and Hill, J. C. 1988. Direct Numerical Simulation of Homogeneous Turbulent Reacting Flow. AIAA-88-3624 paper, USA.

Mathiew, J., and Scott, J. 2000. An Introduction to Turbulent Flow. Cambridge University Press, New York, USA.

Monin, A. S., and Yaglom, A. M. 1971. Statistical Fluid Mechanics, Vol. I. MIT Press, USA. Orszag, S. A., and Kells, L. C. 1980. Transition to Turbulence in Plane Poiseuille and Plane Couette Flow. JFM, Vol. 96, Part 1, pp. 159-205, UK.

Piquet, J. 1999. Turbulent Flows- Models and Physics. Springer-Verlag, New York, USA.

Pope, S. B. 2000. Turbulent Flows. Cambridge University Press, New York, USA.

Rodi, W. 1980. Turbulence Models and Their Application in Hydraulics- A State of the Art Review. Karlsruhe University, Germany.

Schwarz, H. 1971. Multivariable Technical Control Systems. North-Holland & American Elsevier publishing company, USA.

Temam, R., 1981. Navier-Stokes Equations. North Holland Company, Netherlands.

Wilcox, D. C. 1994. The Turbulence Modeling for CFD. DCW Industries, Inc., CA, USA, November 1994.

Abdel K. Mazher

Aerospace Science Engineering Department

Tuskegee University 

Tuskegee, AL 36088,