鄭宇寧

(中國運(yùn)載火箭技術(shù)研究院，空間物理重點(diǎn)實(shí)驗(yàn)室，北京 100076)

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

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

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

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

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

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

1 氣動彈性系統(tǒng)顫振穩(wěn)定性分析方法

根據(jù)結(jié)構(gòu)有限元數(shù)值計(jì)算原理，可以將多自由度氣動彈性系統(tǒng)的時域運(yùn)動方程表示為如下形式：

(1)

質(zhì)量矩陣M可以為如下形式：

(2)

式中：Me為單元的節(jié)點(diǎn)質(zhì)量矩陣；ρ為單元密度；N為單元形函數(shù)；Ve為單元體積。假設(shè)單元阻尼矩陣為Ce，則結(jié)構(gòu)的阻尼矩陣可表示為：

(3)

結(jié)構(gòu)剛度矩陣記為：

(4)

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

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

(5)

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

(6)

在氣動彈性系統(tǒng)顫振穩(wěn)定性分析中，通常對振幅進(jìn)行線性化假設(shè)，即如果振幅足夠小，則認(rèn)為氣動力響應(yīng)與結(jié)構(gòu)位移響應(yīng)之間存在線性關(guān)系。根據(jù)該假設(shè)，Qa(x(t))為：

(7)

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

(8)

在進(jìn)行顫振穩(wěn)定性分析時，可將式(8)動力學(xué)方程的解x(t)假設(shè)為[23]：

x(t)=x0eλt

(9)

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

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

(10)

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

Au=λBu

(11)

這樣，通過求解式(11)所示的廣義特征值問題對氣動彈性系統(tǒng)顫振穩(wěn)定性進(jìn)行判定。給定不同的氣動和結(jié)構(gòu)參數(shù)，利用式(11)求得2n個特征值λi，表示為：

(12)

式中：ωi為第i階振動頻率；ζi為對應(yīng)的阻尼比。根據(jù)系統(tǒng)穩(wěn)定性判定準(zhǔn)則可知，該氣動彈性系統(tǒng)不發(fā)生顫振的條件是所有特征值的實(shí)部都不大于0，表達(dá)式為：

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

(13)

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

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

(14)

因此，氣動彈性系統(tǒng)不發(fā)生顫振的條件表達(dá)式為：

μ<0

(15)

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

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

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

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

(16)

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

圖1 功能函數(shù)的概率密度函數(shù)

(17)

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

(18)

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

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

(19)

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

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

(20)

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

(21)

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

(22)

對氣動彈性系統(tǒng)的非概率顫振可靠度計(jì)算時，可建立如下的功能函數(shù)：

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

(23)

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

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

(24)

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

(25)

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

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

(26)

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

3.1 問題描述與建模

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

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

μcr-μ(αsto,αin)

(27)

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

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

(28)

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

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

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

(29)

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

(30)

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

3.2 含隨機(jī)-區(qū)間參數(shù)的新型混合可靠性分析方法

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

(31)

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

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

pαsto(αsto)

(32)

數(shù)值解為：

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

pαsto(αsto)

(33)

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

(34)

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

q=1,…,Ntotal

(35)

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

q=1,…,Ntotal

(36)

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

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

(37)

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

(38)

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

(39)

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

0,q=1,…,Ntotal

(40)

由于：

R(αsto,q,αin)

(41)

則式(40)可表示為：

q=1,…,Ntotal

(42)

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

q=1,…,Ntotal

(43)

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

(44)

(45)

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

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

(46)

收斂條件為：

(47)

綜上，含隨機(jī)-區(qū)間多源不確定參數(shù)的氣動彈性系統(tǒng)顫振可靠性的分析流程，如圖3所示。具體為：

圖3 顫振可靠性分析流程

(1)構(gòu)造一組代表性點(diǎn)αsto,q,q=1,…,Ntotal，計(jì)算每個代表性點(diǎn)在子域內(nèi)的概率：

(48)

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

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

(49)

(50)

(51)

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

4 數(shù)值算例

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

圖4 二元機(jī)翼模型圖

考慮結(jié)構(gòu)阻尼和立方非線性剛度環(huán)節(jié)，采用三階活塞理論計(jì)算高超聲速流中的氣動力和氣動力矩，運(yùn)用拉格朗日方法建立系統(tǒng)運(yùn)動微分方程[30]：

(52)

(53)

式中：

(54)

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

(55)

式中：

a4=

(56)

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

(57)

Jacobi矩陣的特征值方程為：

(58)

式中，

(59)

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

表1 不確定參數(shù)的均值和中心值

其中隨機(jī)參數(shù)服從正態(tài)分布，變異系數(shù)取為0.01；區(qū)間參數(shù)的區(qū)間半徑和中心值之間滿足：

(60)

其余參數(shù)視為確定性參數(shù)：

a0=0.8,γ=1,b=1

(61)

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

圖5 V∞=1 231.35 m/s時無量綱沉浮運(yùn)動

圖6 V∞=1 231.35 m/s時無量綱俯仰運(yùn)動

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

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

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

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

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

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

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

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

圖11 η隨來流速度的變化

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

表2 混合可靠度η的計(jì)算結(jié)果對比

5 結(jié) 論

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

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

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