An unconditionally stable Runge-Kutta method for unsteady flows

PHILIP JORGENSON, RODRICK CHIMA
1989 27th Aerospace Sciences Meeting   unpublished
A quasi-three-dimensional analysis has been developed for unsteady rotor-stator interaction in turbomachinery. The analysis solves the unsteady Euler or thinlayer NavierStokes equations in a body-fitted coordinate system. It accounts for the effects of rotation, radius change, and stream-surface thickness. The Baldwin-Lomax eddy-viscosity model is used for turbulent flows. The equations are integrated in time using a four-stage Runge-Kutta scheme with a constant time step. Implicit residual
more » ... plicit residual smoothing has been employed to accelerate the solution of the time-accurate computations. The scheme is described and stability and accuracy analyses are given. Results are shown for a supersonic throughflow fan designed at NASA Lewis Research Center. The rotorstator blade ratio was taken as 1:l. Results are also shown for the first stage of the Space Shuttle Main Engine high pressure fuel turbopump. Here the blade ratio is 2:3. Implicit residual smoothing was used to increase the time step limit of the unsmoothed scheme by a factor of six with negligible differences in the unsteady results. We feel that the implicitly smoothed Runge-Kutta scheme is easily competitive with implicit schemes for unsteady flows while retaining the simplicity of an explicit scheme. time steps approaching those used by implicit schemes but retains the speed and simplicity of explicit schemes. Many of the numerical tools used in the analysis of isolated blades can be used for time-accurate rotor-stator interaction analysis. Here the quasi-three-dimensional Euler or thin-layer NavierStokes equations are solved on body fitted C-grids. An explicit four-stage Runge-Kutta algorithm is used with a constant time step. An implicit residual smoothing procedure is used to increase the time step so that a converged periodic solution is reached in fewer iterations. Computer time is reduced by vectoritation. Most turbomachinery analyses assume periodicity from blade to blade. This makes it possible to analyze one blade only. Turbomachinery is designed with unequal number of blades to avoid forced vibration problems. Therefore a full unsteady analysis must include at least a few passages on each wheel. The numerical method is presented with viscous results for a supersonic throughflow (STF) fan designed a t NASA Lewis Research Center. Comparisons are made between viscous results with and without residual smoothing for a scaled 1:l rotor:stator blade configuration of the fan. The solution procedure has also been applied to the first stage of the Space Shuttle Main Engine (SSME) high pressure turbopump. Results are presented on a 2:3 stator:rotor blade configuration of the turbopump. Introduction Governing Equations The major thrust of the computational analysis of turbomachinery to date has been the steady state solution of isolated blades using mass-averaged inlet and exit conditions. Unsteady flows differ from the steady solution due to interaction of pressure waves and wakes between blade rows. To predict the actual complex flow conditions, one must look a t the time-accurate solution of the entire turbomachine. Several solution procedures have been developed to solve unsteady flow in turbomachinery(1-6) including our use of an explicit Runge-Kutta scheme as described in (7). A limitation of an explicit scheme in predicting time accurate flows is the small time step necessary to maintain numerical stability. Residual smoothing was introduced by Lerat(8) for application to the Lax-Wendroff scheme. Jameson(9) and others( 10-11) have applied a similar techique to Runge-Kutta Schemes for steady flows. The residual smoothing technique allows a larger time step to be taken without affecting the accuracy of the scheme. The motivation behind the present work is to use implicit residual smoothing to greatly increase the small time step required for stability of the Runge-Kutta scheme, while maintaining second order accuracy in both time and space. The resulting scheme allows The axisymmetric (m, 9) coordinate system used for the quasi-three-dimensional analysis is shown in Fig. 1 . Here the m-coordinate is defined by (dm)' = ( d~)~ + (dr)' and the 0-coordinate is defined by: where 0' is fixed in space and 9 rotates with the blade row with angular velocity n. The radius r and the stream surface thickness h are taken to be known functions of m. The flow equations are written in the (m, 9) system (see Refs. 12 and 13) and are then transformed to a general body-fitted (e, r]) system using standard methods. The thin-layer approximation is used to eliminate all viscous derivatives in the streamwise (0 direction. The final equations are as follows: (2) 1 A The inviscid flux terms are given by: The relative contravariant velocity components We and W'J along the [ and r] grid lines are given b y W' = [mu, + Eewe; Wv = qmum + Vewe ( 5 ) Here we = component. The energy and pressure are given by: -rR is the relative tangential velocity e = p[C,T + 1 / 2 ( 4 + ugZ)] p = (7 -1)le -1 / 2 ( 4 + u i ) ] (6) (7) The source term Kz arises from the centrifugal force term in the r-component of the m-momentum equation. K2 = (pu~+p-Re-'~~~)r,/r+(p-Re-~a~~)h,/h (8) where 1 dr r dm r,/r = --1 dh h dm h,/h = --Although a similar source term arises from the Coriolis force in the 8-momentum equation, the equation has been made conservative by multiplying through by r. The viscous flux term is r o i where and the shear stress terms are given by: Using the thin-layer approximation, the shear stress terms are evaluated by replacing a, with qmao and (l/r)a, with Teatl. Here a = is the sonic velocity, and the normalized thermal conductivity k equals one. The equations are nondimensionahed by arbitrary reference quantities (here the inlet total density and critical sonic velocity define the reference state), and the Reynolds number Re and the Prandtl number Pr must be specified in terms of that state. These equations assume that the specific heats C, and C, and the Prandtl number are constant, that Stokes' hypothesis X = -2/3p is d i d , and that the effective viscosity may be written as The transformation metrics are found using where the Jacobian is given by Overbars in Eqs. (4-9) denote a rescaling of the metrics: -€0 = €e/r; r e = Ve/r; 7' = rhJ-' For turbulent flows the two-layer eddy-viscosity model developed by Baldwin and Lomax (14) is used. In the (m,6) coordinate system the wall shear rw and vorticity w required by the model are given by r w = ~1 2 w = p(amVe + l/raeumuerm/r)w 1 2 w = -(amwe -I/raeu, + ugrm/r) (14) The eddy viscosity is computed based on the current flow field, but otherwise does not account for unsteady turbulence effects. Computational Grid Body-fitted C-type grids for this work were generated using the GRAPE code developed by Sorenson (15,16). Figure 2 shows typical grids around the first stage stator and rotor of the STF fan. The stator grid was modified by adding a single grid line parallel to the exit to allow approximately one cell of overlap with the rotor grid. The rotor grid was also selectively refined by doubling the number of [=constant grid lines in the inlet region. Both grids were modified by adding one line of dummy points (not shown) corresponding to interior points from neighboring passages. These points are used for imposing the periodic or overlap boundary conditions. Multistage Runge-Kutta Algorithm An explicit multistage Runge-Kutta algorithm based on the work of Jameson (9,17) is used to advance the flow 2 equations in time. Given the residual R from a centraldifference representation of Eq. (3), a k-stage scheme may be written as follows: where D(q) contains the physical and artificial dissipation terms. The time step is given by At = min(Ati,j) where A t , j is calculated from an inviscid stability analysis of where I , and 10 are reciprocal length scales given by In this work we have used a four-stage scheme with the standard coefficients ai=(1/4, 1/3, 1/2, 1). The classical fourth-order Runge-Kutta scheme may be found in most books on numerical methods (see 18 for example) and is given by: The classical scheme matches a Taylor-series expansion of q"+' about q" in time through terms of order (At)', and so is fourth-order accurate. It is stable for Courant numbers C F L 5 2&, and requires one more level of storage than the simplified scheme(l5). The extra storage is needed to accumulate the right hand side for the last stage. For linear problems the simplified scheme with standard coefficients (15) corresponds to the classical scheme (17). For non-linear problems the simplified scheme (15) only matches an expansion of q"+' through terms of order (At)2, and is thus only second-order accurate in time. Furthermore, as shown in (U), we evaluate the artificial and physical viscous terms only at the first stage, so that our scheme is formally second-order accurate in time for the convective terms but only first-order accurate for the viscous terms. Full second-order accuracy can be achieved by another evaluation of the viscous terms. With our use of the thin-layer approximation and with relatively coarse grids between the blades we feel that it is probably not worthwhile to use second-order temporal accuracy for the viscous terms. However we point out that fourth-order accuracy is easily obtained using the classical scheme and will be given further attention in the future. Artificial Dissipation Dissipative terms consisting of fourth and second differences are added to prevent odd-even point decoupling and to allow shock capturing respectively. The dissipative term D in Eq. (15) is given by: The (-direction operator is given by: is a coefficient that partially cancels similar terms in Eq. (15). We have found that using the spatially-varying time step Atij (Eq. (16)) as a coefficient in the dissipation is much less dissipative than using the minimum time step. The terms V2 and V4 are given by: where and P2 = O(1) (23) In smooth regions of the flow the dissipative terms are of third order and thus do not detract from the formal second-order accuracy of the scheme. Near shocks ui,j is large and the seconddifference dissipation becomes locally of first order. Implicit Residual Smoothing The four-stage Runge-Kutta scheme described above has a Courant limit of 2.8, but we typically run the unsteady results at a slightly conservative value of 2.5. For a constant At the Courant number is dominated by the grid spacing, so that the maximum Courant number nearly always occurs in the very fine viscous grid regions around the leading or trailing edge of a blade. Courant numbers in the inviscid parts of the flow are typically two orders of magnitude smaller than the maximum value, 3 i.e. O(0.02). We feel that the Courant numbers in the inviscid part of the flow could be increased to O(1) without loss of accuracy in acoustic wave propogation or wake convection between blade rows. In the viscous part of the flow the maximum Courant number that gives reasonable time-accuracy is questionable, especially with the turbulence model applied in a quasi-steady sense. Overall we suspect that maximum Courant numbers of O(l0-20) in viscous regions would give a good balance between accuracy and efficiency. Residual smoothing was introduced by Lerat (see for example Ref. 8) for use with the Lax-Wendroff scheme and was later applied to Runge-Kutta schemes by Jameson (9). The technique involves running the multistage scheme at a Courant number, C F L , greater than the stability limit of the scheme CFL'. The scheme is stabilised by smoothing the residual using an implicit operator after each stage. If we rewrite stage (k) of the multistage scheme, Eqn. (W, as then the implicit smoothing step is given by where 6" and 6" are seconddifference operators and e is a smoothing parameter. Linear stability analysis has shown that the Runge-Kutta scheme with implicit residual smoothing may be made unconditionaly stable provided that e is sufficiently large (9). In one dimension f [(m)z-l] C F L gives unconditional stabiltiy. In two dimensions different values of e are often used in each grid direction. In one dimension the implicit smoothing step can be written as (27) (1 -eAz26")G = Aq = O(At(JP6 + ...)) and since At = O(Az) = Aq + O(cAz3) Thus implicit residual smoothing modifies the unsteady solution by a term of order eAz3. As long as e is not too large the use of implicit residual smoothing does not upset the formal second-order accuracy in space. Since the smoothing is done in space it has no effect on the temporal accuracy of the scheme. To date implicit residual smoothing has been used in conjunction with a spatially varying time step to accelerate convergence to a steady state. In this case the local Courant number, C F L , is constant so that the smoothing parameter e is constant over the grid. For unsteady flows At is constant and the Courant number varies greatly over the grid, as mentioned previously. If a constant value of c is chosen based on the largest value of CFLICFL' the scheme remains stable, but the inviscid part of the flow, which is well within the stability limit, is grossly oversmoothed and much of the flow physics is lost. If, however, a local smoothing parameter eij is calculated using then the residuals are smoothed no more than needed for stability. In effect, the viscous part of the flow is calculated at a large Courant number using implicit residual smoothing and the inviscid part of the flow is calculated a t a small Courant number without smoothing. Martinelli and Jameson (19) have used a spatially varying smoothing coefficient to account for highlystretched viscous grids. They use an equation similar to (29) in each coordinate direction, with C F L i j replaced by a scaled onedimensional Courant number. For typical viscous flows this minimizes the smoothing in the flow direction when the stability of the scheme is determined by the grid spacing in the normal direction. This type of directionally-biased smoothing could easily be applied to unsteady flow. Implementation of Eqn. (29) results in a scalar tridiagonal matrix equation that must be solved along each grid line in each direction, after each stage of the multistage scheme. The recursive tridiagonal solutions may be vectoriced across one grid direction while solving in the other direction. Implicit residual smoothing adds 26 percent to the CPU time of a solution, but we have used it to raise the time step by a factor of 6.52, for a net decrease in CPU time of about 460 percent. Finally we have a spatially-varying smoothing coefficient with a large constant At (maximum C F L FY 15) to accelerate steady calculations. In general we find that this scheme converges more slowly than one using a constant smoothing coefficient with a spatially-varying At ( C F L = 5.6 =constant). Boundary Conditions Inlet At the inlet to the stator, total pressure, total temperature T', and the whirl rvg are specified. At each time step the upstream-running Riemann invariant R -= urn -2a/(7 -1) is extrapolated to the inlet. The axial velocity component urn is found using urn = (7 -1)R-+ d(7 + 1)(4C"T' -2 4 -2(7l)(R-)z (7-t 1) (28) Density and energy are found using isentropic relations. This is a non-reflective approximation to the axial momentum equation that is first-order accurate in space but 4 sero-order in time. However in the present results we see virtually no unsteady behavior at the inlet boundary. Outlet At the exit from the rotor, static pressure is specified and the other flow quantities are found using secondorder extrapolation of p, pwm, and pue. The interaction between the specified exit pressure and the inlet conditions sets the mass flow through the machine. This exit condition is strictly reflective and may influence the unsteady solution. We hope to sort out the effects of this condition in future work. Solid Wall Blade surface pressures are found from the normal momentulr equation: ( € m~m + €efje)atP+ ( T : + V V~P Surface velocities are found from the tangency or noislip conditions for inviscid or viscous flows, and the blade temperature is specified. Periodic Boundaries The periodic boundaries are solved in the same way as the interior points. The grid for each blade is overlaped one grid cell at its periodic boundary such that the grid points are coincident with the grid points of an adjacent blade. These overlap dummy grid points are updated when an adjacent blade is solved and are then used in the next time step integration. For equal pitch blade rows the flow variables on the spatially periodic boundaries are equal on each blade. This is not true for multi-passage calculations. Here, periodic boundary conditions are enforced only on boundaries of blades that are encompassed by the blade ratio. Flow variables on interior boundaries are updated using quantities from adjacent blade solutions. Interface Formulation For a given displacement between stator and rotor blade rows, the solution on the interface must be determined. Figure 3 shows an enlargement of the interface geometry of the stator and rotor grids. The upstream boundary of the rotor is coincident with the ( = 2 and < = M -1 grid lines in the stator grid. The stator grid extends into the rotor grid such that the exit boundary grid points of the stator lie within the first cell of the rotor domain. The solution is updated at the interface by interpolating the flux variables in the stator computational domain to obtain rotor quantities and interpolating the flux variables in the rotor computational domain to obtain stator quantities. The interface is updated in this way after every Runge-Kutta integration sweep of the stator and rotor domains. Solution Procedure Given a desired leading edge velocity triangle and exit flow angle for each blade, the initial conditions are calculated using an analytic solution of the onedimensional flow equations (see Ref. 7). These initial conditions do not especially speed the convergence to a periodic unsteady solution, but they do provide a smooth transition betmen the inlet and exit boundary conditions. From the initial analytic flow conditione, the solution is integrated using a local time step. This allows the mass flow to stabilice and the wake region to develop quickly. A constant minimum time step is then used to compute the time accurate solution. The solution for each blade row is updated a full time step (four Runge-Kutta stages), then the inflow, exit, blade, and interface boundaries are updated. The solution is considered converged when the loading on the blades becomes periodic in time. The solution for each passage is stored as a separate array in core and accessed as the solver alternates from one blade to the next. The flow solver is relatively unaffected by the number of passages so that the multi-passage problem becomes a problem of data management. Results In 1985 Rai published Euler calculations of the unsteady interaction between a model rotor and stator with supersonic throughflow. The calculations were meant to demonstrate the capabilities of his numerical scheme and did not represent a real Supersonic throughflow machine. Interest has recently increased a t NASA Lewis Research Center in a supersonic throughflow (STF) fan as a key component in an engine designed for sustained supersonic cruise. Such a fan has been designed at Lewis and is now being built for subsequent testing (20,21). The solution procedure has been applied to the STF fan. The actual STF fan rotor:stator blade count ratio is 44:52, but the numerical procedure was applied to a single passage blade row configuration where the upstream rotor was scaled to the stator pitch. The scaling was done such that the blade angles and pitch-to-chord ratio remained unchanged. The grids generated using the GRAPE code have viscous spacings of 0.0008 in. away from the blade. The rotor grid has 1 5 3 x 3 1 points. The stator grid has 1 7 9 x 3 1 grid points. It was necessary to cluster grid points at the inlet of the stator to get the necessary resolution to capture the wake from the upstream rotor. The grids are shown in Fig. 2 with every other point deleted for clarity. The absolute Mach number at the inlet to the rotor is 2.00. The code was run a t a maximum Courant number of 2.3 without residual smoothing. With this time step it takes 2780 steps to move one rotor pitch. A solution is considered converged when the unsteady lift on each blade becomes periodic. For supersonic flow the downstream stator has no effect on the upstream rotor, so that the rotor lift becomes constant after about four 5 blade passings, as shown in Fig. 4 . There is however an unsteady effect on the stator as the rotor wakes pass by, and the stator lift becomes periodic in time after eight blade passings. This requires 55 minutes of computer time on the Cray X-MP. The code was also run at a maximum Courant number of 15 with residual smoothing. With this time step it takes 425 steps to move one rotor pitch. The number of blade passings required to obtain a converged solution remained unchanged from the solution without residual smoothing. However, the computer time was reduced to 12 minutes. Figure 5 shows the static pressure distribution on the rotor. Since the flow is supersonic, the rotor is unaffected by the passing stator and the maximum, minimum, and time-averaged pressures are identical. As shown in Figure 6 , the stator does see fluctuations in static pressure as it passes through the wakes and the weak fishtail shocks that form off the trailing edge of the upstream rotor. In Figure 6 the solid line represents the time-averaged pressure distribution and the dashed lines indicate the maximum and minimum pressures at each point. On both blades the static pressure envelopes are identical whether computed with or without implicit residual smoothing. Mach contours relative to the rotor frame of reference are shown in Fig. 7 . Weak shocks form on the leading edge of both blades. The wake from the rotor convects downstream and interacts with the leading edge shock of the stator. The wake that forms off the trailing edge of the stator is distorted by rotor wakes as they are convected through the downstream boundary. The rotor wakes are more distinguishable in the stator passage when entropy contours are plotted in Fig. 8 . The distortion of the stator wake is shown to coincide with the passing rotor wakes. The solution procedure has also been applied to the firststage of the Space Shuttle Main Engine fuel t u r b o p ump. The blade ratio of the actual machine is 41:63. The stator was rescaled to produce the 2:3 blade ratio shown in Fig. 9 . The pitch-to-chord ratio remained unchanged. The grids generated using the GRAPE code have viscous spacings of 0.0002 in. away from the blade. The stator grid has 115x31 points. The rotor grid has 197x41 grid points. The grid was clustered around the rotor inlet to capture the wake from the upstream stator. The code runs at a maximum Courant number of 2.5 without residual smoothing. With this Courant number it takes 15510 time steps to move one periodic pitch (three passing rotor blades). A converged solution was reached in about seven rotor blade passings. This took about six hours of computer time. Details of results from this calculation were published in (7). A solution was also computed with the use of residual smoothing. A maximum Courant number of 15 was used which decreased the number of interations for the solution to move one periodic pitch to 2585. The solution with residual smoothing was identical to the solution without residual smoothing. The computer time was reduced to about 80 minutes. The unsteady lift diagrams for the scaled first stage of the SSME fuel turbopump is shown in Fig. 10 . The stator blade sees small fluctuations in lift due to deviations in static pressure near the trailing edge on the suction surface. Over one periodic pitch the stator lift shows three flucuations as it is influenced by three rotors passing downstream. The rotor sees larger fluctuations in its lift, with the minimum lift occuring when the rotor's leading edge encounters the wake of the upstream stator. The rotor passes through two stator wakes so two minimums are apparent during its rotation over a periodic pitch. Figure 11 shows the Mach contours for the multipassage solution on the SSME. The absolute flow Mach number varies from .148 at the stator inlet through about .472 between the blade rows to .201 at the rotor exit. The inlet absolute flow angle to the stator is 0.0 degrees and the relative inlet flow angle to the rotor is 31.0 degrees. The asymmetry of the flow is apparent in the rotor passages. Here the stator wakes can be seen at different locations in the flow domain. Entropy contours in Fig. 12 show the migration of the stator wakes toward the rotor suction surface as they convect downstream. Concluding Remarks The quasi-threedimensional Euler or the thin-layer, Navier-Stokes equations are solved for unsteady turbomachinery flows. These equations are written in general coordinates for an axisymmetric stream surface, and account for the effects of blade row rotation, radius change and stream-surface thickness. A four-stage Runge-Kutta scheme based on the work of Jameson is used to predict time-accurate results. Implicit residual smoothing is used to increase the stability limit of the Runge-Kutta scheme and thus allow a larger time step to be used in time-accurate flows. A nonconservative interface formulation and other data management techniques allow the solution of rotor-stator interaction problems in both single and multi-passage machines. Viscous results on the supersonic throughflow fan and Space Shuttle Main Engine fuel turbopump have shown the application of implicit residual smoothing to be viable. The results that have been computed indicate that residual smoothing greatly enhance the Runge-Kutta scheme and does not affect the time or space accuracy of the solution. Acknowledgement
doi:10.2514/6.1989-205 fatcat:yslirx42ibhtlav3yeosikvk34