葉輪機氣彈耦合方程求解的時間推進方法評估

曹琳婷，王丁喜，黃秀全

（西北工業大學動力與能源學院，西安 710129）

0 引言

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

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

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

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

1 氣彈控制方程的求解

1.1 流場控制方程

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

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

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

1.2 結構控制方程

結構控制方程的推導基于牛頓第2運動定律，其矩陣表達式為

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

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

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

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

式中：N為結構動力學方程的維度，也是葉片振動的自由度；qi為標量。

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

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

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

于是式（2）降階為N個獨立的標量方程

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

2 時間推進方法

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

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

1階顯式歐拉（1EEM）

1階隱式歐拉（1IEM）

1階混合歐拉（1HEM）

4階顯式Adams（4EAM）

4階隱式Adams（4IAM）

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

該方法的本質是結構項用4階隱式Adams格式，氣動力項用4階顯式Adams格式離散的混合格式。

基于Adams的預測校正方法（PCAM）

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

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

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

經典龍格-庫塔方法（Classical Runge-Kutta）的推進格式為

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

圖1 經典龍格-庫塔應用于氣彈方程求解的算法

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

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

3 結果比較與分析

3.1 單自由度彈簧系統

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

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

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

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

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

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

圖2 各時間推進方法數值解的相對誤差與時間步長的關系

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

3.2 顫振分析

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

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

圖3 跨聲速轉子NASA Rotor 67顫振計算網格

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

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

式中：W為積累功；Estrain為振動葉片的應變能。

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

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

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

各時間推進方法的最大時間步長如圖6所示，圖中的縱坐標Δt以Δτ的倍數衡量。從圖中可見，經典龍格-庫塔方法的最大時間步長為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的預測校正方法的最大時間步長均為（3.6~4.0）Δτ，即（1.03~1.14）×10-5s；雙時間步長法的最大時間步長為2Δτ，即5.68×10-6s。

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

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

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

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

4 結論

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

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

（3）流固耦合求解過程中流場部分的求解難度和耗時遠大于結構部分，因此選用時間推進方法時應不吝惜耦合難度和內存占用的缺點，選擇穩定性更好的隱式方法。