压力梯度作用下相干结构的演化

摘要

采用共振三波的一个周期作为湍流边界层近壁区相干结构的初值，用直接数值模拟方法对有压力梯度，包括定常压力梯度和变压力梯度作用下相干结构的演化进行了研究，得出其各种特性的变化与实验观测到的结果一致。压力梯度对相干结构初始扰动波的选择以及演化过程都起作用，相比之下对前者的作用更大。条纹结构和流向涡之间有某种内在的紧密的联系，可能是产生流向涡的起因。

正文

压力梯度是工程流动中的普遍现象, 例如沿飞机、船舶、叶轮机叶片等表面, 压力梯度沿流向是不断变化的, 这使得边界层中的湍流经常处在非平衡状态。因此, 弄清楚不同压力梯度对湍流相干结构的影响, 对解决湍流计算的困难有特别重要的意义^[1]。参照文献[2], 笔者采用直接数值模拟方法, 将不稳定波串中的一个周期作为相干结构的初值, 对其在压力梯度作用下的演化进行了计算分析。

《1 理论分析及数值方法》

1 理论分析及数值方法

在近壁区由于粘性作用很强, 湍流的小尺度脉动很弱, 可以将瞬时速度场和压力场分解为基本流项和相干结构项:

$u = U + u_{c}, p = Ρ + p_{c}, (1)$ $u = U + u_{c}, p = Ρ + p_{c}, (1)$

其中U, P表示基本流项, u_c, p_c 表示相干结构项。将它们代入N-S方程可得出相干结构项满足的扰动方程:

$\begin{array}{l} \frac{\partial u_{c}}{\partial t} + (U \cdot Δ) u_{c} + (u_{c} \cdot Δ) U = \\ \frac{1}{R} Δ^{2} u_{c} - Δ p_{c} - (u_{c} \cdot Δ) u_{c}, (2 a) \\ Δ \cdot u_{c} = 0 。 (2 b) \end{array}$ $\begin{array}{l} \frac{\partial u_{c}}{\partial t} + (U \cdot Δ) u_{c} + (u_{c} \cdot Δ) U = \\ \frac{1}{R} Δ^{2} u_{c} - Δ p_{c} - (u_{c} \cdot Δ) u_{c}, (2 a) \\ Δ \cdot u_{c} = 0 。 (2 b) \end{array}$

基本流U ( $\bar{u}$ $\bar{u}$ (x, y) , $\bar{v}$ $\bar{v}$ (x, y) , 0) 由解N-S方程直至稳定而得, 其初始场为经法向压缩的层流剖面一部分, 边界条件为在y⁺=100处与湍流平均速度光滑连接, 关于采用这种剖面的合理性参见文献[2]。

式 (2a) 和式 (2b) 在展向 (Z方向) 采用Fourier谱展开。可得到N个二维方程组:

$\begin{array}{l} \frac{\partial u_{m}}{\partial t} + F_{m} [(u_{c} \cdot Δ) u_{c} + (U \cdot Δ) u_{c} + \\ (u_{c} \cdot Δ) U] + Δ_{m} p_{m} = \frac{1}{R} Δ_{m}^{2} u_{m}, (3 a) \end{array}$ $\begin{array}{l} \frac{\partial u_{m}}{\partial t} + F_{m} [(u_{c} \cdot Δ) u_{c} + (U \cdot Δ) u_{c} + \\ (u_{c} \cdot Δ) U] + Δ_{m} p_{m} = \frac{1}{R} Δ_{m}^{2} u_{m}, (3 a) \end{array}$

$Δ_{m} \cdot u_{m} = 0$ $Δ_{m} \cdot u_{m} = 0$ , (3b)

$\begin{array}{l} m = - \frac{Ν}{2}, \dots, \frac{Ν}{2} - 1 \\ Δ_{m} = {\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, - i m β}, Δ_{m}^{2} = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} - m^{2} β^{2} 。 \end{array}$ $\begin{array}{l} m = - \frac{Ν}{2}, \dots, \frac{Ν}{2} - 1 \\ Δ_{m} = {\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, - i m β}, Δ_{m}^{2} = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} - m^{2} β^{2} 。 \end{array}$

其中:N是截断误差的阶数, β为展向波数, F_m为该项的Fourier谱展开系数。

对于时间离散, 在这里采用混合显-隐的分裂格式, 式 (3a) 的具体离散过程为

$\begin{array}{l} \frac{u^{'}_{m} - \sum_{q = 0}^{J_{i} - 1} α_{q} u_{m}^{n - q}}{Δ t} = - \sum_{q = 0}^{J_{e} - 1} β_{q} F_{m} [(u_{c}^{n - q} \cdot Δ) u_{c}^{n - q} + \\ (U^{n - q} \cdot Δ) u_{c}^{n - q} + (u_{c}^{n - q} \cdot Δ) U^{n - q}], (4 a) \\ \frac{u^{″}_{m} - u^{'}_{m}}{Δ t} = - Δ_{m} p_{m}^{n + 1} ‚ (4 b) \\ \frac{γ_{0} u_{m}^{n + 1} - u^{″}_{m}}{Δ t} = \frac{1}{R} Δ_{m}^{2} u_{m}^{n + 1}, (4 c) \end{array}$ $\begin{array}{l} \frac{u^{'}_{m} - \sum_{q = 0}^{J_{i} - 1} α_{q} u_{m}^{n - q}}{Δ t} = - \sum_{q = 0}^{J_{e} - 1} β_{q} F_{m} [(u_{c}^{n - q} \cdot Δ) u_{c}^{n - q} + \\ (U^{n - q} \cdot Δ) u_{c}^{n - q} + (u_{c}^{n - q} \cdot Δ) U^{n - q}], (4 a) \\ \frac{u^{″}_{m} - u^{'}_{m}}{Δ t} = - Δ_{m} p_{m}^{n + 1} ‚ (4 b) \\ \frac{γ_{0} u_{m}^{n + 1} - u^{″}_{m}}{Δ t} = \frac{1}{R} Δ_{m}^{2} u_{m}^{n + 1}, (4 c) \end{array}$

u′_m, u″_m为中间速度, J_i, J_e为格式精度参数, α_q, β_q, γ₀是适当权数。

对流项的空间离散采用五阶迎风紧致格式计算^[3]。压力梯度项中通过令 $Δ_{m} \cdot u^{″}_{m} = 0$ $Δ_{m} \cdot u^{″}_{m} = 0$ 得到

《图1》

$b = m^{2} β^{2} ‚ p_{m} = Φ, f = \frac{Δ_{m} \cdot u_{m}}{Δ t}$ $b = m^{2} β^{2} ‚ p_{m} = Φ, f = \frac{Δ_{m} \cdot u_{m}}{Δ t}$ 得到压力函数Φ的Helmholtz方程:

$\frac{\partial^{2} Φ}{\partial x^{2}} + \frac{\partial^{2} Φ}{\partial y^{2}} - b Φ = f 。 (5)$ $\frac{\partial^{2} Φ}{\partial x^{2}} + \frac{\partial^{2} Φ}{\partial y^{2}} - b Φ = f 。 (5)$

对于散度, 在内点采用六阶对称紧致有限差分格式 (例如 $\frac{\partial u}{\partial x}$ $\frac{\partial u}{\partial x}$ , 令 $F = \frac{\partial u}{\partial x})$ $F = \frac{\partial u}{\partial x})$

$\begin{array}{l} F_{i + 1} + 3 F_{i} + F_{i - 1} = \\ \frac{1}{12 Δ x} (u_{i + 2} - u_{i - 2}) + \frac{7}{3 Δ x} (u_{i + 1} - u_{i - 1}) 。 (6) \end{array}$ $\begin{array}{l} F_{i + 1} + 3 F_{i} + F_{i - 1} = \\ \frac{1}{12 Δ x} (u_{i + 2} - u_{i - 2}) + \frac{7}{3 Δ x} (u_{i + 1} - u_{i - 1}) 。 (6) \end{array}$

在邻近边界点上, 用四阶对称紧致有限差分格式 ;在边界上, 采用三阶精度的紧致有限差分格式。对于式 (5) , 为同上述方程的差分格式精度相匹配, 采用了9点四阶精度紧致有限差分格式。粘性扩散方程采用四阶精度对称紧致有限差分格式^[4]。

入口速度及压力条件:

$x = 0, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, p_{c} = 0;$ $x = 0, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, p_{c} = 0;$

出口速度及压力条件:

$x = x_{m}, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, \frac{\partial p_{c}}{\partial x} = 0;$ $x = x_{m}, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, \frac{\partial p_{c}}{\partial x} = 0;$

壁面速度压力条件:

$y = 0, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, \frac{\partial p_{c}}{\partial y} = 0;$ $y = 0, \tilde{u}_{c} = \tilde{v}_{c} = \tilde{w}_{c} = 0, \frac{\partial p_{c}}{\partial y} = 0;$

外缘速度压力条件:

$y = y_{n}, \frac{\partial u_{c}}{\partial y} + α u_{c} \approx 0, \frac{\partial p_{c}}{\partial y} = 0;$ $y = y_{n}, \frac{\partial u_{c}}{\partial y} + α u_{c} \approx 0, \frac{\partial p_{c}}{\partial y} = 0;$

展向 (Z) 采用周期性边界条件。

计算域为:在展向取约1个周期的长度, 流向取约为5个周期的长度, 法向取100个粘性长度。计算域 (x, y, z) 上的网格为50×100×32, 时间步长为0.02。

相干结构的初值由一般共振三波模型给出, 即由三个三维波叠加并取一个周期内的值得出, 具体形式参见文献[2]。

《2 计算结果及讨论》

2 计算结果及讨论

《2.1定常压力梯度流动》

2.1定常压力梯度流动

首先分别对逆压、零压和顺压三种定常压力梯度流动进行了计算, 其中压力梯度参数分别为:β_l=-0.1, 0, 0.1, (β_l为压缩层流剖面的压力梯度参数) 。参照文献[2]的结果, 通过大量的计算, 分别找出了三种压力梯度下幅值相对增长最快的共振三波。

图1为三种压力梯度下幅值A的增长曲线, 从图中可见, 逆压梯度下相干结构的幅值增长得最快, 顺压情况下幅值增长的幅度较小, 甚至不增长。另外计算中还发现, 在有逆压梯度作用的情况下, 在很大的范围内都存在能够迅速增长起来的共振三波, 而在顺压梯度下绝大部分的共振三波都是衰减的, 能够增长的波的范围非常小, 说明逆压梯度使得非线性作用加强, 对相干结构具有激励作用;而顺压梯度对相干结构具有抑制作用。

《图2》

图1 三种情况下相干结构幅值随时间的演化

Fig.1 The time evolution curves of the amplitudes of coherent structures

——逆压 ---零压 —··—··—顺压

《图3》

图2 $\tilde{u}$ $\tilde{u}$ , $\tilde{v}$ $\tilde{v}$ 等值线 (实线为正, 虚线为负, 下同)

Fig.2 Contours of $\tilde{u}$ $\tilde{u}$ , $\tilde{v}$ $\tilde{v}$ in X-Y plane

图2分别给出了顺压和逆压两种情况下展向最大雷诺应力位置X-Y平面上扰动速度 $\tilde{u}$ $\tilde{u}$ , $\tilde{v}$ $\tilde{v}$ 的等值线分布。可以看出, 逆压梯度情况下高速和低速相间区域比顺压梯度的多, 单个区域比较小, 速度梯度相对要变大, 这也可以部分反映逆压梯度下相干结构更加活跃。从图2中还可以看出, 逆压梯度作用下相干结构的整体倾向于垂直发展。

图3分别给出了三种压力梯度作用下, 演化后期扰动速度 $\tilde{u}$ $\tilde{u}$ 在X-Z平面上的等值线分布, 可以看出在逆压梯度情况下, 结构变得复杂, 出现密集的高、低速条纹交错排列的现象, 形成了大量的高剪切层, 扰动速度剧烈扭曲。不仅在展向, 而且在流向也出现了大量的高、低速流体交错排列现象, 反映了结构的展向运动很活跃, 而且可以得出在逆压梯度下近壁区不仅流向涡很强, 一定还存在大量的展向涡, 这和文献[5]的实验结果相吻合。而在顺压梯度情况下, 条纹结构比较稳定, 沿展向的脉动非常小, 条纹在展向形成均匀排列, 相互作用变得很弱。从图3中还可以看出, 逆压梯度下结构的流向尺度最短, 而顺压梯度下的流向尺度最长, 而且其流向的传播速度最快。这一结果和文献[6]的实验结果是一致的。

《图4》

图3 $\tilde{u}$ $\tilde{u}$ 在X-Z平面上的等值线分布t=50, y⁺=40

Fig.3 Contours of $\tilde{u}$ $\tilde{u}$ in X-Z plane

(上) 逆压梯度, (中) 零压梯度, (下) 顺压梯度

图4给出了顺压和逆压梯度两种情况下, 流向涡最强的流向位置处的速度矢量图。各图中矢量箭头的长度是由相同的基准量归一化得出的, 由此就可以看出它们之间量的区别来。

《图5》

图4 Y-Z平面的速度矢量图 t=50

Fig.4 Velocity vectors in Y-Z plane

通过比较可以明显看出, 顺压梯度下流向涡比较弱, 并且在展向呈大小相等、方向相反的对称分布, 流向涡很稳定, 在展向的摆动非常小。另外, 涡心在法向的抬升幅度也很小。相比之下, 逆压梯度下流向涡的强度很大, 出现多尺度的、多旋向的多涡结构, 流向涡展向的摆动非常剧烈。在逆压情况中, 一般都存在一个或两个尺度最大的流向涡, 其尺度要远大于顺压梯度下流向涡的尺度, 流向涡的形状变得越来越不规则, 如图4所示, 流向涡在法向拉伸呈椭圆形, 这与前面X-Y平面上扰动速度等值线分布的分析结果是一致的。此外, 相邻的多个涡之间的距离变小, 涡与涡之间的作用加强, 并且涡心都沿法向有明显抬升。

图5给出了三种情况下Y-Z平面的雷诺应力 $(- \bar{\tilde{u} \tilde{v}})$ $(- \bar{\tilde{u} \tilde{v}})$ 分布。从图中可以发现, 顺压梯度情况下的雷诺应力值最小, 约为0.002, 而且始终限制在y⁺=20附近。而在逆压梯度情况下, 随着相干结构的演化, 雷诺应力峰值的位置在法向是逐渐增加的, 到t=50的时候, 其位置已经上升到y⁺=40和50之间, 雷诺应力值也达到0.02。这也说明了逆压梯度激励相干结构并促使猝发产生, 而顺压梯度抑制内部的扰动并将其限制在底层。从图5c还可以看出, 在零压梯度情况下, 到了演化后期, 基本上以双峰值或三峰值居多, 而在逆压梯度下到了演化的后期基本上以一个较强的单峰值为主的居多, 这也和上面所研究的流向涡在X-Z平面上的涡矢分布结果是一致的。

《图6》

图5 三种情况下Y-Z平面的雷诺应力分布 t=50

Fig.5 Reynolds-stress distribution in Y-Z

《2.2变压力梯度流动》

2.2变压力梯度流动

本文对逆压梯度线性递减 (β_l由-0.1→0) (情况一) 和线性递增 (β_l由0→-0.1) (情况二) 两种变压力梯度形式下相干结构的空间演化进行了计算。计算结果表明, 前者幅值增长的速度和范围都略大于后者。图6给出了扰动速度 $\tilde{u}$ $\tilde{u}$ 在X-Z平面上的等值线分布。通过比较可以看出, 图6a中的相干结构要比图6b复杂, 即对于情况一中的相干结构要明显比情况二活跃, 如图所示高速 (实线) 、低速 (虚线) 条纹交错排列的密集, 使得两者的作用面增大, 形成大量的速度剪切层, 使流动更加不稳定。仔细观察不难得出, 图6a中低速流体明显多于图6b中的低速流体, 这也反映了逆压梯度在由大变小的流动下流向涡更强, 将大量的底层低速流体“卷”到上层, 增加了内外层流体的动量交换。

《图7》

图6 X-Z平面上 $\tilde{u}$ $\tilde{u}$ 的等值线分布t=50, y⁺=40

Fig.6 Contours of $\tilde{u}$ $\tilde{u}$ in X-Z plane

图7分别给出了两种情况下流向涡最强的流向位置的速度矢量图。与图4相同, 图中矢量箭头的长度是由相同的基准量归一化得出的。演化初期两种情况下的流向涡都很弱, 尺度也很小, 随着演化过程的发展, 非线性作用的不断增强, 流向涡的强度和尺度都得到了增强。通过比较可以看出情况一的流向涡的强度高于情况二, 虽然在演化的中期出现多涡结构, 但都是相同尺度, 强度都以最大的涡为主。到了演化的后期, 基本上都是单涡结构, 整个流场中仅存在一个强度、尺度都很大的涡。这一点和前面恒定逆压梯度的情况是一样的。

《图8》

图7 Y-Z平面的速度矢量图

Fig.7 Velocity vectors in Y-Z plane

图8给出了两种情况Y-Z平面上雷诺应力 $(- \bar{\tilde{u} \tilde{v}})$ $(- \bar{\tilde{u} \tilde{v}})$ 的分布。通过对比可以看出前者的值要大于后者。另外, 两种情况在演化后期的雷诺应力基本上也都是以单峰值为主。在这里更加明确地体现了雷诺应力强弱和位置与流向涡的强度、位置的对应关系。从上述结果可以得出第一种情况对相干结构的激励作用更大, 与实验结果是一致的^[7]。

《图9》

图8 Y-Z平面上雷诺应力的分布

Fig.8 Reynolds-stress distributions in X-Z plane

在项目的研究中, 压力梯度在两个方面影响相干结构演化的最终结果, 一是利用流动稳定性理论选择共振三波作为相干结构的初始值;二是在共振三波的数值演化过程中, 为了进一步研究压力梯度在这两个方面的作用, 以零压梯度下增长最快的共振三波的初值作为相干结构初始值, 计算其在逆压梯度的平均流场中的演化, 并与定常零压和逆压结果做对比分析, 幅值演化结果见图9。图中同时给出了另两种情况的结果 (与图1相同) 。对比图9中1, 3可以看出, 对于相同的初始扰动 (零压梯度下增长最快的共振三波的初值) , 其幅值在逆压梯度下的流场比在零压梯度下的流场中增长得更快, 这说明逆压梯度对流场中各种扰动的演化具有激励作用。对比图9中1, 2可以看出, 在相同的逆压流场中, 不同的初始扰动增长快慢也不同, 2 最终比1增长得更快。对比三条曲线可以看出, 1与3无论形状还是长幅都更为接近, 说明相干结构初始值对演化结果的影响要大于平均流场的作用。这就意味着压力梯度在初始扰动波的选择上比演化过程中起的作用更大。由此也可以解释前面两种变压力梯度流动的计算结果。

《图10》

图9 相干结构幅值随时间的变化

Fig.9 The time evolution curves of the amplitudes of coherent structures

1 ——○—○—○—零压初值逆压流场 2 ——◇—◇—◇—定常逆压梯度3 ——△—△—△—定常零压梯度

《2.3猝发》

2.3猝发

在研究中发现, 对所有情况都存在一个普遍的现象, 即在绝大多数情况下法向扰动速度 ( $\tilde{v}$ $\tilde{v}$ ) 的负峰值比正峰值大, 并且雷诺应力的峰值更接近于法向扰动速度的负峰值的位置, 这种现象在变压力梯度的流动中更为明显。以β_l=-0.1→0流动为例, 图10给出了雷诺应力最大位置的法向扰动速度的分布, 对照图10、图8a和图7a可以发现, 流向涡左右两侧的强度明显不对称, 卷向壁面的一侧强度要明显强于卷离壁面的一侧, 卷向壁面的一侧在位置上对应着扰动速度的 $\tilde{v}$ $\tilde{v}$ 负峰值, 卷离壁面的一侧对应着扰动速度 $\tilde{v}$ $\tilde{v}$ 的正峰值, 这意味着在湍流边界层近壁区扫掠的强度大于喷射的强度, 而雷诺应力大部分是在扫掠过程中产生的。这与文献[8]中的统计平均结果是一致的。

由前面的研究可以看出, 雷诺应力和条纹结构、流向涡之间不仅在强度、结构上存在对应关系, 而且还在位置和分布上相对应, 这进一步说明了这三者存在着紧密的内在联系。研究中, 初始的扰动值是由三个满足共振关系的三维波给出的, 相干结构的初场同时含有条纹结构和一些准流向涡的元素, 通过仔细分析相干结构的初场, 发现这些准流向涡的元素很弱, 而初始的条纹结构相比要强得多, 到了演化的中后期流向涡结构大大地增强了, 这间接地表明流向涡的增强可能是条纹结构作用的结果, 这一结论与Hussain^[9]的研究结果是一致的。

《图11》

图10 X-Z平面上法向扰动速度 $\tilde{v}$ $\tilde{v}$ 的分布β_l=-0.1→0

Fig.10 Distribution of $\tilde{v}$ $\tilde{v}$ in X-Z plane, β_l=-0.1→0

《3 结论》

3 结论

笔者采用Fourier谱方法和高阶紧致差分格式相结合的计算方法, 并结合流动稳定性的理论对有压力梯度下湍流近壁区的相干结构的演化过程进行直接数值模拟, 得到以下结论:

1) 逆压梯度对相干结构有激励作用, 使得相干结构的流向尺度变小, 流向传播速度变慢;三种压力梯度下慢斑的展向间距都在100个粘性长度左右, 顺压下相对最大。

2) 在逆压梯度作用下, 流向涡很强, 法向的抬升和拉伸都很剧烈, 有明显的不对称性, 而在顺压梯度下, 流向涡则是两个大小近似相等, 在展向呈对称分布的形状规则的弱涡。

3) 和顺压梯度情况相比, 逆压梯度作用下雷诺应力增长很快, 并且以单峰为主, 呈明显的不对称分布, 并且峰值在法向离壁面较远, 与结论2 相对应。

4) 压力梯度对相干结构初始扰动波的选择以及演化过程都起作用, 相比之下对前者的作用更大。

5) 条纹结构和流向涡之间有某种内在的紧密联系, 可能是产生流向涡的起因。

致谢本文的研究工作得到天津大学周恒教授、罗纪生教授的热心指导和帮助, 特此致谢。

展示更多