多源不确定性条件下气动弹性系统颤振可靠性分析方法

2021-02-07郑宇宁

振动与冲击 2021年3期

郑宇宁

(中国运载火箭技术研究院，空间物理重点实验室，北京 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%，在计算效率上具有明显优势。