多源不確定性條件下氣動彈性系統顫振可靠性分析方法

鄭宇寧

(中國運載火箭技術研究院，空間物理重點實驗室，北京 100076)

近年來，氣動彈性系統的顫振可靠性分析問題引起了學術界和工程界的廣泛關注。氣動彈性系統顫振可靠性分析的目的是對系統可靠度或失效可能度進行定量評估，以保證系統具有足夠的安全水平。

現有的氣動彈性系統顫振可靠性分析模型大多為基于隨機理論，將隨機理論中的可靠性分析方法與氣動彈性系統的可靠性準則相結合，建立相應的隨機可靠性分析方法。Liu等[1]考慮了結構剛度和質量分布的隨機性，研究了機翼顫振的概率可靠性。Pourzeynail等[2]利用對數正態分數描述結構幾何、阻尼等隨機參數，研究了結構在隨機風載下的顫振可靠性，分析了不同隨機參數對結構顫振概率可靠度的影響程度。Ge等[3]在考慮氣動力和氣動導數隨機性的基礎上，采用隨機有限元方法研究了顫振可靠性。周崢等[4]考慮了剛度、阻尼、質量等隨機因素對顫振臨界風速的影響，提出了一種結合可靠性理論與有限元分析的隨機有限元方法，分析了大跨度橋梁的顫振失效概率。Song等[5]考慮頻率、質心參數及質量比等隨機參數，提出了基于改進抽樣方法的概率可靠性分析方法，對二元機翼顫振可靠度和靈敏度進行了高精度數值模擬。李毅等[6-7]基于不確定性定量化原理和可靠性設計理論，利用蒙特卡洛模擬方法和核密度估計法得到給定速度下系統發生顫振的概率，對顫振失效風險進行定量評定。戴玉婷等[8]采用標準蒙特卡洛模擬方法，得到了考慮集中質量不確定性的大展弦比雙梁式機翼顫振速度分布和概率臨界速度，通過隨機方法估計了顫振失效概率。Canor等[9]分別采用隨機攝動法、隨機配點法和伽遼金法，基于隨機特征值分析對大跨度橋梁的概率顫振可靠度進行了計算。Pourazarm等[10]將攝動方法和一階、二階及基于加權平均的概率可靠性模型相結合，對渦輪葉片的顫振失效概率進行了數值預測。Wu等[11]考慮氣動和結構不確定性，提出了一種面向氣動彈性系統的改進型蒙特卡洛模擬分析方法。

然而，傳統隨機可靠性方法對樣本數量的要求較高，當樣本數量匱乏時，則其誤差就很大。為解決貧信息、少數據條件下結構可靠性評估這一難題，Ben-Haim[12]提出了非概率可靠性設計思想。在此基礎上，郭書祥等[13]發展了一種基于區間分析的結構非概率可靠性計算方法，結構安全程度可通過非概率可靠性指標進行度量。曹鴻鈞等[14]建立了一種基于凸模型的非概率可靠性指標，適用于評定凸模型與區間模型共存情況下的結構安全可靠程度。王曉軍等[15]利用對原始數據要求低的區間集合來描述影響結構可靠性的不確定性，構建一種基于體積比思想的非概率集合模型。Jiang等[16-17]考慮了多維變量的相關性，發展了基于多維凸模型的結構可靠性分析方法。何新黨[18]基于重要抽樣思想，建立了區間、橢球、超橢球三種不確定性變量描述下的結構非概率可靠性指標求解方法。Wang等[19]將主動控制理論和區間分析方法相結合，用區間過程模型分析閉環系統不確定響應，提出了一種非概率時變可靠性分析方法，能夠有效預估被控結構的動力可靠度。Li等[20]提出了結構振動主動控制系統的非概率可靠性分析方法，利用區間特征值建立表征閉環控制系統的可靠度指標。

非概率可靠性分析方法已經在結構分析中達到了廣泛共識，但是針對氣動彈性系統的非概率可靠性分析還處于起步階段。Wang等[21]考慮結構參數和來流速度的不確定性，利用區間向量對不確定參數進行表征，將區間攝動方法和p-k法相結合，獲取了顫振臨界速度區間，同時將非概率可靠性指標引入顫振安全性度量中，構建了非概率顫振可靠性評定模型，并應用于機翼顫振的可靠性分析中。蔣國慶[22]將魯棒控制理論與非概率可靠性理論相結合，建立了基于非概率可靠性理論的壁板顫振抑制方法，提高了壁板顫振臨界速度區間。

同時，隨著氣動彈性系統的復雜程度日益加劇，氣動彈性系統不可避免地面臨多源不確定性條件，通常情況是隨機性與非概率不確定性交叉影響，共同作用于氣動彈性系統，影響氣動彈性系統的穩定性。然而，考慮混合不確定因素下的氣動彈性系統可靠性評估研究尚未見報道。

鑒于此，針對包含多源不確定性的氣動彈性系統，即綜合考慮隨機性和非概率但有界不確定性，本文首先建立了以系統狀態矩陣特征值作為穩定性判據的非概率顫振可靠度計算方法，在此基礎上，提出了一種新的混合可靠性分析與度量技術，通過分步求解策略計算非概率可靠度的概率密度函數，獲取多源不確定性條件下的顫振可靠度，利用數值算例驗證了本章所提出方法的可信性和實用性。

1 氣動彈性系統顫振穩定性分析方法

根據結構有限元數值計算原理，可以將多自由度氣動彈性系統的時域運動方程表示為如下形式：

(1)

質量矩陣M可以為如下形式：

(2)

式中：Me為單元的節點質量矩陣；ρ為單元密度；N為單元形函數；Ve為單元體積。假設單元阻尼矩陣為Ce，則結構的阻尼矩陣可表示為：

(3)

結構剛度矩陣記為：

(4)

式中：Ke為單元剛度矩陣；Be為單元幾何矩陣；De為單元彈性矩陣。氣動載荷向量Q(x,t)可以為：

Q(x,t)=Qe(t)+Qa(x(t))

(5)

式中：Qa(x(t))為氣動彈性力，代表與結構位移有關的氣動力；Qe(t)為不含氣動彈性效應的外載荷，如控制面載荷和突風載荷。若忽略非氣動彈性力，則可得到以下方程：

(6)

在氣動彈性系統顫振穩定性分析中，通常對振幅進行線性化假設，即如果振幅足夠小，則認為氣動力響應與結構位移響應之間存在線性關系。根據該假設，Qa(x(t))為：

(7)

式中：CΔQa為氣動阻尼矩陣；KΔQa為氣動剛度矩陣。將式(7)代入式(6)可得：

(8)

在進行顫振穩定性分析時，可將式(8)動力學方程的解x(t)假設為[23]：

x(t)=x0eλt

(9)

將式(9)代入式(8)可得：

[λ2M+λ(C+CΔQa)+(K+KΔQa)]x0=0

(10)

假設M,C,CΔQa,K和KΔQa都為n×n維矩陣，則式(10)可改寫為一個2n階的廣義特征值問題：

Au=λBu

(11)

這樣，通過求解式(11)所示的廣義特征值問題對氣動彈性系統顫振穩定性進行判定。給定不同的氣動和結構參數，利用式(11)求得2n個特征值λi，表示為：

(12)

式中：ωi為第i階振動頻率；ζi為對應的阻尼比。根據系統穩定性判定準則可知，該氣動彈性系統不發生顫振的條件是所有特征值的實部都不大于0，表達式為：

Re[λi(A,B)]<0,i=1,2,…,2n

(13)

式中：Re為特征值λi實部。將所有特征值中的最大實部記為μ，即

μ=max(Re[λi(A,B)]),i=1,2,…,2n

(14)

因此，氣動彈性系統不發生顫振的條件表達式為：

μ<0

(15)

2 單源不確定性條件下顫振可靠性分析方法

2.1 基于概率模型的顫振可靠性分析方法

當僅考慮隨機參數時，用αsto=(αsto,1,αsto,2,…,αsto,m)為已知概率密度函數的一組隨機變量，α的下標符號“sto”含義為“隨機”，則系統功能函數表達式為：

L=L(αsto)=L(αsto,1,αsto,2,…,αsto,m)

(16)

功能函數決定了系統的安全狀態，若既定要求系統均能完成，即L(αsto)>0，則稱系統處于安全狀態；若既定要求系統無法完成，即L(αsto)<0，系統處于不可靠或失效狀態；L(αsto)=0為極限狀態曲面或失效平面，代表系統的臨界失效狀態。功能函數的概率密度函數，如圖1所示。則可靠度R的數值大小為圖1中陰影部分面積，可表示為：

圖1 功能函數的概率密度函數

(17)

L=L(αsto)=L(μcr,μ(αsto))=μcr-μ(αsto)

(18)

根據式(15)可知，μcr=0，則式(18)為：

L(αsto)=-μ(αsto)

(19)

因此，將式(19)代入式(17)中，可將氣動彈性系統的概率顫振可靠度表示為：

R=P(-μ(αsto)>0)=P(μ(αsto)<0)=

(20)

2.2 基于非概率區間模型的顫振可靠性分析方法

(21)

R=P{L(αin)>0}

(22)

對氣動彈性系統的非概率顫振可靠度計算時，可建立如下的功能函數：

L=L(αin)=L(μcr,μ(αin))=μcr-μ(αin)

(23)

根據式(15)可知，μcr=0，則式(23)為：

L(αin)=-μ(αin)

(24)

考慮區間不確定性時，μ(αin)在一定區間內變化，可表示為：

(25)

在氣動彈性系統設計過程中，若不發生顫振失效則系統特征值的最大實部須小于0。但是當不確定性存在時，特征值的最大實部有可能大于0。將以上模型定義為氣動彈性系統的非概率可靠性干涉模型。對于圖2中的情況，其對應的非概率可靠度以下式計算[25]：

圖2 μ(αin)的一維數軸表示

(26)

3 多源不確定性條件下顫振可靠性分析方法

3.1 問題描述與建模

考慮隨機變量與區間變量共同作用于氣動彈性系統的功能函數中，則式(23)可改寫為：

L=L(αsto,αin)=L(μcr,μ(αsto,αin))=

μcr-μ(αsto,αin)

(27)

根據式(15)可知，μcr=0，則式(27)可為：

L(αsto,αin)=-μ(αsto,αin)

(28)

若假定隨機向量αsto為常向量，則在此條件下該混合可靠性模型將退化為非概率可靠性模型；同理，若假定αin為確知，則該混合可靠性模型將退化為隨機可靠性模型。

因此，若首先給定隨機向量一個初始值αsto,0=(αsto,10,αsto,20,…,αsto,m0)，根據“2.2”節建立的非概率顫振可靠性分析方法，可求解對應于該隨機樣本點αsto,0的非概率顫振可靠度R，即：

R=R(-μ(αsto,0,αin)>0)

(29)

這樣，借助于αsto的概率分布特性，利用概率密度演化方法[26-27]可求解R的概率密度函數pR(R)，則隨機-區間多源不確定性條件下顫振可靠性度量指標可以定義為：

(30)

式中：η為多源不確定性條件下的顫振可靠度。

3.2 含隨機-區間參數的新型混合可靠性分析方法

根據上節內容可知，在隨機-區間多源不確定性條件下進行顫振可靠性評估的核心是得到非概率顫振可靠度R的概率密度函數p(R)。因此，基于廣義概率密度演化方程[28]，可建立關于R的概率密度演化方程：

(31)

式中：pRαsto為(R,αsto)的聯合概率密度函數，其初始條件可表示為：

pRαsto(R,αsto,αin,t0)=δ(R-R(αsto,αin,t0))×

pαsto(αsto)

(32)

數值解為：

pRαsto(R,αsto,αin,t)=δ(R-R(αsto,αin,t))×

pαsto(αsto)

(33)

在隨機向量αsto的變化域Ωsto內進行積分，可得到R的概率密度函數pR(R,t)：

(34)

然而，直接求解式(34)獲得pR(R,t)的顯式表達式非常困難，這里將采用數值方法進行近似計算。為便于求解，同樣可在Ωsto內均勻地取Ntotal個樣本點，記為αsto，q(q=1,…,Ntotal)，并且將變化域Ωsto分為Ntotal個子域，記為Ωsto,q(q=1,…,Ntotal)。將式(34)在子域Ωsto,q內積分可得：

q=1,…,Ntotal

(35)

根據積分中值定理可將式(35)化簡為：

q=1,…,Ntotal

(36)

將初始條件式(32)在子域Ωsto,q內積分可得：

δ(R-R(αsto，q,αin,t0))Pq,q=1,…,Ntotal

(37)

引入虛擬時間參數τ[29]，可將R表示為如下形式：

(38)

因此，可將式(31)中顫振可靠度R的概率密度演化方程改寫為：

(39)

同理，式(36)中各子域方程為：

0,q=1,…,Ntotal

(40)

由于：

R(αsto,q,αin)

(41)

則式(40)可表示為：

q=1,…,Ntotal

(42)

同樣地，根據式(37)將初始條件改寫為：

q=1,…,Ntotal

(43)

為了避免計算得到的特征值概率密度函數出現失真，本文采用高階無振蕩(Total Variation Diminishing，TVD)差分格式進行求解。引入通量限制器，對式(42)構造TVD差分格式：

(44)

(45)

φ(r)為通量限制器，采用具有較小耗散的Roe-Sweby通量限制器作為構造通量限制器的基礎，表示為：

φ(r)=max{0,min(2r,1),min(r,2)}

(46)

收斂條件為：

(47)

綜上，含隨機-區間多源不確定參數的氣動彈性系統顫振可靠性的分析流程，如圖3所示。具體為：

圖3 顫振可靠性分析流程

(1)構造一組代表性點αsto,q,q=1,…,Ntotal，計算每個代表性點在子域內的概率：

(48)

(2)針對每一個樣本點αsto,q，將其視為確定性參數，同時考慮區間向量αin，利用“2.2”節提出的非概率顫振可靠度分析方法計算R(αsto,q,αin)；

(3)對初始條件式(43)進行離散，表示為：

(49)

(50)

(51)

計算得到pR(R)后，利用式(30)即可計算隨機-區間多源不確定性條件下的顫振可靠度。

4 數值算例

考慮如圖4所示的超音速流中雙楔形二元機翼，具有俯仰和沉浮兩個自由度。

圖4 二元機翼模型圖

考慮結構阻尼和立方非線性剛度環節，采用三階活塞理論計算高超聲速流中的氣動力和氣動力矩，運用拉格朗日方法建立系統運動微分方程[30]：

(52)

(53)

式中：

(54)

式中：b為翼段半弦長，a0b為剛心與坐標原點的距離，V∞為來流速度，M∞為來流馬赫數，ρ∞為大氣密度，mw為單位展長機翼的質量，h為翼段的沉浮運動位移；θ為翼段繞剛心的俯仰角；Sθ為單位展長機翼對彈性軸的靜矩，Iθ為轉動慣量，Ch和Cθ分別為翼段的沉浮和俯仰阻尼系數，Kh和Kθ分別為翼段的沉浮和俯仰剛度，γ為氣體比熱比。

(55)

式中：

a4=

(56)

系統(55)具有平衡點(0,0,0,0)，本算例研究系統在零平衡點處的穩定性問題，由式(55)可得零平衡點處系統的Jacobi矩陣為[30]：

(57)

Jacobi矩陣的特征值方程為：

(58)

式中，

(59)

考慮顫振分析時存在的多源不確定參數，如表1所示。

表1 不確定參數的均值和中心值

其中隨機參數服從正態分布，變異系數取為0.01；區間參數的區間半徑和中心值之間滿足：

(60)

其余參數視為確定性參數：

a0=0.8,γ=1,b=1

(61)

首先，在確定性條件下，計算得到該系統的顫振臨界速度為1 231.35 m/s。在顫振臨界速度下，系統振動響應如圖5～圖6所示。

圖5 V∞=1 231.35 m/s時無量綱沉浮運動

圖6 V∞=1 231.35 m/s時無量綱俯仰運動

利用本文提出的混合可靠性分析法，首先利用概率密度演化方法得到R的概率密度函數，再計算其均值得到混合可靠度η。分別在來流速度V∞=1 060 m/s,1 070 m/s,1 095 m/s條件下，計算得到了R的概率密度函數曲線，并與樣本點數為105的蒙特卡洛模擬法計進行了比較，計算結果如圖7～圖9所示。

圖7 V∞=1 060 m/s時R對比

圖8 V∞=1 070 m/s時R對比

圖9 V∞=1 095 m/s時R對比

觀察圖中曲線可知，采用本文方法獲取的非概率顫振可靠度R的概率密度函數與蒙特卡洛模擬方法吻合較好。同時，還對比了R的概率密度函數隨來流速度的變化情況，如圖10所示。從圖10可知，隨著來流速度的增大，R的分布范圍逐漸向左移動，即表明多源不確定性條件下的顫振可靠度隨來流速度增加而顯著降低。

圖10 R的概率密度函數隨來流速度的變化

在得到R的概率密度函數的基礎上，利用式(30)可得到混合不確定條件下的顫振可靠度η。

圖11所示為混合可靠度η隨來流速度的變化情況，其中來流速度的變化范圍覆蓋了從小于顫振臨界速度到大于顫振臨界速度的整個過程，其結果表明隨著來流速度增大，混合可靠度η在[0,1]范圍內先迅速減小，再平緩下降趨于0。同時可以看出，當考慮不確定因素時，即使來流速度(如V∞=1 060 m/s)小于顫振臨界速度，系統仍然存在顫振失效風險，而對于確定性系統，當來流速度小于顫振臨界速度時，系統不會發生顫振，這即是不確定性系統與確定性系統的區別所在。

圖11 η隨來流速度的變化

表2所示為不同來流速度對應的混合可靠度η計算結果，對于上述四種不同情況下的可靠度結果比較發現，本文方法結果與蒙特卡洛模擬法結果誤差很小，最大誤差的絕對值不超過2%，且本文方法獲得的可靠性計算結果更加保守。同時，在表2的三個工況下，本文方法的平均計算時間約為蒙特卡洛模擬法的40%，在計算效率上具有明顯優勢。

表2 混合可靠度η的計算結果對比

5 結 論

(1)本文首先以概率和非概率區間模型為基礎，建立了單源不確定性條件下顫振可靠性分析方法。在此基礎上，考慮隨機和區間參數共存的多源不確定環境，對多源不確定性條件下混合可靠性建模方法與求解策略進行了研究，構建了多源不確定性條件下氣動彈性系統的顫振可靠度數值計算方法。

(2)結合工程算例，對比了本文獲取的顫振可靠度與蒙特卡洛模擬法計算結果，同時還比較了不同方法所需要的計算時間。

(3)數值算例結果表明本文方法與蒙特卡洛模擬法的最大計算誤差不超過2%，并且在不降低計算精度前提下，能夠減少約60%的計算時間，從而驗證了本文建立的多源不確定性條件氣動彈性系統顫振可靠性分析方法的有效性和可行性。