Real-time solution of nonlinear potential flow equations for lifting rotors

Jianzhe HUANG,David PETERS

aDepartment of Energy and Power Engineering,Harbin Engineering University,Harbin 150001,China

bDepartment of Mechanical Engineering and Materials Science,Washington University in St.Louis,St.Louis,MO 63130,USA

Real-time solution of nonlinear potential flow equations for lifting rotors

Jianzhe HUANGa,*,David PETERSb

aDepartment of Energy and Power Engineering,Harbin Engineering University,Harbin 150001,China

bDepartment of Mechanical Engineering and Materials Science,Washington University in St.Louis,St.Louis,MO 63130,USA

Available online 15 February 2017

*Corresponding author.

E-mail address:jianzhe@harbeu.edu.cn(J.HUANG).

Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

http://dx.doi.org/10.1016/j.cja.2017.02.007

1000-9361©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Analysis of rotorcraft dynamics requires solution of the rotor induced flow field.Often,the appropriate model to be used for induced flow is nonlinear potential flow theory(which is the basis of vortex-lattice methods).These nonlinear potential flow equations sometimes must be solved in real time––such as for real-time flight simulation,when observers are needed for controllers,or in preliminary design computations.In this paper,the major effects of nonlinearities on induced flow are studied for lifting rotors in low-speed flight and hover.The approach is to use a nonlinear statespace model of the induced flow based on a Galerkin treatment of the potential flow equations.

©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Mass-flow;

Nonlinear;

Potential flow;

Real-time flight simulator;

Rotorcraft dynamics;

Swirl;

Wake contraction

1.Introduction

Computation of the induced flow of a lifting rotor for use in flight simulation or preliminary design work requires an efficient treatment of the potential flow equations of a lifting disk with moving pressure singularities(i.e.,the blades).Because a flight simulation must be able to compute flow under all flight conditions,the model must be inherently nonlinear in order to compute flow down to the case of zero flight speed(hover),which is completely nonlinear.In that condition,the momentum terms in the potential flow equations are quadratic and dominate the flow computation.Historically,it has been found that finite-state representations of the flow field are the most advantageous for these real-time flight simulations because they allow easy assembly into the structural dynamics and allow a linearized state-space model from which control system laws may be derived.

The first widely accepted finite-state model was the Pitt-Peters model.1That model represented the flow by three nonlinear states:a uniform flow,a side-to-side gradient,and a fore-to-aft gradient.The nonlinearity,implicit in the equations,meant that the equations gave accurate answers even down under the hover condition,and this was verified by agreement with experimental data.2Within seven years after its inception,this model was utilized in virtually every industrial,government and university code for rotorcraft analysis.3However,this model could only give the normal component of flow at the rotor disk(in terms of these three variables)due to thrust,roll moment and pitch moment.

The next major advancement in finite-state inflow models was an extension to a completely general expansion of the normal flow at the rotor disk due to any arbitrary loading distri-bution.4The model was verified by extensive comparison with laser-Doppler wind-tunnel data taken by the U.S.army.5The method was then expanded to give the correct nonlinear behavior(even in hover)and compared with unsteady data derived from the Georgia Tech hover test facility.6This model also found its way to most industrial and government codes.7,8The Pitt-Peters model is a special case of this newer Peters-He model,and thus the integration was not difficult.

The major drawback of these finite-state codes is that they only give the normal component of flow on the disk.The next advancement was to expand them to include all three components of flow everywhere in the hemisphere above and include the plane of the rotor disk.9,10This model was also extended to include all nonlinearities and then applied to the highly nonlinear condition of the rotor in ground effect11;it contained the Peters-He model as a special case(making it easy to convert codes to the new methodology).However,the Morillo model was ill-conditioned(as opposed to Peters-He which had been robust)due to the lack of certain singular states.Hsieh12was able to find and include the singular states,and Peters et al.were then able to include these states in the nonlinear version of the model.13,14However,the method could not be used in production codes because it did not give the flow below the disk.The flow below the rotor disk is very problematical,especially the flow within the wake(which is not,strictly speaking,potential flow).

The next major advancement for the model was the theoretical development that would allow the computation of the flow below the disk and in the wake by a finite-state model that could be efficient enough for real-time analysis.This advancement was accomplished by the proof of an adjoint theorem that showed that,by adding inexpensive linear adjoint states,the flow below the rotor could be determined.15,16This methodology was then applied to Sikorsky co-axial test data by Nowak et al.17,18The use of co-axial data in hover and forward flight implied that the model must be highly nonlinear and that one had to compute the flow from the lower rotor considering the upper rotor and the flow of the upper rotor considering the lower rotor.Thus,the results of these references validated the model.However,those references showed that a large number of expansion states were necessary to converge,and convergence was seen to deteriorate as one approached edgewise flow,so that the method still lacked the acuity needed to be used in production codes under all flight conditions.

The problem of convergence was finally solved through a series of advancements as follows,which led to a useful model.First,a state-space transformation was introduced(along with blending functions),which could combine the robustness of the Peters-He model with the off-disk accuracy of the Morillo-Duffy model.19,20Next,the Adjoint Theorem was extended so that the flow downstream could be found in perfectly edgewise flow where ‘downstream”can actually be in either the upper or lower hemisphere.21Then,a new set of blending functions had to be determined to find the blending of the new adjoint velocities with existing velocities.22

The final piece of the puzzle is the introduction of all significant nonlinearities into the new formulation(i.e.,the formulation with two Adjoint Theorems and two sets of blending functions),so that this methodology can be used in practical codes.Previous papers have proved the method in highspeed flight.It is the purpose of this paper to introduce those nonlinearities into the model,which are important at low speeds.These include the nonlinearity of wake contraction,the nonlinearity of wake curvature,and the nonlinearity of wake distortion.The work is based on the results in the Ph.D.thesis of the first author.23

2.Nonlinearities

From the basic non-dimensional mass and momentum in three dimensions around the actuator disk,the potential flow equations for incompressible flow are given in Eqs.(1)and(2).

wherevis the induced velocity overV∞;pthe pressure over;τ the reduced time.V∞the free-stream velocity;ξ the non-dimensional coordinate along the free-stream line,positive upstream(Fig.1).

To model the flow fi eld of the rotor disk,the ellipsoidal coordinate system(ν,η,ψ)will be adopted.Such coordinate system can be transformed to the traditional rectangular coordinate system as follows:

Therefore,on the rotor disk,one has η=0.Due to the discontinuity in pressure on the rotor disk,the solution to Laplace’s equations can be obtained in Eq.(6)using ellipsoidal coordinates.

The velocity potentials are defined as the integrals along a streamline of the pressure potentials,

Then the induced velocity in terms of velocity potentials is given in Eq.(9).

Substitution of Eqs.(7)and(9)into Eq.(2),followed by application of the divergence theorem,transforms the Galerkin equations into a set of ordinary differential equations for the components of the velocity potential,which is given in Eq.(10).

whereMis the apparent mass matrix;Lthe influence coefficient matrix;Dthe damping matrix.24The apparent mass matrixMcan be further expressed as

where the letter ‘o’represent odd and ‘e’denotes even.For example,forMoo,the first‘o’meansr+j=odd and the second ‘o’meansm+n=odd.

The closed-form solution for all components of the velocity field above the rotor disk can be obtained from the foregoing differential equations.For the axial velocity in the upper hemisphere,the solution is given in Eq.(11).

The solution in Eq.(11)is terms of Morillo-Duffy variables,and thus it is called as Morillo-Duffy solution for axial velocity which is denoted by

Since the Morillo-Duffy solution does not converge very well on the rotor disk,one can transfer the Morillo-Duffy variables into Nowak-He variablesas

where the subscript ‘odd”means the set of variables withm+n=odd.

The Nowak-He solution then can be expressed as Eq.(15).

For the velocity below the rotor disk,the adjoint theorem can be adopted,as given in Eq.(16)wheresis the distance from the point,at which one desires to compute,along the streamline to the rotor disk atr0and ψ0.For instance,in order to compute the velocity below the rotor disk at pointa(Fig.2),one mustfirst obtain the velocity with time delay at pointb,then add the adjoint velocitywith time delay at pointc,and finally subtract the adjoint velocity at pointd.16

For a helicopter in hover,V∞goes to zero,and the induced velocity computed from linearized potential flow theory will become infinity as the time increases.The reason for this singularity is that the mass flow parameter in the linearized model is based purely on the free-stream velocity.According to nonlinear potential flow theory(as well as momentum theory),the mass flow must also include the induced flow through the rotor.This can be accomplished by attention to momentum theory to give the correct mass flow parameter.When this is done,then the finite-state model in steady,axial flow agrees exactly with classical momentum theory.Therefore,to make the theory compatible with nonlinear momentum,it is necessary to replaceV∞with the total average flow at the rotor.

where χ is the wake skew angle andis the average velocity at the rotor.

The mass-flow parameter is defined by the perturbation quantity:

Eq.(10)with Nowak-He variables is modified by insertion ofVon the left-hand side of the equations.

Eqs.(19)–(22)and Eq.(13)give a nonlinear version of the model that reduces to momentum theory in hover.

Based on the continuity equation,the rotor wake contracts downstream as the pressure expands to atmospheric pressure and the induced velocity increases.The increase in induced flow for the nonlinear model—shown in Eq.(19)—is determined from the Adjoint Theorem for the downstream velocity.What remains,then,is to write a mapping from the linear velocity space of the original model to a warped space that will produce the wake contraction necessary to ensure continuity.The velocity field for the non-deformed geometry is shown in Eq.(23),and the geometric basis for this contraction mapping is illustrated in Fig.3.

From the continuity equation,one obtains

Then the radius of the rotor wake downstream can be found in Eq.(26).

The contraction factor at that plane is defined as

For the hover caseV∞=0,Eq.(27)will be reduced to Eq.(28).

When the contraction factorKis introduced,it gives

3.Swirl velocity

For thexcomponent of the induced velocity with Morillo-Duffy variables on the disk,it can be obtained from Eq.(9)as

When we equate Eqs.(30)and(34)and do some math manipulations,the transformation can be derived as

In a similar manner,theycomponent of the induced velocity based on the present variables(change of variables)can be obtained.Herein,only thexcomponent of the induced velocity will be discussed.In order to improve both the efficiency and accuracy of the solution,only the downstream solution(x＜0)utilizes the special variables developed here.Thexcomponent of the induced velocity in the upper hemisphere(z≤0)is based on our new variables as follows:

For the lower hemisphere(z＞0),similarly to Eq.(16),the adjoint theorem gives

4.Numerical calculations

For the nonlinear model with wake contraction,the Nowak-He variables are used,and only the constant loading case is illustrated.The results have 4 harmonics for the odd terms and no harmonics for the even terms.Thezcoordinate is the axial axis of the rotor disk,andz＞0 represent the location below the disk.The average velocities on the rotor disk,below and above the rotor disk for different free stream velocities are tabulated in Table 1.The contraction factorKcan also be computed by use of Eq.(23).As shown in Table 1,the contraction ratio approaches to unity as the free-stream velocity increases.The wake contracts below the disk,while it expands above the disk.For the free stream velocity is zero,the wake curvature is the largest.ForV∞=0,the real part of axial velocityvzfory=0,withfor ω =0 and χ =0oatz=0,0.4 and-0.4 is shown in Fig.4.The result with wake contraction is represented with black solid line,and the velocity without wake contraction,which is denoted by the green dashed curve,is plotted in the same graph for comparison.For the velocity at the center of the disk(z=0),the two curves coincide.It should be noted that it is impossible to compare these results of Fig.4 with the linear theory because the linear theory gives infinite velocity for this case.The velocity here is completely dominated by the nonlinearity.

To better illustrate the rotor wake shape at the different free stream velocities,thewake contraction streamlinesforV∞=0,0.1 and 0.4 are presented in Fig.5 for the same parameters as in Fig.4.ForV∞=0.4,the wake contraction streamlines are almost flat,and little wake contraction phenomenon can be observed.

Since thexandycomponent of the induced velocity are computed in a similar manner,only thexcomponent of the induced velocity will be illustrated herein.For constant loading,the collective pitch is held steady,while for dynamic loading,the collective pitch oscillates with a constant reduced frequency of ωV∞/RwhereV∞is the free-stream velocity;Rthe radius of the rotor;ω a non-dimensional parameter,and the frequency can be varied by changing ω.For either constant loading or dynamic loading with ω=1,the results are obtained with 12 harmonics for the odd terms and 8 harmonics for the even terms(74 states).The criteria for determination of the number of odd terms or even terms are based on the numerical tests according to the combination which gives the best convergence.In the results to follow,we assume lightly-loaded rotors,so that there is no contraction nonlinearity.These results demonstrate the importance of the newlydefined state variables and the adjoint theorem on the results.

Table 1 Average induced velocity and contraction factor.(y=0,-1 ＜ x ＜ 1 with for ω =0,χ =0°(m-odd=4,m-even=0)).

Table 1 Average induced velocity and contraction factor.(y=0,-1 ＜ x ＜ 1 with for ω =0,χ =0°(m-odd=4,m-even=0)).

z Velocity and contraction factor V∞0 0.1 0.2 0.3 0.4 0¯v0 0.0581 0.0286 0.0172 0.0120 0.0092 0.4¯v10.08630.04250.02560.01790.0137 K 0.8205 0.9500 0.9812 0.9907 0.9945 1.0¯v10.10200.05020.03020.02120.0162 K 0.7547 0.9253 0.9714 0.9856 0.9916-0.2¯v10.04030.01980.01190.00830.0064 K 1.2007 1.0361 1.0124 1.0060 1.0034-0.4¯v10.02990.01470.00890.00620.0047 K 1.3940 1.0589 1.0197 1.0094 1.0055

The pressure distribution on the rotor disk is chosen to be elliptical(i.e.,only τ01is used).In Figs.6 and 7,the results for ω=0(steady loading)are illustrated.The exact solution is shown as open circles.The exact solution is found from a convolution integral from upstream infinity along a streamline down to that particular point.Such exact solutions have been correlated with the measured data.5The solid curve denotes the solution computed based on Huang-He variables.On the same plot of each figure,the Morrilo-Duffy solution is drawn with the dashed-dot line for comparison.For ω=0,the solutions have no imaginary part.Fig.6 gives thexcomponent of induced velocity at the rotor disk(y=0 andz=0)for different skew angles.For small skew angles,the present(Huang-He)solution and Morillo-Duffy solution both perform well compared with the exact solution.When the skew angle grows larger,the Huang-He solution has less oscillation for the ondisk region(-1≤x≤1)and it is also slightly more accurate for computing the flow field upstream(x＞1),while Morillo-Duffy solution has better convergence downstream(x＜-1).In Fig.7(a),thexvelocity with a distance one radii above the disk(z=-1.0)is plotted.Both Huang-He solution and Morillo-Duffy solution correlate with the exact solution well.Fig.7(b)shows the results with distance 1.0 radii below the disk(z=0.1)versusx0,which is the distance from the center of the skewed wake.Therefore,x0=0 is the streamline going through the rotor center.For on-disk region and upstream,Huang-He solution has less oscillation and matches the exact solution better.For the flow field downstream,both of these solutions have large error.The reason for this downstream error is that there is highly concentrated vorticity at the trailing edge of the disk.This makes the convergence slow there.The results here do not have enough harmonics to converge in this region.For a more detailed description of the convergence properties of the method,see Ref.25.

The results forxcomponent of the induced velocity with dynamic loading(ω=1)are presented in Figs.8–11.The results include both real part(in-phase velocity)and imaginary part(out-of-phase velocity).In Figs.8 and 9,the induced velocities on the rotor disk(y=0 andz=0)for χ =30°and 60°are given,respectively.For in-phase velocity,Huang-He solution and Morillo-Duffy solution have almost the same performance.However,Morillo-Duffy solution gives better convergence downstream for the out-of-phase velocity.Forxcomponent of the induced velocity 1.0 radii above and below the rotor disk,they are illustrated in Figs.10 and 11,respectively.The results do not well correlate with the exact solution especially downstream,but they can predict the trend very well.

5.Other comparisons

The above results compare exact solutions to the potential flow equations with our method. Since other numerical approaches26,27such as the vortex-lattice method is also an approximate solution to the exact potential flow equations,there is no point in comparing our method with vortexlattice results,since we have already compared it with the exact solution.However,one might ask how solutions to the potential flow equations are compared with experimental data.This question is actually answered in Ref.25,which shows comparison with experimental data for a co-axial rotor system.Because there are two rotors,one must compute the flow above the lower rotor and below the upper rotor.Figs.12–14,taken from Ref.25,show the comparison.Those results are for hover and include the wake contraction nonlinearities mentioned herein.The results show good correlation for inflow,trust and power.

6.Conclusions

(1)This paper completes the final piece of development that will allow finite-state inflow models(used in flight simulators and preliminary design codes)to compute all three components of induced flow everywhere in the flow field(including inside of the rotor wake).This final piece of development is the addition of all pertinent nonlinearities that will allow the computation to be valid down to hover,which is a completely nonlinear flight condi tion––including the addition of these nonlinearities in such a way that they do not detract from overall computational performance.

(2)This advancement has been accomplished by several contributions:the nonlinear mass flow parameter is added to the mass-flow matrix for the He-like variables;the wake skew parameter is rendered to be nonlinear;the wake curvature parameters are included in the odd Legendre functions;based on conservation of mass,a nonlinear mapping is introduced,which allows the entire wake(both upstream and downstream)to expand or contract according to the conservation laws.

(3)The swirl velocity can also be calculated based on the method of change of variables,which is another breakthrough in the history of the finite-state inflow model.The Huang-He variables and Huang-He solution are derived to compute the swirl velocity.

(4)Only thexcomponent of the induced velocity has been discussed in this paper,and theycomponent of the induced velocity can be easily obtained in a similar manner.With bothxcomponent andycomponent of the induced velocity,the swirl velocity can be constructed.

(5)Some results have been computed for illustration purpose for constant loading and dynamic loading with elliptical pressure distribution.The newly obtained Huang-He solution gives good enough prediction in most of the cases.

(6)Toimprovetheconvergence,Huang-Hesolution,Morillo-Duffy solution and downstream velocity have to be blended,which will be studied in the future.

Acknowledgements

This study was co-supported by the National Natural Science Foundation of China(No.51375104),the Heilongjiang Province Funds for Distinguished Young Scientists(No.JC 201405),the China Postdoctoral Science Foundation(No.2015M581433),and the Postdoctoral Science Foundation of Heilongjiang Province(No.LBH-Z15038).

1.Pitt DM,Peters DA.Theoretical prediction of dynamic-inflow derivatives.Vertica1981;5(1):21–34.

2.Gaonkar GH,Peters DA.Review of dynamic inflow modeling for rotorcraft flight dynamics.Vertica1988;12(3):213–42.

3.Peters DA,HaQuang N.Dynamic inflow for practical applications.J Am Helicopter Soc1988;33(4):64–8.

4.Peters DA,Boyd DD,He CJ.A finite-state induced-flow model for rotors in hover and forward flight.J Am Helicopter Soc1989;34(4):5–17.

5.Peters DA,He CJ.Correlation of measured induced velocities with a finite-state wake model.J Am Helicopter Soc1991;36(3):59–70.

6.Su A,Yoo KM,Peters DA.Extension and validation of an unsteady wake model for rotors.J Aircraft1992;29(3):374–83.

7.Peters DA,Karunamoorthy S,Cao WM.Finite-state induced flow models,part I:Two-dimensional thin airfoil.J Aircraft1995;32(2):313–22.

8.Peters DA,He CJ.Finite-state induced flow models,part II:Three-dimensional rotor disk.J Aircraft1995;32(2):323–33.

9.Morillo JA,Peters DA.Velocity field above a rotor disk by a new dynamic inflow model.J Aircraft2002;39(5):731–8.

10.Peters DA,Morillo JA,Nelson AM.New developments in dynamic wake modeling for dynamics applications.J Am Helicopter Soc2003;48(2):120–7.

11.Yu K,Peters DA.Nonlinear three-dimensional state-space modeling of ground effect with a dynamic flow field.J Am Helicopter Soc2005;50(3):259–68.

12.Hsieh A.A complete finite-state inflow model for rotors in axial flow[dissertation].St.Louis(MO):Washington University in St.Louis;2006.

13.Peters DA,Hsieh A,Garcia-Duffy C.A complete finite-state inflow theory from the potential flow equations.Proceedings of the 3rd international basic research conference on rotorcraft technology;2009.

14.Garcia-Duffy C,Hsieh A,Peters DA.A complete nonlinear induced flow theory for rotors in incompressible flow.Proceedings of the 11th pan American congress of applied mechanics–PACAM XI.Rio de Janeiro:ABCM.

15.Fei Z,Peters DA.Applications of generalized dynamic wake theory to the flow in a rotor wake.Proceedings of the 69th annual forum of the American helicopter society.Alexandria:The AHS International,Inc.;2013.p.371–8.

16.Fei Z,Peters DA.Accurate computation of flow below the rotor disk by a finite-state method.Proceedings of the 31st AIAA applied aerodynamics conference.Reston:AIAA;2013.

17.Nowak M,Prasad JVR,Xin H,Peters D.A potential flow model for a coaxial rotors in forward flight.Proceedings of the 39th European rotorcraft forum;2013.

18.Nowak M,Fei Z,Peters DA,Prasad JVR.Improved finite-State inflow convergence through use of a blended model.Proceedings of the 5th decennial AHS aeromechanics specialists’conference 2014:current challenges and future directions in rotorcraft aeromechanics.Alexandria:The AHS International,Inc.;2014.p.302–11.

19.Huang J,Nowak M,Peters DA,Prasad JVR.Converged velocity field for rotors by a blended potential flow method.Proceedings of the 70th annual forum of the American Helicopter Society.Alexandria:The AHS International,Inc.;2014.p.44–54.

20.Nowak M,Prasad JVR,Xin H,Peters DA.Development of a finite-state model for a co-axial rotor in forward flight.Proceedings of the 70th annual forum of the American Helicopter Society;Alexandria:The AHS International,Inc.;2014.

21.Peters DA,Huang J.Real-time solutions of nonlinear potential flow equations for lifting rotors.Proceedings of the 5th nonlinear science and complexity conference;2006.

22.Huang JZ,Nowak M,Peters D,Prasad JVR.Converged velocity field for rotors by a blended potential flow method.The Journal of Aviation Technology2014;1(2):1–10.

23.Huang JZ.Potential-flow inflow model including wake distortion and contraction[dissertation].St.Louis(MO):Washington University in St.Louis;2015.

24.Morillo JA.A fully three-dimensional unsteady rotor inflow model from a galerkin approach[dissertation].St.Louis(MO):Washington University in St.Louis;2001.

25.Prasad JVR,Nowak M,Xin H.Finite-state inflow models for a coaxial Rotor in hover.Proceedings of the 38th European rotorcraft forum;2012.p.949–58.

26.Xu HY,Ye ZY.Numerical simulation of unsteady flow around forward flight helicopter with coaxial rotors.Chinese J Aeronaut2011;24(1):1–7.

27.Shi YJ,Xu Y,Xu GH,Wei P.A coupling VWM/CFD/CSD method for rotor airload prediction.Chinese Journal of Aeronautics2017;30(1):204–15.

9 May 2016;revised 19 September 2016;accepted 19 October 2016