基于格子波爾茲曼方法的回?zé)崞鲾?shù)值模擬

夏宇棟 陳 曦 馬詩(shī)旻 張 華

(上海理工大學(xué)能源與動(dòng)力學(xué)院 上海 200093)

基于格子波爾茲曼方法的回?zé)崞鲾?shù)值模擬

夏宇棟 陳 曦 馬詩(shī)旻 張 華

(上海理工大學(xué)能源與動(dòng)力學(xué)院 上海 200093)

利用格子玻爾茲曼方法，直接對(duì)蝕刻薄片和層疊絲網(wǎng)回?zé)崞鞯奈⒂^結(jié)構(gòu)流場(chǎng)進(jìn)行了模擬。得到了兩種回?zé)崞魈盍系奈⒂^流場(chǎng)和兩端的壓差。模擬結(jié)果顯示，當(dāng)回?zé)崞鞯闹睆?、水力直徑和填充率相近情況下，不同流速下蝕刻薄片卷裹式回?zé)崞鞯姆€(wěn)態(tài)阻力系數(shù)均比層疊絲網(wǎng)回?zé)崞餍　７€(wěn)態(tài)阻力系數(shù)的模擬變化趨勢(shì)與實(shí)驗(yàn)一致。

格子波爾茲曼方法 回?zé)崞?流阻系數(shù) 數(shù)值模擬

1 引言

回?zé)崞魇堑蜏刂评錂C(jī)的關(guān)鍵部件，也是影響制冷機(jī)性能的最重要因素。目前，對(duì)回?zé)崞鞯难芯恐饕袑?shí)驗(yàn)研究、傳統(tǒng)計(jì)算方法研究和模擬計(jì)算方法研究。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，對(duì)回?zé)崞鞯哪M計(jì)算已經(jīng)成為了回?zé)崞餮芯款I(lǐng)域的一個(gè)新的熱點(diǎn)。張曉青，郭方中等［1］利用基于網(wǎng)絡(luò)模型開(kāi)發(fā)了仿真軟件對(duì)熱聲驅(qū)動(dòng)的制冷機(jī)用回?zé)崞鞯拈L(zhǎng)度、水力半徑以及回?zé)崞髟隈v波聲場(chǎng)中的位置等參數(shù)對(duì)其性能的影響進(jìn)行了模擬分析，并對(duì)其進(jìn)行了優(yōu)化，為更深地理解回?zé)崞鞯臋C(jī)理和實(shí)質(zhì)提供了數(shù)值基礎(chǔ)。章忠敏［2］根據(jù)豪森的數(shù)學(xué)模型，利用有限差分法求解對(duì)稱情況下的理想回?zé)崞鱾鳠嵛⒎址匠?，并編制了?jì)算穩(wěn)流回?zé)崞餍实耐ㄓ贸绦?，分析了理想穩(wěn)流回?zé)崞鞯膫鳠崽攸c(diǎn)和近似的數(shù)學(xué)描述。M.D.Atrey等［3］運(yùn)用有限差分法重點(diǎn)分析了回?zé)崞鲗?duì)整機(jī)性能的影響，分析給出了工質(zhì)質(zhì)量流率、平均壓力、流速、傳熱周期等對(duì)回?zé)崞餍实挠绊?，得出回?zé)崞餍适墙z網(wǎng)孔隙率的函數(shù)。

近年來(lái)，作為一種嶄新的數(shù)值模擬手段格子波爾茲曼法也逐漸被應(yīng)用到對(duì)多孔介質(zhì)的數(shù)值模擬當(dāng)中。格子波爾茲曼方法(Lattice Boltzmann Method(LBM))是一種基于粒子分布函數(shù)演化的數(shù)值方法。早期格子波爾茲曼方法是從格子氣元胞自動(dòng)機(jī)方法發(fā)展而來(lái)［4］。后來(lái)研究表明，通過(guò)離散時(shí)間、物理空間和速度空間，該方法可以直接從連續(xù)波爾茲曼方程嚴(yán)格導(dǎo)出。與傳統(tǒng)CFD方法比較，其控制方程的對(duì)流項(xiàng)是線性的，流體的壓力也無(wú)需求解Possion方程得到，同時(shí)該方法具有二階精度和良好的數(shù)值穩(wěn)定性。

本文就將運(yùn)用格子波爾茲曼方法對(duì)層疊絲網(wǎng)和蝕刻薄片卷裹式回?zé)崞魈盍系奈⒂^機(jī)構(gòu)進(jìn)行數(shù)值模擬，結(jié)合宏觀尺度上的回?zé)崞鲾?shù)值計(jì)算模型，對(duì)回?zé)崞餍阅茏龀鲱A(yù)測(cè)。

2 回?zé)崞魑锢砟Ｐ徒?/h2>

格子玻爾茲曼方法需將計(jì)算空間劃分為離散的流體微團(tuán)。因此首先得建立兩種回?zé)崞鞯奈锢砟Ｐ汀?/p>

2.1 層疊絲網(wǎng)回?zé)崞魑锢砟Ｐ徒?/h3>

層疊絲網(wǎng)是回?zé)崞魈盍现饕遣捎酶吣繑?shù)的不銹鋼或磷青銅層疊而成，其結(jié)構(gòu)形式如圖1所示。

圖1 層疊絲網(wǎng)回?zé)崞魈盍系慕Y(jié)構(gòu)Fig.1 Structure of mesh screen regenerator

其結(jié)構(gòu)特點(diǎn)是徑向?yàn)橐?guī)整的編織網(wǎng)片，軸向?yàn)榫W(wǎng)片的隨機(jī)旋轉(zhuǎn)疊加。本次模擬選用的層疊絲網(wǎng)的詳細(xì)幾何參數(shù)見(jiàn)表1所示。

表1 層疊絲網(wǎng)回?zé)崞鲙缀螀?shù)Table 1 Structure parameters of mesh screen regenerator

根據(jù)層疊絲網(wǎng)的幾何參數(shù)，利用Wolfram Mathematic 8.0軟件，通過(guò)編寫層疊絲網(wǎng)的幾何方程得到層疊絲網(wǎng)的物理計(jì)算模型，如圖2所示，所生成的層疊絲網(wǎng)為7層，每層絲網(wǎng)厚度Lr1=8.5×10-5m。計(jì)算空間為0.6 mm×0.28 mm×0.28 mm。并取格子尺寸δx=5×10-6m。則總共的網(wǎng)絡(luò)計(jì)算節(jié)點(diǎn)約為588萬(wàn)個(gè)。

圖2 三維層疊絲網(wǎng)回?zé)崞魈盍辖Y(jié)構(gòu)Fig.2 Three-dimensional model of meshscreen regenerator

2.2 蝕刻薄片卷裹式回?zé)崞魑锢砟Ｐ偷慕?/h3>

蝕刻薄片卷裹式回?zé)崞魇怯杀菊n題組自主研發(fā)的新型回?zé)崞?，它是在金屬薄片表面采用化學(xué)腐蝕或光刻的方式形成一定的圖案，通過(guò)卷裹形成流道。其幾何尺寸如圖3所示。

圖3 蝕刻薄片的幾何參數(shù)(mm)Fig.3 Geometric parameter of etched foil

同樣利用Wolfram Mathematic 8.0軟件，生成最終的物理計(jì)算模型如圖4所示。

該結(jié)構(gòu)的計(jì)算空間為0.275 mm×1.4 mm×2.6 mm，取網(wǎng)格尺寸δx=5×10-6m，則計(jì)算的網(wǎng)絡(luò)節(jié)點(diǎn)約為800.8萬(wàn)個(gè)。

圖4 層疊絲網(wǎng)回?zé)崞魈盍辖Y(jié)構(gòu)Fig.4 Three-dimensional model of etched foil regenerator

3 D3Q19不可壓等溫格子模型用于填料計(jì)算

運(yùn)用開(kāi)源的CFD框架Palabos可以對(duì)于這兩種回?zé)崞鹘Y(jié)構(gòu)進(jìn)行模擬。Palabos是由日內(nèi)瓦大學(xué)Jonas Latt博士、Orestis博士等成員聯(lián)合開(kāi)發(fā)的開(kāi)源的通用CFD計(jì)算庫(kù)，它基于格子玻爾茲曼法。2011年9月，由FlowKit公司開(kāi)始運(yùn)營(yíng)，為用戶提供解決方案。它基于C++的標(biāo)準(zhǔn)模板庫(kù)(STL)，提供了完整的程序庫(kù)(Library)框架，用戶可以在其框架的基礎(chǔ)上開(kāi)發(fā)基于格子玻爾茲曼的計(jì)算程序。

3.1 回?zé)崞饔?jì)算模型設(shè)定

本次計(jì)算的回?zé)崞鞯妮S向長(zhǎng)度僅為0.6 mm和2.6 mm，其兩端溫差較小，因此忽略溫度變化，即模擬的結(jié)構(gòu)為一個(gè)等溫的模型。所以計(jì)算選用等溫D3Q19-BGK模型，其恢復(fù)的方程與標(biāo)準(zhǔn)的不可壓Navier-Stokes方程組一致，適用于低馬赫數(shù)的流動(dòng)。

對(duì)計(jì)算物理模型，定義其兩端的速度邊界。此處采用了“濕節(jié)點(diǎn)”法中的 Regularized邊界條件［5］，這種邊界條件處理方法的特征是，邊界節(jié)點(diǎn)中的所有局部分布函數(shù)都會(huì)被替換，不管它們是否已知。

3.2 運(yùn)行結(jié)果處理及分析

對(duì)層疊絲網(wǎng)結(jié)構(gòu)回?zé)崞鳎瑸榱四苡行У嘏c實(shí)驗(yàn)進(jìn)行對(duì)比，根據(jù)實(shí)際的回?zé)崞鞴ぷ鲄?shù)，模擬選擇了入口速度為2、5、8、11 m/s 4 個(gè)工況點(diǎn)。

3.2.1 流場(chǎng)分布

圖5所示為入口速度為5 m/s時(shí)，兩種填料結(jié)構(gòu)的速度場(chǎng)向量圖。從圖中可以清晰的看出工質(zhì)在進(jìn)入層疊絲網(wǎng)結(jié)構(gòu)前速度較小，當(dāng)工質(zhì)進(jìn)入填料間的流道，由于流道變窄，工質(zhì)加速，最大速度可達(dá)到入口速度的7—10倍大小。

圖6為兩種結(jié)構(gòu)回?zé)崞魈盍狭鞯纼?nèi)的速度分布向量圖，兩種回?zé)崞魈盍辖Y(jié)構(gòu)的差異使其流場(chǎng)有很大的不同，進(jìn)而影響回?zé)崞鞯幕責(zé)嵝阅堋?/p>

圖5 入口速度為5 m/s的三維速度場(chǎng)分布Fig.5 Velocity field at inlet velocity of 5 m/s

圖6 入口速度為5 m/s的流道內(nèi)速度場(chǎng)分布Fig.6 Velocity field of flow pass at inlet velocity of 5 m/s

從圖中可以分析看出，對(duì)于層疊絲網(wǎng)結(jié)構(gòu)的回?zé)崞鳎z網(wǎng)疊加的隨機(jī)性帶來(lái)了流動(dòng)的隨機(jī)性，使得其性能很難預(yù)測(cè)。即使是相同材料、相同孔隙率、相同大小、相同水力直徑的回?zé)崞?，由于絲網(wǎng)疊加形式的不可控，也會(huì)造成性能的不同。而對(duì)于蝕刻薄片卷裹式回?zé)崞鱽?lái)說(shuō)，其速度分布更為規(guī)則和均勻，沒(méi)有出現(xiàn)圖6(a)中部分填料區(qū)域沒(méi)有發(fā)生換熱的情況。工質(zhì)在平行于流動(dòng)方向的流道只受到邊界粘應(yīng)力的影響，其阻力較小。而斜流道的出現(xiàn)，適時(shí)地打破熱邊界層對(duì)換熱的影響，同時(shí)對(duì)流阻貢獻(xiàn)不大，可以看到，工質(zhì)流過(guò)斜流道后到達(dá)直流道，速度迅速達(dá)到恢復(fù)原來(lái)的大小。這種流場(chǎng)換熱均勻，流動(dòng)阻力小，對(duì)提高回?zé)崞鞯膿Q熱效率，降低流阻教層疊絲網(wǎng)回?zé)崞饔休^大的優(yōu)勢(shì)。

3.2.2 穩(wěn)態(tài)壓降

根據(jù)不同流動(dòng)速度下，模擬的層疊絲網(wǎng)以及蝕刻薄片式回?zé)崞鹘Y(jié)構(gòu)密度云圖，可以得到回?zé)崞鹘Y(jié)構(gòu)進(jìn)出口的格子密度差，參考樣機(jī)的制冷機(jī)的工作壓力pa=1.75 MPa，取平均溫度為Ta=190 K。查氦氣物性后。由式(1)計(jì)算得到穩(wěn)態(tài)的流阻系數(shù)。兩種回?zé)崞鹘Y(jié)構(gòu)的穩(wěn)態(tài)流阻如圖7所示。

圖7 層疊絲網(wǎng)回?zé)崞骱臀g刻薄片卷裹式的流阻系數(shù)的對(duì)比Fig.7 Resistance coefficient vs Reynolds number of mesh screen and etched foil regenerator

式中:ρ為工質(zhì)密度，(kg/m3)，qv為回?zé)崞鲀?nèi)工質(zhì)體積流量，m3/s;dh為回?zé)崞魉χ睆?，m;μ為工質(zhì)動(dòng)力粘度，Pa·S;d為回?zé)崞髦睆?，m;φ為孔隙率;L為回?zé)崞鏖L(zhǎng)，m。

由圖7分析可得:隨著工質(zhì)雷諾數(shù)Re的增加，兩種回?zé)崞鞯墓べ|(zhì)穩(wěn)態(tài)流動(dòng)阻力系數(shù)f均減少，層疊絲網(wǎng)回?zé)崞鞯淖枇ο禂?shù)比蝕刻薄片卷裹式回?zé)崞鞔?，說(shuō)明蝕刻薄片卷裹式回?zé)崞饔辛鲃?dòng)阻力小的特點(diǎn)。

4 回?zé)崞鞣€(wěn)態(tài)流阻測(cè)試及分析

4.1 穩(wěn)態(tài)流動(dòng)的測(cè)試試驗(yàn)系統(tǒng)

如圖8所示為穩(wěn)態(tài)測(cè)試系統(tǒng)中的實(shí)驗(yàn)與測(cè)試系統(tǒng)圖。測(cè)試時(shí)壓力為0.3 MPa，溫度為290 K。測(cè)試對(duì)象為280目的層疊絲網(wǎng)回?zé)崞骱臀g刻薄片卷裹式回?zé)崞鳎べ|(zhì)為氦氣和氮?dú)狻?/p>

圖8 穩(wěn)態(tài)流動(dòng)系統(tǒng)測(cè)試圖1.進(jìn)口混合腔;2.出口混合腔;3.直流電源轉(zhuǎn)換器;4.差壓變送器;5.壓力傳感器;6.Agilent 34901A數(shù)據(jù)采集卡;7.T型熱電偶;8.Agilent 34970A數(shù)據(jù)采集儀。Fig.8 Test bench of stable-state flow system

4.2 LBM模擬數(shù)據(jù)與實(shí)驗(yàn)結(jié)果比較

將運(yùn)用LBM的模擬數(shù)據(jù)與實(shí)驗(yàn)所測(cè)得的數(shù)據(jù)進(jìn)行比較，其比較結(jié)果如圖9所示。

圖9 回?zé)崞髂M值和實(shí)驗(yàn)值的對(duì)比Fig.9 Value of simulation vs value of experiment

從圖中可以看出，模擬值在趨勢(shì)上是符合實(shí)驗(yàn)值的，但出現(xiàn)了整體的偏移。模擬值較實(shí)驗(yàn)值偏小50%左右。

上述模擬結(jié)果與實(shí)驗(yàn)值的對(duì)比可以看出，模擬結(jié)果在趨勢(shì)上符合實(shí)際，但在壓降的精確計(jì)算上產(chǎn)生了偏差。本文認(rèn)為，可能是由以下因素造成的:

(1)流體進(jìn)入填料結(jié)構(gòu)時(shí)，流場(chǎng)會(huì)有一個(gè)發(fā)展過(guò)程。實(shí)際的回?zé)崞饔捎诹黧w發(fā)展所占的空間較回?zé)崞骺傮w積很小，可以忽略不計(jì)。而模擬中的模型體積過(guò)小，由于流場(chǎng)的入口效應(yīng)，流場(chǎng)可能還沒(méi)進(jìn)入充分發(fā)展的階段，因此對(duì)于微觀結(jié)構(gòu)的模擬結(jié)果造成了較大影響。

(2)模擬中使用的回?zé)崞髂Ｐ蜆O小，忽略了溫度對(duì)回?zé)崞餍阅艿挠绊?。而回?zé)崞鲗?shí)際使用中整體溫度、壓力跨度大，工質(zhì)的物性也變化也較大。該模擬無(wú)法完整的模擬回?zé)崞鞔鬁夭畹奶攸c(diǎn)。所以對(duì)模擬的結(jié)果有一定的影響。

(3)本文模擬的模型尺寸過(guò)小。多孔介質(zhì)具有復(fù)雜微觀結(jié)構(gòu)，離散時(shí)需要大量的網(wǎng)格來(lái)描述，很多時(shí)候遠(yuǎn)遠(yuǎn)超出可承受的范圍，這是以往文獻(xiàn)中很少見(jiàn)到關(guān)于對(duì)回?zé)崞魑⒂^尺度模擬的原因。以疊層絲網(wǎng)計(jì)算模型為例，模擬模型長(zhǎng)度為0.6 mm，而實(shí)驗(yàn)采用的模型為29.88 mm，實(shí)際尺寸約為模擬尺寸的50倍。

5 結(jié)論

(1)建立了層疊絲網(wǎng)式和蝕刻薄片卷裹式回?zé)崞鞯奈锢砟Ｐ?，并使用格子波爾茲曼的等溫D3Q19格子模型對(duì)其進(jìn)行模擬計(jì)算。計(jì)算得到了兩種填料結(jié)構(gòu)的速度場(chǎng)分布和兩端的密度差。并在此基礎(chǔ)上分析了兩種回?zé)崞鹘Y(jié)構(gòu)的對(duì)于流場(chǎng)的影響。

(2)通過(guò)模擬得到填料兩端的密度差，計(jì)算了不同雷諾數(shù)下的穩(wěn)態(tài)流阻系數(shù)。結(jié)果顯示:雷諾數(shù)越大，穩(wěn)態(tài)流阻系數(shù)越小且變化幅度區(qū)域平緩。在相同雷諾數(shù)下，水力直徑和填充率相近的卷裹式回?zé)崞鞯姆€(wěn)態(tài)流阻系數(shù)較層疊絲網(wǎng)式回?zé)崞鞯牡汀?/p>

(3)對(duì)模擬結(jié)果與實(shí)驗(yàn)值進(jìn)行了對(duì)比，模擬結(jié)果在趨勢(shì)上符合實(shí)際，但在壓降的精確計(jì)算上產(chǎn)生了偏差。根據(jù)實(shí)驗(yàn)結(jié)果分析了模擬誤差產(chǎn)生的原因。對(duì)進(jìn)一步通過(guò)格子波爾茲曼方法模擬回?zé)崞鹘Y(jié)構(gòu)有一定的指導(dǎo)意義。

(4)下一步工作展望:1)本文計(jì)算的模型尺寸在0.1—1 mm之間，相對(duì)回?zé)崞鱽?lái)說(shuō)，這個(gè)尺度顯得太小了，速度邊界層和充分發(fā)展邊界條件會(huì)干擾流場(chǎng)內(nèi)的流動(dòng)，四周的壁面也帶來(lái)了額外的沿程阻力損失，使得模擬結(jié)果與實(shí)際有很大的偏差。所以下一步應(yīng)該擴(kuò)大模型的尺寸。此外，格子玻爾茲曼計(jì)算過(guò)程中每個(gè)離散微團(tuán)的獨(dú)立性非常適合顯卡的GPU運(yùn)算，因此并行GPU運(yùn)算是今后格子玻爾茲曼方法求解回?zé)崞魑⒂^流場(chǎng)的一個(gè)可以考慮的方向。2)可壓熱格子玻爾茲曼方法的引入:本文所計(jì)算的模型忽略了回?zé)崞餮刂S線的巨大溫差，假設(shè)其為等溫，采用不可壓等溫格子玻爾茲曼方法計(jì)算穩(wěn)態(tài)壓降。若增加計(jì)算規(guī)模，建立大尺寸的回?zé)崞魈盍夏Ｐ?，則可以引入熱格子玻爾茲曼模型，將溫度的影響列入考慮范圍，進(jìn)一步使模擬結(jié)果接近實(shí)際工況。

1 張曉青，郭方中，董凱軍，等.熱聲熱機(jī)回?zé)崞魈匦詤?shù)影響的模擬分析及優(yōu)化［J］.低溫工程，2000，(3):7-11.

2 章忠敏.回?zé)崞餍实牟罘钟?jì)算［J］.制冷學(xué)報(bào)，1982(3):18-28.

3 Atrey M D，Bapat S L，Narayankhedkar K G.Theoretical analysis and performance investigation of Stirling cycle regenerators［J］.Cryogenics，1991，31(12):1044-1052.

4 He Xiaoyi，Luo Lishi.A priori derivation of the lattice Boltzmann equation［J］.Physical Review E，1997，55(6):R6333-6336.

5 Latt J.Hydrodynamic limit of lattice Boltzmann equations［D］.Genève;Univ.Genève，2007.

6 何雅玲，王 勇，李 慶.格子Boltzmann方法的理論及應(yīng)用［M］.北京:科學(xué)出版社，2009.

7 郭永飛.低溫制冷機(jī)用徑軸向混合填充式回?zé)崞鞯牧鲃?dòng)與回?zé)崽匦匝芯浚跠］.上海:上海理工大學(xué)，2011.

8 王 強(qiáng).新型蝕刻金屬薄片式回?zé)崞鞯牧鲃?dòng)特性研究［D］.上海:上海理工大學(xué)，2012.

9 馬詩(shī)旻.低溫制冷機(jī)回?zé)崞鞯臄?shù)值模擬［D］.上海:上海理工大學(xué)，2012.

10 Silva G，Semiao V.First-and second-order forcing expansions in a lattice Boltzmann method reproducing isothermal hydrodynamics in artificial compressibility form［J］.Journal of Fluid Mechanics，2012，698:282-303.

11 Wolff K，Marenduzzo D，Cates M E.Cytoplasmic streaming in plant cells:the role of wall slip［J］.Journal of the Royal Society Interface，2012，9(71):1398-1408.

12 Wu F，Shi W，Liu F.A lattice Boltzmann model for the Fokker-Planck equation ［J］.Communications in Nonlinear Science and Numerical Simulation，2012，17(7):2776-2790.

Numerical simulation of cryocooler regenerator using Lattice Boltzmann Method(LBM)

Xia Yudong Chen XiMa Shimin Zhang Hua

(School of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China)

Based on Lattice Boltzmann Method(LBM)，the flow field in the microstructure of mesh screen and etched foil regenerators were simulated.The velocity and pressure drop were obtained by LBM.The simulation results show that velocity field in etched foil regenerator is better distributed than that of mesh screen.And the etched foil regenerator has less resistance coefficient than that of mesh screen regenerator.The simulation results basically agree with the experiments.

Lattice Boltzmann Method;regenerator;resistance coefficient;numerical simulation

TB651

A

1000-6516(2012)05-0041-05

2012-07-17;

2012-09-25

國(guó)家自然科學(xué)基金(50906054)，教育部博士點(diǎn)新教師基金項(xiàng)目(20093120120006)資助。

夏宇棟，男，24歲，碩士研究生。