View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0002141  OpenFOAM  Bug  public  20160707 11:44  20160707 14:39 
Reporter  sim81  Assigned To  henry  
Priority  normal  Severity  major  Reproducibility  N/A 
Status  resolved  Resolution  fixed  
Platform  GNU/Linux  OS  Ubuntu  OS Version  14.04 
Summary  0002141: mass transfer and continuity error in reactingTwoPhaseEulerFoam  
Description  1) 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.  
Tags  No tags attached.  

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

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) 

What is your proposal to fix the inconsistency? 

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. 

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

Do you have a testcase which demonstrates the need for change to either formulation? 

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. 

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) 

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 cellbased formulation which uses the equivalent 'continutityError_' term. 

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

> 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 testcase to demonstrate the error. 

Resolved in OpenFOAMdev by commit 5571dbf107cd6860aee7d5484001105e57063e12 Resolved in OpenFOAM4.x by commit 3809df37fa0aaeb6fe759bcac3039e9fba318e1b 
Date Modified  Username  Field  Change 

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