[S.Colucci]

Hi all,
I'm looking into the last release of the solver reactingEulerFoam 3.0.x. Mass transfer term (dmdt) associataed with a change of phase is implemented in function massTransfer() in InterfaceCompositionPhaseChangePhaseSystem.C:

dmdtExplicit += dmdtSign*phase.rho()*KD*Yf;
dmdt -= dmdtSign*phase.rho()*KD*eqns[name]->psi();

where dmdt is initialized with dmdtExplicit, thus

dmdt = dmdtSign*rho*KD*( Yf - Yi ),

where Yf is the equilibrium concentration of component i in the phase (given, for example by the Henry's law) and Yi the effective transported concentration of component i in the phase. dmdtSign is "+" for phase 1 and "-" for phase 2, thus for phase 1

dmdt = dmdtSign*rho*KD*( Yf - Yi ),

and for phase 2

dmdt = dmdtSign*rho*KD*( Yi - Yf ).

It means that for phase 1, as well as for phase 2, a mass transfer from 2 to 1 is positive and a mass transfer from 1 to 2 is negative.

In reactingTwoPhaseEulerFoam/pUf/UEqns.H

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

and in U1Eqn

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

while in U2Eqn

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

I think it should be the opposite:

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

because the velocity associated to a mass transfer from phase 1 to phase 2 is U1 and viceversa.

Looking into HeatAndMassTransferPhaseSystem.C, where the momentum transfer term is added to the UEqns for the case when faceMomentum is off, these terms seem to be implemented in the right way. Finally, to verify this, I run a simulation where liquid water is oversaturated in dissolved air in a closed column and I printed on the screen at runtime the value of dmdt. As you can see in the figure, the dissolved air start to exsolve and I it moves from the liquid (phase2 in this case) to the gas phase. I checked that the dmdt, printed on the screen, and associated to a mass transfer from 2 to 1 is positive (it goes in dmdt12, according to the actual implementation).

Should I report this as a bug in the OpenFOAM bug Reporting?

Simone

[S.Colucci]

Sorry, there was a mistake in the correction that I was proposing , the right is

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

in this way the sign of the mass transfer term in the UEqn is right.

[S.Colucci]

Recently, a new version of OpenFOAM has been released (OpenFOAMv1606+). In this new version the mass transfer terms in pUf/UEqns.H have been modified. Now it reads:

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

{
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);
}

The momentum transfer term, dmdt, now is in agreement with the implementation for the case when faceMomentum is off (see HeatAndMassTransferPhaseSystem.C).
But there is still something that I do not understand:

1) in pUf/UEqns.H I read

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

but I expected

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

is there some other dmdt added somewhere in the code? Where is it?

2) in the U1Eqn above, as well as in U2Eqn, where are temporal derivatives, other momentum transfer terms and continuity errors?

Any idea?


Simone