基于本征正交分解技術(shù)的高效氣動(dòng)彈性耦合計(jì)算方法

回慶龍, 曹博超

(1. 中航工業(yè)沈陽(yáng)飛機(jī)工業(yè)(集團(tuán))有限公司 試飛站, 遼寧 沈陽(yáng) 110850；2. 復(fù)旦大學(xué) 航空航天系, 上海 200433)

0 引 言

近些年來(lái)，越來(lái)越多的彈性結(jié)構(gòu)被應(yīng)用到不同的工程領(lǐng)域中。而由于彈性結(jié)構(gòu)大變形的存在，孤立的氣動(dòng)模型已經(jīng)不能對(duì)這些結(jié)構(gòu)的氣動(dòng)性能進(jìn)行準(zhǔn)確評(píng)估了，對(duì)于這些結(jié)構(gòu)的氣動(dòng)分析需要通過(guò)引入耦合的氣動(dòng)彈性模型實(shí)現(xiàn)[1-3]。高精度的氣動(dòng)彈性計(jì)算需要通過(guò)耦合計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)方法和計(jì)算結(jié)構(gòu)力學(xué)(Computational Structural Dynamics, CSD)方法實(shí)現(xiàn)。結(jié)合不同計(jì)算精度的需求，可以在耦合CFD-CSD方法中選擇不同的流體或結(jié)構(gòu)求解器，而在流體和結(jié)構(gòu)求解器之間需要采用合適的策略完成計(jì)算結(jié)點(diǎn)間的數(shù)據(jù)傳遞(氣動(dòng)載荷和結(jié)構(gòu)位移的傳遞)。近年來(lái)，已經(jīng)有很多工作投入到耦合CFD-CSD方法的研究中了[4-7]。

但是無(wú)論是采用無(wú)黏還是黏性流動(dòng)的求解器，CFD方法往往要消耗大量的計(jì)算資源和計(jì)算時(shí)間。尤其在解決一些非定常氣動(dòng)彈性問(wèn)題的時(shí)候需要與CSD方法耦合進(jìn)行時(shí)間推進(jìn)，其帶來(lái)的計(jì)算代價(jià)更是十分巨大的。在一些工程問(wèn)題中，例如在飛機(jī)概念設(shè)計(jì)階段的氣動(dòng)彈性計(jì)算，有時(shí)希望通過(guò)一定程度的計(jì)算精度的犧牲來(lái)?yè)Q取更高的計(jì)算效率。在過(guò)去幾十年中，渦格法作為一種方便高效的定?；蚍嵌ǔ鈩?dòng)計(jì)算方法被廣泛應(yīng)用于工程領(lǐng)域[8-11]，有時(shí)也被用來(lái)與結(jié)構(gòu)動(dòng)力學(xué)方法耦合進(jìn)行氣動(dòng)彈性計(jì)算。Wang等[12]應(yīng)用非定常渦格法與非線性梁模型發(fā)展了一套氣動(dòng)彈性計(jì)算方法，并將該方法在一個(gè)高空長(zhǎng)航時(shí)機(jī)翼的氣彈計(jì)算上進(jìn)行了驗(yàn)證。Palacios等[13]對(duì)于不同氣動(dòng)及結(jié)構(gòu)計(jì)算模型組合成的氣彈耦合計(jì)算方法在低速大展弦比彈性機(jī)翼的飛行模擬中的表現(xiàn)進(jìn)行了評(píng)估。他們的研究結(jié)果表明片條模型在機(jī)翼做小幅運(yùn)動(dòng)時(shí)表現(xiàn)還可以，而對(duì)于機(jī)翼大幅運(yùn)動(dòng)的情況則需要通過(guò)渦格法進(jìn)行模擬。Murua等人[14]綜合利用渦格法、結(jié)構(gòu)動(dòng)力學(xué)和飛行動(dòng)力學(xué)方法對(duì)大展弦比飛行器的尾跡與尾翼干擾問(wèn)題進(jìn)行了研究。

雖然渦格法相對(duì)CFD方法是更高效的計(jì)算方法，但要準(zhǔn)確離散物面上的環(huán)量分布往往仍然需要幾百或幾千個(gè)自由度，而對(duì)于需要進(jìn)行耦合迭代(氣彈計(jì)算)或優(yōu)化迭代(優(yōu)化設(shè)計(jì))的問(wèn)題其計(jì)算量也是很大的。而通過(guò)模型降階可以進(jìn)一步提高渦格法的計(jì)算效率，從而滿足各種問(wèn)題的需要。近年來(lái)，本征正交分解(Proper Orthogonal Decomposition, POD)方法被人們廣泛應(yīng)用于降階模型的構(gòu)建。Romanowski[15]建立了一個(gè)針對(duì)二維翼型的降階氣彈計(jì)算模型。他應(yīng)用降階模型對(duì)Euler解算器求解的二維翼型非線性平衡態(tài)附近的線性動(dòng)力學(xué)行為進(jìn)行了模擬。Kim[16]發(fā)展出了一套針對(duì)線性動(dòng)力系統(tǒng)的頻域POD方法，并將他的方法分別應(yīng)用于帶阻尼的彈簧振子系統(tǒng)和不可壓縮流動(dòng)的非定常渦格法的降階。Hall等[17]應(yīng)用POD技術(shù)對(duì)二維非定常小擾動(dòng)流動(dòng)模型建立了降階模型。他們應(yīng)用該方法分別對(duì)單一的翼型和葉柵陣列進(jìn)行了跨音速的氣動(dòng)和氣彈計(jì)算。Pettit和Beran[18]應(yīng)用基于POD方法的降階Euler解算器對(duì)馬赫數(shù)1.2下的一個(gè)振蕩突起周?chē)牧鲌?chǎng)進(jìn)行了計(jì)算。他們發(fā)現(xiàn)基于POD技術(shù)的降階模型可以在保有較高解算精度的同時(shí)大大降低計(jì)算自由度，同時(shí)降階模型也增大了容許計(jì)算時(shí)間步長(zhǎng)并提高了計(jì)算模型對(duì)不同工況計(jì)算的魯棒性。

最近，Cao[19]提出了一種基于POD技術(shù)的降階非定常渦格法。該方法已經(jīng)在俯仰和拍動(dòng)運(yùn)動(dòng)的剛性機(jī)翼的非定常氣動(dòng)計(jì)算上得到了驗(yàn)證。在運(yùn)動(dòng)剛體的計(jì)算中，非定常渦格法中的影響系數(shù)矩陣并不隨時(shí)間改變，因此在整個(gè)計(jì)算過(guò)程中只需要進(jìn)行一次系統(tǒng)降階。因此，降階帶來(lái)的計(jì)算加速效果在運(yùn)動(dòng)剛體的氣動(dòng)計(jì)算中體現(xiàn)得并不明顯，于是有必要在變形體的氣動(dòng)計(jì)算中對(duì)方法進(jìn)行進(jìn)一步的驗(yàn)證。本文利用降階非定常渦格法耦合梁模型構(gòu)建了一個(gè)高效的氣彈計(jì)算模型。模型中通過(guò)引入偽時(shí)間步迭代實(shí)現(xiàn)了氣動(dòng)和結(jié)構(gòu)計(jì)算的緊耦合。之后，通過(guò)針對(duì)一個(gè)具有展向彈性的振蕩薄翼模型的氣彈計(jì)算對(duì)該降階計(jì)算模型進(jìn)行了驗(yàn)證。

1 理論與方法

1.1 氣動(dòng)與結(jié)構(gòu)建模

在本研究中，選用了非定常單層渦格法對(duì)作用于薄翼(零厚度)的瞬時(shí)氣動(dòng)載荷進(jìn)行計(jì)算。在渦格法中，物面被渦環(huán)構(gòu)成的薄層所取代，而每一個(gè)渦環(huán)的環(huán)量則通過(guò)壁面無(wú)穿透條件所構(gòu)成的方程進(jìn)行求解，如下所示：

AΓ=R(1)

這里A表示影響系數(shù)矩陣，可由Biot-Savart旋渦誘導(dǎo)定律求出，Γ表示由壁面渦環(huán)環(huán)量組成的待求解向量，而R則表示由自由來(lái)流的法向分量、尾跡渦環(huán)誘導(dǎo)速度、剛體運(yùn)動(dòng)速度和壁面變形速度所組成的右端和向量。在每個(gè)時(shí)間步上，假設(shè)壁面渦環(huán)以不變的環(huán)量由薄翼尾緣向下游尾跡中脫落。在本研究中，為了加速氣彈耦合計(jì)算，忽略了尾跡面的自誘導(dǎo)卷起及尾跡中的環(huán)量耗散。在每個(gè)時(shí)間步上，環(huán)量未知向量?？捎墒?1)解出，之后通過(guò)如下的非定常Bernoulli方程可計(jì)算出相應(yīng)的瞬時(shí)氣動(dòng)載荷：

這里ΔFn,ij是作用在每個(gè)渦格上的氣動(dòng)載荷的法向分量；U、V是自由來(lái)流和剛體運(yùn)動(dòng)的和速度在兩個(gè)方向上的分量；uW、vW是由尾跡面誘導(dǎo)的速度分量；Δb和Δc分別是單個(gè)渦格的展向和弦向尺寸。

在當(dāng)前的研究中，假設(shè)薄翼只具有展向彈性并可由懸臂梁模型進(jìn)行模化。在懸臂梁模型中應(yīng)用Euler-Bernoulli假設(shè)，則結(jié)構(gòu)動(dòng)力學(xué)方程可以寫(xiě)為：

這里w表示薄翼的橫向變形量，撇號(hào)和頭上的點(diǎn)分別表示空間和時(shí)間導(dǎo)數(shù)。p(x,t)表示作用于薄翼的瞬時(shí)分布載荷，該載荷由式(2)算出的氣動(dòng)載荷和由薄翼運(yùn)動(dòng)引起的慣性載荷兩部分組成。E、I和μ分別是薄翼的楊氏模量，截面面積矩和單位長(zhǎng)度的質(zhì)量。應(yīng)用模態(tài)展開(kāi)，橫向變形量w(x,t)可表示成以下形式：

這里wi(x)表示薄翼的第i階模態(tài)，而ηi(t)是對(duì)應(yīng)的模態(tài)廣義坐標(biāo)。在這項(xiàng)研究里，結(jié)構(gòu)計(jì)算中包含了梁的前四個(gè)模態(tài)。將模態(tài)展開(kāi)式(4)帶入式(3)，薄翼的結(jié)構(gòu)運(yùn)動(dòng)方程的廣義坐標(biāo)形式可寫(xiě)為：

1.2 流固耦合算法

本文通過(guò)雙時(shí)間步迭代技術(shù)實(shí)現(xiàn)了上述氣動(dòng)和結(jié)構(gòu)計(jì)算的緊耦合。在每個(gè)物理時(shí)間步上，我們引入了偽時(shí)間步并相應(yīng)地在氣動(dòng)和結(jié)構(gòu)控制方程(式(1)和式(6))中都加入了偽時(shí)間導(dǎo)數(shù)項(xiàng)：

這里用τ表示偽時(shí)間。當(dāng)τ不斷增加， 逐漸趨近于0，而式(7)和式(8)也就收斂到式(1)和式(6)。在計(jì)算中，分別應(yīng)用一階和二階隱式格式對(duì)偽時(shí)間步和物理時(shí)間步進(jìn)行推進(jìn)。于是，離散化的氣動(dòng)和結(jié)構(gòu)方程可分別寫(xiě)為：

-BQm+1+Pm(10)

這里Δt和Δτ分別是物理時(shí)間步和偽時(shí)間步，n和m分別表示物理時(shí)間層數(shù)和偽時(shí)間層數(shù)。

1.3 基于POD技術(shù)的模型降階

在氣動(dòng)彈性計(jì)算中，需要在每個(gè)偽時(shí)間迭代步上對(duì)式(9)和式(10)進(jìn)行求解。如方程(9)所示，離散化的氣動(dòng)方程是一個(gè)N自由度的線性問(wèn)題(這里N表示渦格法中渦格的數(shù)量，對(duì)于不同問(wèn)題N可以是幾百或幾千)。根據(jù)Cao[19]之前的研究，式(9)中的未知向量?？梢酝ㄟ^(guò)它的前幾階POD基向量線性疊加而很好的近似出來(lái)，如下所示：

這里φi是POD分析中獲得的第i階基向量，而vi是對(duì)應(yīng)的廣義坐標(biāo)，Nm是疊加近似中所包含的POD特征模態(tài)數(shù)。應(yīng)用這種POD模態(tài)疊加，離散化的氣動(dòng)方程式(9)可變?yōu)椋?/p>

這里I表示單位矩陣，Φ表示由POD特征基向量φi組成的特征矩陣，而v則是由對(duì)應(yīng)的特征廣義坐標(biāo)vi組成的向量。于是，通過(guò)這種對(duì)向量Γ的近似，氣動(dòng)方程的求解自由度從原來(lái)的N降到了Nm。通常來(lái)說(shuō)，Nm是一個(gè)遠(yuǎn)小于N的數(shù)。例如，Cao[19]在之前的研究中顯示，Nm=3時(shí)已經(jīng)能夠很好的近似(誤差小于1%)做俯仰或拍動(dòng)運(yùn)動(dòng)的薄翼的全氣動(dòng)模型了，即使在很高(k=1.2)的振蕩減縮頻率下。因此，這種模態(tài)分解過(guò)程可以顯著的降低問(wèn)題的求解自由度并加速計(jì)算。

在本文中的氣彈計(jì)算中，彈性結(jié)構(gòu)是一個(gè)長(zhǎng)細(xì)的薄翼并且可以通過(guò)Euler-Bernoulli梁的前4個(gè)模態(tài)疊加進(jìn)行很好的?；虼穗x散結(jié)構(gòu)方程式(10)的計(jì)算自由度僅為8，于是也就無(wú)需對(duì)結(jié)構(gòu)方程進(jìn)行降階操作了。否則，如果我們面對(duì)的是一個(gè)復(fù)雜結(jié)構(gòu)的氣彈計(jì)算問(wèn)題，而結(jié)構(gòu)被氣動(dòng)載荷所激勵(lì)起的特征模態(tài)數(shù)量又很多，我們當(dāng)然也可以遵循類(lèi)似的步驟來(lái)降低結(jié)構(gòu)方程的計(jì)算自由度。

2 結(jié)果和討論

2.1 二維剛性薄翼上的驗(yàn)證測(cè)試

(a) k=0.5

(b) k=1.0

2.2 氣動(dòng)彈性計(jì)算和POD分析

在氣彈計(jì)算中，矩形薄翼的幾何尺寸為弦長(zhǎng)0.1 m，展長(zhǎng)0.6 m，厚度0.001 m，薄翼的密度為1000 kg/m3，空氣密度為1 kg/m3。所有氣彈計(jì)算中渦格分布均設(shè)置為30(弦向)×30(展向)。計(jì)算中，在薄翼的展向中心線位置施加了一個(gè)豎向的強(qiáng)迫正弦振蕩(幅度為0.5倍弦長(zhǎng))，而薄翼的其它部分則根據(jù)給定的展向剛度做自由變形，如圖2所示。

圖2 計(jì)算獲得的薄翼變形及尾渦面的示意圖Fig.2 Deformed foil and wake surface obtained from current computation

圖3 前10個(gè)POD模態(tài)的特征值分布Fig.3 Eigenvalue distribution for first 10 POD modes

圖4 前10階POD模態(tài)的累積模態(tài)能量Fig.4 Cumulative energy of first 10 POD modes

1st mode 2nd mode 3rd mode 4th mode 5th mode

圖5POD分析中獲得的特征模態(tài)
Fig.5EigenmodesobtainedfromPODanalysis

2.3 降階模型的計(jì)算結(jié)果

在從計(jì)算樣本中提取出POD模態(tài)之后，獲得了一個(gè)如式(12)所示的渦格法降階模型(reduced order model, ROM)。為了驗(yàn)證降階氣彈模型的計(jì)算精度，在一個(gè)不同于所有樣本計(jì)算工況的條件(k=0.8,E=3×1010Pa)下進(jìn)行了降階計(jì)算，結(jié)果如圖6所示。

圖6 在k=0.8, E=3×1010 Pa工況條件下的ROM計(jì)算結(jié)果和全模型計(jì)算結(jié)果比較Fig.6 Comparison between ROM results and full-model results for the case k=0.8,E=3×1010 Pa

在計(jì)算中，分別采用了包含不同POD模態(tài)數(shù)量的降階模型，并分別將計(jì)算結(jié)果與全模型進(jìn)行了對(duì)比。如圖6所示，只包含3個(gè)模態(tài)的降階模型與全模型有一定的誤差，包含5個(gè)模態(tài)的降階模型結(jié)果則要準(zhǔn)確得多，而包含7個(gè)或9個(gè)模態(tài)的降階模型的計(jì)算結(jié)果就與全模型的結(jié)果幾乎重合了。第一個(gè)測(cè)試算例的參數(shù)完全在樣本參數(shù)的范圍之內(nèi)，即：k=0.8，在0.2～1.0范圍內(nèi)；E=3×1010Pa，在5×109～1×1011Pa范圍內(nèi)。為了進(jìn)一步驗(yàn)證降階模型的適用范圍，又進(jìn)行了第二個(gè)算例驗(yàn)證，并且將這個(gè)算例的參數(shù)設(shè)置在樣本參數(shù)范圍之外，結(jié)果如圖7所示。

圖7 在k=1.2, E=5×1011 Pa工況條件下的ROM計(jì)算結(jié)果和全模型計(jì)算結(jié)果比較Fig.7 Comparison between ROM results and full-model results for the case k=1.2, E=5×1011 Pa

從表1中可以看出，在所包含模態(tài)數(shù)從3個(gè)增加到7個(gè)的過(guò)程中，降階模型的定量化誤差迅速下降，而如果再包含更多的特征模態(tài)則降階模型的誤差下降速度就明顯減慢了。對(duì)于兩個(gè)不同測(cè)試算例，當(dāng)降階模型中包含7個(gè)以上模態(tài)時(shí)，降階模型與全模型之間的誤差均低于2%。因此，對(duì)于目前問(wèn)題的氣彈計(jì)算來(lái)說(shuō)，包含7個(gè)特征模態(tài)的降階模型是一個(gè)不錯(cuò)的選擇。另外，十分重要的是，對(duì)于當(dāng)前的計(jì)算來(lái)說(shuō)包含9個(gè)特征模態(tài)的降階模型所消耗的計(jì)算CPU時(shí)間僅僅為全模型計(jì)算所需要時(shí)間的1/10。由此可以看出，應(yīng)用本文所提出的降階計(jì)算策略，可以在損失很小計(jì)算精度的同時(shí)大大提高計(jì)算速度。

表1 包含不同特征模態(tài)數(shù)的降階模型的計(jì)算百分比誤差Table 1 Percentage RMS errors for ROMs with different numbers of modes

3 結(jié) 論

在本文中，通過(guò)應(yīng)用降階的非定常渦格法，建立了一種高效的氣動(dòng)彈性計(jì)算模型。模型中通過(guò)偽時(shí)間迭代方法實(shí)現(xiàn)了氣動(dòng)和結(jié)構(gòu)計(jì)算的緊耦合。為了降低整個(gè)氣彈計(jì)算模型的計(jì)算自由度，將未知環(huán)量向量近似展開(kāi)成前幾階POD基向量的疊加形式。通過(guò)一個(gè)豎向振蕩的彈性薄翼的氣彈計(jì)算算例對(duì)該降階計(jì)算模型進(jìn)行了驗(yàn)證。計(jì)算結(jié)果顯示僅僅包含前9個(gè)POD特征模態(tài)的降階氣彈模型就已經(jīng)能夠很準(zhǔn)確的近似全模型計(jì)算了。另外，在驗(yàn)證計(jì)算中，通過(guò)降階計(jì)算方法實(shí)現(xiàn)了90%的計(jì)算CPU時(shí)間節(jié)省。

本文中的方法與針對(duì)CFD方法的POD降階模型類(lèi)似，都是針對(duì)計(jì)算方法中最終形成的線性求解系統(tǒng)進(jìn)行特征分析，再應(yīng)用特征基向量對(duì)系統(tǒng)進(jìn)行降階近似。CFD方法不只可以對(duì)模型氣動(dòng)力進(jìn)行計(jì)算，還同時(shí)可以解算模型周?chē)牧鲌?chǎng)，其包含的信息量是遠(yuǎn)大于本文中應(yīng)用的非定常渦格法的。但也由于其包含的信息量太大，因此要對(duì)其最終形成的線性系統(tǒng)進(jìn)行有效的降階重構(gòu)，就需要很多階的POD基向量。并且在構(gòu)建POD基向量時(shí)，要解析系統(tǒng)中包含的大量信息，就需要在所關(guān)心的參數(shù)范圍內(nèi)保證一定的樣本采集密度。而相比之下，由于渦格法所做的僅僅是對(duì)滿足約束條件的壁面渦量分布的求解，因此其構(gòu)成的線性求解系統(tǒng)中有效的信息量就小得多了，或者說(shuō)只應(yīng)用很少的正交特征基就可以對(duì)全系統(tǒng)進(jìn)行較準(zhǔn)確的重構(gòu)，從本文的結(jié)果中也可以看出這點(diǎn)。由此看來(lái)，本文提出的方法更適合針對(duì)大量消耗計(jì)算時(shí)間的復(fù)雜耦合問(wèn)題計(jì)算或優(yōu)化設(shè)計(jì)計(jì)算，尤其對(duì)于需要在很大參數(shù)范圍內(nèi)進(jìn)行大量計(jì)算的工程問(wèn)題。