Open access

A Preconditioned Arbitrary Mach Number Scheme Applied to Rotating Machinery

Written By

Chunhua Sheng

Published: January 1st, 2010

DOI: 10.5772/7111

Chapter metrics overview

2,995 Chapter Downloads

View Full Metrics

1. Introduction

It is well known that the compressible flow equations face difficulties at low Mach numbers due to the large ratio of the acoustic and convective time scales, which leads to an ill-conditioned system when solving low-speed or incompressible flows. The time-dependent system of the Euler and Navier-Stokes equations exhibits stiffness that is strongly dependent on the Mach number and the Reynolds number. In this regard, Briley et al. (Briley et al., 1983) first introduced the preconditioning method using a simple constant preconditioning matrix added to a non-dimensional form of the isoenergetic equations. This generally improved convergence for a test case with reference Mach number Mr = 0.05 using the ADI factorization scheme in primitive variables. However, when applying this preconditioning to rotating flows in either fixed or rotating frame formulation, Sheng, et al. (Sheng & Wang, 2006; Wang & Sheng, 2005) observed the instability of the scheme due to the large variation of rotating speeds across the computational domain. Furthermore, it was found that the convergence of the preconditioned equations is sensitive to the selection of the reference Mach number, especially in rotating flows with a wide range of radial speeds and physical time scales. It was later proved using the Fourier footprint analysis (Wang & Sheng, 2005) that the eigensystem of the compressible governing equations can be significantly affected by both free stream and rotating speeds in rotating flows. A modified preconditioning scheme was thus proposed (Sheng & Wang, 2006; Wang & Sheng, 2005), in which both the global reference Mach number and the rotating Mach number are considered in the formulation of the preconditioning matrix. In general, this modified preconditioning scheme has improved the convergence and accuracy of compressible flows in subsonic, transonic and supersonic Mach number regimes. In this study, the modified preconditioning is further investigated and validated for predicting incompressible viscous flows in rotating machinery, such as a marine propeller P5168.

One of the most important characters for marine propellers is the cavitation observed in high speed flows, which is of vital importance because of the damage of metal surfaces and degradation of performance of lifting surfaces. It is also a source of high-frequency noise in connection with acoustic detection of ships and submarines. Cavitation would take place when the local pressure drops to the vapor pressure. Therefore, accurate prediction of the velocity and pressure field is essential for understanding the process of cavitation inception and improving the hydrodynamic performance of the marine propeller. Since the condition for cavitation inception is related to the tip-vortex location, strength, convection, cavitation position on the propeller surface, and the minimum pressure location. The present study will investigate these important flow characteristics using modified preconditioning scheme, and compare the computed hydrodynamics characters of P5168 with the LDV measurement (Chesnakas & Jessup, 1998).


2. Preconditioned method

2.1. Governing equations

The governing equations for solving the hydrodynamic flow in the present study are based on an unsteady three-dimensional Reynolds-averaged Navier-Stokes equations, coupled with the one-equation turbulent model of Spalart and Allmaras (Spalart & Allmaras, 1991). The unsteady compressible equations are cast in a rotating Cartesian coordinates with a rotational speed of Ω, and are written in terms of absolute velocities. It represents a system of conservation laws for a control volume that relates the rate of change of a vector of node-based average state variables to the flux through the volume surface. The following reference quantities are chosen to normalize the flow variables and equations: density, ρ r ; velocity, U rr ; temperature, T r ; pressure, ρ r U r 2 ; length, L r ; time, L r /U r ; energy and enthalpy, C p T r . The non-dimensional conservative form of the governing equations can be written in integral form as

d d t Ω Q d V + Ω F ( Q ) n ^ d A = Ω S ( Q ) d V E1

The above equation can also be written in terms of primitive variables by introducing a transformation matrix M = Q / q , where Q=(ρ, ρu, ρv, ρw, ρe t ), and q = (ρ, u, v, w, p). The governing equations in primitive variables are

M d d t Ω q d V + Ω F ( q ) n ^ d A = Ω S ( q ) d V E2

The projected normal flux vector to the face of a control volume is given as

F n ^ = [ ρ Θ ρ u Θ + n x p ( n x τ x x + n y τ x y + n z τ x z ) ρ v Θ + n y p ( n x τ y x + n y τ y y + n z τ y z ) ρ w Θ + n z p ( n x τ z x + n y τ z y + n z τ z z ) ρ h t Θ n t E c p ( γ 1 ) M r 2 ( u ( τ n ^ ) ) + Q n ] E3

where E c = (γ−1)M r 2 is an Eckert number. Θ is the relative velocity in the normal direction of the control surface in a rotating frame.Let Ω×r be the speed of the control volume at position r, Θ is given as

Θ = u n ^ × r n ^ = n x u + n y v + n z w + n t E4

The viscous stresses in the flux vector are given as

τ x x = μ 2 3 ( 2 u x v y w z )
τ x y = τ y x = μ 2 3 ( u y + v x ) E5
τ y y = μ 2 3 ( 2 v y u x w z )
τ x z = τ z x = μ 2 3 ( u z + w x ) E6
τ z z = μ 2 3 ( 2 w z u x v y )
τ y z = τ z y = μ 2 3 ( v z + w y ) E7

and the heat flux term is

Q n = μ Re P r ( n ^ T ) E8

where μ is the viscosity, Pr is the Prandt number, and Re is the Reynolds number based on the reference values.

Since the governing equations are cast in the rotating coordinate system, body forces (Coriolis and centrifugal forces) are added in the momentum equations as a source term due to the non-inertial reference frame, which can be expressed as

S = [ 0 ρ w Ω y ρ v Ω z ρ u Ω z ρ w Ω x ρ v Ω x ρ u Ω y 0 ] E9

The formulation of the governing equations in a rotating frame provides advantages for faster convergence in unsteady flows and a more stable solution without the unsteady effects of the grid motion. Furthermore, the above formulation can be easily restored into the ones in the fixed frame by setting the rotational speed Ω to zero in Eq. (6), making the body force source term disappear from the governing equations. All other relations, including the equation of state, remain the same under the current formulation. Using the absolute velocity components to represent the formulation in rotating frame offers great flexibility for solving the governing equations in both unsteady mode (fixed frame) and steady mode (rotating frame).

2.2. Preconditioned equations

In order to solve flows from essentially incompressible regime to supersonic Mach number, a preconditioning matrix Γq -1 is introduced into the time derivative terms of the compressible governing equations, which can be expressed as

M Γ q 1 t Ω q d V + Ω F n ^ d A = Ω S d V E10

where Γ q 1 is a constant diagonal matrix that only depends on the reference Mach number M r

Γ q 1 = d i a g [ 1, 1, 1, 1, 1 / β ( M r ) ] E11

The selection of preconditioning parameter β(M r ) will be discussed in the following section.

The finite–volume scheme needs the approximation of both inviscid and viscous fluxes on the face of a control volume. In order to develop characteristic–based flux formulation, the eigensystem of the preconditioned equations needs to be derived.

The preconditioned Eq. (7) can also be expressed as

q t + a Γ q x = Γ q M 1 S E12

where a Γ = Γ q a = Γ q M 1 A M is the system matrix for preconditioned equations, and A is the flux Jacobian matrix with respect to conservative variables Q. The system matrix for the three–dimensional Euler equations is

a Γ = [ Θ ρ n x ρ n y ρ n z 0 0 Θ 0 0 n x ρ 0 0 Θ 0 n y ρ 0 0 0 Θ n z ρ 0 β ρ c 2 n x β ρ c 2 n y β ρ c 2 n z β Θ ] E13

The eigenvalues of a Γ are

λ ( 1 ) = λ ( 2 ) = λ ( 3 ) = Θ , λ ( 4 ) = Θ β + + σ ,
λ ( 5 ) = Θ β + σ E14

where σ = ( Θ β ) 2 + β c 2 , β ± = ( 1 ± β ) / 2 . A nonsingular set of right eigenvectors for a Γ is given by

R = [ n x n y n z ρ ( Θ β + σ ) β c 2 ρ ( Θ β σ ) β c 2 0 n z n y n x n x n z 0 n x n y n y n y n x 0 n z n z 0 0 0 ρ ( Θ β σ ) ρ ( Θ β + σ ) ] E15

where Ψ = ρ Θ ( 1 β ) / β c 2 , and

c ± = Θ β ± ± σ 2 σ E16

2.3. Preconditioning parameter

From the above equations, it is seen that the eigenvalues and eginevectors of the system matrix have been modified due to the introduction of a preconditioned matrix. If the preconditioning parameter β is selected as one, no preconditioning is applied to the governing equations, and thus the original set of eigenvalues for compressible flowsis restored.

Briley, et al. (Briley et al., 2003) suggested that a global constant preconditioning parameter is evaluated based on the reference Mach number as:

β = min ( 1, M r 2 ) E17

where the reference Mach number should represent the global flow property for the problem being solved. For a steady and non–rotating flow, the above choice provided well–behaved eigenvalues in all characteristic waves, and thus improved the convergence and accuracy for solutions up to nearly incompressible Mach number (Briley et al., 2003). However, it was found that such a preconditioning may cause instability for low Mach number rotating flows, such as the case of a marine propeller P5168 presented here. In rotating machinery, the flow filed is complicated with secondary swirling flow structures due to the high-speed rotating elements. It is thus rational to consider the speed of rotating element in determining the preconditioning parameter. In this work the modified preconditioning parameter is employed for flows aimed at unsteady flows in rotating machinery:

β = min ( 1, M r 2 + α M Ω 2 ) E18

where ΜΩ is evaluated based on the rotational reference speed, and α is a scalar factor to tune the rotational Mach number in the preconditioned scheme. Numerical tests conducted in the previous work (Sheng & Wang, 2006; Wang & Sheng, 2005) indicated that the value of α should be in the range of 0.5~1.0 in order to achieve the optimal convergence and accuracy in preconditioned solutions. In the following, an analysis of the Fourier footprint for the eigensystem of the two-dimensional Euler equations provides evidences for the modification of the preconditioning parameter in low Mach number rotating flows.

2.4. Fourier footprint analysis

Fourier analysis is a useful tool for inspecting analytical and numerical solutions of the partial differential equations of mathematical physics. It is used to investigate the accuracy related to the approximation of hyperbolic equations and characteristics of all frequency modes of the discretized equation system (Vichnevetsky & Bowles, 2005). Euler equation is a pure hyperbolic system in which Fourier analysis finds its full power. The Euler Fourier footprint reveals the characters of error propagation and damping of all modes produced by the numerical discretization. The Fourier footprint of a given spatial discretization scheme is obtained by inserting the discrete Fourier transforms in the semi-discrete equation of the finite-volume governing equations. Detailed analytical formulation of the eigenvalues of Fourier footprint of the preconditioned Euler equation can be found in the reference (Wang, 2005).

In terms of Von Neumann stability analysis (O’Brien et al., 1951), the convective part of numerical solutions of hyperbolic equations corresponds to the extent of the Fourier footprints along the imaginary axis, and the dissipation component corresponds to the extent of the negative real axis. In Fig. 1 (a) for a low Mach number non-rotating flow at Mr = 0.01, there exists large disparity among phase velocities of convective waves (in red and black colors) and acoustic waves (in blue and green colors), which is visually indicated by the extension of the imaginary part of Fourier footprints. The group velocity, which describes the propagation of transient errors, is approximately equal to the amplitude of the sinusoid that defines the imaginary component in the plots. Thus, the damping and propagation of transient errors are much less effective for the convective waves, especially for low frequency modes. For this non-rotating flow, the preconditioned system with β = Μ r 2 equalized the real parts of the complex eigenvalues among all wave families, and rescaled the real and imaginary components to the same order of magnitude as indicated in Fig. 1 (b).

Figure 2 presents the Euler Fourier footprint for the same low Mach number flow, but with a rotating Mach number Μ Ω = 0.1. Without preconditioning, large disparity still remains among the convective and acoustic wave speeds in the same order of non-rotating flow, see Fig.2 (a). When the original preconditioning parameter is used ( β = Μ r 2 ) for Mach number relative to a rotating frame of Μ Ω = 0.1, as shown in Fig. 2(b), the real part extensions of Fourier footprints in all wave families are properly rescaled, and the damping effects of the high frequency error modes are well balanced between acoustic and convective waves. However, there are more footprints collapsed to the origin so that neither damping nor propagation of those low frequency modes is possible, and the system will not converge.

Figure 1.

Fourier Footprint of 1st Order Euler Operator for Non-Rotating Flow at M r = 0.01 .

Figure 2.

Fourier Footprint of 1st Order Euler Operator for Rotating Flow at M r =0.01, M Ω =0.1

This may explain the reason why the original preconditioning showed instability for low Mach number rotating flows. Using a modified a modified preconditioning parameter as β Μ r 2 Ω 2 proposed here, the extensions along the negative real axis for all wave modes are clustered to the same order of magnitude, as indicated in Fig. 2(c). In addition, the imaginary parts of the complex eigenvalues are mostly away from the negative real axis. The propagation properties are thus maintained at the ordinary levels so that the low frequency error modes can rapidly propagate out of the computational domain. The convergence of the Euler system is thus improved.


3. Discretization scheme

The preconditioned governing equations are discretized using a finite volume formulation. This is done by dividing the computational domain into a finite number of elements, from which control volumes are formed that surround each vertex in the mesh. The flow variables are stored at the vertices of the element. The preconditioned governing Eq. (7) or (9) is then numerically integrated over the closed boundaries of control volumes surrounding each node, and the preconditioning matrix Γ q 1 is applied to the convective (inviscid) term only through modifying the system matrix. The diffusive (viscous) term of the governing equations remains unchanged. The discretized form of the Navier-Stoles governing equations over a control volume can be written as:

M Γ q 1 Δ q i Δ t Δ V i + j = 1 n j ( ( F n ) j Δ A j ) = S i Δ V i E19

where F is the numerical flux vector evaluated on the face of the control volume. The subscript i denotes the vertex, and j=(1,… nj) denotes the jth surface of the control volume V surrounding the vertex i.

3.1. Convective flux evaluation

The convective numerical flux of preconditioned governing equations is calculated using in a formula similar to Roe’s flux, as proposed by Briley, et al. (Briley et al., 2003). Consider the surface j of a control volume whose left and right states are denoted by L and R, the numerical flux projected onto the surface j can be written as:

F n = 1 2 ( F ( q L ) + F ( q R ) ) 1 2 M Γ q 1 | a Γ | ( q R q L ) E20

where F ( q L ) and F ( q R ) are the convective fluxes evaluated at the left and right state of the surface. | a Γ | is the preconditioned system matrix with positive eigenvalues. The system matrix is evaluated using the following averaged primitive variables q ^ on the surface of a control volume:

q ^ = 1 2 ( q L + q R ) E21

Since the eigensystem and the flux are evaluated using the above averaged variables instead of Roe variables, the current approach is considered to be an extension of the Roe flux approximation.

In the above Eq. (16) quantities q L and q R are the values of the primitive variables on the left and right side of the face of the control volume. For a first-order accurate differencing, quantities q L and q R are set equal to the value at the vertices lying on either side of the face. For the second-order scheme, these face values (j) are computed with a Taylor series expansion about the central node (i) of the control volume

q j = q i + q i Δ r E22

where r is the vector that extends from the central node to the midpoint of each edge, and q is the first derivatives of the primitive variables at node and is evaluated with an un-weighted least-squares procedure (Anderson et al., 1995).

3.2. Diffusive flux evaluation

Since the local state can be found from the reconstructed solutions at the vertex of each edge, the essence of computing the viscous fluxes is finding the gradients of the velocity vector and the temperature at the midpoints of the interface. A directional derivative approach for viscous gradients has been developed by Hyams (Hyams, 2002). In the present work, a normal derivatibbe approach is presented.

Introducing two unit tangential vectors m ^ , l ^ that are orthogonal to the unit normal vector n ^ on the surface of a control volume

m ^ = ( n y n z n m , n z n x n m , n x n y n m ) T E23
l ^ = ( n y n s 1 n m , n y n s 1 n m , n z n s 1 n m ) T E24


n m = ( n y n z ) 2 + ( n z n x ) 2 + ( n x n y ) 2 E25
n s = n x + n y + n z E26

It is easy to verify that unit vectors are orthogonal to each other. The components of any gradient can be expressed in terms of derivatives in m ^ , l ^ , n ^ directions.

ϕ x = ϕ n n x + ϕ m m x + ϕ l l x E27
ϕ y = ϕ n n y + ϕ m m y + ϕ l l y E28
ϕ z = ϕ n n z + ϕ m m z + ϕ l l z E29

The derivatives in the normal direction n ^ can be approximated by

ϕ n ϕ 2 ϕ 1 d s E30

where ϕ 1 and ϕ 2 are the values of the interest at two vertices of an edge, and ds is the length of an edge. It is seen that this formula guaranties positivity of the operator, and second–order accuracy can be achieved if the direction of the edge is in concurrence with the normal direction of the face. The tangential derivatives are evaluated by projecting the averaging value of the two gradients at the nodes of an edge in m ^ , l ^ directions, respectively,

ϕ m 1 2 ( ϕ 1 + ϕ 2 ) m ^ E31
ϕ l 1 2 ( ϕ 1 + ϕ 2 ) l ^ E32

where ϕ 1 and ϕ 2 are the gradients evaluated at two nodes of the edge using a weighted least square procedure (Anderson, et al., 1995). The key point here is to make sure that these tangential derivatives are positive as well. Because a positive scheme satisfies the local maximum principle such that the solution is bounded by its surrounding neighboring vertices, a limiter is incorporated into the gradients in Eq. (22) to ensure the positivity of the tangential derivatives. A positive scheme for viscous discretization is thus achieved.

3.3. Temporal discretization

The temporal discretization is considered with an implicit Euler method for the system of governing equations. A general expression is available for the temporal discretization, which is given as:

M Γ q 1 Δ q i n θ 1 + θ Δ q i n 1 Δ t + 1 Δ V i j N ( 0 ) ( F n ) j n + 1 Δ A j = S i n + 1 E33

where Δ q i n = q i n + 1 q i n , Δ q i n 1 = q i n q i n 1 . Δt is the time step increment between steps n and n+1, and ΔV is the control volume. F denotes both inviscid and viscous fluxes. Subscript i and j denote the indices of the central node and the dual faces that compose of the control volume. A first order accuracy time Euler implicit scheme is given by the choice θ = 0 . Correspondingly, a second order time accurate Euler implicit scheme is given by θ = 1 .

3.4. Newton’s method

The above nonlinear system of Eq. (26) is solved by Newton’s method, which yields a linear system of equations at each time step. The use of Newton’s method allows relatively large time step to be used in numerical computations, a considerable savings compared to the explicit method, in particular fro the unsteady time-accurate problems. Denote

N ( q i n + 1 ) =M Γ q 1 Δ q i n θ 1 + θ Δ q i n 1 Δ t + 1 Δ V i j N ( 0 ) ( F n ) j n + 1 Δ A j S i n + 1 = 0 E34

The Newton’s method solves the following equation:

N ( q i n + 1, m ) ( q i n + 1, m + 1 q i n + 1, m ) =-N ( q i n + 1, m ) E35

where m = 1,2,3... are the Newton steps, with q i n + 1,0 = q i n . N ( q i n + 1, m ) is the Jacobian matrix of Eq. (27), that is

N ( q i n + 1 ) =M Γ q 1 I Δ t + 1 Δ V i j N ( 0 ) ( F n ) j n + 1 q i n + 1 Δ A j S i n + 1 q i n + 1 E36

In the above Eq. (29), the first term on the right hand side is the contribution from the unsteady time derivative of q, the second term is the contribution from the steady state residual (both inviscid and viscous) of the governing equations, and the third comes from the source term in the rotating reference frame. The flux Jacobian of the steady state residual can be evaluated by an approximate method. Consider the matrix M Γ q 1 | a Γ | in Eq. (16) as a local constant, the approximate flux Jacobian for the inviscid flux can be written as

j N ( 0 ) ( F n ) j n + 1 q n + 1 Δ A j = j N ( 0 ) ( F n ) j n + 1 q L n + 1 Δ A j + j N ( 0 ) ( F n ) j n + 1 q R n + 1 Δ A j E37
( F n ) j n + 1 q L n + 1 = 1 2 ( ( F ( q L ) ) n + 1 q L n + 1 + M ^ Γ q 1 | a Γ | ) E38
( F n ) j n + 1 q R n + 1 = 1 2 ( ( F ( q R ) ) n + 1 q R n + 1 M ^ Γ q 1 | a Γ | ) E39

A similar method is used to evaluate the viscous flux Jacobian by taking the derivatives of discretized viscous flux formulation.

3.5. Gauss-Seidel relaxation

The nonlinear system of equations is linearized by Newton’s method, which results in a sparse system of equations at each time. The solution of the sparse system of equations is obtained by a relaxation scheme in which Δ q n + 1, m is obtained through a sequence of iterations, { Δ q n + 1, m } i , which converge to Δ q n + 1, m . There are several variations of classic relaxation procedures that have been used in the past for solving this linear system of equations (Coirier, 1994; Anderson, 1992; Batina, 1990). In this work, a symmetric implicit Gauss–Seidel procedure as described in (Batina, 1990) is used. To clarify the scheme, the system matrix is first written as a linear combination of matrices representing the diagonal, upper triangular, and lower triangular parts at each time step

[ A ] = [ D + U + L ] E40


[D]=M Γ q 1 I Δ t + 1 Δ V j N ( 0 ) F j n + 1 q n + 1 Δ A j S q n + 1 E41
[U]= 1 Δ V j N U ( 0 ) F j n + 1 q n + 1 Δ A j E42
[L]= 1 Δ V j N L ( 0 ) F j n + 1 q n + 1 Δ A j E43
Let { R n + 1 } be the vector of unsteady residuals, and { Δ q n + 1 } represents the change in the dependent variables, the symmetric Gauss–Seidel relaxation can be written as the following two–step process:
[ L + D ] { Δ q n + 1 , m } i + 1 / 2 + [ U ] { Δ q n + 1 , m } i = { R n + 1 , m } E44
[ D + U ] { Δ q n + 1 , m } i + 1 + [ L ] { Δ q n + 1 , m } i + 1 / 2 = { R n + 1 , m } E45

In the forward pass, { Δ q } i + 1 / 2 are obtained with the previously updated { Δ q } i , which were set to zero at the initial stage. In the backward pass, { Δ q } i + 1 are obtained with the most recent value of { Δ q } i + 1 / 2 from the previous pass. Normally eight symmetric Gauss–Seidel sub-iterations are adequate at each time step.


4. Turbulence model

Several turbulence closure models have been developed for high-Reynolds number flows. One of them is the Spalart-Allmaras one-equation turbulence model (Spalart & Allmaras, 1991), which formulates a transport equation for the turbulence Reynolds number, and is the only one used in the present study. From the original Spalart-Allmaras formulation, a transport equation for the turbulence Reynolds number v ˜ = ρ μ t / f v 2 (working variable) is expressed as

ν ˜ t + U ν ˜ = c b 1 ( f r 1 f t 2 ) S ˜ ν ˜ 1 Re [ c w 1 f w c b 1 k 2 f t 2 ] [ ν ˜ d ] 2 E46
+ 1 σ Re { [ ( ν ˜ + ( 1 + c b 2 ) ν ˜ ) ν ˜ ] c b 2 ν ˜ ( ν ˜ ) } = 0 E47

where the first and second terms on the right hand side of Eq. (34) is the source production, and the third term is the turbulence dissipation. Constants in Eq. (34) can be found in the original paper (Spalart & Allmaras, 1991).


5. Characteristic-based boundary conditions

In order to compute internal and external flows at arbitrary Mach numbers, four commonly used characteristic–based boundary conditions have been developed (Sheng & Wang, 2003), including far field boundary, inflow boundary based on total conditions or mass flow, outflow boundary based on back pressure, and impermeable wall. The characteristic variables are derived based on preconditioned system matrix in order to be consistent with interior point treatment. To solve the governing equations in a rotating frame, a source term is included which should be implicitly integrated into the characteristic boundary conditions. In the following, characteristic-based boundary conditions are provided for completeness.

5.1. Far field

The preconditioned system matrix can be diagonalized as a Γ = R Λ Γ R 1 , where Λ Γ is a diagonal matrix made up of the eigenvalues of a Γ , and R and R -1 are matrices whose columns and rows are the right and left eigenvectors of a Γ , respectively. Multiplying both side of the Eq. (9) by R 0 1 gives the following characteristic relations in one dimensional space

W t + Λ Γ W x = R 1 Γ q M 1 S E48

where W = R 0 1 q is the characteristic variables with the subscript 0 denoting a reference value. The following relation is valid along the characteristic line:

d x / d t = λ Γ E49
d W d t = R 0 1 Γ q M 1 S E50
W b = W r + R 0 1 Γ q M 1 S b Δ t E51

where subscript b denotes a point on the boundary, and subscript r denotes an appropriate reference point depending on the sign of the eigenvalues on the boundary. Five characteristic variables for the three-dimensional preconditioned system are

w 1 = ( ρ p / β c 0 2 ) n x ( Θ n t ) ψ 0 n x w n y + v n z E52
w 2 = ( ρ p / β c 0 2 ) n y ( Θ n t ) ψ 0 n y u n z + w n x E53
w 3 = ( ρ p / β c 0 2 ) n z ( Θ n t ) ψ 0 n z v n x + u n y E54
w 4 = p 2 ρ 0 σ 0 + ( Θ n t ) c 0 + E55
w 5 = p 2 ρ 0 σ 0 + ( Θ n t ) c 0 E56


ψ 0 = ρ 0 Θ 0 ( 1 β ) / β c 0 2 E57
c 0 ± = Θ β ± ± σ 0 2 σ 0 E58

The direction of propagation of five characteristic waves is determined by the sign of each associated eigenvalue. Since all boundary faces have an positive outward pointing normal with the current unstructured grid topology, a negative sign of the eigenvalue means that the related characteristic wave propagates from the free stream to the boundary nodes, and the reference value is the free stream value; while a positive sign of the eigenvalue means that the characteristic wave propagates from the interior of the computational domain to the boundary, so the reference value should be taken from the interior node. The primitive variables at the boundary are then computed by solving the following equations:

q b = R 0 W r + Γ q M 1 S Δ t E59

5.2. Inflow

For internal flows in rotating machineries, total conditions (total pressure Pt, total temperature Tt) are usually given at the inlet, along with the flow angle ( ϕ x , ϕ y , ϕ z ) . Denote the inflow velocity as U b , then the following relations exist

ρ b = γ M r 2 P t T t ( 1 γ 1 2 T t M r 2 U b 2 ) 1 γ 1 E60
p b = P t ( 1 γ 1 2 T t M r 2 U b 2 ) γ γ 1 E61
u b = U b cos ϕ x E62
v b = U b cos ϕ y E63
w b = U b cos ϕ z E64

Since there are six unknowns (ρ b , p b , u b , v b , w b , U b ) out of the above five equations, one more relation is required to solve the equations. For the subsonic inflow, all eigenvalues except λ ( 4 ) are negative. The fourth characteristic relation that propagates from the interior field to the boundary is used

w b ( 4 ) = w i ( 4 ) + [ R 0 1 Γ q M 1 S b Δ t ] ( 4 ) E65
p b 2 ρ 0 σ 0 + ( Θ n t ) b c 0 + = p i 2 ρ 0 σ 0 + ( Θ n t ) i c 0 + + [ R 0 1 Γ q M 1 S b Δ t ] ( 4 ) E66

An alternate way to specify the inflow boundary conditions is to use the total temperature Tt, mass flow rate m ˙ , and the inflow angle α = ( ϕ x , ϕ y , ϕ z ) . Equations (43), (44), and (48) are replaced by the following three equations:

m ˙ = ρ b U b E67
T t = T b + γ 1 2 M r 2 U b 2 E68
p b = ρ b T b γ M r 2 E69

All flow variables at the boundary are obtained by solving the above equations.

5.3. Outflow

Characteristic variable boundary conditions are applied to outflows as well. For the subsonic outflow, the first four eigenvalues λ(1,2,3,4) are positive and the fifth eigenvalue λ(5) is negative. Flow variables at the boundary are connected with the interior values through the first four characteristic relations. The fifth relation is replaced by the static pressure at the exit

( ρ b p b / β c 0 2 ) n x ( Θ n t ) b ψ 0 n x w b n y + v b n z = E70
( ρ i p i / β c 0 2 ) n x ( Θ n t ) i ψ 0 n x w i n y + v i n z + [ R 0 1 Γ q M 1 S b Δ t ] ( 1 ) E71
( ρ b p b / β c 0 2 ) n y ( Θ n t ) b ψ 0 n y u b n z + w b n x = E72
( ρ i p i / β c 0 2 ) n y ( Θ n t ) i ψ 0 n y u i n z + w i n x + [ R 0 1 Γ q M 1 S b Δ t ] ( 2 ) E73
( ρ b p b / β c 0 2 ) n z ( Θ n t ) b ψ 0 n z v b n x + u b n y = E74
( ρ i p i / β c 0 2 ) n z ( Θ n t ) i ψ 0 n z v i n x + u i n y + [ R 0 1 Γ q M 1 S b Δ t ] ( 3 ) E75
p b 2 ρ 0 σ 0 + ( Θ n t ) b c 0 + = p i 2 ρ 0 σ 0 + ( Θ n t ) i c 0 + + [ R 0 1 Γ q M 1 S b Δ t ] ( 4 ) E76
p b = p e E77

The boundary values (ρ b , p b , u b , v b , w b ) are solved by the above equations.

For a supersonic outflow case, the fifth eigenvalue λ(5) is positive, and no characteristic waves can propagate from downstream to the upstream. The last equation (57) is replaced by the fifth characteristic relation (58) where characteristic variables come from the interior field:

p b 2 ρ 0 σ 0 + ( Θ n t ) b c 0 = p i 2 ρ 0 σ 0 + ( Θ n t ) i c 0 + [ R 0 1 Γ q M 1 S b Δ t ] ( 5 ) E78

5.4 Impermeable wall

For impermeable surfaces, the specified condition is that of no flow past perpendicular to the surface, i.e., Θ b = 0 . The remaining boundary conditions are found from the following characteristic relations:

( ρ b p b / β c 0 2 ) n x w b n y + v b n z = ( ρ i p i / β c 0 2 ) n x Θ i ψ 0 n x w i n y + v i n z + [ R 0 1 Γ q M 1 S b Δ t ] ( 1 ) E79
( ρ b p b / β c 0 2 ) n y u b n z + w b n x = ( ρ i p i / β c 0 2 ) n y Θ i ψ 0 n y u i n z + w i n x + [ R 0 1 Γ q M 1 S b Δ t ] ( 2 ) E80
( ρ b p b / β c 0 2 ) n z v b n z + u b n y = ( ρ i p i / β c 0 2 ) n z Θ i ψ 0 n z v i n z + u i n y + [ R 0 1 Γ q M 1 S b Δ t ] ( 3 ) E81
p b 2 ρ 0 σ 0 = p i 2 ρ 0 σ 0 + Θ i c 0 + + [ R 0 1 Γ q M 1 S b Δ t ] ( 4 ) E82
u b n x + v b n y + w b n z = n t E83

It is noticed that the first four characteristic relations (59) through (62) can be obtained by setting Θ b = 0 in Eqs. (53) through (56). It is also worth noting that in all above boundary conditions, the source term contributions are added to the equations implicitly.


6. Results

To validate the new preconditioned algorithm for low Mach number flows in rotating machinery, application of a hydrodynamic flow about the David Taylor marine propeller P5168 is performed in this study. This is a five-bladed, controllable–pitch propeller (Fig. 3) with a design advance coefficient of J = 1.27. The diameter of the propeller is 0.4027m. The 3-D LDV tip-vortex flow field measurement was conducted by Chesnakas and Jessup (Chesnakas & Jessup, 1998) in the CDNSWC 36-inch water tunnel. Four different advance coefficients were measured with detailed velocity data behind the propeller. The numerical simulations were conducted under the four advance ratios corresponding to the experiment. The body-fitted mixed-element anisotropic unstructured mesh was generated using the Advancing-Front Local Reconnection (AFLR) grid generation technique (Marcum, 2001). The far field boundary of the computational domain extends to 3.5 blade diameters in both upstream and downstream, and 1.5 blade diameters in the radial direction. The unstructured mesh contains 3.825 million nodes and 12.182 million elements. For the viscous gird used in this work, the y+ distribution is less than 1.5 over the solid bodies to provide a good viscous sublayer resolution. In addition to the fine grid region over all the viscous surfaces for achieving good boundary layer resolution, a number of interior faces were created downstream in the domain along blade wake paths to refine the volume grid in order to capture the fine details of the flow field in the downstream of the propeller. Those interior faces are transparent to the flow solver with no physic boundary conditions required. As moving further beyond the shaft rear face, the computational mesh becomes coarse to reduce the computational costs. The spatial variation of grid resolution is expected to affect the accuracy of the numerical simulations and its influence will be assessed in the investigation of the simulation results. Computation domain was partitioned into 64 sub-domains for parallel computations.

The velocity measurements were performed by Chesnakas and Jessup (Chesnakas & Jessup, 1998) in the water tunnel, and is used in this work for code validation. Detailed flow conditions for the simulation cases are summarized in Table 1. The characteristic based far-field boundary condition is applied to the far-field boundary with implicit treatment of the source terms. The no-slip boundary conditions are applied to the propeller blade and shaft

Figure 3.

Marine Propeller 5168 and Surface Mesh

Figure 4.

Convergence History of Propeller Thrust

Advance Coefficient J=V/nD Free Stream V(m/s) Rotational Speed RPM Tip Velocity V tip (m/s) Tip Mach Number M tip Reynolds Number Re (106)
0.98 7.893 1200 25.302 0.0744 3.14
1.1 10.705 1450 30.574 0.0899 4.26
1.27 11.081 1300 27.411 0.0806 4.41
1.51 11.730 1150 24.248 0.0713 4.64

Table 1.

Simulation Flow Conditions of P5168 Propeller

surfaces in the rotating frame. All solid surfaces are assumed to be adiabatic, i.e., no heat fluxes across the solid surfaces. The Reynolds number is based on the free stream speed and the diameter of the propeller. The reference Mach number Μ r = 0.1 is selected for all the cases based on the free stream condition. The preconditioning parameter is estimated by the modified scheme as β = Μ r 2 Ω 2 , in which Μ Ω is set equal to the propeller tip Mach number listed in Table 1 for each case. This seems to provide the best accuracy and convergence. Due to the adoption of the relative rotating frame, flows of the rotating propeller can be considered as steady problems. The local time stepping is used in all the simulations with CFL number of 10. One Newton’s iteration and eight Gauss-Seidel sub-iterations are used at each time step. Simulations are initiated from a uniform flow, and converged steady solutions are achieved in less than 1000 iterations. Fig. 4 shows the simulation history of the computed propeller thrust coefficient to demonstrate the convergence at each advance ratio of the propeller.

6.1. Propeller performance

The propeller has been open-water tested in the CDNSWC towing basin, where the cavitation inception conditions are related to either thrust loading (KT) or tunnel velocity (J). The numerical simulations are performed under the same four advance ratios J=0.98, 1.1, 1.27, 1.51 that have been tested in the water tunnel. The computed thrust and torque on the propeller are calculated by integrating the pressure over the propeller surfaces. The thrust and torque coefficients are defined as

K T = t h r u s t ρ n 2 D 4 ,
K Q = t o r q u e ρ n 2 D 5 s E84

where ρ is the flow density, n is the rate of revolutions in the unit of revolution per second, and D is the propeller diameter. Fig. 5 shows that computed thrust and torque coefficients for the P5168 propeller agree well with experiment in all advance ratios. Results obtained by the current preconditioning method seem to provide excellent predictions of the propeller thrust and torque comparing with the experimental data.

6.2. Circumferentially averaged velocity

The measured and computed circumferentially averaged velocities at two downstream locations (x/R=0.2386 and 0.8378) are shown in Fig. 6. The filled symbols are measured values and unfilled symbols are computed values. Predicted axial and tangential velocities match the experiment quietly well, although a slight underestimated radial velocity is seen in Fig. 6 (c). This is related to the prediction of tip-vortex flow in the downstream of the

Figure 5.

Comparison of Predicted and Measured Propeller Thrust and Torque

propeller, which will be discussed in the following section. However, the overall trend in velocity distributions are well captured in the computation, including in the blade tip-vortex region at r/R=1.0. The computed results showed favorable agreement with the experimental data.

6.3. Mean velocity contours

To further investigate the tip-vortex flow generated by the propeller, the computed contours of the axial, radial, and tangential velocity components in the rotating frame at a stream-wise location of x/R = 0.2386 are presented in Fig. 7. The tip vortices generated by the blade tip are clearly visible in both the computation and the experiment, where the vortex core is represented as blue “spokes” in the plot of low velocity. There is a flow circulation around the tip vortex, as well as the radial velocity generated by the blade wake. The strong asymmetry of the vortex can be seen here. It also indicates that the current preconditioned scheme does capture the general feature of the tip-vortex flow. The axial and radial velocity components (Fig. 7(a) and (b)) showed a good agreement regarding the locations of tip vortices between the computation and measurement. However, the numerical diffusion was observed in the current second-order discretization scheme in the vortex core region, as predicted vortex strength appeared not as strong as what was measured in the experiment. To improve the prediction accuracy in the vortex core region, an adaptively refined mesh or a higher-order spatial discretization scheme (Sheng & Wang, 2009) is necessary. Nevertheless, the modified preconditioning method was proved to provide stable and reasonably accurate predictions to the hydrodynamic viscous flows in rotating machinery, which would otherwise impossible to achieve with the original preconditioning method, not to mention the compressible flow solver.

Figure 6.

Comparison of Computed and Measured circumferentially Averaged Velocities

Figure 7.

Comparison of Computed and Measured Mean Velocity Contours at x/R=0.2386

6.4. Tip-vortex trajectory

The tip vortex trajectory was predicted and compared with the experiment measurements. An approximate approach to locate the tip vortex core location is to limit the search for a minimum pressure within the region of the known vortex indicated by a swirl parameter. The physical interpretation for this criterion is given by Jeong and Hussain (Jeong & Hussain, 1995). Within the vortex core, pressure tends to have a local minimum on the axis of a circulatory or swirling motion when the centrifugal force is balanced by the pressure force, which is strictly true only in a steady inviscid planar flow. Using this method, the tip-vortex trajectories were identified from computed flows at the advance ratio of J=1.1. Fig. 8 compared the predicated radial locations of the tip vortex with the experimental observation. In the experiment, the tip vortex was tracked from the blade tip to an axial location of approximately one diameter downstream where x/R=2.0. The current simulation conducted in this study can predict the tip vortex radial trajectory accurately up to a limited distance downstream, due to the grid coarsening in the downstream of the propeller. As shown in Fig. 8, the predicted tip-vortex trajectory is tracked to the location about only one radius downstream at x/R=1.0 before dissipated into the flow field. The largest disagreement between the prediction and the experiment measurement occurs at the last station of x/R=1.0. Beyond that point it appeared to be difficult to trace the tip vortex as it has completely smeared out in the flow field. Examination of the grid resolution of the propeller indicated a finer grid resolution only coving a limited distance in the downstream of the propeller. The grid cells became very coarser as moving beyond the shaft rear surface at station x/R=0.69. Fig. 9 clearly demonstrated the tip vortex convection which is quickly dissipated due to the limited grid resolution. As discussed in the previous section, a refined computational mesh and a higher-order discretization scheme may reduce the numerical dissipation of this vortical flow, although the computational costs would be significantly higher.

Figure 8.

Propeller 5168 Tip Vortex Radial Trajectory

Figure 9.

Computed Tip Vortex Visualized by Swirl Parameter

6.5. Tip vortex convection and decay

Another way to trace the tip vortex of the propeller is to use a non-dimensional parameter termed intrinsic swirl parameter τ (Berdahl & Thompson, 1993), which calculates the velocity gradient tensor over the computational domain. The intrinsic swirl indicates the tendency for the fluid to swirl about a local point, and is more effective to represent the vortical motion in the flow field. In regions where the swirl parameter approaches to zero ( τ 0 ), the fluid convects too rapidly to be captured in swirling motion, while in region where τ 0 , the fluid does not move quickly enough and is trapped in a swirling motion. Fig. 10 compares the tip vortex strength under the four advance coefficients. The intrinsic swirl contours are visualized on six downstream stations of x/R=[0.1, 0.6] with an equal spacing of 0.1. Streamlines are generated from the shaft rear surface at x/R=0.69. The core regions of blade tip vortices are evident by the largest magnitude of swirl parameter. It clearly shows the rapid decay of the tip vortices and the separation of vortices from the wake as the flow moves downstream. The swirl parameter values at the vortex core on selected stream-wise stations are extracted and presented in Fig. 11. The simulation results indicated that vortex strength varied strongly with advance coefficient. At a high loading (lower advance ratio values, J=0.98 and 1.10), strong tip vortices are observed over a longer distance downstream. The tip vortices are very weak at J=1.27 and eventually disappeared at J=1.51. It is noticed that the decreasing of the intrinsic swirl under the advance coefficients of J=1.27 and 1.51 is not monotonic as x/R goes beyond 0.4, where the swirl parameter even increases at further downstream stations. Further investigation found that the locations of high swirl centers at station x/R=0.5, 0.6 under J=1.27 and 1.51 occur in a much finer grid region, which was originally generated for the purpose of capturing the propeller wake. Therefore, the non-physical behaves of the swirl parameter appear in Fig. 11 are due to the variation of grid resolution. The simulation results are qualitatively consistent with the experimental measurements.

Figure 10.

Predicted Wakes and Tip Vortices vs. Advance Coefficient

Figure 11.

Tip Vortex Decay vs. Advance Coefficient J

6.6. Tip-vortex attachment

The water tunnel experiment also observed suction-side tip vortex attachment occurred aft of the blade tip over certain range of tested advance ratio. The locations of the tip vortex attachment on the five blades showed minor variation from blade to blade, and their average radial location has been estimated. In the open-water experiment at an advance ratio of J=0.98, this average location was estimated at the blade trailing edge of 0.99R radius, while at J=1.1, the attachment point was moved aft to 0.998R radius of the trailing edge. From the current computed results, the iso-surface of swirl parameter τ=1.0 is generated to visualize the tip vortex structure as shown in Fig. 12. The iso-surface is shaded by static

Figure 12.

Blade Suction Side Tip Vortex Attachment.

pressure. The propeller geometry is not plotted in the figure, but the tip region of the propeller blade is seen clearly, which is actually the swirl iso-surface in the blade boundary layer. At the blade tip, the attachment occurs aft a local minimum pressure center at the trailing edge. The iso-surface of the swirl parameter indicated that the tip vortex attachment occurs at the trail edge on the blade suction side between 0.994R-0.999R. The variation of the attachment points at J=0.98 and J=1.1 are not distinguishable except that the local minimum pressure center upstream the attachment location has lower pressure at J=0.98 and originates stronger tip vortex. The current simulations correctly captured the tip-vortex attachment observed in the experiment, which is important for the improved understanding of the cavitation inception of the marine propeller.


7. Conclusion

A modified preconditioning method was investigated and validated in the prediction of hydrodynamic viscous flows for marine propeller P5168. The preconditioning parameter is based on the reference Mach number and rotating Mach number to provide stable and accurate solutions for low Mach number flows in rotating machinery, while the original preconditioning method failed to provide converged solutions in the case presented in this study. The predicted overall propeller performance, circumferentially averaged velocities, mean velocity contours, tip-vortex trajectory and decay at blade downstream stations are validated against the experimental data. The comparison between the computation and the experiment indicates that the current preconditioned solver was able to capture the general features of the tip vortex generated by the propeller. In particular, the predicted thrust and torque coefficients at several advance ratios matched well with the open-water experimental data. However, the current numerical simulation showed a quick decay of the tip vortex in the propeller downstream, due to a lack of enough grid resolution. Future studies include using adaptive grid technique to locally refine the computational mesh in the vortex core region, and a higher-order spatial discretization schemes to further reduce the numerical diffusion in the solver.



The author wishes to thank Qiuying Zhao of The University of Toledo and Xiao Wang of High Performance Computing Collaboratory at Mississippi State University for helping processing the simulation data and formatting the manuscript.


  1. 1. Anderson W. K. 1992 Grid Generation and Flow Solution Method for Euler Equations on Unstructured Grids, NASA Technical Report TM- 4295, April 1992.
  2. 2. Anderson W. K. Rausch R. D. Bonhaus D. L. 1995 Implicit/Multigrid Algorithms for Incompressible Turbulent Flows on Unstructured Grids, AIAA Paper 95 1740 , June 1995.
  3. 3. Berdahl C. H. Thompson D. S. 1993 Eduction of Swirling Structure using the Velocity Gradient Tensor, AIAA J., 31 1 (January 1993), 97-103.
  4. 4. Briley W. R. Mc Donald H. Shamroth S. J. 1983 A Low Mach Number Euler Formulation and Application to Time Iterative LBI Schemes, AIAA J., 21 (10), (1983),1467-1469.
  5. 5. Briley W. R. Tayler L. K. Whitfield D. L. 2003 High-Resolution Viscous Flow Simulations at Arbitrary Mach Number, Journal of Computational Physics, 184(1), (2003), 79-105.
  6. 6. Chesnakas C. J. Jessup S. D. 1998 Cavitation and 3 -D LDV Tip-Flowfield Measurements of Propeller 5168, CRDKNSWC/HD-1460-02, (May 1998), Carderock Division, Naval Surface Warfare Center.
  7. 7. Chen J. P. Ghosh A. R. Sreenivas D. Whitfield D. L. 1997 Comparison of Computations Using Navier-Stokes Equations in Rotating and Fixed Coordinates for Flow Through Turbomachinery, AIAA Paper 97-0878 97 0878 , AIAA 35th Aerospace Sciences Meeting and Exhibit, Reno, NV, January 6-10, 1997.
  8. 8. Choi D. Merkle C. L. 1985 Application of Time-Iterative Schemes to Incompressible Flow, AIAA J., 23 10 1518 1524 .
  9. 9. Choi D. Merkle C. L. 1993 Application of Preconditioning to Viscous Flows, J. Comp. Physics, 105 207 223 .
  10. 10. Chorin A. J. 1967 A Numerical Method for Solving Incompressible Viscous Flow Problems, J. Comp. Physics, 2 12 26 .
  11. 11. Coirier W. J. 1994 An Adaptively-Refined, Cartesian, Cell-Based Scheme for the Euler and Navier-Stokes Equations, NASA Technical Memorandum 106754, NASA Lewis Research Center, October 1994.
  12. 12. Hyams D. G. 2000 An Investigation of Parallel Implicit Solution Algorithms for Incompressible Flows on Unstructured Topologies, Ph.D. Dissertation, Mississippi State University, May 2000.
  13. 13. Hageman L. A. Young D. M. 1981 Applied Iterative Methods, 1981, Academic Press, New York.
  14. 14. Jeong J. Hussain F. 1995 On the Indentification of a Vortex, Journal of Fluid Mechanics, 285 (1995), 69-94.
  15. 15. Marcum D. L. 2001 Efficient Generation of High Quality Unstructured Surface and Volume Grids, Engineering with Computers, 17 3 (2001), 211-233.
  16. 16. O’Brien G. G. Hyman M. A. Kaplan S. 1951 A Study of the Numerical Solution of Partial Differential Equations, Journal of Mathematics and Physics, 29 4 (1951), 223-249.
  17. 17. Sheng C. Wang X. 2006 A Global Preconditioning Method for Low Mach Number Viscous Flows in Rotating Machinery, GT2006-91189, Proceedings of ASME Turbo Expo 2006 8 11 May 2006, Barcelona, Spain.
  18. 18. Sheng C. Wang X. 2003 Characteristic Variable Boundary Conditions for Arbitrary Mach Number Algorithm in Rotating Frame, AIAA-2003 3976 , Proceedings of the 16th AIAA Computational Fluid Dynamics Conference, June 23-26, 2003, Orlando, FL.
  19. 19. Sheng C. Wang X. 2009 Aerodynamic Analysis of a Spinning Missile with Dithering Canards Using a High Order Unstructured Grid Scheme, AIAA-2009 1090 th AIAA Aerospace Sciences Meeting, 5-8 January 2009, Orlando, Florida.
  20. 20. Spalart P. Allmaras S. 1991 A One-Equation Turbulence Model for Aerodynamic Flows, AIAA Paper 92 0439 , January 1991.
  21. 21. Vichnevetsky R. Bowles J. B. 1982 Fourier Analysis of Numerical Approximations of Hyperbolic Equations, 1982.
  22. 22. Wang X. Sheng C. 2005 Numerical Study of Preconditioned Algorithm for Rotational Flows, 17th AIAA Computational Fluid Dynamics Conference, June 6-9, 2005, Toronto, Ontario, Canada
  23. 23. Wang X. 2005 A Preconditioned Algorithm for Turbomachinery Viscous Flow Simulation. Ph.D. Dissertation, Mississippi State University, December 2005.

Written By

Chunhua Sheng

Published: January 1st, 2010