[cvr]

Hello,

I am studying a bit of some pressure-coupling procedures used in OpenFOAM. For the moment I am using version 4.x in github.

I've come across some parts of the code which I'm not fully grasping. In several transient solvers (including icoFoam, which I'll use as reference here for simplicity) the phiHbyA field is the flux of HbyA (thus interpolation to faces
and dot product with respective face normal vector) and a correction to avoid an unwanted dependency on the time-step (following Choi (1999) Numerical Heat Transfer A, 36:545-550 and discussed also in Yu et al (2002) Numerical Heat Transfer B, 42:141-166). In solvers/incompressible/icoFoam/icoFoam.C we find:

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",  
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

So we have phiHbyA as the flux of HbyA and a correction related to the times cheme used, ddtCorr(U, phi). Correct me if I'm wrong: in practice as as this is an incompressible solver, phi = Uf·Sf, meaning the mass flux and so it is a quantity defined at each face (i.e. phi = fvc::interpolate(U) & mesh.Sf()).

In finiteVolume/finiteVolume/fvc/fvcDdt.C:

ddtCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
)
{
    return fv::ddtScheme<Type>::New
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    ).ref().fvcDdtUfCorr(U, Uf);
}

Following the source code and going to a simple time scheme such as Euler to better grasp what is going on, in finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C:

EulerDdtScheme<Type>::fvcDdtUfCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
)
{
    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();

    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
    fluxFieldType phiCorr
    (
        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
    );

    return tmp<fluxFieldType>
    (
        new fluxFieldType
        (
            IOobject
            (
                "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
                mesh().time().timeName(),
                mesh()
            ),
            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
           *rDeltaT*phiCorr
        )
    );
}

Keeping the original names of the arguments, in the above expressions we start with (in pseudo-code style):
1. ddtCorr(U, phi), which leads to
2. fvcDdtUfCorr(U, phi), then to
3. Sf & phi.oldTime, and lastly
4. Sf & phi.oldTime - Sf & interpolate(U.oldTime)

Considering phi = Sf & interpolate(U), then the 3rd and 4th step would be substituted by Sf & (Sf & interpolate(U).oldTime), thus phi is being multiplied twice by Sf, which seems incorrect.

Can someone comment on this to try to understand what I'm failing to grasp here?

[cvr]

I have filed two issues regarding this situation which were dismissed (
http://bugs.openfoam.org/view.php?id=3000 and http://bugs.openfoam.org/view.php?id=3001), on the basis of:

> thus phi is being multiplied twice by Sf, which seems incorrect.
This would generate a dimension check error.
> So either there is something I'm failing to grasp
Right, this is a user support request.

and,

For static meshes ddtCorr must be based on phi, there is no need for Uf which is only required for moving meshes. Removing fvcDdtPhiCoeff is not an option as it would cause many case to crash.
Also
Sf & phi.oldTime
is simply not possible, Sf is a surfaceVectorField and phi is a surfaceScalarField, there is no "dot" product defined for these two.

In the meanwhile I have made a test case (in attached file testIcoFoamDdtCorr.tar.gz) where icoFoam is changed to compute Uf (velocity vector at faces) which is then used as argument to ddtCorr(U, Uf) instead of ddtCorr(U, phi).

With this I should get a dimension check error. I did not obtained any error while compiling the solver. I also didn't got any error when running a modified cavity tutorial (also in testIcoFoamDdtCorr.tar.gz). I would like for someone to confirm these results.

This seems to be problematic for the following reasons:
1. If my analysis is incorrect, I should not have been able to run the solver without obtaining a dimension check error.

2. If my analysis is correct (meaning ddtCorr(U, Uf) is correct and ddtCorr(U, phi) not), then I should obtain a dimension check error while running standard icoFoam. This is the same as saying that function ddtCorr(U, phi) is called, in turn calling fvcDdtUfCorr(U, phi) and therein phiUf0 would have dimensions of L^5/T instead of L^3/T. Due to subtraction phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()), this should give an obvious error as dimensionally it doesn't make sense.

3. Again, if my analysis is correct, phiUf0 is computed with. mesh().Sf() & Uf.oldTime(), thus a dot product between a vector and a scalar. It doesn't make sense.

4. The fvcDdtPhiCoeff function returns a factor which is not part of the Choi (1999) method. it was laready pointed out in this forum that this is an ad-hoc correction to avoid instabilities if the full Choi (1999) method is used. If it is confirmed that phiUf0 computation is not correct, I hypotethize this coefficient is compensating for the extra face area being multiplied and is in fact useless if the calculation is made correct.

[mAlletto]

Originally Posted by cvr

Hello,

I am studying a bit of some pressure-coupling procedures used in OpenFOAM. For the moment I am using version 4.x in github.

I've come across some parts of the code which I'm not fully grasping. In several transient solvers (including icoFoam, which I'll use as reference here for simplicity) the phiHbyA field is the flux of HbyA (thus interpolation to faces
and dot product with respective face normal vector) and a correction to avoid an unwanted dependency on the time-step (following Choi (1999) Numerical Heat Transfer A, 36:545-550 and discussed also in Yu et al (2002) Numerical Heat Transfer B, 42:141-166). In solvers/incompressible/icoFoam/icoFoam.C we find:

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",  
    fvc::flux(HbyA)
  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

So we have phiHbyA as the flux of HbyA and a correction related to the times cheme used, ddtCorr(U, phi). Correct me if I'm wrong: in practice as as this is an incompressible solver, phi = Uf·Sf, meaning the mass flux and so it is a quantity defined at each face (i.e. phi = fvc::interpolate(U) & mesh.Sf()).

In finiteVolume/finiteVolume/fvc/fvcDdt.C:

ddtCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
)
{
    return fv::ddtScheme<Type>::New
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    ).ref().fvcDdtUfCorr(U, Uf);
}

Following the source code and going to a simple time scheme such as Euler to better grasp what is going on, in finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C:

EulerDdtScheme<Type>::fvcDdtUfCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
)
{
    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();

    fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
    fluxFieldType phiCorr
    (
        phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
    );

    return tmp<fluxFieldType>
    (
        new fluxFieldType
        (
            IOobject
            (
                "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
                mesh().time().timeName(),
                mesh()
            ),
            this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
           *rDeltaT*phiCorr
        )
    );
}

Keeping the original names of the arguments, in the above expressions we start with (in pseudo-code style):
1. ddtCorr(U, phi), which leads to
2. fvcDdtUfCorr(U, phi), then to
3. Sf & phi.oldTime, and lastly
4. Sf & phi.oldTime - Sf & interpolate(U.oldTime)

Considering phi = Sf & interpolate(U), then the 3rd and 4th step would be substituted by Sf & (Sf & interpolate(U).oldTime), thus phi is being multiplied twice by Sf, which seems incorrect.

Can someone comment on this to try to understand what I'm failing to grasp here?

I think there is a missunderstanging:

there are two different functions defined in https://github.com/OpenFOAM/OpenFOAM...lerDdtScheme.C

with the same name. In C++ if there are two different functions with the same name but different arguments, the arguments count:

template<class Type>tmp<typename localEulerDdtScheme<Type>::fluxFieldType>localEulerDdtScheme<Type>::fvcDdtUfCorr
 
  (const GeometricField<Type, fvPatchField, volMesh>& U,const
 GeometricField<Type, fvsPatchField, surfaceMesh>& Uf)

the above function takes in the case of the piso or icofoam the velocity an the face velocity as argument.

There is another function in the same file which takes the velocity and the mass flux as argument:

[CODE]
template<class Type>tmp<typename
localEulerDdtScheme<Type>::fluxFieldType>localEule rDdtScheme<Type>::fvcDdtPhiCorr

(const GeometricField<Type, fvPatchField, volMesh>& U,const
fluxFieldType& phi)
[\CODE]

which is used actually I think. Correct me if I'm wrong

[cvr]

Thank you for your reply. If I understand you are saying that in fvcDdt.C,

template<class Type>
tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh>>
ddtCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
)
{
    return fv::ddtScheme<Type>::New
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    ).ref().fvcDdtUfCorr(U, Uf);
}


template<class Type>
tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh>>
ddtCorr
(
    const GeometricField<Type, fvPatchField, volMesh>& U,
    const GeometricField
    <
        typename flux<Type>::type,
        fvsPatchField,
        surfaceMesh
    >& phi
)
{
    return fv::ddtScheme<Type>::New
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    ).ref().fvcDdtPhiCorr(U, phi);
}

the solver distinguishes between both ddtCorr formulations based on the type of the argument ("surfaceVectorField" for Uf and "surfaceScalarField" for phi). Is that correct?

[mAlletto]

exactly. see, e.g. http://www.cplusplus.com/doc/tutorial/functions2/


so if you provide a volVectorfield (U) and a surfacescalarfield (phi) as arguments the compiler chooses one function and if you provide U and a surfaceVectorField Uf as argument the compiler chooses another function.

[cvr]

Thank you for your help on this. Always learning.

[Elol]

Hi everybody,

I think this is a very productive discussion however I have some diffuclties to understand:

1- why we use {fvc::ddtCorr(U, phi)} ?

2- what the flux() function do; because what i see here they applied this function to (HbyA)
what I can conclude here that they interpolate the matrix (H/A) to the cell faces!!! is this correct?

3- If my conclusion in question 2 is correct. What is the use of function interpolate (rAU)?
because what i see here is interpolating for the A^-1 (the same function as flux())

Maybe better to say what is the difference between flux() and interpolate() functions?

So; at the end it is something like this (H/A)_f + ( (1/A)_f * ddtCorr(U,phi)

[mAlletto]

If I'm not completely mistaken {fvc::ddtCorr(U, phi)} includes a correction to the rhie and chow interpolation to make it independent of the time step. It increases additionally the pressure velocity coupling. You may find some the derivation in the book of https://www.springer.com/gp/book/9783319168739


For the flux() method you find a description here: Description of flux() method


interpolate() interpolates the variables from the cell centers to the face centers

[Elol]

Thank you mAlletto for reply I understand now the difference between the flux() and interpolation().

I have another question. it would be great if you have the answer ,

Piso algorithm work as follows (correct me if I am wrong) :

1. predicted velocity MU = - grad p

2. seperate the the M matrix to diagonal A and off-diagonal H

3. AU - H = - grad p

4. U = A^-1 H - A^-1 grad (p) this equation at the center interpolate it to the face (U)_f =(A^-1 H)_f - (A^-1 grad (p))_f

5. Enforce contnuity div U =0 ... so equ in 4 will be div (A^-1 grad (p))_f =div (A^-1 H)_f

6. from 5 we calculate the pressure field

7. we update the velocity again from 4

8. start the looop again until the convergence.

my question is what is the meaning of the phiHbyA in this sequence?

[mAlletto]

phiHbyA is the flux computed with Ueqn.H(). you need it when you derive the pressure equation:


U = H/A - grad(p)


If you make the divergence of the above equation and apply a finite volume discretisatoin you get:







You can find a quite detailed description of the pressure correction implemented in OF here: https://openfoamwiki.net/index.php/SimpleFoam. It's only simpleFoam but the principle is the same. For time dependent simulation you have the time step size in the A coefficients and some other terms in the source therm which are incorporated in the H() function.


Best


Michael

[Jesper_Roland]

Originally Posted by Elol

Thank you mAlletto for reply I understand now the difference between the flux() and interpolation().

I have another question. it would be great if you have the answer ,

Piso algorithm work as follows (correct me if I am wrong) :

1. predicted velocity MU = - grad p

2. seperate the the M matrix to diagonal A and off-diagonal H

3. AU - H = - grad p

4. U = A^-1 H - A^-1 grad (p) this equation at the center interpolate it to the face (U)_f =(A^-1 H)_f - (A^-1 grad (p))_f

5. Enforce contnuity div U =0 ... so equ in 4 will be div (A^-1 grad (p))_f =div (A^-1 H)_f

6. from 5 we calculate the pressure field

7. we update the velocity again from 4

8. start the looop again until the convergence.

my question is what is the meaning of the phiHbyA in this sequence?

I think, Tessa Uroic has described the PISO algorithm implemented in OpenFOAM very well in her PhD thesis, chapter 2:


https://www.fsb.unizg.hr/atlantis/we...s_1oct2019.pdf

Cheers,

Jesper