CN110704953A - 一种大展弦比机翼静气弹性能设计敏度的分析方法 - Google Patents

Abstract

本发明利用一次函数与一致紧支径向基函数建立最小二乘拟合曲面函数，从而实现对机翼弹性变形的模拟，解决了传统插值方法进行位移传递时存在的变形曲面波动性问题。在静气弹计算的基础上，提出基于迭代方式的气动载荷设计敏度的模型解析求解方法，其中包括直接方法和伴随方法，该算法避免了对原始刚度矩阵的破坏，同时具有计算的高效性，符合工程实际的需求。

Description

一种大展弦比机翼静气弹性能设计敏度的分析方法

技术领域

本发明适用于飞行器结构静气弹优化设计领域，具体涉及一种大展弦比机翼静气弹性能设计敏度的分析方法。

背景技术

现代飞行器结构轻量化和大柔性的特点，使静气弹效应愈加突出。随着飞行器结构优化设计精细化程度的提高，在优化过程中必须考虑静气弹的影响。静气弹性能设计敏度指结构尺寸或参数的变化所引起的静气弹性能指标的变化程度，它将为梯度优化算法提供必要的信息并指导优化过程的进行，静气弹性能设计敏度的计算是飞机结构优化设计中的重要工作环节。设计敏度分析方法可分为三类：近似方法、离散方法和连续方法。常用的近似方法有有限差分法和复变量法；离散方法采用对离散化的控制方程求导的方式计算设计敏度；连续方法则是在对控制方程进行离散化前对其进行求导运算。离散方法和连续方法是设计敏度的解析计算方法，可通过直接方法或伴随方法实现。在本发明中，为实现结构静气弹性能的高效率计算，分别采用有限元法和基于小扰动势流方程的面元法建立结构模型和气动模型，并使用CFD数据对气动模型进行修正。因此，静气弹计算与设计敏度分析是在离散化的控制方程基础上开展的。

静气弹变形和气动载荷通过松耦合迭代方式获得，结构模型和气动模型之间的数据传递是迭代过程的重要组成部分，包括将结构模型节点位移传递至气动模型，以及将气动模型的气动载荷传递至结构模型。目前常用的数据传递方法是由已知的结构模型节点位移构建插值曲面函数，进而获得其它坐标点的位移值，由此确定不同网格位移之间的关系，气动载荷之间的关系则在此基础上由虚功等效原理确定。但插值曲面通过已知位移点，这将使其在部分区域出现波动性，对于机翼结构模型中缺少内部骨架支撑的蒙皮部分更为明显，最终导致弹性机翼气动载荷计算的误差，进而影响静气弹变形和气动载荷设计敏度的合理性与精确性。因此，有必要寻求一种更好的位移传递方法，以解决插值方法进行位移传递时所存在的问题。

弹性机翼重分布气动载荷的设计敏度是静气弹性能设计敏度分析的重要内容。在弹性机翼变形稳定时，对刚度矩阵形式的平衡方程两端求导，然后根据链式法则将气动载荷的设计敏度列阵写为气动载荷关于位移的导数矩阵与位移的设计敏度列阵的乘积形式。这种情况下，可将刚度矩阵与气动载荷关于位移的导数矩阵之差作为位移的设计敏度列阵的系数矩阵，然后求解气动载荷的设计敏度，文献“Neill D J,Herendeen D L.ASTROSEnhancements,Volume 3:ASTROS Theoretical Manual.WL-TR-93-3006,1995”即采用了类似的思路。该算法存在的技术缺陷为：刚度矩阵是一大型对称稀疏的正定矩阵，在计算中对刚度矩阵进行修改将破坏其原有特性和存储结构，更新刚度矩阵的存储结构对大规模结构优化设计是不利的，它将引起计算上的复杂性。本发明针对该问题，提出一种不改变原始刚度矩阵存储结构的模型解析求解方法。以弹性机翼升力效率为例，对静气弹计算与设计敏度分析的过程进行说明。

发明内容

针对上述问题，本发明提出一种大展弦比机翼静气弹性能设计敏度的分析方法。

主要步骤包括：

1)建立结构模型和气动模型；

2)计算结构模型和气动模型间数据传递矩阵；

3)利用步骤2)的计算结果，通过松耦合迭代方式计算弹性机翼的静气弹变形、气动载荷和升力效率；

4)进行弹性机翼气动载荷和升力效率的设计敏度分析。

进一步，所述步骤2)中，结构模型和气动模型间数据传递矩阵的计算，包括如下步骤：

2-1)选择Wendland C2函数为初始紧支径向基函数，然后在初始紧支径向基函数基础上进行一阶一致性修正，利用一次函数与一致紧支径向基函数建立机翼结构模型上、下表面节点垂向位移的最小二乘拟合曲面函数；

2-2)定义各面元气动单元的扭转角和上反角分别为面元控制点处中弧面的弦向和展向切线与弦平面的夹角，并规定抬头扭转时的扭转角为正，向上弯曲时的上反角为正，然后根据步骤2-1)中所得的最小二乘拟合曲面函数，计算面元扭转角列阵和上反角列阵与结构模型上、下表面节点垂向位移之间的数据传递矩阵；

2-3)根据虚功等效原理，使用基于径向基函数插值的方法计算结构模型与气动模型之间载荷的数据传递矩阵。

进一步，所述步骤4)中是通过直接方法或伴随方法进行弹性机翼气动载荷和升力效率的设计敏度分析。

更进一步，所述气动载荷和升力效率的设计敏度分析的直接方法中采用如下迭代格式计算结构位移列阵u_s关于设计变量b的敏度：

若矩阵的谱半径小于1，则迭代过程将对任意初始值收敛，根据矩阵范数的性质和刚度矩阵的对称性得：

式中，ρ表示矩阵的谱半径，||·||₂表示矩阵的谱范数，相对于

的谱范数，刚度矩阵的谱半径是一个更大的数值，所以使用迭代方式计算

是可行的，最终即可由链式法则计算气动载荷的设计敏度，以及升力效率的设计敏度。

更进一步，所述气动载荷和升力效率的设计敏度分析的伴随方法中，根据链式法则可将结构气动载荷列阵关于设计变量的敏度表示为：

取

作为设计变量，并记为M，则有：

根据刚度矩阵的对称性，最终得到如下迭代格式：

该方程组对任意初始值收敛的充要条件是

的谱半径小于1，根据矩阵范数的性质和刚度矩阵的对称性得：

由上式可知，迭代过程将对任意初始值收敛，在求得伴随变量后，由式(18)即可获得f_s关于设计变量的敏度，进而获得升力效率的设计敏度。

在每一迭代步完成后，采用Aitken迭代加速方法对

和M^T进行修改以提高迭代速率。

本发明的有益效果为：

利用一次函数与一致紧支径向基函数建立最小二乘拟合曲面函数，从而实现对机翼弹性变形的模拟，解决了传统插值方法进行位移传递时存在的变形曲面波动性问题。在静气弹计算的基础上，提出基于迭代方式的气动载荷设计敏度的模型解析求解方法，其中包括直接方法和伴随方法，该算法避免了对原始刚度矩阵的破坏，同时具有计算的高效性。

附图说明

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

图2是机翼有限元结构模型，其中，(a)为结构的整体有限元结构模型，(b)为内部骨架的有限元结构模型；

图3是面元模型图；

图4是机翼前缘变形的展向变化曲线；

图5是机翼前缘变形关于第1个设计变量的敏度的展向变化曲线；

图6是机翼前缘变形关于第2个设计变量的敏度的展向变化曲线；

图7是机翼前缘变形关于第3个设计变量的敏度的展向变化曲线；

图8是机翼前缘变形关于第4个设计变量的敏度的展向变化曲线；

图9是机翼前缘变形关于第5个设计变量的敏度的展向变化曲线；

图10是机翼前缘变形关于第6个设计变量的敏度的展向变化曲线。

图中：1是有限元结构模型中的区域1；2是有限元结构模型中的区域2；3是有限元结构模型中的区域3；4是第1个固支点；5是第2个固支点；6是展向坐标与半展长之比；7是垂向位移，单位为mm；8是扭转角，单位为°；9是上反角，单位为°；10是垂向位移沿展向的变化曲线；11是扭转角沿展向的变化曲线；12是上反角沿展向的变化曲线；13是垂向位移的设计敏度，单位为mm；14是扭转角的设计敏度，单位为°；15是上反角的设计敏度，单位为°；16是垂向位移的设计敏度沿展向的变化曲线；17是扭转角的设计敏度沿展向的变化曲线；18是上反角的设计敏度沿展向的变化曲线。

具体实施方式

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

步骤1：结构模型和气动模型的建立。考虑位于笛卡尔直角坐标系Oxyz内的某一大展弦比机翼，该机翼弦平面位于x-y平面内，x轴沿机翼弦向并由前缘指向后缘，y轴沿机翼展向并由翼根指向翼尖，z轴垂直于机翼弦平面向上。采用有限元法建立结构模型，选择结构优化中的设计变量，给定合适的固支点使模型无刚体运动，然后计算结构总体刚度矩阵。采用基于小扰动势流方程的面元法建立气动模型，计算给定飞行状态下的气动力影响系数矩阵，由刚性机翼的几何模型划分CFD网格并计算CFD数据，然后采用斜率法计算气动力影响系数矩阵的对角修正矩阵。

面元法求解气动力的公式为：

f_a＝q·S·(Cp₀+AIC·ΔDW) (1)

式中，f_a为面元气动载荷列阵；q为飞行动压；S为由面元气动单元面积组成的对角矩阵；Cp₀为0°迎角下刚性机翼面元模型的无量纲压力系数列阵；AIC为修正后的气动力影响系数矩阵；ΔDW为给定迎角下弹性机翼下洗相对于0°迎角下刚性机翼下洗的变化量列阵。下洗为各面元控制点在机翼中弧面的垂向落点处，单位速度来流的法向速度分量，而中弧面表示机翼上、下表面的中间位置曲面，因此ΔDW中各元素是结构模型表面节点垂向位移的函数。

步骤2：结构模型和气动模型间数据传递矩阵的计算。计算弹性机翼气动载荷时，需要将结构模型上、下表面节点的垂向位移传递至各面元控制点，得到变形机翼的中弧面并更新下洗。同时，需要将面元气动载荷传递至结构模型，以求解结构模型节点位移。

本发明中，根据大展弦比机翼的变形特点，假设其弦向只有扭转变形而无弯曲变形，于是设机翼上表面或下表面垂向位移的拟合曲面函数为：

式中，m为选定的不同展向位置个数；a_i和b_i为待求解系数；φ(||y-y_i||)为以展向位置y_i为中心的一致紧支径向基函数，它通过对紧支径向基函数进行修正使其满足一致性条件获得，即

使用最小二乘法计算拟合曲面时，式(2)中的系数需满足：

式中，n为结构模型上表面或下表面用于数据传递的节点个数；w(x_j,y_j)和w_fitting(x_j,y_j)分别为(x_j,y_j)处结构模型节点的真实垂向位移值与近似拟合曲面函数的值。

将面元控制点坐标代入式(2)即可得到结构模型上、下表面节点垂向位移传递至面元模型的结果。事实上，各面元控制点处的位移值可表示为结构模型上、下表面节点垂向位移的线性组合的形式。因此，最终各面元控制点处的中弧面垂向位置列阵可表示为：

w_mean＝G·w_s+w₀ (5)

式中，w_mean为面元控制点处的中弧面垂向位置列阵；w_s为结构模型表面节点垂向位移列阵，其中包含上表面和下表面节点的垂向位移；w₀为刚性机翼模型面元控制点处的中弧面垂向位置列阵；G为位移传递矩阵，其行数为面元气动单元个数，列数为结构模型上、下表面节点个数之和。

定义各面元气动单元的扭转角和上反角分别为面元控制点处中弧面的弦向和展向切线与弦平面之间的夹角，并规定抬头扭转时的扭转角为正，向上弯曲时的上反角为正。由扭转角和上反角即可确定中弧面处的法向量，进而计算下洗。因此，在静气弹计算中实际使用的是扭转角列阵和上反角列阵与结构模型表面节点垂向位移之间的数据传递矩阵，它可通过计算中弧面的曲面函数关于x和y的导数，然后将面元控制点坐标代入获得。该位移传递算法不仅可以解决插值方法引起的变形曲面的波动性，还可使得扭转角和上反角的分布与实际力学规律更为相符。

结构模型和气动模型间载荷的关系由径向基函数插值方法和虚功等效原理获得。气动载荷施加在结构模型的下表面，可表示为：

f_FE＝T·f_a (6)

式中，f_FE为结构模型下表面气动载荷列阵；f_a为面元气动载荷列阵；T为载荷传递矩阵，其行数为结构模型下表面节点个数，列数为面元气动单元个数。

步骤3：松耦合迭代计算静气弹变形和弹性机翼升力效率。由前述可知，气动载荷是机翼结构模型表面节点位移的函数，而结构模型节点位移需要根据结构平衡方程求解。因此气动载荷可通过结构模型和气动模型之间的松耦合迭代获得，当前后两次迭代步计算得到的静气弹变形或气动载荷不再变化时迭代计算结束。

升力效率定义为迎角变化时弹性机翼与刚性机翼升力的变化量之比，即

η＝(dL^f/dα)/(dL^r/dα) (7)

式中，L^f为弹性机翼升力；L^r为刚性机翼升力；α为给定飞行状态的迎角。

本发明中采用复步长方法计算升力关于迎角的导数，将迎角及与迎角有关的变量扩展为复数类型，并给迎角虚部很小的扰动量，通过执行松耦合迭代过程，即可很容易得到升力效率。这时，式(7)中的分子和分母可表示为：

dL^f/dα＝Im[L^f(α+ih)]/h (8)

dL^r/dα＝Im[L^r(α+ih)]/h (9)

式中，h为迎角虚部的扰动量；Im表示复变量的虚部。

步骤4：气动载荷设计敏度分析的基本原理。在步骤3完成后，弹性机翼结构满足如下平衡方程：

K·u_s＝f_s (10)

式中，K为结构总体刚度矩阵；u_s为结构位移列阵；f_s为结构气动载荷列阵；下标s表示结构模型中所有无约束节点的自由度。

气动载荷通过位移间接依赖于设计变量，对式(10)两端关于设计变量b求导，并根据链式法则将f_s关于设计变量的导数展开，得：

式中，

和可以显式地计算得到。若解得

即可根据链式法则获得结构气动载荷列阵f_s关于设计变量的敏度：

气动载荷设计敏度的模型解析解可通过直接方法或伴随方法进行计算。在直接方法中，取位移作为状态变量，计算状态变量关于设计变量的敏度，然后再由式(12)计算气动载荷关于设计变量的敏度，详细过程在步骤5中进行说明；伴随方法可避免求解状态变量关于设计变量的敏度，通过求解与设计变量无关的伴随变量，然后由伴随变量计算气动载荷的设计敏度，详细过程在步骤6中进行说明。在实际应用中，需要根据模型的特征和优化目的的不同对使用直接方法或伴随方法进行合理选择。

在使用直接方法或伴随方法进行敏度分析前，需要完成以下两步操作：

1、计算结构总体刚度矩阵由单元刚度矩阵组装而成，单个设计变量的变化仅对与该设计变量相关的单元刚度有影响。在大型结构中，与单个设计变量相关的单元数目相对于总单元数目很少。所以，

可通过下式进行计算：

式中，NDE表示与设计变量b相关的单元数目；

表示与设计变量b相关的第i个单元的刚度矩阵导数；u_i为与设计变量b有关的单元位移列阵。

2、计算

气动载荷是结构模型节点位移的函数，结构气动载荷列阵f_s关于结构位移列阵u_s的导数由二者之间的关系解析获得，它可表示为

式中，B₁和B₂为布尔矩阵，其作用是按照结构模型上、下表面节点自由度在总体刚度矩阵中的位置，对矩阵

中元素的位置进行重新组织及对矩阵维数进行扩充。

由于将复变量引入了升力效率的计算过程，所以

将变为复数类型，其虚部表示由迎角的扰动量而引起的该矩阵的变化。

步骤5：设计敏度计算的直接方法。为避免计算位移关于设计变量的敏度过程中对刚度矩阵的修改，可将式(11)写为如下迭代格式：

采用迭代方式计算位移关于设计变量的敏度，需要考虑迭代计算的收敛性问题。如果矩阵的谱半径小于1，则迭代过程将对任意初始值收敛。根据矩阵范数的性质和刚度矩阵的对称性得：

式中，ρ表示矩阵的谱半径；||·||₂表示矩阵的谱范数。相对于

的谱范数，刚度矩阵的谱半径是一个更大的数值，所以使用迭代方式计算位移关于设计变量的敏度是可行的。为提高迭代计算速率，可采用Aitken迭代加速方法在每一迭代步结束后对进行修改以提高迭代速率。迭代计算收敛时，由式(12)解得气动载荷关于设计变量的敏度。

在计算升力效率关于设计变量的敏度时，式(15)中的

和

都将为复数类型，在迭代过程中需要同时求解位移的实部和虚部关于设计变量的敏度。最终，由式(12)解得气动载荷的虚部关于设计变量的敏度，进而获得弹性升力效率关于设计变量的敏度。

步骤6：设计敏度计算的伴随方法。根据式(11)，结构位移列阵关于设计变量的敏度可以表示为：

将式(17)其代入式(12)，则结构气动载荷列阵关于设计变量的敏度可表示为：

取式(18)中的

作为伴随变量，并记为M，则有：

根据刚度矩阵的对称性，最终得到如下迭代格式：

该方程组对任意初始值收敛的充要条件是

的谱半径小于1。类似于式(16)，根据矩阵范数的性质和刚度矩阵的对称性得：

由式(21)可知，式(20)迭代计算的收敛性条件与式(15)相同，所以采用迭代方式计算伴随变量是可行的，且可采用与直接方法中相同的方式对迭代速率进行加速。在求得伴随变量后，由式(18)即可获得f_s关于设计变量的敏度。

在计算升力效率关于设计变量的敏度时，需要在式(20)的迭代计算过程中同时求解伴随变量的实部和虚部。在迭代收敛后，由式(18)解得结构气动载荷列阵的虚部关于设计变量的敏度，进而获得升力效率关于设计变量的敏度。

设载荷工况为NL个，设计变量为k个，则在直接方法中，每种载荷工况下对每一设计变量，都需要求解位移关于设计变量的敏度，共需对NL×k个列向量进行迭代求解。对于伴随方法，设优化过程中的约束为NC个，对第i个约束其相应的伴随变量的转置中有NADJ_i个非零列向量，这种情况下，需要对ΣNADJ_i个列向量进行迭代求解。所以，若NL×k<ΣNADJ_i，则使用直接方法进行敏度分析优于伴随方法，反之，使用伴随方法优于直接方法。一般在结构初始设计阶段，约束多但设计变量少，这种情况适合采用直接方法；而在结构后期设计阶段，约束少但设计变量多，适合采用伴随方法。

实施例

本发明算法流程的示意图如图1所示，以下将针对具体实施例，对本发明算法正确性和有效性做进一步说明。

本实施例为一大展弦比后掠机翼，该机翼半展长8.90m，翼根弦长2.14m，翼尖弦长0.44m，前缘后掠角为9°，参考面积为8.92m²。以升力效率及其设计敏度为例，检验本发明算法的基本能力。本实施例的具体步骤为：

步骤1，建立大展弦比机翼的结构模型和气动模型。

使用HYPERMESH建立机翼结构的有限元结构模型，主要结构有蒙皮、梁、肋、梁橼条、肋橼条，包含990个壳单元、400个梁单元、746个节点。蒙皮为[±45/0₄/±45/90₂]_s对称层合板，复合材料单层厚度为0.125mm，纵向弹性模量为140GPa，横向弹性模量为8.6GPa，剪切模量为5.0GPa，泊松比为0.33；其它部分材料为钛合金，其弹性模量为110GPa，泊松比为0.3；蒙皮之外的壳单元厚度为3.0mm。图2为机翼的有限元结构模型，图中标示了设计变量所在的区域和结构固支点位置，机翼结构模型根部节点设置为对称边界条件。本实施例中共有6个设计变量，其中第1个和第2个设计变量分别为区域1和区域2中上表面蒙皮层合板的0°铺层厚度；第3个和第4个设计变量分别为区域1和区域2中下表面蒙皮层合板的0°铺层厚度；第5个和第6个设计变量分别为区域3中前梁和后梁腹板的厚度。

图3为该机翼面元模型，共有面元气动单元802个。设置空气压强为标准大气压，马赫数为0.6。使用ICEM对机翼的几何模型进行CFD网格划分，然后使用FLUENT计算该机翼的气动力，流场计算使用S-A湍流模型。根据0°和2°迎角的CFD数据计算气动力影响系数矩阵的对角修正矩阵，同时根据0°迎角的CFD数据计算式(1)中的Cp₀，由于该机翼翼型不对称，所以0°迎角下刚性机翼的升力不为零。

步骤2，计算结构模型与气动模型之间的数据传递矩阵。

选择Wendland C2函数为初始紧支径向基函数，然后在其基础上进行一阶一致性修正，取紧支半径为1000mm，选择机翼结构模型前缘节点的y坐标分量作为式(1)中的不同展向位置。根据结构模型上、下表面节点坐标和面元模型控制点坐标，计算面元扭转角列阵和上反角列阵与结构模型表面节点的垂向位移之间的数据传递矩阵。根据虚功等效原理，使用基于径向基函数插值的方法计算两模型之间的载荷传递矩阵。保存以上矩阵，用于后续静气弹计算与设计敏度分析。

步骤3，计算静气弹变形、气动载荷和升力效率。

为计算升力效率，将迎角、位移、下洗、气动载荷等与迎角有关的变量扩展为复数类型，取飞行迎角为1°，取迎角的虚部扰动量为1.0E-5°，然后通过松耦合迭代方式进行静气弹计算。迭代初始的结构模型位移列阵取为0，由此计算下洗与面元气动载荷，然后将面元气动载荷传递至结构模型的下表面，求解平衡方程并更新位移列阵，进而可由位移列阵计算新的下洗与面元气动载荷，依此迭代，直至得到最终平衡状态下的气动载荷和弹性变形，最后由刚性机翼和弹性机翼升力的虚部计算得到升力效率。迭代收敛时刚性机翼和弹性机翼的升力分别为79688N和53804N，单位迎角变化引起的刚性机翼和弹性机翼的升力变化量为25858N和23002N，弹性机翼的升力效率为88.96％。由于后掠机翼的弯扭耦合效应，机翼最终呈现低头扭转变形，并且升力效率的值小于1。图4给出了机翼前缘处的位移、扭转角和上反角的展向变化曲线，三条曲线均与实际力学规律相符，且从翼根到翼尖无明显的波动性，说明本发明中的位移传递算法是合理有效的。

步骤4，计算刚度矩阵的设计敏度矩阵与结构位移列阵的乘积。

刚度矩阵的设计敏度列阵与位移列阵的乘积是一个列向量，对每一设计变量根据式(13)计算

并保存，用于后续气动载荷设计敏度的计算。这种方法对于减少计算机内存消耗和提高计算效率的作用非常明显。

步骤5，计算结构模型下表面气动载荷列阵关于结构模型表面节点垂向位移列阵的导数。

根据前述气动载荷列阵与位移列阵之间的关系，需要首先计算扭转角列阵和上反角列阵关于结构模型表面节点垂向位移列阵w_s的导数，接着计算ΔDW关于w_s的导数，然后得到f_a关于w_s的导数矩阵，对其前乘载荷传递矩阵进而得到结构模型下表面气动载荷列阵关于结构模型表面节点垂向位移列阵的导数

实际计算中仅对矩阵

进行保存，敏度分析过程中所使用的

可由式(14)所示形式由

获得。

步骤6，使用直接方法或伴随方法进行敏度分析。

为对本发明所提出的算法进行验证，将分别使用直接方法和伴随方法进行气动载荷和升力效率的设计敏度计算。

若采用直接方法进行升力效率的设计敏度计算，对于6个设计变量，需分别迭代计算u_s的设计敏度列阵。取

的初值为0，由此计算并求解式(15)的右端项，接着通过解线性方程组的方式计算新的

然后判断前后两次迭代步计算结果的误差范数是否满足收敛条件，若否，则进行迭代加速处理并执行下一次迭代计算，若是，则迭代计算结束。最后根据链式法则计算气动载荷关于设计变量的敏度，进而获得升力效率关于设计变量的敏度。

图5～图10为直接方法得到的机翼前缘垂向位移、扭转角和上反角关于设计变量的敏度的展向变化曲线。设计变量的摄动将会引起机翼扭转刚度和弯曲刚度的变化，对于第5个和第6个设计变量，其摄动还会引起机翼弹性轴位置的移动，所以在气动载荷的作用下设计变量的摄动将引起机翼弹性变形的变化；另外，机翼弹性变形的变化又会对气动力产生影响，并反过来影响弹性变形；对于后掠机翼，机翼的弯扭耦合效应也将对结果有明显的影响。图中曲线的变化规律是以上因素综合作用的结果。影响图5～图8曲线变化规律的因素以弯扭耦合效应和弯曲刚度的变化为主；图9曲线的变化规律主要受弯曲刚度变化、弹性轴位置移动、弯扭耦合效应影响，由于弹性轴位置的移动和弯扭耦合效应对扭转角的影响相反，最终扭转角的设计敏度呈现图中变化规律；影响图10曲线变化规律的主要因素为弹性轴位置移动和扭转变形引起的气动力的变化。以上结果进一步说明了本发明中的位移传递算法的合理性和有效性。

若采用伴随方法进行升力效率的设计敏度计算，则需要首先迭代计算获得伴随变量。取伴随变量M的初值为0，由此计算

接着求解式(20)的右端项，进而通过解线性方程组获得M^T的各列，然后判断前后两次迭代步计算得到的伴随变量的误差范数是否满足收敛条件，若否，则进行迭代加速处理并执行下一次迭代计算，若是，则迭代计算结束。最终即可通过伴随变量获得气动载荷的设计敏度，进而计算升力效率关于设计变量的敏度。

为了验证本发明中设计敏度数值算法的正确性，使用中心差分法计算升力效率的设计敏度，并与直接方法和伴随方法的计算结果进行比较。具体操作为，在步骤1中将原始结构模型中设计变量取值由b修改为b+Δb或b-Δb，然后执行松耦合迭代过程，分别计算设计变量取值为b+Δb或b-Δb时弹性机翼的升力效率，进而即可由中心差分法计算得到升力效率关于设计变量的敏度。对于本实施例中的6个设计变量，Δb的值统一取为0.01mm。

表1给出了使用直接方法、伴随方法和中心差分法计算得到的升力效率关于各设计变量的敏度。可以看到，直接方法和伴随方法解得的升力效率的设计敏度完全相同，中心差分法与前两种方法所得的结果十分接近。当设计变量取在区域1时，升力效率设计敏度的绝对值明显大于设计变量取在区域2时的设计敏度绝对值，这是因为相对于区域2内的设计变量，区域1内的设计变量的摄动会对机翼变形产生更大的影响。当以相同区域的上、下表面蒙皮层合板0°铺层厚度作为设计变量时，由于模型并非完全对称，设计敏度的计算结果将存在差异。升力效率设计敏度的正负与扭转角设计敏度的正负相一致。以上结果表明了本发明中所提算法的正确性和有效性，该计算结果可进一步指导机翼结构优化设计的进行。

表1大展弦比机翼升力效率的设计敏度计算结果

	直接方法	伴随方法	中心差分法
				第1个设计变量	2.1738E-2	2.1738E-2	2.1911E-2
第2个设计变量	2.5413E-3	2.5413E-3	2.5746E-3
				第3个设计变量	1.4464E-2	1.4464E-2	1.4628E-2
第4个设计变量	1.7781E-3	1.7781E-3	1.8104E-3
				第5个设计变量	-4.4161E-4	-4.4161E-4	-4.3073E-3
第6个设计变量	1.9547E-3	1.9547E-3	1.9546E-3

Claims (6)

1.一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，包括如下步骤：

1)建立结构模型和气动模型；

2)计算结构模型和气动模型间数据传递矩阵；

3)利用步骤2)的计算结果，通过松耦合迭代方式计算弹性机翼的静气弹变形、气动载荷和升力效率；

4)进行弹性机翼气动载荷和升力效率的设计敏度分析。

2.根据权利要求1所述的一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，所述步骤2)中，计算结构模型和气动模型间数据传递矩阵，包括如下步骤：

2-1)选择Wendland C2函数为初始紧支径向基函数，然后在初始紧支径向基函数基础上进行一阶一致性修正，利用一次函数与一致紧支径向基函数建立机翼结构模型上、下表面节点垂向位移的最小二乘拟合曲面函数；

2-2)定义各面元气动单元的扭转角和上反角分别为面元控制点处中弧面的弦向和展向切线与弦平面的夹角，并规定抬头扭转时的扭转角为正，向上弯曲时的上反角为正，然后根据步骤2-1)中所得的最小二乘拟合曲面函数，计算面元扭转角列阵和上反角列阵与结构模型上、下表面节点垂向位移之间的数据传递矩阵；

2-3)根据虚功等效原理，使用基于径向基函数插值的方法计算结构模型与气动模型之间载荷的数据传递矩阵。

3.根据权利要求1所述的一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，所述步骤4)中是通过直接方法或伴随方法进行弹性机翼气动载荷和升力效率的设计敏度分析。

4.根据权利要求3所述的一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，所述气动载荷和升力效率的设计敏度分析的直接方法中采用如下迭代格式计算结构位移列阵u_s关于设计变量b的敏度：

若矩阵

的谱半径小于1，则迭代过程将对任意初始值收敛，根据矩阵范数的性质和刚度矩阵的对称性得：

式中，ρ表示矩阵的谱半径，||·||₂表示矩阵的谱范数，相对于

的谱范数，刚度矩阵的谱半径是一个更大的数值，使用迭代方式计算

是可行的，最终即可由链式法则计算气动载荷的设计敏度，以及升力效率的设计敏度。

5.根据权利要求3所述的一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，所述气动载荷和升力效率的设计敏度分析的伴随方法中，根据链式法则可将结构气动载荷列阵关于设计变量的敏度表示为：

取

作为伴随变量，并记为M，则有：

根据刚度矩阵的对称性，最终得到如下迭代格式：

该方程组对任意初始值收敛的充要条件是

的谱半径小于1，根据矩阵范数的性质和刚度矩阵的对称性得：

由上式可知，迭代过程将对任意初始值收敛，在求得伴随变量后，由式(18)即可获得f_s关于设计变量的敏度，进而获得升力效率的设计敏度。

6.根据权利要求4或5所述的一种大展弦比机翼静气弹性能设计敏度的分析方法，其特征在于，在每一迭代步完成后，采用Aitken迭代加速方法对

和M^T进行修改以提高迭代速率。

Images