Turbulent Jet Calibration

Isothermal turbulent round jet optimization Dragos Isvoranu, Mircea Marinescu, Dorin Stanciu Dept. of Thermodynamics University Politehnica of Bucharest 313, Spl. Independentei, sect.6 77206 - Bucharest, ROMANIA Tel: 40-1-410.04.00/339 E-mail: [email protected]

Abstract.The paper deals with optimization of k-eps turbulence model constants. The entropy generation rate for viscous isothermal round jet flow is averaged in the terms of the proposed turbulence model. Considering a homogenous turbulence flow it is obvious that the entropy production term will depend on the 5 model constants. We have proved through numerical simulation that these constants could be very well chosen on the basis of Prigogine theorem for steady states for which the entropy generation rate is minimum. Comparison with experiments as well as with the standard k-eps solution for the round jet is available. The proposed method could enable any researcher with a strong tool in order to make numerical turbulent steady flow simulation agree with experiments on a theoretic base.

Nomenclature σ = entropy source term (W/m3.K) T = temperature (K) Π = viscous stress tensor (kg/m2.s2) V = velocity vector (m/s) δ = thickness of layer (m) l = characteristic length (m) p = pressure (Pa) µ = dynamic viscosity (kg/m.s) U = axial velocity (m/s)

k = turbulent kinetic energy (J/kg) ε = turbulent dissipation rate (W/kg) r = jet radius (m) Re = Reynolds number ϕ = generic transport quantity υ = kinematic viscosity (m2/s) η ,ξ = non-dimensional coordinates κ = von Karman constant

K , E , F self-similar representations of k , ε and velocity

1. Introduction The following dissertation will take into account only the processes of steady homogenous turbulent flows. First step is to define the concept of steady state, which is that state far from thermodynamic equilibrium characterized by canceling the time derivatives of all state

parameters. From this point of view it becomes obvious that the steady regime engenders a finite source of entropy that will cancel only at thermodynamic equilibrium. Following Glansdorf and Prigogine demonstration [1], [2], the steady state could be associated with an extreme value of a certain state parameter or function and even more, this regime is stable against local perturbations. As they pointed out, for all steady states the source of entropy passes through a minimum, this principle being valid both for linear and non-linear phenomena. The simplest demonstration of this principle takes into account the case of thermodynamic systems in which the Onsager relations are valid and there is a linear correlation between generalized forces and fluxes, that means constant phenomenological coefficients. Based on this principle, we tried to figure out what might happen for a steady homogenous turbulent flow. Shouldn’t perhaps entropy generation rate be minimum? As it is well known any turbulence model needs a closure model which comprises more or less numerous empirical constants to comply with the necessity for simplicity and accuracy. Most of them are experimentally fitted, but anyway they can not cover all possible cases, not to mention those that are determined by numerical simulations. Taking for example the standard k − ε model there are five characteristic constants commonly involved for any computation. Wouldn’t these constants depend, for example, on the boundary conditions, as usually all steady states do? Or, to put in another words, as these constants have a tremendous influence on the final solution, as I personally found out, shouldn’t they comply to a certain principle intimately connected to reaching the steady state? Of course, this means that these constants are different from one case to another. The case to be analyzed, that is a free round isothermal jet, can easily be circumvented in the domain of linear irreversible thermodynamics. The equation for entropy source for isothermal nonreactive flows as commonly found in literature, Feidt [3], is: 1 σ = − Π : gradV (1) T where: T 2 (2) Π = µ L divV − µ L gradV + gradV 3 is the viscous pressure tensor and V is the velocity vector.

(

(

))

2. TSL approximation of the round jet For a free round isothermal turbulent jet, the system of equations in the thin shear layer hypothesis (TSL), is given in Bradshaw [4] and Wilcox [5]. This system of equations is the result of Reynolds average of turbulent quantities. The diffusion terms implying the laminar viscosity are usually neglected. As Bradshaw points out, the neglected terms especially in the x-direction are of order δ l times the retained ones, whereas in laminar flow the ratio is

(δ l )

2

. On the other hand, the relative thickness of turbulent layer is greater than in the laminar case, which restricts the use of TSL equations only to δ l << 1 situations. Note that pressure is no longer a variable but has been absorbed into the boundary conditions as ∂p ∂y is zero and ∂p ∂x could be replaced by the value at the surface or in the free stream where Bernoulli’s equation applies ( dp dx = − ρU e dU e dx ). All this introspection has the

only purpose to underline that TSL equations are not very accurate but for reason of mathematical and computational simplicity they represent a strong temptation. Therefore, we shall use them, always bearing in mind that possible errors could be generated by the intrinsic nature of the approximation hypothesized. Considering a self-similar representation of the solutions given by: J 0. 5 J J 1. 5 r U (x , r ) = F (η ) ; k (x, r ) = 2 K (η ) ; ε (x, r ) = 4 E (η ) ; η = (3) x x x x where J = mU / ρ is the specific impulse, N = C µ K 2 E and :

ϑ (η ) = −

1

η

η

F (η ′)η ′dη ′

(4)

0

An additional transformation has been used to improve numerical accuracy: dη d d dξ = or =N N dξ dη so that the equations assume the following form: dF 1 d dF ϑ η + NF 2 = dξ η dξ dξ 2

ϑ

dK 1 d dF 1 dK η + NFK + = σ k dξ dξ η dξ dξ

ϑ

E dF dE 1 d 1 dE η + 4 NFE + C1ε = σ ε dξ K dξ dξ η dξ

On another hand, it is easily seen that:

(5)

− Cµ K 2

(6)

2

− C µ C 2ε KE

µ 2 = ρJ 0.5 Re π

From the definition of the instantaneous strain rate tensor: 1 ∂v i ∂v j s ij = + 2 ∂x j ∂x i

(7)

(8)

in combination with equations (1) and (2), the instantaneous entropy generation rate for isothermal flow will look like: µ 2 2 (9) σ = L 2 s ij s ij − (v i ,i ) 3 T Operating a Reynolds average on equation (9), the mean entropy production will be:

σ =

µL

(2S S ) + υ ρ (2s′ s′ )

L (10) ij ij ij ij T T where S ij and s′ij means the Reynolds average and respectively the fluctuation of the strain

rate tensor component. Taking into account the turbulent dissipation definition: ε = υ L s ij′ s ij′

( )

(11)

like in Mohammadi [6], the nullity of velocity divergence for incompressible flows and the observation that the main contribution to the first term in (10), in the case of TSL

approximation, is given by (∂U ∂r ) , we will get the following expression for the global turbulent entropy production rate on unit length for the round jet: 2

Σ′ =

δ

2πρJ 1.5 2 1 dF ηN 2 x 4T 0 N dξ Re π

2

+ E dξ W/m.K

(12)

More details regarding this final equation can be found in the Ph.D thesis of Isvoranu [7].

3. Numerical approach In order to solve system (6) with the following boundary conditions: K ′(0) = E ′(0) = F ′(0) = 0

(13) K = E = F = 0 ;η → ∞ it is best to use a Crank-Nickolson implicit method as presented by Ferziger [8]. This type of discretization scheme is second order accurate in space and first order in time. For example, considering the 1D generic transport equation: ∂ϕ ∂ϕ Γ ∂ 2ϕ = −u + + Sϕ (14) ∂t ∂x ρ ∂x 2 results the following discretized relation: (15) APϕ in +1 + AEϕ in++11 + AWϕ in−+11 = S∆t (ψϕ in + (1 − ψ )ϕ in +1 ) − AEϕ in+1 − AW ϕ in−1 where superscript n is related to time step, subscript i to space step and the coefficients have the following expressions: Γ∆t Γ∆t ρu∆t ρu∆t ; AW = − − AE = − 2 2 4 ∆x 2(∆x ) 4∆x 2(∆x ) (16.a) ρ AP = ρ − ∆t ( AE + AW ) ; S = AE + AW + ∆t

Sϕ =

Sϕ in +1 for S < 0 Sϕ

n i

for S > 0

ψ =

0 for S < 0 1 for S > 0

(16.b)

u∆t is the Courant-Friedrichs-Lewy (CFL) number and Γ is the specific 4 ∆x diffusion coefficient. System (15) is solved by three diagonal matrix algorithm (TDMA), taking the CFL number equal to 1 and 101 spatial nodes of discretization. where λ =

4. Results and comparison with experiments Based on the numerical algorithm described above, a great number of turbulent round jet solutions have been computed, solutions depending on the set of constants C µ , Cε 1 , C ε 2 ,

σ k , σ ε and generating particular values of the entropy production rate. For a quite large range of constants’ values, the turbulent viscous dissipation has been plotted in order to observe the appearance of the alleged minimum value that would have lead us to choose the set of constants related to reaching the steady state. We have encountered several major difficulties in pursuing this purpose as not all sets of constants lead us to viable solutions for the round jet (TSL hypothesis was violated) or simply the system of equations had no

solutions forcing us to discard them. Having to deal with so many constants to optimize, we have undertaken this enterprise gradually, considering at first that σ ε may be determined from the log-layer equations (Wilcox [5]) as:

σ ε = κ 2 (Cε 2 − Cε 1 ) Cµ

(17)

where k is the Karman constant (0.4328 in the standard k − ε model) and that for C µ and

σ k we may retain the standard values of 0.09 and 1.0, respectively. In this case, the nondimensional entropy production rate: Σ′ =

TΣ′x 4 2πρJ 1.5

(18)

will depend only on Cε 1 , C ε 2 and Reynolds number. The distribution of Σ′ for the range C ε 1 ∈ [0.15,2.0] , C ε 2 ∈ [0.6,2.20] and various Reynolds numbers ranging from 15000 and 100000 shown a prominent valley close to the first bisector characterized by 30 to 40 local minimum entropy generation rate values. Figure 1 presents a glimpse of this domain for Re = 15000 , showing a relative minimum localized at Cε 1 = 1.31 and Cε 2 = 1.62 .

Figure 1

Table 1 illustrates the distribution of some of these relative minimum values of the entropy generation rate together with the corresponding pairs of turbulence model constants. Of course, the idea of steady state should be related to a global or absolute minimum for turbulent viscous dissipation, which, unfortunately, is not so obviously displayed by the solutions obtained. Table 1.

Local minimum index 9 13 15 17 19 22 25

Cε 1

Cε 2

Σ′min

0.59 0.76 0.85 0.94 1.03 1.17 1.31

1.06 1.19 1.26 1.33 1.40 1.51 1.62

3.5782 3.6049 3.6125 3.6197 3.6257 3.6233 3.6106

In order to discriminate between this set of minimum values the one value which represents the true steady state, we have checked if the choice made over the other three constants C µ ,

σ k , σ ε , has been properly done, that is if those values minimize the entropy production rate, also.

Figure 2

Hence, considering this time that the values of Cε 1 and C ε 2 are constant and equal with those resulted from the previous minimization, we have analyzed the entropy generation rate dependence on C µ , figure 2, σ k , figure 3 and on σ ε , figure 4. Figure 2 clearly shows that for all points selected (in a Cε 1 , C ε 2 coordinate system), the value 0.09 initially chosen for

C µ is the one that ensures minimum entropy production rate. In the case of the point defined by Cε 1 = 1.44 , Cε 2 = 1.92 , hence corresponding to the standard model, no minimum entropy generation resulted. A somehow different conclusion can be drawn analyzing figure 3. Turbulent viscous dissipation presents an obvious minimum value for all points selected but also, for points outside the valley presented in figure 1, like for example, the one corresponding to the standard model.

Figure 3

In both cases the value for σ k that minimizes the entropy production is 1.01, very close to the standard one. It is also true that in the last case the minimum is not so pronounced. Further analysis has been carried on, regarding the influence of σ ε on entropy production (figure 4). It became clear now that for all points along the specified valley, that is for certain pairs of Cε 1 and C ε 2 , the turbulent viscous dissipation attains an apparent minimum value, the first function drop, around 1.2 to 1.55, and a true minimum value, second drop, around 1.7 to 2.1. The standard value of σ ε is 1.3, which fits quite well with first apparent minimum. On the other hand, for the standard pair of Cε 1 and C ε 2 , the entropy production rate not only has no minimum value but it is a convex function. Next step in our optimization process is to compare the σ ε values that ensure minimum entropy production (figure 4) with those given by equation (17) that sets von Karman parameter κ to the standard theoretical value of 0.4328. We have to admit that upon this quantity there lies a slight controversy as most former experiments give 0.41 while recent ones (Bradshaw, [4]) stand for 0.44, more close to the theoretical value 0.4328.

Figure 4

Figure 5

Figure 5 illustrates clearly that the standard von Karman constant is obtained at the intersection of the two curves which happens only for local minimum number 25, that is for the optimized set of constants: C µ = 0.09 ; C ε 1 = 1.31 ; C ε 2 = 1.62 ; σ k = 1.01 ; σ ε = 2.03 (19)

The final test consists in the comparison between the experimental and optimized nondimensional axial velocity profile U U 0 ( U 0 is the centerline axial velocity) against the non-dimensional coordinate η . The result is shown in figure 6 along with the velocity profile resulted from the standard k − ε model constants, reflecting a better agreement of the optimized model with Bradbury's measurements (Wilcox, [5]).

Figure 6

The same conclusion can be drawn from the analysis of the spreading rate parameter δl δx , (the ratio r x for which U U 0 = 0.5 ), in the three cases mentioned before (table 2). Table 2.

δl δx

Optimized set 0,073

Standard set 0.122

Experiment 0.085-0.095

Relative error

14-23%

28.5-43%

-------

5. Conclusions A new approach of the turbulent incompressible round jet flow has been set forth in this paper. It has been clearly shown the methodology of determining the values of the five constants currently involved in the k − ε turbulence model on the basis of minimum entropy production characterizing steady states. The solution obtained, free from any

Reynolds number influence, is in better agreement with experiments than the one emerging from the standard model. The discrepancies that still remain may easily come from the limitations of the TSL approximation mentioned before. Although only a theoretical exercise, needing to be backed up on a comprehensive survey of other types of flows, we consider it fully expressive in what the Second Law of Thermodynamics may offer related to numerical simulation of flows and in the possibilities of connecting some physical constants to commonly accepted laws of Physics.

Acknowledgements I would like to thank prof. K. Hanjalic, head of the Department of Thermofluids at TU Delft, The Netherlands, for his pertinent remarks during the first stages in setting up this material.

References 1. Glansdorff, P., Prigogine, I. Physica, 20, 773, 1954. 2. Glansdorff, P., Prigogine, I. Thermodynamic theory of structure, stability and fluctuation, Wiley-Interscience, London, New York, Sydney, Toronto, 1976. 3. Feidt, M. Thermodinamique et Optimisation energetique des Systemes et Procedes Technique et Documetation Lavoisier, 1987. 4. Bradshaw, P. Introduction to Turbulence. Lecture Notes for short course. Burgers Certrum / TU Delft, 1996. 5. Wilcox, D.C. Turbulence Models for CFD. DCW Industries, Inc., La Canada, CA, 1993. 6. Mohammadi, B., Pironneau, O. Analysis of the K-Epsilon Turbulence Model. Ed. John Willey&Sons, 1994. 7. Isvoranu, D. Contributions to chemical reacting flows optimization. Ph.D. Thesis, Univ. Politehnica of Bucharest, 1999. 8. Ferziger, J.H., Peric, M. Computational Methods for Fluid Dynamics. Ed. Springer, 1995. 9. Bejan, A. Thermodynamics of an ‘isothermal’ flow: two dimensional turbulent jet’, International Journal of Heat and Mass Transfer, Vol. 34, 1991, pp. 407-413. 10. Bejan, A. Entropy Generation Minimization, CRC Press, Boca Raton, FL, 1996.