引 言

条带和流向涡是壁湍流近壁区的典型相干结构. 在与壁面平行的平面内, 低速流体呈现沿流向拉长并在展向准周期变化的条带形状, 其流向尺度约为1000个黏性尺度, 展向平均间距约为80~120个黏性尺度. 条件统计表明, 低速条带位于一对正负流向涡之间, 流向涡则位于低速流体和高速流体之间, 流向涡的流向平均长度约为200个黏性尺度, 平均直径约为20~30个黏性尺度, 并以(0.6~0.8)U(U为槽道截面平均速度)的速度向下游传播. 在流向涡上洗侧, 低速流体被带离壁面, 称为上抛, 在流向涡下洗侧, 高速流体被带向壁面, 称为下扫, 上抛和下扫合称为猝发, 猝发过程是近壁区Reynolds切应力的来源, 在湍流的产生和维持中起重要作用.

Hamilton等^[1]和Jiménez等^[2]分别对平面Couette湍流和平面Poiseuille湍流进行直接数值模拟, 发现了近壁区条带和流向涡的循环再生过程, 该过程被称为自维持过程(self-sustaining process, SSP). 在此过程中, 流向涡向下游运动, 不断将近壁区流体带离壁面, 形成远大于自身长度的低速条带, 从而使流向速度剖面在法向和展向出现拐点, 变得不稳定, 条带失稳进一步产生流向涡^[3].

湍流具有时空多尺度性以及高维性特点, 其中所包含的大尺度有序运动(相干结构)在剪切湍流的发生和发展中起重要作用, 人们通常利用模态分解的方法来提取湍流场中的相干结构, 利用提取的特征模态可以对流场进行低阶近似. POD方法是较为常用的模态分解法, 该方法最早由Lumley等^[4]提出, 并将其用于获得一组信号的最优基函数. 目前POD方法已被广泛用于湍流场数据的处理和分析, 例如Bernard对流问题^[5]圆管流动^[6]方腔流动^[7]等. Aubry等^[8]将POD方法和Navier-Stokes方程相结合, 得到了描述湍流边界层运动的降阶模型方程, 并且使用动力系统理论解释了湍流边界层中的间歇性现象. Podvin^{[9, 10]}分别对最小槽道和全尺寸槽道湍流的流向平均场进行了POD分析, 获得了POD特征模态, 并构造出4维和10维的动力系统方程. 利用POD方法构造的湍流降阶模型可用于湍流模拟和湍流控制等问题.

湍流的大涡模拟近年来发展迅速, 通过数值求解大尺度运动而对小尺度运动进行模化, 使得该方法比直接数值模拟计算量小而比Reynolds平均数值模拟普适性好. 但对于壁湍流, 其近壁区的模拟成为大涡模拟进一步工程应用的瓶颈问题. 在壁湍流大涡模拟中, 如果要完全分辨近壁区流动, 所需的网格量与直接数值模拟不相上下, 例如对于Re=10⁶的湍流边界层, 90%的网格需要布置在占边界层厚度10%的近壁区^[11]. 因此人们尝试发展了多种模型来模化近壁区流动, 如应力模型两层模型以及RANS-LES混合模拟等. 由于POD模态可以反映近壁区湍流的动力学过程, 人们也尝试用其构造近壁模型. Podvin等^[12]提出利用POD模型构造合成边界作为离面边界条件来代替壁面无滑移条件, 并用DNS进行了检验. 在构造合成边界时最重要的是确定每个POD模态的时间系数, 为此Podvin等^[12]提出两种方法, 一种是在近壁处求解降阶模型方程, 获得时间系数, 但此方法仅使用流向和展向每个波数上一个法向模态, 故获得的湍流统计量误差较大, 会影响合成边界以上的计算结果, 并且Podvin等^[12]在流向和展向各使用了10个波数, 近壁降阶模型方程的维数达到220, 计算量较大; 另一种是利用合成边界以上流场的反馈信息确定模态系数, 但仍然使用较多的模态数, 随着槽道尺寸增加, 所需的模态数也会增大, 不便于计算.

本文通过对最小槽道的直接数值模拟, 获得了近壁区湍流自维持过程的流场数据, 利用POD方法获得了近壁自维持过程的降阶模型, 考察了模型方程的动力学特性, 并利用POD模态构造了槽道湍流直接数值模拟的离面边界条件. 与文献[12]不同的是, 本文得到的是对数区最小槽道的离面边界条件, 并在流向和展向做拼接扩展得到全尺寸槽道的离面边界条件进行计算.

1 POD方法简介

本文的工作是基于槽道湍流的POD分析. 首先对两平行平板间由压力梯度驱动的槽道湍流进行了直接数值模拟, 上下平板分别位于y=0和y=2h处, 流向法向和展向坐标分别为x, y, z(或x₁, x₂, x₃), 相应的速度分量为u, v, w(或u₁, u₂, u₃).

流动的控制方程采用旋度形式的Navier-Stokes方程及不可压缩流体的连续性方程

$ \begin{matrix} \frac{\partial \mathit{\boldsymbol{u}}}{\partial t}=\mathit{\boldsymbol{u}}\times \mathit{\boldsymbol{ }}\!\!\boldsymbol{\omega}\!\!\rm{ }-\nabla \varPi +\frac{1}{\mathit{R}e}{{\nabla }^{2}}\mathit{\boldsymbol{u, }} \\ \nabla \cdot \mathit{\boldsymbol{u}}=0. \\ \end{matrix} $

其中u为速度向量, $ \mathit{\boldsymbol{ }}\!\!\boldsymbol{\omega}\!\!\rm{ }\rm{=}\nabla \times \mathit{\boldsymbol{u}} $为涡量向量, $ \mathit{\boldsymbol{ \boldsymbol{\varPi} }}{\rm{ = }}p{\rm{/}}\rho {\rm{ + }}{\left| \mathit{\boldsymbol{u}} \right|^2}/2 $表示总压, 以上方程采用截面平均速度u_m和半槽宽h为特征速度和特征长度进行了无量纲化, Reynolds数Re=u_mh/ν. 壁面采用不可滑移和不可穿透条件, 在流向和展向采用周期边界条件. 控制方程在流向和展向采用Fourier-Galerkin法^[13]法向采用4阶紧致差分格式进行空间离散, 采用3阶精度的时间分裂法进行时间推进.

根据Navier-Stokes方程的对称变换不变性, 槽道湍流在展向关于z=0以及在法向关于槽道中心平面y=h具有对称性, 即

$ \begin{align} & {{P}_{z}}\centerdot \left[(u, v, w, p)(x, y, z, t) \right]= \\ & \left[(u, v, -w, p)(x, y, -z, t) \right], \\ & {{P}_{y}}\centerdot \left[(u, v, w, p)(x, y, z, t) \right]= \\ & \left[(u, -v, w, p)(x, 2h-y, z, t) \right], \\ & {{P}_{yz}}\centerdot \left[(u, v, w, p)(x, y, z, t) \right]= \\ & \left[(u, -v, -w, p)(x, 2h-y, -z, t) \right]. \\ \end{align} $

(1)

其中, P_z, P_y分 别代表速度场和压力场关于z轴和关于y=h的对称变换算子, $ {{P}_{yz}}={{P}_{y}}\circ {{P}_{z}} $.

POD模态为2点速度相关矩阵的特征向量. 将速度场分解为平均流和脉动两部分, 由于槽道流动在流向和展向统计均匀, 可在流向和展向做Fourier变换.

$ \begin{array}{*{35}{l}} \mathit{\boldsymbol{u}}(x,y,z,t)=\mathit{\boldsymbol{U}}(y,t)+\mathit{\boldsymbol{{u}'}}(x,y,z,t)= \\ \sum\limits_{{{k}_{x}}=-\infty }^{\infty }{\sum\limits_{{{k}_{z}}=-\infty }^{\infty }{{{{\mathit{\boldsymbol{\hat{u}}}}}_{{{k}_{x}}{{k}_{z}}}}}}(y,t){{e}^{\rm{i}(\alpha \mathit{{k}_{x}}\mathit{x}\rm{+}\beta \mathit{{k}_{z}}\mathit{z})}}. \\ \end{array} $

其中, $ \alpha =2\text{ }\!\!\pi\!\!\text{ }/{{L}_{x}}, \beta =2\text{ }\!\!\pi\!\!\text{ }/{{L}_{z}} $, 表示流向和展向的基本波数, L_x, L_z分别代表流向和展向的计算域大小.$ \mathit{\boldsymbol{U}}={{\left[U, V, W \right]}^{\rm{T}}} $代表进行流向和展向平均得到的平均流动, 对于本文研究的槽道湍流，

$ \begin{align} & U(y, t)=\frac{1}{{{L}_{x}}{{L}_{z}}}\int_{0}^{{{L}_{x}}}{\int_{0}^{{{L}_{z}}}{u(x, y, z, t)}}\text{d}x\text{d}z, \\ & V(y, t)=0, W(y, t)=0. \\ \end{align} $

在每个流向和展向波数上定义并求解特征值问题：

$ \int_{{{\Omega }_{y}}}{{{{\hat{R}}}_{ij, {{k}_{x}}{{k}_{z}}}}}(y, y'){{\phi }_{j, {{k}_{x}}{{k}_{z}}}}(y')\mathrm{d}y'={{\lambda }_{{{k}_{x}}{{k}_{z}}}}{{\phi }_{i, {{k}_{x}}{{k}_{z}}}}(y). $

(2)

其中, $ {{\hat{R}}_{ij, {{k}_{x}}{{k}_{z}}}}(y, {y}')=\left\langle {{{\hat{u}}}_{i, {{k}_{x}}{{k}_{z}}}}(y)\hat{u}_{j, {{k}_{x}}{{k}_{z}}}^{*}({y}') \right\rangle $, 为波数(k_x, k_z)上速度的2点相关, 上标“*”代表复共轭, 〈…〉代表系综平均;$ {{\lambda }_{{{k}_{x}}{{k}_{z}}}} $为特征值, 反映了模态的能量大小; $ {{\phi }_{i, {{k}_{x}}{{k}_{z}}}} $为特征模态.

求解以上特征值问题可得到每个波数上的POD模态, 利用POD模态可将脉动速度展开为

$ \begin{array}{l} \mathit{\boldsymbol{u'}}(x,y,z,t) = \\ \sum\limits_{{k_x} = - \infty }^\infty {\sum\limits_{{k_z} = - \infty }^\infty {\sum\limits_{n = 1}^\infty {a_{{k_x}{k_z}}^n(t)} } } \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\rm{ }}_{{k_x}{k_z}}^n(y){e^{{\rm{i}}(\alpha {k_x}x + \beta {k_z}z)}}. \end{array} $

其中, $ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{k_x}{k_z}}^n = {\left[{\phi _{1, {k_x}{k_z}}^n(y), \phi _{2, {k_x}{k_z}}^n(y), \phi _{3, {k_x}{k_z}}^n(y)} \right]^{\rm{T}}} $相应波数下的POD模态, n代表模态阶数, 并且模态满足法向正交性, 即

$ \int_{{{\Omega }_{y}}}{(\phi _{1, {{k}_{x}}{{k}_{z}}}^{*n}\phi _{1, {{k}_{x}}{{k}_{z}}}^{m}+\phi _{2, {{k}_{x}}{{k}_{z}}}^{*n}\phi _{2, {{k}_{x}}{{k}_{z}}}^{m}+\phi _{3, {{k}_{x}}{{k}_{z}}}^{*n}\phi _{3, {{k}_{x}}{{k}_{z}}}^{m})}\text{d}y={{\delta }_{mn}}. $

Ω_y代表POD模态在壁面法向的求解范围. 由于物理流场的速度为实数, 故

$ {{\hat{u}}_{i, -{{k}_{x}}-{{k}_{z}}}}=\hat{u}_{i, {{k}_{x}}, {{k}_{z}}}^{*}, \phi _{i, -{{k}_{x}}-{{k}_{z}}}^{n}=\phi _{i, {{k}_{x}}{{k}_{z}}}^{*n}. $

对于2点速度相关张量$ {{\hat{R}}_{ij, {{k}_{x}}{{k}_{z}}}}(y, {y}') $的计算, 除了采用时间平均, 为了扩大样本数量, 还利用式(1)表示的流动对称性, 将样本集合扩大3倍.

在法向采用数值积分法, 故式(2)可写为

$ \frac{1}{{NT}}{\mathit{\boldsymbol{F}}_{{k_x}{k_z}}}\mathit{\boldsymbol{F}}_{{k_x}{k_z}}^*\mathit{\boldsymbol{D}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{_{{k_x}{k_z}}} = \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{_{{k_x}{k_z}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{_{{k_x}{k_z}}}. $

(3)

其中, NT代表DNS得到的瞬时流场个数, F_kxkz存储每个瞬时流场的每个法向位置在(k_x, k_z)波数上脉动速度的值, D代表法向积分矩阵, 则有$ \mathit{\boldsymbol{D}} = {\rm{diag}}({d_1},{d_2}, \cdot \cdot \cdot {d_N},{d_1},{d_2}, \cdot \cdot \cdot {d_N},{d_1},{d_2}, \cdot \cdot \cdot {d_N}) $,其中, d₁, d₂, …, d_N代表法向每个点的积分权重.

$ {{\mathit{\boldsymbol{F}}}_{{{k}_{x}}{{k}_{z}}}}=\left[\begin{matrix} \begin{matrix} {{{\hat{u}}}_{1, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{1}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{1, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{1}}) \\ \end{matrix} \\ \begin{matrix} \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \end{matrix} \\ \begin{matrix} {{{\hat{u}}}_{1, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{N}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{1, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{N}}) \\ \end{matrix} \\ \begin{matrix} {{{\hat{u}}}_{2, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{1}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{2, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{1}}) \\ \end{matrix} \\ \begin{matrix} \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \end{matrix} \\ \begin{matrix} {{{\hat{u}}}_{2, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{N}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{2, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{N}}) \\ \end{matrix} \\ \begin{matrix} {{{\hat{u}}}_{3, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{1}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{3, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{1}}) \\ \end{matrix} \\ \begin{matrix} \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ \end{matrix} \\ \begin{matrix} {{{\hat{u}}}_{3, {{k}_{x}}, {{k}_{z}}(1)}}({{y}_{N}}) & \cdot \cdot \cdot & {{{\hat{u}}}_{3, {{k}_{x}}, {{k}_{z}}(NT)}}({{y}_{N}}) \\ \end{matrix} \\ \end{matrix} \right] $

使用matlab计算特征值和特征向量, 并且为了使得到的结果具有正交性, 将式(3)写为

$ \frac{1}{NT}\left( \sqrt{\mathit{\boldsymbol{D}}}{{\mathit{\boldsymbol{F}}}_{{{k}_{x}}{{k}_{z}}}}\mathit{\boldsymbol{F}}_{{{k}_{x}}{{k}_{z}}}^{*}\sqrt{\mathit{\boldsymbol{D}}} \right)\sqrt{\mathit{\boldsymbol{D}}}{{\mathit{\boldsymbol{ }}\!\!\mathit{\pmb{\Phi}}\!\!\rm{ }}_{{{k}_{x}}{{k}_{z}}}}={{\mathit{\boldsymbol{ }}\!\!\Lambda\!\!\rm{ }}_{{{k}_{x}}{{k}_{z}}}}\sqrt{\mathit{\boldsymbol{D}}}{{\mathit{\boldsymbol{ }}\!\!\mathit{\pmb{\Phi}}\!\!\rm{ }}_{{{k}_{x}}{{k}_{z}}}}. $

(4)

则Hermite矩阵$ \left( \sqrt{\mathit{\boldsymbol{D}}}{{\mathit{\boldsymbol{F}}}_{{{k}_{x}}{{k}_{z}}}}\mathit{\boldsymbol{F}}_{{{k}_{x}}{{k}_{z}}}^{*}\sqrt{\mathit{\boldsymbol{D}}} \right) $的特征向量满足正交性

$ {{\left( \sqrt{\mathit{\boldsymbol{D}}}{{\phi }_{{{k}_{x}}{{k}_{z}}}} \right)}^{*}}\left( \sqrt{\mathit{\boldsymbol{D}}}{{\phi }_{{{k}_{x}}{{k}_{z}}}} \right)=\phi _{{{k}_{x}}{{k}_{z}}}^{*}\mathit{\boldsymbol{D}}{{\phi }_{{{k}_{x}}{{k}_{z}}}}=\mathit{\boldsymbol{I}}. $

2 降阶模型

对Navier-Stokes方程和连续方程

$ \frac{\partial {{u}_{i}}}{\partial t}+{{u}_{j}}\frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}=-\frac{\partial p}{\partial {{x}_{i}}}+\frac{1}{Re}\frac{{{\partial }^{2}}{{u}_{i}}}{\partial {{x}_{j}}\partial {{x}_{j}}}, \frac{\partial {{u}_{i}}}{\partial {{x}_{i}}}=0. $ 在流向和展向做空间平均, 并结合槽道流动的特点, 得到平均运动方程

$ \frac{\partial U}{\partial t}+\frac{\partial \left\langle {u}'{v}' \right\rangle }{\partial y}=-\frac{\partial P}{\partial x}+\frac{1}{Re}\frac{{{\partial }^{2}}U}{\partial {{y}^{2}}}. $

(5)

平均速度U相对于脉动速度$ {{{u}'}_{i}} $随时间的变化要慢得多, 即$ \left| \partial U/\partial t \right|\ll \left| \partial {{{{u}'}}_{i}}/\partial t \right| $并且随着流向和展向计算域的增大, $ \partial U/\partial t $会进一步减小, 因此在本文模型中, 忽略$ \partial U/\partial t $项, 则得到相应的脉动运动方程

$ \begin{matrix} \frac{\partial {{{{u}'}}_{i}}}{\partial t}+U\frac{\partial {{{{u}'}}_{i}}}{\partial {{x}_{1}}}+{{{{u}'}}_{2}}\frac{\partial U}{\partial {{x}_{2}}}{{\delta }_{i1}}+\left( {{{{u}'}}_{j}}\frac{\partial {{{{u}'}}_{i}}}{\partial {{x}_{j}}}-\left\langle {{{{u}'}}_{j}}\frac{\partial {{{{u}'}}_{i}}}{\partial {{x}_{j}}} \right\rangle \right)= \\ -\frac{\partial {p}'}{\partial {{x}_{i}}}+\frac{1}{Re}\frac{{{\partial }^{2}}{{{{u}'}}_{i}}}{\partial {{x}_{j}}\partial {{x}_{j}}}, \\ \frac{\partial {{{{u}'}}_{i}}}{\partial {{x}_{i}}}=0. \\ \end{matrix} $

在式(5)中忽略平均速度的时间导数项, 并进行分步积分, 可将平均速度表示为

$ U\left( y, t \right)=u_{\tau }^{2}Re\left( y-\frac{{{y}^{2}}}{2h} \right)+Re\int_{0}^{y}{\left\langle uv \right\rangle \left( {y}', t \right)\rm{d}\mathit{y}'}. $

(6)

其中, $ {{u}_{\tau }} $为壁面摩擦速度.

为获得脉动运动的降阶模型方程, 需要对波数和POD模态的阶数进行截断, 设截断集合为 $ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{{K_x}{K_z}}^N = \left\{ {\left( {{k_x},{k_z},n} \right)} \right.{\rm{ }}|{\rm{ }}\left| {{k_x}} \right| \le {K_x},\left| {{k_z}} \right| \le {K_z},\left. {\left| n \right| \le N} \right\} $, 在该集合中, 利用POD模态 $ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{k_x}{k_z}}^n $的正交性, 可以得到模态时间系数$ a_{{{k}_{x}}{{k}_{z}}}^{n}\left( t \right) $的常微分方程. 将截断运算表达成类似于大涡模拟中的过滤运算^[14], 定义过滤算子

$ \begin{align} & {{G}_{ij}}\left( \Delta x, y, {y}', \Delta z \right)=\frac{1}{{{L}_{x}}{{L}_{z}}}\times \\ & \sum\limits_{n=1}^{N}{\sum\limits_{-{{K}_{x}}}^{{{K}_{x}}}{\sum\limits_{-{{K}_{z}}}^{{{K}_{z}}}{\phi _{i, {{k}_{x}}{{k}_{z}}}^{n}}}}\left( y \right)\phi _{j, {{k}_{x}}{{k}_{z}}}^{*}\left( {{y}'} \right)\exp \left( \text{i}\Delta x\alpha {{k}_{x}}+\text{i}\Delta z\beta {{k}_{z}} \right). \\ \end{align} $

则有

$ u_{i}^{<}\left( \mathit{\boldsymbol{x}} \right)=\int\limits_{\Omega }{{{G}_{ij}}}\left( x-{x}', y, {y}', z-{z}' \right){{{u}'}_{j}}\left( {\mathit{\boldsymbol{{x}'}}} \right)\rm{d}\mathit{\boldsymbol{{x}'}}. $

进一步将脉动速度分解为

$ {{{u}'}_{i}}=u_{i}^{<}+u_{i}^{>}. $

其中, $ u_i^ < $代表由截断集合 $ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{{K_x}{K_z}}^N $中的POD模态表示的可解尺度运动, $ u_i^ > $代表不可解尺度运动, 并假设过滤可以和求导运算交换顺序, 则过滤后的可解尺度脉动速度满足的方程为

$ \begin{align} & \frac{\partial u_{i}^{<}}{\partial t}+U\frac{\partial u_{i}^{<}}{\partial {{x}_{1}}}+u_{2}^{<}\frac{\partial U}{\partial {{x}_{2}}}{{\delta }_{i1}}+\frac{\partial }{\partial {{x}_{j}}}{{\left( u_{i}^{<}u_{j}^{<} \right)}^{<}}= \\ & -\frac{\partial {{p}^{<}}}{\partial {{x}_{i}}}+\frac{1}{Re}\frac{{{\partial }^{2}}u_{i}^{<}}{\partial {{x}_{j}}\partial {{x}_{j}}}-\frac{\partial \tau _{ij}^{<}}{\partial {{x}_{j}}}. \\ \end{align} $

(7)

其中, $ \tau _{ij}^{<}=-{{\left( {{{{u}'}}_{i}}{{{{u}'}}_{j}}-u_{i}^{<}u_{j}^{<} \right)}^{<}} $为不可解模态构成的亚格子应力.

将可解尺度运动用截断集合内的POD模态展开

$ \begin{align} & {\boldsymbol{u}^{<}}\left( x, y, z, t \right)= \\ & \sum\limits_{{{k}_{x}}=-{{K}_{x}}}^{{{K}_{x}}}{\sum\limits_{{{k}_{z}}=-{{K}_{z}}}^{{{K}_{2}}}{\sum\limits_{n=1}^{N}{a_{{{k}_{x}}{{k}_{z}}}^{n}}}}\left( t \right)\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\rm{ }}_{{k_x}{k_z}}^n\left( y \right){e^{{\rm{i}}\left( {\alpha {k_x}x + \beta {k_z}z} \right)}}. \\ \end{align} $

(8)

将式(6)和式(8)代入式(7), 利用正交性可得到模态时间系数$ a_{{{k}_{x}}{{k}_{z}}}^{n}\left( t \right) $的演化方程

$ \begin{array}{*{20}{l}} {\;\;\;\;\;\;\frac{{{\text{d}}a{{_{{k_x}}^n}_{{k_z}}}}}{{{\text{d}}t}} = \sum\limits_{m = 1}^N {l_{{k_x}{k_z}}^{mn}} a_{{k_x}{k_z}}^m + } \\ {\sum\limits_{m = 1}^N {\sum _{\begin{array}{*{20}{c}} {{{k'}_x} = - {K_x}} \\ {{{k'}_z} = - {K_z}} \end{array}}^{{K_x},{K_z}}\sum\limits_{p,q = 1}^N {c_{{k_x}{k_z}{{k'}_x}{{k'}_z}}^{mnpq}} } a_{{{k'}_x}{{k'}_z}}^pa_{{{k'}_x}{{k'}_z}}^{q*}a_{{k_x}{k_z}}^m + } \\ {\sum _{\begin{array}{*{20}{c}} {{{k'}_x} = - {K_x}} \\ {{{k'}_z} = - {K_z}} \end{array}}^{{K_x},{K_z}}\sum\limits_{p,q = 1}^N {q_{{k_x}{k_z},{k_x} - {{k'}_x},{k_z} - {{k'}_z}}^{pqn}} a_{{{k'}_x}{{k'}_z}}^pa_{{k_x} - {{k'}_x},{k_z} - {{k'}_z}}^q + } \\ {\sum\limits_{m = 1}^N {\left( {1 + \frac{{\alpha {v_{\text{T}}}}}{v}} \right)} l_{v,{k_x},{k_z}}^{mn}a_{{k_x},{k_z}}^m.} \end{array} $

(9)

其中, 方程右端第1, 2项代表与平均流有关的作用项, 第3项为非线性项, 第4项为黏性和亚格子应力项, 为使表达式简洁, 令$ \mathit{\boldsymbol{k}} = \left( {{k_x}, {k_z}} \right), \mathit{\boldsymbol{k'}} = \left( {{{k'}_x}, {{k'}_z}} \right) $, 则式(9)中相应算子的表达式如下

$ \begin{array}{*{20}{l}} {l_k^{mn} = - \int_{\mathit{\Omega }{_y}} {\phi _{1,k}^{*n}} u_\tau ^2Re\left( {1 - \frac{y}{h}} \right)\phi _{2,k}^m{\rm{d}}\mathit{y}{\rm{ - }}}\\ {\;\;\;{\rm{i}}{\mathit{k}_\mathit{x}}\alpha \int_{\mathit{\Omega }{_\mathit{y}}} {\phi _{{\rm{i}},\mathit{k}}^{{\rm{*}}\mathit{n}}\mathit{u}_\tau ^{\rm{2}}} \mathit{Re}\left( {\mathit{y}{\rm{ - }}\frac{{{\mathit{y}^{\rm{2}}}}}{{{\rm{2}}\mathit{h}}}} \right)\phi _{{\rm{i}},\mathit{k}}^\mathit{m}{\rm{d}}\mathit{y},}\\ {\;\;\;c_{k,k'}^{mnpq} = - Re\int_{\mathit{\Omega }{_\mathit{y}}} {\phi _{1,k}^{*n}\phi _{1,k'}^p\phi _{2,k'}^{*\mathit{q}}} \phi _{2,k}^m{\rm{d}}\mathit{y}{\rm{ - }}}\\ {\;\mathit{i}{\mathit{k}_\mathit{x}}\alpha \mathit{Re}\int_{\mathit{\Omega }{_\mathit{y}}} {\left( {\phi _{\mathit{i},\mathit{k}}^{{\rm{*}}\mathit{n}}\phi _{\mathit{i},\mathit{k}}^\mathit{m}\int_{\rm{0}}^\mathit{y} {\phi _{{\rm{1}},\mathit{k'}}^\mathit{p}\phi _{{\rm{2}},\mathit{k'}}^{{\rm{*}}\mathit{q}}} {\rm{d}}\mathit{y'}} \right)} {\rm{d}}\mathit{y},}\\ {\;q_{k,k - k'}^{pqn} = - \int_{\mathit{\Omega }{_\mathit{y}}} {\phi _{i,k}^{*n}} \left( {{\rm{i}}\alpha \left( {{\mathit{k}_\mathit{x}} - {{\mathit{k'}}_\mathit{x}}} \right)\phi _{\mathit{i,k - k'}}^\mathit{q}\phi _{{\rm{1}},\mathit{k'}}^\mathit{p}} \right. + }\\ {\left. {\;\;\frac{{\rm{d}}}{{{\rm{d}}\mathit{y}}}\phi _{i,k - k'}^q\phi _{2,k'}^p + {\rm{i}}\beta \left( {{\mathit{k}_\mathit{z}} - {{\mathit{k'}}_{\rm{z}}}} \right)\phi _{\mathit{i},\mathit{k - k'}}^\mathit{q}\phi _{{\rm{3}},\mathit{k'}}^\mathit{p}} \right){\rm{d}}\mathit{y},}\\ {l_{v,k}^{mn} = \frac{1}{{Re}}\int_{\mathit{\Omega }{_y}} {\phi _{i,k}^{*n}\left[ {\left( { - {\alpha ^2}k_x^2 - {\beta ^2}k_z^2} \right)\phi _{i,k}^m + \frac{{{{\rm{d}}^2}}}{{{\rm{d}}{\mathit{y}^{\rm{2}}}}}\phi _{i,k}^m} \right]} {\rm{d}}\mathit{y}.} \end{array} $

方程中与压力梯度有关的项相比于其他项为小量, 本文忽略不计. 黏性项中的ν_T称为涡黏系数, 代表了不可解模态的影响.

3 近壁湍流的降阶模型

为了得到近壁湍流的降阶模型, 首先对最小槽道湍流进行直接数值模拟以获得流场数据. 最小槽道是指能维持湍流的流向和展向最小尺寸的槽道, 最早由Jiménez等^[15]针对缓冲区提出, 其展向尺寸只能包含一根条带, 可以维持缓冲区湍流准周期的猝发运动. 后来人们对最小槽道的概念进行了拓展, Flores等^[16]指出, 要在高度y_h以下获得健康的湍流, 则最小槽道尺寸应满足L_z≈3y_h, L_x≈2L_z. 本文希望最小槽道能较准确地预测$ 0\le y_{\rm{h}}^{+}\le 70 $的近壁区湍流性质, 因此将最小槽道尺寸设为$ \left( L_{x}^{+}, L_{z}^{+} \right)=\left( 420, 210 \right) $, 其中“⁺”代表以黏性尺度作为特征量进行无量纲化以后的流动物理量, 例如$ L_{x}^{+}={{L}_{x}}{{u}_{\tau }}/v $, ν为运动黏性系数. 本文采用定流量计算, 摩擦Reynolds数$ R{{e}_{\tau }}={{u}_{\tau }}h/v\approx 180 $, 3个方向的网格数分别为$ {{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}=96\times 96\times 96 $.

图 1显示了最小槽道平均壁面切应力随时间的变化. 壁面切应力随时间有比较大的波动, 反映了条带失稳生成流向涡的准周期的猝发过程, 从图中可以估算出猝发周期T⁺≈400, 这与Jiménez等^[17]的描述相一致. 图 2显示一个猝发过程中摩擦力增长阶段的流向速度等值面的变化, 可见条带发生弯曲失稳, 此时流向涡产生, 壁面摩擦力增大.

从DNS计算结果中选取$ \Delta {{T}^{+}}\approx 3000 $时间内的2000个瞬时场进行统计平均, 法向积分域 $ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_y} $选为0<y⁺<70, 利用式(4)获得POD模态和所对应的特征值. 为了考核预测模型的精度, 首先将DNS数据在POD模态上投影, 获得准确的模态演化系数

$ \begin{array}{l} \hat a_{{k_x}{k_z}}^n\left( t \right) = \\ \int_0^{{y^ + } = 70} {\left( {\phi _{1, {k_x}{k_z}}^{ * n}} \right.{{\hat u}_{1, {k_x}{k_z}}}\left( {y, t} \right)} + \phi _{2, {k_x}{k_z}}^{ * n}{{\hat u}_{2, {k_x}{k_z}}}\left( {y, t} \right) + \\ \phi _{3, {k_x}{k_z}}^{ * n}{{\hat u}_{3, {k_x}{k_z}}}\left( {y, t} \right){\rm{d}}y. \end{array} $

模态系数满足$ \left\langle {\hat a_{{k_x}{k_z}}^n\left( t \right)\hat a_{{k_x}{k_z}}^m\left( t \right)} \right\rangle = \lambda _{{k_x}{k_z}}^n{\delta _{mn}} $, 其中$ \lambda _{{k_x}{k_z}}^n $为POD模态的特征值.

3.1 模态特征值

图 3显示了不同波 数上第1阶POD模态所占能量的比重, 即$ \lambda _{{k_x}{k_z}}^1/\sum\limits_{n = 1}^\infty {\lambda _{{k_x}{k_z}}^n} $, 由 此可以计算出低波数(|k_x|≤2, |k_z|≤2)上第1阶模态所占能量达到该波数能量70%以上, 在|k_x|≤5, |k_z|≤10的波数范围内(大尺度运动), 第1阶POD模态为主要模态, 此结果与文献[18]一致. 选取降阶模型截断波数为|k_x|≤3, |k_z|≤4, 且每个波数上取第1阶模态, 即n=1, 此范围内模态能量占脉动能量的比例可达60%.

3.2 模态的空间分布形式

能量最大的两个模态对应的流向波数k_x=0, 即沿流向模态形式不变. 图 4显示了这两个模态的在y-z截面上流向速度云图及展向-法向速度矢量, 可见这两个能量最大的模态就是条带和流向涡结构.

图 5显示了k_x=0, k_z=1波数上第1阶POD模态的特征函数分布, 并与Podvin等^[9]的结果进行了比较. 本文结果与他们的结果相符.

3.3 降阶模型计算结果

将DNS数据投影到POD模态上作为初始条件, 利用Runge-Kutta方法推进求解式(9), 可得到POD模态系数随时间的演化. 在此之前, 首先需要利用POD模态准确系数$ \hat a_{{k_x}{k_z}}^n $求出涡黏系数ν_T. 由于本文仅取1阶POD模态, n=1, 为书写简洁, 在下文中均用a_kxkz代替$ a_{{k_x}{k_z}}^1 $. 将需要模化的亚格子应力$ \tau _{ij}^ < $变换到波数空间, 并在相应波数的第1阶POD模态上进行投影可得波数(k_x, k_z)上的第1阶模态系数T_kxkz, 根据文献[10], 可设T_kxkz=$ {v_{\rm{T}}}l_{v, {k_x}{k_z}}^{11}{a_{{k_x}{k_z}}} $, 对每个波数做能量加权平均, 并在时间方向上取平均值, 可得

$ {v_{\rm{T}}} = \frac{{\sum\limits_{{k_x}, {k_z} \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{{K_x}{K_z}}^1} {{\lambda _{{k_x}{k_z}}}{\rm{R}}e\left[{{T_{{k_x}{k_z}}}l_{v, {k_x}{k_z}}^{11}\hat a_{{k_x}{k_z}}^ * } \right]} }}{{{{\sum\limits_{{k_x}, {k_z} \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{{K_x}{K_z}}^1} {{\lambda _{{k_x}{k_z}}}\left| {l_{v, {k_x}{k_z}}^{11}{{\hat a}_{{k_x}{k_z}}}} \right|} }^2}}}. $

(10)

式(10)中, Re[…]代表取实部运算, 利用该式可求得ν_T/ν=0.647.

图 6显示了几个波数上的模态系数a_kxkz的实部随时间的演化. 从图中可以看出, 流向波数k_x=0的模态能量较大且频率较低, 随着波数的升高, 特别是流向波数不为0的模态, 其模态系数频率增高, 并且高频运动呈现一定的间歇性. 流向波数不为0的模态以行波的形式向下游传播, 根据Jiménez等^[20]的理论, 行波$ \phi = \hat \phi \exp \left[{{\rm{i}}{k_x}\left( {x-ct} \right)} \right] $的相速度c满足$ c = {\mathop{\rm Im}\nolimits} \left( {\left\langle {\phi \partial {\phi ^ * }/\partial t} \right\rangle /\left( {{k_x}\left\langle {\phi {\phi ^ * }} \right\rangle } \right)} \right) $, 其中Im(…)代表求虚部, 〈…〉代表时间平均. 由此可得到每个流向波数k_x≠0的模态的相速度, 本文所选取模态的最大相速度为$ c_{{k_x}{k_z}}^ + $=13.2, 最小相速度为$ c_{{k_x}{k_z}}^ + $=12.0, 根据Taylor冻结假设, 近壁区的相干结构迁移速度大约为(10~15)u_τ, 故此结果与Taylor冻结假设近似一致.

最小槽道和POD分解均是对槽道湍流进行的简化, 考察降阶模型能否捕捉到近壁区湍流的主要特点, 需要对比降阶模型的近壁区统计量和全尺寸槽道近壁统计量, 图 7显示了Reynolds应力和平均流的分布, 对比了POD降阶模型方程最小槽道DNS以及全尺寸槽道DNS^[19]的结果. 由图可以看出, 在近壁处, 降阶模型给出的流向法向和展向脉动强度比最小槽道DNS的结果偏小, 而Reynolds切应力则相近, 且在更远离壁面的区域降阶模型与DNS会出现较大差异, 这是由于降阶模型方程在流向和展向只选取很少的波数, 并且由于降阶模型方程求解得到的低波数模态, 特别是能量较大的流向不变模态, 其时 间系数a_kxkz的强度比精确时间系数$ {\hat a_{{k_x}{k_z}}} $偏小, 从而造成统计量强度偏小.

4 降阶模型的应用

本文使用最小槽道的近壁POD模态, 利用投影法可以获得合成边界, 并可以通过周期扩展形成较大尺寸槽道的离面边界条件. 将最小流动单元扩展形成较大计算域的边界是一种构造壁模型的新方法, Pascarelli等^[21]将最小槽道进行了展向扩展, Garcia-Mayoral等^[22]则将转捩边界层中的湍斑作为最小流动单元构造离面边界条件, 由于转捩边界层中的主要流动结构为确定频率的T-S波, 故边界条件的时间发展可以由波动的形式确定, 但转捩边界层的无序程度远低于近壁充分发展湍流, 所以这种方法并非合理.

本文将最小槽道DNS的近壁区湍流场投影到POD模态上, 利用投影获得的模态系数构造离面边界条件. 图 8显示了不同截断集合$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{{K_x}{K_z}}^N $下采用投影法获得的近壁湍流统计量的分布, 并和DNS结果进行了对比. 由图可以看出, 当K_x≥10, K_z≥10, N>5时, 再增加波数和法向模态对于结果的改善影响很小, 因此本文计算采用的截断集合为(K_x, K_z, N)=(10, 10, 5), 离面边界条件施加在y⁺=50处, 此高度已经越过Reynolds切应力和流向脉动强度的峰值.

用作离面边界条件的y⁺=50处的速度可由POD模态及其投影系数合成得到, 即

$ \begin{array}{l} \mathit{\boldsymbol{u}}\left( {x, {y^ + } = 50, z, t} \right) = \\ \sum\limits_{{k_x} = - {K_x}}^{{K_x}} {\sum\limits_{{k_z} = - {K_z}}^{{K_z}} {\sum\limits_{n = 1}^N {\hat a_{{k_x}{k_z}}^n\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{k_x}{k_z}}^n} } } \left( {{y^ + } = 50} \right)\exp \left( {{\rm{i}}\alpha {k_x}x + {\rm{i}}\beta {k_z}z} \right). \end{array} $

(11)

该离面边界条件首先在最小槽道的DNS中进行了考核, 计算的下边界位于距下壁面y⁺=50处, 该处的速度由式(11)给出. 图 9显示了湍流统计量的分布, 除了靠近下边界的过渡区以外, 采用离面边界条件的结果和最小槽道DNS的结果相吻合.

若要利用该合成边界计算全尺寸槽道, 需要在展向和法向同时进行周期性扩展, 扩展时要保证流动结构物理尺度不变. Pascarelli等^[21]采用D-scaling方法, 直接对物理空间内最小流动单元做展向拼接, 相当于谱空间内谱系数变换成离散分布, 假设最小槽道流向和展向分别扩展n_x, n_z倍, 将波数(K_x, K_z)上的谱系数平移到波数(n_xk_x, n_zk_z)上, “D”代表discrete. Mizuno等^[23]使用C-scaling构造离面边界条件, 即对谱空间做波数配置后, 新配置的波数对相邻波数平均分配谱系数, 并保证总能量不变, 该方法变换后波数为连续分布, 即“C”代表continuous. 本文采用D-scaling的方法, Pascarelli等^[21]指出虽然构造的合成边界在流向和展向有周期性, 但边界以上流动周期性会很快消失.

图 10显示了在流向和展向分别均扩展2倍(n_x=2, n_z=2)和5倍(n_x=5, n_z=5)的计算结果, 并和全尺寸槽道DNS结果进行了对比, 其中扩展5倍的计算域与全尺寸槽道^[19]相当. 由图可见平均速度Reynolds切应力和流向脉动强度吻合较好, 但法向和展向脉动会在离面边界处有所偏差. 首先因为边界处正好在法向和展向脉动峰值处, 此处对计算条件改变较为敏感；其次, 由于边界的高度周期性, 而内部流动周期性消失, 在边界附近会存在过渡区域.

5 结论

本文对最小槽道湍流进行了直接数值模拟, 获得了近壁自维持过程的数据, 对DNS结果进行了POD分析, 并构造了POD模态系数所满足的降阶模型方程. 本文还利用POD模态构造近壁模型, 这是对全尺寸槽道流动在近壁面做简化, 这样的简化是两个层次. 第1层次是将近壁面的自维持过程简化成若干最小流动单元, 用最小槽道代表. 第2层次是在第1层次简化的基础上再将最小槽道近壁流动用更少的POD模态表示. 实际的计算表明, POD模态可以很好地描述近壁面的统计性质, 并且构造的离面边界条件可以扩展应用于全尺寸的槽道流动计算中. 然而, 本文的计算仍然局限在低Reynolds数的范围内, 在高Reynolds数流动计算中, 由于外区的大尺度运动会对近壁流动产生影响, 使得近壁面的合成边界条件的构造变得更加复杂; 其次, 本文的POD模态的获得需要利用大量DNS或者实验数据, 如何不依赖于预先获得的数据构造模态需要进一步的研究.