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; } //////////////////////////////////////////////////////////////// |
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 |
Thank you for your help.
I see. |
All times are GMT -4. The time now is 16:50. |