黃熙文，陳衍茂，李文龍，劉廣

中山大學(xué)航空航天學(xué)院，廣東 廣州 510006

近幾十年來，已有大量關(guān)于翼型氣動(dòng)彈性響應(yīng)預(yù)測和分析的研究［1-3］。作用在機(jī)翼上的亞音速來流，通常采用Theodorsen理論來進(jìn)行氣動(dòng)力建模［4-5］。在該理論中，升力和俯仰力矩系數(shù)均由卷積積分建立［6］。由于系統(tǒng)方程同時(shí)包含微分和積分算子，因此使用數(shù)值方法直接求解是相當(dāng)困難的［6］。并且從計(jì)算的角度來看，求解積分方程比解微分方程困難得多。

眾所周知，通過Newmark-β法等時(shí)程積分方法［7-8］，可以方便地模擬僅包含微分算子的狀態(tài)空間模型。但倘若同時(shí)存在積分項(xiàng)，此類振動(dòng)方程則很難化簡為狀態(tài)空間方程。對于亞音速來流作用下的氣動(dòng)彈性系統(tǒng)，Lee 等［9］通過引入4 個(gè)輔助變量，將同時(shí)包含卷積積分的系統(tǒng)，化簡為由8 個(gè)常微分方程表示的狀態(tài)空間模型。Trickey［10］則提出了另一種建模方法，通過使用拉普拉斯變換的方式成功處理了卷積，該方法的狀態(tài)空間模型僅包含6 個(gè)常微分方程。此外，Coller 和Chamara［11］也提出了一種僅用6 個(gè)常微分方程表示氣動(dòng)彈性系統(tǒng)的方法。上述三種狀態(tài)空間方程在過去數(shù)十年中都得到了廣泛的應(yīng)用。Lee等建立的模型，被廣泛應(yīng)用于具有結(jié)構(gòu)非線性的氣動(dòng)彈性系統(tǒng)中，包含立方、間隙以及滯回非線性等模型中［12-15］。Trickey等的模型，通常應(yīng)用于帶有控制面的氣動(dòng)彈性系統(tǒng)的顫振抑制問題中［16-18］。Coller等提出的模型多用于研究極限環(huán)顫振和雙顫振的性質(zhì)［19-20］。

上述三種氣動(dòng)彈性系統(tǒng)狀態(tài)空間模型之間是否等價(jià)，對于氣動(dòng)彈性系統(tǒng)的研究至關(guān)重要。但是，目前還沒有解析研究或數(shù)值分析證明它們之間是完全等價(jià)的。在本文中，我們將分別推導(dǎo)翼型氣動(dòng)彈性系統(tǒng)的三種狀態(tài)空間模型，并從數(shù)學(xué)上證明三種化簡模型是等價(jià)的，再通過數(shù)值結(jié)果進(jìn)一步驗(yàn)證三者之間的關(guān)系。

1 翼型氣動(dòng)彈性系統(tǒng)的狀態(tài)空間模型

二元機(jī)翼模型如圖1 所示，其在俯仰方向和沉浮方向上振動(dòng)。定義彈性軸的俯仰角度為α（抬頭為正），沉浮位移為h（向下為正）。彈性軸位于距機(jī)翼中點(diǎn)的ahb處，質(zhì)心位于距彈性軸xαb處。

圖1 俯仰和沉浮方向振動(dòng)的二元機(jī)翼模型［6］Fig.1 Sketch of an airfoil oscillating in pitch(α)direction,and in plunge(h)direction[6]

考慮亞音速來流，機(jī)翼運(yùn)動(dòng)的耦合方程可表示為［6］

根據(jù)Theodorsen空氣動(dòng)力學(xué)模型［1］，在不可壓非定常流作用下，升力CL(t)和力矩CM(t)系數(shù)為

1.1 模型1（Lee等的模型）

為消除積分算子，Lee等［9］引入如下形式的4個(gè)輔助變量

1.2 模型2（Trickey等的模型）

根據(jù)式（2），升力中的積分可展開為

1.3 模型3（Coller和Chamara的模型）

根據(jù)Coller和Chamara的研究［11］，升力中的積分可以表示為

通過簡單的計(jì)算，可以推導(dǎo)出如下常微分方程組

從本質(zhì)上講，這種方法與Coller和Chamara提出的方法是相同的。

2 等價(jià)性證明

式（9）中的最后兩個(gè)方程可表示為

將式（13）代入到式（9）的前4個(gè)方程中可得

將以上4個(gè)方程代入式（14）中可得

其中

另一方面，系統(tǒng)（3）中的前4個(gè)方程可表示為

經(jīng)過計(jì)算，可以發(fā)現(xiàn)B11+B12HX=D11，D12=B12HW，V1=B12J1以及V2=B12J2. 此外，式（16）中的非線性項(xiàng)與式（15）中的完全一致。因此，可以得出結(jié)論：式（15）與式（16）完全相同，系統(tǒng)（9）與系統(tǒng)（3）是等價(jià)的。

對于系統(tǒng)（7），將式（7）中最后兩個(gè)方程寫成為

其中Q=P-1A21. 將式（19）代入Y=PS中，可將系統(tǒng)（7）的前4個(gè)方程寫成為

因此，式（21）與（14）是等價(jià)的，其中A12PT=B12，這說明本質(zhì)上模型2與模型3是等價(jià)的。

3 數(shù)值驗(yàn)證

接下來，我們將對上述三種模型進(jìn)行數(shù)值驗(yàn)證。其中系統(tǒng)參數(shù)設(shè)置為μ= 100，rα= 0.5，ah= -0.5，ζα=ζξ= 0，ωˉ= 0.25，xα= 0.25。根據(jù)Hopf 分岔理論，可以獲得系統(tǒng)的分岔點(diǎn)為Uf= 6.038 5。當(dāng)無量綱來流速度U超過Uf時(shí)，系統(tǒng)的平衡點(diǎn)X= 0將失去穩(wěn)定性，即系統(tǒng)將產(chǎn)生極限環(huán)振蕩［6］。

基于狀態(tài)空間模型，可以通過如龍格-庫塔法等數(shù)值方法，求得系統(tǒng)（3）、（7）、（9）在任意給定初始條件下的數(shù)值解。為了便于對比，使用下標(biāo)1、2和3分別表示系統(tǒng)（3）、（7）和（9）的解?？紤]線性系統(tǒng)，即G(ξ) =ξ，M(α) =α，沉浮方向的位移可表示為ξi(t)(i= 1，2，3)。當(dāng)初始條件都設(shè)置為X(0) =[0.2 0.1 0 0]T時(shí)，三種模型的解完全一致，如圖2 所示。在接下來的所有討論中，初始條件都設(shè)為上述值。此外，圖2 還給出了三個(gè)位移之間的殘差?？梢钥吹?，RK 法獲得的結(jié)果的殘差數(shù)量級在10-6到10-5之間，遠(yuǎn)遠(yuǎn)小于位移解的數(shù)量級10-1。為了進(jìn)一步驗(yàn)證三種模型之間的等價(jià)性，可采用精細(xì)積分法［23-24］（PIM）來求解原系統(tǒng)?？梢钥吹?，精細(xì)積分法的解對應(yīng)的殘差誤差在10-14數(shù)量級。顯然，這種級別的殘差是由計(jì)算機(jī)截?cái)嗾`差造成的。

圖2 龍格-庫塔法（RK）和精細(xì)積分法（PIM）計(jì)算的沉浮方向位移和對應(yīng)的殘差（U/Uf=0.5）Fig.2 Plunge displacements obtained by Runge-Kutta method(RK)and precise integration method(PIM),where U/Uf=0.5

考慮在俯仰方向有立方非線性剛度，即G(ξ) =ξ，M(α) =α+ηα3，其中取η= 80。當(dāng)無量綱風(fēng)速超過Uf時(shí)，氣動(dòng)彈性系統(tǒng)會產(chǎn)生極限環(huán)振蕩，如圖3所示。可以看到，三種模型獲得的結(jié)果完全一致。而當(dāng)U/Uf大約等于2.25時(shí)，系統(tǒng)出現(xiàn)了二次分岔現(xiàn)象。對于這種現(xiàn)象，Liu和Dowell在相關(guān)工作中進(jìn)行了深入的研究［6］。值得注意的是，三種模型的二次分岔結(jié)果也完全相同。此外，無論風(fēng)速U/Uf在二次分岔臨界值之前還是之后，三種模型的相圖都是能夠吻合的，如圖4所示。

圖3 考慮立方俯仰角剛度時(shí)RK法獲得的不同模型的分岔圖Fig.3 Bifurcation charts obtained by RK method as a cubic pitch stiffness is included

圖4 考慮剛度立方非線性，在二次分岔前后的相圖對比Fig.4 Comparison of the phase planes before or after the secondary bifurcation as a cubic pitch stiffness is included

考慮如下所示的分段線性

其中間隙參數(shù)設(shè)置為δ=0.5°. 通過RK 法得到的分岔如圖5 所示?？梢钥吹皆诖蠖鄶?shù)情況下，這三種模型的結(jié)果是相同的，但仍存在一些細(xì)微的差異。例如，當(dāng)U/Uf= 0.31 時(shí)，模型1 與其他兩個(gè)模型的結(jié)果不同。但從圖6中可以看出，模型1的解與模型2的解是對稱的。對于模型2，當(dāng)U/Uf= 0.53時(shí)，以及模型3中U/Uf= 0.4時(shí)也出現(xiàn)了類似的現(xiàn)象。

圖5 在俯仰方向考慮分段線性，由RK法得到的不同模型的分岔圖Fig.5 Bifurcation charts obtained by RK as a piecewise stiffness is included

圖6 考慮分段線性剛度，由RK法得到的不同模型沉浮方向的相圖Fig.6 Phase plane in the plunge direction provided by RK method as a piecewise stiffness is included

事實(shí)上，這種現(xiàn)象是由分段線性系統(tǒng)在數(shù)值積分中，分段線性轉(zhuǎn)換點(diǎn)不準(zhǔn)引起的。早在20世紀(jì)90年代，Conner 等就闡述了這一數(shù)值現(xiàn)象［24］，此現(xiàn)象可以通過預(yù)測-校正的精細(xì)積分法來解決［23］。如圖7 所示，通過精細(xì)積分法求解三種模型的數(shù)值解，分岔圖就完全相同了?？偠灾?，對于含間隙非線性氣動(dòng)彈性系統(tǒng)，這三種模型也是等價(jià)的。

圖7 考慮分段線性剛度時(shí)，精細(xì)積分法得到的不同模型的分岔圖Fig.7 Bifurcation charts provided by PIM as a piecewise stiffness is included

4 結(jié) 論

由于Thordorsen函數(shù)中存在卷積積分，使得求解亞音速流作用下翼型氣動(dòng)彈性系統(tǒng)變得非常麻煩。其原因在于，卷積積分會使氣動(dòng)彈性系統(tǒng)同時(shí)出現(xiàn)積分算子和微分算子，這使得無法通過逐步積分法來直接求解。為了解決這一問題，研究人員陸續(xù)提出了一些化簡方法來處理方程中的積分項(xiàng)。在已有的研究中，有三種廣泛使用的狀態(tài)空間化簡模型。這些模型都可以成功處理系統(tǒng)中的積分項(xiàng)，其中模型1將氣動(dòng)彈性系統(tǒng)表示為8 個(gè)一階常微分方程，其他模型由6 個(gè)常微分方程表示，本文從數(shù)學(xué)上和通過數(shù)值模擬，都證明了這三種狀態(tài)空間模型在本質(zhì)上是等價(jià)的。

此外，數(shù)值結(jié)果表明由于分段非線性的存在，這些模型在某些情況下通過RK法可能得到不同的數(shù)值結(jié)果。分析后發(fā)現(xiàn)，這些差異是由于數(shù)值積分過程中，氣動(dòng)彈性系統(tǒng)轉(zhuǎn)換點(diǎn)不準(zhǔn)確造成的。在本文中，通過將RK 法換成預(yù)估-校正的精細(xì)積分法，精確地確定轉(zhuǎn)換點(diǎn)，從而解決了這個(gè)問題。并且三種模型得到完全相同的結(jié)果。

附錄

參數(shù)ci和di：

式（3）中的矩陣和系數(shù)向量為

式（7）中的矩陣為

式（9）中的矩陣為