CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Add an equation to CompressibleInterFoam (https://www.cfd-online.com/Forums/openfoam-programming-development/240330-add-equation-compressibleinterfoam.html)

AntonioMezzacapo December 28, 2021 03:43

Add an equation to CompressibleInterFoam
 
1 Attachment(s)
I am trying to implement "waveTrasmessive.BC" in compressibleInterFoam. I added this code to createFilds.h

Code:

volScalarField DrhoDp
(
IOobject
(
"DrhoDp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
alpha1*psi1 + alpha2*psi2
);

and in alphaEqns.h I added this:

Code:

DrhoDp = psi1*alpha1 + psi2*alpha2;
but after typing "wmake" I'm getting this error:

Does anyone know by which this error depend on? How can I fix it? Thanks! :confused:

Lorenzo210 December 29, 2021 17:13

Hi,
The errors say that all the variables (alpha1 alpha2 psi1 psi2) were not declared before. Try to look if they are defined after DrhoDp.. If so, move the declaration of DrhoDp after them.
Regards,
Lorenzo

AntonioMezzacapo December 30, 2021 06:30

thanks Lorenzo
 
Quote:

Originally Posted by Lorenzo210 (Post 819317)
Hi,
The errors say that all the variables (alpha1 alpha2 psi1 psi2) were not declared before. Try to look if they are defined after DrhoDp.. If so, move the declaration of DrhoDp after them.
Regards,
Lorenzo


Hi Lorenzo, thank you for your response. As suggested in a previous topic, I added the first equation in createfields.h , I got this:

Code:

  Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
    IOobject
    (
        "p_rgh",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);






#include "createPhi.H"

Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
twoPhaseMixtureThermo mixture(mesh);

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

Info<< "Reading thermophysical properties\n" << endl;

volScalarField& rho1 = mixture.thermo1().rho();
volScalarField& rho2 = mixture.thermo2().rho();

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    alpha1*rho1 + alpha2*rho2
);


dimensionedScalar pMin
(
    "pMin",
    dimPressure,
    mixture
);

mesh.setFluxRequired(p_rgh.name());


#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"


// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
    IOobject
    (
        "rhoPhi",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    fvc::interpolate(rho)*phi
);

volScalarField dgdt
(
    pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);

// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, mixture);

volScalarField DrhoDp
(
IOobject
(
"DrhoDp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
alpha1*psi1 + alpha2*psi2
);


// Construct compressible turbulence model
autoPtr<compressible::turbulenceModel> turbulence
(
    compressible::turbulenceModel::New(rho, U, rhoPhi, mixture)
);

Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));

Later, in alphaEqns.H I added the second equation, getting:


Code:

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

    surfaceScalarField phir(phic*interface.nHatf());

    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
    {
        volScalarField::Internal Sp
        (
            IOobject
            (
                "Sp",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
        );

        volScalarField::Internal Su
        (
            IOobject
            (
                "Su",
                runTime.timeName(),
                mesh
            ),
            // Divergence term is handled explicitly to be
            // consistent with the explicit transport solution
            divU*min(alpha1, scalar(1))
        );

        forAll(dgdt, celli)
        {
            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
            {
                Sp[celli] -= dgdt[celli]*alpha1[celli];
                Su[celli] += dgdt[celli]*alpha1[celli];
            }
            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
            {
                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
            }
        }


        surfaceScalarField alphaPhi1
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

        MULES::explicitSolve
        (
            geometricOneField(),
            alpha1,
            phi,
            alphaPhi1,
            Sp,
            Su,
            1,
            0
        );

        surfaceScalarField rho1f(fvc::interpolate(rho1));
        surfaceScalarField rho2f(fvc::interpolate(rho2));
        rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;

        alpha2 = scalar(1) - alpha1;
    DrhoDp = psi1*alpha1 + psi2*alpha2;
    }

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



May you give me some advice? I am still a beginner , and I'm sure many users might find this useful! thank you :)


All times are GMT -4. The time now is 01:10.