Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://hit-conf.imec.msu.ru/2006/CD/Palymskiy2.pdf Äàòà èçìåíåíèÿ: Wed Feb 22 16:10:20 2006 Äàòà èíäåêñèðîâàíèÿ: Fri Feb 28 20:07:13 2014 Êîäèðîâêà:

International Conference on Hydrodynamic Instability and Turbulence
Proceedings of ICHIT- 06 26 February ­ 5 March 2006, Moscow, Russia

LINEAR AND NONLINEAR ANALYSIS OF NUMERICAL METHOD FOR DNS OF TURBULENT CONVECTION

Igor B. Palymskiy Modern Academy for Humanities, Novosibirsk Branch , Novosibirsk, Russia, 630064 palymsky@hnet.ru

ABSTRACT
We study the spectral characteristics of the numerical method for DNS of turbulent convectional flows. These spectral characteristics compare with spectral characteristics of differential problem. On model nonlinear system equations also the nonlinear analysis of the numerical method has been performed. The results of calculation of turbulent convection with free boundary conditions are given at number Rayleigh up to 34000 times the critical value; these results show the transition from soft turbulence to hard. The derived numerical results are comparing with experimental data and numerical results of other authors.

INTRODUCTION
Many researchers studied thermal Rayleigh-Benard convection using numerical simulation. As rule, they used spectral methods with periodic boundary conditions. In numerical simulations stationary, periodic, quasiperiodic and stochastic regimes have been derived [1]. Some authors performed 2-D and 3-D simulations for high supercriticality with free [2,3] and rigid [4,5] boundary conditions on the horizontal planes. Unfortunately, in these works the spectral analysis of using numerical methods in linear and nonlinear approach is absent. In recent time it is shown that results of 2-D numerical simulations of turbulent convection with free boundary conditions on the horizontal planes at Prandtl number (Pr) is equal to 10 have a good agreement with experimental results on turbulent convection in gaseous He and air and results of numerical simulations with rigid boundary conditions at enough high supercriticality [6].

So far the full numerical simulation of 3-D turbulent convection is very complex problem demanding the large resources. The reasons are: 1. The existence simultaneously of rapidly increasing and rapidly decreasing harmonics in linear approach (at r = Ra/Racr = 1000 and Pr = 1 one of harmonics is increasing as e682·t). 2. The necessity of conformity of spectral characteristics of differential problem and numerical method, it is necessary for correct representation of infinitesimal disturbance development. 3. The necessity of calculations on enough big times of order of several the thermal diffusion times with enough big number of degrees of freedom. The mentioned troubles show that the important question is validation of numerical results. For validation of numerical results the linear analysis is necessary and nonlinear analysis is very desirable. The aim of this work is the linear and nonlinear (on model nonlinear system) analysis of spectral-finite difference method using by author for simulation of 2-D turbulent convection with free boundary conditions and the numerical investigation of transition from soft toward hard turbulence.

PROBLEM FORMULATION
Turbulent convective flow in a horizontal layer numerically is simulated at heating from below. The fluid is viscous and incompressible. The flow is

time-dependent and two-dimensional. Boundaries of a layer are isothermal and free from shearing stresses. The model Boussinesq is used without semiempirical relationships. The


International Conference on Hydrodynamic Instability and Turbulence dimensionless set of equations given in terms of deviations from an equilibrium solution is [1,7]:
x = x = Q = 0 with x = 0, L; 0 y 1 (on the vertical boundaries), L = /.

t +

1 ( y x - x y ) = + Ra Q x , Pr = - , 1 1 1 Q t + ( y Q x - x Q y ) = Q - x , Pr Pr Pr

NUMERICAL METHOD

(1)

where is a stream function, is the vortex, Q is the temperature deviation from equilibrium profile (the total temperature being T = 1 - y + Q), f = fxx +fyy is the Laplace operator, Ra = g3dQ/ is the Rayleigh number, Pr = / is the Prandtl number, g is the gravitational acceleration, , , are the coefficients of thermal expansion, kinematics viscosity and thermal conductivity, respectively, H is the layer height and dQ is the temperature difference on the horizontal boundaries. The required values , and Q are to be sought in the form:

We briefly describe our special spectral-difference numerical algorithm and testing [7]. Following a general ideology of the splitting method, transition from time layer n to time layer n+1 perform in two steps. On the first step, we take into account a linear progress of perturbations neglecting interaction between harmonics. Step 1.

t =

1 + Ra Q x , 2 = - , 1 1 Q - x , 2 Pr Pr

(3)

Qt =

(t , x, y ) = (t , x, y ) =
Q (t , x , y ) =


k ,m


S

k ,m

k cos( kx) sin( my ), k cos( kx) sin( my ),
sin( kx ) sin( my ), ( 2)


k ,m

k ,m

k ,m k ,m


k ,m

Q

where = /L is the wave number, and k = {0.5 (at k = 0, N) and 1 (at 1 k N-1)}, 0 k N, 1 m M-1, Skm = 2k2 + 2m2, here N and M are number of harmonics in x and y directions, respectively. High efficiency of spectral methods is coupled with representation of problem solution through aggregate of the eigenfunctions of the linear stability problem and if N and M are allowed to become infinite the representation is a complete orthogonal set. Also, such representation of the solution allows using the formulas of the linear stability theory; it significantly increases efficiency of numerical method. Solutions is periodic, but we consider this solution only in half of period in x ­ direction, therefore the periodic problem changes on the problem with boundary conditions on the side walls, according to form of solution. Thus, problem is solved in the region G = {(x, y)0 x L, 0 y 1} with the boundary conditions = = Q = 0 for y = 0, 1; 0 x L (on the horizontal boundaries) and

By substituting the solution representation (2) into system (3), we get the system of two ordinary differential equations for two unknown amplitudes km and Qkm in spectral space, which is solved analytically without any approximations on the time. The analytical formulas were derived by program of analytical calculations Maple V Release 4, it is near to the formulas of linear stability theory. Step 2.

t +

1 1 ( y x - x y ) = , Pr 2 = - , 1 1 1 Q - x , ( y Q x - x Q y ) = Pr 2 Pr Pr

(4)

Qt +

Here we use a finite-difference scheme of alternating directions for solving the system equations of nonlinear convective transfer in physical space, earlier used for simulation of turbulent convection [8]. For transition from spectral space into physical space and back, standard programs of FFT were used. The numerical method has the first order of time approximation and the second order of space approximation. The coefficients x and y in (4) are defined: · On the value of stream function on n-th time layer (Scheme 1),


International Conference on Hydrodynamic Instability and Turbulence
· · On the value of stream function after first splitting step (Scheme 2), On the value of arithmetic mean of stream function on n-th and n+1-th time layers (Scheme 3). The realization of Scheme 3 demands introducing of iteration process.

LINEAR ANALYSIS
Here we consider the linear analogs of differential problem (1) and numerical method. The solution (2) comprises only one harmonic and exponential depends from time, for instance (t , x , y ) = e - t c o s ( x ) s i n ( m y ), here is const, harmonics is increasing if <0. By closeness for differential problem (d) and numerical method (sr) we may estimate the precision of numerical method. The correct representation of such solution provides the correct representation of infinitesimal disturbances of equilibrium solution (trivial for system (1)). Fig.1 represents the spectral characteristics for m = 1,2 and 3, Ra = 1000·Ra, Pr = 1, N = 64 and M = 16, the time step is equal to 4·10-4. Here solid line is differential problem, symbol ­ numerical method of present work and dash line ­ finite difference numerical method [8], curves 1,2 and 3 are first, second and third modes (m = 1,2 and 3), respectively.

Figure 2 Instability boundary in spectral space We can see from fig.1 and fig.2 that suggested method has more precision than finite-difference and that both N and M must be proportional to Ra1/4 at correct simulation. Fig.3 shows the spectral characteristics of other spectral methods for most instable mode m = 1, the parameter values are same as fig.1, but time step is equal to 1.6·10-3. Here black line represents the difference problem, symbol ­ suggested numerical method, red line ­ Orszag method (changed for 2-D[9]), blue line ­ [10] and green line ­ [11].

Figure 3 Comparison of spectral characteristics It may be seen in fig.3 that suggested spectral-finite difference method preserved the high accuracy at increasing of time step in forty times. We may derive the analytical formulas [7]:

Figure 1 Spectral characteristics Fig.2 represents the instability boundary in spectral space, here alpha and beta are wave numbers, solid line is a boundary curve for differential problem ( = Ra1/41/2 in polar coordinates), dash line suggested spectral-finite difference method, dadot line ­ finite difference method [8], the parameter values are same as fig.1, but Ra = 10000·Racr now.

H H sr = d + ( + m ) - 1 4 - 2 4 m4 , 96 24 24
6 6 6



2

2

2

where d and sr correspond to differential problem and suggested numerical method, , H1 and H2 are time step and space steps in x and y directions.


International Conference on Hydrodynamic Instability and Turbulence
Unfortunately, the linear analysis doesn't allow to investigate the time approximation of nonlinear terms. Only tools of nonlinear analysis may make it.



n +1 d

-

n +1 sr

=-

A2 ( R a n - S n ) 2

f o r S c h e m e 1,

NONLINEAR ANALYSIS
We perform the nonlinear analysis of our numerical method on model nonlinear system of equations:



n +1 d n +1 d

- -

n +1 sr n +1 sr

A2 Ra n for Scheme 2 and 2 = O( 3 ) for Scheme 3. =

t + y x + x y = + RaQ x ,
= - , Q t + y Q x + x Q y = Q - x .
This nonlinear system has a family of private waveform solutions:

These formulas show that the amplitudes and calculate with same accuracy by Schemes 1-3 and that using of Schemes 1,2 lead to decreasing of calculation accuracy only in phase speed of solution. The test simulations showed that results of simulations are close for Schemes 1-3, therefore for DNS of turbulent convection using of Schemes 1,2 is expedient because of iteration missing.

(t , x , y ) = (t )e

i ( ( t ) + x + y )

, ,

SOFT AND HARD TURBULENCE
We have simulated Rayleigh-Benard convection with = 1, Pr = 10 and supercriticality r = Ra/Racr up to 34000. Fig.4 represents the Nusselt number versus supercriticality r at various resolutions and shows that resolution in 257*63 harmonics is sufficient for correct reproduction of flow development up to r = 34000. In our simulations we used 65*15 (N = 64, M = 16) harmonics at r 1000, 129*31 harmonics at 1000< r < 6000 and 257*63 harmonics at 6000 r. The calculated Nusselt numbers coincide with numerical result [10] at 5 r 40 with graphical accuracy, more detailed comparing is in [6].

Q ( t , x , y ) = - i ( t ) e

i ( ( t ) + x + y )

(t , x , y ) = (t ) =

(t , x , y )
S

,

1 - t - t (C 1e 1 + C 2 e 2 ), 2 1 - t (t ) = (C 1e 1 - C 2 e 2 SRa C C - t (t ) = 0 + A ( 1 e 1 + 2 e

- 2t

), ),

- 2t



1



2

i = S + ( - 1)

i

R a / S , i = 1, 2 .

Here A = 2/S, S = 2+2, C1 and C2 are arbitrary constants. The same solution has and suggested numerical method, thereafter these solutions may compare. In similar way, we investigated earlier finite difference numerical methods for nonlinear equation with oscillating viscosity and for simulation of viscoelastic flows [12]. Let the wave numbers and are small (long waves in x and y directions), let also the values of and on the n-th time layer are well know, we may derive after very cumbersome calculations with Maple program the expressions in form of power series for and on n+1-th time layer.



n +1 d n +1

- -

d

n +1 sr n +1 sr

= O ( 3 ), = O ( 3 ) for Sch em es 1 - 3,

Figure 4 Nusselt number at various resolutions The works of Chicago research group [13,14] show that at supercriticality r = Ra/Racr is equal to


International Conference on Hydrodynamic Instability and Turbulence
approximately 7000 in experiments on turbulent Rayleigh-Benard (R-B) convection the transition toward new mode of turbulent convection take place. The modes of turbulence were named as soft turbulence (at r < 7000) and hard turbulence (at r > 7000). After transition to hard turbulence, the form of probability density for temperature pulsations in the centre of cell changes, the gaussian distribution replaces by exponential, the Nu­r power law also was changed at transition. Fig.5 represents the probability density function for temperature pulsations in the centre of cell at r = 6000 (points), we may see that this distribution is gaussian. Here c = 0.068 is the value of rms temperature pulsations in the middle between the plates.

Figure 6 Temperature probability density at r = 10000

Figure 7 Alpha versus supercriticality Experimental result [13] was received for high values of supercriticality 1.2·105 r 1.2·107. Universality of probability density function (independence of from the supercriticality) was denoted also in experimental work [14]. The result of present simulation also shows approximately universality at 1.6·104 r 3.2·104 and close to both experimental and theoretical values. We note that in numerical simulations [5] (2-D, rigid, Pr = 7) and [4] (3-D, rigid, Pr = 0.7) all profiles of probability density have exponential form. Fig.8 represents the Nu­r power law at free boundary conditions. It is seen that at r 7000 the Nu­r power law was changed.

Figure 5 Temperature probability density at r = 6000 Fig.6 represents also the probability density function for temperature pulsations in the centre of cell (points), but at r = 10000, we may see that this distribution is exponential; here c is equal to 0.067. Red line of fig.6 represents the result of 3-D numerical simulation in air at r = 9800 [2], green line - reestimated experimental result of Chicago group [13] and blue line - theoretical result [15]. We may see the differences only in neighborhood of mean value of temperature pulsations. If the distributions the same shown in fig.6 is fitted by exponential distribution:

p (T ') =


2c

e x p ( - | T ' | / c ) ,

then we obtain results shown on fig. 7, when we show versus supercriticality r, here black points and black solid line ­ present result of 2-D simulation with free boundary conditions, symbol - numerical result [2], symbol ­ numerical result [5], red and blue lines ­ experimental [13] and theoretical [15] results, respectively.

Figure 7 Nu-r power laws


International Conference on Hydrodynamic Instability and Turbulence
CONCLUSION
The suggested numerical method exactly reproduce the spectral characteristics of differential problem, it guarantees the correct reproduction of the infinitesimal disturbance development. It is kept at increasing of time step. On model nonlinear system the nonlinear analysis of suggested numerical method was also performed. It is shown that calculation of nonlinear transfer coefficients on value of stream function with n-th time layer (Scheme1) and on value of stream function after first step of splitting (Scheme2) leads to decreasing of calculation accuracy only in phase speed of solution without dropping of the common accuracy of simulation. For DNS of turbulent convection using of Schemes1,2 is expedient because of missing of iterations. For instance, we represent the simulation results of turbulent convection, it is shown that 2-D simulation with free (stress-free) boundary conditions reflects transition from soft toward hard turbulence. The gaussian distribution was replaced by exponential for temperature pulsations in center of cell, the characteristics of exponential distribution are close to both experimental and theoretical data at 16000 r 32000. The Nu ­ r power law also was changed at transition. 6. Palymskiy, I. B., Direct Numerical Simulation of Turbulent Convection, in book: Progress in Computational Heat and Mass Transfer, R. Bennacer (editor), vol.1, pg. 101-106, Lavoisier, 2005, http://palymsky.narod.ru/Paris.htm 7. Palymskiy, I.B., 2000, Metod Chislennogo Modelirovaniya Konvektivnykh Techeniy, Vychislitel'nye texnologii, 5, pp.53-61; http://palymsky.narod.ru/ 8. Paskonov, V.M., Polezhaev, V.I. and Chudov, L.A., 1984, Chislennoe Modelirovanie Protzessov Teplo- I Massoobmena, Nauka, Moscow. 9. Goldhirsch, I., Pelz, R.,B. and Orszag, S.A., 1989, Numerical Simulation of Thermal Convection in a Two-Dimensional Finite Box, J. Fluid Mech., 199(1), pp.1-28. 10. Babenko, K.I. and Rachmanov, A.I., 1988, Chislennoe Issledovanie Dvumernoj Konvektzii, Preprint of Applied Mathematics Institute, 118, Moscow. 11. Rozhdestvenskiy, B.L. and Stojnov, M.I., 1987, Algoritmy Integrirovaniya Uravneniy Nav'e-Stoksa, Imejuschie Analogi Zakonam Sochraneniya Massy, Impul'sa I Energii, Preprint of Applied Mathematics Institute, 119, Moscow. 12. Palymskiy I.B., 1988, Chislennoe Modelirovanie Konvektivnykh Techenij pri Vysokikh Chislakh Vejssenberga, Preprint of Theoretical and Applied Mechanics Institute, 15, Novosibirsk. 13. Wu, X.-Z. and Libchaber, A., 1992, Scaling Relations in Thermal Turbulence: The Aspect-Ratio Dependence, Physical Review, A 45(2), pp.842-845. 14. Castaing, B., Gunaratne, G., Heslot, F., Kadanoff, L., Libchaber., A., Thomae, S., Wu, X.-Z., Zaleski, S. and Zanetti, G.,1989, Scaling of Hard Thermal Turbulence in Rayleigh-Benard Convection, J. Fluid Mech., 204(1), pp. 1-30. 15. Yakhot, V., 1989, Probability Distributions in High-Rayleigh-Number Benard Convection, Phys. Rev. Letters, 63(18), pp. 1965-1967.
Igor Palymskiy is Professor of Modern University for Humanities, Novosibirsk Branch, Mathematics Department. His main scientific interests are Direct Numerical Simulation of Turbulent Flows and Flows with Hydrodynamical Instabilities.

REFERENCES
1. Palymskiy, I. B., 2003, Determinism and Chaos in the Rayleigh-Benard Convection, Proceeding of the Second International Conference on Applied Mechanics and Materials (ICAMM 2003), Durban, South Africa, pp.139-144; http://palymsky.narod.ru/ 2. Sirovich, L., Balachandar S. and Maxey, M.R., 1989, Numerical Simulation of High Rayleigh Number Convection, J. Scientific Computing, 4(2), pp.219-236. 3. Malevsky, A.V. and Yuen, D.A., 1991, Characteristics Based Methods Applied to Infinite Prandtl Number Thermal Convection in the Hard Turbulent Regime, Phys.Rev., A 3(9), pp. 2105-2115. 4. Kerr, R.M., 1996, Rayleigh Number Scaling in Numerical Convection, J. Fluid Mech., 310, pp.139179. 5. Werne, J., DeLuca, E.E. and Rosner, R., 1990, Numerical Simulation of Soft and Hard Turbulence: Preliminary Results for Two-Dimensional Convection, Phys. Rev. Lett., 64(20), pp.2370-2373.