表面熱流辨識技術在邊界層轉捩位置測量中的應用初步研究

錢煒祺，周 宇，何開鋒，石友安

（中國空氣動力研究與發展中心空氣動力學國家重點實驗室，四川綿陽 621000）

0 引 言

飛行器飛行過程中，邊界層轉捩位置對飛行器的摩阻、壓阻及飛行性能有很大的影響，轉捩位置的確定是飛行器設計的關鍵技術之一。除了航空航天領域，邊界層轉捩位置的確定在船舶、機械制造、降噪、化工等領域也具有實用價值。然而，由于轉捩過程相當復雜，包含著層流失穩、T－S波的出現和放大、三維擾動的發展和流向渦的形成、流動二次失穩、湍流斑生成等一系列復雜的流動現象，使得目前對轉捩機理的認識還不完善，采用流體力學理論方法來計算轉捩位置較為困難。因此，迄今為止，實驗仍是做出轉捩判斷的最可靠的手段。

1 實驗中轉捩測量的主要方法

目前，實驗中轉捩的測量方法主要有升華法、脈動壓力測試技術、紅外成像熱圖技術、磷光粉熱圖顯示技術、液晶涂層法等［1－2］。其中，升華法是利用某些有升華現象的固態物質在邊界層層流區和湍流區的升華速度不同的原理來顯示轉捩位置，但由于需要較長的吹風時間，限制了該方法在高速風洞中的使用。脈動壓力測試技術是利用脈動壓力傳感器測量壓力脈動，根據脈動壓力的變化來確定轉捩位置，但風洞背景噪聲的干擾對測量結果有較大影響。由于層流和湍流區域的壁面溫度不一樣，理論上所有測壁面溫度的方法都可以用來測量轉捩，紅外成像熱圖技術、磷光粉熱圖顯示技術和液晶涂層法正是基于這一思想測量表面溫度的技術，目前在風洞中使用較廣，而在飛行試驗中由于需要獲取飛行器表面的圖像信息，限制了這些技術在飛行試驗中的使用。

在飛行試驗中測量轉捩的主要方法有：近壁熱電偶測溫法、光電二極管指示法、邊界層聲監測法等［3］。由于轉捩前后邊界層內脈動幅值增加，后兩種方法分別對邊界層內的光學信息和壓力脈動進行監測來判斷轉捩的發生，近壁熱電偶測溫法則是利用邊界層轉捩前后壁面溫度的變化來判斷轉捩。由于溫度測量技術相對較為成熟，因而目前在工程實際飛行試驗中主要采用的是近壁熱電偶測溫法。美國20世紀60年代進行的Reentry－F飛行試驗項目［4－5］和Sandia實驗室近年開展的SHARP－B2飛行試驗項目［3］，都采用的是近壁熱電偶測溫法來判斷轉捩。尤其是Reentry－F飛行試驗項目，該飛行試驗的主要目的就是進行飛行器表面的轉捩試驗。試驗飛行器是長3.962m、半錐角5°的尖錐體，最大飛行速度為20馬赫。尖錐體的彈體由鈹合金構成，在彈身上布置了21個溫度測點（如圖1所示），在每個測點沿厚度方向布置了4個熱電偶傳感器，其中距離表面最近的傳感器距離表面0.254mm。飛行試驗結束后，希望利用測點溫度測量數據來得出表面熱流沿飛行器表面的分布，進而確定轉捩位置，文獻［5］中采用了兩種方法來進行溫度測量數據處理，一種是用類似順序函數法的反演方法由測點溫度確定表面熱流［6］；另一種是認為測點溫度近似等于表面溫度，由溫度歷程對時間求導，導數最大的點對應的就是轉捩位置。結果表明，利用兩種方法確定出的轉捩位置基本一致（如圖2所示，圖中“↑”箭頭表示采用溫度歷程求導方法得出的轉捩位置；“↓”箭頭表示采用分析表面熱流辨識結果得出的轉捩位置），同時，這一結果和地面試驗結果也吻合較好。

圖1 Reentry－F飛行器及溫度測點布置Fig.1 Reentry－F spacecraft and its measurement positions

圖2 對飛行試驗數據用兩種方法確定出轉捩位置的比較Fig.2 Comparison of transition location

目前，雖然近壁熱電偶測溫法已在工程上實際應用，但該方法還存在兩個不足：

（1）該方法要求熱電偶安裝位置靠近表面，但靠近表面位置測點的溫升較高，對熱電偶傳感器的量程和精度要求較高。

（2）只有通過在防熱層內鉆孔的方式才能將溫度傳感器安裝在靠近表面的位置，這勢必對防熱層的結構強度產生較大影響。

因此，有必要探索新的處理方法和思路對近壁熱電偶測溫法進行改進與完善，可對熱電偶測點位置進行適當選取，使其既遠離表面，又能夠從測點溫度有效地辨識反演出表面熱流。

2 基于表面熱流辨識技術確定轉捩位置的處理方法

氣動熱參數辨識技術是熱傳導逆問題在氣動熱環境分析中的實際應用，包含防熱材料熱傳導系數辨識和表面熱流辨識兩方面研究內容［7］。轉捩位置的確定屬于表面熱流辨識問題，其基本原理如圖3所示，在轉捩測試平板內壁或材料內部沿流向布置多個溫度測點Pi，通過這些測點的溫度歷程來辨識表面熱流隨流向位置和時間的變化歷程。該方法與近壁熱電偶測溫法的主要差異在于測點離表面相對較遠。但是，溫度傳感器遠離表面后會帶來表面熱流辨識的不適定性問題。體現在兩方面：一是在厚度方向，由于熱傳導的延遲性和耗散性，測點溫度信息可能無法反映表面熱流在時間方向上的某些變化；二是在流動方向，同樣由于熱傳導在x方向的耗散，表面熱流在x方向上的變化可能在測點信息中體現不出來。對于這兩個問題，需要通過仿真計算辨識的方法來確定傳感器安裝位置，即首先根據理論計算或風洞實驗事先確定轉捩位置附近表面熱流的大致形式，然后利用熱傳導正問題計算出測點上的溫度歷程，并在該溫度歷程上疊加隨機誤差來模擬測點溫度的測量值，用來辨識出表面熱流，最后將熱流辨識結果與給定值比較，如果二者符合，則當前的測點位置是合理的，否則需調整測點位置。

圖3 基于表面熱流辨識技術的流動轉捩測試示意圖Fig.3 Sketch of transition determination method based on surface heat flux estimation

上述仿真辨識技術的核心是表面熱流的辨識算法，對于一維傳熱問題的表面熱流辨識算法，文獻［8－9］中有詳細的介紹，下面簡要給出二維傳熱的表面熱流辨識算法。仍以圖3為例，設長度方向（流向）的范圍為x∈［0，H］，厚度方向的范圍為y∈［0，L］，時間方向為t∈［0，tf］。二維傳熱的控制方程為：

初始條件：t＝0：T＝T0。

式中ρ、Cp、k分別為材料的密度、比熱和熱傳導系數，Q為隨時間和空間變化的表面熱流函數，就辨識問題而言，該函數是未知函數，通過如下的觀測信息來確定。測點Pi（i＝1，M）的坐標為（xi，0），對應的溫度測量歷程：

式中，T表示測點溫度的真值，v（t）表示測量誤差噪聲。因此，辨識問題等價于求合適的Q（x，t）使如下目標函數達極小的優化問題：

其中，T（xi，t，Q）表示表面熱流為Q時的測點溫度歷程，通過有限元計算方法來求取；?T為相應的實測值。由于Q（x，t）是空間和時間的函數，因此式（3）實際上是一個泛函表達式，考慮到有限元計算方法對計算空間域和時間域的離散，此時的待優化變量取為各時刻受熱邊界節點上的熱流值，即

其中Nj為受熱邊界離散單元節點數，Np為時間方向離散層數。則式（3）的優化變成了參數優化問題，可以采用梯度類優化算法來進行處理。采用伴隨方程推導的方法來求梯度，首先將目標函數寫為如下的擴展形式：

式中λ稱為伴隨變量。對（5）式進行分部積分再做變分后［8，10］得伴隨變量應滿足：

時域邊界條件：t＝tf，λ（tf）＝0。

同時，可導出目標函數對Q的導數為：

由式（6）、（7）求出梯度值后，即可用最速下降法或共軛梯度法進行優化，注意為克服不適定性的影響，當計算結果滿足如下收斂準則時終止優化迭代計算：

式中σ為溫度測量值的標準差。

3 仿真算例

3.1 X－43飛行器的流動轉捩分析

X－43飛行器馬赫7飛行試驗中，T021點位置表面熱流隨時間變化歷程如圖4所示［11］。從圖中可以看出T021點在哪些時間段處于層流狀態，哪些時間段處于湍流狀態，由此可將多個測點的信息組合起來進行轉捩位置判斷。現考慮采用表面熱流辨識方法，在飛行器殼體內壁測量溫度變化歷程，將殼體內的傳熱簡化為一維傳熱模型，由溫度歷程來辨識出表面熱流，進而對表面流態做出判斷。設防熱層是厚度10mm的不銹鋼，利用圖4中的熱流仿真計算出內壁測點溫度，考慮2%的測量噪聲，得到測點溫度歷程如圖5所示；再由此溫度歷程來辨識表面熱流，得到的辨識結果與圖4中給定熱流的比較如圖6所示，從中可以看到，辨識結果較好地再現了表面熱流的變化情況，能夠用于各時刻表面流態判讀。

3.2 二維平板流動轉捩位置測量仿真分析

二維平板模型及其表面上有轉捩時的表面熱流分布示意圖如圖3所示，設平板模型為長1m厚 10mm的合金，給定時空變化熱流如圖7所示。可以看到，t＝30s后表面熱流逐漸達到穩定狀態。針對計算域生成如圖8所示的有限元網格，利用給定熱流可以計算出內壁單元節點的溫度變化歷程。

圖7 二維平板轉捩分析的給定熱流時空分布Fig.7 Temporal and spatial function of twodimensional plate surface heat flux

圖8 有限元計算網格Fig.8 Computational grid of Finite Element Method

分別考慮內壁x方向均勻安裝10個和20個測點的情況，利用這些測點的溫度測量歷程，用上一節中的方法來辨識表面熱流，在t＝30s時刻得出的x方向熱流辨識結果如圖9中的“Estimated Q”所示，并與給定值“Exact Q”進行了比較。可以看到，在安裝20個測點的情況下，辨識結果能較好地反映出表面熱流的變化情況，對轉捩位置進行較好的預測。

圖9 給定熱流與辨識結果比較（不考慮測量誤差）Fig.9 Comparison of surface heat flux（no measurement noise）

進一步在測量結果中考慮1%的隨機測量誤差，圖10給出了此時的辨識結果。可以看出：在此情況下，安裝20個測點的辨識結果仍能較好地再現出表面熱流的變化情況，對轉捩位置進行較好預測。

此外，針對這一表面熱流辨識問題，也可以考慮在20個測點位置都假設為一維傳熱問題，用一維辨識方法來進行計算，圖11中的虛線給出了這樣處理得出的t＝30s時刻x方向表面熱流分布情況，與給得出的t＝30s時刻x方向表面熱流分布情況，與給定值、二維傳熱辨識結果有較明顯差異，得出的轉捩區域明顯大于真實情況。因此，此時采用基于二維傳熱模型的表面熱流辨識方法是有必要的。

圖10 給定熱流與辨識結果比較（考慮1%測量誤差）Fig.10 Comparison of surface heat flux（with measurement noise of 1%relative error）

圖11 一維和二維傳熱模型辨識結果比較Fig.11 Comparison of estimated heat flux for one－and two－dimensional heat conduction model

4 結 論

介紹了風洞實驗和飛行試驗中測量轉捩的主要方法及其優缺點。提出了基于表面熱流辨識技術確定轉捩位置的基本思想和處理方法，與目前工程上主要采用的近壁熱電偶測溫法相比，表面熱流辨識方法的優勢在于測點可以距離受熱表面相對較遠，從而放寬對溫度傳感器量程和結構強度的要求。但測點的位置也不是沒有限制，因為測點越遠離受熱面，辨識問題的不適定性越強，因此需要采用仿真辨識方法對傳感器安裝位置進行合理選取。在給出二維傳熱模型表面熱流辨識算法的基礎上，對兩個算例進行了仿真辨識分析。結果表明：基于表面熱流辨識技術確定轉捩位置是可行的，能給出較為準確的轉捩區域判斷。

［1］ 卞於中，徐鐵軍，范潔川，等.邊界層綜合診斷技術研究［J］.流體力學實驗與測量，2000，14（4）：7－13.

［2］ 胡成行，黃敘輝，李紅梅，等.應用脈動壓力測試技術探測邊界層轉捩［J］.流體力學實驗與測量，2002，16（2）：67－71.

［3］ KUNTZ D W，POTTER D L.Boundary layer transition and hypersonic flight testing［R］.AIAA－2007－308，2007.

［4］ ILIFF K W，SHAFER M F.A comparison of hypersonic vehicle flight and prediction results［R］.NASA TM－104313，1995.

［5］ WRIGHT R L，ZOBY E V.Flight boundary layer transition measurements on a slender cone at Mach 20［R］.AIAA－77－719，1977.

［6］ HOWARD F G.Single－thermocouple method for determining heat flux to a thermally thick wall［R］.NASATND－4737，1968.

［7］ 錢煒祺，何開鋒，桂業偉，等.飛行器氣動熱參數辨識研究進展［A］.第六屆全國流體力學青年研討會論文集［C］，浙江杭州，2009：79－84.

［8］ ALIFANOV O M.Inverse heat transfer problems［M］.Springer－Verlag，Berlin，1994.

［9］ 錢煒祺，蔡金獅.再入航天飛機表面熱流密度辨識［J］.宇航學報，2000，21（4）：1－6.

［10］錢煒祺，何開鋒，桂業偉，等.非穩態表面熱流反演算法研究［J］.空氣動力學學報，2010，28（2）：155－161.

［11］BERRY S，DARYABEIGI K，WURSTER K，et al.Boundary layer transition on X－43A［R］.AIAA－2008－3736，2008.