王 昊，王運濤，孟德虹，王 毅

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000)

0 引 言

20世紀(jì)80年代，美國、歐洲等發(fā)達國家和地區(qū)已對數(shù)值模擬可信度開展了大規(guī)模、有組織、有計劃的研究工作[1]，在研究過程中提出了驗證與確認(rèn)的概念作為數(shù)值模擬可信度研究的重要內(nèi)容。1998年美國航空航天學(xué)會（AIAA）發(fā)布了CFD驗證與確認(rèn)研究的指南[2]，提出了該領(lǐng)域的一些基本問題和基本概念。在此基礎(chǔ)上，CFD驗證與確認(rèn)工作不斷深入，國內(nèi)外機構(gòu)通過發(fā)布基準(zhǔn)模型構(gòu)建研究和比較平臺，系統(tǒng)性地開展了一系列試驗與計算工作。

可信度研究可以分為驗證與確認(rèn)兩個過程。根據(jù)文獻[2]的定義，驗證過程是評估計算模型的求解精度，確認(rèn)過程則是評估計算模型反映實際物理問題的精度，即本文對計算精度的研究屬于驗證過程。

目前CFD驗證工作主要集中在定常問題，主要內(nèi)容是空間離散以及網(wǎng)格尺度帶來的影響，其中一項重要手段是網(wǎng)格收斂性研究，通過Richardson外插方法來獲得網(wǎng)格無關(guān)解、分析計算空間精度和不確定度。但是對于非定常問題，時間離散帶來的影響同樣需要考慮。時間離散精度目前得到的關(guān)注較少，一般使用算法的理論分析精度來替代實際計算精度，沒有進行相關(guān)的不確定度分析。對于氣動彈性問題的時域求解，涉及到流場和結(jié)構(gòu)的耦合計算，相較于單純CFD問題更為復(fù)雜，其精度不僅取決于CFD以及CSD求解精度還和耦合算法有關(guān)，實際計算精度與理論分析精度難免有所差異。同時作為一種更為精細的氣動彈性問題預(yù)測方法，時域模擬需要耗費大量計算資源，對其結(jié)果進行可信度分析，不僅能夠得到可靠的計算結(jié)果，而且能夠選取更為經(jīng)濟的時間離散步長，提高計算效率。

由于非定常問題的時間離散在數(shù)學(xué)上與空間離散并無本質(zhì)上的差別，因此可以把時間離散看作一維空間離散問題，把空間離散問題的驗證方法加以推廣進行時間離散驗證。本文參照網(wǎng)格收斂性分析方法，對顫振問題的時域模擬進行了時間步長收斂性分析，采用廣義Richardson外插分析計算時間精度并獲得時間步長無關(guān)解，采用網(wǎng)格收斂指數(shù)IGC（Grid Convergence index，GCI）估計時間離散不確定度。建立了顫振問題時域分析時間精度的驗證過程，并檢驗了該方法的可行性。

1 Richardson外插方法及驗證方法

1.1 Richardson外插方法

Richardson外插方法又稱“h2外插法”、“迭代外插法”，于1910年[3]由Richardson首次使用，并于1927年[4]加以完善。該方法使用離散步長h將離散解f與精確解fexact的關(guān)系表示為級數(shù)的形式：

其中g(shù)1、g2等參數(shù)與離散過程無關(guān)。

對于二階精度方法g1= 0，此時只需離散步長分別為h1和h2時的離散解f1、f2即可求得離散無關(guān)解fexact。

根據(jù)離散解f1、f2計算方法的不同，外插fexact為 三階或四階精度。

1.2 廣義Richardson外插方法

采用Roache[5]的方法Richardson外插可推廣至p階方法，可稱之為廣義 Richardson外插方法。此時離散解與精確解之間的關(guān)系可以表示為：

對于非定常問題時間離散，選取時間步長h3＞h2＞h1及其對應(yīng)計算結(jié)果f1、f2、f3，代入公式（3）可得：

由式（4）可得精度p：

其中：

為保證分析結(jié)果準(zhǔn)確性，時間步長選取時需要保證ε32、ε21不能太過接近于0，經(jīng)驗估計r需要大于1.3[6]。

外插得到“精確解”：

在使用Richardson外插方法進行分析的過程中，需要注意所選結(jié)果離散步長不能夠太大，此時計算過程沒有充分收斂，分析過程中忽略的高階項會對分析結(jié)果產(chǎn)生影響，使分析結(jié)果具有很強的不確定性。

1.3 驗證方法

采用上述結(jié)果可繼續(xù)進行誤差及不確定度估計。

近似相對誤差：

外插相對誤差：

網(wǎng)格收斂指數(shù)[7]：

本文選取安全系數(shù)Fs=1.25。

2 氣動彈性分析方法

本文研究工作基于TRIP軟件氣動彈性模塊[8-9]開展。該模塊主要包括流場計算、結(jié)構(gòu)計算、耦合界面插值和動網(wǎng)格四個主要功能部分，依照耦合計算方法的流程依次調(diào)用。

2.1 計算方法

流場求解器采用課題組自主開發(fā)的亞跨超聲速CFD軟件平臺TRIP。經(jīng)過課題組多年的發(fā)展，TRIP已經(jīng)成為一個非常成熟的數(shù)值模擬軟件平臺，大量驗證工作[10-11]表明TRIP軟件具有較高的計算精度和效率。

結(jié)構(gòu)求解模塊采用基于模態(tài)坐標(biāo)的氣動彈性運動方程，可使用龍格庫塔法、線性多步法、雙時間步方法等多種方法求解。為避免CFD/CSD耦合求解中的質(zhì)量不相似問題，本文采用變剛度方法[12-13]確定顫振邊界。

動網(wǎng)格模塊集成了目前工程應(yīng)用中常用的TFI、Delaunay以及RBF插值方法，以及基于這些基礎(chǔ)方法開發(fā)的一些復(fù)合型動態(tài)網(wǎng)格變形方法，如RBF_TFI、RBF_Delaunay等。界面插值模塊集成了無限平板樣條IPS、薄板樣條TPS和徑向基函數(shù)插值RBF三種插值方式。本文二維算例Isogai Wing模型為二元翼型，滿足剛體運動假設(shè)，動網(wǎng)格方法采用剛性旋轉(zhuǎn)法。為簡化計算，該算例中結(jié)構(gòu)求解部分與流場求解采用相同網(wǎng)格，無需使用界面插值技術(shù)。

2.2 耦合方法

目前實際應(yīng)用中的耦合計算方法主要分為松耦合和緊耦合兩類。松耦合方法通過交替求解流場和結(jié)構(gòu)方程進行時間推進，不在單場求解過程中進行流場與結(jié)構(gòu)的數(shù)據(jù)交換。該方法能夠充分利用現(xiàn)有的流場和結(jié)構(gòu)求解器，實現(xiàn)過程簡單，應(yīng)用廣泛。但是松耦合方法為了計算簡便往往造成流場和結(jié)構(gòu)的時間不同步，耦合精度只有一階。為提高松耦合方法計算精度，部分學(xué)者采用預(yù)估校正思想對松耦合方法進行改進，提出了具有二階精度的改進的松耦合方法[14-16]。緊耦合方法則是將流場和結(jié)構(gòu)方程均構(gòu)造為子迭代形式[17-18]，在每個物理時間步內(nèi)，對流場和結(jié)構(gòu)進行多次數(shù)據(jù)交換，來消除流場與結(jié)構(gòu)信息交換不同步的問題，時間精度能夠達到二階。

為充分驗證所建立精度分析方法的可靠性，本文分別選取松耦合方法、改進的松耦合方法和緊耦合方法的結(jié)果進行分析。松耦合方法采用凍結(jié)氣動力的龍格庫塔方法[14]，即在松耦合流程中結(jié)構(gòu)求解使用簡化的龍格庫塔方法。作為對照，使用相同耦合流程，使用雙時間步方法作為結(jié)構(gòu)求解方法，構(gòu)造了使用雙時間步方法的松耦合方法，兩種松耦合方法理論時間精度均為一階。改進的松耦合方法采用二階雜交的線性多步法[15]，理論精度為二階。本文所使用緊耦合方法根據(jù)文獻[17]構(gòu)造，理論精度為二階。同時作為對照，本文將龍格庫塔方法作為結(jié)構(gòu)求解嵌入流場求解子迭代過程中，自行構(gòu)造了一種一階精度的緊耦合方法。本文選取上述多種精度的共五種耦合方法進行計算，并對計算結(jié)果進行精度分析。五種耦合方法的對比如表1所示。

表1 耦合方法對比Table 1 Comparison of coupling methods

3 Isogai Wing模型顫振計算

3.1 計算模型及網(wǎng)格

Isogai Wing[19-21]是由Isogai提出的后掠機翼顫振問題的典型截面，是檢驗氣彈不穩(wěn)定性預(yù)測方法的二維標(biāo)準(zhǔn)算例。翼型截面采用NACA010翼型，屬于二元機翼，翼型結(jié)構(gòu)如圖1所示。圖中來流速度為U∞，機翼半弦長b= 0.5 m，在彈性軸（對應(yīng)于剛心）處固定一個垂直方向的線彈簧以及一個扭轉(zhuǎn)彈簧，彈簧剛度分別為Kn和Kα，彈簧阻尼分別為Ch和Cα。選取機翼后緣一側(cè)為正方向，彈性軸在距翼弦中點（C.L.）ab處，其中a= -2.0，表明彈性軸位置在翼弦中點前方，重心在距彈性軸bxα處，其中xα= 1.8，重心位置在翼弦中點后方。機翼具有繞彈性軸俯仰運動和上下平移的浮沉運動兩個自由度，俯仰運動由轉(zhuǎn)角α表示，前緣向上為正，浮沉運動由彈性軸的上下位移ξ表示，向下為正。機翼無量綱回轉(zhuǎn)半徑rα2= 3.48，浮沉運動固有頻率ωh= 100 rad/s，俯仰運動固有頻率ωα= 100 rad/s，質(zhì)量比μ= 60。計算網(wǎng)格為ECERTA Project網(wǎng)站[22]提供的適用于Euler方程計算的結(jié)構(gòu)網(wǎng)格。

圖1 Isogai Wing結(jié)構(gòu)模型[21]Fig. 1 Structural model of Isogai Wing[21]

3.2 顫振計算結(jié)果

選取Ma= 0.6，顫振臨界狀態(tài)下無量綱來流速度vf= 1.710，不同時間步長下五種耦合方法的結(jié)構(gòu)響應(yīng)對比如圖2所示，橫軸為時間t，縱軸為結(jié)構(gòu)一階廣義位移x1。

由結(jié)果對比可知，時間離散步長對計算結(jié)果具有顯著影響。由圖2（a）可知，時間步長較大時，五種耦合方法的結(jié)構(gòu)響應(yīng)存在較大差異，表明此時不同耦合方法的顫振邊界存在較大差異；隨著時間步長的減?。▓D2（b）），結(jié)果趨向一致，顫振邊界趨向收斂；時間步長進一步減小（圖2（c））五種方法所得結(jié)果基本重合，此時五種方法所得顫振邊界基本一致，即可推斷隨著時間步長減小各耦合方法計算顫振邊界均趨向收斂，時間步長足夠小時五種耦合方法均能夠得到正確結(jié)果。

圖2 不同時間步長下機翼時域結(jié)構(gòu)響應(yīng)對比Fig. 2 Comparison of time domain structural response with different time step sizes

通過圖2可以對五種耦合方法進行初步精度分析：兩種二階精度耦合方法計算所得結(jié)構(gòu)響應(yīng)曲線基本重合，并且不受時間步長的影響，即兩種二階精度方法在計算中表現(xiàn)出的精度基本一致；三種一階精度耦合方法結(jié)構(gòu)響應(yīng)存在一定差異，即實際計算精度存在差異，結(jié)構(gòu)求解采用雙時間步方法的松耦合方法精度略小于凍結(jié)氣動力的龍格庫塔方法、一階精度的緊耦合方法計算精度小于兩種二階精度方法。由于隨時間步長減小結(jié)構(gòu)響應(yīng)曲線趨向收斂的方向不同，這種定性分析方法并不能夠直接對比一階精度的松耦合方法與二階精度方法的計算精度，這也說明了本文所建立的精度分析方法的必要性。

4 Isogai Wing模型數(shù)值精度分析

為系統(tǒng)地研究時間離散對計算結(jié)果的影響，在Ma= 0.6條件下選取一系列時間離散步長計算模型的顫振邊界。顫振邊界的確定采用對數(shù)減幅法，選取顫振臨界速度和顫振臨界頻率作為目標(biāo)變量進行分析。

4.1 時間收斂性

為確定計算結(jié)果的時間收斂性，如圖3所示，對比五種方法在不同時間步長下的顫振臨界速度和顫振頻率。五種方法顫振臨界速度和顫振頻率均隨時間步長的減小而單調(diào)變化，最終收斂于同一點。即五種方法均具有良好的時間收斂性，可以使用Richardson外插方法進行分析。

圖3 不同時間步長下Isogai Wing模型顫振邊界計算結(jié)果Fig. 3 Flutter boundary simulation results of Isogai Wing with different time step sizes

4.2 廣義Richardson外插方法分析

首先選取合適時間步長及其計算結(jié)果。根據(jù)1.2節(jié)分析，所選取時間步長不應(yīng)過大，同時應(yīng)使r大于1.3，另考慮對數(shù)減幅法判斷顫振臨界狀態(tài)帶來的誤差，同樣需要避免選取過小的時間步長，最終選取時間步長0.05 、0.1、0.2 ms，即每周期約800、400、200個時間步。

使用廣義Richardson外插方法分析所選計算結(jié)果實際精度，并計算不確定度?；陬澱衽R界速度分析可得表2。結(jié)果表明，分析所得精度p與理論精度基本相符。兩種松耦合方法和兩種二階精度耦合方法實際精度均略高于理論精度，一階精度緊耦合方法精度略低于理論精度。凍結(jié)氣動力的龍格庫塔方法精度略高于采用雙時間步方法的松耦合方法，改進的松耦合方法和傳統(tǒng)緊耦合方法精度均高于一階精度緊耦合方法，改進的松耦合方法和傳統(tǒng)緊耦合方法精度及誤差基本相同，與3.2節(jié)的定性分析結(jié)果一致，說明該方法得到了可信的分析結(jié)果。五種方法插值所得時間離散無關(guān)解之 間最大誤差為0.15%，在精度允許范圍內(nèi)，五種方法分析得到了一致的“精確解”。三種誤差和的對比則說明相同時間步長下二階精度方法計算誤差更小，具有明顯的精度優(yōu)勢。

表2 基于顫振臨界速度的時間精度分析結(jié)果Table 2 Time accuracy analysis results based on the flutter critical velocity

如表3所示，基于顫振臨界頻率與顫振臨界速度的分析結(jié)果基本一致。但是一階精度緊耦合方法、改進的松耦合方法和傳統(tǒng)緊耦合方法精度與顫振臨界速度分析結(jié)果均有一定差異，表明使用廣義Richardson外插方法進行精度分析時，所選取目標(biāo)量對分析結(jié)果有一定影響。

若計算要求顫振臨界速度誤差小于1%，各耦合方法能夠取的最大時間步長如表4所示。二階精度耦合方法能夠大幅提高計算效率，體現(xiàn)了進行精度分析的重要意義。

表3 基于顫振臨界頻率的時間精度分析結(jié)果Table 3 Time accuracy analysis results based on the flutter critical frequency

表4 外插相對誤差小于1%時的最大時間步長Table 4 Maximum time step for ＜1%

表4 外插相對誤差小于1%時的最大時間步長Table 4 Maximum time step for ＜1%

Coupling method 1st order loosely coupled by RK 2nd order tightly coupled Δtmax /ms 0.16 0.08 0.16 0.5 0.5 Step per period 240 480 240 80 80 1st order loosely coupled by dual 1st order tightly coupled 2nd order loosely coupled

4.3 Richardson外插方法分析結(jié)果驗證

基于廣義Richardson外插方法進行時間精度分析需要選取三個計算時間步長的計算結(jié)果，因此需要考慮特定時間步長是否對分析結(jié)果產(chǎn)生影響。為避免所選時間步長結(jié)果帶來的偶然性，對4.1節(jié)使用的計算結(jié)果進行進一步處理，以確定該方法進行時間精度分析的可行性。

圖4 Isogai Wing模型顫振邊界計算結(jié)果時間精度曲線Fig. 4 Time domain accuracy results for Isogai Wing with flutter boundaries

如1.2節(jié)所述，在時間步長較大時，分析結(jié)果具有很強的不確定性，不適用于采用廣義Richardson外插方法進行精度分析，因此圖4中時間步長大于0.6 ms（每周期約60個時間步）時曲線線性度較差。在時間步長較小時，計算誤差相對較小，使用對數(shù)減幅法判斷顫振邊界引入的誤差對精度分析帶來較大干擾，因此圖4中時間步長小于0.04 ms （每周期約1 000個時間步）時，曲線線性度相對較差。若去除時間步長過大和過小的部分，曲線則具有良好的線性度。為更直觀地說明剩余點擬合曲線良好的線性度，對時間步長在0.04～0.6 ms之間的點進行線性回歸分析。線性回歸所擬合直線如圖5所示，直線斜率k和可決系數(shù)R2在表5和表6中給出。

擬合直線斜率k即耦合計算精度，可決系數(shù)R2則反映了擬合直線與選取數(shù)據(jù)點的接近程度。線性回歸分析所得計算精度與4.2節(jié)分析結(jié)果基本一致，R2非常接近于1，數(shù)據(jù)點與擬合直線相當(dāng)接近，即所選取數(shù)據(jù)點具有良好的線性度。線性回歸結(jié)果表明在恰當(dāng)?shù)臅r間步長區(qū)間內(nèi)，使用廣義Richardson外插方法分析所得結(jié)果與具體時間步長選取無關(guān)，也進一步驗證了本文所提出時間精度分析方法的可行性。

圖5 Isogai Wing模型顫振邊界計算結(jié)果線性回歸擬合直線Fig. 5 Linear regression fitting for Isogai Wing with flutter boundaries

表5 基于顫振臨界速度的線性回歸分析結(jié)果Table 5 Linear regression analysis results based on the flutter critical velocity

表6 基于顫振臨界頻率的線性回歸分析結(jié)果Table 6 Linear regression analysis results based on the flutter critical frequency

5 AGARD 445.6 Wing顫振計算及精度分析

5.1 計算模型

AGARD 445.6 Wing[23-24]被廣泛應(yīng)用于顫振計算程序的驗證工作。機翼材質(zhì)為均勻的桃花心木薄板，1/4弦長處后掠角為45°，機翼截面為NACA65A004對稱翼型。本文計算網(wǎng)格單元數(shù)為321萬，結(jié)構(gòu)計算選取前四階模態(tài)，模態(tài)頻率為試驗所得頻率，計算馬赫數(shù)0.499。

5.2 時間精度分析

選取顫振臨界速度和顫振臨界頻率為分析變量，分析三維顫振標(biāo)模AGARD 445.6的計算精度。由于計算量的限制，相較于二維算例，所選計算時間步長及耦合方法較少。

如圖6所示，該算例具有良好的時間收斂性，可以使用Richardson外插方法進行精度分析。根據(jù)1.2節(jié)及4.2節(jié)所述時間步長選取原則，考慮一階精度與二階精度耦合方法收斂性的差異，一階精度松耦合方法選取時間步長為0.05、0.1、0.2 ms （每周期約900、450、225個時間步），改進的松耦合方法選取時間步長為0.2、0.4、0.8 ms （每周期約225、113、57個時間步）。

采用上述時間步長結(jié)果進行精度分析，結(jié)果如表7、表8所示。改進的松耦合方法基于顫振臨界速度的精度與理論精度差異稍大，其余分析精度與理論精度基本一致。由于三維問題更為復(fù)雜，相較于二維算例分析結(jié)果，兩種耦合方法實際計算精度均有所下降。對于兩個目標(biāo)變量，兩種耦合方法得到的時間離散無關(guān)解誤差分別為0.07%、0.18%，均得到了一致的“精確解”。由以上分析結(jié)果可知，該精度分析方法在該三維算例中得到了合理的分析結(jié)果。

為進一步研究三維算例中時間步長選取對分析結(jié)果的影響。采用4.3節(jié)所述方法，得到如圖7所示精度曲線。由于計算結(jié)果已經(jīng)收斂，計算誤差相對較小，改進的松耦合方法精度曲線在小時間步長處線性度較差。由于精度較低，松耦合方法精度曲線在較大時間步長處線性度較差。兩種耦合方法結(jié)果分別去除較小時間步長或較大時間步長的點進行線性回歸分析，分析結(jié)果如表9、表10所示，擬合直線如圖7虛線所示。

圖6 不同時間步長下AGARD 445.6 Wing模型顫振邊界計算結(jié)果Fig. 6 Flutter boundary simulation results of AGARD 445.6 Wing with different time step sizes

表7 基于顫振臨界速度的時間精度分析結(jié)果Table 7 Time accuracy analysis results based on the flutter critical velocity

表8 基于顫振臨界頻率的時間精度分析結(jié)果Table 8 Time accuracy analysis results based on the flutter critical frequency

圖7 AGARD 445.6 Wing模型顫振邊界計算結(jié)果時間精度曲線Fig. 7 Time domain accuracy results for AGARD 445.6 Wing with flutter boundaries

表9 基于顫振臨界速度的線性回歸分析結(jié)果Table 9 Linear regression analysis results based on the fluttercritical velocity

表10 基于顫振臨界頻率的線性回歸分析結(jié)果Table 10 Linear regression analysis results based on the flutter critical frequency

線性回歸分析所得計算精度k與表7、表8分析結(jié)果基本一致；R2均大于0.99，即數(shù)據(jù)點與擬合直線相當(dāng)接近，具有良好的線性度。線性回歸結(jié)果表明，在恰當(dāng)?shù)臅r間步長區(qū)間內(nèi)，使用廣義Richardson外插方法分析所得結(jié)果與具體時間步長選取無關(guān)，進一步驗證了本文所提出時間精度分析方法針對三維顫振問題的可行性。

6 結(jié) 論

本文參照網(wǎng)格收斂性分析方法，對顫振問題的時域模擬結(jié)果進行了時間收斂性和計算精度分析，建立了基于廣義Richardson外插方法的時間精度分析方法，并對該分析方法進行了驗證。

分析結(jié)果表明：

1）本文建立的顫振問題模擬時間精度的分析方法，相較于傳統(tǒng)的定性分析方法，能夠更為準(zhǔn)確地反映顫振問題時域模擬中的計算精度和誤差，為耦合方法及時間步長的選取提供依據(jù)；

2）相比于一階精度方法，二階精度耦合方法實際計算中表現(xiàn)精度更高，具有明顯的精度優(yōu)勢，能夠減小顫振計算所需時間步長，提高計算效率；

3）基于廣義Richardson外插方法分析所得精度與理論分析基本一致，分析結(jié)果可靠。線性回歸分析證明，在合適的時間步長區(qū)間內(nèi)，分析所得結(jié)果與具體時間步長選取無關(guān)。因此，本文提出的時間精度分析方法能夠應(yīng)用于顫振問題時域模擬的時間精度分析。

致謝：本文得到了課題組洪俊武、孫巖、李偉和楊小川等人的幫助，在此表示衷心的感謝。