基于廣義Richardson外插方法的顫振模擬耦合時間精度研究

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

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

0 引 言

20世紀80年代，美國、歐洲等發達國家和地區已對數值模擬可信度開展了大規模、有組織、有計劃的研究工作[1]，在研究過程中提出了驗證與確認的概念作為數值模擬可信度研究的重要內容。1998年美國航空航天學會（AIAA）發布了CFD驗證與確認研究的指南[2]，提出了該領域的一些基本問題和基本概念。在此基礎上，CFD驗證與確認工作不斷深入，國內外機構通過發布基準模型構建研究和比較平臺，系統性地開展了一系列試驗與計算工作。

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

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

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

1 Richardson外插方法及驗證方法

1.1 Richardson外插方法

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

其中g1、g2等參數與離散過程無關。

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

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

1.2 廣義Richardson外插方法

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

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

由式（4）可得精度p：

其中：

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

外插得到“精確解”：

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

1.3 驗證方法

采用上述結果可繼續進行誤差及不確定度估計。

近似相對誤差：

外插相對誤差：

網格收斂指數[7]：

本文選取安全系數Fs=1.25。

2 氣動彈性分析方法

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

2.1 計算方法

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

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

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

2.2 耦合方法

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

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

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

3 Isogai Wing模型顫振計算

3.1 計算模型及網格

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

圖1 Isogai Wing結構模型[21]Fig. 1 Structural model of Isogai Wing[21]

3.2 顫振計算結果

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

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

圖2 不同時間步長下機翼時域結構響應對比Fig. 2 Comparison of time domain structural response with different time step sizes

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

4 Isogai Wing模型數值精度分析

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

4.1 時間收斂性

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

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

4.2 廣義Richardson外插方法分析

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

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

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

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

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

表3 基于顫振臨界頻率的時間精度分析結果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外插方法分析結果驗證

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

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

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

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

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

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

表6 基于顫振臨界頻率的線性回歸分析結果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]被廣泛應用于顫振計算程序的驗證工作。機翼材質為均勻的桃花心木薄板，1/4弦長處后掠角為45°，機翼截面為NACA65A004對稱翼型。本文計算網格單元數為321萬，結構計算選取前四階模態，模態頻率為試驗所得頻率，計算馬赫數0.499。

5.2 時間精度分析

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

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

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

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

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

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

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

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

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

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

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

6 結 論

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

分析結果表明：

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

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

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

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