Some questions about Laplacian with nonorthogonal correction
 

Hi,
I have taken a look at the nonorthogonal correction for the Laplacian for a scalar variable in Openfoam and compared the treatment of the corrections to the description of the diffusive term in Hrv's thesis.

I will first write what i think is happening in the OF code and then have some questions, because not everthing is clear to me. The code lines in the following refer to the current version of the source guide.

As far as i figured out the interesting things happen in gaussLaplacianSchemes.C after line 39.
In line 56 with fvmLaplacianUncorrected(gammaMagSf, vf) the implicit part of the Laplacian is written into the matrix. I think this is the function defined in gaussLaplacianScheme.C in line 46. It writes into the of-diagonal matrix coeficient
gamma_f*|s|/|d|
with nomenclature according to Hrv's thesis. Then the diagonal matrix entries are filled with the negative sum of the of-diagonal enties (why the negative sum?). The boundaries are also taken care of.
Then we jump back to gaussLaplacianSchemes.C. In line 66 the content of fvm.faceFluxCorrectionPtr() is filled with:
gamma_f*|s|*correction
correction is defined in correctedSnGrads.C in line 42 and just returns the values of fullGradCorrection(vf), which is defined in correctedSnGrad.C in line 44.
fullGradCorrection(vf) a bit confused me, but i think it just returns the inner product of the interpolated gradient at the face with the nonOrthCorrectionVector:
nonOrthCorrectionVector dot nabla(vf)_f

The nonOrthCorrectionVectors and the geometric things needed for them are built in surfaceinterpolation.C as follows:
deltaCoeffs = 1/|d| (line 232)
nonOrthDeltaCoeffs= 1/(n dot d), (stabilised with the max...) (line 292)
nonOrthCorrectionVectors= n - d/(n dot d) (line 340),
which is just
nonOrthCorrectionVectors= k/|s|= (s-(ss/sd)d) /|s|
with k according to Hrv's thesis k=s-(ss/sd)d

So gammaMagSf*this->tsnGradScheme_().correction(vf) in gaussLaplacianSchemes.C line 66 fills the faceFluxCorrectionPtr with:
gamma_f*|s|*correction = gamma_f*|s|*k/|s|*nabla(vf)_f = gamma_f*k*nabla(vf)_f
which is just the second part on the right side of formula 3.28 in Hrv's thesis multiplied by the interpolated diffusion coefficient gamma_f.
The content of the faceFluxCorrectionPtr is summed up with fvc::div for the faces of every cell. It is multiplied with the cell volume (why?) and goes to the right hand side of the equation into the source term with a minus.
If the flux is calculated from the matrix (line 938 in fvMatrix.C) the content of the faceFluxCorrectionPtr is just added.

That's what i think is happening for the corrected laplacian of a scalar equation. If you think i got something wrong, i'll be pleased if you propose your correction.

Now i have some questions.

My first question is as follows: It has often been written in this forum that the approach taken in OF is the overrelaxed approach described in Hrv's PhD Thesis. But if what i wrote above is right, only the explicit part of the overrelaxed approach is used and the implicit part (the matrix coefficients) is taken from the orthogonal correction approach (they should be calculated according to the expression with (phi_N-phi_P) in formula (3.33) and (3.34) in Hrv's thesis, multiplied by gamma_f). But then in gaussLaplacianSchemes.C line 56, fvmLaplacianUncorrected would have to be called with gamma*|s-k|, which is gamma*|Delta| with the Delta defined according to Hrv's thesis, formula (3.32). Does anybody know why the face area is taken instead of the magnitude of the Delta vector?

My second question seems quite simple to answer (and i hope will turn out to be stupid:) : In gaussLaplacianSchemes.C (line 70) the divergence of the corrections is multiplied by the cell volume as far as i can see. But where does this cell volume come from? In my opinion this multiplication by the cell volume should not be there (i hope this question can easily be answered...). I found no way to formulate the discretization of a transport equation in such a way that the explicit flux contribution is multiplied by the volume. Can anybody give me a hint?

Third question: In fvmLaplacianUncorrected there is a coefficient with a positive sign written into the of-diagonal matrix entries and a contribution with negative sign into the diagonal entries. I would have expexted the signs of the diagonal and the of-diagonal entries to be the other way around (of-diagonal entry negative, diagonal positive, as described e.g. in the book by Ferziger and Peric). Can anybody explain?

I hope that my question is understandable. If my notation is not clear, feel free to ask.

Thanks in advance for your answers or comments.

Hi koderer,

Thank you for your clarifying this topic.
I'm also trying to fully understand the discretization of laplacian term.

Below are my answers to your questions:
1) This point is not clear to me. Can somebody help me with it?

2) fvc::div(*fvm.faceFluxCorrectionPtr())().internalF ield() is already divided by
the cell volume, so multiplying this term by the cell volume equals
Σ_f (gamma_f * k & nabla(vf)_f).

3) The orthogonal contribution is evaluated by gamma_f*|s|*((phi_N-phi_P))/|d|,
so I think the off-diagonal coefficient is gamma_f*|s|/|d| and the diagonal
coefficient is -Σ_f (gamma_f*|s|/|d|) and the code is correct.

Please correct the mistakes, if any.

Hope that helps,
fumiya

Thank you for your answer, fumiya.
1) I'm glad i'm not the only one wondering.
2) Can you tell me the code lines where this is done? I am still not sure, because in the implicit counter part fvm::div is always used in a way that makes only sense without a built in division by the cell volume.
3) What you write seems reasonable, but the code lines i refer to seem to say the opposite. Can somebody help?

Hi koderer,

1)
In gaussLaplacianScheme.C

Code:

---

00052    tmp<surfaceScalarField> tdeltaCoeffs =
00053        this->tsnGradScheme_().deltaCoeffs(vf);
00054    const surfaceScalarField& deltaCoeffs = tdeltaCoeffs();

---

and in .H file of each snGradSchemes (e.g. correctedSnGrad.H)

Code:

---

00093        //- Return the interpolation weighting factors for the given field
00094        virtual tmp<surfaceScalarField> deltaCoeffs
00095        (
00096            const GeometricField<Type, fvPatchField, volMesh>&
00097        ) const
00098        {
00099            return this->mesh().nonOrthDeltaCoeffs();
00100        }

---

,so I think the deltaCoeffs used to calculate the off-diagonal coefficients is corrected
one to take into consideration the mesh non-orthogonality. Therefore, the coefficients
calculated in the following line are gamma_f * |Delta|/|d|, where Delta is the vector
defined in the Dr. Jasak's thesis(Eq. (3.32)).

Code:

---

00066    fvm.upper() = deltaCoeffs.internalField()*gammaMagSf.internalField();

---

2)
I think the fvc::div operation to the surfaceScalarField is caluculated using the
surfaceIntegrate function defined in the fvcSurfaceIntegrate.C and in this file

Code:

---

00075    ivf /= mesh.V();

---

When we solve the pressure poisson equation

fvm::laplacian(rAU, p) == fvc::div(phi)

, the fvc::div(phi) term is calculated in the same way. For explicit
source terms, they are incorporated into an equation simply as a field
of values. In this example fvc::div(phi) is volScalarField(not multiplied by the cell volume).