曹琳婷，王丁喜，黃秀全

（西北工業(yè)大學(xué)動力與能源學(xué)院，西安 710129）

0 引言

流體和固體之間的耦合作用（Fluid Structure Interaction，F(xiàn)SI）普遍存在于航空航天、海洋、建筑和生物等工程研究和應(yīng)用領(lǐng)域。在葉輪機研究領(lǐng)域，流固耦合作用產(chǎn)生的氣動彈性振動問題對結(jié)構(gòu)的潛在破壞性極大，嚴重時會造成重大經(jīng)濟損失。在葉輪機流固耦合誘發(fā)的振動問題中，根據(jù)激勵的來源，可以將振動分為自激振動和強迫響應(yīng)[1]。其中顫振是自激振動的典型代表。為了預(yù)防顫振或強迫振動的發(fā)生，高效率和高精度的流固耦合技術(shù)在航空發(fā)動機設(shè)計中的應(yīng)用不可或缺。

定量的壓氣機氣彈性能分析不可避免地需要求解流動控制方程和結(jié)構(gòu)控制方程，其中時間推進的耦合求解方法在多排錯頻設(shè)計以及葉片動響應(yīng)分析等方面仍有解耦分析不具備的優(yōu)勢。基于計算流體力學(xué)（Computational Fluid Dynamics，CFD）和計算結(jié)構(gòu)動力學(xué)（Computational Structural Dynamics，CSD）耦合的時域數(shù)值分析方法可以分為直接求解原始結(jié)構(gòu)動力學(xué)方程的流固全耦合法[2]和求解降階結(jié)構(gòu)動力學(xué)方程的流固部分耦合法[3]。前者非常復(fù)雜，通常需要消耗大量計算資源，其工程應(yīng)用非常少；后者常用于學(xué)術(shù)研究和工程設(shè)計中，因此有大量的研究著重于提高部分耦合法計算的精度和效率。

流固部分耦合法早在2維葉柵的顫振分析中就得到了應(yīng)用，相似的研究[4]都將結(jié)構(gòu)動力學(xué)方程降階為平動和旋轉(zhuǎn)2個自由度，在這個基礎(chǔ)上采用偽時間迭代使得耦合更加緊密，其流場和結(jié)構(gòu)的真實時間推進都采用顯式龍格-庫塔格式。Sadeghi等[6]為解決3維氣動彈性問題給出了一種基于2階向后差分的雙時間步長法的時間推進方法，并應(yīng)用于機翼顫振的預(yù)測；Bendiksen等[7]較早地提出了將Newmark-β、Wilson-θ等一系列常見于結(jié)構(gòu)動力學(xué)方程求解的單步法用于氣動彈性方程的時域推進，同時也提到關(guān)于Adams及其他線性多步法，Adams線性多步法是現(xiàn)代數(shù)值計算中的一類重要方法，適用于求解包含復(fù)雜函數(shù)的常微分方程；Robinson等[8]將基于Adams-Moulton的預(yù)測-校正方法用于氣動彈性方程的時域推進；蔣躍文等[9]利用標(biāo)準(zhǔn)氣動彈性模型校驗了6種Adams相關(guān)格式時間推進的精度和穩(wěn)定性，并提出了一種流場迭代次數(shù)少并且穩(wěn)定性好的半隱式線性多步法；Dokainish等[10]整理并對比了求解單自由度結(jié)構(gòu)動力學(xué)方程的各類時間推進方法，分析了各方法的計算特性。但在多數(shù)涉及葉輪機流固耦合的時間推進方法研究的公開文獻中，研究者側(cè)重于發(fā)展新的時間推進格式，而對不同時間推進格式的效率對比研究較少。

為填補了上述研究領(lǐng)域的空缺，本文通過對比不同時間推進方法的計算效率，力求找到求解流固耦合氣彈方程的最佳方法。

1 氣彈控制方程的求解

1.1 流場控制方程

流場控制方程采用圓柱坐標(biāo)系下的Unsteady Reynolds-averaged Navier-Stokes（URANS）方程

式中：Q為守恒變量；F、G、H分別為軸向、周向和徑向的對流通量；Qvg,x、Qvg,θ、Qvg,r分別為對應(yīng)3個坐標(biāo)方向的由于網(wǎng)格運動引起的通量；Vx、Vθ、Vr分別為對應(yīng)3個坐標(biāo)方向的粘性通量；S為源項。

對該方程的詳細解釋見文獻[11-12]。另外，采用Spalart-Allmaras湍流模型[15]來封閉方程，此處不再贅述。

1.2 結(jié)構(gòu)控制方程

結(jié)構(gòu)控制方程的推導(dǎo)基于牛頓第2運動定律，其矩陣表達式為

式中：從左至右第1項為慣性力；第2項為摩擦力；第3項為彈性力；第4項為氣動力F；M為質(zhì)量矩陣；C為摩擦系數(shù)矩陣；K為剛度矩陣；x為位移向量。

葉輪機葉片小展弦比和實心的特點決定其振動固有頻率高，因此流固耦合對葉片振動振型的影響小，從而使得在流固耦合數(shù)值分析中可以忽略其振型由于流固耦合的影響而引起的變化。葉片固有振型/頻率就是其在真空中振動的振型，因此在無氣動力和阻尼力的真空條件下，式（2）可寫為

根據(jù)模態(tài)分析得到的模態(tài)向量為φi（葉片的振動模態(tài)），模態(tài)向量還滿足正交性

葉片振動的位移向量可以由模態(tài)向量φi的線性組合表示

式中：N為結(jié)構(gòu)動力學(xué)方程的維度，也是葉片振動的自由度；qi為標(biāo)量。

將式（7）的位移表達式代入式（2），可得到

再將上述方程左右同乘振型矩陣的共軛轉(zhuǎn)置ΦH得到

再根據(jù)振型的正交特性，上述方程的左端質(zhì)量矩陣、摩擦系數(shù)矩陣和剛度矩陣分別被對角化為

于是式（2）降階為N個獨立的標(biāo)量方程

上述方程的右端項稱作模態(tài)力，通過URANS方程求解得到。此時，氣彈控制方程（2）的數(shù)值求解簡化成了單自由度方程（13）的求解問題。

2 時間推進方法

流固耦合系統(tǒng)的氣彈控制方程中的結(jié)構(gòu)動力學(xué)方程可以用常微分方程(14)的形式表達，其右邊項可以改寫為

方程(15)的時間導(dǎo)數(shù)項可采用不同的離散格式以達到時間推進的目的，本文選取以下幾種經(jīng)典格式進行比較。

1階顯式歐拉（1EEM）

1階隱式歐拉（1IEM）

1階混合歐拉（1HEM）

4階顯式Adams（4EAM）

4階隱式Adams（4IAM）

4階半隱式Adams[18]（4SAM）

該方法的本質(zhì)是結(jié)構(gòu)項用4階隱式Adams格式，氣動力項用4階顯式Adams格式離散的混合格式。

基于Adams的預(yù)測校正方法（PCAM）

4-1顯式龍格-庫塔方法（MRK）

基于2階向后差分的雙時間步長法

雙時間步長法引入了偽時間導(dǎo)數(shù)項，并對偽時間導(dǎo)數(shù)項差分得到內(nèi)迭代的形式(25)。為了獲得時間精確解，在每個物理時間步內(nèi)需要進行多次內(nèi)迭代，迭代次數(shù)選取為10～20次[16]。

經(jīng)典龍格-庫塔方法（Classical Runge-Kutta）的推進格式為

采用式（26）的時間推進方法求解氣彈方程中的R(tn,Qn)時，每步需要具備3個時刻點tn、tn+Δt/2、tn+Δt的流場求解信息，傳統(tǒng)的簡化處理方式是將氣動力項在每個時間步內(nèi)凍結(jié)。這種處理方式會導(dǎo)致計算精度大大降低。本文發(fā)展了一種算法可以避免這一問題。經(jīng)典龍格-庫塔應(yīng)用于氣彈方程求解的算法如圖1所示。以原先的Δt/2為時間步長，結(jié)構(gòu)控制方程第n步的求解需要的3個流場時刻點為第n步、第n-1步、第n-2步，流場控制方程第n步的求解需要的結(jié)構(gòu)時刻點為第n-2步，推進格式可以重新寫為

圖1 經(jīng)典龍格-庫塔應(yīng)用于氣彈方程求解的算法

Newmark方法本身可以直接求解2階常微分方程（例如結(jié)構(gòu)動力學(xué)方程(2)），相比其他方法省略了將2階方程轉(zhuǎn)換成1階的步驟。該方法的求解迭代過程為

盡管計算過程復(fù)雜，但在γ=0.5、β=0.25時，Newmark方法具有2階精度以及無條件穩(wěn)定的優(yōu)勢，因此也常用于工程計算。

3 結(jié)果比較與分析

3.1 單自由度彈簧系統(tǒng)

為了比較上述時間推進方法的效率和精度，采用2個算例進行驗證。第1個算例是單自由度有阻尼的彈簧-質(zhì)量系統(tǒng)，其受迫振動可以用2階線性常系數(shù)微分方程描述

該方程和降階結(jié)構(gòu)動力學(xué)方程(13)具有相同的形式，不同僅在于右邊項f。假設(shè)f是簡諧力，A為簡諧力的振幅，ωex為簡諧力的頻率，運動方程(30)具有理論解

式中：C1、C2、P1和P2分別為常系數(shù)。

通過上節(jié)提及的時間推進方法求出數(shù)值解，該模型方程的數(shù)值解與理論解的誤差用2-范數(shù)定義

式中：Qnum和Qana分別為方程的數(shù)值解和理論解。

在每個振動周期分別劃分了10、20、50、100、200、500、1000個時間步來考察時間步長與相對誤差的關(guān)系。在總時長45個振動周期內(nèi)各時間推進方法的相對誤差與時間步長的關(guān)系如圖2所示。橫縱坐標(biāo)均采用對數(shù)形式，但該圖橫坐標(biāo)的刻度標(biāo)注為單個振動周期內(nèi)的時間步數(shù)。

圖2 各時間推進方法數(shù)值解的相對誤差與時間步長的關(guān)系

從圖中可見，顯式方法的穩(wěn)定性較差，如1階顯式歐拉每個振動周期需要100個時間步以上才能收斂；1階混合格式、Newmark以及雙時間步長法同是2階精度方法，其準(zhǔn)確性對時間步長大小非常敏感，減小時間步長使誤差減小的效果顯著；基于4階Adams方法的4種方法都有著極為接近的誤差-時間步長關(guān)聯(lián)曲線，其中隱式Adams在較大時間步長（20～50個時間步）略具優(yōu)勢；4-1龍格-庫塔方法與經(jīng)典龍格-庫塔方法的差距很大，原因是伴隨著方程右項外力f的計算次數(shù)減少，4-1龍格-庫塔的準(zhǔn)確度也大大降低；在考查的所有方法之中，經(jīng)典龍格-庫塔方法的精度和效率最高。

3.2 顫振分析

由于壓氣機結(jié)構(gòu)比單自由度彈簧系統(tǒng)復(fù)雜得多，這些時間推進方法在壓氣機葉片顫振耦合分析中的適用性需要專門研究。另外，區(qū)別于單自由度彈簧系統(tǒng)，時域耦合模擬的流場氣動力是未知的。這需要求解URANS方程和結(jié)構(gòu)動力學(xué)方程時交替迭代，因此會對各項時間推進方法的準(zhǔn)確度產(chǎn)生一定影響。

顫振分析的算例選取跨聲速風(fēng)扇轉(zhuǎn)子NASA Rotor 67[17]，此算例存在相對翔實的試驗數(shù)據(jù)，因此被廣泛用于數(shù)值分析結(jié)果的校驗，也常被用于顫振研究[18-20]。采用單通道計算域，計算網(wǎng)格周向41個網(wǎng)格點，徑向89個網(wǎng)格點，軸向161個網(wǎng)格點。跨聲速轉(zhuǎn)子NASA Rotor 67顫振計算網(wǎng)格如圖3所示。進行顫振分析時選取前2階振型作為考查對象：1階彎曲振型（584.04 Hz）和1階扭轉(zhuǎn)振型（1256.95 Hz）。

圖3 跨聲速轉(zhuǎn)子NASA Rotor 67顫振計算網(wǎng)格

以雙時間步長法為例，當(dāng)時間步長取5.68×10-6s時，計算完成了10000步迭代，等效于葉片的1階模態(tài)振動了33個周期，2階模態(tài)振動了71個周期，其模態(tài)位移如圖4所示。展示了1階和2階模態(tài)的模態(tài)位移隨迭代步數(shù)的演化。根據(jù)峰值的衰減速度可以確定葉片振動的對數(shù)衰減率δ（Log-Dec）

圖4 雙時間步長法在時間步長為5.68×10-6 s時迭代10000步的模態(tài)位移

式中：W為積累功；Estrain為振動葉片的應(yīng)變能。

由于模態(tài)位移相鄰2個峰值之間的時間步較少，依據(jù)其計算對數(shù)衰減率時不可避免地帶來較大的數(shù)值誤差，計算得到的對數(shù)衰減率δ出現(xiàn)幅度較大的振蕩。為了減小數(shù)值誤差的影響，計算對數(shù)衰減率δ

本文將δ作為顫振計算時間步長不相關(guān)的考察對象，將ci作為評估時間步長不相關(guān)的指標(biāo)。c1和c2的下標(biāo)分別代表1階和2階模態(tài)。若有這樣的Δt能滿足式(36)，則說明顫振計算滿足時間步長不相關(guān)條件，其中max(Δt)是允許范圍內(nèi)的最大時間步長。在式(36)中，Δτ固定為1個足夠小的時間步長（2.84×10-6s）。在固定時間步長為Δτ時各時間推進方法的對數(shù)衰減率如圖5所示。從圖中可見，在時間步長足夠小時，各方法的對數(shù)衰減率結(jié)果是相吻合的。對所有出現(xiàn)在圖5中的時間推進方法，Δτ對應(yīng)的對數(shù)衰減率計算已經(jīng)達到時間步長不相關(guān)，因此可以作為參考標(biāo)準(zhǔn)。

圖5 在固定時間步長Δτ（2.84×10-6 s）時各時間推進方法的對數(shù)衰減率

各時間推進方法的最大時間步長如圖6所示，圖中的縱坐標(biāo)Δt以Δτ的倍數(shù)衡量。從圖中可見，經(jīng)典龍格-庫塔方法的最大時間步長為6.7倍的Δτ（1.89×10-5s）；其次是4階隱式Adams方法和Newmark方法，最大時間步長為（5.70~6.25）Δτ，即（1.62~1.78）×10-5s；1階混合歐拉、4階 顯式Adams、4階部分隱式Adams以及基于4階Adams的預(yù)測校正方法的最大時間步長均為（3.6~4.0）Δτ，即（1.03~1.14）×10-5s；雙時間步長法的最大時間步長為2Δτ，即5.68×10-6s。

圖6 各時間推進方法的最大時間步長

對于1階顯式歐拉和1階隱式歐拉，前者在時間步長取值很小時（＜Δτ）顫振計算結(jié)果仍發(fā)散；后者在時間步長取值很小時（＜0.4Δτ）計算與時間步長仍相關(guān)。另外，4-1龍格-庫塔方法類似于1階隱式歐拉，時間步長取值在＜0.4Δτ時的顫振計算仍與時間步長相關(guān)。鑒于本文的研究目的是找尋效率最高的時間推進方法，因此上述3種方法均不納入考慮范圍。在所有考慮的時間推進方法中，Newmark方法和基于4階Adams的預(yù)測校正方法增加了全局變量的存儲單元，但憑借如今計算技術(shù)的發(fā)展，足以忽略這些新增存儲單元帶來的負擔(dān)。

分析上述結(jié)果可知，max(Δt)越大，代表計算時需要耗費的迭代次數(shù)越少。因此，經(jīng)典龍格-庫塔方法是所有時間推進方法中最節(jié)省計算時間成本的方法。4階顯式Adams、部分隱式Adams和預(yù)測校正方法的最大時間步長十分接近。在所有納入考慮的方法中，雙時間步長法的最大時間步長是最小的，計算耗時最多。這些結(jié)論也與前述單自由度彈簧系統(tǒng)的相關(guān)數(shù)值計算結(jié)論相符合。

上述結(jié)論和單自由度彈簧系統(tǒng)得出結(jié)論不同之處在于，求解單自由度彈簧系統(tǒng)得出的結(jié)論是Newmark方法略遜于1階混合歐拉格式，但是在顫振耦合計算時Newmark方法表現(xiàn)更佳，并且其優(yōu)越性僅次于經(jīng)典龍格-庫塔的。原因可能在于流固耦合問題的流場方程求解是與結(jié)構(gòu)方程求解交替進行的，1階混合歐拉格式采取顯隱式混合來計算模態(tài)力，而Newmark方法則采取的是隱式，減少了流場求解結(jié)果的數(shù)值損耗。同理也可以解釋4階隱式Adams方法在顫振耦合計算時比其余3種結(jié)果相似的Adams方法更佳。

4 結(jié)論

（1）2項研究的計算結(jié)果有很多共同點，說明時間推進方法的效率更多地與方法本身的數(shù)學(xué)性質(zhì)有關(guān)。因此對于復(fù)雜的流固耦合問題，研究不同時間推進方法的準(zhǔn)確性和耗時成本可以參考簡單模型的。但是流固耦合求解時流場和結(jié)構(gòu)的信息交互可能帶來一定影響，因此還需要校驗；

（2）在所有考查的時間推進方法中，經(jīng)典龍格-庫塔方法的計算效率最高，穩(wěn)定性好；其次是Newmark和4階隱式Adams方法；基于2階向后差分的雙時間步長法效率最低；

（3）流固耦合求解過程中流場部分的求解難度和耗時遠大于結(jié)構(gòu)部分，因此選用時間推進方法時應(yīng)不吝惜耦合難度和內(nèi)存占用的缺點，選擇穩(wěn)定性更好的隱式方法。