夏宇棟 陳 曦 馬詩旻 張 華

(上海理工大學能源與動力學院 上海 200093)

基于格子波爾茲曼方法的回熱器數值模擬

夏宇棟 陳 曦 馬詩旻 張 華

(上海理工大學能源與動力學院 上海 200093)

利用格子玻爾茲曼方法，直接對蝕刻薄片和層疊絲網回熱器的微觀結構流場進行了模擬。得到了兩種回熱器填料的微觀流場和兩端的壓差。模擬結果顯示，當回熱器的直徑、水力直徑和填充率相近情況下，不同流速下蝕刻薄片卷裹式回熱器的穩(wěn)態(tài)阻力系數均比層疊絲網回熱器小。穩(wěn)態(tài)阻力系數的模擬變化趨勢與實驗一致。

格子波爾茲曼方法 回熱器 流阻系數 數值模擬

1 引言

回熱器是低溫制冷機的關鍵部件，也是影響制冷機性能的最重要因素。目前，對回熱器的研究主要有實驗研究、傳統(tǒng)計算方法研究和模擬計算方法研究。隨著計算機技術的迅速發(fā)展，對回熱器的模擬計算已經成為了回熱器研究領域的一個新的熱點。張曉青，郭方中等［1］利用基于網絡模型開發(fā)了仿真軟件對熱聲驅動的制冷機用回熱器的長度、水力半徑以及回熱器在駐波聲場中的位置等參數對其性能的影響進行了模擬分析，并對其進行了優(yōu)化，為更深地理解回熱器的機理和實質提供了數值基礎。章忠敏［2］根據豪森的數學模型，利用有限差分法求解對稱情況下的理想回熱器傳熱微分方程，并編制了計算穩(wěn)流回熱器效率的通用程序，分析了理想穩(wěn)流回熱器的傳熱特點和近似的數學描述。M.D.Atrey等［3］運用有限差分法重點分析了回熱器對整機性能的影響，分析給出了工質質量流率、平均壓力、流速、傳熱周期等對回熱器效率的影響，得出回熱器效率是絲網孔隙率的函數。

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

本文就將運用格子波爾茲曼方法對層疊絲網和蝕刻薄片卷裹式回熱器填料的微觀機構進行數值模擬，結合宏觀尺度上的回熱器數值計算模型，對回熱器性能做出預測。

2 回熱器物理模型建立

格子玻爾茲曼方法需將計算空間劃分為離散的流體微團。因此首先得建立兩種回熱器的物理模型。

2.1 層疊絲網回熱器物理模型建立

層疊絲網是回熱器填料主要是采用高目數的不銹鋼或磷青銅層疊而成，其結構形式如圖1所示。

圖1 層疊絲網回熱器填料的結構Fig.1 Structure of mesh screen regenerator

其結構特點是徑向為規(guī)整的編織網片，軸向為網片的隨機旋轉疊加。本次模擬選用的層疊絲網的詳細幾何參數見表1所示。

表1 層疊絲網回熱器幾何參數Table 1 Structure parameters of mesh screen regenerator

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

圖2 三維層疊絲網回熱器填料結構Fig.2 Three-dimensional model of meshscreen regenerator

2.2 蝕刻薄片卷裹式回熱器物理模型的建立

蝕刻薄片卷裹式回熱器是由本課題組自主研發(fā)的新型回熱器，它是在金屬薄片表面采用化學腐蝕或光刻的方式形成一定的圖案，通過卷裹形成流道。其幾何尺寸如圖3所示。

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

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

該結構的計算空間為0.275 mm×1.4 mm×2.6 mm，取網格尺寸δx=5×10-6m，則計算的網絡節(jié)點約為800.8萬個。

圖4 層疊絲網回熱器填料結構Fig.4 Three-dimensional model of etched foil regenerator

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

運用開源的CFD框架Palabos可以對于這兩種回熱器結構進行模擬。Palabos是由日內瓦大學Jonas Latt博士、Orestis博士等成員聯合開發(fā)的開源的通用CFD計算庫，它基于格子玻爾茲曼法。2011年9月，由FlowKit公司開始運營，為用戶提供解決方案。它基于C++的標準模板庫(STL)，提供了完整的程序庫(Library)框架，用戶可以在其框架的基礎上開發(fā)基于格子玻爾茲曼的計算程序。

3.1 回熱器計算模型設定

本次計算的回熱器的軸向長度僅為0.6 mm和2.6 mm，其兩端溫差較小，因此忽略溫度變化，即模擬的結構為一個等溫的模型。所以計算選用等溫D3Q19-BGK模型，其恢復的方程與標準的不可壓Navier-Stokes方程組一致，適用于低馬赫數的流動。

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

3.2 運行結果處理及分析

對層疊絲網結構回熱器，為了能有效地與實驗進行對比，根據實際的回熱器工作參數，模擬選擇了入口速度為2、5、8、11 m/s 4 個工況點。

3.2.1 流場分布

圖5所示為入口速度為5 m/s時，兩種填料結構的速度場向量圖。從圖中可以清晰的看出工質在進入層疊絲網結構前速度較小，當工質進入填料間的流道，由于流道變窄，工質加速，最大速度可達到入口速度的7—10倍大小。

圖6為兩種結構回熱器填料流道內的速度分布向量圖，兩種回熱器填料結構的差異使其流場有很大的不同，進而影響回熱器的回熱性能。

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

圖6 入口速度為5 m/s的流道內速度場分布Fig.6 Velocity field of flow pass at inlet velocity of 5 m/s

從圖中可以分析看出，對于層疊絲網結構的回熱器，絲網疊加的隨機性帶來了流動的隨機性，使得其性能很難預測。即使是相同材料、相同孔隙率、相同大小、相同水力直徑的回熱器，由于絲網疊加形式的不可控，也會造成性能的不同。而對于蝕刻薄片卷裹式回熱器來說，其速度分布更為規(guī)則和均勻，沒有出現圖6(a)中部分填料區(qū)域沒有發(fā)生換熱的情況。工質在平行于流動方向的流道只受到邊界粘應力的影響，其阻力較小。而斜流道的出現，適時地打破熱邊界層對換熱的影響，同時對流阻貢獻不大，可以看到，工質流過斜流道后到達直流道，速度迅速達到恢復原來的大小。這種流場換熱均勻，流動阻力小，對提高回熱器的換熱效率，降低流阻教層疊絲網回熱器有較大的優(yōu)勢。

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

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

圖7 層疊絲網回熱器和蝕刻薄片卷裹式的流阻系數的對比Fig.7 Resistance coefficient vs Reynolds number of mesh screen and etched foil regenerator

式中:ρ為工質密度，(kg/m3)，qv為回熱器內工質體積流量，m3/s;dh為回熱器水力直徑，m;μ為工質動力粘度，Pa·S;d為回熱器直徑，m;φ為孔隙率;L為回熱器長，m。

由圖7分析可得:隨著工質雷諾數Re的增加，兩種回熱器的工質穩(wěn)態(tài)流動阻力系數f均減少，層疊絲網回熱器的阻力系數比蝕刻薄片卷裹式回熱器大，說明蝕刻薄片卷裹式回熱器有流動阻力小的特點。

4 回熱器穩(wěn)態(tài)流阻測試及分析

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

如圖8所示為穩(wěn)態(tài)測試系統(tǒng)中的實驗與測試系統(tǒng)圖。測試時壓力為0.3 MPa，溫度為290 K。測試對象為280目的層疊絲網回熱器和蝕刻薄片卷裹式回熱器，工質為氦氣和氮氣。

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

4.2 LBM模擬數據與實驗結果比較

將運用LBM的模擬數據與實驗所測得的數據進行比較，其比較結果如圖9所示。

圖9 回熱器模擬值和實驗值的對比Fig.9 Value of simulation vs value of experiment

從圖中可以看出，模擬值在趨勢上是符合實驗值的，但出現了整體的偏移。模擬值較實驗值偏小50%左右。

上述模擬結果與實驗值的對比可以看出，模擬結果在趨勢上符合實際，但在壓降的精確計算上產生了偏差。本文認為，可能是由以下因素造成的:

(1)流體進入填料結構時，流場會有一個發(fā)展過程。實際的回熱器由于流體發(fā)展所占的空間較回熱器總體積很小，可以忽略不計。而模擬中的模型體積過小，由于流場的入口效應，流場可能還沒進入充分發(fā)展的階段，因此對于微觀結構的模擬結果造成了較大影響。

(2)模擬中使用的回熱器模型極小，忽略了溫度對回熱器性能的影響。而回熱器實際使用中整體溫度、壓力跨度大，工質的物性也變化也較大。該模擬無法完整的模擬回熱器大溫差的特點。所以對模擬的結果有一定的影響。

(3)本文模擬的模型尺寸過小。多孔介質具有復雜微觀結構，離散時需要大量的網格來描述，很多時候遠遠超出可承受的范圍，這是以往文獻中很少見到關于對回熱器微觀尺度模擬的原因。以疊層絲網計算模型為例，模擬模型長度為0.6 mm，而實驗采用的模型為29.88 mm，實際尺寸約為模擬尺寸的50倍。

5 結論

(1)建立了層疊絲網式和蝕刻薄片卷裹式回熱器的物理模型，并使用格子波爾茲曼的等溫D3Q19格子模型對其進行模擬計算。計算得到了兩種填料結構的速度場分布和兩端的密度差。并在此基礎上分析了兩種回熱器結構的對于流場的影響。

(2)通過模擬得到填料兩端的密度差，計算了不同雷諾數下的穩(wěn)態(tài)流阻系數。結果顯示:雷諾數越大，穩(wěn)態(tài)流阻系數越小且變化幅度區(qū)域平緩。在相同雷諾數下，水力直徑和填充率相近的卷裹式回熱器的穩(wěn)態(tài)流阻系數較層疊絲網式回熱器的低。

(3)對模擬結果與實驗值進行了對比，模擬結果在趨勢上符合實際，但在壓降的精確計算上產生了偏差。根據實驗結果分析了模擬誤差產生的原因。對進一步通過格子波爾茲曼方法模擬回熱器結構有一定的指導意義。

(4)下一步工作展望:1)本文計算的模型尺寸在0.1—1 mm之間，相對回熱器來說，這個尺度顯得太小了，速度邊界層和充分發(fā)展邊界條件會干擾流場內的流動，四周的壁面也帶來了額外的沿程阻力損失，使得模擬結果與實際有很大的偏差。所以下一步應該擴大模型的尺寸。此外，格子玻爾茲曼計算過程中每個離散微團的獨立性非常適合顯卡的GPU運算，因此并行GPU運算是今后格子玻爾茲曼方法求解回熱器微觀流場的一個可以考慮的方向。2)可壓熱格子玻爾茲曼方法的引入:本文所計算的模型忽略了回熱器沿著軸線的巨大溫差，假設其為等溫，采用不可壓等溫格子玻爾茲曼方法計算穩(wěn)態(tài)壓降。若增加計算規(guī)模，建立大尺寸的回熱器填料模型，則可以引入熱格子玻爾茲曼模型，將溫度的影響列入考慮范圍，進一步使模擬結果接近實際工況。

1 張曉青，郭方中，董凱軍，等.熱聲熱機回熱器特性參數影響的模擬分析及優(yōu)化［J］.低溫工程，2000，(3):7-11.

2 章忠敏.回熱器效率的差分計算［J］.制冷學報，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方法的理論及應用［M］.北京:科學出版社，2009.

7 郭永飛.低溫制冷機用徑軸向混合填充式回熱器的流動與回熱特性研究［D］.上海:上海理工大學，2011.

8 王 強.新型蝕刻金屬薄片式回熱器的流動特性研究［D］.上海:上海理工大學，2012.

9 馬詩旻.低溫制冷機回熱器的數值模擬［D］.上海:上海理工大學，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

國家自然科學基金(50906054)，教育部博士點新教師基金項目(20093120120006)資助。

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