噴動流化床傳熱特性的CFD-DEM模擬

卜昌盛， 王昕曄， 張居兵， 樸桂林

(南京師范大學 能源與機械工程學院， 江蘇省能源系統過程轉化與減排技術工程實驗室， 南京 210042)

噴動流化床傳熱特性的CFD-DEM模擬

卜昌盛， 王昕曄， 張居兵， 樸桂林

(南京師范大學 能源與機械工程學院， 江蘇省能源系統過程轉化與減排技術工程實驗室， 南京 210042)

基于顆粒-顆粒、顆粒-流體間的傳熱機制建立了顆粒尺度下的傳熱模型，并將其與計算流體力學-離散顆粒模型(CFD-DEM)耦合，建立了CFD-DEM傳熱模型，在傳熱計算中采用真實的顆粒接觸剛度修正了顆粒-顆粒間的傳熱。采用典型噴動流化床內的顆粒傳熱實驗數據驗證了 CFD-DEM傳熱模型的準確性，并利用該模型分析了噴動流化床內的傳熱特性。結果表明：噴動流化床內顆粒的傳熱系數受其運動狀態的影響，顆粒在環隙區域外循環的傳熱系數比內循環傳熱系數大；噴動流化床內平均傳熱系數呈對稱分布，流化區域內的平均傳熱系數大于非流化區域，床體底部兩側及氣體入口處的平均傳熱系數最大，床層中央區域的平均傳熱系數較小.

噴動流化床； 傳熱特性； 傳熱模型； CFD-DEM

流化床具有很高的傳熱和傳質速率，被廣泛應用于石油化工、能源、食品及藥品加工等領域. 在這些領域的流化床工藝過程中，氣固及固固間的傳熱特性對流化床內顆粒的物理化學過程至關重要. 深入揭示流化床內氣固兩相的傳熱規律和機理對流化床的設計、運行和控制具有重要意義.

近二十年來，針對稠密氣固兩相流的計算流體力學(CFD)數值模擬發展迅速，主要包括基于歐拉-歐拉框架的雙流體模型(TFM)、基于歐拉-拉格朗日框架的離散顆粒模型(DEM)及直接模擬(DNS)[1]. 在研究稠密氣固流動方面CFD已逐漸成熟，將氣固傳熱機制與CFD耦合，研究和揭示氣固體系內的傳熱特性，已成為當前的研究熱點之一. DNS要求流體網格尺寸小于顆粒尺寸，從而直接求解流體的流動和傳熱控制方程[2]，但該方法計算量巨大，且模擬的顆粒數量極少. TFM將顆粒相假設為類似于流體的連續相. Chang等[3]將顆粒相傳熱與TFM耦合，成功建立了流化床的氣固流動和傳熱模型，但該方法難以在顆粒尺度上對氣固傳熱特性進行細致描述.

Zhou等[4-8]將傳熱模型與CFD-DEM模型耦合建立了CFD-DEM傳熱模型，研究了顆粒尺度下的氣固流動和傳熱特性. Li等[9]將顆粒-顆粒、顆粒-壁面、顆粒-流體傳熱模型與CFD-DEM模型耦合，研究了氣力輸送條件下輸送管道不同位置處的壁面溫度變化。Shimizu[10]將氣固傳熱模型的計算結果與相應的解析解進行對比驗證，并將該模型與CFD-DEM模型耦合,研究了流化床內顆粒與熱氣流間的傳熱規律。Hobbs[11]基于CFD-DEM傳熱模型，模擬了瀝青顆粒在滾筒式干燥機中的傳熱過程.Zhou等[12]建立了顆粒尺度下的氣固、固固傳熱模型，并將其與CFD-DEM模型耦合，研究了固定床及鼓泡流化床條件下的有效傳熱系數和單顆粒的傳熱特性，定量分析了流體-顆粒間傳熱、顆粒-顆粒間傳熱以及輻射傳熱在顆粒傳熱過程中的作用機制.Zhao等[13]采用CFD-DEM傳熱模型考查了流化床埋管周圍的氣固流動及傳熱特性，發現埋管周圍的傳熱系數呈橢圓形分布，即埋管左右兩側的傳熱系數大于埋管上下兩側的傳熱系數. 上述研究進一步揭示了顆粒尺度下流化床內氣固流動及傳熱的作用機制，同時為CFD-DEM傳熱模型與顆粒尺度化學反應模型的耦合提供了便利條件.

然而，上述研究中描述顆粒-顆粒間傳熱過程的模型相對簡單，如未考慮顆粒周圍氣膜及顆粒粗糙表面傳熱、未考慮顆粒-顆粒間傳熱。但已有研究[14-15]顯示，顆粒周圍的氣膜層及顆粒粗糙表面對顆粒間傳熱的影響不可忽略，尤其是在固定床條件下。此外，為選取合理的計算時間步長，在CFD-DEM模型中顆粒接觸剛度的取值遠小于真實值，從而擴大了顆粒碰撞變形量的預測值，產生較大的傳熱計算誤差. 為此，筆者將顆粒-顆粒間傳熱和顆粒-流體間傳熱機制與CFD-DEM模型耦合，修正了模型中因顆粒接觸剛度取值偏小而引起的傳熱計算誤差，建立了改進的CFD-DEM傳熱模型.

將PTV(Particle Tracking Velocimetry)與紅外熱成像法相結合，可同時獲得可視化流化床內顆粒的運動和溫度等實驗數據，已成為驗證CFD-DEM傳熱模型的有效手段[16]. 筆者以Tsuji等[17]在尖晶石二維噴動流化床實驗臺上的傳熱實驗數據為依據，驗證了所建立的CFD-DEM傳熱模型的準確性，并基于該模型詳細分析了噴動流化床內單顆粒的傳熱特性以及不同區域的傳熱特性.

1 CFD-DEM傳熱模型

在CFD-DEM模型中，氣相用Navier-Stokes方程描述，采用DEM方法跟蹤體系內的每個顆粒. 針對每個顆粒，考慮周圍顆粒與其接觸力、顆粒自身的重力以及氣體對顆粒的曳力和壓力梯度力. 氣固兩相間的作用力采用雙向耦合，相間作用力滿足牛頓第三定律. Liu等[18]已采用鼓泡流化床、埋管流化床和循環流化床的實驗數據驗證了冷態CFD-DEM模型預測氣固流動特性的準確性. 在此基礎上，筆者將冷態CFD-DEM模型與相間傳熱機制耦合(主要包括顆粒-顆粒直接導熱、顆粒-流體-顆粒導熱以及顆粒-氣體對流傳熱)，建立了CFD-DEM傳熱模型.

1.1顆粒相控制方程

DEM最早由Cundall等[19]提出，采用牛頓定律描述顆粒的運動、單顆粒的平動及轉動，方程如下：

式中：mi為顆粒i的質量，kg；Vi為顆粒i的體積，m3;Fg、Fij,n、Fij,t和Fd分別為顆粒所受的重力、法向力、切向力以及顆粒與流體間作用力，N；pg為顆粒所在位置的壓力梯度，Pa；vi和wi為顆粒i的線性速度矢量和角速度矢量；Ii為轉動慣量，kg/m2;rj為顆粒j直徑，m；K為顆粒接觸剛度，N/m；η為顆粒碰撞切向恢復系數；θij,n和θij,t為顆粒間法向和切向變形量，m；nij和tij分別為顆粒i和顆粒j間法向單位向量和切向單位向量；N為單個流體網格內的顆粒數目；u為流體的速度矢量，m/s；μ為滑移摩擦因數；φf為流體體積分數，%；β為曳力系數.

1.2氣相控制方程

氣相控制方程主要包括連續性方程、動量方程及能量方程：

(6)

(7)

(8)

式中：ρf為流體密度，kg/m3；pf為流體壓力，Pa;Vcell為單個流體網格的體積，m3；τf為剪切應力，N/m2；g為重力加速度，m/s2；Tf為流體溫度，K；cpf為流體比熱容，J/(kg·K)；λf為流體熱導率，W/(m·K)；Qf,i為顆粒-流體間的傳熱量，W；Qf,wall為流體-壁面間的傳熱量，W.

1.3傳熱模型

單顆粒的能量控制方程為：

(9)

式中：Ti為顆粒i的溫度，K；cpi為顆粒i的比熱容，J/(kg·K-1)；n為同時與顆粒i有熱交換的顆粒數目；Qij,cond為顆粒-顆粒間的導熱傳熱量，W；Qiwall,cond為顆粒-壁面間的導熱傳熱量，W；Qi,conv為顆粒-流體間的對流傳熱量，W.

1.3.1 顆粒-顆粒間導熱傳熱模型

圖1 顆粒-顆粒直接接觸導熱示意圖Fig.1 Particle-particle direct contact heat transfer

圖2 顆粒-流體-顆粒導熱示意圖Fig.2 Particle-fluid-particle heat transfer

由圖2可知，顆粒-流體-顆粒間的總傳熱熱阻為：

Rt,pfp=Rsi+Rsg+Rsi

(10)

當顆粒-顆粒直接接觸時，顆粒表面氣膜受碰撞顆粒擠壓后形成氣隙,存在于顆粒接觸面之間，如圖1中的放大區域所示，此時Rsc、Rsg和Rsgg并聯，顆粒-顆粒間的總傳熱熱阻為：

(11)

表1 顆粒-顆粒傳熱中不同傳熱熱阻的模型

傳熱模型中Rsi、Rsg和Rsgg的大小均與顆粒中心距Lij密切相關，而Lij的大小則由接觸顆粒間的法向變形量θij,n直接決定. 直接接觸顆粒-顆粒間的總傳熱熱阻Rt,pp與θij,n的關系如圖3所示，當θij,n從0增加到2.0×10-5mm時(CFD-DEM模型監測顯示，接觸顆粒間的θij,n在0～2.0×10-5mm內)，Rt,pp急劇減小，因而獲取準確的θij,n值對計算顆粒間的傳熱極為重要.

圖3 顆粒-顆粒間總傳熱熱阻與顆粒間法向變形量的關系

Fig.3 Total heat transfer resistance in particles vs. vertical deformation

DEM依據“軟球模型”計算顆粒間的法向變形量θij,n，計算方程[19]如下：

(12)

至此，顆粒-顆粒間的導熱傳熱量為：

(16)

式中：Tj為顆粒j的溫度，K.

顆粒-壁面間的傳熱采用類似的方法求解.

1.3.2 顆粒-流體間對流傳熱模型

顆粒雷諾數為：

Rei=ρfdiεi|uf-ui|/μf

(17)

式中：di為顆粒i的直徑，m；εi為顆粒i所處網格的空隙率；uf為顆粒周圍流體速度，m/s；ui為顆粒i的速度，m/s;μf為流體的動力黏度，Pa·s.

流化床內氣固流動劇烈，Rei在空間、時間上發生大幅變化. 根據顆粒雷諾數的大小，顆粒-流體間的對流傳熱量采用下式計算：

Qi,conv=hi,convAi(Tf,i-Ti)

(18)

Nui=hi,convdi/λf

(19)

(20)

式中：Nu為努塞爾數；Pr為普朗特數；Ai為顆粒i的表面積，m2；Tf,i為顆粒i所處網格內流體的溫度，K；hi,conv為顆粒i所處網格內顆粒與流體間的傳熱系數，W/(m2·K).

1.4邊界條件和初始條件

初始條件：噴動流化床體壁面初始溫度，T|wall,t=0=Two；顆粒初始溫度，T|particle,t=0=Tpo，其中Two和Tpo分別為壁面和顆粒初始溫度.

1.5顆粒傳熱系數的計算

顆粒傳熱系數是衡量氣固傳熱特性的重要參數之一. 顆粒傳熱系數包括顆粒-流體間對流傳熱系數、顆粒-顆粒間導熱傳熱系數和顆粒-壁面間導熱傳熱系數. 傳熱系數(HTC)的計算方法見表2，其中hi,wall為顆粒i與碰撞壁面的導熱傳熱系數，Twall為壁面溫度.

表2 傳熱系數的計算方程

2 模型的數值求解

采用C語言編寫三維DEM運動及傳熱程序，通過用戶自定義接口編譯進Fluent求解器，建立CFD-DEM數值傳熱模擬平臺. 每個時間步長下，在計算流體流動前，DEM程序提供顆粒的速度、溫度和位置等信息用于計算顆粒-顆粒間的傳熱量、流體網格內空隙率、顆粒與流體間相互作用力及熱量交換等. 氣固兩相流場及溫度場每隔0.025 s保存一次.

3 模型驗證及結果分析

3.1模型驗證

Tsuji等[17]采用PTV及紅外熱成像法對噴動流化床內顆粒的運動行為和溫度變化進行了實驗研究. 在實驗中，噴動流化床床體的下半部為尖晶石，上半部為玻璃。初始溫度為423 K的鋁顆粒從一定高度自由下落，堆積于床體內，溫度為292 K的空氣從床體底部中央狹縫噴嘴進入，狹縫噴嘴尺寸為21 mm×11 mm，鋁顆粒的最小臨界流化風速為1.09 m/s，實驗設定流化風速為1.5 m/s. 采用表3所示的實驗數據[17]對所建立的CFD-DEM傳熱模型進行驗證.

表3 數值計算選用的參數

根據CFD-DEM傳熱模型計算出顆粒溫度，進而繪制出顆粒相的溫度圖，并與實驗結果進行對比. 0.05～7.05 s內顆粒相溫度隨時間的變化如圖4所示，圖4(a)為實驗結果，圖4(b)為采用假設顆粒接觸剛度計算所得的模擬結果，圖4(c)為采用真實顆粒接觸剛度計算所得的模擬結果. 由圖4可以看出，在實驗中，熱顆粒加入噴動流化床后與周圍環境進行傳熱，導致貼近壁面處的顆粒溫度降低. 在模擬中，設置顆粒的初始溫度為423 K. 對比圖4(a)與圖4(b)可知，圖4(b)中顆粒溫度降低的速度明顯比實驗中快，在7.05 s尤為顯著. 這是因為模型在計算傳熱過程中選用的顆粒接觸剛度遠小于真實顆粒接觸剛度，過大地預測了顆粒-顆粒間及顆粒-壁面間的碰撞法向變形量，大幅提高了顆粒-顆粒間及顆粒-壁面間的傳熱量，導致模型預測結果與實驗結果存在較大差異. 當采用真實顆粒接觸剛度修正顆粒間的碰撞傳熱量后，模擬結果與實驗結果基本吻合，過量預測傳熱量的問題得到了有效解決.

(a) 實驗結果

(b) 假設顆粒接觸剛度的模擬結果

(c) 真實顆粒接觸剛度的模擬結果

圖4 實驗和模擬得到的不同時刻噴動流化床內顆粒的位置與溫度

Fig.4 Experimental/simulated position and temperature of particles in a spouted fluidized bed at different moments

3.2顆粒的傳熱特性

3.2.1 單顆粒的傳熱特性

選取初始位置在床層中央的單顆粒作為研究對象，考察其位置、溫度及傳熱系數隨時間的變化規律.圖5(a)給出了采用PTV測量所得不同初始位置顆粒的運動軌跡(其中數字代表不同的顆粒)，圖5(b)～圖5(d)給出了模擬獲得的0～7.05 s內每隔0.025 s顆粒的位置(即運動軌跡，其中X和Y為顆粒坐標位置)、溫度及傳熱系數隨時間的變化. 由圖5(b)可以明顯觀察到顆粒參與了環隙區域的外循環(0～3.5 s)和內循環(3.5～7.05 s)，與PTV的測量結果一致，表明CFD-DEM傳熱模型可準確預測噴動流化床內的氣固流動. 在圖5(c)中t1和t2時間段內，處于噴動流化床外循環下部的顆粒溫度基本不變，但圖5(d)顯示的顆粒瞬時傳熱系數并不為零，主要原因在于該時間段內顆粒與周圍環境的溫差較小，傳熱量較小. 與此相應，圖5(c)中部分時刻的顆粒溫度迅速下降，且持續時間較長，此段時間內顆粒由床層底部經氣流攜帶向上運動，氣固間的相對速度大，對流傳熱較強，熱量由熱顆粒向冷空氣快速傳遞. 由圖5(d)可以看出， 0～3.5 s內顆粒的瞬時傳熱系數較3.5 s后大，表明噴動流化床內的顆粒在環隙區域的外循環傳熱系數比內循環傳熱系數大，因為參與環隙區域外循環的顆粒-顆粒、顆粒-流體的相對速度比內循環大，故而傳熱增強.

(a) 實驗得到的顆粒運動軌跡

(b) 模擬得到的顆粒運動軌跡

(c) 溫度隨時間的變化

(d) 傳熱系數隨時間的變化

圖5 顆粒的運動軌跡以及顆粒溫度和傳熱系數隨時間的變化

Fig.5 Experimental/simulated trajectory of particles and variation of temperature and heat transfer coefficient with time

3.2.2 傳熱系數在噴動流化床內的分布特性

如圖6所示，根據顆粒的運動狀態可將噴動流化床內的流動區域分為流化區域和非流化區域，與非流化區域相比，流化區域內顆粒的運動更活躍，而且其溫度更高.

選取2～7 s內氣固流動及傳熱處于較穩定的狀態，進一步分析噴動流化床內的傳熱特性. 噴動流化床內顆粒相在2～7 s時的平均傳熱系數見圖7. 噴動流化床內左右兩側的平均傳熱系數呈對稱分布，床體底部兩側非流化區域及噴動流化床氣體入口處的平均傳熱系數較大，因此流化區域的平均傳熱系數顯著小于非流化區域，這也解釋了流化區域顆粒溫度高于非流化區域的原因. 盡管非流化區域內的氣流速度較低，但顆粒同時與噴動流化床側面和前后面的壁面接觸，導熱接觸面積較大，平均傳熱系數約為180 W/(m2·K)；氣體入口處的氣固相對速度大，氣固間的對流傳熱顯著，平均傳熱系數約為150 W/(m2·K). 此外，床體上還存在平均傳熱系數相對較小的3個部分區域，分別位于噴動流化床中部的兩側及上部的中間區域，這3個區域內大部分顆粒同時向下或向上運動，顆粒-顆粒間及顆粒-流體間的相對速度較小，熱導率較小，以對流傳熱為主.

圖6 3.2 s時噴動流化床內的流動結構Fig.6 Flow structure in the spouted fluidized bed at 3.2 s

圖7 2～7 s時噴動流化床內的平均傳熱系數分布

Fig.7 Average heat transfer coefficient in the spouted fluidized bed during 2-7 s

4 結 論

(1) 在傳熱模型耦合CFD-DEM模型中，采用真實顆粒接觸剛度對顆粒碰撞傳熱過程進行了修正，計算所得的顆粒溫度與實驗結果吻合較好，所建模型為研究顆粒尺度下的氣固流動及傳熱特性提供了參考。

(2) 噴動流化床內單顆粒的傳熱系數受其運動狀態的影響，顆粒在環隙區域外循環的傳熱系數比內循環傳熱系數大；對噴動流化床整體而言，流化區域內的平均傳熱系數大于非流化區域，床體底部兩側及氣體入口處的平均傳熱系數最大，床層中央區域的平均傳熱系數較小.

[1] DAN C, WACHS A. Direct numerical simulation of particulate flow with heat transfer[J].InternationalJournalofHeatandFluidFlow, 2010, 31(6): 1050-1057.

[2] YU Z S, SHAO X M, WACHS A. A fictitious domain method for particulate flows with heat transfer[J].JournalofComputationalPhysics, 2006, 217(2): 424-452.

[3] CHANG J, YANG S Q, ZHANG K. A particle-to-particle heat transfer model for dense gas-solid fluidized bed of binary mixture[J].ChemicalEngineeringResearchandDesign, 2011, 89(7): 894-903.

[4] ZHOU Z Y, YU A B, ZULLI P. A new computational method for studying heat transfer in fluid bed reactors[J].PowderTechnology, 2010, 197: 102-110.

[5] 武錦濤, 陳紀忠, 陽永榮. 移動床中顆粒接觸傳熱的數學模型[J].化工學報, 2006, 57(4): 719-725.

WU Jintao,CHEN Jizhong,YANG Yongrong. Model of contact heat transfer in granular moving bed[J].CIESCJournal, 2006, 57(4): 719-725.

[6] 劉安源, 劉石. 流化床內顆粒碰撞傳熱的理論研究[J].中國電機工程學報, 2003, 23(3): 161-165.

LIU Anyuan, LIU Shi. Theoretical study on impact heat transfer between particles in fluidized bed[J].ProceedingsoftheCSEE, 2003, 23(3): 161-165.

[7] WAHYUDI H, CHU K W, YU A B. 3D particle-scale modeling of gas-solids flow and heat transfer in fluidized beds with an immersed tube[J].InternationalJournalofHeatandMassTransfer, 2016, 97: 521-537.

[8] HOU Q F, ZHOU Z Y, YU A B. Computational study of heat transfer in a bubbling fluidized bed with a horizontal tube[J].AIChEJournal, 2012, 58(5): 1422-1434.

[9] LI J, MASON D J. A computational investigation of transient heat transfer in pneumatic transport of granular particles[J].PowderTechnology, 2000, 112: 273-282.

[10] SHIMIZU Y.Three-dimensional simulation using fixed coarse-grid thermal-fluid scheme and conduction heat transfer scheme in distinct element method[J].PowderTechnology, 2006, 165(3): 140-152.

[11] HOBBS A. Simulation of an aggregate dryer using coupled CFD and DEM methods[J].InternationalJournalofComputationalFluidDynamics, 2009, 23(2): 199-207.

[12] ZHOU Z Y, YU A B, ZULLI P. Particle scale study of heat transfer in packed and bubbling fluidized beds[J].AIChEJournal, 2009, 55(4): 868-884.

[13] ZHAO Y Z, JIANG M Q, LIU Y L, et al. Particle-scale simulation of the flow and heat transfer behaviors in fluidized bed with immersed tube[J].AIChEJournal, 2009, 55(12): 3109-3124.

[14] TSORY T, BEN-JACOB N, BROSH T, et al. Thermal DEM-CFD modeling and simulation of heat transfer through packed bed[J].PowderTechnology, 2013, 244: 52-60.

[15] BU C S, LIU D L, CHEN X P, et al. Modeling and coupling particle scale heat transfer with DEM through heat transfer mechanisms[J].NumericalHeatTransfer,PartA:Applications, 2013, 64(1): 56-71.

[16] PATIL A V, PETERS E A J F, KUIPERS J A M. Comparison of CFD-DEM heat transfer simulations with infrared/visual measurements[J].ChemicalEngineeringJournal, 2015, 277: 388-401.

[17] TSUJI T, MIYAUCHI T, OH S, et al. Simultaneous measurement of particle motion and temperature in two-dimensional fluidized bed with heat transfer[J].KONAPowderandParticleJournal, 2010, 28: 167-179.

[18] LIU D Y, BU C S, CHEN X P. Development and test of CFD-DEM model for complex geometry: a coupling algorithm for Fluent and DEM[J].Computersamp;ChemicalEngineering, 2013, 58: 260-268.

[19] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J].Geotechnique, 1979, 29(1): 47-65.

[20] RONG D G, HORIO M. DEM simulation of char combustion in a fluidized bed[C]//SecondInternationalConferenceonCFDintheMineralsandProcessIndustriesCSIRO. Melbourne, Australia: [s.n.], 1999: 65-70.

CFD-DEMSimulationonHeatTransferCharacteristicsofaSpoutedFluidizedBed

BUChangsheng,WANGXinye,ZHANGJubing,PIAOGuilin

(Engineering Laboratory for Energy System Process Conversion amp; Emission Control Technology of Jiangsu Province, School of Energy and Mechanical Engineering, Nanjing Normal University, Nanjing 210042, China)

A CFD-DEM heat transfer model was established by coupling the model set up based on the mechanism of both particle-particle and particle-fluid-particle heat transfer with the computation fluid dynamics-discrete particle model (CFD-DEM), which adopts actual contact stiffness to calculate the heat transfer. The CFD-DEM heat transfer model was verified with experimental data of a typical spouted fluidized bed, and the model was then used to further study the heat transfer characteristics of the bed. Results indicate that the heat transfer coefficient of particles is greatly influenced by their flow status, and it is higher in interior cycle than in exterior cycle. The heat transfer coefficient is symmetrically distributed in the spouted fluidized bed, which is higher in fluidized region than in un-fluidized region. The highest heat transfer coefficient appears in both sides of the bed bottom and at the gas entrance, whereas, relatively low heat transfer coefficient exists in the middle region.

spouted fluidized bed; heat transfer characteristics; heat transfer model; CFD-DEM

2017-03-23

2017-04-17

國家自然科學基金青年基金資助項目(51606104)；安徽省科技重大專項計劃資助項目(15czz02045)；江蘇省高校自然科學研究面上資助項目(16KJB470010)

卜昌盛(1987-)，男，江蘇淮安人，講師，博士，研究方向為氣固流動模擬. 電話(Tel.)：025-85481124；E-mail: csbu@njnu.edu.cn.

1674-7607(2017)11-0903-09

TK124

A

470.10