. : http://hit-conf.imec.msu.ru/2006/CD/Palymskiy2.doc
: Wed Feb 22 16:10:18 2006
: Fri Feb 28 20:07:03 2014
:

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 e682t). 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 dimensionless set of equations given in
terms of deviations from an equilibrium solution is [1,7]:

[pic]

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 = g?3dQ/?? 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:

[pic]
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 ?x = ?x = Q = 0 with
x = 0, L; 0 y 1 (on the vertical boundaries), L = ?/?.


NUMERICAL METHOD

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.
[pic]

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.
[pic]
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),
. 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
[pic]
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 = 1000Ra, Pr = 1, N = 64 and

M = 16, the time step ? is equal to 410-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.

[pic]
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/4?1/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 = 10000Racr now.
[pic]
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.610-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].

[pic]
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]:
[pic]

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.
Unfortunately, the linear analysis doesn't allow to investigate the time
approximation of nonlinear terms. Only tools of nonlinear analysis may make
it.


NONLINEAR ANALYSIS

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

[pic]

This nonlinear system has a family of private waveform solutions:

[pic]

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.
[pic]
[pic]
[pic]

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.


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].

[pic]
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 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.

[pic]
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:
[pic]

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.
[pic]
Figure 6
Temperature probability density at r = 10000


[pic]
Figure 7
Alpha versus supercriticality

Experimental result [13] was received for high values of supercriticality
1.2105 r 1.2107. 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.6104 r 3.2104 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.

[pic]
Figure 7
Nu-r power laws
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.


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.139-179.
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.
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.