1、AIAAJOURNALVol.42,No.9,September2004StabilizationofLinearFlowSolverforTurbomachineryAeroelasticityUsingRecursiveProjectionMethodM.S.CampobassoandMichaelB.GilesOxfordUniversity,Oxford,EnglandOX13QD,UnitedKingdomThelinearanalysisofturbomachineryaeroelasticityreliesontheassumptionofsmalllevelofunsteadi

2、nessandrequiresthesolutionofboththenonlinearsteadyandthelinearunsteadyowequations.Theobjectiveoftheanalysisistocomputeacomplexowsolutionthatrepresentstheamplitudeandphaseoftheunsteadyowperturbationforthefrequencyofunsteadinessofinterest.ThesolutionprocedureofthelinearharmonicEuler/NavierStokessolver

3、oftheHYDRAsuiteofcodesconsistsofapreconditionedxed-pointiteration,whichinsomecircumstancesbecomesnumericallyunstable.PreviousworkhadalreadyhighlightedthephysicaloriginofthesenumericalinstabilitiesanddemonstratedthecodestabilizationachievedbywrappingthecorepartofthelinearcodewithaGeneralizedMinimalRe

4、sidual(GMRES)solver.Theimplementationandtheuseofanalternativealgorithm,namely,theRecursiveProjectionMethod,issummarized.Thissolverisshowntobewellsuitedforbothstabilizingthexed-pointiterationandimprovingitsconvergencerateintheabsenceofnumericalinstabilities.Intheframeworkofthelinearanalysisofturbomac

5、hineryaeroelasticity,thismethodcanbecomputationallycompetitivewiththeGMRESapproach.I.IntroductionHEaeroelasticphenomenaofconcernintheturbomachineryindustryarebladeutterandforcedresponse,1astheycanbothleadtodramaticmechanicalfailuresifnotproperlyaccountedforinthedesignoftheengine.Thebladesofanassembl

6、ycanun-dergouttervibrationswhentheaerodynamicdampingassociatedwithcertainowregimesbecomesnegativeandisnotcounterbal-ancedbythemechanicaldamping.Insuchcircumstances,thefreevibrationofthebladestriggeredbyanytemporaryperturbationissustainedthroughtheenergyfedintothestructurebytheunsteadyaerodynamicforc

7、es.2Thehighcyclefatigue(HCF)causedbythesevibrationscanshortenthelifeofthebladesbelowthetargetlifeoftheengine.BladeforcedresponsecanalsoleadtoHCFandiscausedbytherelativemotionofadjacentframesofreference,whichtransformssteadycircumferentialvariationsoftheoweldinoneframeintoperiodictime-varyingforcesac

8、tingonthebladesintheother.Well-knownexamplesincludeforcingcausedbythewakesshedbyanupstreambladerow3andcircumferentialnonuniformitiesproducedbydistressedupstreamvanes.4Theestimationofboththemeanenergyuxbetweenuidandstructureintheuttercaseandtheunsteadyforcesactingonthebladesintheforcedresponseproblem

9、requiresknowledgeoftheunsteadyoweld.Overthepasttwodecades,anumberofap-proacheshaveemergedtocarryouttheanalysisofturbomachin-eryaeroacousticsandaeroelasticity.5Thesemethodsvaryfromun-coupledlinearizedpotentialowsolversinwhichthestructuralequationsaresolvedindependentlyoftheaerodynamics6,7tofullycoupl

10、ednonlinearthree-dimensionalunsteadyviscousmethodsinwhichthestructuralandaerodynamictime-dependentequationsaresolvedsimultaneously.8Withinthisrangetheuncoupledlinearhar-monicEulerandNavierStokes(NS)methods912haveprovedtobeasuccessfulcompromisebetweenaccuracyandcost.Thismethodviewstheaerodynamicunste

11、adinessasasmallperturbationoftheReceived21March2003;revisionreceived13February2004;acceptedc2004byM.S.Campobassoandforpublication17April2004.Copyright󰀅MichaelB.Giles.PublishedbytheAmericanInstituteofAeronauticsandAstronautics,Inc.,withpermission.Copiesofthispapermaybemadeforpersonalorinterna

12、luse,onconditionthatthecopierpaythe$10.00per-copyfeetotheCopyrightClearanceCenter,Inc.,222RosewoodDrive,Danvers,MA01923;includethecode0001-1452/04$10.00incorrespondencewiththeCCC.ResearchOfcer,ComputingLaboratory,ParksRoad.ProfessorofComputationalFluidDynamics,ComputingLaboratory,ParksRoad.LifetimeM

13、emberAIAA.1765Tspace-periodicmeansteadyow.Hencetheunsteadyoweldcanbelinearizedaboutitandbecauseoflinearitycanbedecomposedintoasumofharmonicterms,eachofwhichcanbecomputedin-dependently.Thecyclicperiodicityofboththesteadyandunsteadyowleadstoagreatreductionofcomputationalcostsbecausetheanalysiscanfocus

14、ononebladepassageratherthanthewholebladerowbymakinguseofsuitableperiodicboundaryconditions.Theas-sumptionofsmallamplitudeoftheaerodynamicunsteadinessoftenallowsonetoneglectboththecouplingandvariationsofstructuraleigenmodescausedbytheaerodynamicforces.1Thereforethein-vestigationcanbecarriedoutconside

15、ringonestructuralmodeatatime,determinedbyaniteelementprogramandusedasaninputforcalculatingtheunsteadyaerodynamicforces.Thecompleteaerody-namicanalysisconsistsoftwophases:1)calculationofthenonlinearsteadyoweldaboutwhichthelinearizationisperformedand2)solutionofthelinearharmonicequations.TheHYDRAsuite

16、ofparallelcodes1316includesbothanon-linear(hyd)andalinearharmonic(hydlin)Euler/NSsolver.Thesolutionprocedureforbothhydandhydlincanbeviewedasapre-conditionedxed-pointiteration.Usuallythelinearcodeconvergeswithoutdifculty,butproblemshavebeenencounteredinsituationsinwhichthemeanowcalculationitselffaile

17、dtoconvergetoasteadystatebutinsteadnishedinalow-levellimitcycle,oftenre-latedtosomephysicalphenomenonsuchasvortexsheddingatablunttrailingedge,unsteadyshock/boundarylayer,orshock/wakeinteraction.Inthesecircumstancesthelinearxed-pointiterationonwhichhydlinisbasedbecomesunstable,leadingtoanexponentialg

18、rowthoftheresiduals.TherelationshipbetweentheseinstabilitiesandthephysicalfeaturesoftheunderlyingbaseowisdiscussedinRef.17,whichalsosummarizesthesuccessfulimplementationofaGeneralizedMinimalResidual(GMRES)algorithm18aimedatretrievingthenumericalstabilityofthelinearcode.Forlargethree-dimensionalprobl

19、ems,however,therestartedGMRESsolvercanbecomecomputationallytooexpensiveifthenumberofKrylovvectorsperrestartedcycleneededtopreventtheresidualfromstag-natingbecomeslargerthan30.Toovercomethisproblem,anal-ternativealgorithmhasbeenimplementedinhydlin,namely,theRecursiveProjectionMethod(RPM).19Themainobj

20、ectivesofthispaperareto1)summarizethemainfeaturesofthisalgorithmanditsimplementationinhydlinand2)comparethenumericalperfor-manceoftheRPMandGMRESstabilizediterations.SectionIIpresentsanoverviewofthesteadynonlinearandunsteadylinearequations,whereastheRPMsolverisdiscussedinSec.III.Finally,thenumericalp

21、erformanceoftheRPMandGMRESalgorithmsarecomparedinSec.IV,inwhichthetwomethodsareappliedtothe1766CAMPOBASSOANDGILESutteranalysisofatwo-dimensionalturbinesectionandofacivilenginefan.II.LinearAnalysisofFlowUnsteadinessThetime-dependentEulerandReynolds-averagedNSequationsinconservativeformareapproximated

22、onunstructuredhybridgrids,usinganedge-baseddiscretization.20Consideringthecomputationaldomainconsistingofallofthepassagesofabladerowleadstoasys-temofnonlinearordinarydifferentialequations(ODE)oftheformTdUdt+R(U,Ub,X,X)=0(1)wheretisthephysicaltime,TistheJacobianofthetransformationfromprimitivetoconse

23、rvativevariables,Uisthevectorofowvariables,Risthenodalresidual,andXandXprimitivearethevectorsofnodalcoordinatesandvelocities,respectively.ThevectorUbisusedtoenforcetime-dependentdisturbancesattheinowandout-owboundariessuchaswakesshedbyanupstreambladetheresidualvectorRdependsalsoonthenodalvelocitiesX

24、row,andbecausethegridcandeformconformingtothebladevibration.Therststageofthelinearanalysisrequiresthecomputationofthemeansteadyowaboutwhichthelinearizationoftheunsteadytermswillbecarriedout.Time-averagingthegoverningequation(1)yieldsR(U¯,X¯)=0(2)whereXisthevectorofnodalcoordinatesandthebar

25、overlin-ingUandXowU¯denotestime-averagedquantities.ThemeansteadyisobtainedbysolvingEq.(2)forasinglebladepassage,becausethemeanowiscircumferentiallyperiodic.Thebound-aryconditionstowhichthesystem(2)issubjectcanbeofthreetypes:inow/outow,periodicandinviscid/viscouswall.Thefar-eldboundariesarehandl

26、edthroughuxesthatincorporatepre-scribedowinformation,andthustheybecomepartoftheresidualvectorR.Atmatchingpairsofperiodicnodes,theperiodicitycon-ditionforlinearcascadesisenforcedbysettingtheowstateontheupperboundaryequaltothatonitslowercounterpart.InthecaseofannulardomainsbecauseoftheuseofCartesianco

27、ordinates,theve-locityvectorsontheupperboundaryareobtainedbyrotatingthoseonthelowerone.Combininguxresidualsatthetwoperiodicnodesinasuitablemannertomaintainperiodicity,thisboundarycondi-tioncanalsobeincludedinthedenitionoftheuxresidualvectorR.Thetreatmentofthewallboundariesintroducessomeadditionalter

28、msinEq.(2).Thesetermsarenotreportedhereforbrevity,andtheinterestedreaderisreferredtoRefs.17and21formoredetails.Thediscreteequation(2)isthensolvedusingRungeKuttatime-marchingacceleratedbyJacobipreconditioningandmultigrid.20Thesecondstageoftheanalysisisthelinearizationoftheun-steadyowequations.Assumin

29、gthattheowunsteadinessissmall,thetime-dependentvariablescanbewrittenasthesumofameansteadypartandasmall-amplitudeperturbation:X(t)=X¯+x(t),󰀘x󰀘󰀉󰀘X¯󰀘Ub(t)=U¯b+ub(t),󰀘ub󰀘󰀉󰀘U¯b󰀘U(t)=U¯+u(t),󰀘u

30、󰀘󰀉󰀘U¯󰀘wheretheperturbationsareoverlinedwithatildeLinearizingEq.(1)aboutthemeansteadyconditions(X¯symbol.,U¯)yieldsTdudt+Lu=f1+f2(3)wherethelinearizationmatrixLandthevectorsf1andf2aregivenby󰀃󰀄L=Rf1=RR,x+Xx,f2=RubbBecauseoflinearity,thelin

31、earunsteadyoweldcanbedecom-posedintoasumofcomplexharmonicsoftheformuk(t)=󰀍elements(eiktukof),eachuofnewhichtheamplitudecanbecomputedandphaseseparately.oftheunsteadinessThecomplexkdeatfrequencyk.Analogousexpansionsholdforx(t),InsertingtheminEq.(3)andconsideringonlythexmode(t),andk=ub1(fort).s

32、implicityyieldstheharmonicequation(iT+L)u=f1+f2(4)whicharecomplexandcanbeviewedasthefrequency-domaincounterpartofEq.(3).Theright-hand-sidevectorsf1andf2givethesensitivityoftheresidualstoharmonicdeformationsofthemeshandtoincomingharmonicperturbationsrespectively.BasedonanideaofNiandSisto,22thelineare

33、quationsaresolvedwiththesamepseudo-time-marchingapproachadoptedforthesolutionofthenonlinearsteadyequations,thatis,byintroducingactitioustimederivativedu/dandtimemarchingthesolutionofthesystemoflinearODEs:dud=(iT+L)uf1f2untildu/dvanishes.Discretizingthistimederivativeleadstothelinearxed-pointiteratio

34、ndiscussedingreaterdetailinthefollow-ingsection.Intheuttercase,theobjectoftheanalysisistoassessthestabilityofaparticularstructuralmode.Thefrequencyandtheblademodeshapeareterminecalculatedwithaniteelementprogramandusedtode-f1,whichisnonzerothroughoutthecomputationaldomainbecausethegriddeformstheblade

35、,whereasconformingtotheharmonicvibrationoff2issettozero.Thephasebetweenthemotionofadjacentblades(interbladephaseangleorIBPA)isanadditionalparameteroftheanalysis.Itisgivenbyjindexjusuallycallednodaldiametercantake=any2j/integerNblades,andthevaluebe-tween0and(Nbladesrstfewones,asshown1)in,thoughRef.1.

36、theEquationcritical(4)valuescanarethenusuallybesolvedtheforasinglepassage,introducingthecomplexphaseshifteijbe-tweenthetwoperiodicboundaries.Theoutputofinterestisthenetenergyuxfromthestructuretotheworkinguidoveronecycleofvibration,denedbytheworksumintegral󰀅Tv󰀅W=publade·dSdtSinwh

37、ichTvistheperiodofvibration,pandubladearethetime-dependentbladestaticpressureandvelocityrespectively,dSistheelementalbladesurfacewithoutwardnormal,andSistheover-allbladesurface.Apositivesignindicatesstabilityasenergyistransferredfromthestructuretotheuid,whereasanegativesignindicatestheoccurrenceofut

38、ter.Intheengineeringcommunity,thelogarithmicdecrementisamorefrequentlyusedstabilitypa-rameter,whichdependsontheratiobetweentheamplitudeVoftwoconsecutivecyclesofvibration.Itisdenedas=V(t+Tanditcanbeprovedthatv)/V(t),=W/2Inforcedresponse,theobjectoftheanalysisistodeterminetheunsteadyforcesactingontheb

39、ladeasaresultofanyoftheharmoniccomponents,intowhichtheincomingtime-periodicgustcanbede-composed.TheIBPAdependsonthegeometricpropertiesoftheproblem.Inthecaseofforcingcomingfromcircumferentiallyperi-odicwakes,thebladesandthewakescanhavedifferentpitches,andhencethereisadifferenceinthetimesatwhichneighb

40、oringwakesstrikeneighboringblades.ThereforetheIBPAofthefundamentalharmonicis2Nwakes/singleNblades.Againthelinearharmonicequation(4)canbesolvedforabladeboundaryconditions.Thevectorpassageusingcomplexperiodicf1iszerothroughoutthedomainbecausethemeshisstationary,andthevectorf2isnonzeroonlyattheinletoro

41、utletboundaries,wheretheharmonicperturbationisprescribed.Theunsteadyaerodynamicforceactingonthebladecanbecalculatedinapostprocessingstepforeachstructuralmodeusingtheunsteadypressureelddeterminedwiththeharmonicanalysis.CAMPOBASSOANDGILES1767Thelinearunsteadyanalysisiscompletedbyenforcingsuitablelinea

42、rizedboundaryconditions.Theinow,outow,and(complex)periodicboundaryconditionscanallbesymbolicallyincludedintoEq.(4),whereastheadditionaltermsasaresultofthewallboundaryconditionareomittedhereandreportedinRefs.17and21.Theimplementationofthefar-eldboundaryconditionsisbasedonone-dimensionalnonreectingbou

43、ndaryconditions.23Equation(4)arethensolvedusingthesamepreconditionedpseudo-time-marchingmethodasforthenonlinearequations.III.RPMStabilizationThelinearizedharmonicowequation(4)canbeviewedasasimplelinearsystemoftheformAx=b(5)withA=iT+L,b=f1k=(2×N+f2,andx=u.Thissystemhasdimen-sioneqsis5forinviscid

44、×owsN),andwhere6forNisturbulentthenumberowofanalysesgridnodes,usingNeqsaone-equationturbulencemodel,andthefactor2accountsforrealandimaginarypartofthecomplexoweld.ThoughEq.(4)iscomplex,hydlinhasbeenwrittenusingrealarithmetic,thatis,con-sideringrealvectorsofsizekratherthancomplexvectorsofsizek/2.

45、ThischoicehasbeenmadebecauseoferrorsoftenintroducedbyhighlyoptimizedFORTRANcompilerswhendealingwithcom-plexarithmetic.Thelinearcodeforthesolutionoftheseequationscanberegardedasthexed-pointiteration:xn+1=F(xn)=(IM1A)xn+M1b(6)inwhichM1isapreconditioningmatrixresultingfromtheRungeKuttatime-marchingalgo

46、rithm,theJacobipreconditioner,andonemultigridcycle.LinearstabilityanalysisofEq.(6)showsthatanecessaryconditionforitsconvergenceisthatalloftheeigenvaluesof(IM1A)liewithintheunitcirclecenteredattheorigininthecomplexplaneorequivalentlythatalloftheeigenvaluesofM1Alieintheunitdiskcenteredat(1,0).Formosta

47、eroelasticproblemsofpracticalinterest,thisconditionisfullled,andthelinearcodeconvergeswithoutdifculty.However,anexponentialgrowthoftheresidualhasbeenencounteredinsituationsinwhichthesteadyowcalculationitselffailedtoconvergetoasteadystatebutinsteadnishedinasmall-amplitudelimitcycle,relatedtosomephysi

48、calphenomenonsuchasseparationbubbles,cornerstalls,andvortexsheddingatablunttrailingedge.Thesolutionprocedureofthenon-linearsteadyequation(2)isnottimeaccurate,butitneverthelessreectssometime-dependentphysicalpropertiesoftheoweldbecauseofthepseudo-time-marchingstrategyassociatedwiththeRungeKuttaalgori

49、thm.Physicalsmall-amplitudelimitcyclesdonotpreventthesteadysolverfromconvergingtoanacceptablelevel,andtheireffectissometimesvisibleinsmalloscillationsoftheresid-ual.Howevertheseperiodicinstabilitiesresultinasmallnumberofcomplexconjugatepairsofeigenvaluesofthelinearizationma-trixM1Alyingoutsidetheuni

50、tcircle(outliers)andthuscausingtheexponentialgrowthoftheresidualofthelinearequations.ThisproblemhadbeenpreviouslysolvedbyimplementingarestartedGMRESalgorithminhydlin.17ThedrawbackofthisapproachisthattheGMREScodecanbecomecomputationallyveryexpensivewhendealingwithlargethree-dimensionalproblems.Thisis

51、be-causeeachKrylovvectorhasthesamesizeofthelinearoweldx,andtheextramemoryrequirementwithrespecttothestandardcodegrowslinearlywiththenumberofKrylovvectors(nKr)perrestartedcycle.FurthermorenKrcannotbechosenbelowacase-dependentthresholdtopreventtheresidualfromstagnating.Thememoryre-quirementoftheGMRESc

52、odeisalreadyabouttwicethatofthestandardcodeifnKrTostabilizethelinear=30.codereducingtheadditionalmemoryre-quirement,theRPMrstintroducedbyShroffandKellertostabi-lizeunstableiterativeproceduresfornonlinearparameterdependentproblems19hasbeenimplementedinhydlin.ThisalgorithmisbasedonQoftheRprojectionkof

53、Eq.(6)ontotheorthogonalsubspacesPandassociatedrespectivelywiththesubsetofmoutliersandthatoftheremaining(km)eigenvalueslyingintheunitdisk.AteachQontoisRPMsolvediteration,onlytheprojectionofEq.(6)ontothesubspacethetypicallywiththelow-dimensionalstandardxed-pointsubspaceiteration;Pisinsteadtheprojectio

54、nsolvedwithNewtonsmethod.DenotingbyZanorthonormalbasisofP,theorthogonalprojectorsPandQofthesubspacesPandQaredened,respectively,asP=ZZTandQ=IP.Eachtimethecalculationisdiverging,thebasisZisaugmentedwiththecurrentdominanteigenmode,andtheprojectorsPandQareupdatedac-cordingly.TheprojectionsfandgofEq.(6)o

55、ntoPandQaredened,respectively,asf=PF=P(IM1A)x+M1bg=QF=Q(IM1A)x+M1bandthestabilizediterationcanbewrittenasfollows:p+1=p+(Ifp)1f(p,q)p(7a)q+1=g(p,q)(7b)wherep=Px,q=Qx,fpPFxP,Fx=IM1AItiseasilyveriedthat(Ifp)1=ZIZT(IM1A)Z1ZT=ZIH1ZT(8)where(IH)isasmallmatrixofsizem,whoseinversionrequiresminimumcomputatio

56、naleffort.Thestabilityanalysisofthisalgo-rithmshowsthatitsspectralradiusissmallerthan1,thatis,thestabilizedRPMiterationisstable.19ThebasisZisupdateddirectlyfromtheiteratesqofthemodiediteration(7b),withoutcomputingJacobians.Thisisdonebymoni-toringtherateofconvergenceoftheiteratesq.Iftheresidualstarts

57、growing,itisarguedthatsomeoftheeigenvaluesofgqunitcircle.Thegenericcaseisthateitheran=isolatedQFxQlieoutsidetherealeigenvaluem+1oracomplexconjugatepair(m+1,m+2)causetheinstability.OnehastodecidewhichisthecaseanddeterminetheoneortwovectorstoappendtoZtomakeitspanthelargerinvari-antsubspaceofFxassociat

58、edwiththeaugmentedsetofeigenvalues1,2the,.vectors,m+1orq1,2,.,m+1,m+2matrixgq.Itcanbeshown19that=+1qareapoweriterationwiththeqappliedtothestartingvector󰀍q0.Asymptoticallythesevectorswilltendtobeinthedominanteigenspaceofgq,providedthat󰀍q0hasanonzerocomponentinthisdirection.Iftheiterat

59、ionstartsdiverging,thetwodifferencevectors󰀍qmostrecentpoweriterates)areusedtocomputethe,󰀍Gramq1(i.e.,Schmidtthefactorization24Dq,q1=DT(9)withTR2×2uppertriangularandDRk×2orthogonal.IfT1,1column󰀊T2,2,theofDdominanteigenmodeofgqisreal,andonlytherstisincludedinZ.Otherwisetheinstabilityiscausedbyacomplexconjugatepair,andbothco