劉智侃，劉深深,2,*，劉驍，曾磊，代光月

1.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000 2.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室，綿陽 621000

臨近空間類高超聲速飛行器由于其獨特的飛行環(huán)境及特點，具有十分重要的應(yīng)用價值，已成為當(dāng)前高超聲速飛行器研究的重要發(fā)展方向。該類飛行器需要在稠密大氣中進(jìn)行長時間的高超聲速飛行，因而面臨嚴(yán)酷的氣動加熱問題。長時間氣動加熱帶來的溫升導(dǎo)致結(jié)構(gòu)受熱后產(chǎn)生靜態(tài)變形，而結(jié)構(gòu)變形進(jìn)一步導(dǎo)致流場氣動加熱的部分改變，造成表面熱流、溫度場、結(jié)構(gòu)應(yīng)力/應(yīng)變場之間發(fā)生復(fù)雜的耦合作用，從而引發(fā)綜合熱效應(yīng)[1-3]。獨立分散的熱環(huán)境、熱響應(yīng)和熱應(yīng)力分析方法人為割裂了實際存在的多物理場共同作用的復(fù)雜過程，對存在力/熱/結(jié)構(gòu)耦合嚴(yán)重的情況并不適用[4-5]，因此在當(dāng)前針對高超聲速飛行器的熱防護(hù)設(shè)計和優(yōu)化中必須綜合考慮氣動力/熱/結(jié)構(gòu)的多場耦合計算問題。

當(dāng)前多物理場耦合數(shù)值計算問題的求解主要是松耦合方法[6]，它的基本思想是將緊耦合聯(lián)立求解的問題簡化為遞推求解的過程：在極小的時間步內(nèi)各個模塊分別采用各自適用的求解器進(jìn)行求解，通過不斷地在流體域和固體域交界面上進(jìn)行數(shù)據(jù)傳遞來實現(xiàn)實時耦合。松耦合方法可以復(fù)用現(xiàn)有的發(fā)展比較成熟的氣動力、氣動熱、結(jié)構(gòu)熱響應(yīng)等求解器，同時可以通過控制交互時間步的大小來簡化計算開銷，因而對復(fù)雜的外形具有很好的適用性，應(yīng)用較為廣泛[7-8]。在松耦合方案中，不同物理場之間的數(shù)據(jù)傳遞方法是驅(qū)動和實現(xiàn)耦合計算的關(guān)鍵橋梁，其精度和效率共同決定了耦合計算能否在可承受的時間代價下獲得真實反映物理耦合過程的結(jié)果，因而是實現(xiàn)高效高精度耦合數(shù)值模擬的關(guān)鍵環(huán)節(jié)。

徑向基函數(shù)(Radial Basis Function, RBF)因具有形式簡單，不依賴于求解器離散格式、網(wǎng)格拓?fù)浣Y(jié)構(gòu)，易于編程且精度較高，能保證所傳遞的力載荷的守恒性[9]等特點在多場耦合計算問題中得到了廣泛應(yīng)用。Rendall等[10-12]開展了一系列研究，將其應(yīng)用于高超聲速機(jī)翼氣彈問題分析，表明其對CFD/CSD數(shù)據(jù)插值及大變形位移具有良好的適用性；郭中州等[13]發(fā)展的高效K近鄰-徑向基函數(shù)(KNN-RBF)動網(wǎng)格方法提高了網(wǎng)格變形的計算效率和對復(fù)雜外形的適用性；劉君等[14]利用RBFs-MSA混合算法實現(xiàn)任意幾何形狀的計算網(wǎng)格變形，能夠有效滿足黏性流動模擬需要的異向加密需求；王剛等[15]提出了基于貪心法的點選取方法以消除鄰近點人為選擇對插值精度的影響；為了減少RBF插值過程中支撐點數(shù)目，魏其等[16]提出峰值選擇法以減少迭代次數(shù)提高選點效率。但采用徑向基函數(shù)插值時存在的問題是徑向基本身是各向同性的，更適用于描述各向同性的物理量分布。而具備復(fù)雜外形的高超聲速飛行器物理量分布本身與局部流動特征高度相關(guān)且梯度變化劇烈，例如飛行器頭部、前緣、翼舵干擾區(qū)及機(jī)身大面積區(qū)域壓力、熱流等物理量分布差異巨大；除去上述因素外，由于復(fù)雜外形網(wǎng)格的劃分會在外形曲率或物理量變化劇烈區(qū)域進(jìn)行加密處理因而網(wǎng)格點分布也是各向異性的，人為引入了不同網(wǎng)格點分布的各向異性因素，由此采用各向同性基函數(shù)進(jìn)行高度各向異性物理量空間插值時存在缺陷，造成了精度及守恒性提升瓶頸。

為了解決上述問題，需要對徑向基函數(shù)進(jìn)行各向異性修正?？紤]減少人為網(wǎng)格劃分的影響，劉深深等[17]提出的一種基于網(wǎng)格劃分幾何尺度歸一化方法的各向異性修正方法取得了較好的效果，但該項工作未考慮物理量自身分布特征的影響。為了解決物理量各向異性分布對插值精度帶來的影響，本文考慮將物理量梯度的分布表征因子引入徑向基計算，提出了一種基于物理量梯度進(jìn)行徑向基各方向自適應(yīng)因子修正的新方法。采用高超聲速舵面、翼身組合體算例對該方法在單次數(shù)據(jù)傳遞中的可行性及效果進(jìn)行了驗證分析，同時基于該方法對壓縮拐角外形開展了氣動力/熱/結(jié)構(gòu)多場耦合數(shù)值模擬，并將耦合計算結(jié)果與風(fēng)洞試驗結(jié)果進(jìn)行了對比分析，對該方法在實際耦合問題計算長時間高頻次數(shù)據(jù)傳遞中的適用性及可靠性開展了分析驗證研究。

1 建立基于物理量梯度修正的改進(jìn)模型

1.1 徑向基函數(shù)數(shù)據(jù)傳遞基本原理

在d維歐幾里得空間上給定函數(shù)φ:R+→R，對于空間中n個不同的點{xi,gi|i=1,2,…,n}∈Rd?R，在這些固定點上的標(biāo)量值為{g1,g2,…,gn}。根據(jù)這些點需要確定一個函數(shù)：

(1)

使其滿足插值條件：

(2)

(3)

(4)

其中，添加的多項式為

p(x)=λ1+λ2x+λ3y+λ4z

(5)

需要添加的定解條件[18]為

(6)

依據(jù)基函數(shù)的定義和定解條件便可以推導(dǎo)傳遞矩陣，完成數(shù)據(jù)的傳遞。按照基函數(shù)作用域的范圍可以將RBF分為全局域RBF和緊支RBF，全局域RBF指每個RBF影響區(qū)域為整個定義域，如薄板樣條插值(Thin Plate Spline, TPS)、多重二次函數(shù)插值(Multiquadric, MQ)等基函數(shù)，而緊支RBF的影響區(qū)域為指定的影響半徑，如Wenland定義的C2緊支基函數(shù)[17]。這3種基函數(shù)是目前多場耦合計算中最為常用的3種基函數(shù)，本文研究主要針對上述基函數(shù)進(jìn)行，具體表達(dá)式如表1所示。

表1 3種徑向基函數(shù)

1.2 基于物理量梯度修正RBF方法

吳宗敏[18]指出理論上依據(jù)式(3)計算的歐式距離r是各向同性的，故所述的RBF插值方法更適用于物理上各向同性的問題。而實際流動物理量在各個方向的分布并不是一致的，如圖1所示,在機(jī)翼前緣附近，由于流動特征及網(wǎng)格劃分的影響，熱流變化在x和y方向變化十分劇烈，在z方向變化較為緩慢，在固定空間半徑范圍內(nèi)對于x和y方向的插值所選取的點在待插值點附近越集中插值精度越高，而物理量特征相差較大的D點等的引入會帶來更多空間分布擬合誤差，因此需要通過增大x和y方向的距離權(quán)重來排除物理量相差過大的點，同時相應(yīng)的在物理量分布變化緩慢的z方向，需要通過減小該方向的距離權(quán)重來增加相應(yīng)插值點E，以防止圖中所示的插值點過少的情況。根據(jù)上述思想，可以基于物理量梯度對徑向基計算中的x、y、z3方向進(jìn)行系數(shù)縮比調(diào)控，以求在相同的插值半徑選點范圍內(nèi)用來插值的點能夠更具有代表性和聚集性，從而使插值的結(jié)果更能表征物理量實際分布特征，進(jìn)而提高精度。

圖1 機(jī)翼前緣熱流分布及固定半徑內(nèi)選擇的點

針對徑向基函數(shù)中歐式距離r的計算進(jìn)行如下變換[17]：

(7)

其基本思路是對外形基于物理量梯度分布進(jìn)行適當(dāng)?shù)目s放使得網(wǎng)格點在不同方向上的分布更加均勻，其中c1、c2、c3由物理梯度決定，計算方法如下：

步驟1求解各個已知網(wǎng)格點上的物理量分布梯度。對于離散網(wǎng)格點的梯度求解本文采用的方法主要為最小二乘法。其基本原理是基于網(wǎng)格節(jié)點上物理量值的一階Taylor級數(shù)展開近似得到

(8)

式中:U為物理量的值。將式(8)應(yīng)用到與目標(biāo)節(jié)點i相連接的所有節(jié)點，得到超定約束線性方程組：

(9)

式中：Δ(·)ij=(·)j-(·)i;M表示與目標(biāo)節(jié)點i相連的網(wǎng)格節(jié)點數(shù)目；θj表示權(quán)系數(shù)，如面積、體積或距離等，本文的θj值均為1。通過求解上述超定方程組，可得物理量的梯度分布。

對于網(wǎng)格點物理量梯度的求解，需要用到網(wǎng)格點周圍的連接點信息，由附近連接點的物理量值進(jìn)行梯度的求解，對于網(wǎng)格節(jié)點間間距較大的節(jié)點該方法求得的梯度值精度較低。對于結(jié)構(gòu)/非結(jié)構(gòu)網(wǎng)格在初始插值時均需進(jìn)行一次搜尋并存儲相關(guān)連接信息，耦合計算中多次插值時不需要進(jìn)行重復(fù)搜索，因而不會帶來過多額外的計算開銷。

步驟2求得所有已知節(jié)點在各個方向的梯度總和Gx、Gy、Gz，計算公式為

(10)

式中：N為原始網(wǎng)格點的個數(shù)。

(11)

式中：i=1,2,3分別代表x、y、z3個方向。在此基礎(chǔ)上同時考慮幾何外形的尺度歸一，可以進(jìn)一步提高插值算法的魯棒性：

(12)

步驟4對所有已知網(wǎng)格點和待插值網(wǎng)格點的坐標(biāo)進(jìn)行基于物理量梯度的修正：

步驟5將縮放后的值代入式(3)中，得到更新后的rnew的表達(dá)式為

(13)

對比式(7)與式(13)可知：

(14)

式中:c1、c2、c3為基于物理量梯度分布得到的自適應(yīng)修正各向異性因子。

通過分析式(3)與式(7)可知，對于給定的空間半徑網(wǎng)格點的搜索是在一個球內(nèi)，而梯度修正后的網(wǎng)格點的搜索則是在一個橢球內(nèi)。對網(wǎng)格縮放后外形上某一點來說，距離r相同的所有點組成了一個球面，這等同于在原外形中距離為r的數(shù)據(jù)點組成了一個橢球面，即相同的r所包圍的數(shù)據(jù)點坐標(biāo)范圍根據(jù)不同方向的物理量梯度分布進(jìn)行了拉伸，變化劇烈的方向被拉長，變化平緩的方向被壓縮，進(jìn)而能夠降低物理量梯度和幾何外形的影響。

2 改進(jìn)方法的基礎(chǔ)算例驗證

高超聲速飛行器具有明顯的氣動力/熱/結(jié)構(gòu)的多場耦合現(xiàn)象，涉及到氣動力、熱流、溫度、位移等物理場的傳遞。對于氣動力、熱流等載荷傳遞，不僅需要點對點的精確插值，還需要整個交界面的通量守恒。一般而言，氣動力數(shù)據(jù)傳遞可以通過虛功原理來自動保證力學(xué)能量守恒，而對于熱流的傳遞來說，只能通過更精確的插值來保證交界面上的通量守恒[19]。因而熱流的傳遞是高超聲速飛行器耦合問題插值中的難點，選擇熱流插值作為考核變量。根據(jù)前文對基于物理量梯度修正RBF方法的介紹，采用文獻(xiàn)[17]中高超聲速飛行器中的舵面結(jié)構(gòu)和某一翼身組合體飛行器作為考核算例，對熱流傳遞的準(zhǔn)確性進(jìn)行對比分析。根據(jù)以往的研究[20]，樣條函數(shù)中的TPS方法、MQ方法和緊支C2基函數(shù)方法較為常用，本文將在此基礎(chǔ)上分析物理量梯度修正前后熱流插值效果。

2.1 高超聲速飛行器舵面模型

對于高超聲速飛行器來說，翼舵結(jié)構(gòu)是其一大典型外形特征，且在飛行過程中氣動力/熱/結(jié)構(gòu)耦合較為顯著。因此，本文針對高超聲速飛行器某一舵面進(jìn)行了熱流的插值傳遞計算與分析。舵面翼根長2.67 m，翼展長1.27 m，前緣后掠角為30°，前后緣半徑均為7 mm。對該舵面進(jìn)行網(wǎng)格劃分，得到流體域網(wǎng)格和固體域網(wǎng)格，分別如圖2、圖3所示。其中流場表面網(wǎng)格節(jié)點數(shù)4 073個，單元數(shù)4 016個；結(jié)構(gòu)表面網(wǎng)格節(jié)點數(shù)5 317個，單元數(shù)10 520個。

圖2 舵面流體域網(wǎng)格

圖3 舵面固體域網(wǎng)格

針對該翼面的外形特征，假設(shè)舵面上任意坐標(biāo)數(shù)值分別為(x,y)，則表面熱流分布表達(dá)式為

(15)

根據(jù)該熱流表達(dá)式可以得到舵面的熱流分布如圖4所示。

圖4 舵面流體域網(wǎng)格表面熱流分布云圖

熱流的傳遞由TPS、MQ和緊支C2基函數(shù)6種插值方法(包含原始和梯度修正插值方案)從流體域網(wǎng)格到固體域網(wǎng)格進(jìn)行傳遞，其中MQ的形態(tài)參數(shù)c取0.01，緊支基函數(shù)半徑固定為算例外形3個方向最大尺度的0.3倍。同時依據(jù)表達(dá)式計算固體域網(wǎng)格點的表面熱流作為標(biāo)準(zhǔn)值來衡量插值的精度，誤差的計算公式為

(16)

圖5給出了插值后固體域網(wǎng)格的熱流分布云圖。從圖中可以看出，6種插值方法均能較好地捕捉到流場表面的熱流分布特征，但是從舵面前緣等流動變化劇烈位置來看，原始的3種插值方法會形成間斷性的不規(guī)則波紋；與之對比梯度修正后的3種插值方法在舵面前緣部位均未出現(xiàn)波紋，對于流場局部細(xì)節(jié)的捕捉要優(yōu)于原始插值方法。其原因主要在于在舵面前緣位置，幾何模型曲率較大，且物理量沿z方向梯度較大，導(dǎo)致舵面前緣流場網(wǎng)格呈現(xiàn)長條型，即沿y方向網(wǎng)格節(jié)點要比z方向稀疏得多，而固體域網(wǎng)格在兩個方向的節(jié)點分布則相對均勻得多，因此對于這種從稀疏網(wǎng)格到密網(wǎng)格的插值容易產(chǎn)生數(shù)值振蕩。

圖5 舵面固體域網(wǎng)格熱流分布云圖

圖6為6種插值方法的誤差分布云圖，從結(jié)果比較來看原始方法在前緣部位誤差均較大，相對誤差在10-2左右；梯度修正的插值方法整體插值效果要優(yōu)于原始方法，在前緣部位的改進(jìn)尤為明顯。而從舵面的上表面來看其改進(jìn)效果也相當(dāng)明顯，原始方法的上表面相對誤差在10-3左右，梯度修正后的相對誤差則主要集中于10-4左右，下降了約一個量級。

圖6 舵面表面熱流誤差分布云圖

圖7給出了各方法熱流插值過程中的最大誤差和平均誤差統(tǒng)計結(jié)果，很顯然梯度修正的插值方法的平均誤差、最大誤差均比原始的插值誤差要低，其中TPS方法取對數(shù)后的平均誤差從-3.56 降低到-4.63，緊支C2基函數(shù)方法取對數(shù)后的平均誤差從-3.90降低到-4.80。綜合來看本文采用的改進(jìn)方法要優(yōu)于原始的插值方法。

圖7 舵面6種插值方法平均誤差和最大誤差

MQ插值方法和緊支C2基函數(shù)插值方法中均存在形態(tài)參數(shù)，不同的形態(tài)參數(shù)對插值結(jié)果的影響是巨大的。對于不同算例中不合理的參數(shù)甚至可能導(dǎo)致插值結(jié)果誤差過大不能進(jìn)行下一步的計算。本文以MQ中的形態(tài)參數(shù)為例分析了原始插值方法和梯度修正后插值方法中形態(tài)參數(shù)對結(jié)果的影響。圖8給出了形態(tài)參數(shù)c分別為0.001、0.01、0.1時插值結(jié)果誤差對比云圖。從圖中可以發(fā)現(xiàn)原始的MQ方法在形態(tài)參數(shù)較小時誤差較小，隨著形態(tài)參數(shù)的增大，誤差越來越大；在相同的形態(tài)參數(shù)取值下，采用梯度修正能夠很好地降低插值誤差，并且采用梯度修正后的MQ插值方法可以在較大的形態(tài)參數(shù)變化范圍[0.001，0.05]內(nèi)保持很好的精度，因而提高了在實際工程問題中的適用性。

圖8 不同形態(tài)參數(shù)插值誤差對比云圖

2.2 高超聲速翼身組合體模型

本節(jié)采用類X-37運載器[21]翼身組合體外形考察在復(fù)雜外形和流場分布情況下TPS、MQ和緊支C2基函數(shù)3種徑向基函數(shù)在工程實用外形中的實用性和改進(jìn)插值方法帶來的效果。分別采用結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格生成如圖9、圖10兩套不同的網(wǎng)格，其中流體域網(wǎng)格節(jié)點數(shù)13 216,網(wǎng)格單元數(shù)11 702，固體域網(wǎng)格節(jié)點數(shù)5 704，網(wǎng)格單元數(shù)11 088。采用圖9的網(wǎng)格得到了馬赫數(shù)8.84、攻角0°、壁面溫度300 K條件下的數(shù)值氣動熱結(jié)果，如圖11所示。利用3種不同徑向基函數(shù)分別采用梯度修正方法和原始方法完成熱流的傳遞，其中MQ的形態(tài)參數(shù)c取0.001 1，緊支基函數(shù)半徑固定為算例外形3個方向最大尺度的0.25倍。

圖9 翼身組合體流體域網(wǎng)格

圖10 翼身組合體固體域網(wǎng)格

圖11 翼身組合體流體域網(wǎng)格數(shù)值計算熱流分布

圖12給出了采用原始插值方法和梯度修正插值方法得到的熱流插值結(jié)果云圖，從圖12(a)～圖12(d)中可以發(fā)現(xiàn)TPS和MQ方法能比較準(zhǔn)確地捕捉到模型整體的熱流分布特征，但是在機(jī)翼前緣部位原始MQ方法會出現(xiàn)一些不規(guī)則的條紋，而梯度修正的插值方法可以有效解決該類問題。梯度修正后的緊支C2基函數(shù)插值方法可以有效改善機(jī)身頭部和機(jī)身中段處的插值精度，如圖12(e)和圖12(f)所示。從圖中可以發(fā)現(xiàn)，頭部區(qū)域的等值線不太光滑，部分區(qū)域毛刺較多，與原始熱流等值線相比有所偏離，而梯度修正后的緊支C2基函數(shù)插值結(jié)果較為光滑，更接近原有的熱流分布。經(jīng)過分析發(fā)現(xiàn)，頭部區(qū)域及機(jī)翼前緣是熱流梯度較大區(qū)域，從圖9和圖10對比可知流體域的網(wǎng)格在此處的長細(xì)比較大，但是固體域的網(wǎng)格則相對均勻，因此基于物理量梯度修正的插值方法選擇的點更加均勻插值效果也更好。上述結(jié)果表明基于物理量梯度修正的方法可以在頭部及機(jī)翼前緣等流動變化劇烈的位置更精確捕捉到流動細(xì)節(jié)，實現(xiàn)對復(fù)雜物理量分布的精細(xì)刻畫。

圖12 翼身組合體固體域網(wǎng)格表面熱流分布云圖

由于固體域網(wǎng)格并沒有真實的流場結(jié)果，為了對比每一個數(shù)據(jù)點的插值誤差，將固體域上插值得到的熱流再次采用相同的方法插值回流體域。誤差計算公式與2.1節(jié)的公式一致，其中qinterpolation為從固體域插回流場的數(shù)據(jù)，qreal為原始熱流分布。為了觀察改進(jìn)前后插值結(jié)果的誤差對比，本文對誤差分布取以10為底的對數(shù)，圖13為原始插值方法和梯度修正插值方法的誤差分布云圖，從圖中可以直觀地看出，原始TPS方法的插值誤差主要在機(jī)身中段和機(jī)翼前緣，緊支C2基函數(shù)方法的誤差在機(jī)身頭部區(qū)域也較大，而原始MQ方法在整個機(jī)身表面誤差均較大。經(jīng)過梯度修正后的插值方法不僅改善了前緣部位的插值，還大大改善了機(jī)身中段的插值效果。尤其是MQ方法和緊支C2基函數(shù)方法的改進(jìn)效果均很明顯，其中梯度修正后的緊支C2基函數(shù)插值效果甚至比梯度修正后的TPS、MQ方法效果更優(yōu)，即用更少的點達(dá)到了全域基函數(shù)的效果。此外，經(jīng)過測試發(fā)現(xiàn)當(dāng)梯度修正方法的插值半徑為算例外形3個方向最大尺度的0.03倍時，誤差與原始插值方法的0.25倍得到的結(jié)果相近。經(jīng)過統(tǒng)計分析發(fā)現(xiàn)對于原始網(wǎng)格點數(shù)為5 704，待插值網(wǎng)格點數(shù)為13 216的插值算例，梯度修正方法的平均選取點數(shù)約為181個，而插值半徑為0.25倍的原始方法平均選取點數(shù)為1 859個，選點個數(shù)下降90%以上，插值時間由254.8 s下降到26.4 s。因此本文提出的梯度修正方法實現(xiàn)了精度和效率的雙重兼顧。

圖13 翼身組合體流體域網(wǎng)格表面熱流誤差分布云圖

本文進(jìn)一步對比了在復(fù)雜算例中基于物理量梯度修正的緊支C2基函數(shù)方法與文獻(xiàn)[17]中尺度歸一化的緊支C2基函數(shù)方法插值效果。其表達(dá)式為

(17)

其中：

(18)

式中:Lx、Ly、Lz分別代表3個方向的長度大小。其基本原理是減少飛行器本身3個方向尺度差異過大給插值帶來的影響，主要思想是對飛行器外形在插值前進(jìn)行尺度歸一化處理。

圖14和圖15分別為二者的誤差分布云圖，從圖中可知基于梯度修正的緊支C2基函數(shù)方法在頭部區(qū)域的相對誤差集中于-3.5，尺度歸一化的誤差集中于-2.5，說明前者的插值效果要優(yōu)于后者。定量分析二者的相對誤差可知尺度歸一化的插值方法平均相對誤差為-2.756，最大相對誤差為0.181，而本文基于物理量梯度修正的插值方法平均相對誤差為-3.046，最大相對誤差為0.173，二者均比尺度歸一化的數(shù)值要小。因此可知在復(fù)雜外形算例中梯度修正的插值方法要優(yōu)于文獻(xiàn)[17]中提出的尺度歸一化方法。

圖14 梯度修正方法誤差云圖

圖15 尺度歸一化方法誤差云圖

圖16給出了改進(jìn)前后3種插值方法的最大插值誤差和平均插值誤差。整體來看每種梯度修正后的方法均比原始的方法最大誤差要小。此外，梯度修正后方法的平均誤差均要低于原始的插值方法，其中緊支C2基函數(shù)的平均誤差下降幅度最大，改進(jìn)的插值效果最好。

圖16 翼身組合體6種插值方法平均誤差和最大誤差

為了進(jìn)一步定量分析各種插值方法的誤差變化，圖17給出了改進(jìn)前后各種插值方法的誤差累計百分比曲線圖。橫坐標(biāo)-8代表相對誤差區(qū)間為[-8,-7)，縱坐標(biāo)給出了誤差區(qū)間內(nèi)點數(shù)之和占總點數(shù)的百分比，其他依次類推，如原始緊支C2基函數(shù)的第6個點坐標(biāo)(-3,31.794)表示相對誤差小于-2的點數(shù)占總點數(shù)的31.794%。從統(tǒng)計圖來看，梯度修正后的方法曲線上升的速度要大于原始的方法，說明梯度修正后方法的誤差均集中在誤差較小的區(qū)域，如梯度修正后的緊支C2基函數(shù)方法誤差在[-8,-1]區(qū)間內(nèi)的點數(shù)占比達(dá)到了98%。而原始的緊支C2基函數(shù)方法在該區(qū)間內(nèi)的點數(shù)只占到了86.6%，仍有14.4%的點的誤差在區(qū)間[-1, 0]內(nèi)。綜合來看，對于復(fù)雜外形經(jīng)過物理量梯度修正后的插值方法插值精度均有所提升，緊支C2基函數(shù)方法的改進(jìn)效果最好，原始緊支C2基函數(shù)方法誤差高于TPS、MQ方法，但梯度修正后的緊支C2基函數(shù)方法要優(yōu)于TPS、MQ方法，其結(jié)果與圖13云圖中分析的結(jié)果一致。

圖17 翼身組合體6種插值方法的熱流誤差統(tǒng)計

3 基于風(fēng)洞試驗數(shù)據(jù)的改進(jìn)方法耦合算例驗證

基于所提出的改進(jìn)緊支C2基函數(shù)數(shù)據(jù)傳遞方法，開展了針對壓縮拐角外形氣動力/熱/結(jié)構(gòu)多場耦合風(fēng)洞試驗狀態(tài)的耦合計算數(shù)值模擬，以驗證該方法在長時間耦合計算中的適用性和精度。圖18給出了該壓縮拐角風(fēng)洞試驗?zāi)Ｐ统叽鏪22]，其前緣半徑3 mm，長633.40 mm，寬140 mm，第一壓縮面角度為7.34°，第二壓縮面角度為25.34°，形成18°壓縮拐角。

圖18 試驗?zāi)Ｐ徒Y(jié)構(gòu)尺寸及示意圖

該試驗在中國空氣動力研究與發(fā)展中心?0.6 m 風(fēng)洞進(jìn)行。風(fēng)洞試驗狀態(tài)參數(shù)等效為在高度23.5 km，馬赫數(shù)5.5，攻角α=10°巡航飛行60 s[22]。試驗?zāi)Ｐ筒牧线x取高溫合金IC6，具體參數(shù)見《中國航空材料手冊》[23]。該試驗重點關(guān)注并獲取了氣動力/熱/結(jié)構(gòu)耦合狀態(tài)下的飛行外形變化位移數(shù)據(jù)。圖19為通過特制LED光源照射模型并通過相機(jī)獲得的初始狀態(tài)(t=0 s)、3 s時刻、7 s時刻模型變形特性圖。從圖中可以看到，隨著時間的推進(jìn)，壓縮面前體位置出現(xiàn)了明顯的位移，時間越長，前體位移量越明顯。

圖19 典型時刻前景照明圖

本文耦合計算采用FL-CAPTER軟件平臺內(nèi)置的雙向面/體自適應(yīng)耦合計算策略[24]，結(jié)構(gòu)應(yīng)力/應(yīng)變場計算時，在尾部螺栓孔洞位置施加固定支撐邊界條件。圖20、圖21分別為試驗?zāi)Ｐ土鲌鼍W(wǎng)格示意圖和結(jié)構(gòu)網(wǎng)格示意圖。流場計算湍流模型采用SST模型，得到該狀態(tài)下的流場，圖22 所示為表面熱流分布云圖。

圖20 試驗?zāi)Ｐ土鲌鼍W(wǎng)格

圖21 試驗?zāi)Ｐ徒Y(jié)構(gòu)網(wǎng)格

圖22 試驗?zāi)Ｐ捅砻鏌崃鞣植荚茍D

采用梯度修正后的緊支C2基函數(shù)將氣動力/熱/結(jié)構(gòu)耦合算例中的熱流、溫度進(jìn)行插值以驅(qū)動耦合計算的進(jìn)行，并截取z=0 mm平面內(nèi)的熱流分布曲線得到圖23的熱流曲線分布對比圖，從圖中可以發(fā)現(xiàn)即使在前緣和壓縮拐角等熱流高梯度區(qū)，流場表面熱流與插值后得到的結(jié)構(gòu)表面熱流吻合得很好，因此可以將插值后得到的熱流進(jìn)行熱響應(yīng)分析。由于4 s之前流場還未穩(wěn)定，給模型帶來額外的抖動，熱流和變形數(shù)據(jù)均存在較大波動，有效的試驗數(shù)據(jù)是4 s后的數(shù)據(jù)，因此本文通過耦合計算分析將4 s后的計算結(jié)果與試驗結(jié)果對比，得到沿y方向模型的變形量對比圖，如圖24所示。結(jié)果表明多場耦合仿真結(jié)果與試驗結(jié)果吻合良好，驗證了本文采取的數(shù)據(jù)傳遞策略在長時間多頻次耦合計算狀態(tài)下的有效性，因此能夠?qū)⑵溆糜趯嶋H工程中的多場耦合分析。

圖23 熱流曲線分布對比

圖24 變形量計算和試驗對比

4 結(jié) 論

本文提出了一種基于物理量梯度修正的徑向基函數(shù)數(shù)據(jù)傳遞方法，并采用高超聲速舵面、高超聲速翼身組合體模型和氣動力/熱/結(jié)構(gòu)多場耦合算例測試了該方法對TPS、MQ以及緊支C2基函數(shù)插值效果的影響，得到以下結(jié)論：

1)除模型外形特點和網(wǎng)格分布特征外，目標(biāo)物理量的分布特性同樣會對插值結(jié)果產(chǎn)生影響。舵面前緣附近目標(biāo)物理量變化劇烈，會形成較大的梯度，易造成全局方法的數(shù)值振蕩，而基于物理量梯度修正的TPS、MQ和緊支C2基函數(shù)插值方法能很好地消除這一現(xiàn)象。

2)對于存在形態(tài)參數(shù)的MQ插值方法，在相同的形態(tài)參數(shù)取值下，采用梯度修正能夠很好地降低插值誤差，并且采用梯度修正后的MQ插值方法可以在較大的形態(tài)參數(shù)變化范圍內(nèi)保持很好的精度，因而提高了在實際工程問題中的適用性。

3)基于梯度修正的緊支C2基函數(shù)插值效果要高于原始的TPS、MQ等全域基函數(shù)，達(dá)到甚至超越了基于梯度修正的全域基函數(shù)，尤其在翼身組合體等復(fù)雜外形中，改進(jìn)效果更加明顯。而與原始的緊支C2基函數(shù)插值方法相比，插值時平均選取點數(shù)和插值時間下降幅度較大，實現(xiàn)了精度和效率的兼顧。

4)多場耦合算例的計算結(jié)果表明采用基于梯度修正的插值方法獲得的耦合計算位移變形數(shù)據(jù)與試驗數(shù)據(jù)吻合度很好，這說明基于梯度修正的徑向基函數(shù)數(shù)據(jù)傳遞方法能夠很好地滿足氣動力/熱/結(jié)構(gòu)長時間高頻次耦合計算中的數(shù)據(jù)傳遞需求。

綜上，基于物理量梯度修正的徑向基函數(shù)多場耦合數(shù)據(jù)傳遞方法能夠很好地解決實際復(fù)雜外形物理量分布各項異性特征所帶來的問題，是一種具備應(yīng)用前景的多場耦合數(shù)據(jù)傳遞方法。