基于CFD/CSD技术的跨/超声速飞行器高效高精度颤振时域分析方法

基于cfd/csd技术的跨/超声速飞行器高效高精度颤振时域分析方法
技术领域
1.本发明涉及跨/超声速速域飞行器的颤振分析领域，具体为一种基于cfd/csd技术的跨/超声速飞行器高效高精度颤振时域分析方法。


背景技术：

2.颤振是指由于弹性结构和气流的相互作用，最终导致结构灾难性破坏的不稳定自激振动现象。随着现代飞机飞行包线的不断拓宽，机翼、舵面等重要部件的颤振问题变得日益复杂和突出，直接关系到飞行器的安全。飞机在高空中一旦发生颤振，可能会在极短时间内解体，飞行员处置逃脱的概率也极小。因此，针对颤振问题分析研究已经成为了现代航空航天科学研究与工程应用中的重要一环。然而针对颤振问题的试验研究开展难度较大且成本高昂，飞行试验也存在着相当大的危险性，因此通过数值分析方法来进行颤振问题研究成为了重要的发展方向。
3.在跨声速及超声速飞行速域内，激波-激波干扰、激波-边界层干扰、分离涡脱落等现象将导致气动力具有明显的非线性效应，传统线性颤振分析方法无法获得高置信度结果，因而基于computational fluid dynamics(cfd)/computational structural dynamics(csd)技术的数值模拟方法成为了非线性颤振分析中的关键方法。但是在颤振时域分析时通过高精度cfd方法求解非定常气动力的过程计算量庞大，且耗时冗长，无法满足实际的工程应用需求，因此需要在高精度分析的同时提高分析效率。
4.目前对上述问题的主要处理方式包括如下两类方案：
5.其一是针对非线性非定常气动力的计算方法，处理方案包括降阶模型方法，频域谐波法以及时间谱方法等，但降阶模型方法在生成用于参数训练的样本快照时同样需要花费大量的计算时间，且缺乏对系统参数变化的鲁棒性；频域谐波法和时间谱方法需要对已有成熟cfd程序进行大量改动，且频域谐波法主要适用于周期性流动问题，时间谱方法随着时间步长的减小也存在稳定性和效率显著降低的问题。因此有必要使用高精度cfd方法进行跨/超声速速域的非线性非定常气动力计算，以保证计算精度与稳定性。
6.其次是针对颤振时域推进方法，处理方案包括改进欧拉法、线性多步法等，改进欧拉法精度较低，线性多步法虽然可以达到四阶精度，但是显式推进稳定性较差，若是采用预估-校正的线性多步法进行推进，虽然可以有效提高计算稳定性，但是在每个时间步需要计算两次cfd流场，时间花费同样十分巨大。因此有必要提供一种基于cfd/csd技术的高效高精度颤振时域推进方法以满足工程应用需求。


技术实现要素：

7.为进一步提高颤振时域分析的精度与效率，本发明提供一种基于cfd/csd技术的跨/超声速飞行器高效高精度颤振时域分析方法，以期为跨/超声速速域飞行器机翼、舵面等部件的颤振分析提供的技术支撑。
8.本发明的主要内容为：
9.针对跨/超声速速域由于流场的非线性效应导致活塞理论、牛顿冲击流理论、偶极子网格法等快速气动力计算方法精度较低的问题，选用高精度cfd方法——unsteady-reynolds-averaged-navier-stokes(urans)方法计算非定常气动力；针对改进欧拉法的时域推进方法精度较低以及显式线性多步法稳定性较低的不足，采用四阶milne-simpson预估-校正的线性多步法提高计算精度；针对传统预估-校正法计算在每个时间步需要进行两次cfd计算而导致效率低的弊端，假设广义气动力是时间与广义位移的二元函数，通过采取二元函数的三次lagrange插值多项式对广义气动力进行插值以替代原预估步后所需进行的cfd计算，从而在充分利用预估步计算信息的同时大幅提高计算效率，最终完成了跨/超声速速域飞行器机翼、舵面等部件非线性颤振的高效高精度时域分析。
10.本发明的技术方案为：
11.所述一种基于cfd/csd技术的跨/超声速飞行器高效高精度颤振时域分析方法，包括以下步骤：
12.步骤1：进行流场与结构数据初始化：
13.步骤1.1：建立飞行器的结构有限元模型，并划分有限元网格，施加约束以及载荷边界，且设置材料特性与单元属性；
14.步骤1.2：计算飞行器结构有限元模型的结构模态，包括模态振型与频率信息；
15.步骤1.3：建立飞行器的气动模型，划分网格，设置计算状态；
16.步骤1.4：根据计算状态，采用rans方法计算气动模型的初始定常流场，以作为非定常流场计算的初场，控制方程为n-s方程；
17.步骤2：建立用于时域推进的气动弹性控制方程：
18.步骤2.1：忽略步骤1.2中得到的模态振型中的对结构分析影响很小的高阶模态，仅保留前m阶模态，将对应的结构模态振型记作φ，采用径向基函数插值方法将φ从结构有限元网格点插值到气动物面网格点，记作φa，其中φa＝hφ，h表示将结构有限元网格点插值到气动物面网格点的插值系数矩阵；
19.步骤2.2：建立模态空间结构动力学方程：
[0020][0021]
其中，分别表示结构广义质量矩阵、广义阻尼矩阵和广义刚度矩阵，分别表示结构广义位移向量、广义速度向量和广义加速度向量，表示结构所受的广义力向量；
[0022]
步骤2.3：耦合非定常气动计算与结构动力学计算，将气动物面所受气动力作为结构所受外力，并将模态空间结构动力学方程作为气动弹性控制方程；
[0023]
步骤2.4：引入状态变量x，建立便于时域推进的状态空间形式气动弹性控制方程；
[0024][0025]
其中，a、b为控制方程的系数，x、a、b具体形式为：
[0026][0027]
其中，o表示零矩阵，i表示单位矩阵；
[0028]
步骤3；时域推进求解气动弹性控制方程，进行颤振时域分析：
[0029]
步骤3.1：设定时域推进求解的时间步长δt；
[0030]
步骤3.2：初始化广义位移、广义速度和广义力，获得x0、
[0031]
步骤3.3：设置自由来流速度u
∞
与自由来流动压q
∞
，用于广义气动力的计算；
[0032]
步骤3.4：采用显式方法求解n＝1，2，3时间步的状态变量x1、x2和x3；
[0033]
步骤3.5：当n≥4时，依据计算所得n时刻的广义位移采用计算n时刻的气动物面气动网格点的位移变形向量
[0034]
步骤3.6：通过基于rbf网格变形方法进行空间气动网格变形，获得n时刻变形后的气动网格；
[0035]
步骤3.7：采用urans方法计算n时刻的非定常流场；
[0036]
步骤3.8：计算n时刻的广义气动力
[0037]
步骤3.9：预估步计算，采用显式milnes四步法
[0038][0039]
求解n 1时刻状态变量的预估值
[0040]
步骤3.10：计算状态变量为时对应的广义气动力预估值，假设是关于和t的二元函数，采用二元函数三次lagrange插值多项式
[0041][0042]
插值获得广义气动力的预估值从而避免了采用uans方法计算广义气动力预估值；
[0043]
步骤3.11：校正步计算，采用隐式simpson两步法
[0044][0045]
求解n 1时刻状态变量的校正值x
n 1
；
[0046]
步骤4：判断颤振分析是否满足结束条件，若不满足结束条件，则更新迭代时间步，令n＝n 1，而后返回步骤3.5计算，直至满足结束条件后终止颤振计算；
[0047]
步骤5：结果后处理，具体步骤如下：
[0048]
通过多项式函数拟合自变量来流速度或来流动压与因变量平均对数衰减率或发散率所构成的离散坐标点，所得拟合函数的零点即为颤振速度或颤振动压；通过对计算所
得具有发散特性的时域响应曲线进行功率谱密度分析，获得颤振频率。
[0049]
进一步的，步骤1.4中，控制方程在曲线坐标系下的无量纲形式为
[0050][0051]
其中，t表示时间，ξ，η，ζ表示计算域中的坐标变量，q表示解矢量，fi，gi，hi分别表示ξ，η，ζ三个坐标方向的无粘通量矢量，表示ξ，η，ζ三个方向的粘性通量矢量。
[0052]
进一步的，步骤1.4中，控制方程中无粘项采用roe格式进行离散，粘性项采用二阶中心差分进行离散，湍流模型采用k-ωsst模型，时间推进采取隐式近似因子分解法提高定常流场计算的稳定性，并且为了进一步提高计算效率，采用多重网格加速收敛和并行计算方法。
[0053]
进一步的，步骤2.2中，以及和的具体形式为：
[0054][0055]
其中，m、d、k分别表示结构质量矩阵、阻尼矩阵和刚度矩阵，q表示结构运动的位移向量，f表示结构所受外力向量。
[0056]
进一步的，步骤2.3中，的具体形式为
[0057][0058]
其中，fa表示气动物面所受气动力向量，其具体形式为：
[0059][0060]
其中，fi表示气动物面第i个气动网格点的气动力，na表示气动物面气动网格点的总数，q
∞
表示自由来流的动压，si表示气动物面第i个气动网格点控制单元的面积，表示气动物面第i个气动网格点的当地静压系数，具体形式为：
[0061][0062]
其中，pi表示气动物面第i个气动网格点的当地无量纲静压，p
∞
表示自由来流的无量纲静压，γ表示空气比热比，ma表示自由来流的马赫数。
[0063]
进一步的，步骤2.3中，结合步骤2.1所求模态振型φa，的具体形式改写为更便于计算的形式：
[0064][0065]
进一步的，步骤3.4中，计算状态变量x1、x2和x3的过程为：
[0066]
利用显式euler公式求解n＝1时刻的状态变量x1；
[0067]
依据n＝1时刻的广义位移结合模态振型φa计算n＝1时刻的气动物面气动网格点的位移变形向量
[0068]
[0069]
并通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝1时刻变形后的气动网格；
[0070]
采用urans方法计算n＝1时刻的非定常流场，控制方程为n-s方程；方程中无粘项采用roe格式进行离散，粘性项采用二阶中心差分进行离散，湍流模型采用k-ωsst模型，时间推进采取双时间步法，采用步骤1.4所得定常流场作为初始计算流场，为了进一步提高计算效率，采用多重网格加速收敛和并行计算方法；
[0071]
根据获得n＝1时刻的气动物面气动力向量根据获得n＝1时刻的广义气动力
[0072]
更新迭代时间步，令n＝n 1，采用
[0073][0074]
显式adams两步法求解n＝2时刻的状态变量x2；
[0075]
依据n＝2时刻的广义位移采用计算n＝2时刻的气动物面气动网格点的位移变形向量通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝2时刻变形后的气动网格；继续采用urans方法计算n＝2时刻的非定常流场；根据获得n＝2时刻的气动物面气动力向量根据获得n＝2时刻的广义气动力
[0076]
更新迭代时间步，令n＝n 1，采用显式adams两步法求解n＝3时刻的状态变量x3，进而完成前3时间步的状态变量x的计算。
[0077]
进一步的，步骤4中判断颤振分析是否满足结束条件具体过程如下：
[0078]
监控各阶模态所对应广义位移的时域响应曲线，根据
[0079][0080]
计算曲线的对数衰减率以定量描述曲线的收敛或发散程度；其中，δi表示第i个周期的对数衰减率，ami和am
i 1
表示第i个周期相邻波峰与波谷的幅值，通过读取若干个波峰和波谷的幅值，计算平均对数衰减率，并根据数值的正负判定计算状态的收敛或发散。
[0081]
有益效果
[0082]
本发明提出了一种适用于跨/超声速速域飞行器机翼、舵面等部件非线性颤振的高效高精度时域分析方法，通过采用高精度cfd方法——urans方法计算非定常非线性气动流场，解决了跨/超声速速域由于流场的强非线性效应导致活塞理论、牛顿冲击流理论、偶极子网格法等快速气动力计算方法精度不足的问题，提高了对非定常非线性气动流场的计算精度；通过采用改进预估-校正的四阶线性多步法进行气动弹性控制方程的时域推进，解决了改进欧拉方法精度不足以及高阶显式方法稳定性差的问题，提高了时域推进求解的稳定性；并且在预估-校正计算时，通过假设广义气动力是时间与广义位移的二元函数，采用二元函数三次lagrange插值多项式对广义气动力预估值进行插值以替代原预估步后所需进行的cfd计算，从而克服了传统预估-校正法在一次预估校正计算时需要进行两次cfd求
解导致效率较低的弊端，同时，二元函数假设使得预估步计算的位移信息得到了充分利用，进一步提高了改进预估-校正四阶线性多步法的计算效率与稳定性。
[0083]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0084]
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0085]
图1是本发明基于cfd/csd技术的高效高精度颤振时域分析方法的基本流程示意图。
[0086]
图2是实施例一中，agard445.6机翼标模的几何示意图。
[0087]
图3是实施例一中，agard445.6机翼标模的有限元模型。
[0088]
图4是实施例一中，agard445.6机翼标模的前四阶模态振型云图。
[0089]
图5是实施例一中，agard445.6机翼标模的气动网格。
[0090]
图6是实施例一中，针对马赫数为0.901的计算状态，agard445.6机翼标模的物面压力系数分布云图。
[0091]
图7是实施例一中，针对马赫数为0.901的计算状态，在来流速度为330m/s时，一阶模态的广义位移时域响应曲线。
[0092]
图8是实施例一中，本发明计算所得颤振动压和颤振频率与参考实验值对比结果。
[0093]
图9是实施例二中，超声速全动舵面模型的几何示意图。
具体实施方式
[0094]
本发明为进一步提高颤振时域分析的精度与效率，本发明提供一种基于cfd/csd技术的跨/超声速速域飞行器高效高精度颤振时域分析方法，包括以下步骤：
[0095]
步骤1：流场与结构数据初始化，具体步骤如下：
[0096]
步骤1.1：建立计算模型的结构有限元模型，划分有限元网格，施加约束以及载荷边界，设置材料特性与单元属性；
[0097]
步骤1.2：计算有限元模型的结构模态，包括模态振型与频率信息；
[0098]
步骤1.3：建立计算模型的气动模型，划分多块结构网格，设置计算状态；
[0099]
步骤1.4：根据计算状态，采用reynolds-averaged navier
–
stoke(rans)方法计算气动模型的初始定常流场，以作为非定常流场计算的初场，控制方程为navier
–
stokes(n-s)方程，其在曲线坐标系下的无量纲形式如式(1)所示：
[0100][0101]
其中，t表示时间，ξ，η，ζ表示计算域中的坐标变量，q表示解矢量，fi，gi，hi分别表示ξ，η，ζ三个坐标方向的无粘通量矢量，表示ξ，η，ζ三个方向的粘性通量矢量。
[0102]
方程中无粘项采用roe格式进行离散，粘性项采用二阶中心差分进行离散，湍流模型采用k-ωsst模型，时间推进采取隐式近似因子分解法以提高定常流场计算的稳定性，为了进一步提高计算效率，采用多重网格加速收敛技术和并行计算技术。
[0103]
步骤2：建立用于时域推进的气动弹性控制方程，具体步骤如下：
[0104]
步骤2.1：忽略对结构分析影响很小的步骤1.2所得高阶模态，仅保留前m阶模态，将对应的结构模态振型记作φ，采用径向基函数(rbf)插值方法将φ从结构有限元网格点插值到气动物面网格点，记作φa，且二者满足如式(2)所示的转换关系。
[0105]
φa＝hφ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0106]
其中，h表示将结构有限元网格点变量插值到气动物面网格点的插值系数矩阵；
[0107]
步骤2.2：建立如式(3)所示的模态空间结构动力学方程：
[0108][0109]
其中，分别表示结构广义质量矩阵、广义阻尼矩阵和广义刚度矩阵，分别表示结构广义位移向量、广义速度向量和广义加速度向量，表示结构所受广义力向量，上述各变量的具体形式如式(4)所示；
[0110][0111]
其中，m、d、k分别表示结构质量矩阵、阻尼矩阵和刚度矩阵，q表示结构运动的位移向量，f表示结构所受外力向量；
[0112]
步骤2.3：耦合非定常气动计算与结构动力学计算，将气动物面所受气动力作为结构所受外力，则式(3)可作为气动弹性控制方程，的具体形式如式(5)所示：
[0113][0114]
其中，fa表示气动物面气动力向量，其具体形式如式(6)所示：
[0115][0116]
其中，fi表示第i个气动物面气动网格点的气动力，na表示气动物面气动网格点的总数，q
∞
表示自由来流的动压，si表示第i个气动物面气动网格点控制单元的面积，表示第i个气动物面气动网格点的当地静压系数，其具体形式如式(7)所示：
[0117][0118]
其中，pi表示i个气动物面气动网格点的当地无量纲静压，p
∞
表示自由来流的无量纲静压，γ表示空气比热比，ma表示自由来流的马赫数；
[0119]
进一步的，结合步骤2.1所求模态振型φa，式(5)可改写为如式(8)所示更便于计算的形式：
[0120][0121]
步骤2.4：引入状态变量x，建立如式(9)所示便于时域推进的状态空间形式气动弹性控制方程；
[0122][0123]
其中，a、b为控制方程的系数，x、a、b具体形式如式(10)所示；
[0124][0125]
其中，o表示零矩阵，i表示单位矩阵；
[0126]
步骤3；时域推进求解气动弹性控制方程，进行颤振时域分析，具体步骤如下：
[0127]
步骤3.1：设定时域推进求解的时间步长δt；
[0128]
步骤3.2：初始化广义位移、广义速度和广义力，获得x0、
[0129]
步骤3.3：设置自由来流速度u
∞
与自由来流动压q
∞
，用于式(6)广义气动力的计算；
[0130]
步骤3.4：采用显式方法求解n＝1，2，3时间步的状态变量x，具体步骤如下：
[0131]
采用如式(11)所示的显式euler公式求解n＝1时刻的状态变量x1；
[0132][0133]
依据计算所得n＝1时刻的广义位移结合模态振型φa计算n＝1时刻的气动物面气动网格点的位移变形向量其具体方法如式(12)所示：
[0134][0135]
其中qa表示气动物面位移变形向量；
[0136]
通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝1时刻变形后的气动网格；
[0137]
采用urans方法计算n＝1时刻的非定常流场，控制方程仍为n-s方程，如式(1)所示；方程中无粘项采用roe格式进行离散，粘性项采用二阶中心差分进行离散，湍流模型采用k-ωsst模型，时间推进采取双时间步法，采用步骤1.4所得定常流场作为初始计算流场，为了进一步提高计算效率，采用多重网格加速收敛技术和并行计算技术；
[0138]
依据式(6)获得n＝1时刻的气动物面气动力向量
[0139]
依据式(8)获得n＝1时刻的广义气动力
[0140]
更新迭代时间步，令n＝n 1，采用如式(13)所示的显式adams两步法求解n＝2时刻的状态变量x2；
[0141][0142]
依据计算所得n＝2时刻的广义位移采用式(12)计算n＝2时刻的气动物面气动网格点的位移变形向量
[0143]
通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝2时刻变形后的气动网格；
[0144]
继续采用urans方法计算n＝2时刻的非定常流场；
[0145]
依据式(6)获得n＝2时刻的气动物面气动力
[0146]
依据式(8)获得n＝2时刻的广义气动力
[0147]
更新迭代时间步，令n＝n 1，采用如式(13)所示的显式adams两步法求解n＝3时刻的状态变量x3，进而完成前3时间步的状态变量x的计算；
[0148]
步骤3.5：当n≥4时，依据计算所得n时刻的广义位移采用式(12)计算n时刻的气动物面气动网格点的位移变形向量
[0149]
步骤3.6：通过基于rbf网格变形方法进行空间气动网格变形，获得n时刻变形后的气动网格；
[0150]
步骤3.7：继续采用urans方法计算n时刻的非定常流场；
[0151]
步骤3.8：依据式(6)获得n时刻的广义气动力
[0152]
步骤3.9：预估步计算，采用如式(14)所示的显式milnes四步法求解n 1时刻状态变量的预估值
[0153][0154]
步骤3.10：状态变量为时对应的广义气动力预估值计算，假设是关于和t的二元函数，采用如式(15)所示的二元函数三次lagrange插值多项式插值获得广义气动力的预估值从而避免了采用uans方法计算广义气动力预估值；
[0155][0156]
步骤3.11：校正步计算，采用如式(16)所示的隐式simpson两步法求解n 1时刻状态变量的校正值x
n 1
；
[0157][0158]
步骤4：判断颤振分析是否满足结束条件，具体步骤如下：
[0159]
步骤4.1：监控各阶模态所对应广义位移的时域响应曲线，根据式(17)计算曲
[0160]
线的对数衰减(发散)率以定量描述曲线的收敛或发散程度；
[0161][0162]
其中，δi表示第i个周期的对数衰减(发散)率，ami和am
i 1
表示第i个周期相邻波峰与波谷的幅值，通过读取若干个波峰和波谷的幅值，计算平均对数衰减(发散)率，并根据数值的正负判定计算状态的收敛与发散，进而停止颤振计算；
[0163]
步骤4.2：若不满足结束条件，则更新迭代时间步，令n＝n 1，而后返回步骤3.5计算，直至满足步骤4.1的结束条件后终止颤振计算；
[0164]
步骤5：结果后处理，具体步骤如下：
[0165]
通过多项式函数拟合自变量来流速度(来流动压)与因变量平均对数衰减(发散)率所构成的离散坐标点，所得拟合函数的零点即为颤振速度(颤振动压)；
[0166]
通过对计算所得具有发散特性的时域响应曲线进行功率谱密度分析，获得颤振频率。
[0167]
下面详细描述本发明的实施例，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0168]
实施例一：
[0169]
以agard445.6机翼标模的颤振分析为例，该模型于1961年在nasa兰利中心tdt风洞进行实验并获得了较为完备的实验数，如图2所示，其展弦比为1.65，梢根比为0.66，展长为0.762m，根弦长为0.5587m，机翼面积为0.355m2，四分之一弦线后掠角为45
°
，所采用的翼型为naca65a004。
[0170]
步骤1：流场与结构数据初始化；
[0171]
步骤1.1：建立如图3所示机翼模型的结构有限元模型，首先对其进行有限元网格划分，共包含2475个结构点，单元类型选取为solid，共包含536个solid单元，约束条件为根部固支，材料为桃花心木，密度为381.98kg/m3，弹性模量如式(18)所示。
[0172][0173]
步骤1.2：计算机翼有限元模型的结构模态，前四阶模态振型如图4所示，前四阶模态频率{f1,f2,f3,f4}如表1所示；
[0174]
表1结构模态频率
[0175]
模态一阶振型二阶振型三阶振型四阶振型模态频率(hz)f1＝9.6096f2＝37.5400f3＝47.6890f4＝92.2395
[0176]
步骤1.3：建立机翼模型的气动模型，划分如图5所示的多块结构网格，网格量为223万，计算状态如表2所示：
[0177]
表2实施例1计算状态
[0178][0179]
步骤1.4：采用cfd方法分别求解各计算状态下气动模型的初始定常流场，作为非定常气动计算的初场，如针对马赫数为0.901的计算状态，图6所示为agard445.6机翼定常计算所得的物面压力系数分布；
[0180]
步骤2：建立用于时域推进的气动弹性控制方程；
[0181]
步骤2.1：忽略对本结构影响很小的高阶模态，仅保留前4阶模态，将对应的模态振
型记作φ，采用径向基函数(rbf)插值方法将φ从结构有限元网格点插值到气动物面网格点，记作φa；
[0182]
步骤2.2：建立如式(3)所示的模态空间结构动力学方程，并对进行单位化，假设阻尼为rayleigh阻尼，则的具体形式如式(19)所示；
[0183][0184]
针对本实施例，假设前四阶模态的阻尼比均为0，则具体取值如式(20)所示；
[0185][0186]
步骤2.3：耦合非定常气动计算与结构动力学计算，将气动物面所受气动力作为结构所受外力，则式(3)可作为气动弹性控制方程；
[0187]
步骤2.4：引入状态变量x，根据的取值，可建立如式(9)所示便于时域推进的状态空间形式气动弹性控制方程，其中系数矩阵a、b的取值如式(21)所示；
[0188][0189]
步骤3；时域推进求解气动弹性控制方程，进行颤振时域分析；
[0190]
步骤3.1：设定时域推进求解的时间步长δt＝6
×
10-5
s；
[0191]
步骤3.2：初始化广义气动力、广义位移和广义速度，其取值如式(22)所示；
[0192][0193]
步骤3.3：设置自由来流速度u
∞
与自由来流动压q
∞
，用于式(6)广义气动力的计算，针对不同计算状态下，本实施例预进行每个计算状态下的颤振速度与颤振动压的计算，因此需要针对每一个计算状态，设置一系列来流速度与对应的来流动压并分别进行颤振分析，进而提高后处理时的数据拟合的置信度，如针对马赫数为0.901的计算状态，来流速度序列和来流动压序列的取值如式(23)所示；
[0194][0195]
步骤3.4：针对本实施例不同的计算状态，分别采用显式方法求解n＝1，2，3时刻的状态变量x，；
[0196]
其中n＝1时刻的状态变量x1，采用如式(11)所示的显式euler公式求解；依据计算
所得n＝1时刻的广义位移结合模态振型φa计算n＝1时刻的气动物面气动网格点的位移变形向量
[0197]
通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝1时刻变形后的气动网格；通过采用urans方法求解控制方程式(1)，计算n＝1时刻的非定常流场；依据式(6)获得n＝1时刻的气动物面气动力向量依据式(8)获得n＝1时刻的广义气动力
[0198]
更新迭代时间步，令n＝n 1，为提高时间精度，采用如式(13)所示的显式adams两步法求解n＝2时刻的状态变量x2；
[0199]
依据计算所得n＝2时刻的广义位移采用式(12)计算n＝2时刻的气动物面气动网格点的位移变形向量通过基于rbf网格变形方法进行流场空间气动网格变形，获得n＝2时刻变形后的气动网格；继续采用cfd方法计算n＝2时刻的非定常流场；依据式(6)获得n＝2时刻的气动物面气动力依据式(8)获得n＝2时刻的广义气动力
[0200]
更新迭代时间步，令n＝n 1，采用如式(13)所示的显式adams两步法求解n＝3时刻的状态变量x3，进而完成前3时间步的状态变量x的计算，如针对马赫数为0.901的计算状态，x的取值如式(24)所示；
[0201][0202]
步骤3.5：针对本实施例不同的计算状态，当n≥4时，依据计算所得当前n时刻的广义位移采用式(12)计算n时刻的气动物面气动网格点的位移变形向量
[0203]
步骤3.6：针对本实施例不同的计算状态，通过基于rbf网格变形方法进行空间气动网格变形，获得n时刻变形后的气动网格；
[0204]
步骤3.7：针对本实施例不同的计算状态，继续采用cfd方法计算n时刻的非定常流场；
[0205]
步骤3.8：针对本实施例不同的计算状态，依据式(6)获得n时刻的广义气动力
[0206]
步骤3.9：针对本实施例不同的计算状态，进行预估步计算，采用如式(14)所示的显式milnes四步法求解n 1时刻状态变量的预估值
[0207]
步骤3.10：由于非定常气动力同时受时间与物面变形的影响，因此可以假设广义气动力是关于和t的二元函数，即为了充分利用步骤3.8所得预估值针
对本实施例不同的计算状态，采用如式(15)所示的二元函数三次lagrange插值多项式插值获得广义气动力的预估值避免采用高精度cfd方法——urans方法求解非定常流场，同时保证了广义气动力预估值的预测精度；
[0208]
步骤3.11：针对本实施例不同的计算状态，进行校正步计算，采用如式(16)所示的隐式simpson两步法求解n 1时刻状态变量的校正值x
n 1
；
[0209]
步骤4：判断颤振分析是否满足结束条件；
[0210]
步骤4.1：针对本实施例不同的计算状态，监控各阶模态所对应广义位移的时域响应曲线，通过读取30～50个周期波峰和波谷的幅值来计算各个周期内的对数衰减(发散)率，进而计算平均对数衰减(发散)率，并根据数值的正负判定计算状态的收敛与发散，停止颤振计算，如针对马赫数为0.901的计算状态，其在来流速度为330m/s时，时域响应曲线发散，即来流速度大于颤振速度，计算所得一阶模态广义位移时域响应曲线如图7所示；
[0211]
步骤4.2：若不满足结束条件，则更新迭代时间步，令n＝n 1，而后返回步骤3.5计算，直至满足步骤4.1的结束条件后终止颤振计算；
[0212]
步骤5：结果后处理；
[0213]
针对本实施例不同的计算状态，通过多项式拟合自变量来流速度(来流动压)与因变量平均对数衰减(发散)率所构成的离散坐标点，根据拟合函数零点获得颤振速度(颤振动压)，如针对马赫数为0.901的计算状态，计算所得颤振速度为296m/s；
[0214]
通过对具有发散特性的时域响应曲线进行功率谱密度分析，获得颤振频率；
[0215]
如图8所示为本发明计算所得颤振动压和颤振频率与参考实验值的对比结果，计算颤振动压和颤振频率与实验值吻合良好。
[0216]
表3是本发明的计算结果与改进欧拉方法、四阶显式adams线性多步法、四阶adams预估-校正法的对比结果，可见本发明在计算精度方面相较于改进欧拉方法具有明显优势，在计算稳定性方面相较于显式adams线性多步法具有明显优势，在计算效率方面相较于显式较于四阶显式adams线性多步法以及四阶adams预估-校正法具有明显优势；
[0217]
表3针对马赫数为0.901的计算状态本发明计算方法与传统方法对比
[0218][0219]
实施例二：
[0220]
以超声速全动舵面模型的颤振分析为例，该模型在中国空气动力学研究所的fd-06风洞进行颤振实验并获得了较为完备的实验数，模型具体几何尺寸如图9所示；
[0221]
表4是本发明计算所得颤振速度与参考实验值的对比结果，计算颤振速度与实验值吻合良好。
[0222]
表4超声速全动舵面颤振分析结果对比
[0223][0224]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。