View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003001OpenFOAMBugpublic2018-07-05 18:402018-07-05 19:56
Reportercvr Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilitysometimes
Status closedResolutionno change required 
OSDebianOS Version9 
Summary0003001: Pressure-velocity coupling in transient simulation: ddtCorr(U, phi) and fvcDdtUfCorr
DescriptionThis is a re-submission of http://bugs.openfoam.org/view.php?id=3000. Please, do check the "Steps To Reproduce" as I have put a test case focusing this issue.


I am studying a bit of some pressure-coupling procedures used in OpenFOAM. For the moment I am using version 4.x in github but I have confirmed that version 5.x shares the same formulation.

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:
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.

So either there is something I'm failing to grasp or there is an erroneous dot product of a scalar (phi) with the face-area vector (Sf).
Steps To ReproduceThe manager of this issue pointed out (and very well):

> > 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.

As such I've made a test based on the "icoFoam" solver and the "cavity" tutorial. The problematic part is that I got no "dimension check error", which may indicate further problems.

I made a version of "icoFoam" where I replaced "phi" by the face-velocity vector "Uf" and obtained no "dimension check error", neither while compiling with wmake nor when testing with the "cavity" tutorial.

I submit a file (testIcoFoamDdtCorr.tar.gz) with the solver and tutorial I changed to make this report, for anyone who's interested in reproducing this.


DETAILED STEPS
==============

I copied the icoFoam solver into testIcoFoamDdtCorr, moved icoFoam.C to testIcoFoamDdtCorr.C and changed ./Make/files accordingly to use testIcoFoamDdtCorr.C. I changed/added the following lines:
_______________________________________
    #include "createFields.H"
+ #include "createUf.H" // ADDED TO COMPUTE Uf
    #include "initContinuityErrs.H"
...
+ //+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
+ + fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf) // USES Uf INSTEAD OF phi
...
            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
+ Uf = fvc::interpolate(U); // COMPUTES Uf
_______________________________________

I have compiled it successfully (with no errors nor warnings). My g++ version is (Debian 6.3.0-18+deb9u1) 6.3.0 20170516.

I have copied the cavity tutorial to testIcoFoamDdtCorr_cavity and
changed system/controlDict to have "application testIcoFoamDdtCorr;" instead of "application icoFoam;".

I have run the case with the "testIcoFoamDdtCorr" solver with no error nor warning whatsoever, throughout the whole simulation time. I put an example of the output:q


/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 4.x-d214c8dfd5ba
Exec : testIcoFoamDdtCorr
Date : Jul 05 2018
Time : 18:26:10
Host : "niflheim"
PID : 21793
Case : /home/cvr/OpenFOAM/cvr-4.x/run/testIcoFoamDdtCorr_cavity
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


PISO: Operating solver in PISO mode

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi

Reading/calculating face velocity Uf


Starting time loop

Time = 0.005

Courant Number mean: 0 max: 0
smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 8.90511e-06, No Iterations 19
smoothSolver: Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.0492854, No Iterations 12
time step continuity errors : sum local = 0.000466513, global = -1.79995e-19, cumulative = -1.79995e-19
DICPCG: Solving for p, Initial residual = 0.590864, Final residual = 2.65225e-07, No Iterations 35
time step continuity errors : sum local = 2.74685e-09, global = -2.6445e-19, cumulative = -4.44444e-19
ExecutionTime = 0.01 s ClockTime = 0 s

Time = 0.01

Courant Number mean: 0.0976825 max: 0.585607
smoothSolver: Solving for Ux, Initial residual = 0.160686, Final residual = 6.83031e-06, No Iterations 19
smoothSolver: Solving for Uy, Initial residual = 0.260828, Final residual = 9.65939e-06, No Iterations 18
DICPCG: Solving for p, Initial residual = 0.426283, Final residual = 0.0104012, No Iterations 22
time step continuity errors : sum local = 0.000111592, global = -5.98217e-19, cumulative = -1.04266e-18
DICPCG: Solving for p, Initial residual = 0.29954, Final residual = 5.20644e-07, No Iterations 33
time step continuity errors : sum local = 6.59488e-09, global = 1.40538e-19, cumulative = -9.02123e-19
ExecutionTime = 0.01 s ClockTime = 0 s

...

Time = 0.5

Courant Number mean: 0.222133 max: 0.852126
smoothSolver: Solving for Ux, Initial residual = 2.31861e-07, Final residual = 2.31861e-07, No Iterations 0
smoothSolver: Solving for Uy, Initial residual = 5.06495e-07, Final residual = 5.06495e-07, No Iterations 0
DICPCG: Solving for p, Initial residual = 9.86503e-07, Final residual = 9.86503e-07, No Iterations 0
time step continuity errors : sum local = 1.00346e-08, global = -4.66199e-19, cumulative = -6.50743e-19
DICPCG: Solving for p, Initial residual = 1.0939e-06, Final residual = 2.75801e-07, No Iterations 1
time step continuity errors : sum local = 4.12211e-09, global = 3.25909e-20, cumulative = -6.18152e-19
ExecutionTime = 0.08 s ClockTime = 70 s

End
Additional InformationIf confirmed, this affects all of the time *DdtScheme.C files in finiteVolume/finiteVolume/ddtSchemes/, namely:

CoEulerDdtScheme/CoEulerDdtScheme.C
EulerDdtScheme/EulerDdtScheme.C
backwardDdtScheme/backwardDdtScheme.C
SLTSDdtScheme/SLTSDdtScheme.C
CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
localEulerDdtScheme/localEulerDdtScheme.C

The simplest way to solve this (without changing anything) would be to feed the face-velocity vector field to ddtCorr, so ddtCorr(U, Uf) instead of ddtCorr(U, phi). Even so, the ddtCorr would still be affected by fvcDdtPhiCoeff procedure (which was described as an ad-hoc inclusion due to the fact that fvcDdtPhiCoeff=1 would result in instabilities, vide
www.cfd-online.com/Forums/openfoam-solving/60096-ddtphicorr.html#post516511).

Another way would be to change in each of the *DdtScheme.C file the "phiUf0(mesh().Sf() & Uf.oldTime())" line to "phiUf0(Uf.oldtime())" and to
remove the call to "fvcDdtPhiCoeff".
Tagspressure-velocity coupling

Activities

cvr

2018-07-05 18:40

reporter  

testIcoFoamDdtCorr.tar.gz (3,640 bytes)

henry

2018-07-05 18:54

manager   ~0009844

Last edited: 2018-07-05 19:55

 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.

henry

2018-07-05 19:56

manager   ~0009845

User support request.

Issue History

Date Modified Username Field Change
2018-07-05 18:40 cvr New Issue
2018-07-05 18:40 cvr File Added: testIcoFoamDdtCorr.tar.gz
2018-07-05 18:40 cvr Tag Attached: pressure-velocity coupling
2018-07-05 18:54 henry Note Added: 0009844
2018-07-05 19:55 henry Note Edited: 0009844
2018-07-05 19:56 henry Assigned To => henry
2018-07-05 19:56 henry Status new => closed
2018-07-05 19:56 henry Resolution open => no change required
2018-07-05 19:56 henry Note Added: 0009845