潘 帥，祁影霞，車閆瑾，淡一埔，劉雅麗，張 華

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

0 引言

自1964年Gifford等[1]發(fā)明脈管制冷機（PTR）以來，脈管制冷機以體積小、冷端無運動部件、運行壽命長、可靠性高、振動小等優(yōu)異特性，被廣泛用于低溫電子器件、航空航天技術、紅外超導技術、低溫醫(yī)學等領域。脈管制冷機的發(fā)展也從最開始的基本型逐漸發(fā)展到小孔型脈管制冷機（OPTR）[2]，雙向進氣型脈管制冷機（DPDR）[3]以及慣性管加氣庫調相形式的產生[4]。同時對于脈管制冷機內部機理的研究也在不斷深入。

吳明等[5]采用數(shù)值模擬的方法對脈管的長徑比進行了研究，發(fā)現(xiàn)脈管存在最佳長徑比，在此長徑比下制冷量最大。同時最佳長徑比隨著脈管的結構及運行參數(shù)的不同而變化。Zhi等[6]模擬了脈沖管低溫冷卻器在不同冷端相角下的制冷機理，研究了其關鍵部件氣團的熱力學循環(huán)。通過比較冷端換熱器兩側氣團的傳熱特性，揭示了冷端產生冷卻功率的機理。結果表明，冷端相角決定了脈沖管內氣體工質循環(huán)過程中的傳熱方向和大小。劉少帥等[7]建立阻抗匹配模型，從冷指、壓縮機方面優(yōu)化脈管制冷機的性能。同時，建立了30 K溫區(qū)脈管冷指的一維DeltaEC理論模型，優(yōu)化慣性管脈管制冷機冷指與壓縮機耦合運行參數(shù)，分析回熱器、雙段慣性管長度和氣庫體積對冷指聲學阻抗及壓縮機聲功轉化效率的耦合關系。迄今關于脈管內部的仿真研究主要以計算流體力學（CFD）和數(shù)值模擬為主，均采用節(jié)點分析技術對PTR特定穩(wěn)態(tài)進行研究。將PTR作為封閉體進行單元劃分；每個控制單元采用能量守恒，質量守恒和動量守恒方程，用于模擬某個PTR系統(tǒng)的穩(wěn)態(tài)。因此，不可能模擬脈管中氣體的熱力學不對稱交替流動。

對于非平衡過程，分子動力學（MD）模擬是一種行之有效的方法。隨著計算機技術的快速發(fā)展，MD方法被用于計算復雜的系統(tǒng)。Cao等[8]運用非平衡分子動力學模擬研究表面粗糙度對亞微米鉑通道中氣態(tài)氬氣滑流的影響，得出速度滑移的邊界條件不僅取決于克努森數(shù)，還取決于表面粗糙度。摩擦系數(shù)不僅隨著克努森數(shù)的減少而增加，而且隨著表面粗糙度的增加而增加。Kazemi等[9]對有機納米通道中氣體流動進行了非平衡分子動力學模擬研究。研究表明，氣體分子通過有機納米級管道的傳輸方式是擴散。因此建立一個可以預測回收非常規(guī)流體的實際模型，在有機微納米孔中的傳輸方程應該用擴散傳輸方程代替。本文采用分子動力學模擬，建立了不同幾何條件下的微通道脈管模型，獲得脈管內分子運動的微觀細節(jié)，研究幾何參數(shù)變化對脈管制冷性能的影響。

1 分子動力學模擬方法

分子動力學模擬遵循經典牛頓力學定律。在研究含有大量原子的獨立系統(tǒng)時，考慮系統(tǒng)總能量為系統(tǒng)中所有顆粒的動能（EK）與勢能（U）的總和，勢能分為原子間非鍵合范德瓦爾斯作用（VDW）與分子內勢能（INT）如式（1）所示。

模型中所有粒子均為獨立原子，所以忽略分子內勢能。VDW近似為原子對之間VDW的加成。

式中：rij為i、j兩原子之間的距離。

基于模型的精度要求，選擇采用混合勢力場。Fe-Fe之間采用嵌入原子方法EAM/FS（Embedded Atom Method）描述。該多體式力場的描述方法是把晶體的總勢能分為原子間作用力的短程對勢和原子鑲嵌在電子云中的嵌入能[10]。EAM/FS勢函數(shù)常用于描述金屬離子之間的相互作用力，如式（3）所示：

式中：Ei為原子間的總勢能；Fα為原子嵌入電子云的嵌入能；ραβ為電子密度函數(shù)；φαβ為原子核間的短程對勢；rij為i、j兩原子之間的距離。

文中采用的流體原子為氦原子，壁面原子為鐵原子。固體、氣體分子之間存在相互作用，當兩個粒子間距非常小時產生排斥，當粒子相距一定距離時又會產生吸引力。He-He、He-Fe之間均采用Lennard-Jone勢能，其勢能表達式如式（4）所示。

式中：ε為勢能參數(shù)，其值為0.000 607 098；σ為勢能參數(shù)，其值為2.103；rij為i、j兩原子之間的距離。

運用牛頓經典力學獲得各原子的受力及加速度：

模擬運用Leap-frog方法進行積分預測獲得不同時間i原子的運動位移、速度等參數(shù)：

2 模型建立及模擬方法

建立矩形截面直通道脈沖管模型，通道外為高壓He、通道內為低壓He，模擬脈管的充氣過程；相反，將通道外設置為低壓He，通道內設置為高壓氣體，模擬脈管的放氣過程。模型建立過程：首先建立高低壓He原子為流體介質的穩(wěn)態(tài)模型，同時引入周期性邊界條件，之后分別讓高低壓穩(wěn)態(tài)模型在NTV（粒子數(shù)N、體積V、溫度T）正則系綜下運行，分別使兩個穩(wěn)態(tài)模型內部氣體分子均勻混合；其次分別建立軸向Fe壁面和縱向Fe壁面?；谏鲜龇椒?，建立了基本型脈管的充放氣模型。對于充氣模型，初始通道外高壓1 200 kPa，管內低壓100 kPa?；拘兔}管模型從左向右分為兩部分，分別是微通道脈管及內部低壓He和通道外高壓He。圖1為基本脈管充氣模型是（1B-2L）初始時刻的原子分布，其中1B代表模型寬573.26 nm，2L代表模型長1 145.41 nm。具體參數(shù)如表1所列。

表1 各脈管充氣模型尺寸信息Tab.1

圖1 基本型脈管充氣模型及初始原子分布Fig.1 The basic inflation model and the initial atom distribution

放氣模型是將高低壓He空間對調位置，與充氣模型類似。對系統(tǒng)在NTV正則系綜下運行溫度的標定，將初始溫度設定為300 K，僅對通道內的He氣體粒子進行NVE（原子數(shù)N、氣體體積V、氣體分子總能量E）計算，時間步長選用0.4 fs，每50 000步輸出一次結果，運行至通道內外壓力平衡。

3 模擬結果與分析

3.1 模型軸向溫度、壓力速度分布

基于建立的基本型脈管1B-2L充氣模型，在軸向方向，從左到右進行網格劃分，然后在不同時段內對每一個網格內的He原子壓力、速度和個數(shù)進行統(tǒng)計平均得出模型的軸向速度、壓力以及溫度隨時間的變化分布，如圖2～4所示。將坐標原點定義到模型最左端，從左向右為正方向。

圖2 基本型脈管充氣過程軸向壓力分布曲線Fig.2 Axial pressure distribution during inflation of basic type pulse tube

圖3 基本型脈管充氣過程軸向平均速度分布曲線Fig.3 Axial average velocity distribution of basic vessel during inflation

圖4 基本型脈管充氣過程軸向溫度分布曲線Fig.4 Axial temperature distribution of basic pulse tube during inflation

圖2為基本型脈管的充氣過程軸向壓力分布，初始管外壓力為1 200 kPa，管內壓力為100 kPa，隨時間推移管外壓力逐漸減小，脈管內壓力逐漸增大，在1 200 ps時脈管內壓力高于脈管外壓力。這是由于內外的壓力差，使得通道外He原子流向內部，通道內He原子增多，壓力逐漸增大。通道外情況則相反。當內外壓差為0時，由于慣性作用，He原子繼續(xù)向通道內流動，出現(xiàn)通道內壓力高于通道外的情況。1 200 ps之后，管內壓力下降，管外壓力上升，直至內外壓力再平衡（平衡壓力650 kPa）。圖3為基本型脈管充氣過程軸向平均速度分布，初始狀態(tài)（管外壓力1 200 kPa，管內壓力100 kPa）時He原子在壓差的作用下速度迅速提升，在200 ps出現(xiàn)最大值（434 m/s），之后峰值逐漸減小同時向管內移動。這是因為初始狀態(tài)通道進口截面處（600 nm）壓差最大速度提升最為迅速，同時速度的變化較壓力變化會有一定的延遲，所以速度最大值是在脈管進口靠左（500 nm）處而不是在進口截面（600 nm）處，隨著壓力向管內傳遞，內外壓差逐漸減小，速度在向管內傳遞的同時也在減小。1 400 ps出現(xiàn)速度負值，是因為此時通道內部壓力梯度反向，使得一部分He原子反向流動到通道外直至內外壓力平衡。圖4為基本型脈管充氣過程軸向溫度分布，初始時刻設定系統(tǒng)溫度300 K，從初始時刻到200 ps，在脈管進口截面附近，通道外溫度下降，通道內溫度上升，而模型兩端溫度還基本維持在300 K。之后系統(tǒng)最高溫度逐漸左移，最低溫逐漸右移。從原子做功角度分析，運行初期通道外He原子在壓力差作用下向通道內運動，通道外原子對臨近的通道內原子做功，所以通道內溫度上升，通道外溫度下降，而兩端由于功的作用還沒有傳遞過去，所以還維持在初始溫度300 K。隨著過程進行，通道外原子持續(xù)對通道內做功，功的作用持續(xù)向兩端傳遞，使兩端溫差逐漸擴大，在1 200 ps時熱端溫度達到最高498 K，在遠離封閉端處溫度降低至217 K。之后熱端溫度下降冷端溫度上升，這是由于He原子流動轉向，致使軸向氣體混流加劇，導致通道內外的溫度差減小。微通道脈管的放氣過程與充氣過程類似，只是參數(shù)的變化方向相反。

3.2 微通道長度變化影響

圖5和圖6表示模型中脈管寬度保持不變，長度改變后，軸向溫度分布的變化?？梢钥闯?，對于微通道脈管模型，隨著長度增大，冷端溫度逐漸降低而熱端溫度逐步上升，同時發(fā)現(xiàn)熱端最高溫度隨著長度增大而增大的幅度減小。冷端最低溫度從289.5 K（1B-1L）逐漸降低到265.2 K（1B-3L），熱端的最高溫度從461.8 K（1B-1L）上升到505.8 K（1B-3L），需要說明的是冷端為模型氣體進口，熱端為模型脈管末端。對比冷熱端溫度，取逐時溫差，獲得冷熱端溫度差的逐時變化如圖7所示。

圖5 模型中脈管長度改變后冷端溫度分布曲線Fig.5 Temperature distribution of cold end in the model with the change of length of pulse tube

圖6 模型中脈管長度改變后熱端溫度分布曲線Fig.6 Temperature distribution at the hot end in the model with the change of length of the vessel

圖7 脈管長度改變后冷熱端溫度差的逐時變化曲線Fig.7 Time by time change of temperature difference between hot and cold ends after pulse length change

由圖7可以發(fā)現(xiàn)，隨著長度增加，最大溫差增大，但增大的幅度在減小，表明當壓比一定時（12∶1），存在一最佳長徑比。由于脈管長度的增長，微通道脈管內部氣體分子在壓差作用下，運行時間增長，當脈管內外壓差變?yōu)榱愫?，由慣性力主導的脈管管外分子對脈管封閉端做功的持續(xù)位移也越長，做功增多，所以脈管冷端溫度降低，熱端溫度增高。然而這個過程不是無限的，不能持續(xù)地增加脈管長度，脈管內氣體分子也無法在慣性作用下持續(xù)對脈管封閉端做功。所以適當增加微通道脈管長度可以提高制冷性能，但也存在一個極值。

3.3 微通道寬度變化的影響

模型中脈管長度保持不變，寬度改變后，軸向溫度分布變化如圖8～10所示?？梢钥闯?，隨著模型脈管寬度的增大，冷端溫度逐漸降低而熱端溫度逐步上升，同時可以發(fā)現(xiàn)，熱端最高溫度隨著模型寬度增大而增大的幅度在減小，冷端的最低溫度從265.3 K（1B-3L）逐漸降低到244.2 K（3B-3L），熱端的最高溫度從505.9 K（1B-1L）上升到550.6 K（3B-3L）。對比冷熱端溫度，取逐時溫差，獲得溫差逐時變化如圖10所示?？芍S著寬度增大，冷熱端最大溫差增大但增大的幅度在減小，表明對于當壓力比一定時（12∶1），存在一最佳長徑比。由于保持長度不變增加微通道脈管寬度，同時也增加了脈管中分子數(shù)量，相應的分子對微通道脈管熱端做功越多，在脈管冷端產生的冷量也越多，脈管制冷機性能提高。然而這個促進作用也有限度，由于脈管寬度的增加，也增大了脈管內部分子運動的混亂程度，增加了微通道脈管冷熱端的熱傳導，抑制了微通道脈管制冷機性能的提高。所以增加微通道寬度對脈管制冷機的制冷性能有促進作用，但也存在一個最佳值。

圖8 模型中脈管寬度改變后冷端溫度分布曲線Fig.8 Temperature distribution of cold end in the model with the change of pulse width

圖9 模型中脈管寬度改變后熱端溫度分布Fig.9 Temperature distribution at the hot end in the model with the change of the width of the vessel

圖10 脈管寬度改變后冷熱端溫度差的逐時變化曲線Fig.10 Time by time change of temperature difference between hot and cold ends after pulse width change

3.4 脈管體積變化影響

圖11 、圖12為同一長徑比下改變脈管體積對微通道脈管制冷機冷熱端溫度的影響，即保持長徑比不變，增大脈管體積（同比例放大）。

圖11 同一長徑比下改變脈管體積對微通道脈管制冷機冷端溫度的影響Fig.11 Effect of changing the vessel volume on the cold end temperature of microchannel tube cooler at the same aspect ration

從圖11、圖12中可以看出，隨著脈管體積的增大，脈管冷端溫度逐步降低，熱端溫度逐步升高，同時升高與降低的幅度在減小，說明當壓力比（12∶1）和長徑比一定時，脈管體積存在一最佳值使得冷熱端溫差最大。綜合微通道脈管長度變化與寬度變化的共同影響，在保持微通道脈管長徑比不變的條件下適當提高脈管體積可以提高其制冷性能。這表明在微通道脈管的設計過程中適當?shù)胤糯竺}管體積可以提高脈管制冷機的工作性能。

圖12 同一長徑比下改變脈管體積對微通道脈管制冷機熱端溫度影響Fig.12 Effect of changing the vessel volume on the hot end temperature of microchannel tube cooler at the same aspect ration

4 結論

本文采用分子動力學模擬方法研究了脈管在壓差作用下He充氣過程的原子運動過程，得到的結論為：

（1）脈管充氣過程中原子軸向速度先增大后減小，速度最大值的位置在充氣過程中向脈管內移動；

（2）在保持脈管寬度不變的情況下，增加微通道脈管長度，或者保持長度不變增加寬度均可使脈管熱端最高溫度上升，冷端最低溫度下降，但變化的幅度在減小，說明任何寬度和長度下，都存在一個最佳長徑比；

（3）在保持微通道脈管長徑比不變，同比例放大的情況下，發(fā)現(xiàn)脈管冷端溫度下降，熱端溫度上升，但變化幅度在減小，說明脈管體積存在一個最佳值。當壓力比為12∶1，長徑比為2.86時，該最佳值至少大于0.000 3μm3。