OpenFOAM 中的壁面函数(一)

本系列来看看 OpenFOAM 中的壁面函数。壁面函数的本质,是边界条件。这里主要来看看壁面函数的基本原理,OpenFOAM 中实现了的壁面函数,以及选择壁面函数的一些参考依据。

1. 壁面函数的基本原理

湍流模拟中,需要对近壁区域进行处理。一般来讲,壁面处理方法包含两类,一类是使用很细的网格,使靠近壁面的第一层网格在粘性层内($y^+ <1$),然后里可以直接解析到粘性层的低雷诺湍流模型;另一类,不直接解析粘性层,而是将第一层网格设置在对数区($y^+> 30$),然后用经验公式来将粘性层和对数区关联起来。下图是一个典型的壁面附近的 $U^+ \text{-} y^+$ 关系图。
壁面律
图片来自 Wikipedia:Law of the wall
在粘性层,满足如下关系
$$
u^+ = y^+
$$
而在对数区,则满足
$$
U^+ = \frac{1}{\kappa}\ln(Ey^+)
$$
其中 $U^+ = U/u_\tau$, $y^+ = yu_\tau/\nu$, $u_\tau = \sqrt{\tau_w/\rho}$,$\kappa\approx 0.41$,$E \approx 9.8$,$y$ 表示与壁面的距离。

本篇以标准壁面函数法来讨论一下壁面函数方法的基本原理,以及壁面函数在 OpenFOAM 中的实现。下面的讨论,先局限在 $k-\varepsilon$ 模型,且第一层网格在对数区的情形。
先来看一下壁面函数方法需要解决什么问题。
有限体积方法中,扩散项的离散可以表示如下:
$$
\nabla \cdot (\nu \nabla U) = \sum_f \left [\nu_f \cdot (\nabla U)_f \right]
$$
当 $f$ 表示的是壁面边界单元时,这时就需要知道在壁面上的速度梯度 $(\nabla U)_f$。壁面上一般对速度 $U$ 采用无滑移条件,如何得到正确的壁面速度梯度,这就是一个问题。这个问题有两个解决思路,一是通过实验或者 DNS 模拟等,得到一条连续的 $U-y$ 曲线,然后从这个曲线求壁面上的导数 $dU/dy$ 来得到壁面上的速度梯度;还有一种思路是,由于最终需要得到的是正确的 $\nu_f \cdot (\nabla U)_f$ ,即壁面上的剪应力,虽然
$$
\tau_w = \nu \cdot \frac{\partial U}{\partial n}\left. \right|_w \neq \nu \frac{U_p-U_w}{y}
$$
其中 $U_p$ 表示第一层网格中心的速度,$U_w$ 表示壁面上的速度。
但是,可以构造一个壁面上的有效粘度 $\nu_{eff}$,以使下式成立
$$
\tau_w = \nu \cdot \frac{\partial U}{\partial n} \left. \right |_w = \nu_{eff} \frac{U_p-U_w}{y} = (\nu + \nu_t ) \cdot \frac{U_p-U_w}{y}
$$

后一种解决方法的好处是,不需要修改动量方程,直接使用 $\frac{U_p-U_w}{y}$ 来代替 $\frac{\partial U}{\partial n} \left. \right |_w $,然后通过设置合适的湍流粘度 $\nu_t$ 的边界条件来修正壁面应力 $\tau_w$。

另一方面,在对数区,$k^+$ 是常数
$$
k^+ = \frac{1}{\sqrt{C_\mu}} \\
$$
其中 $k^+ = k/u_\tau^2$。

$$
k^+ = \frac{1}{\sqrt{C_\mu}} = k/u_\tau^2
$$

$$
u_\tau = C_\mu^{1/4}k^{1/2}
$$
于是
$$
\tau_w = \rho u_\tau^2 = \rho u_\tau \cdot \frac{U}{U^+} = \frac{\rho u_\tau (U_p-U_w)}{\frac{1}{\kappa}\ln(Ey^+)}
$$
若令
$$
\nu_{eff} = \frac{u_\tau y}{\frac{1}{\kappa}\ln(Ey^+)}
$$

$$
\tau_w = \rho \nu_{eff}\cdot \frac{U_p-U_w}{y}
$$
这正是上文提到的第二种解决壁面速度问题的形式。

$$
\nu_{eff} = \frac{ u_\tau y}{\frac{1}{\kappa}\ln(Ey^+)} = \frac{ y^+ \nu}{\frac{1}{\kappa}\ln(Ey^+)} = \nu + \nu_{tw}
$$
于是得到壁面上的湍流粘度为
$$
\nu_{tw} = \nu \cdot \left(\frac{\kappa y^+}{\ln(Ey^+)} -1 \right)
$$

$y^+$ 可以通过不同的方式来得到,具体的计算方法,见后文的 nutWallFunctions 部分。

除了得到壁面上的等效湍流粘度,还需要计算靠近壁面第一层网格的湍动能生成和湍动能耗散项。

湍动能生成项计算如下:
$$
G \approx \tau_w\cdot \frac{\partial (U_p -U_w)}{\partial y}
$$
由速度的壁面律
$$
U^+ = \frac{U_p - U_w}{u_\tau} = \frac{1}{\kappa}\ln(Ey^+) = \frac{1}{\kappa} \ln(E\frac{yu_\tau}{\nu})
$$
注意,$G$ 求的是第一层网格内的值,所以,由
$$
U_p -U_w = \frac{u_\tau}{\kappa} \ln(E\frac{yu_\tau}{\nu})
$$
可以求得第一层网格内的梯度
$$
\frac{\partial (U_p -U_w)}{\partial y} \left. \right|_p = \frac{u_\tau}{\kappa y_p}
$$
于是
$$
G = \tau_w \cdot \frac{u_\tau}{\kappa y_p}
$$
注意,这里的 $\frac{U_p-U_w}{d}$,其实是速度在壁面法向方向的梯度的近似值,这一点见上文 $\nu_t$ 的边界条件部分。

再来看 $\varepsilon$,$\varepsilon$ 的计算基于第一层网格内的湍动生成与湍动能耗散项守恒的假设,即
$$
\rho \varepsilon_p = G = \tau_w \cdot \frac{u_\tau}{\kappa y_p} =\rho\cdot \frac{u_\tau^3}{\kappa y_p}
$$
于是得
$$
\varepsilon_p = \frac{u_\tau^3}{\kappa y_p} = \frac{C_\mu^{3/4}k_p^{3/2}}{\kappa y_p}
$$

至于 $k$,一般认为当第一层网格位于对数区时,不需要在壁面上对 $k$ 加任何限制,用零梯度边界条件即可。

2. 在 OpenFOAM 中的实现

在 OpenFOAM 中,$k$,$\varepsilon$ 和 $\nu_t$ 分别有对应的边界条件可以选择,壁面函数的实现是在这些边界条件里进行的。具体地说, k***WallFunction 用于指定 $k$ 的边界条件, epsilon***WallFunction 用于计算 $\varepsilon$ 和 $G$ 在第一层网格内的值, nut***WallFunction 用来计算 $\nu_t$ 在壁面上的值。 还有就是一个要关心的问题是这些边界条件的调用顺序,这需要通过湍流模型的一段代码来说明,以 kEpsilon 为例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
void kEpsilon::correct()
{
RASModel::correct();
if (!turbulence_)
{
return;
}
volScalarField G(GName(), nut_*2*magSqr(symm(fvc::grad(U_))));

// Update epsilon and G at the wall
epsilon_.boundaryField().updateCoeffs();

// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1_*G*epsilon_/k_
- fvm::Sp(C2_*epsilon_/k_, epsilon_)
);

epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());

solve(epsEqn);
bound(epsilon_, epsilonMin_);

// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_/k_, k_)
);
kEqn().relax();
solve(kEqn);
bound(k_, kMin_);

// Re-calculate viscosity
nut_ = Cmu_*sqr(k_)/epsilon_;
nut_.correctBoundaryConditions();
}

从上述代码,可以将湍流模型的具体计算过程归纳如下:

  1. 计算湍动能生成项 $G$,并修正 $G$ 在第一层网格的值。修正是通过 epsilon_.boundaryField().updateCoeffs(); 来实现的,这里调用 epsilon 的边界条件的 updateCoeffs 函数,实现的操作是修正 $G$ 和 $\varepsilon$ 在第一层网格的值。
  2. 利用更新的 $G$ 构建 epsEqn,然后修改 epsEqnepsEqn().boundaryManipulate(epsilon_.boundaryField()); ),这样做的目的是保证在下一步 solve(epsEqn)的时候,epsilonWallFunction 类型的边界所属的网格的值不会变化,而是保持在 epsilon_.boundaryField().updateCoeffs(); 这一步里设置的值(参考 cfd-online 的这个帖子)。
  3. 求解 epsEqn,得到更新的 $\varepsilon$ 场。
  4. 利用更新的 $\varepsilon$ 场构建并求解 kEqn,得到更新的 $k$ 场。
  5. 计算 $\nu_t$,并更新 $\nu_t$ 在边界上的值(nut_.correctBoundaryConditions()

至于具体的 $k$,$\varepsilon$,以及 $\nu_t$ 的边界条件的实现,见后文。

参考

  1. The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM® and Matlab®
  2. http://www.slideshare.net/fumiyanozaki96/openfoam-36426892