- Email: [email protected]
Application of the relaxation method to solve steady difference problems of convection
U.S.S.R. Comput. Maths. Math. Phys. Vol. 21, No. 1, pp. 126-137,198l. Printed in Great Britain
0041-5553/81/010126-12$07.50/O 01982. Pergamon Press Ltd.
APPLICATION OF THE RELAXATION METHOD TO SOLVE STEADY DIFFERENCE PROBLEMS OF CONVECTION*
V.K.F'OLEVIKOV Minsk (Received
26December
1978;revised
16.Iuly 1979)
AN EFFICIENT method based on the relaxation method is presented for improving the convergence of an iterative process of Seidel type for solving steady difference problems of convection, By numerical experiments optimal relaxation parameters are found for the equations of temperature, stream function and vortex in a wide range of Rayleigh numbers. An algorithm for the stabilization of the iterative process at large Rayleigh numbers is constructed.
Introduction
In view of their non-linearity the solution of steady difference convection problems requires the application of iterative procedures. In the majority of papers for increasing the stability of the algorithm the external iterative process is complemented by Internal iterations for each of the equations. The Poisson equation for the stream function is then either iterated, or solved by direct methods, for example by the fast Fourier transform method. The equations for the temperature and vortex are sometimes iterated in an internal iterative cycle independently of each other, but more frequently they are calculated by a scheme of a longitudinal-transverse pivotal condensation. Among the various ways of organizing a computational process used in practice, the method used, for example, in [l-3], which is essentially an analog of Seidel’s method, is distinguished by using the least expenditure of computer memory and time in the computation of one layer of iteration by a simple programmed realization. According to the data in [4] the amount of computer time required for one iteration of Seidel’s method is approximately one sixth of that of the method of longitudinal-transverse pivotal condensation. Its disadvantages are its unsatisfactory properties in respect of computational stability and rate of convergence in conditions of growth of convection, even in the case of monotonic approximations. It has been found to be possible to reduce the destabilizing effect of the Rayleigh number (Reynolds number) by introducing relaxation parameters into this algorithm [5-l 11. However the complex form of the convection equations practically excludes the use of the linear theory of optimization of these parameters. Obviously at the present time only numerical experiments can introduce clarity into this question. Therefore, the main content of the present paper consists of numerical experiments on test problems to investigate the optimal relaxation parameters for the convection equations in a modified Seidel method in the case of a monotonic second-order scheme. Attention is given to formalizing the information obtained with the aim of their possible use in solving other convection problems.
lZh. vjkhisl. Mat. mat. Fir., 21, 1, 127-138,
1981.
126
Steedy difference problems of convection
127
1. Formulation of the test problems Numerical experiments were carried out on two known problems. The first of them concerned the steady convective motion of a viscous incompressible fluid in a square domain with homogeneous heating from the side [ 12 - 141. The object of the second test problem is the same physical situation, but the fluid is heated sinusoidally across the top boundary (14,151. We introduce a rectangular Cartesian system of coordinates with horizontal xl -axis and ordinate x2-axis and origin at the lower left corner of the square. Then the mathematical model for the chosen test problems comprises the system of dimensionless equations [12]
a*
---
aT
ax,ax,
a+ aT =-
I
ax,
PT
_-
ax=,
a2T+ d*T
(ax,= d522)’ -
(1.1)
U-2)
a* a9 a+ acp
a*q
m----z-
ax,
ax2 ax,
ax2
+
ax,*
a*9 -
ax22
+RL aT Pr
zy
(1.3)
Here T, Ji and cpare the dimensionless temperature, stream function and vorticity, Pr is the Prandtl number, and Ra is the Rayleigh number. The stream function is connected with the velocity components by the relations Ui=a$ldx2, u2=-@lax,. It is assumed that the boundaries of the domain are fured and that on them the following adhesion condition holds:
all, [ $1rp= -an
=o.
(1.4)
rP
Boundary conditions of the first kind are specified for the temperature: in the first problem
T(0,x2)=0, T(1,x2)=1,
T(x,, O)=T(x,, i)=zi,
(1.5)
=T(O,x2)=T(1,52) -0.
(1.6)
in the second problem
T(xi, 1) =sin
nx,,
T (x1, 0)
In this paper the Prandtl number was assumed to be fured: Pr = 1.
2. The method of solution and the technique of the numerical experiment Equations (1.1) -( 1.3) were approximated by a monotonic Samarskii scheme [ 16,171 of wcondorder accuracy. The derivatives in the coefficients were replaced by central differences. The resulting system of non-linear difference equations was solved in a unified iterative cycle whose construction was based on Seidel’s method with relaxation. So that at the internal nodes the computations were performed by the formulas
V. K. Polevikov
128
8
T,= ( I-qr)
T,_,+q’
c
(Aa,+T~-‘a)
ad
+
&tJl;ef;)) [i
(Aan+&.,)
]--I,
a-1
where A an=
1 l+,R,,
,
+2G,
1
B an= l+,,R,, , -2&n-,
qJI, qT, 4~ are the relaxation parameters, n is the number of the iteration, and h is the pitch of the square mesh. The remaining notation is that generally accepted [ 171. For convenience the indexes of the mesh nodes are omitted. The boundary conditions for the stream function and temperature were specified in accordance with (1.4), (1 S) or (1.6). A difference formula of accuracy 0 (h4) (see [ 11, 181) was used to approximate cpon the boundary. For example, on the boundaryxz = 0 the vorticity was calculated in the form % = (1-q’)f+&-i+qp{
(--11,+~2) (+il,+fr) -13f$I;+-1:) -cplI-I -cp?l-1
Similarly on the other boundaries. At the corners of the square we put up= 0. The investigations were carried out on a uniform mesh Gh with the number of nodes M = 11 X 11,2 1 X 21, and 41 X 41. As the initial approximation for the iterative process in the first
problem the solution in the heat-conduction mode: To =xl,$0 = cpo = 0 wasused. In the second problem the initial state specified was the zero state: To = $0 = cpo = 0.The calculations were performed on the BESM-6 computer.
Steady difference problems of convection
129
The convergence of the iterations was regarded as achieved if the condition n
c N c1o-3 le,l
En =
m-n--N&i
0
(2.1)
was satisfied for some number n = N. Here
The number N is the least of the numbers n for which condition (2.1) is satisfied. The value of No was specified by the formula No = [2.5 M%] (see section 4). The structure of convergence was controlled by the quantity en? equal to the quantity en T (x), enG (x), 6,~ (x) of maximum modulus on the mesh Gh. To construct the optimal set of relaxation parameters qopt = {$, QI*, p}, ensuring for a given Ra the highest accuracy of convergence of the iterations, the problem of minimizing the function N(qT qJI, qq) was solved. For this the coordinate descent method was used. The computer program was so designed that the i-th approximation {QiT, fi*, ii@} was determined as the result of calculating first zi” for qT =Fiyl and q$ =FF_l, then ajti for qT =;7iTl and the q’ =?i’+’ just obtained, and finally, Fi ’ for qti = ail and 49 = pi 9. The calculations were terminated if the condition
was satisfied for some i = I, and the approximation (Gr’, glQ, g19} was regarded as the solution of the problem. The one-dimensional optimization probiems for computing the coefficients qi’, ‘iilp were solved by numerical experiment. The accuracy 6 was achieved on a sequence of meshes for contracting the domain of existence of the extremum. If there were several minima, then each of them was worked out.
It must be mentioned that if the condition Ien I < 10m3is taken as the convergence criterion, then the dependence of N on the relaxation parameters has a sawtooth form. Figure 1 shows a typical graph of this relation (the first problem, Ra = 104, h = l/20, qT = 4” = 1). This fact may considerably hamper the search for the optimum parameters. As for the integral criterion (2 .l), in the range of Rayleigh numbers investigated it ensured sufficient smoothness of the function NW q$, 49).
130
V. K. Polevikov
3. The principal results Figures 2 -4 illustrate the effect of each of the parameters qT, qJI, qq on the rate of build-up of the iterations (the first problem h = l/20, curves 1,2,3,4,5 correspond to Ra = 10,104,5.104, lOs,106). We see that for large Ra even a small deviation of any of the relaxation parameters from the optimal value rebounds ruinously on the stability of the algorithm. The greatest stabilizing effect on the calculations is exerted by the parameter qv, the least important is the relaxation of the temperature equation. For small Rayleigh numbers (Ra < 104) the relaxations of the temperature and stream function equations can be omitted altogether, since practically optimal convergence can be achieved from qv alone.
025
05
FIG. 3
0.75
?
%Z5
qT=pT, qo=cjo
I.5
Steudy difference problems of convection
ul
025
D.5
FIG. 4
0.75
I
qr-61,
I.25
1.5
131
I.75
qy
q*-g*.
Figure 5 shows curves of the components of the optimal set qopt as functions of the Rayleigh number (u is the first problem, b is the second problem, h = l/20). It is easy to see a number of relationships. In the region Ra < 103, where the principal mechanism of the spread of heat is heatconduction, the most rapid build-up of the iterations is ensured for the upper relaxation of all the functions. The sharp change of the optimal modes of relaxation at Ra - 5.104 corresponds to the beginning of the formation of the boundary layer [ 12,141. In the conditions of the formation of the boundary layer the value of TQ tends steadily to the upper limit 2 as Ra increases, and T@ and ;TT decrease to zero.
FIG The principal numerical results of this investigation are collected in Tables 1,2. The odd rows of Table 1 correspond to the first problem, the even rows to the second problem, and Table 2 refers to the first problem;m is the value of the number N for q*==p, q*=f*, p’-F. The results shown in Table 1 were obtained on the 21X2 1 mesh. It is obvious that even for Ra< 105 the relaxation noticeably reduces the build-up time. Also in conditions of highly intensive convection, when the Seidel iterations no longer converge, the relaxation method yields stationary solutions for small expenditures of computer time. The effect of the mesh on the optimal relaxation parameters is qualitatively of the same nature as in the linear case [4] : their value increases as h decreases. Some deviation from this rule was observed on coarse meshes for large Ra (Table 2). It concerned only the parameters TT and TV and was probably due to the dependence on the mesh of the coefficients of the difference equations for temperature and vortex [8].
132
V, K. Polevikov TABLE
1
-
-I
Ra
10
IE
98
::i
1.05 1.22
102
1.7 1.8
1.05 1.22
1.64 1.41
i03
1.7 1.8
:.t8
1.67 1.55
5.103
1.65 1.8
0.85 0.84
1.6 1.7
109: iti iii
10’
1.63 1.8
0.83 0.75
‘1%
151 97
3.2.10’
1.28 1.8
0.79 0.65
1.4 1.69
320 330
184 112
5.10’
1.22 I.78
0.85 1.0
1.33 1.32
339 397
189 160
1.77 1.7
0.53 0.35
no convergence
1’::
161 128
3.2.10”
1.09 1.12
1.95 1.88
0.23 0.25
no convergence
106
0.93 0.87
1.97 1.94
0.123 0.165
no convergence
280 224
E
0.042 0.085
no convergence
it
105
10’
101
481 1187
211 168
TABLE 2 Ra
I
10
3.2.10‘
lo~
Mesh
-T q
11x11 21x21 41x41
1.5 ::;
0.97 1.05 1.19
IlXll
1.23 1.28 1.46
K ok5
i-99; 11215
1.686 1.97 1.976
FE;: 0:195
2:xxz:
I
IlXll
21x21 41x41 4.
1.43 1.64 1.68
71 195 514
2::
1.06
167
109
:::2
iii
;E
no convergence
52
I73 280 524
An algorithm for stabilizing the iterative process for large Rayleigh numbers
As mentioned, temperature relaxation makes a comparatively small contribution to accelerating the convergence of the iterations even in the high Ra region. Indeed, if we put q T = 1 and find the optimal relaxation parameters of the stream function and vortex (we denote them by 7,” and qz), then, for example for Ra = 106, the number of iterations N = fl. corresponding to the set {1, &*, &o}, only exceeds the value N =8 for the best set {g’, @*, T} by the factor 1.2. The graph of # (Ra) and w (Ra) are qualitatively similar to the corresponding graphs of q* (Ra) and $@@a), shown in Fig. 5. In the range Ra> 105 they pass slightly below the curves of?@ (Ra) and $V (Ra), being situated almost symmetrically about unity. Consequently, for Ra> 105 a nearoptimal rate of convergence may be expected by solving the onedimensional problem of optimizing q(~,where 0
133
Steady difference problems of convection
Then analysis of the structure of the iterative process showed that for large Rayleigh numbers (Ra > 104) the convergence is of an oscillatory nature and is mainly controlled by the vorticity, that is, in the criterion (2.1) e, = e,lp. The period D of the oscillations of the variable-sign quantity c,, is inversely proportional to 4”. For the optimal sets (6, q”, e} and (1, Q.*, QI.O}the oscillations have approximately the same period D = B, whose dependence on Ra in the range 105 G Ra < 107, and also on the boundary conditions was found to be practically negligible. Asql~ varies upwards and downwards from the optimal value the decrease in the growth of D proceeds rapidly. Experiments showed the noticeable effect of the mesh on the oscillation period 25 in the optimal mode of convergence. It is satisfactorily described by empirical formula k-5.
&J&f’“,
(4.1)
Taking these facts into account it is easy to construct a simple algorithm for stabilizing the algorithms in conditions of highly intensive convection. We wiIl explain the essence of this algorithm, 1. We calculate the temperature for qT = 1. 2. At the beginning of the calculation we specify the first approximations qp = 419 and into account the symmetry of the graphs of $f (Ra) and $j (Ra) for large Ra, we choose them as follows:
qb = q1 J, . Taking
where al and bl are the limits of the range of variation of the optimal parameter Tf. If no more precise information about the limits of this range is available, we can put al = 0, bl = 1. 3. We provide for the storage of those numbers n for which en changes sign, and we find the period D. 4. Beginning with n =D, where 5 is defmed by formula (4.1) after every No = [D)2j iterations we compare the mean characteristics of T, and En_No, We regard as a criterion of good convergence the satisfaction of the condition E”<
(l-p,) en-iyo,
Oq3i
(4.2)
If the inequality (4.2) is satisfied, and also if n is not a multiple of No, we make no changes in the computational algorithm in the next No iterations, so that in these cases rt
Q
qn+l=qnsr
qn+1=qn
t 9
an+l=an,
bn+,=bn.
But if (42) is not satisfied, then knowing D, we analyze the structure of the iterative process in the half-interval (n-No, n 1. 5. If D
‘732
(q7l*-GA
1
q:+i=2-q:+,,
O-+-Cl.
134
V. K. Polen’kov
Then for the purpose of contracting the range of variation of these parameters we put for the upper limit b?I+1 =qR 6. The relation D >z testifies that the case Accordingly we correct the parameters: Q:+i=qn’+ps&-qn9,
(a,+i=a,). OCq’ CQev, &‘
is realized.
qn$=2-qn+l, ocpsci
By the same reasoning as before we correct the lower limit of the range of variation of the parameter 49 : 41+i=qnv
(b n+i=bn).
After every No iterations the procedure is repeated. The self-tuning of the iterative algorithm to the optimal solution proceeds in this way. Special investigations into the choice of the corrective coefficients ~1, ~2, p3 have not been made. We merely note that growth of Ra requires their decrease. In the course of the numerical experiments they were assumed to be equal to each other. Good convergence for Ra = 105,106, 107 was achieved by choosing pI = p3 = p3 = 0.6,0.5,0.2, respectively. For lower values, for example for ~1 = p2 = p3 = 0.1, the iterations also converged in all three cases, but more slowly.
FIG. 6 41~~= 1.
FIG.8uisforqlr =O.l,~isforO.O9, c is for 0.08.
FIG. 7qlq = 0.125.
Steady difference problems of convection
135
Figures 6-9 illustrate the efficiency of the proposed method of correcting the relaxation parameters depending on the initial approximation 41~ (the first problem, Ra = 106, h = l/20, p1 = p2 = p3 = 0.5; 1 is without correction of the relaxation parameters, 2 and 3 are with the stabilization algorithm included: 2 is for al = 0, bl = 1,3 is for al = 0, bl = 0.2). The points corresponding to a decrease in 4~ are marked by circles, those corresponding to an increase in q9 by triangles. We note that even in the most unfavourable case 419 = qlti = 1 the algorithm permits the algorithm to be stabilized after a number of iterations not more than double the minimal value.
Conclusion In order to clarify the degree of generality of the results obtained they have been tested for various initial data of the difference problem, although these investigations were not of a systematic and detailed nature. By the example of the second problem with Ra = 10,32.104,106 and h = l/20 the optimal parameters for the other difference schemes and boundary conditions for the vortex were roughly estimated. Schemes of monotonic type were tested: conservative schemes of the first [l] and second [6] orders of approximation, third-order schemes [ 141, and also an improved version of a monotonic scheme of fourth-order accuracy descriied in [19]. For the vortex boundary conditions of second-order accuracy of Woods and Ku&ova and the third-order conditions of Kushkova-Chudov were used. Besides the zero distribution other initial states in the iterative procedure were considered, namely solutions corresponding to smaller Ra. As the first approximation of {ai’, gi*, 4f } in the optimization procedure the FT, FQ, $@ of Table 1 were specified. For the first two Rayleigh numbers the error 6 was about 0.1, and for Ra = 106 it was approximately 0.02. The FT, a$, F9 thus obtained agreed in all cases with those given in Table 1, and the period of oscillation 5 was completely described by the relation (4.1). The results obtained were successfully checked by solving a number of new problems [9 - 11, 20,2 I]. A conservative monotonic second-order scheme was used. In the region of low and moderate values of Ra the calculations were performed on uniform meshes, mainly with step l/20, and the relaxation parameters were everywhere taken as qT = qG = 1,49 = 1.5, thus ensuring an excellent rate of convergence. In conditions of highly intensive convection more precise meshes were also used. To obtain stationary solutions for large values of Ra the stabilization algorithm described in section 4 was used. The value of8 was given by the relation (4.1) or by a similar formula. It sometimes turned out better to consider the structure of the convergence in the interval No = since for No = [z/2] the iterative process could not always be developed to a sufficient degree. Then a build-up of the iterations was obtained which did not hold for qT = qJ, = 49 = 1. As a rule the process converged after a number of iterations commensurate with that achieved in the test problems. If the rate of convergence was unsatisfactory (and this happened on meshes different from those considered in the numerical experiment), then it could be improved by changing (usually increasing) the coefficient k in formula (1 .l), by a factor of not more than 2.
[s],
It must be mentioned that for sufficiently large Ra, corresponding to the domain of the laminar boundary layer, a form of development of the iterations different from the oscillatory form was not encountered. If the iterative algorithm is regarded as some build-up scheme [17 1, then the oscillations of the iterations arising from the growth of Ra may be of a completely physical nature and be caused by turbulence of the convective flow. All this gives a basis for hoping that the basic ideas of the proposed method of stabilizing the iterations will be productive in solving many problems of convection in the mode of the laminar boundary layer.
136
V, K. Polevikov
In passing, convergence conditions weaker than (2.1), sometimes used in practice, were estimated: the build-up of three and four significant figures in the Nusselt number. As a result it became possible to make certain comparisons. Thus, the setting up of three digits in the Nusselt number on solving the fust test problem on a 21X21 mesh for Ra = 5.104 is reached in 25-30 time steps with the time step r = 0.002 (which is close to the maximum permissible), if the implicit build-up procedure proposed in [22 ] is used. As shown in [22], in stability and economy the method surpasses other evolution schemes used to solve the convection problem. On every time layer the temperature and vortex equations are computed by a scheme of longitudinal-transverse pivotal condensation, and the Poisson equation for the stream function, by iterations. The same version by the relaxation method for q r==tjT,q3=QY, qv=cj* is computed in 60-65 iterations. However in view of the approximately 6-fold difference, if [4] is followed, in the time cost of one iteration, we obtain that in the total expenditure of computer time the relaxation method is almost three times more economical. As for the region of high Rayieigh numbers, then already, for example for Ra = 106, the algorithm of [22] d oes not give build-up even for very small values of r. The enhanced stabilizing properties of the proposed relaxation method of solving steady-state convection problems are ensured by an arrangement of the computational process more flexible than that of build-up schemes. The author sincerely thanks B. M. Berkovskii for suggesting the problem, L. A. Chudov, V. I. Polezhaev, Kh. E. Kalis, V. N. Varapaev for fruitful discussions of the results, and V. L. Gryaznov for systematic material provided for comparison. Transkrted by J. Berry
REFERENCES 1. GOSMAN, A. D. et al. Numerical methods of investigating viscousfluid flows (Chislennye metody Issledovaniya techenii vyazkoi zhidkosti), “Mir”, Moscow, 1972. 2. PAU, P. E., KARLI, K. T. and KARRUT, S. L. Numerical solution of the problem of free convection in cylindrical channels of annular cross section. Teploperedacha, No. 2,78-87, 1971. 3. BERKOVSKII, B. M. and IVANOV, L. P. Threshold excitation of photoabsorption NaukSSSR. Mekhanzhidkostiigaza, No. 5,128-135,107l.
convection. Izv. Akad.
4. FEDORENKO, R. P. Iterative methods of solving difference elliptic equations. Uspekhi matem. nauk, 28, 2 (1701, 121-182,1973. 5. CRAWFORD, L. and LEMLICH, R. Natural convection in horizontal concentric cylindrical annuly. Ind. EngngC%em. Funabmentals, 1,4,260-264,1962. 6. BERKOVSKII, B. M. and POLEVIKOV, V. K. The effect of the PrandtI number on structure and heatexchange in natural convection. Inch. -fiz. zh., 24,5, 842-849, 1962. 7. GUSHCHIN, V. A. and SHCHENNIKOV, V. V. A nume-ricaI method of solving the Navier-Stokes Zh. vychisl,Mat. mat. Fiz., 14,2,512-520,1974.
equations.
8. JOHNSON, R. V. and DANAK, A. M. Heat-exchange in a laminar flow about a rectangular cavity with the injection of liquid. Teploperedochn, No. 2,84-90,1976. 9. POLEVIKOV, V. K. and FERTMAN, V. E. Investigation of the heat-exchange through a horizontal annuIar layer of a magnetic fluid when the cylindrical conductors cool with the flow. Mugnifnuyu gidrodinamikrr,No. I, 15-21,1977. 10. NOGOTOV, E. F. and POLEVIKOV, V. K. Convection in a vertical layer of magnetic fluid, situated in the magnetic field of a plate with current. hizgnitnaya gidrodinamika, No. 2,28-34,1977. 11. POLEVIKOV, V. K. Some questions of the numerical investigation of non&mar problems of heat convection by the mesh method. Candidate Dissertation, In-t matem. Akad. Nauk BSSR, 1977.
Steady difference problems of convection
137
12. GERSHUNI, G. Z., ZHUKHOVITSKII and TARUNIN, E. L. Numerical investigation of convective motion in a closed cavity. Izv. Akad. Nauk SSSR. Mekhan. zhidkostigaza, No. 5,56-62,1966. 13. FROMM, J. E. A numerical method for computing the non-linear, time dependent, buoyant circulation of air in a room.IBMJ. Res. Dev., 15,3,451-464,1971. 14. BERKOVSKY, B. M. and POLEVIKOV, V. K. Heat transfer at high-rate free convection. In: Heut Transfer 1974, Vol. 3,85-89, Tokyo, 1974. 15. BERKOVSKII, B. M. and NOGOTOV, E. F. Numerical investigation of free convection with heating from above. Izv. Abd. Nauk SSSR. Mekhan. zhidkosti iguza, No. 2,147-154,197O. 16. SAMARSKII, A. A. On monotonic difference schemes for elliptic and parabolic equations in the case of a nonselfadjoint elliptic operator. Zh. vychlrl. Mat. wt. Fir., 5,3,548-551,196s. 17. SAMARSKII, A. A. Theory of difference schemes (Teoriya raznostnykh &hem). “Nauka”, Moscow, 1977. 18. POLEVIKOV, V. K. A difference scheme of fourth-order accuracy for calculating the vortex function on the boundary in fluid dynamics problems. Dokl. Akud. Nauk BSSR, 23,10,872-875,1979. 19. BERKOVSKY, B. M. and POLEVIKOV, V. K. Numerical study of problems on high-intensive free convection. In Heat nansfer and Turbulent Buoyant Convection, 443-455. Hemisphere Publ., Washington, 1977. 20. BERKOVSKY et al. Specific features of natural convection heat transfer in magnetic fluids. In Heat ?‘Pansfer
1978. Vol. 3, 147-151. Hemisphers Publ., Washington, 1978. 21. BERKOVSKY, B. hi. et al. Heat transfer across vertical ferrofluid layers. Int. J. Heat Mass Pansfer, 19,9, 981-986,1976. 22. GRYAZNOV, V. L. and POLEZHAEV, V. P. Investigation of some difference schemes and approximations of the boundary conditions for the numerical solution of heat convection equations. Reptint In-to probl. mekhan. Akad. Nauk SSSR, 1974, No. 40.
0041-5553/81/010126-12$07.50/O 01982. Pergamon Press Ltd.
APPLICATION OF THE RELAXATION METHOD TO SOLVE STEADY DIFFERENCE PROBLEMS OF CONVECTION*
V.K.F'OLEVIKOV Minsk (Received
26December
1978;revised
16.Iuly 1979)
AN EFFICIENT method based on the relaxation method is presented for improving the convergence of an iterative process of Seidel type for solving steady difference problems of convection, By numerical experiments optimal relaxation parameters are found for the equations of temperature, stream function and vortex in a wide range of Rayleigh numbers. An algorithm for the stabilization of the iterative process at large Rayleigh numbers is constructed.
Introduction
In view of their non-linearity the solution of steady difference convection problems requires the application of iterative procedures. In the majority of papers for increasing the stability of the algorithm the external iterative process is complemented by Internal iterations for each of the equations. The Poisson equation for the stream function is then either iterated, or solved by direct methods, for example by the fast Fourier transform method. The equations for the temperature and vortex are sometimes iterated in an internal iterative cycle independently of each other, but more frequently they are calculated by a scheme of a longitudinal-transverse pivotal condensation. Among the various ways of organizing a computational process used in practice, the method used, for example, in [l-3], which is essentially an analog of Seidel’s method, is distinguished by using the least expenditure of computer memory and time in the computation of one layer of iteration by a simple programmed realization. According to the data in [4] the amount of computer time required for one iteration of Seidel’s method is approximately one sixth of that of the method of longitudinal-transverse pivotal condensation. Its disadvantages are its unsatisfactory properties in respect of computational stability and rate of convergence in conditions of growth of convection, even in the case of monotonic approximations. It has been found to be possible to reduce the destabilizing effect of the Rayleigh number (Reynolds number) by introducing relaxation parameters into this algorithm [5-l 11. However the complex form of the convection equations practically excludes the use of the linear theory of optimization of these parameters. Obviously at the present time only numerical experiments can introduce clarity into this question. Therefore, the main content of the present paper consists of numerical experiments on test problems to investigate the optimal relaxation parameters for the convection equations in a modified Seidel method in the case of a monotonic second-order scheme. Attention is given to formalizing the information obtained with the aim of their possible use in solving other convection problems.
lZh. vjkhisl. Mat. mat. Fir., 21, 1, 127-138,
1981.
126
Steedy difference problems of convection
127
1. Formulation of the test problems Numerical experiments were carried out on two known problems. The first of them concerned the steady convective motion of a viscous incompressible fluid in a square domain with homogeneous heating from the side [ 12 - 141. The object of the second test problem is the same physical situation, but the fluid is heated sinusoidally across the top boundary (14,151. We introduce a rectangular Cartesian system of coordinates with horizontal xl -axis and ordinate x2-axis and origin at the lower left corner of the square. Then the mathematical model for the chosen test problems comprises the system of dimensionless equations [12]
a*
---
aT
ax,ax,
a+ aT =-
I
ax,
PT
_-
ax=,
a2T+ d*T
(ax,= d522)’ -
(1.1)
U-2)
a* a9 a+ acp
a*q
m----z-
ax,
ax2 ax,
ax2
+
ax,*
a*9 -
ax22
+RL aT Pr
zy
(1.3)
Here T, Ji and cpare the dimensionless temperature, stream function and vorticity, Pr is the Prandtl number, and Ra is the Rayleigh number. The stream function is connected with the velocity components by the relations Ui=a$ldx2, u2=-@lax,. It is assumed that the boundaries of the domain are fured and that on them the following adhesion condition holds:
all, [ $1rp= -an
=o.
(1.4)
rP
Boundary conditions of the first kind are specified for the temperature: in the first problem
T(0,x2)=0, T(1,x2)=1,
T(x,, O)=T(x,, i)=zi,
(1.5)
=T(O,x2)=T(1,52) -0.
(1.6)
in the second problem
T(xi, 1) =sin
nx,,
T (x1, 0)
In this paper the Prandtl number was assumed to be fured: Pr = 1.
2. The method of solution and the technique of the numerical experiment Equations (1.1) -( 1.3) were approximated by a monotonic Samarskii scheme [ 16,171 of wcondorder accuracy. The derivatives in the coefficients were replaced by central differences. The resulting system of non-linear difference equations was solved in a unified iterative cycle whose construction was based on Seidel’s method with relaxation. So that at the internal nodes the computations were performed by the formulas
V. K. Polevikov
128
8
T,= ( I-qr)
T,_,+q’
c
(Aa,+T~-‘a)
ad
+
&tJl;ef;)) [i
(Aan+&.,)
]--I,
a-1
where A an=
1 l+,R,,
,
+2G,
1
B an= l+,,R,, , -2&n-,
qJI, qT, 4~ are the relaxation parameters, n is the number of the iteration, and h is the pitch of the square mesh. The remaining notation is that generally accepted [ 171. For convenience the indexes of the mesh nodes are omitted. The boundary conditions for the stream function and temperature were specified in accordance with (1.4), (1 S) or (1.6). A difference formula of accuracy 0 (h4) (see [ 11, 181) was used to approximate cpon the boundary. For example, on the boundaryxz = 0 the vorticity was calculated in the form % = (1-q’)f+&-i+qp{
(--11,+~2) (+il,+fr) -13f$I;+-1:) -cplI-I -cp?l-1
Similarly on the other boundaries. At the corners of the square we put up= 0. The investigations were carried out on a uniform mesh Gh with the number of nodes M = 11 X 11,2 1 X 21, and 41 X 41. As the initial approximation for the iterative process in the first
problem the solution in the heat-conduction mode: To =xl,$0 = cpo = 0 wasused. In the second problem the initial state specified was the zero state: To = $0 = cpo = 0.The calculations were performed on the BESM-6 computer.
Steady difference problems of convection
129
The convergence of the iterations was regarded as achieved if the condition n
c N c1o-3 le,l
En =
m-n--N&i
0
(2.1)
was satisfied for some number n = N. Here
The number N is the least of the numbers n for which condition (2.1) is satisfied. The value of No was specified by the formula No = [2.5 M%] (see section 4). The structure of convergence was controlled by the quantity en? equal to the quantity en T (x), enG (x), 6,~ (x) of maximum modulus on the mesh Gh. To construct the optimal set of relaxation parameters qopt = {$, QI*, p}, ensuring for a given Ra the highest accuracy of convergence of the iterations, the problem of minimizing the function N(qT qJI, qq) was solved. For this the coordinate descent method was used. The computer program was so designed that the i-th approximation {QiT, fi*, ii@} was determined as the result of calculating first zi” for qT =Fiyl and q$ =FF_l, then ajti for qT =;7iTl and the q’ =?i’+’ just obtained, and finally, Fi ’ for qti = ail and 49 = pi 9. The calculations were terminated if the condition
was satisfied for some i = I, and the approximation (Gr’, glQ, g19} was regarded as the solution of the problem. The one-dimensional optimization probiems for computing the coefficients qi’, ‘iilp were solved by numerical experiment. The accuracy 6 was achieved on a sequence of meshes for contracting the domain of existence of the extremum. If there were several minima, then each of them was worked out.
It must be mentioned that if the condition Ien I < 10m3is taken as the convergence criterion, then the dependence of N on the relaxation parameters has a sawtooth form. Figure 1 shows a typical graph of this relation (the first problem, Ra = 104, h = l/20, qT = 4” = 1). This fact may considerably hamper the search for the optimum parameters. As for the integral criterion (2 .l), in the range of Rayleigh numbers investigated it ensured sufficient smoothness of the function NW q$, 49).
130
V. K. Polevikov
3. The principal results Figures 2 -4 illustrate the effect of each of the parameters qT, qJI, qq on the rate of build-up of the iterations (the first problem h = l/20, curves 1,2,3,4,5 correspond to Ra = 10,104,5.104, lOs,106). We see that for large Ra even a small deviation of any of the relaxation parameters from the optimal value rebounds ruinously on the stability of the algorithm. The greatest stabilizing effect on the calculations is exerted by the parameter qv, the least important is the relaxation of the temperature equation. For small Rayleigh numbers (Ra < 104) the relaxations of the temperature and stream function equations can be omitted altogether, since practically optimal convergence can be achieved from qv alone.
025
05
FIG. 3
0.75
?
%Z5
qT=pT, qo=cjo
I.5
Steudy difference problems of convection
ul
025
D.5
FIG. 4
0.75
I
qr-61,
I.25
1.5
131
I.75
qy
q*-g*.
Figure 5 shows curves of the components of the optimal set qopt as functions of the Rayleigh number (u is the first problem, b is the second problem, h = l/20). It is easy to see a number of relationships. In the region Ra < 103, where the principal mechanism of the spread of heat is heatconduction, the most rapid build-up of the iterations is ensured for the upper relaxation of all the functions. The sharp change of the optimal modes of relaxation at Ra - 5.104 corresponds to the beginning of the formation of the boundary layer [ 12,141. In the conditions of the formation of the boundary layer the value of TQ tends steadily to the upper limit 2 as Ra increases, and T@ and ;TT decrease to zero.
FIG The principal numerical results of this investigation are collected in Tables 1,2. The odd rows of Table 1 correspond to the first problem, the even rows to the second problem, and Table 2 refers to the first problem;m is the value of the number N for q*==p, q*=f*, p’-F. The results shown in Table 1 were obtained on the 21X2 1 mesh. It is obvious that even for Ra< 105 the relaxation noticeably reduces the build-up time. Also in conditions of highly intensive convection, when the Seidel iterations no longer converge, the relaxation method yields stationary solutions for small expenditures of computer time. The effect of the mesh on the optimal relaxation parameters is qualitatively of the same nature as in the linear case [4] : their value increases as h decreases. Some deviation from this rule was observed on coarse meshes for large Ra (Table 2). It concerned only the parameters TT and TV and was probably due to the dependence on the mesh of the coefficients of the difference equations for temperature and vortex [8].
132
V, K. Polevikov TABLE
1
-
-I
Ra
10
IE
98
::i
1.05 1.22
102
1.7 1.8
1.05 1.22
1.64 1.41
i03
1.7 1.8
:.t8
1.67 1.55
5.103
1.65 1.8
0.85 0.84
1.6 1.7
109: iti iii
10’
1.63 1.8
0.83 0.75
‘1%
151 97
3.2.10’
1.28 1.8
0.79 0.65
1.4 1.69
320 330
184 112
5.10’
1.22 I.78
0.85 1.0
1.33 1.32
339 397
189 160
1.77 1.7
0.53 0.35
no convergence
1’::
161 128
3.2.10”
1.09 1.12
1.95 1.88
0.23 0.25
no convergence
106
0.93 0.87
1.97 1.94
0.123 0.165
no convergence
280 224
E
0.042 0.085
no convergence
it
105
10’
101
481 1187
211 168
TABLE 2 Ra
I
10
3.2.10‘
lo~
Mesh
-T q
11x11 21x21 41x41
1.5 ::;
0.97 1.05 1.19
IlXll
1.23 1.28 1.46
K ok5
i-99; 11215
1.686 1.97 1.976
FE;: 0:195
2:xxz:
I
IlXll
21x21 41x41 4.
1.43 1.64 1.68
71 195 514
2::
1.06
167
109
:::2
iii
;E
no convergence
52
I73 280 524
An algorithm for stabilizing the iterative process for large Rayleigh numbers
As mentioned, temperature relaxation makes a comparatively small contribution to accelerating the convergence of the iterations even in the high Ra region. Indeed, if we put q T = 1 and find the optimal relaxation parameters of the stream function and vortex (we denote them by 7,” and qz), then, for example for Ra = 106, the number of iterations N = fl. corresponding to the set {1, &*, &o}, only exceeds the value N =8 for the best set {g’, @*, T} by the factor 1.2. The graph of # (Ra) and w (Ra) are qualitatively similar to the corresponding graphs of q* (Ra) and $@@a), shown in Fig. 5. In the range Ra> 105 they pass slightly below the curves of?@ (Ra) and $V (Ra), being situated almost symmetrically about unity. Consequently, for Ra> 105 a nearoptimal rate of convergence may be expected by solving the onedimensional problem of optimizing q(~,where 0
133
Steady difference problems of convection
Then analysis of the structure of the iterative process showed that for large Rayleigh numbers (Ra > 104) the convergence is of an oscillatory nature and is mainly controlled by the vorticity, that is, in the criterion (2.1) e, = e,lp. The period D of the oscillations of the variable-sign quantity c,, is inversely proportional to 4”. For the optimal sets (6, q”, e} and (1, Q.*, QI.O}the oscillations have approximately the same period D = B, whose dependence on Ra in the range 105 G Ra < 107, and also on the boundary conditions was found to be practically negligible. Asql~ varies upwards and downwards from the optimal value the decrease in the growth of D proceeds rapidly. Experiments showed the noticeable effect of the mesh on the oscillation period 25 in the optimal mode of convergence. It is satisfactorily described by empirical formula k-5.
&J&f’“,
(4.1)
Taking these facts into account it is easy to construct a simple algorithm for stabilizing the algorithms in conditions of highly intensive convection. We wiIl explain the essence of this algorithm, 1. We calculate the temperature for qT = 1. 2. At the beginning of the calculation we specify the first approximations qp = 419 and into account the symmetry of the graphs of $f (Ra) and $j (Ra) for large Ra, we choose them as follows:
qb = q1 J, . Taking
where al and bl are the limits of the range of variation of the optimal parameter Tf. If no more precise information about the limits of this range is available, we can put al = 0, bl = 1. 3. We provide for the storage of those numbers n for which en changes sign, and we find the period D. 4. Beginning with n =D, where 5 is defmed by formula (4.1) after every No = [D)2j iterations we compare the mean characteristics of T, and En_No, We regard as a criterion of good convergence the satisfaction of the condition E”<
(l-p,) en-iyo,
Oq3i
(4.2)
If the inequality (4.2) is satisfied, and also if n is not a multiple of No, we make no changes in the computational algorithm in the next No iterations, so that in these cases rt
Q
qn+l=qnsr
qn+1=qn
t 9
an+l=an,
bn+,=bn.
But if (42) is not satisfied, then knowing D, we analyze the structure of the iterative process in the half-interval (n-No, n 1. 5. If D
‘732
(q7l*-GA
1
q:+i=2-q:+,,
O-+-Cl.
134
V. K. Polen’kov
Then for the purpose of contracting the range of variation of these parameters we put for the upper limit b?I+1 =qR 6. The relation D >z testifies that the case Accordingly we correct the parameters: Q:+i=qn’+ps&-qn9,
(a,+i=a,). OCq’ CQev, &‘
is realized.
qn$=2-qn+l, ocpsci
By the same reasoning as before we correct the lower limit of the range of variation of the parameter 49 : 41+i=qnv
(b n+i=bn).
After every No iterations the procedure is repeated. The self-tuning of the iterative algorithm to the optimal solution proceeds in this way. Special investigations into the choice of the corrective coefficients ~1, ~2, p3 have not been made. We merely note that growth of Ra requires their decrease. In the course of the numerical experiments they were assumed to be equal to each other. Good convergence for Ra = 105,106, 107 was achieved by choosing pI = p3 = p3 = 0.6,0.5,0.2, respectively. For lower values, for example for ~1 = p2 = p3 = 0.1, the iterations also converged in all three cases, but more slowly.
FIG. 6 41~~= 1.
FIG.8uisforqlr =O.l,~isforO.O9, c is for 0.08.
FIG. 7qlq = 0.125.
Steady difference problems of convection
135
Figures 6-9 illustrate the efficiency of the proposed method of correcting the relaxation parameters depending on the initial approximation 41~ (the first problem, Ra = 106, h = l/20, p1 = p2 = p3 = 0.5; 1 is without correction of the relaxation parameters, 2 and 3 are with the stabilization algorithm included: 2 is for al = 0, bl = 1,3 is for al = 0, bl = 0.2). The points corresponding to a decrease in 4~ are marked by circles, those corresponding to an increase in q9 by triangles. We note that even in the most unfavourable case 419 = qlti = 1 the algorithm permits the algorithm to be stabilized after a number of iterations not more than double the minimal value.
Conclusion In order to clarify the degree of generality of the results obtained they have been tested for various initial data of the difference problem, although these investigations were not of a systematic and detailed nature. By the example of the second problem with Ra = 10,32.104,106 and h = l/20 the optimal parameters for the other difference schemes and boundary conditions for the vortex were roughly estimated. Schemes of monotonic type were tested: conservative schemes of the first [l] and second [6] orders of approximation, third-order schemes [ 141, and also an improved version of a monotonic scheme of fourth-order accuracy descriied in [19]. For the vortex boundary conditions of second-order accuracy of Woods and Ku&ova and the third-order conditions of Kushkova-Chudov were used. Besides the zero distribution other initial states in the iterative procedure were considered, namely solutions corresponding to smaller Ra. As the first approximation of {ai’, gi*, 4f } in the optimization procedure the FT, FQ, $@ of Table 1 were specified. For the first two Rayleigh numbers the error 6 was about 0.1, and for Ra = 106 it was approximately 0.02. The FT, a$, F9 thus obtained agreed in all cases with those given in Table 1, and the period of oscillation 5 was completely described by the relation (4.1). The results obtained were successfully checked by solving a number of new problems [9 - 11, 20,2 I]. A conservative monotonic second-order scheme was used. In the region of low and moderate values of Ra the calculations were performed on uniform meshes, mainly with step l/20, and the relaxation parameters were everywhere taken as qT = qG = 1,49 = 1.5, thus ensuring an excellent rate of convergence. In conditions of highly intensive convection more precise meshes were also used. To obtain stationary solutions for large values of Ra the stabilization algorithm described in section 4 was used. The value of8 was given by the relation (4.1) or by a similar formula. It sometimes turned out better to consider the structure of the convergence in the interval No = since for No = [z/2] the iterative process could not always be developed to a sufficient degree. Then a build-up of the iterations was obtained which did not hold for qT = qJ, = 49 = 1. As a rule the process converged after a number of iterations commensurate with that achieved in the test problems. If the rate of convergence was unsatisfactory (and this happened on meshes different from those considered in the numerical experiment), then it could be improved by changing (usually increasing) the coefficient k in formula (1 .l), by a factor of not more than 2.
[s],
It must be mentioned that for sufficiently large Ra, corresponding to the domain of the laminar boundary layer, a form of development of the iterations different from the oscillatory form was not encountered. If the iterative algorithm is regarded as some build-up scheme [17 1, then the oscillations of the iterations arising from the growth of Ra may be of a completely physical nature and be caused by turbulence of the convective flow. All this gives a basis for hoping that the basic ideas of the proposed method of stabilizing the iterations will be productive in solving many problems of convection in the mode of the laminar boundary layer.
136
V, K. Polevikov
In passing, convergence conditions weaker than (2.1), sometimes used in practice, were estimated: the build-up of three and four significant figures in the Nusselt number. As a result it became possible to make certain comparisons. Thus, the setting up of three digits in the Nusselt number on solving the fust test problem on a 21X21 mesh for Ra = 5.104 is reached in 25-30 time steps with the time step r = 0.002 (which is close to the maximum permissible), if the implicit build-up procedure proposed in [22 ] is used. As shown in [22], in stability and economy the method surpasses other evolution schemes used to solve the convection problem. On every time layer the temperature and vortex equations are computed by a scheme of longitudinal-transverse pivotal condensation, and the Poisson equation for the stream function, by iterations. The same version by the relaxation method for q r==tjT,q3=QY, qv=cj* is computed in 60-65 iterations. However in view of the approximately 6-fold difference, if [4] is followed, in the time cost of one iteration, we obtain that in the total expenditure of computer time the relaxation method is almost three times more economical. As for the region of high Rayieigh numbers, then already, for example for Ra = 106, the algorithm of [22] d oes not give build-up even for very small values of r. The enhanced stabilizing properties of the proposed relaxation method of solving steady-state convection problems are ensured by an arrangement of the computational process more flexible than that of build-up schemes. The author sincerely thanks B. M. Berkovskii for suggesting the problem, L. A. Chudov, V. I. Polezhaev, Kh. E. Kalis, V. N. Varapaev for fruitful discussions of the results, and V. L. Gryaznov for systematic material provided for comparison. Transkrted by J. Berry
REFERENCES 1. GOSMAN, A. D. et al. Numerical methods of investigating viscousfluid flows (Chislennye metody Issledovaniya techenii vyazkoi zhidkosti), “Mir”, Moscow, 1972. 2. PAU, P. E., KARLI, K. T. and KARRUT, S. L. Numerical solution of the problem of free convection in cylindrical channels of annular cross section. Teploperedacha, No. 2,78-87, 1971. 3. BERKOVSKII, B. M. and IVANOV, L. P. Threshold excitation of photoabsorption NaukSSSR. Mekhanzhidkostiigaza, No. 5,128-135,107l.
convection. Izv. Akad.
4. FEDORENKO, R. P. Iterative methods of solving difference elliptic equations. Uspekhi matem. nauk, 28, 2 (1701, 121-182,1973. 5. CRAWFORD, L. and LEMLICH, R. Natural convection in horizontal concentric cylindrical annuly. Ind. EngngC%em. Funabmentals, 1,4,260-264,1962. 6. BERKOVSKII, B. M. and POLEVIKOV, V. K. The effect of the PrandtI number on structure and heatexchange in natural convection. Inch. -fiz. zh., 24,5, 842-849, 1962. 7. GUSHCHIN, V. A. and SHCHENNIKOV, V. V. A nume-ricaI method of solving the Navier-Stokes Zh. vychisl,Mat. mat. Fiz., 14,2,512-520,1974.
equations.
8. JOHNSON, R. V. and DANAK, A. M. Heat-exchange in a laminar flow about a rectangular cavity with the injection of liquid. Teploperedochn, No. 2,84-90,1976. 9. POLEVIKOV, V. K. and FERTMAN, V. E. Investigation of the heat-exchange through a horizontal annuIar layer of a magnetic fluid when the cylindrical conductors cool with the flow. Mugnifnuyu gidrodinamikrr,No. I, 15-21,1977. 10. NOGOTOV, E. F. and POLEVIKOV, V. K. Convection in a vertical layer of magnetic fluid, situated in the magnetic field of a plate with current. hizgnitnaya gidrodinamika, No. 2,28-34,1977. 11. POLEVIKOV, V. K. Some questions of the numerical investigation of non&mar problems of heat convection by the mesh method. Candidate Dissertation, In-t matem. Akad. Nauk BSSR, 1977.
Steady difference problems of convection
137
12. GERSHUNI, G. Z., ZHUKHOVITSKII and TARUNIN, E. L. Numerical investigation of convective motion in a closed cavity. Izv. Akad. Nauk SSSR. Mekhan. zhidkostigaza, No. 5,56-62,1966. 13. FROMM, J. E. A numerical method for computing the non-linear, time dependent, buoyant circulation of air in a room.IBMJ. Res. Dev., 15,3,451-464,1971. 14. BERKOVSKY, B. M. and POLEVIKOV, V. K. Heat transfer at high-rate free convection. In: Heut Transfer 1974, Vol. 3,85-89, Tokyo, 1974. 15. BERKOVSKII, B. M. and NOGOTOV, E. F. Numerical investigation of free convection with heating from above. Izv. Abd. Nauk SSSR. Mekhan. zhidkosti iguza, No. 2,147-154,197O. 16. SAMARSKII, A. A. On monotonic difference schemes for elliptic and parabolic equations in the case of a nonselfadjoint elliptic operator. Zh. vychlrl. Mat. wt. Fir., 5,3,548-551,196s. 17. SAMARSKII, A. A. Theory of difference schemes (Teoriya raznostnykh &hem). “Nauka”, Moscow, 1977. 18. POLEVIKOV, V. K. A difference scheme of fourth-order accuracy for calculating the vortex function on the boundary in fluid dynamics problems. Dokl. Akud. Nauk BSSR, 23,10,872-875,1979. 19. BERKOVSKY, B. M. and POLEVIKOV, V. K. Numerical study of problems on high-intensive free convection. In Heat nansfer and Turbulent Buoyant Convection, 443-455. Hemisphere Publ., Washington, 1977. 20. BERKOVSKY et al. Specific features of natural convection heat transfer in magnetic fluids. In Heat ?‘Pansfer
1978. Vol. 3, 147-151. Hemisphers Publ., Washington, 1978. 21. BERKOVSKY, B. hi. et al. Heat transfer across vertical ferrofluid layers. Int. J. Heat Mass Pansfer, 19,9, 981-986,1976. 22. GRYAZNOV, V. L. and POLEZHAEV, V. P. Investigation of some difference schemes and approximations of the boundary conditions for the numerical solution of heat convection equations. Reptint In-to probl. mekhan. Akad. Nauk SSSR, 1974, No. 40.