- Email: [email protected]
Coupled solution of the steady compressible Navier-Stokes equations and the k−ε turbulence equations with a multigrid method
MATHEMATICS ELSEVIER
Applied Numerical Mathematics 23 (1997) 49~51
Coupled solution of the steady compressible Navier-Stokes equations and the k-e turbulence equations with a multigrid method E. Dick*,
J. Steelant
Department of Mechanical and Thermal Engineering, Universiteit Gent, Sint-Pietersnieuwstraat 41, B-9000 Gent, Belgium
Received 30 November 1995; revised 13 May 1996
Abstract A relaxation method for the steady turbulent compressible Navier-Stokes equations is developed. The flow equations are fully coupled to the turbulence equations. The principles of the method are illustrated for three different k-e models. The relaxation method can be used in multigrid form. Multigrid results are given for one turbulence model on a fiat plate test case. Keywords: Navier-Stokes equations; Turbulence; Relaxation methods; Multigrid methods
1. Introduction Basically, two classes of methods for the turbulent compressible Navier-Stokes equations are used nowadays. A first class of methods consists of the explicit time stepping schemes. In the multigrid context, we cite as examples the method of Dimitriadis and Leschziner [2] using explicit Lax-Wendroff stepping and the low Reynolds number Chien model and the method of Gerolymos [3] using also Lax-Wendroff stepping but the Launder-Sharma model. These methods suffer from a severe time step restriction due to the presence of the source terms in the turbulence equations. This handicaps the optimization of the multigrid method. Properly defined implicit time stepping schemes overcome this difficulty. Examples of methods of implicit type are the methods of Vandromme and Ha Minh [12] and Morrison [7]. The first method uses the central discretization, the last uses the upwind discretization. These methods are not used in multigrid form. Some methods are of hybrid type like the method of Mavriplis and Martinelli [6], where the Navier-Stokes part is solved explicitly and turbulence part implicitly and where the Navier-Stokes equations are discritized in the central way and the turbulence equations in the upwind way. * Corresponding author. E-mail: [email protected]. 0168-9274/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S0168-9274(96)00061-X
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
50
Implicit time stepping methods generally allow very large time steps. This observation indicates that it is possible to develop methods that solve directly the steady equations without the use of time stepping. Such a method is then equivalent to an implicit time stepping method with infinite time step. In this paper we develop a relaxation method for the steady equations and employ it in multigrid form. To ensure solvability by a relaxation method, the system of discrete equations has to be so-called positive. To obtain positiveness, the different parts of the equations have to be treated carefully: the convective part, the diffusive part and the source part. These three parts have to be split into positive and non-positive contributions. The positive contributions form the left-hand side of the system; the non-positive contributions form the right-hand side. So, with adequate treatment of the different terms in the equations, a solution of the steady equations can be obtained without the need of time stepping. In multigrid methods employing the turbulence equations on the coarse grids [2,3,6], often the coarse grid corrections cause the turbulence quantities k and e to become negative. This leads to breakdown of the calculations. So it is necessary to damp the coarse grid corrections. In [3,6], this is done by limiting the coarse grid corrections to some upper values. Here we follow a more sophisticated approach, explained later, developed by Lien and Leschziner [5] for incompressible Navier-Stokes equations. In this paper, we use the low Reynolds number form of the turbulence equations. We consider two variants of the multigrid formulation. In the first variant the k-e equations are not used on the coarse grids. In the second variant, the k-e equations are solved on the coarse grids, but the coarse grid corrections are damped. We compare the efficiency of both methods and come to the conclusion that it is advantageous to employ the turbulence equations on the coarse grids.
2. The equations We consider the set of steady compressible Navier-Stokes equations for turbulent flow, coupled to a k-¢ turbulence model. In two dimensions, the equations take the form: OF
aG
OFv
OG~
where
F =
pu
pv
puu + p + ~pk
puv
pvv + p + 2 pk
puv
p u e + (p +
G =
pk) u
'
p v E + (p + ~pk)v
puk
pvk
puz
pv~
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
51
o
0
0
7-xx
7-zy
0
Txy
~y
F~ =
,
Gv
S =
-~
0
UTxy + VTyy Jr- qy
urzx + vrxy + qx
0
(].t+ ,t/O'k)]¢x (JZq- ,t/O'e)Ex
Sk \ S~ J
The components of the stress tensor T and the heat flux vector q are given by rxx = (# + #t) [2ux - }(Ux + Vy)],
+ v~),
~-~ = ~-~ = (~, + ~ , ) ( ~
qx=
( Iz Idt) ~+Prt
ryy = (# + #t) [2Vy - ~(uz + vy)],
(~
qv=
hx,
lZt~h
~+PrtJ
y"
The molecular viscosity is denoted by # and the turbulent viscostity by/.z t. Pr is the molecular Prandtl number (taken as 0.71 for air) and Prt the turbulent Prandtl number (taken as 0.91 for air). The static enthalpy is denoted by h. E stands for the total internal energy per unit mass: E-
--1
p+
1 u2 + 1
V2+ k .
Source terms of a k-e model can be written as
Sk = Pk - p e -
D,
S~ = [ C e , f l P k - C ~ 2 f 2 P e ] ~
+ E,
with
{ [ou, ek=
~,LOxj + 7x~
~6ij Oxk J - - ~ i j p k
-
~xj' -
(1)
lZt = Cufupk"]-,
where Cs,, C~2, Ci,, ak and o-~ are the standard k-e model constants; T is the turbulence time scale; fl, f2 and fu are the so-called wall proximity damping functions and D and £ are the low Reynolds number terms. In our study, three models are used: Launder-Sharma (LS) [8], Lam-Bremhorst (LB) with the modification of the f2-function [10] and Yang-Shih (YS) [14]. The five basic constants are: CE, = 1.44, C~- 2 = 1.92, ak = 1, as = 1.3, Cu = 0.09. Table 1 gives the low Reynolds terms Table 1 Low Reynolds terms and boundary conditions for e
k-e LS
T k -
ew - B.C.
79
0
2 . \ i~y ]
2/zut \0--~2 }
C
2
$
(a2u
LB
k e
-- = 0
Oe i~y
0
0
YS
_k+~/-~e
2 # \ Oy ]
o
(02u~ 2 ~, \~-Ty~/
52
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
Table 2 Damping functions. RT = kT-/u, R u = v/-ky/v, aj = 1.5 × 10 4 a3 = 5.0 × 10-7, a5 = 1.0 × 10-I° k-e
f~
fl
fe
LS
exp[-3.4/(1 + RT/50) 2]
1
1 -
0.3 exp(-R~,)
LB
[1 - exp(-0.0165Ru)]2(1 + 20.5/RT)
1 + (0.05/f~) 3
1 -
exp (-R~. - 10-1°)
YS
[1 - exp ( - a , n u - a3R~ - asR~)] '/2
1
1 - 0.22 exp (-R2/36)
together with the boundary conditions for these models and Table 2 gives the damping functions. The distance to the wall is denoted by y. Upstream of a leading edge, the wall distance is taken as the distance to the leading edge.
3.
The
discretization
In the next sections, the discretization of the different parts (convection, diffusion and source) will be discussed and their role in the relaxation method will be detailed. 3.1. Convective terms The convective part is treated by a flux-difference splitting method. Here, the polynomial fluxdifference splitting [1] is employed which is a variant of Roe-splitting. Using the polynomial character of the convective fluxes F and G with respect to the vector of primitive variables W T = {p, u, v , p , k, e}, differences of flux vectors can be expanded as:
p 2fz2 + -~k ~
AF=
~
+ fit ~0
0
0
0
0
0
1
2 -~fi
0
p--~
0
0
2 p---~+ -~pu--
0
5~ + ~k~2
a42
puv
0 "7 "7 _ 1 ~t
f~~
-p k
o
o
p--a
0
gft
pe
0
0
0
pu
!
AW = A1AW,
where (1 = ~ - ~ +
2 vz,
a42 = p q + - ~ f t +
'71p+/~k+ "7-
3P--k.
The bar denotes the algebraic mean value. AG is treated similarly. Most of the terms in the above expression are self-evident. Some ambiguity arises in treating the triad terms. The following choices are made. The term p u u in the momentum equation is seen as the product of the mass flux pu with the velocity component u. This leads to A p u u = - ~ A u + ftApu = ~--~Au + ~(~Ap + fiAu). The other
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
53
possible grouping into p and u 2 being a density and twice a kinetic energy has no physical meaning. The energy flux is ( 2 )
1
1
2
1
2
( 2 )
puE + p + -spk u = 7 - 1pu + -~puu + -~puv + puk + p + -~pk u. The term ½PUU2 is seen as the product of a mass flux pu with a kinetic energy ½u2. Similarly the term puk is a product of a mass flux and the turbulence kinetic energy. The term (p + 2pk)u is the product of the effective pressure (molecular plus turbulent pressure) and the velocity u. This means that the groupings in the term ~pku and the convective term puk are different. The above choice of groupings (i.e., following the physical meaning of the terms) leads to the simplest expressions afterwards. The basic splitting is done with respect to the primitive variables W. A transformation to a splitting with respect to the conservative variables U T = {p, pu, pv, pE, pk, pe} is obtained from .
1
AU= q+/¢
g
0
0
0
0
O
/5
0
0
0
0
0
/5
0
0
0
/sfi
/50
1/(7 - 1)
/5
0
o
o
o
p
0
0
0
0
0
/5
AW = TAW.
The expression A F = A'IAW transforms into AF = A'IT-1AU = A1AU. Similarly AG = A'2T-1AU = A2AU. The flux-differences are split into a positive and a negative part according to the sign of the eigenvalues of the matrices A1 and A2:
AF = A+AU + ATAU,
AG = A+AU + A2AU.
We consider a vertex-centred finite volume formulation as shown in Fig. 1. The flux-difference over the surface Si+l/2 with length ASi+l/2, is .T'i+! -- .:Ki = A~i,i+l ---- (nxAFi,i+l + nyAGi,i+l )A8i+1/2 = (nxAl + nyA2)AUi,i+lASi+l/2
+ + Az,i+l)AUi,i+lASi+l/2" ---- Ai,i+lAUi,i+lASi+l/2 = ( A i,i+l Splitting the matrix A into a positive and negative part allows the definition of the absolute value of the flux-difference by + IA-T'i,i+ll = ( A i,i+l -- A~,i+l)AUi,i+lASi+l/2 •
(2)
Based on (2) a first order accurate upwind definition of the flux is
,~7i+1/2 __-- 21 [-~"i -I- -~i+1 -- ]~i,i+1]]
(3)
= d~i + ASi+l/2AT,,i+lAUi,i+l.
Fluxes on the other surfaces of the control volume can be written analogously. Summation over the surfaces expresses the inviscid flux balance of a control volume, where the index k refers to the surrounding nodes:
Z k
A;Ask(Uk-Uo)=0,
Z(-A;Ask)Uij k
- Z-A;AskUk k
=0.
(4)
POLYMER SCIENCE U.S.S.R. This Journal is published twelve times per year by Pergamon Press and is a cover-to-cover translation of Vysokomolekulyarnye soyedineniya (series A) (Russian original copyright: Nauka, 90 Profsoyuznaya ul., Moscow 117485)
Scientific Translation Editor: Dr. W. COOPER 217 Little Aston Road, Aldridge, Walsall W S 9 0 P A , West Midlands
Translators: C. W. CAPP, R. J. A. HENDRY, G. F. MODLEN, A. CROZY M. KUBIN, J. HORSKA, D. DOSKOCILOV/, and Mrs. E. SEMERE
Chairman of the Honorary Editorial Advisory Board: SIR HARRY MELVILLE,
K.C.B., F . R . S .
Honorary Editorial Advisory Board for this English version: G. R. ALLAN, London P. E. M. ALLEN, Adelaide C. H. BAMFORD, Liverpool J. W. BARRETT, Cheltenham J. C. BEVINGTON, Lancaster S. CLAESSON, Uppsala P. CORRADINI, Naples J. B. DONNET, Mulhouse P. J. FLORY, Stanford J. FURUKAWA, Kyoto G. GEE, Manchester
N. L. P. A. A. C. C.
GRASSIE, Glasgow HOLLIDAY, London MEARES, Aberdeen A. MORTON, Cambridge, Mass. M. NORTH, Glasgow G. OVERBERGER, Michigan V. ScstrLZ, Mainz
G . SMETS,
Louvain
H. M. SPURLIN, Wilmington H . L . WILLIAMS,
Toronto
B. H. ZIMM, California
Publishing and Advertising Offices Pergamon Press Ltd., Headington Hill Hall, Oxford OX3 0BW, England Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford, NY 10523, U.S.A.
1986 Subscription Rates For Libraries, university departments, government laboratories and other multiple-reader institutions: U.S. $650.00 per annum (including postage and insurance). Two year subscription rate (1985/1986): U.S. $1235.00. Subscription enquiries from customers in North America should be sent to: Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford, NY 10523, U.S.A., and for the remainder of the world to: Pergamon Press Ltd., Headington Hill Hall, Oxford OX3 0BW, U.K.
Microform Subscriptions and Back Issues Back issues o f all previously published volumes are available in the regular editions and on microfilm and microfiche. Current sub scriptions are available on microfiche simultaneously with the paper edition and on microfilm on completion o f the annual index.
Copyright Q 1986 Pergamon Press Limited U.S. Copyright Law Applicable to Users in the U.S.A. Ph otocopying information for users in the U.S.A.: T h e Item-Fee Code for this publication indicates that authorization to photocopy items for internal or personal use is granted by the copyright holder for libraries and other users registered with the Copyright Clearance Center (CCC) Transactional Reporting Service provided the stated fee for copying beyond that permitted by Section 107 or 108 of the United States Copyright Law is paid. The appropriate remittance of $10.00 per copy per article is paid directly to the Copyright Clearance Center Inc., 21 Congress Street, Salem, MA 01970. The copyright owner's consent does not extend to copying for general distribution, for promotion, for creating new works, or for resale. Specific written permission must be obtained from the publisher for such copying. In case of doubt please contact your nearest Pergamon office. The Item-Fee Code for this publication is: 0032-3950/84 $10.00+.00
P U B L I S H E D BY
PERGAMON PRESS LTD. OXFORD. NEW YORK B E I J I N G • F R A N K F U R T • SAO P A U L O • S Y D N E Y • T O K Y O • T O R O N T O No part of this publication may be reproduced stored in a retrieval system or transmitted in any form or by any means: electronic, electrostatic, magnetic tape, mechanical, photocopying, recording or otherwise without permission in writing from the publishers.
POLYMER SCIENCE: U.S.S.R. 1984
V O L U M E 26
N U M B E R 12
CONTENTS J. VARGA,J. MENCZEL a n d A. SOLTI
2763
B. V. LEBEDEV,T. G. KULAGINA, V. S. SVISTUNOV, V. S. PAPKOV and A. A. ZHDANOV S. V. BRONNIKOV, V. ]. VETTEGREN', L. N. KORZHAVIN and S. YA. FRENKEL' V. A. PoPov, V. V. GUZEYEV, Yo. A. ZVEREVA, A. N. GRISHIN, T. V. PALAYEVA, A. P. SAVEL'EV and S. ~q. POTEPALOVA V. V. SHILOV, T. E. LIPATOVA, V. V. VORONA, V. N. BLIZNYUK and A. I. SNEGIREV A. I. MAKLAKOV, V. A. SEVRYUGIN, V. D. SKmDA and N. F. FATKULLIN V. V. SHILOV, YE. I. ORANSKAYA and Yu. S. LIPATOV N. V. POGODINA, A. B. MEL'NIKOV, O. 1. MIKRYUKOVA, S. A. DIDENKO and G. N~ MARCHENKO A. V. MAKSIMOV, Yu. YA. GOTLIB and V. G. BARANOV A. G. KORCHAGIN, M. A. MARTYNOV and S. A. TSYGANKOV M. M. KOTON, V. V. KUDRYAVTSEV, V. A. ZUBKOV, A. V. YAKIMANSKII, T. K. MELESHKO and N. N. BOGORAD Yu. I. MITCHENKO, T. I. LEBEDEVA and R. F. TSIPERMAN V. P. BUDTOV, YE. L. PONOMAREVA, N. I. KONDRASHKINA and T. N. ZELENKOVA V. I. SKR1PACHEV, V. N. KUZNETSOV and S. S. IVANCHEV N. A. YURENKO, G. T. YEVTUSHENKO, YU. YE. YERMILOVA, 1. M. SHOLOGON and B. A. ROZENBERG M. V. TSEBRENKO and G. P. DAN~LOVA
2773
A. Yu. BILIBIN,A. V. TEN'KOVTSEV, O. N. PIRANER and S. S. SKOROKHODOV L. YA. MADORSKAYA, V. M. SAMOILOV, G. A. OTRADINA, A. P. AGAPITOV, V. P. BUDTOV, T. G. MAKEYENKO, YE. YU. KHARCHEVA and N. N. LOGINOVA T. I. BORISOVA, L. L. BURSHTEIN, T. P. STEPANOVA and V. P. SHInAYEV
2781 2789
Calorimetric study of the crystallization and melting of polypropylene Calorimetric study of 1,1,3,3,5,5-hexaethylcyclotrisiloxane and polydiethylsiloxane and the process of polymerization of 1,1,3,3,5,5-hexaethylcyclotrisiloxane in the range 13-330 K Low temperature dependence of the strength of polymide fibres Polymerization modification of fillers and its influence on the properties of filled polyvinyl chloride
2797
Supramolecular organization of some block copolyurethanes with saccharide units in the main chain
2804
Self-diffusion of macromolecules in polystyrene solutions
2811
Microphase structure of block copolyurethanes with a different length of the rigid chain blocks Conformational characteristics of low-molecular weight collodions as shown by diffusion-sedimentation analysis
2819
2825 2834 2839
Orientational order and relaxation times in two-dimensional models of polymer systems with inter-chain interactions Structure of high-density polyethylene obtained by hydrostatic extrusion in solid state Experimental and theoretical study of the effect of medium on chemical imidization
2848
Molecular interactions in polymeric amido-salt soluttons
2856
The change in the molecular mass distribution of photosensitized low-density polyethylene during atmospheric aging Th~ mechanism of the modifying action of oligomers in filled compositions based on polyethylene Effect of chemical isomerism of glycidyl diphenolic ethers and esters of dicarboxylic acids on dielectric relaxation in epoxy polymers Features of the fracture of ultrathin fibres in extrudates of polymer blends Synthesis of high-molecular weight liquid crystal polyesters based on a polycondensation mesogenic m o n o m e r Structure of the endgroups of polyvinylidene fluoride synthesized in the presence of fl-hydroxyethyl-tert-butyl peroxide
2862 2867
2872 2882 2891
2897
Influence of the chemical structure of side chains of comb-like polymers on the gel point
(Continued inside back cover} I S S N 0032-3950 Printed in Poland (PWN-DUAM)
257
E. Dick, J. Steelant /Applied Numerical Mathematics 23 (1997) 49-61
56
It is straightforward to verify that all eigenvalues are non-negative. The contribution of the differences in the ~ direction in the diffusive flux balance forms a positive system. The resulting positive system in the left-hand side is expressed in the variables {p, u, v, e, k, e}. The transformation to the conservative variables is done by
1
0
0
0
0
O
~t
p
0
0
0
0
0
~
0
0
0
p~
p~
p
p
o
o
o
o
p
o
0
0
0
0
p
AU=
q + f~ + ~
AW = TAW.
To transform we use AW = T-1AU. The transformation does not change the positive character of the left-hand side.
3.3. Source terms Whereas the construction of the convective and diffusive Jacobians is rather straightforward, this is not the case for the source terms. For the source terms, a proper linearization must be chosen. The Jacobian of the negative source term is then to be brought into the left-hand side to increase the diagonal dominance of the system of equations [13]. Dependent on the k-e-model, a different linearization is necessary. Until now we have
[
~-~(-A k + B~,k)Ask]Uij - ~-~(-A-~ + B~,k)AskUk = R2 + SV, k
k
where the sum extends over the surfaces of the control volume, Ask denotes the length of the surface and k refers to a surrounding node. The term R2 collects the second-order contributions from the inviscid terms and the tangential contributions from the viscous terms. V denotes the volume of the control volume. The source term still has to be treated. A typical relaxation scheme used in the sequel is three Gauss-Seidel steps on the left-hand side (inner iteration) between updates of the right-hand side (outer or defect-correction iteration). Put into ~-formulation, we obtain
~ - ~ ( - A ~ + B~,k)AskSUij ~- F1 = R2 d- S V , k
where F1 is the flux balance based on first-order inviscid fluxes and the normal part of the viscous fluxes. This term is partly on an old iteration level and partly on a new iteration level. The coefficients A~- and B~,k are also partly on new and old iteration levels as these are updated during the inner iterations.
57
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
The source term S is split into positive and negative terms. The negative term is put into the left-hand side and takes part in the inner iteration: S = S + + S - + -0-O--~U~j. So we obtain:
{ Z ( - A - ~ + B~,k)Ask- OS-V}
(6)
k
The negative source terms taken into consideration for linearization are the same for all models: 1
S[ = -[ps 4- D],
S~- = -C~2f:ps -~.
We consider first the LS and LB model. Following Vandromme [13], we write these terms as:
I_ #t J
-- -~ pk,
S~-=-[Ce'f2](Pe)2- pk
In particular, in the expression of S~-, pe has been substituted to pk using the expression of #t (1) in order to increase the diagonal dominance. Considering the quantities in square brackets to be constant, a linearization which guarantees positiveness, is: ~(S~-; S~-)
O(pk; pg)
_
OSk ~pk DS~apk
OS-~ Op~ = OS: Opg
2Cuf~ (pk) - 79_ 0 #t k C~2fe ( pe'~ 2 pg \-~ ] -2Ce2f2-~
For the YS model, the same combinations are kept constant, but due to the different expression of 7-, we have
Sk ~__(pk)24 (V~ 4- ~ ) 2 =
S~- : _[C~2f2](pg)___~, = -- + pg
.
(7)
This results in the Jacobian
) ~(pk;p~)
,'y(V/~ 4- ~ ) 2 ( p / ~ ) Ce2f2~ -
0 1
p
_Ce2f2 ( a T _ ~ V ~ )
~
1
The resulting global matrix formed by the left-hand sides of equations of type (6) still is a generalised M-matrix in the sense that by reduction of all matrix coefficients to positive scalars the global matrix becomes an M-matrix. We observe, without rigorous proof, that this generalised system can be solved by any collective relaxation method. The solvability of the system of equations resulting from the discretization of the Navier-Stokes equations by different relaxation methods is analysed by the authors in another paper [11]. Here we illustrate the practical applicability of relaxation methods in multigrid form with the example of Gauss-Seidel relaxation.
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
58
4. ~ s t c ~ e Flow with a zero-pressure gradient was calculated over a flat plate with a freestream turbulence level of 3% (T3A test case from Savill [9]). A stretched grid of 385 × 97 points was used (Fig. 2). The grid extends upstream of the plate, with the sharp leading edge at station 97. The first grid point in the direction normal to the plate lies at about y+ -- yu~/u -- 1, where u~- is the friction velocity. Stretching was applied normal to the plate and in the flow direction near the leading edge. Uniform inlet profiles for total temperature, total pressure, k and e were specified. At inlet, Mach number was extrapolated from the flow field. The values of k and e at the inlet were calculated with the equations for k and e for uniform flow with velocity U:
U ok U oE Ox Ox where at the leading edge the following values were matched to be in accordance with the experiments (for L = lm): =
ke = 0.03(3Ue2),
=
ee = 0.378m2/s 3,
Ve =
5.4m/s.
The upper and right boundaries are outlet boundaries. There, pressure was imposed. Velocity components, temperature and turbulent quantities were extrapolated. The part of the lower boundary upstream of the leading edge was treated as a symmetry line. At the plate, no-slip and adiabatic boundary conditions were imposed. Density and pressure were obtained by characteristic combinations of the equations [ 1]. Three lexicographic Gauss-Seidel relaxations with underrelaxation factor 0.9 were used as inner iteration. The first relaxation starts from the left bottom point and ends at right upper point. The second relaxation has the reversed ordering. The third relaxation has the same ordering as the first one. Fig. 3 shows the obtained distribution of the skin friction coefficient for the three models. The upper and lower lines correspond with fully laminar and fully turbulent flow fields. Transition point i
i
!!,ll~lll]ll]l IIIII:III:I
1
_
] I I I J
1 mra
Fig. 2. Detail o f the grid at the l e a d i n g edge.
E. Dick, J. Steelant /Applied Numerical Mathematics 23 (1997) 49-61
59
0.010.0090.0080.007. 0.006-
~3
0.005 0.004 0.003 0.002 0.001
160
260
360
460
660
600
Rex (Thousands)
Fig. 3. Skin friction coefficient predicted by the models.
and transition length are not well predicted by all models, when compared to experimental results [9]. This shows that the models still have to be much improved. The convergence results are similar for all models and are discussed in the next section.
5. Multigrid formulation A standard multigrid method using four grids (385 x 97; 193 x 49; 97 x 25; 49 x 13), W-cycle, full weighting as restriction for residuals and bilinear interpolation as prolongation, was employed. The multigrid acts on the left-hand side of the set of equations. The right-hand side is updated in a defect correction cycle. The procedure is the same as the one used for the Euler equations in [1], except for the ordering of the relaxation. Three Gauss-Seidel relaxations are used as prerelaxation and postrelaxation. Two variants of the multigrid formulation are considered. In the first, the turbulence equations are not used on the coarse grids. In the second, the turbulence equations are used on all grids. In this last variant, without damping of the coarse grid corrections for the turbulence quantities, negative values for ]C and e are always encountered leading to a breakdown of the calculation. Therefore, we use the following damping [5]: ]chew --
]Cold + ar]c +
]cold -- aS]c-
]Cold,
(8)
with a similar expression for e, where 5]C+ and 5 k - denote the positive and negative coarse grid correction respectively. The damping factor a is less than 1. For a = 0, the first variant is recovered and the solution of the turbulence equations on the coarse grids can be omitted. A relaxation on the current grid is taken as one local work unit. A residual evaluation plus the associated grid transfer is also taken as one local work unit. Local work units are counted as fractions of work units on the level of the finest grid proportional to the number of cells. The update of the right hand side of the system of equations in the defect correction is also taken as one work unit. With
60
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
1 E-09,1,,~ 1E-1( 1 E-11 1E-12 1E-13~t~ 1E-1, 1E.15. t ....
1E-16]0
1'0
o's
e'0
~s
Work Units (thousands)
Fig. 4. Convergence history for single grid and multigrid calculations for the YS model. Logarithm of the residual as function of work units. these rules, the cost of the cycle is found to be 12.04167 for the first variant and 14.0625 work units for the second. In a single grid calculation, the defect correction cycle is counted as 4 work units. The calculation starts from uniform flow. First a laminar solution is calculated up to a sufficient level of convergence. With this solution, initial values of k and e are calculated according to the boundary layer laws [9]: k = ki
u
,
~ = 0.3k
with ~ ~> ~i,
where the subscript i refers to inlet conditions. Fig. 4 shows the convergence behaviour for the YS model starting from the laminar solution and the initial values of k and ~. The residual shown is the maximum residual over all equations and all nodes on the finest grid at the end of the cycle. The residuals are the residuals of the flux balances. These contain the size of the control volumes resulting in small initial values in Fig. 4. The best multigrid performance is obtained with the second variant for the largest possible value of the damping factor. It was found experimentally that above some critical value of this damping factor, depending on the flow problem, breakdown of the convergence occurs. The best damping value is just below this critical value. For the flat plate problem the critical value is slightly above 0.3. Without solving the turbulence equations on the coarse grids, the convergence is less good, although the cycle is less expensive. For o~ = 0 and 0. l, off-levelling of the convergence is seen. For comparison, in Fig. 4, also the convergence behaviour of the single grid calculation is shown. Although its asymptotic convergence still cannot be seen in the figure, the big gain in convergence speed due to multigrid is clear.
6. Conclusion A relaxation method for the steady turbulent compressible Navier-Stokes equations coupled to the low Reynolds number k-e turbulence equations was developed. Therefore, a careful treatment of the
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49~51
61
convective, diffusive and source terms was necessary. The relaxation method was used in multigrid form. The turbulence equations cannot be used directly in the multigrid formulation. Damping of the coarse grid corrections for the turbulence quantities is necessary. A big gain of convergence speed in comparison with the single grid calculation is obtained.
Acknowledgements The research reported here was granted under contract 9.0001.91 by the Belgian National Science Foundation (N.EW.O.) and under contract IUAP/17 by the Federal Services of Scientific, Technical and Cultural Affairs (D.W.T.C.).
References [ 1] E. Dick, Multigrid solution of steady Euler equations based on polynomial flux-difference splitting, Internat. J. Numer. Methods Heat Fluid Flow 1 (1991) 51-62. [2] K.E Dimitriadis and M.A. Leschziner, Multilevel convergence acceleration for viscous and turbulent transonic flows computed with a cell-vertex method, in: J. Mandel et al., eds., Proceedings 4th Copper Mountain Conference on Multigrid Methods (SIAM, Philadelphia, PA, 1989) 130-148. [3] G.A. Gerolymos, Implicit multiple-grid solution of the compressible Navier-Stokes equations using k-e turbulence closure, AIAA J. 28 (1990) 1707-1717. [4] D.C. Jespersen, A multigrid method for the Euler equations, AIAA-83-0124 (1983). [5] ES. Lien and M.A. Leschziner, Multigrid acceleration for recirculating and turbulent flows computed with a non-orthogonal, collocated finite-volume scheme, Comput. Methods Appl. Mech. Engrg. 118 (1994) 351371. [6] D.J. Mavriplis and L. Martinelli, Multigrid solution of compressible turbulent flow on unstructured meshes using a two-equation model, lnternat. J. Numer. Methods Fluids 18 (1994) 887-914. [7] J.H. Morrison, A compressible Navier-Stokes solver with two-equation and Reynolds stress turbulence closure models, NASA Contractor Report 4440 (1992). [8] V.C. Patel, W. Rodi and G. Scheuerer, Turbulence models for near-wall and low Reynolds number flows: a review, AIAA J. 23 (1984) 1308-1319. [9] A.M. Savill, A synthesis of T3 test case predictions, in: O. Pironneau et al., eds., Numerical Simulation of Unsteady Flows and Transition to Turbulence (Cambridge University Press, Cambridge, 1992) 404442. [10] K. Sieger, A. Schulz, M.E. Crawford and S. Wittig, Comparitive study of low-Reynolds number he turbulence models for predicting heat transfer along turbine blades with transition, in: Proceedings International Symposium on Heat Transfer in Turbomachiner3,, Athens (1992). [11] J. Steelant, S. Pattijn and E. Dick, Smoothing properties of different preconditioners for Navier-Stokes equations, in: Proceedings 3rd Eccomas CFD Conference, Paris (1996). [12] D. Vandromme and H. Ha Minh, About the coupling of turbulence closure models with averaged NavierStokes equations, J. Comput. Phys. 65 (1986) 386-409. [13] D. Vandromme, Turbulence modeling for compressible flows and implementation in Navier-Stokes solvers~ VKI-LS 1991-02. [14] Z. Yang and T.H. Shih, New time scale based k-e model for near-wall turbulence, AIAA J. 31 (1993)
1191-1198.
Applied Numerical Mathematics 23 (1997) 49~51
Coupled solution of the steady compressible Navier-Stokes equations and the k-e turbulence equations with a multigrid method E. Dick*,
J. Steelant
Department of Mechanical and Thermal Engineering, Universiteit Gent, Sint-Pietersnieuwstraat 41, B-9000 Gent, Belgium
Received 30 November 1995; revised 13 May 1996
Abstract A relaxation method for the steady turbulent compressible Navier-Stokes equations is developed. The flow equations are fully coupled to the turbulence equations. The principles of the method are illustrated for three different k-e models. The relaxation method can be used in multigrid form. Multigrid results are given for one turbulence model on a fiat plate test case. Keywords: Navier-Stokes equations; Turbulence; Relaxation methods; Multigrid methods
1. Introduction Basically, two classes of methods for the turbulent compressible Navier-Stokes equations are used nowadays. A first class of methods consists of the explicit time stepping schemes. In the multigrid context, we cite as examples the method of Dimitriadis and Leschziner [2] using explicit Lax-Wendroff stepping and the low Reynolds number Chien model and the method of Gerolymos [3] using also Lax-Wendroff stepping but the Launder-Sharma model. These methods suffer from a severe time step restriction due to the presence of the source terms in the turbulence equations. This handicaps the optimization of the multigrid method. Properly defined implicit time stepping schemes overcome this difficulty. Examples of methods of implicit type are the methods of Vandromme and Ha Minh [12] and Morrison [7]. The first method uses the central discretization, the last uses the upwind discretization. These methods are not used in multigrid form. Some methods are of hybrid type like the method of Mavriplis and Martinelli [6], where the Navier-Stokes part is solved explicitly and turbulence part implicitly and where the Navier-Stokes equations are discritized in the central way and the turbulence equations in the upwind way. * Corresponding author. E-mail: [email protected]. 0168-9274/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S0168-9274(96)00061-X
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
50
Implicit time stepping methods generally allow very large time steps. This observation indicates that it is possible to develop methods that solve directly the steady equations without the use of time stepping. Such a method is then equivalent to an implicit time stepping method with infinite time step. In this paper we develop a relaxation method for the steady equations and employ it in multigrid form. To ensure solvability by a relaxation method, the system of discrete equations has to be so-called positive. To obtain positiveness, the different parts of the equations have to be treated carefully: the convective part, the diffusive part and the source part. These three parts have to be split into positive and non-positive contributions. The positive contributions form the left-hand side of the system; the non-positive contributions form the right-hand side. So, with adequate treatment of the different terms in the equations, a solution of the steady equations can be obtained without the need of time stepping. In multigrid methods employing the turbulence equations on the coarse grids [2,3,6], often the coarse grid corrections cause the turbulence quantities k and e to become negative. This leads to breakdown of the calculations. So it is necessary to damp the coarse grid corrections. In [3,6], this is done by limiting the coarse grid corrections to some upper values. Here we follow a more sophisticated approach, explained later, developed by Lien and Leschziner [5] for incompressible Navier-Stokes equations. In this paper, we use the low Reynolds number form of the turbulence equations. We consider two variants of the multigrid formulation. In the first variant the k-e equations are not used on the coarse grids. In the second variant, the k-e equations are solved on the coarse grids, but the coarse grid corrections are damped. We compare the efficiency of both methods and come to the conclusion that it is advantageous to employ the turbulence equations on the coarse grids.
2. The equations We consider the set of steady compressible Navier-Stokes equations for turbulent flow, coupled to a k-¢ turbulence model. In two dimensions, the equations take the form: OF
aG
OFv
OG~
where
F =
pu
pv
puu + p + ~pk
puv
pvv + p + 2 pk
puv
p u e + (p +
G =
pk) u
'
p v E + (p + ~pk)v
puk
pvk
puz
pv~
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
51
o
0
0
7-xx
7-zy
0
Txy
~y
F~ =
,
Gv
S =
-~
0
UTxy + VTyy Jr- qy
urzx + vrxy + qx
0
(].t+ ,t/O'k)]¢x (JZq- ,t/O'e)Ex
Sk \ S~ J
The components of the stress tensor T and the heat flux vector q are given by rxx = (# + #t) [2ux - }(Ux + Vy)],
+ v~),
~-~ = ~-~ = (~, + ~ , ) ( ~
qx=
( Iz Idt) ~+Prt
ryy = (# + #t) [2Vy - ~(uz + vy)],
(~
qv=
hx,
lZt~h
~+PrtJ
y"
The molecular viscosity is denoted by # and the turbulent viscostity by/.z t. Pr is the molecular Prandtl number (taken as 0.71 for air) and Prt the turbulent Prandtl number (taken as 0.91 for air). The static enthalpy is denoted by h. E stands for the total internal energy per unit mass: E-
--1
p+
1 u2 + 1
V2+ k .
Source terms of a k-e model can be written as
Sk = Pk - p e -
D,
S~ = [ C e , f l P k - C ~ 2 f 2 P e ] ~
+ E,
with
{ [ou, ek=
~,LOxj + 7x~
~6ij Oxk J - - ~ i j p k
-
~xj' -
(1)
lZt = Cufupk"]-,
where Cs,, C~2, Ci,, ak and o-~ are the standard k-e model constants; T is the turbulence time scale; fl, f2 and fu are the so-called wall proximity damping functions and D and £ are the low Reynolds number terms. In our study, three models are used: Launder-Sharma (LS) [8], Lam-Bremhorst (LB) with the modification of the f2-function [10] and Yang-Shih (YS) [14]. The five basic constants are: CE, = 1.44, C~- 2 = 1.92, ak = 1, as = 1.3, Cu = 0.09. Table 1 gives the low Reynolds terms Table 1 Low Reynolds terms and boundary conditions for e
k-e LS
T k -
ew - B.C.
79
0
2 . \ i~y ]
2/zut \0--~2 }
C
2
$
(a2u
LB
k e
-- = 0
Oe i~y
0
0
YS
_k+~/-~e
2 # \ Oy ]
o
(02u~ 2 ~, \~-Ty~/
52
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
Table 2 Damping functions. RT = kT-/u, R u = v/-ky/v, aj = 1.5 × 10 4 a3 = 5.0 × 10-7, a5 = 1.0 × 10-I° k-e
f~
fl
fe
LS
exp[-3.4/(1 + RT/50) 2]
1
1 -
0.3 exp(-R~,)
LB
[1 - exp(-0.0165Ru)]2(1 + 20.5/RT)
1 + (0.05/f~) 3
1 -
exp (-R~. - 10-1°)
YS
[1 - exp ( - a , n u - a3R~ - asR~)] '/2
1
1 - 0.22 exp (-R2/36)
together with the boundary conditions for these models and Table 2 gives the damping functions. The distance to the wall is denoted by y. Upstream of a leading edge, the wall distance is taken as the distance to the leading edge.
3.
The
discretization
In the next sections, the discretization of the different parts (convection, diffusion and source) will be discussed and their role in the relaxation method will be detailed. 3.1. Convective terms The convective part is treated by a flux-difference splitting method. Here, the polynomial fluxdifference splitting [1] is employed which is a variant of Roe-splitting. Using the polynomial character of the convective fluxes F and G with respect to the vector of primitive variables W T = {p, u, v , p , k, e}, differences of flux vectors can be expanded as:
p 2fz2 + -~k ~
AF=
~
+ fit ~0
0
0
0
0
0
1
2 -~fi
0
p--~
0
0
2 p---~+ -~pu--
0
5~ + ~k~2
a42
puv
0 "7 "7 _ 1 ~t
f~~
-p k
o
o
p--a
0
gft
pe
0
0
0
pu
!
AW = A1AW,
where (1 = ~ - ~ +
2 vz,
a42 = p q + - ~ f t +
'71p+/~k+ "7-
3P--k.
The bar denotes the algebraic mean value. AG is treated similarly. Most of the terms in the above expression are self-evident. Some ambiguity arises in treating the triad terms. The following choices are made. The term p u u in the momentum equation is seen as the product of the mass flux pu with the velocity component u. This leads to A p u u = - ~ A u + ftApu = ~--~Au + ~(~Ap + fiAu). The other
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
53
possible grouping into p and u 2 being a density and twice a kinetic energy has no physical meaning. The energy flux is ( 2 )
1
1
2
1
2
( 2 )
puE + p + -spk u = 7 - 1pu + -~puu + -~puv + puk + p + -~pk u. The term ½PUU2 is seen as the product of a mass flux pu with a kinetic energy ½u2. Similarly the term puk is a product of a mass flux and the turbulence kinetic energy. The term (p + 2pk)u is the product of the effective pressure (molecular plus turbulent pressure) and the velocity u. This means that the groupings in the term ~pku and the convective term puk are different. The above choice of groupings (i.e., following the physical meaning of the terms) leads to the simplest expressions afterwards. The basic splitting is done with respect to the primitive variables W. A transformation to a splitting with respect to the conservative variables U T = {p, pu, pv, pE, pk, pe} is obtained from .
1
AU= q+/¢
g
0
0
0
0
O
/5
0
0
0
0
0
/5
0
0
0
/sfi
/50
1/(7 - 1)
/5
0
o
o
o
p
0
0
0
0
0
/5
AW = TAW.
The expression A F = A'IAW transforms into AF = A'IT-1AU = A1AU. Similarly AG = A'2T-1AU = A2AU. The flux-differences are split into a positive and a negative part according to the sign of the eigenvalues of the matrices A1 and A2:
AF = A+AU + ATAU,
AG = A+AU + A2AU.
We consider a vertex-centred finite volume formulation as shown in Fig. 1. The flux-difference over the surface Si+l/2 with length ASi+l/2, is .T'i+! -- .:Ki = A~i,i+l ---- (nxAFi,i+l + nyAGi,i+l )A8i+1/2 = (nxAl + nyA2)AUi,i+lASi+l/2
+ + Az,i+l)AUi,i+lASi+l/2" ---- Ai,i+lAUi,i+lASi+l/2 = ( A i,i+l Splitting the matrix A into a positive and negative part allows the definition of the absolute value of the flux-difference by + IA-T'i,i+ll = ( A i,i+l -- A~,i+l)AUi,i+lASi+l/2 •
(2)
Based on (2) a first order accurate upwind definition of the flux is
,~7i+1/2 __-- 21 [-~"i -I- -~i+1 -- ]~i,i+1]]
(3)
= d~i + ASi+l/2AT,,i+lAUi,i+l.
Fluxes on the other surfaces of the control volume can be written analogously. Summation over the surfaces expresses the inviscid flux balance of a control volume, where the index k refers to the surrounding nodes:
Z k
A;Ask(Uk-Uo)=0,
Z(-A;Ask)Uij k
- Z-A;AskUk k
=0.
(4)
POLYMER SCIENCE U.S.S.R. This Journal is published twelve times per year by Pergamon Press and is a cover-to-cover translation of Vysokomolekulyarnye soyedineniya (series A) (Russian original copyright: Nauka, 90 Profsoyuznaya ul., Moscow 117485)
Scientific Translation Editor: Dr. W. COOPER 217 Little Aston Road, Aldridge, Walsall W S 9 0 P A , West Midlands
Translators: C. W. CAPP, R. J. A. HENDRY, G. F. MODLEN, A. CROZY M. KUBIN, J. HORSKA, D. DOSKOCILOV/, and Mrs. E. SEMERE
Chairman of the Honorary Editorial Advisory Board: SIR HARRY MELVILLE,
K.C.B., F . R . S .
Honorary Editorial Advisory Board for this English version: G. R. ALLAN, London P. E. M. ALLEN, Adelaide C. H. BAMFORD, Liverpool J. W. BARRETT, Cheltenham J. C. BEVINGTON, Lancaster S. CLAESSON, Uppsala P. CORRADINI, Naples J. B. DONNET, Mulhouse P. J. FLORY, Stanford J. FURUKAWA, Kyoto G. GEE, Manchester
N. L. P. A. A. C. C.
GRASSIE, Glasgow HOLLIDAY, London MEARES, Aberdeen A. MORTON, Cambridge, Mass. M. NORTH, Glasgow G. OVERBERGER, Michigan V. ScstrLZ, Mainz
G . SMETS,
Louvain
H. M. SPURLIN, Wilmington H . L . WILLIAMS,
Toronto
B. H. ZIMM, California
Publishing and Advertising Offices Pergamon Press Ltd., Headington Hill Hall, Oxford OX3 0BW, England Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford, NY 10523, U.S.A.
1986 Subscription Rates For Libraries, university departments, government laboratories and other multiple-reader institutions: U.S. $650.00 per annum (including postage and insurance). Two year subscription rate (1985/1986): U.S. $1235.00. Subscription enquiries from customers in North America should be sent to: Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford, NY 10523, U.S.A., and for the remainder of the world to: Pergamon Press Ltd., Headington Hill Hall, Oxford OX3 0BW, U.K.
Microform Subscriptions and Back Issues Back issues o f all previously published volumes are available in the regular editions and on microfilm and microfiche. Current sub scriptions are available on microfiche simultaneously with the paper edition and on microfilm on completion o f the annual index.
Copyright Q 1986 Pergamon Press Limited U.S. Copyright Law Applicable to Users in the U.S.A. Ph otocopying information for users in the U.S.A.: T h e Item-Fee Code for this publication indicates that authorization to photocopy items for internal or personal use is granted by the copyright holder for libraries and other users registered with the Copyright Clearance Center (CCC) Transactional Reporting Service provided the stated fee for copying beyond that permitted by Section 107 or 108 of the United States Copyright Law is paid. The appropriate remittance of $10.00 per copy per article is paid directly to the Copyright Clearance Center Inc., 21 Congress Street, Salem, MA 01970. The copyright owner's consent does not extend to copying for general distribution, for promotion, for creating new works, or for resale. Specific written permission must be obtained from the publisher for such copying. In case of doubt please contact your nearest Pergamon office. The Item-Fee Code for this publication is: 0032-3950/84 $10.00+.00
P U B L I S H E D BY
PERGAMON PRESS LTD. OXFORD. NEW YORK B E I J I N G • F R A N K F U R T • SAO P A U L O • S Y D N E Y • T O K Y O • T O R O N T O No part of this publication may be reproduced stored in a retrieval system or transmitted in any form or by any means: electronic, electrostatic, magnetic tape, mechanical, photocopying, recording or otherwise without permission in writing from the publishers.
POLYMER SCIENCE: U.S.S.R. 1984
V O L U M E 26
N U M B E R 12
CONTENTS J. VARGA,J. MENCZEL a n d A. SOLTI
2763
B. V. LEBEDEV,T. G. KULAGINA, V. S. SVISTUNOV, V. S. PAPKOV and A. A. ZHDANOV S. V. BRONNIKOV, V. ]. VETTEGREN', L. N. KORZHAVIN and S. YA. FRENKEL' V. A. PoPov, V. V. GUZEYEV, Yo. A. ZVEREVA, A. N. GRISHIN, T. V. PALAYEVA, A. P. SAVEL'EV and S. ~q. POTEPALOVA V. V. SHILOV, T. E. LIPATOVA, V. V. VORONA, V. N. BLIZNYUK and A. I. SNEGIREV A. I. MAKLAKOV, V. A. SEVRYUGIN, V. D. SKmDA and N. F. FATKULLIN V. V. SHILOV, YE. I. ORANSKAYA and Yu. S. LIPATOV N. V. POGODINA, A. B. MEL'NIKOV, O. 1. MIKRYUKOVA, S. A. DIDENKO and G. N~ MARCHENKO A. V. MAKSIMOV, Yu. YA. GOTLIB and V. G. BARANOV A. G. KORCHAGIN, M. A. MARTYNOV and S. A. TSYGANKOV M. M. KOTON, V. V. KUDRYAVTSEV, V. A. ZUBKOV, A. V. YAKIMANSKII, T. K. MELESHKO and N. N. BOGORAD Yu. I. MITCHENKO, T. I. LEBEDEVA and R. F. TSIPERMAN V. P. BUDTOV, YE. L. PONOMAREVA, N. I. KONDRASHKINA and T. N. ZELENKOVA V. I. SKR1PACHEV, V. N. KUZNETSOV and S. S. IVANCHEV N. A. YURENKO, G. T. YEVTUSHENKO, YU. YE. YERMILOVA, 1. M. SHOLOGON and B. A. ROZENBERG M. V. TSEBRENKO and G. P. DAN~LOVA
2773
A. Yu. BILIBIN,A. V. TEN'KOVTSEV, O. N. PIRANER and S. S. SKOROKHODOV L. YA. MADORSKAYA, V. M. SAMOILOV, G. A. OTRADINA, A. P. AGAPITOV, V. P. BUDTOV, T. G. MAKEYENKO, YE. YU. KHARCHEVA and N. N. LOGINOVA T. I. BORISOVA, L. L. BURSHTEIN, T. P. STEPANOVA and V. P. SHInAYEV
2781 2789
Calorimetric study of the crystallization and melting of polypropylene Calorimetric study of 1,1,3,3,5,5-hexaethylcyclotrisiloxane and polydiethylsiloxane and the process of polymerization of 1,1,3,3,5,5-hexaethylcyclotrisiloxane in the range 13-330 K Low temperature dependence of the strength of polymide fibres Polymerization modification of fillers and its influence on the properties of filled polyvinyl chloride
2797
Supramolecular organization of some block copolyurethanes with saccharide units in the main chain
2804
Self-diffusion of macromolecules in polystyrene solutions
2811
Microphase structure of block copolyurethanes with a different length of the rigid chain blocks Conformational characteristics of low-molecular weight collodions as shown by diffusion-sedimentation analysis
2819
2825 2834 2839
Orientational order and relaxation times in two-dimensional models of polymer systems with inter-chain interactions Structure of high-density polyethylene obtained by hydrostatic extrusion in solid state Experimental and theoretical study of the effect of medium on chemical imidization
2848
Molecular interactions in polymeric amido-salt soluttons
2856
The change in the molecular mass distribution of photosensitized low-density polyethylene during atmospheric aging Th~ mechanism of the modifying action of oligomers in filled compositions based on polyethylene Effect of chemical isomerism of glycidyl diphenolic ethers and esters of dicarboxylic acids on dielectric relaxation in epoxy polymers Features of the fracture of ultrathin fibres in extrudates of polymer blends Synthesis of high-molecular weight liquid crystal polyesters based on a polycondensation mesogenic m o n o m e r Structure of the endgroups of polyvinylidene fluoride synthesized in the presence of fl-hydroxyethyl-tert-butyl peroxide
2862 2867
2872 2882 2891
2897
Influence of the chemical structure of side chains of comb-like polymers on the gel point
(Continued inside back cover} I S S N 0032-3950 Printed in Poland (PWN-DUAM)
257
E. Dick, J. Steelant /Applied Numerical Mathematics 23 (1997) 49-61
56
It is straightforward to verify that all eigenvalues are non-negative. The contribution of the differences in the ~ direction in the diffusive flux balance forms a positive system. The resulting positive system in the left-hand side is expressed in the variables {p, u, v, e, k, e}. The transformation to the conservative variables is done by
1
0
0
0
0
O
~t
p
0
0
0
0
0
~
0
0
0
p~
p~
p
p
o
o
o
o
p
o
0
0
0
0
p
AU=
q + f~ + ~
AW = TAW.
To transform we use AW = T-1AU. The transformation does not change the positive character of the left-hand side.
3.3. Source terms Whereas the construction of the convective and diffusive Jacobians is rather straightforward, this is not the case for the source terms. For the source terms, a proper linearization must be chosen. The Jacobian of the negative source term is then to be brought into the left-hand side to increase the diagonal dominance of the system of equations [13]. Dependent on the k-e-model, a different linearization is necessary. Until now we have
[
~-~(-A k + B~,k)Ask]Uij - ~-~(-A-~ + B~,k)AskUk = R2 + SV, k
k
where the sum extends over the surfaces of the control volume, Ask denotes the length of the surface and k refers to a surrounding node. The term R2 collects the second-order contributions from the inviscid terms and the tangential contributions from the viscous terms. V denotes the volume of the control volume. The source term still has to be treated. A typical relaxation scheme used in the sequel is three Gauss-Seidel steps on the left-hand side (inner iteration) between updates of the right-hand side (outer or defect-correction iteration). Put into ~-formulation, we obtain
~ - ~ ( - A ~ + B~,k)AskSUij ~- F1 = R2 d- S V , k
where F1 is the flux balance based on first-order inviscid fluxes and the normal part of the viscous fluxes. This term is partly on an old iteration level and partly on a new iteration level. The coefficients A~- and B~,k are also partly on new and old iteration levels as these are updated during the inner iterations.
57
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
The source term S is split into positive and negative terms. The negative term is put into the left-hand side and takes part in the inner iteration: S = S + + S - + -0-O--~U~j. So we obtain:
{ Z ( - A - ~ + B~,k)Ask- OS-V}
(6)
k
The negative source terms taken into consideration for linearization are the same for all models: 1
S[ = -[ps 4- D],
S~- = -C~2f:ps -~.
We consider first the LS and LB model. Following Vandromme [13], we write these terms as:
I_ #t J
-- -~ pk,
S~-=-[Ce'f2](Pe)2- pk
In particular, in the expression of S~-, pe has been substituted to pk using the expression of #t (1) in order to increase the diagonal dominance. Considering the quantities in square brackets to be constant, a linearization which guarantees positiveness, is: ~(S~-; S~-)
O(pk; pg)
_
OSk ~pk DS~apk
OS-~ Op~ = OS: Opg
2Cuf~ (pk) - 79_ 0 #t k C~2fe ( pe'~ 2 pg \-~ ] -2Ce2f2-~
For the YS model, the same combinations are kept constant, but due to the different expression of 7-, we have
Sk ~__(pk)24 (V~ 4- ~ ) 2 =
S~- : _[C~2f2](pg)___~, = -- + pg
.
(7)
This results in the Jacobian
) ~(pk;p~)
,'y(V/~ 4- ~ ) 2 ( p / ~ ) Ce2f2~ -
0 1
p
_Ce2f2 ( a T _ ~ V ~ )
~
1
The resulting global matrix formed by the left-hand sides of equations of type (6) still is a generalised M-matrix in the sense that by reduction of all matrix coefficients to positive scalars the global matrix becomes an M-matrix. We observe, without rigorous proof, that this generalised system can be solved by any collective relaxation method. The solvability of the system of equations resulting from the discretization of the Navier-Stokes equations by different relaxation methods is analysed by the authors in another paper [11]. Here we illustrate the practical applicability of relaxation methods in multigrid form with the example of Gauss-Seidel relaxation.
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
58
4. ~ s t c ~ e Flow with a zero-pressure gradient was calculated over a flat plate with a freestream turbulence level of 3% (T3A test case from Savill [9]). A stretched grid of 385 × 97 points was used (Fig. 2). The grid extends upstream of the plate, with the sharp leading edge at station 97. The first grid point in the direction normal to the plate lies at about y+ -- yu~/u -- 1, where u~- is the friction velocity. Stretching was applied normal to the plate and in the flow direction near the leading edge. Uniform inlet profiles for total temperature, total pressure, k and e were specified. At inlet, Mach number was extrapolated from the flow field. The values of k and e at the inlet were calculated with the equations for k and e for uniform flow with velocity U:
U ok U oE Ox Ox where at the leading edge the following values were matched to be in accordance with the experiments (for L = lm): =
ke = 0.03(3Ue2),
=
ee = 0.378m2/s 3,
Ve =
5.4m/s.
The upper and right boundaries are outlet boundaries. There, pressure was imposed. Velocity components, temperature and turbulent quantities were extrapolated. The part of the lower boundary upstream of the leading edge was treated as a symmetry line. At the plate, no-slip and adiabatic boundary conditions were imposed. Density and pressure were obtained by characteristic combinations of the equations [ 1]. Three lexicographic Gauss-Seidel relaxations with underrelaxation factor 0.9 were used as inner iteration. The first relaxation starts from the left bottom point and ends at right upper point. The second relaxation has the reversed ordering. The third relaxation has the same ordering as the first one. Fig. 3 shows the obtained distribution of the skin friction coefficient for the three models. The upper and lower lines correspond with fully laminar and fully turbulent flow fields. Transition point i
i
!!,ll~lll]ll]l IIIII:III:I
1
_
] I I I J
1 mra
Fig. 2. Detail o f the grid at the l e a d i n g edge.
E. Dick, J. Steelant /Applied Numerical Mathematics 23 (1997) 49-61
59
0.010.0090.0080.007. 0.006-
~3
0.005 0.004 0.003 0.002 0.001
160
260
360
460
660
600
Rex (Thousands)
Fig. 3. Skin friction coefficient predicted by the models.
and transition length are not well predicted by all models, when compared to experimental results [9]. This shows that the models still have to be much improved. The convergence results are similar for all models and are discussed in the next section.
5. Multigrid formulation A standard multigrid method using four grids (385 x 97; 193 x 49; 97 x 25; 49 x 13), W-cycle, full weighting as restriction for residuals and bilinear interpolation as prolongation, was employed. The multigrid acts on the left-hand side of the set of equations. The right-hand side is updated in a defect correction cycle. The procedure is the same as the one used for the Euler equations in [1], except for the ordering of the relaxation. Three Gauss-Seidel relaxations are used as prerelaxation and postrelaxation. Two variants of the multigrid formulation are considered. In the first, the turbulence equations are not used on the coarse grids. In the second, the turbulence equations are used on all grids. In this last variant, without damping of the coarse grid corrections for the turbulence quantities, negative values for ]C and e are always encountered leading to a breakdown of the calculation. Therefore, we use the following damping [5]: ]chew --
]Cold + ar]c +
]cold -- aS]c-
]Cold,
(8)
with a similar expression for e, where 5]C+ and 5 k - denote the positive and negative coarse grid correction respectively. The damping factor a is less than 1. For a = 0, the first variant is recovered and the solution of the turbulence equations on the coarse grids can be omitted. A relaxation on the current grid is taken as one local work unit. A residual evaluation plus the associated grid transfer is also taken as one local work unit. Local work units are counted as fractions of work units on the level of the finest grid proportional to the number of cells. The update of the right hand side of the system of equations in the defect correction is also taken as one work unit. With
60
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49-61
1 E-09,1,,~ 1E-1( 1 E-11 1E-12 1E-13~t~ 1E-1, 1E.15. t ....
1E-16]0
1'0
o's
e'0
~s
Work Units (thousands)
Fig. 4. Convergence history for single grid and multigrid calculations for the YS model. Logarithm of the residual as function of work units. these rules, the cost of the cycle is found to be 12.04167 for the first variant and 14.0625 work units for the second. In a single grid calculation, the defect correction cycle is counted as 4 work units. The calculation starts from uniform flow. First a laminar solution is calculated up to a sufficient level of convergence. With this solution, initial values of k and e are calculated according to the boundary layer laws [9]: k = ki
u
,
~ = 0.3k
with ~ ~> ~i,
where the subscript i refers to inlet conditions. Fig. 4 shows the convergence behaviour for the YS model starting from the laminar solution and the initial values of k and ~. The residual shown is the maximum residual over all equations and all nodes on the finest grid at the end of the cycle. The residuals are the residuals of the flux balances. These contain the size of the control volumes resulting in small initial values in Fig. 4. The best multigrid performance is obtained with the second variant for the largest possible value of the damping factor. It was found experimentally that above some critical value of this damping factor, depending on the flow problem, breakdown of the convergence occurs. The best damping value is just below this critical value. For the flat plate problem the critical value is slightly above 0.3. Without solving the turbulence equations on the coarse grids, the convergence is less good, although the cycle is less expensive. For o~ = 0 and 0. l, off-levelling of the convergence is seen. For comparison, in Fig. 4, also the convergence behaviour of the single grid calculation is shown. Although its asymptotic convergence still cannot be seen in the figure, the big gain in convergence speed due to multigrid is clear.
6. Conclusion A relaxation method for the steady turbulent compressible Navier-Stokes equations coupled to the low Reynolds number k-e turbulence equations was developed. Therefore, a careful treatment of the
E. Dick, J. Steelant / Applied Numerical Mathematics 23 (1997) 49~51
61
convective, diffusive and source terms was necessary. The relaxation method was used in multigrid form. The turbulence equations cannot be used directly in the multigrid formulation. Damping of the coarse grid corrections for the turbulence quantities is necessary. A big gain of convergence speed in comparison with the single grid calculation is obtained.
Acknowledgements The research reported here was granted under contract 9.0001.91 by the Belgian National Science Foundation (N.EW.O.) and under contract IUAP/17 by the Federal Services of Scientific, Technical and Cultural Affairs (D.W.T.C.).
References [ 1] E. Dick, Multigrid solution of steady Euler equations based on polynomial flux-difference splitting, Internat. J. Numer. Methods Heat Fluid Flow 1 (1991) 51-62. [2] K.E Dimitriadis and M.A. Leschziner, Multilevel convergence acceleration for viscous and turbulent transonic flows computed with a cell-vertex method, in: J. Mandel et al., eds., Proceedings 4th Copper Mountain Conference on Multigrid Methods (SIAM, Philadelphia, PA, 1989) 130-148. [3] G.A. Gerolymos, Implicit multiple-grid solution of the compressible Navier-Stokes equations using k-e turbulence closure, AIAA J. 28 (1990) 1707-1717. [4] D.C. Jespersen, A multigrid method for the Euler equations, AIAA-83-0124 (1983). [5] ES. Lien and M.A. Leschziner, Multigrid acceleration for recirculating and turbulent flows computed with a non-orthogonal, collocated finite-volume scheme, Comput. Methods Appl. Mech. Engrg. 118 (1994) 351371. [6] D.J. Mavriplis and L. Martinelli, Multigrid solution of compressible turbulent flow on unstructured meshes using a two-equation model, lnternat. J. Numer. Methods Fluids 18 (1994) 887-914. [7] J.H. Morrison, A compressible Navier-Stokes solver with two-equation and Reynolds stress turbulence closure models, NASA Contractor Report 4440 (1992). [8] V.C. Patel, W. Rodi and G. Scheuerer, Turbulence models for near-wall and low Reynolds number flows: a review, AIAA J. 23 (1984) 1308-1319. [9] A.M. Savill, A synthesis of T3 test case predictions, in: O. Pironneau et al., eds., Numerical Simulation of Unsteady Flows and Transition to Turbulence (Cambridge University Press, Cambridge, 1992) 404442. [10] K. Sieger, A. Schulz, M.E. Crawford and S. Wittig, Comparitive study of low-Reynolds number he turbulence models for predicting heat transfer along turbine blades with transition, in: Proceedings International Symposium on Heat Transfer in Turbomachiner3,, Athens (1992). [11] J. Steelant, S. Pattijn and E. Dick, Smoothing properties of different preconditioners for Navier-Stokes equations, in: Proceedings 3rd Eccomas CFD Conference, Paris (1996). [12] D. Vandromme and H. Ha Minh, About the coupling of turbulence closure models with averaged NavierStokes equations, J. Comput. Phys. 65 (1986) 386-409. [13] D. Vandromme, Turbulence modeling for compressible flows and implementation in Navier-Stokes solvers~ VKI-LS 1991-02. [14] Z. Yang and T.H. Shih, New time scale based k-e model for near-wall turbulence, AIAA J. 31 (1993)
1191-1198.