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

（上海理工大學(xué)能源與動(dòng)力工程學(xué)院，上海 200093）

0 引言

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

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

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

1 分子動(dòng)力學(xué)模擬方法

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

模型中所有粒子均為獨(dú)立原子，所以忽略分子內(nèi)勢(shì)能。VDW近似為原子對(duì)之間VDW的加成。

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

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

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

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

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

運(yùn)用牛頓經(jīng)典力學(xué)獲得各原子的受力及加速度：

模擬運(yùn)用Leap-frog方法進(jìn)行積分預(yù)測(cè)獲得不同時(shí)間i原子的運(yùn)動(dòng)位移、速度等參數(shù)：

2 模型建立及模擬方法

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

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

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

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

3 模擬結(jié)果與分析

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

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

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

3.2 微通道長(zhǎng)度變化影響

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

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

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

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

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

3.3 微通道寬度變化的影響

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

圖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 脈管寬度改變后冷熱端溫度差的逐時(shí)變化曲線Fig.10 Time by time change of temperature difference between hot and cold ends after pulse width change

3.4 脈管體積變化影響

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

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

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

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

4 結(jié)論

本文采用分子動(dòng)力學(xué)模擬方法研究了脈管在壓差作用下He充氣過程的原子運(yùn)動(dòng)過程，得到的結(jié)論為：

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

（2）在保持脈管寬度不變的情況下，增加微通道脈管長(zhǎng)度，或者保持長(zhǎng)度不變?cè)黾訉挾染墒姑}管熱端最高溫度上升，冷端最低溫度下降，但變化的幅度在減小，說明任何寬度和長(zhǎng)度下，都存在一個(gè)最佳長(zhǎng)徑比；

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