CFD中的壁面函数


如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次

1. 引言

近些年来,工业CFD问题的尺度越来越大,虽然计算机的计算能力一直在提高,但对于湍流求解,壁面网格一直不够致密去直接解析。一些湍流模型要求壁面处的网格位于粘性支层内,这导致计算量大大增加。例如,压缩机叶片的研究需要全场的流动结构。如果叶片壁面网格位于粘性支层,会导致计算网格大量增加。壁面函数针对现存问题,可以通过较少的壁面网格结构获得精准的结果。如果应用壁面函数,壁面的网格点可放置于log区,大大减少网格数量并且毫不损失精度。例如,对于非均一分布(延伸率1.15)的未使用壁面函数的壁面网格,对于基于壁面边界层厚度雷诺数为5000的的壁面流动,大约需要40个边界层网格。若使用壁面函数(第一层网格$y^+=100$),只需要约15个壁面网格。大大节省了计算量。同时,壁面函数还可以减少方程的刚性。

壁面函数的提出依赖于壁面法则。壁面法则用于描述普适性流体在壁面附近的速度与壁面距离的规律。下图表明了壁面流动的无量纲速度$u^+$与无量纲距离$y^+$的典型分布(早期,CFD计算都是2维的,$y$方向表示壁面垂直的方向,因此$y^+$的名字顺势而来)。图中的两个蓝线分别为$u^+$与$y^+$的线性分布与log分布,可以看出在粘性支层,$u^+$与$y^+$符合线性关系: \begin{equation}\label{law1} u^+=y^+ \end{equation} 在log区,$u^+$与$y^+$则符合log分布: \begin{equation}\label{law2} u^+=\frac{\ln (Ey^+)}{\kappa} \end{equation} 其中$E$通常取9.8,$\kappa$通常取0.4。图中的红线,则表示壁面附近的边界层内$u^+$与$y^+$的分布。在壁面附近的粘性支层中,可采用低雷诺数的湍流模型或者高雷诺数的湍流模型附加壁面函数来计算。在使用后者的情况下,粘性支层不应放置网格点,而是把第一个网格点放在log区,这样可以减少计算量。

需要注意的是,不同湍流模型的壁面函数略微不同,但主要处理思想一致。本文主要讨论RANS-kEpsilon模型的壁面函数的计算方法。

如果打不开图像,请右键在新标签页打开图像后刷新几次
2. 方程与算法
在湍流核心区域,剪切应力贡献可以区分为粘性应力与雷诺应力,其为一个对称二阶张量。考虑作用于与$y$方向垂直的面,同时施加于$x$方向的剪切应力分量可以表示为 \begin{equation}\label{tau} \tau_{yx}=\tau_{xy}=(\mu_l+\mu_t)\left(\frac{\p u_2}{\p x}+\frac{\p u_1}{\p y}\right) \end{equation} 其中$u_1,u_2$分别表示$x,y$方向的速度,$\mu_l,\mu_t$分别表示层流粘度与湍流粘度。壁面剪切应力表示作用与壁面处的剪切应力,在不可渗透壁面处,$u_2=0$,因此,壁面剪切应力仅仅$\tau_{xy}$不为0。因此壁面剪切应力可以表示为 \begin{equation}\label{tauw} \tau_w=(\mu_l+\mu_t)\frac{\p u_{1,\rP}}{\p y_\rP}\approx (\mu_l+\mu_t)\frac{u_{1,\rP}}{y_\rP} \end{equation} 其中$u_{1,\rP}$表示壁面网格点处的$x$方向速度分量,$y_\rP$表示壁面网格点处网格距离。若网格点布置在粘性支层,则湍流粘度为0。定义壁面摩擦速度$u_\tau$:
\begin{equation}\label{utau} u_\tau=\frac{u_{1,\rP}}{u^+}=\sqrt{\frac{\tau_w}{\rho}}=\sqrt{(\nu_l+\nu_t)\frac{\p u_{1,\rP}}{\p y_\rP}} \end{equation}
有 \begin{equation} \tau=\rho u_\tau^2=(\mu_l+\mu_t)\frac{u_{1,\rP}}{y_\rP} \end{equation} 即 \begin{equation}\label{x} u_\tau^2=(\nu_l+\nu_t)\frac{u_{1,\rP}}{y_\rP} \end{equation} 并定义$u^+$场与$y^+$场:
\begin{equation}\label{uplus} u^+=\frac{u_{1}}{u_\tau} \end{equation}
\begin{equation}\label{yplus} y^+=\frac{u_\tau y}{\nu_l} \end{equation}
将壁面处的$u^+,y^+$代入到方程\eqref{x}有 \begin{equation}\label{x2} u_\tau^2=(\nu_l+\nu_t)\frac{u^+ u_\tau^2}{y^+ \nu_l} \end{equation} 整理有: \begin{equation}\label{nut} \nu_t=\frac{y^+}{u^+}\nu_l-\nu_l \end{equation} 若网格点布置在log区,有: \begin{equation}\label{nutlog} \nu_t=\frac{y^+ \kappa}{\ln (Ey^+)}\nu_l-\nu_l \end{equation} 大量的实验以及DNS数据表明,在壁面处的log区内,湍流动能的产生与耗散相抵。基于这个假设,有:
\begin{equation}\label{utau2} u_\tau=C_\mu^{0.25}k^{0.5} \end{equation}
\begin{equation}\label{utau3} \frac{\p u_{1,\rP}}{\p y_\rP}=\frac{u_\tau}{\kappa y^+}=\frac{C_\mu^{0.25}k^{0.5}}{\kappa y_\rP} \end{equation}
其中$C_\mu$为无经验常数,通常取0.09。因此有壁面的$y^+_\rP$以及$\tau_w$:
\begin{equation}\label{yplus2} y^+_\rP=\frac{C_\mu^{0.25}k_\rP^{0.5} y_\rP}{\nu_l} \end{equation}
\begin{equation}\label{tau3} \tau_w=(\mu_l+\mu_t)\frac{C_\mu^{0.25}k^{0.5}}{\kappa y_\rP} \end{equation}

在这里需要注意壁面的$y^+_\rP$与$y^+$场计算公式的区别。

对于kEpsilon模型,壁面函数的实施主要在于合理的获得壁面附近合理的湍流动能以及湍流动能耗散率,并修正壁面处的湍流粘度。壁面处的湍流动能可以简单的给定零法向梯度。假如壁面处的网格处于粘性支层,则可以给定壁面的湍流动能为0。壁面处的湍流动能耗散率则通过下式进行计算:

\begin{equation}\label{epsi} \varepsilon_\rP=\nu_l\frac{\p ^2 k_\rP}{\p y_\rP^2}, \mathrm{sublayer} \end{equation}
但这种方式如前文所述,需要大量的网格点位于粘性支层。因此,倾向于将网格点布置在log区,在这种情况下,湍流动能可以给定零法向梯度边界条件,湍流动能耗散率则通过下式进行计算:
\begin{equation}\label{epsi2} \varepsilon_\rP=\frac{C_\mu ^{0.75}k^{1.5}_\rP}{\kappa y_\rP}, \mathrm{logRegion} \end{equation}
同时,对于kEpsilon湍流模型,其湍流产生项可以表示为
\begin{equation}\label{G} G=\nabla\bfU:\tau= \tau_{11}\frac{\p u_1}{\p x}+\tau_{12}\frac{\p u_1}{\p y}+\tau_{13}\frac{\p u_1}{\p z}+ \tau_{21}\frac{\p u_2}{\p x}+\tau_{22}\frac{\p u_2}{\p y}+\tau_{23}\frac{\p u_2}{\p z}+ \tau_{31}\frac{\p u_3}{\p x}+\tau_{32}\frac{\p u_3}{\p y}+\tau_{33}\frac{\p u_3}{\p z} \end{equation}
在网格点位于粘性支层的情况下,壁面处产生项必然为0,即$G_\rP=0$。若网格点位于位于log区,$G_\rP$演化为:
\begin{equation}\label{G2} G_\rP=\frac{\p u_{1,\rP}}{\p y_\rP}\tau_{w} \end{equation}
代入方程\eqref{tau3}有
\begin{equation}\label{G3} G_\rP=(\mu_l+\mu_t)\frac{\p u_{1,\rP}}{\p y_\rP}\frac{C_\mu^{0.25}k_\rP^{0.5}}{\kappa y_\rP} \end{equation}
3. 计算过程

在使用壁面函数的情况下,需要进行下述步骤:首先生成计算网格,依据给定的初始值,通过方程\eqref{yplus2}计算壁面处的$y^+_\rP$。通过方程\eqref{epsi2}与\eqref{G3}计算壁面处的湍流动能耗散率$\varepsilon_\rP$以及产生项$G_\rP$并在本时间步下固定。求解$\varepsilon$与$k$的传输方程。随后通过$\nu_t=C_\mu\frac{k^2}{\varepsilon}$更新全场的湍流粘度。壁面处的$\nu_t$通过壁面函数方程\eqref{nutlog}进行计算。

计算后,后处理计算壁面处的$y^+_\rP$,确认$y^+_\rP>30$。若符合,则表示网格适用于高雷诺数壁面函数,若不符合,糙化壁面处网格重新计算,直至$y^+_\rP>30$为止。

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网