A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids

Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

www.elsevier.com/locate/cma

A high-order characteristics upwind FV method for incompressible ¯ow and heat transfer simulation on unstructured grids Yong Zhao *, Baili Zhang School of Mechanical and Production Engineering, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Singapore Received 14 August 1999

Abstract This paper presents a new unstructured-grid upwind ®nite-volume algorithm for accurate numerical simulation of incompressible ¯ows and convection heat transfer on unstructured grids. It is an upwind method at both the di€erential equation level and discretized equation level, based on the method of characteristics. This is made possible with the introduction of Chorin's (J. Comput. Phys. 2 (1967) 12±26), arti®cial compressibility formulation to the incompressible ¯ow equations. Flow variables are calculated along characteristics and their initial values are interpolated based on the signs of the corresponding characteristic speed. In addition, an upwindbiased interpolation method of third-order accuracy is used for interpolating ¯ow variables on unstructured grids. With these inherent upwinding techniques for evaluating convection ¯uxes at control volume surfaces, no arti®cial viscosity is required. And the method is accurate and less sensitive to the grid orientation, thus suitable for unstructured grid application. The use of the ®nite-volume method combined with unstructured grids makes the scheme very ¯exible in dealing with very complex boundary geometries, while remaining to be fully conservative at the cell as well as global level. The discretized equations are solved by an explicit multistage Runge±Kutta time stepping scheme which is found to be ecient in terms of CPU and memory overheads. A computer code has been developed which uses the numerical methods presented, and a highly ecient edge-based data structure for ¯ux calculation and data storage. A number of well-documented test cases, including 2D laminar ¯ow in a channel with a backward-facing step, forced and natural convection ¯ows, have been calculated in order to validate the code and evaluate the performance of the numerical algorithm. Good agreement between the computed and measured results have been obtained for all the cases. To demonstrate the accuracy of the thirdorder scheme, ¯ow in another channel with a backward-facing step studied by Gartling (Int. J. Numer. Meth. Fluids (1990) 11) is simulated with both the third-order and ®rst-order schemes. It is shown that the third-order scheme can produce a velocity pro®le across the channel that has perfect agreement with that of Gartling (local citation), whilst the ®rst-order scheme cannot even predict the separation near the upper wall. The proposed numerical methods are also used to calculate 3D incompressible inviscid ¯ows, and results show that the method is accurate and robust for general 3D simulation. The 3D Navier±Stokes code is later used to study a liddriven cavity ¯ow at a Reynolds number of 400. The unstructured grid used has 97,336 nodes and 455,625 tetrahedral elements with strong grid clustering near the walls. The convergence rate is found to be satisfactory, which shows that the method is also robust for 3D Navier±Stokes simulation. The results obtained agree well with those reported by Jiang et al. (Comput. Meth. Appl. Mech. Eng. 114 (109) (1994)). Results obtained using the second-order central scheme with arti®cial viscosity are also presented for comparison. It is observed that the third-order characteristics scheme is more accurate than the latter, thus con®rming its third-order accuracy. Ó 2000 Elsevier Science S.A. All rights reserved.

1. Introduction In recent years, considerable advancement has been achieved in the development of numerical methods to predict ¯ows around complex geometries on unstructured grids using the ®nite-volume method [4±6]. *

Corresponding author. Tel.: +79-91467; fax: +79-11859. E-mail address: [email protected] (Y. Zhao).

0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 4 4 3 - 0

734

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Many ®nite-volume or ®nite-di€erence incompressible ¯ow solvers are predominantly based on the pressure correction scheme on structured grids [7]. The arti®cial compressibility method of Chorin [1] can be considered as an alternative among others. In this method, a pseudo-time derivative of pressure is added to the continuity equation, thus turning the system of governing equations into hyperbolic ones (if viscous terms are excluded). Therefore the pressure ®eld is closely coupled with the velocity ®eld at the same time level. Using time-marching, the arti®cial waves created are used for propagating information throughout the solution domain and driving the divergence of velocity to zero. This approach has been used successfully by some researchers [8±15]. Among them, Farmer et al. [12] and Schulz et al. [6] have used a central scheme with Jameson's arti®cial viscosity to prevent the pressure±velocity decoupling, with the former group using arti®cial compressibility method and latter group using a pressure-based method. On the other hand, Rogers et al. [9], Liu et al. [15] and Whit®eld et al. [14] used a Roe-type upwinding scheme with a simple averaging of the left and right variables on a body-®tted structured grid system. The Jacobian matrix, left and right eigenvectors as well as the associated metric coecients are relatively complex for 3D unstructured-grid solvers [4]. In this study, a new upwind method for calculating incompressible ¯ows is developed based on the method of characteristics. The upwind bias is introduced at both the di€erential equation level and the discretized equation level. Thus it is less sensitive to grid orientation, making it more suitable for unstructured grid application. Traditionally the method of characteristics is used primarily for compressible ¯ow calculations [16]. However, with the introduction of the pseudo-time derivative in Chorin's formulation, it is possible that the incompressible equations can also be solved by a similar method of characteristics. This idea was tested in [17] for the solution of 2D incompressible ¯ow equations cast in a curvelinear coordinate system on structured grids. Based on this idea a high-order upwind characteristics method has been developed for numerical solution of the incompressible ¯ow Navier±Stokes equations and the energy equation on unstructured triangular (2D) and tetrahedral (3D) grids for convection heat transfer simulation. It is shown that this method is mathematically simpler than the Roe-type method and can be readily applied to arbitrary 3D unstructured grids in the ®nite-volume framework together with high-order schemes such as the upwind-biased MUSCL interpolation. With these inherent upwinding techniques for evaluating convection ¯uxes at control volume surfaces, no arti®cial viscosity is required. And the method is accurate and less sensitive to the grid orientation, thus making it suitable for unstructured grid applications. The use of the ®nite-volume method combined with unstructured grids makes the scheme very ¯exible in dealing with very complex boundary geometries, while maintaining fully conservative properties at both the cell level as well as the global level. The discretized equations are solved by an explicit multistage Runge±Kutta time stepping scheme which is found to be ecient in terms of CPU and memory overheads. A computer code has been developed which uses the numerical methods presented, and a highly ecient edge-based data structure for ¯ux calculation and data storage. A number of well-documented test cases, including 2D laminar ¯ow in a channel with a backward-facing step, forced and natural convection ¯ows, have been calculated in order to validate the code and evaluate the performance of the numerical algorithm. Very good agreement between the computed and measured results has been obtained for all the cases. The proposed numerical method was also used to calculate 3D incompressible inviscid ¯ows, and preliminary results show that the method is accurate and robust for general 3D simulation. Further validation of the 3D Navier±Stokes code with heat transfer is underway and results will be presented in future.

2. Governing equations The Navier±Stokes and energy equations for incompressible ¯ows are derived from the conservation laws of mass, momentum and energy. The equations modi®ed by the arti®cial compressibility method can be expressed as oW ‡r
~ Fc ˆ r
~ Fv ‡ S; ot

…1†

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

where 2

3 p ~ 5; W ˆ 4U T

2

3 ~ bU 6 ~U ~ ‡ p dij 7 ~ Fc ˆ 4U 5; q ~ TU

2

3 0 ~ F v ˆ 4 sij 5; ~ q

735

2

3 0 S ˆ 4 bex DT   T =Fr 5: 0

~ the velocity vector; ~ F v are the convective ¯ux and Here W is the vector of dependent variables and U F c and ~ viscous ¯ux vectors. b is a constant called arti®cial compressibility whose value a€ects the solution convergence to steady state. bex is the thermal volume expansion coecient and DT  is the reference temper 2 † =gl . The viscous stress tensor and heat ¯ux are given in the ature di€erence. Fr is the Froude number …U1 following:   1 oui ouj ‡ ; …2† sij ˆ Re oxj oxi ~ qˆ

1 rT ; Re Pr

…3†

where Re is the Reynolds number and Pr the Prandtl's number. The above equations are actually written in non-dimensional form based on the following scalings:        x y z t u v w ; ; ; …u; v; w† ˆ ; ; ; t ˆ ; …x; y; z† ˆ   U U l l l l =U1 U1 1 1   p ÿ p1 T  ÿ Tref ; T ˆ ; …4† pˆ  ÿ T  †2 T1 q1 …U1 ref where l denotes a reference length; superscript  indicates dimensional quantities and subscript 1 indicates free stream conditions. 3. Solution algorithm 3.1. Finite-volume method The equations in (1) are discretized on an unstructured tetrahedral grid. A cell-vertex scheme is adopted here, i.e., all computed variables in vector W are stored at vertices of the tetrahedral cells. For every vertex, a control volume is constructed using the median dual of the tetrahedral grid. A 2D example is shown in Fig. 1.

Fig. 1. A control volume surrounding a node P .

736

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Spatial discretization is performed by using the integral form of the NS equations over the control volume surrounding node or vertex p Z Z Z Z Z Z Z Z Z Z Z Z o ~ ~ Wp dv ‡ r
Fc dv ˆ r
Fv dv ‡ S dv: …5† ot cv cv cv cv The convection term is transformed in order to introduce the upwind scheme using an edge-based procedure I Z Z Z nedge X ~ Fc
~ r
~ Fc dv ˆ n ds ˆ ‰…~ Fc †ij
~ nDSŠn : cv

Scv

nˆ1

Here DSn is the part of control volume surface associated with edge n (similar to edge PA as shown in Fig. 2). Therefore all the ¯uxes are calculated for all the edges and are then collected at the two ends of each edge for updating of ¯ow variables in time-marching. This cell-vertex approach is more ecient than the cell-centre method, especially in 3D cases, which calculates ¯uxes through cell faces and stores them at cell centres. This is because the number of faces in a 3D tetrahedral grid is about twice that of the edges, whilst the number of cells is about ®ve to six times that of the vertices. The viscous term is calculated by a cell-based method I Z Z Z ncell X ~ ~ Fv
d~ r
Fv dv ˆ sˆ …~ Fv
D~ Sc † i : cv

Scv

i

DSci is the part of control volume surface in cell i (represented by Cij in Fig. 1). By using the following relation I d~ S ˆ 0; …6† s

the total vector surface of the control volume in a cell i (denoted by subscript ci in the following) becomes Spi : D~ Sci ˆ 1=3D~ Thus the calculation of viscous terms can be simpli®ed as I ncell ncell X 1X ~ Fv
d~ sˆ …~ Fv
D~ Sc † i ˆ …~ Fv
D~ Sp † i : 3 i Scv i Here DSpi is the surface vector of the face opposite node p of the tetrahedron under consideration. And ncell is the number of cells surrounding node P . Here the …~ Fv †i is calculated at the centre of the tetrahedron with a vertex p, which can be obtained by using the Green's Theorem based on the variables at the four vertices of

Fig. 2. Control volume surface and the associated edge.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

737

the tetrahedron. Similar to the Galerkin type of formulation, the gradient of a ¯ow variable /, at the centre of a tetrahedron is P4 P Si / 9~ Si 1 4iˆ1 /i~ ˆÿ ; grad/c ˆ ÿ iˆ1 i 27V V 3 where /i is the ¯ow variable at a vertex i of the tetrahedron and ~ Si is the surface vector that is opposite to node i. V is the volume of the tetrahedron. Gradients at vertices are obtained by a volume averaging of the gradients at centres of cells associated with the vertex under consideration, as shown in [18]. 3.2. Upwind-biased interpolation Here an edge-based method of calculating the total inviscid ¯ux is adopted by calculating and storing the ¯ux integrals based on the edges. Such a treatment leads to higher eciency in computation and reduced data storage requirement. First of all, the left and right state vectors WL and WR at a control volume surface are evaluated using a nominally third-order upwind-biased interpolation scheme, as proposed by Lohner et al. for compressible ¯ows [19]. If the left and right state vectors are set to Wi and Wj (i and j corresponding to the two ends of an edge, for example P and A as shown in Fig. 2), it is a ®rst-order upwind scheme. 1 ‡ WL ˆ Wi ‡ ‰…1 ÿ j†Dÿ i ‡ …1 ‡ j†Di Š; 4 1 ÿ WR ˆ Wj ÿ ‰…1 ÿ j†D‡ j ‡ …1 ‡ j†Dj Š: 4 Here D‡ i ˆ Wj ÿ Wi ; ! ! ‡ Dÿ i ˆ Wi ÿ Wiÿ1 ˆ 2 ij
rWi ÿ …Wj ÿ Wi † ˆ 2 ij
rWi ÿ Di ; ‡ Dÿ j ˆ Wj ÿ Wi ˆ Di ;

! ! ÿ D‡ j ˆ Wj‡1 ÿ Wj ˆ 2 ij
rWj ÿ …Wj ÿ Wi † ˆ 2 ij
rWj ÿ Dj : Therefore 1 ! WL ˆ Wi ‡ ‰…1 ÿ j† ij
rWi ‡ jD‡ i Š; 2 1 ! WR ˆ Wj ÿ ‰…1 ÿ j† ij
rWj ‡ jDÿ j Š: 2 Here j is set to 1/3 which corresponds to a nominally third-order accurate scheme. 3.3. Upwind characteristics method The Euler equations modi®ed by Chorin's method are rewritten in partial di€erential form in a Cartesian coordinate system for the derivation of the method of characteristics: op ouj ‡b ˆ 0; ot oxj oui oui ouj op ‡ ui ‡ ˆ 0; ‡ uj ot oxj oxj oxi

738

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

oT oT ‡ uj ˆ 0: ot oxj Here subscripts i and j equal 1, 2 and 3, representing the three coordinates. Suppose that n is a new coordinate normal to the surface of a control volume that surrounds a particular edge (outward normal direction). In order to extend the method of characteristics to the 3D unstructured grid solver, it is assumed that ¯ow in the n direction is approximately 1D and the above equations can then be transformed into op ouj n ˆ 0; ‡b on xj ot

…7†

oui oui ouj op ‡ uj n ‡ ui n ‡ n ˆ 0; ot on xj on xj on xi

…8†

oT oT ‡ uj ˆ 0: ot onnxj

…9†

Here nxi ˆ

on oxi

and

nxj ˆ

on : oxj

In the t±n space as shown in Fig. 3, ¯ow variables W at time level n ‡ 1 can be calculated along characteristics k using a Taylor series expansion and the initial value at time level (n W k ): W ˆ W k ‡ Wn nt Dt ‡ Wt Dt; and Wt ˆ

W ÿWk ÿ Wn nt : Dt

A wave speed kk is introduced p nt ˆ kk nxi nxi ; And the normal vector components are nxj : nxj ˆ p nxi nxi

Fig. 3. t±n coordinates.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

739

Substituting components of Wt into Eqs. (7) and (8), we have: 1 …p ÿ pk † p ÿ pn kk ‡ b…un nx ‡ vn ny ‡ wn nz † ˆ 0; Dt nxi nxi

…10†

1 …u ÿ uk † p ‡ un …k0 ÿ kk † ‡ u…un nx ‡ vn ny ‡ wn nz † ‡ pn nx ˆ 0; Dt nxi nxi

…11†

1 …v ÿ vk † p ‡ vn …k0 ÿ kk † ‡ v…un nx ‡ vn ny ‡ wn nz † ‡ pn ny ˆ 0; Dt nxj nxj

…12†

1 …w ÿ wk † p ‡ wn …k0 ÿ kk † ‡ w…un nx ‡ vn ny ‡ wn nz † ‡ pn nz ˆ 0; Dt nxj nxj

…13†

1 …T ÿ T k † p ‡ Tn …k0 ÿ kk † ˆ 0; Dt nxj nxj

…14†

where k0 is the contra-variant velocity k0 ˆ unx ‡ vny ‡ wnz : In order to derive the compatibility equations, the spatial derivatives, such as un , vn ,wn and pn , have to be eliminated from the above equations. Following the approach of Eberle [16] for compressible ¯ow equations, each of the above four equations is multiplied by an arbitrary variable and all the resulting equations are summed to form a new equation: 1 p A ‡ pn B ‡ un C ‡ vn D ‡ wn E ‡ Tn F ˆ 0; Dt nxj nxj

…15†

where A ˆ a…p ÿ pk † ‡ b…u ÿ uk † ‡ c…v ÿ vk † ‡ d…w ÿ wk † ‡ e…T ÿ T k †; B ˆ ÿakk ‡ bnx ‡ cny ‡ dnz ; C ˆ anx b ‡ b…k0 ÿ kk ‡ unx † ‡ cvnx ‡ dwnx ; D ˆ any b ‡ buny ‡ c…k0 ÿ kk ‡ vny † ‡ dwny ; E ˆ anz b ‡ bunz ‡ cvnz ‡ d…k0 ÿ kk ‡ wnz †; F ˆ e…k0 ÿ kk † and a; b; c; d and e are the arbitrary variables used to multiply the equations. We de®ne the coecients of the partial space derivatives to be zero, i.e., A; B; C; D and E are zero. The following equations are thus obtained: A ˆ a…p ÿ pk † ‡ b…u ÿ uk † ‡ c…v ÿ vk † ‡ d…w ÿ wk † ‡ e…T ÿ T k † ˆ 0; B ˆ 0;

…16† …17†

C ˆ 0; D ˆ 0;

…18† …19†

E ˆ 0;

…20†

F ˆ 0:

…21†

Eqs. (16)±(21) can be represented as a linear system UX ˆ 0 with X ˆ fa; b; c; d; eg.

740

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Variables a, b, c, d and e are generally non-zero, thus the system of equations has non-trivial solution. This means that det(U) ˆ 0 and we can derive the following eigenvalues based on this condition: k1 ˆ k2 ˆ k5 ˆ k0 ; k 3 ˆ k1 ˆ k0 ‡

q 2 …k0 † ‡ b ˆ k0 ‡ c;

k 4 ˆ k2 ˆ k0 ÿ

q 2 …k0 † ‡ b ˆ k0 ÿ c:

For each eigenvalue or characteristics speed, characteristic equations can be derived from Eqs. (16)±(21). For example, for kk ˆ k0 we have aˆ

bnx ‡ cny ‡ dnz : k0

Substituting this into Eq. (16), we obtain bnx ‡ cny ‡ dnz …p ÿ p0 † ‡ b…u ÿ u0 † ‡ c…v ÿ v0 † ‡ d…w ÿ w0 † ‡ e…T ÿ T 0 † ˆ 0; k0 i.e. b‰nx …p ÿ p0 † ‡ k0 …u ÿ u0 †Š ‡ c‰ny …p ÿ p0 † ‡ k0 …v ÿ v0 †Š ‡ d‰nz …p ÿ p0 † ‡ k0 …w ÿ w0 †Š ‡ ek0 ‰T ÿ T 0 Š ˆ 0: For any b; c; d and e the above equation is always satis®ed. Therefore all the terms in square brackets are zero. As a result, we have …u ÿ u0 †ny ÿ …v ÿ v0 †nx ˆ 0; …v ÿ v0 †nz ÿ …w ÿ w0 †ny ˆ 0; …u ÿ u0 †nz ÿ …w ÿ w0 †nx ˆ 0; T ˆ T 0:

…22† …23† …24† …25†

For k ˆ k1 p ÿ p1 ˆ ÿk1 ‰…u ÿ u1 †nx ‡ …v ÿ v1 †ny ‡ …w ÿ w1 †nz Š:

…26†

For k ˆ k2 p ÿ p2 ˆ ÿk2 ‰…u ÿ u2 †nx ‡ …v ÿ v2 †ny ‡ …w ÿ w2 †nz Š: Finally u; v; w and p are determined using the above characteristics Eqs. (22)±(27): u ˆ fnx ‡ u0 …n2y ‡ n2z † ÿ v0 nx ny ÿ w0 nx nz ; v ˆ fny ‡ v0 …n2x ‡ n2z † ÿ w0 ny nz ÿ u0 ny nx ; w ˆ fnz ‡ w0 …n2x ‡ n2y † ÿ u0 nz nx ÿ v0 nz ny ; p ˆ p1 ÿ k1 ‰nx …u ÿ u1 † ‡ ny …v ÿ v1 † ‡ nz …w ÿ w1 †Š or p ˆ p2 ÿ k2 ‰nx …u ÿ u2 † ‡ ny …v ÿ v2 † ‡ nz …w ÿ w2 †Š;

…27†

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

where cˆ

q 2 …k0 † ‡ b;

f ˆ

741

1 ‰…p ÿ p0 † ‡ nx …k1 u1 ÿ k2 u2 † ‡ ny …k1 v1 ÿ k2 v2 † ‡ nz …k1 w1 ÿ k2 w2 †Š: 2c

Flow quantities at …n ‡ 1† time level obtained from the above equations on the characteristics are then used to calculate convection ¯uxes at the control volume interface, whilst those on di€erent characteristics at …n† time level are approximately evaluated by an upwind scheme using the signs of the characteristics: 1 W j ˆ ‰…1 ‡ sign…kj ††WL ‡ …1 ÿ sign…kj ††WR Š; 2 and WL and WR are obtained by the upwind-biased interpolation. 3.4. Time integration Finally for a certain node p, Eq. (5) is discretized as o Wp DVcv ˆ R…Wp †; ot

…28†

and R…Wp † ˆ ÿ

nedge X

1 ‰…~ Fc †ij
~ nDSŠn ‡ 3 nˆ1

nbcell X iˆ1

nDSc †i ‡ …St †p DVcv : …~ Fv~

…29†

The spatially discretized equations form a set of coupled ordinary di€erential equations which are integrated in time by means of an explicit multistage Runge±Kutta scheme. A general m-stage scheme to advance the solution from the nth time step to the next step can be written as W …0† ˆ W n ; W …1† ˆ W …0† ‡ a1




Dt …0† R …W †; DVcv

Dt …mÿ2† R …W †; DVcv Dt …mÿ1† ˆ W …0† ‡ am R …W †; DVcv ˆ W …m† :

W …mÿ1† ˆ W …0† ‡ amÿ1 W …m† W …n‡1†

For a ®ve-stage scheme, the stage coecients are a1 ˆ 1=4;

a2 ˆ 1=6;

a3 ˆ 3=8;

a4 ˆ 1=2;

a5 ˆ 1:

In the present code, the viscous ¯uxes are updated at every time step. As an explicit scheme, the maximum time step is determined by a stability criterion of the numerical method. Although it is impossible to derive an analytical stability condition for general 3D problems, the Courant±Friedrichs±Lewy (CFL) condition for hyperbolic equations can be used to derive such a condition by analysing a linearized model equation. This condition can then be extended approximately to general 3D problems. For Navier±Stokes computations using an explicit multistage Runge±Kutta scheme, the local time step size is estimated as Dt ˆ CFL

Dl ; max…kj †

j ˆ 0; 1; 2:

Here Dl length scale associated with a node under consideration. Normally it is taken as the smallest height of all the cells sharing the node.

742

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

3.5. Implicit residual smoothing The idea behind this technique is to replace the residual at one point of the ¯ow ®eld with a smoothed or weighted average of the residuals at the neighbouring points. The averaged residuals are calculated implicitly in order to increase the maximum CFL number, thus increasing the convergence rate. Normally this procedure allows the CFL number to be increased by a factor of 2 or 3. The smoothing equation for a node k can be expressed as Rk ˆ Rk ‡ d2 Rk ;

…30†

where R is the smoothed residual, R the original residual and  is the smoothing coecient ( " # ) 2 1 CFL ÿ 1 ; 0:0 :  ˆ max 4 CFL Here CFL is the maximum CFL number of the basic scheme. The solution to the above equations can be obtained on an unstructured grid by a Jacobi iterative method: …m†

Rk i.e.,

…m† Rk

…o†

numnod…k† X

ˆ Rk ‡ 

iˆ1

…m†

‰Ri

…m†

ÿ Rk Š;

Pnumnod…k† …mÿ1;m† …o† Rk ‡  iˆ1 Ri : ˆ 1 ‡ 
numnod…k†

Here the numnod…k† is the number of neighbouring nodes of node k. 3.6. Boundary conditions At the solid wall, slip (for inviscid ¯ows) or no-slip (for viscous ¯ows) and no-injection boundary conditions are imposed, i.e., zero normal ¯uxes of mass, momentum and energy are imposed for adiabatic wall. If there is heat transfer between the ¯uid and the wall, a constant wall temperature or constant heat ¯ux is normally used as the boundary condition for the energy equation. In addition, pressure on the wall is directly calculated from the modi®ed continuity equation. For far ®eld, or upstream boundary inlet ¯ow, velocity is given directly while the pressure is also directly calculated. For downstream boundary, pressure is ®xed and other variables have zero gradients.

Fig. 4. Numerical grids used for the two channels.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

3.7. Numerical grid generation The grid generation process can be divided into the following subprocesses: 1. geometry de®nition; 2. surface grid and initial volume grid generation; 3. point insertion;

Fig. 5. Comparison of predicted and measured streamlines for H =h ˆ 1:5.

743

744

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

4. clustering of grid points; 5. smoothing of grid. Here the geometry de®nition is performed by a commercial CAD system and the surface and initial volume grid are generated by using an edge swapping scheme [20]. It is necessary to further modify the initial grids. This will involve point insertion, clustering and smoothing of grids, all based on the edge swapping method. The clustering of grid points for viscous ¯ow simulation is performed based on the distance of a cell to the solid wall (i.e., the shortest distance to all the wall nodes). A function f is de®ned as an exponential function of the distance d f …d† ˆ

1 ; eb…dÿdd=dd †

here dd is a length scale of the ¯ow ®eld and b is a constant less than unit one. We introduce another function for a cell with an area A (or a volume in 3D) U ˆ Af …d†: Then the average of U over the ®eld and its deviation can be calculated as follows: s PN 2 N 1X iˆ1 …U ÿ Ui † : Uˆ Ui ; rˆ N N iˆ1

Fig. 6. Comparison of predicted and measured streamlines for H =h ˆ 2:0.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Fig. 7. Convergence histories for channel ¯ow simulation …H=h ˆ 1:5†.

Fig. 8. Calculated and measured isotherms around a heated cylinder (forced convection): Re ˆ 150.

745

746

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

For any cells whose U is greater than U ‡ r, a point is inserted to the centre of its circumcircle or circumsphere and connected to the vertices of the cell and then the mesh is reconnected based on the edge swapping algorithm [20]. Similar to the implicit smoothing of residual, a Laplacian smoother is employed to redistribute the grid nodes r…x;~ y; z† ~ rk ˆ ~ rk ; rk ‡ d2~

…31†

r is the smoothed node coordinate vector, ~ r the original vector. and  is the smoothing coecient: where ~ The solution to the above equations can be obtained on an unstructured grid by a Jacobi iterative method: …m†

~ rk

…o†

ˆ~ rk ‡ 

numnod…k† X iˆ1

…m†

ri ‰~

…m†

rk Š; ÿ~

i.e., …m† ~ rk

Pnumnod…k† …mÿ1;m† …o† ~ ~ rk ‡  iˆ1 ri : ˆ 1 ‡ 
numnod…k†

Here the numnod…k† is the number of neighbouring nodes of node k.

Fig. 9. Computational grid for free convection over a single cylinder.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

747

4. Results and discussion 4.1. Flows in channels with a backward-facing step This well-studied test case is chosen to validate the base line characteristics scheme and the codes developed. There are subgroups which have an H =h ratio of 1.5 and 2, respectively. Experimental measurements of the ¯ow ®elds in the two were conducted by Kueny and Binder [21]. The Reynolds numbers considered are 50, 150 and 500 for the two subgroups. Fig. 4 shows the ®nal computational grids used for the two cases. There is a recirculation zone behind the step and the re-attachment point is one of the features of the ¯ow that is dicult to calculate. Here the calculated streamlines and locations of the reattachment point for di€erent Reynolds numbers are presented and compared with the corresponding experimental measurements. Fig. 5 shows the comparison of the calculated and measured re-attachment points for H =h ˆ 1:5, whilst Fig. 6 gives a comparison for H =h ˆ 2:0. Good agreement between the measured and calculated streamlines and locations of re-attachment points is found in the ®gures presented. Fig. 7 shows the convergence histories for ¯ows with Reynolds numbers 50, 100 and 500, with a CFL number 0.9 and b ˆ 1. 4.2. Forced convection over a cylinder This is the ®rst test case used to evaluate the performance of the upwind characteristics scheme in calculating convection heat transfer problems. A 1-in.-diameter horizontal cylinder is cooled in an air

Fig. 10. Comparison of calculated and measured isotherms around a heated cylinder (free convection).

748

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

stream with a Reynolds number of 150. Typical of the ¯ow is a thick boundary layer and asymmetrical isotherms due to convection heat transfer. The ¯ow is also stabilized by the buoyancy force. Experimental results by Eckert and Soehngen [22] using a Mach±Zenhder interferometer are compared with the calculated isotherms in Fig. 8. The measured and calculated isotherms in boundary layer and the wake region agree well with each other. 4.3. Free convection over single and multiple cylinders In the ®rst case, a circular cylinder of diameter 6 cm is heated uniformly 9°C above the ambient air with a Grashof number of 30 000. Fig. 9 shows the computational grid used for this case. An interferogram obtained by Grigull and Hauf [23] by experiment shows a steady laminar plume rising over the cylinder, and it is compared with the calculated isotherms in Fig. 10. Again good agreement is obtained between the measured and calculated results. In case two, three heated cylinders are stacked together and this results in a complex interaction of thermal boundary layers generated by the cylinders. A thermal plume generated by one cylinder can a€ect those coming from below and above and it changes their ¯ow directions.The computational grid is shown in Figs. 11 and 12 the computed isotherms are compared with an interferogram obtained by Eckert and Soehngen [24]. It is noted that there is a close qualitative agreement between them even in such a complex ¯ow situation.

Fig. 11. Numerical grid used for calculation of free convection over three cylinders.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Fig. 12. Comparison of calculated and measured isotherms around three heated cylinders (free convection).

Fig. 13. Computational grid in concentric annulus.

749

750

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

4.4. Free convection within two concentric and eccentric cylinders These are internal convection ¯ow test cases where free convection takes place inside a co-centric, upper eccentric and right eccentric annuli. The Rayleigh number is 48 000 and the Prandtl number of 0.706 with air as the working ¯uid. Fig. 13 shows the computational grid used for the concentric annulus with an inner/outer radia ratio of 0.3846. The calculated results, in the form of local equivalent thermal conductivity on cylinder surfaces are given in Fig. 14 whilst the radial temperature distributions are presented in Fig. 15, in comparison with the experimental measurements taken from Kuehn and Goldstein [25]. It is found that the agreement is excellent. In Fig. 16 computed isotherms in annuli with di€erent eccentricities in the upper, lower and right positions are presented and compared with the corresponding interferogram

Fig. 14. Predicted and measured surface equivalent thermal conductivity (concentric annulus: Ra ˆ 48 000†.

Fig. 15. Predicted and measured radial distribution of temperature (concentric annulus: Ra ˆ 48 000†.

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

751

Fig. 16. Comparison of calculated and measured isotherms in various annuli …Ra ˆ 48 000†.

measurements obtained by Kuehn and Goldstein [26]. A satisfactory qualitative agreement is observed. In the upper eccentricity case, no evidence of a well-developed thermal plume above the heated inner cylinder is observed, thus the narrow gap between the cylinders in the upper annulus is conduction dominant. This is di€erent from other cases which show thermal plumes developing above the inner cylinder. Finally, a quantitative comparison of the calculated and measured temperatures in an upper eccentric annulus at di€erent azimuth angles is presented in Fig. 17. At h ˆ 0, the temperature distribution exhibits a near straight line and the agreement is excellent. At other angles, the calculated and measured results generally agree with each other very well except at 120 dec where the calculated result is slightly higher than the measured one in the middle of the annulus.

752

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Fig. 17. Calculated and measured radial temperature distribution at di€erent azimuth angles.

Fig. 18. Comparison of the third-order and ®rst-order schemes in a channel ¯ow with a backward-facing step.

4.5. Comparison of the third-order and ®rst-order schemes In order to evaluate the performance of the third-order upwind scheme, a rigorous investigation is conducted to provide the evidence of high-order accuracy of the proposed scheme. Another channel ¯ow

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

753

with a backward-facing step studied by Gartling [2] is calculated using the two schemes. The ®rst-order scheme is obtained by using the ¯ow properties at the two ends of an edge directly as the left/right states for characteristics calculation. Fig. 18(a) illustrates the geometries of the ¯ow domain and the boundary conditions used. Fig. 18(b) shows a comparison of the velocity distributions across the channel at x=H ˆ 7, predicted by the third-order and ®rst-order schemes, as well as that by Gartling [2]. A very good agreement with that of the third-order and Gartling's is observed, whilst the ®rst-order scheme is completely di€erent from the two results. It is evident that the ®rst-order scheme is so di€usive that it cannot capture the small separation zone existing near the upper wall. 4.6. 3D Flow simulation 4.6.1. Flow around a sphere Inviscid incompressible ¯ow around a sphere is calculated by the solver in order to validate the method developed for 3D ¯ow simulation. Because the ¯ow is symmetrical only one-quarter of the ¯ow ®eld needs to be considered for the sake of simplicity and eciency. The computational grid is shown in Fig. 19 which has 10 452 nodes and 52 729 tetrahedral cells. A comparison of calculated and exact surface pressure coecients along the symmetrical plane is made in Fig. 20 which shows satisfactory agreement. The convergence is found to be very fast with a CFL number of 1.8, as shown in Fig. 21.

Fig. 19. Numerical grid around a sphere for Euler computation.

754

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

Fig. 20. Comparison of calculated and exact wall pressure on symmetrical plane of sphere.

Fig. 21. Convergence history for inviscid ¯ow around a sphere.

4.6.2. 3D Navier±Stokes simulation of lid-driven cavity ¯ow Flow in a cubic cavity driven by a velocity in x direction at its opening (lid) is calculated by the 3D Navier±Stokes code for validation purpose. This case is chosen because the geometry is relatively simple, yet the ¯ow is rather complicated and also well studied over the last three decades. An unstructured mesh with 97 336 nodes and 455 625 cells was generated for the 3D Navier±Stokes simulation. The geometry and grid systems are shown in Fig. 22(a). The convergence history of the lid-driven cavity ¯ow at Re ˆ 400 is

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

755

Fig. 22. Calculation of lid driven cavity ¯ow at Re ˆ 400: nnode ˆ 97 336 and ncell ˆ 455 625.

given in the Fig. 22(b). The convergence rate is found to be satisfactory considering the clustered grid near the walls. The computed velocity of u along the vertical centre line at Re ˆ 400 is compared with calculations by Jiang et al. [3] in Fig. 22(c), which shows that result generated by the present method agrees well with that of Jiang et al. For comparison, result obtained using the second-order central scheme together with Jameson's fourth di€erence arti®cial viscosity (with a coecient of 0.0036) is also presented in the same ®gure. It can be noted that the result of the present third-order upwind scheme is clearly more accurate than that of the central scheme, thus con®rming the higher-order accuracy of the scheme developed. 5. Concluding remarks A new unstructured-grid upwind ®nite-volume algorithm for accurate numerical simulation of incompressible ¯ows and convection heat transfer on unstructured grids has been presented. The use of the ®nite

756

Y. Zhao, B. Zhang / Comput. Methods Appl. Mech. Engrg. 190 (2000) 733±756

volume method combined with unstructured grids makes the scheme very ¯exible in dealing with complex boundary geometries while maintaining mass, momentum and energy conservation at the cell as well as global level. The discretized equations are solved by an explicit multistage Runge±Kutta time stepping scheme which is found to be ecient in terms of CPU and memory overheads. A new upwind characteristics scheme with high-order interpolation is developed as part of the ®nite-volume method for accurate calculation of the convection ¯ux. A general Navier±Stokes code has been developed using the numerical methods presented and has been validated against a number of well-documented 2D and 3D test cases with and without heat transfer. Good agreement between the computed and measured results (or simulation results reported in the literature) has been obtained for all the cases studied. This shows that the numerical method developed is reasonably accurate and robust. Further validation of the scheme using highly complex 3D ¯ow with heat transfer is underway and results will be presented in the future. References [1] A. Chorin, A numberical method for solving incompressible viscous ¯ow problems, J. Comput. Phys. 2 (1967) 12±26. [2] D.K. Gartling, A test problem for out¯ow boundary conditions-¯ow over a backward-facing step, Int. J. Numer. Meth. Fluids (1990) 11. [3] B.N. Jiang, T.L. Lin, L.A. Povinelli, Large-scale computation of incompressible viscous ¯ow by least-squares ®nite element method, Comput. Meth. Appl. Mech. Eng. 114 (109) (1994). [4] W.K. Anderson, R.D. Rausch, D.L. Bonhaus, Implicit/multigrid algorithims for incompressible turbulent ¯ows on unstructured grids, J. Comput. Phys. 128 (1996) 391±408. [5] A. Belov, L. Martinelli, A. Jameson, A new implicit algorithm with multigrid for unsteady incompressible ¯ow calculations, AIAA paper AIAA-95-0049, AIAA, 1995. [6] K.W. Schulz, Y. Kallinderis, Unsteady ¯ow structure interaction for incompressible ¯ows using deformable hybrid grids, J. Comput. Phys. 143 (1998) 569±597. [7] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, 1980. [8] D. Kwak, J.L.C. Chang, S.P. Shanks, S.R. Chakravarthy, A three-dimensional incompressible Navier±Stokes ¯ow solver using primitive variables, AIAA J. 24 (3) (1986) 390±396. [9] S.E. Rogers, D. Kwak, Steady and unsteady solutions of the incompressible Navier±Stokes equations, AIAA J. 29 (1991) 603±610. [10] C. Mercle, M.A. Athavale, Time accurate unsteady incompressible ¯ow algorithims based on arti®cial compressibility, Paper 871137, AIAA, 1987. [11] A. Rizzi, L. Erikson, Computation of inviscid incompressible ¯ow with rotation, J. Fluid Mech. 153 (1985) 275±312. [12] J. Farmer, L. Martinelli, A. Jameson, Fast multigrid method for solving incompressible hydrodynamic problems with free surface, AIAA J. 32 (1994) 1175±1182. [13] E. Dick, J. Linden, A multigrid method for steady incompressible Navier±Stokes equations based on ¯ux di€erence splitting, Int. J. Numer. Meth. Fluids 14 (1992) 1311±1323. [14] D.L. Whit®eld, L.K. Taylor, M. Beddhu, A. Arabshahi, Discretized Newton-Relaxation solution of the three-dimensional unsteady incompressible Naver±Stokes equations, in: D.A. Caughey, M.M. Hafez (Eds.), Frontier of Computational Fluid Dynamics 1994, Wiley, New York, 1994, pp. 576±594. [15] C. Liu, X. Zheng, C.H. Sung, Preconditioned multigrid methods for unsteady incompressible ¯ows, J. Comput. Phys. 139 (1998) 35±57. [16] A. Eberle, 3D Euler calculations using characteristic ¯ux extrapolation, Paper 85-0119. [17] D. Drikakis, P.A. Govatsos, D.E. Papantonis, A characteristic based method for incompressible ¯ows, Int. J. Numer. Meth. Fluids 19 (1994) 667±685. [18] Y. Zhao, D.E. Winterbone, The ®nite volume FLIC method and its stability analysis, Int. J. Mech. Sci. 37 (1995) 1147±1160. [19] R. Lohner, CFD via unstructured grids: trends and applications, in: D.A. Caughey, M.M. Hafez (Eds.), Frontier of Computational Fluid Dynamics 1994, Wiley, New York, 1994, pp. 118±133. [20] T.J. Barth, Aspects of unstructured grids and ®nite-volume solvers for Euler and Navier±Stokes equations. AGARD-VKI Lecture Series AGARD R 787, Advisory Group for Aerospace Research and Development, 1992. [21] K. Morgan, J. Periaux, F. Thomasset (Eds.), Notes on Numerical Fluid Mechanics: Analysis of Laminar Flow Over a Backward Facing Step, vol. 9, Friedr. Vieweg & Sohn, Braunscheweig/Wiesbaden, 1984. [22] V.D. Milton, An Album of Fluid Motion, The Parabolic Press, Stanford, California, 1982, p. 118. [23] V.D. Milton, An Album of Fluid Motion, The Parabolic Press, Stanford, California, 1982, p. 121. [24] V.D. Milton, An Album of Fluid Motion, The Parabolic Press, Stanford, California, 1982, p. 123. [25] T.H. Kuehn, R.J. Goldstein, An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74 (1976) 695±719. [26] T.H. Kuehn, R.J. Goldstein, An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli, ASME J. Heat Transfer 100 (1978) 635±640.