秦英泉， 劉祚秋， 劉濟科， 劉廣

1.中山大學航空航天學院，廣東 深圳 518107

2.深圳市智能微小衛(wèi)星星座技術與應用重點實驗室，廣東 深圳 518107

近幾十年中，對含非線性因素的機翼氣動彈性響應進行了大量的研究（楊超等， 2018；楊智春等， 2016）。機翼系統(tǒng)中的非線性因素通?？梢苑譃闅鈩恿徒Y構非線性，氣動力非線性來自于空氣來流（Gupta et al.， 2019），結構非線性來自于結構的大變形（齊念等， 2013）或者是材料的非線性本構關系（Darabi et al.， 2019）等。此外，機翼在加工和安裝過程中的容限誤差，或者是控制面鉸鏈的松動都會引入間隙非線性（Shi et al.， 2023）。理論分析以及風洞試驗都表明，系統(tǒng)存在上述非線性因素時，若來流速度超過臨界速度機翼就會產(chǎn)生極限環(huán)振蕩（Liu et al.， 2018）。

分析非線性氣動彈性系統(tǒng)周期解的方法主要有數(shù)值法和解析法。因數(shù)值方法的局限性，學者更加關注解析或半解析法。常見的的解析或半解析方法，如諧波平衡法（Miguel et al.， 2020）、同倫分析法（Liao et al.， 2004）和增量諧波平衡法（Liu et al.， 2018）等，都能夠成功獲得非線性系統(tǒng)的周期解。但是由于間隙非線性的非光滑性，上述半解析方法在實施之前都需要先對間隙非線性進行光滑化處理（Liu et al.， 2012）。此類光滑化處理或多或少會引入一些額外的模型誤差。且含間隙非線性的氣動彈性系統(tǒng)還屬于自激系統(tǒng)，這意味著對于處于周期振動狀態(tài)的機翼系統(tǒng)，其振動的頻率是未知的（Liu et al.， 2012；Liu et al.， 2018）。這些進一步增大了求解含間隙非線性的氣動彈性系統(tǒng)半解析周期解的難度。

本文提出了一種新的求解非線性系統(tǒng)的半解析方法，即時域最小殘值法（Liu et al.， 2021a；Liu et al.，2021b；Liu et al.，2022）。該方法通過將含間隙非線性的氣動彈性系統(tǒng)的半解析解求解問題轉化為最小值優(yōu)化問題，然后在時域節(jié)點上展開迭代求解。由于該方法屬于時頻混合方法，并且是直接對含間隙非線性的氣動彈性系統(tǒng)的控制方程進行求解，因此無需對非光滑區(qū)域進行光滑處理，得到的半解析解具有更高的精度。

1 二元機翼模型的運動方程

圖1為大展弦比的機翼模型，僅考慮模型在俯仰和沉浮方向的運動（Liu et al.， 2012）。圖中，b為機翼半跨的長度；h為機翼在沉浮方向的位移，取向下為正；α為機翼俯仰角度，取抬頭為正；E為機翼的彈性軸。ahb為彈性軸E到機翼中心的距離，xαb為機翼質心位置G到彈性軸E的距離。

圖1 二元機翼模型Fig.1 Sketch of the two-dimensional airfoil

假設機翼在亞音速的不可壓縮流中運動，根據(jù)第二類Lagrange 方程，建立如下的氣動彈性系統(tǒng)運動方程（Liu et al.，2018）

式中h，α上的點代表對無量綱時間t'的導數(shù)，t'=Vt/b，V為來流速度，t為真實運動時間。無量綱沉浮位移h1=h/b；rα為機翼繞彈性軸E的回轉半徑；ζh和ζα為機翼在沉浮和俯仰方向的阻尼比。無量綱頻率比-ω=ωh/ωα，ωh和ωα為機翼沉浮和俯仰方向的固有頻率。U=V/bωα為無量綱來流速度；μ為單位長度的機翼與空氣密度的比值；m是單位長度的機翼質量。考慮俯仰方向存在間隙非線性，俯仰方向的非線性剛度M(α)為

M(α)與α的關系如圖2 所示。圖中δ為間隙大小，Mf為圖中間隙對應的斜率。

圖2 M(α)與α的關系Fig.2 Sketch of M(α) versus α

根 據(jù)Theodorsen 氣 動 理 論（Theodorsen，1949），方程（1）中的升力CL(t)為

式中ei，fi(i= 0～7)及gi(i= 0～5)的表達式在附錄中給出。引入x=[h，α，y]T，將方程（5）化為

式中非線性向量F(x) =[0，f7M(α)，0]T，系統(tǒng)質量矩陣為

剛度矩陣為

2 時域最小殘值法

2.1 最小值優(yōu)化問題

時域最小殘值法最早是由劉廣等提出，是一種針對強非線性系統(tǒng)的半數(shù)值半解析方法（Liu et al.， 2021a；Liu et al.，2021b；Liu et al.，2022）。本文也將采用時域最小殘值法來求解方程（5）所示的含間隙非線性氣動彈性系統(tǒng)的周期解。對于方程（5）所示的周期解，可以展開為傅里葉級數(shù)

其中xj(t)表示第j個自由度的位移，aj0是常數(shù)，cjk與sjk是待定的正余弦諧波系數(shù)，ω是和來流速度相關的未知頻率。而求解方程（6）的半解析周期解，本質為求解未知參數(shù)

當確定參數(shù)a后，代入方程（7）中，即獲得了系統(tǒng)的半解析周期解。需要注意的是，在實際求解過程中，方程（7）中的級數(shù)是不可能取無窮項的，因此截取方程（7）的前N項作為系統(tǒng)的近似解，即

方程（8）和（9）的近似級數(shù)解代入系統(tǒng)（6）所獲得的殘差，顯然不能在整個時域上恒等于0。此時的系統(tǒng)殘差為

雖然不能使殘差R在整個時域內都為0，但仍可以通過選取合適的參數(shù)a使得殘差R在一個周期t∈[0，T]內盡可能小。自此，關于含間隙非線性氣動彈性系統(tǒng)的半解析周期解求解問題被轉化為最小值優(yōu)化問題，即

式中?(a，t)為非線性目標函數(shù)，A是未知參數(shù)a的可行域。

2.2 增強的響應靈敏度法

對于方程（11）的最小值優(yōu)化問題，可以通過增強響應靈敏度法來迭代求解。即選定一個合適的迭代初值a()0，令

上述正則化也叫增強的正則化，該響應靈敏度法也被稱為增強的響應靈敏度法。

3 數(shù)值算例

對于方程（6）所描述的系統(tǒng)，根據(jù)Hopf分岔理論（Marsden， 1976；Hassard et al.，1981），當來流速度U超過臨界顫振速度Uf= 6.285 1 時，系統(tǒng)將會失去穩(wěn)定性。本算例首先考慮U= 0.8Uf的情況，此時系統(tǒng)將會產(chǎn)生穩(wěn)定的極限環(huán)，其近似周期解可被展開為如方程（8）的傅里葉級數(shù)，以第一個自由度為例，有

對第二和第三個自由度進行類似的操作，并將所得位移、速度和加速度函數(shù)代入方程（6），則含間隙非線性的氣動彈性系統(tǒng)半解析解求解問題就被轉化為方程（11）所示的最小值優(yōu)化問題。即，尋找一組未知參數(shù)

圖3 給出了N= 30，U/Uf= 0.8 時系統(tǒng)的相圖。圖3 中，（a）為沉浮方向的相圖，（b）為俯仰方向的相圖；黑點表示4 階龍格-庫塔法獲得的數(shù)值解，紅線和藍線為時域最小殘值法獲得的半解析解。從圖中可以看出，數(shù)值解和時域最小殘值法獲得的半解析解吻合非常好，這意味著時域最小殘值法獲得的結果有著較高的精度。此外，當U/Uf=0.8 時，系統(tǒng)兩個方向的相圖均不關于原點對稱，這表明方程（8）的半解析周期解同時存在奇次諧波和偶次諧波。表1 給出了N= 30、U/Uf= 0.8 時該半解析周期解的諧波系數(shù)及振動頻率。圖4給出了含間隙的氣動彈性系統(tǒng)的時程圖。和圖3類似，紅線和藍線分別為時域最小殘值法獲得的沉浮方向和俯仰方向的半解析解。時程圖進一步驗證了時域最小殘值法獲得的半解析解和龍格-庫塔法的結果吻合得非常好。

表1 N = 30，U/Uf = 0.8時的諧波系數(shù)與振動頻率Table 1 The Harmonic coefficient and vibration frequency at N = 30，U/Uf = 0.8

圖3 U/Uf = 0.8時系統(tǒng)的相圖Fig.3 The phase of the system at U/Uf = 0.8

圖4 U = 0.8Uf時系統(tǒng)的時程圖Fig.4 The time histories of the system at U = 0.8Uf

圖5給出了不同諧波數(shù)的俯仰方向相圖。圖中紅點為截取項數(shù)N=1 時對應的極限環(huán)，虛線、綠色實線和藍色實線分別表示N= 5，20，30 時俯仰方向的相圖，而黑點則表示4 階龍格-庫塔法的結果。從圖中可以看出，當僅保留1階諧波時，時域最小殘值法雖然能夠收斂，但結果誤差非常大。隨著保留的諧波數(shù)逐漸增加，半解析解的精度也逐漸升高。在圖5 的局部放大圖中，當保留5 階諧波時，時域最小殘值法獲得的半解析解已經(jīng)可以收斂到正確解上，但和精確解仍有區(qū)別。當保留20 或30 階諧波時，半解析解和數(shù)值解完全重合。這意味著當保留的諧波數(shù)增加到一定數(shù)量時，諧波數(shù)的增加對精度的提升已經(jīng)不大，但更多的諧波數(shù)會顯著影響求解的效率。

圖5 保留不同數(shù)目的諧波獲得的俯仰方向相圖Fig.5 The phase of the system with different order of harmonics

為了定量分析截斷的諧波階次對精度的影響，將獲得的半解析解代到控制方程，繪制系統(tǒng)的殘差曲線。圖6 為N= 5，20，30 時兩個自由度的殘差曲線。其中，紅線和黑線分別為沉浮方向和俯仰方向的殘差曲線。從圖中可以看到，N= 5 時系統(tǒng)的殘差為10-5的量級；N= 20時，系統(tǒng)殘差約為N= 5 時的一半量級；N= 30 時，殘差減小到10-6的量級。這說明保留的諧波越多，對應的半解析解有更高的精度；但當保留的諧波數(shù)達到一定的數(shù)值之后，諧波階次對精度的提升有限。

此外，俯仰方向的殘差曲線中有幾個尖峰，這是由于方程（8）所示的半解析解采用傅里葉級數(shù)作為基函數(shù)。對于系統(tǒng)（6）而言，M(α)存在與α有關的不光滑區(qū)域。當采用截斷的傅里葉級數(shù)作為系統(tǒng)的近似解析解時，將會在非光滑區(qū)域的轉折點附近出現(xiàn)無法避免的Gibbs 現(xiàn)象。Gibbs 現(xiàn)象的存在意味著，以增加截斷諧波的階次來提升近似解析解的精度，需要付出巨大的計算代價。即必須保留足夠多的諧波，才能微量地提升半解析解的精度。這一性質是由系統(tǒng)的非光滑屬性本身決定的，和本文采用的時域最小殘值法無關。

對于來流速度U/Uf= 0.39 時的含間隙的氣動彈性系統(tǒng)（6），當系統(tǒng)的未知參數(shù)ai取不同的初始迭代值時，時域最小殘值法可能會獲得不同的結果。圖8給出了以另一組迭代初值獲得的俯仰方向的相圖。對比圖8 和圖7 可以看到，兩者的相圖是反對稱的。這意味著系統(tǒng)的振動形態(tài)和系統(tǒng)的初始狀態(tài)相關，不同的初始狀態(tài)將會導致系統(tǒng)進入不同的穩(wěn)定極限環(huán)，從另一方面體現(xiàn)了含間隙的氣動彈性系統(tǒng)具有豐富的非線性性質。

圖7 U = 0.39Uf時系統(tǒng)的相圖Fig.7 The phase of the system at U = 0.39Uf

圖8 不同初始迭代值時系統(tǒng)的相圖Fig.8 The phase of the system with different initial iteration values

4 結 論

本文提出了一種新的求解含間隙非線性氣動彈性系統(tǒng)的半解析方法，即時域最小殘值法。通過該方法獲得了來流速度為U/Uf= 0.8 和U/Uf=0.39時系統(tǒng)的半解析解，并和數(shù)值法獲得的結果進行了對比，得到以下結論：

（1）時域最小殘值法獲得的半解析解和數(shù)值法的結果吻合的非常好；

（2）即便只考慮1階諧波，時域最小殘值法的迭代也可以收斂。此外保留的諧波越多，獲得的結果精度越高；

（3）時域最小殘值法從不同的初始迭代值出發(fā)，可能收斂到系統(tǒng)不同的極限環(huán)。