Journal of Computational Physics 183, 165–186 (2002)

doi:10.1006/jcph.2002.7187

A Highly Accurate Technique for the Treatment

of Flow Equations at the Polar Axis
in Cylindrical Coordinates Using
Series Expansions
G. S. Constantinescu and S. K. Lele
Center for Turbulence Research, Department of Mechanical Engineering,
Stanford University, Stanford, California 94305-3030
E-mail: gconstan@ctr.stanford.edu and lele@leland.stanford.edu

Received August 13, 2001; revised August 5, 2002

Numerical methods for solving the flow equations in cylindrical or spherical coor-
dinates should be able to capture the behavior of the exact solution near the regions
where the particular form of the governing equations is singular. In this work we
focus on the treatment of these numerical singularities for finite-difference meth-
ods by reinterpreting the regularity conditions developed in the context of pseudo-
spectral methods. A generally applicable numerical method for treating the singular-
ities present at the polar axis, when nonaxisymmetric flows are solved in cylindrical
coordinates using highly accurate finite-difference schemes (e.g., Pade schemes) on
nonstaggered grids, is presented. Governing equations for the flow at the polar axis
are derived using series expansions near r = 0. The only information needed to cal-
culate the coefficients in these equations are the values of the flow variables and their
radial derivatives at the previous iteration (or time) level. These derivatives, which
are multivalued at the polar axis, are calculated without dropping the accuracy of
the numerical method using a mapping of the flow domain from (0, R) ∗ (0, 2) to
(−R, R) ∗ (0, ), where R is the radius of the computational domain. This allows
the radial derivatives to be evaluated using high-order differencing schemes (e.g.,
compact schemes) at points located on the polar axis. The accuracy of the method
is checked by comparison with the theoretical solution corresponding to a circu-
lar compressible forced jet in the regime of linear growth. The proposed technique
is illustrated by results from simulations of laminar-forced jets and turbulent com-
pressible jets using large eddy simulation (LES) methods. In terms of the general
robustness of the numerical method and smoothness of the solution close to the po-
lar axis, the present results compare very favorably to similar calculations in which
the equations are solved in Cartesian coordinates at the polar axis, or in which the
singularity is removed by employing a staggered mesh in the radial direction with-
out a mesh point at r = 0, following the method proposed recently by Mohseni and

165
0021-9991/02 $35.00

c 2002 Elsevier Science (USA)
All rights reserved.
166 CONSTANTINESCU AND LELE

Colonius [1]. Extension of the method described here for incompressible flows or
for any other set of equations that is solved on a nonstaggered mesh in cylindrical or
spherical coordinates with finite-difference schemes of various levels of accuracy is
immediate. c 2002 Elsevier Science (USA)
Key Words: coordinate system singularity; polar axis; turbulence; large eddy sim-
ulation method; cylindrical system; collocated grids; finite differences method.

1. BACKGROUND

The singularities at the centerline of a cylindrical coordinate system are due to the presence
of terms containing the factor 1/r , where r is the radial distance, in the equations governing
the flow. The flow field itself does not have any singularity at the polar axis. For axisymmetric
flows symmetry conditions may be used to remove these singularities, but in the simulation
of nonaxisymmetric flows this problem must be faced. Furthermore, as the computational
domain is defined as (0, 2) ∗ (0, R), one has to specify numerical boundary conditions at
r = 0, even if physically there is no boundary at the polar axis.
Several numerical methods have been proposed to address the singularity of the flow
equations in cylindrical or spherical coordinates. At first sight they seem to vary greatly
depending on whether a pseudo-spectral, finite-volume, or finite-difference framework is
adopted, but in fact there are many similarities among them. This is especially true for
the new finite-difference method on nonstaggered grids proposed here and the treatment
used with pseudo-spectral methods. There are also common elements with methods that
use l’Hopital’s rule [2] or that use a “shifted” distribution of points in the radial direction
[1] to eliminate the points on the polar axis.
The main idea behind using spectral methods in cylindrical coordinates is to seek an ap-
proximation using polynomial expansions in the radial direction that satisfy some regularity
conditions so as to insure a well-behaved solution near the polar axis. O’Sullivan and Breuer
[3] studied the growth of linear disturbances in pipe flow. They solved the Navier–Stokes
equations using Fourier transforms in both the streamwise and the azimuthal directions,
and Chebyshev polynomials in the radial direction. The singularity at the polar axis was
avoided by using information in wavenumber space related to the form of the series expan-
sions of the velocity and pressure variables at r = 0, and by mapping in the radial direction
to increase resolution near the polar axis. Orszag and Patera [4], in a related study, used a
similar spectral method but they employed parity relations for the Chebyshev coefficients
in the expansions to remove the coordinate singularity of the cylindrical system. Fully
developed pipe-flow DNS calculations in cylindrical coordinates using spectral methods
were reported in [5, 6]. The method in [5] used Jacobi polynomials as the expansion basis
in the radial direction close to the polar axis and Legendre–Lagrangian interpolants away
from the axis, while the latter used B-spline polynomials expansions. An adaptation of the
method proposed in [5] for spherical coordinates, where one is faced with essentially the
same problem at the two singularity axes corresponding to polar angles  = 0 and  = ,
can be found in [7]. Zhang et al. [5] were able to recast the governing system of equations
in the wavenumber space so that the singularities due to coefficients of the form 1/r and
1/r 2 were removed. This was essentially done by taking advantage of the particular form
of the series expansions in wavenumber space of the pressure and velocity components.
However, they still had to use a special form of Jacobi polynomials in the radial direction
 FLOW EQUATIONS AT THE POLAR AXIS 167

close to the origin in conjunction with l’Hopital’s rule to remove potential singularities
arising because of terms where both the numerator and the denominator tend to zero at
approximately the same rate near the polar axis. Loulou [6] followed a similar approach,
but the regularity conditions in terms of having a certain behavior in r near the origin in
the polynomial expansions of the Fourier modes were imposed as constraints on the co-
efficients of the B-splines. These regularity conditions are essentially the same ones we
are going to use in our method, but in the present work everything is formulated in the
physical space (r, , z) as opposed to wavenumber space (r, k , kz ). In the spectral methods
the essential information is the order of the leading term in the polynomial expansions
in r for each mode of the Fourier expansion in . No explicit equations at the polar axis
are derived, and there is no need to calculate the values of the coefficients in these ex-
pansions. Finally, a somewhat different approach was proposed in [8] in the context of
pseudo-spectral (interpolatory) methods. By use of the governing differential equations,
pole conditions were derived in [8] which were subsequently used as numerical boundary
conditions at the coordinate singularity. The advantage of the method is that these pole
conditions are developed in physical space and are easier to implement into a standard
pseudo-spectral method. The method was illustrated for solving a Poisson equation on the
unit disk.
Finite-volume and finite-difference methods are less accurate than spectral methods, but
they can handle complex geometrical configurations more easily. Finite-volume methods
for flows in cylindrical geometries were employed by Eggels et al. [9], Akselvoll and Moin
[10, 11], and Boersma et al. [12], among others. The main difficulty here is related to
the fact that the azimuthal and streamwise fluxes are not defined at the centerline. Due
to the definition of velocity components at the center of the volume cells, the singular-
ity for the azimuthal and streamwise fluxes is removed. However, if the control volume
closest to the centerline is not wedge shaped, the flux at the first control-volume surface
off the polar axis whose normal is oriented in the radial direction must be evaluated. In
one of the implementations of their numerical method, Eggels et al. [9] used first- and
second-order one-sided differences to extrapolate the radial velocity on the faces closest
to the centerline. A similar procedure was used to calculate the diffusive flux at these
locations. Loss of accuracy near r = 0 was found to be not very significant due to the
relatively smooth DNS velocity fields in that region for the pipe flow. They also tried
to use wedge-shaped elements at the centerline, the advantage being that the singular-
ity is automatically removed as the radial flux is identically equal to zero on the faces
“located” at the centerline. In this implementation the method was found to give unreal-
istically high values for the r.m.s. fluctuations near the centerline. However, one should
point out that Boersma et al. [12], using a similar method, did not observed such prob-
lems in their DNS simulations of incompressible jets. Akselvoll and Moin [10] used a
finite-volume method to solve for the turbulent flow in a coaxial jet combustor using LES.
They used regular-shaped elements near the centerline, with the closest faces located at
r/2 from the singularity axis. Provided the radial and azimuthal velocity components
u r and u  are well defined at the centerline for all ’s (they are multivalued at r = 0),
the calculation of the fluxes on these faces presents no special problem. They interpo-
lated the values for the radial velocity at every azimuthal location  on the centerline
using the corresponding values at r/2 for two azimuthal angles,  and  + , and ac-
counting for the change in sign for the radial velocity across the centerline in the limit
r → 0. A similar treatment was followed for the azimuthal velocity at the centerline that
168 CONSTANTINESCU AND LELE

was needed to compute some viscous fluxes in the momentum equation for the azimuthal
fluxes.
Verzicco and Orlandi [13] developed a second-order finite-difference scheme in cylindri-
cal coordinates. The main feature of their method was the introduction of a radial flux (r u r )
on a staggered grid. As for finite-volume methods, this is the only flux to be calculated right
at r = 0, and its value there is obviously equal to zero. The DNS study in [14] used this
algorithm to calculate the fully turbulent flow in pipes with rotating walls. Their method
seems to be very appealing in the context of second-order schemes, but extension to higher
order schemes is not straightforward.
Griffin et al. [2] used finite differences and L’Hopital’s rule to recast the governing
equations at the origin. For all terms involving coordinate singularities they calculated
the radial derivatives using one-sided second-order-accurate finite differences, while the
azimuthal derivatives for the variables that are multivalued at the origin were formally
evaluated by using the values at neighboring azimuthal locations and  = 2/N , where
N is the number of cells in the azimuthal direction. This resulted in a set of new equations
at the origin that do not contain any singularities. Formally this treatment was not very
rigorous, as all these values are physically located at the same point, but the second-order
method was found to produce smooth solutions near the polar axis at least in the framework
of the inviscid calculation of the I.C. engine that they considered. However, its use for DNS
or LES simulations using higher order schemes is doubtful, and in the next section we
discuss in a more rigorous way using series expansions the errors introduced by the use
of L’Hopital’s rule to remove the singularities at the polar axis. This will also answer the
question which they raised in their paper about whether or not a fixed limit exists for terms
like (1/r )∂ p/∂, which using l’Hopital’s rule they recast into ∂ 2 p/∂r ∂. We show that in
this case such a limit exits and can be written using Fourier series expansions of the pressure
in terms of the m = 1 and m = −1 modes (the pressure cross derivative is multivalued at
the origin). For this particular example the Fourier expansions of 1/r ∂ p/∂ and ∂ 2 p/∂r ∂
give the same limit; however, we also show other examples (e.g., the Laplacian operator
that is present in the viscous term) where these kinds of approximations using l’Hopital’s
rule will lead to a wrong result.
Finally, Mitchell et al. [15], Freund [16], and Boersma and Lele [17] used compact
finite differences on nonstaggered meshes in cylindrical coordinates to compute the radi-
ated jet noise of compressible jets. In these studies, the equations at the polar axis were
solved using a Cartesian coordinate system. The directions of the Cartesian system were
arbitrarily chosen, and the multivalued variables (u r , u  ) for the other directions were ob-
tained by rotation of the system at r = 0, where only modes m = 1 and m = −1 may
exist. As the present numerical scheme is very similar to the ones used in these studies,
more details can be found in the validation section. Another option for removing the sin-
gularities at the polar axis is to eliminate the points at r = 0 by distributing the points
in the radial direction starting with r/2 (for a uniform grid spacing in r ) and mapping
the domain (0, 2) ∗ (0, R) into (0, ) ∗ (−R, R) for evaluating the radial derivatives. This
ensures that no numerical boundary conditions have to be specified at the polar axis. The
radial derivative stencils will span the centerline without evaluation at r = 0. This is es-
sentially the method discussed in [1]. We use this idea but recast it for the case in which
we have points at the polar axis. As we derive a set of exact equations which are well
defined on the polar axis, we do not have to avoid the presence of such points at the
centerline.
 FLOW EQUATIONS AT THE POLAR AXIS 169

2. TREATMENT OF THE EQUATIONS AT THE POLAR AXIS

Let us suppose that the system of governing equations may be written as

∂Q
= RHS(Q), (1)
∂t
where, for instance, in the case of compressible three-dimensional flows, Q = ( u x ,  u r ,
 u  ,  , e),  is the density, e is the energy, and the right-hand-side term (RHS) contains
the usual operators in cylindrical coordinates associated with the continuity or the transport
equations of momentum and energy, including the terms associated with the subgrid scale
contributions in the case of LES simulations.
The numerical method is general enough that the particular form of the operators in the
RHS is not really important. The same is true if more equations are included in Eq. (1),
such as transport equations for passive scalars, or for turbulence quantities in the case
of RANS calculations. The only important difference between the different equations in
Eq. (1) is determined by whether the variable in the left hand side of Eq. (1) is single valued
or multivalued at the polar axis. Following Boyd [19], the most general expansion of any
function that is not multivalued at the polar axis (scalars, pressure, streamwise velocity) can
be written as
∞


F(r, ) = ( f m (r ) · cos(m) + gm (r ) · sin(m)), (2)
m=0

where f m and gm are polynomials in r that have mth-order zeroes at r = 0. If m is even, f m

and gm are both symmetric around r = 0 and their power series contain only even powers
of r , while if m is odd, f m and gm are both antisymmetric and accordingly their power series
contain only odd powers of r . A similar statement holds for r u r and r u  , which takes care
of the expansion form for multivalued variables at the origin. These conditions imply that
the most general series expansion of a single-valued quantity (S) at the polar axis can be
written as
∞
∞  ∞
∞ 





S(r, ) = r m
mn r 2n
· cos(m) + r m
mn r 2n
· sin(m), (3)
m=0 n=0 m=0 n=0

while the expressions for multivalued quantities (e.g., u r and u  ) assume the form
∞ ∞
∞ 
1



M(r, ) = A0n r +
2n
r m−1
Amn r 2n
· cos(m)
r n=1 m=1 n=0
∞
∞ 



+ r m−1 Bmn r 2n · sin(m). (4)
m=1 n=0

A sketch of the proof of Eq. (3) can be given, starting with the most general Fourier expansion
of S in the form of Eq. (2) and requiring that all terms be regular at the origin. The Fourier
expansion in terms of complex variables can be written as
∞


F(r, ) = ( fˆm (r )eim ), (5)
m=0
170 CONSTANTINESCU AND LELE

where fˆm (r ) is a polynomial in r with complex coefficients. Each term in the summation
can be written as
∞

∞

∞


fˆm (r )eim = cn (m)r n eim = r m eim cn (m)r n−m = w m cn (m)r n−m , (6)
n=0 n=0 n=0

where w = y + i z = r ei . As w m is regular at origin for all m ≥ 0, the requirement that

each term in Eq. (5) be regular implies cn (m) = 0 for n < m, as 1/r s is not regular at origin
for s > 0. Also r = ( (y 2 + z 2 ))s is regular only for even values of s. For odd values of s
s
√
the derivative of r s with respect to y or z will contain 1/ r , which is not analytic at origin.
The form of the series expansions for multivalued quantities (e.g., for u r ) can be deduced
by observing that u r can be written in terms of the Cartesian velocity components u y and
u z as

u r = u y cos() + u z sin(), (7)

where the form of u y and u z expansions is given by (3). By regrouping the terms and
using algebraic identities of the form 2 cos() cos(m) = cos((m − 1)) + cos((m + 1)),
one can easily recover Eq. (4).
As we previously mentioned, any scalar or Cartesian velocity component is uniquely
defined at the origin, so one can write

∂ S 
= 0. (8)
∂ r =0

This relation holds in particular for u z = u r sin() + u  cos(). So by taking the deriva-
tives with respect to  one can write
   
∂u r  ∂u  
0= − u   sin() + + u r  cos(). (9)
∂ r =0 ∂ r =0

As Eq. (9) should hold for any , one obtains

∂u r ∂u 
= u and = −u r at r = 0. (10)
∂ ∂
There is another important constraint on the general form of the series expansions for u r
and u  . If Ai(rj ) , Bi(rj ) , Ai() ()
j , and Bi j are the coefficients of the series expansions for u r and
u  in Eq. (4), the following relation holds for all i ≥ 1:

() (r ) () (r )
Ai0 = Bi0 and Bi0 = −Ai0 . (11)

This relation can be obtained by using the corresponding series expansions given by Eq. (3)
for u y and by Eq. (4) for u r and u  , respectively. These series expansions are plugged into

u y = u r cos() − u  sin(), (12)

and using algebraic identities similar to the ones employed to deduce the general form of
the series expansions for multivalued variables at the polar axis (Eq. (4)), one can see that
 FLOW EQUATIONS AT THE POLAR AXIS 171

the polynomial multiplying cos(2) in the RHS of Eq. (12) has the form

1  ()
B + A(r10) + r 2 (· · ·) + · · · . (13)
2 10

()
As this polynomial should have a double zero at the origin, one should require B10 = −A(r10) .
The other relations in Eq. (11) follow similarly from analysis of the form of polynomials
multiplying the other modes in the RHS of Eq. (12).
By calculating the derivatives with respect to  and r of the series expansions given by
Eqs. (3) and (4) for all operators present in the RHS of the governing equations Eq. (1),
and taking the limit r → 0, a new form of the governing equations that is valid at r = 0
is obtained. These are a set of exact equations at the polar axis, provided we can calculate
exactly the coefficients Amm , Bmn , mn , mn for all terms present in the RHS of Eq. (1).
However, for a system of PDEs with second-order radial derivatives, as is the case for the
Navier–Stokes equations, it is sufficient to calculate at most the coefficients whose indices
m and n vary between 0 and 2, as we show next. The final expressions in the RHS of Eq. (1)
are dependent on the particular set of equations that is solved (inviscid, laminar, turbulent,
compressible, incompressible, etc.). However, the present method is valid for any system of
equations that can be written in cylindrical (or spherical) coordinates in the form of Eq. (1).
Hence, reference is made to the Appendix for the final expression of the RHS corresponding
to the compressible turbulent flow equations that were used in the present flow calculations
to validate this method. In the following discussion we focus our attention on pointing out
some general features.
At this point it is relevant to look at the expressions for the series expansions of a single-
valued variable at the origin, let us say the streamwise velocity u x (coefficients in series
expansions are labeled mn (x)
and mn
(x)
), and for a multivalued one, the radial velocity u r
(r ) (r )
(coefficients are Amn and Bmn ), as well as their first and second radial derivatives, in the
limit r → 0:

(x)
u x = 00 , (14)
∂u x (x) (x)
= 10 cos() + 10 sin(), (15)
∂r
∂ 2u x (x) (x) (x)
= 201 + 220 cos(2) + 220 sin(2), (16)
∂r 2
u r = A(r10) cos() + B10
(r )
sin(), (17)
∂u r
= A(r10) + A(r20) cos(2) + B20 (r )
sin(2), (18)
∂r
∂ 2 ur
= 2A(r11) cos() + 2B11 (r )
sin() + 2A(r30) cos(3) + 2B30
(r )
sin(3). (19)
∂r 2

As expected the series expansions of scalar quantities (e.g., u x ) contain only the m =
0 mode, while those of u r and u  contain only the m = 1(cos()) and m = −1(sin())
modes. Meanwhile, the first derivative of a scalar quantity is multivalued, as it contains the
m = 1 and −1 modes, while the second derivative contains the m = 0, 2, and −2 modes.
In particular, this applies to the pressure, p, which is a single-valued scalar, but ∂ p/∂r is
multivalued at r = 0. If we go back to Eq. (1) for the u x variable, and if we take the limit
172 CONSTANTINESCU AND LELE

as r → 0, we obtain
(x)
∂u x ∂00
= = RHSu x , (20)
∂t ∂t

where RHSu x contains terms coming from the convective and viscous operators that would
obviously involve modes m = 0, 1, −1, 2, −2, 3, and −3 (see the Appendix). In order for
Eq. (20) to be valid, the coefficients resulting from the series expansions of the different
variables and their derivatives should be such that the coefficients multiplying all but the
m = 0 mode in RHSu x must be equal to zero. Indeed this is what happens when we sum
the contributions coming from all the terms contained in RHSu x and use Eq. (11). However,
individual terms may contain nonzero coefficients multiplying the m = 0 modes. In fact, it
can be proven that this property holds not only for the RHS of the equations describing the
transport of a single-valued variable at the polar axis as a whole, but also for independent
operators such as convective, viscous, dilatation, and viscous dissipation in the energy
equation. For instance, the dilatation operator, which is a scalar quantity and thus should
contain only the m = 0 mode, can be expressed as
  (x)
∂u x ∂u r 1 ∂u  ∂00
+ + ur + = + 2A(r01) , (21)
∂x ∂r r ∂ ∂x

where the streamwise derivative can be calculated with exactly the same method used (sixth-
order compact differences in the present work) for points situated away from the polar axis.
The same observation is true for the continuity equation, recast as an evolution equation for
the density:
   ( ) (x) 
∂ ∂ u x 1 ∂r u r 1 ∂ u  ∂00 00 (r ) ( ) (r ) ( ) (r ) ( )
=− + + =− + A10 10 + B10 10 + 2A01 00 .
∂t ∂x r ∂r r ∂ ∂x

(22)

In Eq. (22) the series coefficients for the density  are mn
( )
and mn
( )
. Another example is the
Laplacian operator (e.g., the viscous term in the u x -momentum equation, or more exactly
the term corresponding to the constant viscosity part):

(x)
∂ 2u x ∂ 2u x 1 ∂u x 1 ∂ 2u x ∂ 200 (x)
+ + + = + 401 . (23)
∂x2 ∂r 2 r ∂r r 2 ∂ 2 ∂x2

For the u r - and u  -momentum equations the RHS will contain only the m = −1 and
m = 1 modes once contributions from all the terms are added. These two equations are
identical in the limit r → 0; thus it is sufficient to calculate the limit for only one of the
equations, e.g.,

∂u r ∂ A(r10) cos() + B10

(r )
sin()
= = RHSur = RHSm=1 cos() + RHSm=−1 sin(). (24)
∂t ∂t

In fact, we will end up with two scalar equations, corresponding to the m = 1 and −1
modes, respectively, as one can see from the previous relation, the reason being that as
r → 0, u r can be obtained from u  by a counterclockwise rotation of /2. This means that
 FLOW EQUATIONS AT THE POLAR AXIS 173

if the radial component at the polar axis is given by Eq. (17), the azimuthal component can
be written as

u  = A() () (r ) (r )
10 cos() + B10 sin() = −A10 sin() + B10 cos(), (25)

where use was made of relation (11). Relation (23) can also be used to show an inherent
problem with methods in which l’Hopital’s rule is applied numerically to remove singular-
ities at the polar axis (e.g., [2]). In these methods the Laplacian will be calculated at the
polar axis by numerically evaluating the expression
  (x)
∂ 2u x ∂ 2u x ∂ 2u x ∂ 2 ∂ 2u x ∂ 200 (x) (x) (x)
+ + + = + 401 − 420 cos(2) − 420 sin(2),
∂x2 ∂r 2 ∂r 2 ∂r 2 ∂ 2 ∂x2
(26)

which introduces a finite error, as the coefficients of the m = 2 and m = −2 modes are now
present.
The last step needed to complete the presentation of the present method is to describe
how the asymptotic series coefficients that are needed to evaluate the RHS in Eq. (1) are
computed. We will show that all that is required to calculate these coefficients accurately
is to be able to estimate numerically the first- and second-order radial derivatives of all the
variables in RHS with the same order of accuracy as points away from the polar axis. In
particular, in our compressible flow solver we use sixth-order Pade schemes to calculate
the radial derivatives. To do this, the following algorithm, similar to the one used in [1], is
adopted. The computational domain is mapped at every x = constant, such as there is no need
to specify numerical boundary conditions at r = 0. The mapping function (r, ) → (r̂ , ) ˆ
is
r̂ = r 0 <  < , r̂ = −r  <  < 2,
ˆ =  for 0 < r < R, and ˆ =  −  for 0 < r < R. (27)

All variables are identical in both systems of coordinates for 0 <  < . The rules con-
cerning the sign changes for the different variables and operators in the subdomain defined
by 0 < r < R and  <  < 2 follow from Eq. (27):

d r̂ dr
u r̂ = =− = −u r ,
d d
r̂ d ˆ r d
u ˆ = =− = −u  , (28)
d d
∂ ∂ ∂ ∂
=− , = , and ŝ = s for any scalar variable s.
∂ r̂ ∂r ∂  ˆ ∂

As all the signs of scalar quantities and streamwise velocities are left unchanged by the
mapping, Eq. (28) is sufficient to determine the sign of all terms involving radial derivatives
in the RHS of Eq. (1) in the mapped domain. The radial derivatives are now taken from
−R to R, with r = 0 being a regular interior point instead of a “numerical” boundary
point. For example, before taking the radial derivative of (r u  u r ), which appears in the
u  -momentum equation, one has to multiply the values for  = (, 2) by −1, because
of the three changes of sign (for r, u  , and u r ) required by the mapping. One should point
out that even if a scalar quantity is single valued at origin, its radial derivative is not. Once
174 CONSTANTINESCU AND LELE

the radial derivatives are calculated, the inverse transformation of Eq. (28) is taken, so the
azimuthal and streamwise derivatives are estimated with the usual coordinate definition
(0, R) ∗ (0, 2).
Next, we expand on the calculation of the required coefficients in the asymptotic expan-
sions for u r , and its radial derivatives. All other variables are treated in a similar fashion
using the appropriate series expansions. Suppose that N + 1 is the number of points in
the azimuthal direction (with modes 0, ±1, ±2, . . . , ±N /2). According to Eq. (17), only
the m = 1 and −1 modes should be present in a series formed by these N values. We
will formally associate the N values with , going from 0(N = 1) to 2(N + 1). A fast
Fourier transform can, in principle, be applied and the coefficients of the m = ±1 modes
will be A(r10) and B10(r )
; the other coefficients should be very small. This can also be done for
the first and second radial derivatives using their asymptotic expansions in Eqs. (18) and
(19) to find A(r01) , A(r20) , B20
(r )
, and A(r11) , B11
(r )
, A(r30) , B30
(r )
, respectively. Once these coefficients
are calculated the limits of terms involving azimuthal and cross derivatives are also known
using Eqs. (4) and (11), while streamwise derivatives of these coefficients are evaluated
with the same scheme used for interior points. However, there is a more economical way
to calculate these coefficients. Supposing that N is divisible by 8, meaning that for every
element in the series there is another situated exactly /4 away, one can take advantage of
the properties of sin() and cos() functions to evaluate the coefficients in Eqs. (17)–(19).
For instance, using Eq. (17) for  and  + /2, one can obtain a system of two linear
equations with two unknowns (A(r10) and B10 (r )
):

u r () = A(r10) cos() + B10

(r )
sin(),
(29)
u r ( + /2) = −A(r10) sin() + B10
(r )
cos().

To eliminate the bias toward a certain direction, one can solve the above system for every
 = (2/N )(n − 1) with n = 1 to N and average the results to get final values for A(r10) and
(r )
B10 . The same kind of treatment can be applied to determine A(r01) , A(r20) , B20
(r )
by first getting
(r )
the coefficient A01 and then using Eq. (18) for  and  + /4 to calculate the remaining two
coefficients. The coefficients involved in the expression for the second derivative of u r can
be determined more easily if we try to obtain at the same time the coefficients appearing in
the expression for the second derivative of u  ; i.e.,

∂ 2u
= 2A() () (r ) (r )
11 cos() + 2B11 sin() + 2B30 cos(3) − 2A30 sin(3), (30)
∂r 2

where use was made of Eq. (11). By writing Eqs. (19) and (30) at ,  + /2,  + /4,  +
3/4 and solving three systems of two linear equations, one can determinate the values of
all coefficients similarly to the procedure employed above.

3. VALIDATION OF PROPOSED NUMERICAL ALGORITHM

The development of this new technique was driven by the interest of the authors in
calculating the turbulence dynamics and acoustic radiation of turbulent compressible round
jets using LES. This problem is especially relevant to the design of more effective jet engines
with reduced noise emissions. The quality of the jet-noise predictions is determined by the
accuracy (resolution) of the numerical method, and by the boundary condition (physical
 FLOW EQUATIONS AT THE POLAR AXIS 175

as well as numerical at the polar axis) treatment and the quality of the mesh. The general
numerical method is described in detail in [18], while the details of the implementation of
the dynamic LES model in the original DNS code can be found in [17].
The numerical method employs compact fourth- or sixth-order Pade schemes [20] for the
spatial derivatives in the nonhomogeneous directions, and Fourier spectral methods in the
homogeneous (azimuthal) direction. The number of modes is dropped near the polar axis so
that the CFL constraint will be determined by the radial (or axial) spacing. No clustering of
the points near r = 0 is needed. The solution is advanced in time using a four-step Runge–
Kutta method. Zonal boundary conditions with artificial damping are employed near the
inlet and outlet to absorb outgoing disturbances and to avoid spurious noise generation at
these boundaries. Nonreflecting boundary conditions are used at the lateral boundary r = R,
as well as damping terms in a layer close to this boundary.
Before discussing the laminar and LES simulations using different treatments at the
centerline, it is important to validate the proposed centerline algorithm by considering a
test case for which a theoretical solution exists and to investigate its effect on the expected
order of accuracy of the numerical solution. An inviscid compressible circular jet at a Mach
number Ma = U/c∞ = 0.9, with the mean velocity, density, and pressure distributions given
by
  
0.25R0
ū x (r ) = U 0.5 − 0.5 tanh (r/R0 − R0 /r ) ,
h
ū r (r ) = 0,
ū  (r ) = 0, (31)
¯ (r ) =  ∞ ,
p̄(r ) = 1/ ,

is forced at the fundamental frequency (temporal wavenumber) = 1.0278c∞ /R0 and

complex streamwise wavenumber k = kr + iki = (1.9275 − 0.6094i)(1/R0 ) corresponding
to m = 1 according to

q(x, r, ) = q̄(r ) + Real (q̂r (r ) + i q̂i (r ))ei(kx+m− t)

, (32)

where q is the instantaneous value of a generic variable, q̄ is its mean value defined in
(31), q̂r and q̂i are the real and imaginary parts of the scaled eigenvector associated with q,
= 10−4 is a chosen coefficient independent of the variable that determines the amplitude
of the disturbances, R0 is the initial radius of the jet in the inlet section, h = 0.1R0 is the
jet momentum thickness, c∞ is the speed of sound at ambient conditions, and = 1.4 is
the ratio of specific heats. All the eigenvectors are scaled with the maximum absolute value
of the eigenvector associated with the streamwise velocity, so the value of plays the
role of amplification factor for the disturbances. The values of , k, and the form of the
eigenvectors associated with u x , u r , u  ,  , and p for m = 1 and the mean flow profiles given
by (31) were determined from the inviscid linear stability theory using a shooting method.
The inlet boundary conditions in our calculations correspond to setting x = 0 in (32). The
calculations were started from conditions identical to the theoretical solution (32) at t = 0.
Given these initial and inflow boundary conditions, it is expected that during the initial phase
of the jet development in the linear growth regime, the numerical solution should follow
176 CONSTANTINESCU AND LELE

the theoretical solution if the maximum amplitude of the oscillations is small. As the radial
and streamwise dimensions of our domain were R = 6R0 and X = 10R0 , the maximum
streamwise velocity disturbance in the exit section will be close to e−ki X ∼ 0.045(U = 0.9).
The subdomain that will be used to estimate the agreement between the simulations, and
the theoretical solution has the dimensions R = 3R0 and X = 5R0 . We verified that in
this region, which does not contain any part of the lateral or exit damping regions, the
growth of the disturbances is indeed governed by linear stability theory. The time step is
t = 0.001R0 /U, corresponding approximately to CFL = 0.1. The small value of the time
step was chosen to eliminate as much as possible the temporal discretization errors for this
model problem. Typically the code is run at CFL ∼ 1.0.
The solutions corresponding to three simulations using three different grids containing
40 ∗ 26 ∗ 16, 80 ∗ 51 ∗ 32, and 160 ∗ 101 ∗ 64 points in the (x, r, ) directions are compared
after t = 12R0 /c∞ (the period is T = 2/ = 6.11R0 /c∞ ). Each subsequent level of coars-
ening is obtained by the elimination of every other point in all three directions, starting from
the finest mesh. The three grids have the characteristic grid size proportional to 0 , 0 /2,
and 0 /4, where 0 corresponds to the coarsest grid. The computed results show that the
agreement between the numerical solution and the theoretical solution (32) improves signif-
icantly as the grid gets finer. This is illustrated in Fig. 1, where the distributions of the radial

FIG. 1. Radial velocity contours in a plan situated at x = 5R0 from the inlet section of the inviscid forced jet
at time t = 12R0 /U. Radial velocity fields are shown with 15 equally spaced contours between −0.0006U and
0.0006U. (a) Theoretical solution via linear stability analysis; (b) fine-grid solution; (c) medium-grid solution;
(d) coarse-grid solution.
 FLOW EQUATIONS AT THE POLAR AXIS 177

FIG. 2. Variation of error with mesh size for an inviscid test case.

velocity component in a plane situated at x = 5R0 are shown for the theoretical solution
and the fine-, medium-, and coarse-grid solutions. In particular, the agreement between the
fine-grid solution and theoretical solution is excellent. A more quantitative assessment of
the level of agreement between the calculated solution on a grid of size  and the theoretical
solution can be obtained by looking at the error defined in the L 2 , L 1 , and L ∞ norms; i.e.,

N N
(2) (q − q)2 (1) |q − q| (∞)
e = 1
, e = 1
, e = max|q − q|, (33)
N N

where q is the calculated value, q is the theoretical value, and N is the number of points
in the test volume (0 < x < 5R0 , 0 < r < 3R0 , 0 <  < 2). The normalized error e (with
respect to the coarse-grid results) for the streamwise velocity component is plotted in Fig. 2
in a log–log scale for the coarse, medium, and fine grids. As in these simulations fourth-order
Pade schemes were used for the discretization of the spatial derivatives in the streamwise
and radial directions, we assume that the error in the numerical solutions obtained on grids
whose local grid size ratio is 2 should ideally scale as 2n = 16, where n = 4 is the expected
order of accuracy. The slope corresponding to n = 4 is shown in Fig. 2 as a continuous
line for comparison. One can see that, as expected, the slope of the error curve (which
corresponds to the order of accuracy of the numerical solution) defined with either norm
(see (33)) is close to 4. Also it is relevant to mention that if we do the same analysis for a
much smaller volume around the centerline (0 < r < 0.3R0 ) we obtain a very close result
for the estimated order of accuracy of the numerical solution near the centerline. So, one
can conclude that the errors near the centerline are not dominant.
In all the remaining calculations discussed here that correspond to laminar and turbulent
jets, the Reynolds number is Re = U(2R0 )/ = 36,000, the Mach number is Ma = 0.9, and
the computational grid consists of 192 ∗ 128 ∗ 64 points in the (x, r, ) directions, which is
about an order of magnitude coarser than the grid used by Freund [16] to calculate a similar
jet at Re = 3600 using DNS. The mean flow distribution at the inlet plane is assumed to
be a rounded top-hat profile with periodic disturbances in the streamwise direction. The
simulations are run at a CFL number of one.
In a second test case (laminar forced jet) no randomized forcing is applied and the LES
model is turned off. A sinusoidal perturbation at a Strouhal number St = 0.5, corresponding
178 CONSTANTINESCU AND LELE

FIG. 3. Total vorticity contours. Vorticity shown using 18 equally spaced contours between 0 and 2U/R0 .
(a) Laminar forced jet; (b) turbulent jet.

to the most amplified wave, is imposed on the streamwise velocity profile. The forced-jet
solution is quasiaxisymmetrical and laminar, as seen from the total vorticity contours in
Fig. 3a. The effects of the different treatments at the polar axis are investigated in Fig. 4,
where snapshots of the dilatation fields in the region x > 20R0 are shown for different
treatments at the polar axis. We choose the dilatation because this quantity is very sensitive
to the centerline treatment and also provides relevant information for the jet acoustics in the
near field. Case L1 corresponds to the method where the equations are solved in Cartesian
coordinates, case L2 corresponds to a grid with no points placed at the polar axis (method
described in [1]), and case L3 corresponds to the treatment using series expansions at the
centerline. In the third test case randomized azimuthal forcing is applied in the input plane
to trigger the three-dimensional instabilities for the LES calculation. All other parameters
of the simulations, as well as the computational mesh, are kept the same as those used in the
second test case. As seen in Fig. 3b the jet undergoes transition to turbulence at the end of

FIG. 4. Dilatation contours in the forced jet. (a) Cartesian coordinates; (b) method of Mohseni and Colonius;
(c) series expansions treatment.
 FLOW EQUATIONS AT THE POLAR AXIS 179

FIG. 5. Dilatation contours in an area situated after the end of the potential core for the turbulent jet. Dilatation
fields are shown with 18 equally spaced contours between −0.08U/R0 and 0.08U/R0 . (a) Method of Mohseni
and Colonius; (b) series expansions treatment.

the potential core situated around x = 15R0 . Results are shown in Fig. 5 for two simulations,
the first (T1) uses the method of Mohseni and Colonius [1], while in the second (T2) the
series expansion treatment is employed.
Considerable success was reported by Freund [16] for DNS simulations of jets at lower
Reynolds numbers using the present numerical method with the equations at the centerline
solved in Cartesian coordinates. As we are interested in simulating jets at very high Reynolds
numbers (two to three orders of magnitude higher than the ones at which DNS is presently
possible) the robustness of the method on coarser meshes is an obvious requirement. Our
simulation showed that even for a laminar forced jet (case L1) strong oscillations in the
dilatation field developed near the centerline in the form of two-delta waves. They are
evident in Fig. 4a in the region near the polar axis, starting with the streamwise (x ∼ 15R0 )
location where ring vortices begin to be shed from the jet. These point-to-point oscillations in
dilatation are on the order of 1.5R0 /U , which is more than one order of magnitude higher
that the values of 0.1R0 /U associated with the maximum dilatation inside the coherent
structures of the jet. These problems were also apparent in the vorticity fields near the
centerline. For case L1 the simulation did not diverge, but obviously the quality of the
solution, especially the sound information that can be collected from these fields, is poor.
When used to calculate a turbulent jet we were not able to get a well-behaved solution on the
relatively coarse mesh used in these simulations without substantial filtering of the solution
to remove the two-delta waves (every couple of time steps), and the quality of the solution
was poor. This shows that the robustness of this polar axis treatment deteriorates rapidly
on coarse meshes and for high Reynolds numbers, when strong nonlinear interactions are
present in the flow. Several features of this method are suspected to be at the origin of
the above-mentioned problems. First, the derivatives in the directions corresponding to the
Cartesian coordinates are evaluated using six-order explicit central differences instead of
180 CONSTANTINESCU AND LELE

sixth-order Pade schemes, while in the evaluation of the radial derivatives with compact
differences, numerical boundary conditions using one-sided differences of lower order have
to be formulated. Another source of errors may be associated with the bias introduced into
the solution by the arbitrary choice (N /2 possibilities) of the two perpendicular directions
after which the flow equations are solved at the centerline.
Both of these approximations are automatically removed when a staggered mesh in the
radial direction is employed in simulation L2. However, our results using the technique
proposed in [1] showed that in the absence of filtering, the solution would diverge after
approximately 1000 time steps. That was due to a very local instability situated at the
first point off the polar axis at a certain streamwise location. The dilatation field remained
fairly smooth at all other streamwise locations close to the polar axis, in contrast to the
results seen in case L1. However, when a sixth-order-accurate explicit filter ( f˜ = ( f i−3 −
6 f i−2 + 15 f i−1 + 44 f i + 15 f i+1 − 6 f i+2 + f i+3 )/64) was applied every 300 time steps
to smooth the velocity field, the dilatation contours shown in Fig. 4b remained smooth
and no unphysical oscillations were observed near the polar axis. The cutoff wavenumbers
corresponding to values of the filtering transfer function of 0.9 and 0.5 are = 0.49 and
= 0.71, respectively. Any of these values can be taken as representative for the cutoff
wavenumber for the filter (see also [20] for more details about the filter properties). The
behavior of the solution for the forced laminar jet seems to suggest that this method is
not very robust for strong nonlinear problems, which would limit their use for turbulent
calculations, especially on fairly coarse grids.
Finally, in case L3, where the new treatment using asymptotic series expansion was used
at the polar axis, the dilatation fields are smooth (Fig. 4c) and the flow features are very
similar to the results obtained with the method of Mohseni and Colonius, but no filtering
was necessary to keep numerical instabilities from forming near the polar axis.
Based on the results from the second (laminar jet) test case, we focus our attention for
turbulent calculations on the method proposed by Mohseni and Colonius [1] and on the
method based on series expansions. Filtering the solution every 30 time steps, despite being
sufficient for keeping the solution from diverging, does not remove the unphysical spurious
oscillations in the dilatation field that are observed very close to the polar axis at all stream-
wise locations for x > 15R0 (Fig. 5a). In fact the maximum values of the dilatation in these
elongated structures very close to the polar axis are around 1.2R0 /U , which is significantly
higher than the levels of dilatation recorded in the rest of the computational domain. That
is true even if the filter is applied more often. One should point out that this amount of
filtering was found to have nonnegligible effects when sound sources were calculated from
these fields. Thus, even though the method proposed in [1] gives comparable results with
the series-expansions based technique for laminar forced jets, in turbulent simulations, the
robustness and accuracy of the method in [1] appears to be inferior. The turbulent case
is much more difficult to handle because the level of nonlinear flow interactions is much
higher and jet structures are passing through the centerline. Convergence of the solution on
relatively coarse meshes can be obtained in our simulations only by filtering the solution
when the method of Mohseni and Colonius [1] is used.
In simulation T2, patches of high dilatation are observed in a region of about 10R0 in the
streamwise direction, starting at the end of the inviscid core. The position of these patches
is not very close to the polar axis (Fig. 5b), except when they are convected by the mean jet
motion through the polar axis. The maximum absolute value of the dilatation field in the
turbulent region is at all times around 0.2R0 /U , including at the polar axis. The short waves
 FLOW EQUATIONS AT THE POLAR AXIS 181

that are seen at the top of Fig. 5b are rather related to the high-aspect ratio x/r of the
present grid and are not a consequence of centerline instabilities. Away from the centerline
the dilatation field appears smoother in the case of T1 compared to that of T2, but this is
just an effect due to the use of the filtering in T1. Even so, spurious waves or unrealistically
high values of the dilatation (Fig. 5a, case T1) or vorticity at the first two to three points off
the centerline are not observed.
The level of the unphysical spurious dilatation oscillations present in the domain close
to the polar axis is greatly reduced compared to the case in which the method of Mohseni
and Colonius [1] was used, while filtering was practically eliminated. The sixth-order
explicit filter was applied every 500 time steps to eliminate the high-frequency waves that
form because of the high grid-stretching ratio. Overall there is a significant increase in the
robustness for the method using asymptotic series expansions compared to the other two
approaches. Moreover, because we are primarily interested in extracting sound information
from the near fields, elimination of filtering that may compromise substantially the quality
of the sound data is a key factor in assessing the different methods.

4. SUMMARY

In this paper a general method for handling the singularities arising at the polar axis of
the governing flow equations in cylindrical coordinates based on power series expansions
in the radial direction and Fourier series in the azimuthal direction was presented. Using the
most general form of these series expansions, we derived a new set of equations at the polar
axis by calculating the limit of the various operators appearing in the governing system
of equations (Eq. 1). The new method is computationally easy to implement and is less
expensive than solving the equations in Cartesian coordinates at the centerline. The method
was proven to be well defined in the sense that for any single-valued variable at the polar axis
(scalar quantity or streamwise velocity) in the RHS of the Navier–Stokes equations only
the m = 0 mode was shown to have a nonzero coefficient in the limit of r → 0, despite the
fact that individual terms had nonzero coefficients for modes ±2 and ±3. The momentum
equations for u r and u  are coupled in the limit r → 0, and only modes m = ±1 are left
in the RHS of the two momentum equations. Again, this is consistent with the asymptotic
behavior of the radial and azimuthal velocities at the origin (u r ( + /2) = u  ()).
The present algorithm avoids the loss in accuracy at the centerline, where most of the
finite-volume and finite-difference methods use some kind of one-sided differences to ap-
proximate the operators at the centerline. This is because in the context of the present
method we were able to calculate the radial derivatives at the origin with the same order of
accuracy as in the rest of the domain (fourth- or sixth-order compact differences) by using
a domain mapping so that the points on the polar axis become regular interior points. As
this is the only information needed to calculate the coefficients of the newly derived polar
equations, there is no loss in the overall accuracy of the method. The coefficients in the
asymptotic expansions are unique, so we do not have to choose any arbitrary directions,
as is the case when the equations are solved in Cartesian coordinates at origin. One of the
advantages is that the first point off the axis, where the various terms with singular behav-
ior have to be evaluated, will be situated at r as opposed to r/2 (as is the case in the
method of Mohseni and Colonius [1]), which may avoid an important source of errors or
instabilities (e.g., due to the nonlinear interactions of eddies near the polar axis in a DNS
or LES simulation).
182 CONSTANTINESCU AND LELE

The accuracy of the present centerline algorithm was tested for a simulation of a forced
circular inviscid jet. The eigenvalues and eigenvectors of the mean jet velocity profile
determined from inviscid linear stability theory were used to excite the inflow boundary
of our jet. During the initial phases of the jet development the numerical solution should
follow the theoretical solution for the forced jet in the linear growth regime. We were able to
show that the numerical solutions converges toward the theoretical solution with increasing
mesh density and that the accuracy of the overall calculation is not affected by the proposed
centerline treatment; in other words, the spatial order of accuracy of the numerical solution
is close to the order of the numerical method.
The robustness of the proposed approach was tested successfully by comparing results
for a deterministically forced jet with similar calculations using the same numerical method
but with the equations solved at the centerline in Cartesian coordinates, or without any
points at r = 0 using the method described in [1]. Finally, the present technique was shown
to give improved results over the method of Mohseni and Colonius [1] for the simulation
of a compressible jet using LES on relatively coarse meshes, in which the flow interactions
taking place near the polar axis are very important, and where a flow solution that is
not contaminated by numerical artifacts originating at the singularity axis is required to
accurately calculate the jet-noise sources.
One should point out that the algorithm presented here to deal with the polar axis sin-
gularities is not restricted to the use of the present numerical method, which uses Pade
schemes to evaluate the derivatives in the radial and streamwise directions, or to solving the
compressible flow equations. Rather it can be adapted in a straightforward way to handle
the system of equations governing magnetohydrodynamics, combustion, and so forth. As
spherical coordinates are locally cylindrical coordinates near the two singularity axes, the
above method is also directly applicable to numerical methods that solve the flow equations
in spherical coordinates. In this case the governing equations are different but the derivation
of the set of equations valid at the two singularity axes is similar.

APPENDIX

The mass conservation, momentum, and equations for the compressible flow equations
in cylindrical coordinates are

 
∂ ∂ u x 1 ∂r u r 1 ∂ u 
=− + + , (A1)
∂t ∂x r ∂r r ∂
 
∂ u x ∂ u x u x 1 ∂r u x u r 1 ∂ u x u  ∂p
=− + + − + Vx , (A2)
∂t ∂x r ∂r r ∂ ∂x
 
∂ u r ∂ u r u x 1 ∂r u r u r 1 ∂ u r u   u 2 ∂p
=− + + − − + Vr ,
∂t ∂x r ∂r r ∂ r ∂r
  (A3)
∂ u  ∂ u  u x 1 ∂r u  u r 1 ∂ u  u   ur u  1 ∂p
=− + + + − + V ,
∂t ∂x r ∂r r ∂ r r ∂
 
∂e ∂(e + p)u x 1 ∂r (e + p)u r 1 ∂(e + p)u 
=− + +
∂t ∂x r ∂r r ∂
+ Ve + u x Vx + u r Vr + u  V + . (A4)
 FLOW EQUATIONS AT THE POLAR AXIS 183

The irreversible viscous dissipation, , is given by

∂u x ∂u r ∂u  ∂u x ∂u r ∂u  x ∂u x
= x x + xr + x + xr + rr + r +
∂x ∂x ∂x ∂r ∂r ∂r r ∂
r  ∂u r  ∂u  u r  u  r 
+ + + − . (A5)
r ∂ r ∂ r r

The viscous terms appearing in the momentum equations are

∂x x 1 ∂r xr 1 ∂x

Vx = + + , (A6)
∂x r ∂r r ∂
∂xr 1 ∂r rr 1 ∂r  
Vr = + + − ,
∂x r ∂r r ∂ r
(A7)
∂x 1 ∂r r 1 ∂ r 
V = + + + ,
∂x r ∂r r ∂ r

while the viscous flux term in the energy equation is

     
∂ ∂T 1 ∂ ∂T 1 ∂ 1 ∂T
Ve = + r + . (A8)
∂ x Pr ∂ x r ∂r Pr ∂r r ∂ r Pr ∂

In Eq. (A8) T is the temperature, Pr is the Prandtl number, and = (T ) is the shear
viscosity. The viscous stresses in cylindrical coordinates are
 
∂u x
x x = 2 + ,
∂x
 
∂u r ∂u x
xr = + ,
∂x ∂r
 
∂u  1 ∂u x
x = + ,
∂x r ∂
  (A9)
∂u r
rr = 2 + ,
∂r
 
∂u  1 ∂u r u
r  = + − ,
∂r r ∂ r
   
1 ∂u  ur
 = 2 + + .
r ∂ r

The coefficient in (A9) is defined as

B 2
= − , (A10)
3

where B is the bulk viscosity ( B/ = 0.6 for air according to [21]). The dilatation of the
velocity field, , is given by
 
∂u x ∂u r 1 ∂u 
= + + ur + . (A11)
∂x ∂r r ∂
184 CONSTANTINESCU AND LELE

In the following the series expansions coefficients introduced in Eq. (3) for a scalar variable
s are i(s) (s)
j , i j (for u x , s = x; for  u x , s =  x), while those introduced in Eq. (4) for multi-
valued variables (u r , u  ) are Ai(s)j , Bi(s)
j , with s = r and s = . The derivative of the viscosity
(T ) is calculated as

∂ ∂ ∂T ∂T
= = (T ) . (A12)
∂ xi ∂ T ∂ xi ∂ xi

If the appropriate series expansions are introduced into Eqs. (A1) through (A8) and A(11)
and the limit r → 0 is taken, a new set of equations valid at the polar axis is obtained; i.e.,

 ( ) (x) 
∂ ∂00 00 (r ) ( ) (r ) ( ) (r ) ( )
=− + A10 10 + B10 10 + 2A01 00 , (B1)
∂t ∂x
 ( x) (x)   ( p) 
∂ u x ∂00 00 (r ) ( x) (r ) ( x) (r ) ( x) ∂00
=− + A10 10 + B10 10 + 2A01 00 − + Vx , (B2)
∂t ∂x ∂x
 ( ) (x) (r )
∂ u r ∂00 00 A10
( )
∂  (x) B (r ) 
=− cos() + 00 00 10 sin() + A(r10) cos() + B10 (r )
sin()
∂t ∂x ∂x
 ( ) ( ) (r ) ( ) 
∗ 10 A(r10) + 10 B10 + 00 3A(r01) ∗ A(r10) cos() + B10 (r )
sin()
 
− A() (r ) (r )
01 ∗ B10 cos() − A10 sin() + cos() ∗ A10 A20 + B10 B20
(r ) (r ) (r ) (r )


 (r ) (r )   ( p) ( p)
+ sin() ∗ A(r10) B20 (r )
− B10 A20 − 10 cos() + 10 sin() + Vr , (B3)
 (x) (e+ p) 
∂e ∂00 00 (e+ p) (r ) (e+ p) (e+ p)
=− + A(r10)10 + B10 10 + 2A(r01)00 + Ve + 00(x)
Vx
∂t ∂x
  (r )    (r ) 
∂ ∂ A10 (r ) ∂ ∂ B10
+ A(r10) + 10 (x)
+ B10 + 10 (x)
∂x ∂x ∂x ∂x
 (x) 
 (T ) (r ) (T ) (r ) ∂10 (r ) ∂10 (x)
(r )
+ 10 A10 + 10 B10 + A10 + B10
∂x ∂x
 () (r ) (r )  (r ) () (r ) 
+ B11 + 3A11 A10 + 3B11 − A11 B10
 (r ) ()  (r ) (r )  (r ) (T ) (r )
+ 5A11 − B11 A(r10) + 5B11 + A() 11 B10 + 2A01 10 A10 + 10 B10
(T ) (r )

 (T ) (r ) (T ) (r )  (T ) (r ) (r ) 
+ 2 10 A20 + 10 B20 A(r10) + 10 B20 − 10 (T ) (r )
A20 B10 + , (B4)

(x) 
∂00 (x)
∂00 
= 2 + +2 A(r01) + 4 A(r01) A(r01) + A(r20) A(r20) + B20(r ) (r )
B20
∂x ∂x
 (r ) 
∂ A10 ∂ A(r10) ∂ B10
(r ) (r )
∂ B10 (x) ∂ A10
(r ) (r )
(x) ∂ B10 (x) (x) (x) (x)
+ + + 210 + 210 +10 10 + 10 10 ,
∂x ∂x ∂x ∂x ∂x ∂x
(B5)
    
∂ ∂ (x) (x) ∂ A(r )
Vx = 2 00 + + 401 + 2 01
∂x ∂x ∂x
 (r ) (r ) 
(T ) ∂ A10 (T ) ∂ B10 (T ) (x) (T ) (x)
+  · 10 + 10 + 10 10 + 10 10 , (B6)
∂x ∂x
 FLOW EQUATIONS AT THE POLAR AXIS 185

     (r ) 
∂ ∂ A(r10) (x) ∂ ∂ B10 (x)
Vr = cos() · + 10 + sin() · + 10
∂x ∂x ∂x ∂x
 (x) 
 (T ) (T ) ∂10 (x)
∂10
+ 10 cos() + 10 sin() + cos() + sin()
∂x ∂x
 ()  (r )   (r )
+ B11 + 3A(r11) cos() + 3B11 − A()11 sin() + 5A11 − B11 ()
cos()
 () (r )   
+ A11 + 5B11 sin() + 2 A(r01) 10 (T )
cos() + 10 (T )
sin()
 (T ) (r ) (T ) (r )  (T ) (r ) (T ) (r ) 
+ 2 10 A20 + 10 B20 cos() + 10 B20 − 10 A20 sin() , (B7)
 (T ) 
∂ ∂00 (T )   (T ) (T ) (T ) (T )
Ve = + 4 01 +   + 10 10 , (B8)
∂ x Pr ∂ x Pr Pr 10 10
(x)
∂00
= + 2A(r01) . (B9)
∂x

In deducing the final form of Eqs. (B1) through (B8) repeated use was made of Eq. (11).
Equations (B3) and (B7) were obtained by using the expressions for the radial momentum
equations and its viscous term (the first equation in (A3) and (A7)). As we pointed out
in the paper the expressions obtained using the azimuthal momentum equation and its
corresponding viscous term (the third equation in (A3) and (A7)) are redundant in the sense
that the two independent equations corresponding to the m = 1 and m = −1 modes will be
exactly the same as the ones deduced using only the expressions of the radial momentum
equation and its viscous term. Nevertheless, if this method is applied to another system of
equations, doing the calculations for both the radial and azimuthal directions may serve as
a check of the calculation. For the LES simulations there are additional terms arising from
the subgrid-scale modeling. For instance, the new term in the energy equation is

     
∂ ∂T 1 ∂ ∂T 1 ∂ 1 ∂T
VSGS =  + r +  , (A13)
∂x ∂x r ∂r ∂r r ∂ r ∂

where  =  t /Prt , t is the eddy viscosity, and Prt is the turbulent Prandtl number. As 
is a function of space that cannot be related directly to the temperature, as was , series
expansions should be used for this coefficient when the limit is taken in Eq. (A13). The
final form will be

 (T ) 
∂ () ∂00 () (T )  () (T ) () (T )
VSGS = 00 + 400 01 + 10 10 + 10 10 . (B10)
∂x ∂x

Note that  is evaluated directly in physical space using Lilly’s contraction in the present
implementation of the dynamical model (as opposed of calculating first Prt ), and its series
expansion is calculated in the same way as for the eddy viscosity t .

ACKNOWLEDGMENT

This work was supported in part by Grant NA G2-1373 from the NASA Ames Research Center.
186 CONSTANTINESCU AND LELE

REFERENCES

1. K. Mohseni and T. Colonius, Numerical treatment of polar coordinate singularities, J. Comput. Phys. 157,
787 (2000).
2. M. D. Griffin, E. Jones, and J. D. Anderson, A computational fluid dynamic technique valid at the centerline
for non-axisymmetric problems in cylindrical coordinates, J. Comput. Phys. 30, 352 (1979).
3. P. L. O’Sullivan and K. S. Breuer, Transient growth in circular pipe flow. I. Linear disturbances, Phys. Fluids
6(11), 3643 (1994).
4. S. A. Orszag and A. T. Patera, Secondary instability of wall-bounded shear flows, J. Fluid Mech. 128, 347
(1983).
5. Y. Zhang, A. Gandgi, A. G. Tomboulides, and S. A. Orszag, Simulation of pipe flow, in Symposium on
Application of Direct and Large Eddy Simulation to Transition and Turbulence (AGARD Conf., Crete,
Greece, 1994). Vol. CP-551, pp. 17.1–17.9.
6. P. Loulou, Direct Numerical Simulation of Incompressible Flow Using a B-Spline Spectral Method, Ph.D.
thesis (Department of Aeronautics and Astronautics, Stanford University, 1996).
7. A. G. Tomboulides, S. A. Orszag, and G. E. Karniadakis, Direct and Large-Eddy Simulations of Axisymmetric
Wakes, AIAA Paper 93-0546 (1993).
8. W. Huang and D. M. Sloan, Pole condition for singular problems: The pseudospectral approximation,
J. Comput. Phys. 107, 254 (1993).
9. J. C. M. Eggels, F. Unger, M. H. Weiss, J. Westerweel, R. J. Adrian, R. Friedrich, and F. T. M. Nieuwstadt,
Fully developed turbulent pipe flow: A comparison between direct numerical simulation and experiment,
J. Fluid Mech. 268, 175 (1994).
10. K. Akselvoll and P. Moin, Large Eddy Simulation of Turbulent Confined Coannular Jets and Turbulent Flow
over a Backward Facing Step, Report TF-63 (Department of Mechanical Engineering, Stanford University,
1995).
11. K. Akselvoll and P. Moin, An efficient method for temporal integration of the Navier–Stokes equations in
confined axisymmetric geometries, J. Comput. Phys. 125, 454 (1996).
12. B. J. Boersma, G. Brethouwer, and F. T. M. Nieuwstadt, A numerical investigation on the effect of the flow
conditions on the self-similar region of a round jet, Phys. Fluids 10(4), 899 (1998).
13. R. Verzicco and P. Orlandi, A finite-difference scheme for three-dimensional incompressible flows in cylin-
drical coordinates, J. Comput. Phys. 123, 402 (1996).
14. P. Orlandi and M. Fatica, Direct simulations of turbulent flow in a pipe rotating about its axis, J. Fluid Mech.
343, 43 (1997).
15. B. E. Mitchell, S. K. Lele, and P. Moin, Direct computation of the sound generated by an axisymmetric jet,
AIAA J. 35(10), 1574 (1997).
16. J. B. Freund, Acoustic Sources in a Turbulent Jet: A Direct Numerical Simulation Study, AIAA Paper 99-1858
(1999).
17. B. J. Boersma and S. K. Lele, Large eddy simulation of compressible turbulent jets, in Center for Turbulence
Research Annual Research Briefs (Dept. of Mechanical Engineering, Stanford Univ., Stanford, CA, 1999),
pp. 365–374.
18. J. B. Freund, P. Moin, and S. K. Lele, Compressibility Effects in a Turbulent Annular Mixing Layer, Report
TF-72 (Department of Mechanical Engineering, Stanford University, 1997).
19. J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springler-Verlag, Berlin/Heidelberg, 1989).
20. S. K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103, 16 (1992).
21. P. A. Thompson, Compressible Fluid Dynamics (McGraw–Hill, New York, 1991).