View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003000OpenFOAMBugpublic2018-07-05 16:482018-07-06 10:50
Reportercvr Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityhave not tried
Status closedResolutionno change required 
PlatformGNU/LinuxOSDebianOS Version4.9.88-1+deb9u1
Summary0003000: Pressure-velocity coupling in transient simulation: ddtCorr(U, phi) and fvcDdtUfCorr
DescriptionI 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).
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".
TagspimpleControl, pressure-velocity coupling

Activities

henry

2018-07-05 17:22

manager   ~0009843

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

Issue History

Date Modified Username Field Change
2018-07-05 16:48 cvr New Issue
2018-07-05 16:48 cvr Tag Attached: pimpleControl
2018-07-05 16:48 cvr Tag Attached: pressure-velocity coupling
2018-07-05 17:22 henry Assigned To => henry
2018-07-05 17:22 henry Status new => closed
2018-07-05 17:22 henry Resolution open => no change required
2018-07-05 17:22 henry Note Added: 0009843