@李东岳 我下面仔细推导了一下：
首先 rAU(1.0/UEqn.A());，这个公式的计算得到的结果是：
对角系数
这是A()函数：
template<class Type>
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
{
tmp<volScalarField> tAphi
(
volScalarField::New
(
"A("+psi_.name()+')',
psi_.mesh(),
dimensions_/psi_.dimensions()/dimVol,
extrapolatedCalculatedFvPatchScalarField::typeName
)
);
tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
tAphi.ref().correctBoundaryConditions();
return tAphi;
}
这是D()对角系数的平均化处理
template<class Type>
Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
{
tmp<scalarField> tdiag(new scalarField(diag()));
addCmptAvBoundaryDiag(tdiag.ref());
return tdiag;
}
边界对对角系数的影响：
template<class Type>
void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
{
forAll(internalCoeffs_, patchi)
{
addToInternalField
(
lduAddr().patchAddr(patchi),
cmptAv(internalCoeffs_[patchi]),
diag
);
}
}
从上面的代码可以得到：
$$
A_p=\frac{\bar{D}}{\Delta V}
$$
其中：
$$
\bar{D}=average(a_p')+diag
$$
式子中average(Ap')表示的是边界对主对角线系数的影响的平均；diag表示的是内部面离散的主对角线系数。
周围系数作为源项：
H()
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
Foam::fvMatrix<Type>::H() const
{
tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi
(
GeometricField<Type, fvPatchField, volMesh>::New
(
"H("+psi_.name()+')',
psi_.mesh(),
dimensions_/dimVol,
extrapolatedCalculatedFvPatchScalarField::typeName
)
);
GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref();
// Loop over field components
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
{
scalarField psiCmpt(psi_.primitiveField().component(cmpt));
scalarField boundaryDiagCmpt(psi_.size(), 0.0);
addBoundaryDiag(boundaryDiagCmpt, cmpt);
boundaryDiagCmpt.negate();
addCmptAvBoundaryDiag(boundaryDiagCmpt);
Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
}
Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
addBoundarySource(Hphi.primitiveFieldRef());
Hphi.primitiveFieldRef() /= psi_.mesh().V();
Hphi.correctBoundaryConditions();
typename Type::labelType validComponents
(
psi_.mesh().template validComponents<Type>()
);
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
{
if (validComponents[cmpt] == -1)
{
Hphi.replace
(
cmpt,
dimensionedScalar(Hphi.dimensions(), 0)
);
}
}
return tHphi;
}
template<class Type>
void Foam::fvMatrix<Type>::addBoundaryDiag
(
scalarField& diag,
const direction solveCmpt
) const
{
forAll(internalCoeffs_, patchi)
{
addToInternalField
(
lduAddr().patchAddr(patchi),
internalCoeffs_[patchi].component(solveCmpt),
diag
);
}
}
这里面的求解过程包含了：
$$
H=\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{\Delta V}
$$
最后：
$$
HbyA= \frac{H}{A}=\frac{\Delta V}{average(a_p')+diag}\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{\Delta V}
$$
$$
HbyA= \frac{H}{A}=\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{average(a_p')+diag}
$$
