叶片单元法

叶片单元法公式基于根据 Rajagoplan [270] 等人提出的方法。它主要用于对直升机旋翼在流场中产生的效应进行建模，但实际并未预先指定旋翼桨叶的确切几何。

此方法包含两个部分。第一步是标识计算域中属于虚拟盘体和应用了源项的网格单元的位置。第二步涉及到计算每个网格单元的动量源项矢量（在第一步中识别和标记）。

步骤 1：基于插值网格的叶片单元/容器的网格单元标记和分配

此步骤的相关信息由用户指定并通过虚拟盘体节点提供：

盘体几何
描述盘体的几何尺寸和方向
虚拟盘体分辨率
在极坐标中指定二维插值网格尺寸，表示实际叶片单元。翼弦、扭转和空气动力特征等叶片属性普遍存在于整个单元中。

源项分布：平均时间

对于“平均时间”法，源项根据插值网格上的每个叶片单元进行计算，然后传输到有限体网格。每个单元与有限体网格上的多个网格单元相关联。对于局部流动条件和建模转子之间的耦合，相关信息从最靠近插值网格上的叶片单元中心的有限体积网格单元获取。因此，指定插值网格尺寸实际上就是计算转子对流场的影响。以下草图说明了有限体网格（红色）与虚拟盘体插值网格（黑色）之间的流场耦合。

标记过程本身分为两步。第一步，将圆盘区域内的所有网格单元标记为接收源项。第二步，为已标记的每个网格单元分配与其相对插值网格的位置对应的容器/单元编号。因此，在第一步中标记的每个网格单元属于插值网格的某一单元/容器。

此方法要求每个单元/容器在其域中至少有一个网格单元。如果有限体网格非常粗糙且插值网格过于精细，则需要降低分辨率。以下两个图显示了两种不同分辨率的插值网格。对于第二个图中的精细插值网格，无有限体积网格单元属于最内环半径处某些单元/容器所在的域。

源项分布：准确时间

“准确时间”法可跟踪叶片的运动，并仅将源项添加到与叶片位置对应的体积网格单元。对于具有四个叶片的转子，将根据在盘体几何节点中指定的叶片形状和位置，对网格单元进行标记，如下所示。

为确保与插值网格成功耦合，有限体网格必须足够精细，在叶片宽度方向上至少标记 5 个网格单元。应向转子的整个扫掠区域应用相近的网格分辨率。在上述示例中，标记了 10 个网格单元。

使用默认插值模式时，Simcenter STAR-CCM+ 支持不同叶片的标记的网格单元的重叠。源项将相应地计算：

步骤 2：计算动量源项强度

叶片单元法的第二个部分计算插值网格单元的动量源项强度。根据转子几何（即，方向和原点），创建了盘体局部笛卡尔坐标系 $(X′, Y′, Z′)$ 。此笛卡尔坐标系的原点位于转子中心和 $Z′$ 轴上，方向对应于盘体的推力方向。 $Z′$ 轴垂直于旋转平面，而 $X′$ 和 $Y′$ 轴则位于盘体平面中。以下草图演示了盘体局部坐标系：

旋转速度 $ω$ 始终沿 $Z′$ 方向。从 + $Z′$ 方向查看时，+ $ω$ 对应于逆时针 (CCW) 旋转，- $ω$ 对应于顺时针 (CW) 旋转。

考虑到叶片俯仰，定义了另外两个坐标系：

基于转子的圆柱坐标系 $(r, ϕ, z)$ ，其相对于基本局部坐标系 $(X′, Y′, Z′)$
笛卡尔坐标系，相对圆柱坐标系的基础。它称为俯仰坐标系。

定义了俯仰坐标系的基础 $(s, ϕ, n)$ ，使 $e_{s}$ 指向叶片的跨度方向。俯仰坐标系的基本矢量 $e_{ϕ}$ 与基于转子的圆柱坐标系的 $e_{ϕ}$ 指向同一方向。圆柱坐标系随叶片旋转，俯仰坐标系相对圆柱坐标系倾斜。给定叶片的方位定位 $ψ$ ，旋转圆柱坐标系，使 $e_{r}$ 与叶展对齐。俯仰坐标系按 $δ (ψ)$ （即，此方位定位的俯仰角）倾斜。所有计算在俯仰坐标系中执行。合力矢量将转换回基准坐标系，随之作为源项添加到动量方程中。

基准坐标系中的局部流场速度为：

v = v_{x} \hat{i} + v_{y} \hat{j} + v_{z} \hat{k}

(5135)

使用变换矩阵序列 $M_{1}, M_{2}, M_{3}$ 将局部流速将传递到俯仰坐标系：

[\begin{matrix} v_{s} \\ v_{ϕ} \\ v_{n} \end{matrix}] = M_{3} M_{2} M_{1} [\begin{matrix} v_{x} \\ v_{y} \\ v_{z} \end{matrix}]

(5136)

其中，这些变换矩阵将各种坐标系的基本矢量相关联，如下所示：

[\begin{matrix} e_{s} \\ e_{ϕ} \\ e_{n} \end{matrix}] = M_{3} [\begin{matrix} e_{r} \\ e_{ϕ} \\ e_{z} \end{matrix}] [\begin{matrix} e_{r} \\ e_{ϕ} \\ e_{z} \end{matrix}] = M_{2} [\begin{matrix} e_{x^{'}} \\ e_{y^{'}} \\ e_{z^{'}} \end{matrix}] [\begin{matrix} e_{x^{'}} \\ e_{y^{'}} \\ e_{z^{'}} \end{matrix}] = M_{1} [\begin{matrix} e_{x} \\ e_{y} \\ e_{z} \end{matrix}]

(5137)

叶片因旋转而产生速度，该速度可写入俯仰坐标系，如下所示：

$v^{b l} = M_{3} M_{2} M_{1} (ω \times r)$

其中， $r$ 为半径矢量。

因此，相对叶片的速度 $v^{r e l} = v_{s} e_{s} + v_{ϕ} e_{ϕ} + v_{n} e_{n}$ 由以下公式给出：

$\begin{matrix} v^{r e l} = v - v^{b l} \\ = M_{3} M_{2} M_{1} v - M_{3} M_{2} M_{1} (ω \times r) \end{matrix}$

要计算有效攻角

α

，需考虑以下叶片角：

诱导攻角（流入角） $β$
叶片俯仰角 $δ (ψ)$
几何俯仰角 $θ (ψ)$
几何扭转角 $θ_{t w} (r)$

叶片俯仰角计算如下：

δ = δ_{0} - δ_{1 c} \cos (ψ) - δ_{1 s} \sin (ψ)

(5138)

其中，

δ_{0}

为锥角，

δ_{1 c}

为纵向俯仰角，

δ_{1 s}

为横向俯仰角。时间导数由以下公式给出：

\dot{δ} = ω \frac{ⅆ δ}{ⅆ ψ} = ω (δ_{1 c} \sin (ψ) - δ_{1 s} \cos (ψ))

(5139)

俯仰速度：

v^{f l a p} = v_{s}^{f l a p} e_{s} + v_{ϕ}^{f l a p} e_{ϕ} + v_{n}^{f l a p} e_{n} = r \dot{δ}

(5140)

现在，诱导攻角计算如下：

β = \arctan [\frac{- (v_{n} + v_{n}^{f l a p})}{s i g n (v_{ϕ}) \sqrt{{(v_{ϕ})}^{2} + {(v_{s})}^{2}}}]

(5141)

其中， $s i g n (v_{ϕ})$ 考虑了高进程比的反流区域中 $v_{ϕ}$ 的方向。

根据 Eqn. (5141)，叶片俯仰角通过叶片俯仰速度影响流入角（诱导攻角），从而影响有效攻角 $α$ 。

俯仰角 $θ$ 定义如下：

θ = θ_{0} - θ_{1 c} \cos ψ - θ_{1 s} \sin ψ

(5142)

其中， $ψ$ 为方位角， $θ_{0}$ 为集合俯仰角， $θ_{1 c}$ 为横向循环俯仰角， $θ_{1 s}$ 为纵向循环俯仰角。

然后，最终有效攻角计算如下：

α = θ + θ_{t w} - β

(5143)

已知翼型遇到的攻角 $α$ 和叶片相对速度 $v^{r e l}$ ，计算马赫数或雷诺数。空气动力系数还用于计算马赫数和雷诺数。在模拟树中以表形式指定这些系数。计算的数取决于如何指定空气动力系数：指定为马赫数 $C_{l}$ , $C_{d} = f (α, M a) or$ 的函数还是雷诺数 $f (α, R e)$ 的函数。

然后，升力和阻力计算如下：

\begin{matrix} L′ = \frac{1}{2} \cdot ρ \cdot {| v^{r e l} |}^{2} \cdot C_{l} \cdot c \cdot d r \\ D′ = \frac{1}{2} \cdot ρ \cdot {| v^{r e l} |}^{2} \cdot C_{d} \cdot c \cdot d r \end{matrix}

(5144)

其中， $c$ 为叶弦。

要考虑转子叶尖因尾涡导致的叶尖损失，需要根据所选的叶尖校正法对升力 $\begin{matrix} L′ \end{matrix}$ 、攻角 $α$ 或升力和阻力 $L' a n d D'$ 进行校正。

然后，将力转换为法向和切向力（从风轴转换为旋翼俯仰坐标系）：

$\begin{matrix} {F′}_{n} = L′ \cos β - D′ \sin β \\ {F′}_{t} = L′ \sin β + D′ \cos β \end{matrix}$

由于叶展方向上不存在空气动力，因此 $F_{s} = 0$ 。在俯仰坐标系中，气动合力矢量可在叶片单元上组合为 $(F_{s}, F_{ϕ}, F_{n})$ 。

最后，力矢量转换为基准坐标系：

(\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}) = M_{1}^{T} M_{2}^{T} M_{3}^{T} (\begin{matrix} F_{s} \\ F_{ϕ} \\ F_{n} \end{matrix})

(5145)

作用于流体单元的瞬时力即为 $- F$ 。

将时间平均源项 $S = (S_{x}, S_{y}, S_{z})$ 添加到离散化动量守恒方程中。由于叶片仅用其总旋转时间的有限部分穿过控制体积，因此时间平均源项按时间分数进行比例缩放：

S = \frac{b Δ φ}{2 π} (- F)

(5146)

其中， $b$ 为叶片数， $Δ φ$ 为叶片在穿过叶片单元时行进的距离。

不需要对准确时间源项进行任何比例缩放。给定叶片单元中的指定源项在有限体网格中的相应网格单元上分布，与网格单元体积成比例。

数值切割算法

以旋翼飞机或直升机为背景，术语切割表示正确调整飞机控件和整体定位旋翼飞机，从而获得所需的稳定飞行条件。例如，对于悬停状态，切割转子的目标在于实现平衡，即，重心周围的净力总和与净力矩总和为零。Simcenter STAR-CCM+ 使用隔离切割（风洞切割）算法。隔离切割算法可计算转子俯仰控件，即：集合俯仰、横向循环俯仰和纵向循环俯仰。

主桨控件

通常，直升机有四个独立控件：三个用于主桨，一个用于尾桨。Simcenter STAR-CCM+ 实施孤立地看待转子，仅考虑转子控件。

主桨控件如下：

集合俯仰：
通过使用此控件，飞行员可以统一更改主桨的所有叶片的俯仰角，即，按相同的量更改。此更改会有效地影响主桨产生的总推力。
循环俯仰：
通过使用此控件，飞行员可以每转更改一次俯仰角。循环俯仰分类如下：
- 横向循环俯仰 $θ_{1 c}$ ：
  随 $\cos (ψ)$ 而变化。叶片在 0° 和 180° 下时，横向循环俯仰最大或最小。转子响应（以 90° 定相）导致横向（y 轴）上产生控制效果并有效地影响旋转力矩。
- 纵向循环俯仰 $θ_{1 s}$
  随 $\sin (ψ)$ 而变化。叶片在 90° 和 270° 下时，纵向循环俯仰最大或最小。转子响应（以 90° 定相）导致纵向（x 轴）上产生控制效果并有效地影响俯仰力矩。

因此，总俯仰角是方位角

ψ

的函数，由以下公式给出：

θ = θ_{0} - θ_{1 c} \cos (ψ) - θ_{1 s} \sin (ψ)

(5147)

算法

转子空气动力响应是其控制输入和空气动力环境的高度耦合函数。尽管转子产生的推力主要受集合俯仰影响，但是也受循环俯仰输入影响。同样，旋转与俯仰力矩与循环俯仰输入直接相关，但也会导致推力发生变化。在虚拟盘体的局部坐标系中定义：

\begin{matrix} T & = & T (θ_{0}, θ_{1 c}, θ_{1 s}) \\ M_{x} & = & M_{x} (θ_{0}, θ_{1 c}, θ_{1 s}) \\ M_{y} & = & M_{y} (θ_{0}, θ_{1 c}, θ_{1 s}) \end{matrix}

(5148)

其中， $T$ 为推力， $M_{x}$ 为旋转力矩， $M_{y}$ 为俯仰力矩。

此耦合方程组可以解释为控制输入矢量 $x$ 与转子响应矢量 $y$ 之间的关系，其中：

x = (\begin{matrix} θ_{0} \\ θ_{1 c} \\ θ_{1 s} \end{matrix}) and y = (\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})

每个飞行条件对应于转子响应矢量 $y$ 的状态，这需要控制输入矢量 $x$ 的特定状态，进而影响直升机的空气动力特征。

要研究直升机的空气动力性能，首先选择飞行条件。这意味着响应矢量 $y$ 已知，而不是控制输入矢量 $x$ 。使用输入矢量 $x$ 的初始猜测值开始迭代过程。然后迭代更新 $x$ ，直至解收敛到目标响应矢量 $y$ 。

要使迭代过程更新

x

，线性化此耦合方程组。然后，应用牛顿-拉夫森方法。要进行线性化，应用

x

y

泰勒级数展开：

y (x + Δ x) = y (x) + J Δ x + {O (Δ x)}^{2}

(5149)

仅考虑一阶项，Eqn. (5149) 重新排列为：

Δ x = J^{- 1} {y (x + Δ x) - y (x)}

(5150)

其中， $J$ 为相关物理量（即，响应矢量 $y$ ）相对控制输入 $x$ 的雅可比矩阵。

Eqn. (5150) 改写为：

{(\begin{matrix} Δ θ_{0} \\ Δ θ_{1 c} \\ Δ θ_{1 s} \end{matrix})}_{(i)} = {(\begin{matrix} \frac{\partial T}{\partial θ_{0}} & \frac{\partial T}{\partial θ_{1 c}} & \frac{\partial T}{\partial θ_{1 s}} \\ \frac{\partial M_{x}}{\partial θ_{0}} & \frac{\partial M_{x}}{\partial θ_{1 c}} & \frac{\partial M_{x}}{\partial θ_{1 s}} \\ \frac{\partial M_{y}}{\partial θ_{0}} & \frac{\partial M_{y}}{\partial θ_{1 c}} & \frac{\partial M_{y}}{\partial θ_{1 s}} \end{matrix})}_{(i)}^{- 1} \times {{(\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})}_{(d)} - {(\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})}_{(i)}}

(5151)

在 Eqn. (5151) 中，下标 (i) 对应于第 i 次切割迭代，d 对应于转子响应矢量 $y$ 的所需目标状态。

$x$ 进的初始猜测值可得出初始值 $y$ ，这通常不等于 $y_{d}$ 的目标状态。要继续进行切割迭代，必须使用冻结流场确定雅可比矩阵。为此，单独扰动控制输入矢量的每个风量，并计算转子响应。然后，使用二阶中心差分计算偏导数。给定 $J$ ，通过求解方程 Eqn. (5150) 确定控制输入矢量 $Δ x$ 校正。计算校正为亚松弛，可实现数值阻尼，然后用于更新控制输入矢量。重复该过程，直至解收敛到指定容差内的目标 $y$ 。

在时间平均值模式下，重复应用控制输入校正。源项的收敛将触发应用校正。

在时间精确值模式下，转子按初始条件旋转几圈，直到流成为周期流。在应用控制输入校正之前，修剪求解器将收集一转的轴数据以建立平均值。 $y$ 应用事件由源项的收敛或由两次校正之间的转数引导，它们在“修剪选项”属性中指定。

源项更新概述

以下流程图概述了叶片单元源项更新过程。当主流动求解器已准备好在动量方程中接受来自可选模型的源项贡献时，将触发此过程。当压力和速度等主流相关场可用时，计算叶片单元源项。

接下来，Simcenter STAR-CCM+ 检查叶片单元源项是否收敛到 Source Tolerance（请参见“用户指南 > 物理建模 > 虚拟盘体模型 > 叶片单元法 > 设置源收敛控制”）。此外，Simcenter STAR-CCM+ 还检查旋转速率跃升是否完成。如果任一检查失败，则叶片单元源项按源亚松弛因子进行亚松弛（请参见“用户指南 > 物理建模 > 虚拟盘体模型 > 叶片单元法 > 设置源收敛控制”），并添加到动量方程的源项中。然后，求解过程继续处理下一个模型。此步骤对应于上述流程图中第一个判定块的 no(否) 分支。

当第一个判定块计算为 yes(是) 时，Simcenter STAR-CCM+ 会检查主桨响应是否等于目标值。此检查对应于上述流程图中的第二个判定块，并解释为切割状态的收敛。如果达到用户指定的切割容差，将收敛切割状态（请参见“用户指南 > 物理建模 > 虚拟盘体模型 > 叶片单元法 > 设置切割收敛控制”）。如果切割状态已收敛，则求解过程将继续处理下一个模型。如果切割状态未收敛，则根据数值切割算法更新转子切割角。计算得出的切割角按切割亚松弛因子进行亚松弛（请参见“用户指南 > 物理建模 > 虚拟盘体模型 > 叶片单元法 > 设置切割收敛控制”）。

每次迭代时，执行源项计算。因此，导致流场发生变化的任何用户设置更改会影响叶片单元法源项矢量。务必注意的是，切割算法仅在叶片单元源项收敛到指定级别时才执行。如果外部流场高度不稳定，建议松弛源项上的收敛容差，使切割角更新可以起作用。对于极端情况，如果主桨在尾流中运行且所需的飞行条件为悬停，建议通过指定容差级别 100% 停用源项收敛检查。此设置可确保在每次迭代时执行该算法。