NASA/TM-2005-213764 

Factorization of the Compressible NavierStokes Equations 

Thomas W. Roberts 

Langley Research Center, Hampton, Virginia 

 

December 2005 
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role. 

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for NASA’s 

scientific and technical information. The NASA STI 

Program Office provides access to the NASA STI 

Database, the largest collection of aeronautical and 

space science STI in the world. The Program Office is 

also NASA’s institutional mechanism for 

disseminating the results of its research and 

development activities. These results are published by 

NASA in the NASA STI Report Series, which 

includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results ... even 

providing videos. 

For more information about the NASA STI Program 

Office, see the following: 

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 

at (301) 621-0134 

• Phone the NASA STI Help Desk at 

(301) 621-0390 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7121 Standard Drive 

 Hanover, MD 21076-1320
NASA/TM-2005-213764 

Factorization of the Compressible NavierStokes Equations 

Thomas W. Roberts 

Langley Research Center, Hampton, Virginia 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

December 2005 
Available from: 

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS) 

7121 Standard Drive 5285 Port Royal Road 

Hanover, MD 21076-1320 Springfield, VA 22161-2171 

(301) 621-0390 (703) 605-6000 
Contents

Abstract iv

1 Introduction 1

2 The Navier-Stokes equations 2

3 Relaxation and principal linearization 2

3.1 Principal linearization of inviscid terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3.2 Principal linearization of the viscous terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.3 The complete principal linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.4 Interpretation of the principal linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4 Factors of the Navier-Stokes equations 6

4.1 Inviscid equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

4.2 Low speed flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4.3 A special case of the Prandtl number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

5 Modes of the principal linearization 9

6 Conclusions 11

iii
Abstract

The Navier-Stokes equations for a Newtonian ideal gas are examined to determine the factorizable form of

the equations relevant to the construction of a factorizable relaxation scheme. The principal linearization

of the equations is found by examining the relative magnitude of the terms for short-wavelength errors.

The principal part of the operator is then found. Comparison of the factors of the Navier-Stokes and Euler

equations differ qualitatively because of the coupling of entropy and pressure through thermal diffusion.

Special cases of the factorization are considered.

iv
1 Introduction

Considerable progress has been made in recent years in the development of fast multigrid methods for the

Euler and Navier-Stokes equations. In particular, the development of factorizable discretizations shows the

promise of attaining ideal, or textbook, multigrid efficiency. A factorizable discretization is one in which the

governing equations are partitioned into component elliptic, hyperbolic and parabolic subsystems. Each of

these subsystems is then solved or relaxed by an appropriate numerical algorithm. By treating each partition

optimally, the full system can be solved with the efficiency of the slowest converging subsystem. These ideas

have been propounded by Achi Brandt for some time (Brandt, 1985) and are the subject of a recent review

by Thomas, Diskin, and Brandt (2003).

In the multigrid literature, factorizability is described in terms of the properties of the determinant of

the system of equations (Sidilkover, 1999; Ta’asan, 1993; Thomas et al., 2003). To relax the system of

equations, the increment in the dependent variables is written in terms of a set of auxiliary variables, or

“ghost” variables in Brandt’s terminology. The transformation of variables is designed such that the system

of equations in terms of the auxiliary variables is in triangular form, with the diagonal elements being the

factors of the determinant. Increments in the auxiliary variables are used to update the original variables

though the transformation. The triangularization of the system by means of the ghost variables is nothing

more than decomposing the system into a set of nearly linearly independent modes. The factors of the

determinant are the operators which act on each mode, and the auxiliary variables are the mode amplitudes.

The key point is that the modes represent the short-wavelength components of the solution error. Only the

terms in the equations that contribute strongly at the short wavelengths need to be considered in factoring

the equations. These terms are called the principal part of the operator.

The factorizable property may be examined in a variety of ways. Thomas et al. (2003) view the construction of a factorizable discrete operator by treating the relaxation process as being an approximation to

a Newton iteration. Starting with a linearization of the full system, they eliminate the subprincipal terms

of the Jacobian and express the solution updates in terms of ghost variables. Ta’asan (1993) works directly

with the full nonlinear Euler equations to derive a set of “canonical variables” corresponding to the different

partitions of the equations. Using these variables he constructs a multigrid solver based on direct relaxation

of the nonlinear system (Ta’asan, 1994). Sidilkover (1999) devises a factorizable scheme for the Euler equations by designing the difference operators such that the discrete potential, entropy and vorticity modes are

preserved exactly in the constant coefficient case. This is generalized to a conservative discretization of the

full nonlinear system in a straightforward way (Roberts, Sidilkover, and Thomas, 2000). As with Ta’asan,

Sidilkover prefers to directly relax the full nonlinear system.

Although factorizable schemes have been constructed and demonstrated textbook multigrid efficiency for

the incompressible Navier-Stokes equations (Brandt & Yavneh, 1992; Thomas, Diskin, and Brandt, 2001;

Swanson, 2001) and for the incompressible and compressible Euler equations (Roberts et al., 2000; Roberts,

Swanson, and Sidilkover, 1999; Ta’asan, 1994) in two dimensions, the development of such schemes for the

compressible Navier-Stokes equations is still in its infancy. An initial attempt at such a scheme has been

reported by Thomas, Diskin, and Brandt (1999). A difficulty arises because the nature of the factors in

the viscous compressible case is qualitatively different than either the viscous incompressible or the inviscid

compressible equations. Thermal conductivity has the effect of coupling the energy and continuity equations

in compressible viscous flow. For certain values of the Prandtl number, described below, the two factors

can be decoupled, and in general the coupling is weak. Thomas et al. (2003) present one possible approach

to factorizing the compressible Navier-Stokes equations. Relaxation of their system will require a block

relaxation of the two thermodynamic variables, which are pressure and internal energy in their case. A

possible disadvantage of their formulation is that the energy and continuity modes do not decouple even for

the appropriate value of the Prandtl number.

In the remainder of this paper, the principal linearization of the compressible Navier-Stokes is developed.

The full Navier-Stokes equations are presented in Section 2. Following Thomas et al. (2003), in Section 3

the motivation is given by viewing relaxation methods as approximate Newton solvers. The principal linearization is derived from a perturbation expansion of the full equations. A discussion of various limits and

1
special cases of the factorization is presented in Section 4. The inviscid and incompressible limits are shown

to lead to the appropriate factorization of the Euler and incompressible Navier-Stokes equation. With this

factorization in mind, a decomposition into vortical, entropy and acoustic modes is developed, which may

form the basis of a fast multigrid solver, in Section 5. These modes transform smoothly to the inviscid case

as well. In addition, this decomposition has a counterpart in the study of compressible turbulent flow.

2 The Navier-Stokes equations

The analysis which follows is for a Newtonian fluid obeying an ideal gas law, with constant specific heats.

Only steady flows are considered. For the purposes of examining the factors of the Navier-Stokes equations,

the primitive variables are the most convenient. Letting (ui, p, s) be the state vector of velocity, pressure,

and entropy, the Navier-Stokes equations may be written

ujui,j +

1

ρ

p,i =

1

ρ

τij,j , (1a)

ρc2uj,j + ujp,j = (γ − 1)(Φ − qj,j ), (1b)

ujs,j =

1

ρT

(Φ − qj,j ), (1c)

where ρ is the density, T is the static temperature and γ is the ratio of the specific heats. The Einstein

summation convention is used, and partial differentiation with respect to the i-th coordinate is indicated by

the comma subscript notation. The viscous stress tensor is defined as τij = μ

0dkkδij +2μdij , where δij is the

Kronecker delta and dij is the deviatoric strain tensor dij = (ui,j + uj,i) /2. The first and second coefficients

of viscosity μ and μ

0 are related through the definition of the bulk viscosity, 3κ = 3μ0 + 2μ ≥ 0. By Stokes’

hypothesis κ = 0, which gives μ

0 = −2μ/3. However, it is not necessary to use this relation to determine the

principal linearization of the equations.

The two terms on the right-hand side of the entropy equation contain the effects of viscous heating and

thermal diffusion. The dissipation function Φ = dij τij ≥ 0 is the viscous work done on the fluid. The heat

flux vector is given by qi = −λT,i = − (μcp/Pr) T,i, where λ is the thermal conductivity, cp is the specific

heat at constant pressure, and Pr is the Prandtl number. The temperature and density may be expressed

in terms of the primary thermodynamic variables s and p:

T

T0

=



p

p0

(γ−1)/γ

e

(s−s0)/cp

,

ρ

ρ0

=



p

p0

1/γ

e

−(s−s0)/cp

where the subscript 0 denotes a reference value. The equation of state can also be written p = ((γ − 1) /γ) ρcpT,

and the speed of sound is given by c

2 = γp/ρ = (γ − 1) cpT.

3 Relaxation and principal linearization

Relaxation methods are a class of iterative methods for the solution of a system of algebraic equations. An

iterative method for solving a system of equations Aφ = b consists of generating a sequence of approximations

{φ0 φ1 φ2 . . .} such that limn→∞ φn = φ. The system is linearized about the n-th iterate to obtain an

equation for the error, Aφ(φ − φn) ≈ b − Aφn where Aφ ≡ ∂(Aφ)/∂φ is the Jacobian of A and b − Aφn is

the residual at the n-th iteration. When the system of equations is linear, then Aφ = A, and the equation

for the error is exact.

A Newton solver is obtained when all the terms of the linearization are kept and the resulting system

for the error is solved exactly. For discretizations of partial differential equations, such as the Navier-Stokes

2
equations, this requires the inversion of a sparse, banded matrix. Newton methods have the outstanding

property of quadratic convergence once the approximate solution is near the exact solution. Unfortunately,

direct inversion is expensive and often impractical. Furthermore, the initial iterates of Newton’s method can

be far from the exact solution—for nonlinear problems, many iterations may be required before the region

of quadratic convergence is reached.

Relaxation methods avoid direct inversion by replacing the Jacobian Aφ with an approximation that

is inexpensive to invert. Generally, this means throwing away all or some of the off-diagonal terms of the

Jacobian. For example, keeping just the diagonal terms leads to Jacobi iteration; the Gauss-Seidel method

results from throwing away the terms above the diagonal. Inverting the matrix involves only the cost of a

matrix multiplication, which because of the sparsity of the matrix is O(n).

Because it involves only local operations, relaxation is effective in eliminating the short-wavelength components of the error, where the short wavelengths are those on the order of h, the grid spacing. For this reason

only those terms of the linearization that contribute strongly to the short-wavelength error components need

be kept. These terms are referred to as the “principal part” of the linearization. Efficient multigrid methods

are based on a relaxation method, or smoother, for reducing the short-wavelength component of the error,

and a coarse-grid operator to reduce the long-wavelength error.

In Thomas et al. (2003) the authors describe the process of determining the principal part by reference

to the operator alone. Taking the determinant, only those terms that contribute to the leading order are

retained. In the following treatment, the principal part is found by expanding the equations in a small

perturbation about a smooth state and keeping the dominant terms. This leads to the same result as

Thomas, et al., but has the advantage of showing directly why certain terms are subprincipal. Furthermore,

the procedure followed here is more familiar to fluid dynamicists and is very similar to a perturbation theory

used in the study of compressible turbulence (Chu & Kov ́asznay, 1958; Kov ́asznay, 1953).

3.1 Principal linearization of inviscid terms

The principal linearization of the Navier-Stokes equations starts by considering a full linearization and

keeping only those terms that are dominant on the scale of the grid spacing. The linearization is found by

replacing the variables (s, ui, p) with

ui ← ui + uˆi, (2a)

p ← p + p, ˆ (2b)

s ← s + sˆ (2c)

where   1. Equations (1) are expanded in  and the linear terms are kept.

The inviscid terms of the Navier-Stokes equations (1) are those terms on the left-hand side. Let U be the

velocity scale and L be the physical length scale of the flow. These are local scales, and can be different in

different regions of the flow. For example, L can be a chord length, leading-edge radius, or boundary-layer

thickness. The grid spacing is h. It is assumed that h  L, or in other words, the grid is sufficiently fine to

resolve the physical features of the flow. Making the substitutions Eqs. (2) into Eqs. (1) and keeping only

the term of O() yields the linearization of the inviscid terms

ujuˆi,j + uˆjui,j +

1

ρ

pˆ,i −

pˆ

ρc2

p,i

ρ

+

sˆ

cp

p,i

ρ

, (3a)

ρc2uˆj,j + γpˆuj,j + ujpˆ,j + p,juˆj , (3b)

ujsˆ,j + uˆjs,j . (3c)

The scaling of the variables is

ui, uˆi ∼ U p, pˆ ∼ ρU2s, sˆ ∼ cpM2.

3
Because only the short-wavelength components of the error are of interest, the derivatives of the perturbation

quantities scale as

uˆi,j ∼

U

h

pˆ,j ∼

ρU2

h

sˆ,j ∼

cpM2

h

but the derivatives of the solution values scale as

ui,j ∼

U

L

p,j ∼

ρU2

L

s,j ∼

cpM2

L

.

Using the assumption that h  L, it is seen that terms involving the derivatives of solution quantities

can be neglected compared to the terms involving derivatives of the perturbation quantities. This gives the

principal linearization of the inviscid terms

ujuˆi,j +

1

ρ

pˆ,i, (4a)

ρc2uˆj,j + ujpˆi,j , (4b)

ujsˆ,j . (4c)

Also, note that the ratio of the second term to the first term of Eq. (4b) is O(M2). Thus for low Mach

number flows, M2  1, the pressure advection term in the continuity equation is subprincipal.

3.2 Principal linearization of the viscous terms

The right-hand sides of Eqs. (1c) and (1b), the entropy and continuity equations, consist of the heat diffusion

and viscous dissipation terms. The principal part of the linearization of each term is considered. The scalings

of the viscosity and temperature are

μ, μ

0 ∼

ρUL

Re

T, Tˆ ∼

U

2

cp

.

The Reynolds number Re is based on the relevant local velocity and length scales. Also, it is assumed

that ∂μ/∂T  μ/T, ∂

2μ/∂T2  μ/T2

, so the variations of the viscosity can be neglected when linearizing.

Examining the scaling of the deviatoric strain tensor and its derivatives gives

dij ∼

U

L

, dij,j ∼

U

L2

,

ˆdij ∼

U

h

,

ˆdij,j ∼

U

h

2

.

As a consequence, the dominant term in the linearization of the dissipation function Φ is

2μ

0

dii ˆdjj + 4μdij ˆdij .

Considering the thermal diffusion term next, note that the temperature derivatives scale as

T,i ∼

U

2

cpL

, T,ii ∼

U

2

cpL2

, Tˆ

,i ∼

U

2

cph

, Tˆ

,ii ∼

U

2

cph

2

,

from which it follows that the dominant term in the heat-flux linearization is

μcp

Pr

Tˆ

,ii.

Using the equation of state to express this in terms of sˆ and pˆ gives

ν

Pr 

ρc2

γ − 1

sˆ,ii

cp

+ pˆ,ii

.

4
Examination of the relative magnitude of the two terms yields

2μ

0

diidˆ

jj + 4μdijdˆij ∼

ρU3L

Re

1

Lh,

μcp

Pr

Tˆ

,ii ∼

ρU3L

RePr

1

h

2

.

It is seen that that thermal diffusion term is larger by O(Pr−1L/h), so the dissipation term is subprincipal

and can be neglected in the linearization.

Linearizing the viscous stress terms τij,j/ρ on the right-hand side of the momentum equation (1a) and

comparing magnitudes of the terms yields

1

ρ



μ

0

dˆ

jj,i + 2μdˆij,j

∼

μ

ρ

U

h

2

,

all other terms of being at least O(h/L) smaller.

In summary, the principal part of the viscous terms in the Navier-Stokes equations are

1

ρ



μ

0 ˆdjj,i + 2μ ˆdij,j

, (5a)

(γ − 1) ν

Pr 

ρc2

γ − 1

sˆ,ii

cp

+ pˆ,ii

, (5b)

(γ − 1) cp

ρc2

ν

Pr 

ρc2

γ − 1

sˆ,ii

cp

+ pˆ,ii. (5c)

in the momentum, continuity, and entropy equations, respectively.

3.3 The complete principal linearization

Having found the principal linearization of the inviscid and viscous terms in Eqs. (4) and (5), the principal

linearization of the full equations is found by comparing the relative magnitude of the terms in each equation.

Starting with the entropy equation, a comparison of the magnitudes of the terms in Eqs. (4c) and (5c)

yields the relative balance of advection and thermal diffusion

U

cpM2

h

∼

cp

ρc2

μ

Pr

U

2

h

2

.

This can be restated as the requirement RehPr = O(1) for both terms to be principal. If the Prandtl number

is O(1), then the cell Reynolds number must be O(1) for the viscous terms in the entropy equation to be

principal.

Examination of Eqs. (4a) and (5a) gives

U

2

h

∼

μ

ρ

U

h

2

as the dominant inertial and shear stress balance, or Reh = O(1).

Finally, the continuity equation terms from Eqs. (4b) and (5b) are comparable when

ρc2U

h

∼

μ

Pr

U

2

h

2

.

The relative magnitude of the viscous terms compared to the inviscid terms is M2/(RehPr), which for

moderate Mach numbers corresponds to the condition that the cell Reynolds number is order unity. However,

5
for low Mach numbers the terms on the right-hand side of the continuity equation are subprincipal. Once

again, this is reflection that in the incompressible limit, Eq. (1b) reduces to the statement that ui,i = 0.

To summarize, the principal linearization of the Navier-Stokes equations is

ujuˆi,j +

1

ρ

pˆ,i =

1

ρ



μ

0 ˆdjj,i + 2μ ˆdij,j

, (6a)

ρc2uˆj,j + ujpˆ,j = (γ − 1) ν

Pr 

ρc2

γ − 1

sˆ,ii

cp

+ pˆ,ii, (6b)

ujsˆ,j = (γ − 1) cp

ρc2

ν

Pr 

ρc2

γ − 1

sˆ,ii

cp

+ pˆ,ii

. (6c)

3.4 Interpretation of the principal linearization

As stated in the beginning of this section, the motivation for finding the principal linearization is that,

because a relaxation method is effective at eliminating short wavelength error components, only those terms

of the linearized operator that contribute strongly at short wavelengths are necessary. This property may be

stated more formally. Let LFL be the full linearization of the operator, and L be the principal linearization.

The principal linearization satisfies the property

kI − L

−1LFLk = O



h

L



. (7)

In other words, the inverse of the principal linearization is a good approximation to the inverse of the full

linear operator at the scale of the grid spacing.

An important point to emphasize here is that the key assumption used to derive the principal linearization

is h  L, which implies that the flow is well-resolved on the grid. It was shown above that formally the

viscous terms are principal when Reh = O(1), which is quite restrictive. For example, a typical calculation

of a laminar boundary layer flow over a flat plate will have on the order of 20–30 grid points across the

boundary layer. The local Reynolds number in the middle of the boundary layer, based on the local flow

speed and the normal grid spacing, is ∼3–10. When solving the system of equations using a full multigrid

(FMG) process, the viscous region will be underresolved on the coarse grids. The principal linearization

will not be strictly valid until the finest grid or two are reached. Also implicit in the condition h  L is

that the region of interest is well-removed from the boundaries and any discontinuities, such as shocks. The

complications that arise near boundaries and discontinuies is discussed at length in Thomas et al. (2003).

For typical aerodynamic flows where the viscous effects are confined to thin shear layers, the variation

in length scales over the domain is quite large. Not all the terms will be principal in the different regions of

the flow. In practice, the variation in physical scales is handled by variation in grid spacing throughout the

computational domain. As a consequence, the principal linearization will differ throughout the computational

domain. In the next section, various limiting cases of the principal linearization for Eq. (6) are considered.

4 Factors of the Navier-Stokes equations

In Section 3 the two components of an efficient multigrid method are given: a relaxation method for smoothing the short-wavelength component of the error, and a coarse-grid operator to reduce the long-wavelength

error. For multigrid to be fully effective, it is necessary that the solution error is well-represented on the

coarse grid. This is the case for elliptic partial differential equations. However, for hyperbolic equations not

all smooth components of the error are visible to the coarse grid operator. Those components normal to

the advection direction are only represented at the level of the truncation error. The coarse-grid correction

for these contains only part of this smooth error component. Multigrid is not effective for these types of

equations.

6
The physics of fluid flow is a combination of advection and diffusion phenomena, and as such the NavierStokes equations are a combination of elliptic and hyperbolic-parabolic operators. Energy and vorticity

are carried along streamlines and transported transversely by diffusion. Conservation of mass represents a

balance between pressure and the divergence of the velocity field, which is elliptic behavior subsonically and

hyperbolic supersonically. To construct an efficient multigrid solver for the Navier-Stokes equations, it is

necessary to split the system into these component parts such that each partition may be treated in the

most appropriate way.

The principal linearization in Eq. (6) can be used to identify and isolate the various subsystems of the

Navier-Stokes equations. Introducing the notation Q = ui∂i and Qν = Q − ν∆ for convenience, the system

may be written in the operator form Lqˆ = r, or





Qν − (ν + ν

0

)∂i∇·

1

ρ

∂i 0

ρc2∇· Q − ν

γ−1

Pr ∆ −ν

ρc2

cpPr∆

0 −ν

cp(γ−1)

ρc2Pr ∆ Q −

ν

Pr∆









uˆi

pˆ

sˆ



 = r. (8)

Because L is a matrix of operators, the formal eigenvalues of L are the “eigen-operators” of the system, in

a generalization of characteristic analysis of hyperbolic equations. These eigenvalues are pseudo-differential

operators. Unfortunately, not all these are true differential operators. However, taking the formal determinant of L gives an expression whose factors are identifiable differential operators. Some of these factors are

the eigenvalues of L, while the remaining factors are products of those eigenvalues which do not correspond

to any real differential operator.

The determinant of Eq. (8) is

detL = Q

d−1

ν

n

Q −

ν

Pr

∆

 h

Q



Q − γ

ν

Pr

∆



− c

2∆

i

+

 ν

Pr

− 2ν − ν

0



Q



Q − γ

ν

Pr

∆



∆

o

, (9)

where d is the number of space dimensions. The repeated advection-diffusion factors correspond to the transport of vorticity in two and three dimensions, and to the conservation of helicity in three dimensions (Ta’asan,

1993). The factor in the curly braces corresponds to energy and mass conservation. This factor is called

the “core operator” by Brandt. Thermal diffusion couples temperature variations to pressure in Eq. (8).

However, there are some special cases in which the core operator can be factorized, which will be examined

below. Furthermore, some of these special cases will occur in regions of many flows of engineering interest.

For this reason, these limits can be used when relaxing the full Navier-Stokes equations in these regions.

4.1 Inviscid equations

As shown in Section 3, the viscous terms are principal only where the cell Reynolds number is O(1). For

most aerodynamic flows, viscous effects are significant only in thin shear layers, and the bulk of the flow is

inviscid. The inviscid limit of Eq. (8) is

Lν→0 =





Q

1

ρ

∂i 0

ρc2∇· Q 0

0 0 Q



 . (10)

The advection-diffusion operator reduces to pure advection, and factors out of the core operator. The

determinant of Eq. (10) is

detLν→0 = Q

d



Q

2 − c2∆



. (11)

Two of the advection factors correspond to vorticity and helicity conservation, as before. The additional

advection factor is the transport of entropy along streamlines, which can be seen by inspection of Eq. (10).

The remaining second-order operator is the full potential operator, which corresponds to conservation of

mass. The full potential operator connects the pressure field to the irrotational part of the velocity field.

7
Note that in the one-dimensional case, the full potential operator factors into two first-order operators,

corresponding to the left- and right-running acoustic characteristics.

Although the discussion concerns the principal linearization of the equations, there is a correspondence of

these terms in the full nonlinear case as well. When the flow is uniform at infinity, the vorticity remains zero

and the entropy is constant. In such a case, all the advection terms are satisfied identically. The solution

of the Euler equations is fully determined by the solution of the continuity equation alone. By writing the

velocity as the gradient of a scalar potential function, the continuity equation reduces to the full potential

operator acting on the potential function. Thus it is seen how the factors of the determinant of L in the

inviscid limit correspond directly to a well-known property of inviscid, compressible flow.

4.2 Low speed flow

Another important limit is at low speed. Specifically, when M2  1, then the pressure advection term in

the continuity equation is subprincipal. With fluctuations in entropy scaling as given in Section 3.1, the

entropy diffusion term in the continuity equation is also subprincipal. In the entropy equation, the diffusion

of pressure and entropy are of comparable magnitude, and remain principal. The resulting reduced form of

the principal linearization is

LM→0 =





Qν − (ν + ν

0

)∂i∇·

1

ρ

∂i 0

ρc2∇· 0 0

0 −ν

cp(γ−1)

ρc2Pr ∆ Q −

ν

Pr∆



 . (12)

This corresponds to a nearly constant-density fluid, and occurs when the boundaries of the domain are

adiabatic.

If there is significant heating at the boundaries, then we can have that s, sˆ  cpM2. In this case, the

temperature (or entropy) fluctuations are decoupled from the pressure fluctuations. If M2  1 the pressure

advection-diffusion term in the continuity equation and the pressure diffusion term in the entropy equation

are subprincipal, leading to the principal linearization

LM→0 =





Qν − (ν + ν

0

)∂i∇·

1

ρ

∂i 0

ρc2∇· 0 −ν

ρc2

cpPr∆

0 0 Q − ν

Pr∆



 . (13)

For both Eqs. (12) and (13) the determinant of the operator is

detLM→0 = −c

2Qd−1

ν

(Q −

ν

Pr

∆)∆. (14)

The decoupling of the temperature and pressure variations is reflected by the core operator factoring into

an advection diffusion term and the Laplacian. The advection-diffusion factor acts on the entropy. The

Laplacian correponds to conservation of mass.

These two low speed limiting cases are distinct—one corresponds to incompressible flow, the other to

low Mach number but compressible flow. An important point to note here is that the determinant of the

operator is the same for both limits discussed above. By the criterion presented in Thomas et al. (2003),

one would neglect altogether the off-diagonal dissipative terms in the continuity and entropy equations in

the incompressible limit, resulting in the principal linearization

LM→0 =





Qν − (ν + ν

0

)∂i∇·

1

ρ

∂i 0

ρc2∇· 0 0

0 0 Q −

ν

Pr∆



 .

This form of the operator does not satisfy the property given by Eq. (7).

8
4.3 A special case of the Prandtl number

Examination of Eq. (9) shows the term in the curly braces can be factored when the Prandtl number is

1/(2 + μ

0/μ). Substituting into the equations yields

LPr=1/(2+μ0/μ) =





Qν − (ν + ν

0

)∂i∇·

1

ρ

∂i 0

ρc2∇· Q − (γ − 1)(2ν + ν

0

)∆ −(2ν + ν

0

)

ρc2

cp

∆

0 −(2ν + ν

0

)

cp(γ−1)

ρc2 ∆ Q − (2ν + ν

0

)∆



 (15)

and this has the determinant

detLPr=1/(2+μ0/μ) = Q

d−1

ν

[Q − (2ν + ν

0

) ∆] Q [Q − γ (2ν + ν

0

) ∆] − c

2∆

	

. (16)

The new advection factor that appears reflects the transport of entropy along a streamline. The remaining

term corresponds to continuity, but the diffusion terms reflect that the pressure is coupled to the entropy

through thermal diffusion.

In general, Stokes’ hypothesis is applicable which gives that μ

0 = −2μ/3, which means that this factorization applies when Pr = 3/4. Note that for air the Prandtl number ranges from 0.75 to 0.70 over the

temperature range −100 to 100 Celsius, and is approximately 0.72 for most atmospheric flight conditions.

Thus this factorization is very nearly exact for aerodynamic flows.

For convenience, Eq. (15) and its determinant Eq. (16) is written below for the case when the Stokes

hypothesis applies and Pr = 3/4:

LPr=3/4 =





Qν − 1

3

ν∂i∇·

1

ρ

∂i 0

ρc2∇· Q −

4

3

(γ − 1)ν∆ −

4

3

ρc2

cp

ν∆

0 −

4

3

cp(γ−1)

ρc2 ν∆ Q −

4

3

ν∆



 ,

detLPr=3/4 = Q

d−1

ν



Q −

4

3

ν∆

 

Q



Q − γ

4

3

ν∆



− c

2∆



.

5 Modes of the principal linearization

To construct an efficient relaxation scheme for the Navier-Stokes equation that makes effective use of the

factors, it is necessary to determine the derived variables upon which the factors are acting. The most

natural decomposition of the perturbation variables into modes corresponding to the factors of the operator

is





uˆi

pˆ

sˆ



 =





Usi

0

0



 +





∂i

−ρ(Q −

ν

Pr∆)

−(γ − 1) cpν

c

2Pr∆



 φ +





ν

cpPr∂i

0

1



 σ (17)

where Usiis a solenoidal velocity field, and σ and φ are scalar potentials. This is related to the classical

decomposition of a velocity field into solenoidal and irrotational parts. Applying the operator L to the above

decomposition yields

L





uˆi

pˆ

sˆ



 = Qν





Usi

0

0



 +





0

−ρQ Q − γ

ν

Pr∆



+ ρc2∆

0



 φ

+



Q −

ν

Pr

∆







ν

cpPr∂i

0

1



 σ

+

 ν

Pr

− 2ν − ν

0



∆





∂i



ν

cpPrσ + φ



0

0



 . (18)

9
Note that the last term vanishes for the special case of Subsection 4.3, i.e., Pr = 1/(2 + μ

0/μ). In this

particular case, three remaining terms correspond exactly to the factors of Eq. (16). Furthermore, in the

inviscid limit it is observed that φ is the velocity potential, and σ becomes the entropy. Most significantly,

the continuity equation depends on φ alone, and the entropy equation depends on σ alone. The coupling

between these modes in the general case appears in the momentum equation, and arises from the effect of

the dilatational viscous stresses.

The quantity Usiis an incompressible vorticity error transported along the mean flow. The quantity σ

is an isobaric temperature error transported by the mean flow. The appearance of the velocity component ν/(cpPr)∂iσ arises from the dilatation of the velocity field induced by heating of the fluid. The quantity φ is essentially a velocity potential. A change in φ corresponds to a change in velocity at a constant

total enthalpy, or total temperature. The continuity equation is the full potential equation modified by the

addition of a term corresponding to the absorption of acoustic energy by thermal diffusion. In the case

of Pr = 3/4, the velocity dilatation due to thermal diffusion is balanced exactly by the viscous action of the

bulk viscosity. For any Pr ∼ 1, the coupling is weak.

In the inviscid limit, these three modes are fully decoupled:

L





uˆi

pˆ

sˆ



 = Q





Usi

0

0



 +





0

−ρ



Q2 − c

2∆



0



 φ + Q





0

0

1



 σ. (19)

It is interesting to note that these modes appear in Kov ́asznay’sstudy of supersonic turbulent flow (Chu &

Kov ́asznay, 1958; Kov ́asznay, 1953). He developed the perturbation equations for small-amplitude, smallscale fluctuations being transported with a uniform mean velocity. In the present case, the short-wavelength

solution errors are analogous to the small-scale turbulent fluctuations, and these errors are transported along

the mean-flow approximate solution. Because the equations are linearized about an unconverged state, the

residuals act as apparent mass sources, body forces, and heat sources on the right-hand side. Kov ́asznay’s

work focuses on disturbances whose Reynolds number is much greater than one, whereas here the viscous

terms are principal when Reh = O(1). For higher Reynolds numbers, the artificial viscosity of the scheme

dominates the physical viscosity.

It is useful to compare the present decomposition to that presented in Thomas et al. (2003). They use

internal energy rather than entropy as their thermodynamic variable. Converting their formulation into the

present notation, and rescaling some terms, yields





uˆi

pˆ

sˆ



 =





I

ρ(ν + ν

0

)∇·

−

1

T

(ν + ν

0

)∇·



 Usi +





∂i

−ρ(Q − (2ν + ν

0

)∆)

1

T

(Q − (2ν + ν

0

)∆)



 φ +





0

0

1



 σ. (20)

Here, the velocity Usiis not assumed to be solenoidal. The quantity σ is an isobaric temperature mode, but

it does not induce any velocity dilatation. Finally, the quantity φ is not the usual velocity potential, but

corresponds to a velocity and pressure change at a constant static temperature, rather than total temperature.

The application of L to this set of variables yields

L





uˆi

pˆ

sˆ



 =





Qν

P∇·

−

ν+ν

0

T Q∇·



 Usi +





0

−ρQ (Q − (2ν + ν

0

)∆) + ρc2∆

1

T Q (Q − (2ν + ν

0

)∆)



 φ +





0

−

ρνc

2

cpPr∆

Q − ν

Pr∆



 σ (21)

where P = ρc2 + (μ + μ

0

)Q. Note that the energy and continuity equations are coupled for all values of Pr.

Also, if the velocity Usiis solenoidal, the tranport equation for this mode becomes fully decoupled from the

other modes. In the inviscid limit, this system becomes

L





uˆi

pˆ

sˆ



 =





Q

ρc2∇·

0



 Usi +





0

−ρ



Q2 − c

2∆



1

T Q2



 φ +





0

0

Q



 σ. (22)

10
6 Conclusions

In this paper, the compressible Navier-Stokes equations have been examined with the objective of identifying

the factors of the operator. The principle linearization of the systems has been developed by examining those

terms of the linearization that are most significant for the short wavelength error components. The factors

of the operator are seen to be identical to the modes identified by Kov ́asznay (1953) in his analysis of

compressible turbulent flow.

It must be emphasized that the decomposition of the Navier-Stokes equations into these factors is only

the first step in the construction of a multigrid algorithm. A factorizable discretization of the system—one

that preserves these factors at the discrete level—must be constructed. Such a discretization must also

be h-elliptic for an efficient relaxation scheme and multigrid solver to be constructed. Creating a scheme

that preserves all these properties is an extremely difficult task (see Sidilkover, 1999). Indeed, there is no a

priori guarantee that such a scheme exists. Nevertheless, the identification of the factors provides a useful

framework for construction of fast converging schemes.

References

Brandt, A. 1985. Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics. GMD-Studie 85,

GMD-FIT

Brandt, A., Yavneh, I. 1992. On multigrid solution of high-Reynolds incompressible entering flow. J. Comput.

Phys. 101:151–64

Chu, B.-T., Kov ́asznay, L. S. G. 1958. Non-linear interactions in a viscous heat-conducting compressible gas.

J. Fluid Mech. 3:494–514

Kov ́asznay, L. S. G. 1953. Turbulence in supersonic flow. J. Aeron. Sci. 20:657–74, 682

Roberts, T. W., Sidilkover, D., Thomas, J. L. 2000. Multigrid relaxation of a factorizable conservative

discretization of the compressible Euler equations. AIAA Paper 2000-2252

Roberts, T. W., Swanson, R. C., Sidilkover, D. 1999. An algorithm for ideal multigrid convergence for the

steady Euler equations. Comput. Fluids 28:427–42

Sidilkover, D. 1999. Factorizable scheme for the equation of fluid flow. ICASE Rep. 99-20

Ta’asan, S. 1993. Canonical forms of multidimensional steady inviscid flows. ICASE Rep. 93-34

Ta’asan, S. 1994. Canonical-variables multigrid for steady-state Euler equations. ICASE Rep. 94-14

Thomas, J. L., Diskin, B., Brandt, A. 1999. Distributed relaxation multigrid and defect correction applied

to the compressible Navier-Stokes equations. AIAA Paper 99-3334

Thomas, J. L., Diskin, B., Brandt, A. 2001. Textbook multigrid efficiency for the incompressible NavierStokes equations: high Reynolds number wakes and boundary layers. Comput. Fluids 30:853–74

Thomas, J. L., Diskin, B., Brandt, A. 2003. Textbook multigrid efficiency for fluid simulations. Annu. Rev.

Fluid Mech. 35:317–40

Swanson, R. C. 2001. Towards optimal multigrid efficiency for the Navier-Stokes equations. AIAA Paper

2001-2574

11
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing

data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or 

any other aspect of this collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate 

for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that 

notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently 

valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR'S ACRONYM(S) 

11. SPONSORING/MONITORING

REPORT NUMBER 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

ABSTRACT 

18. NUMBER 

OF 

PAGES 

19a. NAME OF RESPONSIBLE PERSON 

a. REPORT b. ABSTRACT c. THIS PAGE 

19b. TELEPHONE NUMBER (Include area code) 

Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std. Z39-18