[lihaoyan]

Hi guys,

I made a Solver for three in-compressible, isothermal immiscible fluids with
phase-change based on interPhaseChangeFoam.

I deleted the "interfaceProperties.H" in my main program"CrackingPhaseChangeFoam.C", because I don't need to capture the interface among different phases.
:
These are head files:
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
//#include "interfaceProperties.H"
#include "myphaseChangeThreePhaseMixture.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"

The solver could run normally:

PIMPLE: iteration 1
smoothSolver: Solving for alpha.a, Initial residual = 0.0125493, Final residual = 5.36539e-09, No Iterations 2
Phase-1 volume fraction = 0.654367 Min(alpha.a) = 0.653536 Max(alpha.a) = 0.655
smoothSolver: Solving for alpha.g, Initial residual = 0.193196, Final residual = 5.17285e-09, No Iterations 2
Phase-2 volume fraction = 0.195087 Min(alpha.g) = 0.15 Max(alpha.g) = 0.196981
MULES: Correcting alpha.a
MULES: Correcting alpha.a
MULES: Correcting alpha.g
MULES: Correcting alpha.g
MULES: Solving for alpha.c
MULES: Solving for alpha.c
Phase-1 volume fraction = 0.654367 Min(alpha.a) = 0.653536 Max(alpha.a) = 0.655
Phase-2 volume fraction = 0.195104 Min(alpha.g) = 0.15 Max(alpha.g) = 0.196131
Phase-3 volume fraction = 0.195122 Min(alpha.c) = 0.195 Max(alpha.c) = 0.195283
GAMG: Solving for p_rgh, Initial residual = 0.00316811, Final residual = 0.000128455, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 0.00122, Final residual = 9.3244e-05, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 0.000179942, Final residual = 4.28389e-08, No Iterations 3
ExecutionTime = 4.94 s ClockTime = 11 s

Courant Number mean: 0.0583806 max: 0.175259
deltaT = 0.005
Time = 3

But the results of fractions of three phases weren't output as nonuniform List<scalar> but in such format:
for example:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0.6";
object alpha.a;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 0 0 0 0 0 0];


internalField uniform 0.655;

boundaryField
{
inlet
{
type fixedValue;
value uniform 0.655;
}
outlet
{
type zeroGradient;
}
walls
{
type symmetry;
}
frontAndBack
{
type empty;
}
}


// ************************************************** *********************** //
The calculated fraction of alpha.a is the same as the alpha.a in the "0" file, please point out the reason of this problem, and what should I do to solve it!
Thank you!!!

As for the IOdictionary of alphas, I just modified library files "two mixtures" into "three mixtures", which located in $FOAM_SRC/transportmodels, the codes are:
#include "threePhaseMixture.H"


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

Foam::threePhaseMixture::threePhaseMixture
(
const fvMesh& mesh,
const dictionary& dict
)
:
phase1Name_(dict.get<wordList>("phases")[0]),/
phase2Name_(dict.get<wordList>("phases")[1]),
phase3Name_(dict.get<wordList>("phases")[2]),

alpha1_
(
IOobject
(
IOobject::groupName("alpha", phase1Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),

alpha2_
(
IOobject
(
IOobject::groupName("alpha", phase2Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),

alpha3_
(
IOobject
(
IOobject::groupName("alpha", phase3Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
{}

In alphaEqu.H, the calculating steps of alpha1 are:
if (MULES1Corr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
- fvm::Sp(divU, alpha1)
==

fvm::Sp(-tem,alpha1)
+ em1
);

alpha1Eqn.solve();

Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;

talphaPhi1 = alpha1Eqn.flux();
}

if (MULES1Corr)
{
talphaPhi1Corr.ref() -= talphaPhi1();

volScalarField alpha100("alpha100", alpha10);
alpha10 = alpha1;

MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi1(),
talphaPhi1Corr.ref(),
Sp,
(
divU*(alpha10 - alpha100)
+ tem*alpha10
)(),
oneField(),
zeroField()
);

// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
talphaPhi1.ref() += talphaPhi1Corr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
talphaPhi1.ref() += 0.5*talphaPhi1Corr();
}
}
else
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
talphaPhi1Corr.ref(),
Sp,
(divU*alpha1 + em1)(),
oneField(),
zeroField()
);

talphaPhi1 = talphaPhi1Corr;
}
}
Sp, -tem, em1 are sources.

[Tobermory]

Sounds like your code is simply not writing the fields out, but is simply copying the files from the 0 folder. Check the coding for the phases - how did you define them? As an IO-Object with AUTO_WRITE?

You need to share some of your coding if you want us to help you further.

[lihaoyan]

Originally Posted by Tobermory

Sounds like your code is simply not writing the fields out, but is simply copying the files from the 0 folder. Check the coding for the phases - how did you define them? As an IO-Object with AUTO_WRITE?

You need to share some of your coding if you want us to help you further.

Thank you very much for your reply!
As for the IOdictionary of the alpha, I just modified the library files "two mixtures" in to "three mixtures" in $FOAM_SRC/transportmodels, which are following:

#include "threePhaseMixture.H"


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

Foam::threePhaseMixture::threePhaseMixture
(
const fvMesh& mesh,
const dictionary& dict
)
:
phase1Name_(dict.get<wordList>("phases")[0]),/
phase2Name_(dict.get<wordList>("phases")[1]),
phase3Name_(dict.get<wordList>("phases")[2]),

alpha1_
(
IOobject
(
IOobject::groupName("alpha", phase1Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),

alpha2_
(
IOobject
(
IOobject::groupName("alpha", phase2Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),

alpha3_
(
IOobject
(
IOobject::groupName("alpha", phase3Name_),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
{}

As for the alphaEqu.H:
Sp, em1 and -tem are the sources.
the calculating steps of the first phase is:
if (MULES1Corr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
- fvm::Sp(divU, alpha1)
==

fvm::Sp(-tem,alpha1)
+ em1
);

alpha1Eqn.solve();

Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;

talphaPhi1 = alpha1Eqn.flux();
}

if (MULES1Corr)
{
talphaPhi1Corr.ref() -= talphaPhi1();

volScalarField alpha100("alpha100", alpha10);
alpha10 = alpha1;

MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi1(),
talphaPhi1Corr.ref(),
Sp,
(
divU*(alpha10 - alpha100)
+ tem*alpha10
)(),
oneField(),
zeroField()
);

// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
talphaPhi1.ref() += talphaPhi1Corr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
talphaPhi1.ref() += 0.5*talphaPhi1Corr();
}
}
else
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
talphaPhi1Corr.ref(),
Sp,
(divU*alpha1 + em1)(),
oneField(),
zeroField()
);

talphaPhi1 = talphaPhi1Corr;
}

}

[lihaoyan]

Hi I have already add the key information, please check, thank you very much for your reply.

[Tobermory]

That is very strange ... it is almost as if the version of alpha1 etc that you are solving is different from the one that you have generated in the IOobject. Can you quickly list the coding where you define alpha1, alpha2 & alpha3 (in createFields?). Either that, or something has gone wrong with the MULES implementation - you can try debug that by adding in some write statements into the alphaEqn.H to write out the value of alpha1 say in a few cells, and see how that changes through the alphaEqn.H code. Good luck.

[lihaoyan]

Originally Posted by Tobermory

That is very strange ... it is almost as if the version of alpha1 etc that you are solving is different from the one that you have generated in the IOobject. Can you quickly list the coding where you define alpha1, alpha2 & alpha3 (in createFields?). Either that, or something has gone wrong with the MULES implementation - you can try debug that by adding in some write statements into the alphaEqn.H to write out the value of alpha1 say in a few cells, and see how that changes through the alphaEqn.H code. Good luck.

Thank you very much for your advice, I will try to find out the problem. Here are the definition of alphas:

In createFields.H, I defined the alpha1,alpha2 &alpha3 through:

Info<< "Creating myphaseChangeThreePhaseMixture\n" << endl;
autoPtr<myphaseChangeThreePhaseMixture> mixture =
myphaseChangeThreePhaseMixture::New(U, phi);

volScalarField alpha1(mixture->alpha1());
volScalarField alpha2(mixture->alpha2());
volScalarField alpha3(mixture->alpha3());

"myphaseChangeThreePhaseMixture" is a class in which I defined the sources em1,Sp,tem and which inherits the "three mixture" class.

[lihaoyan]

Originally Posted by Tobermory

That is very strange ... it is almost as if the version of alpha1 etc that you are solving is different from the one that you have generated in the IOobject. Can you quickly list the coding where you define alpha1, alpha2 & alpha3 (in createFields?). Either that, or something has gone wrong with the MULES implementation - you can try debug that by adding in some write statements into the alphaEqn.H to write out the value of alpha1 say in a few cells, and see how that changes through the alphaEqn.H code. Good luck.

These are the codes of threeMixtures.H:

#ifndef threePhaseMixture_H
#define threePhaseMixture_H

#include "volFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
Class twoPhaseMixture Declaration
\*---------------------------------------------------------------------------*/

class threePhaseMixture
{
protected:

// Protected data

word phase1Name_;
word phase2Name_;
word phase3Name_;

volScalarField alpha1_;
volScalarField alpha2_;
volScalarField alpha3_;


public:

// Constructors

//- Construct from components
threePhaseMixture
(
const fvMesh& mesh,
const dictionary& dict
);


//- Destructor
~threePhaseMixture()
{}


// Member Functions

const word& phase1Name() const
{
return phase1Name_;
}

const word& phase2Name() const
{
return phase2Name_;
}

const word& phase3Name() const
{
return phase3Name_;
}

//- Return the phase-fraction of phase 1
const volScalarField& alpha1() const
{
return alpha1_;
}

//- Return the phase-fraction of phase 1
volScalarField& alpha1()
{
return alpha1_;
}

//- Return the phase-fraction of phase 2
const volScalarField& alpha2() const
{
return alpha2_;
}

//- Return the phase-fraction of phase 2
volScalarField& alpha2()
{
return alpha2_;
}

const volScalarField& alpha3() const
{
return alpha3_;
}
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************** *********************** //

[Tobermory]

Originally Posted by lihaoyan

In createFields.H, I defined the alpha1,alpha2 &alpha3 through:

Info<< "Creating myphaseChangeThreePhaseMixture\n" << endl;
autoPtr<myphaseChangeThreePhaseMixture> mixture =
myphaseChangeThreePhaseMixture::New(U, phi);

volScalarField alpha1(mixture->alpha1());
volScalarField alpha2(mixture->alpha2());
volScalarField alpha3(mixture->alpha3());

Aaaah! Is that your problem then? You define alpha1 as a copy of mixture->alpha1(), initialised by the initial value of alpha1_ ... then you update alpha1 ... but since it is a COPY and not a reference to the alpha1_ field, then you are not updating alpha1_ when you run your solver. This means that when alpha1_ is writing, it is still writing its original value.

Look back at the createFields for interPhaseChangeFoam and you will see that it defines alpha1 etc as refs to alpha1_ not as copies:

 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());

[lihaoyan]

Originally Posted by Tobermory

Aaaah! Is that your problem then? You define alpha1 as a copy of mixture->alpha1(), initialised by the initial value of alpha1_ ... then you update alpha1 ... but since it is a COPY and not a reference to the alpha1_ field, then you are not updating alpha1_ when you run your solver. This means that when alpha1_ is writing, it is still writing its original value.

Look back at the createFields for interPhaseChangeFoam and you will see that it defines alpha1 etc as refs to alpha1_ not as copies:

 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());

Hi, thank you very much for your reply! I will change the definitions of alphas in createFields.H by the examples in interPhaseChangeFoam and see if I could get the fields of alphas. Sincerely appreciated for what you did, without your help I couldn't find the problem! I will let you know the following process after the modification. Thank you!!!

[lihaoyan]

Originally Posted by Tobermory

Aaaah! Is that your problem then? You define alpha1 as a copy of mixture->alpha1(), initialised by the initial value of alpha1_ ... then you update alpha1 ... but since it is a COPY and not a reference to the alpha1_ field, then you are not updating alpha1_ when you run your solver. This means that when alpha1_ is writing, it is still writing its original value.

Look back at the createFields for interPhaseChangeFoam and you will see that it defines alpha1 etc as refs to alpha1_ not as copies:

 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());

Hi, your advice is very useful and now my solver could output the field of the alphas.
But again I faced the same problems of my "sources" in the solver, which is only copying the boundary field files from the beginning to the ending of the calculation.
I think I am not wrong with the definitions of sources m1,m2,m3 in "createFields.H", which are:

Info<< "Creating myphaseChangeThreePhaseMixture\n" << endl;
autoPtr<myphaseChangeThreePhaseMixture> mixture =
myphaseChangeThreePhaseMixture::New(U, phi);
const volScalarField& m1 = mixture->m1();
const volScalarField& m2 = mixture->m2();
const volScalarField& m3 = mixture->m3();
const volScalarField tem(m1+m2+m3);

The codes of sources are following:
Foam::volScalarField
Foam::myphaseChangeThreePhaseMixtures::cracking::m 1() const
{
const volScalarField limitedAlpha1
(
min(max(alpha1_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha2
(
min(max(alpha2_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha3
(
min(max(alpha3_, scalar(0)), scalar(1))
);
return
(
-pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho1_o_) *(K1_+K2_)
);
}

Foam::volScalarField
Foam::myphaseChangeThreePhaseMixtures::cracking::m 2() const
{
const volScalarField limitedAlpha1
(
min(max(alpha1_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha2
(
min(max(alpha2_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha3
(
min(max(alpha3_, scalar(0)), scalar(1))
);
return
(
pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho2_o_) *(K1_)
);
}

Foam::volScalarField
Foam::myphaseChangeThreePhaseMixtures::cracking::m 3() const
{
const volScalarField limitedAlpha1
(
min(max(alpha1_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha2
(
min(max(alpha2_, scalar(0)), scalar(1))
);
const volScalarField limitedAlpha3
(
min(max(alpha3_, scalar(0)), scalar(1))
);
return
(
pow(limitedAlpha1*rho1_o_,2)/((limitedAlpha1*rho1_o_+limitedAlpha2*rho2_o_+limi tedAlpha3*rho3_o_)*rho3_o_) *(K2_)
);
}
K1_ and K2_ are constants.
Please help me with it, thank you！！！