祁影霞，張 燦，車閆瑾，王禹賀，張 華

（上海理工大學能源與動力學院，上海 200093）

0 引言

隨著科學技術(shù)的發(fā)展，航天航空技術(shù)、超導技術(shù)、低溫電子學、低溫醫(yī)學等對低溫制冷技術(shù)的需求與日俱增[1-2]。脈管制冷機作為一種小型制冷機，在低溫制冷領(lǐng)域有著獨特的優(yōu)勢，具有無運動部件、結(jié)構(gòu)簡單、機械振動小、無噪音等優(yōu)點[3]，在上述領(lǐng)域應(yīng)用廣泛。

脈管是脈管制冷機的關(guān)鍵部位連接著冷端換熱器和熱端換熱器，氣體在其中交變流動來實現(xiàn)熱量和冷量的傳遞。雖然過程簡單，但內(nèi)部機理至今也沒有統(tǒng)一的定論。Gifford等[4]研究了基本型脈管制冷機，并提出了“泵熱理論”?；拘兔}管缺乏調(diào)相能力，效率低下，Radebaugh等[5]在此基礎(chǔ)上研究了調(diào)相型脈管制冷機中焓流和壓力、速度之間的關(guān)系并提出了“焓流調(diào)相理論”。后來又出現(xiàn)了如向量分析法[6]、熱聲理論[7-8]、特征線法[9]等新的理論來解釋脈管的制冷機理，使得脈管制冷機得到了很好的發(fā)展。

近年來，脈管制冷機的基礎(chǔ)理論不斷發(fā)展。De Waele等[10-11]在假設(shè)氣庫無限大，且壓力為正弦變化時，獲得了小孔型和雙向進氣型脈管制冷機的冷熱端的動態(tài)溫度變化。Jafarian等[12]運用數(shù)值模擬的手段，通過控制容積方法預測了大冷量脈管制冷機的流體動力學和熱力學特性，并且表明脈管制冷機的運行參數(shù)是影響其制冷效率的主要因素之一。Sexl等[13]研究了脈管制冷機脈管內(nèi)部氣體振蕩流動的壓力梯度。

脈管制冷機在運行時涉及熱力學、流體動力學過程，同時伴有復雜的流動過程，因此其機理研究至今仍不完善。為了有效地提升脈管制冷機的制冷效率，需要進一步研究其機理及影響因素?；诜肿觿恿W方法研究了脈管形狀、粗糙度以及壓力對脈管內(nèi)部溫度場的影響，并給出了脈管內(nèi)部交變流動過程溫度、壓力、密度、速度隨時間的變化趨勢。

1 模擬方法及物理模型

1.1 分子動力學

分子動力學模擬方法（Molecular Dynamics Simulation Method）是當前學者們用來計算規(guī)模龐大且復雜，但系統(tǒng)本身尺度微小的一種方法。隨著力學迅速發(fā)展，研究人員建立了適用于生化分子體系、聚合物、金屬與非金屬材料等不同粒子領(lǐng)域的力場，這些力場的應(yīng)用使得計算結(jié)果精準度大幅提升。分子動力學模擬方法是將這些力場和牛頓力學合起來的一種新型的計算方法。

在數(shù)量龐大的原子運動系統(tǒng)中，系統(tǒng)能量是系統(tǒng)中原子的動能與系統(tǒng)總勢能之和。其表達式為：

式中：U為原子間總勢能；UVDW為非鍵結(jié)合范德瓦耳斯作用力；Uint為分子內(nèi)部勢能。

根據(jù)經(jīng)典力學可知，系統(tǒng)中原子i所受的力為勢能梯度：

由牛頓運動定律可得i原子的加速度為：

4、創(chuàng)新考核機制，保證任務(wù)的落實到位。對于各自決策和負責的事權(quán)，各自權(quán)責內(nèi)恰，管理相對容易。對于中央決策、地方執(zhí)行的事權(quán)，與自主權(quán)下放相適應(yīng)，需要創(chuàng)新監(jiān)督機制，既要保障地方自主權(quán)，又要保障地方執(zhí)行與中央大政方針和戰(zhàn)略目標一致。

開始計算時，固定立方體中的分子數(shù)和工質(zhì)氣體原子的初始位置、初始速度是確定的。系統(tǒng)中全部原子運動動能總和為：

式中：KE為系統(tǒng)的總動能；N為總原子數(shù)；kB為玻爾茲曼常數(shù)；T為熱力學溫度。

原子運動的初始速度可由式（2）計算出，得到原子運動的初始位置和初始速度后，之后每一步會產(chǎn)生當前速度和當前位置，通過當前速度可計算固定格子的溫度。計算公式為：

系統(tǒng)運行過程中，在不同時間段對格子內(nèi)工質(zhì)氣體原子數(shù)目進行統(tǒng)計，可以得到不同時間段內(nèi)軸向壓力分布。聯(lián)合上面的公式，可以計算出各個時間段內(nèi)脈管內(nèi)部軸向速度和溫度場的分布。

1.2 物理模型

建立三種不同壁面的脈管通道，模擬帶有氣庫的脈管制冷機在膨脹時系統(tǒng)內(nèi)部工質(zhì)氣體的流動過程。圖1為幾種通道的模擬物理模型，分別為平直光滑通道、彎曲通道、粗糙通道。所有通道高度為200 nm，長度為730 nm，總長（含熱端）1 000 nm。其中，彎曲通道曲率半徑為100 nm；粗糙通道有四種參數(shù)，分別如表1所列。左側(cè)氣體壓力分別取4.05×105Pa、6.08×105Pa、8.11×105Pa，工質(zhì)為Ar氣，模擬過程中，壁面保持固定，氣體與壁面之間沒有熱交換。首先用正則（NVT）系綜（系統(tǒng)中粒子數(shù)N，體積V，溫度T均為定值），使系統(tǒng)平均溫度保持在298 K。之后采用微正則（NVE）系綜（系統(tǒng)中粒子數(shù)N，體積V，總能E均為定值），使系統(tǒng)與外界沒有能量交換。當系統(tǒng)運行時間在3 000 ps以后，通道內(nèi)氣體基本達到平衡，如圖2所示。

表1 粗糙通道粗糙元參數(shù)Table1 Roughness parameter of rough channels

圖2為Ar氣在運行到3 000 ps時，原子位置分布圖。可以看出，直通道的原子分布沿軸向已經(jīng)較均勻，而彎曲通道次之，粗糙通道更次之。

圖1 三種通道模擬物理模型圖Fig.1 Simulation physical model of three channels

圖2 模擬結(jié)束時原子位置分布圖Fig.2 Atomic position distribution status at the end of simulation

2 分子動力學仿真結(jié)果

2.1 軸向速度分布

圖3、圖4和圖5分別是直通道、彎曲通道、粗糙通道原子軸向運動速度分布圖?？梢钥闯?，通道內(nèi)氣體流動速度在1 260 ps時達到最大，這是由于Ar原子剛開始運動時，通道兩端壓差較大，壓差驅(qū)動其快速運動。彎曲通道和粗糙通道軸向速度分布達到最大值在1 380 ps，可知軸向速度分布受壁面的光滑度和彎曲程度影響。在運動一段時間后，原子基本分布在整個通道內(nèi)，兩端壓差減小，原子運動速度降低。從圖中看出，最終在3 180 ps時原子分布均勻且速度接近0，說明系統(tǒng)已達到穩(wěn)定狀態(tài)。

圖3 直通道原子軸向運動速度分布曲線Fig.3 Axial velocity distribution of straight channel atom

圖4 彎曲通道原子軸向運動速度分布曲線Fig.4 Axial velocity distribution of bent channel atom

圖5 粗糙通道原子軸向運動速度分布曲線Fig.5 Axial velocity distribution of rough channel atom

2.2 通道形狀對溫度場的影響

圖6為三種不同形式的通道，在通道內(nèi)冷端溫度達到最低時的軸向溫度分布。圖中虛線左側(cè)為通道，右側(cè)為通道外（熱端）。可以看出，在冷、熱端之間存在一個溫差，這是因為通道內(nèi)氣體初始壓力高，通道外氣體初始壓力低，因而通道內(nèi)氣體向外膨脹溫度降低，通道外氣體被壓縮因而溫度上升。這就是脈管內(nèi)氣體絕熱放氣制冷的熱力學機理。同時，在冷端和熱端都有一個溫度分布相對平緩的部分。

圖6 三種類型通道沿軸向溫度分布對比圖Fig.6 Comparison of axial temperature distribution of three channels

從圖6還可以看出，和直通道相比，彎曲和粗糙通道冷端溫度大幅提升，且整體軸向溫度分布上移，這是由于在彎曲和粗糙通道內(nèi)，氣體與壁面的摩擦轉(zhuǎn)換為摩擦熱而使氣體溫度上升所致。粗糙度和彎曲度的影響程度幾乎一樣。直通道冷端溫度由初始值的298 K降低到138 K，彎曲通道、粗糙通道分別降低到190 K、215 K?？梢姀澢ǖ馈⒋植谕ǖ缹^熱放氣的制冷溫度影響較大。圖7是與圖6相對應(yīng)的軸向壓力分布圖。可以看出，直通道的壓力分布較平緩，彎曲通道、粗糙通道的壓力分布呈上升趨勢，直通道的最大壓差大約是20 kPa，彎曲、粗糙通道的壓差分別為110 kPa、140 kPa，說明彎曲、粗糙通道的流動阻力較大，因而壓力損失也較大。

圖7 三種類型通道沿軸向壓力分布對比圖Fig.7 Comparison of axial pressure distribution of three channels

2.3 壓力對制冷溫度的影響

圖8 三種類型通道冷端最低溫度隨系統(tǒng)壓力的變化圖Fig.8 The change of the minimum temperature of the cold end with the system pressure in three channels

2.4 粗糙度對溫度場的影響

初始壓力在6.08×105Pa時，運行到3 000 ps后，三種類型通道及不同粗糙度通道冷端達到的最低溫度，如圖9所示。可以看出隨著粗糙度的增加，冷端溫度也相應(yīng)提高，相對于平直光滑通道，case4的冷端溫度最大，提高了82 K。同時，也可以看出，粗糙度的影響比彎曲度對冷端溫度的影響要大。由此可見，通道壁面光潔度對脈管制冷的影響還是非常大的。

圖9 不同粗糙度通道冷端最低溫度比較圖Fig.9 Comparison of the minimum temperature of the cold end in different channels

2.5 通道各參數(shù)隨時間的變化

比較在冷端（150 nm）、熱端（790 nm）處，三種類型通道內(nèi)溫度、壓力、密度、速度隨時間的變化。

圖10和圖11分別為150 nm、790 nm處，三種類型通道內(nèi)溫度隨時間的變化。在圖10中，對于直通道，冷端溫度在2 500～2 800 ps時達到最低，之后有明顯回升；彎曲通道在1 500～2 800 ps、粗糙通道在1 700～2 800 ps時達到較低穩(wěn)定波動狀態(tài)。直通道的冷端溫度要比彎曲和粗糙通道的溫度低很多，而彎曲通道又比粗糙通道的溫度低。在圖11中，位置為790 nm處，該位置處于通道出口處，對應(yīng)熱端溫度。直通道在2 400 ps達到最高，之后開始下降；而彎曲、粗糙通道在2 000 ps之后達到較高穩(wěn)定波動狀態(tài)。也就是說，彎曲、粗糙通道的溫度梯度在2 000 ps之后基本不變了，彎曲和粗糙通道的熱端溫度要遠高于直通道。

圖10 三種類型通道在150 nm處溫度隨時間的變化圖Fig.10 The temperature varies of three channels with time at 150 nm

圖11 三種類型通道在790 nm處溫度隨時間的變化圖Fig.11 Temperature varies of three channels with time at 790 nm

圖12～15分別為150 nm、790 nm處，三種類型的通道內(nèi)壓力、密度隨時間的變化。

從圖12、13可知，在150 nm的位置，三種通道的密度和壓力隨著時間的推移都呈下降的趨勢，這是由于靠近冷端的位置，0 ps時，Ar氣密度較大、壓力較大；在壓差作用下，隨著Ar氣向熱端的流動，靠近冷端處密度下降，壓力降低。對比三種通道，直通道的密度、壓力整體低于彎曲、粗糙通道，這是由于直通道壁面光滑，原子運動較快，密度和壓力下降也較快；彎曲、粗糙通道由于曲率半徑和粗糙度的存在，產(chǎn)生阻力，阻力抵消掉了一部分壓差驅(qū)動力，所以彎曲、粗糙通道內(nèi)原子運動速度較慢，密度、壓力下降也較慢，導致密度和壓力明顯高于直通道的壓力。

圖12 三種類型通道在150 nm處軸向密度分布隨時間的變化圖Fig.12 Axial density distribution varies of the three channels with time at 150 nm

圖13 三種類型通道在150 nm處軸向壓力分布隨時間的變化圖Fig.13 Axial pressure distribution varies of the three channels with time at 150 nm

由圖14、15可知，在790 nm的位置，直通道在900 ps時開始出現(xiàn)密度和壓力波動，彎曲通道和粗糙通道在1 400 ps時才出現(xiàn)壓力和密度波動。這是由于彎曲通道和粗糙通道沿程阻力較大，原子運動進程較緩慢，因而到達790 nm位置所花的時間也較長，密度和壓力上升的也較慢。對比三種通道，直通道的密度、壓力整體高于彎曲、粗糙通道，這是由于直通道內(nèi)原子運動的速度更快，因而密度和壓力也上升的更快。

對比圖12～15可知，粗糙通道的阻力比彎曲通道的要大。

圖14 三種類型通道在790 nm處軸向密度分布隨時間的變化圖Fig.14 Axial density distribution varies of the three channels with time at 790 nm

圖15 三種類型通道在790 nm處軸向壓力分布隨時間的變化圖Fig.15 Axial pressure distribution varies of the three channels with time at 790 nm

圖16和圖17為三種類型的通道在150 nm、790 nm處，軸向速度隨時間的變化。由圖16可知，在150 nm處，軸向速度隨時間逐漸下降，在3 000 ps時，基本處于穩(wěn)定，且接近0。彎曲、粗糙通道內(nèi)原子的運動速度整體低于直通道，說明彎曲和粗糙通道對原子運動的阻力較大，原子會比在直通道中更大概率與壁面發(fā)生碰撞，導致速度下降。由圖17可知，在790 nm處，隨著運行時間的增加，軸向速度也逐漸下降，在3 000 ps時，趨于穩(wěn)定，速度基本接近0，說明系統(tǒng)此時達到平衡。在790 nm處，三種通道的軸向速度基本一致，說明通道阻力對出口速度影響不大。

圖16 三種類型通道在150 nm處軸向速度分布隨時間的變化圖Fig.16 Axial velocity distribution varies of the three channels with time at 150 nm

圖17 三種類型通道在790 nm處軸向速度分布隨時間的變化圖Fig.17 Axial velocity distribution varies of the three channels with time at 790 nm

3 結(jié)論

利用分子動力學模擬了脈管形狀、粗糙度以及壓力對脈管內(nèi)部溫度場的影響，給出了脈管內(nèi)部交變流動溫度、壓力、密度、速度隨時間的變化趨勢。得到了以下結(jié)論：

（1）平均壓力的提高，冷端溫度呈現(xiàn)下降趨勢。故應(yīng)提高系統(tǒng)壓力，進而達到降低冷端溫度的目的。

（2）彎曲、粗糙通道由于曲率半徑和粗糙元的存在，造成運動阻力，氣體無法充分膨脹，冷端壓力不能充分降低，制冷溫度有所上升；導致摩擦，造成氣體整體溫度上升。

（3）粗糙度對冷端溫度的影響比彎曲曲率半徑更大。所以在脈管加工過程中，要保證其表面盡可能光潔。

（4）冷端溫度隨運行時間逐漸下降，直至穩(wěn)定；熱端溫度隨運行時間逐漸上升，直至穩(wěn)定。直通道的冷、熱端溫度均顯著低于彎曲、粗糙通道。冷端壓力、密度隨運行時間逐漸減小，直至達到平衡，直通道的冷端壓力、密度低于彎曲、粗糙通道；熱端壓力、密度隨運行時間逐漸上升，直通道的熱端壓力、密度要顯著高于彎曲、粗糙通道。