葉輪機(jī)氣彈耦合方程求解的時(shí)間推進(jìn)方法評(píng)估

曹琳婷，王丁喜，黃秀全

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

0 引言

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

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

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

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

1 氣彈控制方程的求解

1.1 流場(chǎng)控制方程

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

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

對(duì)該方程的詳細(xì)解釋見(jiàn)文獻(xiàn)[11-12]。另外，采用Spalart-Allmaras湍流模型[15]來(lái)封閉方程，此處不再贅述。

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

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

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

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

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

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

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

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

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

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

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

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

2 時(shí)間推進(jìn)方法

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

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

1階顯式歐拉（1EEM）

1階隱式歐拉（1IEM）

1階混合歐拉（1HEM）

4階顯式Adams（4EAM）

4階隱式Adams（4IAM）

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

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

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

4-1顯式龍格-庫(kù)塔方法（MRK）

基于2階向后差分的雙時(shí)間步長(zhǎng)法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 顫振分析

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

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

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

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

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

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

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

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

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

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

圖6 各時(shí)間推進(jìn)方法的最大時(shí)間步長(zhǎng)

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

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

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

4 結(jié)論

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

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

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