圆柱绕流中的雷诺应力弱化现象的讨论

A Review of Weakening of Reynolds Stress in Cylindrical Flow

摘要: 基于对圆柱绕流的直接数值仿真的数据的调查，我们讨论了k-ò模型对雷诺应力的预测与实际雷诺应力的差别，并发现在圆柱后方的雷诺应力湍流时有弱化的趋势，此外在渠道两翼侧存在应力强化的区域，描述这些趋势的比例系数没有表现出和墙壁距离的相关性。

Abstract: With the result of direct numerical simulation on cylindrical flow, it was discussed of the differences of Reynolds stress between that from prediction of k-epsilon model and that from realistic velocity of flow. It was found that the Reynolds stress tends to be weakened in the area behind the cylinder, but tend to be enhanced at the wing-side of the channel; the coefficient for description of the weak-ening didn’t show correlations to the distance to the wall.

1. 概述

在流体的数值仿真中，通常使用湍流模型来解决运算量过大的问题。基于RANS (Reynolds-Averaged Navier-Stokes Equations)的 $k\text{-}ϵ$ 模型是产业中最常用湍流模型。在RANS的方程中，出现了关于浮动量乘积平均的项，它通常被称为雷诺应力，是由小于网格尺度的速度的波动引起的。在 $k\text{-}ϵ$ 模型中，雷诺应力被假定为与牛顿粘性应力的形式相似，由宏观速度剪切率和湍流粘性来表示。

在之前的研究中 [1]，我们通过DNS (Direct Numerical Simulation)调查了圆柱绕流的统计特性，并且发现 $k\text{-}ϵ$ 模型无法再现DNS仿真结果中的带状速度分布和尾流随时间摆动的现象。此外我们还定义了一个湍流粘性期待值，与现 $k\text{-}ϵ$ 模型定义的湍流粘性进行了比较。我们发现，尾流中圆柱附近的地方根据 $k\text{-}ϵ$ 模型定义的湍流粘性被过大预测了。因此，如果我们导入某种方法表现出这种雷诺应力弱化的结果，就可以改善数值仿真的结果真实性。

尽管使用湍流粘性对雷诺应力建模的方法存在某些缺点，并且研究者提出了很多高阶精度的湍流模型，比如EDQNM (Eddy Damped Quasi-Normal Markovian) [2] 和LRA (Lagrangian Renormalization Approximation) [3]，产业中也有LES (Large Eddy Simulation)模型作为替代，但是使用湍流粘性的RANS因为其容易收敛、需要网格数量少等优点，仍旧是产业中最为常用的流体仿真方法。

在本文中，我们将先简要回顾之前的研究中所做的圆柱绕流的数值仿真，并导入两个比例系数用于对湍流应力变化的讨论。然后我们将观察这两个比例系数的分布，以及他们和墙壁、雷诺数之间的关系。最后则给出这项研究的发现和下一步研究的预想。

2. 模型与定义

2.1. 数值仿真

流场为长方体渠道，流场上游为圆柱形障碍物用于激发湍流，通过固定流入口和流出口的压力差进行了雷诺数分别为50，300，800的不可压缩流体的DNS，其控制方程为Navier-Stokes方程：

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

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

数值仿真的具体设置参考之前的研究 [1]，我们导入了拉格朗日系粒子用以统计湍流变量，这个方法的优点是可以避免对流项带来的复杂性，本研究中我们将继续使用这种方法取得的数据。

2.2. 雷诺应力与湍流粘性

通过雷诺分解将速度分解为平均速度和波动速度( $u_i=U_i+{u^{\prime }}_i$ )之后，可以发现平均速度的时间发展方程中出现了波动速度乘积的平均的项

${\tau }_{ij}=-\overline{{u^{\prime }}_i{u^{\prime }}_j}$,(3)

因为该项可以影响平均流的加速度，我们将其称为雷诺应力。我们对湍流模型研究的主要目标之一就是找到雷诺应力的合理的表现形式。 $k\text{-}ϵ$ 模型使用湍流粘性

${\nu }^{\text{turb}}:=C_{\mu }{k}^2/ϵ$, (4)

和宏观速度剪切率 $S_{ij}=\left({\partial }_i{U}_j+{\partial }_j{U}_i\right)/2$ 来表示雷诺应力 ${\tau }_{ij}:=2{\nu }^{\text{turb}}{S}_{ij}$，这里 $k=\overline{{u^{\prime }}_i{u^{\prime }}_i}/2$ 为湍流动能， $ϵ=\nu {\left({\partial }_i{u^{\prime }}_j+{\partial }_j{u^{\prime }}_i\right)}^2/2$ 为湍流动能耗散率，可以发现 ${\nu }^{\text{turb}}$ 为非负的场量。根据湍流的能量串跌理论，湍流中动能主要从大尺度传往小尺度，这意味着雷诺应力整体上起到让动能衰减的效果，而湍流粘性的方法能很好地满足总能量衰减的特性，并且从数值计算的角度考虑，湍流粘性模型具有较好的收敛性。

本研究中，我们定义了两个比例系数 $\alpha$ 和 $\beta$ 用来衡量实际的湍流粘性

${\nu }^{\text{turb}\ast }:=\left|\tau \right|/\left|2S\right|$, (5)

分别表示 $k\text{-}ϵ$ 模型对湍流粘性评价偏差的程度以及高湍流粘性区域的分布：

$\alpha :={\nu }^{\text{turb}\ast }/{\nu }^{\text{turb}}$, (6)

$\beta :=\log \left({\nu }^{\text{turb}\ast }/\nu \right)$. (7)

2.3. 湍流模型的改进方法

Lengani [4] 针对2维问题提出了一种解决方案，他使用湍流粘性 ${\nu }^{\text{turb}\ast }=\left(\begin{array}{ccc}{\nu }^{11}& 0& 0\\ 0& {\nu }^{22}& 0\\ 0& 0& 0\end{array}\right)$ 和旋转矩阵 $X\left(\theta \right)=\left(\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0& 0& 0\end{array}\right)$ 相结合的方法来表示雷诺应力

$\tau :=2{\nu }^{\text{turb}\ast }X\Gamma X^{\text{tr}}$, (8)

这里 $\Gamma$ 为 $S$ 的特征值构成的对角矩阵， ${\nu }^{\text{turb}\ast }$ 和 $\theta$ 由计测得到的大数据复归得到，右上角的tr表示矩阵的转置。其中 $\theta$ 所表示的角度为雷诺应力的特征向量的正主成分(对应的特征值最大的那一个)和平均速度剪切率的特征向量的正主成分的夹角。一个直观的联想是，取零成分(对应的特征值中间的那一个)的夹角为 $\psi$ 的话，就可以做出三维空间的变换矩阵用来解决三维问题。但这么做是没有意义的，我们可以发现Lengani的方法里 ${\nu }^{\text{turb}\ast }$ 和 $\theta$ 只依赖于欧拉坐标，因此他的数据是全时间平均的，尽管Z方向的雷诺应力不可避免的引发能量耗散，但是在他的模型里根据实际平均得到的 $S$ 不包含Z方向的任何成分，那么关于Z方向的角度变换和湍流粘性的扩展都是毫无意义的。

因此我们认为如果要改进模型，历史相关的数据是不可缺少的，也就是大数据应该取自拉格朗日系粒子，这样才能体现尾流随时间摆动的特点。此外Lengani的模型对几何有较强的限制，尽管可以导入几个自由度，比如雷诺数、障碍相对渠道宽度的比例、根据主流速度无量纲化的X位置等，但是一旦几何模型出现了变化，比如圆柱变成方柱、增加圆柱数量等，就必须重新建立大数据数据库。

我们提议的改进方案包括两点：一个是随时间衰减的某种标量，这种标量较大时会引起应力减弱，该势能的减弱的同时将会引起应力增强，该标量的通量则会引起雷诺应力和速度剪切率的夹角；而是这种标量的诱发机制为流域速度的强烈非连续性(比如平均速度的旋转张量)，衰减机制为流体的分子粘性。这些机制与固体力学的“损伤”概念类似但是有着微妙的不同，固体力学里的损伤可以通过能量守恒得到，流体的场景下湍流动能的本质并不是能量，雷诺应力的本质也不是应力。如果我们的模型可以与 $k\text{-}ϵ$ 模型相结合，他的方程组将是：

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

$\nabla \cdot U=0$, (10)

${\nu }^{\text{turb}\ast \ast }={\nu }^{\text{turb}\ast \ast }\left({\nu }^{\text{turb}},\nabla D_{\text{out}}\right)$, (11)

$\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-ϵ+D_{\text{out}}$, (12)

$\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$, (13)

$\partial D/\partial t+U\cdot \nabla D=D_{\text{in}}\left(\Omega \right)-D_{\text{out}}-D_{\text{diss}}\left(\nu \right)$,(14)

如果令 ${\nu }^{\text{turb}\ast \ast }:={\nu }^{\text{turb}}I$，(4)，(9)，(10)，(12)，(13)就是现有的 $k\text{-}ϵ$ 模型解非压缩性流体的湍流问题的方程组。在这个模型中，我们至少要找到三个式子定义式 ${\nu }^{\text{turb}\ast \ast }\left({\nu }^{\text{turb}},\nabla D_{\text{out}}\right)$， $D_{\text{in}}\left(\Omega \right)$ 和 $D_{\text{diss}}\left(\nu \right)$。本研究中我们只调查(6)和(7)作为建立模型的基础。

3. 结果

3.1. 比例系数的空间分布

雷诺应力 ${\tau }_{ij}$ 、宏观速度剪切率 $S_{ij}$ 、湍流动能k和湍流动能耗散率 $ϵ$ 的计算方法参考先前的研究 [1]。

将粒子的数据根据其所在的x，y坐标进行分组平均，由此算出映射在平面位置上的湍流粘性期待值 ${\nu }^{\text{turb}\ast }$ 和 $k\text{-}ϵ$ 模型定义的湍流粘性 ${\nu }^{\text{turb}}$ 的样本平均，并基于这两个量算出系数场 $\alpha \left(x,y\right)$ 和 $\beta \left(x,y\right)$，其空间分布如图1所示。系数 $\alpha$ 对于层流没有实在的物理意义，因此我们首先比较图1(b)和图1(c)。可以发现两者圆柱后部有 $\alpha \in \left[0,2\right]$ 的区域，且随着雷诺数的增加该区域向各个方向扩大。此外，关于 $y~\pm 0.15$ 的强应力区域，可以发现低雷诺数时只出现于上游，而高雷诺数时则是全流域的。对此一个可能的解释是，高雷诺数时相对于墙壁的高速度的主流引起的强速度剪切，有可能引起与圆柱障碍的扰动相类似的效果。此外我们还可以发现，墙壁附近的系数 $\alpha$ 也呈现较低的数值，这与 $k\text{-}ϵ$ 模型在墙壁附近的表现不准确的经验一致，在产业中通常使用单独的壁面函数 [5] 来解决这个问题。但是这种方法也存在着明显的缺点，如果墙壁不是平直的，那基于经验的对数法则的速度分布函数可能就会失效。因此我们提议一种基于边界条件引发的非连续性导致的扰动的统计规律的湍流模型，这种扰动与固体力学里的损伤类似，可以引起应力的弱化并且被流体的粘性摩擦所修复或者抵消。

接下来我们以左右对照的方式去观察图1，左图里的亮色区域代表 $k\text{-}ϵ$ 模型过小评价了湍流粘性，右图里的亮色区域则代表湍流粘性大幅度强于分子粘性。层流时在圆柱后方这种因为圆柱扰动而引发的湍流应力大约10倍于分子粘性应力，但是随着向下游移动则会逐渐减弱。另一方面， $k\text{-}ϵ$ 模型无法正确估算的 $y~\pm 0.15$ 区域，其实际的湍流应力并不比分子粘性应力高太多。对比图2(e)和图2(f)，可以发现随着雷诺数的增加，湍流应力较之于分子粘性应力的比例越来越大，雷诺数从300增加到800后，圆柱后方的倍率从约50提高了约200。这与我们对湍流的可以增强动量输送效果的经验是一致的。对比左侧高雷诺数时，系数 $\alpha$ 没有在 $x>0.3$ 的后方区域呈现递减的变化，系数 $\beta$ 仍旧体现出与x坐标相关的分布。对于这种现象可能的解释是，在扰动没有明显减弱前，雷诺应力的强度与 $k\text{-}ϵ$ 模型有较大差距，但是在远离圆柱障碍以后，这种扰动就快速失去了效果。

(a) (d) (b) (e) (c) (f)

Figure 1. Coefficient α and β Distribution in x-y plane, coefficient on the left α, the actual Reynolds stress intensity and $k\text{-}ϵ$. The difference between the predicted values of the model, the right side is the coefficient β, it represents the specific gravity of Reynolds stress and molecular viscous stress intensity. The Reynolds number distribution of the corresponding DNS from top to bottom is 50, 300, 800

图1. 系数α和β在x-y平面的分布，左侧为系数α，表现实际雷诺应力强度与 $k\text{-}ϵ$ 模型的预测值的差距，右侧为系数β，表现雷诺应力与分子粘性应力强度的比重。自上至下对应的DNS的雷诺数分布为50，300，800

3.2. 墙壁距离的影响

把粒子离最近的墙壁的距离定义为 $Y_{\text{wall}}$，根据各时刻的x坐标和跟墙壁间的距离取样本平均，两个系数 $\alpha$ 和 $\beta$ 的分布如图2所示。可以发现，在层流时，似乎可以发现与墙壁距离的某种相关性，特别是在上游，离墙壁越远， $\alpha$ 的数值越低， $\beta$ 的数值越高。但是在湍流时，我们并没有发现两者呈现出相关性，随机性引发的波动幅度远远大于 $Y_{\text{wall}}$ 的影响。因此我们认为这种由圆柱障碍的扰动引发的雷诺应力的弱化机制基本与墙壁距离无关，虽然在层流中两系数会体现出与墙壁距离的关联，但是对于流场为稳定状态的层流来说动态仿真不是必需的，可以适当加密网格来取得更加精确的结果，这并不需要增加太多的计算量。

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

Figure 2. Relations between coefficient α, β and the distances to the walls. The Reynolds number distribution of the corresponding DNS from top to bottom is 50, 800

图2. 系数α和β与墙壁距离之间的关系，自上至下对应的DNS的雷诺数分布为50，800

3.3. 雷诺数的影响

我们可以选择4个代表性的区域，两翼处应力强化发生的区域的上游A $\left(x\in \left[0.1,0.2\right],y\in \left[-0.16,-0.12\right]\right)$ 和下游B $\left(x\in \left[0.8,0.9\right],y\in \left[-0.16,-0.12\right]\right)$，流场中轴应力弱化发生的区域的上游C $\left(x\in \left[0.1,0.2\right],y\in \left[-0.02,0.02\right]\right)$ 和下游D $\left(x\in \left[0.8,0.9\right],y\in \left[-0.02,0.02\right]\right)$。不同雷诺数四个区域系数 $\alpha$ 的平均值记录在表1里，对比结果我们可以发现在湍流情况下C点的数值相对于层流下发生了明显的变化，D点数值虽然也体现出了和雷诺数的负相关，但是总体上我们可以判断在高雷诺数下，C点的效果的“寿命”更长，因此影响到了D点。A、B点的数值似乎和层流还是湍流没什么直接关系，可能是圆柱绕流这种几何下固有的特征。总体来说，800的雷诺数从湍流的角度来讲并不是特别高的数值，尽管从DNS或者水槽实验的角度来讲，进一步提高雷诺数有很大的难度，但是我们希望能从风洞试验中得到更有代表性的数据，比如雷诺数上万时，C点的 $\alpha$ 能否得到小于1的数值，A、B点的数值是否能仍旧维持在10附近。

Table 1. Averaged values of coefficient alpha with different Reynolds numbers and areas

表1. 系数α在不同雷诺数/区域下的平均值

4. 结论

在先行研究 [1] 中，我们通过直接数值仿真和基于 $k\text{-}ϵ$ 模型的湍流仿真的比较，发现湍流时实际出现的带状结构和各点速度的随时间变化的现象无法通过湍流模型再现。在本研究中，我们考察了湍流粘性基于 $k\text{-}ϵ$ 模型预测的湍流粘性与实际情况的差别。我们推测，在圆柱障碍后方雷诺应力发生了弱化，在计算里大部分区域的湍流粘性期待值都比 $k\text{-}ϵ$ 模型的预测值大10倍左右，而湍流粘性里的常数 $C_{\mu }$ 是根据雷诺数无限大时的渐近值推测出来的，因此我们需要更高雷诺数湍流下的数据来检验这个推测。使用墙壁函数是产业中传统的解决 $k\text{-}ϵ$ 模型在墙壁附近预测结果不准确的方法，但是我们发现描述湍流粘性期待值与预测值比例的系数似乎与墙壁距离没有表现出相关性。另一个发现是在圆柱两翼侧出现了似乎与湍流层流无关的应力强化区域，这可能是圆柱绕流几何条件产生的结果，我们好奇是否可以通过数学方法论证该区域的存在以及应力强化的程度。

如果我们把视角转回Navier-Stokes方程上，雷诺数无限大时的强湍流问题很容易让人联想到理想流体。但事实并非如此，在理想流体中角动量是守恒的，因此再强的宏观速度剪切率也无法引起湍流由涡旋构成的三维结构。可能会有人联想到固体力学的损伤问题中含有偶应力的本构模型是否可以用于描述湍流问题，例如本研究中的应力弱化现象可能描述为强烈非连续性引起的扰动分为湍流动能和湍流角动能，湍流动能可以直接参与雷诺应力引起的大小尺度的能量交换，因此在湍流角动能比例较大的初期(例如本研究中的圆柱后方)，因为雷诺应力的对角和为湍流动能，雷诺应力会表现的比预期值弱；随着扰动的扩散和湍流动能向小尺度各向同性波动的转化，湍流角动能的比例将会变小，此时则会出现像两翼那样雷诺应力比预期值强的效果。此外在研究蠕变流体问题时，有其他研究者使用偶应力本构模型的尝试 [5]，但是分子布朗运动引发的物理粘性与湍流模型中用于描述能量串跌的湍流粘性有着本质的不同，这也意味着我们如果要参考这种做法，需要进一步做更多的工作而非直接套用其数学形式到完全不同的具体问题上。

基金项目

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