引言

时空关联函数是湍流研究中的一个基本统计特性函数，结构函数、能谱都可以由相关函数得到.无论是湍流基础理论研究还是实际工程应用，时空关联性都有广泛的应用.例如，根据时空关联特性，可以改进大涡模拟方法，发展亚格子模型^[1-2]；根据声比拟理论，流动诱导的声音取决于速度时空关联即频率波数能谱^[3-4].如若追其研究历史，早在1959年Kraichnan^[5]根据直接作用近似理论，提出横扫假设理论，给定了速度时空关联应该满足的关系式. 1972年，Orszag等^[6]在采用direct numerical simulation (DNS)研究衰减均匀各向同性湍流时发现，谱空间的速度时空关联函数可以通过去关联时间尺度进行重整化. 1993年，Gotoh等^[7]采用DNS对湍流速度的欧拉关联和拉格朗日关联做了分析和研究，并对"横扫作用"做了详细的讨论. 2006年，何国威等^[8]在Taylor提出的Frozen-Flow假设基础上，根据流场速度时空关联的数据分析，发现湍流时空信号转化的非线性关系，提出了时空关联椭圆模型，解释了泰勒模型缺少涡变形机制的问题，并将其应用于大涡模拟方法、气动噪声、湍流能谱分析等方面研究.而国内外的一系列关于流场时空关联统计性质的实验^[9-11]、数值和理论研究也说明了其对湍流流场分析的意义.

本文用格子Boltzmann方法研究具有多孔介质壁面的槽道剪切湍流的速度时空关联特性.格子Boltzmann模型是近年来迅猛发展的一种模拟流动的新方法^[12]，该方法遵循质量、动量和能量守恒等基本物理规律，用简单的介观模型在Boltzmann方程控制下模拟流动.该方法经过近二十年的研究日渐成熟，并在微纳米尺度流^[13]、多孔介质流^[14-15]、多相多质流^[16-17]、化学反应流^[18]、非牛顿流体^[19-21]、粒子悬浮流^[22]、磁流体^[23]等方面得到充分应用.与传统计算流体力学（computational fluid dynamics, CFD）相比，该方法适用于模拟弱可压流动问题，且算法简单，边界处理方便，常用于高效并行大规模计算.多孔介质的物理特性（孔隙率、渗透率）不同于一般材料，其能够在人为控制下具有诸如减阻、吸附、调节热传输与流动等功能.流体在多孔介质中的流动在自然界中广泛存在，是多种工程及学科的交叉，涉及石油、天然气、煤层气等各种地下流体资源的勘探和开采工程，固定床反应器、膜反应器（燃料电池）等的设计、优化与放大的化学工程，动物体内血液微循环、植物体内水分的输送等生物医学工程，地面沉降、地下水污染、土壤改良等地下水水文学、土壤学以及土力学等^[24].本文意在通过研究多孔介质壁面槽道湍流的时空关联特性, 一方面进一步探索湍流领域的速度时空关联特性，辅助气动声学研究；另一方面促进对多孔介质在流体领域的研究，辅助对多孔介质减噪减阻性能的探索^[25].

1 数理方法和物理模型 1.1 基于通用渗流模型的LBE模型

LBGK（lattice Bhatnagar-Gross-Krook）模型是应用最为广泛的LBE（lattice Boltzmann equation）模型之一，其中D$d$Q$m$（$d$维空间，$m$个离散速度）系列模型最具代表性^[26].本文采用D3Q19模型来数值模拟槽道湍流.

含作用力项的LBGK模型基本方程为

$\begin{array}{l} {f_i}\left( {{\boldsymbol{x}} + {{\boldsymbol{c}}_i}{\delta _t}, t + {\delta _t}} \right) - {f_i}\left( {{\boldsymbol{x}}, t} \right) = \\ \;\;\; - \frac{1}{\tau }\left( {{f_i}\left( {{\boldsymbol{x}}, t} \right) - f_i^{{\rm{eq}}}\left( {{\boldsymbol{x}}, t} \right)} \right) + {\delta _t}{F_i} \end{array}$

(1)

其中，${f}_{ i}$是离散速度空间的分布函数，${\boldsymbol{x}}$为空间位置矢量，$t$为时间，${\boldsymbol{c}}_i $为离散速度，$\delta _t $为单位时间间隔量，${F}_i $为离散速度空间的与外力有关的项, $\tau $为分布函数趋于平衡态的松弛时间. $\tau $与运动黏性系数$\nu $的关系为

$\begin{equation} \nu=c^2_{\rm s}(\tau-0.5)\delta_t \end{equation}$

(2)

式中，$c_{\rm s}$是格子声速. $f_i^{\rm eq}$是离散速度空间的局部平衡态分布函数，可将其看成是密度和速度的泛函

$\begin{equation} \label{eq3} f_i^{\rm eq} ({\boldsymbol{x}}, t) = \omega _i \rho \left[{1 + \frac{{\boldsymbol{c}}_i \cdot {\boldsymbol{u}}}{c_{\rm s}^2 } + \frac{\left( {{\boldsymbol{c}}_i \cdot {\boldsymbol{u}}} \right)^2}{2c_{\rm s}^4 }-\frac{{\boldsymbol{u}}^2}{2c_{\rm s}^2 }} \right] \end{equation}$

(3)

其中，$\omega _i $是权系数，$\rho $是密度，${\boldsymbol{u}}$是速度矢量，对于含作用力项的LBGK模型，对于D3Q19模型其中权系数$\omega _i $和离散速度${\boldsymbol{c}}_i $定义如下

$\begin{eqnarray} \label{eq4} \omega _i = \left\{ {\begin{array}{ll} \dfrac{1}{18}, &i= 1, 2, \cdots, 6 \\ \dfrac{1}{36}, &i= 7, 8, \cdots, 18 \\ \dfrac{1}{3}, &i= 19 \\ \end{array}} \right. \end{eqnarray}$

(4)

$\begin{eqnarray} {\boldsymbol{c}}_i = \left\{ {\begin{array}{lll} \left( {\pm 1, 0, 0} \right)c, \left( {0, \pm 1, 0} \right)c, &\\ \left( {0, 0, \pm 1} \right)c, &i = 1, 2, \cdots, 6 \\ (\pm 1, \pm 1, 0)c, (\pm 1, 0, \pm 1)c,&\\ (0, \pm 1, \pm 1)c, &i = 7, 8, \cdots, 18 \\ (0, 0, 0)c, &i = 19\\ \end{array}} \right. \end{eqnarray}$

(5)

其中，$c$为格子速度.通过D3Q19模型离散速度如图 1所示.

通过该模型，宏观物理量如密度和动量等可以由下式根据守恒性计算得到

$\begin{eqnarray} \sum\limits_i {f_i^{\rm eq} } = \rho \end{eqnarray}$

(6)

$\begin{eqnarray} \sum\limits_i {{\boldsymbol{c}}_i f_i^{\rm eq} } + \frac{\delta _t }{2}\rho {\boldsymbol{ F}} = \rho {\boldsymbol{u}} \end{eqnarray} $

(7)

为了进一步研究多孔介质槽道湍流，本文在含作用力项的LBGK模型中引入Nithiarasu等^[27]近年提出的渗流模型.该模型不但包含线性阻力和黏性项，而且包含了非线性阻力项，并且可以用于瞬态渗流.其中作用力项${\boldsymbol{F}}$的表达式如下

${\boldsymbol{F}} = - \frac{{\varepsilon v}}{K}{\boldsymbol{u}} - \frac{{\varepsilon {F_\varepsilon }}}{{\sqrt K }}\left| {\boldsymbol{u}} \right|{\boldsymbol{u}} + \varepsilon {\boldsymbol{G}}$

(8)

式中，${\boldsymbol{G}}$是体力，结构函数${F}_\varepsilon$与渗透率$K$与孔隙率$\varepsilon$有关，根据Ergun^[28]经验关联式分别给出了

$\begin{eqnarray} F_\varepsilon = \frac{1.75}{\sqrt {150\varepsilon ^3} } \end{eqnarray}$

(9)

$\begin{eqnarray} K = \frac{\varepsilon ^3d_{\rm p} ^2}{150(1 - \varepsilon )^2} \end{eqnarray}$

(10)

其中，$d_{\rm p}$为固体颗粒的直径.

结合上述方程，本文构造可用于求解多孔介质流动的LBE方程，演化方程形式如方程(1)，而作用力采用式(8).由于LBE的速度是离散的，所以$F_i $不能直接通过连续的作用力项进行离散获得.因此前人发展了多种不同的作用力模型^[29-30]，本文采用郭照立等^[31]提出的G-Z-S作用力模型.该模型可以消除作用力对质量方程、动量方程及时间导数导致的离散误差.则平衡态分布函数和作用力分别为

$\begin{eqnarray} f^{\rm eq}_t=\omega_i\rho\left[1+ \dfrac{{\boldsymbol{c}}_i\cdot{\boldsymbol{u}}}{c^2_{\rm s}}+\dfrac{({\boldsymbol{c}}_i\cdot{\boldsymbol{u}})^2}{2\varepsilon c^4_{\rm s}}-\dfrac{{\boldsymbol{u}}^2}{2\varepsilon c^2_{\rm s}}\right] \end{eqnarray}$

(11)

$\begin{eqnarray} F_i=\omega_i\rho\left(1-\dfrac{1}{2\tau}\right)\left[\dfrac{{\boldsymbol{ F}}\cdot{\boldsymbol{c}}_i}{c^2_{\rm s}}+\dfrac{{\boldsymbol{uF}:({\boldsymbol{c}}_i{\boldsymbol{c}}_i-c^2_{\rm s}{\boldsymbol{I}})}}{\varepsilon c^4_{\rm s}}\right]~~~~~~~~~~ \end{eqnarray}$

(12)

其中，${\boldsymbol{I}}$为单位张量.

基于以上模型方法，本文建立研究问题的物理模型(图 2)：流体在压力驱动下流过两个无限延伸的多孔介质壁面，并呈充分发展的湍流流动状态.该模型计算域在流向（$x$轴方向）、法向（$y$轴方向）、展向（$z$轴方向）分别为$2\pi H$，$2H$，$\pi H$. $H$为半槽道高. 图 2中上下表面外层为不可渗透无滑移的理想外壁面，紧贴其加入均匀各向同性材料的多孔介质壁面.

1.2 速度时空关联函数与椭圆理论

两个不同点处物理量乘积的统计平均称之为相关函数.相关函数表征随机物理量之间相互关联作用的程度，是湍流统计理论的基本概念和研究内容.

速度时间相关函数是指同一个空间点不同时刻两个速度脉动分量乘积的系综平均，即

$\begin{equation} R({\boldsymbol{x}}, \tau ) = \left\langle {{u}'_i \left( {{\boldsymbol{x}}, t} \right){u}'_j \left( {{\boldsymbol{x}}, t + \tau } \right)} \right\rangle \end{equation}$

(13)

同理，空间相关函数是指同一时刻不同空间点的两个速度脉动分量乘积的系综平均，即

$\begin{equation} R({\boldsymbol{r}}, t) = \left\langle {{u}'_i \left( {{\boldsymbol{x}}, t} \right){u}'_j \left( {{\boldsymbol{x}} + {\boldsymbol{r}}, t} \right)} \right\rangle \end{equation}$

(14)

其中，${\boldsymbol{r}}$为连接空间两点的矢量（空间分离距离）、$\tau $为两时刻间隔量（时间分离距离）.而速度时空相关函数

$\begin{equation} R({\boldsymbol{r}}, \tau ) = \left\langle {{u}'_i \left( {{\boldsymbol{x}}, t} \right){u}'_j \left( {{\boldsymbol{x}} + {\boldsymbol{r}}, t + \tau } \right)} \right\rangle \end{equation}$

(15)

本文主要研究流向速度时空关联性，因此可以定义，在不同法向位置，流向速度脉动时空关联系数$C$的数学表达式

$\begin{equation} C\left( {r, \tau ;y} \right) = \frac{\left\langle {u_1 ^\prime \left( {x, y, z;t} \right)u_1 ^\prime \left( {x + r, y, z;t + \tau } \right)} \right\rangle }{\left\langle {u_1 ^\prime \left( {x, y, z;t} \right)^2} \right\rangle } \end{equation}$

(16)

由于在流向、展向的脉动量呈统计均匀分布，且沿时间进程流动统计定常, 因此速度时空关联与时间$t$和流向、展向的坐标$x$，$z$无关.

基于二阶近似，何国威等^[8]延伸了泰勒冻结假设与横扫理论，针对剪切湍流提出时空关联椭圆模型，即有

$\begin{eqnarray} C\left( {r, \tau } \right) = C\left( {r_E, 0} \right) \end{eqnarray}$

(17)

$\begin{eqnarray} r_{\rm E} ^2 = \left( {r - U\tau } \right)^2 + V^2\tau ^2 \end{eqnarray}$

(18)

其中，$r_{\rm E} $为椭圆模型的特征长度，传播速度$U$和横扫速度$V$分别为椭圆模型的两个特征速度

$\begin{eqnarray} U = - \frac{\partial ^2C\left( {0, 0} \right)}{\partial r\partial \tau }\left( {\frac{\partial ^2C\left( {0, 0} \right)}{\partial r^2}} \right)^{ - 1} \end{eqnarray}$

(19)

$\begin{eqnarray} V^2 = \frac{\partial ^2C(0, 0)}{\partial \tau ^2}\left( {\frac{\partial ^2C(0, 0)}{\partial r^2}} \right)^{ - 1} - U^2 \end{eqnarray}$

(20)

其中，$U$和$V$也表征时空关联等值图的几何特性.设斜椭圆(图 3)的倾斜角为$\alpha $，椭圆的半短轴$b$与半长轴$a$之比为椭圆率$\lambda $，则有如下关系

$\begin{array}{l} {\tan ^2}\alpha = \\ \;\;\;\;\;\;\frac{{4{U^2}}}{{{{\left[{\sqrt {{{\left( {1 + {U^2}-{V^2}} \right)}^2} + 4{U^2}{V^2}} + \left( {1-{U^2}-{V^2}} \right)} \right]}^2}}} \end{array}$

(21)

$\begin{array}{l} {\lambda ^2} \equiv {\left( {\frac{b}{a}} \right)^2} = \\ \;\;\;\;\;\;\;\frac{{4{V^2}}}{{{{\left[{\sqrt {{{\left( {1 + {U^2}-{V^2}} \right)}^2} + 4{U^2}{V^2}} + \left( {1 + {U^2} + {V^2}} \right)} \right]}^2}}} \end{array}$

(22)

2 数值模拟分析 2.1 计算验证

本文在流向、法向和展向分别采用$280~\times~90~\times~140$的均匀网格，用壁面摩擦速度和半槽道高定义的剪切雷诺数$Re_\tau=180$，$Re_\tau=u_\tau H\mbox{ / }\nu $，其中$u_\tau $为壁面速度.槽道中心剪切雷诺数为$Re_{\rm c}=u_{\rm c}H/\nu$，其中$u_{\rm c} $为中心线速度，计算参数见表 1.

为了验证计算方法和程序的准确性，将计算的速度的相关统计量结果与Kim等^[32]的DNS结果进行对比，如图 4所示流向平均速度和三个方向脉动速度分布，以及如图 5所示的剪切应力和雷诺应力分布.可以看出，本文采用LBM（lattice Boltzmann method）数值计算结果与前人结果符合良好.

2.2 脉动速度两时间关联

为了研究速度时间关联特性，分别计算法向位置$y^ +=16$时，常规壁面剪切湍流在空间固定距离差为$r^ +=0$，$r^ +=36$，$r^ +=72$，$r^ +=108$，$r^ +=144$下的时间关联函数，如图 6所示.随着空间分离距离的增大，其时间关联函数的峰值往右移且峰值从1开始降低.

在不同空间分离距离的时间关联函数基础上，采用数值模拟获得的传播速度$U$和横扫速度$V$定义新的特征尺度$r_{\rm E} $

$\begin{equation} r_{\rm E} = \sqrt {\left( {r^ + - U\tau ^ + } \right)^2 + V^2(\tau^+)^2} \end{equation}$

(23)

并以此作为横坐标对时间关联函数(图 6)进行重整化，得到图 7.从图 7可见，本文所计算的各个不同固定空间分离距离差之下的时间关联在对横坐标处理的情况下，几乎吻合在一起.这说明当以特征尺度$r_{\rm E} $归一化后，不同距离差下的两时间关联函数具有相似性.

为了进一步研究各个法向位置处的解相关时间，计算了不同法向位置处、不同空间分离距离的时间关联函数，得到图 8.可以看出，随着法向位置逐渐远离壁面，在同样的时间分离距离（横坐标）下，其时间关联函数值减小，同时其解相关的时间也越来越短，高值相关性越来越集中.由于壁面附近存在大尺度的拟序结构，随着发卡涡的破裂，会形成强烈的指向壁面的下扫运动并将一部分流体抛往槽道中心方向，因此壁面附近流体流动的相关性是最强的，越远离壁面，受流向流体的影响，大尺度涡逐渐破裂，相关性也就随之减弱.

基于图 8中的不同法向位置的时间关联函数，进一步研究在不同法向位置处的去关联尺度. 图 9的纵坐标为去关联时间$\tau _0 $，横坐标为空间分离距离.从图 9可以看出，在法向方向上，远离壁面位置越远，其流向速度相关性越集中，去关联时间越短.

在多孔介质壁面中，孔隙率和渗透率是多孔介质壁面的特征参数.因此，使用控制变量法，基于速度时间关联函数，通过对比分析不同孔隙率和渗透率，考虑这两个参数对槽道湍流时间关联的影响.其中渗透率以达西数$Da$来表征

$\begin{equation} Da = \frac{K}{H^2} \end{equation}$

(24)

其中，$K$为渗透率，$H$为半槽道高.

首先研究孔隙率$\varepsilon $.计算了在$Da=5\times 10^{-4}$，$y^ +=8$处，孔隙率分别为0.3，0.4，0.5，0.6，0.7和0.8时的时间关联函数，如图 10(a)所示.红色实线为常规壁面槽道剪切湍流计算结果.从图 10(a)可以看出，随着孔隙率的增大，其关联系数也呈现递增趋势.也就是说，对于多孔介质壁面来说，粗糙度可以降低速度关联性.同时，本文研究了距离不同法向位置的关联性问题，因此计算了$y^ +=8$，$y^ +=96$，$y^ +=180$处的不同孔隙率下的时间关联图像，如图 10所示.

从图 10可以发现三个函数特性：其一，在所研究的参数范围内，在不同的法向位置，都保持着随着孔隙率增大而关联性增强的关系.其二，随着法向位置远离壁面，速度关联性有所降低，而且高值关联系数更加集中.其三，随着$y^ + $递增，多孔介质的时间关联函数图线与常规槽道湍流越来越接近，在槽道中央位置的地方，多孔介质与常规槽道湍流的时间关联比较相似.其实说明了壁面多孔介质的影响主要在近壁面区域，而对于远离壁面的位置，其影响能力减弱.

渗透率（$Da$数）对时间关联性的影响如图 11所示，当孔隙率$\varepsilon$分别为$ 0.4$和$0.7$时，随着达西数的增大，关联函数值减小.由于达西数与多孔介质可渗透度成正比，达西数越大渗透率越强.因此在达西数较小的情况下，亦是渗透率较低的情况下，时间关联性较强.随着达西数减小，多孔介质壁面渗透率强度减弱，近壁面的速度脉动增大，而关联性反而增大.从机理上面分析，多孔介质对壁面相关性的主要影响是壁面附近的拟序结构，随着孔隙率的增大和达西数的减小，多孔介质对拟序结构的抑制作用增强，拟序结构的强度减弱，表现为壁面附近的相关性有所减弱.越远离壁面，多孔介质的影响变得越不明显，这是因为越靠近槽道中央，大尺度的湍流结构越稀少，且相关性本来就很小.

2.3 基于LBM的速度时空关联与椭圆模型

根据速度时空关联函数，统计常规槽道剪切湍流在$y^ +=16$处流向速度时空关联性数值，如图 12所示.横坐标为时间分离距离，纵坐标为空间分离距离，等值线为速度关联系数.从图 12可以看出，该函数等值线分布呈现类椭圆特征形状.其次，为了将计算结果与前人研究结果对比，本文在通过LBM方法模拟槽道剪切湍流速度并获得速度时空关联函数基础上，通过测量$\alpha $和$\lambda $，计算$U$和$V$，得到图 13和图 14.并与谱方法DNS获得的研究结果比较，发现基于LBM所获得的壁面剪切湍流与DNS以及椭圆模型的特征速度理论值是较为符合的.

2.4 多孔介质槽道湍流的速度时空关联

基于之前所研究的多孔介质槽道湍流时间关联的特性结果，本文对多孔材料壁面参数取值，对比分析时空关联性在常规壁面与多孔介质壁面间的差别.

为了更加明显地对比分析差别，以常规槽道湍流作为参考值.首先，分析不同孔隙率下的速度时空关联函数.其次，参考时间关联分析的结果，采用对比更加明显的$y^ +=8$法向位置的时空关联数据，其中$Da=5\times 10^{-4}$，如图 15所示（虚线为常规壁面参考值）.发现随着孔隙率不断增大，多孔介质的时空关联等值线越接近常规槽道的时空关联等值线.通过对比参考值，可知在本文所研究范围内的多孔介质，时空关联系数随着孔隙率的增大而增大.

对比常规壁面与多孔壁面可以看出，相对于椭圆率的变化，倾斜角的变化非常小.多孔壁面主要影响槽道湍流速度时空关联的椭圆率，而倾斜角变化不大.由于椭圆率主要与横扫速度有关，推断多孔介质壁面主要影响剪切湍流的横扫速度.

为研究不同渗透率下的速度时空关联函数.统一设置条件为$y^ +=8$, $\varepsilon=0.4$，以常规壁面的时空关联为参考值（虚线为参考值），如图 16所示.可以定性地判断，随着渗透率减小，时空关联性增强.

多孔壁面对槽道剪切湍流的流向速度时空关联函数的影响如图 16所示.从图 15和图 16可以看出，倾斜角相差都不大，而椭圆率相差较大，因此多孔介质主要影响椭圆关联等值线的椭圆率，对其倾斜角影响较小.

3 结论

本文基于Boltzmann方法，获取流场信息，并着重分析了常规槽道湍流与多孔介质槽道湍流的流向速度脉动时空关联特性.得到以下主要结论：

（1）通过格子Boltzmann方法计算得到的槽道剪切湍流流场，研究其速度脉动时空关联函数的统计特性，发现关联系数等值线呈现出近似椭圆的典型特征.

（2）通过分析不同参数的多孔介质壁面剪切湍流，在所研究参数范围内，发现多孔介质壁面的孔隙率、渗透率对近壁区的速度脉动时空关联性有明显影响：随着孔隙率增大，流向速度脉动关联性增强；随着渗透率增大，流向速度脉动关联性降低.

（3）通过研究不同法向位置的流向速度脉动时空关联函数，发现越远离壁面的法向位置（对数律区和中心区）其速度时空关联函数所呈现的椭圆等值线越细长，高值关联值越集中；同时越远离壁面位置，其速度关联性受到多孔介质壁面影响越小.根据关联函数椭圆率的变化，说明多孔介质壁面主要影响剪切湍流的横扫速度.