技术领域
[0001] 本
发明是一种快速获取复材翼面跨音速颤振速度的方法，属于飞机
复合材料翼面
气动弹性剪裁领域，用于飞机复合材料翼面设计中的跨音速颤振速度优化。
背景技术
[0002] 翼面跨音速颤振是翼面设计的一种重要情况，在跨音速颤振分析过程中，求解得到跨音速颤振临界速度是重要步骤。一般的跨音速颤振计算采用非定常气动
力方法，需要很长的计算时间，给翼面优化带来巨大的困难。当前国内外
飞行器复合材料翼面跨音速颤振优化与分析尚未形成成熟的技术方案。过去大都采用实验和工程估算的方法得到翼面的跨音速颤振临界速度；近年来，随着
计算流体动力学(CFD)、计算结构动力学(CSD)的发展和计算机性能的提高，将CFD和CSD结合起来进行复合材料翼面颤振的数值模拟有了相当大的进展，成为当前气动弹性研究的重要方向之一。和以往的方法相比，采用CFD/CSD耦合进行跨音速颤振求解的特点是
精度较高，但其不足之处在于数值模拟过程的计算量很大，耗时过长。在缩短CFD/CSD耦合计算的时间方面，除了依赖计算机
硬件性能的提高、利用大型机参与计算外，还采用了计算机集群并行计算的方法。
[0003] 为设计出满足跨音速颤振要求的飞行器，需要对飞行器复合材料翼面结构进行跨音速颤振优化。优化中若直接采用求解非定常气动力的方法计算跨音速颤振临界速度，计算时间是极其惊人的。造成优化
迭代的单步时间长达数月，这在工程上无法接受。
[0004] 利用响应面进行翼面优化，必须解决复合材料板参数和设计变量可行域描述。这样可以找到响应面的边界，实现高精度插值拟合。而常规复合材料板的描述方法无法满足上述要求。利用三向
刚度方法描述板的力学特性，可以简化复合材料板的多种设计变量参数(包括各种
角度和相应的铺层比例，铺层厚度等)，便于进行翼面优化设计。
[0005] 三向刚度特性是复合材料板元件固有的力学特性，对于一
块给定厚度(层数)的复合材料平板，0°、±45°、90°铺层的比例直接影响到平板的轴向刚度特性EI1、EI2和剪切刚度特性GJ，如
附图1所示。铺层角度、比例各异的复合材料板均可转化为相同厚度的、仅有0°、±45°、90°铺层的复合材料板，且保证板的三向刚度特性相同。由于复合材料板的三向刚度特性参数EI1、EI2、GJ之间存在约束关系，可用其中的两个参数表示另外一个参数，由此得到复合材料板元件刚度可行域的数学表述。具体的公式推导如下：
[0006]
[0007] 其中 δ为板厚，U1、U4、U5为与材料特性相关的不变数，U1、U4、U5的计算方法如下：
[0009] Q66＝G12
[0010] Q12＝μ12Q22＝υ21Q11，
[0011] 其中E11、E22、μ12、μ21、G12为与材料相关的工程弹性常数则[0012]
[0013]
[0014]
[0015] 将C值带入(I)，
整理得到EI1、EI2和GJ的约束关系，如(II)式所示。
[0016] (EI1-EI2)2GJ2+(C4-C6-C2)GJ+(C5+C1-C3)＝0 (II)
[0017] 其中参数C1～C6的表达式如下，
[0018] C1＝2EI1EI2(EI1+EI2)(U1+U5)δ
[0019] C2＝2EI1EI2(EI1+EI2)
[0020] C3＝4EI1EI2(U1+U5)2δ2
[0021] C4＝8EI1EI2(U1+U5)δ
[0022] C5＝(EI1+EI2)2(U5-U4)2δ2
[0023] C6＝2(EI1+EI2)2(U5-U4)δ
[0024] 复合材料平板三向刚度等效转换与刚度特性参数间约束关系式的推导为本方法中建立刚度可行域提供了理论
基础，保证了刚度可行域的正确性与有效性。
发明内容
[0025] 本发明正是针对上述
现有技术存在的问题而设计提供了一种快速获取复材翼面跨音速颤振速度的方法，其目的是解决复合材料翼面跨音速颤振求解、分析与优化中复合材料板参数与设计变量可行域的描述问题，以及计算时间过长的问题。
[0026] 鉴于跨音速颤振优化能力的限制，本方法在亚音速颤振优化的基础上，
选定有限个翼面敏感区称为敏感元件，利用三向刚度来描述敏感元件的力学特性，进而建立敏感元件的三向刚度可行域，在刚度可行域内选择特征点进行非定常气动力颤振计算，利用计算结果分别建立每个敏感元件与整个复合材料翼面的跨音速颤振临界速度响应面。在其后跨音速颤振求解、分析与优化的过程中，利用查询响应面的方式快速获取翼面的跨音速颤振速度，从而代替了长时间的非定常气动力计算，达到减少计算量，缩短计算时间的目的。
[0027] 本发明的目的是通过以下技术方案来实现的：
[0028] 该种快速获取复材翼面跨音速颤振速度的方法，其特征在于：该方法的步骤是：
[0029] (1)对飞机翼面蒙皮进行亚音速颤振优化，优化后上、下翼面蒙皮的有限元模型上，将铺层厚度δ相差5层以内的几个相邻设计单元组成一个敏感元件，单侧上、下翼面蒙皮上共选取2～10个敏感元件，左、右翼面蒙皮上选取的敏感元件对称，敏感元件包含的区域不重叠；
[0030] (2)对于每个敏感元件，计算其所包含的所有设计单元的平均厚度，并将该敏感元件所包含的所有设计单元规整化。规整化的具体步骤为：找出厚度与平均厚度最接近的设计单元，将同一敏感元件中其余设计单元的铺层去除，替换为该设计单元的铺层，并将该设计单元的厚度作为该敏感元件规整后的平均厚度d；
[0031] (3)对每一个敏感元件，建立刚度可行域并选取特征点
[0032] 影响复合材料翼面跨音速颤振特性的变量是翼面蒙皮参数，通常包括铺层厚度(层数)、铺层角度和各种铺层角之间的比例关系。在翼面跨音速颤振分析与优化的过程中，各种铺层之间比例的取值是不连续的，且铺层厚度(层数)、角度、比例的组合情况很多，参数变量的不同组合与跨音速颤振速度Vcp之间很难找到特定的解析关系。鉴于上述原因，难以直接找到参数变量及变量组合的可行域，并建立参数变量与跨音速颤振速度Vcp之间较为精确的响应面(全部为内插值的响应面)。
[0033] 本方法采用三向刚度力学参数代替复杂的铺层组合，建立敏感元件的三向刚度可行域，以可行域中的特征点与对应的颤振速度为基础建立该敏感元件与翼面跨音速颤振速度之间的响应面。根据各敏感元件自身特有的约束条件，建立不同的三向刚度可行域，进一步建立对应于每个敏感元件的跨音速颤振速度响应面。
[0034] 对每一个敏感元件，建立刚度可行域并选取特征点的过程如下：
[0035] 3.1确定在建立刚度可行域时，每层铺层可选用的铺层角θ
[0036] 在优化过程中，每层的铺层角可在0°～180°范围内任意选取。由于其他所有的铺层组合均可在厚度不变的情况下转换为0°、±45°和90°的铺层组合，且转换前后的刚度特性不变。因此在建立刚度可行域时，每层铺层可选用的铺层角θ为0°、±45°和90°之间的一种，以此三种铺层角的组合代表该敏感元件所有可能的铺层状态，这同时也保证了工程上实施的可行性。
[0037] 3.2确定每一个敏感元件的最小许用厚度dmin和最大许用厚度dmax[0038] 本方法选用规整后的设计单元的平均厚度di为最小许用厚度dmin，根据敏感元件许用的结构空间和许用的结构应变
水平，确定该敏感元件的最大许用厚度dmax。在建立刚度可行域时，利用最小许用厚度dmin和最大许用厚度dmax，可以保证可行域是有边界的，进而保证建立的跨音速颤振速度响应面是有边界的。
[0039] 3.3根据每层铺层可选用的铺层角θ、敏感元件的最大许用厚度dmax和最小许用厚度dmin，确定每个敏感元件所有可能的铺层组合
[0040] 在实际应用中，由于工程设计上的要求，铺层板应为复合材料对称均衡板。即铺层角θ为±45°时，对应的铺层为4层的整数倍，铺层角θ为0°或90°时，对应的铺层为偶数层。在这种约束条件下，敏感元件的刚度可行域并非一个连续的曲面，而是由许多离散的点组成，且这些点所对应的三向刚度特性满足(II)式的约束关系，即这些离散点都在(II)式表示的隐式曲面上。
[0041] 3.4分别计算每个敏感元件在每种铺层角与铺层厚度组合下的三向刚度特性值EI1、EI2、GJ，计算步骤如下：
[0042] 首先计算纤维主向刚度系数Qij，
[0043] Q66＝G12
[0044] Q12＝μ12Q22＝υ21Q12，
[0045] 其中E11、E22、μ12、μ21、G12为与材料相关的工程弹性常数然后计算层合板拉伸刚度矩阵，
[0046]
[0047] 其中N为层合板的层数，tk为
单层厚度
[0048] 对于每一层复合材料铺层，铺层角为θ时的刚度系数 计算公式，[0049]
[0050]
[0051]
[0052]
[0053]
[0054]
[0055] 其中m＝cosθ，n＝sinθ，θ为当前层的铺层角
[0056] EI1、EI2、GJ计算式如下
[0057] EI1＝A11-A12*A12/A22；
[0058] EI2＝A22-A12*Al2/A11；
[0059] GJ＝A66；
[0060] 3.5绘制刚度点，分别创建每一个敏感元件的刚度可行域
[0061] 步骤3.4中，每一组EI1、EI2、GJ数据作为敏感元件的一个刚度点，将敏感元件的所有刚度点绘制到EI1、EI2、GJ构成的三维空间中，这些离散的刚度点即构成该敏感元件的刚度可行域，每个敏感元件的可行域所包络的空间均呈现类似三棱台的形状，在同一铺层厚度下，各刚度点在一帆状曲面上，从不同角度观察此帆状曲面如附图2、附图3、附图4；不同厚度的铺层对应的帆状曲面形状相似，彼此为比例缩放关系，缩放比为厚度比的四次方，如附图5。
[0062] 3.6对每一个敏感元件，分别在其刚度可行域内选取特征点
[0063] 在类似三棱台的刚度可行域的顶面、底面和内部等分剖面上选取特征点，选取
位置为上述每面的三个
顶点和三条边的中点，其中内部等分剖面的数量根据刚度可行域的大小选取，一般选取1～4个，特征点数目介于18～36个之间。如计算量可以承受，应考虑在上述各面上适量增选特征点，以提高跨音速颤振速度响应面的插值拟合精度。
[0064] (4)对每个敏感元件，分别建立该敏感元件的跨音速颤振速度响应面[0065] 一般情况下，翼面上会有多个敏感元件，对应每个敏感元件，有不同的刚度可行域和相应的特征点，因此需要分别针对每个敏感元件建立跨音速颤振速度响应面，其建立过程分为如下步骤：
[0066] 4.1创建刚度特征点对应的翼面有限元计算模型文件
[0067] 对于每个敏感元件，去除敏感元件上的原有铺层，用敏感元件的每一个特征点所代表的铺层铺覆在敏感元件上，翼面其他部位的铺层不变，由此建立敏感元件的每一个特征点对应的翼面有限元计算模型文件，每个敏感元件创建与特征点一一对应的18～36个有限元计算模型文件。如增选特征点，则对应增加相应的有限元计算模型文件；
[0068] 4.2分别计算4.1中每个翼面有限元计算模型的非定常气动力，得到相应的跨音速颤振临界速度Vcp，如增选特征点，则Vcp个数会对应增加；
[0069] 4.3对每个敏感元件，分别建立跨音速颤振速度响应面
[0070] 将敏感元件的每个特征点的轴向刚度特性EI1、剪切刚度特性GJ和厚度δ作为自变量(18～36组或更多)，对应的翼面跨音速颤振速度Vcp(18～36个或更多)作为目标值，进行高维曲面拟合与插值，得到敏感元件的跨音速颤振速度响应面。由于另一方向的轴向刚度EI2与EI1和GJ存在(II)式所示的固定的约束关系，则EI2不需要作为自变量参与跨音速颤振速度响应面的建立。
[0071] 在响应面建立后，在敏感元件的刚度可行域内的任意一点，通过查询对应的响应面来代替复杂的非定常气动力计算，即可快速获得该状态下翼面的跨音速颤振速度Vcp。
[0072] 本发明技术方案的优点和积极效果是：
[0073] 1.利用三向刚度来描述敏感元件的铺层和
质量特性，解决了如何建立
覆盖复合材料板各种铺层角度与铺层比例可行域的技术难点；
[0074] 2.通过采用有限个特征点建立响应面的方法，保证了可行域里所有点进行跨音速颤振速度插值时，都落在所建立的响应面边界内，保证了较高的插值精度；
[0075] 3.利用查询响应面的方式代替复杂的非定常气动力计算用以获取翼面的跨音速颤振速度，大幅提高了计算速度，缩短了计算时间，解决了复合材料翼面跨音速颤振优化在工程上应用的难点；
[0076] 4.经过对某飞机复合材料翼面跨音速颤振优化计算的验证，建立的基于三向刚度的跨音速颤振速度响应面误差范围为±1m/s，响应面比较光滑，表明了其鲁棒性好的特点，在颤振分析与优化、以及实际工程上均具实用性。
附图说明
[0077] 图1为复合材料板三向刚度特性示意图
[0078] 图2为某固定厚度下的刚度点构成的帆状曲面
[0079] 图3为图2的另一视角
[0080] 图4为图2的另一视角
[0081] 图5为不同厚度下的刚度点构成的帆状曲面组
[0082] 图6为某个敏感元件对应的刚度可行域与其中的特征点
[0083] 图7为本发明
实施例中上下翼面敏感元件的分布图
[0084] 图8为最大许用厚度面上特征点的选取示意图
[0085] 图9为响应面的误差分析曲线
[0086] 图10为某固定厚度下响应面的剖面图
具体实施方式
[0087] 以下将结合附图和实施例对本发明技术方案作进一步地详述：
[0088] 以下将以某型飞机复合材料
尾翼为例，介绍采用本发明所述快速获取复材翼面跨音速颤振速度的方法的实施过程，该过程的步骤是：
[0089] 1.首先对飞机翼面蒙皮进行亚音速颤振优化，优化方法和过程采用在先
申请的国防
专利技术方案，该专利申请的名称是：飞机翼面蒙皮亚音速颤振优化方法，专利申请号是：201110171379.9，该技术方案的主要内容是：用翼面每块蒙皮元素作为优化对象，同时将蒙皮的铺层角度，铺层厚度增量、铺层位置选作设计变量，保证遗传
算法的搜索是在包括了所有可能解的空间域内进行，采用双循环优化流程，解决了大设计变量(＞1000)的解耦问题，同时采用“小增量填谷法”有效的避免了算法的震荡，确保算法的高重复性。
[0090] 在经过亚音速颤振优化后的水平尾翼上、下翼面蒙皮的有限元模型上，将铺层厚度δ相差5层以内的几个相邻设计单元1组成一个敏感元件2，对于本例来说，共选取8个敏感元件2，敏感元件2位于上下翼面前缘中央略偏翼尖处，以及
后缘升降
舵作动器附近，其分布如附图7所示。平尾左、右翼面蒙皮上选取的敏感元件2对称，敏感元件2之间的区域不重叠。
[0091] 方法中，需针对每一个敏感元件2分别建立跨音速颤振速度的响应面。以下内容以其中一块敏感元件2为例，说明建立刚度可行域，选取特征点以及建立相应的跨音速颤振速度响应面的过程。
[0092] 2.分别计算该敏感元件2所包含的所有设计单元1的平均厚度，找出厚度与平均厚度最接近的设计单元1，将其余设计单元1的铺层去除并替换为该设计单元1的铺层，同时将该设计单元1的厚度作为该敏感元件2规整后的平均厚度d；
[0093] 3.建立该敏感元件2的刚度可行域并选取特征点，其过程如下：
[0094] 3.1按照工程常用的情况，确定在建立刚度可行域3时，敏感元件每层铺层可选用的铺层角θ为0°、±45°和90°之间的一种；
[0095] 3.2确定该敏感元件的最小许用厚度dmin和最大许用厚度dmax
[0096] 3.2.1确定最小许用厚度dmin
[0097] 该敏感元件2的最小许用厚度dmin为8层(单层0.12mm)；
[0098] 3.2.2确定最大许用厚度dmax
[0099] 该敏感元件2的最大许用厚度dmax为56层(单层0.12mm)；
[0100] 3.3根据铺层角θ、敏感元件的最小许用厚度dmin和最大许用厚度dmax，确定所有可能的铺层组合；
[0101] 3.4分别计算敏感元件2在每种铺层组合下的三向刚度EI1、EI2、GJ，计算步骤如下：
[0102] 首先计算纤维主向刚度系数Qij，
[0103] Q66＝G12
[0104] Q12＝μ12Q22＝υ21Q11，
[0105] 其中E11，E22、μ12、μ21、G12为与材料相关的工程弹性常数然后计算层合板拉伸刚度矩阵，
[0106]
[0107] 其中N为层合板的层数，tk为单层厚度
[0108] 对于每一层复合材料铺层，铺层角为θ时的系数 的计算公式，
[0109]
[0110]
[0111]
[0112]
[0113]
[0114]
[0115] 其中m＝cosθ，n＝sinθ，θ为当前层的铺层角EI1、EI2、GJ计算式如下[0116] EI1＝A11-A12*A12/A22；
[0117] EI2＝A22-A12*A12/A11；
[0118] GJ＝A66；
[0119] 对于本实例来说，该尾翼选用的复合材料的工程弹性常数为：E11＝1.29E+11Pa，E22＝9.82E+9Pa，μ12＝0.3056，G12＝5.29E+9Pa，铺层的单层厚度tk＝0.12mm，敏感元件2的最小许用dmin厚度为8层，最大许用厚度dmax为56层。
[0120] 3.5步骤3.4中，每一组EI1、EI2、GJ数据做为敏感元件2的一个刚度点，将敏感元件2的所有刚度点绘制到EI1、EI2、GJ构成的三维空间中，即构成敏感元件2的刚度可行域3，刚度可行域3所包络的空间呈现类似三棱台的形状，参见附图6。
[0121] 3.6在刚度可行域3内选取特征点
[0122] 在类似三棱台的刚度可行域3的顶面4(最小许用厚度面)、底面5(最大许用厚度面)和内部等分剖面6上选取特征点7，本实例中选择两个三等分剖面6，厚度分别为24层和40层，特征点7的选取位置为每个面的三个顶点和三条边的中点，参见附图6所示。
[0123] 在类似三棱台的刚度可行域3的顶面4、底面5和内部等分剖面6上，根据计算能力，可再增选特征点7。如本例中在底面5(最大许用厚度面)的每条边上分别增选两个边等分点71，在对称轴线上增选三个轴等分点72，以及在内部1/4高处，左右对称增选两个内部点73，具体的增选情况参见附图8所示。特征点7总计达到35个。
[0124] 4.建立跨音速颤振速度响应面
[0125] 4.1创建刚度特征点7对应的翼面有限元计算模型文件
[0126] 去除敏感元件2上规整后的铺层，用敏感元件2的每一个刚度特征点7所代表的铺层逐一铺覆在敏感元件2上，建立敏感元件2的每一个刚度特征点7对应的翼面有限元计算模型文件；
[0127] 4.2分别计算4.1中每一个翼面有限元计算模型的非定常气动力，得到翼面跨音速颤振速度Vcp，特征点的轴向刚度EI1、剪切GJ，厚度δ及相应的翼面跨音速颤振速度Vcp如下表：
[0128] 表1.特征点的刚度与跨音速颤振速度
[0129]
[0130] EI1：轴向刚度(MPa)，GJ：剪切刚度(Mpa)，δ：厚度(mm)，Vcp：跨音速颤振速度(m/s)
[0131] 4.3建立敏感元件2跨音速颤振速度的响应面8
[0132] 将表1中每一个刚度特征点7的轴向刚度EI1、剪切刚度GJ和厚度δ作为自变量，对应的翼面跨音速颤振速度Vcp作为目标值，输入至响应面
软件DesParO中进行高维曲面拟合与插值，得到敏感元件2的跨音速颤振速度响应面8。DesParO软件是由德国科学院数学与科学计算研究所开发的，由德国scapos AG公司代销的软件，该软件是用于建立优化设计响应面的
应用软件，是一种用于多目标参数优化设计的工具，软件可以在高维空间建立响应面，并有很高的拟合精度。软件还能够进行响应面的误差分析，保证计算结果具有3西格玛级
置信度。
[0133] 在响应面8建立后，在该敏感元件2的刚度可行域3内的任意一点，通过查询对应的响应面8来代替复杂耗时的非定常气动力计算，即可快速获取该状态下的翼面跨音速颤振速度Vcp。
[0134] 4.4敏感元件2跨音速颤振速度响应面8的误差分析
[0135] 以表1中35个特征点的数据为基础，每次选用34个点建立跨音速颤振速度响应面，另一点作为验证点。采用对比分析的方法，以验证点的刚度EI1、GJ和厚度δ作为输入的自变量，查询响应面得到跨音速颤振速度Vcp1，Vcp1与该点经过非定常气动力求解得到的真实的跨音速颤振速度Vcp相比，即可得到该点的插值误差。此外，响应面软件DesParO有插值误差估算功能，也可大致估算出该点的插值误差范围。分析中获得的数据如下表2所示，误差曲线见附图9。
[0136] 表2.跨音速颤振速度响应面误差分析
[0137]
[0138]
[0139] Vcp：实际跨音速颤振速度，Vcp1：插值跨音速颤振速度，参考误差均为±值Vcp、Vcp1、参考误差、实际误差的单位均为m/s
[0140] 在厚度δ确定时，某个敏感元件1对应的跨音速颤振速度响应面剖面形状如附图10所示。由跨音速颤振速度响应面建立的过程及相关的数据与图表，可以得出如下结论：
[0141] 1.软件估算的插值误差为±0.5m/s，实际验证的误差为±1m/s，实际误差大部分在估算的插值误差范围内，响应面插值拟合的精度较高；
[0142] 2.所建立的基于三向刚度可行域的跨音速颤振速度响应面表面光滑，没有突出的尖点，响应面的鲁棒性好；
[0143] 3.该跨音速颤振速度响应面的建立花费数天时间，计算量在工程上可以接受；
[0144] 4.响应面建立后，在其上查询一次跨音速颤振速度的时间介于0.1s～0.3s之间，远低于完成一次非定常气动力计算所需时间，用查询响应面的方式来代替复杂耗时的非定常气动力计算，可以实现复材翼面跨音速颤振速度的快速获取。
