CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Problem in using new "fvOptions" called "solidificationMeltingSource"

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By alexj
  • 1 Post By alexj

Reply
 
LinkBack Thread Tools Display Modes
Old   October 20, 2015, 13:23
Default Problem in using new "fvOptions" called "solidificationMeltingSource"
  #1
New Member
 
...
Join Date: Jun 2013
Posts: 18
Rep Power: 6
manalis is on a distinguished road
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:
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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:
Code:
/*--------------------------------*- 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
manalis is offline   Reply With Quote

Old   October 20, 2015, 14:45
Default
  #2
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,666
Rep Power: 27
alexeym will become famous soon enoughalexeym will become famous soon enough
Send a message via Skype™ to 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"));
                ...
            }
    ...
Cu and q is a coefficients for calculation of permeability using Kozeny-Carman equation.

Code:
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.
alexeym is offline   Reply With Quote

Old   December 7, 2015, 05:14
Default
  #3
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
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:

Code:
  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

Code:
  // 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.
shipman is offline   Reply With Quote

Old   December 7, 2015, 06:16
Default
  #4
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,666
Rep Power: 27
alexeym will become famous soon enoughalexeym will become famous soon enough
Send a message via Skype™ to 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

Code:
eqn.psi().dimensions() == dimTemperature
checks if source is used in temperature equation, so latent heat of fusion

Code:
L*fvc::ddt(rho, alpha1_)
should be divided by Cp. Else branch does not have this division.
alexeym is offline   Reply With Quote

Old   December 7, 2015, 20:22
Default
  #5
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
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.
shipman is offline   Reply With Quote

Old   December 8, 2015, 04:06
Default
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,666
Rep Power: 27
alexeym will become famous soon enoughalexeym will become famous soon enough
Send a message via Skype™ to 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.
alexeym is offline   Reply With Quote

Old   December 11, 2015, 04:59
Default
  #7
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
Hi Alex,

I set my fvOptions in the compressibleInterFoam as follows

Code:
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
HTML Code:
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
shipman is offline   Reply With Quote

Old   December 11, 2015, 12:25
Default
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,666
Rep Power: 27
alexeym will become famous soon enoughalexeym will become famous soon enough
Send a message via Skype™ to 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.
alexeym is offline   Reply With Quote

Old   December 13, 2015, 20:18
Default
  #9
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
Quote:
Originally Posted by alexeym View Post
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:
Code:
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
Code:
#include "fvIOoptionList.H" // added to use fvOptions
#include "createFvOptions.H" // added to use fvOptions
finally added following line into Make/Options related fvOptions:

Code:
-I$(LIB_SRC)/fvOptions/lnInclude \
 -lfvOptions \
I dont know why it does not work. Do you have any idea?

Thank you.

Baris
shipman is offline   Reply With Quote

Old   December 13, 2015, 21:36
Default
  #10
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
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
shipman is offline   Reply With Quote

Old   January 8, 2016, 13:32
Default
  #11
New Member
 
Paul
Join Date: May 2012
Posts: 23
Rep Power: 7
pmdelgado2 is on a distinguished road
@Baris: Could you post the code for TEqn with the fvOptions that enable the solver to work properly?
pmdelgado2 is offline   Reply With Quote

Old   January 11, 2016, 04:29
Default
  #12
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
Hi Paul,

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

BR
shipman is offline   Reply With Quote

Old   February 22, 2016, 15:19
Default
  #13
New Member
 
Alex Jarosch
Join Date: Dec 2015
Location: Iceland
Posts: 22
Rep Power: 3
alexj is on a distinguished road
@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 is offline   Reply With Quote

Old   February 23, 2016, 12:21
Default
  #14
New Member
 
Alex Jarosch
Join Date: Dec 2015
Location: Iceland
Posts: 22
Rep Power: 3
alexj is on a distinguished road
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
Code:
/*--------------------------------*- 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:
Code:
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
alexj is offline   Reply With Quote

Old   February 23, 2016, 17:00
Default
  #15
New Member
 
Paul
Join Date: May 2012
Posts: 23
Rep Power: 7
pmdelgado2 is on a distinguished road
@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?
pmdelgado2 is offline   Reply With Quote

Old   February 23, 2016, 17:33
Default
  #16
New Member
 
Alex Jarosch
Join Date: Dec 2015
Location: Iceland
Posts: 22
Rep Power: 3
alexj is on a distinguished road
@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:
Code:
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.
Caio Martins likes this.
alexj is offline   Reply With Quote

Old   February 23, 2016, 20:22
Default
  #17
Senior Member
 
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 6
shipman is on a distinguished road
Hi Alex,

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

Baris
shipman is offline   Reply With Quote

Old   February 24, 2016, 06:35
Smile compressibleInterFoam working
  #18
New Member
 
Alex Jarosch
Join Date: Dec 2015
Location: Iceland
Posts: 22
Rep Power: 3
alexj is on a distinguished road
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:
Code:
    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:
Code:
{
    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:
Code:
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:
Code:
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:
Code:
/*--------------------------------*- 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
Xinze likes this.
alexj is offline   Reply With Quote

Old   July 11, 2016, 10:41
Default
  #19
Member
 
Nicklas Linder
Join Date: Jul 2012
Location: Germany
Posts: 35
Rep Power: 7
nlinder is on a distinguished road
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
nlinder is offline   Reply With Quote

Old   May 24, 2017, 11:01
Default SolidificationMeltingSource incompressible
  #20
New Member
 
Farinha
Join Date: Mar 2017
Posts: 1
Rep Power: 0
farinha is on a distinguished road
Quote:
Originally Posted by alexeym View Post
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
farinha is offline   Reply With Quote

Reply

Tags
fvoptions

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
engineFoam new mesh problem ayhan515 OpenFOAM Meshing & Mesh Conversion 5 August 10, 2015 08:45
UDF compiling problem Wouter Fluent UDF and Scheme Programming 6 June 6, 2012 04:43
Gambit - meshing over airfoil wrapping (?) problem JFDC FLUENT 1 July 11, 2011 05:59
natural convection problem for a CHT problem Se-Hee CFX 2 June 10, 2007 06:29
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 19:13


All times are GMT -4. The time now is 23:53.