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

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

（中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng) 621000）

0 引 言

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

1 實(shí)驗(yàn)中轉(zhuǎn)捩測(cè)量的主要方法

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

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

圖1 Reentry－F飛行器及溫度測(cè)點(diǎn)布置Fig.1 Reentry－F spacecraft and its measurement positions

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

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

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

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

因此，有必要探索新的處理方法和思路對(duì)近壁熱電偶測(cè)溫法進(jìn)行改進(jìn)與完善，可對(duì)熱電偶測(cè)點(diǎn)位置進(jìn)行適當(dāng)選取，使其既遠(yuǎn)離表面，又能夠從測(cè)點(diǎn)溫度有效地辨識(shí)反演出表面熱流。

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

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

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

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

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

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

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

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

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

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

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

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

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

式中σ為溫度測(cè)量值的標(biāo)準(zhǔn)差。

3 仿真算例

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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