Table of Contents Author Guidelines Submit a Manuscript
Modelling and Simulation in Engineering
VolumeΒ 2011, Article IDΒ 537464, 17 pages
Research Article

A Preconditioned Method for Rotating Flows at Arbitrary Mach Number

Mechanical, Industrial, and Manufacturing Engineering Department, The University of Toledo, 2801 W. Bancroft Street, MS 312, Toledo, OH 43606, USA

Received 17 April 2011; Accepted 25 July 2011

Academic Editor: AntonioΒ Munjiza

Copyright Β© 2011 Chunhua Sheng. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


An improved preconditioning is proposed for viscous flow computations in rotating and nonrotating frames at arbitrary Mach numbers. The key to the current method is the use of both free stream Mach number and rotating Mach number to construct a preconditioning matrix, which is applied to the compressible governing equations written in terms of primitive variables. A Fourier analysis is conducted that reveals the efficacy of the modified preconditioning. Numerical approximations for the convective and diffusive fluxes are detailed based on the preconditioned system of equations. A set of boundary conditions using characteristic variables are described for internal and external flow computations. Numerical validations are performed on four realistic viscous flows in both fixed and rotating frames. The results indicated that the modified preconditioning has a superior performance compared to the original method to predict flows from extremely low to supersonic Mach numbers.

1. Introduction

It is well known that compressible flow equations face difficulties at low Mach numbers due to a large ratio of acoustic and convective time scales. If the total enthalpy is constant, then the steady solutions approach the incompressible constant-density limit as the Mach number approaches zero. Such problems are often solved using incompressible flow equations such as the artificial compressibility approach of Chorin [1]. However, in many cases compressible equations are a preferred choice even though the Mach number is low. One example is the low speed flow with localized transonic region or shock in the helicopter rotor forward flight. Another example is when the thermal effect is important, and the energy equation needs to be coupled into the governing equations. Over the past decades, a number of researchers have considered preconditioning of compressible governing equations to overcome the algorithmic limitation at low Mach numbers. The use of preconditioning was first introduced by Briley et al. [2]. A simple constant, diagonal preconditioning matrix was added to a nondimensional form of the isoenergetic equations. This generally improved convergence for a test case with a reference Mach number π‘€π‘Ÿ=0.05 using the ADI factorization scheme in primitive variables. Choi and Merkle [3, 4] also considered a similar preconditioning based on local velocities to improve the convergence of an approximate factorization scheme applied to the compressible Euler and artificial compressibility equations in conservative variables. Turkel [5] devised a generalized preconditioning matrix for the compressible equations in terms of primitive and isoenergetic variables. He suggested using a preconditioning matrix based on local velocity with a limiter to avoid singular behaviour near a stagnation point. The VLR (van Leer-Lee-Roe) preconditioner proposed by van Leer et al. [6] is a symmetric method, sometimes referred to as optimal preconditioner. Although in theory it produces an optimal reduction in the condition number for all Mach numbers, in practice it is not very robust (especially at stagnation points) mainly because it requires a well-defined flow angle. Weiss-Smith preconditioner [7, 8], which belongs to the Turkel family of preconditioning, uses the squared Mach number instead of Mach number and two free parameters in the preconditioning matrix. A good review on preconditioning was provided in [9].

Recently, a high-resolution characteristic-based numerical flux approximation for the preconditioned system of hyperbolic equations was developed by Briley et al. [10]. The same single-parameter constant matrix similar to the previous work [2] was adopted for the compressible governing equations written in terms of primitive variables. A number of validation cases from extremely low Mach numbers to various wall temperature ratios on a flat plate were reported [10]. This global preconditioning based on the free stream Mach number generally provides well-behaved solutions with impressive convergence and accuracy if no rotating flow effect is involved. However, if there is a large variation of rotating speeds across the computational domain, the previous single-parameter preconditioning may lead to deteriorated convergence or even numerical instability at low Mach numbers [11, 12]. Numerical analysis using the Fourier footprint analysis [13] also revealed unbalanced wave propagation of the original preconditioning applied to rotating flows [11].

The purpose of this paper is to document the recent development of a modified preconditioning towards viscous flow computations in both rotating and nonrotating flow conditions. In rotating flows, the flow filed is complicated by the secondary swirling structures due to high-speed rotating elements. The speeds of characteristic waves of the physical equations are strongly affected by the free stream velocity as well as the rotational speeds of the fluid. It is thus rational to include the effect of rotating speed in constructing the preconditioning matrix for low Mach number flow computations. The compressible governing equations are written in a generalized form suitable for both rotating and nonrotating flows based on the absolute velocity components [14, 15]. This provides a unified solution algorithm to solve steady and unsteady rotating flows at arbitrary Mach numbers.

This paper is organized as followed. In Section 2 the hyperbolic conservation laws for the unsteady Reynolds-averaged Navier-Stokes governing equations are presented in a general frame. The preconditioned system of equations is described in detail in Section 3, with emphases on the selection of a general preconditioning parameter for both rotating and nonrotating flows. A Fourier footprint analysis is performed based on a two-dimensional Euler flow with rotational speeds. Following that in Section 4 are the numerical discretizations for the convective (inviscid) and diffusive (viscous) fluxes and the iterative scheme for solving the nonlinear system of equations. The Spalart-Allmaras turbulence model is included in Section 5 for completeness. In Section 6, the characteristic-based boundary conditions are developed based on the preconditioned equations for both internal and external flows. Section 7 presents numerical validations of four viscous flows, including two rotating flows and two nonrotating flows at Mach numbers ranging from extremely low to supersonic. Finally, some concluding remarks are given in Section 8.

2. Conservation Laws

2.1. Normalization

It is helpful to express the governing equations in a nondimensional form. Here, the following reference quantities are chosen to normalize the solution variables and governing equations, where the subscript π‘Ÿ represents the state of reference values: density πœŒπ‘Ÿ; velocity π‘ˆπ‘Ÿ; temperature π‘‡π‘Ÿ; pressure πœŒπ‘Ÿπ‘ˆ2π‘Ÿ; length πΏπ‘Ÿ; time πΏπ‘Ÿ/π‘ˆπ‘Ÿ; energy and enthalpy πΆπ‘Ÿπ‘‡π‘Ÿ. Details of the nondimensional equation of state for ideal gas are given in the appendix.

2.2. Governing Equations

In the present study, the integral form of three-dimensional Navier-Stokes governing equations is expressed in a general rotating reference frame. It represents a system of conservation laws for a control volume that relates the rate of change of a vector of average state variables to the flux through the control volume. If the rotating frame has a rotational speed of Ξ©, the integral form of the compressible governing equations can be written asπœ•ξ€œπœ•π‘‘π‘‰βƒ—ξ€Ÿπ‘„π‘‘π‘‰+π‘†ξ‚€π‘ƒξ‚€βƒ—π‘„ξ‚ξ‚ξ€œβ‹…βƒ—π‘›π‘‘π΄=𝑉⃗𝑆⃗𝑄𝑑𝑉,(1) where ⃗𝑄=(𝜌,πœŒπ‘’,πœŒπ‘£,πœŒπ‘€,𝑒)𝑇 is a vector of conservative flow variables based on the absolute velocity components, βƒ—π‘ˆ=(𝑒,𝑣,𝑀)𝑇. ⃗𝑛=(𝑛π‘₯,𝑛𝑦,𝑛𝑧)𝑇 is the normal vector on the face (𝑑𝐴) of a control volume (𝑑𝑉). In the surface integral term, 𝑃 is a tensor consisting of convective and diffusive fluxes. The complete descriptions about the convective and diffusive fluxes are provided in the appendix.

It is noted that the above system of governing equations contains a source term ⃗𝑆, which is the Coriolis and centrifugal forces introduced into the momentum equations by the rotating coordinate system; see the appendix. Since the above governing equations (1) are expressed in terms of absolute velocity components, it is very easy to transform them from the rotating frame to the fixed absolute frame by simply setting the source term to zero. However, a key difference between the governing equations in the rotating and fixed frames is that the grid velocity ξ‹π‘ˆπ‘ =(𝑒𝑠,𝑣𝑠,𝑀𝑠)𝑇 is not included in the former but is included in the latter. More details can be found in [15].

In order to develop a preconditioned numerical algorithm at all Mach numbers, the conservative-variable governing equations are modified and written in terms of primitive variables βƒ—π‘ž=(𝜌,𝑒,𝑣,𝑀,𝑝)𝑇. Denoting that ⃗𝑀=πœ•π‘„/πœ•βƒ—π‘ž is a transformation matrix from the conservative to primitive variables, the governing equation (1) can be written asπ‘€πœ•ξ€œπœ•π‘‘π‘‰ξ€Ÿβƒ—π‘žπ‘‘π‘‰+π‘†ξ‚€π‘ƒξ€·ξ€Έξ‚ξ€œβƒ—π‘žβ‹…βƒ—π‘›π‘‘π΄=π‘‰βƒ—π‘†ξ€·ξ€Έβƒ—π‘žπ‘‘π‘‰.(2)

3. Preconditioned Method

3.1. Preconditioned System of Equations

To simplify the expression of the preconditioned method, consider only the convective part of the three-dimensional compressible governing equations in a differential form:πœ•βƒ—π‘„+πœ•βƒ—πΉξ‚€βƒ—π‘„ξ‚πœ•π‘‘+πœ•βƒ—πΊξ‚€βƒ—π‘„ξ‚πœ•π‘₯+πœ•ξ‹π»ξ‚€βƒ—π‘„ξ‚πœ•π‘¦=⃗𝑆⃗𝑄.πœ•π‘§(3) Equation (3) can be written in terms of primitive variables asπ‘€πœ•βƒ—π‘ž+πœ•π‘‘π΄π‘€πœ•βƒ—π‘ž+πœ•π‘₯π΅π‘€πœ•βƒ—π‘ž+πœ•π‘¦πΆπ‘€πœ•βƒ—π‘ž=βƒ—π‘†ξ€·ξ€Έπœ•π‘§βƒ—π‘ž(4) orπœ•βƒ—π‘ž+πœ•π‘‘π‘Žπœ•βƒ—π‘ž+πœ•π‘₯π‘πœ•βƒ—π‘ž+πœ•π‘¦π‘πœ•βƒ—π‘ž=πœ•π‘§π‘€βˆ’1⃗𝑆,βƒ—π‘ž(5) where π‘Ž=π‘€βˆ’1𝐴𝑀, 𝑏=π‘€βˆ’1𝐡𝑀, and 𝑐=π‘€βˆ’1𝐢𝑀 are the system matrices. Matrices ⃗⃗𝑄𝐴=πœ•πΉ/πœ•, ⃗⃗𝑄𝐡=πœ•πΊ/πœ•, and ⃗𝑄𝐢=πœ•π»/πœ• are the Jacobian matrices of the convective flux vectors (⃗⃗𝐻𝐹,𝐺,and)with respect to conservative variables ⃗𝑄. Introducing a preconditioning matrix Ξ“π‘žβˆ’1 into the time derivative term of (5) yieldsΞ“π‘žβˆ’1πœ•βƒ—π‘ž+πœ•π‘‘π‘Žπœ•βƒ—π‘ž+πœ•π‘₯π‘πœ•βƒ—π‘ž+πœ•π‘¦π‘πœ•βƒ—π‘ž=πœ•π‘§π‘€βˆ’1βƒ—π‘†ξ€·ξ€Έβƒ—π‘ž(6) orπœ•βƒ—π‘žπœ•π‘‘+Ξ“π‘žπ‘Žπœ•βƒ—π‘žπœ•π‘₯+Ξ“π‘žπ‘πœ•βƒ—π‘žπœ•π‘¦+Ξ“π‘žπ‘πœ•βƒ—π‘žπœ•π‘§=Ξ“π‘žπ‘€βˆ’1βƒ—π‘†ξ€·ξ€Έβƒ—π‘ž(7)

Equation (7) is the preconditioned system of governing equations cast in the rotating frame using primitive variables. It is important to notice that (7) is only used to derive the eigenvalues and eigenvectors of the system matrix for calculating the numerical fluxes and boundary conditions. It is not intended for replacing the original governing equations. In other words, the preconditioning technique only modifies the numerical aspects of the flux approximation (such as the convergence and accuracy), not the physics of the governing equations.

3.2. Eigensystem

Comparing to (5), it is noted that system matrices in (7) have been modified by a preconditioning matrix Ξ“π‘ž, which has the following form:Ξ“π‘žξ‚ƒξ‚„=diag1111𝛽,(8) where 𝛽 is the preconditioning parameter which will be discussed in the following section. The preconditioning matrix modifies the eigenvalues and eigenvectors of the original system of equations and thus changes the convergence and dissipation property of characteristics-based numerical fluxes, which are required for the surface integration in finite volume schemes. In the context of a finite volume method, the system matrix of the governing equations (7) can be written in the normal direction ⃗𝑛=(𝑛π‘₯,𝑛𝑦,𝑛𝑧)𝑇on the face of a control volume asπ‘˜Ξ“=Ξ“π‘žξ‚€π‘Žπ‘›π‘₯+𝑏𝑛𝑦+𝑐𝑛𝑧=π‘ŽΞ“π‘›π‘₯+𝑏Γ𝑛𝑦+𝑐Γ𝑛𝑧.(9) The system matrix π‘˜Ξ“ has five real eigenvalues:πœ†(1)=πœ†(2)=πœ†(3)=Θ,πœ†(4)=Ξ˜π›½+πœ†+𝜎,(5)=Ξ˜π›½+βˆ’πœŽ,(10) where Θ=𝑒𝑛π‘₯+𝑣𝑛𝑦+𝑀𝑛𝑧+𝑛𝑑, 𝑛𝑑=βˆ’(𝑒𝑠𝑛π‘₯+𝑣𝑠𝑛𝑦+𝑀𝑠𝑛𝑧) is the grid speed in the normal direction of the face, √𝜎=(Ξ˜π›½βˆ’)2+𝛽𝑐2, and 𝛽±=(1±𝛽)/2. If the preconditioning parameter is set to one, no preconditioning is applied to the governing equations, and the original eigensystem of the governing equations is thus restored. Letting 𝑅 be the matrix composed of the right eigenvectors, the Jacobian matrix, π‘˜Ξ“, can be diagonalized asπ‘…βˆ’1π‘˜Ξ“π‘…=ΛΓ,(11) where ΛΓ is the diagonal matrix containing the real eigenvalues.

3.3. Preconditioning Parameter

As seen in (10), the eigensystem of the three-dimensional Euler equations consists of both convective and acoustic waves. At low Mach numbers, the speed of sound is very large (𝑐2=𝑇/𝑀2π‘Ÿ), causing large disparity between the convective and acoustic wave speeds. For explicit schemes the time step is restricted by the fastest moving wave. If the wave speeds differ significantly, then the slower waves will propagate very slowly and the convergence will also be slow. It has been demonstrated by Briley et al. [2] that this disparity in wave speeds had a strong influence on the approximate factorization errors of the ADI scheme and thus the convergence behaviour at low Mach number in multidimensional problems.

Wang et al. [11] and Sheng and Wang [12] recently conducted a Fourier analysis to investigate the characteristics of various error modes of an implicit Euler scheme analogous to (7). The Fourier footprint reveals the error propagating and damping characters of all frequency modes produced by a numerical discretization. The imaginary axis of the Fourier footprints corresponds to the wave propagation speed, and the real axis (negative) of the Fourier footprints corresponds to the dissipation component of a numerical solution. Due to the fact that an analytic formula of eigenvalues over a fourth order is not possible, the analytic expression of the three-dimensional Fourier footprints is not guaranteed. To this end, the Fourier analysis was performed on a two-dimensional Euler equation based on a uniform grid with a fixed aspect ratio of one [11, 12].

Figure 1(a) shows the Fourier footprints of the original Euler operator (no preconditioning) for a low Mach number flow at π‘€π‘Ÿ=0.01 in the π‘₯ direction, where no rotating velocity component was included. In this low Mach number flow, there exists a large disparity in the wave speeds between the convective (red and black colours) and acoustic (blue and green colours) terms because 𝛽𝑐2=𝛽𝑇/𝑀2π‘Ÿ is overwhelmingly large. In addition, there exists a large numerical dissipation associated with the discretized solution, as indicated by the large magnitude of the real component in the Fourier footprints. To correct the imbalance of the wave speeds and to reduce the condition number of the system matrix, Briley et al. [2] introduced the following constant preconditioning to rescale the acoustic term:𝛽=min1,𝑀2π‘Ÿξ€Έ(12) which scales 𝛽𝑐2 to 𝑇=𝑂(1) and provides well-balanced wave speeds for low Mach number flows as shown in Figure 1(b). The dissipation characteristics associated with the real component of the complex eigenvalues are also improved significantly among all wave modes.

Figure 1: Fourier footprint of Euler operator for nonrotating flow (π‘€π‘Ÿ=0.01).

However, the above single-parameter preconditioning seemed to be effective only in nonrotating flows and became problematic in flows involving a large rotating speed in the computational domain. Figure 2 shows the previous case with the introduction of a rotating speed Ξ© along the π‘₯ direction. The rotating Mach number was 0.1, evaluated based on the reference length scale and rotational speed Ξ©. The Fourier footprints using the original preconditioning parameter indicated more footprints collapsed to the origin for the low wave frequency modes, as shown in Figure 2(a). Neither damping nor propagation of the error waves was effective, and the discretized numerical scheme could not converge to the accurate physical solution.

Figure 2: Fourier footprint of Euler operator for rotating flow (π‘€π‘Ÿ=0.01, 𝑀Ω=0.1).

It is therefore desired to scale the above acoustic term 𝛽𝑐2 to Θ instead of 𝑂(1) as the rotating speed component in Θ may be significantly larger than the free stream speed. A modified preconditioning is then proposed that combines both reference Mach number π‘€π‘Ÿ and rotating Mach number 𝑀Ω into the preconditioning parameter:𝛽=min1,𝑀2π‘Ÿ+𝑀2Ξ©ξ€Έ,(13) where π‘€π‘Ÿ is the global reference Mach number based on the free stream velocity. 𝑀Ω is the rotating Mach number based on the characteristic rotating speed and reference length scale. An alternative way to calculate 𝑀Ω is to use the local rotating speed so that the modified preconditioning is depended on both global free stream Mach number and local rotating speed. In both choices the method restores to the original preconditioning if there is no rotating speed effect involved (Ξ©=0). As conducted by Wang et al. [11] and Shang and Wang [12], Figure 2(b) shows the Fourier footprints for the above rotating flow using the modified preconditioning formula (13). A significant improvement was achieved on all wave modes compared with the original preconditioning method. No Fourier footprints were found collapsed near the origin, and both damping and propagating properties of the error waves were significantly improved.

A validation was performed in the present study to demonstrate the efficacy of the modified preconditioning for a three-dimensional low Mach number viscous flow about a rotating body at a speed Ξ©=2. The test geometry is a simple body of revolution with a free stream velocity at an extremely low Mach number π‘€π‘Ÿ=0.001. The characteristic rotating Mach number is 𝑀Ω=0.08, evaluated based on the reference length scale, free stream velocity, and the speed of sound at a reference condition. Figure 3 shows the convergence of solutions under various preconditioning options: (1) no preconditioning, (2) the original preconditioning, (3) the modified preconditioning with a local rotating Mach number, and (4) the modified preconditioning with a constant global rotating Mach number. The convergence criterion was based on the standard 𝐿2 norm of the unsteady residuals composing of all solution variables on all grid nodes. It is interesting to see that the convergence of the original preconditioning method was even worse than that without any preconditioning. The modified preconditioning using both local and global rotating Mach numbers significantly improved the convergence of the rotating flow at this extremely low Mach number. However, the constant preconditioning based on the characteristic rotating Mach number seemed to perform better than that based on the local rotating Mach number.

Figure 3: Convergence histories of different preconditioning options on a rotating flow.

4. Discretization Method

The governing equations are discretized based on a finite volume formulation. The flow variables are stored at the vertices, and surface integrals are evaluated on the median dual surrounding each of these vertices. The nonoverlapping control volumes formed by the median dual completely cover the domain and form a mesh that is dual to the element grid. Figure 4 shows a diagram of a control volume formed by the median dual faces. The index 𝑖 denotes the vertex of the control volume, and the index𝑗 denotes the edge (or face) number associated with each median dual face.

Figure 4: Control volumes are defined as median duals surrounding each vertex.

The system of governing equations is numerically integrated over the closed boundaries (median dual faces) of a control volume surrounding each node:π‘€Ξ“π‘žβˆ’1πœ•ξ€œπœ•π‘‘π‘‰ξ€Ÿβƒ—π‘žπ‘‘π‘‰+𝑆𝑃Γ⋅⃗𝑛𝑑𝐴+𝑃𝑣=ξ€œβ‹…βƒ—π‘›π‘‘π΄π‘‰βƒ—π‘†π‘‘π‘‰,(14) where 𝑃Γ and 𝑃𝑣 are the tensors of convective and diffusive numerical fluxes evaluated on the face of the control volume, respectively. It is noted that the preconditioning matrix Ξ“π‘žβˆ’1 is only applied to the convective (inviscid) flux term to correct the numerical property of the flux approximation scheme. The diffusive (viscous) term of the governing equations remains unchanged. For unsteady time-accurate flows the preconditioning matrix was not applied to the time derivative term of the governing equations for the finite volume integration, which would otherwise modify the physics of the conservation laws. For this reason, the preconditioning method can be applied to both steady and unsteady time-accurate solutions. The discretized form of the Navier-Stoles governing equations over a control volume can be written asπ‘€Ξ“π‘žβˆ’1Ξ”βƒ—π‘žπ‘–Ξ”π‘‘Ξ”π‘‰π‘–+𝑛𝑗𝑗=1𝑃Γ⋅⃗𝑛𝑗Δ𝐴𝑗+𝑃𝑣⋅⃗𝑛𝑗Δ𝐴𝑗=𝑆𝑖Δ𝑉𝑖,(15) where the subscript 𝑖 denotes the vertex and 𝑗=(1,…,𝑛𝑗) denotes the 𝑗th surface of the control volume surrounding the vertex 𝑖.

4.1. Convective Flux Evaluation

The convective fluxes of the preconditioned governing equations are approximated on the face of a control volume using a formula similar to Roe’s flux scheme, as proposed by Briley et al. [10]. Consider the surface 𝑗 of a control volume whose left and right states are denoted by 𝐿 and 𝑅; the numerical flux projected onto the surface can be written as𝑃Γ1⋅⃗𝑛=2ξ‚€βƒ—πΉξ€·ξ‹π‘žπΏξ€Έ+βƒ—πΉξ€·ξ‹π‘žπ‘…ξ€Έξ‚βˆ’12π‘€Ξ“π‘žβˆ’1|||π‘˜Ξ“|||ξ€·ξ‹π‘žπ‘…βˆ’ξ‹π‘žπΏξ€Έ,(16) where ⃗𝐹(ξ‹π‘žπΏ) and ⃗𝐹(ξ‹π‘žπ‘…) are the physical fluxes evaluated at the left and right states of the surface. The last term on the RHS of (16) is associated with the numerical dissipation for the convective flux approximation, where |π‘˜Ξ“| is the preconditioned system matrix with positive eigenvalues. The following averaged primitive variables Μ‚π‘ž are used on the surface of a control volume:1Μ‚π‘ž=2ξ€·ξ‹π‘žπΏ+ξ‹π‘žπ‘…ξ€Έ.(17) Since the eigensystem and the constructed numerical fluxes are evaluated using the above averaged variables instead of Roe’s variables, the current approach is considered to be a modification of Roe’s flux approximation.

In Equation (18) quantities ξ‹π‘žπΏ and ξ‹π‘žπ‘… are the values of the primitive variables on the left and right sides of the face of a control volume. For a first-order accurate differencing, quantities ξ‹π‘žπΏ and ξ‹π‘žπ‘… are set equal to the value at the vertices lying on either side of the face. For a second-order accurate scheme, these face values (𝑗) are computed with a Taylor series expansion about the central node (i) of the control volumeβƒ—π‘žπ‘—=βƒ—π‘žπ‘–+βˆ‡βƒ—π‘žπ‘–β‹…βƒ—π‘Ÿ,(18) where βƒ—π‘Ÿβ€‰β€‰is the vector that extends from the central node to the midpoint of each edge and βˆ‡βƒ—π‘ž is the gradient of the primitive variables at the vertex and is evaluated with an unweighted least-squares procedure [16–18].

4.2. Diffusive Flux Evaluation

A general diffusive term of the Navier-Stokes equations can be written as 𝑃𝑣𝐹⋅⃗𝑛=βˆ’π‘£,𝐺𝑣,𝐻𝑣𝑇⋅⃗𝑛,(19) where (𝐹𝑣,𝐺𝑣,𝐻𝑣)𝑇 are the viscous fluxes as described in the appendix. The essence to the diffusive flux evaluation is to compute the gradients on the median dual faces of the control volume. Consider an incompressible viscous flow with a constant laminar viscosity; the viscous diffusive term in the momentum equations is identically a Laplace operator acting upon the velocity field. One important property of the Laplace equation is to satisfy the maximum principle, which requires that all stencil weights of a discrete scheme be positive. Positivity is a key property for numerical stability and should be used to guide the viscous flux discretization.

The current finite volume discretization is formulated based on unstructured meshes with mixed elements consisting of tetrahedron, prism, pyramid, and hexahedron, where an edge-based data structured was adopted. If the mesh is composed of tetrahedron elements only, a Galerkin finite element approach [17] may be adopted to accurately evaluate the viscous fluxes. For the mixed element meshes, the calculation of the diffusive fluxes needs to be compatible with the edge-based data structure, where solutions and their gradients are stored at the two vertices of an edge. One way to evaluate the gradients on the face of a control volume is called directional derivative method [19]. The second approach is called a normal directive method, which imposes positivity along the normal direction of a median dual face. As shown in Figure 4, two unit tangential vectors βƒ—π‘™βƒ—π‘š, are introduced that are orthogonal to the unit normal vector ⃗𝑛 on the surface of a control volume:ξ‚΅π‘›βƒ—π‘š=π‘¦βˆ’π‘›π‘§π‘›π‘š,π‘›π‘§βˆ’π‘›π‘₯π‘›π‘š,𝑛π‘₯βˆ’π‘›π‘¦π‘›π‘šξ‚Άπ‘‡,⃗𝑛𝑙=π‘¦π‘›π‘ βˆ’1π‘›π‘š,π‘›π‘¦π‘›π‘ βˆ’1π‘›π‘š,π‘›π‘§π‘›π‘ βˆ’1π‘›π‘šξ‚Άπ‘‡,(20) where π‘›π‘š=ξ”ξ€·π‘›π‘¦βˆ’π‘›π‘§ξ€Έ2+ξ€·π‘›π‘§βˆ’π‘›π‘₯ξ€Έ2+𝑛π‘₯βˆ’π‘›π‘¦ξ€Έ2,𝑛𝑠=𝑛π‘₯+𝑛𝑦+𝑛𝑧.(21) It is easy to verify that the above unit vectors βƒ—π‘š, ⃗𝑙, and ⃗𝑛 are orthogonal to each other. The components of any gradients can be expressed in terms of derivatives in βƒ—π‘š, ⃗𝑙, and ⃗𝑛 directions as:πœ•πœ™=πœ•π‘₯πœ•πœ™π‘›πœ•π‘›π‘₯+πœ•πœ™π‘šπœ•π‘šπ‘₯+πœ•πœ™π‘™πœ•π‘™π‘₯,πœ•πœ™=πœ•π‘¦πœ•πœ™π‘›πœ•π‘›π‘¦+πœ•πœ™π‘šπœ•π‘šπ‘¦+πœ•πœ™π‘™πœ•π‘™π‘¦,πœ•πœ™=πœ•π‘§πœ•πœ™π‘›πœ•π‘›π‘§+πœ•πœ™π‘šπœ•π‘šπ‘§+πœ•πœ™π‘™πœ•π‘™π‘§.(22) The derivative in the normal direction ⃗𝑛 can be approximated byπœ•πœ™β‰ˆπœ™πœ•π‘›2βˆ’πœ™1𝑑𝑠,(23) where πœ™1 and πœ™2 are the values of interest at two vertices of an edge and 𝑑𝑠 is the length of an edge. The above formula guarantees positivity of the operator and results in a second-order accuracy 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 βƒ—π‘š, ⃗𝑙 directions, respectively,πœ•πœ™β‰ˆ1πœ•π‘š2ξ€·βˆ‡πœ™1+βˆ‡πœ™2ξ€Έβ‹…βƒ—π‘š,πœ•πœ™β‰ˆ1πœ•π‘™2ξ€·βˆ‡πœ™1+βˆ‡πœ™2⋅⃗𝑙,(24) where βˆ‡πœ™1 and βˆ‡πœ™2 are the gradients evaluated at two nodes of the edge using a weighted least-squares procedure [18]. To ensure positivity of the tangential derivatives, a limiter is incorporated into the gradients in (22) to ensure the positivity of the tangential derivatives.

4.3. Temporal Discretization

The temporal discretization is considered with an implicit Euler method for the system of governing equations. The use of an implicit scheme allows a larger time step to be used in the computation, which is extremely benficeial in unsteady time-accurate simulations, and results in much faster convergence than any explicit schemes. A general expression is available for the temporal discretization, which is given asMΞ“π‘žβˆ’1Ξ”ξ‹π‘žπ‘›π‘–βˆ’(πœƒ/(1+πœƒ))Ξ”ξ‚”π‘žπ‘–π‘›βˆ’1+1Ξ”π‘‘Ξ”π‘‰π‘–ξ“π‘—βˆˆπ‘(0)𝑃⋅⃗𝑛𝑗𝑛+1Δ𝐴𝑗=𝑆𝑖𝑛+1,(25) where Ξ”ξ‹π‘žπ‘›π‘–=ξ‚”π‘žπ‘–π‘›+1βˆ’ξ‹π‘žπ‘›π‘–, Ξ”ξ‚”π‘žπ‘–π‘›βˆ’1=ξ‹π‘žπ‘›π‘–βˆ’ξ‚”π‘žπ‘–π‘›βˆ’1. Δ𝑑 is the time step increment between steps 𝑛 and 𝑛+1, and Δ𝑉 is the control volume. ⃗𝑃 denotes both inviscid and viscous fluxes. Subscript 𝑖 and 𝑗 denote the indices of the central node and the dual faces that compose of the control volume. A first-order temporal accuracy of the Euler implicit scheme is given by the choice πœƒ=0. Correspondingly, a second-order time-accurate Euler implicit scheme is given by πœƒ=1.

4.4. Newton’s Method

The above nonlinear system of (25) 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 larger time step to be used for numerical computations. Denoteπ‘ξ€·ξ‚”π‘žπ‘–π‘›+1ξ€Έ=MΞ“π‘žβˆ’1Ξ”ξ‹π‘žπ‘›π‘–βˆ’(πœƒ/(1+πœƒ))Ξ”ξ‚”π‘žπ‘–π‘›βˆ’1+1Ξ”π‘‘Ξ”π‘‰π‘–ξ“π‘—βˆˆπ‘(0)𝑃⋅⃗𝑛𝑗𝑛+1Ξ”π΄π‘—βˆ’ξƒ‰π‘†π‘–π‘›+1=0.(26) Newton’s method solves the following equation:π‘ξ…žξ€·ξƒ½π‘žπ‘–π‘›+1,π‘šξ€Έξ€·ξ„¨π‘žπ‘–π‘›+1,π‘š+1βˆ’ξƒ½π‘žπ‘–π‘›+1,π‘šξ€Έξ€·=βˆ’π‘ξƒ½π‘žπ‘–π‘›+1,π‘šξ€Έ,(27) where π‘š = 1,2,3… are the Newton steps, with ξƒ½π‘žπ‘–π‘›+1,0=ξ‹π‘žπ‘›π‘–. π‘ξ…ž(ξƒ½π‘žπ‘–π‘›+1,π‘š) is the Jacobian matrix of (26), that is,π‘ξ…žξ€·ξ‚”π‘žπ‘–π‘›+1ξ€Έ=π‘€Ξ“π‘žβˆ’1𝐼+1Ξ”π‘‘Ξ”π‘‰π‘–ξ“π‘—βˆˆπ‘(0)πœ•ξ‚€ξ‚π‘ƒβ‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπ‘–π‘›+1Ξ”π΄π‘—βˆ’πœ•ξƒ‰π‘†π‘–π‘›+1πœ•ξ‚”π‘žπ‘–π‘›+1.(28)

Equation (28), the first term on the right-hand side is the contribution from the unsteady time derivative of βƒ—π‘ž, 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 π‘€Ξ“π‘žβˆ’1|π‘˜Ξ“| in (16) as a local constant; the approximate flux Jacobian for the inviscid flux can be written as ξ“π‘—βˆˆπ‘(0)πœ•ξ‚€π‘ƒΞ“ξ‚β‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπ‘›+1Δ𝐴𝑗=ξ“π‘—βˆˆπ‘(0)πœ•ξ‚€π‘ƒΞ“ξ‚β‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπΏπ‘›+1Δ𝐴𝑗+ξ“π‘—βˆˆπ‘(0)πœ•ξ‚€π‘ƒΞ“ξ‚β‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπ‘…π‘›+1Δ𝐴𝑗,πœ•ξ‚€π‘ƒΞ“ξ‚β‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπΏπ‘›+1=12βŽ›βŽœβŽœβŽœβŽπœ•ξ‚€βƒ—πΉξ€·ξ‹π‘žπΏξ€Έξ‚π‘›+1πœ•ξ‚”π‘žπΏπ‘›+1+π‘€Ξ“π‘žβˆ’1|||π‘˜Ξ“|||⎞⎟⎟⎟⎠,πœ•ξ‚€π‘ƒΞ“ξ‚β‹…βƒ—π‘›π‘—π‘›+1πœ•ξ‚”π‘žπ‘…π‘›+1=12βŽ›βŽœβŽœβŽœβŽπœ•ξ‚€βƒ—πΉξ€·ξ‹π‘žπ‘…ξ€Έξ‚π‘›+1πœ•ξ‚”π‘žπ‘…π‘›+1βˆ’π‘€Ξ“π‘žβˆ’1|||π‘˜Ξ“|||⎞⎟⎟⎟⎠.(29)

The viscous flux Jacobians are evaluated by taking the derivatives of the discretized viscous flux formulation. Since a constant gradient is assumed for linear distribution of solution variables, the viscous flux Jacobians contain only the local grid information that is constant if no grid motion is involved.

4.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 Ξ”ξƒ½π‘žπ‘›+1,π‘š is obtained through a sequence of iterations, {Ξ”ξƒ½π‘žπ‘›+1,π‘š}𝑖, which converge to Ξ”ξƒ½π‘žπ‘›+1,π‘š. There are several variations of classic relaxation procedures that have been used in the past for solving this linear system of equations [18–20]. In this study, a symmetric implicit Gauss-Seidel procedure [19] 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:[𝐴]=[],𝐷+π‘ˆ+𝐿(30) where[𝐷]=π‘€Ξ“π‘žβˆ’1𝐼+1Ξ”π‘‘ξ“Ξ”π‘‰π‘—βˆˆπ‘(0)πœ•π‘ƒπ‘—π‘›+1πœ•ξ‚”π‘žπ‘›+1Ξ”π΄π‘—βˆ’πœ•βƒ—π‘†πœ•ξ‚”π‘žπ‘›+1,[π‘ˆ]=1ξ“Ξ”π‘‰π‘—βˆˆπ‘π‘ˆ(0)πœ•π‘ƒπ‘—π‘›+1πœ•ξ‚”π‘žπ‘›+1Δ𝐴𝑗,[𝐿]=1ξ“Ξ”π‘‰π‘—βˆˆπ‘πΏ(0)πœ•π‘ƒπ‘—π‘›+1πœ•ξ‚”π‘žπ‘›+1Δ𝐴𝑗.(31)

Letting {𝑅𝑛+1} be the vector of unsteady residuals and {Ξ”ξ‚”π‘žπ‘›+1}the change in the dependent variables, the symmetric Gauss-Seidel relaxation can be written as the following two-step process:[]𝐿+π·Ξ”ξƒ½π‘žπ‘›+1,π‘šξ€Ύπ‘–+1/2+[π‘ˆ]ξ€½Ξ”ξƒ½π‘žπ‘›+1,π‘šξ€Ύπ‘–=𝑅𝑛+1,π‘šξ‚‡,[]𝐷+π‘ˆΞ”ξƒ½π‘žπ‘›+1,π‘šξ€Ύπ‘–+1+[𝐿]ξ€½Ξ”ξƒ½π‘žπ‘›+1,π‘šξ€Ύπ‘–+1/2=𝑅𝑛+1,π‘šξ‚‡.(32)

In the forward pass, {Ξ”βƒ—π‘ž}𝑖+1/2 are obtained with the previously updated {Ξ”βƒ—π‘ž}𝑖, which were set to zero at the initial stage. In the backward pass, {Ξ”βƒ—π‘ž}𝑖+1 are obtained with the most recent value of {Ξ”βƒ—π‘ž}𝑖+1/2 from the previous pass. Normally eight symmetric Gauss-Seidel subiterations are adequate at each time step.

5. Turbulence Model

Several turbulence closure models have been developed for high-Reynolds-number flow computations, such as the one-equation Spalart-Allmaras turbulence model [21, 22] and two-equation π‘˜-πœ€, π‘˜-πœ”, and π‘ž-πœ” turbulence models [23–25]. Spalart-Allmaras turbulence model is a popular one, which is used in the present study. In addition to the original Spalart-Allmaras formulation [21], a variant of this model was also investigated to reduce the smearing of swirl flows [22]. A transport equation for the turbulence Reynolds number ̃𝑣=πœŒπœ‡π‘‘/𝑓𝑣2 (working variable) is expressed asπœ•Μƒπœˆ+βƒ—πœ•π‘‘π‘ˆβ‹…βˆ‡Μƒπœˆ=𝑐𝑏1ξ€·π‘“π‘Ÿ1βˆ’π‘“π‘‘21π‘†Μƒπœˆβˆ’ξ‚Έπ‘Re𝑀1π‘“π‘€βˆ’π‘π‘1π‘˜2𝑓𝑑2ξ‚Ήξ‚ƒΜƒπœˆπ‘‘ξ‚„2+1ξ€½ξ€·πœŽReβˆ‡β‹…ξ€Ίξ€·Μƒπœˆ+1+𝑐𝑏2ξ€Έξ€Έξ€»Μƒπœˆβˆ‡Μƒπœˆβˆ’π‘π‘2ξ€ΎΜƒπœˆβˆ‡β‹…(βˆ‡Μƒπœˆ)=0,(33) where the first and second terms on the right-hand side of (33) are the source production and the third term is the turbulence dissipation. The coefficients in (33) are standard and can be found in [21].

6. Characteristic-Based Boundary Conditions

To solve internal and external rotating flows at arbitrary Mach numbers, four commonly used characteristic-based boundary conditions were developed by Sheng and Wang [26], including the far-field boundary, inflow boundary based on total conditions or mass flow, outflow boundary based on back pressure, and impermeable wall. The characteristic variables were derived based on preconditioned system matrix for consistency 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 briefly outlined for the completeness.

The preconditioned system matrix can be diagonalized as π‘˜Ξ“=π‘…Ξ›Ξ“π‘…βˆ’1,(34) where ΛΓ is a diagonal matrix made up of the eigenvalues of π‘˜Ξ“ and 𝑅and π‘…βˆ’1 are matrices whose columns and rows are the right and left eigenvectors of π‘˜Ξ“, respectively. Multiplying both sides of (7) by 𝑅0βˆ’1 gives the following characteristic relations in one-dimensional spaceπœ•ξ‹π‘Š+πœ•π‘‘Ξ›Ξ“πœ•ξ‹π‘Š=πœ•π‘₯π‘…βˆ’1Ξ“π‘žπ‘€βˆ’1⃗𝑆,(35) where ξ‹π‘Š=𝑅0βˆ’1βƒ—π‘ž are the characteristic variables with the subscript 0 denoting a reference value.

6.1. Far Field

In the far-field boundary, the direction of wave propagation is determined by the sign of each associated eigenvalue. Since all boundary faces have a 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, while a positive sign of the eigenvalue means that the characteristic wave propagates from the interior points of the computational domain to the boundary. It is noted that the following relation is valid along the characteristic line:𝑑π‘₯𝑑𝑑=πœ†Ξ“.(36) For governing equations containing a source term, one hasπ‘‘ξ‹π‘Š=𝑑𝑑𝑅0βˆ’1Ξ“π‘žπ‘€βˆ’1⃗𝑆(37) orξ‹π‘Šπ‘=ξ‹π‘Šπ‘Ÿ+𝑅0βˆ’1Ξ“π‘žπ‘€βˆ’1𝑆𝑏Δ𝑑,(38) where subscript 𝑏 denotes a point on the boundary and subscript π‘Ÿ denotes either an interior point (𝑖) or an outer point (π‘œ), depending on the sign of the characteristic wave propagation on the boundary. The primitive variables at the far-field boundary are then computed by solving the following system of equations:βƒ—π‘žπ‘=𝑅0ξ‹π‘Šπ‘Ÿ+Ξ“π‘žπ‘€βˆ’1⃗𝑆Δ𝑑.(39)

6.2. Inflow

For internal flows in rotating machineries, total conditions (total pressure 𝑃𝑑, total temperature 𝑇𝑑, and flow angle (πœ™π‘₯,πœ™π‘¦,πœ™π‘§)) are often provided at the inlet, and the following relations exist:πœŒπ‘=𝛾𝑀2π‘Ÿπ‘ƒπ‘‘π‘‡π‘‘ξ‚΅1βˆ’π›Ύβˆ’12𝑇𝑑𝑀2π‘Ÿπ‘ˆ2𝑏1/(π›Ύβˆ’1),𝑝𝑏=𝑃𝑑1βˆ’π›Ύβˆ’12𝑇𝑑𝑀2π‘Ÿπ‘ˆ2𝑏𝛾/(π›Ύβˆ’1),𝑒(40)𝑏=π‘ˆπ‘cosπœ™π‘₯,𝑣𝑏=π‘ˆπ‘cosπœ™π‘¦,𝑀𝑏=π‘ˆπ‘cosπœ™π‘§.(41)

For the subsonic inflow, the fourth characteristic equation with the wave propagating from the interior field to the boundary is used to close the above equations.

If the mass flow rate Μ‡π‘š, total temperature 𝑇𝑑, and the inflow angle 𝛼=(πœ™π‘₯,πœ™π‘¦,πœ™π‘§) are specified, (40) and the fourth characteristic equation are replaced by the following three relations:Μ‡π‘š=πœŒπ‘π‘ˆπ‘,𝑇𝑑=𝑇𝑏+π›Ύβˆ’12𝑀2π‘Ÿπ‘ˆ2𝑏,𝑝𝑏=πœŒπ‘π‘‡π‘π›Ύπ‘€2π‘Ÿ.(42)

The total velocity at the boundary point (π‘ˆπ‘) can be calculated based on the above equations.

6.3. Outflow

Characteristic variable boundary conditions are applied to the outflow boundary as well. For a subsonic outflow, the first four eigenvalues are positive and the fifth eigenvalue is negative. Flow variables at the boundary are obtained through the first four characteristic relations propogating from the interior point. The fifth relation is replaced by the static pressure specified at the exit:𝑝𝑏=𝑝exit.(43)

For a supersonic outflow, all characteristic waves propagate from the interior point to the outflow boundary. No flow variables will be specified at the exit.

6.4. Impermeable Wall

For an impermeable surface, the specified condition is that of no flow past through the surface, that is, Ξ˜π‘=0(44) or 𝑒𝑏𝑛π‘₯+𝑣𝑏𝑛𝑦+𝑀𝑏𝑛𝑧=βˆ’π‘›π‘‘.(45) The remaining boundary conditions are obtained based on the first four characteristic relations. It is worth to note that, in all the above boundary conditions, the source term contributions are added to the equations implicitly.

7. Results

Four test cases are chosen in this study to validate the efficacy of the proposed method for viscous flow computations, from extremely low to supersonic Mach numbers. These cases include a SUBOFF bare hull, a marine propeller 5168, a NACA 0012 rotor, and a waisted body. Cases with the SUBOFF bare hull and a waisted body are associated with nonrotating flows in the fixed frame, and cases with marine propeller 5168 and NACA 0012 rotor are associated with rotating flows in the rotating frame. The marine propeller 5168 is an internal flow that is used to validate the characteristic boundary conditions for inflow and outflow, as suggested in Section 6.

7.1. SUBOFF Bare Hull

The first case is concerned with applying the preconditioning method to predict an incompressible flow around SUBOFF bare hull, as shown in Figure 5. The SUBOFF model is a generic submarine model that has been extensively studied by various researchers. It was originally designed at David Taylor Research Center [27] to evaluate the accuracy of CFD tools. The validation data were provided by Roddy [28] and Huang et al. [29] using towing tank and wind tunnel measurements. In this study, a semispan geometric model was considered. The anisotropic grid was generated with boundary layer grids packing near the body. The 𝑦+ value was below 1 on the body surface. The free stream Mach number was chosen at 0.001 to mimic the incompressible flow effect. A wide range of flow angle of attack was chosen for comparison with the experimental data. The Reynolds number was 1.5 Γ— 107 based on the free stream velocity and the body length. The computations were performed with the arbitrary Mach number solver in the fixed frame. The preconditioning parameter for this case was chosen as follows:𝛽=min1,𝑀2π‘Ÿ+𝑀2Ξ©ξ€Έ,π‘€π‘Ÿ=0.001,𝑀Ω=0.0.(46)

Figure 5: SUBOFF bare hull and viscous mesh.

Figure 6 shows the convergence histories for the low Mach number flow computation without and with the preconditioning, measured by the 𝐿2 norm of the unsteady residuals for all flow variables. It is seen that the convergence was severely deteriorated if no preconditioning was applied. However, with the above preconditioning, the convergence was significantly improved, which was comparable to that obtained by the incompressible flow solver, as illustrated in Figure 6.

Figure 6: Convergence histories of SUBOFF bare hull.

The preconditioning not only affects the convergence but also improves the accuracy of the numerical solution at low Mach numbers. To demonstrate this, Figure 7 shows computed and measured axial force, normal force, and pitching moment coefficients. The axial force coefficient represents the drag acting on the body, which was the most difficult one to predict accurately. The error bars in the plots of experiments represent a targeted range of 10% off the measured values for the drag coefficient and 5% off the measured values for the lift and pitching moment coefficients. These plots show that the arbitrary Mach number solutions were compatible to the incompressible solutions, which were within the desired error range compared with the experimental data except at very high angles of attack (>15 deg.). Although not shown here, predicted force and moment coefficients without using preconditioning were considerably larger (10 times) than the measured values, due to the underlying dissipation of the numerical solution at this extremely low Mach number.

Figure 7: Comparison of predicted and measured force and moment coefficients.

The numerical study in this case clearly indicated the importance of the preconditioning to both convergence and accuracy of a numerical solution at low Mach numbers. It also showed that an incompressible flow may be accurately solved by using the preconditioned arbitrary Mach number method as described in the present study.

7.2. Marine Propeller 5168

The previous case considered a low Mach number viscous flow in a fixed frame, where the original preconditioning method was employed. To assess the efficacy of modified preconditioning for low Mach number rotating flows, a David Taylor marine propeller P5168 [30] was investigated. Figure 8 shows the geometry of the five-bladed, controllable-pitch marine propeller with a design advance coefficient of J = 1.27. The diameter of the propeller was 0.4027 meters with the test condition at a rotational speed of 1,450 RPM or tip speed of 30.57 m/s. The free stream velocity was 10.70 m/s. Four advance ratios were considered in the present study which are summarized in Table 1.

Table 1: Flow conditions of marine propeller p516.
Figure 8: Geometric and surface grid for marine propeller 5168.

The characteristic variables-based inflow and outflow boundary conditions were applied on the propeller 5168 grid boundaries. No-slip condition was applied to the propeller blade and shaft surface in the rotating frame. All solid surfaces were assumed to be adiabatic, that is, no heat fluxes across the surfaces. The Reynolds number was based on the free stream speed and the diameter of the propeller. The reference Mach number π‘€π‘Ÿ=0.1 was selected for all cases based on the free stream absolute velocity. The rotational Mach number 𝑀Ω was set to the propeller tip Mach number listed in Table 1 for the constant preconditioning option, while a local rotating speed in the rotating frame was used to determine the local rotating Mach number at each grid point. The constant preconditioning parameter for this case was chosen as follows:𝛽=min1,𝑀2π‘Ÿ+𝑀2Ξ©ξ€Έ,π‘€π‘Ÿ=0.1,𝑀Ω=𝑀tip.(47)

In the rotating frame, the flow can be viewed as steady state in the propeller coordinates. The local time stepping was used in the simulations with a CFL number of 10. One Newton iteration and eight Gauss-Seidel subiterations were used at each time step. Simulations were initiated from a uniform flow, and converged solutions were obtained in less than 1000 iterations. Numerical simulations were conducted using four preconditioning options: (1) no preconditioning, (2) original constant preconditioning, (3) modified local preconditioning, and (4) modified constant preconditioning. The convergence histories of the 𝐿2 norm of the unsteady residual for various preconditioning options were shown in Figure 9. The numerical test indicated that the original preconditioning (𝛽=π‘€π‘Ÿ2) did not produce a converged solution due to numerical instability (unable to obtain a solution). With the modified preconditioning of either local or constant rotating Mach number, the convergence of numerical solutions was noticeably improved. However, the constant preconditioning seemed to perform more superiorly than did the local preconditioning in the current test.

Figure 9: Convergence histories of marine propeller 5168.

Figure 10 shows the convergence histories of the thrust coefficients at four advance ratios. Predicted thrust coefficients were fully converged in less than 500 time steps due to the use of local time stepping in the rotating frame. Predicted thrust and torque coefficients were also compared with the experimental data in Figure 11. Excellent agreements between the computation and experiment were achieved at all four advance ratios using the current modified preconditioning method.

Figure 10: Convergence history of propeller thrust coefficient at different advance ratios.
Figure 11: Comparison of predicted and measured propeller thrust and torque.
7.3. NACA 0012 Rotor

The proposed preconditioning method was designed for not only the low Mach number or incompressible flows but also the transonic or supersonic flows. The following test case was about a Caradonna and Tung rotor [31] hovering at a transonic tip Mach number. Many researchers have used this NACA 0012 rotor and its test data for CFD validations. The rotor geometry was a two-bladed, untwisted, untapered rotor at an aspect ratio of 6 and collective pitch setting of 8 degrees; see Figure 12. The diameter of the rotor was 2.286 meters. The blades rotated at a speed of 2500 RPM with the blade-tip Mach number of 0.877, which was used as the reference Mach number. A Reynolds number was 4.6 Γ— 107 based on the blade-tip speed and the diameter of the blade. The following constant preconditioning parameter was chosen in the computations:𝛽=min1,𝑀2π‘Ÿ+𝑀2Ξ©ξ€Έ,π‘€π‘Ÿ=0.1,𝑀Ω=0.877.(48)

Figure 12: Caradonna and Tung rotor and surface mesh.

Simulations were carried out in the rotating reference frame using a local time stepping of a CFL number 10. A no-slip boundary condition was applied to the viscous surfaces, and the characteristic-variable far-field boundary condition was applied to the outer grid domain. The same four preconditioning options as the previous case were used in the current computation. However, since the reference Mach number was close to one, the effect of preconditioning was not very significant, as shown in convergence histories in Figure 13.

Figure 13: Convergence histories of Caradonna and Tung rotor.

Figure 14 compares computed and measured surface pressure coefficients at four radial span locations. All computed results were almost identical as the preconditioning effect was rather weak in this case. At the inner locations, the flow field was mainly occupied by the low speed flow. Moving further outboard to the blade tip, the flow began to show signs of transonic nature. A shock was developed at approximately 80% of rotor span location. There were slight deviations between the computation and experiment at 0.97% span location, due to a lack of grid resolution in the rotor tip path region. A locally refined grid in both the tip path plane and rotor downwash would significantly improve the overall prediction accuracy.

Figure 14: Computed and measured pressure coefficients at different blade spans.
7.4. Waisted Body

The arbitrary Mach number method was finally validated against a supersonic flow about a waisted body of revolution. The free stream was approaching at a Mach number of 1.4 and at a zero degree angle of incidence. Figure 15 shows the surface and the cutting plane volume grids of the waisted body that was measured by Winter et al. [32]. For supersonic flows, preconditioning was no longer in effect because the preconditioning parameter 𝛽 was always set to one. Nevertheless, the solution obtained by the arbitrary Mach number solver without preconditioning was not always equivalent to that obtained by the compressible flow solver. There were still several major differences between the two methods, including the use of primitive variables instead of conservative variables in the governing equations, a modified Roe scheme for the numerical flux approximation, and a normal derivative method for the viscous flux computation.

Figure 15: Waisted body and volume mesh.

The convergence histories for the compressible and preconditioned solutions (nopreconditioning) are shown in Figure 16. The residual was reduced by 3 orders of magnitude in about a 4000 time step in both cases using a CFL number of 5. Figure 17 shows computed pressure coefficient 𝐢𝑝 compared with the wind tunnel experimental data. Excellent agreement was achieved between the computation and the experiment. Figure 18 shows the skin fraction coefficient 𝐢𝑓 along the body length. The compressible solution showed a slightly deviation from the experiment data in the range from 10% to 50% of the body length. However, the preconditioned solution showed a favourite agreement with the experiment. These validations indicated that the current preconditioned method was equally applicable and robust to both low Mach number and supersonic flows in either fixed or rotating reference frame.

Figure 16: Convergence histories of a waisted body.
Figure 17: Comparison of predicted and measured 𝐢𝑝 distributions.
Figure 18: Comparison of predicted and measured 𝐢𝑓 distributions.

8. Conclusion

An improved preconditioning was presented and validated for solving viscous rotating flows based on unstructured grid topology. The original preconditioning method designed for nonrotating flows was found unstable for computing low Mach number flows involving large rotating speeds. A Fourier analysis revealed that the difficulty of the previous method was due to the collapse of the low wave frequency mode of the system matrix. The modified preconditioning based on both free stream and rotating Mach numbers has been validated against the experimental data, which showed a superior performance for both nonrotating and rotating flows in a wide range of flow regimes from extremely low Mach number to supersonic flows.


Governing Equations

Following the normalization in Section 2.1, the nondimensional relationships for ideal gas are𝑝=πœŒπ‘2𝛾,𝑐2=𝑇𝑀2r,(A.1)𝑒=πœŒπ‘’π‘‘=𝑀2π‘ŸπœŒξ‚΅π‘2𝛾+π›Ύβˆ’12π‘ˆ2ξ‚Άβ„Ž,(A.2)𝑑=𝑀2π‘Ÿξ‚΅π‘2+π›Ύβˆ’12π‘ˆ2ξ‚Ά=𝑒𝑑+πΈπ‘π‘πœŒ,(A.3) where 𝛾=𝐢𝑝/𝐢𝑣 is the constant specific heat, 𝐸𝑐=(π›Ύβˆ’1)𝑀2π‘Ÿ is an Eckert number. 𝑒𝑑 and β„Žπ‘‘ are the total energy and total enthalpy. π‘ˆ is the total velocity, and 𝑐 is the speed of sound.

The three-dimensional compressible governing equations, after introducing a preconditioning matrix Ξ“π‘žβˆ’1 and normalization, can be written in a differential form asπ‘€Ξ“π‘žβˆ’1πœ•βƒ—π‘ž+πœ•ξ‚€βƒ—ξ‹πΉπœ•π‘‘πΉβˆ’π‘£ξ‚+πœ•ξ‚€βƒ—ξ‹πΊπœ•π‘₯πΊβˆ’π‘£ξ‚+πœ•ξ‚€ξ‹ξ‹π»πœ•π‘¦π»βˆ’π‘£ξ‚=βƒ—πœ•π‘§π‘†,(A.4) where π‘ž is the vector of primitive variables and 𝐹, 𝐺, and 𝐻 are the vectors of convective fluxes:βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£πœŒπ‘’π‘£π‘€π‘βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,βƒ—βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£πœŒξ€·βƒ—π‘ž=𝐹=π‘’βˆ’π‘’π‘ ξ€ΈπœŒξ€·π‘’βˆ’π‘’π‘ ξ€ΈπœŒξ€·π‘’+π‘π‘’βˆ’π‘’π‘ ξ€Έπ‘£πœŒξ€·π‘’βˆ’π‘’π‘ ξ€Έπ‘€πœŒξ€·π‘’βˆ’π‘’π‘ ξ€Έβ„Žπ‘‘+π‘’π‘ πΈπ‘π‘βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,βƒ—βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£πœŒξ€·πΊ=π‘£βˆ’π‘£π‘ ξ€ΈπœŒξ€·π‘£βˆ’π‘£π‘ ξ€Έπ‘’πœŒξ€·π‘£βˆ’π‘£π‘ ξ€ΈπœŒξ€·π‘£+π‘π‘£βˆ’π‘£π‘ ξ€Έπ‘€πœŒξ€·π‘£βˆ’π‘£π‘ ξ€Έβ„Žπ‘‘+π‘£π‘ πΈπ‘π‘βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,ξ‹βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£πœŒξ€·π»=π‘€βˆ’π‘€π‘ ξ€ΈπœŒξ€·π‘€βˆ’π‘€π‘ ξ€Έπ‘’πœŒξ€·π‘€βˆ’π‘€π‘ ξ€Έπ‘£πœŒξ€·π‘€βˆ’π‘€π‘ ξ€ΈπœŒξ€·π‘€+π‘π‘€βˆ’π‘€π‘ ξ€Έβ„Žπ‘‘+π‘€π‘ πΈπ‘π‘βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,(A.5) where ⃗𝑛𝑑=(𝑒𝑠,𝑣𝑠,𝑀𝑠)𝑇=βƒ—Ξ©Γ—βƒ—π‘Ÿ is the linear velocity of the rotating coordinates. (𝐹𝑣,𝐺𝑣,𝐻𝑣)𝑇 are the vectors of diffusive fluxes which can be expressed as𝐹𝑣=⎑⎒⎒⎒⎒⎒⎒⎒⎒⎣0𝜏π‘₯π‘₯𝜏π‘₯π‘¦πœπ‘₯π‘§π‘’πœπ‘₯π‘₯+π‘£πœπ‘₯𝑦+π‘€πœπ‘₯π‘§βˆ’π‘žπ‘₯⎀βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,𝐺𝑣=⎑⎒⎒⎒⎒⎒⎒⎒⎒⎣0πœπ‘¦π‘₯πœπ‘¦π‘¦πœπ‘¦π‘§π‘’πœπ‘¦π‘₯+π‘£πœπ‘¦π‘¦+π‘€πœπ‘¦π‘§βˆ’π‘žπ‘¦βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,𝐻𝑣=⎑⎒⎒⎒⎒⎒⎒⎒⎒⎣0πœπ‘§π‘₯πœπ‘§π‘¦πœπ‘§π‘§π‘’πœπ‘§π‘₯+π‘£πœπ‘§π‘¦+π‘€πœπ‘§π‘§βˆ’π‘žπ‘§βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦.(A.6) The Reynolds shear stress and heat flux are given as𝜏π‘₯π‘₯=πœ‡2Re3ξ‚΅2πœ•π‘’βˆ’πœ•π‘₯πœ•π‘£βˆ’πœ•π‘¦πœ•π‘€ξ‚Ά,πœπœ•π‘§π‘₯𝑦=πœπ‘¦π‘₯=πœ‡ξ‚΅Reπœ•π‘’+πœ•π‘¦πœ•π‘£ξ‚Ά,πœπœ•π‘₯𝑦𝑦=πœ‡2Re3ξ‚΅2πœ•π‘£βˆ’πœ•π‘¦πœ•π‘’βˆ’πœ•π‘₯πœ•π‘€ξ‚Ά,πœπœ•π‘§π‘₯𝑧=πœπ‘§π‘₯=πœ‡ξ‚€Reπœ•π‘’+πœ•π‘§πœ•π‘€ξ‚,πœπœ•π‘₯𝑧𝑧=πœ‡2Re3ξ‚΅2πœ•π‘€βˆ’πœ•π‘§πœ•π‘’βˆ’πœ•π‘₯πœ•π‘£ξ‚Ά,πœπœ•π‘¦π‘¦π‘§=πœπ‘§π‘¦=πœ‡ξ‚΅Reπœ•π‘£+πœ•π‘§πœ•π‘€ξ‚Ά,π‘žπœ•π‘¦π‘₯πœ‡=βˆ’Reπ‘ƒπ‘Ÿπœ•π‘‡,π‘žπœ•π‘₯π‘¦πœ‡=βˆ’Reπ‘ƒπ‘Ÿπœ•π‘‡,π‘žπœ•π‘¦π‘§πœ‡=βˆ’Reπ‘ƒπ‘Ÿπœ•π‘‡,πœ•π‘§(A.7) where πœ‡=πœ‡l+πœ‡t is the viscosity, π‘ƒπ‘Ÿ=π‘ƒπ‘Ÿπ‘™+π‘ƒπ‘Ÿπ‘‘ is the Prandt number, and Re is the Reynolds number based on the reference length and velocity. The source term in (A.2) is givenβƒ—βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£0𝑆=πœŒπ‘€Ξ©π‘¦βˆ’πœŒπ‘£Ξ©π‘§πœŒπ‘’Ξ©π‘§βˆ’πœŒπ‘€Ξ©π‘₯πœŒπ‘£Ξ©π‘₯βˆ’πœŒπ‘’Ξ©π‘¦0⎀βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,(A.8) where Ξ©x,Ξ©y,andΞ©z are the components of the rotational speed of the coordinate system.


𝐴:Surface area of control volume
𝐴, 𝐡, 𝐢:Matrices of flux Jacobians with respect to conservative variables
π‘Ž,𝑏,𝑐:Matrices of flux Jacobians with respect to primitive variables
𝑐:Speed of sound
𝐸𝑐:Eckert number
⃗⃗𝐻𝐹,𝐺,:Components of inviscid or viscous flux vector in π‘₯, 𝑦, and 𝑧 directions
β„Žπ‘‘:Total enthalpy
π‘˜:System matrix on the face of control volume
𝑀:Transformation matrix
π‘€π‘Ÿ:Free stream reference Mach number
𝑀Ω:Rotating reference Mach number
⃗𝑛:Normal unit vector
𝑃:General flux vector
𝑃Γ:Preconditioned inviscid flux vector
𝑃𝑣:Viscous flux vector
π‘π‘Ÿ:Prandtl number
⃗𝑄:Conservative flow variables
βƒ—π‘ž:Primitive flow variables
π‘ž:Heat flux
𝑅:Right eigenvector
Re:Reynolds number
⃗𝑆:Source term in rotating frame
βƒ—π‘ˆ:Velocity vector
𝑒. 𝑣. 𝑀:Velocity components
𝑒𝑠. 𝑣𝑠. 𝑀𝑠:Grid velocity components
𝑒. 𝑣. 𝑀:Velocity components
𝑉:Control volume
ξ‹π‘Š:Vector of characteristic variables
𝑀1βˆ’5:Component of characteristic variables
π‘₯, 𝑦, 𝑧:Cartesian coordinates.

Greek Letter

Ξ“π‘ž:Preconditioning matrix
𝛽:Preconditioning parameter
𝛾:Specific heat ratio
Θ:Normal velocity to surface
ΛΓ:Diagonal matrix containing real eigenvalues
πœ†1βˆ’5:Eigenvalues of system matrix π‘˜
𝜏:Viscous stress tensor
Ξ©:Rotational speed of coordinate system.


0:Locally constant value
𝑏:Boundary point of computational domain
𝑖:Interior point of computational domain
π‘œ:Outer point of computational domain
π‘Ÿ:Reference point
𝐿, 𝑅:Left and right states of control volume face
π‘₯, 𝑦, 𝑧:Direction of Cartesian coordinates.


  1. A. J. Chorin, β€œA numerical method for solving incompressible viscous flow problems,” Journal of Computational Physics, vol. 2, no. 1, pp. 12–26, 1967. View at Google Scholar Β· View at Scopus
  2. W. R. Briley, H. McDonald, and S. J. Shamroth, β€œA low Mach number euler formulation and application to time iterative LBI schemes,” AIAA Journal, vol. 21, no. 10, pp. 1467–1469, 1983. View at Google Scholar Β· View at Scopus
  3. D. Choi and C. L. Merkle, β€œApplication of time–iterative schemes to incompressible flow,” AIAA Journal, vol. 23, no. 10, pp. 1518–1524, 1985. View at Google Scholar Β· View at Scopus
  4. D. Choi and C. L. Merkle, β€œThe application of preconditioning in viscous flows,” Journal of Computational Physics, vol. 105, no. 2, pp. 207–223, 1993. View at Publisher Β· View at Google Scholar Β· View at Scopus
  5. E. Turkel, β€œPreconditioned methods for solving the incompressible and low speed compressible equations,” Journal of Computational Physics, vol. 72, no. 2, pp. 277–298, 1987. View at Google Scholar Β· View at Scopus
  6. B. van Leer, W. T. Lee, and P. L. Roe, β€œCharacteristic time-stepping or local preconditioning of the euler equations,” Tech. Rep. 91-1552, AIAA, 1991. View at Google Scholar
  7. J. M. Weiss and W. A. Smith, β€œPreconditioning applied to variable and constant density flows,” AIAA Journal, vol. 33, no. 11, pp. 2050–2057, 1995. View at Google Scholar Β· View at Scopus
  8. J. M. Weiss, J. P. Maruszewski, and W. A. Smith, β€œImplicit solution of preconditioned Navier-Stokes equations using algebraic multigrid,” AIAA Journal, vol. 37, no. 1, pp. 29–36, 1999. View at Google Scholar Β· View at Scopus
  9. E. Turkel, β€œReview of preconditioning methods for fluid dynamics,” Applied Numerical Mathematics, vol. 12, no. 1–3, pp. 257–284, 1993. View at Google Scholar Β· View at Scopus
  10. W. R. Briley, L. K. Taylor, and D. L. Whitfield, β€œHigh-resolution viscous flow simulations at arbitrary Mach number,” Journal of Computational Physics, vol. 184, no. 1, pp. 79–105, 2003. View at Publisher Β· View at Google Scholar Β· View at Scopus
  11. X. Wang, C. Sheng, and J. Chen, β€œNumerical study of preconditioned algorithm for ratational flows,” Tech. Rep. 2005-5126, AIAA, Toronto, Canada, 2005. View at Google Scholar
  12. C. Sheng and X. Wang, β€œA global preconditioning method for low mach number viscous flows in rotating machinery,” in Proceedings of the ASME Turbo Expo, vol. 6, pp. 1555–1564, Barcelona, Spain, May 2006.
  13. R. Vichnevetsky and J. B. Bowles, Fourier Analysis of Numerical Approximations of Hyperbolic Equations, Society for Industrial Mathematics, 1982.
  14. D. Shoe and C. J. Knight, β€œComputation of 3D viscous flows in rotating turbomachinery blades,” Tech. Rep. 89–0323, AIAA, Reno, Nev, USA, 1989. View at Google Scholar
  15. J. P. Chen, A. R. Ghosh, D. Sreenivas, and D. L. Whitfield, β€œComparison of computations using navier–stokes equations in rotating and fixed coordinates for flow through turbomachinery,” Tech. Rep. 97–0878, AIAA, Reno, Nev, USA, 1997. View at Google Scholar
  16. J. T. Batina, β€œImplicit flux–split euler schemes for unsteady aerodynamic analysis involving unstructured dynamic meshes,” Tech. Rep. 90–0936, AIAA, 1990. View at Google Scholar
  17. W. K. Anderson, β€œGrid generation and flow solution method for euler equations on unstructured grids,” NASA Technical Report TM 4295, 1992. View at Google Scholar
  18. W. K. Anderson, R. D. Rausch, and D. L. Bonhaus, β€œImplicit/multigrid algorithms for incompressible turbulent flows on unstructured grids,” Tech. Rep. 95–1740, AIAA, June 1995. View at Google Scholar
  19. D. G. Hyams, An investigation of parallel implicit solution algorithms for incompressible flows on unstructured topologies, Ph.D. Dissertation, Mississippi State University, May 2000.
  20. W. J. Coirier, β€œAn adaptively–refined, cartesian, cell–based scheme for the euler and navier–stokes equations,” NASA Technical Memorandum 106754, NASA Lewis Research Center, 1994. View at Google Scholar
  21. P. Spalart and S. Allmaras, β€œA one–equation turbulence model for aerodynamic flows,” Tech. Rep. 92–0439, AIAA, January 1991. View at Google Scholar
  22. P. R. Spalart and M. Shur, β€œOn the sensitization of turbulence models to rotation and curvature,” Aerospace Science and Technology, vol. 1, no. 5, pp. 297–302, 1997. View at Google Scholar
  23. W. W. Liou and T. H. Shih, β€œTransonic turbulence flow predictions with new two–equation turbulence models,” Tech. Rep., Institute for Computational Mechanics in Propulsion, January 1996. View at Google Scholar
  24. D. C. Wilcon, Turbulence Modeling for CFD, DCW Insdustries, La Cañada, Calif, , USA, 1994.
  25. T. J. Coakley, β€œTurbulence modeling methods for the compressible navier–stokes equations,” Tech. Rep. 83–1693, AIAA, 1983. View at Google Scholar
  26. C. Sheng and X. Wang, β€œCharacteristic variable boundary conditions for arbitrary Mach number algorithm in rotating frame,” in Proceedings of the 16th AIAA Computational Fluid Dynamics Conference, Orlando, Fla, USA, June 2003.
  27. N. C. Groves, T. T. Huang, and M. S. Chang, β€œGeometric characteristics of DARPA SUBOFF models,” Tech. Rep. DTRC/SHD-1298-01, Defense Advanced Research Projects Agency, 1989. View at Google Scholar
  28. R. F. Roddy, β€œInvestigation of the stability and control characteristics of several configurations of the DARPA SUBOFF model,” Tech. Rep. DTRC/SHD-1298-08, Defense Advanced Research Projects Agency, 1990. View at Google Scholar
  29. T. T. Huang, H. L. Liu, N. C. Groves, T. J. Forlini, J. N. Blanton, and S. Gowing, β€œMeasurements of flows over an axisymmetric body with various appendages in a wind tunnel: the DARPA SUBOFF experimental program,” in Proceedings of the 19th Symposium on Naval Hydrodynamics, pp. 321–346, 1992.
  30. C. J. Chesnakas and S. D. Jessup, β€œCavitation and 3–D LDV tip–flowfield measurements of propeller 5168,” Tech. Rep. CRDKNSWC/HD–1460–02, Carderock Division, Naval Surface Warfare Center, May 1998. View at Google Scholar
  31. F. X. Caradonna and C. Tung, β€œExperimental and analytical studies of a model helicopter rotor in hover,” Vertica, vol. 5, no. 2, pp. 149–161, 1981. View at Google Scholar
  32. K. G. Winter, J .C. Rotta, and K. G. Smith, β€œStudies of the turbulent boundary layer on a waisted body of revolution in subsonic and supersonic flow,” Technical Report 68215, Royal Aircraft Establishment, 1968. View at Google Scholar