• 返回主页
  • CFD中的壁面函数

    2017.05.22:经过和邱小平讨论,更正公式$\eqref{28}$,$\eqref{epsi6}$,$\eqref{epsi8}$,$\eqref{epsi88}$,$\eqref{epsi888}$,$\eqref{wallnu}$中的$\nu_t$为$\nu_{eff}$;全文结构重排,统一公式
    2017.05.21:公式$\eqref{G3}$,$\eqref{G}$,$\eqref{GF}$有误,去掉其中的$C_1\varepsilon/k$项

    概述

    在讨论CFD的壁面函数之前。首先介绍CFD中湍流的不同流型。自由剪切流(Free Shear Flow)的典型流型为尾涡(Wake)。顾名思义,自由剪切流距离壁面较远。湍流的起源为平均速度的差异$^{[1]}$。CFD中的壁面流(Wall-bounded Flow)不同于自由剪切流,壁面流通常以壁面为固定的流动边界。其进而可以分为内流和外流,圆管内的流动即为内流,飞机外空气的流动即为外流。

    考虑图$1$中的槽导流。其高度为$2\delta$。这个几何的特点为足够的长:$L \gg \delta$,且足够的宽:$b \gg \delta$。且$z$方向的平均速度为零。即$W=0$。(注:$U$,$V$,$W$表示即时速度分量。$\overline{U}$,$\overline{V}$,$\overline{W}$表示时均后的速度,$U'$,$V'$,$W'$表示速度波动分量)
    槽导流
    图1. 槽导流示意图

    我们定义充分发展的湍流为:$x$方向的时均速度$\overline{U}$不随着$x$的变化而变化。因此,充分发展的槽导流在统计学上是一维的$^{[1]}$,也即速度只和$y$有关。

    在三维问题下我们有雷诺平均后的连续性方程: \begin{equation} \frac{\partial \overline{U}}{\partial x}+\frac{\partial \overline{V}}{\partial y}+\frac{\partial \overline{W}}{\partial z}=0 \label{cont} \end{equation} 依据充分发展的槽导流的定义:$\overline{W}=0$以及$\frac{\partial \overline{U}}{\partial x}=0$。方程($\ref{cont}$)可以简化为: \begin{equation} \frac{\partial \overline{V}}{\partial y}=0 \label{cont2} \end{equation} 对方程($\ref{cont2}$)做积分后(略),并由于壁面处的$\overline{V}=0$,因此在计算域内所有位置处$\overline{V}=0$。

    流体内部的剪切应力

    剪切应力(Shear Stress)$\tau$分为粘性应力(Laminar Shear Stress)贡献$\tau_l$以及雷诺应力(Reynolds Stress)贡献$\tau_t$。其可以表示为: \begin{equation} \tau=\tau_l+\tau_t=\rho\nu\left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right)+\left(-\rho \overline{U'V'}\right) \label{tau} \end{equation} 其中$\nu$为流体的层流粘度,$\overline{U'}$和$\overline{V'}$的乘积总是负值$^{[3]}$,因此雷诺应力为正值。

    壁面剪切应力

    如果考虑壁面处,由于光滑壁面处速度为零,进而湍流应力项$-\rho \overline{U'V'}$在壁面处为零。因此壁面附近处的剪切应力可定义为纯粘性贡献(如图2所示): \begin{equation} \tau_w=\rho\nu\left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right)_{y\rightarrow 0} \label{tauW} \end{equation}
    壁面剪切力
    图2. 壁面剪切力示意图


    湍流粘度

    湍流粘度$\nu_t$可以定义为: \begin{equation} \nu_t\frac{\rd \overline{U}}{\rd y}=-\rho\overline{U'V'} \label{nut} \end{equation} 将方程$\eqref{nut}$带入到方程$\eqref{tau}$,有流体内部的剪切应力为: \begin{equation} \tau=\rho(\nu+\nu_t)\left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right) \label{tauNew} \end{equation} 接下来我们定义四个在边界层理论中比较重要的流体变量:

  • 摩擦速度:
  • \begin{equation} u_\tau=\sqrt{\frac{\tau_w}{\rho}} \label{fu} \end{equation}
  • 摩擦距离:
  • \begin{equation} y_\tau=\frac{\rho\nu}{\sqrt{\rho \tau_w}} \label{fy} \end{equation}
  • $u^+$:
  • \begin{equation} u^+=\frac{\overline{U}_\rP}{u_\tau} \label{uplus} \end{equation}
  • $y^+$(类似与雷诺数,表征粘性和湍流过程的相对重要性):
  • \begin{equation} y^+=\frac{u_\tau y_\rP}{\nu} \label{yplus} \end{equation} 其中$_\rP$表示进壁面网格单元的值。在下文中我们将讨论边界层分区以及壁面函数的实施步骤。有关边界层分区的定义以及$y^+$的基本介绍请参考:湍流神秘y+之旅

    近壁法则

    1. 在粘性支层中,主导的剪切应力为粘性应力,雷诺应力在粘性支层中可以忽略不计($\nu_t=0$)。粘性支层部分非常薄,通常$y^+<5$,且占整个边界层的10-20%$^{[1]}$。依据近壁法则的假定以及粘性支层的特征有: \begin{equation} \tau_w=\tau_l=\rho\nu \left( \frac{\mathrm{d} \overline{U}}{\mathrm{d} y}\right)_{y\rightarrow 0}\approx\rho\nu \left( \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP}\right) \label{sublayer} \end{equation} 对方程($\ref{sublayer}$)做积分有: \begin{equation} \overline{U}_\rP=\frac{\tau_w}{\rho\nu}y_\rP+C \label{sublayer2} \end{equation} 由于壁面处$y=0$,$\overline{U}=0$,因此$C=0$。也即: \begin{equation} \overline{U}_\rP=\frac{\tau_w y_\rP}{\rho\nu} \label{sublayer3} \end{equation} 方程($\ref{sublayer3}$)表示的即为粘性支层内$u^+=y^+$。

    2. 在log区内,进一步假定:

  • 雷诺应力可以用普朗特混合长模型计算;且普朗特混合长$l_m$可以表示为:
  • \begin{equation} l_m=\kappa y_\rP \end{equation}
  • 壁面处的湍流动能耗散和湍流动能的产生相均衡;

  • 由于log区内雷诺应力占主要部分,有: \begin{equation} \tau=\tau_t=-\rho \overline{U'V'} \approx \tau_w \label{logtau} \end{equation} 进一步结合普朗特混合长理论有: \begin{equation} \tau_w=\rho l_m^2 \left| \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP} \right| \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP}=\rho l_m^2 \left( \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP} \right)^2 \label{zero} \end{equation} 其中调用了假定为在图1所示的槽导流中下半部分,$\frac{\mathrm{d} \overline{U}}{\mathrm{d} y}>0$。结合方程($\ref{logtau}$)有: \begin{equation} \left( \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP} \right)^2=\frac{\tau_w}{\rho l_m^2 }=\left(\frac{u_\tau}{l_m}\right)^2 \label{log2} \end{equation} 即: \begin{equation} \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y_\rP} =\frac{u_\tau}{l_m}=\frac{u_\tau}{\kappa y} \label{log3} \end{equation} 依据$y^+$的定义,有$\mathrm{d}y_\rP=\frac{\mu \mathrm{d}y^+}{u_\tau \rho}$。将其带入到方程($\ref{log3}$)有: \begin{equation} \frac{\mathrm{d} \overline{U}_\rP}{\mathrm{d} y^+} =\frac{u_\tau}{l_m}=\frac{u_\tau}{\kappa y^+} \label{log4} \end{equation} 对方程($\ref{log4}$)做积分有: \begin{equation} u^+=\frac{1}{\kappa}\mathrm{ln}\left( y^+ \right)+C=\frac{1}{\kappa}\mathrm{ln}\left( Ey^+ \right) \label{log5} \end{equation} 其中$\kappa=0.4-0.42$,$C=5.0-5.5$,$E\approx 9.8$。

    壁面函数

    首先考虑不使用壁面函数的情况:当壁面处的第一层网格足够致密的时候(处于粘性支层),可以认为: \begin{equation} \tau_w=\mu\left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right)_{y\rightarrow 0}\approx\mu\frac{\overline{U_\rP}}{y_\rP} \label{wallf} \end{equation} 但是通常情况下为了节省计算资源,第一层网格$P$节点并不会布置在粘性支层。因此$ \left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right)_{y=0} \neq \frac{\overline{U_\rP}}{y_\rP}$。进而$\tau_w \neq \mu \frac{\overline{U_P}}{y_\rP}$。图3形象的表示了这一过程。在图3中,B点表示位于粘性支层的网格点。在这部分区域,公式($\ref{wallf}$)是成立的,并且$\tau_w$的值等于图中黑色实线的斜率。A点表示网格节点位于Log区。这时,黑色虚线的斜率和黑色实线的斜率不等,因此其并不真正的等于壁面剪切力的值。因此需要重新计算方程$\eqref{wallf}$中的$\mu$,使得$ \mu \left( \frac{\mathrm{d}\overline{U}}{\mathrm{d} y}\right)_{y=0} =\mu_{eff} \frac{\overline{U_\rP}}{y_\rP}$。可见壁面函数最重要的部分就是计算这个有效粘度$\mu_{eff} $。

    壁面剪切力
    图3. 不同网格节点计算的$\tau_w$示意图


    在壁面函数中又细分为均衡的壁面函数、非均衡的壁面函数、增强型壁面函数、粗糙的壁面函数等。本文只讨论均衡的壁面函数。在使用壁面函数的过程中,$k$以及$\varepsilon$仍然可以通过求解全场的方法来获得。需要补充的是,当网格布置在log区的时候(也即标准的壁面函数),$k_\rP$的获得和不使用壁面函数的情况下相同。仍可以使用法向梯度为零的边界条件即$\frac{\partial k_\rP}{\partial y_\rP}=0$。如果第一层网格点布置在粘性支层内,则可以使用固定值边界条件$k_\rP=0$。复杂的是边界处$\varepsilon_\rP$的确定,依据$\varepsilon_\rP$的定义: \begin{equation} \varepsilon_\rP=\frac{C_\mu^{0.75} k_\rP^{1.5}}{l_m} \label{epsiD} \end{equation} 可以看出在边界处$k_\rP$和$l_m$同时趋向于零。因此,壁面处的$\varepsilon_\rP$的值需要进一步的模化来给定。下文详述如何获得壁面处的$\varepsilon_\rP$的值,以及如何对湍流粘度进行修正。

    壁面处湍流动能耗散率计算

    在log区的均衡的光滑标准壁面函数中,满足上文中提及的“壁面处的湍流动能耗散和湍流动能的产生相均衡”。在三维情况下即: \begin{equation} -\overline{U_i'U_j'}\frac{\partial \overline{U_i}}{\partial x_j}=\varepsilon_\rP \label{epsi} \end{equation} 考虑本文讨论算例,一维情况下有$\frac{\partial \overline{U}_\rP}{\partial x_\rP}=0$以及$\frac{\partial \overline{V}_\rP}{\partial y_\rP}=0$。因此上式可以化为: \begin{equation} -\overline{U_\rP'}\overline{V_\rP'}\frac{\partial \overline{U}_\rP}{\partial y_\rP}=\varepsilon_\rP \label{epsi2} \end{equation} 将公式($\ref{log3}$)代入到公式($\ref{epsi2}$)有: \begin{equation} -\overline{U_\rP'}\overline{V_\rP'} \frac{u_\tau}{\kappa y_\rP}=\varepsilon_\rP \label{epsi3} \end{equation} 依据公式($\ref{fu}$),有: \begin{equation} -\rho\overline{U_\rP'}\overline{V_\rP'}=\tau_w=\rho u_\tau ^2 \label{epsi4} \end{equation} 代入到公式($\ref{epsi2}$),有: \begin{equation} \frac{u_\tau ^3}{\kappa y_\rP}=\varepsilon_\rP \label{epsi5} \end{equation} 回到壁面应力的定义。在这里假定整个边界层内的剪切应力是常量,且和壁面剪切力相同$^{[2]}$,即 \begin{equation} \tau\approx\tau_w=\rho\nu_{eff}\left( \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d} y_\rP}\right) \label{28} \end{equation} 其中$\nu_{eff}=\nu+\nu_t$。结合方程$\eqref{zero}$,有: \begin{equation} \nu_{eff}=\tau_w/\rho \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP}=l_m^2 \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP} \label{epsi6} \end{equation} 同时,方程$\eqref{log2}$可以表示为: \begin{equation} \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP}=\frac{u_\tau}{l_m} \label{epsi7} \end{equation} 结合方程$\eqref{epsi6}$和方程$\eqref{epsi7}$有: \begin{equation} \nu_{eff}=u_\tau l_m \label{epsi8} \end{equation} 依据kEpsilon模型中湍流有效粘度的定义: \begin{equation} \nu_{eff}=\frac{C_\mu k^2}{\varepsilon} \label{epsi88} \end{equation} 并结合$\eqref{epsiD}$,有 \begin{equation} \nu_{eff}=C_\mu^{0.25} k^{0.5}_\rP l_m \label{epsi888} \end{equation} 结合方程($\ref{epsi8}$),有: \begin{equation} u_\tau=C_\mu^{0.25} k^{0.5}_\rP \label{epsi9} \end{equation} \begin{equation} \varepsilon=\frac{u_\tau^3}{\kappa y_\rP}=\frac{C_\mu ^{0.75}k^{1.5}_\rP}{\kappa y_\rP} \label{epsi10} \end{equation} 也即: \begin{equation} y^+=\frac{y_\rP u_\tau}{\nu}=\frac{y C_\mu^{0.25} k^{0.5}_\rP}{\nu} \label{y+} \end{equation} \begin{equation} u^+=\frac{\overline{U}}{u_\tau}=\frac{\overline{U} u_\tau}{u_\tau^2}=\frac{\overline{U} C_\mu^{0.25} k^{0.5}_\rP}{\tau_w/\rho} \label{u+} \end{equation} 当处于log区的情况下,把方程($\ref{u+}$)和方程($\ref{y+}$)带入到方程($\ref{log5}$),并结合方程($\ref{wallf}$)有: \begin{equation} \nu_{eff}=\frac{y^+}{u^+}\nu \label{wallnu} \end{equation} 且壁面剪切应力可以这样计算: \begin{equation} \tau_w=\frac{\rho C_\mu^{0.25} k^{0.5}_\rP \overline{U_P}}{u^+} \label{walltau} \end{equation}

    壁面函数在CFD中的实施过程

    如上所述,壁面函数在CFD中的实施过程主要分为以下步骤:

  • (1)通过求解离散的$k$方程获取全场的$k$值;
  • (2)通过方程($\ref{epsi10}$)确定壁面节点处的$\varepsilon_P$;
  • (3)通过方程($\ref{wallnu}$)确定壁面节点处的$\nu_P$;
  • (4)通过$\nu_P$计算湍流动能耗散率传输方程中的壁面的湍流动能耗散率产生项;
  • (5)求解$\varepsilon$方程;

  • 在这里需要注意的是,三维问题下$\varepsilon$传输方程在壁面处的产生项这样计算: \begin{equation} G_\rP= (\nu_t+\nu) \frac{\partial \overline{U_i}}{\partial x_j}\left( \frac{\partial \overline{U_i}}{\partial x_j} + \frac{\partial \overline{U_j}}{\partial x_i} \right) \label{G3} \end{equation} 在壁面处可以简化为: \begin{equation} G_\rP=(\nu_t+\nu) \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP} \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP} \label{G} \end{equation} 结合方程($\ref{log3}$),其可以化简为: \begin{equation} G_\rP=(\nu_t+\nu) \frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP} \frac{u_\tau}{\kappa y_\rP}=(\nu_t+\nu)\frac{\mathrm{d}\overline{U}_\rP}{\mathrm{d}y_\rP} \frac{C_\mu^{0.25} k^{0.5}_\rP}{\kappa y_\rP} \label{GF} \end{equation}

    其他要点

  • 在使用壁面函数处理分离流的时候,由于在分离点的壁面应力趋向于$0$,这回导致上述的壁面函数失效;
  • 在壁面不光滑的情况下,除了粘性应力、雷诺应力,壁面还会形成一个附加的力,这使得log区内的常数有不同的取值;
  • 均衡的壁面函数假定壁面处的湍流动能耗散和湍流动能的产生相均衡,然而在某些情况下并非如此;
  • 某些情况下,如某些高度三维的流场,壁面处不可以简化为一维模型;
  • 上文讨论的壁面函数模型通过一个网格节点即可实施,没有考虑此节点和壁面直接内部的湍流变量的变化,针对这一特性,可以衍生出其他多层壁面函数模型;


  • 参考文献
    [1]. Pope S B. Turbulent flows[J]. 2001.
    [2]. Versteeg H.K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Pearson Education, 2007.
    [3]. 欧特尔H, 朱自强, 钱翼稷. 普朗特流体力学基础[J]. 2008.
    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国