渦量流函數(shù)法模擬方腔內(nèi)粘性不可壓流動(dòng)

李江飛，李巖芳，謝冬梅，段興華，張 康，胡 凱，逯國(guó)強(qiáng)，李宏磊

（承德石油高等?？茖W(xué)校熱能工程系，河北承德067000）

渦量流函數(shù)法模擬方腔內(nèi)粘性不可壓流動(dòng)

李江飛，李巖芳，謝冬梅，段興華，張 康，胡 凱，逯國(guó)強(qiáng)，李宏磊

（承德石油高等?？茖W(xué)校熱能工程系，河北承德067000）

通過(guò)對(duì)渦量控制方程進(jìn)行無(wú)量綱化推導(dǎo)和離散，用聯(lián)合迭代方法求解二維方腔流動(dòng)這一不可壓縮黏性典型流動(dòng)的解析解.基于Matlab編程，采用渦量流函數(shù)法求解二維方腔流動(dòng)，計(jì)算采用有限體積算法，對(duì)流項(xiàng)采用QUICK格式，擴(kuò)散項(xiàng)采用二階中心差分格式，并采用延遲修正技術(shù)的離散格式對(duì)該問(wèn)題進(jìn)行數(shù)值求解，得到流動(dòng)達(dá)到穩(wěn)定狀態(tài)時(shí)各物理量的分布.

數(shù)值模擬；方腔流動(dòng)；渦量流函數(shù)；有限容積法

Peng等[1]針對(duì)二維不可壓縮黏性流體方腔流動(dòng)邊界條件：流動(dòng)速度u、v采用無(wú)滑移邊界條件，利用動(dòng)量方程推導(dǎo)了壓力p的邊界條件，本文課題組曾采用交錯(cuò)網(wǎng)格MAC算法求解二維方腔流動(dòng)[2]，本文在此基礎(chǔ)上，擬用渦量流函數(shù)法模擬方腔內(nèi)的粘性不可壓流動(dòng)，以得到流動(dòng)達(dá)到穩(wěn)定狀態(tài)時(shí)各物理量的分布.

1 問(wèn)題分析

1.1 渦量控制方程無(wú)量綱化

守恒型渦量控制方程為：

流函數(shù)控制方程為：

用方腔高度H、流體密度ρ,和拖動(dòng)速度utop作為無(wú)量綱標(biāo)尺，將控制方程無(wú)量綱化.無(wú)量綱化后的控制方程：

上式中：

其中：H是方腔長(zhǎng)度和高度；utop為頂蓋拖動(dòng)速度；u指變量?在水平x方向的流速；v為變量?在垂直y方向的流速；μ為黏度；ρ為流體密度；ω是渦量；ψ為流函數(shù).

1.2 邊界條件

邊界條件為：流動(dòng)速度采用無(wú)滑移邊界條件，壁面處法向速度恒為0，切向速度也為0.頂蓋u=1, v=0;其余u=v=0.

流函數(shù)的邊界條件為：

渦量的邊界條件為：

2 方程離散與算法分析

2.1 離散控制方程

渦量控制方程為[3-4]：

對(duì)流項(xiàng)采用QUICK格式[5]：

對(duì)流項(xiàng)渦量?jī)?nèi)邊界點(diǎn)采用一階迎風(fēng)格式：

擴(kuò)散項(xiàng)二階中心差分格式為[6]：

離散代數(shù)方程及系數(shù)表達(dá)式為：

流函數(shù)的離散采用中心差分格式：

2.2 求解步驟及程序框圖

離散方程可采用聯(lián)合迭代方法求解[7-8]：

（1）根據(jù)ψ0i,j計(jì)算u、v,確定離散方程的系數(shù)；

（2）求解渦量代數(shù)方程，得ω0i,j；

（3）由ω0i,j求解流函數(shù)方程，得ψ1i,j；

（4）根據(jù)ψ1i,j和邊界條件確定渦量邊界值ω1B；

（5）根據(jù)ψ1i,j和ω0i,j,迭代計(jì)算直至收斂；

（6）利用壓力的Poisson方程計(jì)算壓力.

具體流程如圖1所示.

渦量方程式與流函數(shù)方程是互相耦合的，渦量以源項(xiàng)形式出現(xiàn)在流函數(shù)方程中，而渦量方程中對(duì)流項(xiàng)的系數(shù)則由流函數(shù)所決定.假設(shè)流函數(shù)的一個(gè)分布，利用該分布計(jì)算u和v，從而確定離散方程的系數(shù)；接下來(lái)計(jì)算渦量和流函數(shù).在獲得了流函數(shù)的收斂解后，可用各節(jié)點(diǎn)上的流函數(shù)離散處理壓力Poisson方程，從而計(jì)算壓力[3].在獲得了流函數(shù)的收斂解后，壓力方程的源項(xiàng)即可用各節(jié)點(diǎn)上的ψ值進(jìn)行離散處理.

3 結(jié)果分析

通過(guò)比較不同雷諾數(shù)（100和1000）分別在網(wǎng)格數(shù)為20、40、80、100時(shí)的流線圖，可以看出網(wǎng)格數(shù)對(duì)結(jié)果的影響并不是很大，在精度要求不高時(shí)可以采用網(wǎng)格數(shù)較少的計(jì)算，節(jié)省時(shí)間.穩(wěn)定狀態(tài)時(shí)，Re=1000，網(wǎng)格數(shù)=20×20、40×40、80×80、100×100的渦量圖和流函數(shù)圖，計(jì)算結(jié)果如圖2所示.

圖1 計(jì)算程序框圖

穩(wěn)定狀態(tài)時(shí)，Re=10、100、1000的流線圖，計(jì)算結(jié)果如圖3所示.

穩(wěn)定狀態(tài)時(shí)，Re=1000，網(wǎng)格數(shù)=20×20，100× 100的壓力、u、v及其中心線速度分布圖，計(jì)算結(jié)果如圖4所示.

對(duì)圖2、圖3、圖4的計(jì)算結(jié)果進(jìn)行分析，可知：

（1）二維不可壓縮黏性流體方腔流動(dòng)的中心大渦并不在中心位置，方腔內(nèi)流動(dòng)也并不對(duì)稱，表明是帶奇性的方腔流動(dòng).

（2）本文處理壓力邊界條件時(shí)，分別采用動(dòng)量方程和NeumannBoundary兩種方法，計(jì)算結(jié)果比較理想.

（3）隨著雷諾數(shù)的增大，在方腔的兩個(gè)底角開(kāi)始出現(xiàn)二次渦流，說(shuō)明有限體積算法QUICK格式計(jì)算精度較高.

（4）通過(guò)分析，對(duì)流擴(kuò)散方程中當(dāng)對(duì)流項(xiàng)占主要作用的時(shí)候，很容易得到與網(wǎng)格數(shù)無(wú)關(guān)的解.

（5）本文計(jì)算分別采用延遲修正處理和不采用延遲修正處理兩種方法.延遲修正處理能保證方程對(duì)角占優(yōu)，減少求解時(shí)間.

圖2 不同網(wǎng)格數(shù)下的渦量圖和流函數(shù)圖（Re=1000）

圖3 不同雷諾數(shù)下的流線圖（自左至右：Re=10，100，1000）

圖4 壓力、u、v及其中心線速度分布圖（Re=1000,網(wǎng)格數(shù)=20×20,100×100）

[1]PENGYF,SHIAUYH,HWANGRR.Transitionina2-Dliddrivencavityflow[J].Computer&Fluids,2003,32(3):337-352.

[2]李江飛，石兆東，段興華，等.MAC算法計(jì)算二維方腔頂蓋流動(dòng)[J].宜賓學(xué)院學(xué)報(bào)，2015，15（6）：28-31.

[3]陶文銓.數(shù)值傳熱學(xué)[M].第二版.西安:西安交通大學(xué)出版社，2001：301-305.

[4]LIUCH,LEUNGDYC.Developmentofafiniteelementsolution fortheunsteadyNavier-Stokesequationsusingprojectionmethod andfractionalθ-scheme[J].Computermethodsinappliedmechan?icsandengineering.2001,190(32-33):4301-4317.

[5]AYDINM,FENNERRT.Boundaryelementanalysisofdrivencav? ityflowforlowandmoderateReynoldsnumbers[J].IntJNumer MethFluids,2001,37(1):45-64.

[6]BRUNEAUCH,SAADM.The2Dlid-drivencavityproblemrevis?ited[J].Computers&Fluids2006,35(3):326-348.

[7]MERCANH,ATALIKK.Vortexformationinlid-drivenarc-shape cavityflowsathighReynoldsnumbers[J].EuropeanJournalofMe?chanicsB/Fluids,2009,28(1):61-71.

[8]WANGJ,LIJF,CHENGWX,etal.ComparisonofFiniteDifference andFiniteVolumeMethodforNumericalSimulationoftheIncom?pressibleViscousDrivenCavityFlow[J].AdvancedMaterialsRe?search,2013.doi:10.4028/www.scientific.net/AMR.732-733.413.

【編校：許潔】

NumericalSimulationofTwo-DimensionalDrivenSquareCavityFlowBasedonStreamFunction-Vor?ticityMethod

LIJiangfei,LIYanfang,XIEDongmei,DUANXinghua,ZHANGKang,HUKai,LUGuoqiang,LIHonglei
(DepartmentofThermalEngineering,ChengdePetroleumCollege,Chengde,Hebei067000,China)

Combinediterationmethodwasusedtogaintheanalyticalsolutionoftwo-dimensionalsquareacvityflowby nondimensionalizedderivationanddiscreteofvorticityequation.BasedonMatlabprogramming,vortexstreamfunction methodsolvedtwo-dimensionalsquarecavityflow.Finitevolumemethodsolvedtheconservationequations.QUICK schemeandtwo-order-central-differenceschemewereappliedintoconvectionanddiffusiontermsseparately.Thedelay correctiontechniqueofdiscreteformatwasappliedintofinitevolumemethodforcalculation.Thedistributionofeach physicalquantitywasobtainedwhentheflowisinasteadystate.

numericalsimulation;cavityflow;vortexstreamfunction;finitevolumemethod

TB126

A

1671-5365(2015)12-0006-04

李江飛，李巖芳，謝冬梅，等.渦量流函數(shù)法模擬方腔內(nèi)粘性不可壓流動(dòng)[J].宜賓學(xué)院學(xué)報(bào),2015,15(12):6-9. LIJF,LIYF,XIEDM,etal.NumericalSimulationofTwo-DimensionalDrivenSquareCavityFlowBasedonStreamFunction-VorticityMethod[J].JournalofYibinUniversity,2015,15(12):6-9.

2015-07-06修回：2015-07-31

河北省科技支撐項(xiàng)目“采用電氣液壓調(diào)節(jié)的蒸汽透平機(jī)運(yùn)行特性和節(jié)能分析”（13211609）

李江飛（1988-），男，講師，碩士，研究方向?yàn)橛蜌鈨?chǔ)運(yùn)

時(shí)間：2015-07-3120:38

http://www.cnki.net/kcms/detail/51.1630.z.20150731.2038.001.html