基于LES的线性分层流尾迹分析

Large-Eddy Simulation Based Linear Stratified Flow Wake Analysis

摘要: 采用LES (Large-eddy Simulation，大涡模拟)研究亚临界雷诺数Re = 3700下，分层效应对球后尾迹的影响。选取了Fr = 0.05、∞，利用UDF (User Defined Function)建立密度线性分层分布，从运动球体的尾迹特征、漩涡结构、速度分布、等分析分层的存在对尾迹的影响。研究表明：分层的存在导致尾迹垂向抑制和速度分布各向异性增强；同时分层对涡结构和涡脱落也有重要的影响。

Abstract: Large-eddy Simulation (LES) was used to study the effect of stratification on the wake of the sphere at subcritical Reynolds number Re = 3700. Fr = 0.05 and ∞ were selected, and the linear hierarchical distribution of density was established by using the User Defined Function (UDF). The influence of the existence of stratification on the wake was analyzed from the wake characteristics, vortex structure, velocity distribution, and so on. The results show that stratification leads to vertical suppression of wake and enhanced anisotropy of velocity distribution. Stratification also has important effects on vortex structure and vortex shedding.

1. 引言

分层尾迹广泛存在于海洋 [1] 、河流、大气尾流 [2] 中，流体受到温度、盐度、压力影响，在垂直方向上形成密度分层结构。当稳定的密度分层结构受到外部扰动或内部流动破坏时，运动体产生的尾迹传播形成内波，由于产生的内波持续时间长，在军事上可通过对海洋潜体产生的内波进行反侦察，进而得知潜体的位置、航速等信息，分层流中尾迹分层流尾迹因此受到学术界的广泛关注。早期对于分层流的研究以实验为主，Lin [3] 、Chomaz [4] 、Boyer [5] 等人分别对圆球、圆柱在分层流形成的尾迹进行可视化分析，发现Fr和Re是影响分层流尾迹的重要因素。

随着计算流体力学的发展，数值模拟被大量运用到分层流尾迹流场的深入探究中，对于分层湍流尾迹的研究越来越多。Gourlay [6] 采用DNS对Re = 10,000、Fr = 10的分层流尾迹涡结构分析，发现在远尾流的煎饼涡流最终在没有相干结构的前提下形成。Pasquetti [7] 采用LES对Re = 10,000，Fr = 25的球体进行空间和时间演化，成功模拟出分层流球体尾迹演化，分别包括3D，NEQ，Q2D。Ortiz Tarin [8] 研究了Re = 104，Fr = 0.5、1、3、∞下椭球绕流，分析分层对长形球体流动的影响，相对于钝体，流过细长体的气流表现出更强的浮力效应。Diamessis [9] 对比了Re = 5000、105下Fr = 4、64分层流尾迹演化，发现尾涡结构演化与Re有关，较高的Re下，三维湍流向准二维、以浮力为主的流动调整的NEQ状态持续的时间明显长于较低的Re。

近年来对于湍流分层流尾迹大多是基于高Re数、高Fr数下的远尾迹研究，对于强分层效应(低Fr)较低Re数下的尾迹涡结构描述较少。Rodriguez [10] 采用DNS对Re = 3700下球体绕流进行数值模拟，发现Re = 3700时，在尾迹不同位置中识别出3种不同的频率，分别对应着：大尺度涡脱落频率、来自剪切层的小尺度KH不稳定性频率、再循环区收缩扩大引起的不稳定性频率。

本文主要研究分层效应对球后尾迹的影响，采用LES对Re = 3700、Fr = 0.05、∞球后尾迹进行数值模拟，分析球后尾迹速度场、漩涡结果，阐述强分层对湍流尾迹的影响以及恢复作用力如何作用尾迹的机理，为分层湍流近尾迹特征的认识提供有益的补充。

2. 数值方法与模型

2.1. 控制方程

分层流中，运动物体的尾迹流场由三维、不可压、定常形式的连续性、动量、能量守恒方程所控制。如下所示。

1) 当流场无源项时，单位时间流体微团的质量改变量等于单位时间内流入或流出的流体微团质量，连续性方程如下所示：

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0$ (1)

式中， $\rho$ 为分层流的背景密度， $v$ 为速度矢量。

2) 单位时间内流体微团的动量变化率等于作用在流体微团上的合外力，动量方程如下：

$\frac{\partial }{\partial t}\left(\rho v\right)+\nabla \cdot \left(\rho vv\right)=-\nabla p+\nabla \cdot \left(\overline{\stackrel{¯}{\tau }}\right)+\rho g$ (2)

式中， $\overline{\stackrel{¯}{\tau }}$ 为切应力。

3)能量守恒方程：

$\frac{\partial }{\partial t}\left(\rho T\right)=-\nabla \cdot \left(\rho vT\right)+\nabla \cdot \left(\frac{k}{c}\nabla \left(T\right)\right)$ (3)

式中，T为温度，k为热导率，c为比热容。

分层流体的背景密度ρ是关于温度T、盐度S、压力P的函数，即 $\rho =f\left(T,S,P\right)$ 。不同海域、不同深度分布中，温度、盐度、压强变化对密度的影响不同。而在实际海洋的开阔海域中，盐度和压力的变化较小，温度分层成为垂向密度分层的主导因素 [11] 。温度和密度随水深(Z)变化的函数如(4)、(5)所示，利用UDF编译实现密度线性分层设置。

$T=T_0+k_t\Delta z$ (4)

$\rho ={\rho }_0+k_{\rho }\Delta z$ (5)

式中， $T_0$ 为平均温度， $k_t$ 为温度梯度，本文中取 $k_t=1$ ， $\Delta z$ 为水深变化值, ${\rho }_0$ 为平均密度， $k_{\rho }$ 为密度梯度。

结合式(2)可知定常z方向压强分布如下：

$\frac{\partial p}{\partial z}=\rho g$ (6)

由(4)、(6)可得

$p={\rho }_{0}g\Delta z+\frac{1}{2}{k}_{\rho }g\Delta z^2$ (7)

$\mathrm{Re}=\frac{UD}{v}$ (8)

$Fr=\frac{U}{ND}$ (9)

$N=\sqrt{\frac{g\Delta \rho }{{\rho }_{0}H}}$ (10)

式中，U为球体运动的速度，D为球体直径，υ为流体运动粘度，N为浮频率，H为流域高度。

2.2. 计算域与网格划分

本文以运动球体为对象，利用CFD (Computational Fluid Dynamic)数值模拟球体在密度线性分层流中的尾迹流场特征。构造了一个长L = 80D，宽W = 20D，高H = 20D的长方体计算域，如图1所示。小球距离入口20D，上、下壁面各10D，两侧壁面各10D。流体域的壁面和中心采用网格加密处理(图2)，为了更好地捕捉球体附近的流场特征，在球体附近划分O-block，并在球体周围加入边界层网格进行加密，如图3所示。

Figure 1. Computational domain of the stratified flow pasta sphere

图1. 流体计算域

Figure 2. Mesh in the vertical median surface

图2. 垂直中心面网格

Figure 3. Boundary layer grid near the sphere

图3. 球体附近的边界层网格

2.3. 模型与初始化条件设置

计算域进口采用速度入口，设置固定的流体速度，代表小球相对于周围流体的匀速运动。出口处给定压力出口，分布形式如式(7)所示。实践表明，给定压力分布能带来更好的计算收敛性。左右壁面采用对称设置，上下壁面设置为滑移壁面，壁面移动速度与进口速度保持一致，并在上下壁面分别赋予对应的上下壁面温度。

采用Large Eddy Simulation (LES)模型，应用Couple算法对离散方程进行求解。空间离散化梯度采用最小二乘单元格法，压力采用Body Force Weighted，动量方程采用Bounded Central Differencing,能量方程采用Second Order Upwind。初始化采用Standard Initialization并通过patch温度，初始化温度场。

分层流的密度线性分布通过式(4)、(5)实现。利用UDF设置温度随深度H的线性分布，关联密度与温度方程，进而实现密度随深度H的线性分布。

本文选择Re = 3700，可以通过Yun [12] 等人的模拟对非分层流过球体的LES模型进行验证，见表1。

Table 1. Flow parameters

表1. 流动参数

3. 速度场

3.1. 瞬态流场

如图4所示的速度等值面图，Fr = 0.05时，饼状垂直相干涡结构在球后交替脱落(图4(a))，垂直涡脱落结构与卡门涡街结构相似，但与卡门涡街不同的是，Fr = 0.05涡结构在垂直方向上受到分层效应的抑制，球后涡结构以片状脱落，Chomaz [4] 将垂直方向上的尾迹结构称为三层结构(上边界层、中心层、下边界层)。图4(b)中，在 $-0.5\le Z/D\le 0.5$ 处发现球后流线与流向方向平行，垂直方向上的流速几乎为0，该区域被称为中心层。呈二维流动(垂向速度几乎为0)的中心层形成原因与浮力作用有关，在较低Fr流体中，浮力主导惯性力，当分层流稳定的结构被破坏时，流体位移被迫发生移动。在浮力的作用下，移动的流体所具有的惯性力被周围流体施加的浮力克服时，则无法从垂直方向上翻越球体，垂向速度几乎为0，流体只能围绕球体水平两侧流动。

Fr = ∞时，球体处于未分层流中，球后尾迹呈现出完全不一致的形态。与Fr = 0.05对比，垂直方向上的中心层消失，取而代之的是球后形成的回流区，回流区内包裹着非对称的涡结构，在 $X/D=1.5$ 时涡结构呈螺旋式脱落。为了进一步分析分层对尾迹流场的影响，对流场进行时均化处理。

3.2. 时均流场

Fr = 0.05、∞时均化后的流场如图5所示。Fr = ∞时尾迹呈各向同性，在垂直、水平方向上分布均匀：水平面、垂直面上的尾迹高度、回流区长度、形态几乎一致。而Fr = 0.05时，尾迹宽度(图4(a))明显大于尾迹高度(图4(b))，速度分布更倾向于水平方向。Fr = 0.05，垂直方向上受到强烈的浮力作用，使得垂直方向上的时均流场与瞬态流场呈现的尾迹结构几乎一致。但在水平方向上，图4(a)脱落涡结构被时均，图5(a)脱落涡结构消失，说明Fr = 0.05时尾迹扰动主要集中在水平方向上，垂直方向上的扰动被抑制。与此同时，对比图5中时均场，发现Fr = 0.05时，水平中心面尾迹两侧存在对称翼型结构(见图5(a))，说明此处的流向时均速度明显高于背景流体和尾迹速度，表示Fr = 0.05时流动更偏向于水平掠过球体。

Figure 4. Instantaneous Velocity contours of Fr = 0.05, ∞. (a) Fr = 0.05, Z = 0; (b) Fr = 0.05, Y = 0; (c) Fr = ∞, Z = 0; (d) Fr = ∞, Y = 0

图4. Fr = 0.05、∞时瞬态速度云图。(a) Fr = 0.05，Z = 0；(b) Fr = 0.05，Y = 0；(c) Fr = ∞，Z = 0；(d) Fr = ∞，Y = 0

Figure 5. Time averaged velocity contours of Fr = 0.05, ∞. (a) Fr = 0.05, Z = 0; (b) Fr = 0.05, Y = 0; (c) Fr = ∞, Z = 0; (d) Fr = ∞, Y = 0

图5. Fr = 0.05、∞时时均速度云图。(a) Fr = 0.05，Z = 0；(b) Fr = 0.05，Y = 0；(c) Fr = ∞，Z = 0；(d) Fr = ∞，Y = 0

3.3. 速度分布

对比时均和瞬态流场在水平、垂直中心面上的云图，发现在分层流尾迹中，水平方向上的扰动比垂直方向上的更大。为了更直观表示，对水平方向和垂直方向处的扰动速度进行定量分析，如图6所示。图6(a)为水平方向的扰动速度 $v^{\prime }/U_{\infty }$ 在流向上的分布，Fr = 0.05在水平方向上扰动速度分布远大于Fr = ∞，且呈周期性逐渐衰减。图6(b)为垂直方向扰动速度 $w^{\prime }/U_{\infty }$ 在流向上的分布，但Fr = ∞垂直方向上扰动速度 $w^{\prime }/U_{\infty }$ 沿流向分布远大于Fr = 0.05，其分量与 $v^{\prime }/U_{\infty }$ 流向分布相近。其中Fr = 0.05， $w^{\prime }/U_{\infty }$ 几乎为0。由此可知，分层效应对 $w^{\prime }/U_{\infty }$ 存在抑制，Fr = 0.05时，尾迹扰动主要在水平方向上，引起尾迹水平涡周期性脱落。

Figure 6. Distribution of horizontal and vertical disturbance velocity in the flow direction (X/D) of Fr = 0.05, ∞. (a) Horizontal disturbance velocity $v^{\prime }/U_{\infty }$ distribution; (b) Horizontal disturbance velocity $w^{\prime }/U_{\infty }$ distribution

图6. Fr = 0.05、∞时，水平方向扰动速度和垂向扰动速度在流向上(X/D)的分布。(a) 水平方向扰动速度分布 $v^{\prime }/U_{\infty }$ ；(b) 垂向扰动速度 $w^{\prime }/U_{\infty }$

3.4. 速度衰减

Fr = ∞时，速度衰减曲线与Spedding [16] 提出的3D湍流曲线一致，沿流向方向以 $X^{-2/3}$ 速度衰减。但在Fr = 0.05时，湍流尾迹衰减速度逐渐减缓，说明分层的存在能够减缓尾迹衰减。

Figure 7. Distribution of Centrel defect velocity $U_o/U_{\infty }$ and $urms/U_{\infty }$ 、 $vrms/U_{\infty }$ 、 $wrms/U_{\infty }$ in the flow direction of (a) Fr = 0.05; (b) Fr = ∞

图7. 中心缺陷速度 $U_o/U_{\infty }$ 与沿流向方向的 $urms/U_{\infty }$ 、 $vrms/U_{\infty }$ 、 $wrms/U_{\infty }$ 分布。(a) Fr = 0.05；(b) Fr = ∞

Fr = ∞，如图7(a)所示，在尾迹回流区 $0\le X/D\le 1.55$ 内，球后中心缺陷速度大于来流速度( $U_o/U_{\infty }\ge 1$ )，此时各向均方根速度也 $0\le X/D\le 1$ 内随 $X/D$ 而增加，流向均方根 $urms/U_{\infty }$ 在 $X/D=1$ 之后开始衰减。 $1.55\le X/D\le 60$ ，中心缺陷速度与各向均方根速度随 $X/D$ 增加而降低，其中，中心缺陷速度衰减曲线速率为 $X^{-2/3}$ ，各向均方根衰减速度快于中心缺陷速度，且衰减呈各向同性，uvw三个方向的衰减速率相近，为 $X^{-1.38}$ 。与均匀流尾迹的速度衰减对比，分层效应的增加引起尾迹各向均方根速度衰减各向异性。在Fr = 0.05时， $vrms/U_{\infty }$ 在 $X/D=8$ 处占主导，流向上分布大于中心缺陷速度、 $urms/U_{\infty }$ 、 $wrms/U_{\infty }$ 。 $wrms/U_{\infty }$ 在 $X/D=0.5$ 处小 $vrms/U_{\infty }$ 一个量级，此时的尾迹以二维结构存在，抑制垂向速度；在 $X/D=8$ 之后 $vrms/U_{\infty }$ 迅速降低，逐渐趋于 $urms/U_{\infty }$ 。

4. 漩涡结构

Fr = ∞时，在X/D = 1.5之前漩涡环在球后保持圆环状，以螺旋式脱落。对比低雷诺数下的非分层流尾迹，Kim & Durbin [13] 观察到均匀流脱落漩涡的尾迹不稳定性与Re有关，Re > 800时处于双频状态，即同时受到K-H不稳定性和螺旋不稳定性的影响。Chomaz [4] 对Kim & Durbin [13] 的结论进一步补充，Chomaz [4] 利用流动可视化和局部探头测量对150 ≤ Re ≤ 30,000下均匀流尾迹空间结构和时间行为进行分析，发现在该雷诺数范围内激发了3种不同的不稳定模式：循环区的轴对称脉动，没有涡脱落；与循环区外围剪切相关的轴对称涡旋脱落；螺旋模式，方位波数为1，与分离线的旋转有关。当Re ≥ 800时(本文Re = 3700)，球后同时受到KH不稳定性和螺旋不稳定性影响。

当在背景流体中加入密度分层时，如图8(a)、图8(b)所示。此时的涡结构呈现饼状涡街，水平方向上的圆饼状涡结构沿流向交替脱落，而垂直方向上则受到分层效应的抑制，依旧保持三层结构，其中上下分离边界层脱落速度较快，牵引着中心层脱落的涡结构，在垂直方向呈“V”状。

Figure 8. Vorticity contours of Fr = 0.05, ∞. (a) Fr = 0.05, Z = 0; (b) Fr = 0.05, Y = 0; (c) Fr = ∞, Z = 0; (d) Fr = ∞, Y = 0

图8. Fr = 0.05、∞涡量云图。(a) Fr = 0.05，Z = 0；(b) Fr = 0.05，Y = 0；(c) Fr = ∞，Z = 0；(d) Fr = ∞，Y = 0

5. 结论

本文采用LES对圆球绕流尾迹从速度场、涡结构、内波展开研究，分析在亚临界雷诺数(Re = 3700)下，Fr = 0.05、∞对尾迹的影响,得出的结论如下：分层效应的存在会引起尾迹速度分布各向异性,主要体现对垂直方向上的抑制。Fr = ∞时，球后尾迹以三维湍流结构存在，在水平和垂直方向上的分布均匀。背景流体中增加密度分层后，垂向方向上会受到抑制，流体逐渐偏向球体水平方向流动，导致尾迹结构由分布均匀的3D结构逐渐向片状大尺度相干涡结构过渡。对比Fr = ∞和Fr = 0.05的速度分布和扰动速度分量发现，Fr = 0.05垂向速度分量明显受到抑制，流体以球体水平两侧流动为主。