引言

均匀各向同性湍流是一种理想的湍流简化形态^[1-2].为了在数值上生成均匀各向同性湍流，目前研究者广泛采用的做法是在三维空间内引入周期性边界条件^[3-10]，将无限大空间近似为无穷多个周期性的立方体.这种做法其实是由于数值分辨率的限制，必须在谱空间进行尺度截断而引入的.小尺度的截断与数值误差、能量积累等问题有关^[11-12]，而周期性条件则是一种最常见的大尺度的截断.这种大尺度截断会引起一系列的问题，包括间歇性^[13]、各向异性等问题.显然，如果由于引入三维周期性条件带来额外的各向异性，势必会影响我们希望生成的“各向同性”湍流场的质量.

文献[14]对此问题已经进行了初步的讨论，将湍流的各向异性问题归结为张量各向异性和标量方向各向异性两种.前者指的是二阶纵向结构函数D_ll (下标l表示结构函数纵向分量)和横向结构函数D_nn (下标n表示结构函数横向分量)在不同矢径长度r下需满足各向同性关系式^[15-18]

$ {D_{{\rm{nn}}}}(r) = {D_{{\rm{ll}}}}(r) + \frac{r}{2}\frac{{{\rm{d}}{D_{{\rm{ll}}}}(r)}}{{{\rm{d}}r}} $

(1)

这在一些数值模型验证工作中曾被作为各向同性检测方法使用^[19-20].后者则指对于标量两点统计量，其统计值应与方向无关，只与两点距离有关.前期工作已经证明^[14]，采用三维周期性条件必定会在大尺度引起棱与对角线之间的标量方向各向异性，其引起的各向异性差异可达10%.这种各向异性可以通过引入六边形谱方法来尽量减小^[21].

本文将关注采用三维周期性条件时初始场3个棱方向的各向异性问题.此各向异性问题不仅在二阶纵向结构函数中体现，而且甚至可在各个速度分量的概率密度分布中体现，这是一种由于样本波动较大而引起的统计问题.本文提出一种改进方法，以减小样本波动，生成各向同性程度更好的初始湍流场.

1 Rogallo初始场生成方法及其棱向各向异性 1.1 Rogallo初始场生成方法

在进行计算机数值模拟时，均匀各向同性湍流场的流场区域一般被设计为一个棱长为2π的立方体形状.这种流场普遍使用Rogallo^[22]在1981年提出的方法进行初始化.该方法的核心思想是在谱空间中构造一个同时满足连续性方程以及给定能谱的流场

$ {k_1}\hat U(\mathit{\boldsymbol{k}}) + {k_2}\hat V(\mathit{\boldsymbol{k}}) + {k_3}\hat W(\mathit{\boldsymbol{k}}) = 0 $

(2)

$ \mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_{A(k)} {\mathit{\boldsymbol{\widehat u}}\left( \mathit{\boldsymbol{k}} \right) \cdot {{\mathit{\boldsymbol{\widehat u}}}^*}\left( \mathit{\boldsymbol{k}} \right)} {\rm{d}}A(k) = E(k) $

(3)

其中，k_i和$ \hat {U} $, $ \hat {V} $, $ \hat {W} $分别为谱空间波数矢量k和速度矢量$ \boldsymbol{\widehat u} $在3个棱方向上的分量，$ \boldsymbol{\widehat u} $^*为速度矢量$ \boldsymbol{\widehat u} $的共轭矢量，A(k)为谱空间中半径为k(波数矢量k的模长)的球面积分域，E(k)为给定能谱.连续性方程要求谱空间中速度矢量与波数矢量垂直^[23].因此，Rogallo方法根据每一个波数矢量k建立局部坐标基(e'₁, e'₂, e'₃)，并令e'₃矢量与k矢量平行，即

$ {{\mathit{\boldsymbol{e'}}}_3}(\mathit{\boldsymbol{k}}) = \frac{{{k_1}}}{k}{\mathit{\boldsymbol{e}}_1} + \frac{{{k_2}}}{k}{\mathit{\boldsymbol{e}}_2} + \frac{{{k_3}}}{k}{\mathit{\boldsymbol{e}}_3} $

(4)

其中，e₁, e₂, e₃分别为坐标基矢量.则速度矢量$ \boldsymbol{\widehat u} $必须处于由e'₁, e'₂张成的平面上，写为

$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) = \alpha (k){{\mathit{\boldsymbol{e'}}}_1}(\mathit{\boldsymbol{k}}) + \beta (k){{\mathit{\boldsymbol{e'}}}_2}(\mathit{\boldsymbol{k}})\; $

(5)

式中，α(k)与β(k)为速度矢量在方向e'₁与e'₂方向上的分量.由于速度场还需要满足能谱方程，因此，两个分量的最终形式为

$ \alpha (k) = \sqrt {\frac{{E\left( k \right)}}{{4{\rm{\pi }}{k^2}}}} \exp \left( {{\rm{i}}{\theta _1}} \right)\cos \phi $

(6)

$ \beta (k) = \sqrt {\frac{{E\left( k \right)}}{{4{\rm{\pi }}{k^2}}}} \exp \left( {{\rm{i}}{\theta _2}} \right)\sin \phi $

(7)

式中，θ₁, θ₂和φ都是在[0, 2π)区间上均匀分布的随机标量.令

$ {{\mathit{\boldsymbol{e'}}}_1}(\mathit{\boldsymbol{k}}) = \frac{{{{\mathit{\boldsymbol{e'}}}_3}(\mathit{\boldsymbol{k}}) \times {\mathit{\boldsymbol{e}}_3}(\mathit{\boldsymbol{k}})}}{{\left| {{{\mathit{\boldsymbol{e'}}}_3}(\mathit{\boldsymbol{k}}) \times {\mathit{\boldsymbol{e}}_3}(\mathit{\boldsymbol{k}})} \right|}} $

(8)

$ {{\mathit{\boldsymbol{e'}}}_2}(\mathit{\boldsymbol{k}}) = \frac{{{{\mathit{\boldsymbol{e'}}}_3}(\mathit{\boldsymbol{k}}) \times {{\mathit{\boldsymbol{e'}}}_1}(\mathit{\boldsymbol{k}})}}{{\left| {{{\mathit{\boldsymbol{e'}}}_3}(\mathit{\boldsymbol{k}}) \times {{\mathit{\boldsymbol{e'}}}_1}(\mathit{\boldsymbol{k}})} \right|}} $

(9)

得到谱空间中的Rogallo初始场

$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) = \frac{1}{{k\sqrt {k_1^2 + k_2^2} }}\left( {\begin{array}{*{20}{c}} {\alpha (k)k{k_2} + \beta (k){k_1}{k_3}}\\ {\beta (k){k_2}{k_3}-\alpha (k)k{k_1}}\\ {-\beta (k)(k_1^2 + k_2^2)} \end{array}} \right) $

(10)

1.2 Rogallo初始场的棱向各向异性

采用一个典型的算例来说明采用Rogallo方法生成的初始场的棱向各向异性问题. 图 1所示为某个空间网格点数为256³的Rogallo初始场的二阶纵向结构函数D_ll在3个棱方向上的曲线(两点距离Y用网格尺度△进行无量纲化；结构函数用脉动速度均方根的平方进行无量纲化).可以观察到，3条曲线在大尺度处差异较大，显示该初始场在棱之间具有明显的各向异性.需要说明的是，虽然在这里只给出了某个初始场的统计结果，但这种现象不是个例，在实际计算中非常普遍.

为了对比，我们也同时计算了采用Rogallo方法生成的初始场各棱向速度分量U, V, W的概率密度函数(possibility density function，PDF)，如图 2所示(横坐标速度用脉动速度的均方根进行无量纲化).同样可以看到，3个棱向速度分量的PDF之间存在明显的差异.由于PDF和一维能谱之间存在对应关系, 这种差异显示了3个棱向湍动能之间分配的各向异性，也正是我们希望在具体应用中尽量减小的.下面我们从统计意义上给这种差异以定量的描述.

2 针对Rogallo方法统计性质的改进方案 2.1 Rogallo初始场的统计性质

从结构函数的各向异性结果中可以看出，棱向各向异性明显的区域位于大尺度区.这可能是由于在单个初始场中大尺度的流动结构数量太少而引起的样本缺乏所导致的.为了证实这一结论，需要增加样本数以研究初始场之间的统计性质.通过改变式(6)和式(7)中随机序列θ₁, θ₂和φ的随机种子，可以得到满足相同能谱的不同初始场.通过计算机程序生成1 024个基于Rogallo方法的、空间网格点数均为256³的初始场样本，每个初始场的棱向D_ll(r)以及速度3个棱向分量的PDF均存在类似图 1和图 2所示的差异.但如果对这1 024个D_ll标量场及PDF曲线做求和平均处理，则可以观察到这些差异几乎被完全消除，结果如图 3和图 4所示.这说明Rogallo初始化方法本身并不会在平均意义上引起棱向各向异性.尽管如此，实际计算中使用的并不是上千个初始场的平均，而只是单个的初始场，仍然不希望每个初始场的波动太大.因此，有必要对初始场之间的波动进行统计计算.需要强调的是，在这里一个“样本”指的是一个完整的初始场，而并非一个离散点，这是由于实际计算中我们总是采用单个初始场，因而对这种单个样本的波动特性具有较高的要求.总之，在这里将单个初始场视为样本是较合理的做法.

为了定量地表征这些初始场样本的波动大小，定义相对标准差

$ {\sigma _r}(a(t)) = \sqrt {\frac{{\frac{1}{N}\sum\limits_{n = 1}^N {{{({a_n}(t)-\mu )}^2}} }}{{{\mu ^2}}}} $

(11)

式中，N=1 024为初始场的总样本数，a_n(t)表示所统计的对象a(t)的第n个样本，t为相应的自变量，μ为a(t)的平均值，即

$ \mu = \frac{1}{N}\sum\limits_{n = 1}^N {{a_n}(t)} $

(12)

如果我们讨论的是用球坐标表示的二阶纵向结构函数D_ll(r, θ, φ)在3个棱向上的曲线各自的样本波动情况，则a(t)分别为D_{ll, θ=π/2, φ=0}(r), D_{ll, θ=π/2, φ=π/2}(r), D_{ll, θ=0, φ=0}(r)；如果讨论的是3个棱向速度分量PDF各自的样本波动情况，则a(t)分别为U, V, W的PDF.这6个对象的相对标准差统计结果分别如图 5和图 6所示.从图 5可以看到，D_ll(r)的样本波动程度在大尺度更明显，达到约5%. 图 6表明当速度分量的值越靠近平均值0时，其相对波动程度越弱；而在脉动速度较大的区域，由于速度在此处概率密度很低，样本之间的差异非常大，其相对标准差的统计值也因此很大，可达2 000%左右.这种异常大的波动可能会导致额外的间歇性，与Brun等^[13]的工作可能有一定联系，此处不做深入讨论.本文中主要关注湍流能量的主要存在区域，即均值附近的区域，相应的可画出此区域U, V, W的PDF相对标准差，如图 7所示.可以观察到此区域的相对标准差大约为3%.

2.2 模量平均法

为了尽量减小初始场样本之间的波动，得到各向同性程度更好的初始场，需要对Rogallo方法进行分析.事实上，式(10)表明谱空间速度矢量$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) $的第1、第2个分量与α(k)和β(k)均显式相关，但第3个分量只与β(k)显式相关，这导致$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) $的3个分量在谱空间是形式上统计不对称的.经过傅里叶逆变换之后，物理空间的棱向速度分量U, V, W之间也因此存在不对称性.

为了使$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) $三个分量之间在谱空间统计对称，可以对每个k矢量分别令k与e'₃, e'₁, e'₂平行，同时生成3个随机速度场

$ {{\mathit{\boldsymbol{\hat u}}}_1}(\mathit{\boldsymbol{k}}) = \frac{1}{{k\sqrt {k_1^2 + k_2^2} }}\left( {\begin{array}{*{20}{c}} {{\alpha _1}(k)k{k_2} + {\beta _1}(k){k_1}{k_3}}\\ {{\beta _1}(k){k_2}{k_3}-{\alpha _1}(k)k{k_1}}\\ {-{\beta _1}(k)(k_1^2 + k_2^2)} \end{array}} \right) $

(13)

$ {{\mathit{\boldsymbol{\hat u}}}_2}(\mathit{\boldsymbol{k}}) = \frac{1}{{k\sqrt {k_2^2 + k_3^2} }}\left( {\begin{array}{*{20}{c}} {-{\beta _2}(k)(k_2^2 + k_3^2)}\\ {{\alpha _2}(k)k{k_3} + {\beta _2}(k){k_2}{k_1}}\\ {{\beta _2}(k){k_3}{k_1}-{\alpha _2}(k)k{k_2}} \end{array}} \right) $

(14)

$ {{\mathit{\boldsymbol{\hat u}}}_3}(\mathit{\boldsymbol{k}})\frac{1}{{k\sqrt {k_1^2 + k_2^2} }}\left( {\begin{array}{*{20}{c}} {{\beta _3}(k){k_1}{k_2}-{\alpha _3}(k)k{k_3}}\\ {-{\beta _3}(k)(k_3^2 + k_1^2)}\\ {{\alpha _3}(k)k{k_1} + {\beta _3}(k){k_3}{k_2}} \end{array}} \right) $

(15)

式中，α_i(k)和β_i(k)均是形如式(6)、式(7)的随机数，但各自采用不同的随机序列.接着对这3个矢量场作如下平均，得到改进的湍流初始场

$ \mathit{\boldsymbol{\hat u}}(\mathit{\boldsymbol{k}}) = \sqrt {\frac{{E(k)}}{{4\pi {k^2}}}} \frac{{({{\mathit{\boldsymbol{\hat u}}}_1}(\mathit{\boldsymbol{k}}) + {{\mathit{\boldsymbol{\hat u}}}_2}(\mathit{\boldsymbol{k}}) + {{\mathit{\boldsymbol{\hat u}}}_3}(\mathit{\boldsymbol{k}}))}}{{|{{\mathit{\boldsymbol{\hat u}}}_1}(\mathit{\boldsymbol{k}}) + {{\mathit{\boldsymbol{\hat u}}}_2}(\mathit{\boldsymbol{k}}) + {{\mathit{\boldsymbol{\hat u}}}_3}(\mathit{\boldsymbol{k}})|}} $

(16)

我们称这种改进方法为模量平均法.容易看出式(16)能同时满足谱空间中的连续方程(2)和能谱方程(3)，且在谱空间3个方向上统计对称.如果把式(16)代入能谱方程(3)的左边部分，给定能谱可自动满足

$ \begin{array}{*{20}{c}} {\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_{A(k)} {\mathit{\boldsymbol{\widehat u}}\left( \mathit{\boldsymbol{k}} \right) \cdot {{\mathit{\boldsymbol{\widehat u}}}^*}\left( \mathit{\boldsymbol{k}} \right)} {\rm{d}}A(k) = \mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_{A(k)} {\frac{{E(k)}}{{4{\rm{ \mathsf{ π} }}{k^2}}}{\rm{d}}A(k)} = }\\ {\frac{{E(k)}}{{4{\rm{ \mathsf{ π} }}{k^2}}}\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_{A(k)} {{\rm{d}}A(k)} = E(k)} \end{array} $

(17)

图 8是实际计算结果，可以看到式(16)生成的初始场能满足给定能谱.

为了观察该方法的统计性质，采用类似1.3节的做法，我们通过修改随机序列的种子生成1 024个基于模量平均法的初始场，并相应地计算3个棱向D_ll(r)的相对标准差以及U, V, W的PDF相对标准差，结果分别如图 9和图 10 (无量纲脉动速度绝对值小于0.4的部分)所示.通过图 5与图 9、图 7与图 10的对比，可以看出各物理量的样本波动程度都有所减弱.定量的比较如图 11和图 12 (无量纲速度脉动值小于0.4的部分)所示.可以看到，在大多数区域，采用模量平均法能将初始场的相对标准差降低约10%.

3 结论

为了生成均匀各向同性湍流场，我们要求初始场应尽可能接近真实的均匀各向同性湍流.尽管某些表征湍流统计状态的统计量(如表征湍流非平衡性的速度梯度扭率^[24-26]和四阶矩^[27]，表征湍流时空特性的时空相关^[28-29])很难在初始场中加以构造，但至少对于结构函数、速度概率密度这种与能谱相关的简单的统计量，初始场应当尽量各向同性，这应当是评价初始场质量的一个基本要求.初始场将能显著影响湍流非均衡发展的时间^[30].

目前被广泛采用的初始场生成方法由Rogallo提出.在实际应用中，Rogallo方法生成的初始场并不能很好保证结构函数和速度概率密度的各向同性，如果将每个初始场视为一个统计样本，则样本波动较大.在Rogallo方法的基础上我们提出一种模量平均法.与Rogallo方法相比，新方法能分别将结构函数和速度概率密度的相对标准差减小约10%，在一定程度上有所改进.因此我们建议，在数值计算中可以采用模量平均法来取代Rogallo方法生成初始场.