曾 磊 ,石友安 ,孔榮宗,賀立新,桂業偉
(中國空氣動力研究與發展中心,四川綿陽 621000)
熱環境地面風洞試驗是考核理論計算方法,為防熱設計提供依據的重要手段。地面試驗的測熱精度直接影響著計算方法和防熱結構的選取,是決定飛行成敗與否的關鍵因素之一。薄膜電阻溫度計作為激波風洞中最常用的熱流測量傳感器,試驗工作人員從傳感器制作、標定、安裝,數據采集系統的采樣頻率等多方面入手以提高測熱試驗精度。實際上,薄膜電阻溫度計在試驗中測量得到的電壓信息(即溫度值變化量),還需要通過理論分析和計算將溫度信息轉化為熱流值。一直以來,薄膜電阻溫度計的理論分析都是從一維兩層介質出發,推導簡化為一維半無限數學模型,再經過數值積分(Cook-Felderman公式)或熱電模擬網絡轉化獲得傳感器表面熱流值。但是,薄膜電阻溫度計是三維的,其表面鉑層內也是有溫度分布的,且隨著鉑層溫度的升高進入鉑層內的熱流會逐漸降低。所以,試驗中測量得到的溫度信息實際上是由一個時變的熱壁熱流引起的,這在后期的數據處理中是必須予以考慮的。
另外,目前常用的由溫度變化量計算熱流值所采用的Cook-Felderman公式或熱電模擬網絡都是基于一維半無限假設的數據處理方法。實際上,熱電模擬網絡的處理方法也是由Cook-Felderman公式轉化而來的,這種數值積分的方法對于測試噪聲有著累積和放大的效應,使得處理得到的熱流值存在較大的波動。
筆者從三維薄膜電阻溫度計計算模型入手,對比研究一維簡化、一維半無限簡化計算方法對鉑層內溫升規律的影響,研究熱壁熱流對表面溫升的影響,并由此建立可行的基于熱傳導反問題計算方法的表面溫升——熱流值的計算方法,為提高測熱精度提供可行的分析和數據處理手段。
(1)一維簡化影響分析
薄膜電阻溫度計的外形如圖1,基體材料為硼硅酸玻璃(ρ=2500kg/m3,Cp=1000J/kg?K,k=0.66W/m?s),表面鍍層的金屬薄膜材料為鉑(ρ=21450kg/m3,Cp=132J/kg?K,k=71.1W/m?s)。

圖1 薄膜傳感器的結構示意圖Fig.1 Structure of thin film resistance thermometer
三維熱傳導計算采用有限元計算方法,直角坐標系下的熱傳導控制方程為

傳感器網格如圖2所示。三維計算網格中取鉑層厚度為 5μ m,玻璃基體長度為 20mm,直徑為3mm,為了簡化計算鉑層外形取為長方體。計算了不同來流熱流條件下,采用三維傳熱計算與一維傳熱計算鉑層內溫度的分布情況,鉑層厚均取5μ m,初始溫度均為300K,三維傳熱計算時除鉑層所在的上表面取恒熱流邊界條件外,其余面均取絕熱邊界條件,這與試驗條件中上表面暴露在模型表面,其余各面均用航空密封泥填充的邊界條件是相適應的。圖3為1MW/m2熱流條件下,鉑層內溫度對比圖,由該圖可見采用三維傳熱算法得到的鉑層內溫度平均值,與一維傳熱計算得到的鉑層內溫度值基本相同。同時,分別采用三維和一維傳熱理論計算了表面熱流由85kW/m2至5MW/m2條件下,鉑層內溫度隨時間變化規律,兩者所得結果基本相同,這說明由于玻璃基體材料導熱系數很小,鉑層附近玻璃基體的橫向傳熱是可以忽略的。由此可以說明薄膜電阻溫度計可以作為一維問題處理。

圖2 傳感器網格圖Fig.2 Grid of sensor

圖3 鉑層內溫度對比圖Fig.3 Comparison of temperature of film
(2)半無限簡化影響分析
目前的處理方法是將薄膜電阻溫度計當作一維半無限體處理,忽略鉑層厚度的影響。圖 4為1MW/m2熱流條件下,鉑層5μ m厚時薄膜傳感器表面溫度分布云圖,可見鉑層的溫度是低于玻璃基體表面的溫度的。圖5所示為表面1MW/m2熱流下,不同鉑膜厚度對鉑層內溫升的影響。當膜厚為5、2和1μ m時,其表面溫升與采用半無限假設得到的溫升相差分別為8%、4%和2%。這說明隨著鉑層厚度的變薄,一維半無限簡化理論逐漸適用于薄膜電阻溫度計的傳熱分析。

圖4 傳感器表面溫度分布Fig.4 Temperature distribution of sensor's surface
(3)熱壁熱流影響分析
之前兩小節中給定的表面熱流條件均為進入傳感器的熱流值,即為熱壁熱流。實際上,由測量得到的電信號(即表面溫度信號)直接轉換得到熱流信息是熱壁熱流信息,需要通過轉換得到冷壁熱流值才能應用于后期的燒蝕試驗和防熱設計。熱壁熱流與冷壁熱流的關系可以用測試風洞的總溫和傳感器被加熱面的溫度來確定,具體如公式(2)

圖5 不同鉑層厚度溫升對比圖Fig.5 Temperature in different films

圖6 熱壁熱流對表面溫升的影響Fig.6 Effect of heat flux on temperature rising

圖7 表面溫升對熱流的影響Fig.7 Effect of temperature rising on heat flux

以中國空氣動力研究與發展中心的φ 2m激波風洞為例,當T∞=58K,Ma=9時,風洞運行總溫約為1000K,若給定1MW/m2的冷壁熱流,考慮表面溫升與熱流的耦合效應,則在10ms之后進入傳感器內部的熱流為冷壁熱流的88%,表面溫升比不考慮耦合效應的情況下低了近10%。
由此可見,在使用測量得到的溫度數據計算表面冷壁熱流時,必須要考慮傳感器結構溫升與熱流的耦合影響。
現有的薄膜電阻溫度計測熱數據處理方法采用的是建立在一維半無限假設基礎上的Cook-Felderman公式(公式3),以及在實際測熱試驗中采用的熱電模擬網絡,具體處理方法見文獻[1]。

熱流測試原理上是測量熱流引起的溫度變化,通過數據處理得到熱流,本質上是一個熱傳導反問題。從數學上講,熱傳導反問題屬于偏微分方程逆問題,是個不適定問題。它的不適定性主要體現在不滿足穩定性上,即解不連續依賴于數據,表現為微小的誤差可能導致計算結果的巨大變化。由于不適定的存在,熱傳導反問題的求解難度比相應的正問題大得多。
對于測試數據而言,誤差是無法避免的,即:Tm=Texact+ε,ε為測量誤差。

分母實際上就變成了一個放大因子,對測溫誤差進行放大。這就是Cook-Federman公式對誤差放大的原因之一。

式中符號意義:觀測點的坐標 xm;測量溫度Tm;辨識溫度 T;誤差ε。
(1)共軛梯度法(CGM:Conjugate Gradient Method)
共軛梯度法也稱為迭代正則化方法,其基本思想則是將優化問題分解為熱傳導正問題、靈敏度問題和伴隨變量問題來進行求解。
梯度的計算與伴隨變量有關,伴隨變量可通過求解伴隨方程來獲取。為得到伴隨方程,首先將目標函數(4)寫為如下的擴展形式

式中P稱為伴隨變量。對(5)式后半部分分部積分后再做變分,整理得到伴隨方程


初值條件 :t=tmax:P=0
1.2.2 對照組干預措施 在接受基礎干預措施的同時,由醫務人員在就醫期間和門診隨診時采取常規健康教育:包括醫學營養指導、運動指導、血糖監測指導等。
整理得到,目標函數對q的梯度為

步長由下式計算

式中ΔT是由Δ q=pn引起的各時刻溫度場的變化值。
在優化過程中,收斂準則根據輸出誤差原則獲得,即

σ為測量誤差的均方差,M為測點數
至此,共軛梯度法中涉及的梯度,步長均已求出,共軛梯度可由梯度、步長構造。具體優化細節可參閱文獻[2]。
(2)順序函數法(FSM:Function Specification Method)
順序函數法是基于熱傳導過程具有阻尼性和延遲性提出的。它不同于共軛梯度法的全時間域估計,而是按時間順序估計表面熱流。用此方法辨識未知熱流q(t)時,問題可以描述成:已知tM時刻的前M-1個時刻處的熱流q1,q2……qM-1和后r個時刻的溫度值TM,TM+1,……RM+r-1,辨識tM時刻的熱流qM。因此,辨識qM時,對應的正問題應該體現在時間域[tM-1,tM+r-1]內的熱傳導過程。選取如下目標泛函

由熱傳導過程的延遲性可知,T(tM+n)不僅取決于qM+n,還取決于qM+1,qM+2……qM+n-1,因此,為了辨識qM,必須假設它們的值。現假設各邊界點上qM至qM+r-1是線性變化的,當時間間隔是常數時,即

則目標函數式(10)僅為qm的函數,關于qm求導,即得到辨識值

至此,辨識邊界熱流時序值的順序函數法已經完全建立?;静襟E是在時間方向上按順序向后推進辨識。根據初始q0和 T0,辨識得到 q1和 T1。它們又作為辨識下一時刻q2的初始條件,重復上述步驟進行辨識。依此類推,即可得到q(t)。
(3)不同方法的對比分析
應用Cook-Felderman公式計算方法、共扼梯度法、順序函數法對中國空氣動力研究與發展中心的φ 2m激波風洞某車次測熱試驗數據進行了處理分析,并考慮第一小節中熱壁熱流的影響。
如圖8所示,采用Cook-Felderman公式計算得到的熱流數據波動相對較大,且在0.02s之后計算熱流值與導熱反問題計算得到熱流值相差較大。且由于薄膜傳感器鍍膜厚度很薄(<1μ m),所以Cook-Felderman公式計算結果并沒有出現明顯的滯后現象。

圖8 不同測熱數據處理方法比較Fig.8 Comparison between different data processing methods
圖9為以數據處理后得到的熱流值為邊界條件,計算得到的傳感器表面溫升情況。由圖可見在全時間域內,共扼梯度法反算得到的溫度值與實測溫度值都有很好的吻合,而Cook-Felderman方法反算得到的溫度值在0.02s之后明顯高于實測值,這是由于Cook-Felderman公式本身有對誤差積累和放大的作用。

圖9 反演熱流再計算溫度對比Fig.9 Comparison between different temperatures
圖10為熱壁熱流效應對采用順序函數法辨識表面熱流值的影響。由圖可見考慮表面溫升對熱流的影響后,在15~20ms時間段內,溫升為 20~30K的情況下,得到的傳感器表面冷壁熱流比熱壁熱流高近5%。由(2)式可知,隨著表面溫度的升高,冷壁熱流與熱壁熱流相差越來越大。同時,公式2也是用于對Cook-Felderman方法和共扼梯度法計算得到的熱流值的修正。

圖10 熱壁熱流效應影響Fig.10 Effect of hot-wall heat flux
(1)通過對比計算,薄膜電阻溫度計的數據處理采用一維半無限假設是合理的。當膜厚小于1μ m時,一維半無限模型和三維計算模型得到的結果相差不到2%;
(2)采用導熱反問題計算得到的熱流結果與采用Cook-Felderman方法得到的結果在趨勢上很相似。雖然Cook-Felderman方法得到的結果波動較大,但對于激波風洞測熱來說有效數據一般取在15~20ms相對平滑的階段,且取該時間段的平均值,其有效熱流值結果與導熱反問題計算結果基本符合。同時,Cook-Felderman方法的計算時間也比反問題的計算略短,可以很快地處理大量的測熱數據;
(3)導熱反問題計算方法具有良好的抗噪性和適定性,可以有效地獲得平滑的熱流數據,有利于反映熱流隨時間的變化規律,對于辨識時變的熱流值有一定的優勢;
總之,采用Cook-Felderman方法處理薄膜電阻溫度計在脈沖式風洞中測量得到的溫度信號是可行的,但測量傳感器必須符合一維半無限假設條件。熱傳導反問題的計算方法計算時間相對較長,適用于處理較長時間的時變熱流值。另外,反問題的計算方法可以推廣至軸對稱或多維的情況,適用于處理結構相對復雜的傳感器測量得到的溫度信號(如同軸熱電偶熱流傳感器),或大面積測熱數據。
[1] 張志成,劉初平,孔榮宗,等.高超聲速氣動熱和熱防護之再入飛行器熱環境地面試驗[M].北京:國防工業出版社,2002.
[2] 解可新,韓健,林友聯.最優化方法[M].天津:天津大學出版社,1997.
[3] 曾磊,桂業偉,賀立新,等.鍍層式同軸熱電偶數據處理方法研究[J].工程熱物理學報,2009,4:661-665.
[4] ALIFNOV O M.Inverse heat transfer problems[M].Springer-Verlag,Berlin,1994.
[5] BECK J V,BLACKWELL B,HAJI-SHEIKH A.Comparison of some inverse heat conduction methods using experimental data[J].Heat Mass Transfer.1996:3649-3657.
[6] JONES D P,PAT EL D K,ALEXANDER C.Numerical determination of surface temperature for in-depth thermocouples[R].AIAA-84-1763.