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

王 昊，王運(yùn)濤，孟德虹，王 毅

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000)

0 引 言

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

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

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

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

1 Richardson外插方法及驗(yàn)證方法

1.1 Richardson外插方法

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

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

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

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

1.2 廣義Richardson外插方法

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

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

由式（4）可得精度p：

其中：

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

外插得到“精確解”：

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

1.3 驗(yàn)證方法

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

近似相對(duì)誤差：

外插相對(duì)誤差：

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

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

2 氣動(dòng)彈性分析方法

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

2.1 計(jì)算方法

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

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

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

2.2 耦合方法

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

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

表1 耦合方法對(duì)比Table 1 Comparison of coupling methods

3 Isogai Wing模型顫振計(jì)算

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

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

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

3.2 顫振計(jì)算結(jié)果

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

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

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

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

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

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

4.1 時(shí)間收斂性

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

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

4.2 廣義Richardson外插方法分析

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

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

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

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

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

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

表4 外插相對(duì)誤差小于1%時(shí)的最大時(shí)間步長(zhǎng)Table 4 Maximum time step for ＜1%

表4 外插相對(duì)誤差小于1%時(shí)的最大時(shí)間步長(zhǎng)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é)果驗(yàn)證

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

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

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

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

圖5 Isogai Wing模型顫振邊界計(jì)算結(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顫振計(jì)算及精度分析

5.1 計(jì)算模型

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

5.2 時(shí)間精度分析

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

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

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

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

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

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

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

圖7 AGARD 445.6 Wing模型顫振邊界計(jì)算結(jié)果時(shí)間精度曲線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

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

6 結(jié) 論

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

分析結(jié)果表明：

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

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

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

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