1. 引言

对旋翼尾迹流场的流动结构与演化规律的物理机理清楚认识是发展旋翼空气动力学必须解决的关键问题之一。它是研究直升机飞行特性、旋翼载荷、振动及噪声、操纵性及稳定性的理论基础，也是深入研究直升机气动干扰问题的重要手段。直升机在飞行过程中，旋翼桨叶会脱出强烈集中的桨尖涡，对旋翼流场具有强烈的影响，而且这些脱出的桨尖涡与桨叶会发生相互干扰，从而在桨叶上产生出强烈的冲击气动载荷和噪声。旋翼的气动特性与其尾迹流场密切相关，准确的尾迹流场信息对旋翼气动计算、性能计算起着决定性的作用；而对旋翼尾迹流场的准确试验测量，又反过来制约着已有理论的可用性和旋翼理论研究的进一步深入与发展。由于桨尖涡是旋翼尾迹流场中的主导结构，测量和研究旋翼尾迹流场的桨尖涡就成为非常重要的研究内容。

多种定性流动显示及定量流场测量方法已经被运用在对旋翼尾迹的研究上。对于旋翼滑流边界的收缩和桨尖涡的运动轨迹，可以通过烟流流场显示及纹影摄像等方法测定。而对于旋翼尾迹速度分布测量，经历了从最开始的皮托静压探头测量、热线测量到非接触式的激光多普勒测速(LDV)和粒子图像测速(PIV)技术的过程。由于传统二维PIV (平面两速度分量)具有测量整个二维平面内瞬时速度场的能力，被认为比LDV更适合用于测量非定常旋翼桨尖涡流场 [1]。近年来，借助于LDV和PIV技术，国内外展开了一系列旋翼尾迹流场的测量试验研究，获得了大量的桨尖涡结构和发展演化试验数据 [2]。袁红刚等人使用二维PIV技术研究了悬停状态下旋翼模型的轴向诱导速度场，确定了桨尖涡位置和尾迹边界随旋翼转速和总距的变化规律 [3]。李春华等人通过二维PIV技术对悬停和前飞状态下模型旋翼流场进行了测量，对比分析了不同前飞速度、总距、转速、方位等条件下的旋翼速度场和桨尖涡运动和演化规律，获得了在悬停和前飞条件下旋翼流动特性 [4]。随着PIV技术的发展，立体PIV (Stereo PIV)技术(平面三速度分量)也广泛应用于旋翼尾迹流场的测量中，获得了桨尖涡轴向速度的变化规律。为了解决PIV测量在全局流场信息获取和局部涡核区域空间解析力不足的矛盾，Rafflel等人使用测量窗口大小不同的两套立体PIV系统对桨尖涡进行了同步测量，得到了精细的涡核轨迹和涡核内部流动特性的实验结果 [5]。

在理论研究中，人们常常把螺旋涡模型用来等效替代桨尖涡。直线涡(点涡)、涡环、螺旋涡是涡旋动力学研究的三大基本对象 [6]，如图1所示。相对于直线涡，对螺旋涡诱导的流动，毕奥萨法尔积分没有一个封闭的形式。对于点涡，它不能表示为一个简单的极点，也不能表示为一个无限薄的涡环的完全椭圆积分。

对于螺旋涡的研究历史十分悠久，可追溯到1880年Kelvin关于柱状涡的螺旋扰动的著名著作，即大螺距螺旋涡 [7]，相对于直线涡，对螺旋涡诱导的流动，Biot-Savart积分没有一个封闭的形式。对于点涡，它不能表示为一个简单的极点，也不能表示为一个无限薄的涡环的完全椭圆积分。因此，Maykopar指出，螺旋涡的诱导速度必须通过由Biot-Savart积分产生的特殊函数来计算，或者通过Kapteyn型贝塞尔函数系列来计算 [8]。Ricca对几种推导螺旋涡速度的方法进行了比较，重点讨论了扭转对速度的影响 [9]。Hardin推导出了Kapteyn级数的完整三角形式，给出了螺旋涡在柱坐标下 $\left(\rho ,\theta ,z\right)$ 三维速度的各个分量 [10]。Kuibin和Okulov对Kapteyn级数主要部分提出了一种有效的改良方法，得到了更精准的解析解 [11]。Joukowsky清楚地证明了螺旋涡自激速度的形式化，表示在忽略毕奥萨伐积分中的某些项后，可以转化为密切涡环的自激运动，这一方法由Moore和Saffman正式提出来 [12] [13]。

正是基于上述要求，本文拟建立涡环模型来验证插值算法对三维流场重构的可行性。在柱坐标 $\left(\rho ,\theta ,z\right)$ 下将涡环沿将涡环沿 $\theta =\pi /2$ 顺时针旋转一定的角度(30度、45度等)后，基于MATLAB重新计算旋转后新的诱导速度场以及在不同方位角下每个平面内速度绝对值大小的等值线和三维等值面。且与之前不同的工作是在于取多个 $\theta$ 下的r-z平面，计算出不同平面下的速度大小，基于MATLAB插值程序从而重构出整个三维空间速度场，为在直升机旋翼尾流用立体PIV测量不同切面下的速度场，并用插值算法重构出整个三维空间速度场的实验提供理论支撑和可行性分析。

2. 涡环模型构建和结果分析

2.1. 正放的涡环模型

如图2考虑半径为 $a=1$ 的涡环，取直角坐标系O-xyz，涡环所在的平面为x-y平面，z轴通过圆心O，同时我们取柱坐标系 $\left(\rho ,\theta ,z\right)$，两坐标系之间的关系为 $x=\rho \cos \theta ,y=\rho \sin \theta$，由于轴对称性，通过Oz轴的所有平面上的运动都是一样的，因此不失普遍性可考察平面 $\theta =0$ 处的流体运动。

2.2. 涡环速度公式

假设给定圆形涡丝的强度是 $\Gamma$，在无散度有涡旋场中，即 $\nabla \cdot u=0$， $\nabla \times u=\Omega$，推出 $u=\nabla \times A$，A称为矢势，将之代入 $\nabla \times u=\Omega$，即 $\nabla \times \left(\nabla \times A\right)=\Omega$，根据场论中基本运算公式可以得到：

$A\left(\rho ,z\right)=\frac{\Gamma }{2\pi }\sqrt{\frac{a}{\rho }}\left\{\left(\frac{2}{\kappa }-k\right)\kappa \left(k\right)-\frac{2}{\kappa }E\left(k\right)\right\}$ (1)

$\kappa \left(k\right)={\displaystyle {\int }_0^{\pi /2}\frac{\text{d}\beta }{\sqrt{1-k^2{\sin }^2\beta }}}$ (2)

$E\left(k\right)={\displaystyle {\int }_0^{\pi /2}\sqrt{1-k^2{\sin }^2\beta }\text{d}\beta }$ (3)

$k^2=\frac{4\rho a}{z^2+{\left(\rho +a\right)}^2}$ (4)

$u_{\rho }=-\frac{\partial A}{\partial z},u_{\theta }=0,u_z=\frac{\partial \left(\rho A\right)}{\rho \partial \rho }$ (5)

令涡环半径 $a=1\text{m}$，在柱坐标系下，取 $\theta =0$ 的平面，如图3所示：

在 $-0.9\text{m}<z<0.9\text{m},0.1\text{m}<r<1.9\text{m}$ 的平面上，引入无量纲量 $r^{\ast }=r/a,z^{\ast }=z/a$ 基于Matlab可以算出 $u_{\rho },u_z$ 的速度等值线图以及矢量图和流线图如下，从图4可以看出涡环诱导的速度场在涡核两端严格对称的。越靠近涡核处速度等值线越密集。

以z轴作为对称轴线，从图中可以看出，涡核附近的流线仍为闭合曲线，这大概是因为，涡核附近旋转运动的速度以 $1/\rho$ 阶次趋于无穷，而涡丝本身的运动则是 $In\rho$ 阶次的无穷大，与 $1/\rho$ 相比是高阶小量，因此起主要作用的仍为旋转运动，流线呈封闭形式。

再引入速度绝对值大小，也就是 $U_{abs}=\sqrt{u_{\rho }^2+u_z^2}$，计算得到如图，为了便于观看，只画了某些特定值的等值线图。

因为涡环诱导的速度场与 $\theta$ 无关，所以只需要速度公式对 $\theta$ 做一次0到 $2\pi$ 的循环，即可得到速度绝对值大小的等值面图。选取 $U_{abs}=8$ 来进行作图，如图5所示，从图中可以清晰地看出，涡环的速度等值面是呈圆环分布，且涡环速度等值面的切面为圆形。

2.3. 坐标轴旋转下的涡环速度场

在正放的涡环模型中，地面坐标系保持不动，将涡环的坐标系沿着y轴，也就是沿着 $\theta =2/\pi$ 顺时针旋转一定的角度 $\alpha$，本文采取旋转30度，如图6(a)所示，旋转之后的涡环模型图如图6(b)所示，图中标示出了地面坐标系和涡环自身坐标系。

坐标轴旋转之后，可以找到新旧坐标轴之间的关系如下：

$y^{\prime }=y$ (6)

$x^{\prime }=x\cos \alpha +z\sin \alpha$ (7)

$z^{\prime }=z\cos \alpha -x\sin \alpha$ (8)

但是不能直接将上式直接带入使用，需要用新坐标系表示出旧坐标系，如下：

$y=y^{\prime }$ (9)

$x=x^{\prime }\cos \alpha -z^{\prime }\sin \alpha$ (10)

$z=z^{\prime }\cos \alpha +x^{\prime }\sin \alpha$ (11)

再将 $\alpha =\pi /6$ 代入上式，再代回涡环速度公式中，在MATLAB里进行计算，即可得到坐标轴旋转后涡环的速度绝对值大小的等值面，如图7所示。

从上面几张视图可以看出，坐标变换后倾斜的角度为 $\pi /6$，且是沿着y轴顺时针旋转的，侧面反映出公式的正确性。且 $U_{abs}=20$ 上的点明显比 $U_{abs}=10$ 的少，所以其圆环的厚度小于 $U_{abs}=10$ 的厚度。

2.4. 坐标轴旋转后的涡环不同切面速度分布

在之前的工作介绍了坐标旋转之后的等值面图，现在研究各个切面上速度绝对值大小的分布。坐标系依旧采用柱坐标，在r-z平面内，得到不同 $\theta$ 下速度的分布。柱坐标与笛卡尔直角坐标系之间的变换关系为： $x=\rho \cos \theta ,y=\rho \sin \theta$，将上述坐标变换公式代入此式，再令 $\theta$ 等于不同的值，即可得到不同切面的速度绝对值分布。

令 $\theta =0$，切面示意图如图8所示：

在 $-2\text{m}<z<2\text{m},0.1\text{m}<r<1.2\text{m}$ 的平面上，引入无量纲量 $r^{\ast }=r/a,z^{\ast }=z/a$，得到在坐标轴旋转后 $\theta =0$ 的速度分布，与未旋转之前的进行对比，如图9所示。

对比前后两次图，可以发现涡核的中心发生了转变，其变化规律满足坐标变换公式。等值线上的点分布也随之发生了改变，旋转之后速度大小的值分布的更广了。

再比较坐标轴旋转之后的 $\theta =0$ 和 $\theta =\pi$ 这两个对称平面的速度等值线图，如图10所示。

从图中可以看出，涡核中心和等值线的分布完全呈对称形式，也证明了旋转之后速度公式的正确性。

令 $\theta =\pi /4$，切面示意图如图11所示。

从图12可以明显的看出，涡核中心随之发生了变化，其变化规律也满足上述的坐标变换公式。在坐标轴倾斜之后涡环速度场通过修改 $\theta$ 的值，都可以在MATLAB中计算得到任意切面下r-z平面内速度绝对值大小的分布。

3. 插值法重构倾斜涡环速度场

假设我们在做PIV实验测量旋翼螺旋涡尾迹的时候，通过测量不同方向角下的r-z平面上速度大小，如果测量全部的方向角，会增加实验成本以及时间。基于上述要求，可以通过测量几个等分方向角平面内的速度大小，通过插值的方法，得到其他平面内速度的信息，从而重构整个三维流场。

3.1. 插值算法简介

分段三次样条插值：分段是把区间 $\left[a,b\right]$ 分成， $\left[\left(x_0,x_1\right),\left(x_1,x_2\right),\cdots ,\left(x_{n-1},x_n\right)\right]$ 共有 $n+1$ 个点，其中两

个端点 $x_0=a,x_n=b$。三次样条就是说每个小区间的曲线是一个三次方程，三次样条方程满足以下条件：

1) 满足插值条件，且在节点处函数值相等；

2) 插值函数的一阶导数、二阶导数连续。

Hermite三次插值：与三次样条函数不同点在于Hermite插值只需要满足节点处一阶导数连续且相等，所以计算效率优于三次样条插值，但曲线不如其光滑。

在MATLAB里直接调用‘spline’，即三次样条插值函数的命令以及‘pchip’，即Hermite三次插值函数命令。

3.2. 插值算法结果与分析

在第二章坐标轴旋转后涡环速度场的讨论后，得到了不同切面下速度大小的分布。现在将 $\theta =0,\pi /12,\pi /6,\pi /4,\pi /3,5\pi /12$ 这几个点上作为等分区间，用三次样条插值和Hermite插值进行插值。以 $\theta =\pi /7$ 为例，从图13可以看出，由三次样条插值函数和Hermite插值出来的等值线的涡核中心和真实情况由一点偏离，这可能是由于网格点的数量和涡核附近诱导的速度极值差异性较大造成的。等值线形状大致一致，也说明了插值方法的可行性，除了在涡核附近有一些差异，在涡核的中远场区域速度大小都有很好的一致性。

从插值出来的数据显示，只有在涡核附近的点会有一定的误差，所以在涡核中心沿着r方向和z方向各取10个网格点来进行误差分析，如图14所示，x表示从涡核的r方向和z方向的网格点的个数，可以看出除了在涡核附近几个点误差相对大之外，其他点都有一个很好的一致性。涡核附近会有较大差异是因为选取网格点的数量和间距有关，因为涡核处的速度是趋于无穷大，这个时候网格的精度对数据大小的影响就十分大，所以在涡核处出现较大误差是可以接受的。

从图中可以看出除了在涡核附近几个点误差相对大之外，其他点都有一个很好的一致性。涡核附近会有较大差异是因为选取网格点的数量和间距有关，因为涡核处的速度是趋于无穷大，这个时候网格的精度对数据大小的影响就十分大，所以在涡核处出现较大误差是可以接受的。

图15是两种插值方法的相对误差分析，同样在涡核中心沿着上下和左右各取10个网格点来进行分析。x表示从涡核的r方向和z方向的网格点的个数。

从相对误差图来看，除了在涡核附近点以外，相对误差大小都在5%-10%左右，有些点的相对误差甚至为0，这也证明了误差是在可以接受的范围内。可以人为通过调节网格点数量和间距，使得涡核附近的点的插值精度更高，从而达到更优化。

由上述例子 $\theta =\pi /7$ 的切面上分别用三次样条插值和Hermite插值得到速度大小的结果和误差分析来看，用已知几个等分方位角的平面上速度大小通过插值的方法重构整个三维流场是切实可行的。为以后的实验提供了新的思路和理论基础以及可行性的分析。

4. 结论

本文利用涡环的速度诱导公式，将坐标轴旋转之后得到新的涡环模型，利用Matlab编程计算出旋转以后的不同速度大小的等值面，以及各个切面上的速度等值线。利用三次样条插值和Hermite插值，得到其他方位角上的速度大小，与公式计算得到的真实值比较有较好的一致性，为立体PIV实验测量旋翼尾迹从而通过插值的方法重构三维流场提供了理论模型基础和可行性分析。

参考文献

参考文献