单向及双向流固耦合方法在无人机翼的应用分析

刘基海，易阅城，谢科，宋雨雨，龙怡行，郭东

四川航天系统工程研究所，四川 成都 610100

大展弦比机翼多应用在高空长航时无人侦察飞行器，民用方面常见于气象预测、环境监控以及中继信息等，军事领域则通常应用于空中情报收集、监视和侦察以及空中打击等[1]。

这类高空长航时无人机需具备在临近空间长时间巡航战备能力，故机翼通常被设计为轻质柔性结构。另一方面，无人机在上升到任务高度前（如气流平稳的平流层）需要经过气流紊乱的对流层，飞行器将面临气动弹性稳定性问题，要求结构具有一定的强度和刚度，使飞行器在飞行过程中拥有足够的安全裕度[2]。大展弦比柔性结构在气动载荷作用下产生大变形导致结构几何非线性，使得传统的线性假设不再适用于这类问题的分析。同时，柔性结构会导致颤振临界点进一步降低，气动弹性失稳不可忽视。要分析与解决上述问题，需结合结构的弹性振动方程与流体运动方程求解，这就是流固耦合分析方法[3]。

针对这一复杂的2 种介质之间相互作用问题，国内外学者开展了一系列理论研究与试验验证。柳兆伟[4]针对大展弦比大挠性机翼流固耦合开展数值分析，研究了流体与结构弱耦合计算的关键技术与难点，分别采用单向耦合和双向耦合的方法，对大展弦比大挠度机翼流固耦合问题进行了静态和动态分析。阮远[5]以复合材料潮流能涡轮叶片为研究对象，进行铺层方式及流固耦合研究，分析了不同铺层方式对叶片固有振动特性的影响并进行单向流固耦合分析。研究认为，叶片变形量较小，对流场中流体的运动的干扰不大，只进行了稳态分析。邓艳波[6]针对风力机翼型研究了其气动性能并利用Ansys-CFX 开展流固耦合分析，结果表明，设计的流场模型可以较为准确地反映流场实际流动情况；叶尖参考点位置振幅不断衰减，说明叶轮在设计工况下是气弹稳定的。

马艳峰等[7]采用流固耦合方法对大展弦比机翼开展了非线性颤振特性分析。结果显示：颤振速度随攻角的增大呈现出先略微增大而后减小的趋势，结构非线性给大展弦比机翼颤振带来不利影响，传统的线性分析方法对于大展弦比机翼不适用。赵林建[8]建立了一种求解直升机旋翼气弹耦合特性的流固耦合算法，进行了在纯离心力、无升力旋转和有升力悬停状态下的流固耦合分析。文中未考虑桨叶弹性变形对流场的影响，得到的应力和变形情况与实际桨叶在流场中变形和应力存在些许差异。

渠晓溪[9]设计了一种复合材料结构机翼，并通过ANSYS-Fluent 软件实现了单向流固耦合分析，结果与全金属材料机翼作比较。但在研究中只做了静力学分析，并未进行机翼的固有振动特性和颤振分析。仲继泽等[10]为准确预测机翼颤振边界，发展了一种基于流固单向耦合的能量方法，结果显示，计算得到的颤振边界与实验值吻合，验证了算法的正确性。

本文基于ANSYS-CFX 数值计算方法完成了复合材料机翼大攻角下单向稳态流固耦合、单向瞬态流固耦合与双向流固弱耦合分析，对机翼的固有振动特性以及特定工况下流场和机翼表面压力分布、结构稳态和瞬态响应进行了系统的研究。通过几种分析方式的结果对比，说明上述分析方法的差异性与必要性。

1 流固耦合基本理论

1.1 流体控制方程

流体在流动时遵循质量守恒定律、动量守恒定律及能量守恒定律，混合不同组分的流体系统还将遵循组分守恒定律。一般的可压缩牛顿流体遵循的守恒定律可以通过以下控制方程描述。

质量守恒方程：

动量守恒方程：

式中：t为时间；ff为流体体积力矢量；ρf为流体密度；v为速度矢量；τf为剪切力张量，

其中：p为流体压力；μ为动力黏度；e为速度应力张量，

能量守恒方程的总焓形式为

式中：λ为导热系数，SE为能量源项，htot为流体部分总焓，ρ为流体密度，T为温度。

1.2 固体控制方程

根据牛顿第二定律，可导出固体部分的守恒方程：

式中：ρs为固体密度，σs为柯西应力张量，fs为固体体积力矢量，为固体域当地加速度矢量。

能量守恒方程与流体相似：

式中fT表示温差引起的热变形项：

其中αT表示与温度相关的热膨胀系数。

1.3 流固耦合方程

同理，流固耦合情况下，基本的守恒原则依然适用[11-13]，因此，流固耦合交界面处，流体与固体应力（τ）、位移（d）、热流量（q）、温度（T）等变量应满足相等或守恒条件，即满足以下方程组（下标f 和s 分别表示流体和固体）：

2 结构建模与单向流固耦合

利用Workbench 平台可实现流场与结构数据交换、模型前后处理以及不同模块之间的数据共享[14-15]。图1 为单向流固耦合分析的模块流程图，单向耦合仅将流场数据传递到结构上去，求解结构静变形。虽然变形并不反馈回流场，不涉及到结构与流场的迭代，但在一定程度上能够反映出结构变形的趋势。这种方法简单快捷，计算周期短，对于设计初期的结构性能计算具有重要的意义。

图1 单向流固耦合分析的模块流程

2.1 流体控制方程

基于某飞行状态（飞行攻角-16°，速度20 m/s），在Workbench 平台下，采用Design 模块建立机翼流场几何模型，并设置边界条件：以无穷远处来流风速作为入口边界条件，压力为标准大气压；采用自由出流作为出口边界条件；两侧面为对称边界条件。运用ICEM 模块划分流场网格，流场设置和翼型附近局部网格如图2 所示，利用CFX软件进行流体分析。

图2 流场设置与翼型附近局部网格

本文采用剪切应力传递（shear stress transport，SST）湍流模型建模。SST 模型翼型在绕流、机翼失速等计算上具有更大的优势，其在流动分离和逆压梯度方面能够非常准确地预测光滑自由壁面上的流动情况。同时，SST 模型也考虑了湍流剪切应力的作用，可以精确地模拟出负压力梯度条件下流体的分离量。

计算完成后，通过后处理可以查看流场的流线和压力分布，如图3 所示。图4 为机翼沿展向中间位置的截面压力分布与涡量云图。可以看出：翼型前缘为压力最大值处，上缘压力大于下缘；翼型下表面附近出现大面积负压区，结合涡量云图判断此时机翼已经失速。

图3 机翼表面压力分布与周围流线

图4 机翼截面附近压力分布与涡量云图

为进一步观察翼型压力面气流分离情况，图5给出了翼型表面速度矢量及速度云图。可以看出，在-16°安装攻角下翼型下表面产生了流动分离，出现湍流区。

图5 翼型表面速度矢量及速度云图

2.2 结构网格划分与铺层设计

在进行复合材料铺层前，首先建立机翼的参考平面，即提取出机翼外表面作为参考面进行网格划分。获得网格节点29 369 个，网格单元8 884个，如图6 所示。根据实际加工工艺，由机翼参考面向内侧铺层，使用叠层单元铺层，共铺6 层，角度依次为0°、-45°、90°、0°、45°和 90°，总厚度为0.19 mm。

图6 机翼结构网格划分

复合材料层合板主方向（normals）及参考方向(reference direction)如图7 所示。图8 是根据铺层信息绘制出的极坐标下层合板力学性能。

图7 机翼铺层参考方向

图8 铺层信息和极坐标下层合板力学性能

2.3 静力学计算与预应力模态分析

在CFD 计算及结构建模完成基础上，按照单向定常流固耦合对复合材料机翼开展静力分析。机翼结构载荷来源于流场作用的翼型表面，因此在结构分析过程中，将机翼与流体相接触面设置为流固交界面，导入施加流体对结构压力载荷，如图9 所示，并在机翼根部设置相应的固定面约束。

图9 气动压力导入机翼表面

在气动力作用下，机翼的变形量在叶尖处最大为23.2 mm；应力分布从翼根到翼尖处逐渐减小，最大应力为111.85 MPa，位置为靠近翼根上表面处。机翼前后缘部分应力较小且分布均匀。应变分布情况和应变类似，具体的变形、应力和应变云图见图10。本研究中，为使模型在双向流固耦合计算中不会因为大变形而出现负网格等报错，设置的动压较小，故变形量也不大。

图10 稳态分析结果

对处于飞行状态的机翼进行模态分析时，需要考虑预应力效应。流固耦合作用引起的不均匀气动力会造成机翼产生弯曲、扭转等变形，其模态与静止环境下的值必然会有差别，因此对预应力条件下的机翼进行模态分析是很有必要的。

为分析翼载荷对机翼固有振动特性的影响，计算了预应力下频率与固有频率（未加载）并作对比分析。前六阶振型如图11，对应模态频率见表1。

表1 预应力工况下前六阶频率与固有频率对比

图11 机翼前六阶振型

根据计算结果显示，由于存在几何非线性，机翼在预应力作用下结构刚度特性发生了微小变化，固有模态频率普遍降低，且对扭转频率的影响较面内弯曲频率更为显著。在本文中，可认为压力对机翼的作用几乎不会对其振动频率产生影响。

2.4 瞬态动力学分析

将求解的流场非定常空气动力作用到结构上，并设置结构动态分析时间步长与流场非稳态分析一致，时间步长为0.05 s，计算到10 s。图12为不同时刻机翼的变形情况，可以看出，机翼结构随时间出现周期性的弯曲振荡并逐渐收敛至稳态，瞬态分析收敛后的结果与稳态结果存在差异。

图12 单向瞬态分析中不同时刻的机翼变形

为研究机翼随时间的振动过程，绘制了翼尖结构响应曲线和对应的功率谱密度曲线，同时给出最大应力应变的时间响应历程，如图13 所示。可以看出，翼尖位移最大值为35.7 mm，振动持续约2 s 收敛至28.4 mm，振动频率7.2 Hz。瞬态稳态结果对比见表2。

图13 机翼翼尖单向瞬态分析

表2 单向稳态与瞬态分析结果对比

通过对比单向瞬态与稳态结果发现，稳态结果相对于瞬态分析结果更加保守，瞬态分析中无论是最大变形与应力应变还是收敛后的结果都较稳态分析大。在力学强度校核中，进行单向瞬态分析是十分必要的。

3 双向流固耦合

3.1 CSD 与CFD 计算模型

与流场单向作用于结构的单向耦合分析不同，双向流固耦合是将结构动态分析与流场非稳态分析相互耦合，结构的变形与流场的相互作用通过流固耦合面进行数据传递。本文通过Workbench实现双向流固耦合分析，算法流程如图14 所示。

图14 双向流固耦合求解流程

机翼根部端面设定为固定约束，其余表面设定成为流固耦合面，如图15。CFD 计算模型、边界条件以及工况等与单向分析一致，合理考虑计算量与计算精度，取时间步长0.05 s，总分析时间10 s。计算方式为顺序求解，即弱耦合。

提交计算后在监控界面设置流固耦合面升力变化为监控曲线，可方便观测计算收敛情况。图16为单向与双向分析中流固耦合升力变化。

图16 计算过程中流固耦合面升力变化情况

由图16 可以看出，双向耦合分析时耦合面升力最大值明显低于单向分析。

3.2 计算结果分析

计算得到了不同时刻机翼变形与翼尖截面处流场流线分布，如图17 所示。绘制了翼尖结构响应曲线和对应的功率谱密度曲线如图18 所示。

图17 不同时刻机翼变形情况和翼尖截面处流场流线分布

图18 机翼翼尖双向流固耦合分析

由图18 可以看出，双向流固耦合中机翼翼尖最大位移为44.2 mm，振动持续约1.5 s 收敛至33.2 mm，振动频率7.2 Hz。图19、图20 给出了机翼应力应变分布情况。

图19 双向流固耦合分析中机翼应力分布情况

图20 双向流固耦合分析中机翼应变分布情况

为与单向瞬态流固耦合分析结果对比，绘制了单向瞬态与双向流固耦合分析中翼尖响应曲线与频响曲线对比图，如图21 所示。

图21 机翼翼尖单向瞬态与双向流固耦合分析对比

由图21 可以看出，无论翼尖最大位移还是振动收敛后平衡位置，双向流固耦合分析结果都比单向瞬态分析结果大。双向耦合震荡收敛得更快，但振动频率几乎一致。表3 和表4 给出位移、应力应变的具体比较结果。

表3 双向与单向瞬态流固耦合分析变形对比

表4 双向与单向瞬态流固耦合分析应力应变对比

通过表4 可以得出，双向流固耦合分析得到的最大应力应变比单向瞬态分析结果小。图22给出了双向与单向瞬态分析中应力分布情况对比，通过对比发现，单向瞬态分析应力最大值发生在靠近翼根的蒙皮上表面，而双向分析应力最大值发生在翼根附近。

图22 应力分布情况对比

由于双向流固耦合分析考虑双向数据传递与互相作用，更加接近真实情况。为获得更近似的数值解，开展双向流固耦合分析是必要的。但双向耦合分析计算周期长、工作量大、鲁棒性也非常差。尤其在计算过程中，由于插值与网格变形容易出现错误而导致计算终止。

同时也可发现，单向流固耦合的应力数据偏大，也就是说，用单向方式求解获得的结果来进行校核是偏于安全的。采用双向流固耦合分析的方法可以对机翼的气弹稳定性进行研究，通过观测振动位移响应发散与否来判断是否发生气弹失稳现象。

4 结论

本文建立了某无人机机翼大攻角下的流场分析模型与复合材料结构有限元模型，计算了非定常气动力与压力分布，分别完成了单向稳态流固耦合分析、预应力下模态分析、单向瞬态流固耦合分析与双向流固耦合分析。根据结构响应分析得到结果如下：

1）由于存在几何非线性，机翼在预应力作用下，结构刚度特性发生了微小变化，固有模态频率普遍降低，且对扭转频率的影响较面内弯曲频率更为显著。

2）单向流固耦合分析中，定常稳态分析结果相对于瞬态分析结果更加保守，瞬态分析中无论是最大变形与应力应变还是收敛后的结果都较稳态分析大。

3）采用双向流固耦合分析的方法可以对机翼的气动弹性稳定性进行研究。单向流固耦合的应力数据较双向耦合分析结果偏大，用单向耦合方法求解获得的结果来校核是偏于安全的。尽管双向耦合计算周期长、工作量大，但为获得更近似的数值解，就有必要开展双向流固耦合分析。