微通道脈沖管中He氣交替振蕩的分子動力學(xué)模擬

蘆 洋， 祁影霞,2， 陸熠璇, 車閆瑾， 劉雅麗， 陳衛(wèi)誼

(1.上海理工大學(xué) 能源與動力工程學(xué)院， 上海 200093; 2.上海市動力工程多相流動與傳熱重點實驗室， 上海 200093)

1 引 言

斯特林脈沖管制冷機具有結(jié)構(gòu)簡單、活動部件少、體積小、重量輕、使用壽命長等特點. 可滿足醫(yī)療器械，紅外探測器，超導(dǎo)量子干涉儀等精密儀器在低溫工作環(huán)境的要求[1, 2]. 近幾年，專家們對微型脈沖管制冷機進(jìn)行了大量的研究[3-7]，Tang用100 Hz的頻率來驅(qū)動直徑為5 mm的脈沖管，該脈沖管的空載溫度為49.5 K，在80 K時的冷卻能力為1.7 W[8]. 研究結(jié)果表明，為了滿足微型脈沖管制冷機的冷卻要求，通常使用高頻線性壓縮機來提供功率[3, 9, 10]，高頻脈沖管制冷機具有能量流密度大、熱聲效率高、振動和噪聲小等特點. 因此，脈沖管制冷系統(tǒng)可以通過提高工作頻率來進(jìn)一步小型化. 對MPT中的微流過程和交替振蕩的氣體特性的研究可以增進(jìn)對制冷機理的了解，為微通道脈沖管制冷機的結(jié)構(gòu)優(yōu)化提供基礎(chǔ).

Gan等人研究了線性壓縮機的阻抗特性，優(yōu)化了線性壓縮機的運行機理和低溫制冷系統(tǒng)中的阻抗匹配[11]. 楊森等人通過CFD仿真研究了脈管預(yù)冷對脈管內(nèi)部溫度場和速度場的影響，對脈管進(jìn)行預(yù)冷會改變脈管工作時的內(nèi)部溫度分布,進(jìn)而提高二級脈管制冷能力[12]. 巢偉等人利用REGEN 設(shè)計了一臺單級大功率斯特林型脈管制冷機的回?zé)崞? 采用功效系數(shù)法對制冷量和COP進(jìn)行了多目標(biāo)優(yōu)化,使得系統(tǒng)的性能最佳[13]. 王海敏等人采用SAGE對制冷機進(jìn)行模擬，引入虛擬的振子阻尼調(diào)相機構(gòu)對無負(fù)荷制冷溫度進(jìn)行模擬，結(jié)果表明30K時出現(xiàn)最佳制冷量[14]. 綜上所述，向量分析方法主要使用正弦壓力波和質(zhì)量流波來演示小壓力比的周期性物理特性，簡化數(shù)值計算過程[11, 15]. CFD模擬方法需要預(yù)先設(shè)置恒定的進(jìn)口邊界值以及經(jīng)驗數(shù)據(jù). 當(dāng)流體被假定為連續(xù)介質(zhì)時，校正系數(shù)用于描述湍流模型，因此與實際的交替振蕩過程相比可能會有誤差[16-20]. REGEN用于計算再生器[21, 22]，SAGE可以基于熱聲理論對整個設(shè)備進(jìn)行一維數(shù)值模擬. 但是二次流和渦流的描述可能不準(zhǔn)確[23, 24]，脈沖管內(nèi)部的溫度難以獲得. 因此，壓縮機聲功率和冷端負(fù)荷的有效耦合方法理論仍然缺乏. 分子動力學(xué)模擬遵循經(jīng)典的牛頓力學(xué)定律，可以重建原子運動演化過程. 非平衡分子動力學(xué)(NEMD)通常用于評估非穩(wěn)態(tài)原子塊，例如施加的力場，電場和磁場. NEMD可以在具有復(fù)雜流和許多分子的微系統(tǒng)中實現(xiàn)[25].

另外，實驗中實際的流動過程和脈管參數(shù)的演變也難以監(jiān)控和獲取，阻礙了MPT數(shù)值模擬. 分子動力學(xué)方法只需要提供初始狀態(tài)，例如時間，充氣壓力和整體溫度，而不是操作過程中的過程點. 因此，在微觀水平上用數(shù)字描述了管內(nèi)溫度梯度的變化. 結(jié)果表明，脈沖管的軸向壓力和密度分布不均勻，這也是沿軸向溫度梯度變化的原因. MPT氣柱的聲振動頻率必須與線性壓縮機的固有頻率結(jié)合在一起，以建立合理的聲場，以確保脈沖管的絕熱膨脹效率. 由于缺乏理論和數(shù)值分析，實驗數(shù)據(jù)被用于基于壓縮機的掃氣容積與MPT容積之比來指導(dǎo)設(shè)計. 因此，本文使用分子動力學(xué)方法建立了線性壓縮機和脈沖管的模型，以便在微觀水平上研究MPT的聲場和溫度梯度.

2 分子動力學(xué)模擬方法

分子的總能量是總勢能U和總動能EK的總和. 總勢能U通常包括分子的范德華力和內(nèi)部勢能. 分子之間的范德華力通常由力場描述. 這項研究中，在He原子之間使用了Lennard-Jones勢能. 其中包括12次冪的排斥項與6次冪的吸引項. 如下所示：

(1)

式中，ε(eV)為原子間作用力，取0.000607098；σ(?)為原子間距離，取2.103. 則系統(tǒng)內(nèi)總勢能可表示為：

(2)

通過原子的初始位置獲得系統(tǒng)的總勢能之后，就可以計算出原子在系統(tǒng)中的力和加速度[27, 28]. 并對時間進(jìn)行積分，獲得t+δt時刻的分子位置與速度,

(3)

(4)

Verlet算法用于求解遞歸方程，并根據(jù)位置的差分獲得速度[10]：

(5)

(6)

當(dāng)系統(tǒng)的平均溫度用于賦予下一次循環(huán)中原子的初始速度，并由該時刻原子的坐標(biāo)得到新的系統(tǒng)勢能與動能，進(jìn)行迭代計算：

(7)

(8)

當(dāng)系統(tǒng)中設(shè)置為(NVE)系綜時，那么系統(tǒng)中的壓強是變量，相對應(yīng)的公式可以計算出其壓強與時均溫度：

(9)

(10)

3 仿真模型分析

為了研究MPT的充氣和放氣過程之間的耦合機制，建立了包括線性壓縮機和MPT的三維絕熱模型. 使用非平衡分子動力學(xué)進(jìn)行模擬計算. 壓縮機內(nèi)使用正弦速度活塞提供He氣動力. 在X維度上使用了周期性邊界，在Y(Yi和Yh)和Z(Zi和Zh)方向上使用了固定邊界. 為了在軸向上區(qū)分結(jié)果，將模擬盒子沿管每100 nm分成一個格子，以獲得實時的物理特性，例如溫度和壓力. Yh，Yi和Zh邊界設(shè)為反彈壁面，使原子與反射壁之間發(fā)生完全彈性碰撞. 活塞直徑和MPT直徑相同. MPT尺寸為5.72 ?*2000 ?*17160 ? (x*y*z)，其中充滿He氣. 將移動的反射正弦速度壁面作為活塞放置在Zi處. 活塞的運動方向定義為Zi→Zh，其位移為572nm，因此其瞬時功能位置如公式11所示，

Z=A[1-cosω(tp-t0)dt]

(11)

圖1 在t=0和時的分子動力學(xué)模型Fig. 1 Molecular dynamics models at t=0 and

首先，在NVT系綜(原子數(shù)、模型體積、溫度不變)中，標(biāo)定10萬步以獲得300 K的均勻溫度分布. 接下來，在NVE系綜中使用絕熱模型來模擬MPT在0.4 fs的時間步長處的充氣和放氣過程. 根據(jù)Qi的研究[7]，1H-3L的脈管在壓比為2時自然振蕩周期為2240 ps. 振蕩時間受充氣壓力[8, 26, 27]影響. 本研究中使用的時間和充氣壓力如表1所示(每個模型的名稱包括兩部分，分別涉及充氣壓力和時間. 例如，在A1中，模型的充氣壓力為 10 bar，時間為4000 ps).

表1 模型參數(shù)

4 結(jié)果與討論

4.1相位角

MPT中壓縮與膨脹過程為等熵絕熱過程. 由公式12可知冷卻能力受壓力波和質(zhì)量流之間的相角影響[28, 29].

(12)

如圖1所示，模型的左端是冷端，右端是熱端. 如圖2、圖3所示，定義壓力波超前質(zhì)量流時相角為正值. 由仿真結(jié)果可知相角在冷端為正、熱端為負(fù). 冷端活塞位移與冷端壓力波之間的相角為54°，仿真數(shù)據(jù)與實驗結(jié)果一致[12,17,24].

圖2 模型C1活塞位移與冷端壓力波之間的相角Fig. 2 The phase shift between the displacement of piston and the pressure curve at cold end of model C1

圖3 模型C1冷端壓力與質(zhì)量流量之間的相角Fig. 3 The phase shift between the pressure curve and mass flow curve at cold end of model C1

4.2 MPT的冷卻機制

圖4 模型C1壓縮過程通道內(nèi)軸向壓力分布Fig. 4 The axial pressure distribution during compression process of model C1 in the tube

圖5 模型C1膨脹過程通道內(nèi)軸向壓力分布Fig. 5 The axial pressure distribution during expansion process of model C1 in the tube

通過對壓縮和膨脹過程作對比可獲知，聲波前壓力梯度分布的動態(tài)傳遞過程是不對稱的. 這是由于封閉的表面和氣體消散聲波的能量引起的. 在膨脹過程中， 活塞的反向運動導(dǎo)致冷端附近產(chǎn)生壓差，產(chǎn)生絕熱膨脹過程，從而導(dǎo)致冷端溫度遠(yuǎn)低于300 K，而熱端溫度遠(yuǎn)高于300K. 在這種動態(tài)過程中，MPT中形成了不均勻的溫度梯度，小孔型脈沖管制冷機(OPTC)和慣性管型脈沖管制冷機(ITPTC)將聲波慣性從熱端通過孔板或慣性管耗散到氣庫，而不是直接與封閉端造成碰撞.

圖6 模型C1壓縮過程通道內(nèi)軸向溫度分布Fig. 6 The axial temperature distribution during compression process of model C1 in thetube

圖7 模型C1膨脹過程通道內(nèi)軸向溫度分布Fig. 7 The axial temperature distribution during expansion process of model C1 in the tube

圖8 模型C1中混流過程通道內(nèi)的軸向壓力分布Fig. 8 The axial pressure distribution during mixed flow process of model C1 in the tube

4.3 瞬時平均參數(shù)分析

如圖10所示，測定循環(huán)中的冷端瞬時平均溫度. 根據(jù)這項研究的模型，當(dāng)τ≥6000 ps時，冷端的瞬時平均溫度變化很小. 這表明MPT的聲場建立過程與強制對流和氣體振蕩過程中壓力梯度引起的自然對流有關(guān). 自然對流強度與MPT的軸向壓力梯度呈正相關(guān). 當(dāng)Po≤15 bar時，瞬時平均溫度受強制對流時間的影響較小. 但是，在高壓條件下(20 bar-30 bar)，當(dāng)強制對流時間短于自然對流時間時，冷端溫度隨時間的降低而降低，如圖11所示，最終進(jìn)一步增加. 將回?zé)崞鞣胖迷贛PT的前面，在反復(fù)的氣體壓縮和膨脹過程中，冷端的溫度不斷降低.

圖10 不同時刻與充氣壓力下冷端的瞬時平均溫度Fig. 10 Instantaneous average temperatures of the cold junction at different times and inflation

圖11 不同時刻和充氣壓力下冷端和熱端之間的時均溫差Fig.11 Time-averaged temperature differences between cold and hot ends under different times and inflation pressures

綜上所述，適當(dāng)提高充氣壓力對冷卻效果是有益的. 當(dāng)運行頻率高于MPT的固有頻率時，冷端溫度會進(jìn)一步降低，但熱端溫度較高，因此對換熱器和調(diào)相元件的散熱量需求更大. 另外，回?zé)崞魍ǔＬ畛浣饘倬W(wǎng)或多孔介質(zhì)，較高的運行頻率會相應(yīng)產(chǎn)生更高的流動阻力，增加回?zé)崞鞯哪Σ翐p失與壓降損失，因此，低流動阻力和高回?zé)嵝实幕責(zé)崞鲬?yīng)成為下一階段研究的重點. 該模型忽略了軸向表面的傳熱損失和表面之間的摩擦損失，從而導(dǎo)致輕微的誤差[30].

5 結(jié) 論

用分子動力學(xué)方法模擬了微通道結(jié)構(gòu)脈管中氦氣的交變振蕩過程，建立了宏觀與微觀狀態(tài)參數(shù)之間的關(guān)聯(lián)式. 仿真結(jié)果揭示了高頻微通道脈沖管制冷機的微觀機理，為設(shè)計MPT制冷機的結(jié)構(gòu)和運行參數(shù)提供了理論依據(jù).

(1) 隨著時間的減少，MPT中出現(xiàn)了壓力波傳輸特性，并且主壓力梯度分布于波前. 在壓縮過程中，從冷端傳遞到熱端的壓力波保持恒定的壓力梯度值，在封閉端形成負(fù)壓梯度且數(shù)值顯著增加. 在膨脹過程中，壓力梯度會隨著時間連續(xù)下降. 隨著時間的增加，平均壓力分布與軸向壓力梯度分布交替出現(xiàn)，梯度值恒定的壓力波消失，因此多重振蕩過程影響了冷卻效果.

(2) 聲波傳輸效應(yīng)是建立冷卻效應(yīng)和軸向溫度梯度的關(guān)鍵. 冷端溫度對強制對流時間和充氣壓力敏感，對高壓更敏感. 當(dāng)強制振蕩周期的時間短于與自然振蕩周期相同的壓力時，冷端的溫度進(jìn)一步降低，而熱端的溫度升高. 這也是為什么高頻斯特林脈沖管制冷機在冷端具有更高的冷卻能力和更低的空載溫度的原因.

(3) MPT相移約為20°至40°. 最佳高充氣壓力周期小于低充氣壓力周期，較短的時間導(dǎo)致聲波傳輸更大的能量流密度，因此聲波兩側(cè)的壓力梯度值都較大. 較大的壓比有利于制冷，但必須具有較高性能的換熱器. 仿真結(jié)果揭示了高頻微通道脈沖管制冷機的微觀機理，為設(shè)計MPT制冷機的結(jié)構(gòu)和運行參數(shù)提供了理論依據(jù).