非定常气动力最小状态有理近似的非线性优化算法
【技术领域】
[0001] 本发明涉及气动弹性力学领域，具体是一种用于频域非定常气动力有理近似的非 线性优化方法。
【背景技术】
[0002] 为便于气动伺服弹性系统的多学科优化设计，需要将弹性飞行器的运动方程转换 为时域的一阶时不变状态空间方程，而结构在任意运动下的非定常气动力模型则是其中一 个关键组成部分。在亚音速或超音速情况下，以偶极子格网法为基础计算得到的广义非定 常气动力系数矩阵（GAF)是一系列给定减缩频率的函数，代表结构在简谐振荡时所受的广 义气动力。为了将广义气动力变换到时域空间，需要将其近似延拓为拉氏域的有理函数形 式，再经过整理，结合弹性飞行器的运动方程即可得到气动弹性系统的状态空间模型。
[0003] 现在工程中最常用的有理函数近似（RFA)方法均以最小二乘法为基础，包括 Roger法、修正矩阵（MMP)法和最小状态（MS)法等。从这些方法出发得到的状态空间方程 中均包含由气动滞后根产生的气动状态扩充项。Roger法产生的气动扩充项数量为结构模 态数与气动滞后根数的乘积，MMP法对应的气动力扩充项数量为广义气动力系数矩阵各列 对应的气动滞后根数量之和，而MS法对应的气动力扩充项数等于气动滞后根的数量。大 量的工程实践表明，在气动力扩充项的数量相同时，MS法的拟合精度最高，但其中的非线性 最小二乘迭代过程计算量很大。文献1 "陈青.一种建立非定常气动力频域模型的简单方 法[J].空气动力学学报，1988年，第6卷第4期"吸取了 Roger近似式中同一矩阵各元素 独立确定时精度较高的特点，对选定模态受到的非定常气动力进行了精确拟合，并保留了 MS法的形式，改善了拟合精度，但对拟合项无合理加权且对气动滞后根的一维优化方法效 果不佳。针对MS法计算量大、初值的选取依赖于工程经验等缺点，文献2"何程.一种改 进的非定常气动力拟合方法[J].航空学报，1993年，第14卷第7期"通过调整拟合式最后 一项的形式并将其中气动滞后根取为大的互异负实数的方式将问题转化为线性拟合问题， 但要求滞后根的数量必须与结构模态数一致且拟合精度不高。文献3 "E. Nissim. On the Formulation of Minimum-State Approximation as a Nonlinear Optimization Problem I. Journal of Aircraft, Vol. 43, No. 4, 2006, pp. 1007 - 1013. " 在保留 MS 法拟合形式的基 础上释放了等式约束，并通过解析法得到了自变量的梯度算法，再结合对广义气动力系数 矩阵的比例缩放使拟合效率和精度都得到了提高，但非线性优化项的数量远大于气动滞后 根的个数，应用难度较大。目前非线性优化算法已经被应用到各类拟合方法中来优化气动 滞后根以减小拟合误差，这就进一步对拟合方法和非线性优化算法的效率提出了要求。
【发明内容】
[0004] 为克服已有方法精度有限或效率不高的缺点，本发明提出了一种非定常气动力最 小状态有理近似的非线性优化算法。
[0005] 本发明的具体过程是：
[0006] 步骤1，建立机翼有限元模型和气动力模型，计算机翼在给定减缩频率下的广义气 动力系数矩阵，具体是：
[0007] I建立机翼有限元模型，通过Nastran软件对建立的机翼有限元模型进行模态分 析，模态分析中取机翼的前9阶弹性模态。所述建立的机翼有限元模型中，机翼下翼面的弦 向剖面共有8个网格节点，分别位于距前缘0. 0%、5. 0%、10. 0%、27. 0%、45. 5%、63. 5%、 82.0%和100%处。沿机翼展向各网格节点数量与翼肋的数量相同，并与各翼肋的展向位置 一致。
[0008] II建立气动力模型，并将得到的模态分析结果导入到Zaero软件中，得到机翼在 给定减缩频率下的广义气动力系数矩阵：
[0009] Q(ik) = F(ik)+iG(ik) (I)
[0010] 其中Q(ik)为给定减缩频率下的广义气动力系数矩阵，/ = ，k = cob/v为减 缩频率，ω为机翼振荡的圆频率，b为机翼参考弦长，V为来流速度。F(ik)和G(ik)分别 表示广义气动力系数矩阵的实部和虚部。
[0011] 对Q (ik)拟合时采用MS法的拟合公式，即
[0013] 其中y =沾/F，S为拉氏变量。表示拟合得到的广义气动力系数矩阵， F即(ik)和Gap(ik)分别表示 的实部和虚部。A0e RnxnJ1G RnXn、A2e RnXn、D e Rnx' R e ΙΓΧ1\ E e Rmxn均为待定系数矩阵，n为结构模态数，m为气动滞后根数量。I为 m阶单 位矩阵。R为由气动滞后根组成的对角矩阵，表示为
[0014] R = diag(x) = diagl^x! X2 …xj) (3)
[0015] 其中X表示气动滞后根组成的向量。
[0016] 所述建立的气动力模型的弦向被均分为15等份，展向被均分为24等份。
[0017] 步骤2,给定气动滞后根向量的初值X。，对MS法的拟合公式进行总体拟合并计算 该总体拟合的误差f (xD)。
[0018] 在对得到的广义气动力系数矩阵进行拟合前，需先给定气动滞后根向量的初值， 初始化方法为：当 m = 1 时，X。= [-0· 3];当 m = 2 时，X。= [-0· 3-0. 5];当 m = 3 时，X。= [-0· 3-0. 5-0. 7];当 m = 4 时，X。= [-0· 3-0. 5-0. 7-0. 9]。
[0019] 根据给定的气动滞后根向量的初值X。得到由气动滞后根组成的对角矩阵R JtMS 法的拟合公式进行总体拟合并计算总体拟合误差f 〇〇，具体过程是：
[0020] I取结构的第r阶模态为关键模态，并令D矩阵的第r行元素为任意非零常数，其 中D矩阵的第r行元素代表结构的第r阶模态对应的广义气动力。此时方程（2)变为
[0022] 其中下标r均表示矩阵的第r行。
[0023] II约束方程在s = 0处的有理逼近值等于气动力系数矩阵列表值，在s = ikf处的 实部有理逼近值等于矩阵的列表值，在s = ikg处的虚部有理逼近值等于矩阵的列表值，其 中M5P k g均为指定的减缩频率值。之后从D矩阵第r行出发拟合E矩阵各列元素，所述E 矩阵第j列的加权最小二乘求解式如下
[0028] 其中L为减缩频率个数。Uik1)表示Q(ik)第r行第j列元素在减缩频率4处 的值，W rj表示Q(ik)的加权矩阵W的第r行第j列元素。巧(^)和$_)分别表示FQk1) 和GQk 1)的第j列。
[0029] 因当前机翼在颤振点处的减缩频率为0.07,故取kf= kg= 0.05以使颤振点附近 的拟合精度更高。
[0030] III求解出E矩阵后，再逐行用加权最小二乘法求解D矩阵的各元素，第r行除外。 D矩阵第i行的加权最小二乘求解式如下
[0034] W12表示加权矩阵W第i行各元素的平方值构成的对角矩阵，$(木)和&(味）分别 表示F Qk1)和G Qk1)的第i行；
[0035] IV计算对MS法的拟合公式进行总体拟合结果的总误差f (X。）
[0037] 其中i和j分别表示各矩阵的行和列。
[0038] 步骤3,开始对滞后根向量进行非线性优化。
[0039] 给定初始Hessian阵的逆阵H。，为m阶单位阵；令迭代下标kk = 0 ;
[0040] 步骤4,计算第一个迭代点X。处的梯度值V/Q <其第i项的算式如下
[0042] 其中O1为第i个元素为1的m阶单位向量，α为一个充分小的实数，此处取为 0. 001。f(xQ+a · ei)和f(xQ_a · ei)的计算方法同步骤2中f(x。）的计算方法。
[0043] 步骤5,确定搜索方向。
[0044] 通过公式（12)确定搜索方向
[0046] 步骤6,沿dkk线性搜索步长因子a kk。具体过程是：
[0047] I给定参数0〈 ξ〈0. 5和0〈 β〈1。取线性搜索时计算目标函数值的最大许可次数 为N1。令mm = 0代表本次循环目标函数值的计算次数，令nn = 0代表本次循环自标量的 越界次数；
[0048] II 令 xkk,st= X kk+P mm+nndkk，xkk,st代表从 X kk 出发得到的试探点；
[0049] m判断^ 或<·〒Τ = 1，…，相）是否满足，其中4 的第i项，I1和f分别表示自变量第i项的下界和上界。若都满足则转至步骤IV ;否则令 nn+1 - nn并转至步骤II，其中"一"表示赋值运算；
[0050] IV计算试探点xkk, st&的目标函数值f (X kk, st);
[0051] V判断mm是否超过&。若不超过，则转到步骤VI ;否则取mm为使f(xkk，st)最小的 数，令步长因子akk= β mm+nn，fkk+1= f(xkk,st)，停止线性搜索；
[0052] VI判断试探点处Armi jo不等式的满足情况
[0055] 若满足，则令步长因子Cikk= β mm+nn，fkk+1= f(xkk，st)，停止线性搜索；否则令 mm+1 - mm，转至步骤V继续搜索计算，直至mm超过N1或者在试探点处满足Armi jo不等式 (13)和（14)。
[0056] 步骤7,求下一个迭代点xkk+