CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Foam::lduMatrix::faceH (http://www.cfd-online.com/Forums/openfoam-programming-development/68487-foam-ldumatrix-faceh.html)

yijin_li September 22, 2009 05:21

Foam::lduMatrix::faceH
 
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 September 23, 2009 03:41

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 September 23, 2009 08:00

Thank you for your help.

I see.


All times are GMT -4. The time now is 06:46.