[manalis]

Hello FOAMers,

I have the following problem when trying to run a case with "buoyantBoussinesqPimpleFoam" with the new "fvOptions" in OF 2.4 called "solidificationMeltingSource" for simulating phase change (melting/solidification). I used a typical setup for this solver and also included the respective "fvOptions" file in the "system" folder. However, I get the following error:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Reading thermophysical properties

Reading field T

Reading field p_rgh

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Creating turbulence model

Selecting RAS turbulence model laminar
Reading field alphat

Calculating field g.h

Radiation model not active: radiationProperties not found
Selecting radiationModel none
Creating finite volume options from "system/fvOptions"

Selecting finite volume options model type solidificationMeltingSource
    Source: solidificationMeltingSource1
    - applying source for all time
    - selecting cells using cellZone PCM
    - selected 10000 cell(s) with volume 0.0001

Courant Number mean: 0 max: 0

PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Time = 1e-05

Courant Number mean: 0 max: 0
deltaT = 1.2e-05
PIMPLE: iteration 1


--> FOAM FATAL ERROR: 

    request for basicThermo thermophysicalProperties from objectRegistry region0 failed
    available objects of type basicThermo are
0()

    From function objectRegistry::lookupObject<Type>(const word&) const
    in file /home/openfoam/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 198.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::basicThermo const& Foam::objectRegistry::lookupObject<Foam::basicThermo>(Foam::word const&) const at ??:?
#3  Foam::fv::solidificationMeltingSource::Cp() const at ??:?
#4  Foam::fv::solidificationMeltingSource::addSup(Foam::fvMatrix<Foam::Vector<double> >&, int) at ??:?
#5  ? at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8  ? at ??:?
Aborted (core dumped)

The problem probably lies in line where is stated:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solidificationMeltingSource1
    {
        type            solidificationMeltingSource;
        active          on;
        selectionMode   cellZone;
        cellZone        PCM;

        solidificationMeltingSourceCoeffs
        {
            Tmelt            505.00;
            L                 60000.0;
            thermoMode       thermo;   // The problem is here //
            relax                   0.7;
            beta             2.67e-4;
            Cu              1.0e+05;    // Default value Cu=1.0e+05 //
            q                 1.0e-03;     // Default value q=1.0e-03 //
            rhoRef          7500.0;    // Solid density //
        }
    }

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

I wasn't able to find much info about this, so I tried changing the "thermo" word but without success. Does someone have any idea how I should handle this entry or generally how does this "fvOptions" work? I searched the respective page (http://foam.sourceforge.net/docs/cpp...6.html#details) but I didn't find anything helpful (or at least I didn't understand it!)

Any help would be very appreciated!

Thank you in advance

[alexeym]

Hi,

In addition to thermo there is lookup thermoMode. Since you are using buoyantBoussinesqPimpleFoam there is really no thermo object.

So you put lookup instead of thermo in your solidificationMeltingSourceCoeffs. But since you do not have thermo object, you have to set CpName to CpRef and provide CpRef (specific heat) in solidificationMeltingSourceCoeffs.

Foam::tmp<Foam::volScalarField>
Foam::fv::solidificationMeltingSource::Cp() const
{
    switch (mode_)
    {
         ...
        case mdLookup:
        {
            if (CpName_ == "CpRef")
            {
                scalar CpRef = readScalar(coeffs_.lookup("CpRef"));
                ...
            }
    ...

Cu and q is a coefficients for calculation of permeability using Kozeny-Carman equation.

scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);

Default value for Cu (1e5) is OK, while q is quite large. It is a value to avoid division by zero, in my simulations I usually use 1e-6 or 1e-8.

[shipman]

Hi Alexey,

thank you for information. Nowadays, I would like to use solidificationMeltingSource option inside compressibleInterFoam to model the nitrogen phase change within Naval Nozzle. At the beginning I tried to understand the source code of solidificationMeltingSource option. Since in .H file, it is stated that model is based on following papers:

  1. V.R. Voller and C. Prakash, A fixed grid numerical modelling         methodology for convection-diffusion mushy phase-change problems,         Int. J. Heat Mass Transfer 30(8):17091719, 1987.      

2. C.R. Swaminathan. and V.R. Voller, A general enthalpy model for         modeling solidification processes, Metallurgical Transactions         23B:651664, 1992.

Both of them are using different equations and when I looked at the source code I cant understand which one is used exactly. For example, I didnt see followingKozeny-Carman equation in the papers.

Code:
scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
So could you help me to understand exactly which parts of papers are used in the code?

Also I would like to ask your help to understand following lines of source codes:

In solidificationMeltingSourceTemplates.C

  // 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_));

both equations (if and else) have same explanation. could let me know differences.


thanks in advance.

[alexeym]

Hi,

Do I get you right, you want me to

1. Download articles,
2. Read them,
3. Point out difference between articles and implemetation of solidificationMeltingSource?

Sounds rather funny. Unfortunately I have never seen two articles you are referencing, yet in general Voller in his papers deals with the part of algorithm responsible for updating liquid fraction. The way you calculate matter permeability depends on the nature of solidification process. "Usually" people use Kozeny-Carman equation but it is not mandatory to go this way.

Concerning the second part of your questions, yes, equations have the same explanation since they are (almost) identical. This line

eqn.psi().dimensions() == dimTemperature

checks if source is used in temperature equation, so latent heat of fusion

L*fvc::ddt(rho, alpha1_)

should be divided by Cp. Else branch does not have this division.

[shipman]

Hi Alexey,

No, you misunderstood. Of course I did not mean that. I just thought that you already had a look these papers since you are using this fvOption and solidificationMeltingSource is based on these papers as written in .H file.

By the way thank you very much for informative explanation.

My last question is that I have gas and liquid nitrogen at the inlet of the naval nozzle. They have reaction so fast and inside of naval nozzle solid phase appears with liquid phase. Do you think that compressibleInterFoam with the solidificationMeltingSource can be a good choice?

Thank you in advance.

[alexeym]

Hi,

In fact I do not use solidificationMeltingSource, to answer thread-starter's question I have just looked at the sources. One of the drawbacks of Voller's liquid fraction update method, it is applicable only to pure substances (i.e. you have melting temperature instead of melting range).

Unfortunately description of the problem is rather vague to choose solver. Maybe you could start with a list of phenomena you have to address in your simulation and then choice of the solver becomes obvious? Right now I do not see the reason for messing with VOF.

[shipman]

Hi Alex,

I set my fvOptions in the compressibleInterFoam as follows

solidificationMeltingSource1
{
    type            solidificationMeltingSource;
    active          yes; //on;
	


    solidificationMeltingSourceCoeffs
    {
		selectionMode	all;
		cellZone		hotplate;
        Tmelt		   	63.15; //K
        L				25702; //Latent heat of fusion [J/kg]
		relax			0.8; // relaxation coefficient [0~1]
		thermoMode		lookup; //thermo;
		CpName			CpRef;
		CpRef			2064; 
		beta			0.00753; //thermal expansion coefficient [1/K]
		rhoRef			860.65;	//solid density
		Cu				1.0e5;	//Model coefficient default value
		q				1.0e-6; //to avoid dividing zero        
    }
}

when I run the case it doesnt give any error but it gives following warning :

PIMPLE: iteration 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0 Min(alpha.N2liquid) = 0 Min(alpha.N2gas) = 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0 Min(alpha.N2liquid) = 0 Min(alpha.N2gas) = 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
--> FOAM Warning :
From function void option::checkApplied() const
in file fvOption/fvOption.C at line 120
Source solidificationMeltingSource1 defined for field T but never used
--> FOAM Warning :
From function void option::checkApplied() const
in file fvOption/fvOption.C at line 120
Source solidificationMeltingSource1 defined for field T but never used


You can see also some steps taken from logfile below

Create time

Create mesh for time = 0


PIMPLE: Operating solver in PISO mode

Reading field p_rgh

Reading field U

Reading/calculating face flux field phi

Constructing twoPhaseMixtureThermo

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectFluid;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

createFields.H::Reading field alpha1

createFields.H::Calculating field alpha1

Reading thermophysical properties


Reading g

Reading hRef
Calculating field g.h

Selecting turbulence model type laminar
Creating field kinetic energy K

Creating finite volume options from "constant/fvOptions"

Selecting finite volume options model type solidificationMeltingSource
    Source: solidificationMeltingSource1
    - selecting all cells
    - selected 618000 cell(s) with volume 3.5095642e-08
Courant Number mean: 7.1233916e-08 max: 0.0006097561

Starting time loop

Courant Number mean: 7.1233916e-08 max: 0.0006097561
deltaT = 1.1904762e-09
Time = 1.19047619e-09

PIMPLE: iteration 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 2.5463107e-14, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 1.1932199e-14, No Iterations 2
smoothSolver:  Solving for T, Initial residual = 1, Final residual = 1.4906003e-09, No Iterations 2
min(T) 300
 max(T) = 300.15601
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 3.616368e-10, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 0.00019165183, Final residual = 4.8057813e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 2.7732851e-07, Final residual = 4.7384965e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMGPCG:  Solving for p_rgh, Initial residual = 4.9004527e-10, Final residual = 4.9004527e-10, No Iterations 0
max(U) 10
min(p_rgh) 100000
ExecutionTime = 7.47 s


Courant Number mean: 8.5337227e-08 max: 0.0007281608
deltaT = 1.4115646e-09
Time = 2.602040816e-09

"Red"]PIMPLE: iteration 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
--> FOAM Warning : 
    From function void option::checkApplied() const
    in file fvOption/fvOption.C at line 120
    Source solidificationMeltingSource1 defined for field T but never used
--> FOAM Warning : 
    From function void option::checkApplied() const
    in file fvOption/fvOption.C at line 120
    Source solidificationMeltingSource1 defined for field T but never used

smoothSolver:  Solving for Ux, Initial residual = 0.13588954, Final residual = 1.0310746e-14, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.58552683, Final residual = 6.4076165e-14, No Iterations 2
smoothSolver:  Solving for T, Initial residual = 0.13594691, Final residual = 1.7317232e-10, No Iterations 2
min(T) 300
 max(T) = 300.33976
GAMG:  Solving for p_rgh, Initial residual = 0.17584259, Final residual = 3.9993902e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 0.00013956399, Final residual = 1.8270438e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 2.3811297e-07, Final residual = 1.8059967e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMGPCG:  Solving for p_rgh, Initial residual = 4.4884261e-10, Final residual = 4.4884261e-10, No Iterations 0
max(U) 10
min(p_rgh) 100000
ExecutionTime = 10.53 s


Courant Number mean: 1.0268911e-07 max: 0.0008693571
deltaT = 1.6792752e-09
Time = 4.281315975e-09

PIMPLE: iteration 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
MULES: Solving for alpha.N2liquid
Liquid phase volume fraction = 0  Min(alpha.N2liquid) = 0  Min(alpha.N2gas) = 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Ux, Initial residual = 0.095160283, Final residual = 1.7536634e-14, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.6744119, Final residual = 1.3559478e-13, No Iterations 2
smoothSolver:  Solving for T, Initial residual = 0.078889661, Final residual = 8.5225492e-11, No Iterations 2
min(T) 300
 max(T) = 300.55591
GAMG:  Solving for p_rgh, Initial residual = 0.088569436, Final residual = 1.7222934e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 0.00011342716, Final residual = 1.0712451e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMG:  Solving for p_rgh, Initial residual = 2.2675389e-07, Final residual = 1.0576804e-11, No Iterations 1
max(U) 10
min(p_rgh) 100000
GAMGPCG:  Solving for p_rgh, Initial residual = 4.8130432e-10, Final residual = 4.8130432e-10, No Iterations 0
max(U) 10
min(p_rgh) 100000
ExecutionTime = 13.62 s


Courant Number mean: 1.2547452e-07 max: 0.0010467447
deltaT = 1.9941393e-09
Time = 6.275455225e-09

Do you have any idea how to fix it?

thank you in advance.

Baris

[alexeym]

Hi,

According to the output you have posted, I guess it is your modified version of compressibleInterFoam, since at least neither in OpenFOAM 2.4.x, nor in OpenFOAM 3.0.x the solver does not have fvOptions functionality implemented (or maybe I did not search thorough enough).

To answer your question I need to see the modifications.

[shipman]

Originally Posted by alexeym

Hi,

According to the output you have posted, I guess it is your modified version of compressibleInterFoam, since at least neither in OpenFOAM 2.4.x, nor in OpenFOAM 3.0.x the solver does not have fvOptions functionality implemented (or maybe I did not search thorough enough).

To answer your question I need to see the modifications.

Hi Alex,

In fact, there is no FvOptions Functionality in the compressibleInterFoam, but I added following lines into UEqn:

fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(U)
      ==
       	fvOptions(rho, U) //added to use solidification model
    );

    UEqn.relax();

    fvOptions.constrain(UEqn); //added to use solidification model

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    interface.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );
	
	fvOptions.correct(U); //added to use solidification model

        K = 0.5*magSqr(U);
    }

Then added dependencies into .C file related to fvOptions

#include "fvIOoptionList.H" // added to use fvOptions
#include "createFvOptions.H" // added to use fvOptions

finally added following line into Make/Options related fvOptions:

-I$(LIB_SRC)/fvOptions/lnInclude \
 -lfvOptions \

I dont know why it does not work. Do you have any idea?

Thank you.

Baris

[shipman]

Hi Alex,

Ok I found the problem that I forgot to add the fvOptions for TEqn. then it works for now OK without that warning message.

thank you for replies.

Baris

[pmdelgado2]

@Baris: Could you post the code for TEqn with the fvOptions that enable the solver to work properly?

[shipman]

Hi Paul,

It is totally same with U equation which i post above. Just change U parameter with T.

BR

[alexj]

@pmdelgado2: Dear Paul,

Did you work out how shipman has modified the Teqn.H? I'm working on the same problem and would appreciate help.

Best regards,
Alex

[alexj]

Dear all,

I have also now modified the compressibleInterFoam solver according to this thread. When it is fully working, I will post the solver here. However I am currently stuggling to run a test case, based on the depthCharge2D case.

Introducing a porous area works without a problem, so the fvOptions functionality must work in principle. However when I set my fvOptions to model with solidificationMeltingSource as such

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ice1
    {
        type            solidificationMeltingSource;
        active          yes;
        solidificationMeltingSourceCoeffs
        {
            selectionMode   cellZone;
            cellZone        ice;
            Tmelt           273.15;
            L               334000;
            thermoMode      thermo;
            beta            50e-6;
            rhoRef          800;
        }
    }

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

the solver complains right at the start with:

Creating finite volume options from "constant/fvOptions"

Selecting finite volume options model type solidificationMeltingSource
    Source: ice1
    - selecting cells using cellZone ice
    - selected 6400 cell(s) with volume 0.2


--> FOAM FATAL ERROR: 
Not implemented

    From function twoPhaseMixtureThermo::he() const
    in file twoPhaseMixtureThermo.H at line 132.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::twoPhaseMixtureThermo::he() const at ??:?
#3  Foam::fv::solidificationMeltingSource::solidificationMeltingSource(Foam::word const&, Foam::word const&, Foam::dictionary const&, Foam::fvMesh const&) at ??:?
#4  Foam::fv::option::adddictionaryConstructorToTable<Foam::fv::solidificationMeltingSource>::New(Foam::word const&, Foam::word const&, Foam::dictionary const&, Foam::fvMesh const&) at ??:?
#5  Foam::fv::option::New(Foam::word const&, Foam::dictionary const&, Foam::fvMesh const&) at ??:?
#6  Foam::fv::optionList::reset(Foam::dictionary const&) at ??:?
#7  Foam::fv::optionList::optionList(Foam::fvMesh const&, Foam::dictionary const&) at ??:?
#8  Foam::fv::IOoptionList::IOoptionList(Foam::fvMesh const&) at ??:?
#9  ? at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  ? at ??:?
Aborted

I have no idea why I get the "Not implemented" error. I work in OpenFoam 3.0.0.

Any help at this point is highly appreciated and please do let me know if you need more info from my side to work on that.

Best regards,
Alex

[pmdelgado2]

@Alexj: I suspect the fvOption is inherently assuming that the flow model is purely incompressible. Most likely, you will have to edit the source code pertaining to solidificaionMeltingSource to get it to work right.

I'm curious though... what exactly do you mean when you say "introducing a porous area works without a problem". How did you come to this conclusion? What verification tests could you have possible run if the solver doesn't run at all?

[alexj]

@pmdelgado2: thanks for the hint.

I came to the conclusion that fvOptions must work in principle as I was able to run the case I have set up with an explicit porosity source in my fvOptions:

porosity1
{
    type            explicitPorositySource;
    active          yes;

    explicitPorositySourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        ice;

        type            DarcyForchheimer;

        DarcyForchheimerCoeffs
        {
            d   (7e5 -1000 -1000);
            f   (0 0 0);

            coordinateSystem
            {
                type    cartesian;
                origin  (0 0 0);
                coordinateRotation
                {
                    type    axesRotation;
                    e1      (0.70710678 0.70710678 0);
                    e3      (0 0 1);
                }
            }
        }
    }
}

However that only tells me that the U field modification works. I suspect that porosity source does not act on the T field, which the solidificationMeltingSource does however.

[shipman]

Hi Alex,

Can you post T and U eqns here. I feel that you didnt add fvOptions into TEqn to use it.

Baris

[alexj]

Hi Baris,

thanks for the hint. You are right, I have forgotten to edit the T.eqn. Darn my mistake. Now the solver runs. But of course as there is not really a thermo object in the compressibleInterFoam, similar to buoyantBoussinesqPimpleFoam, I have to set up my fvOptions to use the thermoMode "lookup" as you have been doing all along up in post #7 and which has been pointed out above in post #2 by Alexey as well.

Now my solver works. Here are the modifications for reference.

UEqn.H:

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(U)
     ==
        fvOptions(rho, U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    interface.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

        fvOptions.correct(U);
    
        K = 0.5*magSqr(U);
    }

TEqn.H:

{
    fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T)
      + fvm::div(rhoPhi, T)
      - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
      + (
            fvc::div(fvc::absolute(phi, U), p)
          + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
        )
       *(
           alpha1/mixture.thermo1().Cv()
         + alpha2/mixture.thermo2().Cv()
        )
      ==
        fvOptions(rho, T)
    );

    TEqn.relax();
    fvOptions.constrain(TEqn);
    TEqn.solve();

    mixture.correct();
    
    fvOptions.correct(T);

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

Make/options:

EXE_INC = \
    -ItwoPhaseMixtureThermo \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude

EXE_LIBS = \
    -ltwoPhaseMixtureThermo \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
    -ltwoPhaseMixture \
    -ltwoPhaseProperties \
    -linterfaceProperties \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lfiniteVolume \
    -lmeshTools \
    -lfvOptions

Make/files:

myCompressibleInterFoam.C

EXE = $(FOAM_USER_APPBIN)/myCompressibleInterFoam

as you can see I called my version creatively "myCompressibleInterFoam" and place it in the USER APPBIN.

Here is my working, but not physically very correct fvOptions file:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ice1
    {
        type            solidificationMeltingSource;
        active          yes;

        solidificationMeltingSourceCoeffs
        {
            selectionMode   cellZone;
            cellZone        ice;
            Tmelt           273.15;
            L               334000;
            thermoMode      lookup;
            beta            207e-6;
            rhoRef          913;
            Cu              1.0e+05;
            q               1.0e-06;
            CpName          CpRef;
            CpRef           4195.0;
        }
    }
}

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

I hope that is useful for others who try to go along the same route. Now I will work on an example case and place it on the forum as well (http://www.cfd-online.com/Forums/ope...ingsource.html)

Again, thanks for your input Baris,

best regards,
Alex

[nlinder]

Hi all,

thanks for sharing your thoughts and modifications here.. I would like to know about your impression of the stability of the solver. The cases I set up are running, but from my point of you, it is kind of unstable. It does not diverge, but sometimes the minimum temperature drops unphysically low and the velocities rise.
Do you restrict the solid.melt-source to the fluid where alpha=1 or do you apply it to the whole field?

Best

[farinha]

Originally Posted by alexeym

Hi,

In addition to thermo there is lookup thermoMode. Since you are using buoyantBoussinesqPimpleFoam there is really no thermo object.

So you put lookup instead of thermo in your solidificationMeltingSourceCoeffs. But since you do not have thermo object, you have to set CpName to CpRef and provide CpRef (specific heat) in solidificationMeltingSourceCoeffs.

[CODE]
Foam::tmp<Foam::volScalarField>
Foam::fv::solidificationMeltingSource::Cp() const
{
switch (mode_)
{
...
case mdLookup:
{
if (CpName_ == "CpRef")
{
scalar CpRef = readScalar(coeffs_.lookup("CpRef"));
...
}
...

Hi all,

I've been trying to setup this source term for the buoyantBoussinesqPimple Foam solver.
I'm doing exactly what alexeym is saying here, but the fvOption stops the simulation due to the lack of a volScalarField for Cp.

I wonder if editing the createFields.H file to create a field for Cp would do the trick...

Does anyone know if this source term was coded solely for compressible? I mean, is there a way to set it up without editing the solver or the SMS.C file?

Thanks