61 不可压缩$k-\epsilon$湍流模型

61.1 文献中的$k-\epsilon$湍流模型

Wilcox[67]文献中给出单向流的$k-\epsilon$的控制方程

涡粘方程: $$ \mu_T=\rho C_{\mu}\frac{k^2}{\epsilon} $$ 湍动能: $$ \rho\frac{\partial k}{\partial t}+\rho U_j\frac{\partial k}{\partial x_j} = \tau_{ij}\frac{\partial U_i}{\partial x_j} - \rho\epsilon+\frac{\partial}{\partial x_j}\left[(\mu+\frac{\mu_T}{\sigma_k})\frac{\partial k}{\partial x_j}\right] $$ 湍耗散： $$ \rho\frac{\partial\epsilon}{\partial t} +\rho U_j\frac{\partial\epsilon}{\partial x_j} = C_{\epsilon 1}\frac{\epsilon}{k}\tau_{ij}\frac{\partial U_i}{\partial x_j}-C_{\epsilon 2}\rho\frac{\epsilon^2}{k}+\frac{\partial}{\partial x_j}\left[(\mu+\frac{\mu_T}{\sigma_k})\frac{\partial\epsilon}{\partial x_j}\right] $$ 封闭系数: $$ C_{\epsilon 1} = 1.44, \quad C_{\epsilon 2} = 1.92, \quad C_{\mu} = 0.09, \quad \sigma_k=1.0, \quad \sigma_{\epsilon} = 1.3 $$ 通过如下基本结构形式对$k$和$\epsilon$的输运方程重新组织: $$ \textrm{local derivative + convection + diffusion = source & sink terms} $$ 湍动能: $$ \rho\frac{\partial k}{\partial t}+\rho U_j\frac{\partial k}{\partial x_j}-\frac{\partial}{\partial x_j}\left[\underbrace{(\mu+\frac{\mu_T}{\sigma_k})}{D_k}\frac{\partial k}{\partial x_j}\right] =\underbrace{\tau{ij}\frac{\partial U_i}{\partial x_j}}{G} - \rho\epsilon $$ 湍耗散: $$ \rho\frac{\partial\epsilon}{\partial t} +\rho U_j\frac{\partial\epsilon}{\partial x_j}-\frac{\partial}{\partial x_j}\left[\underbrace{(\mu+\frac{\mu_T}{\sigma{\epsilon}})}{D{\epsilon}}\frac{\partial\epsilon}{\partial x_j}\right] = C_{\epsilon 1}\frac{\epsilon}{k}\underbrace{\tau_{ij}\frac{\partial U_i}{\partial x_j}}{G}-C{\epsilon 2}\rho\frac{\epsilon^2}{k} $$ 扩散系数: $$ D_k = \mu + \frac{\mu_T}{\sigma_k} $$

$$ D_{\epsilon} = \mu + \frac{u_T}{\sigma_{\epsilon}} $$

扩散项常数表达由扩散常数$D_k$和$D_{\epsilon}$合并组成。湍动能方程的右边第一项为湍动能生成$G$

61.2 OpenFOAM中$k-\epsilon$湍流模型

61.2.1 控制方程

OpenFOAM中的$k-\epsilon$模型基本上与61.1节的控制方程相同。因为OpenFOAM语法以向量为主，所以本节以向量描述。OpenFOAM中某型有一定的修正。

首先，$k$和 $\epsilon$的输运方程除以了密度$\rho$。因此，包含粘性的所有项使用运动粘度$\nu$而不是动力粘度$\mu$

第二，OpenFOAM中的标准$k-\epsilon$模型取消了模型常数$\sigma_k$。使该常数取1，则该常数取消了。这种操作不会改变模型。然而，如果用户尝试改变模型常数，不会有任何效果。可参考32.4.2节的讨论和例子。

最后，对流项根据微分规则分为两项。参考公式222。

涡粘，参考Listing 486 $$ \mu_T = \rho\nu_T $$

$$ \nu_T = C_{\mu}\frac{k^2}{\epsilon} $$

湍动能，参考Listing 487 $$ U_j\frac{\partial k}{\partial x_j} = \mathbf{U}\cdot\frac{\partial k}{\partial \mathbf{x}} = \mathbf{U}\cdot\nabla k $$

$$ \mathbf{U}\cdot\frac{\partial k}{\partial \mathbf{x}} = \nabla\cdot(\mathbf{U}k)-(\nabla\cdot\mathbf{U})k $$

$$ \frac{\partial k}{\partial t} + \nabla\cdot(\mathbf{U}k) - (\nabla\cdot\mathbf{U})k - \nabla\cdot(D_k\nabla k) = G - \epsilon $$

湍耗散： $$ \frac{\partial \epsilon}{\partial t} + \nabla\cdot(\mathbf{U}\epsilon) - (\nabla\cdot\mathbf{U})\epsilon - \nabla\cdot(D_{\epsilon}\nabla\epsilon) = C_1G\frac{1}{k} - C_2\frac{\epsilon^2}{k} $$

扩散常数$\sigma_k$通过下面方程可以消除: $$ D_k = \texttt{DkEFF} = \nu + \nu_T $$

$$ D_{\epsilon} = \texttt{DepsilonEff} = \nu + \frac{\nu_T}{\sigma_\epsilon} $$

封闭系数，默认值: $$ C_{\epsilon 1} = 1.44, \quad C_{\epsilon 2} = 1.92, \quad C_{\mu} = 0.09, \quad \sigma_{\epsilon} = 1.3 $$ 湍流模型类的构造函数中设置了模型常数的默认值。

61.2.2 源代码

Listing 486给出了涡粘的计算代码。简短的代码可能会导致混淆，例如函数$sqr()$可能与平方根函数$sqrt()$的功能类似。

Listing 487给出了涡粘的输运方程。方程右边最后一项进行了扩展。 $$ \epsilon = \underbrace{\frac{\epsilon}{k}{k}}_{\texttt{fvm::Sp(epsilon/k, k)}} $$

nut_ = Cmu_ * sqr(K_)/epsilon_

Listing 486: 涡粘的计算

tmp<fvScalarMatrix> kEqn
(
  fvm::ddt(k_)
  + fvm::div(phi_, k_)
  - fvm::Sp(fvc::div(phi_), k_)
  - fvm::laplacian(DkEff(), k_)
==
  G
  - fvm::Sp(epsilon_/k, k_)
);

Listing 487: 湍动能输运方程

构造函数

Listing488给出了$\texttt{kEpsilon}$类的构造函数第一行代码。该构造函数接收5个参数。分号后，为初始化列表。初始化列表包含模型常数的默认值。更多关于C++的细节可以参考56.5节。模型常数$C_{\mu}$在18行定义

kEpsilon::kEpsilon
(
  const volVectorField& U,
  const surfaceScalarField& phi,
  transportModel& transport,
  const word& turbulenceModelName,
  const word& modelName
)
:
  RASModel(modelName, U, phi, transport, turbulenceModelName),
  
  Cmu_
  (
    dimensioned<scalar>::lookupOrAddToDict
    (
      "Cmu",
      coeffDict_,
      0.09
    )
  )
  /* code continue */

**Listing 488: kEpsilon类的构造函数

61.3 $bubbleFoam$和$twoPhaseEulerFoam$中的$k-\epsilon$模型

$bubbleFoam$和$twoPhaseEulerFoam$两个求解器中已经植入了$k-\epsilon$模型，所以不能像OpenFOAM中其他求解器一样使用广义湍流模型。

关于色散两相流中的湍流模型仍然没有完全解答。但是有以下几种策略可供选择。

Per phase，每一项流体使用各自独立的湍流模型

Mixture，使用基于混合量的湍流模型。

Liquid Phase，使用基于液相数量的湍流模型。对于色散项可以忽略或者使用常数。

61.3.1 控制方程

$bubbleFoam$和$twoPhaseEulerFoam$中的$k-\epsilon$湍流模型在某些方面与OpenFOAM的标准$k-\epsilon$模型有区别

1. 使用有效涡粘计算扩散常数。比较方程219，220和234，235
2. 模型常数$\sigma_k$和$\sigma_{\epsilon}$用它们各自的对偶值代替。
3. 与标准$k-\epsilon$模型不同，模型常数$\sigma_k$没有舍去。通过定义常数$\alpha_{1,k} = 1/\sigma_{k}$，来赋值$\sigma_{k}$

$bubbleFoam$和$twoPhaseEulerFoam$使用的湍流模型基于液相数量。气项湍流模型使用模型常数$C_t$。该常数连接液相湍流粘度和气项。如果该常数设为0，则忽略气相湍流模型。

涡粘 $$ \nu_{2,T} = C_{\mu}\frac{k^2}{\epsilon} $$

$$ \nu_{2,eff} = \nu_2 + \nu_{2,T} $$

$$ \nu_{1,eff} = \nu_{1} + C_t^2\nu_{2,T} $$

湍动能，参考Listing487 $$ \frac{\partial k}{\partial t} + \nabla\cdot(\mathbf{U}2k)-(\nabla\cdot\mathbf{U}2)k-\nabla\cdot(\alpha{1,k}\nu{2,eff}\nabla k) = G - \epsilon $$

$$ \frac{\partial \epsilon}{\partial t} + \nabla\cdot(\mathbf{U}2\epsilon)-(\nabla\cdot\mathbf{U}2)\epsilon-\nabla\cdot(\alpha{1,\epsilon}\nu{2,eff}\nabla \epsilon) = C_1G\frac{1}{k} - C_2\frac{\epsilon^2}{k} $$

扩散系数，注意不同定义 $$ \alpha_{1,k} = \frac{1}{\sigma_k} $$

$$ \alpha_{1,\epsilon} = \frac{1}{\sigma_{\epsilon}} $$

$$ D_k=\alpha_{1,k}\nu_{2,eff}=\frac{\nu_{2,eff}}{\sigma_k} $$

$$ D_{\epsilon} = \alpha_{1,\epsilon}\nu_{2,eff} = \frac{\nu_{2,eff}}{\sigma_{\epsilon}} $$

61.3.2 源代码

$bubbleFoam$和$twoPhaseEulerFoam$的输运方程其他内容可参考$\texttt{kEpsilon.H}$文件。Listing 489给出了$\texttt{kEpsilon.H}$文件中的重要代码。

tmp<volTensorField> tgradU2 = fvc::grad(U2);
vosScalarField G(2*nut2*tgradU2() && dev(symm(tgradU2()))));

// Dissipation equation
fvScalarMatrix epsEqn
(
	fvm::ddt(epsilon)
  + fvm::div(phi2, epsilon)
  - fvm::Sp(fvc::div(phi2), epsilon)
  - fvm::laplacian
    (
    alpha1Eps*nuEff2, epsilon,
    "laplacian(DepsilonEff, epsilon)"
    )
  ==
    C1*G*epsilon/k
  - fvm::Sp(G2*epsilon/k, epsilon)
);

// Turbulence kinetic energy equation
fvScalarMatrix eps
(
	fvm::ddt(k)
  + fvm::div(phi2, k)
  - fvm::Sp(fvc::div(phi2), k)
  - fvm::laplacian
    (
    alpha1k*nuEff2, k,
    "laplacian(DkEff, k)"
    )
  ==
    G
  - fvm::Sp(epsilon/k, k)
);

//- Re-calculate turbulence viscosity
nut2 = Cmu*sqr(k)/epsilon;

Listing 489: $bubbleFoam$和$twoPhaseEulerFoam$求解器的湍流输运方程

61.4 湍动能生成模型

在比较文献和源码中的湍流模型方程时，二者在湍动能生成定义上有着巨大的差别。

61.4.1 文献和源码的定义

湍动能生成的定义不同。

H.Rusche[54]的论文中，$bubbleFoam$和$twoPhaseEulerFoam$基于 $$ P_b = 2\nu_{2,eff}(\nabla\mathbf{U}_b\cdot \textrm{dev}(\nabla\mathbf{U}_b+(\nabla{U}_b)^T)) $$ 源码-bubbleFoam对应的文件$\texttt{kEpsilon.H}$中，参考Listing 489的第二行 $$ G = 2\nu_T(\nabla\mathbf{U}_2:\textrm{dev(sym}(\nabla\mathbf{U}_2))) $$ 源码-标准$k-\epsilon$模型对应的文件$\texttt{kEpsilon.C}$中， $$ G=2\nu_T[\textrm{sym}(\nabla\mathbf{U})]^2 $$ Ferzinger Peric[35] $$ P = \mu_T\nabla\mathbf{U}:(\nabla\mathbf{U}+(\nabla\mathbf{U})^T) $$ Wilcox[67] $$ G=\mu_T\nabla\mathbf{U}:(\nabla\mathbf{U}+(\nabla\mathbf{U})^T) - \frac{2}{3}\rho k\mathbf{I}:\nabla\mathbf{U} $$ 有些定义使用动力粘度，有些使用运动粘度。对于不可压流动，这不是主要的区别。

61.4.2 粘度使用的差别

方程237是唯一使用有效粘度代替湍流粘度的。具体原因此处不表。

然而，在Fluent理论指南[7]中所描述那样，当对于$k-\epsilon$模型在高雷诺数的情况下，使用有效粘度计算生成项。但是指南也没有进一步指明什么是$k-\epsilon$模型的高雷诺数情况。

61.4.3 符号

在61.4.1节定义使用向量形式描述相关内容。但是对于方程237，这样会有一点问题。 $$ P_b = 2\nu_{2,eff}(\nabla\mathbf{U}_b\cdot \textrm{dev}(\nabla\mathbf{U}_b+(\nabla{U}_b)^T)) $$

其中，点积符号并不是向量内积。从量纲角度看，除非点积符号表示缩并。因此，该公式应该为： $$ P_b = 2\nu_{2,eff}(\nabla\mathbf{U}_b : \textrm{dev}(\nabla\mathbf{U}_b+(\nabla{U}_b)^T)) $$

61.4.4 文献中的定义

公式240和241的唯一差别为最后一项 $$ G=\mu_T\nabla\mathbf{U}:(\nabla\mathbf{U}+(\nabla\mathbf{U})^T) - \frac{2}{3}\rho k\mathbf{I}:\nabla\mathbf{U} $$ 使用下面等式，缩并可以用内积代替

$$ \mathbf{I}:\nabla\mathbf{U} = \textrm{tr}(\nabla\mathbf{U}) = \nabla\cdot\mathbf{U} $$ 对于不可压流体，通过连续性方程可以推出速度的散度为0 $$ \nabla \cdot\mathbf{U} = 0 $$

$$ G=\mu_T\nabla\mathbf{U}:(\nabla\mathbf{U}+(\nabla\mathbf{U})^T) - \underbrace{\frac{2}{3}\rho k\mathbf{I}:\nabla\mathbf{U}}_{=0} $$ 所以对于不可压缩流动而言，方程240和241是相等的。那么现在我们可以方程240作为参考检验生成项定义的不同。

61.4.5 Rusche和bubbleFoam的定义

bubbleFoam和twoPhaseEulerFoam求解器基于H.Rusche的论文。但是，它们生成项的定义不同。比较方程237和238 $$ P_b = 2\nu_{2,eff}(\nabla\mathbf{U}_b\cdot \textrm{dev}(\nabla\mathbf{U}_b+(\nabla{U}_b)^T)) $$

$$ G = 2\nu_T(\nabla\mathbf{U}_2:\textrm{dev(sym}(\nabla\mathbf{U}_2))) $$

我们忽视连续项速度符号的不同 $$ \mathbf{U}_2 = \mathbf{U}_b $$ 上述两个公式第二个缩并符号不同。可以问一句，下述公式是否成立 $$ \nabla\mathbf{U}_2:\textrm{dev(sym}(\nabla\mathbf{U}_2)) = \nabla\mathbf{U}_b\cdot \textrm{dev}(\nabla\mathbf{U}_b+(\nabla{U}_b)^T) $$ 如果有下述等式，上面这个问题很好回答 $$ \textrm{dev}(\mathbf{T}) = \mathbf{T} - \frac{1}{3}\textrm{tr}(\mathbf{T}) $$

$$ \textrm{sym}(\mathbf{T}) = \frac{1}{2}(\mathbf{T} + (\mathbf{T})^T) $$

$$ \mathrm{dev(sym(\nabla\mathbf{U}_2))} = \mathrm{dev\bigg(\frac{1}{2}(\nabla\mathbf{U}_2 + (\nabla\mathbf{U}_2)^T)}\bigg) $$

$$ \mathrm{dev(sym(\nabla\mathbf{U}_2))} = \frac{1}{2}\mathrm{dev(\nabla\mathbf{U}_2 + (\nabla\mathbf{U}_2)^T)} $$

$$ \mathrm{dev(sym(\nabla\mathbf{U}_2))} = \underbrace{\frac{1}{2}\bigg(\mathrm{(\nabla\mathbf{U}_2 + (\nabla\mathbf{U}_2)^T)} - \frac{1}{3}\mathrm{tr}\mathrm{(\nabla\mathbf{U}_2 + \mathrm{(\nabla\mathbf{U}2}})^T)\bigg)}{=\mathrm{dev(\nabla\mathbf{U}_2 + (\nabla\mathbf{U}_2)^T)}} $$

$$ \mathrm{dev(sym(\nabla\mathbf{U}_2))} = \frac{1}{2}\mathrm{dev(\nabla\mathbf{U}_2 + (\nabla\mathbf{U}_2)^T)} $$

这样就给出了答案 $$ \nabla\mathbf{U}_2:\textrm{dev(sym}(\nabla\mathbf{U}_2)) = \frac{1}{2}\nabla\mathbf{U}_b\cdot \textrm{dev}(\nabla\mathbf{U}_b+(\nabla\mathbf{U}_b)^T) $$ 从源代码的定义中可以发现有两种不同的生成项定义

1. 使用不同的粘性系数，参考公式237和238
2. 公式247和公式254差了2倍的系数。

这个差异的原因不明。H. Rusche参考的文献不在有效。