壁面函数matlab,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$ 表示与壁面的距离。1$)，然后里可以直接解析到粘性层的低雷诺湍流模型；另一类，不直接解析粘性层，而是将第一层网格设置在对数区($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

47void 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 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 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();

}

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

计算湍动能生成项 $G$，并修正 $G$ 在第一层网格的值。修正是通过 epsilon_.boundaryField().updateCoeffs(); 来实现的，这里调用 epsilon 的边界条件的 updateCoeffs 函数，实现的操作是修正 $G$ 和 $\varepsilon$ 在第一层网格的值。

利用更新的 $G$ 构建 epsEqn，然后修改 epsEqn( epsEqn().boundaryManipulate(epsilon_.boundaryField()); )，这样做的目的是保证在下一步 solve(epsEqn)的时候，epsilonWallFunction 类型的边界所属的网格的值不会变化，而是保持在 epsilon_.boundaryField().updateCoeffs(); 这一步里设置的值(参考 cfd-online 的这个帖子)。

求解 epsEqn，得到更新的 $\varepsilon$ 场。

利用更新的 $\varepsilon$ 场构建并求解 kEqn，得到更新的 $k$ 场。

计算 $\nu_t$，并更新 $\nu_t$ 在边界上的值(nut_.correctBoundaryConditions())

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

参考

The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM® and Matlab®