Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Development and applications of a body-force propulsor model for high-fidelity CFD
(USC Thesis Other)
Development and applications of a body-force propulsor model for high-fidelity CFD
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DevelopmentandApplicationsofaBody-ForcePropulsorModel forHigh-FidelityCFD by TianboXie ADissertationPresentedtothe FACULTYOFTHEGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (AEROSPACEENGINEERING) August2023 Copyright 2023 TianboXie Acknowledgements I wish to express my profound gratitude to my advisor, Dr. Alejandra Uranga, whose support and guidance throughout my doctoral studies have been invaluable. The successful completion of this PhDwouldnothavebeenpossiblewithoutDr. Uranga’sencouragementandbeliefinmyabilities. I am grateful to Arturo Cajal, my officemate for the past six years, and Michael Kruger, who became my officemate for the last two years, and all members of the USC ADRL group. Their emotional and intellectual support played a significant role as I navigated the challenging journey of my PhD project. My deepest thanks go to my fianc´ ee, Fangning Zheng. Her companionship andsupporthavebeenmystrengthinthefaceofthemoststressfultimes. Lastbutnotleast,Iwish toexpressmygratitudetomyparents,YongyiXieandXiaoweiZhao. Theirsteadfastsupportover the last decade, as I pursued my dreams far away from home, has been the bedrock on which my successisbuilt. ii Contents Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx Chapter1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Background: PropulsorModelingMethodologies . . . . . . . . . . . . . . . . . . 5 1.2.1 Propeller,HelicopterRotorandWindTurbineModeling . . . . . . . . . . 6 1.2.2 TurbomachineryModeling . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2.1 WorkInputandFlowTurning . . . . . . . . . . . . . . . . . . . 10 1.2.2.2 Through-Flow(TF)Methods . . . . . . . . . . . . . . . . . . . 11 1.2.2.3 Body-ForceMethods . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.3 ProsandConsofModelingMethodologies . . . . . . . . . . . . . . . . . 16 1.3 ThesisScope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Chapter2: TheBody-ForceModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.1 Body-ForceModelingApproach . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 BladeMetalBlockage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 InviscidLoading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 ImprovedAnalyticalFormulation . . . . . . . . . . . . . . . . . . . . . . 30 2.4 ViscousLosses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4.1 GeneralForm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4.2 ReferenceProfileLosses . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4.3 Off-DesignLosses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.4 ShockLosses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Chapter3: Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.1 ComputationalMethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2 ModelImplementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.1 GeometryExtraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.2 DataAssignmenttoComputationalGrid . . . . . . . . . . . . . . . . . . . 45 3.2.3 SourceTermEvaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Chapter4: ModelAssessmentandValidationinUniformFlow . . . . . . . . . . . . . . . 48 iii 4.1 OverviewofTestCases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 TransonicRotor: NASARotor67 . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1 ComputationalSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.2 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.2.3 QualitativeObservations . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.4 FlowTurning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2.5 MetalBlockage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2.6 Losses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2.7 OverallPerformance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2.8 DeterminationofReferenceDeviationAngle . . . . . . . . . . . . . . . . 66 4.3 TransonicStage: NASAStage35 . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3.1 ComputationalSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3.2 OverallPerformance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.4 Low-SpeedStage: TF8000FanStage . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4.1 ComputationalSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4.2 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4.3 EffectofFlowAngleCorrection . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.4 OverallPerformance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Chapter5: CapturingtheEffectofPropulsorInflowDistortion . . . . . . . . . . . . . . . . 80 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2 Rotor67withRadialDistortion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.3 Rotor67withCircumferentialDistortion . . . . . . . . . . . . . . . . . . . . . . . 86 5.3.1 ProblemSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3.2 UpstreamFlowRedistribution . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3.3 DownstreamFlowFieldandDistortionTransfer . . . . . . . . . . . . . . . 89 5.3.4 DistortionEffectonOverallPerformance . . . . . . . . . . . . . . . . . . 93 5.3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.4 TF8000FanStagewithBLIInflowDistortion . . . . . . . . . . . . . . . . . . . . 95 5.4.1 OverviewofBoundaryLayerIngestion(BLI) . . . . . . . . . . . . . . . . 95 5.4.2 ProblemSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.4.3 PropulsorMeridionalFlowField . . . . . . . . . . . . . . . . . . . . . . . 97 5.4.4 PropulsorUpstreamFlowField . . . . . . . . . . . . . . . . . . . . . . . . 99 5.4.5 FanResponseandDistortionTransfer . . . . . . . . . . . . . . . . . . . . 101 5.4.6 DistortionEffectonOverallPerformance . . . . . . . . . . . . . . . . . . 105 5.4.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Chapter6: ExtendedModelforFree-TipPropulsors . . . . . . . . . . . . . . . . . . . . . 108 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.2 ModificationstoBody-ForceFormulation . . . . . . . . . . . . . . . . . . . . . . 109 6.2.1 WakeInducedCorrections . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.2.2 TipLossCorrection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.2.3 LoadLimit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.2.4 SummaryofChangestoBody-ForceFormulation . . . . . . . . . . . . . . 113 iv 6.3 PropellersinUniformFlow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.3.1 ComputationalSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.3.2 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3.3 QualitativeObservations . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.3.4 EffectofFlowAngleCorrections. . . . . . . . . . . . . . . . . . . . . . . 118 6.3.5 OverallPerformance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.4 PropellersatNon-ZeroAnglesofAttack . . . . . . . . . . . . . . . . . . . . . . . 121 6.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.4.2 PropellerFlowField . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.4.3 OverallPerformance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Chapter7: CapturingPropeller-WingInteractionsinBlown-LiftSystem . . . . . . . . . . . 127 7.1 OverviewoftheBlown-LiftConcept . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.1.1 Blown-LiftWingGeometry . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.1.2 ActuatorDiskModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.1.3 ComparisonbetweenADandBFMinisolation . . . . . . . . . . . . . . . 134 7.1.4 LocalFlowPropertyDistributions . . . . . . . . . . . . . . . . . . . . . . 134 7.1.5 OverallPerformancePredictedbyBody-ForceModel . . . . . . . . . . . . 137 7.2 Wing-PropellerInteractions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 7.2.1 BasicFeaturesoftheFlowOveraBlown-Wing . . . . . . . . . . . . . . . 139 7.2.2 PropellerInflowNon-UniformityandModelResponse . . . . . . . . . . . 141 7.2.3 CoupledWingandPropellerPerformances . . . . . . . . . . . . . . . . . 144 7.2.4 JetMixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 7.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 Chapter8: Conclusions,FutureWorkandPotentialApplications. . . . . . . . . . . . . . . 152 8.1 SummaryandConclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 8.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 8.3 FutureWorks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 8.3.1 ModelImprovements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 8.3.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 A Derivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 A.1 EnergySourceTerm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 A.2 EntropyGeneration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 A.3 BlockageForceTerm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 B CompleteBody-ForceFormulation . . . . . . . . . . . . . . . . . . . . . . . . . . 167 C BlockageModelValidation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 v ListofTables 4.1 Summaryoftestcasesglobalparameters. . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Resultofmeshconvergencestudy . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 Performance parameters at peak efficiency and near stall operating points: mass flow rate ˙ m, mass-averaged total pressure ratio FPR, and isentropic efficiencyh is . Comparisonbetweenexperimentalmeasurementsandbody-forcepredictions. . . . 63 7.1 Geometricalparametersoftheblown-liftwing. . . . . . . . . . . . . . . . . . . . 129 7.2 Operatingconditionoftheexaminedisolatedpropellercase. . . . . . . . . . . . . 134 7.3 Referenceoperatingconditionsfortheblown-liftcase. Thebody-forcemodeloper- atingconditionissetviathewheelspeedW ,andtheactuatordiskmodelcondition issetviatheisolatedthrustmagnitude. . . . . . . . . . . . . . . . . . . . . . . . 139 vi ListofFigures 1.1 Aircraft designs with boundary layer ingestion (BLI): (left) the D8 double bubble aircraft,and(right)theNASASTARC-ABLconceptualaircraft. . . . . . . . . . . 2 1.2 Aircraft designs with blown-lift: (left) the Electra super eSTOL, and (right) the NASAX-57Maxwell. (ImagebyElectra.aeroandNASA) . . . . . . . . . . . . . 3 1.3 GeneralizedpropulsorrepresentationswiththeNASARotor67geometry. . . . . . 5 1.4 Force,velocityandflowanglesonabladeelement. . . . . . . . . . . . . . . . . . 7 2.1 Bladenumbereffects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Body-forcerepresentationofatransonicfanrotor: (left)actualpropellergeometry, and (right) body-force volume in the computational domain showing camberline normals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3 Illustrationoftheinviscidturningmodel. . . . . . . . . . . . . . . . . . . . . . . . 28 2.4 Demonstration of the effect of no induced angle correction on BFM load predic- tions for propulsors with a smaller number of blades: the same amount of loading correspondstoasmallerflowcoefficientor,inotherwords,under-predictionofthe loadingataspecificflowcoefficient. . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Pictorialdemonstrationoflocalangleofattack. . . . . . . . . . . . . . . . . . . . 32 2.6 Illustrationoftheviscousbody-forcemodel. . . . . . . . . . . . . . . . . . . . . . 35 2.7 Boundary layer dissipation coefficient versus Reynolds number based on momen- tumthicknessforlaminarflowandturbulentflow. . . . . . . . . . . . . . . . . . . 39 3.1 Modelimplementationworkflow . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Pictorialdemonstrationofgeometryextraction . . . . . . . . . . . . . . . . . . . . 44 vii 3.3 Contour plots showing distribution of relevant blade geometrical properties for NASA Rotor 67: Blade metal angle k , blade thickness t, blade blockage factor bandaxialgradientofbladeblockagefactor— b. . . . . . . . . . . . . . . . . . . 45 4.1 Distributionoftestcases: diameterversusnumberofblades. . . . . . . . . . . . . 49 4.2 3DRenderingofNASARotor67. . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 MeridionalmeshNASARotor67. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.4 Meridional contours of flow properties through NASA Rotor 67: (a) Normalized total temperature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number,(e)Absoluteflowangle,(f)Relativeflowangle,(g)Normalizedtotalpres- sure,and(h)Normalizedpressure,atpeakefficientoperatingpointobtainedusing theproposedbody-forcemodel. Thinblacklinesshowwheretheleadingandtrail- ing edges of the blades would be located, and enclose the volume where body- forcesareapplied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.5 Normalizedpressurecontours: qualitativecomparisonbetween(left)ouranalytical formulation,and(right)inviscidflowtangencymethod. [-4ex] . . . . . . . . . . . 57 4.6 Meridional contours of flow properties through NASA Rotor 67: (a) Normalized total temperature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number,(e)Absoluteflowangle,(f)Relativeflowangle,(g)Normalizedtotalpres- sure,and(h)Normalizedpressure,atchokedconditionobtainedusingtheproposed body-force model. Thin black lines show where the leading and trailing edges of thebladeswouldbelocated,andenclosethevolumewherebody-forcesareapplied. 58 4.7 Normalized pressure contours showing the captured shock pattern at choked con- dition: comparisonbetween(left)analyticalformulationand(right)flowtangency method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.8 Spanwise distribution of incidence angle (i), outlet deviation angle (d TE ), camber angle (J ), flow turning angle (e ⌘ b 1 b 2 ), relative flow angle at inlet (b 1 ) and outlet(b 2 ),forNASARotor67atnearpeakefficiency(PE)operatingpoint. . . . . 59 viii 4.9 Weinig’s coefficient versus pitch-to-chord ratio at spanwise location varying from 10%to80%,forNASARotor67withvaryingsolidity(B=2to50),withcompar- isontotheoreticaltrends. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.10 Normalized pressure contours showing the captured shock pattern at choked con- dition: comparisonbetween(left)analyticalformulationand(right)flowtangency method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.11 Total pressure ratio versus mass flow at design rotational speed: comparison be- tweencaseswithandwithoutblockagemodelapplied. . . . . . . . . . . . . . . . 60 4.12 Total pressure ratio and isentropic efficiency versus mass flow: comparison of re- sultswithandwithoutlossmodels. . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.13 Spanwisedistributionofflowpropertiesinpeakefficiency(PE)andnearstall(NS) operation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.14 NASA Rotor 67 performance map, with comparision between prediction from the proposedbody-forcemodel,experimentalmeasurementsandCFDresults. . . . . 67 4.15 Diffusion factor versus normalized mass flow for NASA Rotor 67 as represented viathebody-forcemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.16 3DRendering(left)androtorbladesections(right)ofNASAStage35. . . . . . . . 68 4.17 MeridionalmeshNASAStage35. . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.18 Meridional contours of flow properties through NASA Stage35: (a) Normalized total temperature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number,(e)Absoluteflowangle,(f)Relativeflowangle,(g)Normalizedtotalpres- sure, and (h) Normalized pressure, at design operating point obtained using the proposed body-force model. Thin black lines show where the leading and trailing edges of the blades would be located, and enclose the volume where body-forces areapplied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 ix 4.19 Meridional contours of flow properties through NASA Stage35 after choked: (a) Normalized total temperature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number, (e) Absolute flow angle, (f) Relative flow angle, (g) Nor- malized total pressure, and (h) Normalized pressure, at design operating point ob- tained using the proposed body-force model. Thin black lines show where the leading and trailing edges of the blades would be located, and enclose the volume wherebody-forcesareapplied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.20 NASA S35 overall performance as predicted by the body-force model (solid line) withcomparisonwithexperimentalmeasurements. . . . . . . . . . . . . . . . . . 72 4.21 3D rendering of TF8000 fan stage (left) and rotor blade sections at 10%, 50% and 100%spanwiselocations(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.22 MeridionalmeshTF8000fanstage. . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.23 MeridionalcontoursofflowpropertiesthroughTF8000: (a)Normalizedtotaltem- perature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number, (e)Absoluteflowangle,(f)Relativeflowangle,(g)Normalizedtotalpressure,and (h) Normalized pressure, at design operating point obtained using the proposed body-force model. Thin black lines show where the leading and trailing edges of thebladeswouldbelocated,andenclosethevolumewherebody-forcesareapplied. 76 4.24 Meridional contours of flow properties through TF8000 fan stage: total pressure coefficient (left column) and absolute swirl angle (right column), comparison be- tweencaseswithflowanglecorrection(toprow)andwithoutflowanglecorrection (bottomrow)atf =0.37obtainedusingtheproposedbody-forcemodel. . . . . . 77 4.25 TF8000fanstageperformance: comparisonbetweenexperimentaldatafromSiu[62], body-forcemodelresultswithandwithoutflowanglecorrectionD d ,withandwith- outoff-designlosscomponent,atarotationalspeedof13500rpm. . . . . . . . . . 78 5.1 Illustrationoftheparallelcompressormodelconcept. . . . . . . . . . . . . . . . . 81 5.2 Schematicillustratingupstreamflowredistributioneffects. . . . . . . . . . . . . . 82 x 5.3 NASA Rotor 67 distortion transfer: flow properties at different axial stations in cleanflowandwithinletdistortion. . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.4 NASARotor67performanceincleanflowandwithinletdistortionat90%design wheel speed: effect of distortion and comparison with experimental data found in Fidalgo[5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.5 Inlet total pressure contour (left) and circumferential total pressure profile at mid- span(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.6 Upstream (station 1) static pressure contour (left) and flow coefficient contour (right)producedbythebody-forcemodel. . . . . . . . . . . . . . . . . . . . . . . 87 5.7 Upstream(station1)swirlanglecontour(left)andradialanglecontour(right)pro- ducedbythebody-forcemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.8 Upstream(station1)circumferentialdistributionofflowpropertiesatmid-span,as comparedwithURANSresults. . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.9 Downstream(station2)flowcontoursoftotalpressure(leftcolumn)andtotaltem- perature(rightcolumn)producedbythebody-forcemodel(toprow),ascompared withURANSresults(bottomrow). . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.10 Downstream(station2)circumferentialdistributionofflowpropertiesatmid-span, ascomparedwithURANSresults. . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.11 Contoursofentropygenerationupstream(left)anddownstream(right)oftherotor. 92 5.12 Rotor67overallperformanceunderdistortion. . . . . . . . . . . . . . . . . . . . . 92 5.13 MeridionalviewofcomputationaldomainandgridsforTF8000fanstage. . . . . 96 5.14 Verticaldistortionprofile(left)andinletprofileofadvanceratio(right). . . . . . . 97 xi 5.15 MeridionalcontoursofflowpropertiesthroughTF8000: (a)Normalizedtotalpres- sure, (b) Normalized pressure, (c) Normalized axial velocity, (d) Absolute swirl angle, (e) Normalized total temperature, and (f) Entropy generation, at nominal operating point obtained using the proposed body-force model. Thin black lines show where the leading and trailing edges of the blades would be located, and enclosethevolumewherebody-forcesareapplied. . . . . . . . . . . . . . . . . . 98 5.16 Axial contours of flow properties: total pressure coefficient, pressure coefficient, and absolute swirl angle and radial angle, at upstream of the rotor at station 2, viewingfromdownstream(-xdirection). . . . . . . . . . . . . . . . . . . . . . . . 100 5.17 Axial contours of flow properties: normalized axial velocity (first row), pressure coefficient (second row), absolute swirl angle (third row), and radial angle (fourth row)atthreedifferentaxiallocations(2,3,4and6),viewingfromdownstream(-x direction). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.18 Axialcontoursofflowproperties: totalpressurecoefficient(firstrow),normalized totaltemperature(secondrow),andentropyrise(thirdrow),atthreedifferentaxial locations(2,4and6),viewingfromdownstream(-xdirection). . . . . . . . . . . . 103 5.19 ComparisonofC p t contouratstatorexit(station3)betweenbody-forcemodel(left) andexperiment(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.20 TF8000fanstageperformanceincleanflowandwithinletdistortion. . . . . . . . . 106 6.1 IllustrationofcirculationintegrationpathforobtainingrelationinEquation(6.1). . 110 6.2 Schematicofthecomputationaldomain(left)andmeridionalgridviewzoomedin onthemodelregion(right)fortheAPC8b7Spropellertestcase. . . . . . . . . . . 115 6.3 APC8x7Spropellervelocitycontourwithstreamlines,atJ =0.5,W =6800RPM. . 116 6.4 APC8x7SmeridionalflowcontoursofflowpropertiesatJ =0.5,W =6800RPM. . 117 6.5 Q-criterion isosurfaces from the body-force model, overlaid with normalized ve- locitycontours. Indicatingcylindricalvortexsystematthetipandrootoftherotor. 117 xii 6.6 Thrust coefficient versus advance ratio for APC8X7S propeller, predicted via the body-force model at W = 6800RPM, with comparisons between cases with and withoutwakecorrection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.7 Overall performance (thrust coefficient, power coefficient and efficiency) versus advanceratioforvariousAPCpropellers,predictedbythebody-forcemodel,com- paredwithexperiemnetaldata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.8 APC8x7S meridional flow contours ofC p (first row) andC p t (second row) at J = 0.5,andatangleofattacka p of0,30,50,and70degrees(fromlefttorightineach subplot),wherea p isdefinedasflowanglewithrespecttox-axiswithinx-zplane. 123 6.9 APC8x7SaxialflowcontoursofC p t downstreamofthepropelleratJ =0.5,andat angleofattacka p of0,30,50,and70degrees(fromlefttoright). Propellerrotates clockwiseinthisview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.10 APC8x7Spropellerperformanceat6800rpm,includingthrustcoefficients,power coefficients, and propeller efficiency at increasing angles of attack as predicted by proposedbody-forcemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.11 Trends of c T variation with increasing propeller angle of attack a p as found by experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 7.1 Blown-wing geometry used in wind tunnel tests: 2D airfoil geometry wth flap angled f =20 (left);and3Dgeometryofquarterwingspanwithmotorpylonand propeller(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.2 Spanwisecutviewthroughthecenterofthepropelleroftheoversetcomputational gridsfortheblownflappedwingata =10 (left),andmeridionalcutoftheisolated propellergridwiththepropellerhubatthebottomfigureedge(right). . . . . . . . 131 7.3 Actuatordiskradialloaddistributionfordifferentx o constantswithB=4,l =0.1, andr hub =0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.4 Actuatordiskradialloaddistributionforvariousadvanceratiosl withB=4,x 0 =0, andr hub =0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 xiii 7.5 Flowcontourspredictedbythebody-forcemodelforisolatedpropulsoratl =0.085, W =23000RPM: pressurecoefficient,totalpressurecoefficient,normalizedveloc- ity,swirlangle,andentropy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7.6 Contour of body-force model quantities for isolated propulsor at l =0.085,W = 23000RPM: localmomentumboundarylayerthicknessq ,ReynoldsnumbersRe q andRe x ,dissipationcoefficientC D ,andflowdeviationangled . . . . . . . . . . . . 135 7.7 Flowcontourspredictedbytheactuatordiskmodelforisolatedpropulsoratl =0.085, W =23000RPM: pressurecoefficient,totalpressurecoefficient,normalizedveloc- ity,swirlangle,andentropy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 7.8 Radialprofilesatupstream(station1)anddownstream(station2)locationsforiso- lated propulsor at l =0.085,W =23000RPM: pressure coefficient, total pressure coefficient,swirlangle,andaxialvelocity. . . . . . . . . . . . . . . . . . . . . . . 137 7.9 Overallpropellerperformanceforisolatedpropeller: thrustcoefficient(top),aero- dynamic efficiency (bottom left), and power coefficient (bottom right) versus ad- vanceratiol . Comparisonbetweenbody-forcemodel(BFM),andbothwindtun- nelmeasurements(Exp)andQproppredictions. . . . . . . . . . . . . . . . . . . . 138 7.10 Velocity contours of wing sectional flow at the reference condition with the body- forcemodel(top)andtheactuatordiskmodel(bottom),inazoomed-outview(left) andzoomed-inview(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.11 Vertical profiles for blow-lift case reference condition at upstream (station 1) and downstream (station 2) locations: pressure coefficient, total pressure coefficient, normalizedaxialvelocity,andswirlangle. . . . . . . . . . . . . . . . . . . . . . . 141 7.12 Body-force model solution: upstream (station 1) contours of local advance ra- tio (a), and normalized swirl velocity (b); and downstream (station 2) contours of pressurecoefficient(c),totalpressurecoefficient(d),normalizedaxialvelocity(e), andnormalizedswirlvelocity(f). . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xiv 7.13 Actuator disk model solution: upstream (station 1) contours of local advance ra- tio (a), and normalized swirl velocity (b); and downstream (station 2) contours of pressurecoefficient(c),totalpressurecoefficient(d),normalizedaxialvelocity(e), andnormalizedswirlvelocity(f). . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.14 Propellerperformanceparametersasfunctionofwingangleofattack: comparison betweenbody-forcemodel(BFM)andactuatordisk(AD)modelpredictions. . . . 146 7.15 Wing force and moment coefficients as function of angle of attack: comparison betweenbody-forcemodel(BFM)andactuatordisk(AD)modelpredictions. . . . 147 7.16 Jet mixing produced by the body-force model (left column) and by the actuator disk model (right column), as seen from total pressure contours on axial planes at different streamwise stations. Vertical lines mark the propeller location, and propellerrotationisclockwise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 7.17 Total pressure contour on axial plane near the motor pylon at station 3: the swirl flowisblockedbythepylon,inducingaleftwardmotionatthebottomofthepylon. 150 8.1 Converging-divergingnozzle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 8.2 Validationoftheblockagemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . 171 xv Nomenclature b = Blockagefactor B = Numberofblades C p = Pressurecoefficient C p t = Totalpressurecoefficient C T = Propulsorthrustcoefficient C P = Propulsorpowercoefficient e t = Specifictotalenergy f = Body-forcevector f blk = Blockageforcevector f n = Inviscidbody-forcevector f v = Viscousbody-forcevector h t = Specifictotalenthalpy J = Propelleradvanceratio⌘ V • /(2nR) K v = Localviscousdissipationcoefficient K c = Compressibilitycorrectionfactor ˙ m = Massflowratein[kg/s] ˙ m corr = Correctedmassflow⌘ ˙ m p T t in p t in p t ref p T t ref M r = RelativeMachnumber ˆ n = Normalunitvectortobladecamber n = Rotationalspeedin[rev/s] p = Staticpressure p t = Totalpressure r = Radiallocation R = Propulsortipradius Re q = MomentumthicknessReynoldsnumber xvii s = Specificentropy T = Statictemperature T t = Totaltemperature ˆ t = Tangentialunitvectortobladecamber t = localbladethickness U = Bladevelocityvector =W r V = Absoluteflowvelocityvector W = Relativeflowvelocityvector x = Axiallocation a = Absoluteswirlangle;orAngleofattack a loc = Localangleofattack b = Relativeswirlangle d = Deviationangle d ref = Referencedeviationangle e = Flowturningangle⌘ b i b o h is = Isentropicefficiency h aero = Aerodynamicefficiency h prop = Propellerefficiency⌘ C T J/C P G = Bladeboundcirculation k = Localblademetalangle w = Losscoefficient W = Angularspeed[rad/s] q = Circumferentiallocation J = Bladesectioncamberangle f = Flowcoefficient⌘ V • /(W R) y = Totalpressurerisecoefficient p = Totalpressureratio r = Density ¯ ¯ t = Viscousstresstensor xviii () ref = Referencequantity () x = Axialcomponent () q = Circumferentialcomponent () • = Freestreamquantity () TE = Quantityatbladetrailingedge () = Massaveragedquantity ||·|| = Vectormagnitude xix Abstract This thesis presents a body-force based propulsor model that replaces the propulsor blades with a source volume in computational fluid dynamics (CFD) simulations, and that effectively repro- ducestheequivalentflowturning,workinput,bladeblockage,andlosses. Theprimarymotivation behind this work is to establish a propulsor model capable of interacting with local flow non- uniformities and producing physically accurate responses when implemented in aero-propulsive coupled CFD simulations. The model formulation is fundamentally physics-based, thus does not require prior knowledge of propulsor performance or empirical input, and only necessitates the blade’s geometrical information, such as a camber surface, thickness distribution and number of blades. An inviscid formulation for the body-force was previously developed in the literature and demonstrated its potential for predicting the inviscid distortion transfer effects. However, with- outconsideringlosses,blademetalblockage,andbladesolidityeffects,theinviscidmodelcannot accuratelyrepresentapropulsor’sperformance. In this thesis, an improved formulation with additional physical models is proposed and is shown to properly predict both propulsor work input and efficiency. Blockage models are intro- ducedasadditionalsourcetermstoaccountfortheblockageeffectproducedbythebladerowonto the passing flow. Corrections are introduced into the inviscid load formulation to address com- pressibility and solidity effects and ensure compatibility with the blockage model. Loss terms are includedtomodelbladeprofilelossesandshocklosses. Theproposedmodelisimplementedwithintheopen-sourceflowsolverADflowandvalidated against experimental data for two highly loaded turbomachinery fan stages (NASA Rotor 67 and Stage 35) and a low speed low-solidity fan stage (TF8000), operating in clean flow. The total xx pressure ratio shows good agreement and the choking mass flow was captured within a 1% error. Theisentropicefficiencyanditsoff-designtrendswerequalitativelywellpredicted. Thebody-forcemodelwasappliedtotheNASArotor67andTF8000fanstagetestcaseswith varioustypesofinlettotalpressuredistortionprofiles. Theinteractionsbetweenthefanstagesrep- resented via the body-force model and the inlet total pressure distortion were investigated. Com- parisons with full-annulus bladed CFD results demonstrated that the body-force model accurately predicted the distortion transfer behavior and local variations in loading and efficiency. The time- averagedflowfieldofhigh-fidelityURANSresultswereessentiallyreproducedbythebody-force model with much lower computational cost and the 1-2% drop in fan efficiency due to boundary layeringestion(BLI)wereaccuratelypredicted. Additional model formulations were proposed, including wake-induced corrections, tip loss correction,andaloadlimitcorrection,toimprovetherepresentationoffree-tippropulsorssuchas aircraftpropellersusingthebody-forcemodel. Themodelwasappliedtosmall-scalepropellersof various specifications and validated against wind tunnel test data. The body-force model success- fully generated the non-axisymmetric responses of the propellers to skewed inflow. Moreover, it demonstratedperformancetrendsconsistentwithexperimentalresults. Finally, the body-force model was utilized to represent propellers in the RANS simulation of a blown-lift wing test case where strong propeller-airframe interactions are present. Within a blown-lift wing system, distributed propellers blow onto a flapped wing to produce high lift. Comparisons were made between results produced by the body-force model and those generated by an advanced actuator disk model. The body-force model demonstrated superior performance in capturing the wing-propeller interaction effects, highlighting its effectiveness in dealing with complexinteractionscenariosthatarebeyondthecapabilitiesofmoretraditionalmodels. xxi Chapter1 Introduction 1.1 Motivation Due to the growing environmental concerns and decarbonization demands in the aviation indus- try, innovative aircraft concepts that leverage the benefits of aero-propulsive coupling are gaining increased attention. By tightly integrating the airframe and propulsion systems, these otherwise separate subsystems can undergo more complex interactions, which, if carefully designed, can lead to unprecedented improvements such as reduced fuel consumption and enhanced vehicular capabilities. Oneexampleofthisistheadoptionofboundarylayeringesting(BLI)enginesinrecenttrans- port aircraft concepts, as illustrated in Figure 1.1. In these designs, propulsors are integrated with theairframe,allowingpartoftheairframeboundarylayertobeingestedandre-acceleratedbythe propulsors. This approach results in an overall improvement in propulsive efficiency through the reduction of airframe wake dissipation and jet dissipation [1]. The concept has been employed in MIT’sD8DoubleBubbledesign[2]andNASA’sSTARC-ABLadvancedtransportconcept[3]. A total system-level benefit of nearly 19% was reported in the wind tunnel experiment of a powered sub-scaleD8model[4]. Anotherexampleofaero-propulsivecouplingcanbeseeninthenumerouselectricaircraftcon- cepts that have emerged in recent years, driven by the growing trend of aircraft electrification and 1 increasedinterestinurbanmobility,aspromotedbyNASA’sAdvancedAirMobility(AAM)cam- paign. Aero-propulsive coupled designs are often found in these emerging electric aircraft, which utilize electrically-driven propulsors connected only to energy sources via electric wires. This setup allows for greater flexibility in propulsor sizing and arrangement, enabling the full potential ofaero-propulsivecouplingbenefitstobeharnessed,particularlywhencombinedwithdistributed propulsion. In these designs, the thrust stream or jets created by an array of wing-mounted pro- pellersblowontoawingoraflappedwing,inodertoleveragethesynergisticinteractionbetween the propellers and the wing. This allows the aircraft to maintain the necessary lift at much lower speeds,ultimatelyachievingultra-shorttake-offandlanding(STOL)capability. Examplesofsuch configuration, commonly referred to as blown-lift designs (or blown-flap in some applications), are shown in Figure 1.2: NASA’s X-57 Maxwell and the eSTOL aircraft design created by Elec- tra.aero,Inc. Alongwiththeadvantagesandpotentialofaero-propulsivecouplingdesignsmentionedabove, there arise challenges concerning the design processes of such aircraft, specifically the accurate representation of complex aero-propulsive effects in high-fidelity full-aircraft computational fluid dynamics (CFD) simulations. With advancements in computational technology, CFD analyses have become an indispensable part of the aircraft design process, helping refine the aerodynamic design at a significantly lower cost compared to wind tunnel testing. Conventionally, propulsion Figure 1.1: Aircraft designs with boundary layer ingestion (BLI): (left) the D8 double bubble aircraft,and(right)theNASASTARC-ABLconceptualaircraft. 2 systems are installed away from the airframe, such that the propulsion system’s influence on air- frameaerodynamicsisminimal. Ontheotherhand,incoupledpropulsionsystems,thepropulsion effects actively modify the airframe flow field, and the propulsors operate under non-uniform in- flow conditions due to the airframe’s influence. The propulsor’s operating point may vary both radiallyandcircumferentiallyduetonon-uniformityinflowanglesandvelocity. Thus,itisessen- tial to properly represent propulsion effects to accurately predict system performance and enable simultaneousdesignofthecoupledpropulsionsystemandairframe. Fully resolving propulsor geometries in CFD offers the most accurate representation, but the high computational costs and reliance on fully-defined propulsor geometry make this method un- suitable for design. Instead, a simplified propulsor model is often used to represent propulsion effects in CFD. However, commonly used 1D models, such as an actuator plane, cannot capture theeffectofflowdistortion. Foraero-propulsivecouplingcases,thepropulsormodelmustprovide accurate predictions of overall propulsor characteristics under normal operating conditions while representing off-design performance variations under distorted inflow conditions. Therefore, a cost-effective, design-compatible propulsor model capable of accounting for locally off-design flowvariationsisneededtodesignnovelvehicleswithcomplexpropulsor-airframeflowcoupling. Figure1.2: Aircraftdesignswithblown-lift: (left)theElectrasupereSTOL,and(right)theNASA X-57Maxwell. (ImagebyElectra.aeroandNASA) 3 The main goal of this thesis is thus to address the propulsor modeling challenge within high- fidelity CFD. To that end, we propose an improved physics-based body-force model that can rep- resent propulsor work input, flow losses, and blade blockage effects while being able to account fortheinfluenceoflocalflowconditionvariations. 4 1.2 Background: PropulsorModelingMethodologies Various methods, with differing levels of complexity and accuracy, exist to represent propulsors suchasenginefanstagesandpropellerswithinComputationalFluidDynamics(CFD)simulations. A generalized categorization of these methods, demonstrated using NASA Rotor 67 geometry, is presented in Figure 1.3: Category (a) utilizes an infinitesimally thin disk, often referred to as the actuator plane method, where the propulsor’s influence is imposed as an internal boundary condition. Category (b), often known as the actuator disk method, uses a thin disk with a finite volumetoproduceprescribedforces. Thisistypicallyimplementedasasourcevolume,leadingto smoother property changes across the disk rather than a discontinuity. Category (c) depicts an ac- tuatorductmethod,whereadiskaccuratelyresolvestheaxialextentofthepropulsor. Thisresults in a disk with non-constant thickness. Category (d) represents the actuator line model, wherein the blade loads are concentrated into rotating lines that represent the blades. Unlike the previous cases, it accounts for the unsteady effects. Category (e) is the most detailed and accurate method, necessitatingfullyresolvedrotorbladerowgeometrywithinanunsteadysimulation. Thesemeth- ods collectively span a range of trade-offs between computational cost and simulation accuracy, therebycateringtoavarietyofpropulsormodellingrequirements. (a)Actuatorplane (b)Actuatordisk (c) Actuator duct or body-force (d)Actuatorline (e)Fullgeometry Figure1.3: GeneralizedpropulsorrepresentationswiththeNASARotor67geometry. 5 Figure1.3(d)illustratesthatthemostaccuratemethodtorepresentpropulsorsinComputational Fluid Dynamics (CFD) is to fully resolve the blade geometry with 3D unsteady simulations. This methodaccuratelypredictscomplexthree-dimensionaleffectswithintheblade-to-bladeflowfield using high-fidelity full-wheel simulations [5, 6]. While this approach may be suitable for focused propulsordesignduetoitshighlevelofprecision,itiscomputationallyintensive. Thelargenumber ofcomputecellsrequiredtoaccuratelyrepresentthebladewallsurfacesrendersthismethodoverly costlyandimpracticalforfull-aircraftdesignsimulations. Conversely,thepropulsioninfluencecanberepresentedmoresimply,asshowninFigure1.3(a). Here, a thin disk produces a uniform pressure jump, akin to the classic actuator disk momentum theory [7]. This actuator disk model signifies a limiting case where the number of blades ap- proaches infinity, commonly referred to as the infinite number of blade assumption. This uniform pressure jump represents a constant disk loading case, which, while simple, is the least realistic representation. The disk load can be distributed following prescribed distribution as a function of radius, as applied in [8]. More realistic disk load variation can be achieved through various methods, de- pending on the application type. Despite all propulsors operating on the same basic aerodynamic principles (using rotating blades to impart energy onto passing fluid flow), there are mainly two waystomodeltheseeffects,asevidencedintheliterature. 1.2.1 Propeller,HelicopterRotorandWindTurbineModeling In the context of aircraft propellers, helicopter rotors, and wind turbines, which are often referred toas”free-tiprotors”,thenumberofbladesistypicallylow. Thisallowseachbladetobetreatedas an independent rotating wing, as interference between blade-bound circulations is not significant. However,sincetheseapplicationsusuallyfeatureunshroudedbladetips,unlikeinturbomachinery applications, the effects of the wake and tip vortices must be considered when approximating the diskload. 6 Rotation dT V dD dL V U W ↵ W x ✓ Figure1.4: Force,velocityandflowanglesonabladeelement. Blade element theory based methods are the most common method to approximate the blade load for free-tip rotors. As illustrated in Figure 1.4, it relies on pre-computed airfoil polar data to providethesectionalaerodynamicforcesatcertainradiallocation, d L=B 1 2 r W 2 cc ` (a ,Re,M)dr d D=B 1 2 r W 2 cc d (a ,Re,M)dr d T =d Lcos(f ) d Dsin(f ) d Q=(d Lsin(f )+d Dcos(f ))r. (1.1) Note that the angle of attack, a , used to look up the polar data usually needs to be corrected, or iterativelydetermined[9],toobtainthecorrectflowangleaccountingforwakeeffects. Theoverall aerodynamicforcesarethenacquiredbyintegrating2Dsectionalforcesalongthespan. ThecouplingofBladeElementTheory(BET)andMomentumTheory(MT)wasinitiallyintro- ducedbyGlauert[10]in1935. ThisresultedinthedevelopmentoftheBladeElementMomentum (BEM)method,whichisnowwidelyadoptedinindustrydesigntoolsduetoitslowcomputational cost. The Momentum Theory applies the principles of axial and angular momentum conservation to connect the load predicted by the Blade Element Method with the blade-induced flow velocity in both the axial and tangential directions. This allows the flow angle used by the Blade Element 7 Method to be corrected. Since the Momentum Theory is based on an infinite number of blades assumption, tip losses are not inherently predicted and must be accounted for using methods like Prandtl’stiplossfactor[11]. The BEM model has become the standard method in the wind industry [12], with an unsteady BEM implementation included in the Aerodynamic module (AeroDyn) of the open-source wind turbinedesigntool,OpenFAST[13]. Forpropellersandairscrews,theBladeelementtheorywere used in combination with a lift line model by Lock as early as in 1945 [14]. A propeller analysis code,Qprop,wasdevelopedbyDrela[9],whichisbasedonbladeelementsandasimplifiedwake model. TheMomentumTheory(MT),whilewidelyadopted,hasitslimitations,includingtheassump- tions that the flow is axisymmetric and the blades’ wake remains a uniform helical surface. To providemoreaccuratecorrectionsontheflowangle,moreadvancedwakemodelscanbeapplied. Drela [9] integrated the blade element theory (BET) with a simplified vortex wake model, which connects the tangential induced velocity with blade bound circulation via the Helmholtz relation, enabling iterative solving for blade load. Another development in wake modeling was made by Eduardo and Ning [15] with the introduction of a vortex particle wake method. This method dis- cretizes the vorticity field into vortex particles, enabling them to be convected and diffused. In addition, a Free Vortex Wake (FVW) method was introduced in OpenFAST [16], which uses pan- els of vortex doublets to model the wake surface, allowing it to vary its shape in a time-accurate manner. Theseadvancedwakemodelsprovideincreasedaccuracyinpredictingthecomplexinter- actionsbetweenthepropellerorrotoranditswake. Thebladeelementtheory(BET)canalsobeusedwithinComputationalFluidDynamics(CFD) simulations to provide an estimation of local disk loading for an actuator disk model. A compre- hensive summary of actuator disk methods, with a focus on helicopter rotor modeling, can be found in the work by Le Chuiton [17]. Fejtek and Roberts [18] applied the actuator disk model to studywing-rotorinteractionintiltrotoraircraft. Intheirapproach,theloadwasdeterminedbased on the blade element momentum (BEM) method, but this was done in a time-averaged fashion. 8 Similarly, steady actuator disk models paired with the blade element theory have also been used for modeling wind turbines and wind farms [19, 20]. Unsteady implementations of these models were presented by Sørensen [21]. These methods effectively capture the physics of the flow in a morecomputationallyefficientmanner. For a more sophisticated representation of the rotor that models the full unsteady effects, ad- vanced methods like the actuator line model [22], as depicted in Figure 1.3(d), and the actuator surface model [23] have been developed, primarily for wind turbine applications. In the actuator linemodel,thebladeloadestimatedfromthebladeelementmethodisdistributedalongarotating line. The force on the blade is distributed following a Gaussian distribution centered along the span. Because of the individual blade effects fully represented in these models, there is no need for additional wake models in the actuator disk or actuator line models. These methods, which provide a higher level of detail and accuracy, allow for a more nuanced simulation of the rotor’s effectsontheflowfield,withouttheneedtoactuallyresolvebladegeometry. 1.2.2 TurbomachineryModeling In contrast to free-tip rotors, turbomachinery applications typically involve a high number of blades, wherein the blade rows function more like internal flow passages that direct the passing fluid to induce a certain angle of turning, rather than individual blades operating independently. In these cases, the blades are positioned close enough to each other, leading to interference of their bound circulation. As a result, their two-dimensional performance has to be evaluated using cascades rather than isolated airfoil. Furthermore, in the case of turbomachinery, blade tips are typically shrouded by the duct walls. This shrouding prevents the formation of the vortex sys- tem typical for propellers. As a result, the upstream flow angles aren’t affected by the tangential induced flow, as would be the case with propellers. This implies that for turbomachinery rotor models,thereisnonecessityforwakemodeling. 9 1.2.2.1 WorkInputandFlowTurning Assuming circumferential uniformity in the blade-to-blade flow of turbomachinery, the rotor load for such applications can be derived from the change in flow angle across the rotor. This can be demonstratedusingastraightforwardcontrolvolumeanalysisshownbelow. By combining angular momentum and energy conservation, Euler’s turbine equation can be writtenas h t 2 h t 1 =W (r 2 V q 2 r 1 V q 1 ). (1.2) Assumingconstantaxialvelocityandconstantradius,theaboveformcanbere-writtenas h t 2 h t 1 =UV x (tan(a 2 a 1 )), (1.3) or h t 2 h t 1 =UV x (tan(b 1 b 2 )), (1.4) whichindicatesthattotalworkinputisdirectlylinkedtoflowturning. Equation1.2canbebefurtherre-writtenas (h 2 + 1 2 V 2 2 ) (h 1 + 1 2 V 2 1 )=U 2 V q 2 U 1 V q 1 (1.5) Expressabsolutevelocityintermsofcomponentsinthecylindricalcoordinates V 2 =V 2 q +V 2 x +V 2 r (1.6) Alsorewritethecircumferentialcomponentofabsolutevelocityinthefollowingform V q =U W q (1.7) 10 Knowingthat W 2 =W 2 q +W 2 x +W 2 r =W 2 q +V 2 x +V 2 r (1.8) Wecanarriveat h 2 h 1 =W 2 1 /2 W 2 2 /2+U 2 2 /2 U 2 1 /2 (1.9) Assumingconstantradius(i.e.U 2 =U 1 ),staticenthalpyrisecanbesimplyexpressedasthereduc- tionofrelativekineticenergy h 2 h 1 =W 2 1 /2 W 2 2 /2 (1.10) Assumingincompressibleandisentropicflow,asimilarrelationcanbewrittenforpressure p 2 p 1 =r W 2 1 /2 r W 2 2 /2 (1.11) This shows that static pressure rise through a compressor is due to diffusion process in the blade passages. Morespecifically,theslowingdownofrelativevelocity. The derivation above clearly demonstrates that in a turbomachinery blade row, work input is directly associated with flow turning. As a result, a method known as the through-flow method is more typically employed to model load and performance in turbomachinery applications. In through-flow modeling, the influence of the blades on the flow is represented by specifying the flow turning based on a given blade geometry, as opposed to directly specifying the load implied fromairfoilpolardata. 1.2.2.2 Through-Flow(TF)Methods Traditional throughflow methods involve solving the radial equilibrium equations for turboma- chineryflowonanaxisymmetricmeridionalsurface(S2surface)[24]. Thesemethodsaccountfor the radial variation of flow properties. Two major solution methods are extensively adopted: the streamfunction method and the streamline curvature method. Hirsch and Denton [25] provided a comprehensive review of throughflow methods. Although this reference is somewhat dated, its 11 contents remain relevant and applicable to most of the throughflow methods in use today. The upcomingsectionswillbrieflysummarizetheevolutionandapplicationofthroughflowmethods. Wu [24] proposed a quasi-3D theory to solve flow by iterating between one S2 and several S1 surfaces. However,duetolimitedcomputationalpoweratthattime,themethodwasrarelyapplied. The iterative process was deemed too laborious and the source of error overwhelming. Despite these challenges, Wu’s method served as the foundation for later methods. In the mid to late 1960s,Smith[26]andNovak[27]independentlydevelopedthestreamlinecurvaturemethods,each employingauniqueapproach. Thismethodsubsequentlybecamethemostpopularmethodinuse. Aroundthesametime,Marsh[28]simplifiedWu’sgeneraltheory,solvingtheequationofmotion intheformofastreamfunctiononanS2surfaceusingfinitedifferences. Thismethod,knownasthe streamfunction method, was adopted by researchers as an alternative to the Streamline Curvature Method(SCM).Later,HirschandWarzee[29]appliedthefiniteelementmethodtosolvethesame streamfunction, an approach that was capable of handling complex geometries. However, this methodfacedmajorchallengeswhendealingwithtransonicflow,asthetypeofequationchanged. The most widely adopted streamline curvature formulation was described by Denton [30], which was simplified based on a later version of Novak’s formulation [31]. Denton also extended this method to multistage calculations. Adamczyk [32] developed a rigorous mathematical approach for three-dimensional calculations of multiple blade rows. His method incorporated sophisticated averaging techniques, including passage-averaged unsteadiness effects and Reynolds averaging effects. Effectsofotherbladerowsweremodeledasadditionalsourceterms. Thesourcetermsthat Adamczyk developed have been widely adopted by numerous other axisymmetric models. One shortcomingformosttraditionalthroughflowmethodsistheinabilitytopredictrecirculatingflows duetothespatiallymarchingnature. AsDentonpointedoutinhisreview,throughflowcalculations arelittlemorethanvehiclesforinclusionofempericismintheformoflosscoefficients,deviations andblockagecorrections[33]. Forourpurposes,thetraditionalthroughflowisnotsuitablesincethe methodareessentiallyaxisymmetricalandwon’tbeabletohandlecircumferentiallynon-uniform inletcondition. 12 1.2.2.3 Body-ForceMethods In the 1960’s, Marble [34] introduced a forcing term into the equation of motion to represent the local blade force exerted on the fluid, and thus introduced the idea of re-casting the through- flow method as an additional term in Euler’s or Navier-Stokes equations. Since the blade force is distributedwithinthebladeregionasforceperunitvolume,theapproachofaddingasourcetermto thegoverningequationsisoftenreferredtoasabody-forcemethod. AsCFDtechniqueadvances, methods solving more general governing equations, such as axisymmetric Euler’s equation and RANS, became available. One of the advantage of these CFD-based methods is that there’s no limitation when dealing with transonic and supersonic problems. In 1980 Spurr [35] implement the body-force into steady-state Euler’s equation and solved it with a time-marching technique. The force was iteratively updated such that the change of angular momentum produced a relative flow angle that follows the mean blade camberline direction. This way of iteratively solving the body-forcesothattheresultedflowbecometangenttoaprescribedcambersurfaceisreferredtoas flow-tangency method. Later, the representation of blades via a camberline and a body-force was adoptedandimprovedbymanyauthors[36,37]withlossmodelingdoneviacorrelations[38,39], but the implementations largely remained to be axisymmetric. Damle and Dang [37] presented a body-force model in axisymmetric Euler’s equation with both an analysis mode in which flow- tangency is enforced, and an inverse mode in which the swirl distribution is prescribed and the bladecambersurfaceleadingtotherequiredswirldistributionisresulted. Gong[40]implemented the body-force method in the full three-dimensional Euler’s equation. He developed an analytical formulation for the body-force, which represent the blade response to local deviations in flow anglerelativetothebladecamberline. Heshowedthatthemodeliscapableofrespondingtolocal flow variations due to distortion. However, his formulation requires empirical calibration based on data from measurements or bladed CFD computations. Peters [41] improved Gong’s model and applied the model for ultra short fan nacelle design. Thollet [42] proposed a modification of Gong’s model based on lift and drag analogy reducing the need for calibration with multiple operating points. More recently, research efforts have been focusing on body-force formulation 13 without the need for empirical calibration. Hall [43] reformulated the turning part of body-force based on local flow variables and showed that the model can predict the inviscid redistribution effect of a BLI propulsor. However, since neither losses nor blockage effect were included, the model is still un-calibrated and unable to predict propulsor performance accurately. Akaydin and Pandya [44] applied Hall’s model to several fan stages and showed that Hall’s body-force model hasgoodpotentialasanaffordableandaccuratemodelingmethodforon-designcaseswithinflow distortion. Pazireh [45] proposed a loss formulation to use along with Hall’s turning model based on an artificial neural network trained with a dataset of cascade computations and showed good predictionforturbomachineryefficiency. A survey of the flow tangency methods used to solve for the body-force is given below. Early body-forcemethodapplications[37,46,47,39]solvesanadditionalequationtocomputethebody- forcebasedonflowtangencytothebladecambersurface. Theseapproachesarethereforereferred toasflowtangencymethodsinthiscontext. Damle[37]suggeststheinviscidbladeforceinthefollowingform. f b = f b (r,z)—F (1.12) whereF =q + f(r,z)istheexpressionforthemeancambersurfaceoftheblades. VariationofF inzrepresentsbladestaggerandvariationinr representsbladelean. When f =0,F =q meaning thatstreamsurfacesaresurfaceswithconstantq . f b =V ·— V q f v |W q | |W| (1.13) (1.13)suggeststhatwhentherearelossespresent,sameamountofbladebodyforceproducesless flowturning. ThetermV ·— V q isestimatedsothatflowtangencyissatisfied. W ·—F =0! V ⇤ q =r(w +V r ∂ f ∂r +V z ∂ f ∂z ) (1.14) 14 V ⇤ q solvedin(1.14)isusedinthebody-forceexpression(1.13)toestimatedthebodyforcemodulus required. Baralon [46] proposed adding a separate time dependent equation to solve for the modulus of thebladeforce. Thebladeforcewouldeventuallyreachasteadystatewhenflowtangencyismet. ∂ f b ∂t =Cr (V x n x +V r n r +(V q W r)n q ) (1.15) where n is the normal vector of the mean flow surface. The rate of change of f b is proportional to the deviation from tangential condition. (1.15) is changing the magnitude of the body-force such that the flow reach tangential condition to the prescribed streamsurface at steady state. The parameterC controlshowfastistherespondoftheforceforgivenamountofflowdeviation. A penalty formulation was proposed by Persico(2012) [48] to determine the modulus of the bladeforce. L= K(w·n g ) (1.16) Intuitively, the modulus of the force (denoted as Lift force L here) is defined to be proportional to thedotproductofrelativeflowvelocityandthebladenormalvector. NotethatLwillbezerowhen flowtangencyatstreamsurfaceismet. PenaltyfactorK isempericallydefinedby K =K min +(K max K min )exp max(w·n 2 g ) s 2 ! (1.17) so that K is larger when the dot product is small (close to meeting flow tangency). This way flow tangencyconditioncanbeguaranteedwhensolverconverged. Thismethodensuredstabilitywhile eliminatingtheneedtosolveanadditionalequationasBaralon’smethodsuggested. The resulting force distribution would turn the flow so that it exactly match the blade camber. This method over-predicts the turning and thus the loading because the real flow leaves the blade trailing edge with a deviation angle dueto viscous de-cambering effect. The additionof deviation 15 correlationsisthenrequired. Thismethodisreferredtoastheflowtangencymethodintherestof thispaper. Actuator disk methods had also been used to model turbomachinery applications, utilizing the principles of through-flow methods. In the context of turbomachinery modeling, the actuator disk is typically a discontinuity that produce certain amount of flow turning and entropy rise, rather than providing specified amount of load [49]. An application of actuator disk model to simulate turbomachinery blade rows were presented by Joo [50] in which the actuator disk with zero thickness is implemented via boundary condition and the effect of rotor is specified as flow anglechangesandentropyriseacrossthedisk. 1.2.3 ProsandConsofModelingMethodologies While the actuator disk model with blade element-based method is reliable and fairly accurate, it hasseverallimitationspreventingitfromservingasageneralizedmodelforpropulsorsofvarious solidities. These include: (1) It requires extensive input data for utmost accuracy, including blade sectional performance data at different radii, which must account for variations in Mach number and Reynolds number. This requirement becomes particularly demanding for high-solidity tur- bomachinery applications, which necessitate cascade data with various pitch distances. (2) The method encounters difficulties when applied to cases featuring a large number of blades, such as bladerowsinturbomachinery,asitoftenneglectstheeffectofblade-to-bladeinteractions. (3)The modelassumesaspecificloaddistributioninthechordwisedirection,whichmaynotaccuratelyre- flecttheactualchordwisedistribution acrucialfactorinmodelingdistortionpropagationthrough arotor. (4)Itfaceschallengeswhendealingwithcasesfeaturingsubstantialradialflowvariations, as these changes effectively alter the shape of the blade sections, an adjustment not accounted for in the input aerofoil performance data. (5) The model lacks design compatibility. The propulsor geometry needs to be fully defined to apply the model, which might not be available during the designphase. Anygeometricchangesrequirerecomputingtheperformancedata. 16 Ontheotherhand,through-flowbasedmodelsrequireempiricaldatathatmightnotexistforall types of propulsors, which limits their generality. Clearly, there is a gap between free-tip propul- sor models and turbomachinery models. Ideally, a single model could fill this gap, capable of predicting both propeller and turbomachinery types of propulsors. The body-force model with load determined via an analytical method emerges as a potential candidate due to the following advantages: (1) The analytical formulation can be designed to approach thin airfoil theory when the blade number is small and exact turning when the blade number is large, satisfying the needs of both high-solidity and low-solidity applications. (2) It does not require resolving the boundary layers on blade surfaces, thus significantly reducing computational costs. Furthermore, flow turn- ingcanbeproducedwithoutunsteadycomputationsortheuseofarelativeframeofreference. (3) Astatorcanbeeasilyrepresentedwithinthesamedomainandinthesamemannerbyspecifyinga zero-rotational velocity, allowing for the representation of a full rotor-stator stage without chang- ingreferenceframes. (4)Theforceisevaluatedbasedonlocalflowconditions,allowingthemodel torespondtonon-axisymmetricflowvariations. Thismeansthatthemodelcanautomaticallyhan- dle flow distortion transfer across the blade region when faced with a non-uniform inlet. (5) The definitionofthebladescanbelimitedtotheircamberlineandthicknessdistributions,allowingfor amoredesign-friendlyrepresentationwithfewerdegrees-of-freedom. Thegeometryspecification requirementcanbeflexibleintermsoflevelsofdetail,makingitpossibletouseasimplegeometry generatedviadesigncodeasaninput. Afterathoroughcomparison,thebody-forcemodelemergesasastrongcandidateforadesign- compatible model applicable to both high-solidity turbomachinery and low-solidity rotors, while alsocapableofcapturingaero-propulsivecouplingeffects. 17 1.3 ThesisScope The main objective of this thesis is to develop a physics-based propulsor model that leverages the body-force modeling approach for use in the analysis and design of aircraft with aero-propulsive couplingconfigurations. Theoverarchinggoalsforthemodelareasfollows: • Accurately represent the performance of a rotor-only propulsor or a single-stage propulsor undercleanflowoperatingconditions,withinareasonabledegreeofaccuracy. It’sacknowl- edged that achieving perfect accuracy for all types of propulsors using a body-force model is unrealistic due to the inherent simplification of flow physics in the model. However, with simpletuningparameters,thismodelshouldbepossibletocloselymatchtheoverallcharac- teristicswithknownperformancedata. • Becapableofcapturingtheinteractionsbetweentheairframeandpropulsorthatarepertinent forpredictingtheeffectsofinletdistortiononperformance. • Be able to represent the propulsors at a level of fidelity and computational cost that is con- ducive to design, enabling simultaneous optimization of the airframe and propulsion sys- tems. Notably, it should not necessitate the design of an actual propulsor in advance, but instead require only minimal blade geometry specifications such as camber lines and thick- nessdistribution. Otherareasoffocusofthepresentworkinclude: • Verifyingthattheproposedmodelcorrectlyrepresentsphysicalphenomenabycomparingit withtheoreticaltrends. • Evaluating the proposed model’s capability to predict stand-alone propulsor performance in uniform flow across a wide variety of test cases for both turbomachinery and free-tip propulsors. • Assessingtheproposedmodel’seffectivenessinpredictingperformanceofapropulsorwith adistortedinletflow. 18 • Assessing the applicability of the proposed model in scenarios featuring aero-propulsive coupling. This thesis is organized as follows: Chapter 2 introduces the body-force model concept and formulation, which covers the inviscid turning model, viscous losses model, and blockage model. ThisisfollowedbyChapter3,wherethecomputationalapproachandtheimplementationmethod are detailed. The subsequent Chapter 4 provides an evaluation of the model on isolated turboma- chinery propulsors, validated against experimental data and high-fidelity URANS results. Chap- ter 5 then presents propulsor test cases with distorted inflow to further validate the body-force model’s ability to accurately predict distortion responses. In Chapter 6, an extended formulation for free-tip propulsors is provided. Finally, Chapter 7 demonstrates an application to a blown-lift wing test case to highlight the potential of the model in the context of aero-propulsive coupling design. 19 Chapter2 TheBody-ForceModel This chapter presents the development process of a body-force propulsor modeling framework. Model components and formulations representing different physical effects are proposed and in- corporated,withtheoverallgoalofmakingmoreaccuraterepresentationofpropulsorperformance inCFDacrossvariousapplicationswithminimumempiricalinformationrequired. Asastartingpoint,Hall’sinviscidloadingformulation[43]wasadoptedasthebaselineofthe model,providingaphysics-basedanalyticalmethodforproducingflowturningandcontributingto themajorityofworkinput. Sincetheturningmodeltakesaccountofflowconditionlocallyineach of the cells in the model region, it is capable of interact interacting with distorted flow condition and predict physical propagation and transformation of the flow distortion through the propulsor. Thisisamajoradvantageforaero-propulsivecouplingapplications. Hall’s model can be understood as a differential form of thin airfoil theory, which provides a local estimation of inviscid load based on local flow deviation from the blade surface. However, assumptions inherited from thin airfoil theory may limit the model’s accuracy and require careful consideration. Firstly,theincompressibleflowassumptionmustbeaddressed,asturbomachineryflowistypi- callycompressibleandthecompressibilityeffectsmustbeaccountedforwhenestimatingloading. AcompressibilitycorrectionfactorispresentedinSection2.3. Secondly, the implications of varying blade spacing or solidity must be considered. As illus- trated in Figure 2.1, as the gap between blades diminishes, or the solidity increases, the nature of 20 the flow around the blades transitions from external to internal. In applications such as propellers andwindturbines,theinteractionsbetweenadjacentbladesareminimal,allowingeachbladetobe considered as operating independently. In contrast, in high-solidity applications like compressor blade rows, the blade spacing is much smaller, causing the blade row to be seen as an internal flow passage directing the flow through a turn. Studies on 2D potential flow indicate that as the bladespacingdecreases,theblade-boundcirculationandliftalsodecrease. Assolidityapproaches infinity, we encounter the condition of an infinite number of blades, where the flow field is cir- cumferentially uniform, and an uniform flow turning is achieved along the blade camberline. In this condition, there’s no concept of individual blades, and therefore no blade-bound circulation is resolved. The solidity effect as predicted by the inviscid body-force model is investigated in Section4.2. Anothersignificantimplicationtiedtothevariationofblade-boundcirculationwithbladespac- ing pertains to the impact on inflow angles upstream of the blades. The blade-bound circulation generates an upwash in front of the blade, effectively increasing the incidence angle seen by the leading edge. The adoption of the infinite number of blades assumption in the body-force model impliesthatthisupwasheffectcannotbeautomaticallyproducedbythemodel,astheblade-bound circulation isn’t part of the flow field. This upwash effect alters the incidence angle and, conse- quently, the deviation angle perceived by the model, which in turn influences the load prediction. A flow angle correction is developed in Section 2.3 to account for the upwash effect on local flow angle. Infree-tipexternalflowapplicationslikepropellers,thepressuredifferencecannotbesustained atthetipbecausetheflowcanleakaroundthetip,leadingtoareductionofliftnearthetip. Thistip losseffectmustbetakenintoaccount,andproceduressimilartothoseemployedinbladeelement- based methods can be adapted for the body-force model. A wake correction for the body-force modelisdevelopedandshowninSection6.2. Beyond inviscid loading, it’s essential to incorporate flow losses from different sources, such asprofilelossesfrombladeboundarylayersandlossesthatoccurwhenoperatingoff-design. This 21 (C) innite number of blades 1 2 2 2 3 (B) Cascade (A) Individual blade A 1 A 2 (Small numer of blades) (Large numer of blades) L D 1 > 2 > 3 ≥ 0 Figure2.1: Bladenumbereffects inclusion allows for the depiction of efficiency variations across a wide spectrum of operating points. Loss generation can be expressed simply based on local flow conditions and local blade loading, drawing from established boundary layer theories. A loss formulation is developed and showninSection2.4. Moreover,ininternalflowapplicationswithahighnumberofblades,suchasenginefanstages or compressor stages, the thickness of the blades creates a blockage effect on the passing flow. This effect limits the maximum possible flow rate (or choking mass flow rate) through the ma- chine. Representing this effect is crucial to accurately depict the choking behavior at higher mass 22 flow rate. In Section 2.2, we adopted a blockage model method widely used in early body-force methods. In the following sections of this chapter, we will introduce the body-force formulation that takesintoaccounttheabove-mentionedconsiderations. 2.1 Body-ForceModelingApproach The body-force modeling approach inherits the idea of turbomachinery through-flow models that aim to create the necessary flow turning via blade camber shape, instead of directly specifying bladeloadinglikeBEMmethods. AsillustratedinFigure2.2,aforcefieldwhichspanstheregion sweptbythepropellerbladesisimposedonthepassing-throughflowtomimicthepitch-averaged effect of the turning blades without actually including the blade surfaces in the flow domain. As suggestedbyMarble[34],thelocalbody-force,f,isimplementedassourcetermsinthreedimen- sionalNavier-Stokesequationsasfollows ∂r ∂t +— ·(r V)=0, (2.1) ∂ ∂t (r V)+— ·(r VV+p ¯ ¯ I ¯ ¯ t )=r f, (2.2) ∂ ∂t (r e t )+— ·(r Vh t ¯ ¯ t ·V)=r f·U. (2.3) Here r , p, V, e t , and h t are the flow density, pressure, velocity, specific total energy, and specific total enthalpy, respectively, while ¯ ¯ t is the shear stress tensor and ¯ ¯ I the identity tensor. The blade’s rotation is represented via its blade velocityU (zero for a stator stage) and calculated as the cross productofrotationrateandlocalradialdistancetotheengineaxis. Thelocalbody-forcedirection andmagnitudearedeterminedbasedonlocalflowvariablesandbladecambersurfacenormalsvia modelingrelations. 23 Figure2.2: Body-forcerepresentationofatransonicfanrotor: (left)actualpropellergeometry,and (right)body-forcevolumeinthecomputationaldomainshowingcamberlinenormals. The local body-force, i.e., force per unit mass, are determined based on local flow variables andbladegeometryinformation. Itcanbeexpressedingeneralasfollows f = fn(B,ˆ n(x,r),t(x,r),W ,V), (2.4) where ˆ nisthelocalbladecambersurfacenormaldistribution,t isthebladethicknessdistribution, W thewheelrotationalvelocityandVthelocalflowvelocity. Itisnotedthatthereisnoazimuthal variation in the blade geometry distributions, and so the turning blade row is modeled in a pitch- averaged— and hence axisymmetric—way. However, this modeling approach can still produce a non-axisymmetric response when the incoming local flow (V) is non-uniform, and therefore the body-forceapproachiscapableofcapturingtheeffectofinflowdistortion. The body-force contribution to the energy equation, f·U =W (f⇥ r), represents a local shaft power, and can be better understood by writing the blade velocity U in terms of absolute and relative (rotating-frame) velocities, V and W via the velocity triangle (V =U+W). The source termontheRHSoftheenergyEquation(2.3)thusbecomes r f·U = r f·V r f·W. (2.5) 24 The term r f ·V represents the work done by both inviscid and viscous body-forces. To show the physical implication of the term r f·( W), we derived the entropy equation from the energy equation,withthedetailedderivationshowninAppendixA.2, T Ds Dt = T ˙ s = 1 r ( ¯ ¯ t :— V)+f·( W), (2.6) where ˙ s =Ds/Dt is the total rate of entropy change. Thus, the term r f·( W) in Equation (2.5) contributestoentropygenerationsimilarlytohowtheviscousdissipationterm 1 r ( ¯ ¯ t :— V)produces entropy. Insteadystate,therateofentropygenerationbyjustthebody-force, ˙ s b ,canthenbewritten as T ˙ s b =TV·— s b =f·( W). (2.7) ThisconclusionshownbyEquation(2.7)impliesthatonlytheforcecomponentintheopposite direction to the relative velocityW contributes to the production of entropy. Therefore, the body- force vector is split into two parts: f =f n +f v and modeled separately. The inviscid body-force, f n , perpendicular to W and produces deflection of the flow without incurring losses (f n ·W =0). The viscous body-force, f v , in the direction opposite to W, represents viscous dissipation and is solelyresponsibleforentropygeneration. Bothoftheseforcescontributetotheshaftloading. The magnitudeoftheinviscidandtheviscousbody-forcesarecomputedviaseparatedphysicalmodels thatareintroducedinthefollowingsections. 2.2 BladeMetalBlockage The effective reduction in mass flow due to blockage by blades with finite thickness is accounted forbyincorporatingablockagefactorasdefinedbelow, b = 1 t (2p r/B)|n q | , (2.8) 25 into the governing equation, where t is the local thickness perpendicular to the camberline of the bladesectionsand |n q |approximatesthecosineofthelocalstaggeringangle. Thisfactorissmaller than unity inside the body-force volume and equal to one outside. While the inviscid and viscous body-forcecomponents,f n andf v ,onlyrequireknowledgeofthebladecamber,thebladethickness distribution is need to model the blockage effect. Blade metal blockage effects become important only when the number of blades is large, such as in turbomachinery cases, or when the blade thickness is substantial. Without the metal blockage effect included, the mass flow rate or flow coefficientrequiredtoachievecertainamountofpressurerisewillbeover-predicted. Thechoking massflowrateforhigh-speedcaseswillalsobeover-predicted. Tomodeltheblockageeffect,aforcingterm,f blk =p— b/b,needstobeaddedtothemomentum equation,whichrepresentsthelateralpressureforcefromaconvergentordivergentwallinaquasi- 1D way. A control volume analysis also dictates that the blockage factor be applied to each term ofthegoverningequations,whichthenbecome ∂ ∂t (br )+— ·(br V)=0, (2.9) ∂ ∂t (br V)+— · h b(r VV+p ¯ ¯ I ¯ ¯ t ) i =bf blk +br f, (2.10) ∂ ∂t (br e t )+— ·[b(r Vh t ¯ ¯ t ·V)]=br f·U. (2.11) Implementationoftheblockagefactoristypicallydonebyeithermodifyingthesolver’simple- mentation of spatial integration, or artificially reducing cell volumes and surface area in specific directions, but either way requires intrusive modification of the solver’s original source code. An alternative and simpler approach was adopted in this work. By manipulating the governing equa- tions,thecontributionofblockagecanbeisolatedintoseparatesourceterms ∂r ∂t +— ·(r V)= 1 b (r V·— b), (2.12) 26 ∂ ∂t (r V)+— ·(r VV+p ¯ ¯ I ¯ ¯ t )=r f+ 1 b ( ¯ ¯ t r VV)·— b, (2.13) ∂ ∂t (r e t )+— ·(r Vh t ¯ ¯ t ·V)=r f·U+ 1 b ( ¯ ¯ t ·V r Vh t )·— b. (2.14) This way, blockage appears as additional source terms that can be easily implemented. Note the contributionofthestresstensor ¯ ¯ t totheblockagereflectedintheright-hand-sideofEquation(2.13) and(2.14). Several previous works implemented metal blockage effects into pitch-averaged through-flow modelsbasedonEuler’sEquation[46,47],andshowedthataxisymmetricshockscanbeproduced that represent pitch-averaged blade-to-blade normal shocks, with the shock strength being set by the blockage factor. The inclusion of a blockage factor via a source term was previously adopted forinviscidflow[42]. 2.3 InviscidLoading Theinviscidloadingforcegeneratesflowdeflectionthroughaspecificcambershape. Flowturning typically goes hand in hand with pressure gradients resulting from flow curvature, which is akin to the lift generation mechanism due to pressure. Two methods have been employed to determine thedistributionofinviscidbody-forcebasedonbladegeometry: theflowtangencymethodandthe analyticalmethod. As introduced in Section 1.2, the flow tangency methods utilized an additional equation to solve for the body-force based on flow tangency to the blade camber surface. The resulting force distribution would turn the flow so that it exactly match the blade camber. This method over- predictstheturningandthustheloadingbecausetherealflowleavesthebladetrailingedgewitha deviation angle due to viscous de-cambering effect. The addition of deviation correlations is then required. 27 Pitch Average f n Inviscid loading body-force field Figure2.3: Illustrationoftheinviscidturningmodel. An alternative approach is to model the inviscid loading with an analytical formulation. It providesaphysicalestimationofthebladeloadingbasedonlocalrelativeflowdeviationfromthe blade camber surface. The core of the turning force formulation is based on a form developed by Hall[43]withtheforcemagnitudecomputedasbelow f n ⌘ (2pd ) 1 2 W 2 (2p r/B) |n q | , (2.15) inwhichn q isthecircumferentialcomponentofthecambernormalvector,rislocalradialcoordi- nateandBisnumberofblades. d islocaldeviationanglewhichmeasuresthelocalflowdeviation from a fictitious blade camber surface. The force magnitude scales with d with a factor of 2p , analogous to thin airfoil theory’s “c ` =2pa ”, since this force is the result of a local pressure dif- ferenceproducedbytheflowcurvaturecreatedbytheblade. Theforceperbladeisthendistributed acrossthebladepitch(2p r/B)resultinginaforceperunitmass. Due to circumferentially averaged modeling approach, only the axial extent of the propulsor bladesisrepresentedintheflowdomain,hencetheactualchordlength,orliftingsurfacepermodel volume,needstobecorrectlyreflectedusingthelocalstaggerangleofthebladesurface. Knowing that |n q |⇠| cos(k )|, (2.16) wherek isthebladestaggeringangle,afactorof1/|n q |wasincludedinthedenominatortoaccount forliftsurfacescaling. Thiswasexplainedin[51]inmoredetails. 28 Theinviscidbody-forceworkstoreducethedeviationangleandcreatesanaturalflowcurvature following the camber curves. At the beginning of the computation, the initialized flow has very largedeviationangleandthusverylargeturningforcesthatinitiatestheflowturning. Thedeviation angle reduces in the following iterations when the flow curvature gets generated by the body- force. When the model converges, the deviation angle no longer changes and this is the result of the equilibrium between the turning force and the natural de-cambering effects of the flow. The deviation angle d distribution represents the flow’s natural response to certain amount of loading via given camber surface, in the sense that it is resulted from the flow simulation instead of a quantitythatcanbeknownapriori. A fundamental difference between the body-force model and blade element based methods is that the former method provides a localized estimation of blade loads, rather than distributing a pre-computed integrated force. The body-force model provides a more physical chordwise load distributionandthusphysicalpropagationofinflowdistortiondownstream. The deviation angled plays a crucial role in the modeling approach. Not only does it provide a way of estimating the local inviscid loading, it is essential in representing the load variation due toeitherdifferentbladegeometries(i.e. numberofblades,solidityandcamberangle)oroperating points(i.e. incidenceangle). Ingeneral,thedeviationanglereduceswithhigheramountofloadingandisafunctionofblade geometryandflowangles. d =fn(i,k ,s/c,q ). (2.17) Higherincidenceangleiandhighercamberangleq willcauselargeramountofflowdeviation and thus results in more turning and larger blade loading. Increasing number of blade or reducing blade spacing (s/c) increases the inviscid loading and in turn reduces the deviation angle d . This trendistrueas,inreality,thereducingbladespacingcausesblade-to-bladeinteractionsthatreduce the lift and circulation from each blade. Hence the model possess the mechanism to reflect the effect of blade numbers. However, the actual trend predicted by the body-force model has yet to 29 be validated. Validation efforts are made in the result sections where the predicted blade loadings areshownwithdifferentnumberofblades. In the limit of small number of blades (B! 1), such as the case of propellers, the distributed local loading is so small that chordwise change of deviation angle can be assumed to be small, or inotherwords, thecreatedflowswirlissmall. Foranairfoilwithangleofattacka , thefollowing relationcanbewrittenforthedeviationangle Z c 0 2pd dx⇠ 2p (a +a 0 )c, (2.18) whichimpliesthattheloadpredictionfrombody-forcemodelissimilartothatpredictedusingthin airfoil theory. This argument serves as a first order justification for the load predicting capability of the body-force model when applied to small number of blade applications. As the number of blades increases and approaches infinity, a significantly larger amount of flow turning is created thatcloselyfollowsthecambersurface,whichisappropriateforhigh-solidityapplications. 2.3.1 ImprovedAnalyticalFormulation Inwhatfollows,weproposeanimprovedinviscidloadingformulationbasedonHall’sformulation. Acompleteformoftheupdatedinviscidbody-forceiswrittenas f n ⌘ K c 2pd 1 2 W 2 (2p r/B)b|n q | ˆ e n , (2.19) Toaccountforthechangeinpitchdistancecausedbyblockage,theblockagefactorbisintroduced inthedenominator. Inaddition,compressibilityeffectsareaccountedforbyincludingacorrection factor,K c ,withform K c ⌘ min 1 p |1 M 2 r | , 1.66 ! . (2.20) 30 1 corr target BFM prediction with no induced angle correction actual propulsor performance V 1 V corr shift w w⌘ upwash Figure2.4: DemonstrationoftheeffectofnoinducedanglecorrectiononBFMloadpredictionsfor propulsorswithasmallernumberofblades: thesameamountofloadingcorrespondstoasmaller flowcoefficientor,inotherwords,under-predictionoftheloadingataspecificflowcoefficient. Forsimplicity,K c representsamodifiedPrandtl-GlauertconstantbasedonthelocalrelativeMach number M r , namely by setting and using this same factor in the presence of transonic and super- sonic flow. The maximum value is caped at 1.66, which is the value of the Prandtl-Glauert factor atMach0.8,inordertoavoidthenon-physicalbehaviornearsonicconditions. Thelocaldeviationangleisfoundvia d = sin 1 ✓ W·ˆ n ||W|| ◆ +D d . (2.21) ThefirstterminEquation(2.21)representsthelocalrelativeflowdeviationfromthebladesurface, andD d wasintroducedasacorrectiontoaccountfortheeffectofbladeboundcirculationsonflow angles. In two-dimensional flow over an isolated airfoil, the bound circulation creates an upwash in front of the airfoil, as such, effectively increasing the incidence angle felt by the leading edge. In blade-element based methods (BEM), this bound circulation effect is automatically taken into accountbyusing2Dbladesectionperformancedata. 31 ↵ loc = +0.5( TE ) W 0.5( TE ) Figure2.5: Pictorialdemonstrationoflocalangleofattack. However,inthecaseofBFM,thebladesarerepresentedinacircumferentiallyaveragedfashion undertheinfinitenumberofbladeassumptionandthere’snonotionofindividualblades,thus,the boundcirculationeffectscannotbereflected. Asaresult,whenappliedtolow-blade-countpropul- sors such as propellers, the BFM under-predicts loading due to smaller incidence angles, which effectively shifted the loading characteristic to the left, as illustrated in Figure 2.4. Therefore, a flow angle correction is needed to increase the incidence angle simulating the bound circulation effects. Note that this bound circulation correction is not to be confused with wake-induced flow corrections used in BEM, which is a three-dimensional effect due to downwash caused by the trailing vortex system of the propulsor. It is also worth noting that the bound circulation correc- tion is typically not needed for turbomachinery that has large number of blades, considering that the bound circulation is significantly reduced by the influences from the closely-spaced adjacent blades; As number of blades increases, the baldes function transits from producing lift force to producing flow turning, diffusion, and pressure. For this reason, such correction is rarely seen in turbomachinerythrough-flowmethods. ThelocalizedfashionofmodelinginBFMposesasignif- icantchallengeforformulatingacorrection,asintegralquantitiessuchastotalliftanddragforces arenotreadilyavailableduringruntime. Inwhatfollows,wedevicealocalizedcorrectionformula basedonlocalcirculationestimatedfromlocalbladeloadingforce. Recognizing the fact that a cambered airfoil should generate positive lift at a zero angle of attackdefinedwithrespecttoitschordline,weinferthatthecamberangleshouldbeconsideredin 32 thecorrection. Assuch,wedefinealocalangleofattackthatincorporatesthelocalflowdeviation andcamberangle,asfollows: a loc =d + 1 2 J rem (2.22) where d is the local deviation angle, J rem is the remaining camber angle at a certain chordwise location,definedas J rem =k k TE , (2.23) and k is the local blade metal angle. The angles d , J rem , and a loc are schematically shown in Figure 2.5. Using the Kutta-Joukowski theorem, a local circulation can be related to a local lift forcecontributionapproximatedusingthethinairfoiltheoryas L loc =r WG loc =2p K c a loc ✓ 1 2 r W 2 l ◆ , (2.24) wherel isacellcharacteristiclength. Thus,thelocalcirculationG loc canbedefinedas G loc =p K c a loc Wl . (2.25) Thetangentialvelocityinducedbyapointvortexisgivenas v q = G loc 2p r . (2.26) where r is the distance from the center of the point vortex within a 2D blade sectional plane. CombiningEquation(2.24),(2.25),and(2.26),thechangeofflowvelocityduetolocalcirculation canbewrittenas w=v q = p K c a loc W 2p (r/l) . (2.27) Thechangeofdeviationangleisthus D d ⇡ w W = K c a loc 2(r/l) . (2.28) 33 Assumethat r l =1.0, (2.29) thenwecanwrite D d ⇡ w W =0.5K c a loc . (2.30) It’s recognized that this local form of correction only provides an estimation of the integral effect and is not expected to be able to fit every cases perfectly, given that each propulsor has its unique design that might also deforms and flexes in different ways. Therefore a tuning factor K t , with a default value of 1, is included in order to help achieve better matching of a certain know performance characteristic. The solidity effect also has to be accounted, because as the number of blade increases, the induced effect from each blade reduces due to influences from adjacent blades. Assume that the induced effect disappears at small s/c and takes full effect at s/c =4.0 thecorrectionbecomes D d ⌘ 0.5K t K c a loc ✓ min(s/c,4.0) 4.0 ◆ (2.31) Thedirectionoftheinviscidforceisfoundby ˆ e n ⌘ f n ||f n || = sign(W·ˆ n) W⇥ (W⇥ ˆ n) ||W⇥ (W⇥ ˆ n)|| , (2.32) suchthattheforceisnormaltotherelativeflowandturningonlyoccurswithinthesurfaceformed by the blade normal vector and the relative flow direction. This consideration helps account for theradialcomponent’seffectofthebladenormalvector,essentiallytheleaningoftheblades. The sign function is utilized to ensure that the normal force always points in the direction that reduces thelocalflowdeviation. 34 Losses f v Body-forces representing pitch averaged losses Pitch Average Figure2.6: Illustrationoftheviscousbody-forcemodel. 2.4 ViscousLosses As introduced in Section 2.1, the viscous body-force, f v , is a dissipative force responsible for entropygenerationsimulatingvariouslossmechanisms. Theforcewasderivedtobeintheopposite directiontothatofthebladerelativeflowdirection. 2.4.1 GeneralForm Following a similar approach to Hall’s inviscid loading force formulation, the viscous body-force field is obtained from distributing the dissipation losses that occurs in the vicinity of the blade surfacesacrossthepitchdistance,likeFigure2.6illustrates. Theforcemagnitudeisfoundvia f v = K v ( 1 2 W 2 ) (2p r/B)b|n q | , (2.33) which represents the loss generated from a single blade after evenly distributed across the blade pitch. K v represents a local dissipation factor which in general depends on the local blade loading andboundarylayercondition K v = f(d ,Re q ) . (2.34) It was shown in Section 2.1 that entropy generation per unit mass from viscous body-force in steadyflowis ˙ s b =W ·— s= Wf v T . (2.35) 35 Thelocalentropygenerationperunitmassfromtheviscousbody-forcecanthenbeexpressedas ˙ s b = K v ( 1 2 W 3 ) (2p r/B)b|n q |T (2.36) K v can be linked to various loss mechanisms as will be introduced in the following sections. If theentropygenerationviacertainlossmechanism,D s spec ,isknownfromempiricalsources,using Equation (2.36), the resulted viscous body-force can be solved by evenly distributing the entropy generationacrossthechordas f v spec = TD s spec l m ||V m || ||W|| . (2.37) wherel m ismeridionaldistanceacrosswhichthelossisspread. Asintroducedin[52],lossmech- anisms in flow machines are quite complicated. To keep simplicity and generality, the viscous lossesaremodeledinonlytwopartsinthiswork,namely, K v ⌘ K v o + K od , (2.38) in which K v o accounts for profile losses without any pressure gradient effects, and are due to dissipation in boundary layer regions on both sides of the fictitious propulsor blade. To some degree,thislossrepresentstheon-designlosses. TheamountofK v o dissipationisbasedonknown flat plate correlations for laminar and turbulent flow. K od accounts for the extra off-design losses estimated based on a second order loss bucket model. We recognize that there are other loss sources that might be important for certain applications, and can be implemented to help achieve betterpredictions. 36 2.4.2 ReferenceProfileLosses Denton [52] provided an in-depth review of loss mechanisms in turbomachines where he men- tionedthatrateofentropygenerationfromaboundarylayerpersurfacearea, ˙ S a ⌘ d dx Z d 0 (r V x (s s d ))dy (2.39) canberelatedtoadissipationcoefficientC D via ˙ S a =C D r V 3 d /T , (2.40) in whichV d is the boundary layer edge velocity andC D is a dissipation factor that is function of momentumthicknessReynoldsnumber,Re q ,asplottedinFigure2.7. Entropygenerationspecified fromEquation(2.40)canberelatedtoK v o via ˙ S a =r ˙ s b (2p r/B)b|n q |= K v o ( 1 2 r W 3 ) T , (2.41) andresulting, K v o = C D (V 3 d u +V 3 d l ) 1 2 W 3 (2.42) IfweassumethelocalrelativevelocityWistheaveragedvelocitybetweenupperandlowerblade surfacethentheboundarylayeredgevelocityontheupperandlowersurfacescanbeassumedas V d u =W +D V (2.43) and V d l =W D V (2.44) 37 respectively. Assuming incompressible flow,D V can then be related to the blade inviscid loading force f n viaBernoulli’sequation,resulting D V = f n 2W (2.45) CombiningEquation(2.45),(2.43),(2.44)and(2.42),thereferencelosscoefficientcanbederived asfollows K v o = 4.0C D + 3.0C D (K c pd ) 2 , (2.46) Asimilarformoflosscoefficientisadoptedin[53]withaconstantC D . Alternatively,asimilarderivationcanbecarriedoutusingthefrictioncoefficientC f insteadof the dissipation coefficientC D . However, aC D based formulation was chosen in this work due to that(1)C D islesssensitivetothestateoftheboundarylayerthanC f ,and(2)wecantakeaccount ofthelocalV 3 d /V 3 • ratiowhichleadstomoreaccurateprofiledragapproximation[54]. Following the formulation in [52], the dissipation coefficientC D is taken as a function of the momentum thickness Reynolds number, Re q , as plotted in Figure 2.7. The momentum thickness is estimated locally using classical Blasius and turbulent 1/7 power boundary layer expressions. Thus,forlaminarflow,thedissipationcoefficientcanbecomputedas q = 0.664x (Re x ) 1/2 C D =0.173Re 1 q (2.47) andforturbulentflowas q = 0.016x (Re x ) 1/7 C D =0.0056Re 1/6 q (2.48) 38 Figure 2.7: Boundary layer dissipation coefficient versus Reynolds number based on momentum thicknessforlaminarflowandturbulentflow. Figureadoptedfrom[52]. whereRe x isevaluatedbasedonthelocalrelativevelocityW,localdensityr ,andlocaldistancex. xiszeroattheleadingedgelocationandequalstothechordattrailingedgelocation. Localxvalue isrequiredasaninputtothemodel. AcriticalReynoldsnumberneedstobeassumedforswitching betweenEquation(2.47)and(2.48). Alternatively, theflowcanbeassumedtobeentirelylaminar or turbulent in which case only one formulation will be applied. Additional 20% of profile losses is included to account for the trailing edge mixing losses [52], resulting the following adjustment totheon-designlosscoefficientK v o =1.2K v o . 2.4.3 Off-DesignLosses In traditional throughflow models, the losses are typically modeled via emperical correlations in form of an on-design reference value plus an off-deisign contribution [39]. As shown in Equa- tion2.49,theoff-designlossesareoftenexpressedwithasimplequadraticlossbucketasfunction ofincidenceangle, ¯ w = ¯ w ref +C(i i ref ) 2 , (2.49) 39 inwhichthei ref representsanon-designreferenceincidenceangleconstantCcanbefoundviacor- relations[55]. Inthiswork,thereferenceprofilelosscoefficientK v o derivedpreviouslyrepresents thew ref partofthelossesandtheoff-designcontributionfollowsaquadraticfunction, K od = K c K v 1 (d d ref ) 2 , (2.50) where d ref is a reference deviation angle distribution analogous to i ref in Equation (2.49), d is the local deviation angle, K c is the compressibility correction factor included here to account for the Mach number related effects and K v 1 is a constant that controls the quadratic curves shape similar to the function of constant C in Equation (2.49). The value of K v 1 varies with different types of applications experiencing different flow conditions. In this work, the K v 1 value is chosen to fit the given performance trends, with magnitude ranging from 5.0 for a lightly loaded fan rotor to 12.0 foraheavilyloadedhigh-speedrotor. 2.4.4 ShockLosses For high-speed turbomachines, inlet relative mach numbers are often greater than unity, resulting intheemergenceofshockwavesthatincurlosses. Asshownin[52],entropyriseperunitmassof theflowacrossweaknormalshockscanbeapproximatedas D s shock ⇡ c v 2g (g 1) 3(g +1) 2 (M 2 1) 3 , (2.51) whereM isMachnumberupstreamoftheshock,g istheheatcapacityratio,andc v istheconstant volume heat capacity. The rate of entropy generation via viscous body-force as shown in Equa- tion (2.7) can be used in combined with Equation (2.51) to derive a form of viscous body-force generatingshocklosses f v,shock ⌘ Tc v c 2g (g 1) 3(g +1) 2 (max(M r ,1) 2 1 3 , (2.52) 40 whereT islocaltemperature,cisthechordatcertainspanwiselocationandM r isthelocalrelative Mach number . This form of body-force spreads the losses generated by the shock across the entirechord. Under-predictionofshocklossispossiblesincediffusionhappenswhichreducesthe relative Mach number along the streamline, which in turn leads to an under-prediction of shock losses. 41 Chapter3 Methodology 3.1 ComputationalMethod Theopen-sourceCFDcodeADflow,amulti-blockstructuredandoversetfinitevolumesolverde- veloped at the University of Michigan [56], is used to carry out compressible Reynolds-Averaged Navier-Stokes (RANS) simulations with the Spalart-Allmaras turbulence model. The code uses second-order central difference discretization with matrix artificial dissipation. We use the fully implicitapproximateNewton-Krylov(ANK)method[57]whichreliesonbackwardEuler’smethod to iteratively solve the system of equations. This allows the use of very large CFL numbers, an efficient startup for overset meshes, and thus faster convergence. The body-force model is imple- mentedinADflowcodeasasourcetermroutine. All meshes are constructed using the commercial meshing software Pointwise. For blown-lift wing test cases, we took advantage of the ADflow’s overset grid capability to mesh around the complex geometry of the flapped wing and motor pylon, and added a dedicated region for the propulsormodel. EachtestcasesolutionwascomputedonsupercomputersattheCenterforAdvancedResearch Computing (CARC) of the University of Southern California. Convergence reduction of 8 orders in total residual were generally achieved. Except for cases with more inherent unsteadiness, con- vergencemaystallearlier,inwhichcasecarewastakentomakesurethatthequantitiesofinterest, suchaslift,dragorthrust,hadconverged. 42 3.2 ModelImplementation Theimplementationofthebody-forcemodelconsistsofthreemajorsteps,namelytheextractionof geometrical information from a full blade geometry, assigning the blade geometry information to eachofthecomputingcellsinsidethemodelsourceregion,andthecomputationofthelocalbody- forcescontributionatruntime. AsimilarprocedureusedforHall’sinviscidmodelwasintroduced indetailbyDogusandShishier[44]andusedasareferenceforthecurrentimplementation. Geometry extraction (Python + Rhino) Data preparation (Python) Assign geo info to computational grid (Solver pre-processing, ADflow) ✓ x camber surface normal vectors Axisymmetric Volume Computational grid Model region Source evaluation in runtime (ADflow) Solution pressure contour Body-Force model Solver Input r y z x Figure3.1: Modelimplementationworkflow 3.2.1 GeometryExtraction To provide geometry input for the body-force model, the camber surface normals and thickness distributionareextractedfromtheactualpropulsorbladegeometry. Inthiswork,commercialCAD software Rhinoceros was used to construct the blade surface and to extract geometry information from the surfaces, due to its good capability to handle complex surfaces and to be scripted with python. ThestartingpointisusuallyaCADmodelofthebladesurfaceseitherconstructedwithknown manufacture coordinates or 3D scanned from the actual blade geometry. We extract blade infor- mationusingastructuredgridofpointsfollowingtheU,Vcoordinatesofthecambersurface. The extraction points are clustered toward leading edges and tailing edges of the blade sections where geometry changes more drastically. The resulted geometry info are listed below and pictorially showninFigure3.2. 43 Camberpointcoordinates ~ x Surfacenormalvectors ˆ n(~ x) Localthickness(normaltocamber) t(~ x) Localmetalangle k (~ x) Chord c(r) U V ˆ n ~ x x y z c t x r flatten ~ x Figure3.2: Pictorialdemonstrationofgeometryextraction The next step is to flatten the camber surface and transform form the data collected in the last step into an axisymmetric format, which allows for easier direct assignment or more accurate interpolationtocomputationalgrids. TheCartesiancoordinatesofthecamberpoints,~ x=(x,y,z), arefirstconvertedtoacylindricalcoordinatesystemthatisalignedwiththepropulsoraxis(x-axis in this work) resulting in cylindrical coordinates,~ x=(x,r,q ). The camber surface is flatten by settingq =0 and only keeping meridional coordinates,~ x=(x,r). The vector geometry variables associated with each points, specifically the surface normal vectors, need to be rotated by q so thatthecorrectsurfacenormalsareretainedafterflattening. Theoutputfromthisstepistheblade geometry information stored as 2D structured-grid data. Additional data requried are evaluated includingthegradientofblockagefactor. Theresultedgeometryvariablesarethen ˆ n(x,r),b(x,r),— b(x,r),k (x,r),c(r) , (3.1) asplottedinFigure3.3. 44 Figure3.3: ContourplotsshowingdistributionofrelevantbladegeometricalpropertiesforNASA Rotor 67: Blade metal angle k , blade thickness t, blade blockage factor b and axial gradient of bladeblockagefactor— b. We recognizes that this geometry extraction step is the most laborious procedure for applying themodeltoanewtestcase. Inthefuture,workcanbedonetodevelopanalternativeroutinethat generates the blade geometry following ideal load distributions with only simple inputs such as propulsordimension,requiredthrust,androtationalspeeds. Thiswillallowmuchfasterandeasier applicationofthebody-forcemodelwhenthepropulsorgeometryisnotknown. 3.2.2 DataAssignmenttoComputationalGrid The geometrical data prepared in the previous step are used as input volumetric data to the solver aspartofapre-processingstep. Thecomputationalcellslyingwithinthemodelregionarelabeled as source cells. This is achieved by utilizing a geometric searching routine that already exists in the solver for the purpose of grid-to-grid interpolation. The searching algorithm takes the cell center of each grid cell and determines if this cell lies within the source volume. If it does, the 45 algorithm interpolates geometrical data from the input volume data to the current cell and stores thedatawithinadatamoduleforthebody-forcemodel. Thisstephelpsavoidtheneedtore-extract geometricaldataeachtimethecomputationalgridischanged. Thesamegeometricalinputcanbe usedfordifferentgrids,improvingflexibilityofmodelimplementation. 3.2.3 SourceTermEvaluation The body-force model is implemented in ADflow as a source term routine. Source contribute in each source cell are computed in runtime utilizing the input geometrical data and following formulation introudced in Chapter 2. Assume that surface normal points to the suction side of the blade camber. Given the following known info: current cell center~ x, propulsor axial vector ~ x a , axis origin ~ x 0 , rotational speedW , number of blades B. The following computations are repeated for each cell in the source volume. First, we need to find the coordinate in the cylindrical system centeredalongthepropulsoraxis. Radialcoordinatesfromthepropulsoraxisiscomputedas r = ||~ x a ⇥ (~ x ~ x 0 )||. (3.2) Unitvectorofthecircumferentialdirectioniscomputedas ~ t= ~ x a ⇥ (~ x ~ x 0 )/r (3.3) Rotationalvelocityvectoristhencomputedusinginputrotatonalspeed,W . U=W r ~ t (3.4) KnowingbladevelocityUandacquiringVfromthesolver,relativevelocitycanbecomputedas W=V U (3.5) 46 Surfacenormalvectorisprojectedontothetangentialdirection ~ n q =~ n· ~ t (3.6) Nextwefindtherelativevelocity’scomponentinthenormaldirection W n =W·ˆ n (3.7) Themagnitudeofdeviationanglecanbefoundby d =sin 1 (|W n |/||W||) (3.8) The body-force magnitude and directions can then be evaluated following formulation introduced in Chapter 2. The body-force contribution is added to mass, momentum and energy residual for celli: R m,i =R m,i +r i f i V i , (3.9) and R e,i =R e,i +r i (f i ·U i )V i . (3.10) whereRistheresidualvectorforeachequationinthegoveningequation,andV isthecellvolume. 47 Chapter4 ModelAssessmentandValidationinUniformFlow Thischapterfocusesonassessingtheproposedmodelforitscapabilitytoaccuratelyrepresentthe performance of stand-alone propulsors with uniform inflow condition. Test cases are shown for bothhigh-solidityfanstages,lowsoliditylowspeedfanstages. Investigationsonphysicalaspects ofthemodelsuchasflowturning,lossgeneration,bladesolidityeffectsarealsoshown. 4.1 OverviewofTestCases A variety of propulsors, spanning a wide range of sizes, numbers of blade and flow regimes, are selectedastestcasesinordertothoroughlyvalidatetheincludedphysicsandlimitsoftheproposed body-force model. As shown in Figure 4.1, the selected test cases mainly fall into two categories: propellers and turbomachinery blade rows. While all the test cases shown in Chapter 4 are under the turbomachinery category, the propeller test cases are presented seprately in Chapter 6. The adoptedpropellersareallAPCpropellersdesignedforR/Cairplanesanddrones. Propellergeom- etry information are found from APC website [58] and wind tunnel test results for a list of APC propellersusedforvalidationwerefoundin[59]. Twoturbomachineryapplicationswereadopted: NASA rotor 67 and NASA stage 35. Both rotors operate with transonic tip speeds at design con- dition. Their geometries and test data were found in several NASA reports [60, 61]. TF8000 fan stageisalow-speedR/Cductedfanthathasfewerbladesthantypicalturbomachinery. Italsohas more blades with lower aspect ratio than most propellers. It was used as model propulsors in the 48 powered D8 wind tunnel tests [4, 62, 63]. The geometry was acquired by directly scanning the rotorandstatorblades. Amongallthetestcases,thebladecountsrangefromtheminimumnumberof2bladesforpro- pellers to 36 blades in the case of turbomachinery. The tip mach number covers both incompress- ible and compressible ranges, and both subsonic and supersonic speeds. The propulsor diameter range from 0.127m (5in) as seen in propeller cases to 0.5m as observed in turbomachinery cases. Two of the test cases, NASA stage 35 and TF8000 fan stage, include a full stage, i.e., both a rotor andastator. DetailedparametercomparisonaresummarizedinTable4.1. 5 10 15 20 25 30 35 Number of blades, B 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 Diameter [m] NASA Rotor 67 NASA Rotor 35 Turbomachinery RC Propeller/Ducted fan TF8000 APC5x4-4 APC10x5E APC8x7 APC7x6 APC19x12E Figure4.1: Distributionoftestcases: diameterversusnumberofblades. 49 Table4.1: Summaryoftestcasesglobalparameters. TestCase B rotor (B stator ) D[m] AR W M tip Type NASARotor67 22 0.510 1.56 16040 1.38 ducted NASAStage35 36(46) 0.487 1.19 17220 1.33 ducted TF8000FanStage 5(4) 0.144 1.51 13500 0.3 ducted APC19x12E 2 0.4826 3000 open APC14x13S 2 0.3556 3000 open APC8x7S 2 0.2032 6000 open APC5x4E-4 4 0.127 10000 open 50 4.2 TransonicRotor: NASARotor67 NASA Rotor 67 is a 1.56 aspect-ratio transonic rotor with a tip relative Mach number of 1.38 and pressureratioof1.63atdesignmass flowof33.25kg/sanddesignrotationalspeedof16043rpm. Thisrotorisoftenusedasvalidationcaseforturbomachinerycomputationalworks. Ameridional view with the rotor geometry and axial measurements is shown in Figure 4.3. A detailed descrip- tionoftherotordesignandexperimentalmeasurementsatdesignrotationalspeedwasreportedby Strazisar[60],whoprovidedspan-wisemeasurementsattwooperatingpoints: peakefficiencyand near stall. Detailed design information was found in [64] and measured performance at different rotational speeds was found in [65]. The flow and rotor performance quantities are normalized usingthevaluesfrom[60],namelystandardsealevelpressureandtemperature: p ref =101325Pa, T ref =288.15K. 0.00 0.02 0.04 0.06 0.08 x[m] 0.00 0.02 0.04 0.06 0.08 y[m] Figure4.2: 3DRenderingofNASARotor67. 4.2.1 ComputationalSetup A meridional view of the computational grid created with Pointiwise over the domain is shown in Figure4.3. Withtheuniformflowinletconditionsusedforthistestcase,thethroughflowremains axisymmetric, so a 24 degree wedge shaped flow domain with periodic boundaries, as shown in Figure4.3, isutilizedinsteadofafullrevolution. Meshconvergencestudywascarriedoutforthe 51 NASA Rotor 67 test case with four levels of systematic refinement and shown in Table 4.2. The quantities of interest including total pressure ratio and isentropic efficiencies are shown to reach convergencewithincreasingmeshsizes. Thedomainhas230⇥ 150⇥ 12cellsintheaxial, radial, and circumferential directions. The grid in the body-force volume region is shown in Figure 4.3: ithas100cellsinboththeaxialandradialdirections. r x centerline 1 2 y x z Figure4.3: MeridionalmeshNASARotor67. Theboundaryconditionattheinletendofthedomainsetsthetotalpressure,totaltemperature, and flow direction, while the static pressure boundary condition (or back pressure) is imposed at thedomainoutlet. Aradialpressuregradientispresentattheoutletduetotheswirlingflowcreated 52 Totalcell Ncellacrossmodelregion p h isen 12,480 10 1.546 0.887 37,765 25 1.651 0.938 102,245 100 1.674 0.944 809,796 200 1.674 0.944 Table4.2: Resultofmeshconvergencestudy bytherotor,andthereforethebackpressureateachoutletcomputationalcellisdeterminedbased onsimpleradialequilibrium,namely p j = p BC , if j =1 (hub) p j = p j 1 +D p j 1 , if j>1 (4.1) with D p j 1 = r j 1 V 2 q ,j 1 r j 1 (r j r j 1 ), where j isthecellindexintheradialdirectionand p BC istheimposedbackpressureatthehub. 4.2.2 Metrics For turbomachinery cases, mass flow rate, ˙ m, are commonly used to specify the operating point becausetherotoroftenoperatesincompressibleflow, wheredensityvariationsaresignificantand needstobetakenaccountof. IntheNASARotor67andNASAStage35testcases,themassflow isnormalizedbythechokingmassflowratereportedfromexperiments. Theloadingcanbecharacterizedbythetotalpressureratioacrossthepropulsor p = ¯ p t 2 /¯ p t 1 , (4.2) whichismoreconvenientwhenlargetotalpressureriselevelsareexpected. Theoverbar,¯ ·,denotes massaveragedquantities. 53 Rotor or stage efficiency is measured by the isentropic efficiency, defined as the ratio of isen- tropicenthalpyrisetoactualenthalpyrise. andcomputedas h is = ✓ ¯ p 2 ¯ p 1 ◆ g 1 g 1 T 2 T 1 1 . (4.3) The isentropic efficiency quantifies the efficiency of the diffusion process that happens through the rotor in the presence of thermodynamic irreversibilities, and is more appropriate for use in turbomachinery applications for which the outlet dynamic pressure does not amount to useful energy. EntropygenerationasameasureoflossesisdefinedviaGibbsequationas s s • = R g 1 ln ✓ T T • ◆ +R ln ✓ r • r ◆ . (4.4) wherethereferencestateistakentobesetbyfreestreamtemperatureanddensity,R =287J/kg·K isthegasconstantforstandardair,andg =1.4theairspecificheatratio. Notethatinthisequation onlyT denotesatemperature. 4.2.3 QualitativeObservations The meridional contour plots of various flow properties obtained with the body-force model for NASA rotor 67 at peak efficiency are shown in Figure 4.4. As the flow passes through the rotor blade(body-force)region,diffusionoccursalongthemeridionalstreamline,andstaticpressureand temperature rise gradually. Local minima of static pressure and temperature can be observed near the root and are the result the metal blockage. Total pressure and total temperature also increase indicating work done on the flow, and swirling momentum was created as can be seen in absolute flow angle contours. Relative flow angle reduces across the blade region indicating flow turning and diffusion. Entropy is generated due to the presence of end-wall boundary layer (viscous flow 54 along the tip and hub walls resolved by the RANS simulations), and also across the blade via the body-force loss models. The entropy generated by body-force shows an increasing trend towards the tip due to higher relative blade velocities and higher shock losses. Since the blade losses are modeled in a pitch averaged fashion, its magnitude appears smaller than that in the end wall regions. Flow property variations exhibit unique contour pattern that is solely resulted from the utilizationofbladecambershape. Thesmoothvariationoftheflowpropertiesthroughthebody-forceregiondemonstrateanother advantage of the analytical formulation over the flow tangency method mentioned in Section 2.3. With the latter, Sturmayr et al.[47] showed that the appearance of the body-force region abruptly forces the flow to be tangent to the blades, which leads to non-physical discontinuities in flow an- gles—andthusinotherflowproperties—attheleadingedgeofthebladeregion. Variousmethods wereproposedtoremovetheleadingedgediscontinuityincludingmodifyingthebladeshape[46]. However, with the analytical formulation adopted in this work, flow turning is created gradually andthusnoleadingedgediscontinuitiesoccur. Pressurecontoursfrombody-forcesimulationsus- ingouranalyticalformulationandtheflowtangencymethodarecomparedinFigure4.5,showing theleadingedgediscontinuityresultingfromthelatter. Theresultswiththeflowtangencymethod areobtainedbyupdatingtheinviscidbody-forceuntilthedeviationangleiszero. ContoursoftheflowpropertiesatthechokedconditionareshowninFig.4.6. Thepresenceof anaxisymmetricshockiscapturedbyourapproach,withtheshockpatternevidentinthepressure contours. The shock location is close to that observed by Sturmayr et al. [47], but its strengths appearstopredictalowerstrength. Figure4.10showsacomparisonbetweentheshockpatterncapturedbytheanalyticalformula- tionofthepresentworkandthatcapturedbytheflowtangencymethod. Thedifferenceisduetothe factthatthepresenceofashockdeflectstheflowdirection,asshownintheabsoluteswirlingangle ofFig.4.6. Theflowtangencymethodiscapableofreactingtoflowdirectionchangesquicklyand thusformalocalpeakforceattheshocklocationthatbringsbacktheflowtowardstheblade. The 55 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.4: Meridional contours of flow properties through NASA Rotor 67: (a) Normalized total temperature,(b)Entropyrise,(c)AbsoluteMachnumber,(d)RelativeMachnumber,(e)Absolute flow angle, (f) Relative flow angle, (g) Normalized total pressure, and (h) Normalized pressure, at peak efficient operating point obtained using the proposed body-force model. Thin black lines showwheretheleadingandtrailingedgesofthebladeswouldbelocated,andenclosethevolume wherebody-forcesareapplied. analyticalformulationappliedinthisstudycannotrespondtoflowdirectionchangesinstantly,and thussmears-outthediscontinuity. 56 0.779 0.806 0.832 0.832 0.859 0.859 0.885 0.885 0.912 0.938 0.965 0.991 1.071 1.203 1.256 1.282 0.779 0.806 0.832 0.859 0.885 0.885 0.912 0.912 0.938 0.938 0.965 0.965 0.991 0.991 0.991 1.018 1.018 1.044 1.044 1.071 1.097 1.124 1.124 1.150 1.150 1.176 1.203 1.229 1.256 1.282 1.309 1.335 1.362 1.388 1.388 1.415 1.415 1.441 1.441 1.468 1.468 1.494 1.521 1.547 Figure 4.5: Normalized pressure contours: qualitative comparison between (left) our analytical formulation,and(right)inviscidflowtangencymethod. 4.2.4 FlowTurning Radialdistributionofrelativeflowangleattheinletandoutletoftherotor,incidenceangle,camber angle, and deviation angle are plotted in Figure 4.8 at the peak efficiency operating point. The relativeflowangleatinlet,b 1 ,varieswithradius,asthisrotorisafirst-stagerotorthatingestsonly axial flow. On the other hand, the outlet flow angle, b 2 , varies more significantly with radius due to the varying stagger angle and camber angle of the blades. It should be noted that the net flow turning angle |b 2 b 1 | closely follows the camber angle distribution, with its magnitude slightly smallerduetodeviation,indicatingthatflowturningiscreatedbasedontheprovidedcamber. The resulted deviation angle is on the same order with incidence angle. Theses effects demonstrate that the body-force model produces flow turning effects representative of those would otherwise becreatedbythespecifiedrotor. ItwasshowninaninviscidcascadestudybyWeinig[66]thatasblade-to-bladespacingreduces inabladerow,theliftcoefficientforeachindividualbladereduceslinearlywithspace-chord-ratio andscaleswith1/cos(k ),resultinginWeinig’scoefficient, K ⌘ c ` c ` isolated ⇡ p 2 (s/c) 1 cos(k ) . (4.5) Inanefforttoverifythatthebody-forcemodelproducesthecorrecttrendofinviscidloadvariation with respect to blade spacing, a series of body-force test cases with increasing number of blades, from B = 2 to 50, were carried out using the NASA Rotor 67 geometry. The flow coefficients for these cases are kept the same at f =0.35, so that the incidence angle distribution is the same 57 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.6: Meridional contours of flow properties through NASA Rotor 67: (a) Normalized total temperature,(b)Entropyrise,(c)AbsoluteMachnumber,(d)RelativeMachnumber,(e)Absolute flow angle, (f) Relative flow angle, (g) Normalized total pressure, and (h) Normalized pressure, at choked condition obtained using the proposed body-force model. Thin black lines show where theleadingandtrailingedgesofthebladeswouldbelocated,andenclosethevolumewherebody- forcesareapplied. across all cases. To conform the potential flow assumptions, Hall’s original turning model was usedwithoutmodifications. TheWeinig’scoefficient,K,isapproximatedvia K⇡ R c 0 2pd ds/c 2pa = d ave i+J /2 , (4.6) 58 0.510 0.540 0.570 0.600 0.630 0.660 0.690 0.720 0.750 0.780 0.780 0.810 0.810 0.840 0.840 0.840 0.870 0.870 0.900 1.080 0.500 0.550 0.600 0.600 0.650 0.650 0.700 0.700 0.700 0.750 0.750 0.750 0.800 0.800 0.850 0.850 0.850 0.900 0.900 0.950 1.300 1.350 1.400 Figure4.7: Normalizedpressurecontoursshowingthecapturedshockpatternatchokedcondition: comparisonbetween(left)analyticalformulationand(right)flowtangencymethod. and plotted against the blade pitch-to-chord ratio, s/c, in Figure 4.9 with comparison to the theo- retical trends. Consistent with theory, the body-force model predicts a linear variation of inviscid loading with s/c for smaller spacings (s/c<1). As blade spacing increases, the inviscid loading trendflattensandapproachesaconstantvalue,whichisindependentofbladespacing. Thisbehav- ior can be attributed to fewer interactions between blades when the blade number is small. Thus, Hall’s turning model effectively captures the correct inviscid loading variation concerning blade spacing for both high-solidity and low-solidity propulsors. When comparing the trends at various spanwise locations with varying staggering angles, it becomes apparent that the load scaling with 1/cos(k ) is accurately represented for smaller spacings, approximately s/c<0.5. This is mostly 0 20 40 60 80 Angle [deg] 0.0 0.2 0.4 0.6 0.8 1.0 Normalized span incidence, i deviation, camber, ✓ flow turning, ✏ 1 2 Figure 4.8: Spanwise distribution of incidence angle (i), outlet deviation angle (d TE ), camber angle(J ),flowturningangle(e ⌘ b 1 b 2 ),rel- ative flow angle at inlet (b 1 ) and outlet (b 2 ), for NASARotor67atnearpeakefficiency(PE)op- eratingpoint. 0 2 4 6 8 s/c (2 r)/(Bc) 0.0 0.2 0.4 0.6 0.8 1.0 K ave /(i+ /2) 10%, = 20.2 , = 49.8 20%, = 27.1 , = 36.5 40%, = 39.9 , = 21.0 50%, = 45.1 , = 15.7 60%, = 49.2 , = 12.1 80%, = 55.6 , = 9.0 Theory =50 =0 B=50 B=22 B=12 B=8 B=4 B=2 B=30 # Figure4.9: Weinig’scoefficientversuspitch-to- chord ratio at spanwise location varying from 10% to 80%, for NASA Rotor 67 with varying solidity (B = 2 to 50), with comparison to theo- reticaltrends. 59 attributed to the incorporation of |n q | in the denominator of f n . However, the body-force model’s trends deviate from theoretical trends for approximately s/c > 0.5, specifically under-predicting therateofincreaseofinviscidloadingwithspacing. Thisdeviationresultsfromtheturningmodel’s failure to account for bound circulation effects, leading to an under-prediction of local deviation angles. Incorporating the flow angle correction as specified in Equation (2.31) is anticipated to compensate for this under-prediction, which also justifies the necessity of a flow angle correction intheinviscidloadingformulation. 4.2.5 MetalBlockage Pressure contour at the choked condition are shown in Figure 4.10 with comparison to that pro- duced using the flow tangency method. The presence of shock is captured by our approach, with theshockpatternevidentinthepressurecontoursrepresentingapitch-averagedblade-to-bladenor- mal shock. The shock location is close to that observed by Sturmayr et al. [47], but its strengths appears to predict a lower strength. Figure 4.10 shows a comparison between the shock pattern captured by the analytical formulation of the present work and that captured by the flow tangency method. Thedifferenceisduetothefactthatthepresenceofanaxisymmetricshockreducesaxial 0.510 0.540 0.570 0.600 0.630 0.660 0.690 0.720 0.750 0.780 0.780 0.810 0.810 0.840 0.840 0.840 0.870 0.870 0.900 1.080 0.500 0.550 0.600 0.600 0.650 0.650 0.700 0.700 0.700 0.750 0.750 0.750 0.800 0.800 0.850 0.850 0.850 0.900 0.900 0.950 1.300 1.350 1.400 Figure 4.10: Normalized pressure con- tours showing the captured shock pattern at choked condition: comparison between (left)analyticalformulationand(right)flow tangencymethod. 0.90 0.95 1.00 1.05 Normalized mass flow, ˙ m/ ˙ m ref 1.4 1.5 1.6 1.7 1.8 Total pressure ratio, Body-force Hall Exp Fidalgo Exp Strazisar Hall with blockage Turning only Hall with blockage With blockage and on-design losses Figure 4.11: Total pressure ratio versus mass flow at design rotational speed: comparison be- tween cases with and without blockage model applied. 60 velocity but do not affect circumferential velocity thus deflects the flow direction. The flow tan- gency method is capable of reacting to flow direction changes quickly and thus form a local peak forceattheshocklocationthatbringsbacktheflowtowardstheblade. Theanalyticalformulation applied in this study cannot respond to flow direction changes instantly, and thus smears-out the discontinuity. Todemonstratetheeffectoftheblockagemodelingonthepropulsorperformance,theproposed model is compared with the original inviscid turning model by Hall in Figure 4.11. Comparisons are shown between (i) Hall’s original model, (ii) Hall’s original model with added blockage, (iii) the proposed inviscid loading model, (iv) the inviscid loading model with blockage, and (v) the on-design model (turning, blockage, on-design losses, but no off-design losses). It is evident that the addition of blockage effectively reduces the choking mass flow rate and the work done by the body-forcemodel,aswellasproducingthecorrectchokingbehaviorcharacterizedbyasharpdrop of pressure rise near choking. The reduction of work input is expected as work input is depen- dent on mass flow. Without blockage, Hall’s original turning model over-predicts the pressure ratioby5%andthechokingmassflowbyover10%,failingtoproducephysicallycorrectchoking behavior. However, applying blockage without further modification causes the original model to under-predict the pressure ratio by up to 9%. By employing the proposed inviscid formulation, whichtakesintoaccountcompressibilitycorrection,flowanglecorrection,andtheblockagefactor introduced in the f n denominator, notable improvements in the predictions are observed. Specifi- cally, the total pressure ratio is predicted with an overestimation of approximately 5%, while the physicallyaccuratechokingbehavioriscapturedwithachokingmassflowrateconsiderablycloser to the experimental value. Incorporating the on-design loss model further decreases the pressure ratioandchokingmassflow. Strazisarmeasuredachokingmassflowof34.9kg/s,whilethepro- posedon-designbody-forcemodelpredicts34.6kg/s,resultinginachokingpointcapturedwitha 1%error. Aspreviouslydiscussed,onedisadvantageoftheblockagemodelisthataddingasource termtothemassequationcanleadtoerrorsinmassconservation. Wesawonaveragea0.08%loss inmassflowacrossthebody-forcevolumewhenusingthemodelwithblockage. 61 4.2.6 Losses The impact of incorporating losses on the predicted performance, including pressure rise and effi- ciency,isdemonstratedinFigure4.12. Thisfiguredisplayscurvesforresultswithonlyon-design losses, results with all losses, as well as predictions from the inviscid loading model without loss application. It is observed that, without loss implementation, the isentropic efficiency resulted fromtheinviscidloadingmodelisnearly100%,withminimallossesattributedtoend-wallbound- ary layers and a higher amount of losses generated by shock near choking, as discussed in Sec- tion 4.2.5. The adoption of on-design losses, in conjunction with inviscid loading and blockage, effectively reduces the previously overestimated total pressure ratio, resulting in better agreement with experimental data at the on-design operating point. This validates that the correct amount of inviscid loading and losses are predicted by the body-force model for the on-design condition. However, as anticipated, the off-design behavior is inadequately captured when solely including theon-designcomponentofthelossmodel,leadingtoanover-predictionofboththetotalpressure ratio and isentropic efficiency under off-design conditions. Upon applying off-design losses, with a K v1 of 10.0 and d ref extracted at the design operating point, the off-design trends in total pres- sure ratio and isentropic efficiency are accurately captured. It should be noted that the off-design prediction is sensitive to the choice of operating point to extract d ref and the choice of off-design loss parameter K v1 . Ideally, d ref should be extracted at an operating point representative of the 0.88 0.90 0.92 0.94 0.96 0.98 1.00 Normalized mass flow, ˙ m/ ˙ m ref 1.4 1.5 1.6 1.7 1.8 Total pressure ratio, Body-force Hall Exp Fidalgo Exp Strazisar No loss On-design loss only All losses Near stall Peak efficiency (Design) 0.88 0.90 0.92 0.94 0.96 0.98 1.00 Normalized mass flow, 0.80 0.85 0.90 0.95 1.00 Isentropic efficiency, isen Body-force Hall Exp Fidalgo Exp Strazisar No loss On-design loss only All losses Near stall Peak efficiency (Design) ˙ m/ ˙ m ref Figure4.12: Totalpressureratioandisentropicefficiencyversusmassflow: comparisonofresults withandwithoutlossmodels. 62 PeakEfficiencyOperation NearStallOperation ˙ m[kg/s] FPR h is ˙ m[kg/s] FPR h is Experiments 34.57 1.64 0.93 32.31 1.73 0.90 Body-Force 34.53 1.62 0.91 32.42 1.73 0.90 Table 4.3: Performance parameters at peak efficiency and near stall operating points: mass flow rate ˙ m,mass-averagedtotalpressureratioFPR,andisentropicefficiencyh is . Comparisonbetween experimentalmeasurementsandbody-forcepredictions. minimum loss condition, which can be inferred either by the diffusion number or the incidence angle. Furthermore, while it is possible to introduce physical models to give the K v1 coefficient a more physical interpretation, such as the influence from streamwise adverse pressure gradients, it is still treated as a turning constant in this work. Finally, It is noteworthy that losses generated by shock at the choking condition (right-most points in the curves of Figure 4.12 are excessively high;althoughapressureratiosimilartoexperimentalmeasurementsispredicted,theefficiencyis underestimatedby10%atthechoking(maximummassflowpoints). 4.2.7 OverallPerformance Tovalidatethemodel’slocalpredictions,thespanwisedistributionofflowpropertiesiscompared withexperimentalmeasurementsattwooperatingpoints: peakefficiency(PE)andnearstall(NS). Theseoperatingpointsaresetinthesimulationsbymatchingtheinlettotalconditionandadjusting theoutletboundaryconditionstocloselymatchtheexperimentaldataatstation1(seeFigure4.13). Forthenearstalloperation,theupstreamaxialvelocityandpressureprofileatstation1exhibitgood agreementwithexperimentaldata,exceptnearthetip. However,weobservedlargerdiscrepancies intheupstreamprofilesforthepeakefficiencyoperatingpoint,likelyduetotheeffectofblockage modelbeingstrongeratthepeakefficientoperatingpointwhichisclosetothechokingpoint. The blockageeffectisstrongerneartheroot,asaresulttheaxialvelocityprofileshowslowervaluenear the root. Discrepancies in upstream conditions implies differences in local operating points and eventually leads to differences in the predicted performance downstream. As we can observe that downstream flow properties for the near stall operating condition shows overall good agreement 63 Figure 4.13: Spanwise distribution of flow properties in peak efficiency (PE) and near stall (NS) operation. with experiments, whereas larger deviation from the experiments are seen for the peak efficiency operation results. However, these small local differences don’t have significant impact on overall performance. The good agreement at the near stall operating point validates that the local flow turning and work input are well predicted by the body-force model, given the correct upstream 64 condition. Conversely,theover-predictionoftotalpressureandtemperatureforthepeakefficiency operatingpointcanbeattributedtotheunder-predictionofupstreamflowcoefficients,resultingin higherincidenceangles. Across 40% of the span in the tip region where the blade operates with supersonic flow, large discrepanciesareobservedforbothoperatingpoints,suggestingthatsomeloss-inducingflowfea- turesarenotadequatelycapturedbythebody-forcemodel. Althoughshocklossesareincludedin thelossmodeling,theresultingeffectontheflowfielddiffersfromthatproducedbyactualoblique shocks, specifically the shock compression and heating effects are not captured. This explains the under-prediction of static pressure rise and total temperature ratio in the tip region. Additionally, an actual shock produces a discontinuity that reduces flow speed instantly, allowing for higher flow turning, as indicated in the absolute flow angle from the experiment. The body-force model distributestheloss,thusitisunabletoproducethesamediscontinuityeffect. Local over-prediction of total pressure, total temperature, and absolute flow angle can be ob- served near the end-walls where the flow is significantly off-design. This suggests the presence ofanon-physicaldiffusioninthebody-forceresults,resultingfromtheover-turningoftheslower flow near the boundary-layer, which contributes to the over-prediction of static pressure. Tradi- tional through-flow models address this diffusion problem by adjusting the deviation angle based on correlations. However, the analytical loading model employed in this work does not directly controlthedeviationangle. Analternativeapproachtopreventnon-physicalflowturninginvolvesapplyingaload-limiting factor to the inviscid body-force. For example, Gallimore[38] introduced a correction procedure to account for end-wall loading limits in a viscous through-flow method. Although this method imposes a correction on an integrated blade force, it is incompatible with the local nature of our body-force model. A localized limiting procedure is required for the current loading model and willbeconsideredinfuturework. 65 OverallperformanceatthetwooperatingpointsarecomparedwiththeexperimentinTable4.3 which shows that the body-force model captures the peak efficiency performance within 2% error andthenearchokeperformancewithin1%error. Varying the rotational speeds as an input to the body-force model influences the evaluation of relative velocity W and, in turn affects both the force magnitude and work input. The predicted performance maps at different rotational speeds are shown in Figure 4.14, and are compared with the experimental measurements and the CFD data reported by Fidalgo [5]. The overall agreement oftotalpressureratiowithexperimentsforallrotationalspeedsindicatesthatthebody-forcemodel effectivelycapturesthevariationofworkinputduetochangingoperatingspeeds. Aroughly1.5% over-predictionoftotalpressureratioisobservedatrotationalspeedslowerthandesign. Efficiency trends are well predicted across all speeds and align closely with full-wheel CFD results, except at near choking operating conditions (higher mass flow), where the blockage model introduces excessive losses that cause the efficiency to under-predict compared to full-wheel CFD. When compared with experiments, the efficiency under-predicts by as much as 3% for lower RPMs, a result that is consistent with full-wheel CFD predictions. The aforementioned discrepancies in performance predictions highlight the primary limitation of the body-force model. Physical effects are introduced through simple physics models, which cannot guarantee an exact match withactualperformancesincetherewillalwaysbeagapbetweenphysicalmodelingandtheactual flowbehavior. Nevertheless,thecurrentmodelcontainsallessentialcomponentsthatcanbefurther improvedbyintroducingmoredetailedmodelingormanuallycalibratingthemodeltobettermatch knownperformanceiftheexactreproductionofaspecificcharacteristicisdesired. 4.2.8 DeterminationofReferenceDeviationAngle Sinced ref isalocalizedvariableanditsdistributioncanonlybeobtainedduringsimulationruntime, in this work, we acquire the d ref distribution by extracting the deviation angle distribution from a test case (without the application of off-design loss) at a reference operating point, which is typicallyclosetothedesignpointorminimumlossoperatingpoint. 66 Figure 4.14: NASA Rotor 67 performance map, with comparision between prediction from the proposedbody-forcemodel,experimentalmeasurementsandCFDresultsreportedbyFidalgo[5]. However,inmanycases,thedesignoperatingpointisunknownandhastobeinferredthrough certainflowproperties. Thereareprimarilytwoindicatorsthatcanbeusedtodeterminetherefer- enceoperatingpoint-thediffusionfactorandtheincidenceangle. Inturbomachinery,lossgenera- tionistypicallycorrelatedwiththediffusionfactor[67]. Therefore,theoperatingpointassociated withtheminimumdiffusionfactorcanbeusedasthereferencepointfromwhichthed ref profileis extracted. The diffusion factor as a function of flow coefficients for Rotor 67 at different RPM, as plotted in Figure 4.15, indicates that the minimum diffusion factor occurs at a specific mass flow rate. Alternatively,theoperatingpointassociatedwiththeminimumlossincidenceanglefoundin empiricalcorrelations[68]mayalsobeusedasthereferencepoint. 0.90 0.92 0.94 0.96 0.98 ˙ m/ ˙ m ref 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 D Figure 4.15: Diffusion factor versus normalized mass flow for NASA Rotor 67 as represented via thebody-forcemodel. 67 4.3 TransonicStage: NASAStage35 NASAStage35ispresentedasthesecondtestcasewithafullstageconsistsofarotorandastator, thus can demonstrate the body-force model’s capability to represent a full stage within the same frame of reference. . It is one of the four high-load, high-speed axial flow fan stages presented by Reid and Moore [61]. As shown in Figure 4.16, it has a 1.19 aspect-ratio rotor with a tip relative Machnumberof1.33anddesignpressureratioof1.8atdesignrotationalspeedof17188RPMand mass flow of 20.188 kg/s. The rotor has 36 blades and the stator has 46 blade. Performance map at various rotational speeds are found in [61]. The rotor in this case has nearly twice the number ofbladesandloweraspectratiocomparedtoNASARotor67,resultinginahigherloading. Figure4.16: 3DRendering(left)androtorbladesections(right)ofNASAStage35. 4.3.1 ComputationalSetup A meridional view of the computational grid is shown in Figure 4.17. Total conditions are set at theinletandstaticpressureissetattheoutletwithsimpleradialequilibriumappliedsimilartothat shown in Equation 4.1. The same metrics as adopted in the Rotor 67 test case is used in this test case. 68 This test case demonstrates the body-force model’s capability to represent both a rotor and a statorwithoutusingdifferentframesofreference. Toapplybody-forcemodeltorepresentastator, wesimplysetthemodelinputrotationalspeedstozero. Figure4.17: MeridionalmeshNASAStage35. InFigure4.18,themeridionalcontourplotsofvariousflowpropertiesobtainedwiththebody- forcemodelforNASAstage35arepresented. Boththerotorandthestatorregionsexhibitarisein staticpressure. However,thetotalpressureriseisonlycreatedbytherotorandisreducedthrough the stator due to losses. In the rotor region, there is a total temperature rise, while there is no changeinthestatorregion,asthebody-forcehasnocontributiontotheenergyequationinthecase ofastator(f q W r =0). Theabsoluteswirlangleisobservedtobecreatedbytherotorandreduced bythestator. 4.3.2 OverallPerformance The prediction of overall performance by the body-force model at different rotational speeds is shown in Fig. 4.20 compared with experimental data presented in [61]. General good agreement with the experiment was achieved with the physically correct choking behavior produced. The predicted choking mass flow was captured within 1% of the experimental value. For this rotor at 100%designRPM,therelativeMachnumberreachesupto1.5nearthetip,thussignificantshock lossesareexpected. Experimentaldatarevealsthatthemaximumefficiencyislowerat100%RPM comparedto90%RPM.Byincorporatingtheshocklossmodel,thebody-forceisabletoaccurately reflectthisobservedtrend. SimilareffectsareobservedinNASARotor67resultsaswell. 69 (a) (b) (c) (d) (e) (f) (g) (h) Figure4.18: MeridionalcontoursofflowpropertiesthroughNASAStage35: (a)Normalizedtotal temperature,(b)Entropyrise,(c)AbsoluteMachnumber,(d)RelativeMachnumber,(e)Absolute flow angle, (f) Relative flow angle, (g) Normalized total pressure, and (h) Normalized pressure, at design operating point obtained using the proposed body-force model. Thin black lines show wheretheleadingandtrailingedgesofthebladeswouldbelocated,andenclosethevolumewhere body-forcesareapplied. 70 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.19: Meridional contours of flow properties through NASA Stage35 after choked: (a) Normalized total temperature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number, (e) Absolute flow angle, (f) Relative flow angle, (g) Normalized total pressure, and (h) Normalized pressure, at design operating point obtained using the proposed body-force model. Thin black lines show where the leading and trailing edges of the blades would be located, and enclosethevolumewherebody-forcesareapplied. 71 Figure 4.20: NASA S35 overall performance as predicted by the body-force model (solid line) withcomparisonwithexperimentalmeasurements[61](squaremarker). 72 4.4 Low-SpeedStage: TF8000FanStage TF8000 is a low-speed ducted fan with a tip diameter of 0.144m. It consists of a rotor and a stator. Therotorhadfiveblades,anaveragechordof0.039m,andahub-to-tipratioof0.178atthe leading edge and 0.383 at the trailing edge. The stator has four blades and a constant hub-to-tip ratio of 0.453. The operating wheel speeds range from 8000RPM to 13500RPM with a design flow coefficient of around 0.37. The meridional view of the geometry is shown in Figure 4.22. Wind tunnel experimental tests of the TF8000 fan stage with and without inlet distortion were carriedoutbySiu[62]andareusedforcomparisonandvalidationinthepresentwork. Figure 4.21: 3D rendering of TF8000 fan stage (left) and rotor blade sections at 10%, 50% and 100%spanwiselocations(right). 4.4.1 ComputationalSetup The computational domain extends 3.9 fan chords upstream of the fan face and 4.4 chords down- stream as shwon in Figure 4.22. For the simulations of the TF8000 in clean flow operation, the computationaldomainspans30 withperiodicboundaryconditions. 73 Figure4.22: MeridionalmeshTF8000fanstage. 4.4.2 Metrics Theoperatingpointofapropulsorisusuallysetbytheflowcoefficientdefinedas f ⌘ V x W R , (4.7) whereV x isfreestreamaxialvelocityandW iswheelspeed. Bladeloadingisquantifiedbythenon-dimensionaltotalpressurerisecoefficient y p t = ¯ p t 2 ¯ p t 1 r U 2 tip . (4.8) Aerodynamicefficiencyofthefanisdefinedas h ⌘ 1 P loss P tot (4.9) via the ratio of non-loss-generating power to total power produced by the propulsor model. Note that efficiency is only computed for the body-force model since its loss-generating force and non- loss-generatingforcearedefinedseparately. Thetotalshaftpowerisfoundbyintegratingthetotal body-forcecontributiontotheenergyequation, P S = ZZ Z V r f·UdV (4.10) 74 andtheloss-generatingpowerisfoundvia, P loss = ZZ Z V r f·WdV . (4.11) Pressureandtotalpressurearenon-dimensionalizedaspressurecoefficient C p ⌘ p p • 1 2 r V 2 • , (4.12) andtotalpressurecoefficient C p t ⌘ p t p t • 1 2 r V 2 • , (4.13) respectively. 4.4.3 EffectofFlowAngleCorrection Figure 4.23 shows meridional flow contours of total pressure coefficient and absolute swirl an- gles, comparing cases with and without the flow angle correctionD d . In the uncorrected case, a counter-swirlisgeneratedattherotorleadingedgeduetothenegativeincidenceangle,leadingto a reduction in total pressure. Total pressure starts increasing through the rotor when blade cam- ber takes effect further downstream. However, after applying the flow angle correction,D d , the negative incidence angle and the resulting counter-swirl are effectively reduced. This leads to the creation of more flow swirl and a larger total pressure rise by the rotor. Additionally, the stator is found to be more compatible as the stator inflow angle matches better with the stator blade angle, asindicatedbytheswirlanglecontourneardesignflowcoefficient. 4.4.4 OverallPerformance Figure4.25presentsthetotalpressureriseandfanefficiencycharacteristicsoftheTF8000fanstage aspredictedbythebody-forcemodel. Comparisonsaremadebetweencaseswithandwithoutflow 75 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.23: Meridional contours of flow properties through TF8000: (a) Normalized total tem- perature, (b) Entropy rise, (c) Absolute Mach number, (d) Relative Mach number, (e) Absolute flow angle, (f) Relative flow angle, (g) Normalized total pressure, and (h) Normalized pressure, at design operating point obtained using the proposed body-force model. Thin black lines show wheretheleadingandtrailingedgesofthebladeswouldbelocated,andenclosethevolumewhere body-forcesareapplied. anglecorrectionapplied, aswellaswithandwithoutoff-designlossapplied. Theproposedbody- force model, with blockage and losses applied but excludes flow angle correction, underestimates thetotalpressurerisebyroughly30%. However,uponapplyingflowanglecorrection,themodel’s predictions significantly improve and demonstrate good overall agreement with the experimental measurements of Siu [62] in terms of pressure rise characteristics. Due to fewer blades, this fan 76 Figure4.24: MeridionalcontoursofflowpropertiesthroughTF8000fanstage: totalpressurecoef- ficient(leftcolumn)andabsoluteswirlangle(rightcolumn),comparisonbetweencaseswithflow angle correction (top row) and without flow angle correction (bottom row) at f =0.37 obtained usingtheproposedbody-forcemodel. stagehaslowersoliditycomparedtopreviousturbomachinerycases,renderingtheboundcircula- tion correction more critical and explaining the substantial impact of flow angle correction. With considerablylessloadingandpressurerisethanpreviousturbomachinerycases,thefanbladessuf- fer less from off-design losses due to adverse pressure gradients and thus exhibits a linear trend in pressure rise characteristics as shown by the experimental data. As indicated in Figure 4.25, the body-force model was able to accurately captures the slope of this linear trend without the need to apply off-design loss model. Predicted efficiency without off-design loss applied aligns withexperimentalresultsonlyinthevicinityofthedesignoperatingpoint,butdeviatesunderoff- designconditions. Uponapplyingoff-designloss,withad ref extractedatflowcoefficientf =0.37 and K v1 = 5, the efficiency trend exhibits an overall closer agreement with experimental data at higher flow coefficients, while the linear trend in total pressure rise characteristic remains largely unaffected. 4.5 Summary Three test cases ranging from high-solidty high-speed turbomachinery (NASA Rotor 67, NASA Stage 35) and low-solidity low speed fan rotor (TF8000) are presented For high-solidity turbo- machinery cases, the blade metal blockage plays a crucial role in performance. The body-force 77 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.04 0.06 0.08 0.10 0.12 BFM Exp RPM8000 Exp RPM10600 Exp RPM12250 Exp RPM 13500 No flow angle correction With flow angle correction With off-design losses 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.80 0.85 0.90 0.95 No flow angle correction With flow angle correction With off-design losses Figure4.25: TF8000fanstageperformance: comparisonbetweenexperimentaldatafromSiu[62], body-force model results with and without flow angle correctionD d , with and without off-design losscomponent,atarotationalspeedof13500rpm. model, thanks to its blockage component, captures the physically correct choking behavior and predicts the choking mass flow within a 1% error compared to experimental measurements. The proposed model achieves reasonably accurate predictions near the on-design operating condition for both pressure rise and efficiency, indicating the correct amount of inviscid load and losses. To capture off-design trends due to increased losses, the off-design loss component must be applied, requiringareferencedesignoperatingpointtobechosen. Withallmodelcomponentsapplied,the totalpressureandefficiencyarewellpredictedforvariousrotationalspeeds. Theeffectofreduced 78 efficiency at higher RPMs due to higher tip Mach numbers was captured because of the inclusion ofashocklossformulation. InthecaseofTF8000fanstage,theflowanglecorrectionsignificantly improved prediction, fixing the previously under-predicted pressure rise. The linear trend of pres- sure rise characteristics was accurately captured without the need for off-design losses, whereas thecorrectefficiencytrendwasproducedonlyafterapplyingtheoff-designcomponent. 79 Chapter5 CapturingtheEffectofPropulsorInflowDistortion With inflow distortion, some portion of the propulsor operates under more unfavorable conditions thantheuniformflowcondition,resultinginvariationsinperformance. Theinfluenceonpropulsor performance due to various types of total pressure distortion as predicted by the proposed body- force model is investigated in this chapter. In the first two section, NASA rotor 67 test cases was set up with radial and circumferential total pressure distortion. The results produced by the body-force model are compared with that predicted from a bladed high-fidelity CFD study using the same rotor geometry. In the following section, the TF8000 fan stage is given a vertically stratified inlet profile that resembles a fuselage boundary layer profile as ingested by integrated propulsors. This distortion profile represents inlet distortion in both radial and circumferential directions. The body-force predictions are compared with experiments and the distortion transfer effectsarediscussed. 5.1 Introduction Inthischapter,wefocusonthreetypesofinlettotalpressuredistortions,namelyradialdistortion, circumferentialdistortion,andboundarylayeringestion(BLI)typedistortion,whichencompasses acombinationofbothradialandcircumferentialdistortions. Thesespecificvariationswerechosen to represent the diverse types of distortions a propulsor may encounter during operation, and to evaluate the body-force model’s prediction of distortion response under different flow conditions. 80 The chosen distortion length scales are much larger than the blade spacing, so that the modeled distortionrepresentlow-frequencysquarewaveratherthansmallperturbations. The effect of distorted inflow of various kind on turbomachinery performance has been ex- tensively studied in early works [69, 70, 71]. These studies often utilized simplified modeling techniques, such as the parallel compressor model, and focused on idealized distortion profiles, suchascircumferentialtotalpressuredistortion. AsillustratedbyFigure5.1,theparallelcompres- sormodelsuggeststhatthetwoinletstreamsoperatingwithdifferenttotalpressurecanbeviewed as two compressors operating in parallel at different operating points. Assumptions are made that there’s no cross-flow between the two streams and each stream operates with the same uniform flow propulsor characteristic. The propulsor performance under distortion can be implied by tak- ing the weighted average between the performance taken at the lower flow coefficient operating pointandhigherflowcoefficientoperatingpoint. Assumingthatthepropulsorperformancecurves are concave curves, losses in both pressure rise and efficiency can be resulted from this simple analysis. ˙ m p T t p t ⌘ TS Loss in⌘ 180 distortion 120 distortion No distortion Loss in Highp t Lowp t High total pressure Low total pressure Figure 5.1: Illustration of the parallel compressor model concept. Figure reproduced based on schematicsfrom[70]. However,it’simportanttonotethattheassumptionsoversimplifythetruebehaviorofdistorted flowwithinaturbomachinerystageandmaynotaccuratelycapturecomplexflowinteractions[71]. Specifically, the parallel compressor model does not take into account of the complicated flow 81 migration upstream of the rotor due to its response to flow distortion, which effectively changes the operating mass flow in both the clean and the distorted sectors. Experiments [72, 73, 74] and high-fidelityfull-annulusunsteadysimulationresults[75,5,73]revealedthatthefaninteractswith flow distortion in a highly three-dimensional manner and causing upstream flow redistribution in bothradialandcircumferentialdirections. Theeffectofinflowdistortiononpropulsoroperationcanbeconceptuallydescribedasfollows. The total pressure distortion creates corresponding non-uniformities in axial velocity that leads to non-uniformpropulsoroperatingconditionsaroundtheannulus. Consequently,thepropulsorpro- duces non-uniform pressure rise that modifies the pressure distribution in front of the propulsor. Specifically,thepartofthepropulsorthatoperateswithlowertotalpressure,thuslowerflowcoef- ficients, produces more work, creates larger pressure difference, and causes lower upstream pres- sure along with higher upstream axial velocity. The influence on upstream axial velocity works to reduce the amount of non-uniformity. As illustrated in Figure 5.2, the upstream pressure non- uniformitycausesflowredistributionfromthehightotalpressureregiontothelowertotalpressure region. The flow redistribution entails swirl and radial flow that further changes the local propul- sor operating points. The amplitude and shape of the distortion gets modified by upstream flow redistributionandthroughthepropulsor. Figure5.2: Schematicillustratingupstreamflowredistributioneffects. Figureadoptedfrom[70]. 82 5.2 Rotor67withRadialDistortion To test the capability of the body-force model to predict propulsor performance under distorted flow, we start by applying simple axisymmetric inlet total pressure distortion profiles with two different magnitudes, referred to as low and high distortion levels, to NASA Rotor 67. The inlet radial profiles for total pressure and axial velocity can be seen in the first column (Station 1) of Figure5.3. Forsimplicity,theapplieddistortionvarieslinearlyontheouterhalf-spanoftheblade between freestream total pressure and either 87% or 90% of p t • at the tip for the low or high distortion cases, respectively. The applied flow distortion simulates a scenario where a propulsor seesslowerflowinthetipregionduetoboundarylayersdevelopingoveralongduct. Itisimportant to note that the results presented in this section are generated using an early version of the body- forcemodel. Asaresult,thenon-distortionpredictionmightnotbeconsistentwiththefindingsin Section4.2. However,theprimaryfocusshouldbeonthequalitativetrendsobserved. Spanwise plots of predicted flow properties upstream and downstream of the rotor (Stations 2 and 3) are plotted in Figure 5.3 with and without distortion at the same non-dimensional mass flow of around 0.875, which illustrates how the distortion specified at the domain inlet transfers downstreamacrosstherotor. Thetotalpressuredistortionisseentobeattenuatedthroughtherotor, butthetipregionstillexperiencesalowertotalpressureinthepresenceofdistortionwhiletheroot sees a larger velocity due to the higher loading and the streamlines get deflected radially inward. The total pressure ratio across the rotor at the tip is much higher than at the root implying that the blades do more work at the tip due to slower flow. This is also reflected in the total enthalpy profileatStation3,whichshowsthatthecenter-bladeregionseesalowerworkinput—wherethe velocity is higher—in the presence of distortion compared to the clean flow case. The absolute swirl angle at Station 3 also shows that more flow turning is produced in the tip region where the flowisslower. The predicted overall performance under flow distortion is plotted in Figure 5.4. Both effi- cienciesandtotalpressureratioaregenerallyreducedinthepresenceofdistortioncomparedwith operation under clean flow. The presence of inlet distortion is seen to shift the pressure ratio and 83 0.8 0.9 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Station 1 Clean flow Low distortion High distortion 0.8 0.9 1.0 Station 2 1.3 1.4 1.5 1.6 Station 3 0.0 0.2 0.4 0.6 0.8 1.0 P t /P t 0.0 0.2 0.4 0.6 0.8 1.0 Span, percent (a)Totalpressure p t /p t • 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 0.4 0.6 0.8 1.0 1.2 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 V x /V x 0.0 0.2 0.4 0.6 0.8 1.0 Span, percent (b)AxialvelocityV x /V x • 0.99 1.00 1.01 0.0 0.2 0.4 0.6 0.8 1.0 0.99 1.00 1.01 1.100 1.125 1.150 1.175 1.200 0.0 0.2 0.4 0.6 0.8 1.0 H t /H t 0.0 0.2 0.4 0.6 0.8 1.0 Span, percent (c)Specifictotalenthalpyh t /h t • 1.0 0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.0 0.5 0.0 0.5 1.0 40 50 60 0.0 0.2 0.4 0.6 0.8 1.0 Absolute swirl angle, 0.0 0.2 0.4 0.6 0.8 1.0 Span, percent (d)Absoluteswirlanglea Figure 5.3: NASA Rotor 67 distortion transfer: flow properties at different axial stations in clean flowandwithinletdistortion. 84 efficiency curves towards lower mass flows. These results also indicate that stronger inlet distor- tionresultsinhighefficiencylevelsbeingachievedoveranarrowerrangeofoperatingmassflows: forinstanceanisentropicefficiencyhigherthan92%canbereachedforvaluesof ˙ m/˙ m ref between roughly 0.86 and 0.92 in clean flow but only between 0.86 and 0.89 in heavy distortion. Further- more,themaximumefficiencyexperiencessmalldropsfrom92.9%incleanflowto92.7%inlow distortion and 92.3% in high distortion. Thus, the body-force model gives the desired trends in relationtothedistortioneffectsonwork,efficiency,anddistortiontransfer. Figure 5.4: NASA Rotor 67 performance in clean flow and with inlet distortion at 90% design wheelspeed: effectofdistortionandcomparisonwithexperimentaldatafoundinFidalgo[5]. It is worth noting that the operating point matching process might be different when assessing the influence of inlet distortion from the perspective of an aircraft than it is for isolated propul- sors. Comparing performance between clean and distorted inflow cases at a given mass flow is appropriate from the perspective of the propulsor. However, when considering an aircraft system with and without integrated BLI propulsors, the same aircraft flight condition will likely result in a different propulsor mass flow due to the mass displacement effect of boundary layers. Thus, as discussed in the context of the D8 studies [63, 1] there is more than one way to define equivalent propulsoroperationforaircraftwithandwithoutdistortion,andtheoptimalpropulsoroperationin aconventionalengineinstallationmightbeverydifferentfromtheoptimaloperationwithBLI. 85 5.3 Rotor67withCircumferentialDistortion 5.3.1 ProblemSetup Inthissection,weinvestigatefullwheelRotor67testcasesusingthebody-forcemodel,focusing on circumferential distortion that spans a 120-degree sector of the inlet flow field. Our results are compared with full-annulus blade-resolved unsteady RANS results by Fidalgo [5] which used the same distortion profile and the same fan geometry. Experimental investigations by Gunn [72] and Mistry [74], focusing on fan distortion response with circumferential distortions of 60 and 90-degreesectors,alsoprovidevaluableinsightsintotheproblem. Figure 5.5 shows theC p t contours of the distorted inlet flow field, as well as the circumferen- tial distortion profile at mid-span. The comparison with full-wheel CFD results demonstrates the consistency of the applied distortion. The intensity of inlet distortion is measured by a indicator, DC, defined by Reid [71], the distortion applied in the current test case correspond to a distortion level of DC 120 = 0.83. In the following sections, we present flow contours both upstream and downstream of the rotor to illustrate the effects of distortion. These are taken at the same station locationsasusedinthecleanflowtestcasesdiscussedinSection4.2. ✓ =0 120 Figure5.5: Inlettotalpressurecontour(left)andcircumferentialtotalpressureprofileatmid-span (right). 86 5.3.2 UpstreamFlowRedistribution Whenexposedtocircumferentialinletflowdistortion,therotoroperatesundervaryingconditions alongtheannulus. ThisisillustratedintheflowcoefficientcontourinFigure5.6,whichresultsin non-uniform work input from the rotor around the annulus. Specifically, the sector operating with lowertotalpressureoperateswithalowermassflowrate,therebyproducingalargerpressurerise. Thisnon-uniformoperationinducesvariousupstreaminfluencesthatcanbeobservedbyanalyzing the upstream flow field as shown in Figure 5.6 and Figure 5.7. Axial contours are shown below which are taken upstream of the rotor with the view direction aligned with the rotor axis and the rotationdirectionbeingcounter-clockwise. Figure 5.6: Upstream (station 1) static pressure contour (left) and flow coefficient contour (right) producedbythebody-forcemodel. ⇠ [deg] Figure5.7: Upstream(station1)swirlanglecontour(left)andradialanglecontour(right)produced bythebody-forcemodel. 87 Figure 5.8: Upstream (station 1) circumferential distribution of flow properties at mid-span, as comparedwithURANSresults. As shown in Figure 5.6, a non-uniform pressure field can be observed upstream of the rotor, whereinaregionoflowerpressureformsupstreamofthedistortionsector. Thisoccursastherotor exertsagreaterpullontheslowermassflowinthedistortionsector,thuscreatingalargerpressure difference. The non-uniform upstream pressure distribution induces a flow in the circumferential direction and radial direction, moving from the high-pressure region to the low-pressure region. ThisphenomenoncanbeobservedintheupstreamswirlangledistributiondepictedinFigure5.7, which shows both a counter-swirl and a co-swirl region due to the circumferential redistribution. A radial flow migration can also be observed in the radial angle contour shown in Figure 5.7, in whichastrongradialflowacrosstheannulusisindicatednearthehub. Thesameupstreamredistributioneffectswerereportedinthehigh-fidelityURANSresults[5], with similar magnitudes of swirl and radial angles. Figure 5.8 provides a qualitative comparison 88 of the mid-span circumferential distribution of swirl angle, normalized mass flux, and pressure at three axial locations upstream of the rotor, against URANS results. The comparison shows strong agreement, validating the body-force model’s ability to reproduce upstream effects akin to those predicted by high-fidelity URANS. An additional observation is that distortion-induced non-uniformity in pressure and swirl angle intensifies as the flow approaches the rotor, while the imposed distortion in mass flow is attenuated near the rotor. This suggests that the variation in operationalconditionscausedbytheimposedmassflowdistortionisattenuatedbyupstreamredis- tributioneffects,whichisanimportantmechanismthatenablesenginefantooperatewithstability under distortion. A similar effect is also observed in an experiment [72] with fan operating under 60degreeinletdistortion. 5.3.3 DownstreamFlowFieldandDistortionTransfer Figure 5.9 shows the total pressure and total temperature contours downstream of the rotor, as generated by the body-force model, compared with time-averaged flow contours from URANS resultsatthesameaxiallocation. The total temperature contour reveals that the maximum work input occurs in a region where upstream counter-swirl flow combines with lower axial velocity in the distortion region, leading to the highest local incidence angle. In contrast, the lowest work input is observed in a region where upstream co-swirl results in a locally lower incidence angle. The total pressure contour is more complex due to its origination from a non-uniform profile upstream of the rotor and the subsequent non-uniform work input, which generates a larger pressure rise within the distorted region. The resulted total pressure contour exhibits a maxima just outside the distortion sector on the side where counter-swirl occurs, leading to a high incidence angle. Lower total pressure is observedintheco-swirlregion. Acomparisonbetweenthebody-force-producedflowfieldandthetime-averagedURANSflow field, shown in Figure 5.9, clearly illustrates that the body-force model is capable of qualitatively 89 reproducing the time-averaged distortion transfer effects as predicted by bladed URANS simula- tions. The comparisonalsoshowsthatthemaximaandminimaoftotalpressureandtotaltemper- ature from the body-force predictions are higher and lower respectively, when compared with the URANScontours. Figure 5.9: Downstream (station 2) flow contours of total pressure (left column) and total tem- perature (right column) produced by the body-force model (top row), as compared with URANS results(bottomrow). Figure 5.10 presents a quantitative comparison of the mid-span circumferential distribution of flow swirl angle, total pressure, and total temperature downstream of the rotor, against URANS results and experimental measurements from [5]. The results generated by the body-force model alignwellwiththeURANSresults,validatingthemodel’sabilitytoaccuratelyreproducedistortion transfer. 90 We observed that the body-force model tends to overestimate the peaks of total temperature and total pressure, as previously noted in the flow contours. Higher flow turning is produced in thedistortedregion,asindicatedbytheabsoluteswirlangleplots. Thebody-forcemodel’sresults exceedthosepredictedbyURANSinthisregard,butthebody-forcemodelshowsbetteragreement withexperimentaldata. Figure5.10: Downstream(station2)circumferentialdistributionofflowpropertiesatmid-span,as comparedwithURANSresults. ComparingentropycontoursupstreamanddownstreamoftherotorinFigure5.11,weobserve higher entropy production in areas where the work input is highest, as indicated by the total tem- peraturecontours. Thisresultisexpected, asareasofhigherworkinputarealsoregionsofhigher energy conversion and hence, greater entropy production. The maxima of entropy observed near thehubisassociatedwithradialflowmigrationacrosstherotor. 91 It’s worth noting that the body-force model test case required significantly less computational resources than the URANS test case. Specifically, the former utilized 5 million cells and required 5 hours on 32 processes, while the latter involved 42 million cells and took between 1-2 months of runtime on 128 processes. The previously shown comparisons demonstrate that that the body- forcemodelcandeliversimilarpredictionsofdistortiontransfereffectswithasignificantlyreduced computationalcost,whichmakesitafeasibleoptionforuseinaircraftdesignsimulations. Figure5.11: Contoursofentropygenerationupstream(left)anddownstream(right)oftherotor. Figure5.12: Rotor67overallperformanceunderdistortion. 92 5.3.4 DistortionEffectonOverallPerformance Thenon-axisymmetricupstreamflowconditioncausestherotortooperateundervariousoff-design conditions around the annulus, leading to a decrease in both pressure rise and efficiency. A com- parison presented in Figure 5.12, involving cases with and without a 120-degree circumferential distortion,revealsareductioninpressureriseandfanisentropicefficiencythroughouttheoperating massflowrange. Our body-force model predicted a decrease of approximately 2% in the pressure ratio and a roughly 5% drop in isentropic efficiency due to the distortion. In contrast, the URANS results from Fidalgo [5] showed a less substantial impact of distortion on overall performance, with al- mostnegligiblechangeinthepressureratioandonlyapproximatelya1.5%decreaseinefficiency. The discrepancy in these predictions could be partially attributed to differences in the averaging processesemployed. Ourperformancemetricswerecomputedbasedonmass-averagedquantities across the annulus at upstream and downstream measurement stations. Meanwhile, Fidalgo [5] utilized a more complex sector-based method that tracks sector streamtubes from upstream to downstreamoftherotorandevaluatesperformancebasedonthesesectors. However, upon comparing our results with an experiment conducted by Gunn [72], which featuredasimilardistortionlevel(DC 60 =0.83),wefindtheimpactonperformancetobeconsistent. Gunn’s experiment reported a 4.6% efficiency drop and a 1.6% decrease in pressure rise, which alignswellwithourpredictionsusingthebody-forcemodel. 5.3.5 Summary Inthissection,weutilizedthebody-forcemodeltoexamineRotor67testcaseswithcircumferen- tial distortion in the inlet flow field. The results were compared with a high-fidelity full-annulus unsteady RANS simulation that fully resolved the rotor blades. The body-force model success- fully reproduced upstream flow redistribution effects, including non-uniform pressure field, swirl angle, and radial angle distribution, which qualitatively agree with the URANS results. The in- let distortion and the upstream redistribution caused the propulsor to operate at various off-design 93 conditionsaroundtheannulus,leadingtonon-uniformworkinputfromtherotor. Therotordown- streamflowfieldpredictedbythebody-forcemodelshowedstrongcorrelationwiththatproduced by URANS. A 2% drop in pressure rise and a 5% decrease in efficiency were predicted by the body-force model due to the distortion, in line with findings from relevant experimental works. Moreimportantly,thebody-forcemodelachievedtheseresultsusingsignificantlyfewercomputa- tionalresourcesthantheURANSmodel,demonstratingitspotentialapplicabilityinaircraftdesign simulations. 94 5.4 TF8000FanStagewithBLIInflowDistortion 5.4.1 OverviewofBoundaryLayerIngestion(BLI) Boundary layer ingestion (BLI) is a novel concept in aircraft design that aims to improve propul- sion efficiency by strategically positioning the propulsion system to ingest the boundary layer generatedbytheairframe. Thisapproachhasgainedincreasingattentioninrecentyearsduetoits potential to enhance overall aircraft performance, reduce fuel consumption, and contribute to the decarbonizationoftheaviationindustry. In conventional aircraft designs, the propulsion system is typically installed relatively away fromtheairframetominimizetheinteractionbetweenthetwo. However,byintegratingthepropul- sionsystemwithinorneartheairframe,BLIallowstheaircrafttotakeadvantageofthesynergistic effects between the airframe and the propulsion system. Ingesting the boundary layer reduces the wake dissipation behind the aircraft, which in turn reduces the overall drag experienced by the airframe. Simultaneously,thepropulsionsystemexperiencesareductioninjetdissipation,further enhancingpropulsiveefficiency[1]. Novel aircraft designs incorporating BLI concepts have emerged in recent years, spurred by the growing demand for more efficient and environmentally friendly aircraft. These designs of- ten feature tighter integration between the airframe and the propulsion system, resulting in more complex interactions that can be exploited for improved performance. Examples of such designs includetheD8DoubleBubblebyMIT[2]andNASA’sSTARC-ABL[3]advancedtransportcon- cept, which have demonstrated significant potential of improvements in overall efficiency. As the aviation industry continues to pursue more sustainable and efficient solutions, the application of BLIinnovelaircraftdesignsisexpectedtoplayacrucialroleinshapingthefutureofairtravel. 5.4.2 ProblemSetup In this section, we employ the body-force model to examine full wheel TF8000 fan stage test cases, with an emphasis on boundary layer ingestion type distortion. This distortion type entails 95 centerline 1 2 3 4 5 6 Rotor Stator r x Figure5.13: MeridionalviewofcomputationaldomainandgridsforTF8000fanstage. total pressure distortion in both radial and circumferential directions. While earlier works primar- ily focused on idealized distortion profiles, the fan distortion response with boundary layer types of total pressure inlet distortion has been thoroughly investigated in a recent work by Gunn and Hall [73], which provides valuable insights into our current problem. We compare our results withwindtunneltestresultsontheTF8000fanoperatinginbothcleanflowanddistortedflow,as reportedbySiu[62]. A meridional view of the computational grid is shown in Figure 5.13 with axial location la- beled. Notethatoversetgridisutilizedtoincludethespinner. Thefinalgridisafullrevolutionof themeridionalgridwith180cellsinthecircumferentialdirection. Six axial locations are considered for the TF8000 stage. The inlet conditions are shown at Station1(x= 0.1475m). Therotorfanfaceislocatedatx=0(fanrootleadingedge),designated asStation2. Station3islocatedroughlyhalfwayalongthefanchordatx=0.0161m,andStation 4 at x=0.0341m is between the rotor and the stator. Midway along the stator at x=0.0689m is Station 5. The most downstream Station 6 is placed roughly one rotor chord after the stator at x=0.1624m. One difference between the simulations performed in this work and the experimental data of Siu [62] is that in the latter, the propulsor includes a nozzle and therefore the static pressure will not be uniform there, while at that same x location the domain is straight and the pressure will be forced to be uniform. This is however not expected to affect the comparison of total pressure and efficiency. 96 The TF8000 propulsor performance under flow distortion is evaluated using a stratified inlet total pressure profile that is representative of that ingested by an integrated BLI propulsor. Fig- ure 5.14 shows the two inlet distortion profiles applied at the domain inlet (Station 1), plotted togetherwiththetwoinletdistortionprofilestestedintheexperimentsandsampledalongthecen- terline(z-axis). ThecontourplotoftheadvanceratiocontourisalsoshowninFigure5.14forthe lowdistortioncase. Forsimplicity,theexperimentaldistortionprofilesareapproximatedaspiece- wise linear functions that reduce the total pressure from its freestream value to aC p t of 0.8 and 1.0 for the low and high distortion levels, respectively. The distortion extends from the 60% of thediametertothebottomofthepropulsor. Itcanbeseenthatthelowerdistortionleveliscloseto thetwodistortionlevelsappliedintheexperiments[62],whilethehighdistortionlevelrepresents astrongerdistortion. Figure5.14: Verticaldistortionprofile(left)andinletprofileofadvanceratio(right). 5.4.3 PropulsorMeridionalFlowField Theresultsataflowcoefficientof0.37andwheelspeedof13500RPMwithhighdistortionlevel are chosen to illustrate the distortion transfer effects through the fan stage predicted by the body- forcemodel. The contour plots of various flow properties on the x-z meridional plane of Figure 5.15 and 5.15 show the difference between the top half of the propulsor, where less distortion is ingested, 97 (a) (b) (c) (d) (e) (f) Figure 5.15: Meridional contours of flow properties through TF8000: (a) Normalized total pres- sure, (b) Normalized pressure, (c) Normalized axial velocity, (d) Absolute swirl angle, (e) Nor- malized total temperature, and (f) Entropy generation, at nominal operating point obtained using theproposedbody-forcemodel. Thinblacklinesshowwheretheleadingandtrailingedgesofthe bladeswouldbelocated,andenclosethevolumewherebody-forcesareapplied. 98 and the bottom half where distortion is heavier, as well as how distortion evolves along the axial direction. AsobservedintheNASARotor67tests,theeffectsofdistortiononthepropulsorareevident. Notably, increased turning and work input at the bottom of the propulsor, where distortion causes slower flow speeds, can be seen in the total temperature and absolute swirl angle contours. Fol- lowing the rotor, total pressure is lower, but the total pressure increment is higher on the bottom half. Under distorted conditions, the highest work input occurs near the blade tip, where the ro- tor operates with a locally lower flow coefficient. By contrast, in undistorted sectors, the highest work input is found near the mid-span. It’s also noteworthy that the radial velocity gradient due to distortion experiences some attenuation upstream of the fan, indicated by an increasing axial velocity and decreasing static pressure at the fan face on the lower half of the propulsor. The un- balanced upstream pressure field due to distortion causes a radial flow migration upstream of the rotor. This, in turn, moves the stagnation point upward towards the clean flow region. Similar movements of the stagnation point on the spinner are also observed in the study by Fidalgo [5]. Significantly higher absolute swirl angles are generated at the bottom of the propulsor, leading to higher incidence angles on stator blades and greater stator loading, which results in larger axial flow acceleration through the stator. The propagation of entropy stratification downstream is evi- dent, with some entropy generation visible across the fan stage due to the model’s introduction of losses. Particularly, significant entropy generation is seen at the bottom of the propulsor near the tipregionwhereworkinputishigh. 5.4.4 PropulsorUpstreamFlowField Axialcontourplotsony-zplanesatthefan’supstreamstation2aredisplayedinFigure5.16. This reveals a similar effect of upstream flow redistribution to what was observed in the Rotor 67 test case. Comparison of the fan upstream flow coefficient with the inlet flow coefficient profile shown in Figure 5.14 indicates an attenuation of axial velocity distortion upstream of the fan. Here, the 99 Figure 5.16: Axial contours of flow properties: total pressure coefficient, pressure coefficient, and absolute swirl angle and radial angle, at upstream of the rotor at station 2, viewing from downstream(-xdirection). lowestflowcoefficientinthedistortedregionincreasesfrom0toaround0.2,andadecreaseofthe flowcoefficientisalsonotedatthetopofthefan. The pressure field displays higher pressure at the top and lower pressure at the bottom. This is because the fan pulls harder on the bottom sector due to the slower mass flow. Near the top of the spinner, higher pressure and lower axial flow are observed, corresponding to the upward- shiftedstagnationpoint,asshowninthemeridionalpressureplot. Thesepatternsofupstreamflow distributionwerealsoreportedintheexperimentalresultspresentedbyGunnandHall[73]. The pressure difference induces a flow migration, both radially downward (as illustrated in theradialcontour)andcircumferentiallytowardthebottom(asevidentintheabsoluteswirlangle plot). This results in a counter-swirl region on the left and a co-swirl region on the right, with a higher swirl angle magnitude near the spinner, where the pressure difference is more pronounced duetoflowimpingingonthespinner. 100 Asmootherradialangledistributionisseeninthiscase,comparedtothatobservedintheRotor 67 test case with circumferential distortion. This is likely due to the smoother distortion profile in thecurrentcase. Insummary,thedistortionprofileinitiallyimposedatthedomaininletundergoesmodifications as it approaches the fan face. The axial velocity distortion is notably attenuated, and additional distortions in flow angles are introduced in both the circumferential and radial directions. This further increases the complexity of the non-uniformity of local operating conditions at the fan face. 5.4.5 FanResponseandDistortionTransfer The contour plots on axial y-z planes at three axial stations, just upstream of the rotor (station 2), in the middle of the rotor (station 3), downstream rotor the rotor (station 4), and downstream of thestator(station6)areshowninFigures5.17and5.17inwhichwearelookingfromdownstream with the freestream directed out of the page. They reveal distortion transfer effects along the circumferenceofthepropulsor. The fan responses to the distortion are characterized by two main flow features: the co- swirl/counter-swirl effects upstream of the rotor and the slower flow region at the bottom of the propulsor. The fan face absolute swirl angle a contour in Figure 5.17(b) at Station 2 reveals the presence of regions of counter-swirl on the left of the spinner hub and co-swirl on the right at the fan face, both with maximum magnitudes of about 20 , which result from the interaction of the incoming distortion and the rotor’s upstream influence, The counter-swirl increases the incidence anglelocallyandresultsinahigherworkinput,withtheoppositeeffectintheco-swirlregion. The low mass flow at the bottom of the fan also experience higher work input due to high incidence angles. Theseeffectcanbeseenintheevolutionofvariousflowproperties. The stratified velocity profile develops in such a way that a higher axial velocity increment is producedinthecounter-swirlregionduetoincreasedloading. Thefinalaxialvelocitydistribution at station 6 displays higher velocity at the top, originating from the higher velocity region, and 101 Station2 Station3 Station4 Station6 Figure 5.17: Axial contours of flow properties: normalized axial velocity (first row), pressure coefficient (second row), absolute swirl angle (third row), and radial angle (fourth row) at three differentaxiallocations(2,3,4and6),viewingfromdownstream(-xdirection). 102 Station2 Station3 Station4 Station6 Figure 5.18: Axial contours of flow properties: total pressure coefficient (first row), normalized totaltemperature(secondrow),andentropyrise(thirdrow),atthreedifferentaxiallocations(2,4 and6),viewingfromdownstream(-xdirection). 103 higher velocity on the left, a result of the higher fan loading. The lower part of the fan also experiencesaconsiderableincreaseinaxialvelocity,butitsmagnituderemainslowerthantherest ofthefan,asitoriginatesfromaregionoflowervelocity. Regarding pressure, a significant rise is observed from the bottom of the fan, with a greater pressure increase coming from the counter-swirl region than the co-swirl region. The circumfer- entialnon-uniformityinpressureisalmostentirelyattenuatedasitpassesthroughthestator. Similartoourpreviousobservationinthemeridionalcontours,significantlargerswirlangleis createdatthebottomofthefanthanotherparts. Atrotordownstream,wesawsimilarlevelofflow swirlattheco-swirlregionandthebottomlowmassflowregion,withtheformerregionoriginated fromahigh-swirlangleregion. Theupstreamradialmigrationremainnoticeablethroughtherotor withitscircumferentialreducedandeventuallyfullyattenuatedacrossthestator. Inthetotalpressurecontours,highertotalpressuregenerationisobservedclosertothecounter- swirl region, as anticipated. The non-uniformity evolves similarly to that shown in the axial ve- locity, with higher total pressure observed at the top-left side of the fan. This is because there is greaterworkinputinthecounter-swirlregion,andthetopofthefanoriginatesfromaregionwith highertotalpressure. Thetotaltemperaturecontoursprovideaclearerdepictionoftotalworkinput, consideringthe upstream total temperature is uniform. The highest work input is evident in the low mass flow region at the bottom, while the second highest work input comes from the counter-swirl region. Theco-swirlregioncontributestheleastamountofworkinput. The entropy generation is indeed higher on the lower half where the work input is higher, indicatingalocallylowerefficiency. Thisisanexpectedoutcomeashighworkinputgenerallyac- companieshigherlosses. Thehigherentropyregionisthenshiftedclockwiseduetotherotational movementoftheflow. Figure5.19presentsthestageexit(Station3)totalpressurecontourfromthebody-forcemodel alongsideflowsurveyresultsfromexperiments[62]. Thecomparisondemonstratesthatthebody- force model qualitatively replicated the pitch-averaged distortion transfer effects observed in the 104 Figure 5.19: Comparison ofC p t contour at stator exit (station 3) between body-force model (left) andexperiment(right). experiments, with a similar magnitude. The contour indicates that a higher total pressure was generated from the top left portion of the fan, which is attributable to the combined influences of upstreamswirlandmassflowdistortion. 5.4.6 DistortionEffectonOverallPerformance Thepredictedoverallperformanceofthefanstageundertheinfluenceofflowdistortionsisplotted in Figure 5.20. The inlet distortion has a very small (less than 1%) effect on loading, as indicated bythetotalpressureincreaseshownbyboththebody-forcepredictionandtheexperimentaldata. Thereisaqualitativediscrepancybetweensimulationsandexperiments: thebody-forcemodel predictsthathigherdistortionlevelgivesslightlyhigher(⇠ 1%)totalpressureriseatthesameflow coefficient, while the experiment shows that higher distortion reduces the total pressure rise (still by less than 1%). Given that in the presence of distortion part of the rotor operates at lower flow coefficients, the discrepancy observed is likely due to the over-prediction of total pressure rise at lowerflowcoefficientsareas. Withregardstotheefficiencyprediction,thebody-forcemodelcorrectlypredictedthe1%–2% reduction of efficiency near design operating condition (f ⇠ 0.37) observed in the experimental data at the lower level of inlet distortion. It is also noted that the degradation in fan efficiency is muchsmallerathigherflowcoefficientsandismorepronouncedatlowerflowcoefficients. 105 Figure5.20: TF8000fanstageperformanceincleanflowandwithinletdistortion. 106 5.4.7 Summary The body-force model predicts the physically correct distortion transfer behavior and local varia- tionsofloadingandefficiency. Higherworkinputandflowturningisobservedforthelowertotal pressure and lower velocity (flow coefficient) part of the flow. Upstream co-swirl/counter-swirl effects are observed on the TF8000 due to the interaction between the spinner and the incoming stratified distortion. The counter-swirl/co-swirl creates locally higher/lower incidence angles and thus higher/lower work input and flow turning. The parts of the propulsor that operate at lower flowcoefficientscreatealowerefficiencyregionthatpropagatesintheclockwiseduetotherotor- induced turning. The performance reduction is more pronounced for the NASA Rotor 67 which operatesatmuchhigherspeedsthanfortheTF8000fanstage,andwhichisthusmoresusceptible to local operating point variations caused by the distortion. The 1%-2% reduction of overall fan efficiency due to BLI type of distortion as indicated in the experiments [62, 73] was captured by thebody-forcemodel. 107 Chapter6 ExtendedModelforFree-TipPropulsors 6.1 Introduction Theapplicationofbody-forcebasedmethodsonfree-tiprotatingrotors,suchaspropellersorwind turbines, remains to be an unexplored field. These rotors, characterized by their open tips, exhibit unique flow characteristics and challenges that are distinct from those found in closed turboma- chinery components. Free-tip rotors experience complex tip vortex interactions, induced effects, and wake dynamics that influence their performance and efficiency. The accurate modeling of thesephenomenaiscrucialforthedesignandoptimizationofsuchrotors. Whilebody-forcebased methods have shown promise in modeling closed turbomachinery components, their applicability andadaptabilitytofree-tiprotorstillrequirefurtherinvestigationanddevelopment. BladeElementMomentum(BEM)methodhaslongbeenthemainstreammethodforanalyzing and designing propellers and wind turbines due to its unique combination of simplicity, accuracy, and ease of implementation. BFM can offer several advantages when compared to BEM. First, BFMisinherentlydesignedtohandlenon-uniforminflowconditions,asitaccountsforlocalflow variations and interactions with the rotor. This is in contrast to the BEM method, which assumes a uniforminflow across therotor plane andmay require additionalcorrections or modificationsto handle distorted flow cases accurately. Moreover, BFM provides flexibility in geometry represen- tation, necessitating only minimal knowledge of the blade geometry. This makes it particularly 108 well-suited for design studies and optimization tasks, where the propulsor geometry may be un- knownorsubjecttochange. Thisadaptabilityenablestheexplorationofabroaderrangeofdesign possibilities. Inthischapter,weintroducemodificationstothebody-forceformulationinanefforttoextend its performance prediction capabilities for free-tip propulsors. We validate the model using a va- riety of small-scale propeller test cases with diverse specifications, and compare the predictions to wind tunnel test results. Furthermore, we analyze the propeller performance predicted by the body-force model at various angles of attack relative to the free-stream. In the next chapter, we exploreablown-liftwingtestcase,wheresignificantwing-propellerinteractionsarepresent. 6.2 ModificationstoBody-ForceFormulation 6.2.1 WakeInducedCorrections Propellers, similar to finite wings, experience induced effects due to self-induced velocities in the tangential and axial directions, which effectively changes the local flow angles and affects the bladeloading. Inthecontextofthebody-forcemodel,theaxialinductioncanbeaccountedforvia streamtube contraction created by the propulsor model, from far upstream to far downstream of the propulsor. However, the tangential induced velocity can not be readily captured due to BFM’s infinite number of blade assumption. A flow angle correction method is required to avoid load over-prediction. In what follows, we derive a local flow angle correction formulation to account forthetangentialinductioneffects. InsteadofafullBiot-Savartintegrationovertheentirespan,we estimate a pitch-averaged tangential velocity v t based on blade bound circulation via Helmholtzs Theorematacertainradius v t = BG 4p r (6.1) where B is number of blades. The extra factor of 1/2 in Equation (6.1) is due to the circular integration path includes only semi-infinite trailing vortices [9] as illustrated in Figure 6.1. Next, 109 Figure 6.1: Illustration of circulation integration path for obtaining relation in Equation (6.1). Figureadaptedfrom[9]. inordertodevelopalocalformulation,theboundcirculationG isreplacedbyalocalestimationof circulationbasedona loc similartothederivationofEquation(2.25),resulting G loc =p K c a loc Wc, (6.2) andalocalestimationofv t ,writtenas v t = Bp K c a loc Wc 4p r (6.3) wherea loc as defined in Equation (2.22) is a local angle of attack. A correction of flow angle due tov t canthusbederivedas D d i ⇡ v t W = BK c a loc 4(r/c) , (6.4) Equation (6.4) provides a local estimation of the induced angle and accounts for the tangential inductionwhenappliedtocorrectthelocaldeviationangle. The self-induced effects on inviscid loading are introduced via reducing the local deviation angle byD d i in Equation (6.4), similar to the effect of an induced angle for a wing. This angle reducesthelocalflowangle,causingpartoftheusefulloadtobecomeasourceofdragorloss. By 110 projecting inviscid loading force based on the induced angle correctionD d i , an additional viscous body-forcecomponentcanbederivedtoaccountforinducedloss: f v,i = K c (2pd )sin(D d i ) 1 2 W 2 (2p r/B)b|n q | (6.5) 6.2.2 TipLossCorrection Due to pressure equalization at the tip of a propeller, differences arise between an ideal rotor with an infinite number of blades and one with a finite number of blades. Tip loss correction factors are commonly applied in Blade Element Momentum (BEM) methods, which incorporate the assumption of an infinite number of blades to account for such difference. Similarly, in the Body-ForceModel(BFM),atip-losscorrectionisneededduetothesameinfinitenumberofblade assumption. Prandtl derived a correction factor to adapt the optimal circulation theory to account for a finite number of blades. A modified form of Prandtl’s tip loss factor F is adopted to adjust the deviation angle near the tips and account for the tip losses. As introduced in [9], a modified Prandtl’stiplossfactorisdefinedas, F⌘ 2 p cos 1 (e f ), (6.6) where f = B 2 ✓ 1 r r tip ◆ r 1+ 1 l 2 (6.7) andl isalocaladvanceratiodefinedasl ⌘ V ax /(W r tip ). NotethatF tendstounityasBapproches to infinity and tends to zero as r approches to r tip . Tip loss factor, F, is directly applied to the inviscidloadingforce,resultinginaflowanglecorrection D d tip =(F 1)d . (6.8) 111 6.2.3 LoadLimit For any airfoil, there exists a range of incidence angles within which the airfoil can maintain lift production without experiencing flow separation. In BFM, however, there is no mechanism to predict flow separation from propulsor blades. Consequently, the inviscid loading forces continue to scale with local deviation angles, even at operating points where flow separation might occur, resulting in over-predictions of load at lower advance ratios. A generic mechanism is required to limit blade loading or, equivalently, the flow turning created by the body-force. One option is to directlysetalimitontheincidenceangle,butthelocalizedmodelingapproachusedinBFMmakes it difficult to determine the incidence angle when evaluating source terms away from the leading edge. Alternatively, we can leverage the observation that axial velocity does not experience rapid changes across the propulsor, allowing us to set a limit based on advance ratio evaluated based on local velocities. To do this, we first define the lower bound on the tip advance ratio, which acts as aproxyfortheincidenceangle, l min = V ax min W r tip , (6.9) beyond which the blade loading will be reduced. The value of l min can be determined based on known design information or selected based on modeling needs. Then, a local reduction of the deviationanglerequiredtolimittheloadingcanbedevisedas D d lim =tan 1 ✓ W r V ax min ◆ tan 1 ✓ W r V ax ◆ (6.10) whichisonlyappliedwhenthefollowingistrue V ax <V ax min . (6.11) This correction method effectively changes the blade shape by reducing the deviation angle thus reducingtheflowturningthatwouldotherwisebeproduced. 112 6.2.4 SummaryofChangestoBody-ForceFormulation The proposed formulation to extend body-force model for freetip propellers is applied as follows. The additional flow angle corrections proposed, as seen in Equations (6.4), (6.2.3), and (6.8), are directlyappliedtothelocaldeviationangle d =d +D d i +D d lim +D d tip , (6.12) and the induced losses specified in Equation (6.5) is applied by introducing a new component to theviscousbody-forcevia f v = f v + f v,i . (6.13) It is important to note that the application of body-force models to free-tip propulsors has not been extensively explored in the literature, and the theoretical foundation remains an open area of research. We propose this preliminary model in an effort to make body-force prediction more accurate without complex calibration steps. However, without properly accounting for nuanced flowphenomenasuchasflowseparationandlaminarseparationbubbles(LSB),discrepanciesfrom actual performance, especially in off-design conditions, are expected. Manual calibration can be performedbyadjustingthestrengthofboundcirculationcorrectionspecifiedinEquation(2.31). 113 6.3 PropellersinUniformFlow A selection of APC propellers with different diameters, pitches, and numbers of blades were adopted as test cases, and the results were compared to experimental data from both the UIUC propeller database [59] and experimental measurements conducted at the USC Dryden wind tun- nel. It is worth noting that the propeller blade geometry data from APC website [58] does not cover the entire span down to the hub; it only provides information on the thrust-producing part of the span, starting from the tip to a specific radial location. To address this issue, we visually approximated the missing portion near the hub to ensure numerical stability. It is important to emphasizethatthisapproximationshouldnotaffecttheoverallthrustprediction. 6.3.1 ComputationalSetup Forpropellertestcases, far-fieldboundariesareappliedasthedomainboundarycondition. These boundaries are located at least 60 times the propeller radius away from the model region, as illus- trated by the schematics in Figure 6.2. Two types of domains are used: full wheel domain, which includes the full circumferential extent of the propulsor, and periodic domain, which spans 30 de- grees in the circumferential direction with periodic boundary conditions applied on the domain boundaries in the azimuthal direction. Periodic domains are used for test cases with uniform up- stream flow, where the flow field is expected to be axisymmetric. Full wheel domains are applied fortestcaseswithnon-uniformupstreamflow,suchasthoseshowninSection6.4. One major difference between the two domains is the treatment of grid singularity at zero radius. For the full wheel domain, the center core domain is meshed with a multi-block strategy. For the periodic domain, the center core cannot be meshed due to singularity; instead, a hub wall surface is included to enable axisymmetric grids, and a slip wall condition is applied along the hubwall. Theoperatingpointofthepropellerissetbycontrollingthefreestreamvelocityandthe rotationalRPMwithinthemodelsetting. 114 Far-field Model 60R Figure 6.2: Schematic of the computational domain (left) and meridional grid view zoomed in on themodelregion(right)fortheAPC8b7Spropellertestcase. 6.3.2 Metrics A set of performance metrics commonly used for propellers, as defined by Glauert [7], is adopted forthistestcase. Thepropelleroperatingpointismeasuredbytheadvanceratio,definedas: J⌘ V nD , (6.14) whereV isfreestreamvelocity,nisthepropellerrotationspeedinrevolutionspersecond,andDis thepropellerdiameter. Thethrustcoefficientisdefinedas: c T ⌘ T r n 2 D 4 , (6.15) WhereT ispropellerthrust. Thepowercoefficientisgivenby: c P ⌘ P r n 3 D 5 , (6.16) Where P is shaft power. we compute thrust and shaft power by integrating the source field region as T = ZZ Z V br ||f|| x dV , (6.17) 115 and P = ZZ Z V br f·UdV = W ZZ Z V br ||f|| q rdV . (6.18) wherefisbody-forcevectorandV isthesourcevolume. Lastly,thepropellerefficiencyisdefined astheratiobetweentheusefulpropulsivepowerandthetotalshaftpowerrequried: h ⌘ c T J c P (6.19) 6.3.3 QualitativeObservations ThemeridionalflowcontoursofafewflowpropertiesthroughthepropellerareshowninFigure6.3 and6.4. Physicallycorrecteffectsontheflowfieldareproducedbythebody-forcemodel: Pressure andtotalpressurerisesareproducedacrossthepropellerregion,exhibitinguniquegradientpatterns linkedtothebladegeometry. Ajet(momentumrise)isproducedasshownbythevelocitycontour, with stream tube contraction predicted. Swirl flow is generated with higher swirl angles near the root. It can be seen that less than 10% of the span near the hub has less effectiveness compare to the rest of the span due to the lower blade tangential velocity and smaller chord. This part of the span produces much smaller pressure rise and work input, and blocks the flow by reducing axial velocityandproducingacounter-swirl. Figure 6.5 shows visulization of Q-criterion isosurfaces produced by the body-force model. The isosurfaces are taken at Q crit =10,100 and 1000. As number of blades tends to infinity, The Figure6.3: APC8x7Spropellervelocitycontourwithstreamlines,atJ =0.5,W =6800RPM. 116 (a)C p (b)C p t (c)Swirlangle( ) (d)Entropygeneration Figure6.4: APC8x7SmeridionalflowcontoursofflowpropertiesatJ =0.5,W =6800RPM. Figure 6.5: Q-criterion isosurfaces from the body-force model, overlaid with normalized velocity contours. Indicatingcylindricalvortexsystematthetipandrootoftherotor. 117 helicalvortexsystemofapropellerreducestoacylindricalvortexsheetattherimofthepropeller diskandarootvortexsheetatthehub. 6.3.4 EffectofFlowAngleCorrections In order to illustrate the effect of the wake corrections introduced in Section 6.2, a comparison hasbeenmadeinFigure6.6betweencaseswithandwithoutwakecorrection,theoriginalturning modelwithnoflowanglecorrections,andexperimentaldatafromUIUC[59]atarotationalspeed of 6800RPM. Without any flow angle corrections (blue line), the model underestimates C T by roughly 30% in comparison with the experimental results. When the bound circulation upwash correctionisapplied(redline),C T increases,exceedingtheexperimentaltrendwithahigherdegree of over-prediction observed at lower advance ratios, where larger blade loading implies a larger amount of circulation. With the wake correction applied (black line), the body-force prediction closelyalignswiththeexperiment. Itshouldbenotedthattheblackcurvealsoincorporatesaload limitcorrectionataroundJ =0.4,effectivelycausingthecurvetoflattenatsmalladvanceratios. 0.0 0.2 0.4 0.6 0.8 1.0 J 0.025 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 T No wake correction Exp Qprop Hall With all correc- tions applied Figure6.6: ThrustcoefficientversusadvanceratioforAPC8X7Spropeller,predictedviathebody- forcemodelatW =6800RPM,withcomparisonsbetweencaseswithandwithoutwakecorrection. 118 6.3.5 OverallPerformance Figure 6.7 presents the overall performance trends, including thrust coefficient, power coefficient, and efficiency versus advance ratio for four APC propellers of varying sizes, pitches, and blade counts. Comparisons are made with experimental data from UIUC [59] for all two-bladed pro- pellers, and with experimental data from Long [76] for the four-bladed APC5x4-4 propeller, for which only theC T -J data is available. Despite the diverse propeller specifications, the body-force modelprovidesaccuratepredictionsfortheC T -J trendintermsofbothslopeandmagnitude. Atsmalleradvanceratioranges,thepropellerbladesencounterhighincidenceanglesandlarger loadings. This leads to significant local axial induction and structural deflection of the propeller, collectively causing the C T -J trend to flatten at smaller advance ratios instead of maintaining a Figure 6.7: Overall performance (thrust coefficient, power coefficient and efficiency) versus ad- vance ratio for various APC propellers, predicted by the body-force model, compared with ex- periemnetaldatafromUIUC[59]. 119 linear trend. This effect is partially replicated by the body-force model due to the application of a loadlimitcorrectionspecifiedinSection6.2. However,largerdiscrepanciesfromtheexperimental dataareobservedatsmalleradvanceratioranges. The C P -J trend is also well represented by the body-force model, with under-prediction ob- served near the higher advance ratio range of operating points where loading is small. It should be noted that without the wake-induced losses term shown in Equation (6.5), a consistent under- predictionofthepowercoefficientacrossalladvanceratiorangeswasobserved. The under-prediction ofC P subsequently results in an over-prediction of efficiency or a post- poned drop in efficiency at higher advance ratios. This is prominently evident in the efficiency plot, which indicates that the maximum efficiency is predicted to occur at a larger advance ratio withaslightlyelevatedmagnitude. 6.3.6 Summary The body-force model formulation is extended to include a wake corrections along with a tip loss correction in order to provide better performance prediction of free-tip propulsors. Using this preliminarybody-forceformulationdevelopedforfree-tippropulsors,themodel’spredictivecapa- bilitieswereassessedacrossaseriesofAPCpropellersofvaryingsize, pitch, andbladenumbers. A qualitative inspection of the flow field verified that the model successfully generates physically correctfloweffects,suchaspressureriseandmomentumriseacrossthepropeller. Regardingoverallperformance,theoriginalturningmodel,withoutanyflowanglecorrection, was found to under-predict the thrust coefficient by 30% for APC 8X7 two bladed propeller. The proposed flow correction improved the prediction, leading to results that closely align with exper- imentalfindings. Agenerallysatisfactorymatchwasachievedforthethrustcoefficientandpower coefficient trends versus the advance ratio for four different APC propellers. An under-prediction of the power coefficient was noted near the higher end of the advance ratio range, resulting in an over-prediction of maximum efficiency and the advance ratio at which the maximum efficiency occurs. 120 6.4 PropellersatNon-ZeroAnglesofAttack 6.4.1 Introduction Recently, there has been a growing interest in eVTOL vehicles due to the increasing demand for air-mobility solutions. Novel eVTOL aircraft concepts are raised to offer efficient, sustainable, andflexibletransportationsolutionswithinandbetweenurbanareas,addressingthegrowingneed for an alternative to traditional ground-based transportation methods. As congestion in urban en- vironments continues to worsen, eVTOL vehicles have the potential to provide a faster, cleaner, and more convenient mode of transportation for both passengers and cargo. eVTOL aircraft are characterized by their ability to take off and land vertically, similar to helicopters, but typically incorporate advanced distributed propulsion systems and innovative airframe designs that enable more efficient forward flight. This unique combination of capabilities allows eVTOL vehicles to operate in a wide range of environments, including dense urban areas with limited infrastructure, whileminimizingnoisepollutionandreducingtheircarbonfootprint. Propellersplayacrucialrole in the performance of eVTOL aircraft, as they are responsible for providing both lift and propul- sion. Depending on the specific eVTOL configuration, the propellers may be required to operate under varying inflow conditions, which can significantly impact their performance and efficiency. For example, during the transition between vertical take-off and forward flight, propellers may experience non-uniform inflow due to changes in the relative orientation of the aircraft and the airflow around it. This can result in complex aerodynamic interactions and operating conditions thatmustbecarefullyconsideredduringthedesignprocess. Inthissection,wesetupanAPC8x7Spropellerbody-forcemodeltestcasewithvariousinflow angles of attack, aiming to investigate the non-axisymmetric responses by the body-force model to skewed inflow. The purpose of this analysis is to evaluate the model’s ability to predict perfor- mancevariationsofthepropellerundernon-uniforminflowconditions. Wecomparethepredicted performance trends from the body-force model against experimental data to assess the model’s 121 capatbility. The angle of attack,a p , is defined as the angle between freestream flow direction and axialdirection(x-axis)withinthex-zplane. 6.4.2 PropellerFlowField We investigated flow angles of attack ranging from 0 degrees to 70 degrees for an isolated pro- peller. Flowcontoursfora p =0,30,50,70aredisplayedinFigure6.8,whichshowcasespressure coefficientandtotalpressurecoefficientflowcontoursasviewedfromtheside(x-zplane)andthe top (x-y plane). The model exhibits clear non-axisymmetric responses with non-zero angles of attack. In Figure 6.8 (a) and (c), variations in load are compared between the top and bottom of the propeller in response to the skewed inflow. Due to the angled inflow, the bottom of the propeller experiences locally smaller axial velocity, leading to a lower local advance ratio and increased loading. This is evident from the flow contours showing a larger total pressure rise and pressure differencecreatedatthebottomofthepropellerwithhigheranglesofattack. Figure 6.8 (b) and (d) demonstrate load variation responses compared between the left and right sides of the propeller. The upward-traveling side (left) exhibits a reduced load while the downward-traveling side (right) shows an increased load. This is due to the vertical component of the freestream velocity compounding with the blade rotation speed, resulting in an increased incidenceangleforthebladetravelingagainsttheverticalfreestreamvelocity(rightside),andvice versa for the blade traveling in the same direction as the vertical velocity. This effect is sometime referredtoP-factorinthecontextofpropellerdrivenaircraftperformance. Figure6.9showcasestotalpressurecoefficientcontoursdownstreamofthepropelleratvarious anlge of attacks, viewed from the axial direction. The propeller rotates clockwise in this view. It’s evident that higher inflow angles of attack yield more significant non-uniform load responses. Moreover, the P-factor effect is more pronounced at lower angles of attack, while the load asym- metry caused by locally lower axial velocity gains prominence at higher angles of attack. This is 122 illustratedbytheshiftinglocationofthemaximumC p t ,whichtransitionsfromtherightsideofthe propellertowardsthebottomastheangleofattackescalates. (a)C p onx-zplane (b)C p onx-yplane (c)C p t onx-zplane (d)C p t onx-yplane Figure 6.8: APC8x7S meridional flow contours ofC p (first row) andC p t (second row) at J =0.5, andatangleofattacka p of0,30,50,and70degrees(fromlefttorightineachsubplot),wherea p isdefinedasflowanglewithrespecttox-axiswithinx-zplane. Figure 6.9: APC8x7S axial flow contours of C p t downstream of the propeller at J = 0.5, and at angle of attack a p of 0, 30, 50, and 70 degrees (from left to right). Propeller rotates clockwise in thisview. 123 6.4.3 OverallPerformance Figure 6.10 plots the overall performance of the propeller, including the thrust coefficient, power coefficient, and efficiency, against the advance ratio for various angles of attack. In general, a net increase in thrust is observed across the studied range of advance ratios, resulting from an overall decreaseinvelocityperpendiculartothepropellerdiskastheangleofattackincreases. Greater performance variations are seen at higher advance ratios due to the higher freestream velocity or lower blade rotation speed at these points, causing changes in the angle of attack to havemoresignificantimpactsonlocaloperatingconditions. For zero or small angles of attack, c T and cP decrease with the advance ratio. However, for a p > 60 , cT and c P start to increase with the advance ratio J. This trend is consistant with experimentalresultsreportedbyStokkermans[77],asshowninFigure6.11. The performance variations observed originate from the locally lower advance ratios expe- rienced by propellers at larger angles of attack. This lower advance ratio results in a delayed efficiency drop at larger a p . It’s worth noting that when a propeller encounters non-uniform up- stream axial velocity, the efficiency definition provided in Equation (6.19), which is premised on a uniform advance ratio, is no longer accurate. This discrepancy may result in efficiencies greater than one. In such cases, it becomes inappropriate to use the freestream advance ratio to calculate efficiency. Wehavethereforeintroducedaslightlymodifiedformofpropellerefficiencytoaccount for the reduced velocity perpendicular to the propeller disk due to skewed inflow. This adjusted formulationispresentedinEquation(6.20). h ⌘ c T Jcos(a p ) c P (6.20) 124 Figure 6.10: APC 8x7S propeller performance at 6800 rpm, including thrust coefficients, power coefficients,andpropellerefficiencyatincreasinganglesofattackaspredictedbyproposedbody- forcemodel. 6.4.4 Summary This section evaluates the influence of inflow angle of attack variations on propeller performance as predicted by the body-force model. The flow field of the propeller revealed that the body-force modelgeneratesphysicallyaccurateloadvariationsinresponsetoskewedinflow. Increased load is observed at the bottom of the propeller due to locally lower advance ratios, while higher load is generated on the right side of the propeller owing to the increased blade incidenceanglecausedbytheverticalcomponentofthefreestreamvelocity. Thebody-forcemodelhasdemonstrateditscapacitytoreproducetrendsofc T andc P variations with flow angle of attack, as shown by the experiments. This success is attributed to the model’s abilitytoadapttolocalflowvariations. 125 increasing↵ p Advance ratio J 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Experimental trends Figure6.11: Trendsofc T variationwithincreasingpropellerangleofattacka p asfoundbyexper- iments. Figureadaptedfrom[77]. 126 Chapter7 CapturingPropeller-WingInteractionsinBlown-LiftSystem This chapter presents the application of a body-force propulsor model to represent propellers in RANSsimulationsandaccuratelycaptureaero-propulsivecouplingeffects. Thebody-forcemodel isappliedtothesimulationofacase withstrongpropeller-airframeinteractions, namelyablown- liftwinginwhichdistributedpropellersblowontoaflapped-wingtoproducehighlift. Resultsare compared with an advanced actuator disk model and with experimental measurements. We look at the load distribution, propeller performance, jet properties, and wing performance for a range ofoperatingconditions. Foranisolatedpropeller,bothbody-forceandactuatordiskmodelsmake comparable flow-field predictions. However, the former is able to predict propeller performance (thrust,power,efficiency)andmatchexperimentaldatawithonlythewheelspeedspecified,while thelatterrequiresthethrustasaninputanddoesnotpredictefficiency. Intheblown-liftapplication, thebody-forcemodelisfoundtosuccessfullypredictwing-propellerinteractions,andrespondsto thelocalflowinawaythatbettermimicsactualpropelleroperation. 7.1 OverviewoftheBlown-LiftConcept Propeller-wing interactions can be taken advantage of to augment lift generation from a wing. Specifically, with the propeller thrust stream blowing onto the wing, higher amount of lift can be generated. This effect is sometimes called blown-lift and is receiving more attention recent years due to advancement in aircraft electrification and the increasing usage of distributed electrical 127 propulsion(DEP).Comparedtoconventionalpropulsionsystemwithtwoengines,DEPenablesa muchlargerportionofthewingbeingtotakeadvantageofblown-lift,thusprovidingmuchhigher lift augmentation. NASA X57 utilizes an array of wing mounted propellers that blows a jet sheet ontothewingtogeneratetherequiredliftfortake-offwithanarrowerwing. Similar ideas were adopted in another concept called blown-flap, in which the DEP jet sheet is deflected downward by a multi-element flap, and this results in even higher lift coefficient, as muchas15.0. Theblown-flapconceptenablesultrashorttake-offandlandingcapabilitiesandwas adoptedbyeSTOLconceptfromElectra.aerothatcantakeoffinultrashortdistance. There are two main mechanisms behind the lift augmentation brought by blown-lift wings. The downward deflection of high-momentum jets via the flap results in larger flow curvature and amountsto higherflowcirculationaroundtheairfoil, hencehigherlift. Thehigh-momentumflow throughtheflapslotsenergizetheflapboundarylayer, whicheffectivelydelaysflowseparationat higheranglesofattack,andhelpachievehighermaximumliftcoefficient. Experimentalworkswereconductedtostudytheaerodynamiccharacteristicsandperformance ofblown-liftwings[78,79]. TwodimensionalinviscidtheoriesbySpence[80]utilizesvortexsheet tomodeltheairfoilandtreatedthedeflectedjetstreamasanextensiontotheflapthathelpedcreate extra pressure loading. The theory provides reasonable approximation for blown-lift performance whenpropellerssizesapproachestothesmallerlimit,representingathinjet. Analyzingtheperformanceofblown-liftwing-propellersystemswithCFDrequirespropulsor modelsthatcanaccuratelyrepresentpropulsorperformancesundertheinfluenceofthenear-wing flowfield. In this section, a blown-lift wing test case is setup in both quasi-3D conditions, with two dif- ferentpropulsormodelsapplied,namely,3Dactuatordiskwithoptimumloaddistribution,andthe body-forcemodel. OurinterestsarefocusedoninvestigatingInthissection,ablown-liftwingtest case is established under quasi-3D conditions, employing two different propulsor models: a 3D actuator disk with optimal load distribution, and the body-force model. Our focus is on exploring theinteractionsbetweenthewingandthepropeller,understandinghowthesecomplexinteractions 128 Table7.1: Geometricalparametersoftheblown-liftwing. Chord(atd f =0 ) 0.2286m Span(perpropeller) 0.1524m Motoraxisangle 10 Propellertipdiameter 0.127m Propellerhubdiameter 0.0108m Propellerreferencearea 0.0126m 2 influence the system performance, and assessing the ability of the proposed body-force model to accurately capture these intricate interactions. This is done through comparison with results ob- tainedfromanadvancedactuatordiskmodel. 7.1.1 Blown-LiftWingGeometry Asatestbedforstudyingthepropulsormodels’capabilities,weconsiderablown-liftwingadopted from a sub-scale wind tunnel experiment conducted at MIT [76]. This wing features distributed propellers that produce a jet system blowing over the flapped wing for lift augmentation. In turn, the high-lift wing creates significant upstream influences that affects propeller performance. This interactionisamajorcharacteristicofblown-liftdesigns. Thewing’sgeometricalspecificationsarelistedinTable7.1andits2Dsectionalgeometrywith propeller location shown in Figure 7.1. The original wing tested in the MIT wind tunnel has four propellers distributed evenly along the span and has walls on both ends of the span (end walls), thus represents a quasi-2D configuration with a 2D (infinite span) wing performance but fully 3D jetdevelopment. For the sake of reducing computational cost, we only simulate a quarter of the span blown by a single propeller as shown on the right of Figure 7.1, with periodic boundary condition applied on the span-wise domain boundary faces to represent an infinite span. Only a small number of full-widthfour-propellercaseswereconductedforthepurposeofinvestigatingjetmixingeffects. The propeller camber surface geometry and thickness distribution needed for the body-force modelwereextractedfromtheAPC5x4E-4propellerthattheMITwindtunnelexperimentsused[76]. 129 10 Figure 7.1: Blown-wing geometry used in wind tunnel tests by [76]: 2D airfoil geometry wth flap angle d f =20 (left); and 3D geometry of quarter wing span with motor pylon and propeller (right). This propeller has 4 blades with a diameter of 5 inches (0.127m), and was tested experimentally atamaximumof29000RPMandproducedupto10.84Nofthrust. All meshes are constructed using the commercial meshing software Pointwise. For blown-lift wing test cases, we took advantage of the ADflow’s overset grid capability to mesh around the complex geometry of the flapped wing and motor pylon, and added a dedicated region for the propulsormodel. Meshconvergencestudieswereperformedseparatelyforthewingmeshandthe propulsormodelregion. Thesewerecarriedoutwith2Dsimulationsofthewingairfoilsduetothe 2D nature of the wing meshes. Three mesh levels were created, and the coarsest of the three was chosenforthepresentstudysinceitproducesalessthan2%differenceinforcecoefficientswhen comparedwiththefinestmeshlevel. Formeshinthebody-forcemodelregion,duetotherelianceontheactualbladegeometry,grid finesse may influence the model’s prediction of propeller performance. There needs to be enough cellsinboththeradialandaxialdirectionstoproperlyresolvetheloaddistribution. Gridpointsare clusteredneartheleadingandtrailingedgesofthemodelregionwhereflowundergoesmorerapid changes. The detailed model region mesh convergence is shown in the Appendix for the isolate propulsor study. The thickness of the actuator disk region was chosen such that it is similar to the averagethicknessofthebody-forcemodelregion. 130 1 2 3 4 5 Figure 7.2: Spanwise cut view through the center of the propeller of the overset computational gridsfortheblownflappedwingata =10 (left),andmeridionalcutoftheisolatedpropellergrid withthepropellerhubatthebottomfigureedge(right). A side view of the assembled blown-lift wing computational grid close to the wing surfaces is shown in Figure 7.2. Axial measuring locations are labeled and marked withvertical dashed lines for reference in later sections. After the overset cutting, the mesh has a total of approximately 5.6 million compute cells with a y + of about unity on solid walls. The far-field domain boundaries extendmorethan100chordsawayfromthewingsurfaceinalldirections. Each single-propeller blown-wing solution was computed in around 600 to 800 CPU-hours with 32 processes on supercomputers at the Center for Advanced Research Computing (CARC) oftheUniversityofSouthernCalifornia. Convergencereductionof8ordersintotalresidualwere generally achieved. Except for cases with more inherent unsteadiness, convergence may stall ear- lier, in which case care was taken to make sure that the lift and drag coefficients had converged. The computational cost of the body-force model is slightly higher than that to the actuator disk modelduetoitsrequirementoffinermodelregiongridsandsomeadditionalcomputations,butin practicetheauthorsfoundthecostofsolvingwithbothmodelscomparableduetothemuchhigher costofsolvingtherestoftheproblem. Notethatfortheblown-wingcases,thepropellersourceregiondoesnotgobeyondthehubof thepropellerwherethebladeends. Inotherwords,thesourceregiononlycoversthevolumeswept bytheblades,andnotthecenterhubandmotor-mountingareas. 131 7.1.2 ActuatorDiskModel Forcomparison,weimplementedanadvancedactuatordiskmodel,whichprovidesmorerealistic load distribution and flow swirling than just an uniform pressure jump. The actuator disk model specifiesatotalthrust,T spec ,thatisradiallydistributedfollowingBetzoptimalcirculationdistribu- tionwithPrandtl’stip correction. Its localforcepervolume canbegenerallydescribedasfollows f AD = f AD (r,R,B,T spec ,l ), (7.1) where r is the local radius (distance to propeller axis), R the propeller tip radius, B the number of blades and l the tip advance ratio. Note that the application of this model needs as input the thrust versus advance ratio relation for the specific propeller being represented. For this, we used the c T -l information from the wind tunnel experiments (shown in the results section) in order to settheadvanceratiowhenthrustisspecifiedastheinputparameter. Theforcetermf AD hasaxialandtangentialcomponentsexpressedas f a = T spec V K(r) (7.2) and f q = T spec V K(r)l w (r), (7.3) respectively, where l w is the local wake advance ratio which varies with radial location, K is a scaling factor responsible for distributing the loading, and V is the total volume of the model region. By assumeing Betz optimal circulation distribution and applying Prandtl’s tip correction, K takestheform K(r)= 1 K c F(r)C G Be (r), (7.4) withthemodifiednon-dimensionalformofBetzoptimaldistribution C G Be (r)= x 2 0 +1/l 2 w 1+1/l 2 w +x 2 0 , (7.5) 132 Figure 7.3: Actuator disk radial load distri- bution for different x o constants with B=4, l =0.1,andr hub =0. Figure 7.4: Actuator disk radial load distribu- tion for various advance ratios l with B=4, x 0 =0,andr hub =0. ascalingfactorforapplyingPrandtl’stipcorrectionsetto F(r)= 2.0 p cos(exp( f)) (7.6) where f = 1 2 B ⇣ 1.0 r R ⌘ s 1.0+ 1.0 l 2 w tip , (7.7) andthescalingfactor K c = 1 R r hub Z R r hub F(r)C G Be,o (r)dr (7.8) thatensuresthetotalintegratedthrustisthesameasthespecifiedvalueT spec . The constant x 0 in Equation (7.5) provides a way to adjust the loading near the center of the propeller. As shown in Figure 7.3, lower x 0 gives lower loading near the center and vise versa. In this work, we choose x 0 =0 to better compare with the body-force model which generates zero loadingatthecenter. Figure7.4showstheeffectofadvanceratioonradialloaddistribution. Thecomputationalcostofthebody-forcemodelisslightlyhigherthanthattotheactuatordisk modelduetoitsrequirementoffinermodelregiongridsandsomeadditionalcomputations,butin 133 practicetheauthorsfoundthecostofsolvingthepropulsormodelissmallcomparedtothecostof solvingtherestoftheprobleminapplicationstofullaircraftsimulations. 7.1.3 ComparisonbetweenADandBFMinisolation Inthissection,thetwopropulsormodelsareappliedtorepresentastand-aloneAPC5x4-4propeller withanuniforminflowthatisalignedwithitsaxis,implyingthatallsolutionsareaxisymmetrical. We examine the local flow property for a nominal operating point, as specified in Table 7.2. The overall performance predicted by the body-force model is also shown for a range of operating pointsandcomparedwithwindtunnelmeasurements. Table7.2: Operatingconditionoftheexaminedisolatedpropellercase. l 0.085 V • 13.0m/s W 23000RPM c T 0.035 7.1.4 LocalFlowPropertyDistributions Figure7.5showsmeridionalcontourplotsofvariousflowpropertieswithinthebody-forceregion that is marked by thin black lines. Pressure and total pressure rises are gradually produced across the propeller region, exhibiting unique gradient patterns linked to the blade geometry. A jet is produced as shown by the velocity contour, and swirl flow is generated with higher swirl angles neartheroot. Entropyisproducedbythelossmodelingwithinthebody-forcemodel. Higherloss isobservednearthetipduetothehigherblade-relativevelocityandhigherloading. Itcanbeseenthatapproximately10%ofthespannearthehubhaslesseffectivenesscompare to the rest of span due to the lower blade tangential velocity and smaller chord. This part of the span produces much smaller pressure rise and work input, and blocks the flow by reducing axial velocityandproducingacounter-swirl. However,theimpactonperformanceisnegligibleandthat hubregionwillbeblockedbythemotoronceinstalled. 134 The distribution of various quantities computed by the body-force model are shown in Fig- ure 7.6. They are related to fictitious blade relative flow and are used in the model formulation. ThesecontoursthusshowtheactualdistributionsusedinthemodeldescriptionequationsofChap- ter2. Figure7.7showsthesamesetofflowcontoursproducedbyusingtheactuatordiskmodelatthe samecondition. Overall,similarchangesofflowpropertiesacrossthebladeregionwereobserved when compared to that produced by the body-force model in Figure 7.5. The tip region flow propertiesrevealasmoothertransition(smallergradients)fromthemodelregiontothesurrounding flow than that produced by the body-force model, since the latter’s thickness reduces near the tip Figure 7.5: Flow contours predicted by the body-force model for isolated propulsor at l =0.085, W =23000RPM: pressure coefficient, total pressure coefficient, normalized velocity, swirl angle, andentropy. Figure 7.6: Contour of body-force model quantities for isolated propulsor at l = 0.085, W = 23000RPM: local momentum boundary layer thicknessq , Reynolds numbers Re q and Re x , dissi- pationcoefficientC D ,andflowdeviationangled . 135 Figure7.7: Flowcontourspredictedbytheactuatordiskmodelforisolatedpropulsoratl =0.085, W =23000RPM: pressure coefficient, total pressure coefficient, normalized velocity, swirl angle, andentropy. and thus causes sharper changes. Major differences between the models are found in entropy production, with the actuator disk producing an effect that is deemed non-physical and due to numerical errors. One of the shortcomings of the actuator disk model is that it does not include a mechanismtorepresentlossgeneration. Quantitative comparison of the radial profiles of several flow properties taken at locations up- stream (station 1) and downstream (station 2) of the propulsor model region are shown in Fig- ure 7.8. Despite fundamental differences in the two modeling approaches (one using the blade geometry and the other a prescribed load distribution), the predicted radial profiles agrees very well for C p , C p t , and axial velocity for this isolated propeller case. This is likely because this particular propeller (and its geometry used by the body-force) does produce a close-to-optimal loading (which is well-described by the distribution imposed in the actuator disk model). On the otherhand,largedifferencesareseenintheabsoluteswirlangleawayfromthetip,withtheactua- tordiskmodelcreatingnearlytwiceasmuchflowswirlat30%ofthespan. Thisislikelyduetothe simple estimations of circumferential force as scaling with 1/r, which thus becoming inaccurate atverysmallradii. 136 2 1 0 1 2 Cp 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 r/R 0 2 4 6 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 r/R BFMupstream BFMdownstream ADupstream ADdownstream Cpt 0 10 20 30 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 r/R ↵ 0.0 0.5 1.0 1.5 2.0 Vx/V 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 r/R Figure 7.8: Radial profiles at upstream (station 1) and downstream (station 2) locations for iso- latedpropulsoratl =0.085,W =23000RPM: pressurecoefficient,totalpressurecoefficient,swirl angle,andaxialvelocity. 7.1.5 OverallPerformancePredictedbyBody-ForceModel Figure 7.9 shows the characteristic of the isolated APC5x4-4 propeller over a range of operating points as predicted by the body-force model with comparisons to wind tunnel measurements [76] and Qprop prediction. The body-force model predicts the experimental propeller c T l charac- teristicslopewell,andisslightlyoffinmagnitudewithanapproximately5%over-predictionover most of the operating range when compared with the wind tunnel measurements. The body-force prediction also closely matches results from Qprop 1 , even though the latter could be expected to over-predict the performance of small propellers with low Reynolds numbers like this one. The operating point where thrust reduces to zero (aroundl =0.2) is well predicted by the body-force model. Performancepredictionswiththebody-forcemodelatsmalladvanceratios(l ¡0.03)could be improved by adding the off-design term, which was not applied in this work due to its lower importanceinthisapplication. With the advantage of profile loss modeling, the body-force model predicts the physically correct trend in propeller aerodynamic efficiency. Its predictions show lower efficiency at higher advanceratiountilthethrustreducestozero. Wewouldseetheefficiencydropatverylowadvance ratiosiftheoff-designlossmodelwasincluded. Physicallycorrectvariationsofpowercoefficient wasalsoobserved,withahigherrequiredshaftpoweratloweradvanceratio. 1 aclassicalblade-element/vortexcodedevelopedbyMarkDrelaatMIT 137 Figure 7.9: Overall propeller performance for isolated propeller: thrust coefficient (top), aerody- namicefficiency(bottomleft),andpowercoefficient(bottomright)versusadvanceratiol . Com- parison between body-force model (BFM), and both wind tunnel measurements (Exp) from [76] andQproppredictions. 7.2 Wing-PropellerInteractions Inthissection,weconsidertheblown-wingtestcaseandfocusonunderstandinghowtheinterac- tionsbetweenwingandpropelleraffectseachother’sperformance,andonhowwellthisinteraction iscapturedbyeachofthepropulsormodels. The flow fields at a reference operating condition, defined in Table 7.3, is used to describe the basicfeaturesfoundinblown-liftwingsflowfieldsandunderstandthemodelresponsetopropeller inflow. Thischosenoperatingpointrepresentsahigh-liftconditionrepresentativeofwhatisdesir- ableforultrashort-distancetake-offorlanding. Wethenanalyzetheaero-propulsiveperformance 138 and the interaction between wing and propeller, before looking more closely at the propeller jet developmentandmixing. Table7.3: Referenceoperatingconditionsfortheblown-liftcase. Thebody-forcemodeloperating conditionissetviathewheelspeedW ,andtheactuatordiskmodelconditionissetviatheisolated thrustmagnitude. a 8 d f 30.0 V • 12.7m/s l 0.08 W 23900RPM c T (isolated) 0.036 7.2.1 BasicFeaturesoftheFlowOveraBlown-Wing Figure 7.10 shows velocity contour on a wing cross-section through the propeller center, as pre- dicted by the body-force model and the actuator disk model at the reference condition. The basic Figure7.10: Velocitycontoursofwingsectionalflowatthereferenceconditionwiththebody-force model(top)andtheactuatordiskmodel(bottom),inazoomed-outview(left)andzoomed-inview (right). 139 flow features of a blown-lift wing are evident in this view. A jet produced by the propeller blows over the wing and energizes the flap boundary layers through the slot. This allows the flow over the flap to remain attached even at high angle of attack and with high flap deflection angles, thus generatinghighlift. Thejetisalsodeflecteddownwardsbytheflapped,whichalsoprovideslift. A consequence of the high lift and corresponding high circulation is that the bending of the stream- linesandtheirverylargecurvaturecausesthepropellertooperatewithitsaxismisalignedwiththe inflow,ascanbeseenfromthestreamlinesinthefigure. Some differences between the results from the two models can be observed. The velocity magnitude contours indicated that the body-force model produces stronger blowing at the bottom of the propeller than the actuator disk model and, as result, more flow deflection is produced far downstream of the wing as indicated by the streamline directions (left side of Figure 7.10). This contributestopredictingahigherliftandpitchingmomentfromtheblown-wingaswillbeshown inthewingperformanceresultsofSection7.2.3. The sampled flow properties, or profiles, along vertical lines upstream (station 1) and down- stream(station2)ofthepropulsormodelregionareshowninFigure7.11. Thesearetakenalonga linethroughpropulsorcenterthatisperpendiculartothepropelleraxisasillustratedinFigure7.2. These profiles confirm that overall the body-force model predicts higher static and total pressure rises than the actuator disk model. Furthermore, whereas the actuator disk model produces rel- atively similar changes at the top and bottom of the propeller, the body-force results in higher pressure (total and static) rises and higher jet velocity being generated at the propeller’s bottom than at its top. This results in a nearly 15% higher maximumC p rise and a 42% higher maximum C p t riseatthepropellerbottomthanatitstopinthesimulationswiththebody-forcemodel. On the axial velocity profile, we can see that the upstream location has higher V x at the top of the propeller than at the bottom, irrespective of which model is used. The models however differindownstreamaxialvelocityprediction: thebody-forcemodelrisesthevelocitymoreatthe bottom enough to invert the top-bottom ratio of the upstream station, and thus results in a higher axialvelocityat thebottomthanatthe top. The actuatordiskmodel, on theotherhand, maintains 140 2 0 2 4 1.5 1.0 0.5 0.0 0.5 1.0 1.5 y/R Cp BFMupstream BFMdownstream ADupstream ADdownstream 0.0 2.5 5.0 7.5 10.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 y/R Cp t 0 1 2 1.5 1.0 0.5 0.0 0.5 1.0 1.5 y/R Vx/V1 20 0 20 40 1.5 1.0 0.5 0.0 0.5 1.0 1.5 y/R ↵ Figure 7.11: Vertical profiles for blow-lift case reference condition at upstream (station 1) and downstream (station 2) locations: pressure coefficient, total pressure coefficient, normalized axial velocity,andswirlangle. similar relative magnitudes between top and bottom at the downstream station (i.e. higher at the top). Finally, both models produce similar amounts of overall swirl, but the actuator disk induces higherswirllevelsnearthehub. The above results clearly show that the propulsor models react differently to the misaligned inflow, while they behaved similarly for the isolated propulsor. Since propeller jets are crucial componentsoftheblown-liftsystemwhichimpacttheoverallperformancebyaffectingtheforces experienced by the wing, the observed differences emphasize the importance of choosing the cor- rectpropellermodelsforsuchapplications. Inordertobetterexplaintheobserveddifferences,we willexaminethenon-uniformpropellerflowonaxialplanesinthenextsection. 7.2.2 PropellerInflowNon-UniformityandModelResponse Thehighliftfromtheblown-wingproducesasignificantupstreaminfluencethatchangesthepro- peller’s operating condition. Figure 7.12 shows the flow contours on axial planes for various flow properties overlaid with arrows that indicate flow velocity directions measured upstream (station 1) and downstream (station 2) of the propeller. Axial planes in this paper are perpendicular to the propeller axis and the contours plotted looking downstream. Also note that the rotation direction ofthepropellerisclockwiseinthisview. 141 Propeller inflow non-uniformity is primarily characterized by two mechanisms, both of which are the result of large upstream flow curvatures that causes the propulsor’s inflow to be at large and non-uniform incidence angles. First, the propeller inflow shows lower axial velocity at the bottom due to locally higher inflow misalignment, as seen in Figure 7.12(a) and as indicated by the direction of the streamlines in Figure 7.10. Higher axial velocity at the top is also contributed by the suction (and hence higher acceleration) on the wing leading edge. This non-uniform axial velocitydistributionimpliesthatthepropellerisoperatingwithlocallysmalleradvanceratioatthe bottomthanatthetop,thuswithhigherloadingatthebottom. Thebody-forcemodelproducedthe expectedresponseasshownbythedownstreamC p t contoursofFigure7.12(d),whichconfirmthe higherworkinputfromthepropellerbottom. (a) (b) (c) (d) (e) (f) Figure 7.12: Body-force model solution: upstream (station 1) contours of local advance ratio (a), andnormalizedswirlvelocity(b);anddownstream(station2)contoursofpressurecoefficient(c), totalpressurecoefficient(d),normalizedaxialvelocity(e),andnormalizedswirlvelocity(f). 142 Thesecondeffectisswirlflownon-uniformitycausedbytheoverallup-washexperiencedatthe propellerlocation. ThiscanbebetterseenintheupstreamswirlvelocitycontoursofFigure7.12(b), which indicate a co-swirl (red) region on the left and a counter-swirl (blue) region on the right. Thus, the up-wash compounds with the propeller rotation on the co-swirl (left) side and counter- acts it on the counter-swirl (right) side. In other words, the counter-swirl side of the propeller bladeisoperatingwithhigherbladeincidenceanglesandthusbearshigherloading,andviseversa for the co-swirl side. The body-force model produced the expected response by creating slightly higher work input (as reveled in downstreamC p t ), higher axial velocity rise, and higher pressure riseontherightsideofthepropeller. Inadditiontoaffectingtheloading,theincomingswirlnon-uniformityimpactsthedownstream swirlfield. Theresulteddownstream(station2)swirlvelocitydistributionofFigure7.12(f)shows that significantly higher swirl velocity is produced by the co-swirl (left) side of the propeller than the counter-swirl (right) side, with on the right side the swirl produced by the propeller almost cancelingouttheupstreamcounter-swirl. Thisunbalancedswirlvelocitydistributionatthedown- streamlocationinfluenceshowthejetthenpropagatesandmixes. Intermsofrelativeimportanceofthesetwoeffects(V x versusV q non-uniformities),thedown- stream total pressure contours show stronger gradients (more work input) near the bottom than on the right, thus indicating that the impact from the inflow axial velocity non-uniformity (more pronouncedatthebottom)ismoresignificantthanthatfromtheswirldistributionnon-uniformity (moreimportantontheright). The same set of axial flow contours as predicted by the actuator disk model are shown in Figure 7.13. Similar upstream non-uniform conditions can be seen from Figure 7.13(a) and (b) exceptthattheadvanceratiodistributionappearstobemoreleft-rightsymmetricalthanitwasfor the body-force model, due to a more symmetrical axial velocity profile (streamtube contraction). ThedownstreamC p t contourisrelativelyuniformwithslightlyhigherworkinputfromtheleft. As such, upstream non-uniformity in axial velocity and swirl distribution is not correctly reflected by theactuatordiskworkinput,duetoitsprescribedandfixedloading. Downstreamswirldistribution 143 (a) (b) (c) (d) (e) (f) Figure7.13: Actuatordiskmodelsolution: upstream(station1)contoursoflocaladvanceratio(a), andnormalizedswirlvelocity(b);anddownstream(station2)contoursofpressurecoefficient(c), totalpressurecoefficient(d),normalizedaxialvelocity(e),andnormalizedswirlvelocity(f). ismoresimilartothatpredictedbythebody-forcemodelduetoacircumferentiallyuniformswirl flow added by the advanced actuator disk model we are using, except for the higher downstream flowswirlnearthehubthatwasalreadyobservedintheisolatedpropellercase. 7.2.3 CoupledWingandPropellerPerformances Propellerperformanceparametersincludingthrustcoefficient,powercoefficientandaerodynamic efficiencyareplottedasafunctionofangleofattacksinFigure7.14forthetwopropulsormodels. These performance variations illustrate the impact of inflow misalignment and non-uniformity produced by the wing’s influence on the predicted propeller performance. Note that the thrust 144 value shown in the figure for the actuator disk results is computed by integrating the source field (andmatchthespecifiedthrustmodelvalue). Thisblown-liftwingtestcaseisanexampleofahighlycoupledpropulsor-airframesystem. It is worth noting that although the propeller is aligned with the freestream flow direction (i.e. far upstream outside of wing’s influence) when the wing is at an angle of attack of a =10 (due to the propeller being installed 10 pitching-down), the actual propeller inflow was observed to be aligned with the propeller axis at an a between 8 and 10 instead. Overall, there is thus a nearly 20 difference between the freestream direction and the propeller axis when the wing is at a = 10 , i.e., a nearly 20 up-wash at the propeller location due to the very high lift generated bythewing. Thus,itisclearthatthepresenceofthewingstronglyinfluence(viathelift/up-wash) the propulsor flow. Conversely, the propulsor jet is responsible for some lift-augmentation effect onthewingflap,andthereforethepropulsorflowinturninfluencesthewinganditslift. Ata = 10 wingangleofattack,thepropellerthrustcoefficientwhenpredictedbythebody- force model is undistinguishable (within 0.8%) from that for the isolated propulsor, a reflection of the fact that the propeller inflow is aligned with the propeller axis at this condition and that the major wing influence on the propeller is exerted via flow direction changes more than non- uniformities. As the angle of attack increases, thrust output from the actuator disk model remains constant despite flowfield changes, whereas the body-force model produces increasing thrust lev- els. A11%thrustincreaseisseenacrosstherangeofangleofattacksstudied(a 2 [ 10 ,15 ])as predictedbythebody-forcemodel. Thisistheeffectofinflowmisalignmentandwinginfluenceas previousexplained. Itwasshowninarelevantstudy[77]thatanisolatedpropellerproduceshigher thrust at a larger angle with respect to freestream: the presence of the high-lift wing exaggerates thiseffect. Thisincreaseofthrustoutputisnotcapturedbytheactuatordiskmodelsinceitsthrust valueisprescribedasaninputparameter. Propeller power and aerodynamic efficiency were only computed for the body-force model. The power coefficient is seen to initially increase with angle of attacks due to higher required powertomaintaintheincreasingloading. Thecompetingfactorofincreasingaverageaxialinflow 145 velocity(thusadvanceratio)takeseffectathigheranglesofattackwherethepowercoefficientstart to decrease. The maximum predicted shaft power occurs near a = 5 . The maximum propeller efficiencyisfoundtobe89%(whileitwas93%intheisolatedpropulsorcase), andtheefficiency curve has a relatively flat optimal region over a values between 0 and 5 , where the propeller experiencesmoderatemisalignmentandproduceshigherthrust. Figure 7.15 shows the lift, drag and pitching moment coefficients of the blown-lift wing as function of angle of attack, to illustrate how the blown-wing performance gets influenced by the propulsor models. Note that these coefficients are computed (defined) by integrating the surface stresses on the wing-and-pylon only (i.e. the propeller force contribution was not included), and then projecting along the freestream and normal directions. Similar trends are seen from the two modelsbutdiscrepanciesinmagnitudesexist. Byusingthebody-forcemodel,theblownwingpro- ducedonaverage5%largerc ` ,18%largerc d ,and18%largerc m magnitudethanwhattheactuator disk model predicts, with slightly higher slopes in lift and drag curves. Significant differences in pitchingmomentslopesareobservedbetween 10 and5 angleofattack,withmagnitudediffer- encesoverall,andthesepredictiondifferencescouldpotentiallyhaveanimpactonaircraftcontrol considerations. The difference between models diminishes as the angle of attack is reduced due tobetteralignmentofthepropellerplanewiththelocalflowdirection. Theobserveddifferencein forcesandmomentsaretheproductofthedifferentresponsesfromeachpropulsormodelatvarious anglesofattack,andtheresultofslightlydifferentjetsblowingoverthewingandbeingdeflected Figure 7.14: Propeller performance parameters as function of wing angle of attack: comparison betweenbody-forcemodel(BFM)andactuatordisk(AD)modelpredictions. 146 Figure 7.15: Wing force and moment coefficients as function of angle of attack: comparison be- tweenbody-forcemodel(BFM)andactuatordisk(AD)modelpredictions. bytheflap. Therefore,thepropulsormodel’sresponsetowinginfluencedoesmakealargeenough difference in wing performance (in addition to the difference already shown for propulsor perfor- mance) to justify the need for an advanced propulsor model that can adapt to non-uniform flow variations. From the comparison between the models for the blown-lift simulations, it is evident that the body-force model is able to adapt to the inflow variation caused by wing influence and produces aphysicalresponsethatbettermimicsrealpropelleroperation,whiletheactuatordiskmodelfails tocapturetheinteractionsandrequiresasetthrustlevel. 7.2.4 JetMixing AnimportantflowfeatureoftheDEPblown-liftwingisthemixingofthepropellerjets. Testcases with the full wing span and four propellers and zero angle of attack were created to visualize the jetmixingprocess. Figure 7.16 shows contours of C p t on axial planes placed at the streamwise stations when looking downstream, with vertical lines indicating the propeller center spanwise positions. The closely spaced propellers produce jets that gradually mix and coalesce into a jet sheet with an irregular pattern that moves slightly to the left as it travels downstream. This irregular jet profile is shaped by two mechanisms: the unbalanced swirl distribution (as shown in Figure 7.12(f)) and the blockage by the motor pylon as revealed in Figure 7.17. The stronger swirl region seen on 147 the propeller’s left side is deflected towards the top left by the motor pylon and impinges onto the main wing element. This in turn generates a leftward motion at the bottom of the motor pylon near station 3, clearly visible in Figure 7.17, which pushes the jet from that propulsor into the neighboring jet. This pattern is very similar to the jet-wake pattern reported in the wind tunnel experiments [76]. The leftward jet motion could be much more pronounced in a 3D full aircraft situation where the wing tip is unbounded. If all propellers on a wing rotate in this clockwise direction (as seen from the front), then the leftward motion induced by the outer-most propeller would compound with the motion induced by the tip vortex (from lower surface to upper surface) nearontheleftwingtip,butopposethetipvortexmotionneartherightwingtip. Figure 7.16 also shows the jet mixing process produced by the actuator disk model. Compar- ing the prediction made by the two models, minor difference in jet patterns can be seen due to differenceinswirldistributionsandpressurerise,butoverallasimilarmixingprocessisobserved. In addition, the body-force produced jets that appear to be more non-uniform probably due to its responsetoinflownon-uniformity. 7.3 Summary Body-force propulsor model has been implemented to represent propellers in high-fidelity CFD simulationsofablown-liftwing,atestcasewithstrongpropulsor-airframeinteractions,andcom- pared with results using an advancedactuator diskmodel andexperimental data. The overallgoal wastovalidatethebody-forcemodelforpropellerapplications,andinvestigatetheextenttowhich the propulsor models are able to capture the complex aero-propulsive interactions between the wingandthepropellers. Thebody-forcemodel’spredictioncapabilitieswasassessedfirstforanisolatedpropellerwith a uniform inflow aligned with propeller axis. The results show that the body-force model predicts the propeller’s c T l characteristic within 5% of experiment measurements and produces a real- istic trend for efficiency and power variation. The actuator disk model uses the thrust value as an 148 Figure 7.16: Jet mixing produced by the body-force model (left column) and by the actuator disk model(rightcolumn), asseenfromtotalpressurecontoursonaxialplanesatdifferentstreamwise stations. Verticallinesmarkthepropellerlocation,andpropellerrotationisclockwise. 149 Figure7.17: Totalpressurecontouronaxialplanenearthemotorpylonatstation3: theswirlflow isblockedbythepylon,inducingaleftwardmotionatthebottomofthepylon. input,andproducedsimilarflowfieldsasthebody-forcemodelincludingpressurerise,totalpres- surerise,andaxialvelocity. Radialsamplingoftheflowfieldshowedoverallsimilarpressureand velocity profiles, but higher swirl angle was observed near the hub from the actuator disk model duetoitssimpleestimationofcircumferentialforceswhichbecomesinaccurateatsmallerradius. The wing-propeller flow interaction was investigated via a blown-lift wing test cases. It was found that the propeller operates under inflow conditions with non-uniformities in both axial ve- locity (thus advance ratio) and swirl velocity, due to the large upwash at propeller location that is produced by the high-lift wing. The body-force model predicts more realistic response by pro- ducing higher pressure (total and static) rise and higher jet velocity in areas with locally lower advance ratio and upstream counter-swirl. An unbalanced downstream swirl distribution was ob- served as result of upstream non-uniformities which impacts how the propeller jets develop and mix. Acrossasweepoverarangeofangleofattacks,thebody-forcemodelproduceshigherthrust at higher angles of attack thanks to its ability to adjust to the varying inflow conditions caused by misalignment,whiletheactuatordiskwasunabletoreactasitmaintainsaconstantthrust. Theblown-liftwingperformancewasalsoaffectedbythedifferentpropulsormodelresponses. Onaverage,5%higherlift,18%higherdragand18%higherpitchingmomentwasobservedwhen usingthebody-forcemodel. Significantdifferencesinpitchingmomentslopesandmagnitudesare observed, an important consideration when selecting a model for estimating aircraft stability and 150 control. More jet flow deflection by the flapped wing was observed from cases using the body- force model, an effect that partially contributes to the higher lift and pitching moment. Similar jet mixingprocesseswereproducedbybothmodels,andthejet-wakepatternatflaptrailingedgewas foundtobesimilartothatreportedinwindtunnelexperiments. Overall, the body-force model was shown to be able to capture the wing-propeller interaction effectsbetterthantheactuatordiskmodel. Tobeabletomoreaccuratelypredictdistributedelectric propulsion(DEP)andblown-wingperformancewithstrongaero-propulsivecoupling,anadvanced propulsor model that is able to adapt to the varying local flow conditions is required, such as the body-forcemodelpresentedinthiswork. 151 Chapter8 Conclusions,FutureWorkandPotentialApplications 8.1 SummaryandConclusions With the goal of capturing the effects of inlet distortion on propulsor performance in numerical simulations of aircraft with propulsion system-airframe integration, a propulsor model was devel- opedusingabody-forcemodelingapproach. Theactualpropulsorbladesarereplacedbyavolume withasourcefieldthatproducestheequivalentpitch-averagedflowturning,blademetalblockage andlosses. Thecomputationofthebody-forcefieldutilizesaphysics-basedanalyticalformulation,which takes into account the blade camber and thickness distributions as well as local flow conditions. Thisallowsforanon-axisymmetricresponsetonon-uniformincomingflow. The proposed formulation behaves similar to thin airfoil theory when the number of blades is small, and as the number of blades approaches infinity, the formulation trends towards providing exact flow turning along the camber line. Thus, it is suitable for a wide range of applications, includingbothhigh-solidityturbomachineryandlow-solidityrotors. Thelocalnatureofthemodel makesitcompatiblewithgeneralCFDcodes,andsinceitrequiresonlyminimalknowledgeofthe bladegeometryitisalsowell-suitedfordesignstudies. The present model formulation builds upon the inviscid formulation introduced by Hall et al.[43], which demonstrated the capability to predict inviscid inflow distortion transfer. However, 152 Hall’s model exclusively captures the inviscid turning of the flow and cannot accurately represent propulsorperformancesuchasefficiency,asitdoesnotaccountforanylosses. Toaddresstheseshortcomings,weproposeanenhancedbody-forceformulationthatalsomod- elsthefollowing: (i)thebladeupwashinfluenceontheflowdeviationangle;(ii)themetalblockage effects;(iii)compressibilityeffects;and(iv)bladeprofilelossesandshocklosses. We implemented the propulsor model in the open-source flow solver ADflow and validated against experimental data for high-solidity cases including transonic NASA Rotor 67 and NASA Stage 35, and low-speed low solidity fan stage TF8000. The body-force model, thanks to its blockage component, captures the physically correct choking behavior and predicts the choking massflowwithina1%errorcomparedtoexperimentalmeasurements. The proposed model achieves accurate predictions for stand-alone rotors near the on-design operating condition for both pressure rise and efficiency, indicating the correct amount of inviscid loadandlosses. Tocaptureoff-designtrendsduetoincreasedlosses,theoff-designlosscomponent mustbeapplied,requiringareferencedesignoperatingpointtobechosen. Withallmodelcompo- nentsapplied,thetotalpressureandefficiencyarewellpredictedforvariousrotationalspeeds. The effectofreducedefficiencyathigherRPMsduetohighertipMachnumberswascapturedbecause of the inclusion of a shock loss. The flow angle correction significantly improved prediction for low-solidityfanstage,fixingthepreviouslyunder-predictedpressurerise. Next, The proposed model were tested with various total pressure inlet distortion including a simpleradialdistortion,a120degreesectorofcircumferentialdistortion,andaboundarylayerin- gestiontypedistortion. Three-dimensionalflowredistributionakintothatobservedinfull-annulus high-fidelity bladed simulations and experiments were accurately reproduced by the body-force model. This upstream redistribution results in a non-uniform operating condition at propulsor in- let,thebody-forcemodelwasabletoprovidenon-uniformresponseandprovideaccuratedistortion transfer through the rotor, resulting in downstream flow property distribution well correlated with thatshowninbothURANSandexperimentalresults. Furthermore,thebody-forcemodelcorrectly predictedthedecreaseinfanefficiencyresultingfromtheinletdistortion. 153 Inordertoextendthemodeltoapplyforfreetiprotorssuchaspropellers,addtionalcorrections including a wake correction and a tip loss correction are proposed for the body-force model. Test casesusingthegeometryoffoursmall-scalepropellersyieldedpromisingresultsthatcloselyalign with experimental data for thrust coefficient, power coefficient, and efficiency. This proves that body-force based model can be used for low-solidity rotors that relies on blade element based models pre-dominantly. This demonstrates that the body-force-based model can be effectively appliedtorepresentlow-solidityrotors,suchaspropellers,whichhavehistoricallypredominantly reliedonbladeelement-basedmodels. Last, thebody-forcepropulsormodelwasimplementedtorepresentpropellersinhigh-fidelity CFD simulations of a blown-lift wing, a test case with strong propulsor-airframe interactions. Overall, the body-force model was shown to be able to capture the wing-propeller interaction effectsbetterthantheactuatordiskmodel. In conclusion, the proposed body-force model produces encouraging results. In its current state, it should be capable of predicting the physically-correct trend of local efficiency variations whenusedforfull-aircraftCFDwithengineinletdistortion,andprovideinsightsintotheeffectof inletflownon-uniformityonsystemlevelperformance. 8.2 Contributions The present work improves the state-of-art in body-force propulsor modeling for CFD via the followingcontributions: • Introduced formulation for a body-force model that represents inviscid loading, efficiency variation, and blockage effects concurrently while using an analytical formulation and car- riedoutvalidationsagainstexperiments. • Validatedtheabilityoftheinviscidformulationtocapturethemechanismtoproperlyrepre- sentinviscidloadingtrendsinrelationtobladespacing(solidity),whileconcurrentlyidenti- fyingtheneedforaflowanglecorrectiontoaccountforincreasedupwasheffectswithblade 154 spacing, particularly in bridging the gap between high solidity rotors (turbomachinery) and lowsolidityapplications. • Demonstratedthattheapplicationofblockagemodelingwithouranalyticalbody-forcefor- mulationleadstocorrectchokingbehaviorandpredictsthechokingmassflowrateconsistent withexperimentalfindings. • Applied the proposed body-force model to propulsor test cases with various type of inlet distortion and confirmed its ability to accurately predict distortion induced upstream redis- tribution effects, distortion transfer across fan stage, and flow distortion’s impact on overall performance,consistentwithresultsfrombothhigh-fidelityblade-resolvedURANSsimula- tionsandexperiments. • Confirmed the applicability of the body-force model to low solidity external flow (free-tip) applications,suchaspropellers,providinganaccuraterepresentationofperformancetrends, including thrust coefficients, power coefficients, and efficiency. Necessary adjustments to themodelformulationareproposedtoaccountforwakeeffects. 8.3 FutureWorks 8.3.1 ModelImprovements ReducedModelInput Theminimumrequiredinputofabladecambersurfaceposesadifficulty for applications when such information is not known or available. It is possible to derive An alternativeroutinecouldbedevelopedtogeneratebladegeometrybasedonidealloaddistributions, requiring only simple design inputs such as propulsor dimensions, required thrust, and rotational speeds. This approach would enable much faster and easier application of the body-force model whenthepropulsorgeometryisnotknown. 155 Loss Model Our model is based on a simplified loss model that relies on a known relationship of momentum thickness as a function of Reynolds number for flat plate boundary layers, coupled withalossbucketmodel. ThislossmodelprovidesreasonablepredictionsacrossvaryingReynolds numbers. However, to further improve the model’s reliability, a more rigorous loss model might beneeded. Theboundarylayeronapropulsorbladeexperiencesvaryinglevelsofadversepressuregradi- ents,leadingtodifferentstatesofboundarylayerthickness. Furthermore,differencesexistbetween highly loaded turbomachines that experience larger adverse pressure gradients, and lightly loaded propellers. One potential improvement could involve modifying the momentum thickness formulation to account for pressure gradients or the inviscid loading generated by the body-force, leading to a more comprehensive representation of loss. Currently, the extra losses are accounted for via the off-designorlossbucketmodel. Alternatively,wecouldseektoleverageexistinglossdatafromempiricalsourcesorlow-order aerotoolssuchasXfoilandMESIS.Pazireh[45]proposedalossformulationbasedonanartificial neural network trained with a dataset of cascade computations, which resulted in very good loss prediction for turbomachines. However, it may not be generalizable to low solidity applications likepropellers. Thebody-forcemodeliscapableofproducingtheentireperformancetrendsacrosstheoperat- ingrange,providedoneoftheoperatingpointsiscorrectlycalibrated. Simplecalibrationstepsfor thelossmodelcanbeincorporatedtocalibratelossatasingleon-designoperatingpoint. 8.3.2 Applications DesignandOptimization The body-force model can serve an important role in the design of a propulsor. Hall [43] demonstrated the use of the body-force model in suggesting possible design features for a distortion-tolerant stator. Given that the body-force model generates performance predictionsbasedontherotorandstatorgeometry,thesegeometriescanbeeasilymodifiedduring 156 runtime without necessitating any grid modifications. This feature enables the simultaneous opti- mization of airframe and propulsor geometries, which can offer novel insights into the design of aero-propulsiveaircraft. Turbines There’s no technical difficulty that prevents the application of body-force model to represent turbines. This is because the energy source term can change its sign based on the flow’s reaction to the specified camber surface at a given rotation speed and direction. As a result, a turbine, which acts as an energy sink, can be easily represented by the body-force model, given a camber surface and its rotational speeds. The wake models suggested in Chapter 6 can also be applied to wind turbines. As discussed in Section 1.2, the current modeling of wind turbines primarily relies on models based on the blade element theory. The body-force model provides an alternative approach with the ability to respond to non-uniform inflow conditions, making it a valuabletoolforthedesignandoptimizationofwindfarms. Full Engine The body-force model offers a way to do full engine modeling with low compu- tational cost, with thermodynamic model required to represent the combustion chamber. A sig- nificant value proposition of the body-force model lies in its ability to investigate the propagation of distortion in multi-stage turbomachines. This model’s adaptability to local flow variations en- ablesittoeffectivelyaddressthecomplexitiesinvolvedinmodelingmulti-stagesequences,thereby providinginvaluableinsights. 157 References [1] A. Uranga, M. Drela, D. K. Hall, and E. M. Greitzer, “Analysis of the aerodynamic benefit fromboundarylayeringestionfortransportaircraft,”AIAAJournal,vol.56,no.11,pp.4271– 4281,2018. [2] M. Drela, “Development of the D8 transport configuration,” in 29th AIAA Applied Aerody- namicsConference,p.3970,2011. [3] J. Welstead and J. L. Felder, “Conceptual design of a single-aisle turboelectric commercial transport with fuselage boundary layer ingestion,” in AIAA 2016-1027, AIAA SciTech, 54 nd AerospaceSciencesMeeting,4–8Jan2016,SanDiego,CA,2016. [4] A. Uranga, M. Drela, E. Greitzer, N. Titchener, M. Lieu, N. Siu, A. Huang, G. M. Gatlin, andJ.Hannon,“Preliminaryexperimentalassessmentoftheboundarylayeringestionbenefit for the D8 aircraft,” in AIAA SciTech, 52 nd Aerospace Sciences Meeting, 13–17 Jan 2014, NationalHarbor,MD,2014. [5] V. J. Fidalgo, C. A. Hall, and Y. Colin, “A study of fan-distortion interaction within the nasa rotor67transonicstage,”JournalofTurbomachinery,vol.134,052012. [6] C. P. Arroyo, J. Dombard, F. Duchaine, L. Gicquel, B. Martin, N. Odier, and G. Staffelbach, “Towards the large-eddy simulation of a full engine: Integration of a 360 azimuthal degrees fan,compressorandcombustionchamber.parti: Methodologyandinitialisation,”Journalof theGlobalPowerandPropulsionSociety,no.Specialissue,p.133115,2021. [7] H. Glauert, The Elements of Aerofoil and Airscrew Theory. Cambridge university press, 1983. [8] T. Xie and A. Uranga, “Propeller representation in full-vehicle cfd: Actuator disk versus body-force modeling approches,” in Proceedings of the Eleventh International Conference onComputationalFluidDynamics,ICCFD11,Hawaii,UnitedStates,on11-15July,2022. [9] M.Drela,“QPropformulation,”pp.1–14,2006. [10] H.Glauert,“Airplanepropellers,”Aerodynamictheory,1935. [11] L.Prandtl,ApplicationsofModernHydrodynamicstoAeronautics. USGovernmentPrinting Office,1925. [12] M.O.Hansen,AerodynamicsofWindTurbines. Routledge,2015. 158 [13] “OpenFAST.”https://github.com/OpenFAST/openfast. Accessed: 2023-04-15. [14] C. N. H. Lock and R. Pankhurst, “Strip theory method of calculation for airscrews on high- speedaeroplanes,”tech.rep.,HMStationeryOffice,1945. [15] E. J. Alvarez and A. Ning, “High-fidelity modeling of multirotor aerodynamic interactions foraircraftdesign,”AIAAJournal,vol.58,no.10,pp.4385–4400,2020. [16] K. Shaler, E. Branlard, A. Platt, and J. Jonkman, “Preliminary introduction of a free vortex wakemethodintoopenfast,”inJournalofPhysics: ConferenceSeries,vol.1452,p.012064, IOPPublishing,2020. [17] F. Le Chuiton, “Actuator disc modelling for helicopter rotors,” Aerospace Science and Tech- nology,vol.8,no.4,pp.285–297,2004. [18] I.FejtekandL.Roberts,“Navier-stokescomputationofwing/rotorinteractionforatiltrotor inhover,”AIAAjournal,vol.30,no.11,pp.2595–2603,1992. [19] I. Ammara, C. Leclerc, and C. Masson, “A viscous three-dimensional differential/actuator- diskmethodfortheaerodynamicanalysisofwindfarms,”J.Sol.EnergyEng.,vol.124,no.4, pp.345–356,2002. [20] R.G.RajagopalanandJ.B.Fanucci,“Finitedifferencemodelforverticalaxiswindturbines,” JournalofPropulsionandPower,vol.1,no.6,pp.432–436,1985. [21] J. N. Sørensen and A. Myken, “Unsteady actuator disc model for horizontal axis wind tur- bines,”JournalofWindEngineeringandIndustrialAerodynamics,vol.39,no.1-3,pp.139– 149,1992. [22] N. Troldborg, Actuator Line Modeling of Wind Turbine Wakes. PhD thesis, Technical Uni- versityofDenmark,2009. [23] W. Z. Shen, J. H. Zhang, and J. N. Sørensen, “The actuator surface model: A new navier– stokes based model for rotor computations,” Journal of solar energy engineering, vol. 131, no.1,2009. [24] C.-H. Wu, “A General Though-Flow Theory of Fluid Flow With Subsonic or Supersonic Velocity in Turbomachines of Arbitrary Hub and Casing Shapes,” Tech. Rep. March 1951, 1951. [25] C. Hirsch and J. Denton, “Through Flow Calculations in Axial Turbomachines,” in NATO AdvisoryGroupforAerospaceResearchandDevelopmentReport,no.175,1982. [26] J.L.H.Smith,“TheRadial-EquilibriumEquationofTurbomachinery,”Journalofengineer- ingforpower,1966. [27] R.A.Novak,“StreamlineCurvatureComputingProceduresforFluid-FlowProblems,”Jour- nalofEngineeringforGasTurbinesandPower,vol.89,p.478,oct1967. 159 [28] H. Marsh, “A Digital Computer Program for the Through-flow Fluid Mechanics in an Arbi- trary Turbomachine using a Matrix Method A Digital Computer Program for the Through- flowFluidMechanicsinanArbitraryTurbomachineusingaMatrixMethod,”no.3509,1968. [29] C. H. Hirsch and G. Warzee, “A Finite-Element Method for Through Flow Calculations in Turbomachines,”no.September1976,pp.403–414,1976. [30] J. D. Denton, “Throughflow Calculations for Transonic Axial Flow Turbines,” Journal of EngineeringforPower,vol.100,no.2,p.212,1978. [31] R.A.NovakandR.M.Hearsey,“ANearlyThree-DimensionalIntrabladeComputingSystem forTurhomaohinery,”no.March1977,pp.154–166,1977. [32] J.J.Adamczyk,“ModelEquationforSimulatingFlowsinMultistageTurbomachines,”Nasa, 1985. [33] J. D. Denton and W. N. Dawes, “Computational Fluid Dynamics for Turbomachinery De- sign,”ProceedingsoftheInstitutionofMechanicalEngineers,PartC:JournalofMechanical EngineeringScience,vol.213,no.2,pp.107–124,1998. [34] F. E. Marble, Three-Dimensional Flow in Turbomachines, pp. 83–166. Princeton University Press,1964. [35] A.Spurr,“Thepredictionof3dtransonicflowinturbomachineryusingacombinedthrough- flow and blade-to-blade time marching method,” International Journal of Heat and Fluid Flow,vol.2,no.4,pp.189–199,1980. [36] J.D.Denton,“TheCalculationofThree-DimensionalViscousFlowThroughMultistageTur- bomachines,”JournalofTurbomachinery,vol.114,no.1,p.18,1990. [37] S. V. Damle, T. Q. Dang, and D. R. Reddy, “Throughflow method for turbomachines appli- cableforallflowregimes,”JournalofTurbomachinery,vol.119,p.256262,041997. [38] S. J. Gallimore, “Viscous Throughflow Modelung of Axial Compressor Bladerows Using a TangentialBladeForceHypothesis,”ProceedingsoftheASMETurboExpo,vol.1,no.Octo- ber,pp.662–670,1997. [39] J. F. Simon, Contribution to Throughflow Modelling for Axial Flow Turbomachines. PhD thesis,UniversityofLiege,2007. [40] Y.Gong, A Computational Model for Rotating Stall and Inlet Distortion in Multistage Com- pressors. PhDthesis,MassachusettsInstituteofTechnology,1999. [41] A. Peters, Z. S. Spakovszky, W. K. Lord, and B. Rose, “Ultrashort nacelles for low fan pres- sureratiopropulsors,”Journalofturbomachinery,vol.137,no.2,p.021001,2015. [42] W. Thollet, G. Dufour, X. Carbonneau, and F. Blanc, “Body-Force Modeling for Aerody- namic Analysis of air Intake-Fan Interactions,” International Journal of Numerical Methods forHeatandFluidFlow,vol.26,no.7,pp.2048–2065,2016. 160 [43] D.K.Hall,E.M.Greiter,andC.S.Tan,“Analysisoffanstagedesignattributesforboundary layeringestion,”JournalofTurbomachinery,vol.139,pp.071012–1–10,July2017. [44] H. D. Akaydin and S. A. Pandya, “Implementation of a Body Force Model in OVERFLOW forPropulsorSimulations,”35thAIAAAppliedAerodynamicsConference,2017. [45] S. Pazireh and J. Defoe, “A new loss generation body force model for fan/compressor blade rows: Application to uniform and non-uniform inflow in rotor 67,” Journal of Turbomachin- ery,vol.144,no.6,p.061005,2022. [46] S.BaralonandU.Hall, “ViscousThroughflowModelingofTransonicCompressorsUsinga Time-Marching Finite Volume Solver,” Thirteenth International Symposium on Air Breath- ingEngines(ISABE),pp.502–510,1997. [47] A.SturmayrandC.Hirsch,“ShockRepresentationbyEulerThroughflowModelsandCom- parison with Pitch-Averaged Navier-Stokes Solutions.” VRIJE UNIV BRUSSELS (BEL- GIUM)DEPTOFFLUIDMECHANICS,1999. [48] G. Persico and S. Rebay, “A Penalty Formulation for the Throughflow Modeling of Turbo- machinery,”ComputersandFluids,vol.60,pp.86–98,2012. [49] J. H. Horlock, “Actuator disk theory-discontinuities in thermo-fluid dynamics,” New York, 1978. [50] W. Joo and T. Hynes, “The simulation of turbomachinery blade rows in asymmetric flow usingactuatordisks,”1997. [51] D. Hall, Analysis of Civil Aircraft Propulsors with Boundary Layer Ingestion. PhD thesis, MassachusettsInstituteofTechnology,February2015. [52] J.D.Denton, “The1993IGTIScholarLecture: LossMechanismsinTurbomachines,” Jour- nalofTurbomachinery,vol.115,no.4,p.621,1993. [53] Y.Chen,“TopicsonConceptualDesignoftheD8Aircraft: FuelBurnUncertaintyEstimate and BLI Propulsor Loss Modeling,” Master’s thesis, Massachusetts Institute of Technology, 2017. [54] M.Drela,FlightVehicleAerodynamics. MITpress,2014. [55] H.CrevelingandR.Carmody,“Computerprogramforcalculatingoff-designperformanceof multistageaxial-flowcompressors,”NASACR-72427,1967. [56] C.A.Mader,G.K.Kenway,A.Yildirim,andJ.R.Martins,“Adflow: Anopen-sourcecompu- tational fluid dynamics solver for aerodynamic and multidisciplinary optimization,” Journal ofAerospaceInformationSystems,vol.17,no.9,pp.508–527,2020. [57] A. Yildirim, G. K. W. Kenway, C. A. Mader, and J. R. R. A. Martins, “A jacobian-free approximate newtonkrylov startup strategy for rans simulations,” Journal of Computational Physics,p.108741,Nov.2019. 161 [58] “APC propeller technical information.” http://www.apcprop.com/ technical-information. Accessed: 2023-04-15. [59] J. B. Brandt, “Small-scale propeller performance at low speeds,” Master’s thesis, University ofIllinoisatUrbana-Champaign,2005. [60] A. Strazisar, J. R. Wood, M. D. Hathaway, and K. L. Suder, “Laser Anemometer Measure- mentsinaTransonicAxialFlowFanRotor,”Tech.Rep.July,NASA,1989. [61] L.ReidandR.D.Moore,“Designandoverallperformanceoffourhighlyloaded,highspeed inletstagesforanadvancedhigh-pressure-ratiocorecompressor,”tech.rep.,1978. [62] N. M. Siu, “Evaluation of Propulsor Aerodynamic Performance for Powered Aircraft Wind TunnelExperiments,”Master’sthesis,MassachusettsInstituteofTechnology,2015. [63] A.Uranga,M.Drela,E.Greitzer,D.Hall,N.Titchener,M.Lieu,N.Siu,A.Huang,G.Gatlin, and J. Hannon, “Boundary layer ingestion benefit of the D8 aircraft,” AIAA Journal, vol. 55, no.11,pp.3693–3708,2017. [64] W. S. Cunnan, W. Stevans, and D. C. Urasek, “Design and performance of a 427-meter-per- second-tip-speedtwo-stagefanhavinga2.40pressureratio,”tech.rep.,1978. [65] D. C. Urasek, W. T. Gorrell, and W. S. Cunnan, “Performance of two-stage fan having low- aspect-ratiofirst-stagerotorblading,”tech.rep.,1979. [66] F. S. Weinig, Theory Of Two-Dimensional Flow Through Cascades, pp. 13–82. Princeton UniversityPress,1964. [67] N.A.Cumpsty,CompressorAerodynamics. LongmanScientific&Technical,1989. [68] S. Lieblein, “Incidence and deviation-angle correlations for compressor cascades,” Journal ofBasicEngineering,1960. [69] A.H.Stenning,“Inletdistortioneffectsinaxialcompressors,”JournalofFluidsEngineering, vol.102,pp.7–13,031980. [70] J.LongleyandE.Greitzer,“Inletdistortioneffectsinaircraftpropulsionsystemintegration,” in Steady and transient performance prediction of gas turbine engines, Advisory Group for AerospaceResearchandDevelopment(AGARD),1992. [71] C. Reid, “The response of axial flow compressors to intake flow distortion,” in Turbo Expo: PowerforLand,Sea,andAir,031969. [72] E.J.Gunn,S.E.Tooze,C.A.Hall,andY.Colin,“Anexperimentalstudyoflosssourcesina fanoperating withcontinuousinletstagnation pressuredistortion,” Journal of Turbomachin- ery,vol.135,no.5,p.051002,2013. [73] E. Gunn and C. Hall, “Aerodynamics of boundary layer ingesting fans,” in Turbo Expo: PowerforLand,Sea,andAir,vol.45578,p.V01AT01A024,AmericanSocietyofMechani- calEngineers,2014. 162 [74] C. S. Mistry and A. Pradeep, “Investigations on the effect of inflow distortion on the perfor- manceofahighaspectratio,lowspeedcontrarotatingfanstage,” inTurboExpo: Powerfor Land, Sea, and Air, vol. 55225, p. V06AT35A005, American Society of Mechanical Engi- neers,2013. [75] C.Hah, D. C.Rabe, T. J. Sullivan, and A. R.Wadia, “Effects of InletDistortion on theFlow FieldinaTransonicCompressorRotor,”JournalofTurbomachinery,vol.120,pp.233–246, 041998. [76] T. Long, “An Experimental Investigation of Propeller Blown-Flap Airfoils,” Master’s thesis, MassachusettsInstituteofTechnology,2021. [77] T.C.A.StokkermansandL.L.M.Veldhuis,“Propellerperformanceatlargeangleofattack applicabletocompoundhelicopters,”AIAAJournal,vol.59,no.6,pp.2183–2199,2021. [78] C.B.Courtin,AnAssessmentofElectricSTOLAircraft. PhDthesis,MassachusettsInstitute ofTechnology,2019. [79] T.Long,M.Drela,J.Hansman,andC.Courtin,“Parametricdesignstudyforablownflapped wing,”inAIAAAviation2021Forum,p.2871,2021. [80] D.Spence,“Theliftcoefficientofathin,jet-flappedwing,”ProceedingsoftheRoyalSociety of London. Series A. Mathematical and Physical Sciences, vol. 238, no. 1212, pp. 46–68, 1956. 163 Appendices A Derivations A.1 EnergySourceTerm Given that the blade row is driven by the shaft, local shaft power per unit mass can be computed fromthelocalbladebody-forceas ˙ w=W (f⇥ r)=(W r)f q =f·U (A.1) Euler’s turbine equation dictates that contribution of the energy equation should equal to the shaft power. Thus Equation (A.1) need to be added to the right hand side of the energy equation to representtheshaftworkinput. A.2 EntropyGeneration Next,wewillshowthederivationoftheentropyequationwithbody-forcecontribution. Assuming that the flow is steady and adiabatic. The equation for total energy, total enthalpy and momentum withthebody-forcecontributionarewrittenas r V·— (e+ 1 2 V 2 )= — ·(pV)+— ·(V·t )+r f·U, (A.2) 164 V·— h t = 1 r — ·(V·t )+f·U (A.3) and V·— V= 1 r ( — p+— ·t )+f (A.4) ThekineticenergyequationcanbederivedbyV·Equation(A.4) V·— ( 1 2 V 2 )= 1 r ( V·— p+V·(— ·t ))+f·V (A.5) Subtract the kenetic energy Equation (A.5) from total energy Equation (A.2) gives the internal energyequation V·— e= p r — ·V+ 1 r t :— V f·W (A.6) Duetomassconservation,wecandothefollowingsubstitutioninEquation(A.6) — ·V= 1 r V·— r =r V·— v (A.7) TheGibbsequationiswrittenas T— s=— e+p— v (A.8) V·Equation(A.8)andthensubstituetheinternalenergyEquation(A.6)givestheentropyequation TV·— s= 1 r t :— V f·W (A.9) A.3 BlockageForceTerm Pressureforceonlateralwallintheaxialdirection F p x = 2p~ n·~ e x D lD z= 2psinq D lD z= 2p ∂h ∂x |{z} tanq cosq D l | {z } D x D z= 2p ∂h ∂x D xD z (A.10) 165 Expressblockageanditsgradientintermsofwallsurfacegradient b= L h(x) L =1 h(x) L (A.11) db dx = 1 L dh dx (A.12) Thus, F p x =(2L)p db dx D xD z (A.13) Distributetheforceuniformlyinthepitch F p x =(2L)p db dx 1 2Lb D xD zD y= p b db dx D xD zD y (A.14) Due to pitch-averaging, blockage term only has a gradient in the axial direction thus the force per volumetobeaddedintothedifferentialformofgoverningequationis: f p x = p b — b (A.15) V ` h(x) x y ~ n L o ✓ Figure8.1: Converging-divergingnozzle 166 B CompleteBody-ForceFormulation Thegoverningequationswiththecompletebody-forcemodelwithblockageeffectsarewrittenas follows ∂r ∂t +— ·(r V)= 1 b (r V·— b), (B.1) ∂ ∂t (r V)+— ·(r VV+p ¯ ¯ I ¯ ¯ t )=r f+ 1 b ( ¯ ¯ t r VV)·— b, (B.2) ∂ ∂t (r e t )+— ·(r Vh t ¯ ¯ t ·V)=r f·U+ 1 b ( ¯ ¯ t ·V r Vh t )·— b (B.3) wherebisablockagefactor b = 1 t (2p r/B)|n q | , (B.4) t is local blade thickness, |n q | is the circumferential component of the local blade normal vector, r is local radius from propulsor axis, B is number of blades. The source terms involving blockage factor b contribute to the blockage effects and the sources terms with body-force f representing inviscidloadingandviscouslossesthataremodeledseparatelyas f=f n +f v . (B.5) Theinviscidbody-forcemagnitudeisevaluatedas f n ⌘ K c 2pd 1 2 W 2 (2p r/B)b|n q | ˆ e n , (B.6) where blockage factor b is included in the denominator to account for the change of the pitch distancesduetoblockage,K c isacompressibilitycorrectionfactorthattakesthefollowingform K c ⌘ min 1 p |1 M 2 r | , 1.66 ! , (B.7) 167 localdeviationangleisfoundvia d = sin 1 ✓ W·ˆ n ||W|| ◆ +D d , (B.8) whereW isvelocityrelativetotheblade, ˆ nisthelocalbladenormalvector. EquationB.8represents thelocalrelativeflowdeviationfromthebladesurface,andtheflowanglecorrectionsarereflected viatermD d whichtakesdifferentformsforductandfree-tiprotors. In the case of ducted rotors where the effects of the tip vortex system can be neglected, only theimpactoftheboundcirculationneedstobeconsidered,thus D d =D d G ⌘ 0.5K t K c a loc ✓ min(s/c,4.0) 4.0 ◆ (B.9) where a loc =d + 1 2 (k k TE ), (B.10) inwhichd isthelocaldeviationangle,k islocalbladeangle,andk TE isbladeangleatthetrailing edge. K t isacalibrationfactorwithadefaultvalueof1. Inthecaseoffreetiprotors(propeller): D d =D d G +D d i +D d tip +D d lim (B.11) whereD d G istheflowanglecorrectionduetoboundcirculationasshowninEquation(B.9),D d i is acorrectionaccountingforwakeeffects D d i ⌘ BK c a loc 4(r/c) , (B.12) D d tip accountfortiplossesthatreducestheeffectiveloading D d tip =(F 1)d . (B.13) 168 where F⌘ 2 p cos 1 (e f ), (B.14) and f = B 2 ✓ 1 r r tip ◆ r 1+ 1 l 2 (B.15) and l is a local advance ratio defined as l ⌘ V ax /(W r tip ), and lastlyD d lim is a load limiting cor- rectionforpropellersatsmalladvanceratioorstaticthrustcondition,i.e.,V ax <V ax min D d lim =tan 1 ✓ W r V ax min ◆ tan 1 ✓ W r V ax ◆ (B.16) Thedirectionoftheinviscidforceisfoundby ˆ e n ⌘ f n ||f n || = sign(W·ˆ n) W⇥ (W⇥ ˆ n) ||W⇥ (W⇥ ˆ n)|| , (B.17) such that the force is normal to the relative flow and the radial component of the blade surface normal is taken account of. The sign function ensures that the normal force always points in the directionthatreduceslocalflowdeviation. Theviscousbody-forceaccountingfortheprofilelosses iscomputedas f v ⌘ K v ( 1 2 W 2 ) (2p r/B)b|n q | ˆ e v . (B.18) where K v is a dissipation factor accounting for profile losses for both on-design condition and off-designdeviation K v =K v o +K v 1 (d d ref ) 2 (B.19) where K v o = 4.0C D + 3.0C D (K c pd ) 2 , (B.20) 169 and the boundary dissipation factor, C D , is functions of Reynolds number based on momentum thickness. ForlaminarflowC D isevaluatedas q = 0.664x (Re x ) 1/2 C D =0.173Re 1 q (B.21) andforturbulentflowas q = 0.016x (Re x ) 1/7 C D =0.0056Re 1/6 q (B.22) Anadditionalviscousbody-forcecomponentaccountingforshocklossesneedtobeaddedfor transonicapplication: f v,shock ⌘ Tc v c 2g (g 1) 3(g +1) 2 (max(M r ,1) 2 1 3 , (B.23) Anadditionallosscomponentneedtobeaddedforfree-tiprotorsaccountforinducedcompo- nentofdrag: f v,i = K c (2pd )sin(D d i ) 1 2 W 2 (2p r/B)b|n q | (B.24) Thedirectionˆ e v oftheviscousbody-forceisalwaysoppositetothelocalrelativeflowdirection, i.e., ˆ e v ⌘ f v ||f v || = W ||W|| , (B.25) sothatlossesareintroducedwithoutaffectingtheflowturning. 170 C BlockageModelValidation As shown in Figure 8.2, a choked 2D converging-diverging nozzle was utilized as a validation test case for the blockage model. A case with a straight duct was created, with the blockage model applied to simulate a converging-diverging nozzle. The pressure variation along the duct’s centerlinewasthencomparedwithanothercasethatfeaturedtheactualductgeometrybutwithout theblockagemodel. Both cases had identical total pressure at the inlet and the same static pressure at the outlet. The comparison showed a good agreement between the pressure variations in the two cases. This verifiesthattheblockagemodelisfunctioningasanticipatedinsimulatingaconverging-diverging nozzle. Discrepancieswereobservedattheinlet,whichcanbeattributedtotheinteractionbetween thedifferentductgeometryandtheinletboundaryconditionthatsetsthevelocitydirection. Blockage Model Blockage model Actual Nozzle Pressure contour Pressure distribution along centerline x pressure (Pa) Acutal Nozzle Figure8.2: Validationoftheblockagemodel. 171
Abstract (if available)
Abstract
This thesis presents a body-force based propulsor model that replaces the propulsor blades with a source volume in computational fluid dynamics (CFD) simulations, and that effectively repro- duces the equivalent flow turning, work input, blade blockage, and losses. The primary motivation behind this work is to establish a propulsor model capable of interacting with local flow non- uniformities and producing physically accurate responses when implemented in aero-propulsive coupled CFD simulations. The model formulation is fundamentally physics-based, thus does not require prior knowledge of propulsor performance or empirical input, and only necessitates the blade’s geometrical information, such as a camber surface, thickness distribution and number of blades. An inviscid formulation for the body-force was previously developed in the literature and demonstrated its potential for predicting the inviscid distortion transfer effects. However, with- out considering losses, blade metal blockage, and blade solidity effects, the inviscid model cannot accurately represent a propulsor’s performance.
In this thesis, an improved formulation with additional physical models is proposed and is shown to properly predict both propulsor work input and efficiency. Blockage models are intro- duced as additional source terms to account for the blockage effect produced by the blade row onto the passing flow. Corrections are introduced into the inviscid load formulation to address com- pressibility and solidity effects and ensure compatibility with the blockage model. Loss terms are included to model blade profile losses and shock losses.
The proposed model is implemented within the open-source flow solver ADflow and validated against experimental data for two highly loaded turbomachinery fan stages (NASA Rotor 67 and Stage 35) and a low speed low-solidity fan stage (TF8000), operating in clean flow. The total pressure ratio shows good agreement and the choking mass flow was captured within a 1% error. The isentropic efficiency and its off-design trends were qualitatively well predicted.
The body-force model was applied to the NASA rotor 67 and TF8000 fan stage test cases with various types of inlet total pressure distortion profiles. The interactions between the fan stages rep- resented via the body-force model and the inlet total pressure distortion were investigated. Comparisons with full-annulus bladed CFD results demonstrated that the body-force model accurately predicted the distortion transfer behavior and local variations in loading and efficiency. The time-averaged flow field of high-fidelity URANS results were essentially reproduced by the body-force model with much lower computational cost and the 1-2% drop in fan efficiency due to boundary layer ingestion (BLI) were accurately predicted.
Additional model formulations were proposed, including wake-induced corrections, tip loss correction, and a load limit correction, to improve the representation of free-tip propulsors such as aircraft propellers using the body-force model. The model was applied to small-scale propellers of various specifications and validated against wind tunnel test data. The body-force model success- fully generated the non-axisymmetric responses of the propellers to skewed inflow. Moreover, it demonstrated performance trends consistent with experimental results.
Finally, the body-force model was utilized to represent propellers in the RANS simulation of a blown-lift wing test case where strong propeller-airframe interactions are present. Within a blown-lift wing system, distributed propellers blow onto a flapped wing to produce high lift. Comparisons were made between results produced by the body-force model and those generated by an advanced actuator disk model. The body-force model demonstrated superior performance in capturing the wing-propeller interaction effects, highlighting its effectiveness in dealing with complex interaction scenarios that are beyond the capabilities of more traditional models.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Aerodynamics at low Re: separation, reattachment, and control
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
Tracking and evolution of compressible turbulent flow structures
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Large eddy simulations of turbulent flows without use of the eddy viscosity concept
PDF
Modeling and analysis of propulsion systems and components for electrified commercial aircraft
PDF
Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows
PDF
Model based design of porous and patterned surfaces for passive turbulence control
PDF
Numerical and experimental investigations of ionic electrospray thruster plume
PDF
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Multiscale modeling of cilia mechanics and function
PDF
Passive rolling and flapping dynamics of a heaving Λ flyer
PDF
Biomimetics and bio-inspiration for moderate Reynolds number airfoils and aircraft
PDF
Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection
PDF
Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves
PDF
Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders
PDF
Large eddy simulations of laminar separation bubble flows
PDF
Accuracy and feasibility of combustion studies under engine relevant conditions
PDF
Thermal modeling and control in mobile and server systems
Asset Metadata
Creator
Xie, Tianbo
(author)
Core Title
Development and applications of a body-force propulsor model for high-fidelity CFD
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Aerospace Engineering
Degree Conferral Date
2023-08
Publication Date
07/19/2023
Defense Date
06/26/2023
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
aero-propulsive coupling,body-force model,computational fluid dynamics,fan distortion response,OAI-PMH Harvest,propulsor modeling
Format
theses
(aat)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Uranga, Alejandra (
committee chair
), Bermejo-Moreno, Ivan (
committee member
), Nakano, Aiichiro (
committee member
)
Creator Email
rayevision77@gmail.com,tianboxi@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC113281399
Unique identifier
UC113281399
Identifier
etd-XieTianbo-12112.pdf (filename)
Legacy Identifier
etd-XieTianbo-12112
Document Type
Dissertation
Format
theses (aat)
Rights
Xie, Tianbo
Internet Media Type
application/pdf
Type
texts
Source
20230719-usctheses-batch-1070
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
aero-propulsive coupling
body-force model
computational fluid dynamics
fan distortion response
propulsor modeling