CN109446557A - 一种基于概率密度演化的随机气动弹性系统稳定性分析方法 - Google Patents

Abstract

本发明公开了一种基于概率密度演化的随机气动弹性系统稳定性分析方法。充分考虑气动弹性系统中存在的不确定因素，利用随机方法对不确定参数进行定量化表征。建立气动弹性系统的有限元方程，并将其转化为广义特征值方程。在此基础上，建立气动弹性系统稳定性分析的概率密度演化方程，通过引入虚拟参数，将概率密度演化方程转化为标准形式，采用有限差分方法和总变差减小格式，求解系统特征值最大实部的概率密度函数，根据最大实部的分布范围进行稳定性分析。数值结果表明，本发明方法得到的气动弹性系统特征值最大实部的概率密度函数与蒙特卡洛方法吻合较好，并且能够大幅度减小计算时间，为随机气动弹性系统的稳定性分析提供了新思路。

Description

一种基于概率密度演化的随机气动弹性系统稳定性分析方法

技术领域

本发明涉及气动弹性设计领域，特别涉及一种基于概率密度演化的随机气动弹性系统稳定性分析方法。

背景技术

气动弹性力学具有三个显著特点，一是其本质上是一个流固耦合问题，结构在气动力的作用下发生弹性变形，同时结构变形反过来又改变流场的边界；二是涉及的非线性因素多，包括结构和气动方面的非线性因素，如控制面铰链处的间隙非线性和大攻角飞行引起的气动非线性；此外，由于气动弹性系统是一个复杂的多学科耦合系统，涉及气动、结构、热等多个学科，实际气动弹性系统中不可避免地存在不确定因素。

对于实际气动弹性系统而言，不确定性的来源是多种多样的，体现在以下四个方面：(1)模型的不确定性，建模过程中对相关因素进行了简化或忽略次要因素导致所建立的气动弹性分析模型与实际对象之间存在模型误差；(2)材料参数的不确定性，由于制造环境、技术条件、材料的多相特征等因素影响，使工程材料的弹性模量、泊松比、质量密度具有不确定性；(3)几何尺寸的不确定性，由于制造及安装误差使结构几何尺寸如厚度、横截面积等具有不确定性；(4)载荷的不确定性，由于测量条件、外部环境等因素影响，使作用在结构上的气动和热载荷具有不确定性。

对于由材料和几何参数引起的不确定性，一般采用随机方法对不确定因素进行量化表征，利用随机变量描述参数的分布特征。对含随机参数的气动弹性系统进行稳定性分析，可将其转化为广义特征值方程，借助特征值分析实现稳定性判定。现有处理随机特征值分析的途径可分为三大类：(1)蒙特卡洛模拟方法；(2)响应面展开方法；(3)矩方程法。然而，蒙特卡洛模拟方法需要进行大量的样本点分析，耗费计算资源较大；响应面展开方法和矩方程法只能处理小不确定性问题。当需要获取特征值最大实部概率密度函数时，目前还没有有效的分析方法能够解决，在一定程度上限制了气动弹性系统稳定性分析技术的发展。综上所述，亟需发展一种能够快速、准确求解气动弹性系统广义特征值概率密度函数的新方法，以克服传统方法计算时间久、精度低的弊端，从而为稳定性分析提供技术支撑。

发明内容

本发明要解决的技术问题为：针对传统处理随机气动弹性系统稳定性分析方法计算效率低、特征值概率密度函数难以获取等问题，提出一种基于概率密度演化的随机气动弹性系统稳定性分析方法。该方法充分考虑气动弹性系统中存在的不确定因素，利用随机方法对不确定参数进行定量化表征。在确定性条件下，建立气动弹性系统的有限元方程，并将其转化为广义特征值问题。在此基础上，建立气动弹性系统稳定性分析的概率密度演化方程，通过引入虚拟参数将概率密度演化方程转化为标准形式，采用有限差分方法和总变差减小格式，求解系统特征值最大实部的概率密度函数，根据最大实部的分布范围进行稳定性分析。

本发明解决上述技术问题采用的技术方案为：一种基于概率密度演化的随机气动弹性系统稳定性分析方法，包括以下步骤：

步骤(1)、建立气动弹性系统的有限元方程：

式中，M为质量矩阵，C为结构阻尼矩阵，为气动阻尼矩阵，K为结构刚度矩阵，为气动刚度矩阵，x(t)为广义坐标，为广义速度，为广义加速度，t为时间，下角标ΔQ_a表示气动；

步骤(2)、令x(t)＝x₀e^λt，可以将气动弹性系统的有限元方程转化为广义特征值方程：

Au＝λBu (2)

式中，λ为广义特征值；

步骤(3)、在确定性条件下，气动弹性系统特征值最大实部μ可以通过下式得到：

μ＝max(Re[λ(A,B)]) (3)

式中，Re表示特征值实部，λ(A,B)通过步骤(2)中的广义特征值方程得到；

步骤(4)、建立气动弹性系统稳定性分析的概率密度演化方程，表示为如下形式：

式中，α＝(α₁,...,α_s)为s维随机不确定参数，p_μα(μ,α,t)为(μ,α)的联合概率密度函数，

步骤(5)、在不确定参数α的变化域Ω内，均匀地取N_total个样本点，记为α_q(q＝1,...,N_total)，并且将变化域Ω分为N_total个子域，记为Ω_q(q＝1,...,N_total)；

步骤(6)、将方程(4)在子域Ω_q内积分，可以得到：

步骤(7)、通过交换积分和求导次序，可以将方程(5)化简为：

式中，为对应于第q个样本点的概率密度函数；

步骤(8)、引入虚拟参数τ，令代入方程(6)中可以得到：

步骤(9)、确定初始条件为：

式中，δ为狄拉克函数，

步骤(10)、采用有限差分方法和总变差减小格式可以得到如下差分格式：

式中， 为差分网格比，为限流器，和可表示为：

步骤(11)、将在N_total个样本点处计算出求和，可以得到：

步骤(12)、取τ_k＝1，则可得到特征值最大实部μ的的表达式：

步骤(13)、根据p_μ(μ)，若存在μ＞0，则该气动弹性系统具有颤振失效风险。

其中，所述步骤(7)中，交换积分和求导次序的具体方式为：

其中，所述步骤(9)中，可将初始条件离散为：

式中， 为方向的网格尺寸；

其中，所述步骤(10)中，差分格式要满足的收敛条件为：

本发明的有益效果是：

本发明提出了一种基于概率密度演化的随机气动弹性系统稳定性分析方法，可以对随机气动弹性系统的特征值进行分析，获取特征值最大实部的概率密度函数，从而对系统稳定性进行判定。本发明方法得到的特征值实部概率密度函数与蒙特卡洛方法得到的特征值实部概率密度函数吻合较好，并且能够大幅度减小计算时间，为随机气动弹性系统稳定性分析提供了新思路。

附图说明

图1为悬臂翼板模型示意图；

图2为V_∞＝2330m/s时μ的概率密度函数；

图3为V_∞＝2354m/s时μ的概率密度函数；

图4为V_∞＝2360m/s时μ的概率密度函数；

图5为本发明的方法实现流程。

具体实施方式

以下将参照附图，对本发明的设计实例进行详细描述。应当理解，所选实例仅为了说明本发明，而不是限制本发明的保护范围。

(1)以悬臂翼板结构为对象，其几何模型如图1所示；

(2)给定翼板结构单元材料属性参数和来流参数，如表1所示；

表1翼板结构材料属性和来流参数

表1中，E为弹性模量，μ_s为为泊松比，ρ_s为翼板密度，ρ_∞为来流密度；

(3)选取壁板弹性模量E和大气密度ρ_∞作为随机参数，服从正态分布，其均值和标准差如表2所示；

表2机翼结构材料属性参数

(4)将E和ρ_∞的变化区间[θ6σ,θ+6σ]等分为20个子区间，则子区间边界处的E和ρ_∞为：

这样，形成的样本点(E_i,ρ_∞j)共有441个；

(5)取第q个样本点(E_i,ρ_∞j)_q，根据(E_i,ρ_∞j)_q设置材料和来流参数，建立该气动弹性系统的有限元方程，通过广义特征值分析获取对应用于样本点(E_i,ρ_∞j)_q的特征值最大实部μ；

(6)建立针对特征值最大实部μ的概率密度演化方程：

(7)初始条件可以离散为：

式中，

(8)有限差分格式设定为：

式中，

(9)限流器设定为：

(10)利用通过有限差分方法，在满足收敛条件的前提下，可以得到对应于样本点(E_i,ρ_∞j)_q的概率密度函数

(11)重复步骤(5)～(10)，计算所有441个样本点对应的概率密度函数将其求和可以得到：

(12)取τ_k＝1，则能够计算出特征值最大实部的概率密度函数表达式为：

(13)利用蒙特卡洛模拟方法，取N＝10000个服从正态分布的样本点，通过对每个样本点对应特征值最大实部μ的计算，得到μ的概率密度函数；

(14)在来流速度为2330m/s、2354m/s和2360m/s的条件下，利用以上两种方法得到μ的概率密度函数如图2-4所示，图中结果说明本发明方法得到的结果与蒙特卡洛方法结果吻合较好；

(15)根据μ的概率密度函数，分别计算均值(θ)和标准差(σ)，结果如表3所示：

表3 μ的均值和标准差对比

从表3可以看出，两种方法的最大计算误差不超过1％，说明本发明方法精度较好；

(16)两种方法计算总耗时分别为：T_{本发明方法}＝3610s,T_{蒙特卡洛模拟方法}＝12120s。时间对比结果表明本发明方法可以减小计算耗时，从而显著提高了机翼结构特征值概率密度函数的计算效率；

(17)根据图2-4可以进行稳定性分析，当V_∞＝2330m/s时，由于μ＜0，系统不存在颤振失效风险；当V_∞＝2354m/s或V_∞＝2360m/s时，存在μ＞0的情况，即系统具有颤振失效风险，且V_∞＝2360m/s时颤振失效风险将大于V_∞＝2354m/s时的颤振失效风险。

综上所述，本发明提出了一种基于概率密度演化的随机气动弹性系统稳定性分析方法。首先建立气动弹性系统的有限元方程，并将其转化为广义特征值方程。在此基础上，建立气动弹性系统稳定性分析的概率密度演化方程，通过引入虚拟参数将概率密度演化方程转化为标准形式，采用有限差分方法和总变差减小格式，求解系统特征值最大实部的概率密度函数，根据最大实部的分布范围进行稳定性分析。数值结果表明，本发明方法得到的系统特征值最大实部的概率密度函数与蒙特卡洛方法吻合较好，并且能够大幅度减小计算时间，为随机气动弹性系统稳定性分析提供了新思路。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制，其可扩展应用于飞行器结构设计领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。

Claims (4)

1.一种基于概率密度演化的随机气动弹性系统稳定性分析方法，其特征在于实现步骤如下：

步骤(1)、建立气动弹性系统的有限元方程：

式中，M为质量矩阵，C为结构阻尼矩阵，为气动阻尼矩阵，K为结构刚度矩阵，为气动刚度矩阵，x(t)为广义坐标，为广义速度，为广义加速度，t为时间，下角标ΔQ_a表示气动；

步骤(2)、令x(t)＝x₀e^λt，可以将气动弹性系统的有限元方程转化为广义特征值方程：

Au＝λBu (2)

式中，λ为广义特征值；

步骤(3)、在确定性条件下，气动弹性系统特征值最大实部μ可以通过下式得到：

μ＝max(Re[λ(A,B)]) (3)

式中，Re表示特征值实部，λ(A,B)通过步骤(2)中的广义特征值方程得到；

步骤(4)、建立气动弹性系统稳定性分析的概率密度演化方程，表示为如下形式：

式中，α＝(α₁,...,α_s)为s维随机不确定参数，p_μα(μ,α,t)为(μ,α)的联合概率密度函数，

步骤(5)、在不确定参数α的变化域Ω内，均匀地取N_total个样本点，记为α_q(q＝1,...,N_to_tal)，并且将变化域Ω分为N_total个子域，记为Ω_q(q＝1,...,N_total)；

步骤(6)、将方程(4)在子域Ω_q内积分，可以得到：

步骤(7)、通过交换积分和求导次序，可以将方程(5)化简为：

式中，为对应于第q个样本点的概率密度函数；

步骤(8)、引入虚拟参数τ，令代入方程(6)中可以得到：

步骤(9)、确定初始条件为：

式中，δ为狄拉克函数，

步骤(10)、采用有限差分方法和总变差减小格式可以得到如下差分格式：

式中，τ_k＝kΔτ(k＝0,1,…)，为差分网格比，为限流器，和可表示为：

步骤(11)、将在N_total个样本点处计算出求和，可以得到：

步骤(12)、取τ_k＝1，则可得到特征值最大实部μ的概率密度函数的表达式：

步骤(13)、根据p_μ(μ)，若存在μ＞0，则该气动弹性系统具有颤振失效风险。

2.根据权利要求1所述的一种基于概率密度演化的随机气动弹性系统稳定性分析方法，其特征在于：所述步骤(7)中，交换积分和求导次序的具体方式为：

3.根据权利要求1所述的一种基于概率密度演化的随机气动弹性系统稳定性分析方法，其特征在于：所述步骤(9)中，可将初始条件离散为：

式中， 为方向的网格尺寸。

4.根据权利要求1所述的一种基于概率密度演化的随机气动弹性系统稳定性分析方法，其特征在于：所述步骤(10)中，差分格式要满足的收敛条件为：