2004-Fez conference on Differential Equations and Mechanics Electronic Journal of Differential Equations, Conference 11, 2004, pp. 95–102. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

GREEN’S FUNCTIONS AND CLOSING IN PRESSURE IN PARTIALLY NONHOMOGENEOUS TURBULENCE

MOHAMMED OUÇAMAH CHERKAOUI MALKI, OMARI ALAOUI ELKEBIR

Abstract. In this work, we are interested with the numerical solution of the equations of the correlations -or moment of order two - associated with the Navier-Stokes equations. The method of closing in pressure which is based on the elimination of the terms of pressure present in these equations by using the functions of Green, is completely re-examined. We underline the Green’s functions divergence problems which occurred with the traditional method of closing. Then we establish a new formalism which makes it possible to circumvent these problems. We present and confront in the course of our presentation two methods of construction of Green’s functions according to the choice of the boundary conditions, namely the method of images and the spectral method.

1. Introduction

When one is interested in non homogeneous turbulence, the Navier-Stokes equations and the models based on a description in only one point such as turbulent viscosity [1], “the k − ε” [8] or other models, provide only one partial description of all of the phenomena associated with these flows. One thus expects the development of new writings in two(or more) points. (i.e. descriptions taking into account the interaction between the various structures of turbulence)

\frac{∂}{∂t}V_{i} + V_{j} \frac{∂}{∂x_{j}} V_{i} − ηΔV_{i} + \frac{∂}{∂x_{i}} P = 0  (1.1)

where V_{i} is the velocity component along x_{i}, P is the pressure, and Δ is the Laplacian The models in two points are based on the resolution of the equations of the correlations [2], Q_{ij}(x, x′), given by

Q_{ij}(x, x′) = ⟨v_{i}(x)v_{j}(x′)⟩  (1.2)

where ⟨·⟩ stands for the statistical average and where v_{i} (i=1,2,3) are the velocity fluctuations, which are defined by V_{i} = ⟨V_{i}⟩ + v_{i} one must then solve the tensor of

2000 Mathematics Subject Classification. 65Z05, 65R10, 65R32, 76F55. Key words and phrases. Navier-Stokes; nonlinearity; Green’s function; quasi spectral analysis; mathematical formalism; data construction. \begin{matrix}c \\ ◯ \end{matrix}2004 Texas State University - San Marcos. Published October 15, 2004.

Reynolds in which each term is described by a nonlinear equation which contains other terms of the tensor.

(\begin{matrix}Q_{11} & Q_{12} & Q_{13} \\ Q_{21} & Q_{22} & Q_{23} \\ Q_{31} & Q_{32} & Q_{33} \end{matrix})

Obviously, development of such models is not easy and the numerical resolution of the new equations obtained starting from the equations of Navier-stokes is even more complicated; we must deal with problems in terms of mathematical formalism, numerical processing, and physical modelling [3, 4, 5, 7]. Besides, The complexity of the double correlation’s tensor and the non linearity make very difficult the control of the evolution of any inaccuracy allowed at the beginning of calculations or at the representation of pressure as a function of velocities.

2. Conventional closing in pressure of the Navier-Stokes equation

To simplify the search for a model of turbulence and the corresponding algorithm, it is adequate from our point of view to restrict the equations to the velocity terms. We thus understand by “closing in pressure” the representation of the terms of pressure present in the Navier-Stokes equations as functions of the velocity components.

The general approach of closing consists in introducing the Green’s kernel during a formal calculation and to use the equations relative to the incompressible fluids: ∑ \frac{∂}{∂xi}V i = 0.

In general, Green’s functions are presented in the form of series of functions; their introduction systematically generates truncation errors in the algorithm of the numerical resolution. Considering the complexity of the tensorial writing of the double correlations and the non linearity, the control of the evolution of these errors is then very difficult.

The conventional procedure consists of a first step, in writing the pressure in an integral form such as

P (x) =∫_{Ω} G(x, x”)ΔP (x”)dx”+∫_{Γ} \frac{∂}{∂n}G(x, x”)P (x”)dx”−∫_{Γ} G(x, x”)\frac{∂}{∂n}P (x”)dx”

(2.1)

where Ω is the flow domain, Γ its border and \frac{∂}{∂n} the normal derivative to Γ. G(x, x′) is the Green’s function which on Ω satisfies

Δ_{x}G(x, x′) = δ(x − x′) .  (2.2)

Note that here no more boundary conditions are imposed for G(x, x′) on Γ. One obtains (2.1) by two formal applications of Gauss theorem from:

P (x) = ∫_{Ω} ΔG(x, x”)P (x”)dx”  (2.3)

which is the starting point of all the attempts of closing of the conventional procedure.

The aim is now to determine ΔP in the domain P (x”), and \frac{∂}{∂n}P (x”) on Γ as functions of the velocities; see equation (2.1).

- In general ΔP is obtained from the Poisson’s equation (2.4); in the case of incompressible fluids, the writing of the Laplacian of pressure in terms of velocities is simplified and takes the form:

ΔP (x) = f(∂V i/∂xk)i, k = \frac{∂}{∂xi}(V j \frac{∂}{∂xj} V i)  (2.4)

- Moreover\frac{∂P}{∂n} is deduced by making degenerate the Navier-Stokes equations at Γ which only needs to know the geometry of the flow.

- P (x): for the calculation of the remaining integral at Γ, one introduces the following two limit problems: Dirichlet type problem:

\begin{matrix}Δ_{x}G(x, x′) = δ(x − x′) in Ω , \\ G(X, X′) = 0 on Γ ; \end{matrix}  (2.5)

Neumann type problem:

\begin{matrix}Δ_{x}G(x, x′) = δ(x − x′) in Ω , \\ \frac{∂G}{∂n}(x, x′) = 0 on Γ \end{matrix}  (2.6)

Now the assumptions can be taken on the Green’s function at Γ. In [6] it has been established a judicious procedure which consists in choosing boundary conditions on the Green’s kernel in such a manner to eliminate the quantities which one does not manage to translate. So, he imposes on the Green’s kernel boundary conditions of the Neumann type (2.6) of such kind to eliminate the integral on Γ containing the terms in P (x) in (2.1). One obtains then:

P (x) = ∫_{Ω} G(x, x”)f(\frac{∂V i}{∂x_{k}})_{i,k}dx” − ∫_{Γ} G(x, x”)\frac{∂P (x”)}{dn}dx  (2.7)

Thus, this simplified equation is of great interest and does not introduce any additional assumption on the pressure.

However, the problem which remains is the convergence of G(x, x′) given by (2.6) and which we are dealing with in the next section.

3. Green’s functions and their convergence

We present here two techniques of construction of the Green’s functions. We will try by the same occasion to clear up the origin of the divergence of some of these functions.

Method of images. The method of images has the advantage of providing an analytical solution. This method takes as a starting point that known in electrostatics and more precisely the method of coverings. It consists in constructing a charge distribution which would produce the equipotential fixed beforehand.

More generally, when in electrostatics one imposes equipotentials and that one wants to know the potential in an other given point of space, one adopts the process known as the process of covering. This process consists in determining the charge distribution, which would produce these equipotentials. Then, one identifies the desired potential with that which would produce the charge distribution thus determined.

To understand the concept of charge and image let us point out the following solutions of the problems of Dirichlet and Neumann for the common cases:

The entire ℝ^{3} space. When the domain represents the entire ℝ^{3} space, the two systems (2.5) and (2.6) have the same solution G∞ which is written as:

G_{∞}(x, x′) = \frac{−1}{4π∥x − x′∥}  (3.1)

This function is identified (except for a constant) with the potential associated with a charge q = −1 placed in x. One will thus speak about “turbulence charge” by associating the potential G∞ to it (In electrostatics, one associate a potential to each charge).

Semi infinite domain. When the domain is semi infinite (or semi closed) the Green’s functions are respectively:

G_{1}(x, x′) = \frac{−1}{4π∥x − x′∥} + \frac{1}{4π∥x^{∗} − x′∥}  (3.2)

For system (2.5), and

G_{2}(x, x′) = \frac{−1}{4π∥x − x′∥} + \frac{−1}{4π∥x^{∗} − x′∥}  (3.3)

for system (2.6). where x∗ is the symmetric of x with respect to the wall.

We remark that G_{1} can be obtained in placing a charge q′ = −q in a symmetrical point x∗ of x with respect to the plane(x_{2} = 0). In the same way, if we changes the sign of q’ then we find the analytical solution of (2.6):

The determination of the Green’s function is thus equivalent to that of a “turbulence charges” distribution whose equipotential zero coincides with the frontier of the domain. One will retain that in the first case, one has an alternation of sign whereas, in the second, one cumulates quantities of the same sign.

Case of the plane channel. Let us now consider the plane channel represented below:

[figure: ^{0}q^{1}; Γ_{2}; (Ω); x_{2}; x_{1}; x_{3}; ^{0}q; ^{0}q^{0}; Γ_{1}]

Figure 1. Case of plane channel

According to what we have just observed in the case of half spaces and, by considering initially only the plane Π_{1}, we start by placing a symmetrical charge q′ of q with respect to Π_{1} in order to satisfy the condition at the edge of this plane. We will place according to this same point of view another charge q” symmetrical of q with respect to Π_{2} satisfying the condition at the edge of this second plane. Obviously, these two new charges produce secondary effects which we must neutralize thereafter. These effects are due to the influence of q”on the Π_{1} plane and that of q′ onΠ_{2}.

We then place two new charges, the first one symmetrical of q” with respect to Π_{1} and the second one symmetrical of q′ with respect to Π_{2} and so on. Gradually, we set up an (infinite) charge distribution which satisfies at the same time the two conditions.

It is then verified easily that the contribution to the potential of each charge at the frontiers is neutralized respectively on each wall by that of the two charges which are respectively symmetrical for them.

The solution that we then obtain for the problem of Neumann is an ensemble of contributions resulting from charges of the same sign placed at distances in “1/r”:

G_{N} = \frac{−1}{4π} ∑_{n=−∞}^{n=+∞} \frac{1}{∥x − x′ − 4ane_{2}∥} + \frac{1}{∥x − x′ − 2a(2n + 1)e_{2} + 2x2^{′}e_{2}∥} (3.4)

and that relative to the problem of Dirichlet is presented in the form of alternate series (resulting from charges of opposite signs):

G_{D} = \frac{−1}{4π} ∑_{n=−∞}^{n=+∞} \frac{(−1)^{n}}{∥x − x′ − 4ane_{2}∥} + \frac{(−1)^{n}}{∥x − x′ − 2a(2n + 1)e_{2} + 2x2^{′}e_{2}∥} (3.5)

Note that the function G_{N} diverges for any couple of points (x, x′) of Ω × Ω and that the functionG_{D} is convergent in any point of Ω × Ω. Thus when one uses the method of images, the study of the coherence of the final solution obtained constitutes the final step. So nothing prevents from pushing a little further the analogy and of saying that, as in electrostatics, where one cannot conceive infinite potential in a point, one cannot preserve, as for as we are concerned , only the Green’s functions which take finished values inside the domain.

More complex geometry. When the geometry of the field is more complex the method of the images can still be applied. It is enough for that to realize that the process that we applied for two planes can be done for a polyhedric field delimited by N planes. Moreover, in the case of curvilinear frontiers, we can always approach the desired Green’s function by that associated to a close polyhedric geometry. The problem which remains is the divergence of the functions associated to Neumann’s system: the solution is provided by cumulating contributions of the same sign.

Quasi spectral method. This method which serves here to validate the method of images consists in transporting the differential systems (2.5) and (2.6) in a spectral space by a Fourier transform (on one or more space variables) and to solve the new system thus obtained.

It is hoped whereas that the solution of the spectral problem admits an opposite Fourier transform and that this transform is also solution of the initial problem. In the case of the entire space or of the half space, one can see easily that this process leads to the same solutions referred to above.

With regard to the plane channel and taking into account the geometry of the domain, it is necessary to take the spectral transform:

f̂(k_{1}, x_{2}, k_{3}) = ∫ ∫ f(x_{1}, x_{2}, x_{3}) exp(−i(k_{1}x_{1} + k_{3}x_{3}))dx_{1}dx_{3}  (3.6)

The differential system associated with the boundary conditions of the Dirichlet type can be written

\begin{matrix}(\frac{∂^{2}}{∂x_{2}^{2}} − q^{2})W^{13}(x_{2}, x”_{2}) = \frac{1}{4π^{2}} δ(x_{2} − x”_{2}), \\ W^{13} = 0 in Γ × Γ \end{matrix}  (3.7)

where q^{2} = k_{1}^{2} + k_{3}^{2} and

W^{13}(x_{2}, x”_{2}) = 4π^{2}Ĝ(k_{1}, x_{2}, x”_{2}, k_{3}) exp(i(k_{1}x_{1} + k_{3}x_{3}))  (3.8)

The solution of this differential system can be obtained by taking into account the following remarks:

(1) On any interval not containing x′_{2}, W is the solution of

\begin{matrix}(\frac{∂^{2}}{∂x_{2}} − q^{2})W^{13}(x_{2}, x′_{2}) = 0, \\ W^{13}(x_{2}, x′_{2}) = 0 for x_{2} = a and x_{2} = −a \end{matrix}  (3.9)

(2) W is continuous on [−a; a]. Indeed, a discontinuity in x_{2} = x′_{2} would lead to δ (undesirable in the second derivative).

(3) The jump: w^{13}(x′_{2}+) − w^{13}(x′_{2}−) = 1 (Theorem of derivation of the distributions represented by discontinuous functions).

The procedure then consists in finding W +, solution of (3.9) on the right of x′_{2}, and W, solution of (3.9) on the left of x′_{2}, and to connect them using conditions 2 and 3.

One finds that

W_{+}^{13}(x_{2}, x”_{2}) = \frac{C_{1}(x”_{2})}{q}sh[q(x_{2} − a)] for x_{2} > x”_{2}

W_{−}^{13}(x_{2}, x”_{2}) = \frac{C_{2}(x”_{2})}{q}sh[q(x_{2} + a)] for x_{2} < x”_{2}

with

C_{1}(x”_{2}) = −\frac{sh[q(x”_{2} + a)]}{sh(2qa)}, C_{2}(x”_{2}) = −\frac{sh[q(x”_{2} − a)]}{sh(2qa)}

One can then transform this writing in the following condensed form:

W^{13}(x_{2}, x”_{2}) = −\frac{e^{−q|x_{2}−x”_{2}|}}{2q} + \frac{ch[q(x_{2} + x”_{2})] − e^{−2qa}ch[(q(x_{2} − x”_{2})]}{2qsh(2qa)}

In the same way, the solution of the problem with condition of the Neumann type is written:

Z^{13}(x_{2}, x”_{2})) = −\frac{e^{−q|x_{2}−x”_{2}|}}{2q} + \frac{−ch[q(x_{2} + x”_{2})] − e^{−2qa}ch[(q(x_{2} − x”_{2})]}{2qsh(2qa)}

The fractional terms of W and Z behave respectively in: - a continuous function inW in the vicinity of zero (sh(q) is equivalent to q in the vicinity of zero) - a function in [1/q2] in the vicinity of zero in the expression of Z.

We remark that contrary to W (the Dirichlet solution), the function Z (Neumann solution) does not admit a transform of opposite Fourier in the physical space.

In addition, the Green’s function which we obtained in the physical space by the method of images can be easily obtained from the functions Z and W. It is enough, for that, to use the developments into exponential functions sh and ch, and to apply then the inverse Fourier transforms.

Conclusion. When constructing the solution relative to the plane channel, the recourse to infinite charge distributions, caused by the principle of compensation, is due to the presence of two parallel plans located on both sides of the charge q. In other words it is the confinement of the charge (of the flow) in its domain which gives solutions in the form of infinite series of functions. As the contribution of each charge to the potential is presented in the form of a general term of a series in 1/n, the convergence of this series of function cannot be carried out only if one utilizes charges of contrary signs. Consequently, the choice of boundary condition of the Neumann type must be ignored.

Let us notice that the method described above can be extended to more complex geometries even if it means to check the validity of the solution that it provides.

Table 1 gathers some results obtained for different geometries concerning the two classes of limit problems which interest us.

Table 1. Green’s functions nature for different geometries

Neumann Dirichlet
ℝ^{3} space  Convergent Convergent
Half plane  Convergent Convergent
Channel flow  Divergent Convergent
Polyhedric section Divergent Convergent

4. New closing

To summarize what proceeds we can say that the problem is to choose a model of closing which satisfies two fundamental conditions:

- The series of Green’s function should not be divergent.

- The resulting integral should not lead to introducing additional assumptions on the pressure.

And we can say that neither the condition of Dirichlet nor that of Neumann applied to the pressure are appropriate. The former necessitates additional physical conditions, the latter lead to a divergence of the Green’s function. Our governing idea consists in finding models of closing for each gradient components of pressure

With the help of an identical procedure to that which we provided for the pressure (see section 2.), we write each component of the gradient of pressure in an integral form and establish the relations:

\begin{aligned}\frac{∂P}{∂x_{i}} (x) & = ∫_{Ω} G(x, x”)Δ_{x} \frac{∂P}{∂x”_{i}} (x”)dx” + ∫_{Γ} \frac{∂}{∂n}G(x, x”)\frac{∂}{∂x”_{i}} P (x”)dx” \\ & − ∫_{Γ} G(x, x”)\frac{∂}{∂n}(\frac{∂}{∂x”_{i}} P (x”))dx” \end{aligned} (4.1)

starting from

\frac{∂P}{∂x_{i}} (x) = ∫_{Ω} ΔG(x, x”)\frac{∂}{∂x”_{i}} P (x”)dx”  (4.2)

Now, the choice of the limit conditions is then determined in a natural way on the basis of the constraints referred to above. The conditions of Dirichlet type (2.5) imposed on the limits of the domain on the Green’s function eliminate the second integral at G of the equation (4.1). And it remains:

\frac{∂P}{∂x_{i}} (x) = ∫_{Ω} G(x, x”)Δ_{x} \frac{∂P}{∂x”_{i}} (x”)dx” + ∫_{Γ} \frac{∂}{∂n}G(x, x”)\frac{∂}{∂x”_{i}} P (x”)dx” (4.3)

To calculate the remaining integral at Γ, we must know the gradients of pressure on it. However these quantities can directly be written starting from the Navier-Stokes equations as we explained previously in section 2.

In summary, we can say that this new formulation has a double advantage, On the one hand, concerning the determination of the values at the edge, we find the same physical considerations as in the judicious choice of the conventional procedure and, on the other hand, we utilize a new Green’s function which converges. Finally, formal calculation described by the new procedure of closing can be exploited numerically.

Conclusion. We have solved the problem of closing in pressure of the Navier -Stokes equations and we have provided a new mathematical formalism which extends the field of application of the models of closing using Green’s function. This formalism makes it possible to circumvent the problems of divergence of the Green’s functions and offers an alternative to the conventional formalisms.

References

[1] J. Boussinesq, Théorie de l’écoulement tourbillonnant et tumultueux des liquides dans les lits rectilignes à grande section, Gautier-Villard, Paris (1897).

[2] R. Comolet, Mécanique expérimentale des fluides Dynamique des fluides réels, Masson 1994.

[3] M. M. O. Cherkaoui, Contribution à l’étude de la turbulence non homogène par méthode intégrale et spectrale, Thèse de doctorat de l’Ecole Centrale de Lyon. (1992)

[4] Cherkaoui M. M. O., Omari Alaoui El Kébir, Hamadiche M.; Modélisation des équations de Navier-stokes pour un écoulement non homogène” 5ème, Congrès international de mécanique. Mohamedia (Avril 1999)

[5] S. F. Edwards and W. D. McComb, Local transport equation for turbulent shear flows, Proc. R. Soc. A. 330, pp 495-516 (1972)

[6] Hamadiche M., Analyse spectrale des mécanismes linéaires des écoulements turbulent s non homogènes, Thèse de doctorat d’état Université Claude Bernard Lyon I (1985)

[7] Kim J. , Moin P., & Moser R., Turbulence statistics in fully developed channel flow at low Reynolds number, Journal of Fluid Mechanic vol 177 pp 133 - 166. (1987)

[8] B. E. Launder, D. B. Spalding Mathematical models of turbulence, Academic Presse London (1972).

M. O. Cherkaoui Malki Laboratoire d’informatique, Département de Mathématiques et Informatique, Faculté des Sciences Dhar Mehraz, B.P. 1796 Fès-Atlas Morocco

   E-mail address: cherkaouimmo@hotmail.com

E. Omari Alaoui Laboratoire de Mécanique des fluides, Département de Physique, Faculté des Sciences Dhar Mehraz, B.P. 1796 Fès-Atlas Morocco

   E-mail address: cherkaouimmo@hotmail.com