组合式无人机兼具大型无人机长航时与小型无人机高灵活性的特点，备受世界各国关注^[1]。组合体无人机对长航时性能、机身轻量化及结构优化设计具有更高要求。随着结构优化理论的不断完善，在满足机身结构强度的前提下，机身轻量化设计得到长足发展，其中，拓扑优化具有多自由度设计空间特点，是结构优化最具前景的一个方向。主流拓扑优化方法包括均匀化方法^[2]、变密度方法^[3]和水平集方法^[4]等。在拓扑优化理论研究中，Guillaume等^[5]提出一种改进罚函数——渐进展开优化算法，对材料单元的去除进行优化分析；Novotny等^[6]建立了一组正则化公式，满足任何条件下的可行体积约束；刘宏亮等^[7]以固体各向同性惩罚微结构(solid isotropic microstructures with penalization，SIMP)模型为基础，建立了基于等几何分析-固体各向同性惩罚微结构(isogeometric analysis-solid isotropic microstructures with penalization, IGA-SIMP) 模型框架下的连续体拓扑优化方法。拓扑优化在飞行器领域的应用中，沈浩杰等^[8]利用Hypermesh仿真平台对机翼结构进行了拓扑优化分析，机翼重量降低26.9%；邱福生等^[9]基于气动-结构耦合作用，对三维拓扑优化方法进行了初步研究；兰剑英等^[10]利用拓扑优化方法对翼肋结构进行优化设计，机翼减重30%。在机翼结构的疲劳寿命分析方面；赵玉龙等^[11]对某型号机翼过渡梁和副油箱挂架开展了疲劳寿命分析，通过对后接头的耳片宽度调整优化，使该型号机翼满足寿命目标；高鹏飞等^[12]对比分析不同装配工艺对机盒疲劳耐久性能的影响特性。综上所述，机身细节疲劳寿命均基于已有型号对关键部位开展了详细研究，研究结果反馈至局部乃至整机设计的初始环节，这会增加重新优化设计的迭代次数，甚至会影响机身整体设计参数的重新设计，目前研究工作较少将疲劳寿命分析工作置于结构设计初期阶段。本文重新构建组合体机翼承力结构初始设计流程，首先结合飞行域进行全工况气动载荷分析，以极限气动载荷对机翼承力框架进行静力学仿真计算，结合Fe-safe平台对应力集中区域进行疲劳寿命特性分析，在机身满足疲劳寿命约束的前提下，基于拓扑优化理论开展机身结构优化设计。相较于传统设计方法，该设计可有效减少针对疲劳寿命而开展的全局迭代优化次数，在保证结构安全疲劳特性的同时，实现组合式无人机单机的轻量化设计，为组合式飞行器的初步设计工作提供技术支撑。

1 组合体无人机整机方案与单体构型设计 1.1 组合体无人机整机方案设计

1) 组合体整机构型设计。

组合体无人机采用“主体+单体”的总体构型方案，如图 1所示。其中: 主体为组合体的动力单元与飞行控制单元，单体通过连接分离机构并联于主体翼端，为组合体的主要作战单元。组合体无人机的飞行工况如表 1所示。

2) 作业模式设计。

组合体远程飞行至作战空域，在自主分离机构作用下实现有序分离后，单体通过自身气动外形与内部操控系统实现无动力滑翔，平稳着陆后实施后续作战任务，主体在单体全部分离结束后自主返航。

1.2 单机构型设计

通过单机构型初始设计，得到单机构型原始模型，再根据拓扑优化理论对原始模型开展承力结构轻量化研究工作。单机为平直翼构型设计，机翼承力框架材料选取航空铝合金，材料性能如表 2所示。

对单机进行初步结构设计时，在满足极限工况载荷的要求下，对机翼外部结构布局进行参数设定，构型参数如表 3所示。

单机内部选择由左、中、右3组肋板配合前、中、后3组翼梁的组合方案，前、中、后3组连接销固连于左侧肋板，各单机的升阻力载荷主要由连接销依次传递至机身，故左侧肋板进行加强设计，连接销、承力肋板与三组翼梁组成机翼的整体承力架构如图 2所示，后续优化区域主要针对单机承力结构，分别为前、中、后翼梁与中间肋板。

2 单体机翼拓扑优化模型

通过对单体机翼的气动特性分析获取许用飞行工况下的最大气动升阻合力，以最大气动升阻合力为输入进行单体机翼有限元静力分析，获取机翼框架在极限载荷工况时关键部位的应力与应变响应，分析过程的输出变量将作为疲劳寿命分析环节的输入变量，以Fe-safe疲劳分析软件为基础，通过仿真计算获取结构的疲劳分析结果，得到关键部位的疲劳耐久特性。然后将有限元分析结果联合Tosca拓扑优化软件进行拓扑优化，将最小应变与结构体积作为优化约束，设定体积最小为优化目标函数。用所保留的子单元重新描述结构的最优拓扑优化结果^[13]，最后通过2个环节的计算分析，获取机翼框架最优设计，机翼设计优化流程如图 3所示。

2.1 机翼气动载荷分析

本文气动载荷计算采用不可压定常三维连续方程和动量N-S(Navier-Stokes)方程^[9]。Descartes坐标系下的表达式分别为:

$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0, $

(1)

$ u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial x}+\frac{\mu}{\rho} \nabla^2 u, $

(2)

$ u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial y}+\frac{\mu}{\rho} \nabla^2 v, $

(3)

$ u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}=-\frac{1}{\rho} \frac{\partial p}{\partial z}+\frac{\mu}{\rho} \nabla^2 w . $

(4)

其中：u、v、w分别为x、y、z 3个坐标方向的流速，ρ为流体密度，μ为流体黏性系数，p为流体压力，∇²为Laplace算子。

气动升力L和气动阻力D计算方程可分别表示为：

$ L=\frac{1}{2} \rho V^2 S C_L, $

(5)

$ D=\frac{1}{2} \rho V^2 S C_D. $

(6)

其中：L为气动升力，D为气动阻力，ρ为空气密度，V为来流速度，S为机翼面积，C_L为升力系数，C_D为阻力系数。

无人机正常飞行工况下的空气Re为1.428×10⁶，机翼表面为湍流流动，应用Fluent仿真平台构建机翼有限单元网格模型，补充模型采用Menter-剪切应力传输(Menter's shear stress transfer, Menter's SST)模型，计算网格为非结构网格，机翼表面采用无滑移绝热固体壁面。飞行攻角取值范围为0°~32°，飞行高度H为10 km，空气密度ρ为0.414 kg/m³，动力黏度系数k为1.446×10^-5 Pa/s，飞行速度V为100 m/s，攻角变化引起升阻比变化曲线如图 4所示。

升阻比直接影响飞行器的气动特性^[14]，本文选取临界升阻比值为8，所对应的最大许用攻角为10°，机翼随攻角变化而产生的升力/阻力变化曲线如图 5所示。通过对不同攻角下机翼气动合力曲线分析，当飞行攻角为3°时，弹体机翼所产生的升力为355 N，略大于自身重力350 N，因此3°飞行攻角为最小巡航攻角，即组合体无人机的攻角范围为3°~10°。

由图 5可知，在最小飞行攻角时，机翼所受气动阻力F_D为36.22 N，气动升力F_L为355 N；在最大飞行攻角时，机翼所受气动阻力F_D为88.96 N，气动升力F_L为712.57 N，如表 4所示。通过仿真计算获得的气动合力，将作为后续机翼框架疲劳寿命分析与拓扑优化分析的输入载荷。

2.2 单体机翼载荷环境与疲劳寿命

1) 单体机翼载荷环境分析。

组合体无人机的任务剖面分为滑翔起飞、稳定巡航、自主分离、单体滑翔和主体返航5个阶段。滑翔起飞阶段与稳定巡航阶段为单体机翼承受气动载荷较大的2个主要阶段，其中滑翔起飞阶段的气动载荷最为复杂。单体机翼受到主体发动机所施加的推力T、垂直地面的重力G、来流方向的F_D与垂直于来流方向的F_L合力作用，并在以F_L为主要载荷的举升作用下稳定爬升至H后，开始过渡到巡航阶段，此时组合体无人机为相对稳定飞行状态，即T与F_D相互平衡，G与F_L相互平衡，并持续到下一飞行阶段。

由以上分析可知，在滑翔起飞与稳定巡航任务剖面内，机翼受到4个方向载荷的联合作用，其中滑翔起飞阶段为合力最大阶段，为验证机翼关键部位的疲劳寿命，需要针对上述2个阶段进行疲劳寿命分析。

2) 单体机翼细节疲劳分析方法。

细节疲劳额定值(detail fatigue rating，DFR)是民机机身结构细节本身所固有的疲劳性能特征值，因其适用范围广(循环寿命达到10⁵量级)且稳定可靠，被认为是飞机结构疲劳耐久性设计的主要方法之一^[15]。采用细节疲劳额定值对飞机结构进行疲劳分析，步骤如下：

(1) 确定高疲劳特征细节。首先确定所优化部件的高疲劳特征的细节部位，通常对所建立的三维模型进行有限元分析，根据计算结果确定危险细节部位，然后针对该部位进行疲劳寿命分析。

(2) 确定细节部位可靠性寿命与细节疲劳额定值。结构可靠性寿命表示如下：

$ N_{95 / 95}=\frac{\beta}{S_{\mathrm{T}} S_{\mathrm{C}} S_{\mathrm{R}}} . $

(7)

其中：β为材料特征寿命；S_T为部件系数，铝合金通常取值1.3；S_C为置信系数；S_R为可靠度系数。

细节疲劳额定值与N_95/95的数值关系为

$ \left\{\begin{array}{l} X=S^{\left(5-\lg N_{95 / 95}\right)}, \\ D=\frac{\sigma_{m_0}(1-R)}{0.94 \frac{\sigma_{m_0}}{\sigma_{\max }} X-(0.47 X-0.53)-R(0.47 X+0.53)} . \end{array}\right. $

(8)

其中：X为计算参数，S为材料S-N曲线斜率参数，σ_m0为材料特征参数，σ_max为材料最大应力参数，R为应力比，D为结构细节疲劳额定值。

(3) 确定当量载荷下的许用应力[σ_max]。

当结构细节疲劳额定值已确定且在应力比不变的前提下，结构的疲劳寿命能达到目标寿命所允许的最大载荷循环次数时的最大应力值，被称作当量载荷下的许用应力值，具体表示如下：

$ \left[\sigma_{\max }\right]=\frac{0.9 \sigma_{m_0} \cdot D}{(1-R) \cdot \sigma_{m_0} \cdot\left(\frac{N_{\mathrm{e}}}{N_0}\right)^{\frac{1}{S}}+(R-1) \cdot D} . $

(9)

其中：N₀为最大载荷循环数，取值1×10⁵；N_e为当量等辐载荷循环数，取值5×10⁴。

(4) 计算最大应力下的疲劳裕度MF，具体表示如下：

$ \mathrm{MF}=\frac{\left[\sigma_{\max }\right]}{\sigma_{\max }}-1. $

(10)

2.3 基于SIMP单体机翼拓扑优化模型

1) SIMP算法连续体拓扑优化。

SIMP算法是针对各同性材料，基于幂律法则的具有罚因子的计算方法^[16]，使优化区域内的每个单元与1个密度相关的连续函数关联，关联函数的取值范围限于0~1，将单元的密度与弹性模量以某种关系式相关联，表示如下：

$E_i=x^j E_0 . $

(11)

其中：E_i为第i个单元的惩罚函数；x为单元密度，j为罚因子；E₀为材料弹性模量。

以连续体体积为约束条件的拓扑优化数学模型可表示如下：

$\left\{\begin{array}{l} \min c(x, \boldsymbol{u})=\boldsymbol{f}^{\mathrm{T}} \boldsymbol{u}=\sum\limits_e \boldsymbol{u}_e^{\mathrm{T}}\left(x_e \boldsymbol{k}_0\right) \boldsymbol{u}_e \\ \text { s. t. } \sum\limits_e x_e v_e \leqslant V \\ 0 \leqslant x_{\min } \leqslant x_e \leqslant x_{\max } \end{array}\right. . $

(12)

其中：c为柔顺度(目标函数)；u为位移向量；f为载荷向量；u_e为单元e位移矩阵；x_e为单元e设计变量；k₀为单元刚度矩阵；v_e为单元e对应的体积变量；V为体积约束变量；x_max为单元密度最大值，取值1；x_min为单元密度最小值，为防止出现刚度奇异现象，取值0.000 1。

2) 优化过程灵敏度分析模型。

灵敏度分析是实现最优求解的关键步骤，式(12)中c对单元设计变量x_i的灵敏度表示如下：

$ \frac{\partial c}{\partial x_i}=\frac{\partial \boldsymbol{f}^{\mathrm{T}}}{\partial x_i} \boldsymbol{u}+\boldsymbol{f}^{\mathrm{T}} \frac{\partial \boldsymbol{u}}{\partial x_i} . $

(13)

其中 x_i即第i个单元的单元密度。引入平衡方程K·u=f，即约束条件λ(f-K·u)=0，其中：K为刚度矩阵，λ取值无限制。为求解方便，取值λ=u后，式(13)可化简为

$ \frac{\partial c}{\partial x_i}=-\boldsymbol{u}^{\mathrm{T}} \frac{\partial \boldsymbol{K}}{\partial x_i} \boldsymbol{u}. $

(14)

优化区域总体积对单元的灵敏度模型，则有

$ \frac{\partial V}{\partial x_i}=\int \mathrm{d} V=v_e. $

(15)

3) 基于最优准则法的拓扑优化求解。

式(12)为典型体积约束下的拓扑优化模型，由于最优准则法(optimality criteria method，OC)具有高效、快速的特点，通常被认为是解决此类优化问题的主要方法之一。本文基于Kuhn-Tucker最优求解准则^[17]，将上述带有体积约束条件拓扑优化模型通过Lagrange乘子构造成Lagrange函数，表示如下：

$\begin{gathered} L=c+\lambda\left(\sum\limits_{i=1}^n v_i x_i-V^*\right)+ \\ \lambda_1\left(x_{\min }-x\right)+\lambda_2\left(x-x_{\max }\right) . \end{gathered} $

(16)

其中：v_i为第i个单元的微元体积，V^*为最大体积约束变量，λ、λ₁和λ₂分别为各设计变量约束条件下的Lagrange乘子。当式(16)取极值时，需满足Kuhn-Tucker必要条件，具体表示如下：

$ \left\{\begin{array}{l} \frac{\partial L}{\partial x_i}=\frac{\partial c}{\partial x_i}+\lambda v_i-\lambda_1+\lambda_2=0 ; \\ \sum\limits_{i=1}^n v_i x_i-V^* \leqslant 0 ; \\ x_{\min } \leqslant x_i \leqslant 1 ; \\ \lambda\left(\sum\limits_{i=1}^n v_i x_i-V^*\right)=0, \lambda \geqslant 0 ; \\ \lambda_1\left(x_{\min }-x_i\right)=0, \lambda_1 \geqslant 0 ; \\ \lambda_2\left(x_i-1\right)=0, \lambda_2 \geqslant 0 . \end{array}\right. $

(17)

将设计变量进行变换，可表示为

$B_i^k=\frac{\frac{\partial c}{\partial x_i}}{-\lambda v_i}=1. $

(18)

其中：B_i为第i个单元的驻值条件，k为第k个迭代步。当x_i满足x_i≠x_min≠1时，λ₁与λ₂均等于0，则迭代方程可表示为

$ x_i^{k+1}=\left(B_i^k\right)^\eta x_i^k. $

(19)

最优解迭代方程可表示为

$ x_i^{\text {new }}=\left\{\begin{array}{l} \max \left(x_{\min }, x_i-m\right), \\ x_i\left(B_i^k\right)^\eta, \\ \min \left(1, x_i+m\right) . \end{array}\right. $

(20)

其中：m为移动界限，η为数值阻尼系数。当x_i(B_i^k)^η ≤max(x_min, x_i-m)时，x_i^new取值max(x_min, x_i-m)；当x_i(B_i^k)^η≤min(x_min, x_i-m)时，x_i^new取值x_i(B_i^k)^η；当max(x_min, x_i-m) < x_i(B_i^k)^η时，x_i^new取值min(1, x_i+m)。

2.4 单体机翼模型

步骤1  建立有限元模型。

为方便有限元模型计算以及在不影响计算精度的前提下提高计算效率，对所建立模型进行如下简化：

1) 蒙皮主要在气动特性分析时承担升阻力，但在机翼框架拓扑分析环节，不作为主要传力构件，故不考虑蒙皮在拓扑优化过程中的影响作用。

2) 忽略蒙皮与机翼框架铆接处的影响。

3) 忽略翼梁和翼肋之间螺栓连接的影响。

4) 忽略加工误差导致的机翼之间销-孔连接的缝隙影响。

通过对所建立的机翼框架模型进行合理简化，在Abaqus前处理模块中对机翼框架进行离散化处理，所选网格类型为六面体扫掠网格，得到机翼框架有限元模型如图 6所示。

步骤2   设定机翼框架的载荷及边界条件。

实际飞行工况中，气体通过流过上翼面与下翼面时的气动压差产生升力，本文以极限工况(升阻比10)与巡航状态(升阻比3)下的气动力分别作为输入载荷，在F_L与F_D加载环节，极限载荷定义为飞行攻角10°下的气动合力，在下翼均匀加载气动升力载荷362.57 N，在机翼前缘加载均布阻力载荷88.96°；常值载荷定义为飞行攻角3°下的气动合力，在下翼面加载均布净升力载荷5 N，在机翼前缘加载均布阻力载荷36.22°，单体机翼与机身通过连接销进行固连，且在机身坐标系下机翼相对机身为静止状态，故本文针对连接销进行边界条件设置时，默认机身为刚体固定状态，固定销选择完全固定约束。最大应力位置即为疲劳寿命敏感区域，极限载荷下最大应力为71.99 MPa，如图 7a所示；常值载荷下最大应力为0.63 MPa，如图 7b所示。

步骤3   分析关键部位疲劳寿命。

疲劳寿命是指在循环加载情况下，材料产生疲劳破坏所需的应力或应变的循环次数^[18]，结合2.2节中的机翼载荷环境分析，分别对极限载荷与常值载荷作用下的关键部位开展疲劳寿命分析，设循环载荷加载频率为1 Hz的正弦波函数，基于循环载荷作用下的材料极限寿命，针对步骤2分析得出的敏感区域开展疲劳寿命计算，循环载荷谱如图 8所示。

根据图 7有限元计算结果确定疲劳寿命敏感区域，在最大载荷力作用下，敏感区域处的应力、疲劳寿命与细节疲劳额定值之间的关系^[15]可表示如下：

$ \left\{\begin{array}{l} Z=\frac{(1-R)\left(\sigma_0-0.53 D\right) \sigma_{\max }}{D\left[0.94 \sigma_0-0.47(1+R) \sigma_{\max }\right]}, \\ N_{95 / 95}=10^{\left(5-\frac{\mathrm{lg} Z}{\mathrm{lg}S}\right)} . \end{array}\right. $

(21)

其中：各参数根据2.2节仿真计算结果与文[19]进行取值，计算结果如表 5所示。

基于上述有限元仿真计算结果，联合Fe-safe疲劳寿命分析平台，对机翼框架进行全域疲劳寿命分析，计算结果如图 9所示。其中，在极限载荷循环作用下，敏感区域最低疲劳寿命N_95/95为132 434次，常值载荷作用下，最低疲劳寿命N_95/95为15 848 931次，均达到材料许用循环寿命10⁵受次循环^[20]。组合体无人机设计安全标准为安全起飞1 000次，其中正常起飞(攻角5°~8°)占比90%，极限起飞(攻角10°)占比10%，爬升角等价于最大攻角，爬升时间为1 000 s，由极限载荷敏感区寿命计算可得极限起飞为132架次，即安全标准起飞1 320架次，满足组合体无人机设计标准。

步骤4   分析机翼框架拓扑优化。

根据SIMP算法对机翼框架开展拓扑优化分析，在分析设置中，由于左翼肋与右翼肋为机翼内部装置的固连肋板(见图 2右肋板)，对其结构参数具有特殊要求，以结构布局最优为目标的拓扑优化会对实际工程指标造成影响，因此左翼肋与右翼肋不在优化范围内。同时，由于前、中、后梁，以及中翼肋受力权重比例不同，本文选择分别对上述4组优化单元进行拓扑分析。通过有限元分析的结果，联合Tosca拓扑优化模块并建立优化模型，以最小应变能作为目标函数，选择各优化单元的体积优化指标作为约束条件，优化设置迭代步数为30步，优化过程及结果数据如图 10所示。

通过对机翼框架初始模型的拓扑优化分析，根据翼尖升力方向最大变形量D_tip应小于5% b的刚度要求^[21]，优化后的机翼框架翼尖最大变形量为0.27 mm，远小于许用刚度标准。翼端在循环载荷作用下的最大位移点的位移量如图 11所示，优化后的机翼框架体积比如表 6所示。

3 结论

本文根据未来组合体无人机的应用需求，进行组合体无人机机翼框架的初步结构设计，针对空间构型冗余的问题，综合考虑其气动载荷加载与应力集中区域的疲劳寿命问题，以SIMP算法下的连续体拓扑优化方法为技术支撑，建立机翼框架的有限元模型，以该模型提供的应力响应作为疲劳寿命分析与拓扑优化的桥梁，建立了基于气动载荷下机翼敏感区域疲劳寿命的拓扑优化方法，并通过最小应变能作为目标函数，机翼体积比作为约束条件，经计算，全工况下单体机翼最低疲劳裕度达到3.31，优化后体积减少35%，翼端最大位移仅增大0.01 mm，最终在考虑极端气动载荷作用下，同时率先分析敏感区域疲劳寿命，达到机翼框架的轻量化设计要求。

本文所采取的基于实际飞行工况极端气动载荷作用下的机翼敏感区域疲劳寿命分析与拓扑优化设计流程，为组合体无人机设计提供了更全面、高效的构型设计思路与方案，未来研究可在本文设计流程基础上，增加多目标约束优化设计，建立多学科优化设计方法，实现机翼满足更多约束条件下的轻量化设计目标。