[cfd@kgp]

Hi Foam experts,

I was trying to modify the solidificationMeltingSource (Isothermal ohase change source) for the mushy ploblems in OpenFoam version 3.0.

I was successful till defining the liquid fraction using Tsolidus and Tliquidus.

But some how I am unable to get the "Phi field" in order to add the term " fvc::div(phi, alpha1_))" in solidificationMeltingSourceTemplates.C.

\*---------------------------------------------------------------------------*/

#include "fvMatrices.H"
#include "fvcDdt.H"
#include "fvcDiv.H"


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

template<class RhoFieldType>
void Foam::fv::mushysolidificationMeltingSource::apply
(

    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn
)
{
    if (debug)
    {
        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
    }

    const volScalarField Cp(this->Cp());

    update(Cp);

    dimensionedScalar L("L", dimEnergy/dimMass, L_);

    // contributions added to rhs of solver equation
    if (eqn.psi().dimensions() == dimTemperature)
    {
        // isothermal phase change - only include time derivative
         eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change  
      //  eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
    }
    else
    {
        // isothermal phase change - only include time derivative
         eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change  
       // eqn -= L*(fvc::ddt(rho, alpha1_));
    }
}


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

I get following compilation error,

In file included from sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.H:269:0,
                 from sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C:26:
sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C: In member function ‘void Foam::fv::mushysolidificationMeltingSource::apply(const RhoFieldType&, Foam::fvMatrix<double>&)’:
sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C:56:57: error: ‘phi’ was not declared in this scope
          eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change 
                                                         ^
sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSourceTemplates.C:62:54: error: ‘phi’ was not declared in this scope
          eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_)); //mushy phase change 
                                                      ^
sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C: In member function ‘void Foam::fv::mushysolidificationMeltingSource::update(const volScalarField&)’:
sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.C:170:16: warning: unused variable ‘Cpc’ [-Wunused-variable]
         scalar Cpc = Cp[cellI];
                ^
make: *** [Make/linuxGccDPInt32Opt/sources/derived/mushysolidificationMeltingSource/mushysolidificationMeltingSource.o] Error 1

I have included the header file for the divergence term fvcDiv.H. should I include any other header file?

Thanks in advance.
kindly help

[cnsidero]

The flux, phi, is not stored in that class (or parent classes) so you need to get a reference to it through the registry. Look in one of the other fvOptions sources that requires the flux for an example. The mesh is accessible so this should work:

const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

-Chris

[cfd@kgp]

Originally Posted by cnsidero

The flux, phi, is not stored in that class (or parent classes) so you need to get a reference to it through the registry. Look in one of the other fvOptions sources that requires the flux for an example. The mesh is accessible so this should work:

const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

-Chris

Thanks a lot cnsidero!

[janghel]

hi cfd@kgp
were you able to run make the solver for solidificationMelting using fvOptions for mushy zone??

[janghel]

Anyways I got it working for mushy zone ( both linear and schiel law case)

[MauricioReis]

Originally Posted by cnsidero

The flux, phi, is not stored in that class (or parent classes) so you need to get a reference to it through the registry. Look in one of the other fvOptions sources that requires the flux for an example. The mesh is accessible so this should work:

const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

-Chris

Hello,

I'm using OpenFOAM 4.1, I have made the proposed modifications, but the following error occurs during compilation:

error: ‘div’ is not a member of ‘Foam::fvc’
eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));


can anybody help me?

[blebon]

Did you include fvcDiv.H?

My mushyZoneSourceTemplates.C file looks like this and it compiles fine:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvMatrices.H"
#include "fvcDdt.H"
#include "fvcDiv.H"

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

template<class RhoFieldType>
void Foam::fv::mushyZoneSource::apply
(
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn
)
{
    if (debug)
    {
        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
    }

    const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

    const volScalarField Cp(this->Cp());

    update();

    dimensionedScalar L("L", dimEnergy/dimMass, L_);

    // contributions added to rhs of solver equation
    if (eqn.psi().dimensions() == dimTemperature)
    {
        // isothermal phase change - only include time derivative
        eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
    }
    else
    {
        // isothermal phase change - only include time derivative
        eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L*(fvc::ddt(rho, alpha1_));
    }
}


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

[janghel]

Originally Posted by blebon

Did you include fvcDiv.H?

My mushyZoneSourceTemplates.C file looks like this and it compiles fine:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvMatrices.H"
#include "fvcDdt.H"
#include "fvcDiv.H"

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

template<class RhoFieldType>
void Foam::fv::mushyZoneSource::apply
(
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn
)
{
    if (debug)
    {
        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
    }

    const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

    const volScalarField Cp(this->Cp());

    update();

    dimensionedScalar L("L", dimEnergy/dimMass, L_);

    // contributions added to rhs of solver equation
    if (eqn.psi().dimensions() == dimTemperature)
    {
        // isothermal phase change - only include time derivative
        eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
    }
    else
    {
        // isothermal phase change - only include time derivative
        eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L*(fvc::ddt(rho, alpha1_));
    }
}


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

Yes you have to include fvcDiv.H. your mushy template looks good to me. i think it will work.

-Diwakar Janghel

[blebon]

Thanks, I was replying to MauricioReis's error ... but since you are here, did you ever try to use a tabulated fraction liquid:

void Foam::fv::mushyZoneSource::update()
{
    if (curTimeIndex_ == mesh_.time().timeIndex())
    {
        return;
    }

    if (debug)
    {
        Info<< type() << ": " << name_ << " - updating phase indicator" << endl;
    }

    // update old time alpha1 field
    alpha1_.oldTime();

    const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);

    interpolationTable<scalar> fraction_curve("constant/fL");

    forAll(cells_, i)
    {
        label celli = cells_[i];

        scalar Tc = T[celli];
        scalar alpha1New;
        if (Tc < Tsolidus_)
        {
            alpha1New = 0.0;
        } else if (Tc < Tliquidus_) {
            alpha1New = fraction_curve(Tc);
        } else {
            alpha1New = 1.0;
        }

        alpha1_[celli] = max(0, min(alpha1New, 1));
    }

    alpha1_.correctBoundaryConditions();

    curTimeIndex_ = mesh_.time().timeIndex();
}

[janghel]

Originally Posted by blebon

Thanks, I was replying to MauricioReis's error ... but since you are here, did you ever try to use a tabulated fraction liquid:

void Foam::fv::mushyZoneSource::update()
{
    if (curTimeIndex_ == mesh_.time().timeIndex())
    {
        return;
    }

    if (debug)
    {
        Info<< type() << ": " << name_ << " - updating phase indicator" << endl;
    }

    // update old time alpha1 field
    alpha1_.oldTime();

    const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);

    interpolationTable<scalar> fraction_curve("constant/fL");

    forAll(cells_, i)
    {
        label celli = cells_[i];

        scalar Tc = T[celli];
        scalar alpha1New;
        if (Tc < Tsolidus_)
        {
            alpha1New = 0.0;
        } else if (Tc < Tliquidus_) {
            alpha1New = fraction_curve(Tc);
        } else {
            alpha1New = 1.0;
        }

        alpha1_[celli] = max(0, min(alpha1New, 1));
    }

    alpha1_.correctBoundaryConditions();

    curTimeIndex_ = mesh_.time().timeIndex();
}

hii bruno

i did not get your question? are you talking about schiel law?

[blebon]

Yes, I implemented the fraction solid dependence on temperature using a temperature table ... which was itself calculated assuming Scheil solidification. Did you ever try this before? It converges fine ... to the wrong temperature profile.

[janghel]

Originally Posted by blebon

Yes, I implemented the fraction solid dependence on temperature using a temperature table ... which was itself calculated assuming Scheil solidification. Did you ever try this before? It converges fine ... to the wrong temperature profile.

Can you share the formula which you are using to relate solid fraction with temperature..?

[blebon]

Not a formula, but a table:

(
    (  0   0)
    (873   0)
    (883.3 0.05)
    (884.7 0.1)
    (885.4 0.2)
    (886.8 0.3)
    (887.5 0.4)
    (891.6 0.5)
    (900.5 0.6)
    (907.4 0.7)
    (912.2 0.8)
    (913.6 0.9)
    (915.0 1.0)
   (9999   1.0)
)

[janghel]

Originally Posted by blebon

Not a formula, but a table:

(
    (  0   0)
    (873   0)
    (883.3 0.05)
    (884.7 0.1)
    (885.4 0.2)
    (886.8 0.3)
    (887.5 0.4)
    (891.6 0.5)
    (900.5 0.6)
    (907.4 0.7)
    (912.2 0.8)
    (913.6 0.9)
    (915.0 1.0)
   (9999   1.0)
)

Sorry i did not try this. i just used the formula.

[MauricioReis]

Originally Posted by blebon

Did you include fvcDiv.H?

My mushyZoneSourceTemplates.C file looks like this and it compiles fine:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvMatrices.H"
#include "fvcDdt.H"
#include "fvcDiv.H"

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

template<class RhoFieldType>
void Foam::fv::mushyZoneSource::apply
(
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn
)
{
    if (debug)
    {
        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
    }

    const surfaceScalarField& phi =
        mesh_.lookupObject<surfaceScalarField>(phiName_);

    const volScalarField Cp(this->Cp());

    update();

    dimensionedScalar L("L", dimEnergy/dimMass, L_);

    // contributions added to rhs of solver equation
    if (eqn.psi().dimensions() == dimTemperature)
    {
        // isothermal phase change - only include time derivative
        eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
    }
    else
    {
        // isothermal phase change - only include time derivative
        eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
        // eqn -= L*(fvc::ddt(rho, alpha1_));
    }
}


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

Thank you.

[Bloerb]

Can someone explain to me why alpha is updated in the following fashion in the original solidificationMeltingSource? I understand that C*(Tc - Tmelt_) is basically the energy above the melt temperature. We devide by L to get a fraction....why is this needed though? Why not simple if T<Tmelt alpha=0 if T>Tmelt alpha=1...

        scalar alpha1New = alpha1_[celli] + relax_*Cpc*(Tc - Tmelt_)/L_;

And how would you update it for mushy problems? Something like this?

    alpha1 = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5);

And again this would set it directly unlike the above method

[cfd@kgp]

Originally Posted by Bloerb

Can someone explain to me why alpha is updated in the following fashion in the original solidificationMeltingSource? I understand that C*(Tc - Tmelt_) is basically the energy above the melt temperature. We devide by L to get a fraction....why is this needed though? Why not simple if T<Tmelt alpha=0 if T>Tmelt alpha=1...

        scalar alpha1New = alpha1_[celli] + relax_*Cpc*(Tc - Tmelt_)/L_;

And how would you update it for mushy problems? Something like this?

    alpha1 = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5);

And again this would set it directly unlike the above method

Bloerb/Stefan I am curious how did you implemented it using functions, please let us know your progress. meanwhile, I had send you a messege too, but not sure about its delivery so thought to get your response here.

Thanks in advance

[blebon]

Originally Posted by Bloerb

Can someone explain to me why alpha is updated in the following fashion in the original solidificationMeltingSource? I understand that C*(Tc - Tmelt_) is basically the energy above the melt temperature. We devide by L to get a fraction....why is this needed though? Why not simple if T<Tmelt alpha=0 if T>Tmelt alpha=1...

        scalar alpha1New = alpha1_[celli] + relax_*Cpc*(Tc - Tmelt_)/L_;

And how would you update it for mushy problems? Something like this?

    alpha1 = 0.5*Foam::erf(4.0*(T-Tmelt)/(Tl-Ts))+scalar(0.5);

And again this would set it directly unlike the above method

If you update melt_fraction using the if statement, energy is not conserved, because solution of the energy balance equation is based on previous values of melt fraction and hence inconsistent with the current temperature. If I am not mistaken, the correction process (i.e. what is implemented in solidificationMeltingSource) is based under the assumption that enthalpies is the same both for solving the energy balance equation and updating the melt fraction. Upon convergence of the iterative process, both energy is conserved and melt_fraction residuals should be very small.

When you have a mushy zone, you need to replace Tmelt_ but the temperature corresponding to alpha1_[celli] (either using a formula -- lever rule, Scheil or whatever --) or a lookup table.

[gtarang]

Hi,
I have made some changes in solidificationMeltingSource and I use it with tinkered buoyantPimpleSolver which I named as meltFoam. I basically solve energy equation before NS. Further, I solve Energy equation 'n' (nEnergyCorrector in PIMPLE dict) times so that energy equation converges reasonably before I move to next PIMPLE corrector. I am attaching both things here, feel free to comment or suggest improvements.

[mikulo]

Originally Posted by gtarang

Hi,
I have made some changes in solidificationMeltingSource and I use it with tinkered buoyantPimpleSolver which I named as meltFoam. I basically solve energy equation before NS. Further, I solve Energy equation 'n' (nEnergyCorrector in PIMPLE dict) times so that energy equation converges reasonably before I move to next PIMPLE corrector. I am attaching both things here, feel free to comment or suggest improvements.

Hi, thank you for your shared case files and solver. May I ask in which solver will it work? It would have been nice if it will work on OF 7. I will try it later..