鱼骨柔性机翼提供了一个替代传统离散襟翼和其他主动变弯度机翼的设计方案，是由Woods和Friswell[2]受自然界鱼类构造所启发设计的，只需要很小的驱动载荷就能实现后缘大角度变形，是一款优秀的变弯度机翼结构。由于后缘变形大，并且存在机翼结构主动变形，其存在显著的气动弹性问题。传统线性气动弹性分析方法并不适用于鱼骨柔性机翼的动力学分析，而需要与之匹配的气动弹性分析方法。
因此，本文充分考虑鱼骨柔性翼段大变形时气动弹性的问题，运用非线性有限元(FEM)和曲面涡格法(VLM)为基础的非线性静气动弹性分析方法，完成大变形下的鱼骨柔性翼段非线性静气动弹性分析，并将结果与传统线性静气动弹性结果进行对比。
1 相关工作 早在1920年，工程师们就设计了一种刚性铰链结构，这种结构通过各种机构运动和柔性蒙皮的组合，可以显著地改变机翼弯曲度。最近主动变弯度机翼设计已经成为一个非常活跃的研究领域，仅在过去20年里就有60多个设计方案被提出。Barbarino等[3]对所做工作进行了全面概述。这些方案中，基于柔性结构的设计方案比基于刚性铰链的方法更常见，主要是考虑到基于刚性铰链的设计具有结构复杂、质量笨重、结构可靠性差等缺点。图 1显示了主动变弯度机翼的几个典型设计[2]。可以从中看出, 随着时间的推移变形结构设计从复杂刚性铰链结构向简单柔性结构变化的发展规律。
 |
| 图 1 主动变弯度机翼设计[2] Fig. 1 Design of active camber morphing airfoil[2] |
| 图选项 |
鱼骨柔性机翼只要很小的驱动能量就能实现翼剖面中弦线大幅度连续变化，如图 2所示[3]。其主结构是一个脊柱主梁和连接在梁上的一系列纵梁，结构外侧是一层柔性蒙皮。驱动器通过钢索输出扭矩给脊柱主梁。这种结构具有较高的各向异性。弦向的弯曲刚度相对较低，以减少驱动载荷，翼型厚度方向刚度较高，以维持承载能力。这种设计减少了非变形区域承受的载荷和随之产生的应力。Woods提出鱼骨柔性机翼概念后，跟进了大量的研究。首先，使用XFOIL低精度的面元法与OpenFORM高精度的计算流体力学方法分别对鱼骨变形结构的气动性能进行研究，发现鱼骨柔性机翼相比传统的后缘襟翼设计具有良好的气动外形[4]。然后，制作简单的模型，并对模型进行了一系列地面变形实验，结果表明，鱼骨柔性机翼具有变形能力优越、驱动载荷小的优点[5]。最后，进行了吹风实验，通过风洞试验发现，与传统的尾缘襟翼相比，鱼骨柔性机翼气动效率有所提高，在等效升力条件下升阻比提高25%[6]。在风速为20 m/s和迎角为0°的条件下，鱼骨柔性机翼对比传统固定机翼升力系数提高0.72。
 |
| 图 2 鱼骨结构示意图[3] Fig. 2 Schematic diagram of fishbone structure[3] |
| 图选项 |
鱼骨结构的详细设计是一个困难的、多学科耦合的问题。主要是由于鱼骨的特殊结构：柔性结构弦向弯曲刚度低，变形时气动载荷变化大，气动弹性问题显著。因此，Woods和Friswell[7]提出了引入弛豫参数提高收敛性的静气动弹性分析方法，用来研究鱼骨机翼的流固耦合问题。结构部分以欧拉-伯努利梁理论的结构分析模型，将钢索建模为线性刚度单元，其驱动力简化为在翼尖对脊柱主梁上的弯矩。气动部分用XFOIL软件求解得到了气动载荷。基于这种方法，Woods和Friswell[8]又在后续研究中用多目标优化的遗传算法对鱼骨结构进行了详细的优化设计。
可以看出，Woods对于鱼骨结构的分析设计和制作这一套过程都是详细并且完备的，可以论证这样的鱼骨结构是真实有效并且可靠的。但是值得注意的是，文献[5]进行的结构实验中，鱼骨结构后缘位移已经超过整体弦长的15%，从材料力学角度来说，这是属于结构大变形的范畴，其结构变形具有很强的几何非线性。而Woods在结构建模、有限元分析及后续的气动弹性研究中，结构部分的理论基础都是采用了线性的欧拉梁理论和有限元理论。虽然后面通过引入鱼骨结构脊柱主梁的长度不变这一非线性条件对模型进行修正，但仍然没有考虑到大变形带来的随动载荷的影响。这就跟其实际变形存在一定误差。因此，从准确反映实际鱼骨变形的角度上来说，鱼骨结构应该采用非线性大变形结构分析方法来研究其结构特性和气动弹性特性。
另外，针对大变形结构的气动弹性问题，现阶段研究大多以面向理论研究为目标的简单梁式模型结合二维气动理论为主，比较成熟的有：美国密西根大学基于应变梁方程和入流理论开发的非线性气动弹性仿真工具UM/NAST[9]，能够进行简单梁或梁组合模型的非线性气动弹性仿真分析，应用于高空长航时无人机设计；麻省理工学院基于大变形梁理论与升力线模型结合，发展的ASWING综合设计分析方法，与UM/NAST框架类似[10]。但是非线性梁理论作为结构求解器限制了复杂模型的适用性，二维气动力计算方法作为气动求解器也限制了其计算精度。为了提高计算精度，CFD/CSD松耦合迭代方法也被应用于计算大变形结构非线性气动弹性分析。党会学等[11]采用CFD/CSD方法对大展现比柔性机翼进行非线性气动弹性分析，并用黄金分割加速松耦合计算的收敛，结果表明，几何非线性对静气动弹性大变形起着重要作用。基于CFD/CSD计算气动弹性计算精度高，但是带来了计算耗费巨大的缺点。涡格法作为基于势流方程的气动力计算方法，建模简单，计算效率高，在静气动弹性分析中有较为明显的优势。特别的是，涡格法在许多方面被证明是非常准确的，并在大柔性构型的气动弹性分析中得到了广泛的应用[12]。基于涡格法的计算具有耗费少、精度相对较高的优势。北京航空航天大学的谢长川等基于曲面涡格法和几何非线性有限元法，建立了一种大柔性机翼静气动弹性分析方法[13]，并进一步设计了地面振动和风洞试验来验证该方法[14]，计算值与风洞实测值趋于一致；在此基础上，又进行了大柔性飞机全机的静气动弹性分析[15]。可以看出，基于非线性有限元和曲面涡格法的大变形结构静气动弹性分析方法适合工程实践中快速、高效的分析要求，其既可以用于飞机设计的初级阶段，又可以用于飞机设计的详细阶段。
2 理论分析 2.1 结构几何非线性 传统线性气动弹性分析中，结构分析均是基于小变形假设的线性结构分析。鱼骨主动变弯度结构弦向刚度小，在气动载荷作用下会发生较大的弯曲变形，其变形往往不再符合线性小变形的假设。这样一来，几何关系和平衡方程都是非线性的，构成了典型的几何非线性问题[16-17]。
有许多方法能解决几何非线性结构分析，如基于位移的有限元法、基于应变的有限元法和几何精确梁理论。这里使用位移基础有限元的更新拉格朗日公式，即
 | (1) |
式中：单元矩阵由2部分构成，即线性刚度项tkL和非线性刚度项tkNL；u为结构位移；t+ΔtQ为新增量步的外载荷；tF为等效结构内力。线性刚度项只与结构本身的刚度和质量有关，而非线性刚度与结构变形构型、载荷形式等因素相关。
另一个在几何非线性问题中必须引起重视的问题就是“随动载荷”效应，即外载荷是变形的函数，尤其表现在其作用方向会随着变形一起改变。图 3以悬臂梁为例，给出了受力方向的变化情况。梁长度为L，忽略轴力引起的长度变化，引入弧长坐标s，自由端作用始终沿着y轴正方向的集中力P。可以看出，随着P的增加，挠度w也会增加，而轴向变形u也会增加。鱼骨柔性机翼作为典型的非常柔性结构，在对其进行气动弹性研究时，考虑气动压力作为随动载荷更为合理，因为气动载荷的一部分将因结构的变形成为侧向力，当变形足够大时，这部分力将不可忽略，会影响最终的气动弹性结果。
 |
| 图 3 一端受垂直随动载荷的悬臂梁 Fig. 3 Cantilever beam with one end subjected to vertical follower load |
| 图选项 |
2.2 曲面涡格法 在传统气动弹性分析中，气动力是由平面涡格法进行计算的。平面涡格法求解过程中，机翼被视为一个平面，马蹄涡网格被布置在气动面上。利用Biot-Savart定理计算气动网格控制点上每个马蹄涡的诱导速度。控制点处满足物面不穿透条件得到有关马蹄涡涡强的线性代数方程，进而求解气动面压强差得到升力分布。
定常曲面涡格法将机翼中弧面作为气动面，在曲面气动面上布置涡环基本解进行气动力的求解。曲面涡格法涡格划分如图 4所示[15]。在曲面气动面上布置空间涡环，能够较为准确地反映其真实气动面形状和流场状态，满足大变形鱼骨柔性翼段气动力建模的需求。
 |
| 图 4 薄翼型机翼曲面涡格模型[15] Fig. 4 Non-planar vortex lattice model of thin airfoil[15] |
| 图选项 |
涡环在任意位置引起的速度由Biot-Savart定律计算。
 | (2) |
式中：Γ为涡环强度；W为气动力影响系数；V为当地流速度；r为任意点到涡环控制点的位移；C为涡环长度。
将其代入物面不穿透边界条件，得到线性代数方程组为
 | (3) |
式中：A为气动力影响系数矩阵；RHS为列向量，其元素由式(4)得到
 | (4) |
其中：元素包含无穷远来流速度V∞和控制点法向ni的信息，含义是物面不穿透。涡格法细节参见文献[18]中详细描述。
得到各涡环的强度后，利用库塔-儒科夫斯基定理计算气动载荷:
 | (5) |
式中：ΔLij为第i行第j列气动网格升力；ρ为空气密度；Δyij为第i行第j列气动网格等效弦长。
在几何非线性气动弹性分析中，需要在迭代过程中对曲面涡格法的气动模型进行更新，使其与结构变形保持一致。虽然气动力的计算是线性的，但是几何精确的边界条件将非线性引入了气动弹性问题。
2.3 曲面大变形插值方法 传统线性气动弹性分析所用的二维插值方法无限平板样条(IPS)及其三维推广的薄板样条(TPS)方法运用广泛，是商用软件如NASTRAN[19]与ZAERO[20]的标准插值方法。鱼骨柔性机翼曲面气动网格和节点需随机翼变形成空间曲面，传统气动弹性分析所用的二维插值不再适用，需要三维曲面的样条插值方法确定其曲面变形位置信息。谢长川和杨超将传统二维插值IPS方法推广到多维空间[21]。
设在N维空间的n个给定向量Xi={xi1, xi2, …, xiN}(i=1, 2, …, n)，及其对应的多元函数Wi={wi1, wi2, …, wiM}(i=1, 2, …, n)。对W的每个分量按照IPS方法建立如下插值函数：
 | (6) |
式中：
将上述插值理论应用于气动弹性特性分析的结构/气动界面耦合，就能够将结构界面的运动传递到气动界面上。
令n个结构节点坐标为XS，其变形矢量为US，需得到对应m个气动网格节点XA处的变形矢量UA。根据结构节点坐标写出插值函数系数方程为
 | (7) |
式中：AS为结构节点坐标构成的系数矩阵；WS为结构位移坐标构成的系数矩阵；C为中间变形量系数矩阵。
根据式(7)求得C为
 | (8) |
则气动网格节点处的位移矢量可表示为
 | (9) |
式中：AA为AS去掉前4列元素得到。因为WS前4行元素均为0，式(9)可化为
 | (10) |
式中：G为结构节点与气动节点之间的位移插值矩阵，由矩阵AAAS-1去掉前4列元素得到。
气动弹性分析中另一个重要问题是要将气动界面的载荷信息转化到结构界面之上，这一过程称为力插值。气动弹性分析中的力插值不采用一般结构分析中的静力学等效或变形等效原理，而是要满足所谓结构等效性，也就是在结构任意可能的虚位移情况下2套等效的载荷系统虚功相等。根据虚位移的任意性，可以推出气动模型与结构模型之间的力插值关系式为
 | (11) |
式中：FS为节点力；FA为气动力。
2.4 静气动弹性计算流程 鱼骨柔性机翼作为典型的柔性结构，相比传统机翼气动弹性问题更为显著。鱼骨柔性机翼受气流作用的影响发生变形，而气动力又与机翼的形状有关，两者是相互影响、相互制约的非线性关系。本文采用定常曲面涡格法与非线性结构有限元，利用分析迭代计算的方法得到鱼骨柔性机翼的静气动弹性非线性平衡状态，具体流程如图 5所示。传统线性静气动弹性分析方法与本文提供的非线性静气动弹性分析方法流固耦合的迭代流程是相同的，同样是先求解气动载荷下的结构静力分析，再通过插值方法更新气动载荷，最终判断收敛条件得出计算结果。不同之处在于：传统线性静气动弹性分析方法的气动力求解是基于平面涡格法，结构求解时基于线性有限元分析，对应的插值方法为薄板样条插值。
 |
| 图 5 鱼骨结构非线性静气动弹性计算流程 Fig. 5 Calculation process of nonlinear static aeroelasticity analysis of fishbone structure |
| 图选项 |
3 算例计算 3.1 模型验证 根据文献[7]提供的模型数据，建立鱼骨柔性翼段的有限元模型，具体参数如表 1所示。
表 1 鱼骨结构模型参数 Table 1 Structural parameters of fishbone structure model
| 参数 | 数值 |
| 基础翼型 | NACA0012 |
| 弦长c/mm | 305 |
| 展长b/mm | 150 |
| 变形起始位置/mm | 107 |
| 变形结束位置/mm | 260 |
| 纵墙数量 | 14 |
| 纵墙厚度/mm | 0.8 |
| 蒙皮厚度/mm | 1.5 |
| 脊柱主梁厚度/mm | 2 |
| 纵梁弹性模量/GPa | 2.14 |
| 脊柱主梁弹性模量/GPa | 2.14 |
| 蒙皮弹性模量/MPa | 4.56 |
| 钢索弹性模量/GPa | 131 |
表选项
采用商用软件Pratran2018对鱼骨柔性翼段进行有限元建模。前缘翼盒和后缘翼尖采用HEX8单元，材料选择为铝合金，弹性模量为207 GPa。鱼骨部分和蒙皮采用板单元QUAT4建模。钢索布置于翼段中间位置。边界条件为前缘翼盒固支，有限元建模如图 6所示。
 |
| 图 6 鱼骨结构机翼有限元模型 Fig. 6 Finite element model of wing with fishbone structure |
| 图选项 |
有限元建模完成后，对其进行模态分析，得到前6阶模态，如图 7所示。
 |
| 图 7 鱼骨柔性翼段模态 Fig. 7 Modes of fishbone flexible wing |
| 图选项 |
模态结果如表 2所示，结果显示，鱼骨柔性翼段模态变形符合预期。
表 2 鱼骨柔性翼段模态频率 Table 2 Modal frequency of fishbone flexible wing
| 阶数 | 模态名称 | 频率/Hz |
| 1 | 一阶弯曲 | 8.314 1 |
| 2 | 一阶扭转 | 15.359 |
| 3 | 二阶弯曲 | 22.495 |
| 4 | 二阶扭转 | 44.228 |
| 5 | 三阶弯曲 | 45.852 |
| 6 | 面内模态 | 61.048 |
表选项
为了进一步验证鱼骨柔性翼段有限元建模，在不同加载条件下对模型进行线性有限元分析和非线性有限元分析，并将分析结果与文献[7]中的实验结果进行对比。
第1种加载方式为：鱼骨柔性翼段前端固定，在后缘中点处，垂直向下施加外力，如图 8所示。加载后鱼骨柔性翼段变形仿真结果如图 9和图 10所示。
 |
| 图 8 鱼骨柔性翼段后缘外力加载示意图 Fig. 8 Schematic diagram of external force on the trailing edge of fishbone flexible wing |
| 图选项 |
 |
| 图 9 鱼骨柔性翼段后缘z方向位移与外力的关系 Fig. 9 Relationship between z-direction displacement of the trailing edge of fishbone flexible wing and external force |
| 图选项 |
 |
| 图 10 鱼骨柔性翼段外力加载变形结果 Fig. 10 Deformation of fishbone flexible wing underexternal force |
| 图选项 |
图 9中将后缘加载有限元计算结果和文献实验值[7]进行了比较。计算结果表明，对于后缘z方向位移，线性有限元结果和非线性有限元结果变化趋势是一致的，随着加载力的增大而增大。非线性有限元仿真结果和线性有限元结果差距小于0.5%，并且与文献实验值[7]符合得很好。此时后缘z方向最大位移为翼型弦长的7.51%，鱼骨主梁偏转角为12.52°。可以看出，鱼骨变形比较小的情况，看出后缘z方向位移和加载外力之间的关系基本呈线性的，有限元非线性结果和线性结果差距不大。随着加载力的进一步增加，鱼骨柔性翼段进一步发生弯曲，鱼骨柔性翼段有限元变形分析的非线性结果如图 10所示。2种分析方法得到的结果差别主要体现在鱼骨柔性机翼的x方向。线性分析结果在x方向的位移为0，随着加载力增加，鱼骨结构只向z方向发生变形，实际上拉长了鱼骨柔性翼段中弦线的长度，而非线性分析则保持了中弦线的长度不变。随着加载力的增加，这2种分析结果差距会越来越大。外载荷为10 N时，2种分析方法得到的后缘点位移差为2.455 mm，是弦长的0.805%，两者几乎重合。外载荷为20 N时，2种分析方法得到的后缘点位移差为9.38 mm，是弦长的3.07%，此时鱼骨柔性翼段后缘z方向位移为-53.18 mm。可以看出，当结构变形比较大时，非线性有限元结果才能准确反映鱼骨柔性翼段真实变形结果。
第2种加载方式为：模拟鱼骨柔性翼段真实作动情况。通过给转轴施加扭矩，拉动钢索，使得鱼骨结构发生变形，并且同样通过线性和非线性2种有限元方法研究鱼骨的翼尖变形情况。加载示意图如图 11所示，鱼骨结构变形仿真结果如图 12和图 13所示。
 |
| 图 11 鱼骨柔性翼段外力矩加载示意图 Fig. 11 Schematic diagram of external moment of force on fishbone flexible wing |
| 图选项 |
 |
| 图 12 鱼骨柔性翼段后缘z方向位移与外力矩的关系 Fig. 12 Relationship between z-direction displacement of the trailing edge of fishbone flexible wing and external moment of force |
| 图选项 |
 |
| 图 13 鱼骨柔性翼段外力矩加载变形结果 Fig. 13 Deformation of fishbone flexible wing under external moment of force |
| 图选项 |
如图 12所示，考虑真实作动情况下鱼骨柔性翼段的线性有限元和非线性有限元结果最大相差0.45%，并且都与文献实验值[7]符合得很好，最大变形为翼型弦长的5.77%。这与在后缘加载力的情况得到的结果趋势一致，进一步验证了有限元的准确性。从图 13进一步看出，不管是后缘加载还是真实作动情况，线性分析依然不考虑x方向的变形。当驱动力矩为2.0 N·m时，后缘z方向位移为-45.33 mm，是弦长的14.86%，此时鱼骨主梁偏转角为21.04°，仍然在鱼骨柔性翼段合理的变形范围中。此时，非线性分析结果与线性分析结果得到后缘位移差为6.83 mm，是弦长的2.24%，不能忽略。
通过2种不同加载形式的计算验证，本文所建立的模型与文献[7]模型具有很好的一致性。在后缘挠度比较大的情况下，非线性分析更符合翼段实际的变形情况。真实作动情况下，当后缘z方向位移到达弦长14.86%时，弦向x方向变形为弦长2.24%，其结果值得关注。
3.2 气动力分析 本节运用曲面涡格法分别计算3.1节中真实作动力矩加载下的鱼骨柔性翼段线性变形和非线性变形对应的气动载荷，并将计算结果与XFOIL软件进行对比。对鱼骨柔性翼段气动面进行网格划分，如图 14所示。气动网格数量为10×20，其中弦向20个网格，展向10个网格。气动网格前7×10个网格为前缘固定端，后13×10个网格为鱼骨柔性结构气动面。
 |
| 图 14 鱼骨柔性翼段气动网格 Fig. 14 Aerodynamic mesh of fishbone flexible wing |
| 图选项 |
假设无穷远处来流速度为10 m/s，鱼骨柔性翼段迎角为5°。本节计算仅考虑作动变形情况下鱼骨柔性翼段外形对应的气动力结果，不考虑气动载荷对结构的影响，相关计算结果如图 15所示。
 |
| 图 15 鱼骨柔性翼段气动面压力系数分布 Fig. 15 Pressure coefficient distribution of aerodynamic surface of fishbone flexible wing |
| 图选项 |
图 15(a)为最大驱动力矩2.0 N·m的加载情况下鱼骨柔性翼段的压力系数分布。由于考虑到三维效应，翼段上不同翼剖面上压力系数分布并不完全相同，取翼段中剖面压力分布，如图 15(b)所示，实线对应鱼骨翼段非线性变形气动力分布，虚线对应线性变形气动力分布。可以看出，随着驱动力矩的增加，前缘固定端压力系数分布变化较小，而鱼骨柔性翼段压力系数显著增加。相同驱动载荷下，线性变形柔性段压力系数要大于非线性变形的结果，并且2个结果之间的差距随着驱动力矩的增加而增大。这表明鱼骨柔性翼段如果用线性有限元分析计算变形，那么在后续气动计算中得到的气动载荷将大于真实变形对应的气动载荷，并且误差会随着变形的增加而增大。
为了验证气动力计算的准确性，将曲面涡格法得到的三维气动力结果换算成二维气动力系数，并与XFOIL软件计算结果进行对比。XFOIL代码是基于势流理论的二维气动力计算软件，已被证明在计算二维翼型时与高精度的CFD计算结果具有可信度，结果如图 16所示。
 |
| 图 16 曲面涡格法与XFOIL计算结果对比 Fig. 16 Comparison of calculation results between non-planar vortex lattice method and XFOIL |
| 图选项 |
从图 16(a)可以看出，不管是非线性变形结果还是线性变形结果，曲面涡格法和XFOIL计算结果符合得很好。随着驱动力矩增加，气动力系数增加。2种方法得到升力系数最大误差为-1.81%。随着驱动力矩增加，2种方法得到升力系数误差增加，当驱动力矩为2.0 N·m时误差最大。这是因为曲面涡格法计算的是薄板中弦线的升力系数，而XFOIL计算的是二维翼型的升力系数，随着驱动力矩的增加，上下翼面弯度变化和中弦线弯度变化差距增加，从而引起升力系数的差别，但其误差仍满足工程分析的精度要求。另外，则看到对于不同气动力程序，同样驱动力矩下线性变形结果计算得到的升力系数要大于非线性结果，在驱动力矩为2.0 N·m时，曲面涡格法得到的非线性结果升力系数相对于线性结果减少4.71%。这是因为线性有限结果在大变形时不考虑x方向变形，相当于增加鱼骨柔性变形部分的相对弦长，而本算例中柔性变形部分比前缘固定端对升力系数影响要大，这一点也可以从图 15(b)翼段压力系数分布中可以看出，引起线性变形后鱼骨翼段升力系数要大于非线性变形。
图 16(b)为不同驱动力矩下，曲面涡格法与XFOIL气动力矩系数计算结果对比。气动力矩取矩点为1/4弦线处，负值代表低头。可以看出，对于同一气动外形，2种气动力计算方法得到的气矩系数符合得很好，相对误差小于1.6%。与图 16(a)结果相似，驱动力矩下线性变形结果计算得到的低头力矩系数要大于非线性结果，在驱动力矩为2.0 N·m时，曲面涡格法得到的非线性结果力矩系数相对于线性结果差异为-7.81%。原因与升力系数差别类似，线性结果不考虑x方向变形，增大了柔性段气动力分布(见图 15(b))，同时也增加了气动力臂，最终引起线性变形的鱼骨柔性翼段力矩系数要大于非线性变形。
结合以上分析，对于鱼骨柔性翼段不同的气动外形曲面涡格法和XFOIL结果误差小于1.81%，变形趋势保持一致。因此，用曲面涡格计算鱼骨柔性翼段气动载荷是可靠的。当驱动力矩较小时，曲面涡格法计算得到鱼骨柔性翼段非线性变形对应的气动力系数和线性变形气动力系数误差不大，但是驱动力矩较大时，两者误差显著增加。其中，反映真实变形情况的非线性结果计算得到的升力系数要小于线性结果，差异最大为-4.71%。非线性结果的力矩系数同样小于线性结果，差异最大为-7.81%。
3.3 气动弹性分析 3.1节和3.2节分别从结构和气动的角度对鱼骨柔性翼段进行了分析。本节将运用基于曲面涡格法和非线性有限元分析的非线性静气动弹性分析方法对鱼骨结构进行静气动弹性分析，并将结果与基于线性有限元分析的传统线性静气动弹性结果进行对比。气动力条件与3.2节一样，假设无穷远处来流速度为10 m/s，鱼骨柔性翼段迎角为5°。加载方式为真实作动模式，参考点为后缘中点，不同驱动力矩下鱼骨柔性翼段收敛过程如图 17所示。可以看出，整个气动结构耦合过程只需要4步到5步迭代就能收敛，计算效率比较高。
 |
| 图 17 鱼骨柔性翼段静气动弹性计算收敛过程 Fig. 17 Convergence process of static aeroelasticity calculation of fishbone flexible wing |
| 图选项 |
图 18(a)为驱动力矩与鱼骨柔性翼段后缘z方向最大位移的关系，可以看出，线性气动弹性分析结果和非线性气动弹性分析结果在后缘z方向上变形差距不大，两者相差小于1.33%，而两者主要区别在于翼段x方向的位移变形，如图 18(b)、(c)所示。可以看出，不管是线性气动弹性分析方法还是非线性分析方法，鱼骨柔性翼段考虑气动载荷作用后翼段都向上偏转。表明在此气动条件下，气动弹性变形对于鱼骨柔性翼段的主动向下弯曲起到了抑制作用。此外，气动载荷引起的鱼骨柔性机翼额外变形随着驱动载荷的增加而增大。这是因为驱动力矩增加，鱼骨柔性翼段柔性段压力系数增加，气动载荷增加，造成载荷作用下附加的位移也增加了。在驱动力矩为2.0 N·m时，气动载荷引起的后缘z方向附加位移最大，线性分析结果为-3.20 mm，是后缘变形的7.11%，非线性分析结果为-2.60 mm，是后缘变形的5.764%，变形抑制程度相当。
 |
| 图 18 不同驱动力矩作用下鱼骨柔性翼段变形 Fig. 18 Deformation of fishbone flexible wing under different driving torque |
| 图选项 |
图 19(a)为考虑气动载荷作用下驱动力矩与鱼骨柔性翼段升力系数的关系。不管是线性分析还是非线性分析，同样驱动力矩作用下，考虑气动载荷作用下鱼骨柔性翼段的升力系数要小于不考虑气动载荷作用时。这是因为此时气动载荷作用下于鱼骨柔性翼段向下变形量减少，引起了升力系数的下降。2种气动弹性分析方法在小变形时结果差距不大，随着翼段主动变形差距显著增加。线性气动弹性分析得到平衡状态升力系数要大于非线性气动弹性分析结果，2种结果最大偏差为8.28%。此外，当驱动力矩为0时，非线性分析得到升力系数损失最多，为-11.78%。随着驱动力矩增加，升力系数下降比例减少，当驱动力矩为2.0 N·m时，升力系数下降-3.71%。这反映出鱼骨柔性翼段升力系数在小变形时更为敏感。
 |
| 图 19 驱动力矩与气动系数关系 Fig. 19 Relationship between aerodynamic coefficients and driving torque |
| 图选项 |
图 19(b)反映了气动载荷对于平衡状态下鱼骨柔性翼段力矩系数的影响。同样因为鱼骨柔性翼段向下变形量减少，同样驱动力矩作用下，考虑气动载荷作用的鱼骨柔性翼段力矩系数要小于不考虑气动载荷作用的情况。非线性气动弹性分析得到，当驱动力矩为2.0 N·m时，力矩系数下降最大，为-8.77%。随着驱动力矩增加，力矩系数下降量增加。线性气动弹性分析得到平衡状态力矩系数也要大于非线性气动弹性分析结果，2种结果最大偏差为6.86%。值得注意的是，传统气动弹性分析中更看重由气动力矩产生的机翼附加迎角的影响。本文计算中鱼骨柔性翼段前缘固定，迎角不随气动力矩变化，也就忽略了气动力矩对于翼段的气动弹性影响，实际的气动弹性效应会更加显著。
通过以上算例可以看出，基于非线性有限元分析和曲面涡格法的大变形静气动弹性分析方法能快速地对鱼骨柔性翼段进行静气动弹性分析，相比传统线性静气动弹性分析方法，气动弹性平衡状态升力系数最多减少8.28%，力矩系数最多减少6.86%。由于非线性有限元分析能反映鱼骨柔性翼段真实变形情况，本文提供的非线性静气动弹性分析方法对于鱼骨柔性翼段分析结果更为准确。鱼骨柔性翼段气动弹性效应较为明显，考虑气动载荷作用下，后缘z方向最大位移减少5.764%，翼剖面升力系数最大下降11.78%，力矩系数最大下降8.77%。
4 结论 本文针对鱼骨柔性翼段在大变形状态下的气动弹性问题，首先运用非线性有限元分析和曲面涡格法分别从结构和气动角度揭示了鱼骨柔性翼段大变形对于气动和结构的影响及线性分析方法的不足，然后运用非线性静气动弹性分析方法从气弹角度对鱼骨柔性翼段进行静气动弹性分析，并与传统线性静气动弹性分析结果进行对比，得到如下结论：
1) 线性有限元分析结果和非线性有限元分析结果在小变形的情况下误差小于0.5%，并且都与文献[7]中的实验果符合很好。但是线性有限元分析方法不考虑鱼骨柔性翼段在x方向位移，在大变形时与非线性分析得到的真实变形结果会产生较大误差。当驱动力矩达到2.0 N·m时，误差达到弦长的2.24%，进而影响气动力计算和气动弹性分析结果。
2) 基于曲面涡格法和几何非线性的大变形静气动弹性分析方法对鱼骨柔性翼段进行静态分析时只需要4步到5步迭代就能收敛，计算效率高。相比线性静气动弹性分析方法，气动弹性平衡状态升力系数最多减少8.28%，力矩系数最多减少6.86%。由于非线性静气动弹性分析方法能反映鱼骨柔性翼段真实变形情况，分析结果更为准确。该方法是三维的静气动弹性分析方法，可以为鱼骨柔性机翼的后续优化设计研究提供理论工具。
3) 鱼骨柔性翼段大变形下气动弹性效应显著，在气动载荷作用下，后缘z方向最大位移减少5.764%，翼剖面升力系数最大下降11.78%，力矩系数最大下降8.77%，对精确的气动外形控制带来不利影响。需要搭建反馈控制系统，对驱动力矩进行补偿，以抵消或者利用气动弹性效应对于鱼骨柔性机翼的影响。
参考文献
| [1] | RAYMER D. Aircraft design: A conceptual approach[M]. Reston: AIAA, 2018. |
| [2] | WOODS B K, FRISWELL M I. Preliminary investigation of a fishbone active camber concept[C]//Proceedings of the ASME 2012 Conference on Smart Materials, Adaptive Structures and Intelligent Systems. New York: ASME, 2012: 1-8. |
| [3] | BARBARINO S, BILGEN O, AJAJ R M, et al. A review of morphing aircraft[J]. Journal of Intelligent Material Systems and Structures, 2011, 22(9): 823-877. DOI:10.1177/1045389X11414084 |
| [4] | WOODS B K, FINCHAM J H, FRISWELL M I. Aerodynamic modelling of the fish bone active camber morphing concept[C]//Proceedings of the RAeS Applied Aerodynamics Conference, 2014: 737-748. |
| [5] | WOODS B K, FRISWELL M I. Structural characterization of the fish bone active camber morphing airfoil: AIAA 2014-1122[R]. Reston: AIAA, 2014. |
| [6] | RIVERO A E, WEAVER P M, COOPER J E, et al. Progress on the design, analysis and experimental testing of a composite fish bone active camber morphing wing[C]//ICAST 2017: 28th International Conference on Adaptive StructUres and Technologies, 2017: 1-11. |
| [7] | WOODS B K, FRISWELL M I. Fluid-structure interaction analysis of the fish bone active camber mechanism: AIAA 2013-1908[R]. Reston: AIAA, 2013. |
| [8] | WOODS B K, FRISWELL M I. Multi-objective geometry optimization of the fish bone active camber morphing airfoil[J]. Journal of Intelligent Material Systems and Structures, 2016, 27(6): 808-819. DOI:10.1177/1045389X15604231 |
| [9] | DRELA M. Integrated simulation model for preliminary aerodynamic, structural, and control-law design of aircraft: AIAA 99-1394[R]. Reston: AIAA, 1999. |
| [10] | MURUA J, PALACIOS R, GRAHAM J M R. Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics[J]. Progress in Aerospace Sciences, 2012, 55: 46-72. |
| [11] | DANG H X, YANG Z, LI Y. Accelerated loosely-coupled CFD/CSD method for nonlinear static aeroelasticity analysis[J]. Aerospace Science and Technology, 2010, 14(4): 250-258. DOI:10.1016/j.ast.2010.01.004 |
| [12] | WERTER N P M, DE BREUKER R, ABDALLA M M. Continuous-time state-space unsteady aerodynamic modeling for efficient loads analysis[J]. AIAA Journal, 2017, 56(3): 905-916. |
| [13] | XIE C C, WANG L B, YANG C, et al. Static aeroelastic analysis of very flexible wings based on non-planar vortex lattice method[J]. Chinese Journal of Aeronautics, 2013, 26(3): 514-521. DOI:10.1016/j.cja.2013.04.048 |
| [14] | 谢长川, 胡锐, 王斐, 等. 大展弦比柔性机翼气动弹性风洞模型设计与试验验证[J]. 工程力学, 2016, 33(11): 249-256. XIE C C, HU R, WANG F, et al. Aeroelastic wind tunnel test model design and experiment on very flexible high-aspect-ratio wings[J]. Engineering Mechanics, 2016, 33(11): 249-256. DOI:10.6052/j.issn.1000-4750.2015.04.0254 (in Chinese) |
| [15] | 刘燚, 杨澜, 谢长川. 基于曲面涡格法的柔性飞机静气动弹性分析[J]. 工程力学, 2018, 35(2): 249-256. LIU Y, YANG L, XIE C C. Study on the static aeroelasticity for flexible aircraft based on non-planar vortex lattice method[J]. Engineering Mechanics, 2018, 35(2): 249-256. (in Chinese) |
| [16] | 王勖成. 有限单元法基本原理和数值方法[M]. 北京: 清华大学出版社, 1997. WANG X C. Basic principle and numerical method of finite element method[M]. Beijing: Tsinghua University Press, 1997. (in Chinese) |
| [17] | 宋天霞, 郭建生, 杨元明. 非线性固体计算力学[M]. 武汉: 华中科技大学出版社, 2002. SONG T X, GUO J S, YANG Y M. Nonlinear solid computational mechanics[M]. Wuhan: Huazhong University of Science and Technology Press, 2002. (in Chinese) |
| [18] | KATZ J, PLOTKIN A. Low-speed aerodynamics-from wing theory to panel methods[M]. Cambridge: Cambridge University Press, 2001: 400-404. |
| [19] | RODDEN W P, JOHNSON E H. MSC/NASTRAN aeroelastic analysis user's guide[Z]. Los Angeles: MSC, 1994. |
| [20] | CHEN P C, LIU D D, KARPEL M. ZAERO user's manual(Version 6.2)[Z]. Scottsdale: ZONA Technology, 2006. |
| [21] | XIE C C, YANG C. Surface splines generalization and large deflection interpolation[J]. Journal of Aircraft, 2007, 44(3): 1024-1026. DOI:10.2514/1.24571 |
