An Artificial Compressibility Method for Incompressible Flows

This article was downloaded by: [Aalto-yliopiston kirjasto] On: 23 September 2014, At: 03:52 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/unhb20

AN ARTIFICIAL COMPRESSIBILITY METHOD FOR INCOMPRESSIBLE FLOWS M. M. Rahman, T. Siikonen Published online: 30 Nov 2010.

To cite this article: M. M. Rahman, T. Siikonen (2001) AN ARTIFICIAL COMPRESSIBILITY METHOD FOR INCOMPRESSIBLE FLOWS, Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology, 40:5, 391-409 To link to this article: http://dx.doi.org/10.1080/104077901753243188

PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever

or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions

Numerical Heat Transfer, Part B, 40: 391± 409, 2001 Copyright # 2001 Taylor & Francis 1040-7790 /01 $12.00 + .00

AN ARTIFICIAL COM PRESSIBILITY M ETHOD FOR INCOM PRESSIBLE FLOWS M. M. Rahman and T. Siikonen

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

Helsinki University of Technology, Department of Mechanica l Engineering, Laboratory of Applied Thermodynamics , SaÈhkoÈmiehentie, Finland An artiŽ cial compressibility method characterized by the pressure-based algorithm is developed on a nonorthogonal collocated grid for incompressible  uid  ow problems, using a cell-centered Ž nite-volume approximation . Unlike the traditional pseudo-compressibilit y concept, the continuity constraint is perturbed by the material derivative of pressure, the physical relevance of which is to invoke matrix preconditionings. The approach provokes density perturbations, assisting the transformation between primitive and conservative variables. To account for the  ow directionality in the upwinding, a rotational matrix is introduced to evaluate the convective  ux. A rational means of reducing excessive numerical dissipation inherent in the pressure–velocity coupling is contrived which has the expedience of greater  exibility and increased accuracy in a way similar to the MUSCL approach . Numerical experiments in reference to a few laminar  ows demonstrate that the overall artifacts expedite enhanced robustness and anticipated oscillation damping properties adhering to the factored pseudo-time integration procedure.

1.

INTRODUCTION

Considerable research has been devoted to emulating the compressibility ideas around density-base d phenomena for low-Mach-number and incompressible ¯ ows [1± 6], and thereby circumventing the computational complexity through a separate computation of the velocity and pressure ® elds. Preconditioning is an eminent means of simulating nearly incompressible ¯ ows with numerical algorithms designed for the compressible code. The fundamental issue herein is to modify the transient behavior of the equations without altering the steady-state solution achieved [3]. Another alternative for extending a compressible ¯ ow code for use at a signi® cantly low Mach number springs from a ® ctitious equation of state or arti® cial compressibility introduced by Chorin [6], entailing essentially the introduction of a ``small’ ’ perturbation in the incompressibility constraint. Customarily, the arti® cial compressibility approach incorporates the addition of a time-dependent pressure term to the continuity equation together with a multiplicative variable designated as arti® cial compressibility. With this arti® cial

Received 29 November 2000; accepted l5 June 2001. Address correspondence to Prof. M. M. Rahman, Helsinki University of Technology, Department of Mechanical Engineering, Laboratory of Applied Thermodynamics, SaÈhkoÈmiehentie 4, FIN-02015 HUT, Finland. E-mail: mizranur.rahman@hut.® 391

392

M. M. RAHMAN AND T. SIIKONEN

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

NOMENCLATURE A; B c CFL F; G Gr L p Pr Q Re S t T u; v U; V

Jacobian matrices in x and y directions arti® cial sound speed Courant number ¯ ux vectors in x and y directions Grashof number characteristic length static pressure Prandtl number source term Reynolds number cell face area time rotational matrix velocity components in x and y directions contravariant or dimensionless velocity components

W W~ x; y b G y l m n r f 8

conservative variable vector primitive variable vector Cartesian coordinates compressibility parameter di€ usion coe cient Viscous stability factor or dimensionless temperature eigenvalue dynamic viscosity kinematic viscosity density scalar transport variable cell volume

Subscript ref reference value

quantity, the resultant system of equations is symmetric hyperbolic for the inviscid terms. Consequently, the system is well posed and e cient numerical algorithms developed for compressible ¯ ows can be used to advance the system in time. However, the analogy is not always a strict equivalence, since the arti® cial term plays a crucial role in the determination of a credible success [7]. In principle, the solution methodology resorting to arti® cial compressibility for viscous incompressible ¯ ows can deliberately be ascribed to the pressure-based method, requiring physically a precluded decoupling of the pressure± velocity ® elds. When using collocated grids, this approach confronts some critical issues that provoke arti® cial damping terms or a special cell-face interpolation strategy to eliminate the occurrence of spurious oscillations. A more rational appliance is to ® nd a mechanism for introducing just enough cell face dissipation to extenuate the destabilizing e€ ect arising from the pressure checkerboarding. Several formulations based on primitive variables have been contrived to surmount the di culties associated with collocated grid arrangements [7± 13]. Nevertheless, a prerequisite to evaluating the convective ¯ ux on a cell face with the approximative Riemann solver [8] is a numerically convenient choice for the arti® cial compressibility factor on which the steady-state computation envisages a progressive dependency. The present solution method advances with recourse to the cell-centered, nonorthogonal , fully collocated ® nite-volume method devised for two-dimensional incompressible ¯ uid ¯ ow problems. Unlike Chorin’s arti® cial density-base d algorithm, the mass continuity is perturbed by a quantity amenable to the material derivative of pressure, remunerating vitally matrix preconditionings accompanied by an implicit coupling of the primitive variables. The scheme approves density preconditionings=perturbations whose participation is exclusively con® ned to the transformation between primitive and conservative variables. To counteract the additional complexity comprising the scaling of numerical dissipation components embedded with the Riemann-based approach, appropriate ¯ ux-di€ erence scheme is

ARTIFICIAL COMPRESSIBILIT Y METHOD

393

pursued without the occurrence of nonphysical pressure modes. The equations of motion are solved with an improved alternating direction implicit (ADI) factorized time integration method [14]. Detailed presentations of the spatial discretization and associated features are provided.

2.

GOVERNING EQUATIONS

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

The two-dimensional Navier-Stokes equations, including an equation for a scalar f, can be written in the following form: qW q(Finv ± Fvis ) q(Ginv ± Gvis ) ‡ ‡ ˆQ qt qx qy

(1)

where W ˆ (r; ru; rv; rf) T and Q represents the source term. The inviscid ¯ uxes are Finv

0

1 ru B ru2 ‡ p C C ˆB @ ruv A ruf

Ginv

0

1 rv B rvu C C ˆB @ rv2 ‡ p A rvf

(2)

Here r is the density, u and v are the Cartesian velocity components, and p represents the pressure. The viscous ¯ uxes are

Fvis

0

1 0 B G (qu=qx) C C ˆB @ G (qv=qx) A G (qf=qx)

Gvis

0

1 0 B G (qu=qy) C C ˆB @ G (qv=qy) A G (qf=qy)

(3)

where G signi® es the exchange coe cient, i.e., either viscosity or conductivity of the ¯ uid.

3.

SPATIAL DISCRETIZATION

A cell-centered ® nite-volume scheme is applied to solve the ¯ ow equations having an integral form d dt

Z

8

W d8 ‡

Z

S

F(W ) ¢ dS ˆ

Z

8

Qd8

(4)

for an arbitrary ® xed region 8 with a boundary S. Performing the integration for a computational cell i provides 8i

dW i X ˆ ± S F^ ‡ 8i Q i dt faces

(5)

M. M. RAHMAN AND T. SIIKONEN

394

where the sum is taken over the faces of the computational cell. Each face has a unit normal vector n de® ned by n ˆ nx i ‡ ny j ˆ

Sx Sy i‡ j S S

(6)

and the corresponding expression for the cell-face ¯ ux becomes

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

F^ ˆ nx F ‡ ny G

(7)

Here F and G are the ¯ uxes de® ned by Eqs. (2) and (3) in the x and y directions, respectively. To account for the directional in¯ uences in the upwinding process, the inviscid ¯ ux on the cell face (i ‡ 21 ) is calculated as F^ i‡1=2

0

1 rU B C ± 1 B rUu ‡ p C ˆ T i‡1=2 @ rUv A rUf i‡1=2

0

1 B0 T ˆB @0 0

0 nx ± ny 0

0 ny nx 0

1 0 0C C 0A 1

(8)

where the contravariant velocity U ˆ unx ‡ vny , and u ˆ unx ‡ vny and v ˆ vnx ± uny are the velocity components, normal and parallel to the cell face. A fully upwinded second-order (FUS) scheme is adopted to approximate the scalar and velocity components. The rotational matrix T [13] transforms the dependent variables to a local coordinate system with the principal direction being perpendicular to the cell face. After multiplying by the rotational matrix, the ¯ ux has a functional form akin to the Cartesian coordinate system and can be split into its corresponding contributions. A second-order central di€ erencing is applied while evaluating the derivatives in the viscous ¯ uxes. The relevant aspects of arti® cial compressibility and the construction of dissipation models for the pressure± velocity coupling are discussed in some detail in subsequent sections.

3.1. Art ificial Com pressibilit y M odel The formulation revives a structure resembling the compressible equations by inserting arti® cial compressibility in the derivative of density with respect to the pressure, eventuating in a hyperbolic system. Within the framework of a two-dimensional analysis, the nonconservativ e form of the mass continuity can be introduced as qr qr qr qu qv ‡u ‡v ‡r ‡ qt qx qy qx qy

ˆ0

(9)

Having relaxed the incompressibility constraint, the continuity equation is perturbed as

ARTIFICIAL COMPRESSIBILIT Y METHOD

qr qp qp qp qu qv ‡u ‡v ‡r ‡ qp qt qx qx qx qy

395

ˆ0

(10)

where r ˆ r( p). Upon rearranging, Eq. (10) yields qp qp qp qu qv ‡u ‡ v ‡ rc2 ‡ qt qx qy qx qy

ˆ0

(11)

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

where c is denominated herein as the arti® cial sound speed=perturbation parameter recognized by 1 qr ˆ c2 qp

(12)

Reverting to the justi® cation of continuity modi® cation, it can be immediately seen that the arti® cial sound speed must be su ciently large to have a signi® cant regularizing e€ ect, and at the same time must be as small as possible to minimizing perturbations on the incompressibility equation. To reconcile this criterion, c is estimated as cˆb

q                                           max (u2 ‡ v2 ); 12 U2ref

(13)

where Uref represents a reference velocity and b is a compressibility parameter. Evidently, c depends on b, in¯ uencing the convergence rate and stability of the solution method. Values of b in the range of 1± 10 are recommended for better convergence to the steady state at which the mass conservation is enforced. Using Eq. (11), the time-dependent linearized inviscid ¯ ow equations can be written in condensed notation as qW~ qW~ qW~ ‡A ‡B ˆ0 qt qx qy

(14)

with the primitive variable vector, W~ ˆ ( p; u; v; f) T and the Jacobian matrices 0

u B1 B B A ˆ Br B @0 0

rc2

0

u

0

0

u

0

0

0

1

C 0C C C C 0A u

0

v

B B0 B B ˆ B1 B @r 0

0

rc2

v

0

0

v

0

0

0

1

C 0C C C 0C A v

(15)

The eigenvalues=characteristic speeds l of A and B are readily obtained as la ˆ u ± c; u ‡ c; u; u

lb ˆ v ± c; v; v ‡ c; v

(16)

The inviscid matrix A can be decomposed in compliance with the primitive variables as

M. M. RAHMAN AND T. SIIKONEN

396

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

A ˆ L la L ± 1

0

rc B±1 L ˆB @ 0 0

rc 1 0 0

0 0 1 0

1

0 0C C 0A 1

L ±1

0 1 B 2rc B B 1 ˆB B 2rc B @ 0 0

1 2 1 2 0 0

An analogous eigensystem for the matrix B can be devised as 0 1 0 1 B 2rc 0 rc 0 rc 0 B B 0 B 0 1 1 0 0C B b ±1 ±1 B C B ˆ Ml M Mˆ@ M ˆ B A 1 ±1 0 1 0 B B 0 0 0 0 1 @ 2rc 0

0

±

0

0 1 0 1 2 0 1 2 0

±

0

1

C C C 0C C (17) C 0A 1 1 0C C 0C C C C 0C A 1

(18)

The characteristic variables are obtained from L ± 1 DW~ ˆ L ± 1 (Dp; Du; Dv; Df) T for the x direction as 1 ( Dp ± rc Du) 2rc 1 w ~2 ˆ ( Dp ‡ rc Du) 2rc 3 ~ ˆ Dv w ~ 4 ˆ Df w w ~1 ˆ

(19)

Conformable characteristic variables can be de® ned for the y direction. A di€ erential=arti® cial change is incorporated with the density to allow for a transformation between conservative and primitive variables, qr Dp Dr ˆ Dp ˆ 2 (20) qp c Using Eq. (20), the changes in the conservative variables can be expressed explicitly in terms of the primitive variables as follows: Dp Dru ˆ r Du ‡ u Dr ˆ r Du ‡ u 2 (21) c Dp Drv ˆ r Dv ‡ v Dr ˆ r Dv ‡ v 2 (22) c Dp Drf ˆ r Df ‡ f Dr ˆ r Df ‡ f 2 (23) c Obviously, Eqs. (20)± (23) provide explicit de® nitions of the ¯ uid constitutive relations by which a transformation between conservative and primitive variables can be made. Therefore, the transformation relation from conservative to primitive variables can be cast in matrix form as 0 1 0 2 10 1 Dp Dr rc 0 0 0 B Du C 1 B ± u 1 0 0 CB Dru C B C B CB C (24) @ Dv A ˆ r @ ± v 0 1 0 A@ Drv A Df Drf ±f 0 0 1

ARTIFICIAL COMPRESSIBILIT Y METHOD

397

where the conservative explicit change DW is calculated from

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

DW ˆ (Dr; Dru; Drv; Drf) T ˆ

Dti Ri 8i

(25)

where the residual Ri is de® ned by the right-hand side of Eq. (5) and the spatially varying time step Dti is evaluated from [13]. This transformation is applied before the implicit treatment of the solution procedure. The implicit stage features the explicit density residual Dr as a consequence of the use of primitive rather than conservative variables. In principle, Dr is convincingly identical to the mass residual induced by the tentative velocity ® elds. The realization that Dr is explored in order to linearize the residuals as well as to predict a compatibility pressure ® eld in collaboration with the perturbation parameter c. After the implicit stage, the conservative variables are retrieved by allowing an arti® cial change in density as referred to Eq. (20), leading to develop the following linearized relations: 01 c2 0 1 B B Dr Bu B Dru C B c2 B C B B CˆB B @ Drv A B v B c2 B Drf @ f c2

0

0

r

0

0

r

0

0

The primitive variables are updated from: W~

n‡1

0

1

C0 Dp 1 C CB C 0C Du C CB B C CB C CB Dv C 0C A C@ C A Df r

1 ˆ W~ n ‡ (r Dp; Dru; Drv; Drf) T r

(26)

(27)

where the superscript (n ‡ 1) denotes the time level (n ‡ 1) Dt. Unlike Eq. (27), the updating of primitive variables from DW~ ˆ (Dp; Du; Dv; Df) T may end up with an exclusion of the convergence acceleration contributed by the arti® cial density change Dp=c2 to the solution method. To this end, it can appropriately be stressed that the continuity equation is perturbed by the quantity 1 qp qp qp ‡u ‡v c2 qt qx qy which is physically conducive to diagonalizing the inviscid Jacobian matrices, subjected presumably to the matrix preconditionings. Alternatively, the preconditioned continuity equation, containing the pressure-convectiv e term V ¢ Hp assists in forming convenient matrix inversions relating to the eigenstructure. However, Eqs. (24) and (25) show that residuals are calculated using the conservative ¯ ow Eq. (1), indicating no involvement of V ¢ Hp in the mass residual, based on which Dp is evaluated. Consequently, the steady-state solution remains una€ ected by the

M. M. RAHMAN AND T. SIIKONEN

398

arti® cial compressibility modi® cations. Essentially, both Dr’ s in Eqs. (24) and (26) can be interpreted as density perturbations=preconditionings to the incompressible limit. As the solution converges to the steady state, perturbations disappear. In particular, the system is pressure based and the arti® cial density change can be construed as mass imbalance which is driven to zero with the iterative convergence.

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

3.2. M om ent um -Based Dissipat ion Schem e This section addresses the interface mass ¯ ux representation related to the suppression of checkerboard pressure distribution re¯ ecting the pressure± velocity decoupling. To prevent the destabilizing in¯ uence emerging from signi® cant pressure variation, the simpli® ed momentum equations are solved for the cell face velocities. Referring to [13], the interpolation technique to curvilinear coordinates can be provided for the cell face (i ‡ 12) as Ui‡1=2 ˆ Ui‡1=2 ‡ DUi‡1=2

ˆ Ui‡1=2 ± D i‡1=2 qp ‡ rU qU ± Dn q

i‡1=2

(28)

with Dt¤ r Dn

D i‡1=2 ˆ Ri‡1=2 ˆ

Dt¤ r Dn

qp ˆ pR ± pL ¤

1:0 ‡ CFL¤

i‡1=2

|R| |U| ‡ e

qp ‡ rU qU ± Dn q

i‡1=2

i‡1=2

qU ˆ UR ± UL

CFL ˆ max(1:0; CFL)

q ˆ qu nx ‡ qv ny

where qu and qv are the appropriate source terms with the pressure gradients extracted, and Dni‡1=2 is the distance between the cell centers in the i direction. The small number e (¹ 10 ± 6 ) is included to prevent the singularity. The overbar indicates the regular results obtained by averaging the values at the two adjacent cells. The state variables (p; U) L and (p; U) R , respectively, can be evaluated on the left and right sides of the cell face using the MUSCL (monotone upstream-centered schemes for conservation laws) approach of Van Leer [15]. Seemingly, the dissipative mechanism herein has its potential ability to a€ ord a further evaluation of the time step Dt¤ from the relation Dn Dt¤

i‡1=2

ˆ max |U| ‡ c;

2n Dn

(29) i‡1=2

where n ˆ m=r. Equation (28) reveals that the cell face velocity is linked directly to the velocity, pressure, and source terms at the adjacent nodal points, leading to a nonlinear cell-face interpolation scheme. The framework guarantees that the damping term DUi‡1=2 does not introduce local extrema into the cell face velocity and could be reducible with the convergence acceleration. Analogous expressions follow for the remaining cell faces.

ARTIFICIAL COMPRESSIBILIT Y METHOD

3.3.

399

Dual Dissipat ion Schem e

The elimination of decoupling for the odd and even points can be accelerated by introducing dissipations in both the cell face velocity and pressure. To facilitate the subsequent development, the calculation commences with resorting to one-dimensional linearizations of (u; p) in x:

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

u~ i‡1=2 ˆ

u ‡ Dx

qu qx

p~ i‡1=2 ˆ

i‡1=2

p ‡ Dx

qp qx

(30) i‡1=2

Ostensibly, the interface velocity and pressure can be regarded as consisting of central di€ erence and di€ usive contributions . Restricting attention to the steadystate and ascribing u=c º 1 (i.e., the maximum crucial limit being considered to support subsonic ¯ ows), Eq. (11) in one dimension reduces to rc

qu qx

i‡1=2

ˆ±

qp qx

(31) i‡1=2

Combining Eqs. (30) and (31) together with an upwind de® nition conventionally demands that 1 u~ i‡1=2 ˆ (ui ‡ ui‡1 ) ± 2

1 2rc

1 p~ i‡1=2 ˆ ( pi ‡ pi‡1 ) ± 2

rc 2

( pR ± pL )

(32)

( uR ± uL )

(33)

i‡1=2

i‡1=2

where the factor (1=2) is included with the dissipation terms so that a full upwinding is achieved. In particular, the framework of Eq. (31) is chosen such as to preserve positivity in the cell face dissipations to evaluate the mass ¯ ux and pressure gradient. Remarkably, as the continuity constraint (rqu=qx ˆ 0) is evoked, the cell face dissipations approach zero with the iterative convergence. Equations (32) and (33) warrant that a strong coupling between the pressure and velocity ® elds is retained. Extension of the dissipative mechanism to curvilinear coordinates results in the following expressions for the velocity and pressure at the computational cell face (i ‡ 12): ~ i‡1=2 ˆ 1 (Ui ‡ Ui‡1 ) ± U 2 1 p~ i‡1=2 ˆ ( pi ‡ pi‡1 ) ± 2

1 2rCu rCp 2

i‡1=2

( pR ± pL )

i‡1=2

(UR ± UL )

2n Dn

(34)

p              u2 ‡ v2

(35)

Cu ˆ max c; Cp ˆ b

where the restriction on c ˆ (Cu ; Cp ) renders the physical viscosity e€ ective. Additional limitation on c in this way is well suited to viscous ¯ ows containing extensive low-speed portions such as thick boundary and shear layers. To recover particularly the desirable attributes pertaining to the cell-face interpolation method that minimize the continuity error and contribute a physically meaningful smooth pressure ® eld, relation (34) is reconstructed as

400

M. M. RAHMAN AND T. SIIKONEN

1 1 Ui‡1=2 ˆ (Ui ‡ Ui‡1 ) ± Cui‡1=2 (pR ± pL ) 2 2rCu i‡1=2 " # 4(|pi| ‡ |pi‡1 |) u Ci‡1=2 ˆ max 0; 1 ± r( |Ui | ‡ |Ui‡1 | ) 2 ‡E

(36)

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

and the corresponding reasonable expression for the cell face pressure pi‡1=2 becomes 1 rCp p pi‡1=2 ˆ (pi ‡ pi‡1 ) ± C i‡1=2 ( UR ± UL ) 2 2 i‡1=2 Á ! r( |Ui | ‡ |Ui‡1 | ) 2 Cpi‡1=2 ˆ max 0; 1 ± 4(|pi | ‡ |pi‡1 |) ‡ E

(37)

where the sign (~:) is dropped out for notational convenience. The monitor function C ˆ (Cu ; Cp ), dealing physically with the local pressure coe cient, provides insight into the manner by which the scheme is susceptible to controlling the degree of biasing the dissipation. The expressions for the remaining cell faces can be derived identically. 4.

BOUNDARY CONDITIONS

The boundary values are given in ghost cells such that the actual boundary conditions are obtained at the cell face when applying the central di€ erencing. Consequently, referring to a boundary cell face (i ± 12), the Dirichlet and Neumann conditions are formulated as W~ i± 1 ˆ 2W~ i± 1=2 ± W~ i ;

W~ i± 1 ˆ W~ i ±

Dn

qW~ qn

(38) i± 1=2

where (i ± 1) indicates the ghost cell nodal point and W~ i± 1=2 or (qW~ =qn) is the speci® ed boundary value. On the solid surface, the central expression of the viscous ¯ ux is replaced by a second-order one-sided formula. A second-order extrapolation from the computational domain is applied to evaluate the wall pressure. The boundary condition for the change in W~ is prosecuted by setting DW~ i± 1 ˆ 0 at the beginning of each iteration. 5.

SOLUTION ALGORITHM

The discretized equations are integrated in time using an improved ADI factorization. This method is based on the approximate factorization and the splitting of the Jacobians of the ¯ ux terms. The implicit stage is composed of a backward and a forward sweep in every coordinate direction. The sweeps are based on a ® rst-order upwind di€ erencing. The boundary conditions are treated explicitly. The implicit stage can be written after factorization as follows:

ARTIFICIAL COMPRESSIBILIT Y METHOD

Dti ± ‡ ± q S i‡1=2 A‡ i ± qi S i± 1=2 Ai 8i i Dti ± ‡ ± ~ ~ ¤¤¤ £ I‡ (q S j‡1=2 B‡ j ± qj S j± 1=2 Bj DW i ˆ DW i 8i j

401

I‡

(39)

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

where I is the identity matrix, qi;± j and q‡ i; j are ® rst-order backward and forward spatial di€ erence operators in the i and j directions, respectively, and A and B are the ¯ ux Jacobian matrices. The quantity DW~ ¤¤¤ ˆ (Dp; Du; Dv; Df) T represents the prii mitive explicit change in the dependent variables as de® ned by Eq. (24). In order to ensure stability of the viscous term, the Jacobian of the ¯ ux difference is reconstructed in the i direction as A§ ˆ L |l§ | ‡ yI L ± 1

(40)

where l§ are the diagonal matrices containing the positive and negative eigenvalues of the Jacobian pertaining to the convective term. The right and left eigenvector matrices L and L ± 1 , respectively, are given by Eq. (17). The factor y ensures the stability of the viscous term evaluated from yˆ

2G r Dn

(41)

where Dn is the cell height. In the i direction the tridiagonal equation set evolving from Eq. (39) is replaced by two bidiagonal sweeps and a matrix multiplication as ± ~ ¤¤ ~ ¤¤¤ (8i ‡ DtS i± 1=2 |A| i ) DW~ ¤¤ i ± DtS i‡1=2 Ai‡1 DW i‡1 ˆ 8i DW i ± (8i ‡ DtS i± 1=2 A‡ i ‡ DtS i‡1=2 Ai )

±1

8i DW~ ¤i ˆ DW~ ¤¤ i

(42)

~ ~¤ (8i ‡ DtS i‡1=2 |A| i ) DW~ i ± DtS i± 1=2 A‡ i± 1 DW i± 1 ˆ 8i DW i where Dt ˆ ZDti and Z ˆ 1:5, a factor governing the stability of implicitness. A similar procedure is attributed to the j direction. The basis herein for the splitting is the sequence of discrete transformation between primitive and characteristic representations. Accordingly, each of the equations (42) is solved in an e€ ectively explicit manner. For instance, the backward sweep can be inverted as ~ ¤¤¤ ‡ DtS i‡1=2 A ± DW~ ¤¤ DU¤¤ i ˆ 8i D W i i‡1 i‡1 ±1 DW~ ¤¤ L i ˆ T

L ± 1 T DU¤¤ i 8i ‡ DtS i± 1=2 (|l| ‡ yI)

(43)

starting with a boundary value for DW~ ¤¤ i‡1 ˆ 0 in the ghost cell. The interior nodal ± point value for Ai‡1 DW~ ¤¤ can be determined as i‡1 ± ±1 Ai‡1 DW~ ¤¤ LL i‡1 ˆ T

L ± 1 T DU¤¤ i 8i ‡ DtS i± 1=2 (|l| ‡ yI)

L ˆ max(0; ± l) ‡ yI

(44)

M. M. RAHMAN AND T. SIIKONEN

402

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

For completeness in simplicity and clarity of interpretation, L , T , (|l| ‡ yI), and L are evaluated at (i ± 12). The succession of matrix multiplication is executed as demonstrated. The other two steps can be performed analogously. The solution procedure consists of the following sequences: 1. Guess the velocity, pressure, and scalar ® elds. 2. Calculate the cell-face mass ¯ uxes using relations like Eq. (28). Alternatively, estimate the cell-face mass and pressure ¯ uxes from relations like Eqs. (36) and (37), respectively. 3. Evaluate primitive explicit changes using Eq. (24). 4. Solve the system of equations in accordance with the above-mentioned procedure. 5. Construct conservative correction-type vectors using Eq. (26). 6. Update the primitive variables from Eq. (27) to provide physical ® elds. 7. Repeat steps 2± 6 until convergence is achieved.

6.

TEST COM PUTATIONS

To demonstrate the e€ ectiveness and generality of the proposed formulation, a few applications to laminar ¯ ows are addressed. Tested ¯ ows consist of a viscositydriven ¯ ow in a square cavity, a buoyancy-drive n ¯ ow in a square cavity and a buoyancy-drive n ¯ ow in a half-concentric annulus. These cases are probably the most common computational experiments to validate the predictive performance of a solver. Since no analytically exact solutions to these problems are available, the results of the SIMPLE-based pressure correction method (PCM), details of which can be found in [13], are included for comparison. In the calculations that follow, CFL ˆ 10 and the compressibility parameter b ˆ 2. The state variables ( p; U) L and ( p; U) R associated with Eqs. (28), (36), and (37) are determined using the FUS scheme that induces con® dingly third-di€ erence dissipations. Within the proposed algorithmic framework, a comparative assessment of the momentum-based dissipation scheme (MDS) and the dual dissipation scheme (DDS) is made based on the convergence behaviors. The root mean square of the arti® cial density residual is used as a measure for the convergence, which can be de® ned as s                       PNP 2 |Dr | i iˆ1 Dr¤ ˆ µ Zt NP

Dr ˆ

Dp c2

(45)

where NP refers to the number of control volumes of the computational domain, and Zt is the user-de® ned tolerance limit.

7.

VISCOSITY-DRIVEN CAVITY FLOW

A schematic of the test problem con® guration is shown in [13] with the top wall moving to the right at a velocity Uref ˆ uw while the three sides are at rest. The

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

ARTIFICIAL COMPRESSIBILIT Y METHOD

403

problem can be converted to a dimensionless form using the characteristic length L , Uref , and rU2ref to nondimensionaliz e distances, velocities, and the pressure, respectively. The main dimensionless group evolving is the Reynolds number, Re ˆ Uref L =n. The boundary conditions are simple inasmuch as u ˆ v ˆ 0 on the bottom and two vertical walls, and u ˆ uw , v ˆ 0 on the top wall. A nonuniform 40 £ 40 grid, containing ® ner grid points near the walls than in the core to resolve the sharp velocity gradients, is employed for the computation. Plots of the horizontal velocity pro® le at the vertical centerline of the cavity are presented for Re ˆ 5; 000 in Figure 1, where U ˆ u=Uref and Y ˆ y=L . As is observed, the agreement is quite striking for all three formulations. Figure 2 displays the convergence histories of the arti® cial density residuals with the same parameters mentioned above, starting from the same initial conditions. The trends regarding the relative performances indicate that the DDS appears to have slight superiority over the MDS.

Figure 1. Horizontal velocity pro¢le at the vertical centerline for viscosity-drive n cavity £ow.

M. M. RAHMAN AND T. SIIKONEN

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

404

Figure 2. Convergence of arti¢cial density residuals for viscosity-drive n cavity £ow.

8.

BUOYANCY-DRIVEN CAVITY FLOW

The problem of natural convection is shown schematically in [13] with a characteristic dimension L , insulated at the top and bottom walls. The hot and cold vertical walls are both isothermal, having temperatures T h and T c , respectively. The overall system of equations in nondimensional form for the two-dimensional case is given in [11]. The problem is governed by two dimensionless parameters resulting from the normalization procedure: Gr ˆ gg(T h ± T c )L 3 =n2 and Pr ˆ n=a. The parameters g, a, and g are the acceleration due to gravity, thermal di€ usivity, and coe cient of thermal expansion, respectively. In the present computations, Pr is ® xed at unity. Comparisons are made by plotting the results in the form of V ˆ vL =n and y ˆ (T ± T p c )=(T     h ± T c ) versus X ˆ x=L . The reference velocity considered herein is Uref ˆ Gr : Computations involving a 40 £ 40 nonuniform grid for Gr ˆ 107 are conducted for the purpose of comparison. Shown in Figure 3 are the velocity and temperature pro® les at the horizontal midplane, which shows encouraging qualitative agreement. Figure 4 deals with the convergence histories with the same condition mentioned

405

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

ARTIFICIAL COMPRESSIBILIT Y METHOD

Figure 3. Velocity and temperature pro¢les on the horizontal midplane for buoyancy-drive n cavity £ow.

M. M. RAHMAN AND T. SIIKONEN

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

406

Figure 4. Convergence of arti¢cial density residuals for buoyancy-drive n cavity £ow.

above. Apparently, both the MDS and DDS correspond to the same approval about the arti® cial density convergence. 9.

BUOYANCY-DRIVEN FLOW IN AN ANNULUS

Attention is drawn to the half-concentric annulus shown schematically in [13]. The radii of the inner and outer surfaces are presented by Ri and Ro , respectively, with (Ro =Ri) ˆ 3. The origin of the Cartesian coordinate system is located at the center of the circles. The system has a characteristic length L set equal to (Ro ± Ri ) and is insulated at the top and bottom vertical walls. The inner curved surface is held at a temperature T h and the outer at T c , with T h > T c . The governing equations and boundary conditions are similar to those implemented for the buoyancy-drive n ¯ ow in a square cavity. The nonuniform grid employed herein is composed of 38 radial and 28 circumferential line seqments. Figure 5 exhibits the vertical velocity and temperature pro® les, respectively, at the horizontal midplane for Gr ˆ 106 , where X is measured

407

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

ARTIFICIAL COMPRESSIBILIT Y METHOD

Figure 5. Velocity and temperature pro¢les on the horizontal midplane for buoyancy-driven £ow in an annulus.

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

408

M. M. RAHMAN AND T. SIIKONEN

Figure 6. Convergence of arti¢cial density residuals for buoyancy-driven cavity £ow in an annulus.

exactly from the inner surface. Remarkably, for this particular grid re® nement, both the MDS and DDS solutions excel in producing impressive agreement with PCM data. Figure 6 shows the convergence histories of the arti® cial density residuals on the same grid with the same parameters. The plot is self-explanatory. 10.

CONCLUSIONS

A pressure-based method with an adherence to the arti® cial compressibility is developed for two-dimensional incompressible ¯ ows on a curvilinear collocated grid. The methodology uses a cell-centered ® nite-volume approximation with a factored pseudo-time integration approach. The arti® ce of compressibility persuades matrix/ density preconditionings, contributing assertively an implicit coupling among the primitive variables followed by an e€ ective linearization to the residual. The associated cell-face interpolation contrivance for ensuring the pressure± velocity coupling has strong motivation toward reducing the numerical dissipation being produced. In particular, the matrix inversion technique relating to the eigenstructure is quite simple and straightforward.

ARTIFICIAL COMPRESSIBILIT Y METHOD

409

To substantiate the validity of the proposed method, computations have been carried out for both forced- and natural-convectio n ¯ ows inside some selected geometries. Results demonstrate that the present algorithm compares favorably with the PCM. The convergence studies essentially dictate that the primitive formulation, utilizing either the momentum-base d dissipation scheme or the dual dissipation scheme at the cell face to prevent nonphysical oscillations, has competency in producing satisfactory stabilization for the iteration process. Obviously, the present procedure can readily be extended to the three-dimensional case.

Downloaded by [Aalto-yliopiston kirjasto] at 03:52 23 September 2014

REFERENCES 1. E. Turkel, Preconditioned Methods for Solving the Incompressible and Low Speed Compressible Equations, J. Comput. Phys., vol. 72, pp. 277± 298, 1987. 2. Y.-H. Choi and C. L. Merkle, The Application of Preconditioning to Viscous Flows, J. Comput. Phys., vol. 105, pp. 207± 223, 1993. 3. E. Turkel, R. Radespiel, and N. Kroll, Assessment of Preconditioning Methods for Multidimensional Aerodynamics, Comput. Fluids, vol. 26, pp. 613± 634, 1997. 4. J. M. Weiss and W. A. Smith, Preconditioning Applied to Variable and Constant Density Flows, AIAA J., vol. 33, pp. 2050± 2057, 1995. 5. M. M. Rahman, P. Rautaheimo, and T. Siikonen, Numerical Study of Turbulent Heat Transfer from a Con® ned Impinging Jet Using a Pseudo-Compressibility Method, 2nd Int. Sym. on T urbulence, Heat and Mass Transfer, Delft, The Netherlands, pp. 511± 520, 1997. 6. A. J. Chorin, A Numerical Method for Solving Incompressible Viscous Flow Problems, J. Comput. Phys., vol. 2, pp. 12± 26, 1967. 7. W. K. Anderson, R. D. Rausch, and D. L. Bonhaus, Implicit=Multigrid Algorithms for Incompressible Turbulent Flows on Unstructured Grids, J. Comput. Phys., vol. 128, pp. 391± 408, 1996. 8. P. L. Roe, Approximate Riemann Solvers, Parameter Vectors, and Di€ erence Schemes, J. Comput. Phys., vol. 43, pp. 357± 372, 1981. 9. C. M. Rhie and W. L. Chow, Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation, AIAA J., vol. 21, pp. 1525± 1532, 1983. 10. S. W. Arm® eld, Finite Di€ erence Solutions of the Navier± Stokes Equations on Staggered and Nonstaggered Grids, Comput. Fluids, vol. 20, pp. 1± 17, 1991. 11. M. M. Rahman, A. Miettinen, and T. Siikonen, Modi® ed SIMPLE Formulation on a Collocated Grid with an Assessment of the Simpli® ed Quick Scheme, Numer. Heat T ransfer B, vol. 30, pp. 291± 314, 1996. 12. M. M. Rahman, T. Siikonen, and A. Miettinen, A Pressure Correction Method for Solving Fluid Flow Problems on a Collocated Grid, Numer. Heat T ransfer B, vol. 32, pp. 63± 84, 1997. 13. M. M. Rahman and T. Siikonen, An Improved SIMPLE Method on a Collocated Grid, Numer. Heat T ransfer B, vol. 38, pp. 177± 201, 2000. 14. C. Lombard, J. Bardina, E. Venkatapathy, and J. Oliger, Multi-Dimensional Formulation of CSCMÐ An Upwind Flux Di€ erence Eigenvector Split Method for the Compressible Navier± Stokes Equations, 6th AIAA Computational Fluid Dynamics Conference, Danvers, MA, AIAA Paper 83-1895-CP, pp. 1453± 1460, 1986. 15. W. K. Anderson, J. L. Thomas, and B. Van Leer, Comparison of Finite Volume Flux Vector Splittings for the Euler Equations, AIAA J., vol. 24, pp. 649± 664, 1983.