Abstract
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 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 ; 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 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
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: Equation (3) can be written in terms of primitive variables as or where , , and are the system matrices. Matrices , , and are the Jacobian matrices of the convective flux vectors ()with respect to conservative variables . Introducing a preconditioning matrix into the time derivative term of (5) yields or
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: 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 The system matrix has five real eigenvalues: where , is the grid speed in the normal direction of the face, , and . 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 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 (), 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 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 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: which scales 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.
(a) β = 1
| (b) β = |
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.
| (a) β = |
| (b) |
It is therefore desired to scale the above acoustic term 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: 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 (). 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 . The test geometry is a simple body of revolution with a free stream velocity at an extremely low Mach number . The characteristic rotating Mach number is , 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 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.
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.
The system of governing equations is numerically integrated over the closed boundaries (median dual faces) of a control volume surrounding each node: 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 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 where the subscript denotes the vertex and 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 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: 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 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 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: where 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: The derivative in the normal direction can be approximated by where and 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, where and 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 as where , . is the time step increment between steps and , 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 . Correspondingly, a second-order time-accurate Euler implicit scheme is given by .
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 Newton’s method solves the following equation: where = 1,2,3 are the Newton steps, with . is the Jacobian matrix of (26), that is,
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 in (16) as a local constant; the approximate flux Jacobian for the inviscid flux can be written as
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 is obtained through a sequence of iterations, , which converge to . 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: where
Letting be the vector of unsteady residuals and the change in the dependent variables, the symmetric Gauss-Seidel relaxation can be written as the following two-step process:
In the forward pass, are obtained with the previously updated , which were set to zero at the initial stage. In the backward pass, are obtained with the most recent value of 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 (working variable) is expressed as 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 where is a diagonal matrix made up of the eigenvalues of and and are matrices whose columns and rows are the right and left eigenvectors of , respectively. Multiplying both sides of (7) by gives the following characteristic relations in one-dimensional space where 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: For governing equations containing a source term, one has or 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:
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:
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:
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:
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, or 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:
Figure 6 shows the convergence histories for the low Mach number flow computation without and with the preconditioning, measured by the 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.
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.
(a) Axial force coefficients
(b) Normal force coefficients
(c) Pitching 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.
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 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:
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 norm of the unsteady residual for various preconditioning options were shown in Figure 9. The numerical test indicated that the original preconditioning () 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 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.
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:
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 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.
| (a) distribution at 68% span |
| (b) distribution at 80% span |
| (c) distribution at 89% span |
| (d) distribution at 96% span |
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.
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.
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.
Appendix
Governing Equations
Following the normalization in Section 2.1, the nondimensional relationships for ideal gas are where is the constant specific heat, 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 and normalization, can be written in a differential form as where is the vector of primitive variables and , , and are the vectors of convective fluxes: where is the linear velocity of the rotating coordinates. are the vectors of diffusive fluxes which can be expressed as The Reynolds shear stress and heat flux are given as where 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 where are the components of the rotational speed of the coordinate system.
Nomenclature
| : | 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 |
| : | Pressure |
| : | Prandtl number |
| : | Conservative flow variables |
| : | Primitive flow variables |
| : | Heat flux |
| : | Right eigenvector |
| Re: | Reynolds number |
| : | Source term in rotating frame |
| : | Temperature |
| : | Time |
| : | Velocity vector |
| . . : | Velocity components |
| . . : | Grid velocity components |
| . . : | Velocity components |
| : | Control volume |
| : | Vector of characteristic variables |
| : | Component of characteristic variables |
| , , : | Cartesian coordinates. |
Greek Letter
| : | Preconditioning matrix |
| : | Preconditioning parameter |
| : | Specific heat ratio |
| : | Normal velocity to surface |
| : | Diagonal matrix containing real eigenvalues |
| : | Eigenvalues of system matrix |
| : | Viscosity |
| : | Density |
| : | Viscous stress tensor |
| : | Rotational speed of coordinate system. |
Subscript
| 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. |