[wavefunction]

Hello experts.

I am trying to incorporate source terms in the alpha equations of interMixingFoam in order to model phase change and thereby build my own solver.

According to this paper, https://pdfs.semanticscholar.org/001...78aa4ea653.pdf, I want to accomodate eqns (4), (10) and (11) into my solver. Eqn (4) describes the calculation of the source term to be added into the alpha governing eqns (10) and (11).

On trying to compile the solver, I am getting an error that says

no matching function for call to ‘Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricField()’
rho3_("rho", dimDensity, nuModel3_->viscosityProperties())

It seems that the data type of variable rho3 is changing to something else.

I have attached relevant portions of my incompressibleThreePhaseMixture.C, incompressibleThreePhaseMixture.H and alphaEqn.H files. (most relevant files) The changes made to existing files are in bold.

incompressibleThreePhaseMixture.C:

#include "incompressibleThreePhaseMixture.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"
#include "fvc.H"

// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

void Foam::incompressibleThreePhaseMixture::calcNu()
{
nuModel1_->correct();
nuModel2_->correct();
nuModel3_->correct();

// Average kinematic viscosity calculated from dynamic viscosity
nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
}

void Foam::incompressibleThreePhaseMixture::calcrhogp()
{
rhogp_ = (alpha2_*rho2_ + alpha3_*rho3_)/(alpha2_ + alpha3_);
}

void Foam::incompressibleThreePhaseMixture::calcmDot()
{
// Cell gradient of alpha
volVectorField gradAlpha(fvc::grad(alpha1_));

volScalarField Y_ = (alpha2_*rho2_)/(scalar(1) - alpha1_);

volVectorField gradY(fvc::grad(Y_));

mDot_ = (gradAlpha&gradY)/(scalar(1) - Y_);
}

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

Foam::incompressibleThreePhaseMixture::incompressi bleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
)
:
IOdictionary
(
IOobject
(
"transportProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),

phase1Name_(wordList(lookup("phases"))[0]),
phase2Name_(wordList(lookup("phases"))[1]),
phase3Name_(wordList(lookup("phases"))[2]),

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

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

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

U_(U),
phi_(phi),

nu_
(
IOobject
(
"nu",
U.time().timeName(),
U.db()
),
U.mesh(),
dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
calculatedFvPatchScalarField::typeName
),

nuModel1_
(
viscosityModel::New
(
"nu1",
subDict(phase1Name_),
U,
phi
)
),
nuModel2_
(
viscosityModel::New
(
"nu2",
subDict(phase2Name_),
U,
phi
)
),
nuModel3_
(
viscosityModel::New
(
"nu3",
subDict(phase3Name_),
U,
phi
)
),


rho1_("rho", dimDensity, nuModel1_->viscosityProperties()),
rho2_("rho", dimDensity, nuModel2_->viscosityProperties()),
rho3_("rho", dimDensity, nuModel3_->viscosityProperties())
{
alpha3_ == 1.0 - alpha1_ - alpha2_;

calcNu();
calcrhogp();
calcmDot();
}

incompressibleThreePhaseMixture.H:

namespace Foam
{

/*---------------------------------------------------------------------------*\
Class incompressibleThreePhaseMixture Declaration
\*---------------------------------------------------------------------------*/

class incompressibleThreePhaseMixture
:
public IOdictionary,
public transportModel
{
// Private data

word phase1Name_;
word phase2Name_;
word phase3Name_;

volScalarField alpha1_;
volScalarField alpha2_;
volScalarField alpha3_;

const volVectorField& U_;
const surfaceScalarField& phi_;

volScalarField nu_;

autoPtr<viscosityModel> nuModel1_;
autoPtr<viscosityModel> nuModel2_;
autoPtr<viscosityModel> nuModel3_;

dimensionedScalar rho1_;
dimensionedScalar rho2_;
dimensionedScalar rho3_;

volScalarField rhogp_;
volScalarField mDot_;

// Private Member Functions

//- Calculate and return the laminar viscosity
void calcNu();
void calcrhogp();
void calcmDot();


public:

// Constructors

//- Construct from components
incompressibleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
);


//- Destructor
~incompressibleThreePhaseMixture()
{}


// Member Functions

const word phase1Name() const
{
return phase1Name_;
}

const word phase2Name() const
{
return phase2Name_;
}

const word phase3Name() const
{
return phase3Name_;
}

const volScalarField& alpha1() const
{
return alpha1_;
}

volScalarField& alpha1()
{
return alpha1_;
}

const volScalarField& alpha2() const
{
return alpha2_;
}

volScalarField& alpha2()
{
return alpha2_;
}

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

volScalarField& alpha3()
{
return alpha3_;
}

//- Return const-access to phase1 density
const dimensionedScalar& rho1() const
{
return rho1_;
}

//- Return const-access to phase2 density
const dimensionedScalar& rho2() const
{
return rho2_;
};

//- Return const-access to phase3 density
const dimensionedScalar& rho3() const
{
return rho3_;
};

//- Return the velocity
const volVectorField& U() const
{
return U_;
}

//- Return the flux
const surfaceScalarField& phi() const
{
return phi_;
}

//- Return const-access to phase1 viscosityModel
const viscosityModel& nuModel1() const
{
return nuModel1_();
}

//- Return const-access to phase2 viscosityModel
const viscosityModel& nuModel2() const
{
return nuModel2_();
}

//- Return const-access to phase3 viscosityModel
const viscosityModel& nuModel3() const
{
return nuModel3_();
}

//- Return the dynamic laminar viscosity
tmp<volScalarField> mu() const;

//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> muf() const;

//- Return the kinematic laminar viscosity
tmp<volScalarField> nu() const
{
return nu_;
}

//- Return the laminar viscosity for patch
tmp<scalarField> nu(const label patchi) const
{
return nu_.boundaryField()[patchi];
}

//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> nuf() const;

//- Correct the laminar viscosity
void correct()
{
calcNu();
calcrhogp();
calcmDot();
}

tmp<volScalarField> rhogp()
{
return rhogp_;
}

tmp<volScalarField> mDot()
{
return mDot_;
}

//- Read base transportProperties dictionary
bool read();
};


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

} // End namespace Foam

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

#endif

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

alphaEqn.H:

{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");

tmp<volScalarField> rhogp = mixture->rhogp();
tmp<volScalarField> mDot = mixture->mDot();

mDot = (Dc23 + Dc32)*rhogp*mDot;


surfaceScalarField phir
(
IOobject
(
"phir",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
);

for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
// Create the limiter to be used for all phase-fractions
scalarField allLambda(mesh.nFaces(), 1.0);

// Split the limiter into a surfaceScalarField
slicedSurfaceScalarField lambda
(
IOobject
(
"lambda",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimless,
allLambda,
false // Use slices for the couples
);


// Create the complete convection flux for alpha1
surfaceScalarField alphaPhi1
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha3, alpharScheme),
alpha1,
alpharScheme
)
);

// Create the bounded (upwind) flux for alpha1
surfaceScalarField alphaPhi1BD
(
upwind<scalar>(mesh, phi).flux(alpha1)
);

// Calculate the flux correction for alpha1
alphaPhi1 -= alphaPhi1BD;

// Calculate the limiter for alpha1
if (LTS)
{
const volScalarField& rDeltaT =
fv::localEulerDdt::localRDeltaT(mesh);

MULES::limiter
(
allLambda,
rDeltaT,
geometricOneField(),
alpha1,
alphaPhi1BD,
alphaPhi1,
zeroField(),
zeroField(),
1,
0
);
}
else
{
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha1,
alphaPhi1BD,
alphaPhi1,
zeroField(),
zeroField(),
1,
0
);
}

// Create the complete flux for alpha2
surfaceScalarField alphaPhi2
(
fvc::flux
(
phi,
alpha2,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(phir, alpha1, alpharScheme),
alpha2,
alpharScheme
)
);

// Create the bounded (upwind) flux for alpha2
surfaceScalarField alphaPhi2BD
(
upwind<scalar>(mesh, phi).flux(alpha2)
);

// Calculate the flux correction for alpha2
alphaPhi2 -= alphaPhi2BD;

// Further limit the limiter for alpha2
if (LTS)
{
const volScalarField& rDeltaT =
fv::localEulerDdt::localRDeltaT(mesh);

MULES::limiter
(
allLambda,
rDeltaT,
geometricOneField(),
alpha2,
alphaPhi2BD,
alphaPhi2,
zeroField(),
zeroField(),
1,
0
);
}
else
{
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha2,
alphaPhi2BD,
alphaPhi2,
zeroField(),
zeroField(),
1,
0
);
}

// Construct the limited fluxes
alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;


// Solve for alpha1
solve(
fvm::ddt(alpha1)
+ fvc::div(alphaPhi1)
==
- fvm::Sp(mDot, recRho1)
);


// Add the diffusive flux for alpha3->alpha2
alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(al pha1);

// Solve for alpha2
fvScalarMatrix alpha2Eqn
(
fvm::ddt(alpha2)
+ fvc::div(alphaPhi2)
- fvm::laplacian(Dc23 + Dc32, alpha2)
==
+ fvm::Sp(mDot, recRho2)
);
alpha2Eqn.solve();

// Construct the complete mass flux
rhoPhi =
alphaPhi1*(rho1 - rho3)
+ (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
+ phi*rho3;

alpha3 = 1.0 - alpha1 - alpha2;
}

Info<< "Air phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;

Info<< "Liquid phase volume fraction = "
<< alpha2.weightedAverage(mesh.V()).value()
<< " Min(" << alpha2.name() << ") = " << min(alpha2).value()
<< " Max(" << alpha2.name() << ") = " << max(alpha2).value()
<< endl;
}

Kindly tell me how to proceed solving this problem. Thanks in advance.

[nw_ds]

Hello,
I am trying to do the same , so did you have any luck!?

[Pavithra]

Hi,

Any update on this ? I am stuck at this problem.

Thank You.