低溫恒溫箱的非穩(wěn)態(tài)數(shù)值模擬及驗(yàn)證

杜子崢 謝 晶 朱進(jìn)林

（1.上海海洋大學(xué)食品學(xué)院，上海 201306；2.上海水產(chǎn)品加工與貯藏工程技術(shù)研究中心，上海 201306）

低溫恒溫箱已廣泛應(yīng)用到醫(yī)療、科研、軍事等多個(gè)領(lǐng)域，其在食品貯藏研究中使用更是廣泛，試驗(yàn)樣品往往需要在恒定溫度下進(jìn)行貯藏，低溫恒溫箱內(nèi)部不合理溫度分布將導(dǎo)致試驗(yàn)精度下降，從而影響試驗(yàn)結(jié)果，因此有必要對低溫恒溫箱內(nèi)部流場和溫度分布進(jìn)行研究。

計(jì)算流體力學(xué)（computational fluid dynamics，CFD）相對于高成本試驗(yàn)，因其低花費(fèi)、高可靠性，已經(jīng)廣泛應(yīng)用在航空航天、氣象、水文等多個(gè)領(lǐng)域，其在改進(jìn)現(xiàn)有工藝、設(shè)計(jì)和開發(fā)新產(chǎn)品上有明顯的優(yōu)勢，而且可作為試驗(yàn)數(shù)據(jù)的補(bǔ)充和參照，研究人員可以修改多個(gè)參數(shù)，獲得不同試驗(yàn)結(jié)果，這往往是單一試驗(yàn)所無法完成的。目前，計(jì)算流體力學(xué)在冷庫和冰箱流場分布模擬中應(yīng)用較多，Gonalves等［1，2］、Belleghem等［3］和Jaramillo等［4］對冷庫風(fēng)幕進(jìn)行數(shù)值模擬，在不同風(fēng)幕吹風(fēng)角度、速度下確定了冷庫風(fēng)幕最優(yōu)運(yùn)行參數(shù)。Chourasia等［5］和Delele等［6］對有貨物貯藏情況冷庫內(nèi)部的溫度和氣流組織分布進(jìn)行數(shù)值模擬，確定了最優(yōu)貨物擺放方式和堆垛尺寸以改善內(nèi)部溫度分布情況。Amara等［7］對家用冰箱進(jìn)行數(shù)值模擬并利用粒子圖像測速技術(shù)（PIV）對冰箱內(nèi)部氣流流速分布進(jìn)行試驗(yàn)檢測，模擬值和測量值具有良好一致性。António等［8］以冷藏柜為研究對象，分別采用CFD技術(shù)和ANN模型對冷藏柜內(nèi)部溫度分布進(jìn)行模擬，發(fā)現(xiàn)采用ANN預(yù)測溫度分布絕對誤差（0.8K）小于CFD預(yù)測誤差（1K），ANN模型預(yù)測溫度分布更準(zhǔn)確。Kumluta等［9］在家用冰箱的數(shù)值模擬中也采用了ANN（人工神經(jīng)網(wǎng)絡(luò)）模型，試驗(yàn)最大相對誤差僅有2.32%，驗(yàn)證了ANN模型預(yù)測冰箱溫度分布的準(zhǔn)確性。張敏等［10］對兩種不同送風(fēng)形式（上出風(fēng)型和下出風(fēng)型）低溫恒溫箱氣流組織進(jìn)行模擬和實(shí)驗(yàn)驗(yàn)證，提出上出風(fēng)型恒溫箱降溫效果優(yōu)于下出風(fēng)型恒溫箱，但下出風(fēng)型恒溫箱氣流分布更均勻。Guo-Liang Ding等［11］利用Star-CD軟件對采用自然對流形式的冰箱進(jìn)行了CFD模擬，指出貨架與冰箱前后壁的間距是影響其內(nèi)部溫度分布的關(guān)鍵因素，縮小間距有助于冰箱內(nèi)部溫度均勻分布，最后設(shè)計(jì)了一種強(qiáng)制通風(fēng)裝置改善冰箱內(nèi)部氣流組織分布。目前，有關(guān)低溫恒溫箱降溫過程非穩(wěn)態(tài)模擬文獻(xiàn)鮮見，而為提高低溫恒溫箱控溫精度，有必要對其內(nèi)部降溫過程和氣流組織分布進(jìn)行研究，提出改進(jìn)低溫恒溫箱設(shè)計(jì)的建議，從而提高其溫度分布的均勻性。

1 模型的建立

1.1 物理模型

本試驗(yàn)研究對象是一個(gè)尺寸為55cm長×64cm寬×124cm高，采用強(qiáng)制通風(fēng)形式的低溫恒溫箱，其圓形出風(fēng)口位于低溫恒溫箱壁面中上部，距離頂壁15cm，直徑為1 0cm；矩形回風(fēng)口位于出風(fēng)口同側(cè)底部，尺寸為9cm寬×6.5cm高，低溫恒溫箱隔熱材料為0.06m厚聚氨酯夾芯板，其熱導(dǎo)率為0.022W/（m·K），模型結(jié)構(gòu)見圖1。

圖1 低溫恒溫箱結(jié)構(gòu)示意圖Figure 1 Mathematical model of the cryostat

為簡化計(jì)算，本次模擬做了以下幾點(diǎn)假設(shè)：

（1）低溫恒溫箱內(nèi)氣體為不可壓縮理想氣體且符合Boussinesq假設(shè)；

（2）不考慮恒溫箱箱門處的冷風(fēng)滲透率；

（3）忽略恒溫箱內(nèi)部的輻射傳熱；

（4）不考慮恒溫箱內(nèi)部降溫過程中因相對濕度變化導(dǎo)致系統(tǒng)吸收的潛熱；

（5）壁面具有相同的熱導(dǎo)率且外界溫度為恒定測量值。

1.2 網(wǎng)格劃分及邊界條件設(shè)定

首先用GAMBIT軟件對計(jì)算域進(jìn)行非結(jié)構(gòu)網(wǎng)格劃分，非結(jié)構(gòu)網(wǎng)格對幾何模型的適應(yīng)性好，多用于對復(fù)雜區(qū)域進(jìn)行網(wǎng)格劃分。計(jì)算域網(wǎng)格尺寸采用1.5cm，共生成網(wǎng)格數(shù)量為152 928個(gè)。

低溫恒溫箱出風(fēng)口平均流速經(jīng)過實(shí)際測量為0.6m/s，經(jīng)計(jì)算雷諾數(shù)Re＝24 500為湍流，因此采用標(biāo)準(zhǔn)k－ε模型進(jìn)行求解。近壁區(qū)域需要使用壁面函數(shù)進(jìn)行處理，大多數(shù)高雷諾數(shù)情況可使用標(biāo)準(zhǔn)或非平衡的壁面函數(shù)（Re＞106），本次模擬使用增強(qiáng)壁面處理，因其對復(fù)雜流動(dòng)尤其是低雷諾數(shù)流動(dòng)很適合。

表1 k－ε方程中的系數(shù)Table 1 Coefficients in k－εmodel

1.3 計(jì)算條件

謝晶等［12］利用PISO算法在對小型冷庫開關(guān)門過程進(jìn)行非穩(wěn)態(tài)模擬，提出PISO算法適用于瞬態(tài)過程；繆晨等［13］同樣利用PISO算法對有空氣幕冷庫開關(guān)門過程溫度和流場分布進(jìn)行預(yù)測；楊磊等［14］在冷藏庫預(yù)冷降溫過程非穩(wěn)態(tài)模擬中同樣使用PISO算法，PISO算法對于瞬態(tài)問題收斂效率較高，因此本次模擬采用PISO算法，微分離散格式選項(xiàng)中，壓力空間離散因引入Boussinesq假設(shè)采用PRESTO格式，其余項(xiàng)采用二階迎風(fēng)格式以提高計(jì)算精度。

2 出風(fēng)口溫度函數(shù)

降溫過程出風(fēng)口溫度隨著時(shí)間增加逐漸降低，其過程呈現(xiàn)非穩(wěn)態(tài)變化，若在Fluent中將出風(fēng)溫度設(shè)定為一恒定值，模擬結(jié)果與測量結(jié)果之間將存在極大誤差，從而影響整個(gè)試驗(yàn)結(jié)論。因此有必要對低溫恒溫箱出風(fēng)口溫度進(jìn)行測量，繪制出風(fēng)溫度變化曲線，選擇準(zhǔn)確的函數(shù)對溫度變化曲線進(jìn)行擬合，利用UDF模塊將擬合函數(shù)編入Fluent中來定義出風(fēng)溫度。

2.1 測量方法

為驗(yàn)證模擬精度和實(shí)現(xiàn)對出風(fēng)口溫度準(zhǔn)確測量，需對恒溫箱空箱降溫過程進(jìn)行測試，試驗(yàn)共在恒溫箱內(nèi)部布置16個(gè)T型熱電偶（見圖2），利用Fluke-NetDAQ32多點(diǎn)測溫儀（美國Fluke公司）對數(shù)據(jù)進(jìn)行采集，T型熱電偶精度為±0.5℃，熱點(diǎn)偶主要布置在恒溫箱6個(gè)壁面中心，回風(fēng)口以及出風(fēng)口?？紤]到蒸發(fā)器的安放位置，測試將4個(gè)T型熱電偶以十字形式布置在恒溫箱圓形出風(fēng)口上，每個(gè)T型熱電偶間距為5cm。將恒溫箱貯藏溫度設(shè)定為－4℃，啟動(dòng)低溫恒溫箱同時(shí)開始記錄溫度數(shù)據(jù)，每隔20s對各測試點(diǎn)溫度數(shù)據(jù)進(jìn)行采集，直到監(jiān)測點(diǎn)溫度開始出現(xiàn)周期性波動(dòng)后停止測量，對出風(fēng)口4個(gè)監(jiān)測點(diǎn)溫度取平均值并繪制溫度變化曲線（見圖3）。

圖2 恒溫箱內(nèi)部測溫點(diǎn)布置示意圖Figure 2 Arrangement of temperature measuring

2.2 擬合函數(shù)

出風(fēng)溫度變化曲線有多種函數(shù)可對其進(jìn)行擬合。楊磊等［14］在對冷庫預(yù)冷降溫過程數(shù)值模擬和試驗(yàn)研究中，使用二次多項(xiàng)式對出風(fēng)溫度曲線進(jìn)行線性擬合；楊昭等［15］利用指數(shù)函數(shù)對保鮮庫實(shí)測送風(fēng)溫度進(jìn)行線性回歸。為實(shí)現(xiàn)出風(fēng)溫度變化的準(zhǔn)確預(yù)測，本試驗(yàn)建立了二維恒溫箱模型，使用線性回歸方程（三次多項(xiàng)式、五次多項(xiàng)式、指數(shù)函數(shù)）對出風(fēng)溫度測量值進(jìn)行擬合，擬合優(yōu)度分別為0.992 69，0.998 00和0.991 64，擬合優(yōu)度是指回歸直線對觀測值的擬合程度，擬合優(yōu)度越接近1，說明回歸直線對觀測值擬合程度越好。從擬合優(yōu)度指標(biāo)判斷，五次項(xiàng)擬合函數(shù)對出風(fēng)溫度擬合度最高。

圖4是將3種不同擬合函數(shù)編入Fluent中并對降溫過程進(jìn)行模擬得到的出風(fēng)口溫度線性回歸曲線。當(dāng)降溫過程達(dá)到第26分鐘（第1 600秒）時(shí)，出風(fēng)口溫度開始周期性變化，由于出風(fēng)溫度開始發(fā)生周期變化時(shí)，庫內(nèi)溫度已經(jīng)達(dá)到相對穩(wěn)定，因此波動(dòng)過程不計(jì)入擬合過程，只對降溫過程曲線進(jìn)行擬合。

圖5是基于二維CFD技術(shù)預(yù)測不同出風(fēng)溫度函數(shù)下的恒溫箱平均降溫曲線，3種出風(fēng)溫度函數(shù)溫度模擬值均低于實(shí)際測量值，這可能是由于二維模擬無法預(yù)測水平方向上的傳熱，并不能夠準(zhǔn)確反映真實(shí)的降溫過程。三次多項(xiàng)式、五次多項(xiàng)式和指數(shù)函數(shù)在平均溫度預(yù)測上最大誤差分別僅為1.12%，1.13%，1.26%，三者較接近，因而僅從這一指標(biāo)無法判斷這3個(gè)函數(shù)的優(yōu)劣。模擬結(jié)果顯示在降溫前40s，五次項(xiàng)出風(fēng)溫度函數(shù)預(yù)測的平均溫度有一個(gè)明顯升溫過程，在第20秒溫度上升最高到29.5℃，這可能導(dǎo)致錯(cuò)誤預(yù)測降溫前期恒溫箱內(nèi)部溫度分布，延長整個(gè)模擬降溫的時(shí)間。其他兩種函數(shù)均能正確反映恒溫箱降溫過程，與測量值誤差分別僅為0.1℃和0.2℃。三次項(xiàng)函數(shù)在降溫后期過度預(yù)測了恒溫箱內(nèi)部的降溫效果，在第1 600秒模擬降溫值達(dá)到－6℃，而恒溫箱內(nèi)部平均溫度測量值為－4.1℃，指數(shù)函數(shù)在整個(gè)恒溫箱降溫過程對內(nèi)部溫降預(yù)測準(zhǔn)確，且最終模擬值與測量值誤差僅為0.08℃，遠(yuǎn)遠(yuǎn)低于三次項(xiàng)函數(shù)1.95℃和五次項(xiàng)函數(shù)的1.35℃，因此本次模擬出風(fēng)口溫度選擇指數(shù)擬合函數(shù)，指數(shù)方程為：

式中：

T——出風(fēng)口的平均溫度，K；

t——降溫時(shí)間，s。

該控制方程使用UDF編入Fluent中（udf＿temperature.c），對三維恒溫箱降溫過程進(jìn)行預(yù)測，3種出風(fēng)擬合函數(shù)見表2。

圖3 恒溫箱平均出風(fēng)溫度Figure 3 Curve of cryostat average temperature

圖4 不同函數(shù)形式對出風(fēng)溫度擬合曲線Figure 4 Fit curves of outlet temperature utilizing different function

圖5 不同擬合函數(shù)下箱內(nèi)平均溫度的模擬值與測量值對比Figure 5 The comparison between simulation temperature and experiment temperature by different functions

表2 3種函數(shù)對出風(fēng)溫度的擬合方程Table 2 Fitting equation of outlet temperature with three kinds of function

3 非穩(wěn)態(tài)模擬結(jié)果與分析

3.1 溫度分布

圖6為三維恒溫箱降溫過程測量值和模擬值降溫曲線，模擬選取指數(shù)函數(shù)對出風(fēng)口溫度進(jìn)行線性回歸后，實(shí)際測量值和模擬值之間的誤差進(jìn)一步縮小，最大誤差僅為0.82%，低于二維恒溫箱模型最大誤差1.26%，說明三維CFD模擬技術(shù)與二維相比更具有可靠性，模擬的整個(gè)精度得到顯著提高，測點(diǎn)擬合溫度和測量溫度對比見圖7。

圖6 三維恒溫箱平均溫度模擬值和測量值對比Figure 6 The comparison between simulation temperature and experiment temperature

圖7 恒溫箱測點(diǎn)溫度試驗(yàn)值與模擬值的比較Figure 7 Comparison of predicted temperatures with experimental temperatures in cryostat

圖8 恒溫箱降溫1min時(shí)溫度場和渦旋場分布示意圖Figure 8 Temperature fields and vorticity distribution of cryostat after running 1minute

圖8（a）和圖8（b）顯示恒溫箱降溫1min時(shí)溫度分布剖面圖，中心主流沿射流方向溫度逐漸降低，到達(dá)壁面發(fā)生彎折并貼近壁面流動(dòng)，造成近壁面處降溫迅速，中心主流由遠(yuǎn)離出風(fēng)口一側(cè)壁面逐漸向內(nèi)部進(jìn)行對流換熱，恒溫箱內(nèi)回風(fēng)口處形成一個(gè)回風(fēng)主流區(qū)域，回風(fēng)主流兩側(cè)有兩個(gè)對稱的溫度較高區(qū)域，熱氣流在密度差作用下向恒溫箱內(nèi)部和上部區(qū)域擴(kuò)散，使得這兩個(gè)區(qū)域降溫速率較緩慢，荊有印等［16］在冷庫CFD研究中指出渦旋區(qū)域往往是低溫貯藏貨物最差區(qū)域，圖8（c）為降溫1min時(shí)Y＝0.325m的渦旋分布截面圖，在恒溫箱頂部風(fēng)機(jī)遠(yuǎn)端形成一個(gè)渦流區(qū)，在恒溫箱中心區(qū)域同樣發(fā)現(xiàn)兩個(gè)渦旋區(qū)域，而渦旋中心也是降溫過程中溫度最高區(qū)域；圖8（d）顯示在中心主流兩側(cè)也形成多個(gè)渦旋區(qū)域，上部形成的較大渦旋迫使擴(kuò)散主流發(fā)生偏移，而渦旋的產(chǎn)生往往會(huì)伴隨機(jī)械能損耗，圖8（c）和圖8（d）在近壁區(qū)域形成多個(gè)渦旋，導(dǎo)致氣流動(dòng)能減少，從而降低了恒溫箱內(nèi)部氣流流速，影響內(nèi)部降溫過程。圖9顯示隨著降溫過程繼續(xù)，恒溫箱底部溫度均勻度得到改善，但中心和頂部區(qū)域仍是整個(gè)結(jié)構(gòu)溫度分布最高區(qū)域，在恒溫箱頂端和中心左側(cè)仍發(fā)現(xiàn)存在兩個(gè)渦旋，渦旋處溫度逐漸降低，這可能是由于紊流能量級聯(lián)過程［17］，時(shí)均紊流向流場最大尺度渦旋提供其運(yùn)動(dòng)所需能量，隨過程發(fā)展向能量較小尺度渦旋傳遞，最終將機(jī)械能消耗轉(zhuǎn)化成內(nèi)能，造成大渦處溫度降低，渦旋周邊溫度升高，圖9（c）顯示恒溫箱風(fēng)口遠(yuǎn)端近壁處渦旋上移并且渦量增大迫使出風(fēng)擴(kuò)散主流向上偏移。圖10為26min時(shí)恒溫箱內(nèi)部平均溫度達(dá)到設(shè)定值－4℃數(shù)值模擬的溫度和渦旋分布剖面圖，中心和頂端仍是降溫最差區(qū)域，近壁面左側(cè)兩個(gè)渦流仍未消失。整個(gè)降溫過程，恒溫箱中心主流兩側(cè)及恒溫箱中心近壁區(qū)域始終存在渦旋，造成氣流分布不均勻，影響恒溫箱頂部右側(cè)和中心區(qū)域降溫。

圖9 恒溫箱降溫10min時(shí)溫度場和渦旋場分布示意圖Figure 9 Temperature fields and vorticity distribution of cryostat after running 10minute

圖10 恒溫箱降溫26min時(shí)溫度場和渦旋場分布示意圖Figure 10 Temperature fields and vorticity distribution of cryostat after running 26minute

3.2 速度分布

圖11為降溫過程Y＝0.325m處XZ軸速度矢量分布剖面圖，隨著降溫過程進(jìn)行，恒溫箱內(nèi)速度分布得到明顯改善，出風(fēng)主流沿程流速逐漸降低，在末端壁面上下兩處形成風(fēng)速較高區(qū)域，恒溫箱內(nèi)部平均風(fēng)速從降溫初期0.09m/s逐漸降低到0.065m/s，內(nèi)部主要貯藏區(qū)域風(fēng)速始終處在0.05～0.10m/s區(qū)間，出風(fēng)主流遠(yuǎn)端壁面上部和下部始終存在速度較低區(qū)域，這也恰好存在渦旋，流體合理的溫度和速度是貨品冷卻過程至關(guān)重要的兩個(gè)因素，恒溫箱內(nèi)部渦旋的產(chǎn)生不但使溫度分布不均勻而且導(dǎo)致渦旋處和恒溫箱內(nèi)部流速降低，影響樣品冷卻過程，因此有必要對恒溫箱進(jìn)行改進(jìn)以消除渦旋的產(chǎn)生。

圖11 恒溫箱截面Y＝0.325速度分布示意圖Figure 11 Cross-section of velocity distribution（Y＝0.325）

4 結(jié)論與展望

本研究引入Boussinesq假設(shè)，采用PISO算法并對出風(fēng)溫度進(jìn)行修正，模擬了恒溫箱內(nèi)部降溫過程，從模擬和試驗(yàn)結(jié)果，得到以下結(jié)論：

（1）結(jié)合前人研究，本試驗(yàn)對3種不同出風(fēng)溫度擬合函數(shù)進(jìn)行討論，確定指數(shù)函數(shù)更適合預(yù)測恒溫箱降溫過程，模擬值和試驗(yàn)值最大誤差僅為0.82%。

（2）整個(gè)降溫過程恒溫箱頂端和中心區(qū)域左側(cè)始終存在兩個(gè)渦旋區(qū)域，其對氣流的動(dòng)能耗散效應(yīng)導(dǎo)致內(nèi)部氣流速度減慢，影響內(nèi)部降溫過程。

本研究對恒溫箱空箱降溫過程進(jìn)行了CFD模擬，未來研究還可以從以下幾方面展開：

（1）利用CFD技術(shù)對有熱源情況下恒溫箱降溫過程進(jìn)行數(shù)值模擬，對貨物在不同呼吸強(qiáng)度、幾何形狀、擺放方式下溫度分布進(jìn)行研究，尋找最優(yōu)化方案。

（2）研究開關(guān)門過程對恒溫箱內(nèi)部降溫及溫度分布影響，對采用強(qiáng)制通風(fēng)形式的恒溫箱，選取合適的氣流速度及設(shè)備抵御開關(guān)門過程中外部熱質(zhì)侵入。

3 Belleghem M，Verhaeghe G，T＇Joen C，et al.Heat transfer through vertically downward-blowing single-jet air curtains for cold rooms［J］.Heat Transfer Engineering，2012，33（14）：1 196～1 206.

4 Jaramillo J，Pérez-Segarra C，Oliva A，et al.Analysis of the dynamic behavior of refrigerated spaces using air curtains［J］.Numerical Heat Transfer：Part A – Applications，2009，55（6）：553～573.

5 Chourasia M K，Goswami T K.Simulation of effect of stack dimensions and stacking arrangement on cool-down characteristics of potato in a cold store by computational fluid dynamics［J］.Biosystems Engineering，2007，96（4）：503～515.

6 Delele M A，Schenk A，Tijskens E，et al.Optimization of the humidification of cold stores by pressurized water atomizers based on a multiscale CFD model［J］.Journal of Food Engineering，2009，91（2）：228～239.

7 Ben Amara S，Laguerre O，Charrier-Mojtabi M C，et al.PIV measurement of the flow field in a domestic refrigerator model：Comparison with 3Dsimulations［J］.International Journal of Refrigeration，2008，31（8）：1 328～1 340.

8 António C C，Afonso C F.Air temperature fields inside refrigeration cabins：a comparison of results from CFD and ANN modelling［J］.Applied Thermal Engineering，2011，31（6）：1 244～1 251.

9 Kumluta爧D，Karadeniz Z H，Avci H，et al.Investigation of design parameters of a domestic refrigerator by artificial neural networks and numerical simulations［J］.International Journal of Re-frigeration，2012，35（6）：1 678～1 689.

10 張雷杰，張敏，楊樂，等.兩款恒溫箱的氣流組織模擬與實(shí)驗(yàn)驗(yàn)證［J］.安徽農(nóng)業(yè)科學(xué)，2009（18）：8 752～8 754，8 764.

11 Ding G L，Qiao H T，Lu Z L.Ways to improve thermal uniformity inside a refrigerator［J］.Applied Thermal Engineering，2004，24（13）：1 827～1 840.

12 謝晶，吳天.小型冷庫開關(guān)門過程溫度場的數(shù)值模擬［J］.上海水產(chǎn)大學(xué)學(xué)報(bào)，2006（3）：3 333～3 339.

13 繆晨，謝晶.冷庫空氣幕流場的非穩(wěn)態(tài)數(shù)值模擬及驗(yàn)證［J］.農(nóng)業(yè)工程學(xué)報(bào)，2013（7）：246～253.

14 楊磊，汪小旵.冷藏庫預(yù)冷降溫過程中溫度場的數(shù)值模擬與試驗(yàn)研究［J］.西北農(nóng)林科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2008（9）：219～223.

15 楊昭，徐曉麗，李喜宏.保鮮庫三維非穩(wěn)態(tài)流場模擬及實(shí)驗(yàn)［J］.天津大學(xué)學(xué)報(bào)，2007（2）：157～162.

16 荊有印，王長海，丁桂艷，等.旋轉(zhuǎn)氣流及貨物對冷庫流場影響的研究［J］.冷藏技術(shù)，2008（4）：34～37，41.

17 葉飛.非穩(wěn)態(tài)漩渦運(yùn)動(dòng)及其產(chǎn)生機(jī)理［D］.天津：天津大學(xué)，2010.