低溫可視化力學裝置傳熱分析與計算研究

吳姍姍,王玨,張恒成,劉輝明,黃傳軍,黃榮進,李來風,周遠

(1. 中國科學院理化技術研究所航天低溫推進劑技術國家重點實驗室,100190,北京;2. 中國科學院大學,100049,北京; 3. 航天科工空間工程發(fā)展有限公司,100854,北京)

材料的低溫力學性能與實際工程應用息息相關,ITER計劃通過全超導Tokamak產(chǎn)生的強磁場來實現(xiàn)大規(guī)模的核聚變反應[1-2]。然而,超導線纜產(chǎn)生強磁場時在電磁力的作用下容易發(fā)生變形甚至破壞,從而影響超導磁體正常運行[3-4]。超導磁體在極端溫度(4.2 K)下運行,此過程中材料的力學性能發(fā)生改變[5]。因此,超導磁體材料以及線圈外部的保護夾套在工作環(huán)境下的力學性能是保證超導磁體安全運行的重要因素[6]。此外,材料的低溫性能在航空航天液體燃料儲存系統(tǒng)、材料改性等應用中也引起了眾多學者的關注。材料的力學性能與材料的成分、組織結構、結構缺陷以及材料所處的環(huán)境溫度有關,特別是在低溫下應用時,其力學性能與室溫下截然不同[11-13]。以金屬材料為例,面心立方的金屬在低溫下仍具有較好的韌性與強度,體心立方的金屬在低溫下易發(fā)生低溫脆性,而316LN不銹鋼材料在低溫環(huán)境中則會發(fā)生應力誘導馬氏體相變,因此需要了解材料隨溫度變化的綜合力學行為[14-15]。

在材料的變形測試中,位移和應變是保證測量精度的重要數(shù)據(jù),而機械引伸計測量和光學測量是應變測量中常用的方法。其中機械引伸計受測量環(huán)境(溫度、震動)的影響大,且測量范圍有限,無法測量位移場及應變場,在使用時受到諸多限制[16]。數(shù)字圖像相關(DIC)光學測量方式是通過對比試樣變形前后的數(shù)字圖像來獲取力學性能的測試方法,可以對試樣表面進行高精度全場動態(tài)測量[17-22],本文中的低溫力學性能測試裝置采用后者對試樣進行測量。

為了將3D-DIC技術應用到材料的低溫變形測試中,需要為試樣提供便于采集數(shù)字圖像的低溫環(huán)境。目前,常用的低溫冷卻方式包括低溫液體冷卻和制冷機冷卻兩種方法。低溫液體冷卻包括液氦(4.2 K)、液氫(20.3 K)、液氮(77.4 K)和液化天然氣(111 K)等,這種冷卻方式是通過低溫液體的汽化潛熱和氣體顯熱將物體進行冷卻[23-26]。然而,低溫液體會使采集到的數(shù)字圖像出現(xiàn)扭曲與變形,低溫氦氣的熱對流也會引起折射率的變化,從而使相機拍到的數(shù)字圖像模糊,影響DIC方法測量的準確性。制冷機冷卻方式通過有效的傳熱方式將物體溫度在一定溫度范圍內連續(xù)、可調、可控地降低,且無需采用其他低溫工質[27-29],本文采用制冷機冷卻的方式來提供連續(xù)可調可控的低溫環(huán)境。由于低溫裝置需要配置光學窗口,其發(fā)射率較大,無疑增加了被測樣品與外界的輻射換熱量,且制冷機在低溫下提供的冷量有限,為了在更低溫度下測得材料的力學性能,需可靠估計低溫裝置內的漏熱量,選取直觀的傳熱分析方法尤為重要。

本文首先介紹了4.2～300 K溫區(qū)內材料力學性能測試的可視化低溫裝置的幾何結構,并對裝置進行了傳熱過程分析;其次,對此研究對象中各部件之間的傳熱關系進行了介紹,并分別采用了COMSOL、Sage軟件對試樣進行傳熱計算;最后,分析對比了兩種計算方法的溫度分布,并基于Sage研究了試樣尺寸、試樣發(fā)射率、窗口尺寸、窗口表面發(fā)射率對試樣溫度分布的影響,為低溫裝置中傳熱分析與計算提供理論參考。

1 研究對象及傳熱過程分析

本文中低溫DIC力學性能測試系統(tǒng)的結構具有不同的換熱方式,需要對其進行理論傳熱分析。

1.1 研究對象

本裝置中使用的3D-DIC光學成像系統(tǒng)主要由兩臺高速攝像機和普通白光組成,其中DIC測量系統(tǒng)主要包括光學成像系統(tǒng)、光電轉換傳感器和數(shù)字圖像處理系統(tǒng),光學成像系統(tǒng)將試樣表面的散斑圖像通過攝像機讀取至圖像卡中,通過計算機中的數(shù)字圖像處理系統(tǒng)將圖像進行相關計算轉換。低溫-DIC力學性能測試裝置結構及其傳熱關系示意圖如圖1所示,整個裝置為偏置式,真空罩和輻射屏上分別開設了兩個窗口,一個為50 mm×230 mm的矩形石英玻璃窗口,另一個為Ф50 mm的圓形石英玻璃窗口。真空罩內真空度可達1.3×10-5Pa。試樣測試前需均勻噴涂一層啞光漆便于圖像識別,兩臺相機在豎直方向上通過兩層矩形窗口呈15°～45°對試樣表面進行圖像采集,白光光源通過兩層圓形窗口在試樣表面形成均勻光場。 低溫力學測試裝置內主要部件如圖1(b)所示,真空罩、支撐管與夾具均為不銹鋼材質,拉伸桿為鈦合金材質,輻射屏、熱橋以及制冷機冷頭為銅質,矩形窗口和圓形窗口均為石英玻璃,熱橋與連接端面均勻涂覆Apiezon N脂并壓接。真空罩內部在實驗過程中始終保持真空,兩臺住友RDK-415D制冷機作為冷源安裝在真空罩法蘭上,輻射屏通過制冷機一級冷頭連接并冷卻,制冷機二級冷頭通過銅質柔性熱橋連接至力學拉伸部件使其冷卻。力學拉伸組件主要由拉伸桿、支撐管、拉伸夾具組成,拉伸桿直接連接真空罩法蘭(室溫端,300 K)和低溫夾具(低溫端,4.2 K),支撐管分別連接真空罩法蘭、輻射屏以及低溫夾具,夾具通過熱橋與冷頭相連,且直接與試樣接觸。

(a)低溫-DIC裝置外部示意圖

(b)裝置內部幾何結構及傳熱關系示意圖

1.2 傳熱過程分析

熱量通過拉伸桿和支撐管從高溫端向低溫端傳導,材料的導熱換熱量為

(1)

式中:A為橫截面積;L為導熱長度;ks為熱導率,是固體材料溫度依賴的特性;Qcond為凈熱流。

除熱傳導外,部件間的輻射傳熱亦不可忽略。裝置內存在溫差的部件間皆存在輻射換熱,如真空罩與輻射屏間的輻射換熱Qrad1、輻射屏窗口與試樣間的輻射換熱Qrad2-1、輻射屏窗口與拉伸組件表面的輻射換熱Qrad2-2、輻射屏對制冷機組件的輻射換熱Qrad3、輻射屏與拉伸組件表面的輻射換熱Qrad4以及輻射屏與試樣間的輻射換熱Qrad5等。

兩表面間的輻射換熱量為

Qrad=A1X1,2J1-A2X2,1J2

(2)

(3)

式中:Qrad為兩表面間的輻射換熱量;A1、A2分別為兩輻射表面的輻射面積;X1,2、X2,1分別為從表面1到表面2的角系數(shù)和從表面2到表面1的角系數(shù);Ji為i表面的有效輻射力;εi為i表面的表面發(fā)射率;C0為黑體輻射系數(shù),其值為5.67 W/(m2·K4);Ti為i表面的溫度。

在實驗過程中,真空罩內恒為高真空(10-4Pa)環(huán)境,因此忽略對流換熱。若Th、Tl、Ts在此裝置內分別表示室溫、低溫以及輻射屏溫度,根據(jù)以上裝置結構和傳熱分析結果,與試樣有關的導熱和輻射換熱關系可表示為熱阻關系圖,如圖2所示。熱阻R可表示為

(4)

式中:ΔT為物體兩端的溫差;Q為輻射換熱量或導熱換熱量。

(a)低溫端Tl向室溫Th及輻射屏Ts間的傳熱

(b)輻射屏Ts向室溫Th的傳熱

2 模型簡化及計算方法

為了探索更準確簡潔的傳熱計算方法,針對上節(jié)中的幾何模型和傳熱關系,首先使用COMSOL三維求解方法計算了各個傳熱結構的傳熱量,然而此方法所需的計算模型復雜,且計算時間較長。因此,使用Sage一維求解方法分別開發(fā)出導熱模型和輻射模型,此方法的計算成本顯著降低。

2.1 COMSOL三維求解

COMSOL Multiphysics可將單一類型物理場擴展到耦合的物理場模型,且可以同時進行計算。對整個低溫力學拉伸裝置建立涵蓋兩臺制冷機、真空罩、輻射屏、石英窗口、被測試樣、拉伸夾具、拉伸桿、支撐管、熱橋等結構的三維輻射場與導熱場的耦合求解域物理模型。根據(jù)以上結構分析與傳熱分析結果,在COMSOL軟件仿真中做出如下假設:①除與熱橋接觸的界面外,其余所有界面間為理想接觸;②所有參與導熱的結構溫度呈線性分布;③玻璃窗口與真空罩、輻射屏的溫度始終保持一致。

首先建立幾何模型,為節(jié)約軟件計算成本,三維耦合求解時只關注輻射屏以內的傳熱情況。用于COMSOL軟件的三維簡化幾何模型如圖3所示,輻射屏內部為制冷機、熱橋、夾具、試樣、拉伸桿、支撐管等組件,各個配合結構間均簡化為簡單的面接觸,忽略法蘭結構。本仿真中采用固體傳熱與表面對表面輻射耦合的多物理場,以穩(wěn)態(tài)求解方式計算溫度場的分布情況。

在穩(wěn)態(tài)傳熱場中,定義在固體區(qū)域的溫度方程為傅里葉定律的微分形式,即

(5)

(6)

式中:ρ為密度;cp為比定壓熱容;u為速度矢量;T為絕對溫度;q為熱傳導通量;qr為輻射熱流;Q為外加熱源;k為導熱系數(shù)。

(a) 輻射屏

(b)輻射屏內部組件

輻射屏和二級冷頭處的溫度初始值分別為50、4.2 K。邊界條件設置如下:考慮實際裝置中支撐管上端面、內層石英玻璃、制冷機一級冷頭均與輻射屏緊密接觸,將其溫度統(tǒng)一設置為50 K;二級冷頭及其接觸表面溫度為4.2 K;制冷機冷頭與熱橋界面、夾具與熱橋界面的熱接觸均設為等效熱阻。計算中需要的材料的導熱系數(shù)、密度、恒壓熱容等相關特性數(shù)據(jù)均來源于COMSOL軟件中內置材料庫,具體參數(shù)如表1所示。

初始溫度為軟件計算的多場耦合結果。邊界條件分別設置如下:將試樣表面、制冷機冷頭表面、回熱器表面、拉伸結構表面、玻璃表面以及輻射屏內表面分別設置為漫反射表面,均勻地向各個方向反射輻射強度,將輻射屬性的波長相關性簡化為常數(shù),各表面的表面輻射率如表2所示。由于輻射屏內為真空環(huán)境,環(huán)境輻射率和散射輻照度均定義為0。輻射傳熱的計算方法為Hemicube(半立方體法),模型為穩(wěn)態(tài)研究。通過表面對表面輻射傳熱耦合接口建立了傳熱-表面對表面輻射模型,分析了固體傳熱和輻射換熱對溫度場的影響。由于計算網(wǎng)格較密時,軟件所需內存超出計算機內存能力,在滿足網(wǎng)格獨立性的基礎上,本計算中將網(wǎng)格劃分為24 264個。

表1 COMSOL模型中材料的導熱參數(shù)

表2 表面對表面輻射場中各個表面的輻射率[21-22]

2.2 Sage一維求解

對于輻射換熱,Sage中兩表面之間的輻射通過將它們連接至自帶的輻射交換模型組件來實現(xiàn)建模,每個輻射面均以這種方式連接至其他輻射面并與之交換輻射。由于實際裝置的模型的復雜度較高,因此對于低溫力學拉伸裝置中的局部結構特征,需要將其近似為Sage中可提供的模塊。本裝置結構轉化為Sage模型的簡化如下:裝置中偏置的試樣及拉伸組件在Sage中將采用與真空罩同軸的方式表達;在輻射換熱計算中將拉伸組件統(tǒng)一視為一個整體,由于其形狀復雜,簡化為具有固定輻射面積的表面;初步假設室溫端溫度為300 K,試樣、夾具、制冷機二級冷頭等低溫部件溫度為4.2 K,輻射屏及其連接界面溫度為50 K。

首先采用Sage計算導熱傳熱量,根據(jù)上述裝置的傳熱關系分析可知,參與導熱的部件包括拉伸桿和支撐管,其材料、尺寸及傳熱參數(shù)如表3所示,Sage模型如圖4所示。該模型為簡單的固體傳導,主要包括熱端、冷端和導熱體,導熱體設置為一維導熱。

表3 導熱部件尺寸及溫度參數(shù)

圖4 低溫拉伸裝置中導熱部件的Sage模型圖Fig.4 Sage model of thermal conductive components in cryogenic mechanical testing device

各表面之間輻射換熱的Sage模型均由熱端、冷端和視圖配置組件組成,視圖配置組件根據(jù)兩個相鄰輻射表面的幾何形狀以及它們在空間中的相對方向在Sage自帶的組件中選取,輻射換熱的Sage模型如圖5所示。視圖配置組件匯總在子模塊中。以輻射屏窗口與試樣表面的輻射傳熱為例,圖5(b)展示了子模塊中輻射換熱關系的基本結構,可以看出矩形窗口表面和圓形窗口表面為兩個連接至熱端點熱源的輻射表面,試樣表面連接至冷端點熱源。假設試樣表面與兩窗口表面為平行放置,選取Sage中的平行表面組件,通過設置輻射流將輻射表面與組件連接。通過輸入各表面計算參數(shù),可計算出自熱端表面至冷端表面的輻射換熱量。計算過程中各表面的計算參數(shù)如表4所示。采用Sage計算時,共劃分為420個節(jié)點。

(a)輻射屏窗口與試樣的Sage輻射模型

(b)輻射模型子模塊的基本結構

表4 Sage輻射換熱模型中采用的參數(shù)

根據(jù)上述計算的低溫部件與外部環(huán)境的總換熱量計算每個部件的熱阻,還原部件自身因熱阻產(chǎn)生的溫差。

Apiezon N脂兩側的界面溫差、導熱部件兩端溫差為

(7)

(8)

式中:φ為換熱量;A1-2為界面換熱面積;A1為導熱截面積;hAp為界面換熱系數(shù);λ為導熱系數(shù)。

裝置中的熱阻如圖6所示,圖中C為冷頭,tb1為熱橋冷端,tb2為熱橋熱端,jig1為夾具冷端,jig2為夾具熱端,Sam1為試樣端部,Sam2為試樣中部。熱橋端部與其他部件的接觸面之間界面換熱系數(shù)hAp取1 835 W/(m2·K)[12]。部件的導熱參數(shù)取平均溫度對應的值,其結構參數(shù)與三維仿真中模型的取值相同。

圖6 低溫力學拉伸裝置內熱阻簡化圖Fig.6 Simplified diagram of thermal resistance in the device

3 計算結果及討論

對比分析了三維和一維求解方法對試樣溫度的計算結果,在保證計算精度的前提下,選擇計算成本相對較低的一維求解方法分別計算了試樣尺寸、試樣表面發(fā)射率、窗口表面發(fā)射率和窗口尺寸對試樣溫度影響,為后續(xù)實驗提供了參考。

3.1 三維求解方法與一維求解方法的計算結果對比

三維物理場計算結果如圖7所示,導熱部件的溫度均呈線性分布,試樣的兩端溫度為4.219 4 K,試樣中部溫度(最高溫度)為9.676 5 K。

一維計算中,導熱和輻射換熱量如表5所示。從圖1(b)可以看出,熱阻計算中從試樣至冷頭的換熱量為Qrad2-1+Qrad5,從夾具至冷頭的換熱量為Qrad2-1+Qrad5Qrad2-2+Qrad4。各部件或接觸面間的溫差如表6所示,試樣自身的最低溫度(即兩端溫度)、最高溫度(即中部溫度)分別為TSam1=4.211 5 K、TSam2=9.612 9 K。對比三維計算與一維計算的結果,可以看出兩種計算結果幾乎一致,認為在三維傳熱模型中可以將其解耦為各部件的一維傳熱單獨計算,節(jié)約了計算時間,降低了操作難度。此外,由于一維計算中各部件被拆分為單獨的模型,因此各部件的傳熱路徑更加清晰,且傳熱參數(shù)更易提取。

(a)各部件的溫度分布情況

(b)試樣溫度分布

表5 Sage一維方法計算的制冷機二級冷頭負載

表6 各溫度節(jié)點間的溫差

3.2 試樣尺寸對其溫度分布的影響

為進一步驗證兩種計算方法的一致性,保持其他結構參數(shù)不變,比較了3種不同尺寸的試樣溫度分布情況,其中試樣1、2、3的尺寸分別取25 mm×230 mm×2 mm、12.5 mm×100 mm×2.5 mm和3 mm×30 mm×2.5 mm。試樣尺寸對其溫度的影響如圖8所示,分別使用COMSOL、Sage計算的試樣最低溫度、最高溫度,可以看出兩種計算方法的結果基本吻合。3種尺寸的試樣最低溫度始終保持在4.20～4.22 K之間;最高溫度卻因尺寸不同有較大差別,試樣1的最高溫度高達約9.7 K,而試樣3的最高溫度最高約為4.4 K,說明試樣尺寸是影響試樣溫度是否均勻的重要因素,且試樣尺寸越小,其溫度分布越趨于均勻。

根據(jù)圖8中兩種方法的結果,將其存在差異的原因總結如下:COMSOL通過三維幾何模型可以將計算模型與實際裝置中各部件的關系完好的呈現(xiàn)與匹配;而Sage的默認模型中只能根據(jù)各部件的關系選擇較為合適的固有模型,在此基礎上設置各參數(shù)的值。然而此裝置中的很多部件與輻射屏為偏心放置,無法在Sage中匹配到完全適用的模型,故只能做近似處理與計算。COMSOL中將所有參與輻射換熱的表面設置為漫反射表面,最終可獲得綜合輻射換熱的結果;Sage中只考慮了各部件與輻射屏或光學窗口之間單獨的輻射關系,而未將輻射屏內各部件之間的輻射關系考慮在內。

圖8 試樣尺寸對其溫度的影響Fig.8 Effect of specimen size on specimen temperature

3.3 試樣表面發(fā)射率對其溫度分布的影響

由一維和三維計算方法的對比可知,兩種方法的計算結果幾乎一致,而Sage的計算時間和操作難度遠小于COMSOL。因此,使用Sage研究了上述3種尺寸的試樣在表面發(fā)射率分別為0.01、0.25、0.5、0.75以及0.95時試樣的溫度分布情況。試樣尺寸與試樣表面發(fā)射率對其溫度的影響如圖9所示,以試樣1為例,隨著試樣發(fā)射率由0.95減小至0.01,試樣最低溫度由4.213 K降至4.202 K,降低幅度較小;而試樣最高溫度降低幅度較大,由9.799 K降至4.447 K。隨著試樣尺寸減小,表面發(fā)射率對溫度的影響顯著減小,其溫度分布逐漸趨于均勻。試樣3在發(fā)射率為0.01時,端部、中部溫度分別為4.201、4.206 K,由此可見,在低溫力學拉伸裝置的實際使用中,減小試樣尺寸和試樣的表面發(fā)射率是使試樣溫度均勻分布的有效方法。

圖9 窗口發(fā)射率為0.94時試樣尺寸與試樣表面發(fā)射率對其溫度的影響Fig.9 The effect of sample size and specimen surface emissivity on its temperature for a window emissivity

3.4 窗口的表面發(fā)射率及其尺寸對試樣溫度分布的影響

使用Sage計算了窗口參數(shù)對3種試樣溫度分布的影響,其中試樣的表面發(fā)射率均為0.95。玻璃窗口的表面發(fā)射率對試樣溫度的影響如圖10所示,由圖10可知:隨著表面發(fā)射率由0.94降至0.01,試樣1的自身溫差由5.586 K降至3.195 K,尤其在窗口表面發(fā)射率從0.25降至0.01時,試樣自身溫差降幅顯著;對于小尺寸試樣3,試樣自身溫差始終保持最低,當窗口發(fā)射率為0.94時,試樣自身最大溫差僅為0.295 K。由此可見,降低窗口的表面發(fā)射率對大尺寸試樣具有較為明顯的效果,但小試樣的溫度分布幾乎不受窗口表面發(fā)射率的影響。

圖10 試樣發(fā)射率為0.95時窗口表面發(fā)射率對試樣溫度的影響Fig.10 Effect of window surface emissivity on sample temperature at sample emissivity of 0.95

窗口表面發(fā)射率為0.94、試樣表面發(fā)射率為0.95時,窗口的有效輻射面積對3種試樣溫度的影響如圖11所示。由圖11可知,隨著窗口的有效輻射面積增大,試樣溫度近似呈線性增長。隨著窗口有效輻射面積由初始面積減小至初始面積的25%,試樣1的自身溫差由5.586 K減小至3.796 K,試樣3的自身溫差由0.294 K減小至0.264 K。因此,窗口的有效輻射面積對大尺寸試樣的溫度分布影響較大,對小尺寸試樣幾乎沒有影響。

圖11 試樣發(fā)射率0.95、窗口發(fā)射率0.94時窗口有效輻射面積對試樣溫度的影響Fig.11 Effect of window effective radiation area on sample temperature for sample emissivity of 0.95 and window emissivity of 0.94

4 結 論

本文基于低溫力學拉伸裝置,分別使用COMSOL三維計算方法和Sage一維計算方法對裝置中部件的傳熱情況進行了詳細的數(shù)值研究,在此基礎上對不同尺寸的試樣溫度進行了對比,最后研究了試樣參數(shù)與窗口參數(shù)對試樣溫度分布的影響。得到如下主要結論。

(1)在相同的簡化條件下,Sage一維計算方法與COMSOL三維仿真方法在傳熱計算與試樣的溫度場分析中得到的最低溫度僅相差0.007 9 K,結果基本吻合。使用Sage軟件建立的導熱和輻射模型有助于對裝置內的傳熱路徑和熱阻進行更為直觀的分析。在傳熱場分析中,可以利用各部件間的空間關系選取Sage模塊來代替復雜的三維耦合場計算,為低溫裝置的傳熱計算提供了新型計算方法與參考依據(jù)。

(2)對于大尺寸試樣,將其表面發(fā)射率由0.95降至0.01后才可得到較為均勻的溫度分布,自身溫差由5.586 K降至0.245 K;隨著試樣尺寸減小,試樣自身溫度分布逐漸趨于均勻,對于本文中的小尺寸試樣3(尺寸為3 mm×30 mm×2.5 mm),其自身溫差幾乎不受自身表面發(fā)射率的影響。

(3)隨著窗口的表面發(fā)射率由0.94降至0.01,大尺寸試樣1的自身溫差由5.586 K降至3.195 K,而小尺寸試樣3始終保持較為均勻的溫度分布,其溫差由0.295 K降至0.286 K。因此,小尺寸試樣的溫度分布幾乎不受材料表面發(fā)射率和窗口參數(shù)的影響,且自身溫差始終保持最小,滿足實驗要求。