高保真CFD/CSD方法在旋翼叶片优化中的应用

《Aerospace Science and Technology》：High Fidelity CFD/CSD Method for Rotor Blade Optimization

　　A. Woodgate | George N. Barakos
CFD实验室，格拉斯哥大学工程学院，格拉斯哥，G12 8QQ，苏格兰，英国

**摘要**
本文展示了一个使用高保真方法设计旋翼叶片的框架。首先对一个初始的矩形叶片进行了优化，优化考虑了叶片表面参数化的各种变化，包括扭转、弦长、后掠角和扭转角的变化。特别关注了在期望飞行条件下将气动弹性与伴随优化和修剪方法相结合的问题。本文得出了关于新设计效果以及该框架在实现设计方面的效率的结论。所获得的设计具有平衡悬停和前飞性能的扭转分布，并且需要一种叶片尖端设计来减轻前飞过程中的不利影响。显著的特点包括前后尖端的后掠角以及弦长分布的轻微减小。

**1. 引言**
改进的直升机旋翼叶片设计旨在减少燃油消耗、扩大飞行范围，并最小化噪音和振动。要实现这些目标，需要仔细考虑叶片形状，同时考虑实际约束条件。近年来，有证据表明在实际叶片设计中至少使用了一些优化技术。诸如BERP设计[1]、Blue-Edge叶片[2]或Advanced Chinook Rotor Blade (ACRB)[3]等先进叶片形状表明，旋翼设计仍在不断进步，在为每种直升机提供最优叶片形状之前还有许多工作要做。观察这些工作可以发现，BERP叶片的开发使用了缩比旋翼进行实验，而计算流体动力学（CFD）仅在开发后期才被应用。BlueEdge的最终设计是在多个研究项目和广泛实验的支持下经过多个阶段迭代后完成的。ACRB叶片是较新的发展，其设计更多依赖于数值模拟和CFD。然而，在初始飞行测试阶段，ACRB遇到了两个问题：振动问题和后旋翼叶片的失速现象。

展望未来，旋翼设计的方向是将高保真模型与优化方法相结合。相关文献主要分为两类研究方向：无梯度优化方法和基于梯度的优化方法。总体而言，无梯度方法需要更多的CFD计算才能得到最优解；相比之下，基于梯度的方法需要更多的开发工作，如果不加以注意，可能会导致接近初始叶片设计的局部最优解。无梯度方法之所以普遍使用，是因为它们相对容易使用现有的计算机代码进行开发和实现。例如，任何能够计算旋翼流动的CFD求解器都可以与某种形式的元模型结合使用，以“拟合”预测的气动特性，并将其输入到网格生成工具和优化方法中。最流行的元模型包括克里金插值、有理基函数、神经网络和简单的查找表。如果元模型中包含了大量不同设计的CFD结果，还可以计算出气动特性相对于设计参数的梯度，从而形成基于梯度的方法。随着设计参数数量的增加，所需的CFD计算量也会增加，这使得这些方法在实际应用中变得不太实用。

另一方面，基于梯度的方法需要的CFD计算较少，并且可能更好地适应设计变量的数量。与参数化研究相比，基于梯度的优化方法允许探索更全面的设计空间，更好地帮助设计师选择叶片的最佳形状。它们还允许探索更多新颖的叶片设计，并提供关于设计特性对最终气动特性的影响的信息。然而，这些先进方法在解决方案的保真度、准确性、计算成本以及多学科分析（包括气动、气动弹性、结构和振动）的耦合方面依赖于算法的进步。捕捉旋翼设计的多学科方面对于工业应用至关重要，因为它能确保获得现实的设计结果。最后，选择优化方法和设计参数也将驱动设计变更，从而影响最终结果。未来的旋翼设计可能会基于工程判断自适应地进行，而不需要预设参数。

本文的主要贡献在于观察到在叶片优化研究中并不总是考虑叶片的弹性。关于旋翼叶片设计的综合文献可以在[4]中找到。显然，虽然使用CFD和伴随方法优化叶片的气动形状很重要，但大多数文献中的工作都忽略了问题的结构和气动弹性方面。本文解决了这一问题，将叶片弹性的影响加入到形状优化所需的计算中。此外，随着优化的进行，结构模型的属性也会更新，以便在优化循环中反映截面积、扭转等的变化对叶片载荷和性能的影响。

**2. 方法论**
本节描述了本研究中使用的CFD-CSD方法论。这包括HMB3求解器、谐波平衡方法、伴随方法、优化框架以及将HMB3与NASTRAN接口的中间件工具。下一节将描述在本项目中实施的发展，包括叶片参数化和设计框架、将谐波平衡方法与非稳态气动计算的模态方程相结合，以及中间件和模态方法与HMB3的伴随求解器的功能增强。

**2.1. 计算流体动力学方法**
Helicopter Multi-Block 3 (HMB3)工具集是本研究中所有模拟所使用的内部开发工具。该计算流体动力学工具是在结构化多块网格[5]上使用的有限体积求解器，采用重叠网格方法[6]。HMB3使用任意拉格朗日欧拉（ALE）公式求解时变域中的非稳态雷诺平均纳维-斯托克斯（URANS）方程，包括移动边界。本工作中应用的两种不同求解公式如下：
(1) \(dW/dt + R(W) = S\)
(2) \(dW_b/dt + \omega D_W_b + R_hb(W) = 0\)
其中，\(W\)是守恒变量向量，\(R\)是残差向量，\(S\)是源项，\(t\)是实时步长，\(\omega\)是归一化频率，\(D\)是傅里叶配置导数算子，\(HB\)（和\(hb\)）指的是谐波平衡方法。

对于这两种方法，为了评估对流通量，使用了Osher[7]近似黎曼求解器，而粘性项则使用二阶中心差分空间方案离散化。Leer[8]开发的MUSCL方法在大梯度区域提供了高空间精度，并采用了Albada限制器[9]的替代形式。Jameson采用了隐式双时间步进方法。线性化方程组使用广义共轭梯度方法求解，并以BILU分解作为预处理器。HMB3求解器提供了一到四方程的湍流模型，本文的计算使用了Menter[10]的1994 k-ω SST模型。

HMB3中的悬停离地效应（OGE）模拟使用了弗劳德边界条件。旋翼网格包含\(2\pi/N_b\)个旋翼域，其中\(N_b\)是叶片总数。因此，该网格只包含一个叶片，这对于使用稳态公式进行悬停计算是足够的。在离散化的纳维-斯托克斯方程中添加了一个源项，以考虑旋翼中的旋转速度分量。

**2.2. 谐波平衡方法**
为了降低非稳态伴随设计相关的计算成本，对于前飞情况使用了谐波平衡方法[11][12][13][14]。该方法通过假设时间的周期性将解决方案视为一个大的稳态问题，并增加了一个非稳态源项。解决方案的稳态特性使其适用于伴随优化研究，因为可以使用稳态伴随求解器（经过某些修改）。

谐波平衡方法将流动解和残差向量表示为截断的傅里叶级数，假设时间具有频率\(\omega\)的周期性。流动由\(NH\)个谐波平衡模式表示，并分为\(NT = 2NH + 1\)个子区间，这些子区间通过傅里叶配置导数算子\(D\)耦合。与时间推进模拟相比，这显著降低了计算成本，因为流动场可以作为一个大的（\(NT \times NT\)）稳态问题来解决。

从半离散形式的常微分方程系统开始：
(3) \(I(t) = V dW(t)/dt + R(t) = 0\)
（其中\(V\)是单元体积）。假设\(W(t)\)和\(R(t)\)是时间上的周期函数，并且是频率\(\omega\)的函数，它们可以表示为傅里叶级数：
(4) \(W(t) = \sum_{k=-\infty}^{\infty} W^k e^{ik\omega t}\)
(5) \(R(t) = \sum_{k=-\infty}^{\infty} R^k e^{ik\omega t}\)
这个级数可以截断为指定数量的谐波\(NH\)：
(6) \(W(t) = \sum_{k=-N_H}^{N_H} H^k e^{ik\omega t}\)
(7) \(R(t) = \sum_{k=-N_H}^{N_H} H^k e^{ik\omega t}\)
方程(3)也可以使用傅里叶级数展开来截断：
(8) \(I(t) = \sum_{k=-N_H}^{N_H} H^k e^{ik\omega t}\)
将离散傅里叶变换代入半离散形式的方程(3)，并将状态变量的时间导数移到级数求和内部。利用傅里叶项的正交性，可以得到每个波数\(k\)的方程：
(9) \(\omega_k V W^k + R^k = 0\)
需要注意的是，由于\(R(W(t))\)是\(W(t)\)的非线性函数，因此每个系数\(R^k\)都依赖于所有系数\(W^k\)。因此，方程(9)表示一个必须迭代求解的非线性方程组：
(10) \(I^k = i\omega_k V W^k + R^k = 0\)
该系统有\(NT = 2NH + 1\)个方程，对于\(NH\)个谐波，可以写为：
(11) \(\omega A W^k + R^k = 0\)
其中矩阵\(A\)的大小为\(NT \times NT\)。该系统使用伪谱方法求解，通过将方程(10)转换回时间域并将解分成周期\(T = 2\pi/\omega\)内的\(NT\)个等间隔子区间：
(12) \(W_b = (W(t_0 + \Delta t) W(t_0 + 2\Delta t) \cdots W(t_0 + T)\)
\(R_b = (R(t_0 + \Delta t) R(t_0 + 2\Delta t) \cdots R(t_0 + T)\)
其中\(\Delta t = 2\pi/(NT\omega)\)。

通过使用变换矩阵\(E\)，使得\(W^k = E W_b\)和\(R^k = E R_b\)，我们可以写出：
(13) \(0 = A E W_b + E R_b = \omega E^{-1} A E W_b + E^{-1} R_b = \omega D W_b + R_b\)
其中\(D = E^{-1} A\)，定义为：
(14) \(D_{ij} = \frac{2NT}{\sum_{k=1}^{N_H} k \sin(2\pi k(j - i)/NT}\)
然后对方程(13)应用伪时间推进：
(15) \(dW_b/dt + \omega D W_b + R_b = 0\)
该方程使用隐式方法求解，对源项进行隐式处理（以提高CFL数）：
(16) \(\omega D W_b^{n+1} = \omega D W_b^n + \omega D(\Delta W_b)\)
完整的谐波平衡方程变为：
(17) \(W_b^{n+1} - W_b^n \Delta t^* = -[\omega D W_b^{n+1} + R_b(W_b^{n+1})\)
这种非稳态源项的实现将所有变量在所有\(NT\)个解快照中耦合起来。在线性系统中使用基于残差函数一阶离散化的近似雅可比矩阵来改善系统的条件数：
(18) \(J = \left[ \frac{\partial R}{\partial W|t_0 + \Delta t} \omega D_1,2 \cdots \omega D_{1,NT} \frac{\partial R}{\partial W|t_0 + 2\Delta t} \cdots \omega D_{NT,1} \frac{\partial R}{\partial W|t_0 + T} \right]\)
然后使用Krylov子空间方法和BILU分解求解该系统：
(19) \(P[V\Delta t^* + J] \Delta W_b = -P R_hb\)
谐波平衡方法的主要缺点是需要求解2\(NT + 1\)个耦合系统方程，对于大量的谐波平衡模式来说，这可能会占用大量内存。由于需要在雅可比矩阵和预处理器中存储额外的非对角线块，隐式方法的内存占用会更快地增加。为了减少这种开销并获得与谐波数量线性相关的存储需求，不显式存储雅可比矩阵的非对角线项，而是在矩阵-向量乘法过程中将其加到右侧，并在方程(18)中将\(\omega\)设置为0来构建预处理器。这种对预处理器的更改对求解线性系统的Krylov迭代次数影响很小，但也有助于减少线性系统的整体CPU时间，因为预处理器向量的乘法操作减少了[15]。

这种实现与其他常用的频域方法不同。McMullen等人的非线性频域方法[16]直接在频域求解流动方程，但残差是在时域中制定的。这需要在迭代过程中进行快速傅里叶变换和逆快速傅里叶变换操作。本研究采用了Hall[17]的方法，在谐波平衡方程中使用离散傅里叶变换矩阵（在伪谱导数算子内部），允许直接在时域中获得傅里叶级数系数和残差的解。Gopinath和Jameson的时间-谱方法[18]与我们的方法类似，但在时间导数项的实现上略有不同。

**2.3. 伴随方法**
伴随方法提供了一种高效的方法来生成设计变量的梯度，这些梯度与基于梯度的优化器一起使用。与有限差分法相比，后者中梯度的计算依赖于变量的数量，而伴随方法可以在获得基流解的情况下生成所有设计变量的梯度。基本思想是将目标函数表示为流量和设计变量的函数。我们假设流量变量满足纳维-斯托克斯方程，因此有：
(20)R[W(α),α]=0
正式地，对I关于α求导得到：
(21)dIdα=?I?α+?I?W?W?α
通过引入伴随变量λ作为以下线性系统的解：
(22)(?R?W)Tλ=?(?I?W)T
方程(21)可以重写为：
(23)dIdα=?I?α+λT?R?α
对偶灵敏度问题(22)-(23)的计算成本随输出数量的增长而增加，因为(22)的右侧依赖于I，但它与输入参数无关。当（输出）成本函数的数量较少而（输入）设计变量的数量较多时，伴随形式的灵敏度方程特别有效。对于相反的情况，即（输出）成本函数的数量较多时，可以使用 tangent 求解器。因此，对于转子设计优化问题，伴随方法更有意义，并且消除了对设计变量数量的限制（受优化器算法的约束）。方程22中显示的线性系统使用基于嵌套广义最小残差（GMRES）求解器的Krylov子空间求解器[19]来解决。GMRES是一种用于求解非对称线性系统的稳健方法，它基于使用Arnoldi方法[20]生成的Krylov子空间上的方程残差最小化。一个众所周知的理论结果是，GMRES保证最多在n步内收敛，其中n是系统矩阵的维度[21]。该方法的缺点是内存需求较大，因为它需要存储Krylov子空间的所有基向量，而每个Arnoldi迭代后基向量的大小会增加一个。为了限制内存需求同时确保良好的收敛性，采用了Deflated Restarting（DR）技术（表示为GMRES-DR），即在预定义的m次Arnoldi迭代后丢弃Krylov子空间，并从最后一个解开始重新开始。为了提高收敛性，在重新开始过程中，从Krylov子空间中提取与最小特征值相关的k个近似特征向量，并在下一个循环中使用它们来扩展子空间。GMRES求解器以嵌套方式实现，其中GMRES-DR求解器作为（外部）GMRES-DR循环的预处理器。为此，需要一个外部求解器的灵活变体[22]，尽管这会导致更高的内存需求，因为在每次Arnoldi迭代时需要存储两个向量。内部的GMRES-DR求解器反转预处理矩阵J^，并反过来使用解耦的ILU(0)预处理器。整个方案可以表示为FGMRES-DR(m, k)-GMRES-DR(mp, kp)。这种嵌套Krylov子空间求解器被证明能够解决复杂的伴随问题，这对于刚性谐波平衡伴随解是必需的。
伴随方法是通过源代码微分在HMB3中开发的，使用了TAPENADE工具[23]。求解器的各个函数被微分并用于计算灵敏度方程23所需的偏导数。嵌套Krylov子空间求解器的实现不需要雅可比矩阵J，而只需要矩阵-向量积JTv。因此，J从未被显式存储，JTv是通过残差函数的自动微分代码计算得到的。关于计算空气动力灵敏度所需的偏导数和矩阵-向量积JTv的更多细节在[24]中给出。对于设计变量的灵敏度，项?R?α的计算如下：
(24)?R?α=?R?N?N?Xv?Xv?Xs?α
这表示相对于设计变量的残差向量灵敏度的计算，它依赖于每个设计变量的表面网格灵敏度?Xs?α（从叶片参数化工具获得）、相对于表面网格的体积网格灵敏度?Xv?Xs、相对于体积网格的网格度量灵敏度?N?Xv，以及最终相对于网格度量的残差灵敏度?R?N。总之，所采用的伴随方法是一种非常灵活的方式来计算梯度，无需分析伴随方法的复杂性和限制，也不需要为大量设计变量进行成本高昂且难以获得的导数数值评估。当前伴随方法的验证和确认已在早期出版物中详细说明，例如[24]。

2.4. SLSQP优化器
优化问题可以如下表述：
(25){Findminαi,i∈{1,?,n}I(α1,?,αn)subjecttoαi,min≤αi≤αi,max,i∈{1,?,n}gj(α1,?,αn)≤0,j∈{1,?,m}hk(α1,?,αn)=0,k∈{1,?,p}
这旨在在设计空间αi,min到αi,max的范围内最小化或最大化目标函数I，同时满足一系列等式（h）和不等式（g）约束。当前的优化问题使用序列最小二乘二次规划（SLSQP）算法[25] [26]来解决，其实现基于NLopt软件库[27]。优化器代码是HMB3求解器套件中的一个独立工具，因此可以轻松使用其他优化算法。尽管如此，SLSQP的性能对于当前类别的设计问题来说已经足够，且优化时间只是设计过程总时间的很小一部分。
任何基于梯度的算法的主要目标都是更新设计空间中的变量，以将目标（I）驱动到最优值。为此，优化器根据解的梯度和步长评估一个最优搜索方向。SLSQP优化算法解决一系列子问题，每个子问题都在约束线性化的情况下优化目标函数的二次模型。SLSQP使用Hans-Powell拟牛顿（基于近似二阶信息）方法和Broyden-Fletcher-Goldfarb-Shanno（BFGS）[28]方法来更新Hessian矩阵。近似Hessian矩阵是基于前几次迭代的函数和梯度值构建的。SQP方法也可以看作是将牛顿方法应用于Karush-Kuhn-Tucker（KKT）最优性条件[29] [30]，它找到以下二次规划子问题的最优搜索方向d：
(26)mindf(xk)+?f(xk)Td+12?xx2L(xk,λk,σk)ds.t.g(xk)+?g(xk)Td≥0h(xk)+?h(xk)Td≥0
其中f(x)是目标函数I的值，g(x)和h(x)是不等式和等式约束的值，L是优化问题的拉格朗日函数，形式为L(x,λ,σ)=f(x)?λg(x)?σh(x)，λ和σ是拉格朗日乘数。也可以通过提供其他来源的值来改变优化计算Hessian的方式。可以通过改进的伴随方法来计算Hessian，从而进一步减少工具链优化步骤所需的努力。

2.5. 中间件工具
中间件是一种计算机代码，为HMB3工具集中的应用程序提供服务，并允许它们交换数据。由于每个工具只与中间件通信，而不是为应用程序之间的交换编写专用代码，因此实现通信和输入/输出变得更加容易。它还提供了所有必要的构建块，使开发人员能够快速使用此框架在HMB3工具集中添加未来的应用程序。
图1显示了HMB3和MSC/NASTRAN如何通过中间件进行通信。HMB3将当前的表面网格和流动解发送给中间件。中间件处理数据以生成MSC/NASTRAN使用的空气动力载荷，然后使用系统调用启动分析。一旦MSC/NASTRAN完成处理，结果会被插值到HMB3的表面网格上，并将其传递给HMB3。
下载：下载高分辨率图像（356KB）
下载：下载全尺寸图像
图1. 使用中间件将HMB3与MSC/NASTRAN耦合。
由于中间件是一个单核进程，它具有CPU和内存限制，这意味着它不能处理完整的网格或流动解。例如，解决方案和网格的存储需求大约是每1000万个单元1GB的RAM。目前，中间件具有以下功能：
a) 读取并初始化数据结构。
b) 与HMB3和MSC/NASTRAN协作。
c) 根据给定的边界条件标志计算截面的和完整的叶片载荷。
d) 为MSC/NASTRAN写入结构载荷数据。
e) 处理来自MSC/NASTRAN的位移和正常模式。
f) 解决结构位移的模态位移方程。
g) 在结构网格和CFD网格之间插值位移数据。
该功能允许使用结构位移的模态方程或使用MSC/NASTRAN中的完整FEM和HMB3计算的叶片载荷来进行悬停和前飞状态下机翼和转子的空气弹性计算。
已经实现了一维梁、二维板和三维有限元（FE）模型的加载交换。移动最小二乘法[31]用于平滑和插值CFD和CSD网格之间的数据。该方法从空间中任意固定点的加权最小二乘公式开始。然后，将该点移动到整个参数域上，对每个点分别计算和评估加权最小二乘拟合。然后使用Luke等人提出的逆距离加权（IDW）方法[32]变形体积网格。IDW是一种多变量插值方法，它使用一组已知散布点的值的加权平均值来计算未知点的值。分配给已知点值的权重与已知点和未知点之间的逆距离成正比。

2.6. 模态方法
空气弹性系统的模态分析从描述其行为的方程组开始。这些方程的一般形式为：
(27)[M]U¨(t)+[C]U˙(t)+[K]U(t)=F(t)
其中[K]、[C]和[M]分别是系统的刚度、阻尼和质量矩阵，U(t)和F(t)是时间t的广义位移和力向量。在模态方法中，位移向量U(t)被限制在所有可能位移的一个子集内，允许前m个模态形状，每个模态形状都有一个缩放因子。
(28)U(t)=sumi=1m?i·αi(t)=[Φ]A(t)
其中?i是第i个模态形状，αi(t)是时间t的第i个缩放因子。由于模态矩阵[Φ]与时间无关，方程28可以代入27并乘以[Φ]T，得到：
(29)[Mm]A¨(t)+[Cm]A˙(t)+[Km]A(t)=[Φ]TF(t)
其中[Mm]=[Φ]TM[Φ]，[Cm]=[Φ]TC[Φ]，[Km]=[Φ]TK[Φ]。然后对模态形状进行缩放，使得模态质量矩阵变为单位矩阵[I]，模态刚度矩阵变为对角矩阵，其条目等于模态频率ωi的平方。假设每个模态的阻尼与其临界阻尼成正比，则方程29变为：
(30)αi¨(t)+2ζiωiαi˙(t)+ωi2αi(t)=[Φ]TF(t)
其中ζi是结构阻尼系数。然后可以独立计算αi，并使用方程28通过线性叠加组合它们以获得位移向量。
在这项工作中，模态方程30也通过谐波平衡公式求解，假设αi(t)缩放因子也与流动变量具有相同的频率。这导致方程30的时间导数被替换为它们的谐波平衡对应物，从而得到方程：
(31)(ωHB2[D]2+ωHB[D]ζ˙iωi+ωi2)αi=[Φ]TF
模态方程的解是线性的，其规模取决于谐波平衡快照的数量。线性系统针对每个固体表面和结构模态求解，以获得模态振幅。线性系统使用高斯消元法求解。由于方程以维度单位（模态频率以rad/s为单位）求解，因此还需要将谐波平衡频率从归一化频率转换为转子的角频率。大约每2000次谐波平衡迭代进行一次空气弹性更新。对于修剪后的前飞计算，在每个修剪状态需要2到3次弹性更新，在整个计算过程中需要进行3到4次重新修剪。这允许在结合伴随方法的情况下优化非稳态转子流动和空气弹性效应。

2.7. 参数化
转子叶片使用表面参数化来改变其形状。在这种情况下，基线叶片的参数化由扭转、弦长 sweep 和 anhedral 设计变量组成。不同的设计变量组合会产生一系列叶片设计。一个具有七个系数的Bernstein多项式函数改变叶片的扭转，而弦长、sweep 和 anhedral 在叶片上的任意数量的预定径向站点上进行变化。在这个例子中，第一个径向站和叶片根部不允许有任何变化。所有叶片部分都允许设计参数的线性变化，除了倒数第二个和最后一个径向站之间的叶片部分，这部分可以是线性的或抛物线的。当前示例使用了八个弦长站、七个后掠角站和两个扭转角站（见图2）。增加更多的站和更多的设计参数可以扩展产生的形状家族，鉴于当前工作使用了伴随公式，主要的代价在于形状的复杂性和分析每个设计变化的原因，而不是CPU时间。可以在每个设计阶段应用表面平滑处理，以保证参数化站之间的高阶连续性。即使使用这个示例中的29个设计变量，也有足够的灵活性来捕捉几种叶片尖端形状（见表1）。

下载：下载高分辨率图像（36KB）
下载：下载全尺寸图像

图2. 沿叶片半径的控制站分布示例。

表1. AH-64A叶片优化研究的设计参数上下限。

设计参数 边界值
TWIST1 - TWIST7 ± 6度
CHORD at r/R = 0.5 [0.8c,1.2c]
CHORD at r/R = 0.70 [0.8c,1.4c]
CHORD at r/R = 0.75 [0.8c,1.4c]
CHORD at r/R = 0.80 [0.8c,1.4c]
CHORD at r/R = 0.85 [0.8c,1.4c]
CHORD at r/R = 0.90 [0.8c,1.4c]
CHORD at r/R = 0.943 [0.8c,1.4c]
CHORD at r/R = 1.0 [0.3c,1.4c]
SWEEP at r/R = 0.70 [-0.4c,0.4c]
SWEEP at r/R = 0.75 [-0.4c,0.4c]
SWEEP at r/R = 0.80 [-0.4c,0.4c]
SWEEP at r/R = 0.85 [-0.4c,0.4c]
SWEEP at r/R = 0.90 [-0.4c,0.4c]
SWEEP at r/R = 0.943 [-0.4c,0.4c]
SWEEP at r/R = 1.0 [-0.7c,0.4c]
ANHEDRAL at r/R = 0.70 [-0.1c,0.1c]
ANHEDRAL at r/R = 0.75 [-0.1c,0.1c]
ANHEDRAL at r/R = 0.80 [-0.1c,0.1c]
ANHEDRAL at r/R = 0.85 [-0.1c,0.1c]
ANHEDRAL at r/R = 0.90 [-0.1c,0.1c]
ANHEDRAL at r/R = 0.943 [-0.1c,0.1c]
ANHEDRAL at r/R = 1.0 [-0.2c,0.1c]

2.8. 优化框架

图3显示了前飞优化的流程。该图表明，气动弹性、叶片修剪和形状优化是主要协同工作的模块。然而，优化可以改变的是叶片的参数化和设计方面的选择，这体现了负责设计的工程师的设计理念。首先，使用稳态求解器或谐波平衡方法计算基线流动解。然后从基线流动解评估基线目标函数和约束条件。接着解决伴随问题，并使用伴随方法计算设计灵敏度。通常，伴随问题被近似求解以最小化计算成本。伴随方法的主要目的是提供足够精确的梯度，以正确地指导设计方向。成本函数、约束条件和梯度被输入到基于梯度的优化器中，以获得新的设计变量向量。根据新的设计变量向量，生成新表面的坐标和用于伴随方法的表面灵敏度。即使使用这个示例中的29个设计变量，也有足够的灵活性来捕捉几种叶片尖端形状（见表1）。

下载：下载高分辨率图像（259KB）
下载：下载全尺寸图像

图3. HMB3优化流程。

3. 结果与讨论

3.1. 基线叶片

基线叶片结合了Bell UH-1 Model 205的几何形状和Apache AH-64A旋翼的翼型部分。所得到的矩形叶片的纵横比为12.2，线性名义叶片扭转角为11度。叶片包含两个翼型部分：HH-02部分直到0.943R，以及叶片尖端的NACA64A006部分。NACA64A006部分经过修改，以考虑HH-02翼型稍厚的后缘。这两个翼型部分在0.04R的距离上进行了混合，考虑到不同零升力角度的部分。最后，根部切口的位置为0.15R。

3.2. 叶片结构模型

在这项工作中，使用了梁有限元模型。由于难以获得结构设计和属性，因此将使用基于AH-64A叶片已发布参数的近似动态模型。该模型在NASTRAN中开发，其得到的模态频率与其他计算工作和实验数据进行了比较。对于1D梁模型，这里采用了与[33]类似的方法。叶片使用CBEAM元素离散化，并使用刚性杆元素（RBAR）将梁的扭转可视化并映射到CFD网格中。对于基线叶片，叶片结构沿叶片跨度使用39个CBEAM元素离散化，如图4所示。

下载：下载高分辨率图像（123KB）
下载：下载全尺寸图像

图4. 基线旋翼叶片的结构1D梁模型。

在叶片优化中包含气动弹性效应需要调整叶片的结构模型以考虑形状的变化。某些优化研究包括气动弹性效应，但在优化研究过程中不更新结构模型[34]，[35]，这可能导致例如叶片沿跨度有较大变化时的扭转变形不准确。因此，即使结构属性没有优化，它们也会根据叶片形状的变化进行调整。在这项工作中，我们遵循了Stanger [36]的方法，Wilke [37]也使用了这种方法。调整结构属性是基于结构缩放定律，随着旋翼叶片弦长和后掠角的变化而变化。这里使用了单壳缩放，假设叶片蒙皮厚度保持不变，只有内部结构被缩放。这种缩放导致惯性矩线性缩放，刚度三阶缩放。后掠角的变化改变了沿叶片跨度的弹性轴位置，而扭转角改变了局部梁元素轴的位置。忽略了扭转角/反扭转角的变化。表2显示了应用的缩放因子。

表2. 随叶片形状变化而应用的结构属性缩放因子。α=cnew/cref [36]。

属性 缩放因子
质量 α
面积 α
极惯性矩 α
扭转刚度 α
展向刚度 α
弦向刚度 α

3.3. 悬停结果

首先在悬停状态下对叶片进行了优化。这展示了即使计算是在假设空间和时间周期性的稳态流动下进行的，气动弹性也与伴随方法的耦合。这种优化是针对气动弹性旋翼进行的，涉及功率最小化问题，同时约束推力。名义旋翼速度设置为270 RPM，导致尖端马赫数为0.606，尖端雷诺数为874万。为了比较目的，表3显示了优化前基线叶片的推力、扭矩和优值。选择集体角度为7.7度时CT=0.0085作为优化的约束推力，此时优值开始下降。

表3. 集体角度对基线叶片推力、扭矩和优值的影响。

集体角度 推力（N） 扭矩（Nm） 优值

4.0 105 229 30.0
0.6 25.0 131 0.63
1.0 150 0.67 16.0
1.6 203 90 0.68
2.7 98 67 0.67
5.0 80 60 0.65
2.8 230 0.65

在优化的每一步中，都会根据新的翼型形状更新结构模型，然后进行气动弹性悬停计算。最后，启动一个单独的伴随计算，将气动弹性变形内置到网格中。优化过程持续进行，直到目标函数或设计变量变化达到设定的容忍度。优化框架的整体性能以及各个组件的性能在之前的工作中已经进行了详细介绍[4]，[14]，[24]。然而，包含气动弹性是新颖的，这里首次提出。

下载：下载高分辨率图像（259KB）
下载：下载全尺寸图像

图3. HMB3优化流程。

3. 结果与讨论

3.1. 基线叶片

基线叶片结合了Bell UH-1 Model 205的几何形状和Apache AH-64A旋翼的翼型部分。所得到的矩形叶片的纵横比为12.2，线性名义叶片扭转角为11度。叶片包含两个翼型部分：HH-02部分直到0.943R，以及叶片尖端的NACA64A006部分。NACA64A006部分经过修改，以考虑HH-02翼型稍厚的后缘。两个翼型部分在0.04R的距离上进行了混合，考虑到部分不同的零升力角度。最后，根部切口的位置为0.15R。

3.2. 叶片结构模型

在这项工作中，使用了梁有限元模型。由于难以获得结构设计和属性，因此将使用基于AH-64A叶片已发布参数的近似动态模型。该模型在NASTRAN中开发，其得到的模态频率与其他计算工作和实验数据进行了比较。对于1D梁模型，这里采用了与[33]类似的方法。叶片使用CBEAM元素离散化，并使用刚性杆元素（RBAR）来可视化和将梁的扭转映射到CFD网格中。对于基线叶片，叶片结构沿叶片跨度使用39个CBEAM元素离散化，如图4所示。

下载：下载高分辨率图像（123KB）
下载：下载全尺寸图像

图4. 基线旋翼叶片的结构1D梁模型。

在叶片优化中包含气动弹性效应需要调整叶片的结构模型以考虑形状的变化。某些优化研究包括气动弹性效应，但在优化研究过程中不更新结构模型[34]，[35]，这可能导致例如叶片沿跨度有较大变化时的扭转变形不准确。因此，即使结构属性没有优化，它们也会根据叶片形状的变化进行调整。在这项工作中，我们遵循了Stanger [36]的方法，Wilke [37]也使用了这种方法。调整结构属性是基于结构缩放定律，随着旋翼叶片弦长和后掠角的变化而变化。这里使用了单壳缩放，假设叶片蒙皮厚度保持不变，只有内部结构被缩放。这种缩放导致惯性矩线性缩放，刚度三阶缩放。后掠角的变化改变了沿叶片跨度的弹性轴位置，而扭转角改变了局部梁元素轴的位置。忽略了反扭转角/正扭转角的变化。表2显示了应用的缩放因子。

表2. 随叶片形状变化而应用的结构属性缩放因子。α=cnew/cref [36]。

属性 缩放因子
质量 α
面积 α
极惯性矩 α
扭转刚度 α
展向刚度 α
弦向刚度 α

3.3. 悬停结果

首先在悬停状态下对叶片进行了优化。这展示了即使计算是在假设空间和时间周期性的稳态流动下进行的，气动弹性也与伴随方法的耦合。这类似于零模优化。优化是针对气动弹性旋翼进行的，涉及功率最小化问题，同时约束推力。名义旋翼速度设置为270 RPM，导致尖端马赫数为0.606，尖端雷诺数为874万。为了比较目的，表3显示了优化前基线叶片的推力、扭矩和优值。选择集体角度为7.7度时CT=0.0085作为优化的约束推力，此时优值开始下降。

表3. 集体角度对基线叶片推力、扭矩和优值的影响。

集体角度 推力（N） 扭矩（Nm） 优值

4.0 105 229 30.0
0.6 25.0 131 0.63
1.0 150 0.67 16.0
1.6 203 90 0.68
2.7 98 67 0.67
5.0 80 60 0.67
8.0 60 0.65

在优化的每一步中，都会根据新的翼型形状更新结构模型，然后进行气动弹性悬停计算。最后，启动一个单独的伴随计算，将气动弹性变形内置到网格中。优化过程显示出显著的好处，见表4，气动弹性悬停优化叶片的优值增加到0.789。

表4. 集体角度对气动弹性悬停优化叶片推力、扭矩和优值的影响。

集体角度 推力（N） 扭矩（Nm） 优值

5.0 254 82 480 0.61
7.0 074 60 52
1.0 130 0.78 68
1.6 203 90 0.78
9.0 136 30 0.72

基线和气动弹性悬停优化叶片的翼型形状如图5所示。在叶片优化中包含气动弹性效应需要调整叶片的结构模型。叶片扭转角的增加，见图6，以及先进的尖端形状非常明显。在叶片优化中包含气动弹性效应需要调整叶片的结构模型。优化收敛到最大可用扭转角，并在跨度的最后10%显示出强烈的反扭转角/正扭转角分布。高扭转角和非平面翼型形状是性能提升的主要驱动力。优化叶片的表面压力分布显示出更有利的行为，如图7所示。基线叶片在叶片尖端附近的高载荷对于优化叶片完全消失了。优化叶片显著减轻了尖端的载荷，仅在叶片前缘附近有轻微的吸力。优化叶片设计的尖端涡流场显示出较弱的尖端涡流。然而，由于尖锐的反扭转角/正扭转角分布，叶片的脱落移动了复杂的涡流结构。

下载：下载高分辨率图像（426KB）
下载：下载全尺寸图像

图5. 基线和气动弹性悬停优化叶片的形状。

下载：下载高分辨率图像（115KB）
下载：下载全尺寸图像

图6. 基线和气动弹性悬停优化叶片的扭转角。

下载：下载高分辨率图像（474KB）
下载：下载全尺寸图像

图7. 在相同集体角度下原始叶片和优化叶片的比较。

图8显示了原始叶片和优化叶片的优值，显示出显著的改进。一些改进是因为扭转角现在增加了。然而，还有另一个非常重要的贡献：叶片尖端与尖端涡流之间的间隙，如图9所示，集体角度为9.5度的情况。

下载：下载高分辨率图像（158KB）
下载：下载全尺寸图像

图8. 原始叶片和优化叶片之间的悬停优值比较。

下载：下载高分辨率图像（332KB）
下载：下载全尺寸图像

图9. 使用涡度可视化显示悬停优化叶片后的流动（θ75=9.5度）。

3.4. 前飞结果

使用模态谐波平衡方法对基线旋翼叶片进行了计算。模拟的飞行条件设置为CT=0.009，三个不同的前掠角μ=0.25、μ=0.3和μ=0.35。通过观察叶片拍动模式的位移和速度，如图10所示，四个谐波平衡模式足以充分解决模态振幅随时间的变化。

下载：下载高分辨率图像（132KB）
下载：下载全尺寸图像

图10. 随谐波平衡模式数量增加拍动系数的收敛。

悬停优化叶片也在相同的飞行条件CT=0.009下进行了计算，三个不同的前掠角μ=0.25、μ=0.3和μ=0.35。结果如图11所示。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像

图11. 在CT=0.009和三个不同前掠角μ=0.25、μ=0.3和μ=0.35下，悬停优化叶片的压力系数比较。

使用悬停参数化，计算了不同前掠角下的设计变量灵敏度，如图12所示。在所有前掠角下，灵敏度都表明减少叶片扭转角是有益的，且这种好处随着前掠角的增加而增加。对于弦长、后掠角和反扭转角，不同前掠角下的灵敏度显示出类似的趋势。对于弦长，最大的好处可以在尖端区域减少弦长时获得。对于后掠角，先向后然后向前再向后似乎都有好处。最后对于反扭转角，从尖端区域去除它可以获得最大的好处。这些变化都会导致悬停性能的下降。

下载：下载高分辨率图像（295KB）
下载：下载全尺寸图像

图12. 在CT=0.009和三个不同前掠角μ=0.25、μ=0.3和μ=0.35下，悬停优化叶片的设计灵敏度比较。

前飞优化从悬停优化叶片开始，而不是基线叶片，其扭转角分布固定为悬停优化叶片获得的扭转角分布。做出这一决定的原因是为了在减轻由于高叶片扭转角导致的前飞性能 penalty 的同时，保持悬停性能的大幅提升。最终叶片设计的压力系数如图13所示。最终设计的尖端弦长减小，尖端区域的后掠角更加明显，使得CQ从0.001260降低到0.001231，减少了2.5%。在前飞优化之后，还获得了最终形状的悬停性能结果。图8还显示，最终叶片的悬停性能没有受到影响。这一点从图中的虚线和实线可以明显看出，这些线条分别代表了在刚性或弹性假设下，经过悬停优化处理的叶片以及重新优化后适用于前飞状态的叶片的悬停性能。总体而言，由于采用了受悬停优化结果限制的叶片扭转设计，悬停性能得到了保持。另一方面，前飞性能的提升得益于在叶片尖端区域允许进行的弦长、后掠角和扭转角的变化，这些变化是通过所采用的参数化方法实现的。

下载：下载高分辨率图片（426KB）
下载：下载全尺寸图片

图13. 前飞优化叶片在CT=0.009和迎角μ=0.3.4时的压力系数

结论
本文展示了一种基于HMB3流动斜面的耦合伴随 Aeroelastic 叶片优化方法，该方法适用于悬停和前飞状态。考虑到问题的复杂性，该方法具有鲁棒性和高效性，其中谐波平衡方法有助于避免时间步进计算，并且不依赖于设计参数化中的变量数量。从矩形叶片开始，悬停效率提高了最大FoM的15%，而这并没有对前飞性能产生负面影响。该方法是一种比基于元模型的解耦技术更为现实的选择，为探索叶片设计空间提供了良好的框架，无需创建结果数据库或进行曲线拟合。未来，该方法将应用于不同的初始叶片设计，以防止初始叶片形状与最终叶片形状之间存在强烈的依赖性。

作者贡献声明
A. Woodgate：撰写——原始草稿、可视化、方法论、研究、形式分析、数据整理、概念化。
George N. Barakos：撰写——审阅与编辑、撰写——原始草稿、可视化、监督、软件开发、资源管理、项目协调、方法论、资金获取、概念化。