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