剪切应力张量divDevReff 与 divDevRhoReff解析

在OpenFOAM中，我们用函数divDevReff 与 divDevRhoReff来计算剪切应力。

1.函数名称含义

下面，解释这两个函数名称的由来：
$$\nabla \cdot \tau =\nabla \cdot {\sigma }^{dev}$$

关于这个公示的推导可查看Cauchy应力张量、剪切力张量与压力的关系。

• div:散度（divergence）；
• dev:每个二阶张量都可分解为 dev 部分（deviatoric）和 hyd 部分（hydrostatic）；
• R:雷诺时均(Reynolds-Average Simulation)；
• eff:有效的(effective)流动方式，包括 层流（laminar）或湍流(turbulent)；
• rho:取决于流体密度是否发生变化。

由此，不难判断divDevReff 用于不可压缩流体，divDevRhoReff用于可压缩流体 。

2.不可压缩流

2.1代码解读

对于turbulence->divDevReff(U)，即操作对象turbulence通过指针调用名为divDevReff(U)的函数。该函数声明如下：
(文件位置：)

tmp<fvVectorMatrix>laminar::divDevReff(volVectorField &U) const 
{
	Return
	(
		-fvm::laplacian(nuEff(), U)
		-fvc::div(nuEff()*dev(T(fvc::grad(U))))
	)
}

注意：计算过程中，采用的是运动粘度。方程一部分用了隐式(fvm)，而另一部分用了显式(fvc)。此外，还调用了一个名为dev()的函数:

//- Return the deviatoric part of a tensor
template<class Cmpt>
inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::oneThirdI*tr(t);
}

（文件位置：$FOAM_SRC/turbulenceModel/incompressible/RAS/laminar/laminar.C）

2.2函数方程

关于函数dev()，可以简单地理解为一个矩阵的dev(deviatoric part)部分：
$${\mathbf{A}}^{dev}=\mathbf{A}-{\mathbf{A}}^{hyd}=\mathbf{A}-\frac{1}{3}tr(\mathbf{A})\mathbf{I}$$

函数dev()的参数为速度场的转置梯度矩阵:
$$\mathbf{A}=[\nabla \mathbf{U}{]}^T$$

我们对代码进行重写：

-fvm::laplacian(nuEff(), U)=$-\nabla \cdot \left({\nu }^{eff}(\nabla \mathbf{U})\right)$

-fvc::div(nuEff()*dev(T(fvc::grad(U))))=$-\nabla \cdot \left({\nu }^{eff}\ast dev((\nabla \mathbf{U}{)}^T)\right)$

$$dev((\nabla \mathbf{U}{)}^T)=(\nabla \mathbf{U}{)}^T-\frac{1}{3}tr((\nabla \mathbf{U}{)}^T)\mathbf{I}$$

将上述代码整合在一起，有：
$$-\nabla \cdot {\mathbf{R}}^{eff}=-\nabla \cdot \left({\nu }^{eff}(\nabla \mathbf{U})\right)-\nabla \cdot \left({\nu }^{eff}\left[(\nabla \mathbf{U}{)}^T-\frac{1}{3}tr((\nabla \mathbf{U}{)}^T)\mathbf{I}\right]\right)$$
整理，可得：
$$-\nabla \cdot {\mathbf{R}}^{eff}=-\nabla \cdot \left({\nu }^{eff}\left[\nabla \mathbf{U}+(\nabla \mathbf{U}{)}^T-\frac{1}{3}tr((\nabla \mathbf{U}{)}^T)\mathbf{I}\right]\right)$$
根据CFD张量公式
$$\mathrm{t}\mathrm{r}\left(\nabla \mathbf{U}\right)\mathbf{I}=\mathrm{t}\mathrm{r}\left((\nabla \mathbf{U}{)}^{\mathrm{T}}\right)\mathbf{I}=\left(\nabla \cdot \mathbf{U}\right)\mathbf{I}$$
对于不可压缩流连续性方程 $\nabla \cdot \mathbf{U}=0$，则：
$$-\nabla \cdot {\mathbf{R}}^{eff}=-\nabla \cdot \left({\nu }^{eff}\left[\nabla \mathbf{U}+(\nabla \mathbf{U}{)}^T\right]\right)=-\nabla \cdot \left(2{\nu }^{eff}\underset{‘形变率\mathbf{D}‘}{\underbrace{\left[\frac{1}{2}\left(\nabla \mathbf{U}+(\nabla \mathbf{U}{)}^T\right)\right]}}\right)=\nabla \cdot \tau$$

对于不可压缩流，$\frac{1}{3}tr((\nabla \mathbf{U}{)}^T)$ 应该恒等于零，是可以忽略的，但考虑到连续性不可能百分百为零，这里我们将其作为误差引入进去。

3.可压缩流

对于可压缩流，调用turbulence->divDevRhoReff(U)函数。该函数包含关键字Rho，这说明在计算可压缩流的剪切应力时，需要考虑由密度变化引起的膨胀粘度项。

3.1代码解读

函数divDevRhoReff(U)

tmp<fvVectorMatrix>laminar::divDevRhoReff(volVectorField &U) const
{
	 Return
	(
		 -fvm::laplacian(muEff(), U) 
		 -fvc::div(muEff()*dev2(T(fvc::grad(U)))) 
	)
}

注意，这里我们应用的是dev2()函数，这种新函数的代码如下：

//- Return the deviatoric part of a tensor
template<class Cmpt>
inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::twoThirdI*tr(t);
}

3.2 函数方程

关于函数dev2()，可以理解为一个矩阵减去两次hyd部分

$${\mathbf{A}}^{dev2}=\mathbf{A}-2{\mathbf{A}}^{hyd}=\mathbf{A}-\frac{2}{3}tr(\mathbf{A})\mathbf{I}$$
经整理（方法同上），其方程为：
$$-\nabla \cdot {\mathbf{R}}^{eff}=-\nabla \cdot \left(2{\nu }^{eff}\underset{‘形变率\mathbf{D}‘}{\underbrace{\left[\frac{1}{2}\left(\nabla \mathbf{U}+(\nabla \mathbf{U}{)}^T\right)\right]}}-\frac{2}{3}{\mu }^{eff}(\nabla \cdot \mathbf{U})\mathbf{I}\right)=\nabla \cdot \tau$$