1D Stefan problem validation of new solver
 

Hello everyone,

I developed a new evaporation solver based on some excellent solvers (interFoam, interPhaseChangeFoam and evapVOFHardt), and when I validated the 1D Stefan problem, the velocity field is not correct (see attached), and some key code can be found as follows:
1. Source terms in alphaEqn (the same as the source terms in interPhaseChangeFoam)

// explicit source term
volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
- dotM/rho1
);


//zeroField Sp;

// implicit source term
volScalarField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
fvc::div(phi) - dotM*(1.0/rho2-1.0/rho1)
);

2. UEqn with new source terms
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) // new
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);

3. TEqn

Info << "Solve TEqn" << endl;
{
volScalarField k = alpha1*k1 + alpha2*k2; // one-field formulation
volScalarField rhoCp = alpha1*rho1*cp1 + alpha2*rho2*cp2;
// one-field formulation; rho:density; cp: specific heat capacity

surfaceScalarField alphaPhi = (rhoPhi-phi*rho2)/(rho1-rho2);
// hint: rhoPhi = alpha1*rho1*phi + (1-alpha1)*rho2*phi
// rho1 must not be equal to rho2
surfaceScalarField rhoCpPhi = alphaPhi*(rho1*cp1-rho2*cp2)+phi*rho2*cp2;
// rhoCpPhi = alpha1*rho1*cp1*phi + (1-alpha1)*rho2*cp2*phi

fvScalarMatrix TEqn
(
fvm::ddt(rhoCp,T)
+ fvm::div(rhoCpPhi,T)
- fvm::laplacian(k,T)
==
fvm::Sp(fvc::ddt(rhoCp) + fvc::div(rhoCpPhi),T)
// new source term 01 (Coupled to mass source term)
+ fvm::Sp(TEqnSou,T)
- TEqnSou*Tsat
// new source term 02 (cooling due to evaporation)
);

TEqn.relax();
TEqn.solve();

Info<< "min/max(T) = " << min(T).value() << "; "
<< max(T).value() <<endl;
}

4. pEqn
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) - dotM*(scalar(1.0)/rho2 - scalar(1.0)/rho1) //new
);

5. source terms (modified from evapVOFHardt)

volScalarField jLoc = (T-Tsat)/(Rint*hEvap);

dotM = jLoc*magGradAlpha;

forAll(alpha1.ref(),iCell)
{
if (magGradAlpha.ref()[iCell] > 0.0)
{
dotM[iCell] *= 1.0;
// only set source terms in interface cells
}

else
{
dotM[iCell] *= 0.0;
}
}

// the second term of the TEqn
TEqnSou = -magGradAlpha/Rint;

6. main program
#include "sourceTerms.H"
evapRate = fvc::domainIntegrate(dotM);

#include "TEqn.H"

Info<< "Evaporation rate is = " << evapRate.value() << "kg/s at t = " << runTime.timeName()
<< "s" << endl;

If anyone know how to fix the problem, it will be best, and if I have something wrong in my code, please let me know.

Many thanks,
Huihui

Hi,

You work on 1D twho-phase-flow probleme, , Me too, i work on 1D two-phase-flow but compressible configuration . I hope that you resolve your problem



Please, I would like to know how you chose your discretizations (fvSchemes and fvSolution) for your problem ??

best regret,