表面熱流辨識技術(shù)在邊界層轉(zhuǎn)捩位置測量中的應(yīng)用初步研究

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

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

0 引 言

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

1 實驗中轉(zhuǎn)捩測量的主要方法

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

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

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

圖2 對飛行試驗數(shù)據(jù)用兩種方法確定出轉(zhuǎn)捩位置的比較Fig.2 Comparison of transition location

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

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

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

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

2 基于表面熱流辨識技術(shù)確定轉(zhuǎn)捩位置的處理方法

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

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

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

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

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

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

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

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

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

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

同時，可導(dǎo)出目標函數(shù)對Q的導(dǎo)數(shù)為：

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

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

3 仿真算例

3.1 X－43飛行器的流動轉(zhuǎn)捩分析

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

3.2 二維平板流動轉(zhuǎn)捩位置測量仿真分析

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

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

圖8 有限元計算網(wǎng)格Fig.8 Computational grid of Finite Element Method

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

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

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

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

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

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

4 結(jié) 論

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

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

［2］ 胡成行，黃敘輝，李紅梅，等.應(yīng)用脈動壓力測試技術(shù)探測邊界層轉(zhuǎn)捩［J］.流體力學(xué)實驗與測量，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］ 錢煒祺，何開鋒，桂業(yè)偉，等.飛行器氣動熱參數(shù)辨識研究進展［A］.第六屆全國流體力學(xué)青年研討會論文集［C］，浙江杭州，2009：79－84.

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

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

［10］錢煒祺，何開鋒，桂業(yè)偉，等.非穩(wěn)態(tài)表面熱流反演算法研究［J］.空氣動力學(xué)學(xué)報，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.