CN102880804B - 一种基于振型叠加法的飞行器翼肋动力学确定方法 - Google Patents

Abstract

一种基于振型叠加法的飞行器翼肋动力学确定方法，包括步骤:导入飞行器翼肋初始几何模型与动力载荷；将翼肋的几何模型化为有限元模型；初始化有限元模型需优化部分的各个单元的伪密度值；将赋予伪密度值的翼肋有限元模型整体代人求解器中进行模态分析；基于模态分析结果采用振型叠加法进行翼肋动力响应的近似求解；基于动力分析结果求解翼肋的变形能，采用的纽曼展开差分方法求解变形能和约束对伪密度值的灵敏度；将这些值传递给优化器，得到更新值后的伪密度值；如果满足收敛准则，停止迭代否则返回到有限元求解部分；生成翼肋最终的动力学设计方案。本发明在考虑动力学的工况条件下，降低了飞行器质量，提高了飞行器性能。

Description

一种基于振型叠加法的飞行器翼肋动力学确定方法

技术领域

本发明涉及一种基于振型叠加法的飞行器翼肋动力学确定方法，属于飞行器翼肋结构设计领域。

背景技术

飞行器结构不可避免地受到气动、风载等来自于外界环境的其它振动或冲击干扰作用。飞机、火箭、导弹等飞行器往往要求结构的圆频率远离发动机振动频率以避免结构共振现象的发生。通常要求载体颠簸产生的振动对仪器、仪表等装置造成的危害降至最低限度，而且还要增加乘员乘坐的舒适感。有效地降低振动强度，可以驱除乘员的疲劳与不适，确保装置正常工作。显然欲降低飞行器重量的同时提高机械系统的动态特性，达到控制振动目的或者确保它们在动力环境下能够安全可靠的工作，最为行之有效的办法就是进行结构的动力拓扑优化设计。

在飞行器结构拓扑优化设计研究中，大部分研究都集中在静载荷的研究中，在动力结构拓扑优化上所做的工作非常有限。而飞行器中存在大量翼肋需要进行考虑动力学的拓扑优化设计。由于振动和噪音的主要原因是机械装置承受周期的载荷作用，以往讨论的对周期载荷下飞行器结构进行拓扑优化设计中，构造动柔度这一概念来衡量物体振动的强弱。但是以往研究中对动柔度的定义不统一，并且现有的动柔度的定义无法与静柔度的定义相容。由于相容的动柔度求解与灵敏度分析存在很多困难，同时计算过程中存在局部模态问题，现有方法可以实现对一些与静柔度不相容的目标进行拓扑优化，并且绕过局部模态的处理，最终得到的设计方案没有实际的工程意义。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种基于振型叠加法的飞行器翼肋动力学确定方法，采用相容的变形能作为动柔度，变形能为设计目标对飞行器翼肋进行拓扑优化设计，采用振型叠加法克服了灵敏度求解困难与局部模态处理问题。

本发明技术解决方案：一种基于振型叠加法的飞行器翼肋动力学确定方法，其特征在于实现步骤如下：

(1)导入飞行器翼肋初始几何模型或有限元模型，以及飞行条件下的周期激励动力载荷；

(2)将飞行器翼肋的初始几何模型离散为由有限元单元组成的有限元模型，如果已有有限元模型则直接进行下一步；

(3)确定飞行器翼肋结构中需要拓扑优化的部分，赋予所述步骤(1)中有限元模型中的有限元单元不同属性以区分需进行拓扑优化部分及不需要进行优化部分，再初始化飞行器翼肋有限元模型中需要进行优化部分的各个单元的伪密度值ρ_n，n为需要进行优化部分有限元单元编号，所述伪密度值ρ_n为0到1间连续的数值，0代表此处单元为空，1代表此处有单元，其它数值代表中间状态；

(4)基于伪密度值对有限元模型中需进行优化部分的有限元单元的刚度矩阵进行插值：k_n＝ρ_n ^mk_n ⁰m>1，其中n为需要进行优化部分有限元单元编号，m为惩罚因子，k_n ⁰为实体材料单元刚度矩阵,不需要优化部分的有限元单元为实体，单元刚度矩阵保持不变；然后将基于伪密度值插值的飞行器翼肋有限元模型整体代人有限元求解器中进行模态分析，同时避免数值模态的发生有限元分析时不考虑的低伪密度值单元，只考虑高伪密度值单元，得到飞行器翼肋阵型与圆频率，当伪密度值大于0.1定义为高伪密度值，伪密度值不大于0.1定义为低伪密度值；

(5)基于步骤(4)模态分析得到的飞行器翼肋阵型与圆频率与飞行条件下的周期激励动力载荷，采用振型叠加法进行飞行器翼肋动力响应的近似求解；

(6)利用步骤(5)得到的飞行器翼肋的动力学的响应，来求解飞行器翼肋在动力载荷下的设计目标，即飞行器翼肋的变形能、约束函数的值；

(7)进行飞行器翼肋优化过程中需要灵敏度信息，采用纽曼展开差分方法求解步骤(6)中变形能对伪密度值的灵敏度，以及约束函数的灵敏度，同时低伪密度值单元的灵敏度由周围高伪密度值单元插值得到；当伪密度值大于0.1定义为高伪密度值，伪密度值不大于0.1定义为低伪密度值；

(8)将步骤(7)计算得到的飞行器翼肋的变形能，约束函数的值以及灵敏度传递给优化器，以变形能最小为优化目标得到更新的伪密度值；

(9)判断更新后得到的有限单元的当前伪密度值与前一次更新前的有限单元的伪密度值之间的差异，如果所有有限单元的伪密度值都小于设定的小量，则停止计算转到步骤(10)，否则返回到步骤(4)继续进行飞行器翼肋的动力学更新设计；

(10)设定阈值，伪密度值大于此阈值为实体，伪密度值小于阈值为空洞，将得到的拓扑结构的有限元模型导回为几何模型，确定飞行器翼肋最终的动力学方案。

所述步骤(5)用振型叠加法进行翼肋动力响应的近似求解公式：

当结构受到周期激励时，M为飞行器翼肋结构质量阵，K为飞行器翼肋结构刚度阵，F为激励载荷向量，X为系统随着时间t的响应，阻尼C为瑞利阻尼，C＝aM+bK，a,b为绝对值小于1的正实数，振动控制方程为：假设激励F为多个强迫振动载荷组合而成i代表载荷的编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量。

满足飞行器翼肋近似的动力学响应即用下式来表示：

其中φ_j，ω_j为K，M所组成翼肋的阵型与圆频率，j为圆频率与阵型的编号，对应的第i个载荷对应于第j阶阵型的相位角，m_i，n_i代表选取p_i附近的圆频率对应的阵型编号的上界与下界；

$B_{\mathrm{ij}}=\frac{H_{\mathrm{ij}}}{\sqrt{{({\mathrm{\ω}}_j^2-p_j^2)}^2+{(a+b{\mathrm{\ω}}_j^2)}^2{p}_i^2}},$ $H_{\mathrm{ij}}=f_i^T{\mathrm{\φ}}_j.$

B_ij表示第i个载荷对应于第j阶阵型的引起的响应的振幅，H_ij表示第i个载荷的向量与第j阶阵型的内积；

所述步骤(6)的飞行器翼肋的变形能的求解过程如下：采用振型叠加法求得翼肋振动响应X，由飞行器翼肋结构质量阵M，飞行器翼肋结构刚度阵K与φ_j相互正交，通过对时间t的单变量在一个周期内进行优化即求出以下公式的最大值即翼肋的变形能：

其中φ_j，ω_j为K，M所组成翼肋的阵型与圆频率，ⁱ代表载荷的编号，_j为圆频率编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量，m_i，n_i代表附近的圆频率对应的阵型编号的上界与下界，B_ij表示第i个载荷对应于第j阶阵型的引起的响应的振幅，对应的第i个载荷对应于第j阶阵型的相位角；

所述步骤(7)中采用纽曼展开差分方法求解变形能对伪密度值的灵敏度的步骤如下：

对不同单元的伪密度值ρ_n进行摄动△ρ_n，计算刚度矩阵和质量矩阵的变化小量△K(ρ_n)、△M(ρ_n)，将动力学方程组转换为复数空间中的方程，利用复变函数求解；采用振型叠加法求解确定部分，采用纽曼展开法求解摄动部分，由于伪密度值ρ_n变化，振动系统的稳态响应有变化量△X(ρ_n)，则变形能有变化量△(XKX)＝2XK△X(ρ_n)+X△K(ρ_n)X，第一项采用正交关系求得，第二项采用稀疏计算求得；因此变形能对密度值的灵敏度为：

$\frac{\∂\max \left|X^T\mathrm{KX}\right|}{\∂{\mathrm{\ρ}}_n}=\frac{\max |2\mathrm{XK\ΔX}\left({\mathrm{\ρ}}_n\right)+\mathrm{X\ΔK}\left({\mathrm{\ρ}}_n\right)X|}{\mathrm{\Δ}{\mathrm{\ρ}}_n}.$

本发明与现有技术相比的优点在于：虽然现有方法可以实现对一些与静力学不相容的目标进行翼肋的优化设计，并且绕过局部模态问题的处理，但是最终得到的设计方案没有实际的工程意义。本发明考虑在周期载荷下的作用下对翼肋进行优化设计，以与静力学相容的最小化翼肋的变形能为设计目标，此设计目标有实际工程意义，可以保证翼肋的动态刚度最大化。采用振型叠加法与纽曼展开法花费少量计算量近似求解设计目标与相关灵敏度。提出局部模态问题处理方案，防止在优化设计中发生计算中断。降低飞行器结构质量，提高飞行器性能，而且得到设计方案所用计算量不大。

附图说明`

图1为本发明中飞行器翼肋有限元模型示意图；

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

具体实施方式

如图2所示，本发明实现步骤具体如下：

(1)导入飞行器翼肋初始几何模型或有限元模型，以及飞行条件下的周期激励动力载荷；

(2)将飞行器翼肋的初始几何模型离散为如示意图1所示的由有限元单元组成的有限元模型，如果已有有限元模型则直接进行下一步；

(3)确定飞行器翼肋结构中需要拓扑优化的部分，赋予所述步骤(1)中有限元模型中的有限元单元不同属性以区分需进行拓扑优化部分及不需要进行优化部分，再初始化飞行器翼肋有限元模型中需要进行优化部分的各个单元的伪密度值ρ_n，n为需要进行优化部分有限元单元编号，所述伪密度值ρ_n为0到1间连续的数值，0代表此处单元为空，1代表此处有单元，其它数值代表中间状态；

(4)基于伪密度值对有限元模型中需进行优化部分的有限元单元的刚度矩阵进行插值：k_n＝ρ_n ^mk_n ⁰，m>1，k_n ⁰为实体材料单元刚度矩阵，n为进行优化部分有限元单元编号，ρ_n代表进行优化部分单元伪密度值，m为惩罚因子，不需要优化部分的有限元单元为实体，单元刚度矩阵保持不变；选择惩罚因子的目的在于：对中间大小的伪密度值单元进行惩罚，尽量减少结构中间密度单元数目，使结构单元的伪密度值尽可能趋于0或1；当ρ_n＝1时表示这个单元是体材料ρ_n＝0时表示这个单元没有材料；当处于中间状态时可以解释为复合材料；m通常取为3，为避免矩阵奇异性选取ρ_min通常取ρ_min＝10^-3；由于在设计过程中，过多的低伪密度值会引起的局部模态问题，低伪密度值单元即最终不存在部分先发生振动，并且计算出的应变能与只计算高伪密度值的情况下相差很大，这使下面需要求解的变形能出现很大偏差，在插值完成后将基于伪密度值插值的飞行器翼肋有限元模型整体代入有限元求解器中进行模态分析，同时避免数值模态的发生有限元分析时不考虑的低伪密度值单元，只考虑高伪密度值单元，得到飞行器翼肋阵型与圆频率，当伪密度值大于0.1定义为高伪密度值，伪密度值不大于0.1定义为低伪密度值；

(5)基于步骤(4)模态分析得到的飞行器翼肋阵型与圆频率与飞行条件下的周期激励动力载荷，采用振型叠加法进行飞行器翼肋动力响应的近似求解；

当结构受到周期激励时，M为翼肋结构质量阵，K为翼肋结构刚度阵，F为激励载荷向量，X为系统随着时间t的响应；

当考虑含阻尼的振动系统时，振动控制方程为：

$\mathrm{KX}+C\stackrel{\·}{X}+M\stackrel{\·\·}{X}=F$

设阻尼C为瑞利阻尼，a,b一般为绝对值小于1的正实数：

C＝aM+bK

假设激励F为多个强迫振动载荷组合而成，i代表载荷的编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量；

$F=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_i\sin (p_{i}t+{\mathrm{\α}}_i)$

方程解也为对应载荷的解叠加而成，即:

${\mathrm{KX}}_i+C{\stackrel{\·}{X}}_i+M{\stackrel{\·\·}{X}}_i=f_i\sin (p_{i}t+{\mathrm{\α}}_i)$

$K\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{X}_i+C\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\stackrel{\·}{X}}_i+M\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\stackrel{\·\·}{X}}_i=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_i\sin (p_{i}t+{\mathrm{\α}}_i)$

设 $X_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{x}_{\mathrm{ij}}{\mathrm{\φ}}_j$

其中φ_j，ω_j为K，M所组成翼肋的阵型与圆频率，j为圆频率与阵型的编号，满足

$K{\mathrm{\φ}}_j={\mathrm{\ω}}_j^{2}M{\mathrm{\φ}}_j$

${\mathrm{\φ}}_j^{T}M{\mathrm{\φ}}_j=1,{\mathrm{\φ}}_j^{T}K{\mathrm{\φ}}_j={\mathrm{\ω}}_j^2$

${\mathrm{\φ}}_j^{T}M{\mathrm{\φ}}_j=O,{\mathrm{\φ}}_j^{T}K{\mathrm{\φ}}_j=O$

${\mathrm{KX}}_i+C{\stackrel{\·}{X}}_i+M{\stackrel{\·\·}{X}}_i=f_i\sin (p_{i}t+{\mathrm{\α}}_i)---\left(1\right)$

在方程(1)左乘

${{\mathrm{\φ}}_j^T\mathrm{KX}}_i+{\mathrm{\φ}}_j^{T}C{\stackrel{\·}{X}}_i+{\mathrm{\φ}}_j^{T}M{\stackrel{\·\·}{X}}_i={{\mathrm{\φ}}_j^{T}f}_i\sin (p_{i}t+{\mathrm{\α}}_i)$

${\stackrel{\·\·}{x}}_{\mathrm{ij}}+(a+b{\mathrm{\ω}}_j^2){\stackrel{\·}{x}}_{\mathrm{ij}}+{\mathrm{\ω}}_j^2{x}_{\mathrm{ij}}=H_{\mathrm{ij}}\sin (p_{i}t+{\mathrm{\α}}_i)$

$H_{\mathrm{ij}}={{\mathrm{\φ}}_j^{T}f}_i$

$B_{\mathrm{ij}}=\frac{H_{\mathrm{ij}}}{\sqrt{{({\mathrm{\ω}}_j^2-p_j^2)}^2+{(a+b{\mathrm{\ω}}_j^2)}^2{p}_i^2}},$

B_ij表示第i个载荷对应于第j阶阵型的引起的响应的振幅，H_ij表示第i个载荷的向量与第j阶阵型的内积，对应的第i个载荷对应于第j阶阵型的相位角；

现在，计算方法上可以用凝聚迭代的方法求解在激励频率p_i附近的圆频率与特征向量，而不用全部都同时求解出来，所以系统响应可以用振型叠加法近似表示为：

系统的响应可以用下式子表示：

m_i，n_i代表p_i附近的圆频率对应的阵型编号的上界与下界；m_i，n_i选取原则为：选取矩阵范数m_i，n_i使S₁＝m_i-1，S₂＝n_i+1，同时满足：

S₁，S₂代表m_i，n_i区间外的临近阵型编号，表示第i个载荷对应于第S₁，S₂阶阵型的引起的响应的振幅，对应的第i个载荷对应于第S₁，S₂阶阵型的相位角， 为K，M所组成翼肋的第S₁，S₂阶阵型；

(6)利用步骤(5)得到的飞行器翼肋的动力学的响应，来求解飞行器翼肋在动力载荷下的设计目标，即飞行器翼肋的变形能、约束函数的值；

采用振型叠加法求得翼肋振动响应X，由飞行器翼肋结构质量阵M，飞行器翼肋结构刚度阵K与φ_j相互正交，通过对时间t的单变量在一个周期内进行优化即求出以下公式的最大值即翼肋的变形能：

其中φ_j，ω_j为K，M所组成翼肋的阵型与圆频率，i代表载荷的编号，j为圆频率编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量，m_i，n_i代表附近的圆频率对应的阵型编号的上界与下界，B_ij表示第i个载荷对应于第j阶阵型的引起的响应的振幅，对应的第i个载荷对应于第j阶阵型的相位角；

(7)由于进行翼肋优化设计过程中需要灵敏度信息，采用纽曼展开差分方法求解步骤(6)中最小变形能对伪密度值的灵敏度，以及约束函数的灵敏度，同时低伪密度值单元的灵敏度由周围高伪密度值单元插值得到；

当载荷或者材料为摄动量时假设：

$M=M_0+\mathrm{\ΔM}=M_0+\underset{p=1}{\overset{q}{\mathrm{\Σ}}}{\mathrm{\δ}}_p{M}_p,{\mathrm{\φ}}_j^T{M}_0{\mathrm{\φ}}_j=1$

将系统引入人复数域：

$Z_i=f_i{e}^{i(p_{i}t+{\mathrm{\α}}_i)}$

${\mathrm{KY}}_i+C_i\stackrel{\·}{Y}+M{\stackrel{\·\·}{Y}}_i=Z_i$

$(K+\mathrm{iC}{p}_i-M{p}_i^2)Y_i=Z_i$

$\begin{array}{c}{Y}_i={(K+\mathrm{iC}{p}_i-M{p}_i^2)}^{-1}{Z}_i\\ ={(K_0+\mathrm{\ΔK}+i(a(M_0+\mathrm{\ΔM})+b(K_0+\mathrm{\ΔK}))p_i-(M_0+\mathrm{\ΔM})p_i^2)}^{-1}{Z}_i\\ ={(K^{\′}+\mathrm{\Δ}{K}^{\′})}^{-1}{Z}_i\end{array}$

$K^{\′}=K_0+i(a{M}_0+b{K}_0)p_i-M_0{p}_i^2$

$\mathrm{\Δ}{K}^{\′}=\mathrm{\ΔK}+i(\mathrm{a\ΔM}+\mathrm{b\ΔK})p_i-\mathrm{\Δ}{\mathrm{Mp}}_i^2$

Y_i＝(K'+△K')^-1Z_i＝(I+K'^-1△K')^-1K'^-1Z_i

纽曼展开后：

P＝K'^-1△K',(I+K'^-1△K')^-1＝(I-P+P²-P³....)

Y_i＝(I-P+P²-P³....)K'^-1Z_i

＝Y_i0-PY_i0+P²Y_i0-P³Y_i0

＝Y_i0-Y_i1+Y_i2-Y_i3...

PY_i(s-1)＝Y_is

Y_i0可由采用振型叠加法求解得到，同理可以按照阵型叠加方法求解低阶分量：

PY_i(s-1)＝K'^-1△K'Y_i(s-1)＝Y_is

K'Y_is＝△K'Y_i(s-1)

Y_i＝Y_i0-Y_i1+Y_i2-Y_i3..Y_is取Y_i复数部分即为原方程的解；

对伪密度值ρ_n进行摄动△ρ_n，刚度矩阵和质量矩阵都有△K(ρ_n)、△M(ρ_n)的小量变化，按照上述纽曼展开方法求出，由于伪密度值ρ_n变化，振动系统的稳态响应有变化量△X(ρ_n)，则变形能有变化量△(XKX)＝2XK△X(ρ_n)+X△K(ρ_n)X，第一项可以采用正交关系求得，第二项采用稀疏计算求得，因此，机翼应变能对密度值的灵敏度为：

$\frac{\∂\max \left|X^T\mathrm{KX}\right|}{\∂{\mathrm{\ρ}}_n}=\frac{\max |2\mathrm{XK\ΔX}\left({\mathrm{\ρ}}_n\right)+\mathrm{X\ΔK}\left({\mathrm{\ρ}}_n\right)X|}{\mathrm{\Δ}{\mathrm{\ρ}}_n}$

式中，△X(ρ_n)可以有一中的摄动公式快速求解，X、△X(ρ_n)、K都正交，△K(ρ_n)为稀疏矩阵，此步骤计算量非常小；

同时低伪密度值单元的灵敏度由周围高伪密度值单元插值得到；当伪密度值大于0.1定义为高伪密度值，伪密度值不大于0.1定义为低伪密度值；

(8)将步骤(7)计算得到的翼肋的变形能，约束函数的值以及灵敏度传递给优化器，以变形能最小为优化目标得到更新的伪密度值；

使用有限元离散后可以得到：

$\underset{u\∈U,E}{\mathrm{Min}}\max \left|X^T\mathrm{KX}\right|$

${\mathrm{st}.\mathrm{KX}}_i+C{\stackrel{\·}{X}}_i+M{\stackrel{\·\·}{X}}_i=f_i\sin (p_{i}t+{\mathrm{\α}}_i)$

$\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_i{v}_i<V_0$

$M=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_i$

(9)判断更新后得到的有限单元的当前伪密度值与前一次更新前的有限单元的伪密度值之间的差异，如果所有有限单元的伪密度值都小于设定的小量，则停止计算转到步骤(10)，否则返回到步骤(4)继续进行飞行器翼肋的动力学更新设计；

(10)设定阈值，伪密度值大于此阈值为实体，伪密度值小于阈值为空洞，将得到的拓扑结构的有限元模型导回为几何模型，确定飞行器翼肋最终的动力学方案。

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

本发明未详细阐述的部分属于本领域公知技术。

Claims (3)

1.一种基于振型叠加法的飞行器翼肋动力学确定方法，其特征在于实现步骤如下：

(1)导入飞行器翼肋初始几何模型或有限元模型，以及飞行条件下的周期激励动力载荷；

(2)将飞行器翼肋的初始几何模型离散为由有限元单元组成的有限元模型，如果已有有限元模型则直接进行下一步；

(3)确定飞行器翼肋结构中需要拓扑优化的部分，赋予所述步骤(2)中有限元模型中的有限元单元不同属性以区分需进行拓扑优化部分及不需要进行优化部分，再初始化飞行器翼肋有限元模型中需要进行优化部分的各个单元的伪密度值ρ_n，n为需要进行优化部分有限元单元编号，所述伪密度值ρ_n为0到1间连续的数值，0代表此处单元为空，1代表此处有单元，其它数值代表中间状态；

(4)基于伪密度值对有限元模型中需进行优化部分的有限元单元的刚度矩阵进行插值：k_n＝ρ_n ^mk_n ⁰m>1，其中n为需要进行优化部分有限元单元编号，m为惩罚因子，k_n ⁰为实体材料单元刚度矩阵,不需要优化部分的有限元单元为实体，单元刚度矩阵保持不变；然后将基于伪密度值插值的飞行器翼肋有限元模型整体代入有限元求解器中进行模态分析，同时避免数值模态的发生有限元分析时不考虑的低伪密度值单元，只考虑高伪密度值单元，得到飞行器翼肋振型与圆频率，当伪密度值大于0.1定义为高伪密度值，伪密度值不大于0.1定义为低伪密度值；

(5)基于步骤(4)模态分析得到的飞行器翼肋振型与圆频率与飞行条件下的周期激励动力载荷，采用振型叠加法进行飞行器翼肋动力响应的近似求解；

(6)利用步骤(5)得到的飞行器翼肋的动力学的响应，来求解飞行器翼肋在动力载荷下的设计目标，即飞行器翼肋的变形能、约束函数的值；

(7)进行飞行器翼肋优化过程中需要灵敏度信息，采用纽曼展开差分方法求解步骤(6)中变形能对伪密度值的灵敏度，以及约束函数的灵敏度，同时低伪密度值单元的灵敏度由周围高伪密度值单元插值得到；

(8)将步骤(7)计算得到的飞行器翼肋的变形能，约束函数的值以及灵敏度传递给优化器，以变形能最小为优化目标得到更新的伪密度值；

(9)判断更新后得到的有限单元的当前伪密度值与前一次更新前的有限单元的伪密度值之间的差异，如果所有有限单元的伪密度值的差异都小于设定的小量，则停止计算转到步骤(10)，否则返回到步骤(4)继续进行飞行器翼肋的动力学更新设计；

(10)设定阈值，伪密度值大于此阈值为实体，伪密度值小于阈值为空洞，将得到的拓扑结构的有限元模型导回为几何模型，确定飞行器翼肋最终的动力学方案；

所述步骤(5)用振型叠加法进行翼肋动力响应的近似求解公式：

当结构受到周期激励时，M为飞行器翼肋结构质量阵，K为飞行器翼肋结构刚度阵，F为激励载荷向量，X为系统随着时间t的响应，阻尼C为瑞利阻尼，C＝aM+bK，a,b为绝对值小于1的正实数，振动控制方程为：假设激励F为多个强迫振动载荷组合而成i代表载荷的编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量；

满足飞行器翼肋近似的动力学响应即用下式来表示：

其中j为圆频率与振型的编号，对应的第i个载荷对应于第j阶振型的相位角，m_i，n_i代表选取p_i附近的圆频率对应的振型编号的下界与上界；

$B_{\mathrm{ij}}=\frac{H_{\mathrm{ij}}}{\sqrt{{({\mathrm{\&omega;}}_j^2-p_i^2)}^2+{(a+{\mathrm{b\&omega;}}_j^2)}^2{p}_i^1}},$ H_ij＝f_i ^Tφ_j

B_ij表示第i个载荷对应于第j阶振型的引起的响应的振幅，H_ij表示第i个载荷的向量与第j阶振型的内积，φ_j，ω_j为K，M所组成翼肋的振型与圆频率。

2.根据权利要求1所述的一种基于振型叠加法的飞行器翼肋动力学确定方法，其特征在于：所述步骤(6)的飞行器翼肋的变形能的求解过程如下：采用振型叠加法求得翼肋振动响应X，由飞行器翼肋结构质量阵M，飞行器翼肋结构刚度阵K与φ_j相互正交，通过对时间t的单变量在一个周期内进行优化即求出以下公式的最大值即翼肋的变形能：

其中φ_j，ω_j为K，M所组成翼肋的振型与圆频率，i代表载荷的编号，j为圆频率编号，k代表载荷数目，p_i，α_i，f_i为对于第i个载荷的激励圆频率，相位角与载荷的向量，m_i，n_i代表选取p_i附近的圆频率对应的振型编号的下界与上界，B_ij表示第i个载荷对应于第j阶振型的引起的响应的振幅，对应的第i个载荷对应于第j阶振型的相位角。

3.根据权利要求1所述的一种基于振型叠加法的飞行器翼肋动力学确定方法，其特征在于：所述步骤(7)中采用纽曼展开差分方法求解变形能对伪密度值的灵敏度的步骤如下：

对不同单元的伪密度值ρ_n进行摄动Δρ_n，计算刚度矩阵和质量矩阵的变化小量ΔK(ρ_n)、ΔM(ρ_n)，将动力学方程组转换为复数空间中的方程，利用复变函数求解；采用振型叠加法求解确定部分，采用纽曼展开法求解摄动部分，由于伪密度值ρ_n变化，振动系统的稳态响应有变化量ΔX(ρ_n)，则变形能有变化量Δ(X^TKX)＝2X^TKΔX(ρ_n)+X^TΔK(ρ_n)X，第一项采用正交关系求得，第二项采用稀疏计算求得；因此变形能对伪密度值的灵敏度为：

$\frac{\&part;max\left|X^{T}KX\right|}{\&part;{\mathrm{\&rho;}}_n}=\frac{max|2X^{T}K\mathrm{\&Delta;}X\left({\mathrm{\&rho;}}_n\right)+X^T\mathrm{\&Delta;}K\left({\mathrm{\&rho;}}_n\right)X|}{{\mathrm{\&Delta;\&rho;}}_n}.$