纳维斯托克斯方程（N-S）方程数值解

1UmeåUniversity2002-02-07DepartmentofPhysicsStaffanGrundbergNumericalsolutionofNavier-StokesequationusingFEMLAB231.IntroductionThepurposeofthiscomputerlaboratoryistostudysomeexamplesofincompressibleflows:flowsinsidechannelsandpipes,flowsoutsidecircularcylindersandspheresandflowsoutsideotherobjects,e.g.,aerofoils.InallthesecasestheNavier-Stokesequationwillbesolvedtogetherwiththecontinuityequationforanincompressiblefluidv-)vvv2(t∇+∇=∇⋅+∂∂ηρρp(1a)0=⋅∇v,(1b)wherevdenotesthevelocity,ρthedensity,pthepressure,andηtheviscosity.ThenumericalequationsolverFEMLABwillbeusedtosolvethesystemofequationsconsistingoftheNavier-Stokesequationandthecontinuityequation.FemlabisafiniteelementprogramwhichrunsontopofMatlab.Thespecificationofthegeometry,thepartialdifferentialequationsandtheboundaryconditionscanbedonefromtheMatlabcommandline,inaMatlabscript(M-file),orusingtheFEMLABgraphicaluserinterface(GUI).FEMLABsolvesnon-linearpartialdifferentialequationsonthegeneralformFtud=Γ⋅∇+∂∂0,xΩ∈(2a)R=0,x∈Ω∂,(2b)whereΩisaboundeddomainandΩ∂itsboundary.Whenthedependentvariableisavectoruwithncomponents,theequationandtheboundaryconditionstakethegeneralform()ijjijnjjijFxtud=∂Γ∂+∂∂∑∑==3110,i=1,…,n,x∈Ω(2c)iR=0,i=1,…,n,x∈Ω∂.(2d)Inthecaseofatwo-dimensionalstationary(00=d)incompressibleflow,thesystemofequationsisxpyuvxuuyuyxux∂∂−∂∂−∂∂−=⎟⎟⎠⎞⎜⎜⎝⎛∂∂∂∂−⎟⎠⎞⎜⎝⎛∂∂∂∂−ρρηη(3a)ypyvvxvuyvyxvx∂∂−∂∂−∂∂−=⎟⎟⎠⎞⎜⎜⎝⎛∂∂∂∂−⎟⎠⎞⎜⎝⎛∂∂∂∂−ρρηη(3b)0=∂∂+∂∂yvxu,(3c)andinthiscase⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂−=Γ00yvxvyuxuijη(4)4and⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂+∂∂∂∂+∂∂+∂∂∂∂+∂∂+∂∂−=yvxuypyvvxvuxpyuvxuuFiρρρρ.(5)InaMatlabscriptthearraysΓandFareenteredintothearraysfem.equ.gaandfem.equ.f.TheboundaryconditionR=0ontheboundaryisaDirichletboundarycondition,i.e.,thevaluesofthetwovelocitycomponentsandthepressureisspecifiedbythisboundarycondition.HowthevectorRisspecifiedisillustratedinthefollowingexample:inthecaseofaninflowboundarywherethevelocityattheinletis(u,v)=(u0,0)theboundaryconditionis⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−−==⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛00000vuuRi.(6)Atwo-dimensionalflowinachannel,i.e.,aflowbetweentwoparallelwalls,isconsideredinsection2.Aflowinsideacircularpipeisstudiedinsection3.Inthefourthandfifthsectionstheflowsoutsideacircularcylinderandoutsideasphereareinvestigated.Thedependenceoftheliftanddragonanaerofoilontheangleofattackisdeterminedinsection6.Finally,thedependenceofthedragonastreamlinedobjectonitsshapeisstudiedinsection7.2.ChannelflowConsiderastationaryincompressibleflowinachannelofwidthdformedbytwoparallelplates,seefigure1.Theplatesareassumedtobeinfinitelylonginthez-directionsothevelocityonlydependsonthexandycoordinates.Figure1.Flowbetweentwoparallelplates.Thevelocityattheinlet(x=0)isu0inthex-direction.Downstreamtheinletthevelocitywillapproachthevelocityprofileofafullydevelopedchannelflow,i.e.,vxyzd5()()2max21)(dyuyu−=.(7)Theentrylength,i.e.,thedistancefromtheinlettothepointonthex-axiswherethevelocityreachese.g.90%ofthefullydevelopedvalue,umax,dependsontheReynoldsnumberηρdu0Re=,(8)whereρisthedensityandηisthedynamicviscosity.ItisconvenienttorewritetheNavier-Stokesequationandthecontinuityequationondimensionlessformmeasuringlengthsinunitsofd,velocitiesinunitsofu0,andpressureinunitsof20uρ.Ifthedimensionlessquantitiesaredenotedbyprimes(dxx/=′,dyy/=′,0/uuu=′and)/(20uppρ=′)thedimensionlessequationsreadv-v)v′∇′+′∇′=′∇′⋅′2Re1(p(9a)0=′⋅∇′v.(9b)Attheinletthevelocityis1'=uandatthewallstheno-slipboundaryconditionisused,i.e.,0)2/1()2/1(=±=′′=±=′′yvyu.(10)WearenowreadytoimplementourmodelinFEMLABusingthegraphicaluserinterface(GUI).Thelistofinstructionsbelowdescribesstepbystephowtospecifythegeometry,theboundaryconditionsandtheReynoldsnumber.1.StartMatlabandtype:femlab.2.ChoosePhysicsModesintheModelNavigator.ChooseIncompressibleNavier-StokesandStationary.3.SetAxesEqualoffunderoptions.Setxlimitsto[-0.2,5.2]andylimitsto[-1,1]inAxes/GridSettingsunderoptions.Setx-spacingto0.5andy-spacingto0.1inAxes/GridSettingsunderoptions.4.Drawarectanglewithlength5andheight1withthelowerleftcornerat(0,-0.5).Usetherectanglebuttoninthetopofthedrawtoolbar.MakesurethesnapmodemodeisactivebydoubleclickingonSNAP.Whenthesnapmodeisactivethedrawingissnappedtothegridpoints.5.SpecifyboundaryconditionsunderBoundary.Settheinletvelocityto(u,v)=(1,0).Theinletisboundary1.Usetheno-slipconditionforthewallsofthechannelandsettheoutletpressuretozero.PressApplyandOKwhenyouhavespecifiedtheboundaryconditions.6.SpecifytheReynoldsnumberbychoosingAdd/EditConstantsunderOptions.Givethenameofthevariable(Re)andavalue,e.g.,1.PressApplyandOK.ChooseSubdomainSettingsunderSubdomain.Theequationswhichwillbesolvednowappearonthescreen.Theseequationshavethesameformasequation(1).Inordertosolvethedimensionlessequations(9)weformallysetηto1/Reandρto1.Observethatthisdoesnotimplythatthephysicaldensityis1andthatthedynamicviscosityis1/Re,itisjustawayofchoosingthecoefficientsinthedisplayedequationinordertosolvetheproblemspecifiedbytheequations(9).PressApplyandOK.7.Beforetheproblemissolvedatriangulationofthechannelmustbeperformed.ChooseInitializeMeshunderMesh.ChooseRefineMeshtwiceinordertoobtainasufficientlyfinegrid.68.Solvetheproblembypressingthebuttonwithanequalitysign.9.ThesolutioncanbepresentedgraphicallyinvariouswaysbyselectingPlotParametersinthePostmenu.Task:DeterminethedependenceoftheentrylengthontheReynoldsnumber.Thevaluesof,e.g.,thevelocitycomponentuatthenodescanbeinterpolatedinMatlabusingthecommand:[uu,pe]=postinterp(fem,’u’,xx)wherexxisa2×nimatrixwherethefirstandsecondrowcontainthexandycoordinatesoftheniinterpolationpointsanduuisavectoroflengthnicontainingthevaluesofuattheinterpolationpoints.TheFEMstructurehasfirsttobeexportedfromFEMLABtoMatlab,whichcanbedoneunderFileintheGUI.Thespecificationofthegeometry,thepartialdifferentialequations,boundaryconditionsandthemeshcanbesavedasaM-file(inMatlab)oraMAT-file(inGUI).3.PipeflowThenextmatterofconcernisthedependenceoftheentrylengthontheReynoldsnumberforastationaryincompressibleflowinacylindricalpipeofinternalradiusR.TheReynoldsnumberinthiscaseisdefinedasηρRu0Re=,(11)where0uisthevelocityattheinlet.Duetothegeometryoftheproblemitisconvenienttousecylindricalcoordinates.Weintroducedimensiolessvariablesbymeasuringlengthinunitsofd,velocitiyinunitsofu0,andpressureinunitsof20uρ.ThedimensionlessNavier-Stokesequationsincylindricalcoordinatesforastationaryflowwithoutswirl,i.e.,noflowintheazimuthaldirection,are⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂+∂∂+∂∂−=∂∂+∂∂rvrrvzvzprvvzvvzzzzrzz1Re12222(12a)⎟⎟⎠⎞⎜⎜⎝⎛−∂∂+∂∂+∂∂+∂∂−=∂∂+∂∂222221Re1rvrvrrvzvrprvvzvvrrrrrrrz(12b)wherethez-axisisalignedwiththeaxisofthepipe.Thecontinuityequationincylindricalcoordinatescanbewrittenas0=+∂∂+∂∂rvrvzvrrz.(12c)Equation(12)canberewrittenonaformmoresuitableforFEMLAB7zprrvvzvvrrvrrzvrzzrzzzz∂∂−⎟⎠⎞⎜⎝⎛∂∂+∂∂−=⎟⎠⎞⎜⎝⎛∂∂∂∂−⎟⎠⎞⎜⎝⎛∂∂∂∂−Re1Re1(13a)rprrvvzvvrrvrvrrzvrzrrrzrrr∂∂−⎟⎠⎞⎜⎝⎛∂∂+∂∂−=⎥⎦⎤⎢⎣⎡⎟⎠⎞⎜⎝⎛−∂∂∂∂−⎥⎦⎤⎢⎣⎡∂∂∂∂−2222Re1Re1(13b)rrzvrvrzvr−∂∂−∂∂−=0(13c)Theboundaryconditionsattheboundariesoftheregioninfigure2arespecifiedinthefollowingway:Attheinlet(z=0)thevelocityinthez-directionis1.No-slipboundaryconditionsareappliedatthewallofthepipe,i.e.,r=R(orr=1indimensionlessunits).Thecentralaxis,i.e.,r=0,constitutesasymmetryboundarywherethevelocityisparalleltotheboundary.Observethatthevelocityisnotzeroatthisboundary.Attheoutletaconstantpressureisprescribed.Figure2.Boundariesandboundaryconditionsforpipeflow.Thespecificationofthegeometry,boundaryconditionsandequationsisperformedintheMatlabscriptinappendixA.YoumayrunthematlabscriptdirectlyorusetheGUI.IfyouchoosetousetheGUI,youcanusethespecificationoftheequations(ΓandF)andtheboundaryconditionsofthescriptinappendixA.Task:DeterminetheentrylengthasafunctionoftheReynoldsnumber.E.g.,ifyouwanttoplotthepressuredistributioninthepipe,usethecommandpostplot(fem,’tridata’,’u3’,’axisequal’,’on’),ifyouwantanarrowplotofthevelocityfieldyoutypepostplot(fem,’arrowdata’,{’u1’,’u2’},’axisequal’,’on’),andifyouwanttoplotstreamlinesyouusethecommandpostplot(fem,'flowdata',{{'u1','cont','on'},{'u2','cont','on'}}).4.DragonacircularcylinderInthisexampletheflowaroundacircularcylinderofradiusRwillbeconsidered.ThedependenceofthedragcoefficientontheReynoldsnumberisthemainconcernofthisexample.ThedragcoefficientisgivenbyWall,no-slipOutletp=constSymmetryboundary,0=rvInlet0,1==rzvv8RuDCD2021ρ=,(14)whereDisthedragperunitlength.TheReynoldsnumberinthiscaseisdefinedasηρRu0Re=.(15)Task:UsethematlabscriptinappendixBortheGUItodeterminethedragcoefficientforReynoldsnumbersrangingfrom210−to200.InthecaseoflowReynoldsnumbersnumericalconvergenceisreachedquickerifbothsidesoftheNavier-StokesequationaremultipliedbytheReynoldsnumberv-v)v′∇′+′∇′=′∇′⋅′2ReRe(p(16)correspondingtomultiplyingthetwofirstrowsofbothΓandFbyRe.5.DragonasphereInthisexamplethesametypeofinvestigationasinthepreviousexampleisrepeated,butinsteadofacylinderasphereofradiusRisconsidered.Thedragcoefficientisgivenby22021RuDCDρ=,(17)whereDisthedragandRistheradiusofthesphere.TheReynoldsnumberinthisexampleisdefinedasηρRu0Re=.(18)Task:UsethematlabscriptofappendixCortheGUItodeterminethedragcoefficientforReynoldsnumbersrangingfrom210−to310.6.FlowaroundanaerofoilInthisexamplethematterofinvestigationwillbehowthedragandliftonanaerofoildependsontheangleofattack,α,seefigure3.Dimensiolessunitsinwhichthechord,l,oftheaerofoil,theflowvelocity,0u,farfromtheaerofoilandthedensityallareequaltoonewillbeused.TheReynoldsnumberischosentobeηρlu0Re=.(19)Thedragandliftcoefficientsoftheaerofoilaredefinedas9luDCD2021ρ=(19)andluLCL2021ρ=,(20)whereDandLarethedragandliftperunitlength.Figure3.Flowaroundanaerofoil.Task:Determinehowthedragandliftcoefficientsforanaerofoildependsontheangleofattack,α,at1000Re=.UsethematlabscriptinappendixDortheGUI.Theangleofattackoftheaerofoilisspecifiedinthelinesalpha=12;winga=rotate(wing,-pi/180*alpha);where’alpha’isgivenindegrees,inthiscase12°.TheNavier-StokesequationsarerewrittenondimensionlessformandtheReynoldsnumberisspecifiedinthelinefem.variables={'Re',1000};7.FlowaroundastreamlinedbodyInthisexamplethedependenceofthedragonastreamlinedbodyonitsshapewillbeconsidered.Atwo-dimensionalflowwillbeconsidered,assumingthebodytobeinfinitelylonginthedirectionperpendiculartotheflow.Themaximumwidthofthebody,d,iskeptconstantwhilethelength,L,ofthetailisvaried,seefigure4.Figure4.Flowaroundastreamlinedbody.u0Ldαlu010TheReynoldsnumberηρdu0Re=(21)iskeptatthevalueRe=200.Task:DeterminehowthedragvariesasafunctionofthelengthofthetailofthebodyatRe=200.UsethematlabscriptinappendixEortheGUI.Lengthsaremeasuredinunitsofd,velocitiesaremeasuredinunitsofu0,andpressuresaremeasuredinunitsof20uρ.ThebodyconsideredisacylinderwhenL=0.5andbecomeselongatedinthedownstreamsidewhenLisincreasedfromthisvalue.Sincethebodyissymmetricaroundaplaneparalleltotheflowvelocityfarawayfromthebody,onlytheflowononesideofthesymmetryplaneiscomputedinordertosavecomputertime.11AppendixA%Pipeflowclearfem;%Geometryfem.geom=rect2(0,14,0,1.0);%PDEspecificationfem.dim=3;fem.form='general';%PDEcoefficientsfem.equ.f={{...'-y.*(u1.*u1x+u2.*u1y+u3x)';...'-y.*y.*(u1.*u2x+u2.*u2y+u3y)';...'-y.*(u1x+u2y)-u2'}};fem.equ.ga={{...{'-1/Re*y.*u1x';'-1/Re*y.*u1y'};...{'-1/Re*y.*y.*u2x';'-1/Re*y.*(y.*u2y-u2)'};...{'0';'0'}}};%Boundaryconditionsfem.bnd.r={...{'0';'-u2';'0'}...%Symmetryboundaryconditionatr=0{'0';'0';'-u3'}...%Constantpressureattheoutlet{'-u1';'-u2';'0'}...%No-slipbcatthewall{'1-u1';'-u2';'0'}};%Inletvelocity1alongthepipe%Symbolicderivativesfem=femdiff(fem);%Reynoldsno.fem.variables={'Re',100};%Initializationfem.init={{1;0;0}};%Meshfem.mesh=meshinit(fem,'hmax',0.20);%Solvenonlinearproblemfem.sol=femnlin(fem,'report','on',...'sd','on',...'norm','lpnorm',...'toln',{21e-7},...'jacobian','equ');12AppendixB%Dragonacircularcylinderclearfem;%Geometryfem.geom=rect2(0,40,0,20)-circ2(10,10,1)%geomplot(fem,'edgelabel','numeric','sublabel','numeric');%Numberofvariables(u(u1),v(u2)andp(u3))fem.dim=3;%Equationsfem.form='general';fem.equ.f={{...'-(u1.*u1x+u2.*u1y+u3x)';...'-(u1.*u2x+u2.*u2y+u3y)';...'-(u1x+u2y)'}};fem.equ.ga={{...{'-1/Re*u1x';'-1/Re*u1y'};...{'-1/Re*u2x';'-1/Re*u2y'};...{'0';'0'}}};%Boundaryconditions%4typesofboundaryconditionsand8boundariesfem.bnd.r={...{'1-u1';'-u2'}...%1:inlet:u1=1{0,'-u2'}...%2:symmetrywalls:u2=0{0,0,'-u3'}...%3:outlet:u3=0(constantpressure){'-u1','-u2'}};%4:cylinder:no-slipfem.bnd.ind=[12234444];%Symbolicderivativesfem=femdiff(fem);%Reynoldsnumberfem.variables={'Re',50};%Initializationofvelocityandpressurefem.init={{1;0;0}};%Meshfem.mesh=meshinit(fem,'hmax',0.5);%Solvenon-linearproblemfem.sol=femnlin(fem,'report','on',...'sd','on',...'norm','lpnorm',...'toln',{21e-7},...'jacobian','equ');%Calculatedragfx=posteint(fem,{'0''0''0''-(-u3+2/Re*u1x1).*nx1-1/Re*(u1x2+u2x1).*nx2'});fy=posteint(fem,{'0''0''0''-1/Re*(u1x2+u2x1).*nx1-(-u3+2/Re*u2x2).*nx2'});13AppendixC%Dragonasphereclearfem;%Geometryfem.geom=rect2(0,20,0,10)-circ2(10,0,1);%geomplot(fem,'edgelabel','numeric','sublabel','numeric')%Numberofvariables(u(u1),v(u2)andp(u3))fem.dim=3;%Equationsfem.form='general';%PDEcoefficientsfem.equ.f={{...'-Re*y.*(u1.*u1x+u2.*u1y+u3x)';...'-Re*y.*y.*(u1.*u2x+u2.*u2y+u3y)';...'-y.*(u1x+u2y)-u2'}};fem.equ.ga={{...{'-y.*u1x';'-y.*u1y'};...{'-y.*y.*u2x';'-y.*(y.*u2y-u2)'};...{'0';'0'}}};%Boundaryconditions%inlet:u1=1%walls:u2=0%outlet:u3=0(constantpressure)%sphere:no-slipfem.bnd.r={...{'1-u1';'-u2'}...{0,'-u2'}...{0,0,'-u3'}...{'-u1','-u2'}};fem.bnd.ind=[1222344];%Symbolicderivativesfem=femdiff(fem);%Reynoldsnumberfem.variables={'Re',0.01};%Initializationofvelocityandpressurefem.init={{1;0;0}};%Meshfem.mesh=meshinit(fem,'hmax',0.4);%Solvenon-linearproblemfem.sol=femnlin(fem,'report','on',...'sd','on',...'norm','lpnorm',...'toln',{21e-7},...'jacobian','equ');%Calculatedragdrag=posteint(fem,{'0''0''0''-2*pi*y.*((-u3+2/Re*u1x1).*nx1+1/Re*(u1x2+u2x1).*nx2)'});14AppendixD%Flowaroundanaerofoilclearfem;%Geometryclearpp=[-0.79497098646034836-1.3017408123791103-0.82205029013539666...0.0;-0.00193423597678921140.01740812379110257...0.237911025145067790.0];rb={[14],[1;4],zeros(3,0),(1:4)'};wt={zeros(1,0),[1;1],zeros(3,0),ones(4,1)};lr={[NaNNaN],[1;0],zeros(2,0),[0;1]};wing=solid2(p,rb,wt,lr);alpha=12;winga=rotate(wing,-pi/180*alpha);fem.geom=rect2(0,10,0,5)-move(winga,4,2);%geomplot(fem,'edgelabel','numeric','sublabel','numeric')%Numberofvariables(u(u1),v(u2)andp(u3))fem.dim=3;%Equationsfem.form='general';fem.equ.f={{...'-(u1.*u1x+u2.*u1y+u3x)';...'-(u1.*u2x+u2.*u2y+u3y)';...'-(u1x+u2y)'}};fem.equ.ga={{...{'-1/Re*u1x';'-1/Re*u1y'};...{'-1/Re*u2x';'-1/Re*u2y'};...{'0';'0'}}};%Boundaryconditions%inlet:u1=1%walls:u2=0%outlet:u3=0(constantpressure)%aerofoil:no-slipfem.bnd.r={...{'1-u1';'-u2'}...{0,'-u2'}...{0,0,'-u3'}...{'-u1','-u2'}};fem.bnd.ind=[122434];%Symbolicderivativesfem=femdiff(fem);%Reynoldsnumberfem.variables={'Re',1000};%Initializationofvelocityandpressurefem.init={{1;0;0}};%Meshfem.mesh=meshinit(fem,'hmax',0.20);15%Solvenon-linearproblemfem.sol=femnlin(fem,'report','on',...'sd','on',...'norm','lpnorm',...'toln',{21e-7},...'jacobian','equ');%CalculatedragandliftD=posteint(fem,{'0''0''0''-(-u3+2/Re*u1x1).*nx1-1/Re*(u1x2+u2x1).*nx2'});L=posteint(fem,{'0''0''0''-1/Re*(u1x2+u2x1).*nx1-(-u3+2/Re*u2x2).*nx2'});16AppendixE%Dragonastreamlinedobjectclearfem;%Geometryclearscpobjsnamesltail=1.0;p=[000.50.5+ltail1;00.50.500.5];rb={[134],[1;4],[13;25;34],zeros(4,0)};wt={zeros(1,0),[1;1],[11;0.707106781186547680.70710678118654768;11],zeros(4,0)};lr={[NaNNaNNaN],[1;0],[00;11],zeros(2,0)};obstacle=solid2(p,rb,wt,lr);fem.geom=rect2(0,20,0,4)-move(obstacle,4,0);%geomplot(fem,'edgelabel','numeric','sublabel','numeric');%Numberofvariables(u(u1),v(u2)andp(u3))fem.dim=3;%Equationsfem.form='general';fem.equ.f={{...'-(u1.*u1x+u2.*u1y+u3x)';...'-(u1.*u2x+u2.*u2y+u3y)';...'-(u1x+u2y)'}};fem.equ.ga={{...{'-1/Re*u1x';'-1/Re*u1y'};...{'-1/Re*u2x';'-1/Re*u2y'};...{'0';'0'}}};%Boundaryconditions%inlet:u1=1%walls:u2=0%outlet:u3=0(constantpressure)%body:no-slipfem.bnd.r={...{'1-u1';'-u2'}...{0,'-u2'}...{0,0,'-u3'}...{'-u1','-u2'}};fem.bnd.ind=[1222344];%Symbolicderivativesfem=femdiff(fem);%Reynoldsnumberfem.variables={'Re',200};%Initializationofvelocityandpressurefem.init={{1;0;0}};%Meshfem.mesh=meshinit(fem,'hmax',0.20);%Solvenon-linearproblemfem.sol=femnlin(fem,'report','on',...'sd','on',...17'norm','lpnorm',...'toln',{21e-7},...'jacobian','equ');%CalculatedragandliftD=posteint(fem,{'0''0''0''-(-u3+2/Re*u1x1).*nx1-1/Re*(u1x2+u2x1).*nx2'});