CN110866353B - 基于应变邻域的飞机复合材料结构优化方法 - Google Patents

Abstract

本发明涉及飞机复合材料结构优化设计领域，具体涉及基于应变邻域的飞机复合材料结构优化方法，本发明方法在划分“应变邻域”之前利用最小二乘法对应变进行曲面拟合，用拟合后的应变值来作为划分邻域的依据，在结构优化之前，对结构蒙皮传力过程进行应变敏度分析，将应变相近的单元分成一个“应变邻域”，在结构优化时以“应变邻域”为对象进行优化，并在结构优化完成之后对优化后的蒙皮厚度进行二次曲面拟合；本发明解决了超大变量的飞机复合材料结构精细化设计中优化结果不满足传力连续性和不满足工艺制造要求的问题；本方法使优化后铺层结果体现力的传递连续性，解决了蒙皮铺层的厚度连续性和蒙皮表面光顺度问题，利于工艺制造。

Description

基于应变邻域的飞机复合材料结构优化方法

技术领域

本发明属于飞机复合材料结构优化设计领域，具体涉及基于应变邻域的飞机复合材料结构优化方法。

背景技术

随着航空技术的发展，结构精细化设计成为结构减重的重要手段。结构优化使结构设计可以实现结构精细化设计，成熟可靠的优化技术与工具可以有效的提高结构利用效率，实现结构重量的大幅降低。随着先进复合材料在航空飞机结构中的使用比例越来越高，复合材料飞机结构的优化设计成为了人们研究的重点方向。当前世界大型民用飞机被美国波音公司和欧洲空中客车公司垄断，双方对飞机的结构设计与优化开发了自己的软件系统。设计变量数成为优化技术发展的重要指标，空客公司研发的LAGRANGE系统已经能够处理一万五千设计变量的多学科优化设计，成为该系统领先的标志。

中国航空研究院开发的SAMOS优化系统实现了十万以上设计变量的优化。超大设计变量具有独立性和离散性，使得目前超大设计变量的优化结果无法满足工艺制造的要求，也不能做到传力的连续性。然而现有的基于敏度分析的遗传算法SAMOS系统，在复合材料结构超大设计变量优化中，虽然能够在满足设计要求的前提下，有效地减轻结构的重量，但是优化后很难得到一个理想的铺层图，如图1所示。主要原因有如下三点：

1、由于超大设计变量(10万设计变量)的相互独立性，设计变量单元离散地分布在结构的各个部位，不同敏度的细小单元交叉相邻，因而，单元铺层的增加也会随之呈现出离散化；

2、遗传算法本身具有较大的随机性，随机投点造成每一步单元厚度变化不相同，形成了相邻单元高度突变，导致优化后的机翼结构厚度凹凸不平；

3、发动机等机翼上的集中载荷造成结构力的突变。

目前基于遗传算法的十万设计变量SAMOS系统结构优化后，出现复合材料蒙皮铺层不光顺的问题，给制造工艺带来极大的困难，增加生产成本，而且从传力角度来讲，结构有明显的厚度突变传力也不合理。工程上常利用数学方法如曲面拟合等进行光顺后处理遇到以下问题：

1、对于单极值的厚度变化可以得到较好的结果；

2、采用数学方法拟合得到的蒙皮厚度改变了最初的优化结果，此时的结构需要进一步进行优化。

发明内容

本发明的目的是：设计基于应变邻域的飞机复合材料结构优化方法，以解决在超大变量(10万以上变量)的飞机复合材料结构精细化设计中，由于变量的离散性和独立性而使优化结果不满足传力连续性和不满足工艺制造要求的技术问题。

为解决此技术问题，本发明的技术方案是：。

提供一种基于应变邻域的飞机复合材料结构优化方法，所述的飞机复合材料结构优化方法在结构优化之前，对结构蒙皮传力过程进行应变敏度分析，将应变相近且空间位置相邻的单元分成一个“应变邻域”，在结构优化时以“应变邻域”为对象进行优化；所述空间位置相邻的单元是指根据四联通原理有公共边的单元。

在划分“应变邻域”之前利用最小二乘法对计算应变值进行曲面拟合，用拟合后的计算应变值来作为划分邻域的依据，并在结构优化完成之后对优化后的蒙皮厚度进行二次曲面拟合。

所述的飞机复合材料结构优化方法包含以下步骤：

步骤一、建立蒙皮的有限元计算模型；模型变量包括：蒙皮单元个数、蒙皮单元厚度、蒙皮单元铺层角度、蒙皮初始厚度、复合材料单层厚度、蒙皮单元铺层层数；

复合材料单层厚度就是每次外循环增加一层铺层的厚度，蒙皮单元铺层层数是指这个单元的总厚度，是多个单层厚度的厚度叠加；

步骤二、针对模型进行蒙皮受力的有限元分析；

步骤三、进行应变邻域划分：读取每个有限元单元对应的应变大小，利用最小二乘法对每个单元的计算应变值进行曲面拟合，

根据各个邻域所包含有限元单元的计算应变值的最大值，按照计算应变值从大到小对各应变邻域排序；

步骤四、利用遗传算法以“应变邻域”为对象进行结构优化：

以邻域为初始种群的遗传算法来计算最佳铺层块数、铺层角度和蒙皮单元铺层层数，形成初始计算n个个体的种群；种群中的每一个个体代表着某一种添加重量的邻域块数和相应添加的铺层角度的组合；计算种群中个体的计算应变值大小，进行遗传算法操作，基于最优保存策略保存最优个体，确定结构重量最佳的添加块数和相应的铺层角度；

步骤五、对每个应变邻域进行整体增重迭代：每块应变邻域内包含的每个单元增重相同，每块应变邻域内包含的每个单元的铺层角度也相同；

步骤六、当所有单元的拉/压/剪切应变都满足许用应变时进入下一步；不满足，返回步骤一，更改模型变量，继续迭代；

步骤七、将上一步的优化后的蒙皮厚度结果进行二次曲面拟合；输出拟合后的厚度云图。

步骤一中：蒙皮初始厚度为0.8mm，铺层角度0°∶90°∶±45°＝1∶1∶1。

步骤三中有限元单元对应的计算应变值ε_i的计算公式为：

其中，[ε_X]为复合材料纵向/横向拉伸许用应变，[ε_Y]为复合材料纵向/横向压缩许用应变，[ε_XY]为复合材料剪切许用应变，ε_x、ε_y、ε_xy分别是纵向/横向拉伸应变、纵向/横向压缩应变、剪切应变。

步骤三有限元分析中复合材料单层的强度失效准则采用最大应变失效准则；对不满足失效准则的蒙皮单元进行优化，以蒙皮单元形状形心的位置作为蒙皮单元的样本点，对n个样本点(X_i，,Y_i,Z_i)，(i＝0,1,2...n)，根据如下曲面方程进行拟合：

f(x,y)＝P^T(x,y)A，其中A^T＝(a₀,a₁,…a_m)是曲面方程的系数，X_i，Y_i是单元的坐标，Z_i是单元的计算应变值。

优选地，在步骤三中还包含对划分后的邻域进行优化的操作，具体为：

在划分出所有邻域后，设定一个面积阈值，若某一邻域的面积小于这个阈值，则再判定该邻域的周围哪块邻域的计算应变值最大，将该邻域归并到周围计算应变值最大的邻域中。

步骤四中具体遗传算法步骤如下：

每块应变邻域的增重通过如下方式确定：每一步铺层选取总重量的5％-10％迭代，按照邻域的面积大小和应变范围确定每个邻域铺层重量分配和增加的铺层厚度：

若遗传算法每次输出的结果是给n块邻域铺层，那么应变最大的n块邻域被选中，染色体解码的角度就是最佳铺层角；

选中的第i个邻域的平均计算应变公式如下：

N是第i个邻域内包含的单元个数；

邻域面积S_i(i＝1,2…,n)为邻域内包含的各个单元的面积之和；

被选中的第i个邻域总计算应变值ε_sum-i＝ε_avg-i*S_i，

第i个邻域铺层重量分配：

第i个邻域增加的铺层厚度：

其中，Δw为一次给所选中的所有应变邻域的增重的总和，Δw_i为一次给选中的第i个应变邻域的增重，i是领域编号，i＝1,2…n，j为第i个邻域内包含的第j个单元编号，ρ是所选用的复合材料的密度。

步骤四中第i个邻域铺层角度确定方法具体如下：

根据如下公式

N是第i个邻域内包含的单元个数，i＝1,2…,n；α_j是领域内第j个单元的主应变角；

计算出每个邻域的平均主应变角，在主应变角的±128°范围内搜索最佳铺层角，并按照最佳铺层角进行铺层；最后用刚度等效原理把角度转换成工程用的铺层角；所述工程用的铺层角为0°∶90°∶±45°。

角度转换成工程用的铺层角具体操作为：

按照实际的主应变角作为铺层角度，优化完成后再转换成0°∶90°∶±45°，等效后层合板的总厚度以及刚度特性保持不变。

优选地，步骤四中种群初始化采用随机投点方法。

本发明的有益效果是：本发明的基于应变邻域的飞机复合材料结构优化方法加入基于应变邻域的优化步骤能解决上述问题，能使优化后铺层结果体现力的传递连续性。也能解决蒙皮铺层的厚度连续性和蒙皮表面光顺度问题，利于工艺制造。

本发明方法在一定载荷下(本方法是采用极限载荷，2.5g，安全系数1.5)按照计算应变值的大小划分邻域，为使应变划分结果有较少的噪点且邻域边界光顺，对实际计算应变值采用最小二乘法进行高阶曲面拟合，用拟合后的应变值来作为划分邻域的依据。

把各个邻域的铺层角度作为设计变量考虑，而不是采用工程常用的铺层角度：0°∶90°∶±45°。计算出每个邻域的平均主应变角，在主应变角的±128°范围内搜索最佳铺层角，并按照最佳铺层角进行铺层。最后用刚度等效原理再把角度转换成工程常用的铺层角。

由于超大设计变量的独立性和离散性，采用超大设计变量需要去噪点和边界光顺，除了每一步邻域划分中需要采用到最小二乘法进行曲面拟合外，对多步优化后的单元厚度也要进行最小二乘法进行曲面拟合。经试验此方法可以有效去除噪点和光顺边界，优化处理后的结果显示的传力路径更合理。

本方法采用利用遗传算法的全局搜索能力，并且拥有较大的设计变量空间(变量可达10万以上)，在所有可行域空间进行搜索，能够更准确的找到影响应变大小的关键因素，在多个型号应用中得到了成功的验证。

附图说明

为了更清楚地说明本发明实施的技术方案，下面将对本发明的实例中需要使用的附图作简单的解释。显而易见，下面所描述的附图仅仅是本发明的一些实施例，对于本领域的技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为目前基于敏度分析的遗传算法SAMOS系统优化后的铺层示意图；

图2为四连通原理示意图；

图3为基于应变邻域的示意图；图中，C1、C2、C3、C4分别代表机翼上的四块邻域；

图4为初始优化模型单元划分示意图；

图5为本发明的方法流程图；

图6为本发明的遗传算法的流程图；

图7为本发明优化完成后对单元厚度的拟合流程图；

图8为按照等分计算应变值范围的计算应变值范围-单元数目示意图；

图9为去除计算应变值范围最大的前5％单元后等分计算应变值范围的计算应变值范围-单元数目示意图；

图10为最小二乘法对计算应变值大小进行曲面拟合后的拟合图；

图11为染色体编码范围示意图；

图12为最小二乘法对单元厚度大小进行曲面拟合后的拟合图；

图13为应变邻域划分结果图；

图14为优化一定步骤后的厚度云图显示。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域的普通技术人员在没有做出创造性劳动的前提下，所获得的所有其他实施例，都属于本发明保护的范围。

下面将详细描述本发明实施例的各个方面的特征。在下面的详细描述中，提出了许多具体的细节，以便对本发明的全面理解。但是，对于本领域的普通技术人员来说，很明显的是，本发明也可以在不需要这些具体细节的情况下就可以实施。下面对实施例的描述仅仅是为了通过示出本发明的示例对本发明更好的理解。本发明不限于下面所提供的任何具体设置和方法，而是覆盖了不脱离本发明精神的前提下所覆盖的所有的产品结构、方法的任何改进、替换等。

在各个附图和下面的描述中，没有示出公知的结构和技术，以避免对本发明造成不必要的模糊。

下面以某客机飞机机翼复合材料结构优化实例为例说明本发明的基于应变邻域的飞机复合材料结构优化方法：

·优化初始条件

如图4为客机飞机机翼初始优化模型，模型共有33476个节点，47128个蒙皮单元，模型大小为15M。选择客机机翼上下翼面蒙皮作为优化对象。蒙皮初始厚度为0.8mm，复合材料单层厚度0.1mm，铺层角度0°∶90°∶±45°＝1∶1∶1，上下翼面均为对称。初始机翼中央翼盒蒙皮壁板重量质量66.7kg。

·优化变量选择

选择单侧机翼上下翼面共25004个蒙皮单元，每次划分邻域上下翼面总共不超过35块。优化铺层角度范围-90°至90°。

·优化参数设定

外循环最大代数设定200步，内循环代数为10代，种群大小为100，交叉概率为0.7，变异概率为0.1。每次添加初始重量的1％。

·优化结果

如图13为第一步的应变划分结果示意图，图14为优化一定代数后的复合材料蒙皮单元铺层层数拟合曲面。

通过机群进行计算，每步计算时间需耗时10min，优化至160步，共需时26个小时，约1天。

在满足结构强度的要求下，未经优化的机翼中央翼盒蒙皮壁板质量240kg，未划分邻域只进行遗传算法优化后的重量为151.2kg,采用本方法优化后的重量变为173.4kg。与原结构减重达27.75％。当然这是理想的结果，但有效减重是可以实现的。结合流程图5图6，具体实现步骤描述如下：

步骤一：在有限元前处理软件Patran里建立机翼的有限元模型，划分好网格单元，模型的初始厚度0.8mm，铺层角度0°∶90°∶±45°＝1∶1∶1，上下翼面均为对称铺层，单层厚度均为0.1mm，模型共有33476个节点，47128个单元，模型大小为15M。优化的单元是上下翼面的单元共25004个。

步骤二：将初始模型的.bdf文件提交至Nastran计算得到.f06文件。

前处理主要是对利用Patran建好的有限元模型的信息进行提取，.bdf文件是由MSC.Patran生成的、供MSC.Nastran读取的文件，保存着在MSC.Patran中所建立的有限元模型的所有信息，MSC.Nastran就是根据.bdf文件来进行运算的。

将.bdf文件提交给Nastran计算，得到.f06文件，.f06文件是分析运算信息记录文件。

步骤三：编写C++程序，从.f06文件里读取单元的各项应变信息和角度信息，计算出每个单元的计算应变值ε_i，材料牌号T300/5208，碳纤维材料的力学性能参数如下表1所示。材料密度ρ＝1.6g/cm³ε_X＝3500，ε_Y＝-3200，ε_XY＝5000。用最小二乘法对计算应变值进行曲面拟合，

表1

记ε_i为单元计算的计算应变值，计算方法为根据选用的飞机蒙皮复合材料层合板的材料牌号来确定纵向/横向拉伸许用应变[ε_X]，纵向/横向压缩许用应变[ε_Y]和剪切许用应变[ε_XY]。一般所说的应变值是公知的物理概念，这里的应变值ε_i是用公式算出来的，所以称之为“计算应变”。公式里用到的ε_x、ε_y、ε_xy分别是纵向/横向拉伸应变、纵向/横向压缩应变、剪切应变，这个公式算出来的ε_i是距离三项许用应变的距离，所以用“计算应变”。复合材料单层的强度是计算层合板强度的基础，单层的强度分析包括单层的强度失效判据，单层的失效准则是用以判别单层在偏轴向应力作用或平面应力状态下是否失效的准则。本发明采用的单层失效准则是最大应变失效准则，该准则认为当单层在平面应力的任何应力状态下，单层正轴向的任何一个应变分量达到极限应变时，单层就失效。这个极限应变在单轴应力或纯剪应力状态下即是相应的基本强度所对应的应变。也由于单层的基本强度在纵向、横向、面内剪切向是不同的，故其对应的极限应变也是不同的。所以最大应变失效准则是由三个互不影响、各自独立的表达式组成，如下式：若只要满足下式中任何一个，单层即失效。

若.f06文件里提取的单元的拉、压、剪应变都满足该方向的拉/压和面内剪切许用应变，则该单元为基本层单元，令ε_i＝0。否则该单元为非基本层单元。对于非基本层的单元，计算其计算应变值ε_i：ε_X和ε_Y取该方向最大拉/压应变，ε_XY取面内最大剪切应变，[ε_X]和[ε_Y]根据ε_X和ε_Y是拉应变或是压应变取相应的许用应变，得到每个单元的计算应变值。为了使单元的计算应变值随空间位置连续变化，去除少数不连续的噪点，过滤面积较小的邻域，采用最小二乘法用曲面对计算应变值进行拟合。对计算应变值的曲面拟合重点是拟合后的曲面能反映出应变的大小相对关系，而不是着重于拟合出精确值。蒙皮单元分为四边形单元和三角形单元，为了忽略单元形状的影响，取三角形/四边形单元的形心位置来代替这个单元的位置作为单元样本点。

曲面的拟合次数是超参数，多项式次数在小于某个值时，次数越高，拟合效果越好，计算量也会相应增大。但是当拟合次数超越某个临界值再增加时，系数矩阵是病态的，误差可能不是单调递减而会变大，会出现不稳定的情况。经试验，选用8次，而且此时误差是满足需求的，既可以去噪点，拟合值也和真实值的变化也是符合的。即拟合出了单元计算应变值的相对大小，所以拟合次数取8次是目前试验中得到拟合效果较好而且计算量合适的。

对n个样本点(X_i，,Y_i,Z_i)，(i＝0,1,2...n)，将曲面方程记为f(x,y)＝P^T(x,y)A，其中A^T＝(a₀,a₁,…a_m)是曲面方程的系数。此处X_i，Y_i是单元的坐标，Z_i是单元的计算应变。

由最小二乘法的定义，误差平方和J(A)是所有单元已知特征向量的期望输出与实际输出的误差平方和的累加值，

表示为为使J(A)最小。

对于8次多项式拟合，

系数矩阵A为45*45的可逆矩阵。

求解如下所示45*45阶可逆矩阵，解出拟合曲面方程的系数A，把所求的系数A代入曲面方f(x,y)＝P^T(x,y)A即可求出拟合曲面方程。可计算每个单元计算应变值的拟合值：

下文中是按照每个单元的计算应变拟合值进行邻域划分。

如图10所示是在Pycharm软件里用Python编写程序调用Anaconda库绘制的三维立体图，x和y坐标表示单元形心的实际位置，z坐标表示单元的计算应变值。黑色的点为单元的实际计算应变值，灰色点为拟合后单元的计算应变值，可以看出拟合趋势和实际相符，翼根部分的计算应变值相对更大。以下计算都采用的拟合后的计算应变值。

按照拟合后的计算应变值来划分邻域，将单元按照拟合后的计算应变值从大到小排序。邻域划分的原则是四连通原理：如图2所示，和核心单元直接四连通的单元为核心单元的一阶邻域单元R1，由一阶邻域单元R1向外延伸得到和一阶邻域单元四连通的单元称为核心单元的二阶邻域单元R2，由二阶邻域单元R2向往延伸得到和二阶邻域单元四联通的单元称为核心单元的三阶邻域单元R3，以此类推。若两个单元有一条公共边则称这两个单元是连通的；相应的还有八连通原理：两个单元共顶点则称这两个单元是连通的，这里为了得到形状均匀的邻域，采用四连通原理。

把单元的拟合后的计算应变值的范围分成6-8层。先从拟合后计算应变值最大的单元开始遍历循环，若和其四联通的单元的拟合后的计算应变值范围相同，则划分为同一邻域，再以该邻域内新的单元作为起始单元重新开始遍历，直到该邻域内所有单元都遍历循环完成后。开始下一邻域的划分。

探究单元计算应变值范围的分布规律，做出统计图下图8所示，纵坐标为拟合后计算应变值，横坐标为大于这个计算应变值的单元数目。可以发现曲线的开始段斜率很大，如果按照等分拟合后计算应变值的标准即按照图8中5条虚线标示的方法来划分邻域范围会使单元数目分布不均匀，拟合后高计算应变值的单元数目过少，仅为总数目的1％，而有99％的单元拟合后计算应变值都小于最大应变的17％，显然不合理。

先提取按拟合后计算应变值降序排列后数目占单元总数的前5％的单元，再对剩下的单元按照拟合后计算应变值等分划分层，如图9所示，即为对水平线分割的曲线下半段部分单元按照拟合后计算应变值范围进行等分。可以得到较均匀的邻域范围。经多次试验，包括基本层一共有8层是可以得到比较合理的邻域分布的。中间值通过公式(i＝0,1,2,…,Slice)计算，Min为除去基本层外的单元的拟合后计算应变值最小值，Max为除去拟合后计算应变值最大的前5％单元后剩下单元中的最大拟合后计算应变值，经试验此处取Slice＝6。

上述步骤中已经将所有单元分到应变范围不同的8层里，下一步是在每一层中划分出不同的应变域。

从计算应变值最高的第8层开始即从i＝6开始，遍历所有单元，y为单元的计算应变值，当y>＝Y时取出该单元，并遍历和这个单元四连通的单元，若这些相邻单元的计算应变值y也满足y>＝Y，则这些单元为一个邻域的，给这个邻域编号赋值为1，重复遍历四连通单元这一步，直到周围单元的计算应变值y<Y时停止就找到了所有邻域编号为1的单元，取邻域1里所有单元拟合后计算应变值的最大值作为该邻域的计算应变值。重复上述操作，可以得到i＝6时对应得所有邻域。

减小i的值重复上述操作至i＝0。可以得到所有邻域，编号为0的为基本层，1为拟合后计算应变值最大的邻域，2为第二大，依次类推，保证邻域的最大拟合后计算应变值是随着邻域的编号从1开始递减。划分好邻域后，取该邻域内单元的最大拟合后计算应变值作为该邻域的计算应变值，邻域的平均主应变角取该邻域内所有单元主应变角的平均值即可。

由于在邻域的划分过程中采用了最小二乘法进行计算应变值的曲面拟合，所以划分出来的邻域没有噪点、边界也比较光顺，按照计算应变值范围划分出来的每块邻域的面积大体也都比较均匀。但是也有极少邻域的面积过小，从工艺方面考虑，所以在划分出所有邻域后，会设定一个面积阈值：S_min＝fliterratio*S_total，fliterratio＝0.01。若单个邻域的面积小于这个阈值，则再判定这块邻域的周围哪块邻域的计算应变值最大，将这块面积较小的邻域归并到周围计算应变值最大的邻域中。

每次邻域划分后上下翼面的邻域总块数和分层数以及面积阈值大小有关，考虑到制造工艺，邻域块数不宜过多，邻域面积也不宜过小。Slice＝6和fliterratio＝0.01是在多次试验中得到的比较合适的组合，每次划分后邻域总块数不超过50块。

划分好邻域后，根据之前从.f06文件里读取单元的角度信息，按照公式N是第i个邻域内包含的单元个数，i＝1,2…,n计算每块邻域的铺层角度。

步骤四：每个邻域作为遗传算法的操作对象，

本方法采用二进制编码方式，染色体由两部分组成，第一部分是每次选中的邻域块数，是可变长度的。根据每次划分邻域统计出的上下翼面邻域总块数确定，每次取邻域总块数的1/3，用N表示，表示块数的染色体长度用L₁表示，则N要满足下式。

2^L-1≤N<2^L (2)

由于二进制解码转换成十进制的数是离散的，为了得到和邻域块数N最近的编码块数，L₁由下式决定。

第二部分是铺层角度的变化范围，角度范围确定后，长度即可固定，本方法选取的角度是在每块邻域的平均主应变角的±128°范围内按1°为单位变化。考虑到染色体的长度不宜过长，目前采取的方式是随机投点方法，(即遗传算法的第一步种群初始化。)随机生成-128°至128°之间的整数角，加上每块邻域的平均主应变角作为角度搜索范围，所以表示角度的染色体长度是8位。其中第一位是符号位，后7位是角度值，如图11所示。

在遗传算法的优化迭代过程中，需要一个评价标准来评价一条染色体的优劣，即适应度，适应度可以是一个函数也可以是实验结果，是进行遗传操作的依据。合适的适应度可以给出更好的进化方向，适应度高的个体更有可能复制到下一代而适应度小的个体会被淘汰。提高优化效率。

根据优化的静强度约束条件，对机翼蒙皮的优化是在增重最小的情况下满足许用应变，适应度f和目标函数是正相关的。每次更新计算模型后调用Nastran计算得到应变值，需要一个能表示整体应变变化的值作为个体的评判标准。如下式所示，n为单元数目。

因为个体的选择是采用轮盘赌的方式，适应度大的个体被选中的概率更大，所以适应度函数应该取正值，且是求解最大值问题，需要对上述评价标准作转换，其转换关系如公式(5)所示，F即为适应度函数，C是实验中可能f可能取到的最大值。

本方法编写的遗传算法是标准遗传算法，只有三个遗传算子：选择、交叉和变异。遗传算法的进化过程中，基于最优保存策略保存最优个体，保留每一代的最优个体，为了保留最优个体不被破坏，因此最优个体不进行交叉变异等遗传操作而是直接保存到下一代种群中。其中交叉和变异算子在迭代过程中保持不变，交叉概率P_C＝0.7，变异概率P_m＝0.1。初始种群规模P＝100，进化代数为10。

步骤五：

优化中外循环给出邻域按计算应变值大小排序结果，优化没有将重量增量全部添加到计算应变值最大的邻域处，而是进入到内循环中，寻找最适宜添加的邻域块数，以达到计算应变值接近0，因为没有考虑每次添加的邻域之间的影响，所以要按照初始重量的一定比例逐步增加重量到选定的几块邻域上，而不能一次加上所有重量，避免算法的震荡，并可以保证算法的鲁棒性。

如前所述，每次增加机翼蒙皮质量Δw固定，占初始重量的10％。机翼蒙皮的初始厚度0.8mm，铺层角度0°:90°:±45°＝1:1:1，单层厚度0.1mm，上下翼面均是对称铺层。初始机翼中央翼盒蒙皮质量66.7kg，即每次增加0.667kg的材料重量。

若遗传算法每次输出的结果是给n块邻域铺层，那么计算应变值最大的n块邻域被选中，染色体解码的角度就是最佳铺层角。在主应变角的±128°范围内搜索最佳铺层角，并按照最佳铺层角进行铺层；最后用刚度等效原理把角度转换成工程用的铺层角：纵/横和剪切三向刚度是层合板的基本刚度特性，对称均衡板的纵/横和剪切三向刚度确定后，无论选择其铺层的各种角度，板的总厚度相同，因此在设计时可以按照实际的主应变角作为铺层角度，优化完成后再转换成0°∶90°∶±45°，等效后层合板的总厚度以及刚度特性保持不变，所述工程用的铺层角为0°∶90°∶⊥45°。设选中的第i个邻域的平均计算应变值是：N是第i个邻域内包含的单元个数，i＝1,2…n；

，邻域面积S_i(i＝1,2…,n)，被选中的每个邻域总应变ε_sum-i＝ε_avg-i*S_i，第i个被选中的邻域铺层重量分配：增加的铺层厚度：

步骤六：当所有单元的计算应变值都为0时，即所有单元的拉/压/剪切应变都满足许用应变时，即时优化结束，否则返回步骤一开始新一轮的外循环。

步骤七：当优化结束时，从优化模型的.bdf文件里提取单元的蒙皮厚度信息，按照厚度值递减给单元排序，从最厚的单元开始遍历，取这个单元为中心单元，厚度为这个邻域的厚度，若这个单元的四连通单元的厚度值和这个邻域的厚度值差小于单层材料的厚度0.125mm，则判定这个单元也是这个邻域的，直到邻域里所有单元的四连通单元厚度和邻域厚度差值超过0.125mm(两个单元的厚度差值如果不到一层的厚度就可以认为是同一厚度范围)停止遍历。重复上述操作。邻域编号从1开始为最厚，依次随邻域编号增大厚度递减。实施例中初始化厚度0.1mm是没有优化之前模型的单元的单层材料的厚度，是建模者给定的模型的参数，在做优化时采用工程上常用的复合材料比较规范的的单层厚度0.125mm，因此后续优化时以0.125mm为判定标准。

划分好厚度邻域后，按照和应变拟合相同的最小二乘法原理对厚度进行拟合，流程图如图7所示，拟合次数选8次。得到拟合后的厚度云图如图12所示。

本发明的优化方法中采用双循环优化流程，有效的解决了大设计变量的解耦问题。优化中采用“小增量填谷法”，避免了一般优化算法中出现的震荡，并且能够找到增加一个小重量最有利的铺层位置，以及铺层状态，因此这样优化出来的结果，能够保证重量最小，且满足结构强度要求。通过多次优化对比，发现修改每次重量增量，得到优化结果基本接近，说明应用该法进行优化，重复性高，具有较好的鲁棒性。

最后应该说明的是：以上实施例仅用以说明本发明的技术方案，但本发明的保护范围并不局限于此，任何熟悉本领域的技术人员在本发明揭露的技术范围内，可以轻易想到各种等效的修改或者替换，这些修改或者替换都应该涵盖在本发明的保护范围之内。

Claims (9)

1.一种基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述的飞机复合材料结构优化方法在结构优化之前，对结构蒙皮传力过程进行应变敏度分析，将计算应变值相近且空间位置相邻的单元分成一个“应变邻域”，在结构优化时以“应变邻域”为对象进行优化；所述空间位置相邻的单元是指根据四联通原理有公共边的单元；包含以下步骤：

步骤一、建立蒙皮的有限元计算模型；模型变量包括：蒙皮单元个数、蒙皮单元厚度、蒙皮单元铺层角度、蒙皮初始厚度、复合材料单层厚度、蒙皮单元铺层层数；

所述复合材料单层厚度就是每次外循环增加一层铺层的厚度，蒙皮单元铺层层数是指这个单元的总厚度，是多个单层厚度的厚度叠加；

步骤二、针对模型进行蒙皮受力的有限元分析；

步骤三、进行应变邻域划分：读取每个有限元单元对应的应变值大小，利用最小二乘法对每个单元的计算应变值进行曲面拟合，根据各个邻域所包含有限元单元的计算应变值的最大值，按照计算应变值从大到小对各应变邻域排序；

步骤四、利用遗传算法以“应变邻域”为对象进行结构优化：

以邻域为初始种群的遗传算法来计算最佳铺层块数、铺层角度和蒙皮单元铺层层数，形成初始计算n个个体的种群；种群中的每一个个体代表着某一种添加重量的邻域块数和相应添加的铺层角度的组合；计算种群中个体的计算应变值大小，进行遗传算法操作，基于最优保存策略保存最优个体，确定结构重量最佳的添加块数和相应的铺层角度；

步骤五、对每个应变邻域进行整体增重迭代：每块应变邻域内包含的每个单元增重相同，每块应变邻域内包含的每个单元的铺层角度也相同；

步骤六、当所有单元的拉伸/压缩/剪切应变都满足许用应变时模型优化结束，进入下一步；不满足，返回步骤一，更改模型变量，继续迭代；

步骤七、将上一步的优化后的模型中的蒙皮厚度结果进行二次曲面拟合；输出拟合后的厚度云图。

2.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述的飞机复合材料结构优化方法在划分“应变邻域”之前利用最小二乘法对计算应变值进行曲面拟合，用拟合后的计算应变值来作为划分邻域的依据，并在步骤六之后对优化后的蒙皮厚度进行二次曲面拟合。

3.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述步骤三中有限元单元对应的计算应变值的计算公式为：

，其中，/>为复合材料纵向/横向拉伸许用应变，/>为复合材料纵向/横向压缩许用应变，/>为复合材料剪切许用应变，/>、/>、分别是纵向/横向拉伸应变、纵向/横向压缩应变、剪切应变。

4.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：步骤三中每个单元的计算应变值是取单元上所有铺层的最大应变值中的最大值，该值大小表明单元是否失效，有限元分析中复合材料单层的强度失效准则采用最大应变失效准则；最大应变失效准则是根据复合材料单层拉伸/压缩/剪切三个方向的最大值来判断复合材料是否失效；对不满足失效准则的蒙皮单元进行优化，以蒙皮单元形状形心的位置作为蒙皮单元的样本点，对n个样本点，i =0,1,2...n，根据如下曲面方程进行拟合：

，其中/>是曲面方程的系数，/>，/>是单元的坐标，/>是单元的计算应变值，/>为系数矩阵，p ^T=(1, x, y, …, x ^m y ^m-n, x ^m-n y ⁿ), m=1,2,…,N, n=1,2,…N, m≥n, N为拟合次数。

5.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述步骤四中具体遗传算法步骤如下：

每块应变邻域的增重通过如下方式确定：每一步铺层选取总重量的5%-10%迭代，按照邻域的面积大小和应变范围确定每个邻域铺层重量分配和增加的铺层厚度：

若遗传算法每次输出的结果是给n块邻域铺层，那么应变最大的n块邻域被选中，染色体解码的角度就是最佳铺层角；

选中的第i个邻域的平均计算应变值公式如下：

，N是第i个邻域内包含的单元个数；

邻域面积为邻域内包含的各个单元的面积之和；

被选中的第i个邻域总应变，

第i个邻域铺层重量分配：，

第i个邻域增加的铺层厚度：，

其中，Δw为一次给所选中的所有应变邻域的增重的总和，为一次给选中的第i个应变邻域的增重，i是邻域编号，i = 1,2…n，j为第i个邻域内包含的第j个单元编号，/>为第j个单元计算的计算应变值，/>是所选用的复合材料的密度。

6.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：步骤四中铺层角度确定方法具体如下：

根据如下公式

N是第i个邻域内包含的单元个数，i=1,2…,n；是邻域内第j个单元的主应变角；为第i个邻域铺层角度；

计算出每个邻域的平均主应变角，在主应变角的范围内搜索最佳铺层角，并按照最佳铺层角进行铺层；最后用刚度等效原理将角度转换成工程用的铺层角；所述工程用的铺层角为/>。

7.根据权利要求1所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述步骤四中种群初始化采用随机投点方法。

8.根据权利要求6所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：所述将角度转换成工程用的铺层角具体操作如下：

按照实际的主应变角作为铺层角度，优化完成后再转换成，等效后层合板的总厚度以及刚度特性保持不变。

9.根据权利要求1至8任一项所述的基于应变邻域的飞机复合材料结构优化方法，其特征在于：在步骤三中还包含对划分后的邻域进行优化的操作，具体为：

在划分出所有邻域后，设定一个面积阈值，若某一邻域的面积小于这个阈值，则再判定该邻域的周围哪块邻域的应变值最大，将该邻域归并到周围应变值最大的邻域中。