[发明专利]非定常气动力最小状态有理近似的非线性优化算法

摘要

一种非定常气动力最小状态有理近似的非线性优化算法。本发明对D矩阵的关键模态行初始化再顺序求解E矩阵的各列和D矩阵的其余行，并充分利用当前迭代点处的函数值信息，有效提高了逆Hessian阵的近似精度，提高了外层非线性优化算法的效率。本发明有效提高了广义气动力矩阵中关键模态行元素的拟合精度，有利于提高时域颤振分析结果精度和突风响应分析结果精度。将频域下离散的非定常气动力系数矩阵以MS法的形式延拓至拉氏域，且能有效提高计算效率和关键模态项的精度。

主权项

一种非定常气动力最小状态有理近似的非线性优化算法，其特征在于，具体过程是：步骤1，建立机翼有限元模型和气动力模型，计算机翼在给定减缩频率下的广义气动力系数矩阵，具体是：Ⅰ建立机翼有限元模型，通过Nastran软件对建立的机翼有限元模型进行模态分析，模态分析中取机翼的前9阶弹性模态；Ⅱ建立气动力模型，并将得到的模态分析结果导入到Zaero软件中，得到机翼在给定减缩频率下的广义气动力系数矩阵：Q(ik)＝F(ik)+iG(ik)                         (1)其中Q(ik)为给定减缩频率下的广义气动力系数矩阵，k＝ωb/V为减缩频率，ω为机翼振荡的圆频率，b为机翼参考弦长，V为来流速度；F(ik)和G(ik)分别表示广义气动力系数矩阵的实部和虚部；对Q(ik)拟合时采用MS法的拟合公式，即Qap(s‾)=Fap(ik)+iGap(ik)=A0+A1s‾+A2s‾2+D(s‾I-R)-1E·s‾---(2)]]>其中s为拉氏变量；表示拟合得到的广义气动力系数矩阵，Fap(ik)和Gap(ik)分别表示的实部和虚部；A0∈Rn×n、A1∈Rn×n、A2∈Rn×n、D∈Rn×m、R∈Rm×m、E∈Rm×n均为待定系数矩阵，n为结构模态数，m为气动滞后根数量；I为m阶单位矩阵；R为由气动滞后根组成的对角矩阵，表示为R＝diag(x)＝diag([x1 x2…xm])                  (3)其中x表示气动滞后根组成的向量；步骤2，给定气动滞后根向量的初值x0，对MS法的拟合公式进行总体拟合并计算该总体拟合的误差f(x0)；给定气动滞后根向量的初值x0；根据给定的气动滞后根向量的初值x0得到由气动滞后根组成的对角矩阵R；对MS法的拟合公式进行总体拟合并计算总体拟合误差f(x0)，具体过程是：Ⅰ取结构的第r阶模态为关键模态，并令D矩阵的第r行元素为任意非零常数，其中D矩阵的第r行元素代表结构的第r阶模态对应的广义气动力；此时方程(2)变为Qap,r(s‾)=A0,r+A1,rs‾+A2,rs‾2+Dr(s‾I-R)-1E·s‾---(4)]]>其中下标r均表示矩阵的第r行；Ⅱ约束方程在s＝0处的有理逼近值等于气动力系数矩阵列表值，在s＝ikf处的实部有理逼近值等于矩阵的列表值，在s＝ikg处的虚部有理逼近值等于矩阵的列表值，其中kf和kg均为指定的减缩频率值；之后从D矩阵第r行出发拟合E矩阵各列元素，所述E矩阵第j列的加权最小二乘求解式如下Σl=1L(PlTWrj2Pl+QlTWrj2Ql)Ej=Σl=1L{PlTWrj2F‾j(ikl)+QlTWrj2G‾j(ikl)}---(5)]]>其中Ej表示E矩阵第j列，Pl=kl2Dr[(kl2I+R2)-1-(kf2I+R2)-1]Ql=klDr[(kl2I+R2)-1-(kg2I+R2)-1]R---(6)]]>Wrj=1maxl{|Qrj(ikl)|,1}---(7)]]>其中L为减缩频率个数；Qrj(ikl)表示Q(ik)第r行第j列元素在减缩频率kl处的值，Wrj表示Q(ik)的加权矩阵W的第r行第j列元素；和分别表示F(ikl)和G(ikl)的第j列；因当前机翼在颤振点处的减缩频率为0.07，故取kf＝kg＝0.05以使颤振点附近的拟合精度更高；Ⅲ求解出E矩阵后，再逐行用加权最小二乘法求解D矩阵的各元素，第r行除外；D矩阵第i行的加权最小二乘求解式如下Σl=1L(Pl*TWi2Pl*+Ql*TWi2Ql*)DiT=Σl=1L(Pl*TWi2F^i(ikl)+Ql*TWi2G^i(ikl))---(8)]]>其中Di表示D矩阵第i行，Pl*=kl2ET[(kl2I+R2)-1-(kf2I+R2)-1]Ql*=klETR[(kl2I+R2)-1-(kg2I+R2)-1]---(9)]]>Wi2表示加权矩阵W第i行各元素的平方值构成的对角矩阵，和分别表示F(ikl)和G(ikl)的第i行；Ⅳ计算对MS法的拟合公式进行总体拟合结果的总误差f(x0)f(x0)=Σl=1LΣj=1nΣi=1nWij2{[Fap,ij(kl)-Fij(kl)]2+[Gap,ij(kl)-Gij(kl)]2}---(10)]]>其中i和j分别表示各矩阵的行和列；步骤3，开始对滞后根向量进行非线性优化；给定初始Hessian阵的逆阵H0，为m阶单位阵；令迭代下标kk＝0；步骤4，计算第一个迭代点x0处的梯度值其第i项的算式如下(▿f0)i=f(x0+α·ei)-f(x0-α·ei)2α,(i=1,2,...m)---(11)]]>其中ei为第i个元素为1的m阶单位向量，α为一个充分小的实数，此处取为0.001；f(x0+α·ei)和f(x0‑α·ei)的计算方法同步骤2中f(x0)的计算方法；步骤5，确定搜索方向；通过公式(12)确定搜索方向dkk=-Hkk▿fkk---(12)]]>步骤6，沿dkk线性搜索步长因子αkk；具体过程是：Ⅰ给定参数0<ξ<0.5和0<β<1；取线性搜索时计算目标函数值的最大许可次数为N1；令mm＝0代表本次循环目标函数值的计算次数，令nn＝0代表本次循环自标量的越界次数；Ⅱ令xkk,st＝xkk+βmm+nndkk，xkk,st代表从xkk出发得到的试探点；Ⅲ判断或是否满足，其中表示xkk,st的第i项，xi和分别表示自变量第i项的下界和上界；若都满足则进行本步骤的第Ⅳ步；否则令nn+1→nn并转至本步骤第Ⅱ步,其中“→”表示赋值运算；Ⅳ计算试探点xkk,st处的目标函数值f(xkk,st)；Ⅴ判断mm是否超过N1；若不超过，则进行本步骤的第Ⅵ步；否则取mm为使f(xkk,st)最小的数，令步长因子αkk＝βmm+nn，fkk+1＝f(xkk,st)，停止线性搜索；Ⅵ判断试探点处Armijo不等式的满足情况f(xkk,st)-f(xkk)≤ξβmm+nn▿fkkTdkk---(13)]]>f(xkk,st)-f(xkk)≥(1-ξ)βmm+nn▿fkkTdkk---(14)]]>若满足，则令步长因子αkk＝βmm+nn，fkk+1＝f(xkk,st)，停止线性搜索；否则令mm+1→mm，转至本步骤的第Ⅴ步继续搜索计算，直至mm超过N1或者在试探点处满足Armijo不等式(13)和(14)；步骤7，求下一个迭代点xkk+1：xkk+1＝xkk+αkkdkk                         (15)步骤8，计算点xkk+1处的梯度值计算方法同步骤4；步骤9，判断非线性优化过程是否收敛，具体是：Ⅰ判断||xkk||>ε6是否成立，若不成立则进行本步骤的第Ⅱ步，否则再判断||xkk+1‑xkk||/||xkk||≤ε1是否成立，若成立则停止优化迭代，否则转本步骤的第Ⅲ步；Ⅱ判断||xkk+1‑xkk||≤ε2是否成立，若成立则停止优化迭代，否则进行本步骤的第Ⅲ步；Ⅲ判断|fkk|>ε7是否成立，若不成立则转本步骤的第Ⅳ步，否则再判断|fkk‑fkk+1|/|fkk|≤ε3是否成立，若成立则停止优化迭代，否则转本步骤的第Ⅴ步；Ⅳ判断|fkk‑fkk+1|≤ε4是否成立，若成立则停止优化迭代，否则进行本步骤的第Ⅴ步；Ⅴ判断是否成立，若成立则停止优化迭代，否则转步骤10，继续进行优化迭代；在上述的收敛准则中，ε1、ε2、ε3、ε4、ε5为误差限，取为10‑6；ε6、ε7分别用于判断||xk||和|fk|的数量级，取为10‑4；步骤10，计算Hkk+1；具体是：Ⅰ令自变量变化为skk＝xkk+1‑xkk，梯度差为Ⅱ通过改进的BFGS校正公式(16)和(17)计算逆Hessian阵Hkk+1：Hkk+1=(I-skky‾kkTy‾kkTskk)Hkk(I-y‾kkskkTy‾kkTskk)+skkskkTy‾kkTskk,y‾kkTskk>0Hkk,y‾kkTskk≤0---(16)]]>其中y‾kk=ykk+12(fkk-fkk+1)+5▿fkk+1Tskk+(7-αkk)▿fkkTskkskkTθkkθkk---(17)]]>式中，表示修正的梯度差，θkk∈Rm是满足的任意向量，此处取为skk；步骤11，令迭代下标kk→kk+1，转到步骤5继续迭代，直至满足步骤9中的收敛条件；至此，完成了对NACA0012对称翼型的M6机翼的非定常气动力最小状态有理近似的非线性优化算法。