塞塊量熱計的熱流計算與修正方法研究

朱新新, 王 輝, 楊 凱, 朱 濤, 楊慶濤, 劉進博

(中國空氣動力研究與發展中心 超高速空氣動力研究所, 四川 綿陽 621000)

0 引 言

地面防熱試驗是高超聲速飛行器設計與考核的必要環節[1-2]。隨著航天技術的飛速發展和高超聲速飛行器精細化設計熱考核需求的提升[3-4]，地面防熱試驗中對熱流等測量參數的精準度要求越來越高[5]。在以電弧風洞為代表的地面防熱試驗設備的流場校測中，塞塊量熱計常被選作測量熱流的主要傳感器[6]。塞塊量熱計可分為空氣隙結構和隔熱套結構兩類[7-8]，空氣隙結構的裝配工藝較難控制，更適合駐點熱流測量；隔熱套結構的工藝易于實現且不受流場條件限制。兩種結構的熱流計算均采用絕熱假設的近似計算方法，這種近似計算會帶來較大誤差，國內外研究者力圖對其進行優化和修正：

一方面，通過改進塞塊量熱計的隔熱結構使其盡量滿足絕熱條件，如楊慶濤等[9]提出的尖楔狀隔熱套結構、許考等[10]提出的線接觸隔熱結構以及朱新新等[11]提出的中空隔熱套封裝結構，并給出了這種隔熱套結構的相對優越性。即便如此，這些改進結構仍然未能達到較理想的絕熱條件。

另一方面，是對熱流計算方法進行修正。陳連忠[12]考慮到塞塊熱物性會隨著溫度的升高而改變，對其進行了一定程度的修正，降低了一部分計算誤差。Mark Hightower等[13]考慮了空氣隙結構下的熱損失和變熱物性兩個因素，建立了一個近似的指數衰減理論模型，并基于該模型給出了修正公式，在一定程度上降低了熱流測量誤差。但正如文獻[13]所述，將基于測試數據得到的熱流修正參數代入仿真模型進行熱流驗證會存在一定偏差。出現偏差的原因有兩方面：一方面，衰減理論模型中假定加熱時隔熱結構的溫度始終不變，而實際情況是一直升高且不能忽略；另一方面，整個修正方法始終沒有獲得可溯源的熱流真值，僅是一個基于流場測試數據的相對修正。有鑒于此，本文擬通過熱流標定系統[14-15]獲取可溯源的準確熱流，并以此為依據建立塞塊量熱計的仿真模型，對熱流測量結果進行計算修正。

1 塞塊量熱計的傳熱模型和熱流計算方法

塞塊量熱計測量熱流是基于能量守恒原理，假定塞塊除受熱面外，其余均為絕熱面。但在實際研制和使用中，受材料和裝配工藝等限制，很難達到理想絕熱條件，熱損失始終存在。為減少熱損失，朱新新等[11]設計了一種中空型隔熱套結構的塞塊量熱計，其結構如圖1所示，從左上至右下依次為：耐熱填充劑、封裝殼、焊有熱電偶的塞塊、中空隔熱套、螺釘。這種塞塊量熱計的優勢在于塞塊熱損失較少且能將塞塊量熱計熱流標定時內部的結構和熱阻等傳熱特性保留到熱流測量試驗中。

圖1 塞塊量熱計結構示意圖

塞塊量熱計的傳熱模型如圖2所示。其前表面為受熱面，熱流(qin)入射塞塊后，塞塊溫度快速上升，可以通過測量塞塊的溫升率計算得到熱流值。熱流計算誤差主要來自4個方面：一是塞塊前端與耐熱填充劑的熱傳導，熱導率很低的耐熱填充劑前端暴露于熱環境中，局部溫度會較高，表現為對塞塊放熱，但其后端溫度較低，表現為對塞塊吸熱，總的表現為吸熱(qcond1)；二是隔熱套前端內側面和后表面的內側面對塞塊的吸熱(qcond2)；三是塞塊對外熱輻射帶來的熱損失(qrad)；四是塞塊熱物性發生了變化，若按常溫時的熱物性計算熱流，會帶來一定熱損失(qvar)。

圖2 塞塊量熱計的傳熱模型

假設圖2中圓柱狀塞塊量熱計的長度為l，橫截面積為A，則根據能量守恒原理可得，任意時刻單位時間內的入射熱能等于儲存至塞塊里的熱能與損失熱能之和：

qinA=Estore+Eloss

(1)

在某時刻t，塞塊儲存的熱能Estore可根據其溫升率計算得出：

(2)

其中，Tave為塞塊平均溫度，ρ(Tave)和cp(Tave)為塞塊在Tave時的密度和質量定壓熱容。式(2)取實時溫度下的熱物性，不會產生如圖2所示的熱損失qvar。則總的熱損失包括前文所述的熱傳導損失和輻射熱損失：

Eloss=qradS+qcond1S1+qcond2S2

(3)

其中，S為塞塊表面積，S1為塞塊與填充劑的接觸面積，S2為塞塊與隔熱套的接觸面積。由于塞塊體積和長度很小，且具有較高的熱導率，則有以下近似關系成立：

(4)

再設接觸面上耐熱填充劑的某接觸位置z的接觸熱阻為Rcond1(z)，溫度為T1(t,z)，則由接觸熱阻定義可推導出t時刻該接觸面上的平均熱流qcond1，同理可得qcond2：

(5)

(6)

其中，Z1為塞塊和填充劑的接觸面長度，Z2為塞塊和隔熱套的接觸面長度，Tb(t,z)取塞塊后表面的測點溫度Tb(t)。

另外，塞塊量熱計測量熱流屬于瞬態測量，塞塊溫度并不會很高(一般不超過600 K)，其輻射損失的熱能與入射熱能相差至少5個數量級，所以輻射熱損失可以忽略。綜合式(1)～(6)，可得塞塊量熱計的熱流計算模型如下：

(7)

假設隔熱結構完全絕熱(Rcond1(z)和Rcond2(z)無窮大)，且塞塊溫升較小(熱物性變化小)，則可得到常用的簡化公式(8)：

(8)

其中，ρ、cp為常溫時的塞塊密度和質量定壓熱容。

事實上，由于忽略了較多的熱損失(這在后文的仿真計算和標定試驗中會進一步說明)，這種假設對于隔熱套結構的塞塊量熱計很難成立，按照式(8)計算熱流會產生較大誤差；而由于無法獲得塞塊與填充劑和隔熱套之間的熱流qcond1和qcond2，采用直接解析計算方法求解qin相當困難。因此，本文將基于熱流標定系統和塞塊量熱計仿真模型對塞塊量熱計的熱流測量結果進行計算和修正。

2 塞塊量熱計的仿真模型和影響因素分析

依據圖1的塞塊量熱計結構建立如圖3所示的塞塊量熱計軸對稱仿真模型。塞塊直徑為4 mm，長度為7 mm。在圖3中，上表面為熱流加載面，除了8個不同材料間的接觸面(圖中編號1～8)，其余邊界面設為絕熱條件。整個仿真模型建立的難點是8個接觸面熱阻值的設置，因此，首先需要分析各個接觸熱阻對熱流測量結果的影響。

塞塊量熱計包括4種材料：無氧銅(塞塊和封裝殼)、氧化鋁(填充劑)、玻璃鋼(隔熱套)和螺釘(不銹鋼)，根據查表后擬合得到的溫度函數設置各材料的熱物性參數。接觸熱阻主要跟表面粗糙度和受力情況有關，一般受力越大、接觸面越光滑，接觸熱阻越小。接觸面1和4的受力情況與材料均相同，其接觸熱阻相同，定義為R1,4；接觸面2和3材料相同，且兩個接觸面受力都控制為較小力(為減小熱損失)，也可定義接觸熱阻為同值，以R2,3表示；接觸面6盡管材料與接觸面2和3一樣，但其裝配時受緊固螺釘擠壓，受力明顯大于后兩者，故將其定義為R6；接觸面5、7和8的熱阻分別定義為R5、R7和R8。

結合各接觸面熱阻特點，參考文獻[16]對接觸熱阻的研究，對各接觸熱阻值作表1所示的定義和說明。

圖3 塞塊式量熱計的仿真模型

表1 不同接觸面的熱阻值Table 1 Thermal resistance of different contact surfaces

按表中第二列接觸熱阻估計值設定塞塊量熱計各接觸面的熱阻。不失一般性，設入射熱流qstd為1 MW/m2，加熱時間3 s，加熱前所有部件溫度初始值設為300 K。通過仿真計算，可以得到塞塊后表面溫度(Tb)、填充劑中靠近接觸面位置的平均溫度(T1)、隔熱套中靠近接觸面位置的平均溫度(T2)的溫升曲線，如圖4所示。隨著測試時間推移，塞塊溫升最快、填充劑其次、隔熱套最慢，必然導致塞塊與填充劑、塞塊與隔熱套的溫差持續升高，其結果就是：在被測熱流不變的前提下，測試時間越長，熱損失越多，在測試時間段內提取溫升數據的時間越晚，計算得到的熱流誤差就越大。

圖4 塞塊式量熱計不同部位的溫升曲線

為定量描述誤差大小，定義一個計算誤差ε：

(9)

其中，qcal為采用式(8)計算得到的熱流。

將圖4中塞塊的溫升數據進行分段擬合，將擬合得到的溫升率分別代入式(8)計算熱流。如表2所示，在測試時間段內提取溫升數據的時間越晚、測得熱流越低，誤差越大。即使進行了熱物性修正(表2最后一列)，也僅稍有改善，這也進一步驗證了式(7)中熱損失的存在和不可忽略性。因此，在標定試驗和實際熱流測量中應盡可能早地提取溫升數據。本文結合送進系統到達指定位置的時間(0.20 s)和塞塊量熱計的響應時間(0.25 s)，統一設定提取時間為0.5～0.7 s。

表2 不同時間段的熱流及誤差Table 2 The heat flux and error in different time periods

為分析接觸熱阻對熱流測量結果的影響，表3給出了各個接觸熱阻單獨變化時引起的計算誤差變化情況。由表3可知，當接觸熱阻縮小或放大1個數量級時，對計算誤差影響最大的是R2,3(塞塊與隔熱套之間的接觸熱阻)，其余5種接觸熱阻的變化對計算誤差影響很小。這是因為與塞塊直接接觸的僅有填充劑和隔熱套，而塞塊與填充劑的接觸面積小于塞塊與隔熱套的接觸面積，且塞塊與填充劑的溫差也小于塞塊與隔熱套的溫差(見圖4)，因此，R2,3對塞塊溫升率的影響最大，即整個塞塊量熱計仿真模型對接觸熱阻R2,3的值比較敏感，需精確設置。

表3 不同熱阻下的計算誤差Table 3 Calculation error under different thermal resistance condition

3 塞塊量熱計熱流計算的修正方法

塞塊量熱計的結構尺寸確定后，每加載一個不同的熱流qin，都會在固定時間段內輸出唯一的dTb/dt。將式(7)變形為：

(10)

其中，X表示某具體傳感器的固有屬性，可視作一未知常量。基于式(10)，就可以在使用塞塊量熱計之前，先在一個熱流已知的熱源系統中對其進行標定，獲得不同熱流下對應的dTb/dt。在塞塊量熱計仿真模型中加載不同熱流(如表4所示，其他邊界條件和初始值與圖4的仿真情況相同)，可得到對應的dTb/dt和熱流計算誤差。

表4 不同熱流下的計算誤差Table 4 Calculation error with different heat flux

仿真結果顯示，加載熱流越大，計算誤差也越大；但這種增大趨勢并非嚴格線性，這就給修正帶來一定困難。比較極端的做法是標定足夠多的試驗數據，甚至做成一張熱流分度表的形式，測試時通過溫升率查熱流分度表直接獲得對應熱流。但這并不現實：一方面，試驗中使用塞塊量熱計的數量較大，且不同塞塊量熱計的固有屬性X也不同；另一方面，會花費大量時間和成本在熱流標定環節上，且目前熱流標定系統的熱流范圍還比較窄，并未完全覆蓋測試區間。基于此，結合設備實際及應用需求，提供兩種方法：一種是直接比對標定修正方法；另一種是基于標定的數值計算修正方法。

3.1 直接比對標定修正方法

通過進一步開展仿真和標定試驗發現，在較小的熱流區間內，不同的加載熱流引起的計算誤差差異很小。因此，可以在一定熱流測量區間內均勻選取n個(n≥3)不同大小的熱流狀態，在每個熱流狀態下用標準熱流計測量得到參考熱流qstd，將其值作為縱軸，將同樣熱流狀態下由式(8)計算得到的塞塊量熱計熱流qcal的值作為橫軸，然后進行線性擬合，可得到該塞塊量熱計的標定曲線：

qstd=a+bqcal

(11)

其中，a和b為修正參數。在熱流測量中，就可以用獲得的修正參數a和b對計算熱流進行修正，得到修正熱流如下：

qcor=a+bqcal

(12)

需要說明的是，式(11)和(12)中的qcal不同。式(11)中的qcal是標定試驗中的計算熱流，其對應的真實入射熱流(qstd)是已知的，用于求解修正參數a和b；而式(12)中的qcal是測量試驗中計算得到的熱流，其對應的真實入射熱流未知，需要通過式(12)修正得到。

為評估上述修正方法的有效性，定義一個修正誤差ε′，用于定量評估修正后的誤差。

(13)

研制加工了5個塞塊量熱計(結構如圖1所示)。圖5為1#塞塊量熱計在5個不同熱流下獲得的標定曲線，其線性度極好(R2=0.9999)，通過線性擬合得到其修正參數a和b分別為-0.039和1.068。表5給出了1#塞塊量熱計的標定結果，包括5個標定狀態及每個狀態對應的計算誤差和修正誤差。

圖5 塞塊量熱計的熱流標定曲線

表5 不同狀態修正前后的誤差對比

表6給出了5個塞塊量熱計的標定曲線參數和對應的最大誤差。由表6的標定結果可知，直接比對標定修正方法能夠有效減小塞塊量熱計的熱流計算誤差，5個塞塊量熱計修正前的誤差為5.0%～8.9%，修正后的誤差降為0.4%～1.7%。

表6 不同塞塊量熱計修正前后的誤差對比

3.2 基于標定的數值計算修正方法

對于熱流標定系統不能覆蓋的較高或較低的熱流區間，擬借助求解塞塊量熱計仿真模型的途徑進行修正。如前文所述，塞塊量熱計仿真模型是否準確主要受接觸熱阻R2,3影響，因此，求解塞塊量熱計仿真模型主要就是求解接觸熱阻R2,3的值。

在熱流標定試驗中，能夠準確知道不同熱流狀態下進入塞塊量熱計的熱流，且能測得塞塊后表面相應的溫升率，此時將R2,3視作唯一未知量(相當于式(10)中的X)，并賦予一個預估初值；然后以實際測得的熱流值作為輸入條件，加載到該塞塊量熱計的仿真模型中，就能得到仿真模型輸出的溫升率；再將此溫升率與標定試驗中測得的溫升率進行比較，根據比較結果，改變R2,3的值；如此反復迭代，直至仿真溫升率與標定試驗溫升率的偏差小于誤差要求(如1%)，此時對應的R2,3就近似為該塞塊量熱計塞塊與隔熱套的最佳接觸熱阻。

采用上述方法求解接觸熱阻R2,3，有以下說明：

(1) 不同熱流條件下算出的接觸熱阻R2,3可能會略有不同，因為試驗測量結果并不僅受R2,3的影響，還受到熱電偶焊點工藝、標定系統精度以及數采誤差等影響，所以僅能得到一個較為準確的熱阻；此熱阻除了包含實際的熱阻影響，還包含了其他因素的誤差影響，可將其視作一個唯一影響該塞塊量熱計仿真模型的待求變量；

(2) 求解該熱阻理論上僅需一個標定狀態，但考慮到影響因素較多，應考慮使用2～3個熱流狀態相互驗證。

使用上述方法求解表5中1#塞塊量熱計的接觸熱阻R2,3，求解過程如表7所示。熱阻初始值設為1×10-3m2·K/W，第二輪迭代結果的溫升率偏差已經小于1%；繼續以相鄰兩個高熱流狀態進行驗證優化，最終確定接觸熱阻約為1.7×10-3m2·K/W。

表7 接觸熱阻的求解過程Table 7 Solution process of contact thermal resistance

受熱流標定系統能力限制，未能開展更高熱流的驗證，僅以表5中的狀態4和5這兩個相對較高的熱流狀態驗證上述方法的可行性，即假設狀態4和5的標準入射熱流值未知。對狀態4和5的入射熱流的求解過程如下(見表8)：

(1) 求仿真模型中的R2,3。表7中已求得其值為1.7×10-3m2·K/W。

表8 入射熱流的求解過程Table 8 Solution process of incident heat flux

(2) 求初始入射熱流。試驗中已測得狀態4和5的溫升率分別為77.70和91.15 K/s，根據式(8)可計算得到兩狀態的熱流初始值為1.869和2.192 MW/m2。

(3) 迭代優化入射熱流。不斷更新仿真模型中入射熱流的值，當仿真溫升率與已知的試驗溫升率偏差小于1%時，對應的仿真輸入熱流值即為較優入射熱流。按此迭代方法，得到狀態4和5的熱流修正值分別為1.930和2.300 MW/ m2，與實際標準熱流的誤差分別為1.02%和0.30%，小于修正前的4.2%和5.0%，可見該方法能夠有效減小計算熱流的誤差。

4 結 論

本文建立了塞塊量熱計的仿真模型，并給出了仿真模型的求解方法，仿真和試驗結果表明了該方法的可行性，這種求解思路也可應用于類似的傳熱仿真模型中。提出的直接比對標定方法和基于標定的數值計算修正方法均能降低塞塊量熱計的熱流測量誤差。考慮熱流標定系統自身誤差(3%)后，塞塊量熱計的熱流測量誤差小于5%；直接比對標定修正方法的誤差相對更小，但要求熱流標定系統能覆蓋被測熱流范圍；基于標定的數值計算修正方法對熱流標定系統要求較低，適用范圍更廣。

本文的兩種修正方法應用于中空型隔熱結構的封裝塞塊量熱計中，可以提高塞塊量熱計的熱流測量精準度，已在電弧風洞試驗中得到了驗證(與戈登計的風洞試驗比對測試結果見文獻[11])和推廣應用。