[yijin_li]

Based on the pressure solution, assemble conservatie face flux F
F=s_f*H(u)-a_p_n(p_n-p_p)
Correct face flux phi -= pEqn.flux();
I found that in the code
faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];

why use this expression to correct face flux?

////////////////////////////////////////////////////////////////
Foam::lduMatrix::faceH
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::lduMatrix::faceH(const Field<Type>& psi) const
{
tmp<Field<Type> > tfaceHpsi(new Field<Type> (lduAddr().lowerAddr().size()));
Field<Type>& faceHpsi = tfaceHpsi();

if (lowerPtr_ || upperPtr_)
{
const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();

// Take refereces to addressing
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();

for (register label face=0; face<l.size(); face++)
{
faceHpsi[face] = Upper[face]*psi[u[face]]
- Lower[face]*psi[l[face]];
}
}
else
{
// No off-diagonal. Bug fix for conjugate matrices. HJ, 27/Nov/2008
faceHpsi = pTraits<Type>::zero;
}

return tfaceHpsi;
}
////////////////////////////////////////////////////////////////

[hjasak]

Because of consistency: since the pressure equation IS the continuity equation, fuxes must be built in exactly the same way as the pressure matrix coefficients and sources.

Work out the maths by hand and you will see that what you need is the off-diagonal pressure coefficient (+ explicit corrections, eg. non-orthogonality, which are treated separately).

Hope this helps,

Hrv

[yijin_li]

Thank you for your help.

I see.