View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002141OpenFOAMBugpublic2016-07-07 11:442016-07-07 14:39
Reportersim81 Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityN/A
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Summary0002141: mass transfer and continuity error in reactingTwoPhaseEulerFoam
Description1) When faceMomentum is off, dmdt12 and dmdt21 are defined in HeatAndMassTransferPhaseSystem.C as

        const volScalarField dmdt21(posPart(dmdt));
        const volScalarField dmdt12(negPart(dmdt));

When faceMomentum is on, they are defined in the opposite way (UEqns.H)

    const volScalarField dmdt12(posPart(fluid.dmdt()));
    const volScalarField dmdt21(negPart(fluid.dmdt()));

2) When faceMomentum is off, the continuity error is added to UEqn in MovingPhaseModel.C

template<class BasePhaseModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
{
    return
    (
        fvm::ddt(*this, this->thermo().rho(), U_)
      + fvm::div(alphaRhoPhi_, U_)
      - fvm::Sp(continuityError_, U_)
      + this->fluid().MRF().DDt(*this*this->thermo().rho(), U_)
      + turbulence_->divDevRhoReff(U_)
    );
}

When faceMomentum is on, I do not see any continuity error in the U1Eqn and U2Eqn

    {
        U1Eqn =
        (
            fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
          + fvm::Sp(dmdt12, U1) - dmdt12*U2
          + MRF.DDt(alpha1*rho1, U1)
          + phase1.turbulence().divDevRhoReff(U1)
          + Vm*(UgradU1 - (UgradU2 & U2))
        );
        U1Eqn.relax();
        fvOptions.constrain(U1Eqn);
        U1.correctBoundaryConditions();
        fvOptions.correct(U1);
    }

Without continuity error I suppose that mass transfer terms (+ fvm::Sp(dmdt12, U1) - dmdt12*U2) are not implemented in the right way.
TagsNo tags attached.

Activities

henry

2016-07-07 11:49

manager   ~0006494

> Without continuity error I suppose that mass transfer terms (+ fvm::Sp(dmdt12, U1) - dmdt12*U2) are not implemented in the right way.

Why do you suppose it is not implemented in the right way?

sim81

2016-07-07 12:15

reporter   ~0006495

Because

ContinuityError = dmdt

that is the continuity equation with phase change, so when you subtract the ContinuityError in the UEqn (as done for the case when faceMomentum is off in the file MovingPhaseModel.C) you have (for phase 1)

dmdt12*U1 - dmdt12*U2 - ContinuityError*U1 = dmdt12*U1 - dmdt12*U2 - dmdt*U1 = dmdt12*U1 - dmdt12*U2 - (dmdt12 + dmdt21)*U1 = dmdt21*U1- dmdt12*U2

thus the mass transferred from 2 -> 1 is multiplied by U2 and gives a positive contribution (on the rhs) to phase and the mass transferred from 1 -> 2 is multiplied by U1 and gives a negative contribution (on the rhs)

henry

2016-07-07 12:30

manager   ~0006496

What is your proposal to fix the inconsistency?

sim81

2016-07-07 13:25

reporter   ~0006497

in pUf/UEqn.H we can add the ContinuityError in U1Eqn and U2Eqn, in analogy to the faceMomentum off case (MovingPhaseModel.C)

        U1Eqn =
        (
            fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
          + fvm::Sp(dmdt12, U1) - dmdt12*U2
          - fvm::Sp(phase1.continuityError(), U1)
          + MRF.DDt(alpha1*rho1, U1)
          + phase1.turbulence().divDevRhoReff(U1)
          + Vm*(UgradU1 - (UgradU2 & U2))
        );

or we can write instead of

+ fvm::Sp(dmdt12, U1) - dmdt12*U2

directly this

 fvm::Sp(dmdt21, U1) - dmdt12*U2.

henry

2016-07-07 13:53

manager   ~0006498

The current pUf/UEqn.H implementation has been checked and tested, see commit 1b9d8d295464450bde5f83d85320ad80936c3a6d

henry

2016-07-07 13:56

manager   ~0006499

Do you have a test-case which demonstrates the need for change to either formulation?

sim81

2016-07-07 14:14

reporter   ~0006500

Sorry I don't know how to provide a valid test case to show this, but checking the mathematics something is missing. If you add the continuity error in the pUf/UEqns.H, in analogy to the same equation implemented when faceMomentum is off (in MovingPhaseModel.C), the formulation seems to be correct.

sim81

2016-07-07 14:25

reporter   ~0006501

In fact, if you subtract the continuity error on the left side of the U1Eqn, as I shown here in the second post, you obtain

-dmdt21*U1- dmdt12*U2

thus the mass transferred from 2 -> 1 is multiplied by U2 and gives a positive contribution (on the rhs) to phase 1 and the mass transferred from 1 -> 2 is multiplied by U1 and gives to phase 1 a negative contribution (on the rhs)

henry

2016-07-07 14:27

manager   ~0006502

The current implementations have been tested and demonstrated to be correct. The only difference I see is in the naming convension for the dmdt terms and I propose to change pUf/UEqns.H to

    const volScalarField dmdt21(posPart(fluid.dmdt()));
    const volScalarField dmdt12(negPart(fluid.dmdt()));

    {
        U1Eqn =
        (
            fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
          + fvm::Sp(dmdt21, U1) - dmdt21*U2
          + MRF.DDt(alpha1*rho1, U1)
          + phase1.turbulence().divDevRhoReff(U1)
          + Vm*(UgradU1 - (UgradU2 & U2))
        );
        U1Eqn.relax();
        fvOptions.constrain(U1Eqn);
        U1.correctBoundaryConditions();
        fvOptions.correct(U1);
    }

    {
        U2Eqn =
        (
            fvm::div(alphaRhoPhi2, U2) - fvm::Sp(fvc::div(alphaRhoPhi2), U2)
          - fvm::Sp(dmdt12, U2) + dmdt12*U1
          + MRF.DDt(alpha2*rho2, U2)
          + phase2.turbulence().divDevRhoReff(U2)
          + Vm*(UgradU2 - (UgradU1 & U1))
        );
        U2Eqn.relax();
        fvOptions.constrain(U2Eqn);
        U2.correctBoundaryConditions();
        fvOptions.correct(U2);
    }

In the above note carefully the ' - fvm::Sp(fvc::div(alphaRhoPhi1), U1)' terms
which are not present in the cell-based formulation which uses the equivalent 'continutityError_' term.

sim81

2016-07-07 14:33

reporter   ~0006503

Does it mean that there is non continuity error correction in the pUf/UEqns formulation? Or it is elsewhere in the code?

henry

2016-07-07 14:35

manager   ~0006505

> Does it mean that there is non continuity error correction in the pUf/UEqns formulation?

Yes, I believe the formulations are consistent, apart from the names which I will fix. If you still believe that either of the formulations are not correct you will need to provide a test-case to demonstrate the error.

henry

2016-07-07 14:39

manager   ~0006506

Resolved in OpenFOAM-dev by commit 5571dbf107cd6860aee7d5484001105e57063e12
Resolved in OpenFOAM-4.x by commit 3809df37fa0aaeb6fe759bcac3039e9fba318e1b

Issue History

Date Modified Username Field Change
2016-07-07 11:44 sim81 New Issue
2016-07-07 11:49 henry Note Added: 0006494
2016-07-07 12:15 sim81 Note Added: 0006495
2016-07-07 12:30 henry Note Added: 0006496
2016-07-07 13:25 sim81 Note Added: 0006497
2016-07-07 13:53 henry Note Added: 0006498
2016-07-07 13:56 henry Note Added: 0006499
2016-07-07 14:14 sim81 Note Added: 0006500
2016-07-07 14:25 sim81 Note Added: 0006501
2016-07-07 14:27 henry Note Added: 0006502
2016-07-07 14:33 sim81 Note Added: 0006503
2016-07-07 14:35 henry Note Added: 0006505
2016-07-07 14:39 henry Note Added: 0006506
2016-07-07 14:39 henry Status new => resolved
2016-07-07 14:39 henry Fixed in Version => 4.x
2016-07-07 14:39 henry Resolution open => fixed
2016-07-07 14:39 henry Assigned To => henry