CN112464372B - 一种弹性机翼副翼舵面效率的设计敏度工程数值方法 - Google Patents

Abstract

一种弹性机翼副翼舵面效率的设计敏度工程数值方法，包括如下步骤：建立气动面元模型和结构有限元模型，并计算气动面元模型和结构有限元模型间数据传递矩阵；计算刚性机翼单位副翼偏转角产生的滚转速率；计算弹性机翼单位副翼偏转角产生的驱动力矩和单位滚转速率产生的阻尼力矩；进行弹性机翼副翼舵面效率及其设计敏度的计算。本发明将弹性机翼单位副翼偏转角产生的滚转速率表示为单位副翼偏转角产生的驱动力矩和单位滚转速率产生的阻尼力矩的组合的形式，解决了滚转速率形式定义的副翼舵面效率的设计敏度难以解析求解的问题，在此基础上，通过直接法或伴随法对副翼舵面效率设计敏度的工程计算，从而辅助后续飞机结构的精细数值优化设计。

Description

一种弹性机翼副翼舵面效率的设计敏度工程数值方法

技术领域

本发明适用于飞机结构的数值优化设计领域，具体涉及一种弹性机翼副翼舵面效率关于结构设计变量敏度的工程数值方法。

背景技术

施加于弹性机翼上的气动载荷会改变机翼的变形，并反过来改变机翼的气流特征，从而产生气动弹性现象并影响最终稳定变形状态下机翼的气动载荷。副翼用于控制飞机的飞行状态，随着飞行速度的增大，副翼舵面效率会不断降低，当达到反效速度时，副翼偏转将对飞行状态不产生任何影响。副翼舵面效率的大小反映了飞机的机动性能，常作为约束条件用于飞机结构的数值优化设计，弹性机翼副翼舵面效率的设计敏度可为梯度优化算法提供必要的指导。在实际工程中一般有两种副翼舵面效率的定义形式：机翼根部固支情况通过弹性机翼和刚性机翼的翼根弯矩之比进行定义，机翼稳定滚转情况通过弹性机翼和刚性机翼单位副翼偏转角产生的滚转速率之比进行定义。前者是一种较简单的形式，未能考虑副翼对飞机滚转运动的实际作用效果，因此，本发明中副翼舵面效率采用第二种定义形式。如文献“Wright JR,Cooper JE.Introduction to AircraftAeroelasticity andLoads.2nd ed.NewYork:Wiley,2015”所述，目前常用的副翼舵面效率计算方法是依据机翼气动弹性原理，利用简单的工程算法公式以及粗定量技术参数进行求解。其计算结果常作为一种定性的参考，并且该方法也不利于后续设计敏度的解析计算，这些缺陷使其难以应用于飞机结构的精细数值优化设计。

发明内容

为克服滚转速率形式定义的副翼舵面效率设计敏度难以解析求解的缺点，本发明中，分别采用势流理论面元法和有限元法建立机翼的气动模型和结构模型，在此基础上，进行弹性机翼副翼舵面效率及其设计敏度的计算。由前面的叙述可知，副翼舵面效率设计敏度求解的关键是弹性机翼单位副翼偏转角产生的滚转速率的设计敏度，但滚转速率和结构设计变量之间的关系难以直接显式确定，因此需要对副翼舵面效率的形式进一步转化。当飞行速度小于副翼反效速度时，副翼偏转产生的驱动力矩将使机翼绕其纵轴滚转，同时机翼的滚转运动将产生阻尼力矩，且阻尼力矩随着滚转速率的增大而增大。另外，阻尼力矩有使机翼滚转速率降低的趋势。当滚转速率达到某一值时，副翼偏转产生的驱动力矩和滚转运动产生的阻尼力矩之和为零，这时机翼达到稳定滚转状态。特别地，对于弹性机翼，弹性变形也将引起驱动力矩和滚转力矩的变化，从而使弹性机翼具有和刚性机翼不同的滚转速率。

基于以上思路，本发明将弹性机翼单位副翼偏转角产生的滚转速率转化为单位副翼偏转角产生的驱动力矩和单位滚转速率产生的阻尼力矩的组合形式。设计敏度工程数值方法是基于离散化的气动模型和结构模型实现的。驱动力矩或阻尼力矩为气动载荷和模型节点坐标的函数，而弹性机翼的气动载荷通过位移间接依赖于结构设计变量，由此即可通过解析法方式计算副翼舵面效率的设计敏度。具体说明如下，

本发明提出一种弹性机翼副翼舵面效率的设计敏度工程数值方法，包括如下步骤：

1)建立气动面元模型和结构有限元模型，并计算气动面元模型和结构有限元模型间数据传递矩阵；

2)利用步骤1)的计算结果，计算刚性机翼单位副翼偏转角产生的滚转速率；

3)利用步骤1)的计算结果，计算弹性机翼单位副翼偏转角产生的驱动力矩；

4)利用步骤1)的计算结果，计算弹性机翼单位滚转速率产生的阻尼力矩；

5)根据步骤2)、3)、4)的计算结果，进行弹性机翼副翼舵面效率及其设计敏度的计算。

优选地，所述步骤3)，具体包括如下步骤：

31)在某一副翼偏转角下计算刚性机翼副翼偏转角产生的驱动力，将其作为初始气动载荷施加于结构有限元模型；

32)通过松耦合迭代的方式计算最终稳定的弹性变形和气动载荷，以及弹性变形状态下的驱动力矩；

33)获得弹性机翼单位副翼偏转角产生的驱动力矩。

优选地，所述步骤4)，具体包括如下步骤：

41)在选定的滚转速率下计算刚性机翼滚转运动产生的阻尼力，将其作为初始气动载荷施加于结构有限元模型；

42)通过松耦合迭代的方式计算最终稳定的弹性变形和气动载荷，以及弹性变形状态下的阻尼力矩；

43)获得弹性机翼单位滚转速率产生的阻尼力矩。

优选地，所述步骤5)中弹性机翼副翼舵面效率的计算，通过以下公式实现：

式中，η为副翼舵面效率；ω^r,β为刚性机翼单位副翼偏转角产生的滚转速率；β为副翼偏转角；f_sezβ为弹性机翼由β产生的施加于结构模型的气动载荷；Δω为选定的滚转速率；f_sezΔω为弹性机翼由Δω产生的施加于结构模型的气动载荷；y^T为有限元模型节点坐标y分量组成的行向量。

优选地，所述步骤5)中是通过直接法或伴随法进行弹性机翼副翼舵面效率的设计敏度计算。

优选地，所述通过直接法进行弹性机翼副翼舵面效率的设计敏度计算中，首先通过解线性方程组的方式计算弹性机翼在副翼偏转角β下对应的平衡状态的结构位移列阵的设计敏度du_sβ/db，以及弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构位移列阵的设计敏度du_sΔω/db，然后得：

式中，b为设计变量；u_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构位移列阵；f_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构气动载荷列阵；u_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构位移列阵；f_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构气动载荷列阵；B₁和B₂是布尔矩阵，其功能是分别按照f_sezβ和f_sezΔω中的元素在f_sβ和f_sΔω中的位置由df_sβ/db和df_sΔω/db得到df_sezβ/db和df_sezΔω/db；最终副翼舵面效率关于设计变量b的敏度可表示为：

优选地，所述通过伴随法进行弹性机翼副翼舵面效率的设计敏度计算中，首先计算以下两个伴随变量：

式中，

为y^T·df_sezβ/du_sβ对应的伴随变量；

为y^T·df_sezΔω/du_sΔω对应的伴随变量；K为结构总体刚度矩阵，最终副翼舵面效率关于设计变量b的敏度可表示为：

本发明的有益效果为：将弹性机翼单位副翼偏转角产生的滚转速率表示为单位副翼偏转角产生的驱动力矩和单位滚转速率产生的阻尼力矩的组合的形式，解决了滚转速率形式定义的副翼舵面效率的设计敏度难以解析求解的问题。在此基础上，实现了直接法或伴随法对副翼舵面效率设计敏度的工程计算，从而辅助后续飞机结构的精细数值优化设计。

附图说明

图1是本发明的算法流程示意图；

图2是机翼面元模型图；

图3是机翼整体有限元模型；

图4是机翼内部骨架有限元模型；

图5是副翼舵面效率随飞行速度变化的曲线；

图6是副翼舵面效率设计敏度随飞行速度变化的曲线；

图7是直接法的计算时间随设计变量数目变化的曲线。

图中标注说明如下：

1是机翼面元模型中副翼对应的子块；2是结构有限元模型第1个区域；3是结构有限元模型第2个区域；4是第1个固定点；5是第2个固定点；6是飞行速度，单位为m/s；7是副翼舵面效率；8是副翼舵面效率设计敏度；9是第1个设计变量对应的曲线；10是第2个设计变量对应的曲线；11是第3个设计变量对应的曲线；12是第4个设计变量对应的曲线；13是设计变量数目；14是计算时间，单位为s。

具体实施方式

设机翼所在的直角坐标系为Oxyz，其中机翼的弦平面位于x-y平面内，x轴沿机翼的弦向，y轴沿机翼的展向，z轴垂直于x-y平面向上。副翼舵面效率的具体定义形式为：

式中，η为副翼舵面效率；ω^e,β为弹性机翼单位副翼偏转角产生的滚转速率；ω^r,β为刚性机翼单位副翼偏转角产生的滚转速率。

式(1)中的副翼舵面效率是直接通过滚转速率形式进行定义的，当结构设计变量的值变化时，弹性机翼的滚转速率也将随之发生变化。为能计算滚转速率关于设计变量的敏度，需按照力学原理将式(1)转化为可计算表达的形式。

在稳定滚转状态下，机翼绕x轴的总滚转力矩为零，即：

(M^r,β+ΔM^e,β)·β+(M^r,ω+ΔM^e,ω)·ω＝0 (2)

式中，β为副翼偏转角；ω为副翼偏转角β引起的弹性机翼滚转速率；M^r,β为刚性机翼单位副翼偏转角产生的滚转驱动力矩；M^r,ω为刚性机翼单位滚转速率产生的滚转阻尼力矩；ΔM^e,β为弹性机翼单位副翼偏转角产生的驱动力矩相对于M^r,β的变化量；ΔM^e,ω为弹性机翼单位滚转速率产生的阻尼力矩相对于M^r,ω的变化量。对于弹性机翼，可认为其弹性变形包括由副翼偏转所产生的弹性变形部分和由滚转运动所产生的弹性变形部分，ΔM^e,β反映了单位副翼偏转角引起的弹性变形部分所产生的滚转驱动力矩变化量，ΔM^e,ω反映了单位滚转速率变化引起的弹性变形部分所产生的滚转阻尼力矩的变化量。

根据式(2)，可将式(1)中的ω^e,β写为：

式中，f_sezβ为弹性机翼由β产生的施加于结构模型的气动载荷；f_sezω为弹性机翼由ω产生的施加于结构模型的气动载荷；y^T为有限元模型节点坐标y分量组成的行向量。当机翼的弹性变形较小时，可认为弹性机翼由滚转速率引起的施加于结构模型的气动载荷随滚转速率的增大而线性变化。所以，可将弹性机翼副翼舵面效率写为以下形式：

式中，Δω为选定的滚转速率；f_sezΔω为弹性机翼由Δω产生的施加于结构模型的气动载荷。由于原始机翼中弧面和迎角产生的气动载荷左右对称，其气动载荷对机翼的滚转力矩贡献为零，因此在计算弹性机翼的气动载荷时直接取机翼的迎角为零，并忽略原始机翼中弧面对下洗的影响。本发明中的副翼舵面效率设计敏度的工程分析方法流程示意图如图1所示。

本发明的具体实施步骤是：

步骤1：建立机翼的气动模型和结构模型，并计算气动模型和结构模型之间的数据传递矩阵。通过求解定常线性小扰动速度势方程可获得气动力影响系数矩阵，进而得到机翼的气动载荷，其一般形式的公式为：

f_a＝q·S·AIC·DW (5)

式中，f_a为面元气动载荷列阵；q为飞行动压；S为各面元单元的面积组成的对角矩阵；AIC为气动力影响系数矩阵；DW为机翼下洗列阵。

采用有限元法建立机翼的结构模型，并选择合适的根部固支节点使机翼无刚体运动，其平衡方程可表示为：

K·u_s＝f_s (6)

式中，K为结构总体刚度矩阵；u_s为结构位移列阵；f_s为结构气动载荷列阵。

式(5)中，机翼下洗列阵DW的元素为单位速度来流在相应的面元控制点处的法向速度分量，即由中弧面处的零法向流动边界条件确定。中弧面为机翼结构模型上、下表面节点垂向位移得到的中间位置曲面，显然f_a是结构有限元模型节点位移的函数。此外，f_s由f_a传递至结构模型获得，结构位移列阵u_s需要通过式(6)解线性方程组获得。所以，弹性机翼的气动载荷由松耦合迭代的方式确定。在该过程中，需要气动模型和结构模型之间的数据传递，一方面需要将结构模型上、下表面节点垂向位移传递至面元控制点以计算下洗列阵，另一方面则需要将面元气动载荷传递至结构有限元模型。在本发明中，两模型之间的数据传递通过径向基函数插值的方法实现，首先进行位移传递矩阵的计算，然后根据虚功原理即可获得载荷传递矩阵，数据传递矩阵中的元素是气动模型和结构模型节点坐标的函数，与设计变量无关。

步骤2：计算刚性机翼单位副翼偏转角产生的滚转速率。副翼偏转角β产生的施加在面元模型上的气动载荷可表示为：

f_arβ＝q·S·AIC·DW_rβ (7)

式中，f_arβ为刚性机翼在副翼偏转角β下产生的面元气动载荷列阵；DW_rβ为刚性机翼在副翼偏转角β下对应的下洗列阵。

设待求解的刚性机翼滚转速率为ω^r，则ω^r对应的下洗列阵DW_rω可表示为：

式中，y_a为面元控制点的展向坐标组成的列阵；V为来流速度。

进而刚性机翼在滚转速率ω^r下产生的面元气动载荷列阵f_aωr可表示为：

f_arω＝q·S·AIC·DW_rω (9)

在稳定滚转状态下，f_arβ产生的驱动力矩和f_arω产生的阻尼力矩之和为0，此时仅ω^r未知，据此即可解得刚性机翼的滚转速率。进一步，得

步骤3：计算弹性机翼单位副翼偏转角产生的滚转驱动力矩。该过程中不考虑副翼偏转产生的滚转运动，对于弹性机翼，以式(7)给出的刚性机翼在副翼偏转角β下产生的气动载荷作为初始载荷，将其传递至结构有限元模型计算弹性变形，然后求解因弹性变形产生的气动载荷，并与f_arβ一同施加于有限元模型，如此反复迭代计算直至达到平衡状态，这时有：

K·u_sβ＝f_sβ (11)

式中，K为结构总体刚度矩阵；u_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构位移列阵；f_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构气动载荷列阵。进而可得到弹性机翼单位副翼偏转角所产生的滚转驱动力矩，即式(4)中的y^T·f_sezβ/β，其中f_sezβ可通过提取f_sβ中相应的元素获得。

步骤4：计算弹性机翼单位滚转速率产生的滚转阻尼力矩。对于选定的滚转速率Δω，则初始的刚性机翼滚转运动产生的面元气动载荷为：

f_aΔω＝q·S·AIC·DW_Δω (12)

式中，f_aΔω为刚性机翼在选定的滚转速率Δω下产生的面元气动载荷列阵；DW_Δω为Δω对应的下洗列阵。以f_aΔω作为初始气动载荷，将其传递至结构有限元模型计算弹性变形，然后求解因弹性变形产生的气动载荷，并与f_aΔω一同作为外力施加于有限元模型，如此反复迭代计算直至达到平衡状态，则该情况下的结构平衡方程可表示为：

K·u_sΔω＝f_sΔω (13)

式中，K为结构总体刚度矩阵；u_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构位移列阵；f_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构气动载荷列阵。进而可得到弹性机翼单位滚转速率引起的滚转阻尼力矩，即式(4)中的y^T·f_sezΔω/Δω，其中f_sezΔω可通过提取f_sΔω中相应的元素获得。

步骤5：副翼舵面效率设计敏度计算的基本原理。根据式(4)，副翼舵面效率关于设计变量b的敏度可表示为：

由此可知，副翼舵面效率设计敏度求解的关键是计算f_sezβ和f_sezΔω的设计敏度，而f_sezβ和f_sezΔω分别与f_sβ和f_sΔω直接关联，所以，需要首先求解f_sβ和f_sΔω的设计敏度。式(11)和式(13)关于设计变量b求导可得：

dK/db·u_sβ和dK/db·u_sΔω为列阵，由于设计变量b对刚度矩阵的影响仅限于与该设计变量相关的单元，因此在计算中不必存储矩阵dK/db，只需要将相关联的单元刚度矩阵的敏度与其相应的节点位移相乘再叠加即可。这种方式可以减少计算机内存的消耗并提高计算效率。f_sβ和f_sΔω与结构模型上、下表面节点垂向位移具有直接显式关系，所以df_sβ/du_sβ和df_sΔω/du_sΔω可通过气动载荷和位移之间的关系解析推导获得。f_sβ和f_sΔω通过结构模型节点位移而间接依赖于设计变量，若由式(15)和式(16)解得du_sβ/db和du_sΔω/db，即可根据链式法则得到f_sβ和f_sΔω的设计敏度，进而得到副翼舵面效率的设计敏度。副翼舵面效率设计敏度的计算有两种方法，即直接法和伴随法，详细过程将在步骤6和步骤7进行说明。

步骤6：副翼舵面效率设计敏度的直接法。在利用直接法进行设计敏度的解析求解时，首先由式(15)和式(16)通过解线性方程组得到du_sβ/db和du_sΔω/db，进而得

在式(17)和式(18)中，B₁和B₂是布尔矩阵，其功能是分别按照f_sezβ和f_sezΔω中的元素在f_sβ和f_sΔω中的位置由df_sβ/db和df_sΔω/db得到df_sezβ/db和df_sezΔω/db。将式(17)和式(18)的计算结果代入式(14)即可获得副翼舵面效率的设计敏度。若结构设计变量数目为N_dv，则在直接法中，为了求解位移关于设计变量的敏度，需要进行N_dv×2次线性方程组的求解。显然，随着设计变量数目的增大，该部分的计算量也将越来越大，从而影响设计敏度计算的效率。

步骤7：副翼舵面效率设计敏度的伴随法。在伴随法中，需要引入两个伴随变量，即：

式中，

为y^T·df_sezβ/du_sβ对应的伴随变量；

为y^T·df_sezΔω/du_sΔω对应的伴随变量。以[K-df_sβ/du_sβ]^T作为线性方程组的系数矩阵，以-[df_sezβ/du_sβ]·y作为右端项，可解得

的转置；同理，以[K-df_sΔω/du_sΔω]^T作为系数矩阵，以-[df_sezΔω/du_sΔω]·y作为右端项，可解得

的转置。最终，使用伴随法计算副翼舵面效率设计敏度的公式可表示为：

相对于直接法，因为伴随变量与设计变量无关，所以对不同数目的设计变量，

和

仅需要求解一次，所以当设计变量数目多的时候，采用伴随法进行副翼舵面效率设计敏度的分析是一种更高效的方式。

实施例

本发明以某一大展弦比机翼作为具体实施例，对副翼舵面效率设计敏度的工程数值方法进行验证和说明。该机翼模型的几何参数如下：翼根弦长为2140mm，翼尖弦长为440mm，机翼的半展长为8900mm，前缘后掠角为9°。

步骤1：建立大展弦比机翼的气动模型和结构模型，并计算气动模型和结构模型之间的数据传递矩阵。面元模型的几何外形是刚性机翼在x-y平面的投影，如图2所示，该面元模型由6个面元子块组成，共包含气动单元802个。大气压强设置为101325Pa，飞行马赫数设置为0.6。

采用Hypermesh和Patran建立结构模型，如图3和图4所示，有限元模型共包含有限元节点746个、四边形壳单元904个、三角形壳单元86个、梁单元400个。机翼蒙皮为复合材料对称层合板[±45/0₄/±45/90₂]_s，其单层属性如下：纵向弹性模量为1.40×10¹¹Pa，横向弹性模量为8.60×10⁹Pa，剪切模量为5.00×10⁹Pa，泊松比为0.33，单层厚度为0.125mm。其它壳单元为钛合金材料，其属性如下：弹性模量为1.10×10¹¹Pa，泊松比为0.3，壳单元厚度为3.0mm。该模型中有4个设计变量，均表示机翼蒙皮复合材料层合板0°铺层的厚度，其中第1个设计变量在结构有限元模型第1个区域的上表面，第2个设计变量在结构有限元模型第2个区域的上表面，第3个设计变量在结构有限元模型第1个区域的下表面，第4个设计变量在结构有限元模型第2个区域的下表面。

使用径向基函数插值方法计算数据传递矩阵时，以Wendland C2函数作为径向基函数，并取紧支半径为1500mm。在弹性机翼气动载荷计算过程中，为建立机翼中弧面，需要分别将结构模型上、下表面节点垂向位移传递至面元控制点，进而计算机翼下洗列阵。由于气动面元模型是二维的，而机翼结构模型是三维的，所以在实施例中将气动载荷施加于有限元模型的下表面节点。

步骤2：计算刚性机翼单位副翼偏转角产生的滚转速率。取副翼偏转角为1°，得刚性机翼的滚转驱动力矩为6904.97N·m，最终解得刚性机翼单位副翼偏转角(单位为°)产生的滚转速率为3.61°/s。

步骤3：计算弹性机翼单位副翼偏转角产生的驱动力矩。取副翼偏转角为1°，当弹性机翼前后两次迭代步最大翼尖垂向位移的相对误差小于0.1％时，则认为迭代收敛，最终得弹性机翼单位副翼偏转角(°)产生的滚转驱动力矩为2954.20N·m。

步骤4：计算弹性机翼单位滚转速率产生的阻尼力矩。选定的滚转速率为1°/s，当弹性机翼前后两次迭代步最大翼尖垂向位移的相对误差小于0.1％时，则认为迭代收敛，最终得弹性机翼单位滚转速率(单位为°/s)产生的滚转阻尼力矩为-1519.88N·m。

步骤5：根据式(4)得马赫数0.6时的副翼舵面效率为：1/3.61×2954.20/1519.88＝53.84％。如图5所示，给出了副翼舵面效率随飞行速度变化的曲线，副翼舵面效率随着飞行速度的增大而减小，直至为零，此时发生副翼反效。由图可知，副翼反效速度约为315m/s。

步骤6，副翼舵面效率的设计敏度计算。分别采用直接法和伴随法进行副翼舵面效率的设计敏度计算，然后使用中心差分法对其计算结果的正确性进行验证。将设计变量b增大或减小较小的值Δb，可得到修改后的设计变量b+Δb或b-Δb。基于修改后的设计变量所对应的结构有限元模型进行副翼舵面效率的计算，之后即可得到副翼舵面效率设计敏度的中心差分结果，在本实施例中，Δb取为0.01mm。表1给出了直接法、伴随法和中心差分法的副翼舵面效率设计敏度的计算结果，其中直接法和伴随法的结果完全相同，中心差分法与直接法或伴随法的结果高度吻合。

由于后掠机翼的弯扭耦合效应，弹性机翼将产生低头扭转变形，从而使气动载荷减小，所以最终副翼舵面效率小于1。设计变量值的增大将使机翼的刚度增大，进一步使弹性机翼的低头扭转变形程度减小，最终引起副翼舵面效率的增大。所以本实施例中4个设计变量所对应的副翼舵面效率设计敏度均为正值。

图6给出了副翼舵面效率设计敏度随飞行速度变化的曲线。相对于靠近翼尖的设计变量，靠近翼根的设计变量的摄动将会对弹性变形产生更大的影响，所以在相同的飞行速度下，副翼舵面效率关于第1个设计变量的敏度大于关于第2个设计变量的敏度，同理，副翼舵面效率关于第3个设计变量的敏度大于关于第4个设计变量的敏度。结构有限元模型的第1个区域或第2个区域的两个设计变量对应的副翼舵面效率设计敏度不同，这是由结构模型的不对称性引起的。

直接法中位移的设计敏度和伴随法中的伴随变量都通过解线性方程组的方式获得，在该过程中，需要首先对线性方程组的系数矩阵进行LU分解，然后前后回代得到待求解项。相对于矩阵分解，前后回代的时间非常小。虽然分解的系数矩阵可以重复使用，但是随着设计变量数目的增多，直接法中由前后回代引起的时间花费也将越来越大。为了对直接法的计算时间进行测试，以本实施例中的4个设计变量作为一个设计变量组，然后不断地将这4个设计变量添加到该设计变量组并进行计算时间的测试。这里仅考虑前后回代的计算时间，图7给出了直接法的计算时间随设计变量数目变化的曲线。所以，在设计变量数目多的情况下，伴随法比直接法具有更高的计算效率。

以上实施例的计算结果说明了本发明中弹性机翼副翼舵面效率设计敏度的工程数值方法的有效性。

表1副翼舵面效率设计敏度的计算结果

	直接法	伴随法	中心差分法
				第1个设计变量	3.1399E-2	3.1399E-2	3.1531E-2
第2个设计变量	2.4340E-2	2.4340E-2	2.4287E-2
				第3个设计变量	2.1975E-2	2.1975E-2	2.2104E-2
第4个设计变量	1.8574E-2	1.8574E-2	1.8516E-2

Claims (4)

1.一种弹性机翼副翼舵面效率的设计敏度工程数值方法，其特征在于，包括如下步骤：

1)建立气动面元模型和结构有限元模型，并计算气动面元模型和结构有限元模型间数据传递矩阵；

2)利用步骤1)的计算结果，计算刚性机翼单位副翼偏转角产生的滚转速率，通过以下公式实现：

式中，刚性机翼滚转速率为ω^r，β为副翼偏转角；

3)利用步骤1)的计算结果，计算弹性机翼单位副翼偏转角产生的驱动力矩，具体包括如下步骤：

31)在某一副翼偏转角下计算刚性机翼副翼偏转角产生的驱动力，将其作为初始气动载荷施加于结构有限元模型；

32)通过松耦合迭代的方式计算最终稳定的弹性变形和气动载荷，以及弹性变形状态下的驱动力矩；

33)获得弹性机翼单位副翼偏转角产生的驱动力矩；

4)利用步骤1)的计算结果，计算弹性机翼单位滚转速率产生的阻尼力矩，具体包括如下步骤：

41)在选定的滚转速率下计算刚性机翼滚转运动产生的阻尼力，将其作为初始气动载荷施加于结构有限元模型；

42)通过松耦合迭代的方式计算最终稳定的弹性变形和气动载荷，以及弹性变形状态下的阻尼力矩；

43)获得弹性机翼单位滚转速率产生的阻尼力矩；

5)根据步骤2)、3)、4)的计算结果，进行弹性机翼副翼舵面效率及其设计敏度的计算，通过以下公式实现：

式中，η为副翼舵面效率；ω^r,β为刚性机翼单位副翼偏转角产生的滚转速率；β为副翼偏转角；f_sezβ为弹性机翼由β产生的施加于结构模型的气动载荷；Δω为选定的滚转速率；f_sezΔω为弹性机翼由Δω产生的施加于结构模型的气动载荷；y^T为有限元模型节点坐标y分量组成的行向量。

2.根据权利要求1所述的一种弹性机翼副翼舵面效率的设计敏度工程数值方法，其特征在于，所述步骤5)中是通过直接法或伴随法进行弹性机翼副翼舵面效率的设计敏度计算。

3.根据权利要求2所述的一种弹性机翼副翼舵面效率的设计敏度工程数值方法，其特征在于，所述通过直接法进行弹性机翼副翼舵面效率的设计敏度计算中，首先通过解线性方程组的方式计算弹性机翼在副翼偏转角β下对应的平衡状态的结构位移列阵的设计敏度du_sβ/db，以及弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构位移列阵的设计敏度du_sΔω/db，然后得：

式中，b为设计变量；u_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构位移列阵；f_sβ为弹性机翼在副翼偏转角β下对应的平衡状态的结构气动载荷列阵；u_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构位移列阵；f_sΔω为弹性机翼在选定的滚转速率Δω下对应的平衡状态的结构气动载荷列阵；B₁和B₂是布尔矩阵，其功能是分别按照f_sezβ和f_sezΔω中的元素在f_sβ和f_sΔω中的位置由df_sβ/db和df_sΔω/db得到df_sezβ/db和df_sezΔω/db；最终副翼舵面效率关于设计变量b的敏度表示为：

4.根据权利要求2所述的一种弹性机翼副翼舵面效率的设计敏度工程数值方法，其特征在于，所述通过伴随法进行弹性机翼副翼舵面效率的设计敏度计算中，首先计算以下两个伴随变量：

式中，

为y^T·df_sezβ/du_sβ对应的伴随变量，

为y^T·df_sezΔω/du_sΔω对应的伴随变量，K为结构总体刚度矩阵，最终副翼舵面效率关于设计变量b的敏度表示为：

Images