基于时域气声伴随算法的直升机转子叶片形状优化分析

何 磊，魏海军

(1.安徽农业大学 经济技术学院，安徽 合肥 230011；2.上海海事大学 商船学院，上海 201306)

直升机机动灵活，可以在低空飞行，战略打击，空中救援，航拍侦查，被广泛应用于军事领域以及民用航空。旋翼作为直升机的核心动力部件，研究提升其性能对直升机的动力性能、飞行性能、飞行稳定性有着重要作用。

伴随方程已被广泛用于计算流体力学的气动外形优化问题，是一个非常强大的工具，该方法将流场控制方程作为等式约束引入目标函数，将物体边界作为控制函数，创新性地把优化设计问题转化为控制问题，这为复杂气动外形优化设计提供了可能，因为它允许以基本参数输入，来计算一组给定输入的目标函数的灵敏度导数。单次梯度信息求解计算量只相当于两倍的流场计算量，与设计变量的数目无关。基于伴随的稳态问题形状优化被广泛应用于复杂气动外形优化设计的二维和三维问题[1]。Mavriplis使用非定常伴随方法，开展了空气动力学中三维问题的应用研究[2]。Nielsen等人采用离散型伴随方法，使用ANSYS CFX开展直升机旋翼的外形气动优化设计的研究[3]。李亚等研究了在噪声计算中采用了FW-H方程、结合旋转偶极子辐射、涡声方程三种方法，实现了大涡模拟预报螺旋桨辐射噪声的三种声学方法[4]。南京航空航天大学的吴琪建立了一套适合于先进旋翼气动外形综合优化设计分析的黏性伴随方法和代理模型方法相结合的混合设计方法，实现了基于黏性伴随方法的旋翼先进气动外形优化设计分析[5]。伴随方法在多学科设计优化问题的应用主要集中在固定翼和旋转翼飞机的流固耦合上[6]。本文将三维多学科伴随方法推广到旋翼机的降噪问题。由于对旋翼机噪声的要求越来越严格，旋翼气动声学在直升机的设计过程中起着主导作用。通常，转子噪声已通过选择转子叶片配置、翼型、平面形状和叶尖形状来解决。虽然这些方法主要依靠风洞实验和飞机设计工程师的经验，但很少有人尝试在降噪背景下设计优化翼型的形状。在之前开展的项目研究中，我们推导了二维翼型的耦合伴随气动声学问题，并证明了非定常气动声学优化的可行性，这里我们研究了一种三维非定常伴随方法，使用基于梯度优化用于直升机转子气动声学应用的方法，以便通过应用最佳形状修正，在不产生显著性能损失的情况下使得直升机转子噪声最小化。

1 时域气声伴随算法的计算模型的构建

1.1 基于CFD方法的直升机转子叶片流场计算模型

关于流体的数值模拟分析公式，本文中使用的是基本三维非结构网格雷诺平均Navier-Stokes流体数值模拟方法。ANSYS CFX已被广泛用于验证稳态和时变流，并包含离散切线和伴随灵敏度能力，这已被证明用于稳态和时变流问题的优化。因此，本文将仅对这些公式进行简要描述，其他细节可在以前的参考文献中获得[7]。流体的数值模拟基于Navier-Stokes方程的保守形式，可写成：

(1)

对于移动网格问题，这些问题以任意拉格朗日-欧拉形式表示为：

(2)

使用以顶点为中心的中值双控制体积公式，该公式为二阶精度，其中闭合控制体积周围的无粘通量积分S离散为：

(3)

(4)

以简化的形式展示了不同时间水平面上U和x的函数关系，如下所示：

Rn(Un,Un-1,Un-2,xn,xn-1,xn-2)=0

(5)

在每个时间步长n，隐式残差相对于未知解向量Un线性化，并使用牛顿法求解如下：

(6)

其中，向量U表示所有时间步的流量值，并且该方程中的每个(块)行对应于方程(5)中给出的特定时间步的解。方程(6)表示所有时间步的联立解，在实践中，由于每个新的时间步取决于前两个时间水平，因此通过使用前向块替换(即时间上的前向积分)的牛顿法求解。

表1 直升机叶片转子参数

1.2 推导气动灵敏度分析公式

基本敏感性分析的实施遵循参考文献中开发的策略[7]。考虑一个任意的目标函数L，该函数使用非定常流动解集U和非定常网格解集x进行计算，表达式为：

L=L(U,x)

(7)

假设状态变量(即U，x)依赖于一些输入设计参数D，目标函数L对设计输入的总灵敏度可以表示为设计输入的状态灵敏度向量和状态变量的目标灵敏度向量之间的内积，如下所示：

(8)

如前所述的非线性流动残差算子和线性弹性网格残差算子提供了约束，这些约束可以在整个空间和时间域上以一般形式表示为：

(9)

当对设计输入进行线性化时，推导得到：

(10)

这些构成了正向灵敏度或切线灵敏度方程。然后，网格和流量灵敏度矢量可以代入方程(8)，以获得目标相对于设计变量D的完全灵敏度。

对于每个设计参数D，正向灵敏度方法需要方程(10)的新解。另一方面，伴随方法可以以近似独立于设计变量数目的代价获得任意数目的设计输入D的灵敏度。伴随问题可通过将方程(10)预乘大耦合矩阵的逆，并将灵敏度的结果表达式代入方程(8)，并将伴随变量定义为系统的解来获得：

(11)

其中ΛU和Λx分别是流动和网格伴随变量。最终目标灵敏度可通过以下方式获得：

(12)

回顾方程(11)和(12)适用于整个时域，方程(13)表示正向时间积分，而方程(12)中的回代过程导致时间上的反向积分，从最后一个物理时间步，然后进入初始时间步。

1.3 前飞旋翼桨距公式

由于本文的目标是优化配平前飞状态的旋翼，因此必须将旋翼桨叶桨距角作为控制输入，以便旋翼能够获得所需的推力和侧向力矩[8]。俯仰角可以表示为平均俯仰角和几个谐波分量的组合：

(13)

其中ψ是转子方位角。在本文中，我们只使用一个谐波，因此i=1，设计变量D=[θ0，θci，θsi]的选择应使转子得到调整，其中θ0是总参数，θci，θsi是两个周期螺距参数。为了将旋翼桨距驱动纳入计算流体动力学分析中，需要一个叶片表面桨距驱动方程(剩余形式)：

Sθ=xSθ+[Tθ(θ(D))]x=0

(14)

和网格运动方程

G′(xθ,xsθ)=0

(15)

其中，xsθ表示俯仰后的叶片表面网格坐标，xθ表示俯仰和方位旋转后网格运动的内部网格坐标，xp表示俯仰和方位旋转后的表面和内部网格坐标，而U表示的是非定常流动解。

2 转子叶片噪声优化分析

在本文中建立的模型是前飞状态的 NACA0012 四叶旋翼，其自由流马赫数为0.095，翼尖马赫数为0.638，雷诺数为2×106。在此条件下，非线性流动效应可以忽略不计，因此可以使用机身和机翼组成的不可渗透的集成表面。声学表面与ANSYS CFX转子表面网格重合，表面上的每个节点都是一个声源，通过 FW-H 集成过程在观察者位置产生声压[9-10]。将声学模块的线性化与流动求解器的线性化相结合，可以计算声学目标函数相对于在基于梯度的形状优化过程中使用的设计变量的完整向量的伴随灵敏度。

2.1 声学分析公式：FW-H方程

FW-H方程可以用微分形式表示为：

(16)

其中，

其中，公式(17)表示流体气动弹性非稳定力关系函数，公式(18)表示流体质量变化函数。f表示所定义的求解积分面，是时间和位移的函数。ρ表示空气总密度，p表示总压力。ui为流体速度，而vi为转子叶面速度，c0为自由流体产生的声速。ni是声学表面的单位法线。然后通过方程(22)给出积分解，忽略四阶无穷小量。

其中y是观察者位置，t是观察者时间，r是源x和观察者y之间的距离，M是表面马赫数。方程(22)是基于时域气声伴随算法的FW-H方程，p′T(y,t)代表由厚度产生的声压在叶面面壁边界上的积分，p′L(y,t)代表由载荷产生的声压在叶面面壁边界上的积分，其中含1/r的项为远场噪声，含1/r2的项为近场噪声。

表2 声学目标函数相对于直升机转子叶片的观察者位置(R是转子叶片半径)

2.2 基于时域气声伴随算法的FW-H灵敏度公式

计算耦合气动声学问题的灵敏度需要将声学模型相对于设计变量进行线性化[11]。如前一节所述FW-H积分方程的灵敏度公式与其他参考文献中详述的CFD数值模拟曲线的线性化非常相似。首先，声码的正向线性化是通过离散FW-H积分方程的精确微分设计的。然后转置正向线性化，并从相反的方向推导求解得到关于设计变量完整向量的伴随灵敏度公式。针对声码的复杂微分形式使用正向灵敏度进行验证。最后，通过伴随灵敏度公式来验证切线的双重性问题。

声学目标函数的切线灵敏度公式，是计算关于单个设计变量的目标函数的正向线性化灵敏度。本文将观察者位置和时间、声压关系表示为：

p′(y,t,D)=F(U(D),x(D))

(23)

其中D是设计变量的向量表示，F(U(D)，x(D))表示以数值方式评估方程 (23) 所需的所有离散运算。在这个表达式中，U(D)和x(D)是时间积分过程中每个时间步的流动和网格解。声压切线灵敏度时间方程可表示为：

(24)

所以推导得到定义的声学目标函数的切线线性化方程表示为：

(25)

伴随灵敏度公式，伴随灵敏度可以通过转换正切灵敏度公式得到。转换方程(25)得到

(26)

3 转子叶片算例验证

相应的转子转速为1100转每分钟(旋翼轴向进给比μ=0.15)。ANSYS CFX数值模拟由2.32×106个节点所组成的网格，进行了两次旋翼叶片旋转试验，网格由2度时间步长、每个时间步长6个流固耦合迭代和位于距转子轮毂两个半径的静止平面内观测器组成，方位角ψ= 180度。进行了优化设计计算。如图1所示的结果，在这两种情况下，对比发现本文所建立的FW-H模型实现与实验值之间一致性非常好[11]。表明本文计算方法和建模可以有效地预测前飞状态的旋翼转子叶片的噪声。

图1 某型直升机叶片转子结构图

基于时域气声伴随算法的直升机转子叶片FW-H积分方程的灵敏度优化收敛结果如图2所示，经过12次非线性迭代和20次设计循环之后，优化问题的可行性和最优性伴随灵敏度提高了3个数量级。对直升机四叶NACA0012旋翼进行配平之后，其推力和侧向力矩关于时间历程的数值变化进一步表明了本文计算方法和建模可以有效地预测前飞状态的旋翼转子叶片的噪声。图3是直升机转子旋翼叶片的厚度噪声、载荷和总声压随时间的变化曲线，厚度噪声仅受到主要作用于加载压力的优化的轻微影响，从而得到最小化的总声压。噪音最小化是通过改变厚的翼型实现的，特别是在内侧部分。由于转子扭矩不包括在优化范围内，因此尽管代表了当前数值模拟的优化方法可实现最大的降噪量，但这里实现的降噪是以牺牲直升机旋翼转子性能为代价的，如图4所示。最后得到了气动声学优化直升机转子旋翼叶片旋翼的优化翼型形状，如图5所示。

图2 直升机转子旋翼计算噪声声压与实验值对比曲线

图3 直升机转子旋翼配平非线性迭代函数优化的可行解和最优解对比曲线

图4 直升机转子旋翼叶片的厚度噪声、载荷和总声压随时间的变化曲线

图5 气动声学优化直升机转子旋翼叶片旋翼的优化翼型形状

4 结论

本文建立了一种直升机气动声学时变问题的灵敏度分析公式。采用ANSYS CFX数值模拟方法与FW-H声学模块耦合的混合气动声学方法，将直升机转子噪声信号传播给远场观测者。得到以下结论：

(1)推导了与气动声学分析问题类似的切向灵敏度公式和伴随灵敏度公式。首次验证了全耦合气动声学分析公式能有效地预测四叶NACA0012旋翼在前飞时的噪声。在验证成功的基础上，采用伴随公式对旋翼噪声信号进行优化。

(2)前飞状态的直升机旋翼旋转总噪声、厚度噪声、载荷和总声压优化后数值有明显降低的趋势，厚度噪声仅受到主载荷压力的优化轻微影响，从而得到最小化的总声压。噪音最小化是通过改变厚的翼型实现的，特别是在内侧部分。

(3)本方法对其他汽车噪声降噪处理、航空固定翼机型机翼气动声学优化都有较强的参考意义。