A New Coupled Free Wake-CFD Method for Calculation of Helicopter Rotor Flow-Field in Hover

ABSTRACT: This paper concerns a new coupled free wake-CFD method for proper calculation of aerodynamic loads on a two bladed helicopter rotor in hovering flight. Loading is computed by solving the three dimensional Euler equations in a rotating coordinate system. However, since direct simulation of the tip vortices and wake requires a very fine grid, the rotor’s wake effects are modeled by a free wake approach and included into the CFD calculation by a transpiration boundary condition at the rotor surface. An influence coefficient solution method is used to find the rotor’s wake shape, being steady in a rotating frame. Euler equations are also considered in the form of absolute flow variables and solved by a multi grid Jameson’s finite-volume method. The accuracy of the proposed method is illustrated by comparing numerical results to the available experimental results for the pressure distribution on a blade of a model helicopter rotor at different tip Mach numbers.


INTRODUCTION
The flow field around rotors is inherently complex.In the vicinity of the blade tip, it is highly three dimensional and vortical and it is characterized by non-uniform inflow with blade-vortex interactions and transonic conditions.In hover, the complex' three dimensional rotor wake has a dominant influence on the blade loading (Jenny et al., 1968;Clark and Langrebe, 1971;Landgrebe, 1972).
A proper evaluation of the rotor aerodynamic loading has an important role in the analysis and design of helicopter rotor blades.There are two general approaches for simulating the flow field around a rotor and determining its loads.In the first approach, the inviscid or viscous set of governing equations is solved by a proper Computational Fluid Dynamic (CFD) method (Roberts and Murman, 1985;Chen and McCroskey, 1988;Kramer et al., 1988).Obviously, in this approach, the vortex system near the blade tips and the wake region below the rotor disk must be captured accurately for realistic pressure and load estimations.The computational cost is extremely increased for this approach (Roberts and Murman, 1985;Chen and McCroskey, 1988;Kramer et al., 1988;Ramachandran et al., 1989;Lee and Kwon, 2006;Chen et al., 1987).In the second approach, the tip vortices and the helical wake system are not solved by the CFD calculation and instead, the influence of strong tip vortices and the geometry of the wake are included by a vortex-wake model.
Solution scheme which use the later idea are often grouped under methods using wake models encompassing the potential flow (Strawn and Caradonna, 1987;Strawn and Tung, 1986;Chang and Tung, 1985;Egolf and Sparks, 1987), the Euler (Sankar et al., 1986, Agarwal andDeese, 1987;Chang and Tung, 1987) and the Navier-Stockes methods (Agarwal and Deese, 1988;Narramore and Vermeland, 1989;Wake and Sankar, 1989;Srinivasan and MacCroskey, 1988).Several approaches have been applied to model the vortical field of rotors.The most famous of these models is the prescribed wake model, which uses experimental data to drive empirical formulas relating the rotor blade geometric parameters and thrust coefficient (C T ) to the vortex' wake geometry (Kocurek and Tangler, 1976;Landgrebe, 1972;Landgrebe and Cheney, 1972;Landgrebe et al., 1977).These schemes do not correctly model the flow characteristics and their reliance on experimental data causes to not appropriately predict blade loading for a new rotor configuration.Free wake models have been developed in an attempt to overcome these limitations by simulating in details the actual shape and motion of the wake (Crimi, 1965;Landgrebe, 1969;Clark and Leiper, 1970;Sadler, 1971;Johnson, 1980;Johnson, 1981;Summa and Clark, 1979;Summa and Maskew, 1981;Saberi, 1983;Rosen and Graber, 1988).This can be done by finding the geometry of the wake and the strength of tip and sheet vortices, which are modeled by a finite number of vortex line segments.Then, the effect of the vortex-wake is brought into CFD calculation by modifying either the boundary conditions at the blade surface (Sankar et al., 1986, Agarwal andDeese, 1987;Chang and Tung, 1987) or explicitly modifying the blade angle of attack (Wake and Sankar, 1989;Srinivasan and McCroskey, 1988).
There are some advantages for calculating the rotor aerodynamics loads by a coupled wake-CFD method when compared to the direct method.The first benefit is that accurate enough results can be obtained with significant reduction in the computational expense.This is a key factor when a low fidelity approach is required for design or multidisciplinary optimization purposes.The second benefit is to relieve the difficulties arisen in meshing the 3-D domains around the rotor.Since the tip region is not computed directly in the CFD calculation, the 3-D grid can be obtained from an assembly of two-dimensional spanwise sectional grids (Caradonna et al., 1984).
The shape of the rotor wake has an intense effect on the aerodynamic performance of a helicopter rotor in hover flight.The present research explains a new interactive wake-CFD solver for calculating the inviscid transonic flow on a two-bladed rotor in hover flight.In this method, the CFD flow solver is coupled with a free wake model by using the surface transpiration approach, i.e., effects of rotor vortex-wake are included into CFD calculation by using the estimated downwash velocity to modify boundary condition along the blades.
In most of the previously reported works, the induced velocity is computed prior to the CFD analysis by a separate program such as Hover, B-TRIM of McDonnell Douglas Helicopter Company (Agarwal and Deese, 1987), or Host (Benoit et al., 2000).In this work, however, the induced downwash of rotor blades are computed along with the CFD solution and one can interact with it through an iterative procedure.The model is based on a lifting line theory, which provides estimation for induced downwash along the blade.In this method, the blade is substituted by a row of square vortex panels in radial direction and, at some location of blade span, a finite number of spiral vortices are shed downwards the disc rotor.Then, the positions of helices are corrected by induced velocities, which are computed by employing the Biot-Savart law and the use of the Basic Curved Vortex Element (BCVE) or Self-Induction Vortex Element (SIVE) elements (Quackenbush et al., 1989).The procedure continues until the shape of the wake converges to a steady state and therefore, the self-induced velocities become from zero up to a tolerance.
Then, the transonic flow field about the rotor is calculated by the solution of the three-dimensional Euler equations on a multi-block structured grid about the blades.However, as mentioned earlier, the three-dimensional grid is constructed from an assembly of 2D grids about the blade sections which is extended only to the blade tip and consists of two blocks that are symmetric in relation to the planform.In this study the finite-volume scheme of Jameson, Schmidt, and Turkel (Jameson et al., 1981) is used to discrete governing equations.The method is reformulated for a rotating coordinate system which is attached to the rotor blades, but equations are finally rewritten in terms of absolute-flow variables.This approach has been used previously by Holmes and Tong (1985) in order to develop an Euler solver for some applications in turbomachinery.
The iterative procedure starts by determining the spanwise circulation and the shape of the wake from the blade geometry.Then the flowfield around the blade is computed by the CFD solver, using the calculated downwash velocity along the blade as a boundary condition.Now, a new bound circulation can be obtained from the lift of each section and therefore, it is possible to improve the wake model accordingly.The procedure is repeated until the blade bound circulation or downwash velocity remains unchanged in two subsequent iterations.
The paper is organized as follows: the methodology used for modeling the rotor wake and tip vortices is explained.Governing equations for CFD calculations are presented and the numerical procedure for discretizing these equations is described.The necessary boundary conditions and the way through which the rotor wake model interacts with the CFD solver are explained.Finally, some illustrative numerical results are presented.Calculations are performed for the flow field of a model helicopter rotor in hover at various collective pitch angles.Comparisons with the experimental data of Caradonna and Tung (1981) are also presented.

ROTOR WAKE MODELING
The rotor wake model is based on the proposition that the free wake hover problem possesses a self-preserving state solution when viewed in a rotating frame that is fixed to blades (Bliss et al., 1985).The fundamental idea in the model is that at each point on the rotor wake, when viewed in coordinates rotating with the blade, the velocity contributions due to all effects must be such that there is no convection to change the wake shape.Any velocity component in a plane, normal to the vortex filament curve, will convect the vortex to a new position, thereby changing the filament shape (Quackenbush et al., 1989).
Applying the idea, two features of the wake are properly determined: the contraction of the helical wake as it extends in the axial direction below the rotor, and the interaction of the tip vortex of two blades.For modeling the wake, the two curved vortex elements introduced in Bliss et al. (1985) and Baskin et al. (1967) are used.These elements include the BCVE and SIVE.The BCVE is used in order to compute the velocity induced at any point in the flow field, not on the element itself, and the SIVE is used to evaluate the velocity induced by a vortex filament on itself, including the effect of distributed vorticity in the bent vortex core.Figure 1 shows two types of vortex elements on a typical helical vortex filament.The BCVE is a parabolic arc element positioned between control points whereas the SIVE is a circular arc element passing through each three adjacent control points.
The solution procedure involves a way of reducing the cross flow velocity to zero at every collocation point.Figure 2 shows the idea.In order to achieve this goal, for a given arrangement of the free wake collocation points, the net velocity in each cross flow plane is computed first.Next, the effect of small displacements of each collocation point on the velocity at every other collocation point, and the displaced point itself, are determined.These effects are known as influence coefficient.
In practice, the solution must be obtained numerically with the wake modeled by using discrete elements.The blade is represented by a set of vortex quadrilaterals in a way that . BCVE and SIVE elements for discretizing wake vortex filaments.

Cross-ow Planes Displaced Plane
Initial Plane Method of displacing the vortex filament to eliminate cross flow.
their forward edges are located at the quarter chord line of the blade and the distance of their aft edges from the blade trailing edge is also one-fourth of the local chord.The vortex filaments of sheet vortex and the tip vortex are shed from the points on the aft legs of these panels, and the panels including the shed points are in the form of truncated vortices, as shown in Fig. 3.
Figure 4 shows a schematic of a typical load distribution on a rotor blade.It is assumed that the corresponding wake consists of four filaments, whose strengths are defined as shown in the figure.The locations of the release point are also computed from the Batchelor principle (Batchelor, 1967) as: (1) where is the distance of the release point of the vortex filament from the center of rotation.
The velocity contribution of each individual vortex element is determined analytically using the Biot-Savart law.The total velocity vector at a collocation point includes the local selfinduction effect (SIVE), the contributions from the rest of the wake (BCVE's), the velocity induced by the bound circulation on the rotor blades, and the swirl velocity component associated with the rotating coordinate system.
The influence matrix relates cross flow velocity perturbations to small displacements for a given wake configuration.This matrix can be used to predict the set of collocation point displacements, which will eliminate the cross flow velocity.If the free wake location is described in terms of the positions of N collocation points, and the normal and binormal cross flow velocity components at the ith point are denoted by and , respectively; then the N-dimensional vectors based on these velocity components are defined as: (2) The free wake solution is achieved when a set of collocation point locations is found, such as .The amplitudes of small displacements on the ith collocation point in the local normal and binormal directions are denoted by and .From these, changes of displacement vectors are defined as: (3) and the normal and binormal velocity component changes, associated with the above displacement vectors, are:  The velocity change vectors can be related to the change of displacement vectors through the influence coefficient sub matrices, namely: (5) For instance, the sub matrix gives the changes in normal velocity due to binormal displacements.These matrices must be evaluated numerically.To obtain zero net cross flow velocity, the collocation point displacements must be chosen so that the velocity changes cancel the net cross flow velocity, namely, = and = , which means solving the following matrix equation for and : (6) Since the wake problem is actually nonlinear, this process must be repeated several times.
A free wake analysis involves only a finite number of turns of free vortex.Beyond the free wake region is an adaptive mid wake which completes the contraction process.The final part is also a semi-infinite helical wake having constant radius and constant pitch based on the momentum theory.The location and geometry of the adaptive mid wake depend on the position of the last points in the near section of the wake.In other words, the last points of near wake section are the beginning points of the mid-section, and after the mid-section the wake model is continued with the far section as shown in Fig. 5.
From the conservation of mass flow within the contracting mid wake and by applying the momentum theory, it is possible to drive the boundary shape of the mid wake (Bliss et al., 1985) as: (7) where and are the axial and radial locations of the last free wake points, and is the final contraction radius.In hover flight, the radial contraction exponent is given by: (8) where is the azimuth location of the last free wake point.
When calculating two-bladed rotors, the solution is symmetric from one blade to the next blade.This symmetry allows the induced effect of the second blade on the primary wake to be expressed in terms of the corresponding effect of the primary wake on the second blade.
Once all the different effects have been computed and summed, the resultant velocity vectors at each collocation point are expressed in the blade fixed reference frame.These resultant velocity vectors are solved into the local normal and binormal directions in order to obtain the cross flow velocity components and their changes, associated with point displacements.The velocity changes are normalized by the unit displacement magnitudes to obtain the influence coefficients.These coefficients, as described earlier, are used to predict how the collocation points should be displaced in order to reduce the velocity components in the cross flow planes.The obtained result at this stage is not the final solution due to the nonlinear nature of the problem.Therefore, starting from an initial wake shape, the solution procedure is repeated with a proper relaxation factor in quasi linear steps to improve the convergence rate (Bliss et al., 1984).
The Initial wake shape can be estimated from a contracted wake model with a contraction ratio of = 0.85, where is the radius of the contracted far section of the wake.The circulation distribution is also obtained from an approach explained by Baskin et al. (1967); in other words, equating the 2-D lift coefficient of blade section to the lift coefficient obtained from the Kutta-Joukowski theorem.As a result, a system of linear algebraic equation is obtained by defining the sought circulation.By knowing the circulation, the convergent shape of the wake is computed as explained in the previous section.Finally, at the end, the distribution of the downwash on blade surface is computed.This downwash is then fed to the CFD solver to include the compressibility effect and to correct the circulation.The circulation at this stage is computed from the integration of the pressure field around the blade sections and by using the Kutta-Joukowski theorem.With this new distribution of the circulation, the free wake solver provides a new downwash distribution and this process is iterated until a convergent solution for the downwash is attained.

GOVERNING EQUATIONS
Let (u,v,w) denote the absolute velocity components in the rotating Cartesian coordinate system (x,y,z), as shown in Fig. 6.Compressible Euler equations can be formulated as (Agarwal and Deese, 1987): (9) where: (10) and ( , , ) denotes the rotational velocity components: (11) Defining the flux vector and rotational speed vector as: and: where ( , , ) are the unit vectors in the (x, y, z) coordinate system, Eq. ( 9) can be written as: (12) Equation ( 12) can be solved by the finite-volume method described in the next section.

NUMERICAL METHOD
In order to apply the finite-volume method, Eq. ( 12) is written in its integral form as: The computational domain is meshed by a multiblock structured hexahedral grid.Applying Eq. ( 13) to a typical (i,j,k) cell, for a cell center algorithm, we obtain a system of ordinary differential equations as: (14) where is the cell volume, represents the net convectional flux out of the cell, and and are rotational fluxes out of each cell.We rewrite Eq. ( 14) as:  The convective fluxes can be determined by a second order central scheme.However, the resultant scheme is not dissipative, and therefore, undamped oscillations at odd and even mesh points can be developed during the computation.In order to suppress the tendency for odd-even decoupling, and to prevent the appearance of oscillations in regions containing severe pressure gradients near shock waves and stagnation points, the finite-volume scheme is modified by the addition of artificial dissipative terms as below: Typical values of the constants κ (2) and κ (4) are κ (2) = 0.5, κ (4) = 1/64.Also, is defined as: (23 where CFL is the courant number and ( , , ) are the spectral radii for each cell: (24) a is the speed of sound and (U, V, W) are the contavariant velocities: (25) (U r , V r , W r ) are the components of relative velocity obtained from: (26) ( , , ),( , , ) and ( , , ) are the metrics of transformation for computational coordinates ( , , ).Other dissipation terms, D y W y,j,k and D z W y,j,k , in Eq. ( 17), are calculated in an analogous manner.Equation ( 16) can be integrated in time by a multistage Runge-Kutta scheme.The conserved variables are updated to the new time level n+1 as: (27) where α m = 1, and D(w (k) )=β (k) D(W (n+1,k) )+(1-β (k) )D(W (k-1) ).For a fifth-stage scheme the coefficients are given as: (28) In this work a local time stepping and residual smoothing procedure is also used in order to improve the convergence rate.

BOUNDARY CONDITIONS
Special care must be given to the boundary conditions when solving Euler equations for rotor blades of a helicopter by a coupled free wake-CFD method.The computational domain, as shown in Fig. 7, is a cylinder with radius of 1.5 times the blade span.A far-field boundary condition is applied at the outer surface of the cylinder and also at the blade tip plane.In addition to that, at the root of the blade, the symmetry boundary condition is applied where z axis is along the blade span.At the solid surface of the blade, the transpiration condition is applied to take into account the effects of the wake in the CFD calculations.In this method, the usual velocity to the blade surface is calculated so that it exactly cancels out the effect of the wake normal downwash to the surface.
The boundary condition at the blade surface must be modified in order to include the wake effects.For invicid calculations, this becomes: (29) where is the wake induced velocity and is the normal unit vector of the blade surface.Thus, the induced velocity is used to modify the conventional tangency of velocity vector at the solid surface.Also, the fluxes of density and the enthalpy at the solid surface are modified accordingly to meet the new boundary condition.
It is interesting to note that, for a computational region covering only a portion of the rotor disk, the effect caused by the other rotor blades and the flow outside the finite computational region are incorporated into the calculation by employing the above transpiration velocity technique if the effects of other blades are included in the induced downwash velocity.
The treatment of the far-field boundary condition is based on the concept of Riemann invariants for a one-dimensional flow normal to the boundary (Jameson and Baker, 1983).Let subscripts ∞ and e denote far-field values and values extrapolated from the interior cells adjacent to the boundary, respectively, and let V n and a be the velocity component normal to boundary and to the speed of sound (a).If the flow in the far field is subsonic, fixed and extrapolated Riemann invariants are introduced as: (30) Which are corresponding to incoming and outgoing waves.These invariants may be added and subtracted to obtain: The procedure is completed by extrapolating the entropy at an outflow boundary S = S e or setting it equal to its free stream value at inflow boundary S = S ∞ .The density, energy and pressure can then be calculated from a and s.

RESULTS
In this section some illustrative examples are presented as to demonstrate the capability of the proposed method for predicting the aerodynamic load and the flow field of a model helicopter rotor in hover state at various tip speeds.The rotor is a two bladed rotor with an aspect ratio of 6, which are untwisted and have rectangular planform with the NACA-0012 airfoil section.The experiments on this rotor have been presented by Caradonna and Tung (1981) at the NASA Ames Research Center.Computations are performed on a two-block mesh with 141570 nodes, on a personal computer.Figure 7 shows the computational domain which is symmetric in relation to the blade planform.
Five orders of magnitude reduction in the density residual are considered as convergence criteria.However, in most simulations for a Courant number of 2.0, the convergence criteria is not met in less than 1,600 time steps for the first iteration of the free wake model, and it is reduced to around 800 steps for the next iterations.Some of the results are discussed in the following sub-sections.

NON-LIFTING CASE
For this case, the free wake model is switched off and only the CFD calculation is performed.Pressure distribution on the blade upper and lower surfaces is shown in Fig. 8 for a tip Mach number of M t = 0.52 and zero collective angle.Except near the tip, the agreement between the Euler calculation and the experimental data is obvious.The discrepancy near the tip can be caused by boundary conditions.However, the difference between numerical and excremental results is not significant.

LIFTING CASES
As stated before, for this case, the rotor downwash is first calculated by modeling the blade vortices with the BCVE vortex elements.Then, the effect of the wake and tip vortices is fed into the CFD solver and a new bound circulation is calculated from the sectional lift.Finally, a new iteration is started by correcting the shape of the wake and its strength according to the new bound circulation.The procedure is continued until the induced downwash velocity converges to a steady form.For the lifting cases, calculations are performed for five blade tip Mach numbers and in four different collective angles.For these simulations, the bound vortex is modeled by 32 vortex rings.In the near field part, the tip vortex is modeled by 3 turns of a helical vortex filament, each turn composed by six BCVE elements.The pitch of helix is initially set to in this region, where the initial C T is estimated from a collective angle by using the rotor momentum theory.The sheet vortices are also modeled by 4 filaments, each helical filament has 2 turns and each turn has six basic vortex elements.The mid part of the wake consists of 8 turns of helical filaments with the same axial pitch as the near wake filaments, i.e., . Also the far part of the wake consists of 50 turns with an axial pitch, which is twice the value used in the mid part.
For these simulations the number of iterations required for the convergence of induced downwash velocity is less than 300 iterations for a known circulation distribution in the free wake solver, and the number of iterations for achieving a steady distribution of downwash velocity from the combined free wake-CFD solver at low collective angles is six while it is reduced to less than 4 iterations for higher angles.
The sectional lift coefficient and the circulation distribution obtained from these calculations are shown in Fig. 9.The nondimensional circulation (experimental or numerical) is defined as: where C 1 is the section lift coefficient at z radial station.The circulation shown in these figures is non-dimensioned byR 2 Ω, where Ω is the rotational velocity of the rotor and R is the rotor radius.As can be seen from these figures, the agreement with the experimental data is satisfactory.The thrust and inviscid power coefficients have also been computed and presented in Table 1.Comparison reveals that coefficients have good coincidence with experimental data (Caradonna and Tung, 1981) for all collective angles.
As a typical result, non-dimensional downwash and circulation distribution at intermediate iterations between the free wake model and the CFD solver is depicted in Fig. 10.It should be pointed out that through this iterative process, the compressibility effects in the tip region of the blades are included into the free wake model.Figure 10 shows that the downwash, or circulation distribution, is converged to a steady form in a few iterations.The flow solver can also be run in a full multi-grid (FMG) mode, which accelerates the process of convergence.Figure 16 shows the results obtained by applying the FMG approach in three levels and in a two-block cylindrical mesh.The time of computation can be decreased one order of magnitude and approximately two hours on a personal computer.

CONCLUSION
A new coupled free wake-CFD methodology for evaluation of aerodynamic loads on a helicopter rotor in hovering condition was introduced.In this approach, the Euler equations are solved in a rotating coordinate system and the rotor wake effects are modeled by a free wake approach and included into the CFD solver by a transpiration boundary condition at the blade surface.
The obtained results from this approach for the pressure distribution on the blade surface are promising; the performance parameters such as thrust and induced torque coefficients are particularly estimated with a rather good accuracy.The main advantage of this method is that satisfactory results can be achieved by minimum computational efforts, particularly for the grid generation problem and the required CPU time.This feature becomes important where a low fidelity accurate rotor blade flow solver is required for primarily design or optimization purposes.However, the proposed method ignores important phenomena such as the Blade-Vortex Interaction (BVI) or shock vortex interactions.

Figure 3 .
Figure 3. Schematic of the quadrilateral vortex panels of the blade.

Figure 4 .
Figure 4. Defining the strength of trailing vortex filaments from the bound circulation.

Figure 5 .
Figure 5. Three part wake model for the hover analysis.
fluxes.Jameson et al. (1981) have established that an effective form of dissipative terms for flows with discontinuities is a blend of second and fourth differences, with coefficients which depend on the local pressure gradient.Therefore, dissipative terms are written as follows:

Figure 7 .
Figure 7. Computational region is covered by two blocks of structured grids.

Figure 9 .Figure 10 .
Figure 9. Lift coefficient and circulation distribution at various tip Mach numbers and collective angles.

Table 1 .
The computed thrust and power coefficients and comparison with experimental data.