[liam_montgomery]

Hi Fabian,

Thank you for the effort which you have put into the meltFoam solver. I am interested in using your solver in the development of a thermal storage model.

In your paper you indicate the choice of constants C and b in the porosity function for the meltFoam test case based on the work of Gau and Viskanta. In your own experiment, you follow the work of Shmueli et al. in guide the selection of the large constant C.

I intend to model the location of the phase front in a melting Silicon body. Could you please point me in the right direction regarding the selection of the constants C and b? Brent and Voller et at. (1988) state "Further work needs to be carried out in the area of mushy region phase change to establish guiding principles in assigning appropriate values for both the constants C and b". ...I am hoping there has been some development in this area since 1988!

[fabian_roesler]

Hi Liam

Unfortunately, there are not much more references I can recommend to you. Shmueli et al. performed a parametric study on the large constant C. However, they modeled a VOF gas phase above the melt. This may lead to wrong results. In my present study I found out, that the large constant is not the only parameter to vary. Next to the permeability (which is the large constant C), the viscosity of the melt in the transition region has to be taken into account. Moreover, by using the liquid fraction as porosity function, porosity is changing linear. This is not necessarily the case for all PCM.
As you use a pure PCM (silicon) you should choose a large C for the Darcy source term. You could also use the switch-of technique instead of the Darcy term. Soon I will publish a paper on all this stuff concerning the Darcy type source term. Also keep in mind that my solver is for non-isothermal phase change of mixtures.

Regards

Fabian

[liam_montgomery]

Hi Fabian,

Thanks for getting back to me on this. Your advice has been helpful so far - I'll be sure to share my findings with you. Just to keep you in the loop, I also emailed Vaughan Voller with a similar question and this was his response:

With a pure material where one would physically expect a sharp interface between the solid and liquid the choice of the parameters in the porosity model are not critical--provided they are sufficient to "switch of" the liquid velocity in the computational cell undergoing the phase change

So in a first case I would use the Gallium melting values. But there would no harm in seeing what happens when you adjust these values--I would not expect adjustment up to 1/2 and order to have much impact--but worth and easy to check.

Thanks again for your help.

Liam Montgomery

[fabian_roesler]

Hi Liam

Thanks that you share your knew findings with us. Vaughan Voller is for sure wright when he says, that the C has no big impact when it comes to isothermal melting. A higher C does not change much. Is C to low, you will see your solid PCM sinking downwards. I would suggest a parametric study varying C from a low value lime 100 to say 10e8. This will give you a feeling. Apropos, when you use different densities for solid and liquid, C will become more and more important.

Regards and thanks again for the discussion

Fabian

[einatlev]

I just downloaded the erfMelt piece and the meltFOAM piece. Will try to use this to model lava!
anyone has any kind of documentation for these? no worries if not. Just figured it will save me some time if someone already wrote something.

[einatlev]

I just downloaded the erfMelt piece and the meltFOAM piece. Will try to use this to model lava!
anyone has any kind of documentation for these? no worries if not. Just figured it will save me some time if someone already wrote something.

[fabian_roesler]

Hi

I did not write any documentation for the solver, since it is not 100 percent finished. The erfMeltSolver reduces the non-linearity effects in the energy conservation equation. However, it does not entirely get rid of it. For simple problems, the results converge with a very small error for the energy conservation. With increasing convective transport, the error increases and some iteration to account for the non-linearity have to be performed.

F. Rösler, D. Brüggemann (2011): Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments. Heat and Mass Transfer, Vol. 47 Issue 8 , 1027-1033, DOI: 10.1007/s00231-011-0866-9
http://www.springerlink.com/content/b1tp01k2u7q8j432/

In my previous work, I use a linear liquid fraction function and the linear source method proposed by Voller.

Regards

Fabian

[Tobi]

Originally Posted by Yahoo

Hi Guys
Have you tried to compile erfConvectiveMeltingPimpleFoam.C (Fabian's code for melting) on OpenFOAM 2.2.0? I am getting this compile error:

make: *** No rule to make target `/opt/OpenFOAM-2.2.0/src/finiteVolume/lnInclude/newFvPatchField.C', needed by `erfConvectiveMeltingPimpleFoam.dep'. Stop.

Hi guys,

if I have time I ll rebuild the solver for OpenFOAM-2.2.x
At the Moment I rebuild the flamelet solver and build a transient solver for combustion based on the flamelet thermodynamic model.

After that (and please write me a message if I dont Keep you updated ) I try to rebuild the solver for 2.2.x

Tobi

[ahmmedshakil]

Hi,
I'm working with the meltIcoFoam solver which is basically cloned from the meltFoam. In my cases, I have seen the liquid fraction at the boundary gives some unusual value. Is there anyone who has worked with the same code of meltFoam? It will be great if anyone can upload validated case for me for the meltFoam.

Cheers!

Originally Posted by AnjaMiehe

Hello David,

this I not trivial for me either as the solver was coded in 1.4.1, or at least that is the tag in the test case. Also, the test case proposes icoFoam, a PISO based solver, to solve the test case. I tried to rewrite meltFoam in 2.1.0 using the SIMPLE Algorithm but I found far more similarities in coding with the icoFoam solver in 2.1.0.

Therefore, and as my try for the meltSimpleFoam does compile but not solve the test case, here is my idea for the meltIcoFoam in 2.1.0 . For me, it compiles only with warnings and the rough test case given earlier in this thread works alright. I can't guarantee for anything.


I hope this is a help for you,
Regards, Anja

[AnjaMiehe]

Hello Mohammad,
could you post your solver and test case? Without more information, helping you is a difficult task.
Regards, Anja

[ahmmedshakil]

The solver I'm using is exactly the meltIcoFoam "http://www.cfd-online.com/Forums/openfoam/93620-melting-problem.html", and the test case is as attached.

[AnjaMiehe]

Hello Mohammad,
I do not see anything in your set up for the test case that might force the solver to crash. As I do not know, how you changed the solver to work with a newer version of OpenFOAM than 1.4 (which I guess is the one posted in #2 of this thread) I can only guess.

• Does your problem occur in single or parallel computation?
• To prepare for parallel processing, you might consider to change the max command in the TEqn to Foam::gMax to account for all processors.
• In the UEqn: the additional buoyancy term is now situated in the pEqn.
• To bound the lf, in the sources file you could try lf=max(0,min(1,DH/L)) instead of only lf=DH/L

I hope, that will help you.
Regards, Anja

[fabian_roesler]

Hi

I did not see any boundary condition for the liquid fraction in your test case. It should be zeroGradient everywhere except for inlets. Check the erfMeltFoamSolver in the thread. You do not necessarily need the boundary condition as the liquid fraction (alpha) is a function of temperature but then it has to be handled a little bit different in the solver code. All fields depending on the liquid fraction, so basically all thermo physical properties, should use the boundary condition of the liquid fraction.

      volScalarField xyz
      (
          IOobject
          (
              "xyz",
              runTime.timeName(),
              mesh,
              IOobject::NO_READ,
              IOobject::AUTO_WRITE
          ),
          alpha*x*y*z,
          alpha.boundaryField().types()
      );

Regards

Fabian

[AnjaMiehe]

Hello Fabian,
thanks for the reply. Again, I learned something new.
Regards, Anja

[ahmmedshakil]

Hi,
Has anyone already implemented the linear source method proposed by Voller ? Will it be possible to post code for that (Just the iteration loop between temperature and liquid fraction)?

Cheers
shakil

[AnjaMiehe]

Hello Shakil,

here is the frame of the way I do it.
For the hEqn.h file, it is:

    int iter=0;
    double diff=0.0;
      
  do
  {
    
    fl.storePrevIter();
    
    iter++;
    
    fvScalarMatrix hEqn
    (
      ...
    );
    
    hEqn.relax();
    hEqn.solve();

    fvScalarMatrix flEqn
    (
      or what way you want to calculate fl out of T
    );

    flEqn.relax(); OR fl.relax(); depending on the way of calculation as well
    flEqn.solve(); // if flEqn in a fvScalarMatrix is used

    diff = Foam::gMax(mag(fl.internalField()-fl.prevIter().internalField()));

  }while ((iter < iterMax) && (diff1 > convergence));

gMax is a version of seeking out the maximum taking into account the parallel distribution if used.
The additional solving parameters iterMax and convergence are declared as follows in a file like readAdditional.h.

dictionary flUpdate = mesh.solutionDict().subDict("flUpdate");   
    scalar iterMax = flUpdate.lookupOrDefault("iterMax", 100);
    doubleScalar convergence = flUpdate.lookupOrDefault("convergence", 1e-6);

For every test case, your fvSolution file is expanded with something like:

flUpdate
{
    iterMax 150;
    convergence 1e-8;
}

I hope, that helps you, it took me quite some time, to make it up myself.

Have a nice day,
Regards,
Anja

[liam_montgomery]

Hi Fabian,

I'm back onto using your solver to model a melting problem. Essentially I am trying to model the melting of a converging prism, with a heat flux applied to the top
(eg: http://imgur.com/RPsaWHa)

The melting occurs as expected in the x-z plane with convection visibly effecting the melting front as shown here: http://imgur.com/iEl6kfu

However, since your solver only appears to work in 2D (one of the assumptions in shell-and-tube paper), the melting is not correctly modelled in the x-y plane: http://imgur.com/9J6bAbM

How much effort would be involved in extending the solver to 3D if at all possible? Would I be better off using 2D results and trying to extrapolate to 3D?

Thanks again for your help.

Kind regards,

Liam Montgomery

[fabian_roesler]

Hi

the solver is fully capable of solving 3D cases . The assumption in the paper is only stated to reduce cell number and so the calculation time. However, I have to state again, that the energy conservation in my error function solver is only valid for small time steps as one would actually need an iterative correction procedure like Voller is using it.

Regards

Fabian

[liam_montgomery]

Ok, thanks - that's great news. Do you have any suggestions as to why convection is not appearing the same in both the x-z and x-y planes as shown in the images?

[ThomasV]

Hi everyone!

I'm going to enter the fray of writing a melting/solidification solver too...
That's why I'd like to ask for some heads-up on the solver for OF2.1 which Fabian kindly shared with us. I understand that it's based upon the enthalpy-porosity technique by Voller et al and the fluid flow is modelled using a PIMPLE solution. I'd like to ask though what general idea the simulation of the flow is based upon. Is it a VOF approach, an Euler-Euler approach or even something different?