湍流现象出现在与流体力学相关的各种学科问题中，如大坝排水、大气流动、机翼尾流等. 湍流流动在工业设备中较为常见，如何有效实现湍流流动的控制，改善设备性能，一直以来备受相关领域的学者关注. 利用仿生学思想，设计工业元件来控制流动的方法已经成为趋势^[1]. 在仿生学中Mandelbrot^[2]将分形定义为具有某种自相似性结构、功能和性质的集合. 因此分形流动(fractal-generated turbulent flow)泛指流体（液体或气体）绕/经过具有分形外形的物体时引发的紊乱、不规则流动. 流体流经分形结构或规则结构网格的现象被称为分形网格湍流(fractal grid turbulence, FGT)或规则网格湍流(regular grid turbulence, RGT).

国内外学者开展了一系列有关分形结构流场的实验和数值模拟研究. 在实验方面，Ferko等^[3]利用压电俘能器装置回收分形流场中的能量，研究表明与常规结构相比，利用分形结构可以使其能量转化的效率显著增加. Steiros等^[4]研究发现在搅拌器中利用分形刀刃，可以有效降低阻力系数. Sponfeldner等^[5]将分形结构应用在燃烧实验中，研究表明与规则湍流场相比，分形湍流场中的火焰燃烧更稳定以及燃烧速率更大. Comte-Bellot等^[6]利用风洞实验来探索改善网格湍流各向同性的方法. Antonia等^[7]对网格湍流场的涡量变化进行了深入研究. Seoud等^[8]和Mazellier等^[9]开展了分形流场机理研究，结果表明利用分形流场，将能量直接注入不同大小尺度的湍流结构中，可以改变流场中能量传递机制和湍动能耗散特性. Lavoie等^[10-11]的实验表明FGT的统计特性受到初始条件的影响较大. 严磊等^[12]利用风洞实验发现风参数变化主要与测试点和格栅断面间距、格栅板条厚度和单元格栅边长有关. 周蓉^[13]在实验中发现风洞流向脉动风速在不同水平间距和湍流积分尺度下的相关性与两者比值相关，比值越大，相关性越弱.

研究人员也利用直接数值模拟(direct numerical simulation, DNS)方法进行了有关方面的积极探索. Mazzi等^[14]研究了三维不可压缩湍流受幂律多尺度作用力影响的周期性和统计稳定的行为. Laizet等^[15-18]采用DNS仿真开展了有关FGT和RGT的研究，结果表明分形湍流可以显著提高流体混合效率、增大标量的传递，并在此基础上提出了一种空间尺度展开(space-scale unfolding, SSU)机制. Nagata等^[19]对分形流场空间演化的DNS模拟结果，证实了在靠近网格的近场区域存在湍流产生区，远离网格的区域存在湍流衰减区. Suzuki等^[20]模拟并分析了规则网格湍流特性，发现即便是在流场远下游处( $ X/M \approx 100 $，X为规则网格湍流空间演化时沿流向的距离，M为规则网格的尺寸大小)，大尺度脉动速度场也并非是完全各向同性的，各向异性约为17%.

据文献统计，目前关于FGT的研究方法主要是实验与数值模拟，并且绝大部分都集中在空间演化流动(spatially evolving turbulent flow). 网格湍流空间演化的数值模拟在计算成本上存在较大的劣势，而时间演化湍流(temporally evolving turbulent flow)指的是，尽管流动在流向上是非均匀的，但仍然采用周期性边界条件，这样可以大幅缩小流向长度，节约计算成本. 近些年也有部分学者将此方法应用于湍流各领域的研究，如射流^[21-22]、尾流^[23]、混合层^[24-25]和边界层^[26]等流动. Watanabe等^[27]也对RGT开展过时间演化的DNS研究，但至今缺乏对分形网格湍流进行时间演化模拟研究. 因此本研究基于伪谱Fourier-Galerkin方法，布置三阶正方形分形网格流场，开展随时间演化的DNS仿真，探讨FGT相关统计特性随时间演化的特征. 并且开展了具有相同阻塞率的RGT的数值模拟，以进行对比分析.

DNS仿真是基于伪谱-伽辽金方法^[28]进行空间离散实现的. 控制方程为连续性方程和不可压缩Navier-Stokes (NS)方程组，表达式如下：

(1) $ \nabla \cdot {\boldsymbol{u}} = 0 , $

(2) $ \frac{{\partial {\boldsymbol{u}}}}{{\partial t}} - {\boldsymbol{u}} \times {\boldsymbol{\omega}} = \nu {\nabla ^2}{\boldsymbol{u}} - \nabla P , $

(3) $ {\boldsymbol{u}}({\boldsymbol{x}}+2{{\text{π}}} {{\boldsymbol{e}}^i},t) = {\boldsymbol{u}}({\boldsymbol{x}},t);\;i = 1,2,3 , $

(4) $ {\boldsymbol{u}}({\boldsymbol{x}},0) = {{\boldsymbol{u}}_0}({\boldsymbol{x}}) . $

式中： ${\boldsymbol{u}}({\boldsymbol{x}},t)$为瞬时速度矢量， ${\boldsymbol{\omega}} = \nabla \times {\boldsymbol{u}}$为涡量矢量， ${{\boldsymbol{e}}^i}$为笛卡尔坐标单位矢量， $P = p+{\boldsymbol{u}} \cdot {\boldsymbol{u}}/2$为无量纲化后的压力， $ \nu $为运动黏度. X、Y、Z分别表示计算域的3个方向. 流场Y-Z平面网格效果如图1所示，规则网格条的横向厚度 $ d = 0.2M $，M为规则正方形网格条之间的长度尺寸，同时等于分形网格第3阶正方形长度尺寸M₂.

流场计算域利用均匀网格进行空间离散，其网格节点在物理空间表达式如下：

(5) $ \begin{split} {\boldsymbol{x}} = [x,y,z] =\;& \left\{ [{x_i},{y_j},{z_k}] = \left[\frac{{2{\text{π}} i}}{{{N_X}}},\frac{{2{\text{π}} j}}{{{N_Y}}},\frac{{2{\text{π}} k}}{{{N_Z}}}\right];\right. \\ \;&i,j,k = 0,\cdots ,{N_i} - 1 \Bigg\}. \end{split} $

在伪谱-伽辽金方法中，要求将所有的变量从物理网格转换到离散的傅立叶波数网格进行计算. 在傅立叶空间（谱空间或波数空间）中，假设在流场物理空间中，沿流向、法向及展向均有N个网格节点，即N_X=N_Y=N_Z，则波数网格表达式如下：

(6) $ {\boldsymbol{k}} = [{k_x},{k_y},{k_z}] = \left\{ {[l,m,n];l,m,n = - \frac{N}{2}+1,\cdots ,\frac{N}{2}} \right\} .$

对速度场进行傅立叶变换(Fourier transform)，用 ${F}$和 ${{F}^{ - 1}}$分别表示傅立叶正向变换和其逆向变换过程，则速度在傅立叶空间和物理空间的具体形式分别如下：

(7) $ {{\hat {\boldsymbol{u}}}_k}(t) = \sum\limits_{\boldsymbol{x}} {{{\boldsymbol{u}}}({\boldsymbol{x}},t){{\rm{e}}^{ - l{\boldsymbol{k}} \cdot {\boldsymbol{x}}}} = } {F}({{\boldsymbol{u}}}({\boldsymbol{x}},t)) , $

(8) $ {\boldsymbol{u}}({\boldsymbol{x}},t) = \sum\limits_{\boldsymbol{k}} {{{{\hat {\boldsymbol{u}}}}_k}(t){{\rm{e}}^{l{\boldsymbol{k}}\cdot {\boldsymbol{x}}}}/{N^3} = } {{F}^{ - 1}}({{\hat {\boldsymbol{u}}}_k}(t)) . $

式中： ${{\hat {\boldsymbol{u}}}_k}(t)$为傅立叶系数， $\mathrm{e}^{l{\boldsymbol{k}}\cdot {\boldsymbol{x}}}$为伪谱-伽辽金方法的基础函数，k为波数矢量， $ l = \sqrt { - 1} $为虚数单位. 其中傅立叶正变换和逆变换都涉及3个连续的变换，每个周期方向1个. 则傅立叶空间的控制方程组表达式如下：

(9) $ \frac{{{\rm{d}}{{\hat {\boldsymbol{u}}}_k}}}{{{\rm{d}}t}} - {(\widehat {{\boldsymbol{u}} \times {\boldsymbol{\omega}} })_k} = - \nu {\left| {\boldsymbol{k}} \right|^2}{\hat {\boldsymbol{u}}_k} - l{\boldsymbol{k}}{\hat P_k} , $

(10) $ l{\boldsymbol{k}} \cdot {\hat {\boldsymbol{u}}_k} = 0 . $

压力项通对式(9)中的压力点乘 $l{\boldsymbol{k}}$重新整理，其表达形式为 ${\hat P_k} = - l{\boldsymbol{k}} \cdot {(\widehat {{\boldsymbol{u}} \times {\boldsymbol{\omega}} })_k}/|{\boldsymbol{k}}{|^2}$. 最终不可压缩NS方程组在傅立叶空间表达式如下：

(11) $ \frac{{{\rm{d}}{{\hat {\boldsymbol{u}}}_k}}}{{{\rm{d}}t}} = {(\widehat {{\boldsymbol{u}} \times {\boldsymbol{\omega}} })_k} - \nu {\left| {\boldsymbol{k}} \right|^2}{\hat {\boldsymbol{u}}_k} - {\boldsymbol{k}}\frac{{{\boldsymbol{k}} \cdot {{(\widehat {{\boldsymbol{u}} \times {\boldsymbol{\omega}} })}_k}}}{{{{\left| {\boldsymbol{k}} \right|}^2}}} . $

考虑分形网格(fractal grid)流场和规则网格(regular grid)流场2组工况进行对比研究. 具体参数设置如表1、2所示. T_end为数值模拟结束的时间. 计算域3个方向的长度为 ${L_{X}} = {L_{Y}} = {L_{Z}} = 8M{\text{ = }} 8{M_2}$， $ {N_{\text{f}}} $表示分形网格在Y-Z平面的几何维度(N_f= $3 $， $ i = 0,1,2 $). 阻塞率的定义为 $ \sigma = {T_{{\text{bar}}}}/({L_{{Y}}}{L_{{Z}}}) $， $ {T_{{\text{bar}}}} $为网格条的总面积（减去重复区域）. 雷诺数基于初始平均速度 $ {U_0} $和网格尺寸M， $ R{e_{\text{M}}} = {U_0}M/\nu $. 由于须与规则格栅进行对比分析，同时为了研究小尺度结构的运动规律，须保证分形格栅具有比规则格栅更小的网格. 在相同阻塞率 $ \mathrm{\sigma }=0.36 $的条件下，通过计算得到分形流场采用三阶正方形( ${{N}}_{\mathrm{f}}= 3$)分形网格较为合适，而正方形网格也是目前研究最多，流动特性最为特殊的分形模式. 另外，Nagata等^[29]和Laizet等^[30]的研究表明，在法向/展向的最小网格条上最少须保证3个节点才能够精确模拟出合适的小尺度尾流迹象，本研究分形格栅的最小网格条上有5个节点 ${N}=5$，已经可以满足精度要求. 因此，本研究对象设置为三阶正方形分形网格. 其中每一阶网格的长度尺寸为 $ {M_i} = R_{\text{M}}^i{M_0} $， $ {M_0} $、 $ {M_1} $和 $ {M_2} $分别表示一阶、二阶和三阶分形网格尺寸，其中一阶网格尺寸 ${M_0} = 0.5{L_{Y}}$，长度比系数 $ {R_{\text{M}}} = 0.5 $. 网格条横向厚度 $ {d_i} = R_{\text{d}}^i{d_0} $， $ {d_0} $、 $ {d_1} $和 $ {d_2} $分别表示一阶、二阶和三阶分形网格横向厚度，其中厚度比系数 $ {R_{\text{d}}} \approx 0.320 $. 网格Y-Z平面的具体效果可在图1中观察. 分形网格的横向厚度比 $ {D_{\text{r}}} $定义为：构成最大正方形网格的横向厚度与最小正方形网格横向厚度的比值，即 $ {D_{\text{r}}} = {d_0}/{d_2} $. 计算域的网格节点数量为 ${N_{X}} {N_{Y}} {N_{Z}} = 800 \times 400 \times 400$.

DNS仿真采用PYTHON语言编写程序实现，数据后处理部分采用FORTRAN、MATLAB语言编写程序完成. 由于DNS数据信息量大，利用更加高效的HDF5格式存储和处理数据，以节约计算和存储成本. 计算域的3个方向均采用周期性边界条件，时间递进采用四阶Runge-Kutta方法求解，计算时间步长为 $ \Delta t/(M/{U_0}){\text{ = }}0.005 $，根据库朗数(Courant Friedrichs Lewy, $ {\text{CFL = }}\Delta tU/\Delta X $)定义可以确保数值方法稳定性，以及时间递进精度满足要求.

在式(11)中，采用谱方法求解对流项、黏性项和压力项. 其中求解对流项须先将速度和涡量转换到物理空间进行叉乘，然后将 $ {\boldsymbol{u}} \times {\boldsymbol{\omega}} $转换回傅立叶空间进行计算. 因此傅立叶空间的涡量矢量和对流项计算式分别如下：

(12) $ {F}(\nabla \times {\boldsymbol{u}}) = {\hat {\boldsymbol{\omega}} _k} = l{\boldsymbol{k}} \times {\hat {\boldsymbol{u}}_k} , $

(13) $ {(\widehat {{\boldsymbol{u}} \times {\boldsymbol{\omega}} })_k} ={F}({\boldsymbol{u}} \times {\boldsymbol{\omega}} ) = {F}({\boldsymbol{u}} \times {{F}^{ - 1}}(l{\boldsymbol{k}} \times {\hat {\boldsymbol{u}}_k})) . $

设置网格湍流时间模拟的初始条件为

(14) $ u = {\langle U\rangle _x}+u' , v = v' , w = w' . $

式中：u、v和w分别为X、Y、Z这3个方向的瞬时速度， $ u' $、 $ v' $和 $ w' $为对应的脉动速度，符号“ $ {\langle \cdot \rangle _x} $”表示对统计均匀方向（X方向）进行平均. 初始速度沿X方向均匀分布，其表达式如下：

(15) $ {\langle U\rangle }_{x}=\left\{\begin{array}{cc}{U}_{0},& 网格条的位置;\\ 0,& 其他位置.\end{array}\right. $

时间模拟的网格湍流相当于以初速度 $ {U_0} $拖曳网格条，由网格条产生的尾流湍流合并而来. 其时间演化的距离 $ X = t{U_0} $，相当于风洞实验中空间演化在流向上远离网格的距离 $ X $. 因此本研究无量纲时间采用 $ T/{T_{\text{r}}} = X/M $(其中 $ {T_{\text{r}}} = M/{U_0} $)，即风洞实验中无量纲的流向距离 $ X/M $. 网格湍流具体效果如图2所示.

由于时间模拟的网格湍流初始是层流状态，利用一定强度的初始脉动速度可以有效促进网格条产生的尾流从层流转变为湍流状态. 本研究采用Kempf等^[31]提出的扩散方法获得式(14)中的初始脉动场，其均方根(root-mean-square, RMS)速度约为 $ 0.5\text{%} {U_0} $，特征长度尺度约为0.2M. 该脉动速度产生的效果相当于风洞实验中网格尾流湍流，脉动速度强度与部分风洞实验^[8,32]报道的尾流湍流强度基本一样.

为了验证数值算法的精度可靠性，利用经典的二维Taylor-Green Vortex (TGV)数值解与解析解的误差大小作为验证方案. 二维TGV控制方程组如下：

(16) $ \frac{{\partial {\boldsymbol{u}}}}{{\partial t}}+ \nabla \cdot ({\boldsymbol{u}}{\boldsymbol{u}}) = - \nabla P+\nabla \cdot (\nu {\nabla }{\boldsymbol{u}}) , $

(17) $ \nabla \cdot {\boldsymbol{u}} = 0 . $

由于不存在外力加持，系统的能量会逐渐耗散衰减. X和Y方向均采用周期性边界条件，计算域范围为 $0 \leqslant X,Y \leqslant 2{\text{π}}$. 雷诺数设置为 $Re = 1\;600$，无量纲时间步长 $ \Delta t = 0.001 $. 模拟2种不同节点数的工况Run1: $ {N_{2{\text{D}}}}{\text{ = }}{32^2} $和Run2: $ {N_{2{\text{D}}}}{\text{ = }}{64^2} $. 二维TGV流动数值解的初始条件如下：

(18) $ \left. \begin{array}{l} u(x,y,t = 0) =\; \sin \; (2{\text{π}} x)\cos\;(2{\text{π}} y), \\ v(x,y,t = 0) =\; - \cos\;(2{\text{π}} x)\sin\;(2{\text{π}} y) . \end{array} \right\} $

与之对应的无外力二维不可压缩TGV的解析解表达式如下：

(19) $ \left. \begin{array}{l} u(x,y,t) =\; \sin\; (2{\text{π}} x)\cos\;(2{\text{π}} y){{\text{e}}^{ - 2\nu t}} , \\ v(x,y,t) = \; - \cos\;(2{\text{π}} x)\sin\;(2{\text{π}} y){{\text{e}}^{ - 2\nu t}}. \end{array} \right\} $

二维TGV数值模拟结果与解析解的均方根误差表达式^[33]为

(20) $ \left. \begin{aligned} {e_{{{\rm{RMS}}}}} = \;& {\left[ {\frac{1}{{{N_X}}}\frac{1}{{{N_Y}}}\sum\limits_{i = 1}^{{N_X}} {\sum\limits_{j = 1}^{{N_Y}} {{{({u_{\rm{n}}}({x_i},{y_j},t) - {u_{\rm{a}}}({x_i},{y_j},t))}^2}} } } \right]^{1/2}} ; \\ \;&{x_i} = (i - 1)/{N_{X}},\;{y_j} = (j - 1)/{N_{Y}} . \end{aligned} \right\}$

式中：下标n和a分别代表数值解和解析解.

为了明确数值解与解析解之间的误差大小，两者的均方根误差如图3(a)所示. 如图3(b)所示为最大误差e_max，计算式参考Sengupta等^[34]的方法，即 ${\text{max}}\;(\left| {{{\boldsymbol{u}}_{\text{n}}} - {{\boldsymbol{u}}_{\text{a}}}} \right|)$. 在无量纲时间 $ \tau \leqslant 150 $时，Run1和Run2的误差大小基本保持一致，但随着时间的增加，Run2的误差要小于Run1的. 当计算时间步数达到200 000时，Run2的最大误差仍能保证在约10⁻⁹. 由此可见本研究的数值方法满足计算精度要求.

网格湍流的时间模拟须在高分辨率下才能更好地捕捉流场主要特征的变化. 通常主要集中观察网格中心线上统计信息的变化. 如图4所示为分形网格和规则网格中心线上空间分辨率 $ {(\Delta X\Delta Y\Delta Z)^{1/3}}/\eta $随时间演化结果. 其中 $ \eta {\text{ = }}{({\nu ^3}/\varepsilon )^{1/4}} $为Kolmogorov长度尺度，是湍流场中最小的尺度空间分辨率，计算位置为虚拟中心(Y, Z)=(0, 0)， $ \varepsilon {\text{ = }}2\nu {\langle {s_{ij}}{s_{ij}}\rangle _x} $为沿网格中心线发展的瞬时耗散率. $ {s_{ij}} = (\partial {u_j}/\partial {x_i}+ \partial {u_i}/\partial {x_j})/2 $表示瞬时应变率张量. 分形流场中心线上最大空间分辨率为0.65，表现在 $ T \approx 50{T_{\text{r}}} $处. 规则流场中心线上最大空间分辨率为1.2，该值出现在 $ T \approx {\text{6}}{T_{\text{r}}} $处. 在网格湍流大部分的衰减阶段，网格尺寸 $ {(\Delta X\Delta Y\Delta Z)^{1/3}} $小于 $ \eta $. Laizet等^[35]的研究表明，当流场空间分辨率最差为7时，数值模拟结果与实验值相比误差仍能保证在约10%. Moin等^[36]也指出，当网格尺寸大小达到 $ 1.5\eta $时，利用谱方法进行DNS仿真可以精确捕捉到大部分小尺度耗散，同时可以将数值计算误差控制在5%以内. 由此可以认为本研究计算精度满足空间分辨率要求.

如图5、6所示分别展示了FGT和RGT不同时刻下流向(X方向)瞬时速度三维可视化结果. 如图5(a)、6(a)所示分别为2种网格湍流的初始状态( $ T = 0 $). 初始平均速度的影响在图5(b)~(f)和图6(b)~(e)中表现明显，平均速度产生的切应力导致网格条的尾流随时间从层流加速转变为湍流状态. 从图5(g)开始，FGT中大尺度网格产生的尾流相互作用逐渐增强，流场湍流混合度增大. Laizet等^[17-18]关于网格湍流空间演化研究表明随着分形网格横向厚度比 $ {D_{\text{r}}} $的增大，FGT需要演化更长的空间距离 $ X/M $（时间模拟的 $ T/{T_{\text{r}}} $），流向瞬时速度才能达到与RGT中相同的状态. 本研究分形网格的 $ {D_{\text{r}}}{\text{ = }}9.5 $，由图5、6对比也可以发现，FGT要演化更长的时间才能达到与RGT相同的湍流状态. 由图6(f)开始，规则网格条的尾流相互作用逐渐合并形成网格湍流，并随时间逐渐发展为近似统计均匀状态.

如图7、8所示分别为FGT和RGT的流向平均速度 $ {\langle U\rangle _x} $(本研究网格中心线上 $ {\langle U\rangle _x} $采用U_c表示，下标c表示中心线)随时间演化的可视化结果. 可以看出，由于初始平均速度剖面的不均匀性，统计量在初始阶段表现为高度不均匀状态. 与RGT相比，由于FGT受到初始条件影响更久，其流向平均速度分布与空间的相关性减弱的更慢，但随着时间演化，2类网格湍流终将会达到统计均匀状态. Zhou等^[37]的研究也表明，随着网格湍流空间演化的流向距离 $ X/M $（时间演化的 $ T/{T_{\text{r}}} $）增加，网格湍流中的小尺度将逐渐演化为类似于各向均匀的状态.

如图9所示为FGT在 $ T = 25{T_{\text{r}}} $时，某流向平面位置(− $4M \leqslant Y,\;Z \leqslant 4M$， $ X = 4M $)的瞬时涡量图，并给出了速度梯度张量第二不变量Q. 如图10所示为RGT对应的可视化结果. 速度梯度的张量 ${{\boldsymbol{A}}_{ij}} = \partial {{\boldsymbol{u}}_i}/\partial {{\boldsymbol{x}}_j}$包括2部分，应变率张量 ${{\boldsymbol{s}}_{ij}} = (\partial {{\boldsymbol{u}}_i}/\partial {{\boldsymbol{x}}_j}+ \partial {{\boldsymbol{u}}_j}/\partial {{\boldsymbol{x}}_i})/2$和涡张量 ${{\boldsymbol{\omega}} _{ij}} = (\partial {{\boldsymbol{u}}_i}/\partial {{\boldsymbol{x}}_j} - \partial {{\boldsymbol{u}}_j}/\partial {{\boldsymbol{x}}_i})/2$，速度梯度张量第二不变量 ${\boldsymbol{Q}} = {{\boldsymbol{Q}}_\omega }+{{\boldsymbol{Q}}_s}{\text{ = }}({{\boldsymbol{\omega}} _{ij}}{{\boldsymbol{\omega}} _{ij}}/2 - {{\boldsymbol{s}}_{ij}}{{\boldsymbol{s}}_{ij}}/2)$. 正的Q值与涡应变占优势的区域有关，而负的Q值出现在与应力应变占优势的区域. 由图9、10可以看出，在此计算域中充斥着类似于网格湍流空间演化的各种尺度结构. 须注意的是，尽管都处于相同的初始时刻( $ T{\text{ = }}25{T_{\text{r}}} $)，FGT依然能识别出最大尺度湍流结构的形状，它起到了强化湍流的作用，导致FGT比RGT产生了更高的涡度值，而RGT此时明显已经进入衰减耗散阶段. 从本质上讲，这2种不同类型的网格湍流是以不同的方式产生的. 如图10所示的规则网格湍流是相同大小的尾流在几个网格内相互作用，并在网格附近以均匀的方式混合在一起而形成的. 图9分形网格的情况下，网格上最小的网格条产生最小的尾流，这些尾流在初始阶段相遇并混合在一起, 而较大的网格条产生较大的尾流，这些尾流将在随时间演化的过程中相遇并混合. 这一现象在最小到最大的湍流尺度上重复进行. Mazellier等^[9]与Laizet等^[38]在分形网格湍流空间演化研究中也得到了类似的结论.

由于缺乏分形网格湍流时间演化的参考案例，选取前人^[39-42]关于规则网格湍流的实验数据和数值模拟结果与本研究RGT工况进行对比，以此验证计算结果的准确性. 规则网格是由多个长度为M的正方形格子重复组成，取其中的 $ {P_1} $(虚拟中心Y/M=0，Z/M=0)、 $ {P_2} $(Y/M=0.5，Z/M=0)和 $ {P_3} $(Y/M=0.5，Z/M=0.5)点分析统计特性，平面位置如图1所示. 如图11所示为RGT中Y-Z平面 $ {P_1} $、 $ {P_2} $和 $ {P_3} $位置的流向平均速度 $ {\langle U\rangle _x} $和均方根速度(又称湍流强度) $ {u_{{\text{rms}}}}{\text{ = }}\langle u{'^2}\rangle _x^{1/2} $随时间演化结果. 图11(a)中，本研究设置的雷诺数 $ R{e_{\text{M}}} $更小，因此 $ {\langle U\rangle _x} $在初始阶段发展要落后于Watanabe等^[27]的数值模拟结果，但从 $ T \approx 25{T_{\text{r}}} $开始，两者 $ {\langle U\rangle _x} $的演化保持高度吻合. 由于采用了周期性边界条件以及系统质量守恒的原因，流向平均速度最终会发展至统计均匀，并且 $ {\langle U\rangle _x}/{U_0} $等于网格的阻塞率0.36. 如图11(b)所示为分别与Melina等^[39]、Nagata等^[40]在风洞实验近场区测的实验数据和Dickey等^[41]、Liu^[42]的拖曳网格湍流实验数据进行对比. 其中利用 $ X/{U_\infty } $将风洞实验的流向距离转化为时间， $ {U_\infty } $为风洞入口速度. 由于和文献工况的阻塞率及雷诺数不同，湍流强度在初始阶段与文献数据存在轻微误差，但同样在 $ T \approx 25{T_{\text{r}}} $之后保持高度吻合. 因此根据规则网格湍流统计量分析结果与参考文献结论吻合良好，进一步证明了本研究数值算法的精度可靠性. 同时，也保证了后文在分析分形网格湍流的湍流强度、能量耗散以及能谱演化等特性时，所采用计算方法的合理性与可靠性.

如图12所示为2类网格Y-Z面中心线上U_c随时间演化结果. 相比于RGT，FGT流向平均速度受到初始条件的影响更大，持续的时间更长. 直到 $ T \approx 300{T_{\text{r}}} $时，U_c才开始趋向统计均匀.

由于分形网格和规则网格Y-Z平面上下左右对称分布，因此截取2类网格的1/4的Y-Z平面进行观察，具体示意图如图13所示. 在2类网格平面上选取相对应的位置进行统计特性对比分析，其中下标fg和rg分别表示分形网格和规则网格，点A_fg和点A_rg可以看作2类网格平面的虚拟中心.

如图14所示为FGT和RGT不同(Y, Z)位置上流向平均速度 $ {\langle U\rangle _x} $和均方根速度 $ {u_{{\text{rms}}}} $随时间演化结果，其中(Y, Z)点对应图13. 在图14(a)中，FGT各位置的流向平均速度在初始阶段表现出强烈的不均匀性，相比于图14(c)，初始条件对FGT的影响持续到 $ T \approx 300{T_{\text{r}}} $才逐渐消失. 其中 $ {E_{{\text{fg}}}} $点位于二阶分形网格上，由于初始阶段主要受到周围三阶网格条产生的小尺度尾流影响，其 $ {\langle U\rangle _x} $迅速降低，但是一阶与二阶网格条产生的尾流大尺度运动随着时间推移而加剧，使得 $ {E_{{\text{fg}}}} $点的 $ {\langle U\rangle _x} $逐渐增大并发展至统计均匀. 其他统计点因自身所处网格位置不同，在初始阶段受到不同大小尾流尺度的影响，会产生不规则波动，但所有位置的 $ {\langle U\rangle _x} $最终趋向统计均匀. 可见，分形正方形网格的较小分形迭代部分的影响仅在初始阶段范围内持续，大部分衰减阶段的湍流特性主要由较大分形迭代部分决定. 因此可以确定时间模拟的FGT在初始阶段产生的多尺度尾流相互作用引起的附加湍流，有利于初始阶段不均匀性的快速恢复. 在Zhou等^[43]模拟的四阶FGT空间演化结果中，也提出了与上述类似的观点. 如图14(b)、(d)所示为FGT与RGT的流向湍流强度. 可以看出，RGT所有位置的特性几乎呈现相同变化趋势，并在 $ T \approx {\text{6}}{T_{\text{r}}} $开始显著衰减. 而FGT由于受到不同大小的多尺度运动影响， $ {u_{{\text{rms}}}} $表现出剧烈的波动行为，其湍流强度峰值时刻并不统一，但A_fg点的中心线上峰值时刻最晚出现在 $ T \approx 50{T_{\text{r}}} $附近. 在大部分衰减阶段，FGT比RGT产生了更高的湍流强度，Laizet等^[17]关于FGT与RGT空间演化的对比研究结果中也存在该现象.

如图15(a)、(b)所示分别为网格中心线上湍流强度 $ {u_{{\text{rms}},{\text{c }}}} $、 $ {v_{{\text{rms}},{\text{c}}}} $和 $ {w_{{\text{rms}},{\text{c}}}} $与泰勒雷诺数 $ R{e_{\lambda ,{\text{c}}}} = {u_{{\text{rms}},{\text{c}}}}{\lambda _{\text{c}}}/\nu $随时间演化结果，其中 $ {\lambda _{\text{c}}} = {u_{{\text{rms}},{\text{c}}}}/\langle {(\partial u'/\partial x)^2}\rangle _x^{1/2} $为中心线上的泰勒尺度. 在图15(a)中，RGT中心线上 $ {u_{{\text{rms}},{\text{c }}}} $、 $ {v_{{\text{rms,c}}}} $和 $ {w_{{\text{rms}},{\text{c}}}} $表现趋势类似，同样在 $ T \approx {\text{6}}{T_{\text{r}}} $时达到峰值. 由于FGT中心线上湍流强度主要受最大正方形网格条产生的大尺度尾流作用影响，其湍流强度峰值出现的时刻( $ T \approx 50{T_{\text{r}}} $)比RGT晚，这一点从图7和图8中时间在 $ T{\text{ = }}50{T_{\text{r}}} $时，两者的流向平均速度可视化结果对比也能得到证明. 前期Zhou等^[44]对RGT空间演化进行数值模拟研究，也验证了中心线上湍流强度峰值的位置 $ {X_{{\text{peak}}}} $很大程度上取决于网格条产生的尾流尺度大小. 尽管本研究的分形网格阶数( $ {N_{\text{f}}} = 3 $)及网格横向厚度比 $ {D_{\text{r}}}{\text{ = }}9.5 $的取值与Hurst等^[45]的空间演化实验不同，但中心线上湍流强度大小刚好介于其实验中 $ {D_{\text{r}}}{\text{ = }}8.5 $和 ${D_{\text{r}}}{\text{ = 13.0}}$的结果之间. 同时图15(a)与Laizet等^[18]的空间演化数值模拟结果相似的是，在本研究的时间演化初始阶段，RGT产生的尾流尺度要大于FGT产生的小尺度，导致其湍流强度大于FGT，但随着时间演化，FGT的大尺度运动加剧，此后( $ T > 22{T_{\text{r}}} $)大部分衰减阶段FGT的湍流强度大小至少为RGT的2倍. 另外，在 $ T = 300{T_{\text{r}}} $之后，根据数据拟合结果显示，RGT的湍流动能衰减规律满足Saffman湍流的特征 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {(T/{T_{\text{r}}})^{ - 6/5}} $，其中 $ {{\boldsymbol{u}}_{{\text{rms}}}} $为均方根速度矢量. 在图15(b)中，RGT中心线上 $ R{e_{\lambda ,{\text{c}}}} $大约在 $ T \approx {\text{6}}{T_{\text{r}}} $达到最大值30.6，该位置与前文分析的湍流强度最大时刻相同. 分形网格不同的几何结构导致了多尺度运动，FGT中心线上的 $ R{e_{\lambda ,{\text{c}}}} $在整个时间演化过程中都呈现出不规则波动行为，并在 $ T \approx 50{T_{\text{r}}} $附近达到最大值43，此后FGT的 $ R{e_{\lambda ,{\text{c}}}} $开始衰减，其在衰减阶段的大小大约为RGT的2倍.

为了观察网格中心线上流向脉动速度偏导数的分布形态，利用其流向脉动速度偏导数的偏斜度 $ {S_{\partial u/\partial x}}{\text{ = }}{\langle {(\partial u'/\partial x)^3}\rangle _x}/\langle {(\partial u'/\partial x)^2}\rangle _x^{3/2} $和平坦度 $ {F_{\partial u/\partial x}}{\text{ = }}{\langle {(\partial u'/\partial x)^4}\rangle _x}/\langle {(\partial u'/\partial x)^2}\rangle _x^2 $2类统计量进行分析. 2类统计量随时间演化结果如图16所示. 可以看出，分形网格的多尺度特性导致了FGT的 $ {S_{\partial u/\partial x}} $和 $ {F_{\partial u/\partial x}} $随时间演化产生较大波动，其中 $ {S_{\partial u/\partial x}} $基本趋势与RGT的保持一致，随着时间增加偏斜度趋于稳定 $ {S_{\partial u/\partial x}} \approx - 0.5 $. 并且Antonia等^[46]和Burattini等^[47]关于均匀各向同性湍流研究中也发现，当 $ R{e_\lambda } < 60 $时， $ - {S_{\partial u/\partial x}} \sim 0.5 $. 此外，尽管FGT的 $ {F_{\partial u/\partial x}} $波动极为剧烈，但其数值大小随时间推移始终分布在2~6，而RGT的 $ {F_{\partial u/\partial x}} $后期稳定在约2.5. Kitamura等^[32]之前在网格湍流空间演化实验中，发现当 $ R{e_\lambda } < 100 $时，平坦度范围满足 $ 2 \leqslant {F_{\partial u/\partial x}} \leqslant 5 $，这与本研究低雷诺数下计算结果高度相似. 最终，两者的流向脉动速度偏导数的偏斜度 $ {S_{\partial u/\partial x}} $为负值，表示其流向脉动速度偏导数分布形态与正态分布相比为负偏. 而平坦度数值与3相比，偏离越大表示流向脉动速度偏导数分布形态的陡缓程度与正态分布的差异程度越大.

流场湍流动能(turbulent kinetic energy)的表达式为 $ K(t){\text{ = }}\langle u_{{\text{rms}}}^2+v_{{\text{rms}}}^2+w_{{\text{rms}}}^2\rangle /2 $，其中符号“ $ \langle \cdot \rangle $”表示对X、Y和Z这3个方向进行空间平均. 如图17所示为网格中心线上 $ {K_{\text{c}}}(t) $、 $ {\varepsilon _{\text{c}}}(t) $、 $ {\lambda _{\text{c}}}(t) $和 $ {\eta _{\text{c}}}(t) $随时间演化结果. 其中 $ {\varepsilon _{\text{c}}}{\text{ = }}2\nu {\langle {S_{ij}}{S_{ij}}\rangle _x} $表示中心线上的脉动耗散率， $ {\eta _{\text{c}}} = {\varepsilon _{\text{c}}}^{ - 1/4}{\nu ^{3/4}} $为Kolmogorov长度尺度. 在图17(a)、(b)中，FGT和RGT网格中心线上湍流动能始终大于能量耗散，2类网格湍流 $ {K_{\text{c}}}(t) $与 $ {\varepsilon _{\text{c}}}(t) $的峰值时刻与图14(b)、(d)湍流强度峰值时刻一致. 在图17(c)、(d)中须注意的是，根据Kolmogorovs长度尺度的定义，在 $ T < 22{T_{\text{r}}} $(见图15(a))时，RGT中心线上湍流强度大于FGT的，尾流相互作用使得RGT的耗散率 $ {\varepsilon _{\text{c}}}(t) $更大，导致其尺度 $ {\eta _{\text{c}}}(t) $小于FGT的. 在 $ T > 22{T_{\text{r}}} $之后，FGT的大尺度运动起主要作用，结合图17(a)、(b)可知其耗散率 $ {\varepsilon _{\text{c}}}(t) $远大于RGT的，因此在衰减阶段FGT的 $ {\eta _{\text{c}}}(t) $大小表现为近乎RGT的1/2. 此外，由图17(c)、(d)可以看出，2类湍流中心线上 $ {\lambda _{\text{c}}} $发展趋势基本相似，并且随着时间的演化，两者 $ {\lambda _{\text{c}}} $的水平大小也越发接近. 在能量级串过程中，由于FGT存在多尺度运动的影响，其泰勒微尺度变化产生不规则波动，但能量会逐渐转移到小尺度上进行耗散. Zhou等^[48]在对双方柱绕流的能量逐尺度转移的研究中发现，在流场下游区域方柱产生的尾流近乎完全湍流，变得更加各向同性和均匀. 因此可以预知的是，本研究FGT和RGT经过较长的时间演化之后也会产生类似的效果.

为了分析波数空间湍动能的变化特征，对2类网格湍流在波数空间的三维能谱 $ E(k) $随时间演化规律展开讨论. 能谱函数的计算式如下：

(21) $ E(k,t) = \int {\int\limits_{ - \infty }^\infty {\int {\frac{1}{2}} } } {{\boldsymbol{\varPhi}} _{ij}}({\boldsymbol{k}},t)\delta (\left| {\boldsymbol{k}} \right| - k){\rm{d}}{\boldsymbol{k}} .$

式中：k为波数空间的标量波数； $ {{\boldsymbol{\varPhi}} _{ij}}({\boldsymbol{k}},t) $为速度谱张量，由2点速度相关函数进行傅立叶变换而来，其中速度采用脉动速度 $ {\boldsymbol{u}} - {\langle {\boldsymbol{U}}\rangle _x} $. 将能谱函数在波数空间进行积分可以得到湍流动能：

(22) $ K(t) = \int_0^\infty {E(k,t)} {\rm{d}}k .$

具体结果如图18所示. 在图18(a)中，当 $ T \leqslant 50{T_{\text{r}}} $时，FGT的能量主要集中在小尺度结构上，根据图5(b)~(f)也可以发现此阶段分形网格较小，分形迭代部分产生的小尺度尾流占主导作用，随着时间推移，结合图14(b)可知湍动能逐渐增大. 当 $ T > 50{T_{\text{r}}} $后，FGT基本上进入衰减阶段，此时大尺度运动开始起到强化湍流的作用，能量集中在低波数范围内，并随时间演化逐渐减小. 在图18(b)中，湍动能在 $ T \approx 6{T_{\text{r}}} $时刻达到最大值，这与图14(d)湍流强度最大时刻保持一致. 同理，RGT的湍动能随时间推移也呈现出先增大后减小的趋势，并且在衰减后期能量也主要集中在低波数范围内. 值得注意的是，由图15(a)可知，在时间 $ T < 22{T_{\text{r}}} $的初始阶段，RGT产生的湍流强度要大于FGT的，这一点由图18(a)和(b)中 $ T/{T_{\text{r}}}{\text{ = }}2 $、6的能谱最大值对比也能得到证明. 这是因为在时间演化初期，FGT中第3阶尾流小尺度要比RGT的尾流尺度更小，从而产生的湍动能更低.

众所周知，泰勒耗散定律^[49]$ \varepsilon (t) \propto {{\boldsymbol{u}}_{{\text{rms}}}}{{\text{(}}t{\text{)}}^3}/{L_{\text{u}}}(t) $和Kolmogorov惯性子区内的标度律是众多湍流统计理论的基础，具有重要的研究意义. 它要求即使在低波数 $ kM \approx L_{\text{u}}^{ - 1}( \approx {k_{\text{f}}}) $下，能量局部平衡假说也是成立的，其中 $ {L_{\text{u}}} $为积分尺度， $ {k_{\text{f}}} = 1 $为截断波数. 但Goto等^[50]进行了高雷诺数下外力强迫周期性盒子湍流的DNS仿真，结果表明该假说在 $ kM \geqslant 20{k_{\text{f}}} $的范围内才成立. 波数在 $ [k,\infty ] $内的单位时间湍流动能 $ {K > }(k,t) $以及其耗散率 $ {\varepsilon > }(k,t) $表达式如下：

(23) $ {K > }(k,t) = \int_k^\infty {E(k',t){\rm{d}}k'}, $

(24) $ {\varepsilon > }(k,t) = 2\nu {\int_k^\infty {k'} ^2}E(k',t){\rm{d}}k'. $

如图19(a)、(b)所示分别为FGT和RGT在不同截断波数下的湍流动能和耗散率随时间演化的结果. 可以看出，2类网格湍流几乎都是当截断波数 $ kM \approx 20{k_{\text{f}}} $时，满足能量与耗散水平相等，但并不等同于能量级串^[51]过程中的能量传输与耗散相平衡的特点. 当波数 $ kM > 20{k_{\text{f}}} $时，随着截断波数的增大，流场中的含能尺度减少，小尺度结构能量耗散加剧，FGT和RGT均表现出耗散大于能量. 在时间演化的初始阶段，FGT和RGT整体的动能在 $ T \approx 6{T_{\text{r}}} $附近达到最大值，这一点也体现在图14(b)和(d)的湍流强度分析中. 结合图15(a)和图18，在时间演化的初期，RGT在尾流运动产生的能量和耗散都大于FGT的，说明可以利用分形结构实现对流动的被动控制，进而有针对性的减缓能量的衰减. 在进入衰减期后，FGT产生的能量及耗散率比RGT更大，足以表明FGT大尺度运动逐渐加剧，起到了强化湍流的作用.

据3.3节分析可知，在时间演化流场的初始阶段，湍动能主要集中在高波数范围内. 因此为了探究初始条件对时间演化流场的具体影响，本研究利用不同的无量纲化方式放大三维能谱在高波数下的主要特征，分析初始阶段小尺度的运动规律. 具体结果如图20、21所示，其中计算三维能谱采用脉动速度 $ {\boldsymbol{u}} - \langle {\boldsymbol{U}}\rangle $，与图18能谱计算所采用的 $ {\boldsymbol{u}} - {\langle {\boldsymbol{U}}\rangle _x} $有所区别. 由图20可以看出，FGT的三维能谱 $ E(k) $峰值在 $ T/{T_{\text{r}}} = {\text{2}} ,4,\cdots ,$50处较为明显，能谱峰值对应的波数与三阶分形网格几何尺寸相关，也包括正方形网格的对角线尺寸. 其中 ${k_{{\text{M}}1}} = 2{\text{π}} /(4M)$、 ${k_{{\text{M}}2}} = 2{\text{π}} /(2M)$和 ${k_{{\text{M}}3}} = 2{\text{π}} /M$分别对应于最大正方形网格尺寸M₀、第2阶网格尺寸M₁和最小网格尺寸M₂. 而且可以看出FGT三维能谱出现的峰值直到 $ T \approx 300{T_{\text{r}}} $才逐渐消失. 同理由图21可以看出，RGT的三维能谱 $ E(k) $同样在初始阶段出现了峰值，对应的波数为规则正方形网格长度尺寸 ${k_{\text{M}}} = 2{\text{π}} /M$，并且峰值在 $ T \approx 25{T_{\text{r}}} $开始消失. 能谱峰值的出现是因为受到初始条件和网格尾流多尺度相互作用的影响，随着时间的演化，能谱峰值逐渐从高波数转移到低波数下，并且2类网格湍流尾流混合更加均匀，初始条件的干扰逐渐减弱，能谱峰值便会随之消失. FGT能谱峰值消失时刻要远远晚于RGT，也可以说明初始条件对FGT的影响持续时间更久，结合图7、8的可视化结果，再次印证了图14(a)、(c)关于2类流场达到统计均匀状态的时间分析. 此外，对比图20(a)、21(a)的结果，纵坐标量纲为速度的平方，其与流场的湍流强度相关. 可以发现，在 $ T < 10{T_{\text{r}}} $内，RGT的纵坐标取值范围为 $ 0 \sim 0.3 $，而FGT的纵坐标取值范围为 $ 0 \sim 0.06 $，可见初始阶段RGT在小尺度范围内产生的尾流强度要远大于FGT的. 这一点与图15(a)和图18(a)、(b)的分析结论一致.

长期以来，学界对于网格湍流中的大尺度结构到底是属于Saffman类型还是Batchelor类型存在很大争议，尚未有明确的结论. 这2种湍流的衰减规律完全不同，其中Saffman湍流^[52]的湍动能衰减规律遵循 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {t^{ - 6/5}} $，Batchelor湍流^[53]则遵循 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {t^{ - 10/7}} $. 结合图15(a)的分析可知，在时间达到 $ T = 300{T_{\text{r}}} $之后，RGT湍动能衰减规律基本满足 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {t^{ - 6/5}} $. 因此对比分析2类湍流场在低波数下的Saffman能谱^[54]${E_{\text{S}}}({{k}}){\text{ = }} L{k^2}/(4{{\text{π}} ^2})$与三维能谱 $ E(k) $，旨在进一步观察时间演化的FGT和RGT中的大尺度结构是否属于Saffman类型. Saffman积分定义如下：

(25) $ L = \int {\langle {\boldsymbol{u}}'({\boldsymbol{x}}){\boldsymbol{u}}'({\boldsymbol{x}}+{\boldsymbol{r}})\rangle } {\rm{d}}{\boldsymbol{r}} .$

式中： $ {\boldsymbol{u}}'({\boldsymbol{x}}) $为脉动速度矢量， $ {\boldsymbol{r}} $为方向矢量.

如图22所示为Saffman积分L的时间演化结果，积分范围为 $0 \leqslant \left| {\boldsymbol{r}} \right| \leqslant {r_0}$， $ {r_0} $的定义为 $\langle {\boldsymbol{u}}'({\boldsymbol{x}}){\boldsymbol{u}}'({\boldsymbol{x}}+ {r_0}{{\boldsymbol{e}}_x})\rangle = 0$. 随着时间的推移，FGT和RGT的Saffman积分L持续减小，尽管FGT减小的更加缓慢，但两者趋势基本相同. 北村拓也等^[55]的规则网格湍流空间演化实验证明了在远离网格的下游区域( $ X/M > 100 $)，L会保持不变. 在图22中，RGT的L在 $ T \approx 300{T_{\text{r}}} $时，整体开始趋于恒定，而此时FGT的L依然在大幅减小，说明RGT中的大尺度结构要比FGT更早显现出Saffman类型特征.

进一步分析FGT和RGT在 $ T = 300{T_{\text{r}}} $时，流场三维能谱 $ E(k) $与Saffman能谱在波数空间的分布情况，结果如图23所示. 可以看出，在能量级串过程中，该时刻FGT能量主要集中在高波数范围内，FGT的小尺度结构含有的湍动能大于RGT的，由图14(b)、(d)的衰减后期湍流强度对比也可以看出这一点. 值得注意的是，在低波数范围内，RGT三维能谱更加符合 $ E(k) \approx {E_{\text{S}}}(k) $的特征，且波数范围为 $1.0 \leqslant kM \leqslant 2.5$. 在该时刻下，FGT的大尺度结构未表现出明显的Saffman特征，可能是因为本研究模拟的FGT的演化时间还不够长，流场须进一步演化至完全各向均匀状态后，再考虑其大尺度结构的特征. 因此可以认为当演化时间 $ T > 300{T_r} $时，规则网格湍流符合Saffman湍流，湍动能衰减规律满足 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {t^{ - 6/5}} $，并且在低波数下( $ 1.0 \leqslant kM \leqslant 2.5 $)的大尺度结构符合Saffman类型.

（1）时间演化的FGT受到初始条件的影响比RGT更持久，直到 $ T \approx 300{T_{\text{r}}} $才逐渐消失. 这导致FGT的流向平均速度之类的统计量在 $ T \approx 300{T_{\text{r}}} $才开始趋近统计均匀，而RGT的各统计量则在 $ T \approx 25{T_{\text{r}}} $已达到统计均匀状态.

（2）在FGT的时间演化过程中，诸如泰勒雷诺数、湍流强度、偏斜度、平坦度等统计特性出现较为剧烈且不规则的波动现象，这是由于分形结构不同尺度的网格条产生的多尺度尾流在相互作用，从而形成了多尺度运动. 与此同时，在相同阻塞率 $ \sigma $和雷诺数 $ R{e_{\text{M}}} $的条件下，FGT产生了比RGT更高的湍流强度 $ {u_{{\text{rms}}}} $.

（3）初始条件对流场的影响主要集中在高波数范围内，且与网格几何尺寸相关，并随着时间演化而逐渐消失. 其次，在网格湍流时间演化的初期，RGT的动能和耗散率均大于FGT的，而在衰减期，由于FGT大尺度运动加剧，FGT的动能及耗散更大. 另外，当截断波数 $ kM \approx 20{k_{\text{f}}} $时，2类流场高波数下的能量与耗散水平保持相等，但有别于能量级串过程中的能量传输与耗散的平衡. 在 $ T > 300{T_{\text{r}}} $后，RGT的大尺度结构符合Saffman类型，其湍动能衰减规律基本满足 $ \langle {\boldsymbol{u}}_{{\text{rms}}}^2\rangle \sim {t^{ - 6/5}} $.

（4）开展时间演化的FGT，设置了周期性边界条件，进而在谱空间中很容易实现探索湍流动能的变化以及流场的能量传输机理，能为后续开展FGT的相关研究奠定理论基础. 同时根据FGT的能量衰减变化特征，发现利用分形结构能够实现流场的被动控制. 但本研究没有考虑分形网格几何形状的变化产生的影响，如不同分形结构或阻塞率，这是今后值得继续探索的研究课题.