展向旋转平板Couette湍流中的多态现象

夏振华

(浙江大学 航空航天学院工程力学系, 浙江 杭州 310027)

0 引 言

旋转条件下的壁湍流问题是自然界和工业工程中常见的湍流现象之一。考虑地球自转作用，那么地球上的壁湍流问题几乎都可以视为旋转壁湍流问题。而航空发动机的转子、涡轮间的流动也是典型的旋转壁湍流问题。为了设计性能更好的航空发动机等旋转机械，我们急需了解旋转条件下的壁湍流的统计规律，以及对应的流动结构随旋转强度的变化规律。另一方面，航空发动机等旋转机械中关键湍流问题有很多，也非常复杂。为了更好地独立出旋转因素的影响，人们提出了一系列简化的标准模型问题，其中展向旋转平板Couette湍流(Spanwise rotating plane Couette flow, RPCF)就是其中之一[1-2]。

在旋转系统下，由于科氏力的存在，基本的稳定性分析告诉我们，旋转会根据其方向的选择性使得剪切流动变得更稳定或者更不稳定。如果旋转方向和平均流涡量同号(cyclonic rotation)，那么旋转会使得流动更稳定；如果两者异号(anticyclonic roation)，那么旋转会让流动变得更不稳定。在展向旋转的槽道湍流中，由于其平均速度存在极值，平均流涡量会在槽道中变号，因此旋转会在一部分槽道中使得流动变得不稳定；而在另一部分槽道中，其会使流动变得更稳定。而在RPCF中，其平均流涡量基本是相同符号的，因此旋转在整个槽道保持相同的效应，即或全场让流动更不稳定；或全场让流动更稳定[3-5]。

图1为RPCF流动的示意图。相距为2h的两块平板以相对速度2Uw沿x方向运动，而整个系统沿展向(z方向)以常角速度Ωz旋转。该系统有两个控制参数：雷诺数Rew=Uwh/ν和旋转数Ro=2Ωzh/Uw。

图1 RPCF简图[5]Fig.1 Sketch of the flow geometry in RPCF[5]

在过去的几十年中，很多研究组针对RPCF开展了系列的实验和数值模拟研究。在实验方面，Alfredsson教授和其团队做出了一系列工作[6-9]。特别需要指出的是，Tsukahara、Tillmark和Alfredsson[7]在2010年时总结并扩充了之前的实验工作，系统地研究了参数范围在0

在数值研究上，Andersson教授及其合作者首先针对Rew=1300时的RPCF开展了直接数值模拟(Direct numerical simulation, DNS)研究[13-17]。他们研究了弱旋转时的流动统计以及流动结构，中高旋转数时的旋转效应以及极强旋转时的流动统计。他们发现在一定旋转数范围里，湍流和二次流随着Ro增加而增强，但是超过一定的旋转数范围，其湍流和二次流则会随着Ro增加而减弱。在极高旋转数条件下，二次流消失，这个时候的湍流特征与无旋转时的流动特征差别很大。Salewski和Eckhardt[18]用DNS研究了四个不同雷诺数(Rew= 650、1300、2600和5200)下，以壁面剪切为表征的动量流随Ro的变化行为。他们考察了-0.04≤Ro≤0.6之间的算例，发现变化曲线在两个较低雷诺数情况下只有一个极值点；而在两个较高雷诺数下，变化曲线有两个极值点，其中一个极值点位于Ro=0.02附近(较小旋转数)。Gai等[5]也针对Rew=1300时的RPCF流动开展了一系列DNS模拟研究，他们分析了中心线处的平均剪切率、壁面摩擦雷诺数、体平均的总湍动能、体平均的二次流湍动能、体平均的剩余场湍动能等随着旋转数的变化规律。他们还分析了不同旋转数下的流动结构的异同。

近些年，本文作者与其合作者还发现了RPCF中存在多湍态现象[19-22]。在经典的湍流理论中，人们都假设雷诺数足够高的时候，湍流是各态遍历的，也就是说所有的湍流态之间都是等价的，或者说湍流态是唯一的。在这个假设下，虽然瞬时场跟初始条件有关，但是湍流的统计规律是与其初始流场无关的[23]。湍流的唯一统计态假设是经典湍流理论的基石，也是很多湍流结论和假设能够从低雷诺数外推到更高雷诺数的原因之一[24]。然而，近些年来，一些实验和数值模拟的结果显示，部分湍流中可能存在不唯一的湍流态[25]。在多态条件下，相同的控制参数，可能会出现不同的湍流统计和流动结构。例如，在Taylor-Couette流动中，Huismann等[24]就通过实验发现在Re为106量级下，通过不同途径达到相同的内外筒转速比a时，流动的Taylor涡对数和扭矩等结果在一定的参数范围里(0.17

本文将就RPCF中的多态现象的研究现状进行综述介绍。

1 RPCF的控制方程和数值解法

1.1 控制方程

RPCF的几何构型如图1所示。在通常的设置中，我们假设流向为x，法向为y，展向为z(或者x1、x2、x3)，其对应的速度为u、v、w(或者u1、u2、u3)。用Uw，h无量纲化之后，系统的控制方程为：

(1)

(2)

这里p是包含了离心力的有效压力(effective pressure)。其中雷诺数Rew和旋转数Ro的定义前文已给出。

1.2 数值解法介绍

在实际计算中，为了保证程序的精度，我们采用高精度Fourier-Chebyshev伪谱法对空间项进行离散。具体而言，在流向和展向周期方向，我们采用Fourier级数展开；而在法向采用Chebyshev多项式展开。在时间推进上，我们采用了半隐式离散格式，即线性项采用2阶精度的Crank-Nicolson根式而非线性项采用的是2阶的Adams-Bashforth格式。此外，我们用3/2法则消除混淆误差。值得说明的是，我们的方法比Bech和Andersson[14-15]和Barri和Andersson[17]的数值解法精度更高。Bech和Andersson用的是2阶精度的有限差分程序ECCLES[14-15]，而Barri和Andersson则用的是2阶精度的有限体积程序MGLET[17]。

在实际计算中，我们在大小为Lx×Ly×Lz的计算域上开展计算，离散的网格为Nx×Ny×Nz。具体的计算网格，在后续会做介绍。

1.3 雷诺分解和三重分解

(3)

由于在RPCF中存在二次流，因此我们在分析的时候，还可以开展三重分解[5,14,19]：

(4)

(5)

(6)

即湍动能和雷诺应力都可以分解为由二次流贡献的部分和剩余场贡献的部分。

2 结果与讨论

本部分将着重介绍相关的结果与讨论部分。

2.1 旋转数对湍流统计和流动结构的影响

Gai等人[5]针对RPCF开展了一系列的DNS模拟，基本控制参数如表1所示。其主要研究目标是旋转数对湍流统计和流动结构的影响。在分析该问题时，研究者们一般重点关注两个量：平均速度在中心线处的法向梯度Ψ和壁面摩擦雷诺数Reτ。

表1 计算参数Table 1 Computational parameters

图2(a、b)分别给出了计算中20个算例的中心线速度的法向梯度的平均值及其标准差以及壁面摩擦雷诺数Reτ随旋转数Ro的变化规律。可以看到，在Ro从0增加到0.02附近时，Ψ逐渐减小；而当Ro继续从0.02增加到0.9时，Ψ则随着Ro增大而增加。非常有意思的是，在Ro=0.02时，Ψ是负值，即在该旋转数条件下，平均流动存在局部回流现象。Salewski和Eckhardt[18]在Ro=0.04时也发现了平均速度的局部回流现象。不过需要说明的是，他们所采用的流向计算域过小，对流动结构和统计有一定的影响(在下文，我们将讨论这种差别)。而Kawata和Alfredsson[10]他们也通过实验测到过这种平均速度局部回流现象。相反地，壁面摩擦雷诺数Reτ随旋转数Ro的变化则呈现先增加后减小的趋势。其最大值约在Ro=0.15附近。在该问题中，Rew保持固定，因此Reτ可以表征壁面摩擦的大小。由图可以看出，在00.611时，旋转可以起到减阻的作用。

(b)

图3给出了体平均的湍动能[k]y以及二次流和剩余场贡献的部分[ks]y和[k″]y随旋转数Ro的变化曲线。其中，算子“[]y”定义如下

(7)

图3 体平均的湍动能及其分量随旋转数Ro的变化规律。相关量都用Ro=0时的无量纲化[5]Fig.3 Volume-averaged kinetic energy and its two shears (the contributions from the secondary flows and the residual field) as a function of Ro. All quantities are normalized by at Ro=0[5]

由于湍动能等项已经是时间、流向和展向平均之后的量，因此其法向平均可以认为是体平均的结果。由图中可以看出，总湍动能随着Ro是先增加后减小的，在Ro=0.25附近达到最大值。即当Ro<0.25时，旋转会增强湍流；而当Ro>0.25时，旋转相反会抑制湍流。二次流部分的湍动能的结果同总湍动能基本保持类似的趋势，在Ro>0.6时，二次流基本消失。剩余场部分的湍动能则体现出同总动能相反的趋势。它是先减小，后增大，最后又减小的过程。在Ro=0.1附近其达到极小值，而在Ro=0.5时，其已经非常接近总湍动能的值。

图4是六个不同旋转数下二次流的统计图。该图中，横截面速度向量由(vs,ws)表征；流向速度分量us则由等值面表征。从图中可以看出，在Ro<0.25时随着Ro增加，二次流逐渐增强；而当Ro>0.25时，二次流则随着Ro增加而减弱，在Ro=0.7时已经几乎无法看到明显的二次流了。这和前面的图3的结论基本一致。值得注意的是，当Ro=0.25和0.32时，二次流的大涡两侧存在次生的二次流。其中Ro=0.25时，次生二次流强于Ro=0.32时的结果。

(a) Ro=0

(b) Ro=0.01

(c) Ro=0.02

(d) Ro=0.25

(e) Ro=0.32

(f) Ro=0.7

2.2 弱旋转条件下的流动结构和统计

在上一节中，我们提到Rew=1300、Ro=0.02时平均流在槽道中心位置有局部回流现象。本部分将介绍在弱旋转Ro=0.02时的流动结构和统计分析。在Huang等人[20]的文章中，他们开展了一系列DNS研究，分析了网格分辨率、流向计算域、展向计算域和计算时间的影响。在计算结果中，他们发现，在Ro=0.02时，流动结构可能会随着计算时长而改变，不过这种计算结果还跟流向计算域有关。当流向计算域较小，例如Lx=4π时，展向Lz=4π的槽道中二次流的涡对(Roll-cells)数恒定为3对；而当Lx≥8π时，展向Lz=4π的槽道中会先出现3对涡对而后流动会演变为2对涡对。在流动处在3对涡对时，平均速度在中心区存在局部回流；而当流动处于2对涡对时，该局部回流消失。如果Lx=8π，Lz=6π时，该流态转换依旧存在，流动会从4对涡对转换成3对涡对。图5给出了该计算域下x=0位置处的中心线瞬时流向速度在展向和时间上的变化规律。可以看出，在长达5000个无量纲时间中，流动从层流加扰动初始解首先演化出4对流向涡对，该4对涡对的状态持续了大约1800无量纲时间(200～2000)。之后该流动的涡对开始出现合并，最后在3000个无量纲时间的时候，流动变成3对涡对结构。值得说明的是，我们尝试过计算超过10 000个无量纲时间，流动后期依旧保持在3对涡对的状态。

图5 在x=0处中心线流向瞬时速度信号u(0,0,z,t)随时间的变化关系[20]Fig.5 Contours of u(0,0,z,t) as time evolves[20]

图6给出了4对涡和3对涡结构下，流场的湍流统计规律。可以看出，两者之间是存在明显差别的。特别的，在4对涡时，平均速度在中心线处存在局部回流(从图6(d)中雷诺剪应力大于1可以推出来)；而当流动变成3对涡的时候，该平均速度的局部回流消失，雷诺剪应力的最大值小于1。

(a) 平均速度

和

(d) -〈u′v′〉+

需要说明的是，在前人的工作中，Lz=4π的槽道宽度中只发现了3对涡的状态[5,14-15]。我们认为这主要原因在于过去的计算时间不够长。在过去的计算中，最多只算了几百无量纲时间[5]，图5显示了计算时间如果不够长，容易忽略掉后面稳定的流态。因此，我们建议算湍流问题时，要用更长的计算时间。

2.3 多态现象及其统计结果、鲁棒性分析

事实上，RPCF中的多湍流态现象在文献[19]中已有论述。在该文章中，本文作者和合作者针对Rew=1300，Ro=0.2的RPCF流动开展了DNS模拟。由于在文献[5]的基础上，我们知道随着Ro从0.02增加到0.2时，二次流是会变得更强更稳定的；而在前期研究Ro=0.02时，我们也发现了Ro=0.02时存在流态转换过程。因此，在设计计算的时候，我们选取了两个不同的初始场，一个初始场是Ro=0.02时较前期的瞬时场，此时该流场拥有3对涡对结构(Lz=4π)；另一个初始场则选取Ro=0.02时后期的瞬时场，该流场拥有2对涡对结构。在计算中，两个初始场分别按照相同的程序，在相同的控制参数，相同的计算域和相同的网格分辨率条件下进行时间推进。经过超过3000个无量纲时间的计算，我们发现流动的最终统计结果和流动结构不一样。另外，我们还在流向和展向进行了网格加密，结果依旧出现多态现象。

图7给出了稳态条件下，两种不同流态下的二次流的横流速度矢量图和流向分量的等值线图。可以看出，两组不同的初始条件会等到不同的流动结构：一个会得到2对涡对结构；另一个则会保持3对涡对的状态。

(a) 2对涡对的结果

(b) 3对涡对的结果

两种不同流态下，其平均速度存在差异。为了说明问题，表2中给出了两组不同网格下，2对涡和3对涡的一些统计结果。从表中的数据可以看出：相同流态下，两组不同的网格下壁面摩擦雷诺数Reτ相差很小，误差都小于1‰。而不同流态之间的Reτ之间的差别则明显增大，差别约为6%。类似地，平均速度在中心线处的法向梯度也呈现类似的行为，只是误差相对较大。在3对涡流态下，两组不同网格的计算给出的值约为0.2；而在2对涡流态下对应的值减小到约0.14。从稳定性分析的角度看，前者在中心区趋向于中性稳定；而后者则处于稳定区域，这也可以解释为何2对涡条件下湍动能在中心区会相对较小。

表2 不同流态下的湍流特征统计[19]Table 2 Turbulent statistics at different states[19]

图8给出了2对涡和3对涡时的湍流二阶统计结果。可以看出，两种不同流态之间的统计存在非常明显的差异，而同一种流态在不同网格下的统计差异几乎可以忽略。在靠近壁面的区域，3对涡条件下的流向和展向速度脉动更强；同时，其对应的法向速度脉动和雷诺剪应力则在槽道中心区更强。

为了验证多态现象的存在是物理的，而不是计算域造成的虚假结果，我们在文献[19]中扩大了流向和展向的计算域，计算结果基本还能保持多态的结果。为了进一步验证多态现象的鲁棒性，我们还考察了多态现象对初始流场的依赖性。我们构造了一个新的流场ui：

(8)

更进一步，我们还模拟实验的方式[25]，在计算中逐渐增加旋转数或者逐渐减小旋转数，进而来研究多态现象的存在区间。在文献[21]中，Huang等人开展了两组DNS模拟。在相同的计算域和计算网格下，两组计算的初始场都是层流解加扰动的无散场。其中一组的旋转数Ro逐渐从0.01增加到0.02，0.03，0.04，0.05，0.1，0.2，0.3，0.5，0.6；而另一组则按照相同的Ro从0.6逐渐减小到0.01。在每个Ro上，都计算1250个无量纲时间，后一个算例(Ro)是以前一个算例的最后一个时刻的瞬时流场作为初始场。

(a) 〈u′u′〉

(b) 〈v′v′〉

(c) 〈w′w′〉

(d) -〈u′v′〉

图9 不同参数α下初场的最终流态示意图[22]Fig.9 A sketch of final states with different initial conditions defined by α[22]

图10给出了法向瞬时速度的前四个模态的幅值随计算时间的变化曲线。其中，模态幅值定义如下：

(9)

其中，Fz(f(z))表征的是函数f(z)的离散傅里叶变换，“〈〉x”表征流向平均。通常情况下，幅值最大的kz值可以视为涡对数。由图10可以看出，在0.03≤Ro≤0.3范围里，随着旋转数Ro逐渐增加，kz=2的模态逐渐占主导；而当Ro逐渐减小时，kz=3的模态逐渐占主导。

(a) 随着Ro增加的结果

(b) 随着Ro减小的结果

图11给出了在Ro在0.01到0.5之间的流向瞬时速度在x=0、y=0处的时间历程。其中，图11(a)给出的是随着Ro增加的情况，图11(b)给出的是随着Ro减小时的情况。可以看出在0.03≤Ro≤0.3范围里，随着Ro增加时，流场呈现2对涡；而随着Ro减小时，流场则有3对涡。图11的结果与Huisman等人[25]在Taylor-Couette流动的实验中的结果非常相似。

(a) Ro逐渐增加

(b) Ro逐渐减小

我们还分析了壁面摩擦雷诺数Reτ、体平均的湍动能以及体平均的二次流湍动能随着Ro增加和减小的结果。图12中给出了体平均的湍动能随着Ro增加和减小时的变化规律。由图可以看出，在0.03≤Ro≤0.3范围里，两组计算出来的体平均湍动能在相同的Ro下存在明显差别，3对涡时(Ro减小方向)的湍动能明显大于2对涡时的湍动能。在文章中，我们还发现Reτ和体平均的二次流的湍动能都存在类似于体平均湍动能一样的迟滞环[21]。

需要说明的是，在这两组计算中，10个Ro的具体数值并非预先有意设计的，而是在计算中随机选取的。计算时间间隔1250个无量纲时间也不是特意选取的。在所有计算结果出来之前，我们也不能保证我们能观测到迟滞环现象。最终结果我们发现了迟滞环，这至少说明了多态现象可以在很大的Ro范围里存在。我们承认如果改变Ro的具体值选取或者改变时间间隔，结果可能不一样。事实上，在Huisman等人的实验研究中，他们也是在一些情况下遇到了多态现象[25]。

图12 体平均的湍动能随着Ro增加和减小的变化规律[21]Fig.12 Changes of the volume-averaged turbulent kinetic energy as Ro increases or decreases[21]

2.4 多态条件下的能量流动分析

在Taylor-Couette湍流中多态现象的工作中，Huisman等人曾经猜测大尺度相干结构的选择性作用在多态现象中占据重要地位[24]。在RPCF中，二次流就是这种大尺度的流动结构。为了验证Huisman等人的假设，我们研究了二次流和剩余场之间的能量传输过程。参考Bech和Andersson[14]和Cai等人[26]的工作，我们将速度场分成四部分：

(u,v,w)=(〈u〉,〈v〉,〈w〉)+(us,0,0)+

(0,vs,ws)+(u″,v″,w″)

(10)

(11)

(12)

(13)

(a) 3对涡的流态

(b) 2对涡的流态

3 结论与展望

我们通过一系列直接数值模拟在Rew=1300的RPCF中发现了多湍流态现象。在不同的湍流态下其流动结构和湍流统计差别明显。该多态现象比较鲁棒，可以在非常大的Ro范围里被观察到。

我们还分析了不同湍流态下，能量在平均场、二次流场和剩余场之间的传递关系。通过数据分析我们发现，3对涡时(Lz=4π)，二次流较强，此时二次流相当于能量传输的一个桥梁：平均场将大部分能量传输给二次流，再由二次流传给剩余场；而当流动处在2对涡的流态时，此时二次流相对较弱，平均流虽然传输了很多能量给二次流，但是大部分都用于维持二次流本身，只有少部分能量能传输给剩余场。而剩余场的能量主要是通过平均场直接传输过来的。这种能量传输的不同表现，侧面验证了大尺度结构在多态现象中占据重要地位。

经典的湍流理论假设湍流是各态遍历的，是跟初始场无关的。这是湍流建模和模拟的基础之一。多态现象的出现将对湍流建模和模拟提出新的挑战。如果在模型的构建(生成项的大小都不一样)和验证中，将多态现象考虑进去是值得研究的课题。此外，多湍态现象在工程问题中是否存在？如果存在，它会对工程设计和安全有什么影响？这些也是非常值得研究的问题。