二維橢圓型流動(dòng)傳熱程序

PAGEPAGE17二維橢圓型流動(dòng)傳熱通用程序變量表及算例說明（本材料僅供教學(xué)參考）西安交通大學(xué)CFD&NHT/EHT研究中心陶文銓教授西安目錄………………2一、FORTRAN變量表…………3二、關(guān)于程序的主要說明………………………6二維橢圓型流動(dòng)和傳熱問題通用計(jì)算機(jī)算法方面的特點(diǎn)………………6各程序的主要功能…………7三種坐標(biāo)系統(tǒng)………………8網(wǎng)格系統(tǒng)與節(jié)點(diǎn)命名方法…………………9十一個(gè)例題的已知條件與求解內(nèi)容…………12例題1直角坐標(biāo)中二維穩(wěn)態(tài)無內(nèi)熱源的導(dǎo)熱……………………12例題2空心圓柱內(nèi)的穩(wěn)態(tài)熱傳導(dǎo)…………………12例題3正方形管道內(nèi)的充分發(fā)展對流換熱………13例題4內(nèi)壁上有直肋的環(huán)形通道內(nèi)充分發(fā)展對流換熱…………13例題5給定流場條件下溫度場的計(jì)算……………14例題6二維突擴(kuò)通道中的流動(dòng)與換熱……………14例題7方形通道內(nèi)的復(fù)雜充分流動(dòng)………………15例題8旋轉(zhuǎn)圓盤上的沖擊流動(dòng)……………………15例題9軸對稱燃燒內(nèi)的瞬間燃燒過程……………16例題10帶中心射流的通道內(nèi)的紊流換熱………16例題11有質(zhì)量源的流動(dòng)問題……………………17一、FORTRAN變量表ListofFORTRANVariablesACOFQuantitycalculatedbysubroutineDIFLOWtogivethecombinedconvectionanddiffusioneffect.AIM(I,J)Thecoefficient.AIP(I,J)Thecoefficient.AJM(I,J)Thecoefficient.AJP(I,J)Thecoefficient.AP(I,J)tThecoefficient;alsoinGAMSOR.APTTheunsteadyterm.AREALocalvariable,usuallytheareaofaC.V.face.ARHOLocalvariable,(area)().ARX(J)TheareaofthemainC.V.facenormaltothexdirection.ARXJ(J)ThepartofARX(J)thatoverlapsontheC.V.forV(I,J).ARXJP(J)ThepartofARX(J)thatoverlapsontheC.V.forV(I,J+1).BLBLCBLMBLPCoefficientsusedintheblockcorrection.CON(I,J)Theconstanttermbinthediscriminationequation;alsostandsforinGAMSOR.DENOMTemporarystorage.DIFFDiffusionconductanceD.DTThetimestep.DU(I,J)InfluencingU(I,J).DV(I,J)InfluencingV(I,J).F(I,J,NF)Various.FLTemporarystorageleadingtoFLOW.FLMTemporarystorageleadingtoFLOW.FLOWFLPMassflowratethroughaC.V.face.TemporarystorageleadingtoFLOW.FV(J)FVP(J)InterpolationfactorswhichgivethemassflowAtamaingridpoint,I,JasFV(J)(I,J)+FVP(I,J)(I,J+1)FX(I)Interpolationfactors,whichgivetheinterface.FXM(I)DensityRHOM(atthelocationofU(I,J))asFX(I)RHO(I,J)+FXM(I)RHO(I-1,J).FY(J)Interpolationfactors,whichgivetheinterface.FYM(J)DensityRHOM(atthelocationofV(I,J))asFY(J)RHO(I,J)+FYM(J)RHO(I,J-1).GAM(I,J)Thediffusioncoefficient.IIndexdenotingthepositioninx.IBEGIENDTemporaryvaluesusedinPRINT.IFSTThefirstvalueofIforwhichtheprint-outisarranged;usedinPRINT.IITemporaryindex.IPREFThevalueofIforthegridpoint,whichisusedasareferenceforpressure.ISTThefirstinternal-pointvalueofI.ISTFIST-1;usedinSOLVE.ITERAcounterforiterations.IT1IT2TemporaryvaluesusedinSOLVE.JIndexdenotingthepositioniny.JFLTemporaryindexusedinPRINT.JFSTSimilartoILST.JJTemporaryindexJLSTSimilartoILSTJPREFSimilartoIPREF.JSTThefirstinternal-pointvalueofJ.JSTFJST-1;usedinSOLVE.JT1JT2TemporaryvaluesusedinSOLVE.LASTThemaximumnumberofiterationsallowedbytheuser.LBLK(NF)When.TRUE.TheblockcorrectionforF(I,J,NF)isused.LPRINT(NF)When.TRUE.,F(I,J,NF)isprinted.LSOLVE(NF)When.TRUE.,WesolveforF(I,J,NF).LSTOPWhen.TRUE.Computationstops.L1ThevalueofIforthelastgridlocationinthexdirection.L2(L1-1).L3(L1-2).MODEIndexforthecoordinatesystem;=1for,=2for,=3for.M1ThevalueofIforthelastgridlocationintheydirection.M2(M1-1).M3(M1-2).NTemporarystorageforNF.NFIndexdenotingaparticular.NFMAXThelargestvalueofNFforwhichstorageisassigned.NGAMNFMAX+3;GAM(I,J)canbeconsideredasF(I,J,NGAM).NPNFMAX+1;P(I,J)canbeconsideredasF(I,J,NP).NRHONFMAX+2;RHO(I,J)canbeconsideredasF(I,J,NRHO).NTIMES(NF)ThenumberofrepetitionsofthesweepsinSOLVEforthevariableF(I,J,NF).P(I,J)Thepressurep.PC(I,J)Thepressurecorrection.PREFThepressureatthereferencepoint.PT(I)orPT(J)QT(I)orQT(J)TransformedcoefficientsintheTDMA.R(J)TheradiusrforamaingridpointI,J.REL1.0-RELAX(NF).RELAX(NF)RelaxationfactorforF(I,J,NF).RHO(I,J)Thedensity.RHOCONThevalueofforaconstant-densityproblem.RMN(J)ThevalueofradiusrforthelocationtowhichV(I,J)refers.SMAXThelargestabsolutevalueofthe“masssource”usedintheequation.SSUMThealgebraicsumofallthe“masssources”intheequation.SX(J)ScalefactorforthexdirectionatthemaingridlocationsY(J).SXMN(J)ScalefactorforthexdirectionatinterfacelocationsYV(J).TEMPTemporarystorage.TIMETimetforunsteadyproblems.TITLE(NF)AlphabetictitleforF(I,J,NF).U(I,J)Thex-directionvelocityu.V(I,J)They-directionvelocityv.VOLVolumeoftheC.V.X(I)Thevaluesofxatgridpoints.XCV(I)Thex-directionwidthsofmainC.V..XCVI(I)ThepartofXCV(I)thatoverlapsontheC.V.forU(I,J).XCVIP(I)ThepartofXCV(I)thatoverlapsontheC.V.forU(I+1,J).XCVS(I)Thex-directionwidthofthestaggeredC.V.forU(I,J).XDIF(I)ThedifferenceX(I)-X(I-1).XLThex-directionlengthofthecalculationdomain.XU(I)ThelocationsoftheC.V.faces;i.e.thelocationofU(I,J).Y(J)Thevaluesofyatgridpoints.YCV(J)They-directionwidthsofmainC.V.YCVR(J)TheareaforamainC.V.YCVRS(J)TheareafortheC.V.forv(I,J).YCVS(J)They-directionwidthofthestaggeredC.V.forV(I,J).YDIF(J)ThedifferenceY(J)-Y(J-1).YLThey-directionlengthofthecalculationdomain.YV(J)ThelocationsoftheC.V.faces;i.e.thelocationofV(I,J).二關(guān)于程序的主要說明二維橢圓型流動(dòng)和傳熱問題通用計(jì)算機(jī)程序算法方面的特點(diǎn)采用原始變量法，即以速度U、V及壓力P作為直接求解的變量守恒型的差分格式，離散方程系對守恒型的控制方程通過對控制容積作積分而得出的，無論網(wǎng)格疏密程度如何，均滿足在計(jì)算區(qū)域內(nèi)守恒的條件；采用區(qū)域離散化方法B，即先定控制體界面、再定節(jié)點(diǎn)位置采用交叉網(wǎng)格，速度U、V與其他變量分別存儲(chǔ)于三套網(wǎng)格系統(tǒng)中；不同的項(xiàng)在空間離散化過程中去不同的型線假設(shè)，源項(xiàng)采用局部線性化方法；擴(kuò)散——對流項(xiàng)采用乘方格式（但很容易轉(zhuǎn)化為中心差分、迎風(fēng)差分或混合格式）；街面上的擴(kuò)散系數(shù)采用調(diào)和平均法，而密度與流速則用線性插值；不穩(wěn)態(tài)問題采用全隱格式，以保證在任何時(shí)間步長下均可獲得具有物理意義的解；邊界條件采用附加源項(xiàng)法處理；耦合的流速與壓力采用SIMPLE算法來求解；迭代式的求解方法，對非線性問題，整個(gè)求解過程具有迭代性質(zhì)；對于代數(shù)方程也采用迭代法求解；采用交替方向先迭代法求解代數(shù)方程并補(bǔ)以塊修正技術(shù)以促進(jìn)收斂。各方程塊的主要功能主程序MAIN:規(guī)定整個(gè)計(jì)算過程的流程，決定是否停止計(jì)算。子程序SETUP1:設(shè)置與網(wǎng)格系統(tǒng)有關(guān)的，在計(jì)算過程中保持不變的幾何參數(shù)及設(shè)置U、V、P’、P、ρ、CON(SC)、AP(SP)的初值。子程序SETUP2：1.建立每一輪迭代中變量U、V及P’和其他變量離散方程的系數(shù)aE、aW、aN、aS、aP及b項(xiàng)。所有變量公用一套數(shù)組逐個(gè)求解，逐個(gè)確定各自目標(biāo)的系數(shù)，以節(jié)省內(nèi)存。在建立系數(shù)過程中調(diào)用DIFLOW子程序；2調(diào)用SOLVE子程序求解代數(shù)方程；3迭代次數(shù)指標(biāo)及時(shí)間步長增值；子程序DIFLOW：計(jì)算乘方格式中的系數(shù)A(|p|);子程序SOLVE:采用交替方向線迭代法并輔以塊修正技術(shù)求解代數(shù)方程；子程序GRID：設(shè)置為建立網(wǎng)格系統(tǒng)所必須的量，包括設(shè)置x,y方向求解區(qū)域的寬度XL及YL設(shè)置x,y方向上各自的節(jié)點(diǎn)數(shù)L1及M1設(shè)置x,y方向上控制體積的界面位置，XU(I),I=2,L1;YV(J),J=2,M1。如為均分網(wǎng)格可確定子程序UGRID,如為非均分網(wǎng)格需由用戶在GRID中把XU(I),YV(J)一一設(shè)置好。規(guī)定坐標(biāo)系，MODE=1，2，3分別相應(yīng)于直角、圓柱軸對稱及極坐標(biāo)當(dāng)MODE不等于1時(shí)，規(guī)定徑向起始點(diǎn)半徑R(I)子程序START:設(shè)置初值，包括（1）、對不穩(wěn)態(tài)問題給出初始條件（2）、對穩(wěn)態(tài)問題給出迭代求解的假定值（3）、已知的邊界值也可在此塊中，在設(shè)置初值或假定值時(shí)一并送入，但對邊界條件隨時(shí)間而異的不穩(wěn)態(tài)不能在此塊中賦值，因在整個(gè)計(jì)算過程中此會(huì)僅執(zhí)行一次。子程序DENSE:規(guī)定密度場，對于密度為常數(shù)的問題，此塊中可不設(shè)任何語句，但必須保留ENTRY及RETUIRN語句子程序BOUND:規(guī)定邊界條件，包括設(shè)置各變量的邊界條件對各特殊問題所需的量，如Nu,fRe等等可在此塊中計(jì)算，但也可在OUTPUT中設(shè)置子程序OUTPUT：打印輸出，包括每做一輪迭代，輸出一行信息，以觀察收斂情形；調(diào)用PRINT子程序，實(shí)現(xiàn)二維物理量場的輸出；用戶所需輸出的其他特殊變量子程序GAMSOR：設(shè)置擴(kuò)散系數(shù)及源項(xiàng)，包括(1)設(shè)置控制方程中規(guī)定的源項(xiàng)Sc,Sp(2)對第二，三類邊界條件設(shè)置與邊界相鄰的控制體中的附加源項(xiàng)，Sc,ad及sp,ad內(nèi)節(jié)點(diǎn)及邊界節(jié)點(diǎn)規(guī)定擴(kuò)散系數(shù)，擴(kuò)散采用附加源項(xiàng)法時(shí)令邊界擴(kuò)散系數(shù)為零。12、子程序UGRID:為均分網(wǎng)格設(shè)置界面位置13、子程序PRINT:輸出二維物理場(3)三種坐標(biāo)系統(tǒng)x,iyjx,iyj軸對稱圓柱坐標(biāo)，計(jì)算對θ=1弧度的區(qū)域進(jìn)行；MODE=2；Y(J)可以從任何起點(diǎn)算起，但R(J)須從對稱軸開始；R(1)維計(jì)算區(qū)的最小半徑。兩個(gè)區(qū)域XL,YL相等，但R(1)不同R(1)=0R(1)=0YLXLYLXLR(1)極坐標(biāo)系，Z方向厚度為1；MODE=3;計(jì)算區(qū)域的θ角須小于2π；Y(J)可從任何點(diǎn)起算，但R(J)比從中心點(diǎn)開始；R(1)為計(jì)算區(qū)域邊界的最小半徑。yyθrMODEX(I)Y(J)R(J)SX(J)1xy1.01.02xyr1.03θyrr(4)網(wǎng)格系統(tǒng)與節(jié)點(diǎn)命名方法V(i,j)U(i,j)V(i,j)U(i,j)P(i,j)P(L1,M1)U(2,M2)U(2,2)V(2,M1)V(2,2)U(L1,M2)V(L1,2)(a)總體網(wǎng)格系統(tǒng)V(L2,M1)U(L1,2)U(L1,2)U(I,J)U(I,J)I,JI-1,JV(I,J)I,JI,J-1(b)速度編號的方法XDIF(3)XDIF(3)XU(3)X(I)XDIF(2)U(3,J)U(2,J)L2L3I321XU(2)X(1)X(2)X(3)XU(L1)X(L1)XU(I)XDIF(I)XCVS(I)XCVS(3)XCV(I)XCVIP(I)XCVI(I)XCVS(L2)(c)x方向的主要幾何參數(shù)U(L1,J)U(L1,J)L1L1ARXJ(J)ARXJP(J)ARXJ(J)ARXJP(J)ARX(J)ycv(J-1)1y(1)yv(2)2y(2)yv(3)3y(3)M1y(M1)yv(M1)M2M3V(I,J)I，JYDIF(J)ycvs(j)(e)y方向的主要幾何參數(shù)U(2,J)XCVS(3)(2,J)U(3,J)XDIF(3)(d)速度U(3,J)的控制容積(3,J)(3,J)P(1,2)P(1,1)P(1,2)P(1,1)P(2,2)P(2,1)(f)壓力參考點(diǎn)（1，1）上壓力的確定P(1,1)=P(1,2)+P(2,1)-P（2,2）變量內(nèi)點(diǎn)下標(biāo)起始值變量ISTJSTФ,P’,P22U32V23三、十一個(gè)例題的已知條件與求解內(nèi)容例題1直角坐標(biāo)中二維穩(wěn)態(tài)無內(nèi)熱源的導(dǎo)熱已知：已知：求解區(qū)域如圖1所示，四個(gè)邊界上的溫度由下式?jīng)Q定；T=x+y+xy求：該區(qū)域中的溫度分布。yyx21圖1例題2空心圓柱內(nèi)的穩(wěn)態(tài)熱傳導(dǎo)已知：求解區(qū)域如圖2所示。左邊界為給定溫度，右邊界為對流，上邊界為絕熱，下邊界為給定熱流。整個(gè)計(jì)算區(qū)域內(nèi)都有源項(xiàng)。求：該區(qū)域中的溫度分布。計(jì)算方法：用附加源項(xiàng)法處理第2，3類邊界條件圖2圖2x211S=100-0.5TK=1.0K=0.2(1+T/100)T=100(1+y)絕熱yq=50例題3正方形管道內(nèi)的充分發(fā)展對流換熱已知：常物性流體在均勻壁溫的正方形截面管道內(nèi)作充分發(fā)展的對流換熱（如圖3）。求：截面上的速度分布、溫度分布、fRe及Nu。計(jì)算方法：充分發(fā)展對流換熱問題的處理特點(diǎn)圖3圖311TwTwyxz例題4內(nèi)壁上有直肋的環(huán)形通道內(nèi)充分發(fā)展對流換熱已知：常物性流體在帶直肋的環(huán)形通道內(nèi)作充分發(fā)展層流對流換熱，內(nèi)表面為均勻壁溫（周向），軸向則呈線性變化，外表面絕熱。R1=1,R2=2,а=15o。求：截面上速度分布、溫度分布及fRe、Nu。計(jì)算方法：耦合問題的一種處理方法。ααR1R2圖4例題5給定流場條件下溫度場的計(jì)算已知：流體流過一個(gè)直角的兩表面，速度場為u=Ax,V=-AY,a=10,Tin=500,Tw=100求計(jì)算區(qū)域中溫度的分布。圖5圖5ToutTin11Tw絕熱xy例題6二維突擴(kuò)通道中的流動(dòng)與換熱已知：由兩平行平