bug in constrainPressure

[李东岳]

constrainPressure主要是更新压力第二类边界条件，公式如下：

\begin{equation}
\left( \nabla p_{rgh} \right)_f \cdot\bfn_f=
\frac{\left(\mathbf{HbyA}_f^{*}- \frac{1}{{{A^n_{\mathrm{P},f}}}}(\bfg\cdot\bfh\nabla\rho)_f - \mathbf{U}_f \right)\cdot\bfS_f}
{
|\bfS_f|
\frac{1}{{{A^n_{\mathrm{P},f}}}}
}
\end{equation}
上述公式与代码并不一致。在OpenFOAM中，constrainPressure为

forAll(pBf, patchi)
    {
        if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
        {
            refCast<fixedFluxPressureFvPatchScalarField>
            (
                pBf[patchi]
            ).updateCoeffs
            (
                (
                    phiHbyABf[patchi]
                  - rho.boundaryField()[patchi]
                   *MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
                )
               /(magSfBf[patchi]*rhorAUBf[patchi])
            );
        }
    }

多乘了一个密度。应该改为：

forAll(pBf, patchi)
    {
        if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
        {
            refCast<fixedFluxPressureFvPatchScalarField>
            (
                pBf[patchi]
            ).updateCoeffs
            (
                (
                    phiHbyABf[patchi]
                  - MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
                )
               /(magSfBf[patchi]*rhorAUBf[patchi]/rho.boundaryField()[patchi])
            );
        }
    }

在非常老的OpenFOAM版本中，看起来是正确的，与公式一致

setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p_rgh.boundaryField(),
        (
            phiHbyA.boundaryField()
          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
    );