Code of the Week #19: fvc::reconstruct

#201806

在 interFoam 的 pEqn.H 里面,有这样一段:

    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
    U.correctBoundaryConditions();

这一段程序是出现在 pEqn.H 里面求解 p 方程,并进行通量非正交修正后,对 U 变量实行的一种修正方式。
那么为什么要这么做?

相信大家对下面这个方程也不再陌生了

\begin{equation} \mathbf U_\mathrm{P}^{n+1} = \mathbf{U}^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla p^{n+1} \end{equation} \tag{5}

这个方程是利用 U 的预测值和 p 的求解值来获得 U 的真实值,在上一期 code of the week 我们提到了通量的修正。但是,当通量修正后,速度并没有自动修正,而要重新处理。

而速度处理有两种方式,一种是直接利用上面这个式子获得,在 incompressible 的一些求解器里面,使用的是这样一种方式,比如 icoFoam:

            U -= rUA*fvc::grad(p);

那开头那一种方式是怎么回事呢?
进行修正了的通量是 phi,未进行修正的通量是 phiU,他们的差值,自然就是

\text{phi} - \text{phiU} = -\frac{1}{A_p} \mathbf S_f \cdot \nabla p \\

式中, \dfrac{1}{A_p} 就是 rAUf。当前, phi 和 phiU 都已求得,那如何从守恒的通量 \text{phi} 获得守恒的 U 呢?

我们来看看 reconstruct 都干了啥? 文件是 fvcReconstruct.C,

template<class Type>
tmp
<
    GeometricField
    <
        typename outerProduct<vector,Type>::type, fvPatchField, volMesh
    >
>
reconstruct
(
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    typedef typename outerProduct<vector, Type>::type GradType;

    const fvMesh& mesh = ssf.mesh();

    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
    (
        new GeometricField<GradType, fvPatchField, volMesh>
        (
            IOobject
            (
                "volIntegrate(" + ssf.name() + ')',
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
          & surfaceSum((mesh.Sf()/mesh.magSf())*ssf),
            zeroGradientFvPatchField<GradType>::typeName
        )
    );
    GeometricField<GradType, fvPatchField, volMesh>& reconField =
        treconField();

    // Note:
    // 1) Reconstruction is only available in cell centres: there is no need
    //    to invert the tensor on the boundary
    // 2) For boundaries, the only reconstructed data is the flux times
    //    normal.  Based on this guess, boundary conditions can adjust
    //    patch values
    // HJ, 12/Aug/2011

    GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
        surfaceSum((mesh.Sf()/mesh.magSf())*ssf);

    reconField.internalField() =
    (
        inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField())
      & fluxTimesNormal.internalField()
    );

    reconField.boundaryField() = fluxTimesNormal.boundaryField();

    treconField().correctBoundaryConditions();

    return treconField;
}

其中,核心的一段是

            inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
          & surfaceSum((mesh.Sf()/mesh.magSf())*ssf),

这一段直观地、返回的量是

\begin{align} {\Delta \boldsymbol X} =& {\left(\sum_f \frac{\boldsymbol S_f\otimes \boldsymbol S_f}{|\boldsymbol S_f|}\right)^{-1} \left(\sum_f \frac{\boldsymbol S_f}{|\boldsymbol S_f|} \underbrace{\boldsymbol A_p(\text{phi}-\text{phiU})}_{\text{(phi-phiU)/rAUf}}\right)}\\ = &{\left(\sum_f \frac{\boldsymbol S_f\otimes \boldsymbol S_f}{|\boldsymbol S_f|}\right)^{-1} \left(\sum_f \frac{\boldsymbol S_f}{|\boldsymbol S_f|} \boldsymbol S_f \nabla p \right)}\\ =& - \nabla p \end{align}

:sweat:为什么可以这样?设一对称张量 \mathcal D

\mathcal D=\sum_f \frac{1}{|\boldsymbol S_f|} \boldsymbol S_f\otimes \boldsymbol S_f

\begin{align} \mathcal D\cdot \boldsymbol U =& \sum_f \frac{\boldsymbol S_f}{|\boldsymbol S_f|}\boldsymbol S_f \boldsymbol U\\ \sum_f \left(\boldsymbol n_f \otimes \boldsymbol S_f \right)\cdot \boldsymbol U =& \sum_f \boldsymbol n_f (\boldsymbol S_f \boldsymbol U_f)\\ \sum_f \boldsymbol n_f \left( \boldsymbol S_f \cdot \boldsymbol U \right) =& \sum_f \boldsymbol n_f (\boldsymbol S_f \boldsymbol U_f)\tag{7} \end{align}

上式的意思可以表述为:对称张量将速度向量 \boldsymbol U 映射到与 \boldsymbol n_f 平行的向量上,其幅值为 |\boldsymbol n_f| \left( \boldsymbol S_f \cdot \boldsymbol U \right) \approx |\text{phi}| ,即 \boldsymbol n_f \text{phi} 。OpenFOAM 认为式 (7) 是成立的,因此,运用了这一方法来重构速度。

有些人对此两种方法进行了对比,认为差别并不大。Henry Weller 也表示过 fvc::reconstruct() 可能会导致较大的误差。openfoamwiki 认为使用这个的原因是消除棋盘压力[2]。因此,暂时未有定论。对于哪一种方法更好,可能是 depends on cases。

参考文献

  1. https://www.cfd-online.com/Forums/openfoam-programming-development/77943-fvc-reconstruct-algorithm.html
  2. http://openfoamwiki.net/index.php/OpenFOAM_guide/Reconstruction

似乎这么说不完全.因为这个reconstruct在布辛涅斯克近似中也使用,而且是在Ueqn.H中编写,那时候尚未需要求解压力泊松方程

我的意思是,这个东西不完全是为了求解压力泊松方程后更新速度用的