非光滑氣動彈性系統的修正增量諧波平衡法

2019-10-22 06:04:02倪迎鴿

兵器裝備工程學報 2019年9期

倪迎鴿,呂毅,張偉

(1.西安航空學院飛行器學院, 西安 710077； 2.西北工業大學無人機特種技術重點實驗室, 西安 710072)

Transform

非線性現象廣泛存在于航空航天領域，而對于二元翼的非線性氣動彈性的預測尤為重要[1]。因為許多飛機在飛行包線內服役時會持續經歷極限環振蕩(limit cycle oscillation，LCO)，此時會導致結構性能下降[2]。因此，對二元翼進行非線性顫振分析很有意義。

Woolston[3]首次利用諧波平衡法(harmonic balance method,HBM)分析了非線性氣動彈性響應，之后便被廣泛地應用于非線性系統。經典的諧波平衡方法，假設解只包含一個諧波，不能預測高階諧波響應[2]。然而，對于大多數工程問題，例如氣動彈性問題，運用經典的諧波平衡方法獲得的解的精度較差。因此學者們提出了改進的諧波平衡方法[4-8]。對于具有立方非線性的二元翼氣動彈性分析表明，較高的諧波次數可以獲得二次分叉，但是隨著諧波個數的增加，非線性項的推導變得復雜，并且隨著系統階數的增加，需要較多的諧波項滿足解的精度[4]。為了獲得更準確的解，文獻[6-8]均對諧波平衡法進行了改進，但是這些改進均會引起更復雜的推導過程。

由Lau提出的增量諧波平衡方法(incremental harmonic balance,IHB)[9]，結合了增量法和諧波平衡方法，已被廣泛地應用于多種非線性方程[10-13]。例如Liu[10]利用增量諧波平衡法研究了俯仰自由度上具有遲滯非線性的二元翼氣動彈性響應，利用不確定系數法解決了增量諧波平衡迭代中的遲滯非線性項的推導問題。Niu[13]結合諧波平衡法和增量諧波平衡方法獲得了強非線性系統的高階穩態近似解，該方法可以簡化高階非線性項的計算。

以上的研究中，大多數是將諧波平衡方法和增量諧波平衡方法應用于光滑非線性系統，例如平方非線性或者立方非線性。一方面，與非光滑系統相比，光滑非線性容易處理。另一方面，隨著系統自由度的增加，分析代價也在增加，十分不利于問題求解。而非光滑系統，例如可以用非光滑函數表達的間隙非線性和遲滯非線性，廣泛存在于工程領域。目前利用諧波平衡方法和增量諧波平衡方法求解非光滑系統的文獻不多。主要是因為對非光滑函數進行傅里葉級數展開很復雜，同時隨著系統階數的增加，需要較多的諧波系數滿足解的精度，求解代價很高[10]。因此，對具有非光滑系統，需要提出有效且通用的增量諧波平衡法，推廣其在工程領域的應用。

本文以具有非光滑氣動彈性系統為對象，結合快速傅里葉變換，修正了增量諧波平衡方法。研究了非光滑項的傅里葉級數展開方法，并獲得了顯式的表達的形式。分別對遲滯非線性和間隙非線性的氣動彈性系統，采用本文的方法獲得了系統的響應，與Runge-Kutta 法獲得的數值解進行對比，對本文方法的有效性進行了驗證，并討論了諧波個數對解的精度影響以及非線性中的參數對響應幅值的影響。

1 二元翼模型與動力學方程

具有俯仰和沉浮兩個自由度的二元翼的示意圖如圖1，可以參考文獻[14]。沉浮用h表示，向下為正方向。關于彈性軸的俯仰用α表示，往上仰為正。彈性軸距翼型中心的距離為ab，質心離彈性軸的距離為xab。

圖1 二元翼示意圖

為了方便地應用于增量諧波平衡方法，對文獻[15]中的非線性氣動彈性方程進行改寫可得：

(1)

如果非線性來源于控制面處鉸鏈的磨損，那么非線性表現為間隙。如果同時考慮摩擦和間隙，則非線性表現為遲滯[10]，可以用非光滑函數來表達。圖2為示意圖，其數學表達式為可以參考文獻[15]。

圖2 非光滑函數示意圖

2 改進增量諧波平衡

第一步為增量步，引入τ=ωt，對式(1)進行變換，可得：

ω2MX″+ωDX′+KX+F(X)=0

(2)

式(2)中“″”和“′”分別代表對τ的一階導數和二階導數。令(X0，ω0)為式(1)的解，則其鄰近點可以表示為：

X=X0+ΔX，ω=ω0+Δω

(3)

將式(3)代入式(2)，并略去高階小量后可以得到以ΔX和Δω為未知量的增量方程：

(4)

(5)

(6)

其中，

將式(5)和式(6)改寫成矩陣的形式：

X0=SA，ΔX=SΔA

(7)

將式(7)代入式(4)，進行諧波平衡過程，得到線性化的代數方程組：

KmcΔA=R+RmcΔω

(8)

(9)

式中，Cstmp=[1 cosθcos2θ… cos(Nc-1)θsinθsin2θ… sinNsθ]，M為|α(τ)|=δ0和|α(τ)|=δ在區間(0,2π)內解的個數。令θ0=0，θM+1=2π，設|α(τ)|=δ0和|α(τ)|=δ的解為θi(i=1,…,M)，且按升序排列。在式(10)中，sgna和sgnb為α(τ)與δ0，δ在[θ0θ1],[θ1θ2],…,[θMθM+1]內的關系，具體的表達式為

(10)

在上述過程中，增量諧波平衡的假設解依賴于諧波個數和頻率。如果諧波項數太多，則式(8)中的矩陣維數很高，意味著計算代價增加。如果初始頻率距離真實值較遠，會導致求解不收斂或者錯誤的解。因此諧波個數和頻率的確定對增量諧波平衡法獲得的解的精度影響較大。

本文改進的基本思想是利用快速傅里葉變換提取數值解中的主導頻率，確定假設解中的諧波個數和近似解的頻率。然后利用獲得的頻率應用到增量諧波平衡方法中，求解線性化方程，有效地確定近似解，同時提高解的精度。

在諧波平衡過程中，令：

(11)

(12)

其中，ck由快速傅里葉變換確定。ajk,bjk,Δajk和Δbjk為幅值。根據前面的求解步驟，可獲得非光滑系統的解。

3 方法驗證與結果分析

這里運用的二元翼模型和相應的參數可以參考文獻[15]。當保留兩個氣動力模態時，對于線性氣動彈性系統(即Mn(α(t))=0)，無量綱顫振速度Vnon=2.0，該結果與文獻[15]中的結果一致。因此，在后面的分析中，均采用兩個氣動力模態進行計算。由于文獻[15]中的氣動彈性模型是無量綱化的，所以文中沉浮和俯仰響應結果沒有標注單位。

3.1 模型驗證

增量諧波平衡法主要用于求解非線性系統的周期解。因此，為了驗證修正的增量諧波平衡法可以求解非光滑系統的響應。分別通過具有遲滯非線性和間隙非線性的二元翼的氣動彈性模型進行驗證。

對于遲滯非線性，分析中使用的參數：Vnon=1.5，2δ=0.017 4 rad,δ0=1.5δ(即k0=0.2kα)，初始條件為[0 0.052 0 0]。利用Runge-Kutta法求解響應的數值解，通過快速傅里葉變換提取響應的主導頻率，確定修正增量諧波平衡法中近似解的諧波個數和頻率。圖3給出了數值解和頻譜圖。

從圖3中可以看出：沉浮響應的主導頻率為0.452 4，俯仰響應的主導頻率為0.452 4和1.357。顯然1.357為0.452 4的3倍。首先假定解中包含1,cos(τ),sin(τ)，然后假定解中包含1,cos(τ),cos(3τ),sin(τ),sin(3τ)，根據前一節的步驟，獲得系統的響應。圖4和圖5給出了本文獲得的解和數值解，可以看出：當1,cos(τ),sin(τ)參與計算時獲得的響應與數值解存在一定的差異。但是當1,cos(τ),cos(3τ),sin(τ),sin(3τ)參與計算時獲得的響應與數值解吻合較好。

圖3 遲滯非線性系統的數值解和頻譜圖

圖4 遲滯非線性系統的俯仰結果(改進的IHB包含1,cos(τ),sin(τ))

圖5 遲滯非線性系統的俯仰結果(改進的IHB包含1,cos(τ),cos(3τ),sin(τ)sin(3τ))

對于間隙非線性，分析中使用的參數：Vnon=1.1，δ=0.017 4 rad，δ0=δ(即k0=0)。利用Runge-Kutta法求解響應的數值解，通過快速傅里葉變換提取響應的主導頻率，確定修正增量諧波平衡法中近似解的諧波個數和頻率。圖6給出了數值解和頻譜圖。

圖6中響應主要包含0.439 8和1.307兩個頻率分量，并且1.307是0.439 8的3倍。分別利用1,cos(τ),sin(τ)和1,cos(τ),cos(3τ),sin(τ),sin(3τ)參與計算。注意到當響應中包含較高階的諧波時，利用本文方法獲得的線性化代數方程組比利用原始的增量諧波平衡法獲得的方程組個數少。例如，用改進的增量諧波平衡法中的xj0=1+aj1cos(τ)+aj3cos(3τ)+bj1sin(τ)+bj3sin(3τ)代替原始的增量諧波平衡法中xj0=1+aj1cos(τ)+aj2cos(2τ)+aj3cos(3τ)+bj1sin(τ)+bj2sin(2τ)+bj3sin(3τ)，獲得的線性化代數方程組，即式(8)，利用原始的增量諧波平衡法建立的代數方程組個數為28，而本文方法獲得的方程組個數為20。圖7和圖8給出了本文獲得的解和數值解，可以看出：當1,cos(τ),sin(τ)參與計算時獲得的響應與數值解存在一定的差異。但是1,cos(τ),cos(3τ),sin(τ),sin(3τ)參與計算時獲得的響應與數值解有很好的一致性。