Amano R.S., Sunden B. (Eds.) Computational Fluid Dynamics and Heat Transfer: Emerging Topics

Подождите немного. Документ загружается.

Sunden CH011.tex 10/9/2010 15: 22 Page 456

456 Computational Fluid Dynamics and Heat Transfer

[51] Oldham, M. J., Mannix,R.C., and Phalen, R.F. Deposition ofmonodisperseparticles

in hollow models representing adult and child-size tracheobroncial airways, Health

Physics, 72(6), pp. 827–833, 1997.

[52] Broday,D.M. Deposition of ultraﬁneparticles at carinal ridges oftheupperbronchial

airways,Aerosol Science andTechnology, 38, pp. 991–1000, 2004.

[53] Zhang, Z., and Kleinstreuer, C. Airﬂow structures and nano-particle deposition in a

human upper airway model, Journal of Computational Physics, 198(1), pp. 178–210,

2004.

[54] Chen,F.,andLai,A.C.K.AnEulerianmodelforparticledepositionunderelectrostatic

and turbulent conditions, Journal of Aerosol Science, 35, pp. 47–62, 2004.

[55] Cohen, B. S., Sussman, R. G., and Lippmann, M. Ultraﬁne particle deposition in a

human tracheobronchial cast, Aerosol Science and Technology, 12, pp. 1082–1093,

1990.

[56] van Ertbruggen, C., Hirsch, C., and Paiva, M. Anatomically based three-dimensional

model of airways to simulate ﬂow and particle transport using computational ﬂuid

dynamics, Journal of Applied Physiology, 98(3), pp. 970–980, 2005.

[57] Pedley, T. J. Pulmonary ﬂuid dynamics, Annual Review of Fluid Mechanics,9,

pp. 229–274, 1977.

[58] Chan,T. L., Schreck, R. M., and Lippmann, M. Effect of the laryngeal jet on particle

deposition in the human trachea and upper bronchial airways, Journal of Aerosol

Science, 11, pp. 447–459, 1980.

[59] Gurman, J. L., Lippmann, M., and Schlesinger, R. B. Particle deposition in replicate

casts of the humanuppertrancheobronchialtree under constant andcyclicinspiratory

ﬂow. I. Experimental,Aerosol Science and Technology, 3, pp. 245–252, 1984.

[60] Martonen,T.B.,Zhang,Z.,andLessmann,R.Fluiddynamicsofthehumanlarynxand

upper tracheobronchial airways, Aerosol Science and Technology, 19, pp. 133–144,

1993.

[61] Corcoran,T.E.,andChigier,N.Inertialdepositioneffects:Astudyofaerosolmechan-

ics in the trachea using laser Doppler velocimetry and ﬂuorescent dye, Journal of

Biomechanical Engineering, 124, pp. 629–637, 2002.

[62] Kabilan, S., Lin, C. L., and Hoffman, E. A. Characteristics of airﬂow in a CT-based

ovine lung: A numerical study, Journal of Applied Physiology, 102, pp. 1469–1482,

2007.

[63] Li, Z., Kleinstreuer, C., and Zhang, Z. Particle deposition in the human tracheo-

bronchialairwaysduetotransientinspiratoryﬂowpatterns,JournalofAerosolScience,

38, pp. 625–644, 2007.

[64] Lin, C. L., Tawhai, M. H., McLennan, G., and Hoffman, E.A. Characteristics of the

turbulent laryngeal jet and its effect on airﬂow in the human intra-thoracic airways,

Respiratory Physiology and Neurobiology, 157, pp. 295–309, 2007.

[65] Nowak, N., Kakade, P. P., and Annapragada, A. V. Computational ﬂuid dynamics

simulation of airﬂow and aerosol deposition in human lungs, Annals of Biomedical

Engineering, 31, pp. 374–390, 2003.

[66] Zhang, Z., Kleinstreuer, C., and Kim, C. S. Cyclic micron-size particle inhalation

and deposition in a triple bifurcation lung airway model, Journal of Aerosol Science,

33(2), pp. 257–281, 2002.

[67] Cheng, K. H., Cheng, Y. S., Yeh, H. C., and Swift, D. L. Deposition of ultraﬁne

aerosols in the head airways during natural breathing and during simulated breath-

holding using replicate human upper airway casts, Aerosol Science and Technology,

23(3), pp. 465–474, 1995.

Sunden CH011.tex 10/9/2010 15: 22 Page 457

Evaluation of continuous and discrete phase models 457

[68] Cheng,Y. S.,Yeh, H. C., Guilmette, R.A., Simpson, S. Q., Cheng, K. H., and Swift,

D. L. Nasal deposition of ultraﬁne particles in human volunteers and its relationship

to airway geometry,Aerosol Science and Technology, 25(3), pp. 274–291, 1996.

[69] Shi, H., Kleinstreuer, C., and Zhang, Z. Laminar airﬂow and nanoparticle or vapor

deposition in a human nasal cavity model, Journal of Biomechanical Engineering,

128, pp. 697–706, 2006.

This page intentionally left blank

Sunden CH012.tex 17/8/2010 20: 31 Page 459

12 Algorithm stabilization and acceleration in

computational ﬂuid dynamics: exploiting

recursive properties of ﬁxed point

algorithms

Aleksandar Jemcov and Joseph P. Maruszewski

ANSYS, Inc., Canonsburg, PA, USA

Abstract

Convergence acceleration of nonlinear ﬂow solvers through use of techniques that

exploit recursive properties ﬁxed-point methods of computational ﬂuid dynam-

ics (CFD) algorithms is presented. It is shown that nonlinear iterative algorithms

create Krylov subspace that can be considered linear within the second order

of accuracy in the vicinity of the local solutions. Moreover, repeated applica-

tion of the ﬁxed-point iterations give rise to various vector extrapolation methods

that are used to accelerate underlying iterative method. In particular, suitability

of the reduced rank extrapolation (RRE) algorithm for the use of convergence

acceleration of nonlinear ﬂow solvers is examined. In the RRE algorithm, the

solution is obtained through a linear combination of Krylov vectors with weight-

ing coefﬁcients obtained by minimizing L

norm of error in this space with

properly chosen constraint conditions. This process effectively deﬁnes vector

sequenceextrapolationprocessinKrylovsubspacethatcorrespondstotheGMRES

method applied to nonlinear problems. Moreover, when the RRE algorithm is

used to solve nonlinear problems, the ﬂow solver plays the role of the precon-

ditioner for the nonlinear GMRES method. Beneﬁts of the application of the

RRE algorithm include better convergence rates, removal of residual stalling, and

improved coupling between equations in numerical models. Proposed algorithm

is independent of the type of ﬂow solver and it is equally applicable to explicit-,

implicit-, pressure-, and density-, based algorithms. RRE algorithm can also be

considered as an generalization of matrix-free algorithms that are used extensively

in CFD.

Keywords: Fixed-point methods, Krylov subspace, Recursive properties, Reduced

rank extrapolation,Vector extrapolation

Sunden CH012.tex 17/8/2010 20: 31 Page 460

460 Computational Fluid Dynamics and Heat Transfer

12.1 Introduction

Computational ﬂuid dynamics (CFD) in modern engineering applications is char-

acterized by complex physical models, complex shapes of computational domain,

and large mesh sizes. Governing equations of ﬂuid ﬂow and associated equations

for the transport of conserved variables are seemingly simple and yet they contain

very prominentnonlinearitiesof whichthe most famousoneis manifested through

the phenomena of turbulence. The need to model turbulent ﬂows coupled with

radiation, heatandmass transfer, andmultiphaseinteractions leadto very complex

physical models that contain many nonlinear terms and wide range of spatial and

temporal scales. This complexity directly translates into difﬁculties in numerical

applications and they are loosely referred to here as computational stiffness (for

the current purposes of the qualitative discussion, we will use the term stiffness

to signify computational difﬁculties). Another source of stiffness in CFD is asso-

ciated with the use of large mesh sizes in order to capture relevant spatial scales

associated with physical models. It is well-known fact that the condition number

of linear system increases with the increase of the mesh size as it can be veriﬁed

by the uniform reﬁnement of the computational mesh. Furthermore, topological

complexities of the computational domains lead CFD practitioners to the use of

unstructured meshes with nearly arbitrary shapes of control volumes. The main

issues of unstructured meshes such as presence of highly skewed volumes, high

aspect ratios close to the boundaries, loss of orthogonality, and lack of convexity

directly contribute to the stiffness of the computational problem. This is usually

manifested in the lossof the diagonal dominance of the matrix inthe linear system

of discrete equations. Iterative solvers for linear systems of equations are not par-

ticularly effective when the matrix of the linear system departs strongly from the

diagonal dominance expressed through M-matrix condition.

Complexity is not just restricted to physical models and related issues as

described above; ﬂow simulation is also characterized by a great number of dif-

fering algorithms ranging from explicit to fractional step to implicit methods [1].

Furthermore, dependingon thechoice ofthesolution variablesand chosendiscrete

formulation, pressure [2–4] and density-based solvers [5,6] can be distinguished.

Pressure-anddensity-based solvers canbe formulated as implicit,explicit, orfrac-

tionalstep methods. Complexity ofthe landscapeis further increasedbythe choice

of linearsolver usedin implicit algorithms.Additionally, variousﬂux formulations

are in the use today including Riemann solver, ﬂux vector splitting, central dif-

ference with dissipation ﬂux formulation [6], to name just a few. Each of these

algorithms have control parameters that have optimal values for a given applica-

tion area thus producing an optimal convergence rate during the iterative process.

Examples of control parameters are CFL numbers in density-based and relaxation

factors in pressure-based algorithms.

Optimal range of control parameters are deﬁned by the stability constraints

of the numerical algorithm, underlying physics, computational mesh, geometrical

complexityofthecomputationaldomain,andboundaryconditions. Duetothenon-

linear nature of equations describing the ﬂow phenomena, selection of the values

Sunden CH012.tex 17/8/2010 20: 31 Page 461

exploiting recursive properties of ﬁxed point algorithms 461

of control parameters that are outside of the optimal range results in either slow

convergence rates or divergence of the residual norm. Linear stability analysis [5]

provides some guidelines for the selection of the optimal values of control param-

eters but in the industrial applications of ﬂow solvers, problem complexity due

to physics, domain geometry, boundary conditions, etc. dramatically changes the

optimal range of control parameters. Furthermore, the nonlinear nature of govern-

ing equations changesthe optimal valuesof control parameters during theiterative

process. Duringinitial iterations, dependingon theinitial guess,optimal values are

very different from the optimalvalues ofthecontrol parameters during laterstages

in iterative process. Usually, initially conservative values of control parameters are

used to maintain the stability of the iterative process and as iterations approach the

solution, more aggressive values can be used.

The need to dynamically change the values of control parameters during itera-

tions is addressed in algorithms that compute optimal or nearly optimal values by

utilizing the knowledge of the state of computed variables and some criterion that

measures stability of the computed state. Continuation methods [7] provide a pos-

sible solution to this problem.Another approach to the problem of optimal control

parametersisintheuseofbacktrackingorlinesearchalgorithminconjunctionwith

Newtonmethods [8]. Both continuation and Newton with backtracking algorithms

are effective tools for the solution of the problem of dynamically changing control

parametersbut they requirealgorithmicchanges totheﬂowsolver. Manytimesitis

not feasible to implement new algorithms in large codes, and even though the new

algorithms provide added stability and increase convergence rates, they cannot be

implemented in timely manner. Less intrusive ways of improving the convergence

rates and overall stability of the ﬂow solver are required.

Anothersourceofinstabilityandlowconvergenceratesinﬂowsolversisrelated

to artiﬁcial segregation of discrete equations in the numerical model. Example of

the algorithms that use segregated approach are the SIMPLE family of algorithms

[2,4]. Another example of segregated algorithms is found in the treatment of tur-

bulence models [9] in density-based algorithms where they are lagged behindﬂow

quantities. In multi-equation turbulence models, each turbulent transport equation

may be solved separately, thus artiﬁcially decoupling equations in the numerical

model.Decouplingofequationsinsegregatedﬂowsolvershasadetrimentalimpact

on convergence rates. The solution to this problem is usually found in creation of

coupled solvers at the expense of increased memory and code complexity.

Vectorsequenceextrapolationmethods[10] offeranalternativesolution tocon-

vergenceaccelerationandstabilizationof ﬂowsolvers.The ideaofvectorsequence

extrapolation is based on the generalization of scalar sequence convergence accel-

erationmethods,suchasRichardsonextrapolation[10],tovectorsequences.Vector

sequences are produced byiterative procedures of solver algorithms while seeking

the solution of the given ﬂow problem. In addition to accelerating properties of

vectorsequence extrapolation methods, theyalso exhibit couplingproperties since

only one set of extrapolation coefﬁcients is calculated for all ﬁelds. This has a

positive effect on convergence rates, very similar to the action of Krylov subspace

solvers acting on systems with block coefﬁcients. Moreover, formal equivalence

Sunden CH012.tex 17/8/2010 20: 31 Page 462

462 Computational Fluid Dynamics and Heat Transfer

between vector sequence extrapolation methods to Krylov subspace solvers such

as GMRES and FOM [11] exists.

Early examples of the application [13] of vector sequence methods in CFD

was the application of reduced rank extrapolation (RRE), minimal polynomial

extrapolation (MPE), and modiﬁed minimal polynomial extrapolation (MMPE)

toexplicitﬂowsolversbasedonJSTscheme [12]andimplicitﬂowsolverwithﬂux

vector splitting formulation. Acceleration of convergence rates of inviscid ﬂow

solverswasexamined[13] andmarginalimprovements wereobtained withexplicit

solvers, whereas implicit solvers showed better results. Another application [14]

of an extrapolation-related technique to several types of ﬂow demonstrated the

suitability ofextrapolation toaccelerateinviscid compressibleand viscous laminar

incompressible ﬂows.The method [14] used in the studies relied on residual mini-

mizationandrequiredmanyevaluationsofresidualsmakingthismethodexpensive.

Other applications of vector sequence extrapolation to nonlinear problems include

Chandrasekhar’s H-equation [15]. Vector sequence extrapolation is also used for

accelerationofstationaryiterativemethods[16]andthealgebraicmultigrid(AMG)

method [17].

Traditional view of extrapolation methods applied to vector sequences is that

they perform the transformation of one sequence into another sequence result-

ing with faster convergence rates. Since a ﬂow solver will be generating vector

sequences that will be used to compute extrapolation coefﬁcients, an alternative

viewofthisprocessisthattheﬂowsolverplaystheroleofanonlinearpreconditioner

for an extrapolation solver. Due to equivalence of RRE and GMRES algorithms,

the proposed acceleration technique corresponds to a nonlinearly preconditioned

GMRES solver.

The outlineof the paper isas follows. CFD algorithms were examined from the

ﬁxed-pointalgorithmpointofviewandtherecursivepropertyofthesealgorithmsis

shown. Recursive property ofthe ﬁxed-point algorithmsrepresent, to second-order

approximation, a power iteration with the matrix of Jacobian of the ﬁxed-point

function.Thisproperty playsan important rolein thederivationof vector sequence

extrapolation algorithms. This is followed by a description of the vector sequence

extrapolation method where the recursive property of the ﬁxed-point function was

usedtoobtaintheextrapolationcoefﬁcients.Thisisdonethroughtheformulationof

an optimization problem in the Krylov subspace and a proper choice ofconstraints

leadtotheRREalgorithm. Inthesectiononnumericalexperiments,theapplication

ofnonlinearaccelerationisdemonstratedonvariousCFDproblemsutilizingseveral

ﬂow solvers. Finally, a summary and conclusion is presented, following results of

numerical experiments.

12.2 Iterative Methods for Flow Equations

In this section we will consider the discrete form of the governing equations and

iterative methods of ﬁnding the solution for nonlinear problems in CFD. Iterative

methods are considered from the point of view of the ﬁxed-point algorithms that

Sunden CH012.tex 17/8/2010 20: 31 Page 463

exploiting recursive properties of ﬁxed point algorithms 463

have a recursive property related to powers of the Jacobian of the ﬁxed-point func-

tion.RecursivepropertyplaysimportantroleinestablishingtheexistenceofKrylov

subspace in which the solution ofthe problem isapproximated. In fact, we demon-

strate how the recursive property generates Krylov subspace for any ﬁxed-point

functionthusenablingacceleration andstabilization ofarbitraryiterativeproblems

under the certain conditions.

12.2.1 Discrete form of governing equations

Consider Navier–Stokes system of equations in integral form

∂

∂t



Qd +

∂

− F

)dA −



Sd = 0 (1)

Here Q represents vector of conserved quantities, F

and F

are convective and

viscous ﬂux vectors, and S is the source term vector

Q =

⎛

⎜

⎝

ρu

ρv

ρw

ρE

⎞

⎟

⎠

, F

⎛

⎜

⎝

ρV

ρuV + n

ρvV + n

ρwV + n

ρHV

⎞

⎟

⎠

, F

⎛

⎜

⎝

+ n

⎞

⎟

⎠

S =

⎛

⎜

⎝

ρf

e,x

ρf

e,y

ρf

e,z

ρ(f

e,x

u +f

e,y

v +f

e,z

w) +˙q

⎞

⎟

⎠

and V and

, and

are deﬁned as:

V = n

u +n

v +n

= uτ

+ vτ

+ wτ

+ k

∂T

∂x

= uτ

+ vτ

+ wτ

+ k

∂T

∂y

= uτ

+ vτ

+ wτ

+ k

∂T

∂z

The system of equations in equation (1) may include additional scalar or vector

transport equations

∂

∂t



d +

∂



,c

− F

,v



dA −





d = 0 (2)

Sunden CH012.tex 17/8/2010 20: 31 Page 464

464 Computational Fluid Dynamics and Heat Transfer

Depending on the choice of the transported quantity, equation (2) can represent

species, passive scalar, turbulent transport, etc. Appropriate boundary conditions

[1] must be also supplied to deﬁne a well-posed problem that can be solved

numerically.

The ﬁrst step in obtaining the discretized equations is to deﬁne a residual R(Q)

R(Q) =

∂

(

− F

)

dA −



Sd (3)

Equation(3)deﬁnesacontinuousresidualthatcanbediscretizedonacomputational

meshusingacell-centeredﬁnitevolumemethod. Forithﬁnitevolumecell,discrete

residual is given by:

(Q) =



− F

)

− S

(4)

Here symbol



( · ) represents summation over all faces of a given ﬁnite volume

cell i. If the discrete value of the conserved variables vector Q

is deﬁned as:



Qd



(5)

The discrete representation of the Navier–Stokes system in equation (1) becomes

(

Q

)

=−R

(Q) (6)

The discrete system of equations given by equation (6) represents the nonlinear

discrete system due to nonlinear nature of the residual R

. An explicit algorithm

canusethesysteminequation(6)directlybyapplyingthetimemarchingprocedure

such as the multistage Runge–Kutta [1] procedure and full approximation scheme

(FAS) [18] to obtainthe solution of theproblem.The explicitalgorithm is given by

the following expression:

Q

(ν+1)

=−

t

(



)

(ν)

) (7)

(ν+1)

= Q

(ν)

+Q

(ν)

In the case of an implicit algorithm, a Newton or Picard linearization of the resid-

ual is applied to form an iterative scheme. If Newton linearization is used, the

discretized system of equations is given by the following expression:



()

t

+ (∂



Q

(ν)

=−R

(ν)

) (8)

(ν+1)

= Q

(ν)

+Q

(ν)

Sunden CH012.tex 17/8/2010 20: 31 Page 465

exploiting recursive properties of ﬁxed point algorithms 465

Since the solution of the system in equation (8) requires inversion of matrix on

the left-hand side, iterative methods for solving large sparse systems of linear

equationsare used.Popularchoicesfor linearsolversareKrylov subspacemethods

[11] and AMG method [18]. There are also number of choices for the convective

and viscous ﬂux evaluations at cell faces, popular choices in density-based solvers

are ﬂux difference splitting and ﬂux vector splitting [1], whereas in the case of

pressure-based solvers, collocated schemes with Rhie–Chow dissipation [4] are

often used.

12.2.2 Recursive property of iterative methods

Both explicit and implicit algorithms are particular instances of the ﬁxed-point

method that can be generally expressed as:

(ν+1)

= Q

(ν)

− M

(ν)

(R(Q

(ν)

)) (9)

The general form of any ﬁxed-point method is

(ν+1)

= F(Q

(ν)

) (10)

with ﬁxed-point function F(Q

(ν)

)

F(Q

(ν)

) = Q

(ν)

− M

(ν)

(R(Q

(ν)

)) (11)

Operator M

(ν)

(R(Q

(ν)

)) is, in principle, a nonlinear operator that may not even

be known explicitly. In the case of stationary linear iteration schemes such as

Gauss–Seidel, Jacobi, and ILU(n), this operator may be constructed explicitly.

However, in the case of theAMG method with multiple coarse levels, it is usually

knownimplicitlythroughimplementationintheﬂowsolver.Inthecaseofnonlinear

ﬂow solvers, this operator is almost always deﬁned algorithmically through code

implementation except in trivial cases. However, as it will be shown later, this

knowledge is not required for the vector sequence extrapolation method to work.

Animportant propertyof any ﬁxed-pointalgorithmisthe recursiverelation that

canbeestablishedbetweentheinitialsolutioncorrectionQ

(0)

andlatercorrections

Q

(ν)

through powers of the ﬁxed-point function Jacobian ∂

F. This relation can

be established by forming the following difference:

Q

(ν)

= Q

(ν+1)

− Q

(ν)

(12)

From the deﬁnition of the ﬁxed-point algorithm, equation (10), equation (12)

becomes

Q

(ν)

= F(Q

(ν)

) −F(Q

(ν−1)

) (13)

By the deﬁnition, Q

(ν)

is given by

(ν)

= Q

(ν−1)

+Q

(ν−1)

(14)