[alexj]

Dear OpenFoamers,

I want to start simulating the melting of a solid wall bounding a flow. Ultimately the flow will be turbulent and transport heat towards the wall, which will result in the phase change. To get me started I have been looking for a tutorial case that uses the fvoption "solidificationMeltingSource" but so far I only found this thread:

http://www.cfd-online.com/Forums/ope...ingsource.html

which modifies the compressibleInterFoam solver to use the "solidificationMeltingSource".

If somebody has a working example case using "solidificationMeltingSource" that would be of great help to get me started. Else I guess I have to start by modifying the compressibleInterFoam case for myself, but I would prefer a working example to start with.

Thanks in advance

Alex

[alexj]

@Baris:

I am new to using the "solidificationMeltingSource" within OpenFoam (not new to OpenFoam as such) and have asked recently in another part of the forum for a workable test case (http://www.cfd-online.com/Forums/ope...ingsource.html).

I was wondering if you would be interested in putting together some information on how to use that feature. So far I have the following list of questions which people like me ask when they get started:

1. Is "compressibleInterFoam" the solver to start with or is there a better introduction solver to the melting concept? Obviously this depends on the application, but I guess a fairly generic assumption would be that people want to simulate solid/fluid mixture and phase changes therein.
2. How is the solid phase represented. The source code of "solidificationMeltingSource" hints towards a porous solid in the fluid. Is this best done with topoSet? Where are examples for doing this?

Any help at this point would be highly appreciated as I have a hard time compiling information to get started.


Thanks in advance,
Alex

[Moderator note: Moved from here: http://www.cfd-online.com/Forums/ope...ingsource.html]

[shipman]

Hi Alex,

Truely speaking, when I decided to use this solidification fvOption, I have also searched a lot to find any previous application or example. However, I could not find. Then, what I did is that as you mentioned above I have modified the compressibleInterFoam and I thought that solid phase is included inside liquid alpha (I have defined one more volume fraction for solid) and when the Temp is lower/upper than melting point set to use this option for both melting/solidification among only solid and liquid phases. Since I used nitrogen as fluid, under some special conditions also there can be seen direct solidification from nitrogen gas to Nitrogen solid which was ignored in my case.

Hope these hints give you some idea.

BR.

Baris

[alexj]

Hi Baris,

thanks for the answer, this helps along the way. I am currently working on implementing a demo case with a buoyantPimpleFoam solver. When I have a working demo case set up, I will post it here along with some info.

In the meantime, everybody is of course welcome and invited to contribute.

Best regards,
Alex

[alexj]

Hi Baris,

would you mind posting the TEqn.H modification you did in compressibleInterFoam to get the fvOptions working? I also want to try your approach but struggle to get the modified compressibleInterFoam working.

Best regards,
Alex

[shipman]

Hi Alex,

I checked the code of the buoyantPimpleFoam Solver of OF v3.0.x. Since it already has been using fvOptions for both U and T equations, you dont need to make any modificiation. You can use directly.

BR

Baris

[alexj]

Hi Baris,

you are right, I can use the buoyantPimpleFoam Solver right away. I am still working on that example case. Thanks for the input.

I will report back as soon as I have a fully working example going.

BR
Alex

[alexj]

Hi phase-change Foamers,

I have completed a demo case for the solidificationMeltingSource function in OpenFOAM 3.0.x. Attached is the case file. I use the buoyantPimpleFoam solver and set up a case where gallium is melted, similar to this post:
http://www.cfd-online.com/Forums/ope...tml#post329995

You can also head over to vimeo and check a video of the simulation:


Hope this is useful for people starting on phase change models.

Best regards,
Alex

[zfaraday]

Hi Alex,

Many thanks for your contribution! I have been thinking of study this fvOption for a while but I haven't had the time nor the need by now... Now it will be easier to start thanks to your test case!

Btw, could you provide an Allrun script in order to ease and automate the process for the beginners? I say so because I saw your case needs to run setFields and it may cause some confusion...

Best regards and many thanks again,

Alex

[alexj]

Hi Alex,

Thanks for the hint on the Allrun script and glad to hear you find the test case useful.

I think it's always the easiest way to understand a setup when one can study a test case which works.

I am out of town the next few weeks and have no access to my computer, but when I return I will update the test case to include a fully functional "Allrun" and "Allclean" script.

Best regards
Alex

[alexj]

Hi all,

I have added the 'Allrun' and 'Allclean' scripts to the gallium melting demo. Also I have removed the currently not needed 'decomposeParDict' and 'setFieldsDict' files. They were left overs from another attempt of setting up the demo.

This demo defines all BCs directly in the 0/ folder and does actually not require to run 'setFields'. However for convenience I have added the run and clean script to be certain people know instantly how to run the demo.

Cheers,
Alex

[neilpaulx]

Hi Alex,

Many thanks for the sample.

I tried running your shared files as is in OF ver 3.x. It runs fine except that it won't generate aMS1_alpha1 in the succeeding time step results. Any idea what went wrong?

Thanks,
Neil

[alexj]

Originally Posted by neilpaulx

Hi Alex,

Many thanks for the sample.

I tried running your shared files as is in OF ver 3.x. It runs fine except that it won't generate aMS1_alpha1 in the succeeding time step results. Any idea what went wrong?

Thanks,
Neil

Hi Neil,

the demo has been tested on OF 3.0.0 and worked fine. It also works perfectly on OF 4.0.0, I just tested it.

Could you maybe give more info on what is wrong in your case. Do you get an error message?

Cheers,
Alex

[cfd@kgp]

Dear Alexj,

Thanks for the case file.

I was interested to calculate global liquid fraction. sum(alpha1*vol)/sum(vol)
I tried to do it with adding a function to controldic

functions
(
GFl
{
type swakExpression;
valueType cellZone;
zoneName solid; // or whatever is your zoneName
accumulations (
min
);
expression "sum(mag(alpha1)*vol())/sum(vol())";
verbose true;
}
);

As "alpha1" is not a IOobject for buoyantPimpleFoam (it is defined in src/....../solidificationMeltingSource....), the swak4Foam utility the function in controlDict gives error.

Any ideas for getting global fluid/solid fraction per time step.

please help

[janghel]

Dear alex,
I am using your solver for metling of gallium but i am getting an error. i have copy pasted the error. can you give me some solutions.
Selecting finite volume options model type solidificationMeltingSource
Source: sMS1
- selecting cells using cellZone solid
- selected 7200 cell(s) with volume 7.2e-06

Starting time loop

Courant Number mean: 0 max: 0
deltaT = 1
Time = 1

PIMPLE: iteration 1


--> FOAM FATAL ERROR:

request for volScalarField Cp from objectRegistry region0 failed
available objects of type volScalarField are

3
(
nu
sMS1:alpha1
p
)


From function const Type& Foam:bjectRegistry::lookupObject(const Foam::word&) const [with Type = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]
in file /home/ubuntu/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 193.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#2 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const& Foam:bjectRegistry::lookupObject<Foam::Geometric Field<double, Foam::fvPatchField, Foam::volMesh> >(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)

[cfd@kgp]

Originally Posted by alexj

Hi phase-change Foamers,

I have completed a demo case for the solidificationMeltingSource function in OpenFOAM 3.0.x. Attached is the case file. I use the buoyantPimpleFoam solver and set up a case where gallium is melted, similar to this post:
http://www.cfd-online.com/Forums/ope...tml#post329995

You can also head over to vimeo and check a video of the simulation:


Hope this is useful for people starting on phase change models.

Best regards,
Alex

Hi Alex and all members on this Thread,

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

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


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

My compilation error is as given below,

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

[cfd@kgp]

Originally Posted by alexj

Hi all,

I have added the 'Allrun' and 'Allclean' scripts to the gallium melting demo. Also I have removed the currently not needed 'decomposeParDict' and 'setFieldsDict' files. They were left overs from another attempt of setting up the demo.

This demo defines all BCs directly in the 0/ folder and does actually not require to run 'setFields'. However for convenience I have added the run and clean script to be certain people know instantly how to run the demo.

Cheers,
Alex

Hi Alex and all phase-change foamers,

The solidification and melting sources (in OF 3.0 and 4.x ) gives accurate transient results. I have checked both the transient meltfront evolution and the time evolution of Nusselt numbers for melting.

However, the definition of fluid fraction do not seems to convince me, especially values of fluid fraction cannot be correlated with the temperature,

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

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

alpha1_[cellI] = max(0, min(alpha1New, 1));
deltaT_[i] = Tc - Tmelt_;
}

can any body explain me the rationale of this fluid fraction defination---> Cpc*(Tc - Tmelt_)/L_;

Any Research paper with this definition?

Thanks in advance!

[gtarang]

Originally Posted by cfd@kgp

Hi Alex and all phase-change foamers,

The solidification and melting sources (in OF 3.0 and 4.x ) gives accurate transient results. I have checked both the transient meltfront evolution and the time evolution of Nusselt numbers for melting.

However, the definition of fluid fraction do not seems to convince me, especially values of fluid fraction cannot be correlated with the temperature,


can any body explain me the rationale of this fluid fraction defination---> Cpc*(Tc - Tmelt_)/L_;

Any Research paper with this definition?

Thanks in advance!

Hi,
I have answered your question melting problem: looking for appropriate solvers
-
Tarang

[gtarang]

Hi,
I have been using solidificationMeltingSource as prescribed by you. Other than simple problems the solver fails and typically gives floating point exception error and runs for CFL no 0.1 or less for fine meshes.
I dug into the code and found that energy equation needs to be solved multiple times to actually correct the step function of liquid fraction as can be seen from the algorithm of Voller and Prakash (1987) http://dx.doi.org/10.1016/0017-9310(87)90317-6. I tried doing it by putting nEnergyCorrectors in PIMPLE dictionary and solving the energy equation n times. While I was debugging the code, I found this in solidificationMeltingSource.C line no 145.

if (curTimeIndex_ == mesh_.time().timeIndex())
    {
        return;
    }

Now according to my understanding, the above code ensures that the liquid fraction is updated only once and found it also. Can anyone suggest what can be the consequences of removing above bit of code.

Thanks
-
Tarang

[rahul62]

Hi,
I am running the case (solidificationMeltingSource) in of4 but i am not getting alpha1 printed. Can anyone help me with this?