基本型脈管內氣體振蕩制冷機理的分子動力學模擬

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

在科學技術迅猛發展的時代，低溫制冷技術的發展具有重大意義，而脈管制冷機在低溫制冷技術中有著不可替代的位置。脈管制冷機與斯特林、G-M制冷機不同，不存在運動部件,所以不會在可靠性、振動、電磁干擾(EMI)等方面影響其制冷性能[1-3]。脈管制冷機被廣泛應用于航空航天、軍事、低溫電子學、低溫醫學、氣象、氣體液化等領域。

雖然脈管制冷機結構較為簡單，但由于脈管制冷機在運行時涉及許多熱力學知識，如熱力學、流體動力學、傳熱學等，其內部流動過程復雜，所以其內部機理的研究至今仍不完善，還有待于進一步研究。當前，脈管制冷機的熱力學原理存在以下幾種不同的理論和解釋。W. E. Gifford等[4]提出脈管制冷機的表面泵熱理論，解釋了脈管制冷機制冷機理，是氣體微團在壓縮膨脹過程中沿脈管壁面從冷端向熱端的泵熱，但其理論只適用于基本型脈管制冷機的分析。E. I. Mikulin等[5]提出經典分析法，此方法是在理想情況下，利用熱平衡方程和已知的熱力學關系式，計算出循環完成后有效制冷量的表達式，最后求得改進型脈管制冷機的總制冷量和效率。R. Radebaugh等[6]提出焓流調相理論，該理論運用能量守恒和質量守恒定律進行研究，揭示了脈管制冷機中焓流和壓力、速度之間的關系。焓流調相理論闡釋了小孔脈管制冷機理,尤其是作用在小孔和氣庫之間的調相。P. J. Storch等[7]提出向量分析法，在壓比較小時，可通過旋轉向量來表示正弦向量。但此種向量分析法較為理想化，表達式僅在計算結果與實驗結果定性時，具有較好的符合度。通過向量分析法可揭示脈管制冷機內部參數變化和相互關系。N. Rott[8]建立了熱聲理論，脈管制冷機依賴于工作流體的振蕩與固體介質熱相互作用產生時均冷量效應，但局限于駐波聲場的熱聲效應。肖家華[9]進一步研究了熱聲理論，將其推廣到一般聲場，并可利用熱聲理論對脈管制冷機進行分析和計算。Liang Jingtao等[10-12]提出關于熱力學非對稱理論。通過建立模型，分別將脈管內的氣體分為靠近管壁的熱黏層內的氣體和脈管中心空間內的氣體，制冷機理分別為表面泵熱效應和熱力學非對稱效應。脈管的制冷量為這兩部分氣體產生的制冷量總和。Wu P. 等[13]提出特征線法，此理論利用氣體動力學中數學描述方法，揭示脈管制冷機內部流動過程。Zhu Shaowei等[14]通過數值模擬方法，得出脈管制冷機整機性能及脈管內部流動和傳熱動態特性等，且對實驗中產生的各種不可逆因素進行數值模擬，可以合理地解釋實驗結論。除了上述研究脈管制冷機理的方法，目前用到的理論方法還有電工學類比分析、向量法、網絡分析法、脈管制冷氣體對外做功分析法[15]等。

本文基于以上脈管制冷機的理論研究現狀，采用分子動力學模擬的方法，通過對微通道的充氣和放氣過程，研究基本型脈管制冷機中脈管內部氦氣被壓縮和膨脹的熱力學過程，可得到脈管內軸向分布的密度、壓力、速度和溫度隨時間變化的結果，這是首次對脈管內部的瞬時溫度變化進行研究，從而進一步探討脈管制冷機中脈管內部的溫度梯度形成的機理及影響因素；也為運用分子動力學模擬方法對目前廣泛應用的調相型脈管制冷機機理研究提供了可能。

由于CFD等模擬方法存在不足，如計算模型復雜、計算過程中的嚴重耦合及為減小計算量所作的一些簡化，使計算結果與實驗結論存在一定差距。本文采用分子動力學模擬方法進行研究，模型簡單，省略了復雜計算，且實驗中無法獲得的微觀細節可較容易得到。

1 脈管制冷機機理

1.1 脈管制冷機熱力學性質

圖1所示為基本型脈管制冷機循環過程氣體溫度分布。高壓氣體在經過回熱器時，制冷溫度被冷卻至Tc，進入冷端交換器制冷。經過層流化原件后使氣體變為層流狀態，再進入到脈管中，即壓縮過程。該過程中氣流以層流形式運動，并受到后續氣體的壓縮，繼續進入氣體，溫度升高，在封閉端溫度升高至Tm，由于熱交換器作用，封閉端處被水冷卻，溫度降至Ta。之后為膨脹過程，管內氣體分布如圖中的虛線所示，高壓氣體流出脈管，因降壓膨脹導致降溫。這時在冷端交換器的氣體溫度降至Ti，Ti

圖1 基本型脈管制冷機中壓縮和膨脹過程的溫度分布Fig.1 Temperature distribution of compression and expansion in a basic pulse tube cryocooler

1.2 焓流調相理論

焓流調相理論解釋了小孔和氣庫的工作原理，通過調節脈管冷端壓力波和質量流之間的相位差而提高制冷效率。脈管制冷機焓流理論公式如下：

(H)∝|p||u|cosθ

(1)

式中：p為管內的氣流壓力波，N；u為速度波，m/s；θ為相位角。

當θ=0°時,即氣流壓力波和速度波同相時,制冷量最大；當θ=90°時,制冷效果最差。故為了提高脈管制冷機制冷能力，可減小脈管內氣流壓力波和速度波之間的相位差。通過對脈管內氣流壓力波和速度波的對比，最佳長徑比可增加制冷量和制冷效率。

2 分子動力學模擬

2.1 分子動力學模擬計算

分子動力學模擬[16]是目前被廣泛采用的計算龐大復雜系統的方式。自1970年起，由于力學發展迅速，建立了許多適用于生化分子體系、聚合物、金屬與非金屬材料的力場，使計算結果精準度大大提升。分子動力學模擬是應用這些力場與牛頓力學原理所發展的計算方法。

對于含有許多原子的運動系統，系統能量為系統中原子的動能與系統總勢能的總和。總勢能為：

U=UVDW+Uint

(2)

式中：U為原子間總勢能，J；UVDW為非鍵結范德瓦爾斯作用力,J；Uint為分子內部勢能,J。

根據經典力學可知，系統中i原子所受的力為勢能梯度：

(3)

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

(4)

(5)

(6)

(7)

自行計算時，已知固定立方體中的分子數和立方體中的氦原子的初始位置及速度。系統中所有原子運動動能總和為：

(8)

式中：KE為系統的總動能，J；N為總原子數；kB為玻爾茲曼常數；T為熱力學溫度， K。

原子運動初速度可由式(2)得出，產生原子的起始位置與速度后，之后每一步產生新的速度和位置，通過新產生的速度可計算固定格子的溫度為：

(9)

式中：Tcal為固定格子內的平均溫度，K。

系統運行中，在不同時間段統計格子內He原子數目，可得不同時間段內軸向壓力分布，并計算出脈管制冷機中各時間段內脈管內部軸向速度與溫度的分布。

2.2 分子動力學模擬方法

本文通過建立微通道外高壓He氣向含有低壓He氣的微通道內流動的模型，模擬脈管制冷機在壓縮時脈管內部的He氣流動過程；建立微通道內高壓He氣向含有低壓He氣的微通道外流動的模型，模擬脈管制冷機在膨脹時脈管內部的He氣流動過程。通過壓縮模型和膨脹模型研究脈管內部的熱力學性質。建立模型步驟：首先建立以低壓(100 kPa)He原子(半徑為0.122 nm)為流體介質的穩態模型和以高壓(1 300 kPa)He原子流體介質的穩態模型。然后分別讓高壓穩態模型和低壓穩態模型在NTV(粒子數N、體積V、溫度T)正則系綜下運行，分別使兩個穩態模型內部氣體均勻混合；其次分別建立軸向Fe壁面和縱向Fe壁面。最終模型如圖2所示，分別為微通道壓縮初始物理模型和微通道膨脹初始物理模型，模型尺寸均為5732.921?(X)×14.303?(Y)×1017.442?(Z)。

圖2 微通道壓縮、膨脹初始物理模型 (左側為微通道)Fig.2 Initial physical model of channel compression and expansion (the left is the channel)

模型建好后，需要設置其勢能函數，一般情況在Fe-Fe原子之間設置Johnson勢能函數[17]，在Fe-He原子之間和He-He原子之間設置UFF勢能函數[18]。

圖3 系統運行最終時氣體分布狀態Fig.3 Gas distribution status at the end of system operation

系統首先在NTV正則系綜下運行，賦予該系統初始溫度，運行結束后，系統內部溫度穩定在300 K。然后在巨正則系綜NVE(E為系統的能量)下運行，最終使系統內部氣體分布達到均勻，如圖3所示。運行時間步長為0.8 fs，每1 000步輸出一次結果。

該模擬采用了兩個限定條件：Fe原子固定，使鐵壁面絕熱，氣體與鐵壁面之間無熱量交換；氦氣體為理想氣體。

3 模擬結果及分析

3.1 軸向密度與壓力分布隨時間的變化

對所建的壓縮和膨脹模型在軸向方向上進行劃分，從左向右，每200 ?劃為一格。然后統計不同時間段內每個格子內的He原子數量，通過He原子數量來反映每個格子內在不同時間段的密度和壓力。本文壓縮過程運行了832 ps，膨脹過程運行了960 ps。圖4所示為在壓縮過程He原子在系統內軸向的原子數密度分布和壓力分布隨時間的變化；圖5所示為在膨脹過程He原子在系統內軸向的原子數密度分布和壓力分布隨時間的變化。

圖4壓縮過程軸向密度分布和壓力分布隨時間的變化Fig. 4 Axial density distribution and pressure distribution vary with time in the compression process

由圖4可知，壓縮過程中，在0 ps時，微通道外的He原子數目較多，密度較大，壓力較高。微通道內的He原子數目較少，密度較低，壓力較小，此時系統中存在較大壓差。在系統運行過程中，由于壓差的作用，微通道外的高壓He氣會向含有低壓He氣的微通道內流動。此時微通道外的He原子不斷減少，密度降低，導致壓力減小。微通道內側He原子數目增加，密度升高，導致壓力增加。運行至640 ps時，在軸向方向內，微通道內壓力顯著高于微通道外，產生這種現象主要原因是He原子隨著流動方向運動，在慣性作用下，分子會繼續向左運動。在運行至832 ps時，封閉端處的密度和壓力降低，因為此過程中壓力占主導作用，壓差作用使氣體流動方向逆轉。系統內軸向密度分布梯度和軸向壓力分布梯度均隨著時間的變化逐漸減小，最終密度與壓力梯度趨于平緩，最終壓力穩定在0.65 MPa，密度為每個格子內470個粒子。

由圖5可知，在膨脹過程中，密度與壓力分布隨時間變化與壓縮過程相似，但膨脹過程中He氣流動方向與壓縮過程相反，所以密度梯度方向和壓力梯度方向均相反。608 ps時，在軸向方向內，微通道內密度、壓力低于微通道外的密度、壓力；960 ps時，微通道內密度、壓力升高，微通道外密度、壓力降低。

圖5 膨脹過程軸向密度分布和壓力分布隨時間的變化Fig. 5 Axial density distribution and pressure distribution vary with time in the expansion process

3.2 軸向速度分布隨時間的變化

關于壓縮過程和膨脹過程中軸向速度分布隨時間變化的研究與上述密度研究方式相同，也是將不同時間段各個格子內的分子速度進行統計，并計算其平均值，得出每個格子內的平均速度。圖6所示為在壓縮過程和膨脹過程軸向速度分布隨時間的變化。

圖6 壓縮過程和膨脹過程軸向速度分布隨時間的變化Fig. 6 Axial velocity distribution varies with time in the compression process and expansion process

如圖6(a)所示，研究壓縮過程系統內He原子的流動過程，并對軸向速度分布隨時間的變化進行分析(定義速度向右移動為正)。由圖6(a)可知，0 ps時，速度分布總體趨勢在0 m/s處波動，即系統處于靜止過程。在系統中，兩側由于界面固定的原因，速度始終在0附近波動，最大速度則最先(64 ps)發生在微通道出口附近，分別為775 m/s和864 m/s，并向微通道內部逐漸傳遞，這是由于微通道出口附近為高壓與低壓交界處，起始壓差最大，且壓差傳遞需要時間，導致此處的速度也最先達到最大值；速度在向微通道內部傳遞的同時，速度值也在逐漸下降，640 ps時基本回到初始過程，832 ps時開始反向流動，因為此時微通道內壓力高于微通道外，加速度方向改變。

由圖6(b)可知，在膨脹過程中，速度分布隨時間變化與壓縮過程相似，但速度梯度方向相反。64 ps時，軸向最大速度發生在微通道出口處；608 ps時，軸向速度分布基本回到初始位置；960 ps時，速度方向改變，反向流動。

3.3 速度與壓力波的相位關系

圖7所示為壓縮過程時微通道熱端、冷端和氣源端壓力波與速度波相位關系。由圖7可知，熱端(40 nm)壓力波相位超前速度波相位；冷端(300 nm)壓力波相位與速度波相位同相；氣源端壓力波滯后速度波。由此可知微通道內各點壓力波與速度波相位差隨位置而變，在高溫處相位差較大，在低溫處相位差較小或為負值。

圖7 熱端、冷端和氣源端壓力波與速度波相位關系Fig.7 Phase relationship between pressure wave and velocity wave at at hot end, cold end and gas source end

3.4 軸向溫度分布隨時間的變化

本模擬不僅計算每個格子內平均速度，還計算了每個格子內在不同時間的溫度。圖8所示為在壓縮過程和膨脹過程時軸向溫度分布隨時間的變化。

圖8 壓縮過程和膨脹過程軸向溫度分布隨時間的變化Fig. 8 Axial temperature distribution varies with time in the compression process and expansion process

如圖8(a)所示，在0 ps時，系統內部初始軸向溫度分布穩定在300 K。當系統進行至64 ps時，在微通道出口處，首先出現較低的溫度，而在系統兩端，溫度還維持在約300 K，這是因為此時系統運行時間還較短，壓力波還未傳到兩端。隨著壓縮過程的進行，微通道內部左側底端(脈管熱端)溫度越來越高，右側出口端(脈管冷端)附近溫度越來越低，沿著系統軸向溫度梯度越來越大，且最低溫度由微通道出口處逐漸向系統右端移動。這是由于系統左側的氣體不斷被壓縮，吸收了越來越多壓縮功的能量，因而溫度越來越高；而右側的氣體對左側氣體持續做壓縮功，輸出了越來越多的壓縮功能量，因而溫度越來越低。320 ps時，左側熱端溫度達到最大，可高達500 K，之后由于壓縮過程減弱，氣體的軸向導熱逐漸處于優勢，熱端溫度開始下降；640 ps時，右側溫度達到最低。溫度梯度變化開始發生逆轉現象；832 ps時，左側溫度降低，右側溫度上升，溫度梯度明顯減小，這是由于此時系統內部氣體的流動發生了逆轉。

可以發現，微通道外部高壓氣體部分，溫度梯度分布一直處于較平緩的過程，且均處于較低溫度，有利于脈管制冷。這是由于沒有微通道的限制，氣體容易發生混合過程，所以溫度分布較均勻。在微通道內，接近底部封閉端處溫度達到最高，是由于此處氣體速度為零，存在“滯止溫度”。由“滯止溫度”可知，在流場中，速度為零的點，靜壓最大，此時溫度最大。

如圖8(b)所示，在膨脹過程中，溫度分布隨時間的變化與壓縮過程相似，但溫度梯度方向相反。沿軸向溫度分布的最低溫度開始也發生在微通道出口處附近，在288 ps時達到最低，可以降至223 K，此時系統左側溫度有所降低，右端溫度達到最高。隨后最低溫度逐漸向左側微通道內部移動，同時微通道出口處溫度逐漸升高，在608 ps時最低溫度移動到微通道最左側底端，但右端最高溫度有所下降，說明氣體軸向導熱開始占優勢。

3.5 一次循環溫度場疊加分布

圖9所示為微通道在經歷了一次壓縮和膨脹后的疊加氣體軸向溫度分布。壓縮與膨脹時間均取為576 ps，此時刻壓縮與膨脹過程中溫度變化達到最佳值。雖然溫度梯度與單獨壓縮、膨脹過程相比有大幅下降，但仍有80 K的平均溫差，可以在脈管兩端實現溫度差。圖10所示為冷、熱端壓縮過程、膨脹過程、壓縮與膨脹過程疊加后溫度隨時間的變化。疊加溫度是對相同時間內的壓縮和膨脹過程瞬時溫度疊加再減去設定的環境溫度300 K，表示此脈管能夠制冷。由圖10(a)中疊加溫度線可知，熱端(微通道左側，20～40 nm)溫度始終維持在300 K以上，時間積分的平均值為375 K，有利于熱端向環境散熱。由圖10(b)中疊加溫度線可知，冷端(微通道出口處，320～340 nm)溫度大部分時間在300 K以下，按照時間積分平均值計算為244 K，因此，脈管在經歷了一次壓縮、膨脹循環后，冷端溫度低于初始溫度300 K。因此，如果連接回熱器，經過多次循環后，冷端溫度會持續下降。

圖9 膨脹與壓縮過程疊加后瞬時軸向溫度分布Fig.9 Instantaneous axial temperature distribution after superposition of expansion and compression process

圖10 溫度隨時間的變化Fig.10 Temperature varies with time

3.6 脈管微通道的優勢

層流化原件工作原理為：增加內部流通阻力，減緩流速，形成層流。層流化原件是由微小通道組合而成，由于所做的脈管通道管徑較小，達到微細尺度，所以此脈管可以起到層流化原件作用，代替脈管制冷機中層流化原件。通過壓縮與膨脹過程所得的溫度變化結果，在脈管管徑達到微細尺度時，對微通道內進行充放氣，管道內部會出現熱端和冷端，可以制冷，突破了焓流理論對脈管長徑比的要求。

4 結論

本文運用分子動力學模擬的方法，通過對微通道的充氣和放氣過程，研究了基本型脈管制冷機中脈管內氣體的壓縮和膨脹過程，并與相關實驗結論吻合[4]，得出如下結論：

1) 脈管內部沿軸向首先建立起壓力、密度梯度，隨著過程的進行，密度、壓力分布梯度逐漸減小，直至達到平衡，此時壓力穩定在650 kPa，密度為每個格子內470個粒子；然后過程繼續進行，因而出現了微小的逆向梯度。

2) 微通道內氣體軸向速度分布不均勻，速度隨著時間的增加，先增大后減小；膨脹過程和壓縮過程在64 ps時，在微通道出口附近，出現最大速度，分別為775 m/s和864 m/s，隨著反應的進行，最大速度處逐漸向壓力低的方向移動。

3) 微通道內各點壓力波與速度波不同相，相位差隨位置而變。

4) 微通道內軸向溫度梯度隨著壓縮、膨脹過程的進行逐漸增大，在壓縮過程中，微通道內靠近封閉端溫度較高，可高達500 K;遠離封閉端溫度較低，可降至223 K；膨脹過程則相反。微通道外部氣體溫度分布較平緩。一次循環溫度場疊加時，在熱端處，時間積分的平均值為375 K；在冷端處，時間積分的平均值為244 K；有利于熱端向環境散熱和冷端制冷。