圆柱绕流的湍流统计

Statistics of Turbulence in Cylindrical Flow

摘要: 为了调查现有的湍流模型无法再现圆柱绕流湍流时的尾流处周期性摆动的问题，我们对三维不可压缩流体的圆柱绕流进行了直接数值仿真，并使用无惯性拉格朗日系粒子来观察湍流统计量的规律，此外我们还提取圆柱壁面的时序数据用以讨论扰动的扩散问题。数值实验的结果显示：1) 高雷诺数情况下湍流动能和湍流动能耗散率会随时间发生剧烈震荡，这个震荡可能是圆柱壁面附近的扰动引起的；2) 湍流粘性模型要求的雷诺应力方向和速度剪切方向一致或者完全无关，但这个条件无法在圆柱绕流问题中得到满足。这些发现对我们将来改进湍流模型以得到更贴近实际的仿真结果可能有所帮助。

Abstract: To investigate why the current turbulent models cannot represent shaking behaviors in the wake of turbulent cylindrical flow, the three-dimensional direct numerical simulations of incompressible cylindrical flow were performed. We used Lagrangian inertialess particles to collect the information of turbulent fields, and discussed the features of diffusion of disturbs near the cylindrical walls. It was found via the numerical experiments that: 1) Fluctuations of turbulent kinetic energy and dissipation rate of kinetic energy occurred in the cases of higher Reynolds numbers; this behavior may be caused by the disturbs from the walls. 2) It was expected that directions of Reynolds stress and macro-scale strain of fluid velocity are aligned or independent, but this precondition is not satisfied in cylindrical flow. These results may help us to improve the turbulent models in the future.

1. 概述

圆柱绕流可以用来描述很多实际中的场景，比如海上平台、桥墩、风力发电机的塔筒和废弃原子能燃料棒的储存池。实验和数值模拟被用于研究其规律，特别是高雷诺数时的尾流的性质。通过风洞实验、数值仿真和真实环境测量 [1] [2] [3] [4]，进行了很多与之相关的研究。随着雷诺数的增加，流体的运动从层流迁移到湍流，对于圆柱尾流，研究 [5] 显示，临界值为230~260，Re > 260时会产生不连续性，流场构造从2维迁移到3维。

层流与湍流的一个巨大的差别可以体现在高次力矩上，例如对于扰动为高斯分布的低频带剪切力的各向同性湍流来说，研究 [6] [7] [8] 显示，在向湍流迁移前 $\langle {\phi }^{2n}\rangle /{\langle {\phi }^n\rangle }^2$ 为常数(高斯分布的性质)，这里 $\phi \left(x,t\right)$ 为某一规格化的常量，例如 $ϵ/{ϵ}_{rms}$，随着雷诺数的增加，迁移为湍流后将趋近于雷诺数的幂函数。这是湍流的间歇性的表现之一，间歇性在时序数据上则体现为尖刺形状的极端值的出现。

因为没有找到Navier-Stokes方程在一般场景下的解析解，湍流场景的力学仿真是一个尚未被完全解决的问题。因为DNS (Direct Numerical Simulation，直接数值仿真)需要极细的网格，产业中的流体仿真一般使用湍流模型，比如RANS (Reynolds Averaged Navier-Stokes，雷诺平均纳维尔斯托克方程)模型和SGS-LES (SubGrid Scaled Large Eddy Simulation，子网格尺度大涡模拟)模型。其中 $k\text{-}ϵ$ 模型 [9] [10] 是一种最常用的RANS模型，但是如果我们尝试用 $k\text{-}ϵ$ 模型对圆柱绕流进行数值模拟，如图3(d)，会发现湍流时本应存在的尾流的随时间摆动(在瞬时的流场构造上，则呈现为压力和速度大小的带状的纹理结构，例如图3(b)和图3(c)现象并不能呈现出来。因为这种差别出现在宏观尺度上，我们认为现有的湍流模型仍有改进的必要。

我们将在第二节介绍物理模型和控制方程以及一些常用的湍流统计中常用的量的定义；第三节会介绍数值仿真的设定和格式以及算法的选择；第四节我们会介绍包括基础物理量和自定义的统计量在内的仿真结果；第五节是总结和展望。

2. 模型和控制方程

计算区域为如图1所示的长方体，对称中轴上游有圆柱形的障碍，长方体的前端为流入口，后端为流出口，其他侧面和圆柱表面为无滑动墙壁，坐标原点为圆柱中心的底部。流入口和流出口的边界条件为固定压力，且流入口的压强比流出口高ΔP以提供动力，流入流出口的切向速度梯度为零；墙壁上的速度为零，压力梯度为零。初始值为各点具有一定的速度 $U_0$。

不可压缩牛顿流体的力学现象可以用Navier-Stokes方程描述，

$\partial u/\partial t+u\cdot \nabla u=-\nabla p+\nu {\nabla }^{2}u$ (1)

Figure 1. Schematic diagram of the calculation area, where X1 = 0.2, X2 = 1.2, R = 0.02, Y = 0.4, Z = 0.2, and the fluid is water

图1. 计算区域的示意图，其中X1 = 0.2，X2 = 1.2，R = 0.02，Y = 0.4，Z = 0.2，流质为水

$\nabla \cdot u=0$ (2)

这里 $u$ 为流体的运动速度， $p=\left(P-\bar{P}\right)/\rho$ 为使用密度规格化后的压强， $\rho$ 为流体的密度， $\nu$ 为动力粘性系数。

流体速度可以通过雷诺分解表示为平均量和变动量之和 $u=\bar{u}+u^{\prime }$，有 ${\bar{u}}^{\prime }=0$， $\overline{\stackrel{¯}{u}}=\bar{u}$。为了记述方便我们通常写成 $U=u$，可以得到雷诺平均速度的控制方程。如果导入湍流粘性，使得 ${\tau }^{\ast }=-\overline{u^{\prime }{u}^{\prime }}={\nu }^{\text{turb}}\nabla U$，并且应用 $k\text{-}ϵ$ 模型湍流动能 $k=\left(u^{\prime }\cdot u^{\prime }\right)/2$，和湍流动能耗散率 $ϵ=\nu {\left(\nabla u^{\prime }+{\left(\nabla u^{\prime }\right)}^{\text{T}}\right)}^2/2$ 来定义湍流粘性 ${\nu }^{\text{turb}}=C_{\mu }{k}^2/ϵ$， $C_{\mu }=0.09$，可以得到工业应用里经常被使用的湍流仿真用的方程组：

$\partial U/\partial t+U\cdot \nabla U=-\nabla \bar{p}+\nabla \cdot \left[\left(\nu +{\nu }^{\text{turb}}\right)\nabla U\right]$ (3)

$\nabla \cdot U=0$ (4)

$\partial k/\partial t+U\cdot \nabla k=\nabla \cdot \left[\left(\nu +{\nu }^{\text{turb}}/{\sigma }_k\right)\nabla k\right]+{\nu }^{\text{turb}}{S}^2-ϵ$ (5)

$\partial ϵ/\partial t+U\cdot \nabla ϵ=\nabla \cdot \left[\left(\nu +{\nu }^{\text{turb}}/{\sigma }_{ϵ}\right)\nabla ϵ\right]+C_1ϵ{\nu }^{\text{turb}}{S}^2/k-C_2{ϵ}^2/k$ (6)

$S=\left(\nabla U+{\left(\nabla U\right)}^{\text{T}}\right)/2$ (7)

$S$ 通常被称为宏观速度剪切率，根据经验取得的常数的值为 $C_1=1.44$， $C_2=1.92$， ${\sigma }_k=1.0$， ${\sigma }_{ϵ}=1.3$。

因为本研究的目的是寻找 $k\text{-}ϵ$ 模型无法再现尾流摆动现象的原因和解决办法，所以我们将使用DNS (方程组(1) (2))进行流场的数值仿真，然后观察与 $k\text{-}ϵ$ 模型相对应的变量的特征。

此外雷诺应力 ${\tau }^{\text{*}}={\tau }_{ij}^{\text{*}}$ 其中的迹成分可以被连续性方程相关联的压力项抵消，因此，我们只计算除此之外的部分 ${\tau }_{ij}={\tau }_{ij}^{\text{*}}-{\tau }_{ii}^{\text{*}}{\delta }_{ij}/3$，这里我们使用了爱因斯坦记号 ${\tau }_{ii}^{\text{*}}={\displaystyle {\sum }_{i=1}^3{\tau }_{ii}^{\text{*}}}$。此外对于对称2阶张量 $S$ 和 $\tau$，都可以使用特征值分解： $A=P{\left(\begin{array}{ccc}{\lambda }_1& {\lambda }_2& {\lambda }_3\end{array}\right)}^{\text{T}}I$，这里 $P$ 为由特征向量构成的旋转矩阵，可以发现模量满足 $\left|A\right|=\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}$，特征向量分解将协助我们讨论湍流粘性模型是否妥当以及如何改进该模型的工作。

3. 数值模拟

我们运行了三个数值计算算例，其雷诺数分别为约50，300，800。其中Run 1和Run 2使用相同的网格，包含288，640个控制体，Run 3使用进一步加密的网格，包含1,600,000个控制体。Run 1和Run 2所用网格的x-y断面如图2所示。此外，我们还使用包含10,800个控制体的较粗的网格，基于 $k\text{-}ϵ$模型进行了数值仿真，压力边界条件与Run 3相同，湍流变量流入条件为湍流强度等于0.05，其结果表示在图3(d)里。

Figure 2. Mesh used for numerical simulation, containing 288,640 control volumes, 309,246 vertices and 886,232 mesh faces

图2. 数值模拟所用的网格，包含288,640个控制体，309,246个顶点和886,232个网格面

(a) (b) (c) (d)

Figure 3. Velocity and pressure distribution of the flow field, the Reynolds numbers from top to bottom are about 50, 300, 800, 800, respectively, the first three are direct numerical simulation results, and the fourth is RANS simulation results. Red represents areas of higher pressure; blue represents areas of lower pressure; arrows represent the (x, y) component of fluid velocity

图3. 流场的速度和压力分布，自上至下的雷诺数分别为约50，300，800，800，前三个为直接数值仿真结果，第四个为RANS仿真结果。红色代表压力较高的区域，蓝色为压力较低的区域，箭头代表流体速度的(x，y)成分

仿真程序为修改版的OpenFOAM，在v2112基础上加上了对无惯性粒子进行处理和统计，以及在特定网格输出每个时间步数据的代码。时间积分的格式为欧拉前进迭代，通量计算格式为MUSCL (Monotonic Upwind Scheme for Conservation Laws) [11]，扩散项计算为二阶线性差分，压力-速度耦合求解器为PISO (Pressure Implicit with Split Operator)，泊松方程的求解器为GAMG (Geometric-Algebraic Multi-Grid) [12]。计算设备为联想ThinkStation P620工作站(CPU: AMD Ryzen™ Threadripper™ PRO 3955WX, RAM:256GiB)，操作系统为Fedora 35，计算中基于区域分割的16进程MPI并行运算。雷诺数通过流入口和流出口的压力差ΔP控制，分别为1000 Pa，5000 Pa和20,000 Pa。流质为水，密度1000 kg∙m⁻³并不会出现在程序计算中，通过动粘性和密度规格化压力的方法体现数值仿真里。在计算到达稳定状态之前的数据则会被去掉，统计计算只使用稳定状态的数据。可视化处理则使用了Paraview和matplotlib。

4. 计算结果

4.1. 流场整体

$z=Z/2$ 平面的某瞬间的流体速度和压力的分布在图3中显示，可以发现湍流里存在着由类周期性的波动所形成的带状结构，而层流的速度和压力分布几乎不随时间发生变化。

4.2. 拉格朗日系粒子

比较方程(3)，(5)和(6)可以发现，导入流场自变量 $\theta \left(x,t\right)$ 用于湍流建模时，其支配方程的左边一般是 $\left(\partial /\partial t+u\cdot \nabla \right)\theta$ 的形式，一般这被称为物质导数，又写作 $\text{D}\theta /\text{D}t$，可以被理解为跟随某个流体小微团，对其各物理量的时间变化描述。这种描述最大的好处是只需要关注其右端项，即扩散项和源项的影响，而忽视对流运动的影响。因此在我们的研究中导入了拉格朗日系粒子，计算该粒子的各时刻的湍流统计量用以后面的分析和研究。

假设存在随着流体运动的无惯性的粒子，其瞬时速度为所在位置的流场速度，即 $v\left(x_{\xi },t\right)={\displaystyle \int u\left(x^{\prime },t\right)\delta \left(x^{\prime }-x_{\xi }\right)\text{d}{x}^{\prime }\text{d}t}$，编号为 $\xi$ 的粒子的位置为速度的积分 $x_{\xi }\left(t\right)={\displaystyle {\int }_0^t{v}_{\xi }\left(t^{\prime }\right)\text{d}{t}^{\prime }}$。计算中我们使用高斯过滤器作为插值函数 $\delta \left(r\right)={\left(r_0\sqrt{\pi }\right)}^{-1}{\text{e}}^{-r^2/r_0^2}$，使用2次Runge-Kutta法进行位移的时间积分。

对于粒子的湍流变量，我们使用粒子附近半径为R的球体内网格点数值的样本平均求得，即 $\Theta \left(\xi ,t,R\right)=N^{-1}{\displaystyle {\sum }_{n^{\prime }}\theta \left(x_{n^{\prime }},t\right)H\left(\left|x_{n^{\prime }}-x_{\xi }\right|-R\right)}$，这里 $n^{\prime }$ 为网格的编号， $N={\displaystyle {\sum }_{n^{\prime }}H\left(\left|x_{n^{\prime }}-x_{\xi }\right|-R\right)}$ 为总样本数，由此可以算得平均速度 $U\left(R,\xi ,t\right)$ 和扰动速度 $u^{\prime }\left(\xi ,t,R,n^{\prime }\right)$。由此我们可以定义单个粒子周围半径R的球体内的湍流动能、雷诺应力、宏观速度剪切率和动能耗散率为

${\tau }_{ij}^{*}\left(\xi ,t,R\right)=N^{-1}{\displaystyle {\sum }_{n^{\prime }}{u^{\prime }}_i\left(\xi ,t,R,n^{\prime }\right){u^{\prime }}_j\left(\xi ,t,R,n^{\prime }\right)H\left(\left|x_{n^{\prime }}-x_{\xi }\right|-R\right)}$ (8)

$k\left(\xi ,t,R\right)={\displaystyle {\sum }_{i=1}^3{\tau }_{ii}^{*}/2}$ (9)

${\tau }_{ij}\left(\xi ,t,R\right)={\tau }_{ij}^{*}-2k{\delta }_{ij}/3$ (10)

$S_{ij}\left(\xi ,t,R\right)=\left({\partial }_j{U}_i+{\partial }_i{U}_j\right)/2$ (11)

$ϵ\left(\xi ,t,R\right)={\left(2N\right)}^{-1}{\displaystyle {\sum }_{n^{\prime }}\nu {\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)}^{2}H\left(\left|x_{n^{\prime }}-x_{\xi }\right|-R\right)}$. (12)

湍流动能的数据在图4里表示。可以发现，湍流对比于层流，不同起始时刻生产的粒子的湍流动能表现出较大范围的变动，这可以解释为湍流动能在流场里的分布并不是均匀的。但是另一方面，变动的幅度在上下游是接近的，因此可以推测在湍流耗散区域，随机性的效果远小于湍流增幅的区域。另一个有趣的发现是近墙壁的湍流动能最强的区域，可以发现雷诺数越高，其位置越处于下游，这可能是来自流场中轴的扰动的影响传播到了墙壁侧。

Figure 4. The turbulent kinetic energy of particles moving to each X coordinate. The starting position of the left image is the near-axis (0.07, 0.02, 0), and the starting position of the right image is (0.07, 0.18, 0). Curves of different colors represent randomly selected 20 particles generated at different times. The top-down examples are laminar flow (Re~50), low Reynolds number (Re~300), and high Reynolds number (Re~800). The average radius is R = 0.05 meters, and the average time T is 0.005 seconds

图4. 粒子移动到各X坐标时的湍流动能，左图起始位置为近中轴(0.07，0.02，0)，右图起始位置(0.07，0.18，0)。不同颜色的曲线表示随机选择的20个不同时刻生成的粒子。自上而下的算例分别为层流(Re~50)，低雷诺数(Re~300)和高雷诺数(Re~800)。平均半径为R = 0.05米，平均时间T为0.005秒

湍流动能耗散率的数据在图5里表示。与湍流动能的分布最大的差别是在高雷诺数的情况下，中心流域的粒子的湍流动能耗散率先是随着流动快速减小到约 ${ϵ}_0/5$ 之后开始缓慢减小(这里我们把初始位置的湍流动能耗散率称作 ${ϵ}_0$ )。而墙壁附近的粒子的 $ϵ$ 是随着流动平缓减小的。

Figure 5. Turbulent kinetic energy dissipation rates when particles move to each X coordinate; the layout of the picture is the same as in Figure 4

图5. 粒子移动到各X坐标时的湍流动能耗散率，图片的布局和图4相同

如果使用湍流粘性模型，那么 $S_{ij}$ 和 ${\tau }_{ij}$ 将被期待具有共同的主方向。这在墙壁附近是很明显不成立的，通常的经验显示两者呈45度左右的夹角，因为 $S_{12}$ 和 ${\tau }_{11}$ 远大于其他成分。因此我们也对两者夹角进行了统计，观察是否两者存在对齐的关系。通过图6可以发现两者的夹角在45˚/135˚附近的概率较大，特别是近墙壁的场景，平均值大约为90度，两者几乎不会重合或呈现相反的方向。雷诺数的影响主要在于中心区域产生的粒子的夹角分布范围更广，但是仍然并非完全随机而是在90度±45度附近。

Figure 6. The angle between turbulent stress and velocity shear rate when particles move to each X coordinate; the layout of the picture is the same as that in Figure 4

图6. 粒子移动到各X坐标时的湍流应力和速度剪切率的夹角，图片的布局和图4相同

一个简单的扩展粘性模型的方法是引入旋转矩阵，即 $\tau =2{\nu }^{\text{turb}}RS$，其中 $\left|R\right|=1$。因此我们也对RANS模型里常用的湍流粘性模型 ${\nu }^{\text{turb}}=C_{\mu }{k}^2{ϵ}^{-1}$ 和两者的模的比值 ${\nu }^{\text{turb}\ast }=0.5\left|\tau \right|/\left|S\right|$ 进行了统计，后者可以视为改进模型后的湍流粘性的期待值。

由图7和图8的比较可以看出，随着雷诺数的增加，两者的大小趋近，但是仍旧有一定的差别：1) 期待值的近墙壁的湍流粘性随流动变化很小，但是 $k\text{-}ϵ$ 模型的结果和湍流动能类似变化幅度很大。2) 期待值随着粒子不同生成时间而产生的变化幅度要比湍流动能的变化小很多，这可能是因为 $k\text{-}ϵ$ 模型的结果基于自变量k和 $ϵ$， ${\nu }^{\text{turb}}\propto k^2$ 且k的变化较为剧烈，导致 ${\nu }^{\text{turb}}$ 表现出很大的变化范围。

Figure 7. The turbulent viscosity calculated based on the $k\text{-}ϵ$ model when the particles move to each X coordinate; the layout of the picture is the same as that in Figure 4

图7. 粒子移动到各X坐标时基于 $k\text{-}ϵ$模型算出的湍流粘性，图片的布局和图4相同

Figure 8. Turbulent viscosity inferred based on turbulent stress and velocity shear when particles move to each X coordinate; The layout of the picture is the same as in Figure 4

图8. 粒子移动到各X坐标时基于湍流应力和速度剪切推断出的湍流粘性，图片的布局和图4相同

通过图4~8可以发现，层流与湍流最大的差别在于：层流下的各种场量在不同时刻的相同位置的数值几乎是相同的，而湍流的场合下，同一个统计量在相同位置的则会随着时间有着大范围的波动。通过图3可以看出压力的分布在湍流的场景下不是随着X方向单调减少，而是高低交替的带状构造的性质也与这一点保持一致。但是如果我们使用现有的湍流模型比如 $k\text{-}ϵ$ 模型进行高雷诺数仿真，并不能再现这种带状构造，而是和低雷诺数相似呈现出维持稳定状态的结果(图3(d))。根据不可压缩流体的压力与速度梯度的关系 $p=-{\nabla }^{-2}\left(\nabla u:\nabla u\right)$，我们可以推断出宏观的速度剪切率也是呈现这样的带状构造。如果从拉格朗日坐标系的视角，我们可能会猜想在湍流状态下，扰动的发生时，因为时刻不同而对每个流体微团赋予了不同的湍流强度，这意味着在图4和图5中，我们应该发现很多平行线。然而中轴附近的粒子并没有表现出这种形态，而是杂乱无章的交错在一起，墙壁附近的粒子则符合这种猜想，强者多半恒强，弱者多半恒弱。中间区域出现这种情况的可能原因之一是，对于中间区域，发生于圆柱两端(图4中的深蓝色区域)的扰动的后续区域会随机地产生涡旋的碰撞导致了湍流动能的交换，而不是严格的遵循轴对称由较近的其中一方的扰动源主导。

我们着重观察的一点是宏观速度剪切率S和雷诺应力 $\tau$ 是否存在某种对齐关系或者完全独立，这对于湍流粘性模型是否需要扩展以及如何扩展至关重要。在低Re时，两者的主方向很明显的呈现45度或者135度夹角，并随机切换，在高雷诺数时，这个稳定状态下的夹角分布于更广泛的范围，但是无论如何，并不是均匀的分布，而是更大概率在30~60度或者120~150度才能维持稳定。这并不是一个意外的结果，因为我们可以很明显的发现宏观剪切率的集中于 $\partial U_x/\partial y$ 这个分量，这是由渠道流的无滑动墙壁的边界条件决定的，但是雷诺应力中最大的分量应该是 ${\tau }_{xx}$，连续性条件使得速度扰动集中在了X方向。因此其中一种可能湍流雷诺应力的本构模型的改进方案，是导入边界条件(流域几何形状)的影响。

此外我们还调查了计算平均量所用的球体半径的影响(图9)，这里我们只比较Run 3 (Re ~ 800)的结果。可以发现墙壁附近和流域中轴附近的曲线的性质有很大区别，总体来说对于湍流建模重要的参考量 $\left|\tau \right|/\left|2S\right|$ 来说，流域中轴附近其大小大约与R²成反比，墙壁附近则是随着 增加趋近于某个极限值。

Figure 9. The relationship between the turbulence variable and the average radius when the particles move to each X coordinate, from top to bottom are the turbulent kinetic energy k, the turbulent kinetic energy dissipation rate $ϵ$, the model-based turbulent viscosity ${\nu }^{\text{turb}}$ and the modulus ratio-based turbulent viscosity ${\nu }^{\text{turb}\ast }$, the different colored lines represent the average radius. The left side is near the axis; the right side is near the wall

图9. 粒子移动到各X坐标时湍流变量和平均半径的关系，自上至下为湍流动能k，湍流动能耗散率 $ϵ$，基于模型的湍流粘性 ${\nu }^{\text{turb}}$ 和基于模量比的湍流粘性 ${\nu }^{\text{turb}\ast }$，不同颜色的线表示平均半径。左侧为近轴心，右侧为近墙壁

4.3. 圆柱附近的时序数据

假定湍流的产生是因为流体微观上的不连续，那么圆柱附件的点可能含有描述这种不连续产生的时候的信息，为此我们选择了8个网格(图10)，尝试根据他们的时序信息进行分析。其中2号点为圆柱附近的点，流体向X方向流动时会因为边界条件产生扰动，并流向3号点；4号点则受扩散项的影响获得扰动，1号项则可能通过压力的影响获得扰动。作为对比，6号点因为受对流影响极小，5号点、7号和8号点被期待与4号点具有共同的机制。

首先我们对不同雷诺数下的2号和6号点的主流方向速度和截面方向速度进行比较，为了方便观察，使用了绝对值的平均数进行了归一化。由图11可知，速度的变化幅度和变化频率随着雷诺数的变化增加而增加。特别是对于2号点来说，中等雷诺数时速度的变化体现出了明显的间歇性。

对于2号点来说，由 $\partial U_x/\partial y$ 所激发 ${\omega }_z$ 很容易被联想为尾流的周期型摆动来源，连同附近3点1、3、4号的数据如图12，可以发现对于较高雷诺数的流场，1号点和3号点产生了很明显的类周期的波动，他们的频率与2号点类似，但4号点几乎看不到这种波动。由此可以推测这种类周期性的行为主要通过压力而非粘性来扩散到四周，并且会被边界条件所限制。速度剪切率的性质与其类似，在图13里表示。

Figure 10. Top view of the center position of the grid for obtaining timing information, the position in the Z direction is the positive cut surface z = Z/2

图10. 获取时序信息的网格的中心位置的俯视图，Z方向的位置为正中断面z = Z/2

Figure 11. Variation of fluid velocity with time at a point near the cylinder. The left picture shows the wing (position 2) and the right picture shows the tail (position 6). The upper figure shows the main flow direction component of the velocity (u_x(t)), and the lower figure shows the cross-sectional direction of the velocity as (u_y(t)), both of which are normalized using the average of absolute values

图11. 圆柱附近点流体速度随时间的变化。左侧图为翼侧(位置2)，右侧图为尾侧(位置6)。上图为速度的主流方向成分(u_x(t))，下图为速度的截面方向成为(u_y(t))，两者均使用绝对值的平均进行规格化

Figure 12. Z-component of the vorticity of the cylinder flanks as a function of time. The upper left is the comparison of time series changes of different Reynolds numbers at position 2; the upper right, lower left and lower right are the changes of three Reynolds numbers (50, 300, 800) at position 2 and nearby points (1~4) in tiny time segments, respectively

图12. 圆柱侧翼的涡度的Z成分随时间的变化。左上为位置2不同雷诺数时序变化的比较；右上，左下和右下分别为三种雷诺数(50，300，800)位置2及附近几点(1~4)在微小时间片段的变化

Figure 13. X-Y components of cylindrical flank velocity shear rate as a function of time. The upper left is the comparison of the time series changes of different Reynolds numbers at position 2; the upper right, lower left and lower right are the changes of three Reynolds numbers (50, 300, 800) at position 2 and nearby points (1~4) in tiny time segments, respectively

图13. 圆柱侧翼速度剪切率的X-Y成分随时间的变化。左上为位置2不同雷诺数时序变化的比较；右上，左下和右下分别为三种雷诺数(50，300，800)位置2及附近几点(1~4)在微小时间片段的变化

对比2号点，6号点似乎更缺少由边界条件导致的抑制扰动的能力，5号和7号点的角速度似乎是背景角速度与类周期波动效果的叠加，而这种类周期性在速度剪切率里似乎也存在(图14)。对比2号点和6号点，我们认为2号点的波动更有规律性，6号点可能不是这类周期型波动的主要来源。

Figure 14. X-Y components of vortex and velocity shear rate at the tail of the cylinder as a function of time. The left side is the comparison of time series changes of different Reynolds numbers at position 6; the right side is the changes of Run 3 (Re~800) at position 6 and nearby points (5~8) in a small time segment

图14. 圆柱尾侧涡旋和速度剪切率的X-Y成分随时间的变化。左侧为位置6不同雷诺数时序变化的比较；右侧为Run 3 (Re~800)位置6及附近几点(5~8)在微小时间片段的变化

通过对对扰动源附近的流动状态进行调查，首先是我们确认了湍流状态下间歇性的存在(图11)，可以看到层流下是没有任何尖刺形状的波动的，而在湍流状态它会不定期的出现，振幅远大于其他时候的波动。此外我们发现扰动几乎不向Y方向传播，但是扰动可以向流动的上游传播，位置1，2，3 (图10)都可以观察到类似于周期性的波动，位置4则保持平稳。这意味着压力可能是湍流扰动传播的关键，因为如果是通过粘性传播的话4号点的状态应该和1、3号点相似，而如果是对流传播的话，则上游的位置1不应该体现出波动。

5. 总结和展望

通过数值仿真和一系列的统计方法，我们比较了圆柱绕流场景下层流与湍流以及近中轴与近墙壁流域的湍流统计变量的性质。我们认为湍流时出现的类似周期性的波动是圆柱绕流不可忽略的特性，现有的湍流的模型需要改进使得数值仿真时可以再现这种随时间变化的流场以及带状的空间结构。为了调查这种差别的原因，我们使用了拉格朗日系粒子和单点时序信息两种数据进行分析。数据显示湍流与层流最大的差别在于：湍流时不同时刻同一位置的湍流变量会在较大范围内波动，而层流时他们的数值只与位置有关。此外我们还对比了等效于湍流粘性的湍流应力与速度剪切率的模量比值和RANS模型算出的湍流粘性的差异，结果显示RANS模型算出的湍流粘性呈现出更大的波动范围。在时序数据方面，湍流速度的间歇性得以确认，这个现象没有出现在层流数据中。此外根据时序数据采集点位置间的关系，我们推测波动的扩散主要通过压力传播，而非通过对流效果或者分子粘性传播。根据这些结果，我们认为以参考固体力学的损伤模型来改进湍流应力的本构模型可能会改善仿真的贴近真实度(比如再现压力和速度场的带状的空间结构)，这将是我们今后的研究方向。

基金项目

中国博士后科学基金，课题编号2021CA020302。