CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   melting problem: looking for appropriate solvers (https://www.cfd-online.com/Forums/openfoam-solving/93620-melting-problem-looking-appropriate-solvers.html)

kuechenrole May 8, 2015 05:12

casting processes
 
2 Attachment(s)
Good morning SSSS,

I'll try to realize all your instructions over the day. The greates challenge seems to be the new convMeltFoam approach.

Attached you'll find the hole code. I also extended the immiscibleIncompressibleTwoPhaseMixture and incompressibleTwoPhaseMixture in order to store the thermophysical properties. You'll find the new transportModels attached as well.

Thanks,

ole

kuechenrole May 8, 2015 08:27

simple test case for casting processes
 
1 Attachment(s)
Here is a simple test case for my solver attempt. It simulates a simple cube which is open to the atmosphere at the top and has an inlet at the bottom.

kuechenrole May 8, 2015 08:44

A good introduction to the ideas, discretisation and solving strategies of interFoam seems to be the following pdf:
http://infofich.unl.edu.ar/upload/3b...7523c8ea52.pdf

kuechenrole May 9, 2015 04:10

casting processes
 
3 Attachment(s)
Here is a revised version of my 3 phase attempt. Following the hints of SSSS the main changes are:
  • reduced size, 2D domain inkl. symmetryPlane
  • rho and cp is equal for liquid and solid
  • rho becomes a volScalarField featuring the Boussinesq approach
  • added TEqn.relax()

The TEqn still doesn't work. Just cooling the cube results in higher temperatures than initiated and filling the isothermal cube results in floating Point error after raising and dropping T and liquid fraction.

Two ideas are still open:
  • Do anybody know an alternative approach to convMeltFoam, which doesn't use the convective term fvc::div( *** , fl )?
  • How can terms in the TEqn, which depend on the melt fraction (alpha1 = liquid + solid) be limited, other than just multiplying it with alpha1?

Thanks,

ole

kuechenrole May 11, 2015 11:16

2 Attachment(s)
New try with the approach from file:///home/konsole/Downloads/PhD_Thesis_ZSSALDI.pdf

kuechenrole May 22, 2015 03:41

Hello there,

does anybody know a common test case, which consider a phase change due to melting or solidification and the propagation of the PCM (phaseChangeMaterial). I want to benchmark my solver.

The frame is a casting process.

ahmmedshakil May 28, 2015 11:41

Quote:

Originally Posted by kuechenrole (Post 547276)
Hello there,

does anybody know a common test case, which consider a phase change due to melting or solidification and the propagation of the PCM (phaseChangeMaterial). I want to benchmark my solver.

The frame is a casting process.

Hi Ole,
I did some validations for some extreme conditions of melting/solidification here: http://people.eng.unimelb.edu.au/ima...ngs/19/142.pdf

Cheers
shakil

high_er July 16, 2015 20:53

Species conservation equation in PimpleFoam
 
1 Attachment(s)
Hi there, I am a new former, which is described in attachment is species conservation equation . Now here is my problem, I don't know how to implement the species conservation equation in PimpleFoam accurately, if there were someone give me a guidance or show me how to do that, I will be very appreciate!Attachment 40892

manalis September 2, 2015 11:39

Regarding different versions of "convMeltFoam" solver
 
Hello phase-change foamers,

I downloaded the Anja's modified version of Fabian's "convMeltFoam" solver for OF2.3 (post #119) and also the latest (if I am not wrong) parallel version from Fabian's solver posted by himself (post #139). I tried the first one (let's call it "Anja's modified version" for clarity) and it runs normally but when I tried the other in serial case (let's call it "Fabian's latest parallel code") I get an error in Transport properties. Specifically, it asks for "mu0" and "muk" and from what I've seen they are declared this way in the respective "ReadTransportProperties.H"; when I put some values in the dictionary I receive an error about dimensions. I suppose that "mu" stands for dynamic viscosity with respective units [Pa*s]-->[kg/(m*s)]. On the contrary, "Anja's modified version" just asks for "nu" (kinematic viscosity) and the case runs without any problem (in OF 2.3.0). What could be the problem?

Thank you in advance!

fabian_roesler September 3, 2015 03:18

Hi manalis

It's a long time ago since I posted this solver here. However, I remember that the viscosity was implemented temperature dependant.

Code:

mu = mu0*(scalar(1)-max(muk*(T-Tl),scalar(0)));
So mu0 has the same units as dynamic viscosity [kg/(m*s)] and muk has reciprocal temperature units [1/K].

Cheers

Fabian

manalis September 3, 2015 06:26

Hi Fabian and thanks for the immediate response!

After your suggestion and some minor changes in the "system" folder dictionaries I made it work! Up to now I have been using Fluent for my melting/solidification cases and very recently (actually last week) I decided to try the approach in OpenFOAM. So I am really new to this. Could you explain a little bit the expression for "mu" that you posted? What is the meaning of "scalar(0)", "scalar(1)" and "muk" in your code regarding the way that viscosity is changing with temperature? Have you tried to implement different properties for the 2 different phases with a high viscosity value for the solid phase?

Best regards

manalis September 10, 2015 09:34

convMeltFoam
 
Any thoughts/suggestions regarding viscosity equation? Apart from that, I would be very interested to know if someone has dealt with contact melting (e.g. melting inside a closed capsule) with this specific solver.

Best regards

fabian_roesler October 7, 2015 03:53

solidificationMeltingSource
 
Hi folks. During a small solver research within OpenFOAM 2.4.x I found a new fvOption for solidification and melting:
Quote:

solidificationMeltingSource
Basically it uses the same approach as most solvers in this thread - the enthalpy porosity method by Voller and Prakash.
So just in case somebody is new to this thread or it's topic. This is one way to go besides the proposed solvers in this enlightening thread.

Cheers

Fabian

pmdelgado2 January 12, 2016 16:13

Ole,
interTempFoam is created for what version of OF? The case files seem to imply version 2.3.0, but when I compile it, I get errors stating that the class 'immiscibleIncompressibleTwoPhaseMixture' (in alphaCourantNo.H) lacks a member named 'nearInterface', while UEqn.H and pEqn.H lack members named surfaceTensionForce. Do these members exist in some other version of OF's immiscibleIncompressibleTwoPhaseMixture?

pmdelgado2 January 14, 2016 11:00

Setting up melting problem with interTempFoam for 3 phases (ice, water, and air)
 
Quote:

Originally Posted by kuechenrole (Post 545883)

Hi Ole,
I successfully compiled your solver based on Z. Saldi's thesis (see post #205). I'm trying to setup a problem with a block of ice immersed in heated air. I'm trying to use the 'setFieldsDict' to set the alpha1=0 in the air region and alpha1=1 in the water region.

Within the water region, I want to make part of it solid (ice) and the other part fluid (water). Should I set alphas=0 for the ice and alphas=1 for the liquid water region (only in the water region)? Or should I set alphas=0 for the ice region and alphas=1 everywhere else?

Please let me know how best to setup this problem with your solver.
Thanks,
Paul

pmdelgado2 January 14, 2016 13:45

Hi Fabian,

Could you post your modified convMeltFoam solver that incorporates the third phase (air), as suggested by your thesis?

Quote:

Originally Posted by fabian_roesler (Post 513359)
Hi Folks

My PhD thesis on modeling and simulation of phase change processes in macro encapsulated latent heat thermal energy storage is now published. It is written in German and was published by the Logos Verlag Berlin.

Modellierung und Simulation der Phasenwechselvorgänge in makroverkapselten latenten thermischen Speichern
Thermodynamik: Energie - Umwelt - Technik, Bd. 24
Fabian Rösler
ISBN 978-3-8325-3787-6

http://www.logos-verlag.de/cgi-bin/e...87&lng=deu&id=

The convMeltFoam solver as well as the solver which takes unconstrained close contact melting and an additional gas phase into account are explained in detail. Experimental validation in an rectangular cavity is performed.
Hope you like it.

Cheers

Fabian


alexj February 19, 2016 09:39

convMeltFoam OpenFoam-3.0.0
 
Dear phase-change Foamers,

did anyone convert the convMeltFoam solver from the post http://www.cfd-online.com/Forums/ope...tml#post484860 above to OpenFoam-3.0.0??

Cheers,
Alex

alexj February 25, 2016 08:59

Quote:

Originally Posted by fabian_roesler (Post 566920)
Hi folks. During a small solver research within OpenFOAM 2.4.x I found a new fvOption for solidification and melting:

Basically it uses the same approach as most solvers in this thread - the enthalpy porosity method by Voller and Prakash.
So just in case somebody is new to this thread or it's topic. This is one way to go besides the proposed solvers in this enlightening thread.

Cheers

Fabian

I have uploaded a functional demo case of gallium melting here on the forum which might be helpful for people starting on phase change modeling using the fvOption "solidificationMeltingSource"

http://www.cfd-online.com/Forums/ope...tml#post586823

Best regards,
Alex

tetra-eder March 16, 2016 06:12

Melting in a body-fixed reference frame
 
Hi,
first I would like to thank for this very helpful thread.

I want to simulate a melting problem and I just thought that someone of you might help me.
A heat source is embedded into a phase change material and the heat source moves with a constant and given velocity. To solve for the phase change with convection, I use a solver based on Fabian Rösler's solver which was posted in this thread. I would like to use a fixed grid with a heat source fixed reference frame, so that the phase change material moves relative to the heat source. When using just an inlet and outlet condition for the velocity, the velocity changes within the solid phase, which is unphysical. I think the reason for this is mass conservation. Another possibility would be to add the velocity in the darcy term but this seems also not to work. Do you have any ideas how i can treat this relative motion phase change problem?

Thanks in advance!

alexj March 16, 2016 08:41

@tetra-eder. Unfortunately I have not worked much with moving reference frame simulations yet in OpenFOAM so I don't have a answer to your question readily available for you.

Maybe if you would detail a bit more the simulation you want to achieve it will be easier for the forum to answer you.

If I understand you correctly, your heat source is moving in a static fluid? So you should be able to simulate that with an inlet, outlet boundary condition and let the fluid move instead of the heat source. However where is the solid phase in in your simulation? The term "heat source" implies that you want to melt something right? Or do you want to freeze the liquid onto the moving heat sink?

As you can see, a bit more detail might help.

One option which comes to mind is that you can also use a cyclic boundary condition for the inlet and outlet and use the fvOption "momentumSource" with a "meanVelocityForce" type to prescribe a fixed flow rate in your simulation for the fluid. However if you do that and keep the temperature field also cyclic, then you will effectively heat the fluid with your heat source as no thermal energy will escape the system.

Hope this helps somewhat.

Cheers,
Alex

tetra-eder March 16, 2016 10:31

Thanks for the answer Alex. You are right, the description was not very clear. Here are more details about the problem.
Let's assume a cylinder with a constant temperature (higher than melting temperature) at its boundary. The cylinder is surrounded by a phase change material, which is initially solid. The simulation starts at a time, when there is already some melt around the cylinder, so that motion through the melt is possible (technically I use the setFieldsDict to do this before I start the solver). When the heated cylinder now moves with a small velocity (<= the velocity of the phase interface), it will melt the solid phase change material in the direction of movement.
In a global reference frame the velocity of the solid phase change material is zero, the velocity of the liquid phase change material (melt) is influenced by the boussinesq term + its interaction with the moving cylinder wall and the velocity of the cylinder is given. In a body-fixed reference frame (which I want to use) the velocity of the cylinder is zero, the velocity of the solid phase change material is spatially constant and has the negative velocity that the cylinder would have in a global reference frame and the velocity of the melt is some value + the velocity of the solid phase change material.

I think I've found a solution in this thread (#169) for the problem.

ahmmedshakil March 16, 2016 22:00

Quote:

Originally Posted by tetra-eder (Post 589997)
Hi,
first I would like to thank for this very helpful thread.

I want to simulate a melting problem and I just thought that someone of you might help me.
A heat source is embedded into a phase change material and the heat source moves with a constant and given velocity. To solve for the phase change with convection, I use a solver based on Fabian Rösler's solver which was posted in this thread. I would like to use a fixed grid with a heat source fixed reference frame, so that the phase change material moves relative to the heat source. When using just an inlet and outlet condition for the velocity, the velocity changes within the solid phase, which is unphysical. I think the reason for this is mass conservation. Another possibility would be to add the velocity in the darcy term but this seems also not to work. Do you have any ideas how i can treat this relative motion phase change problem?

Thanks in advance!

Hi Tetra-eder,

Can you please post the governing equations ? I guess it should not be that tough. I have solved a similar problem, but not considering the melting. In my melting problem, I considered the fixed heat source, which was due to the problem I was solving.

Cheers,
shakil

tetra-eder March 17, 2016 03:03

Thanks for the answer, shakil. I also think that the problem is very similar. Could explain how you define the boundary conditions for U and p_rgh (or p)? And if you also used the solid velocity in the momentum equation, how you define this velocity?

The governing equations I use are similar or even the same that are already posted in this thread. Therefore I will not describe the variables, if its okay for you.

Conservation of energy:
\frac{\partial \left( \rho c T \right)}{\partial t}+\nabla\cdot\left(  \rho \mathbf{u} c T\right)=\nabla\cdot\left( \lambda\nabla T \right)-L\left( \frac{\partial\left( \rho \gamma_L \right)}{\partial t} +\nabla\cdot\left( \rho \mathbf{u} \gamma_L\right) \right)

Conservation of momentum:
\frac{\partial\left(\rho \mathbf{u}\right)}{\partial t}+\nabla\cdot\left( \rho \mathbf{u}\mathbf{u} \right)=-\nabla p + \nabla\cdot\left( \eta\nabla\mathbf{u}\right)-\rho \mathbf{g}\beta(T-T_0)-C\frac{\left( 1-\gamma_L \right)^2}{\gamma_L^3+\epsilon}\left( \mathbf{u}-\mathbf{u}_{\text{heatsource}} \right)

Conservation of mass:
\frac{\partial\rho}{\partial t}+\nabla\cdot \left( \rho\mathbf{u} \right)=0

This does of course not include a velocity superposition in the liquid phase but in a first step I neglect it due to the fact that the expected velocity is small compared to the convection in the melt.

One idea that came into my mind when I wrote down these equation was that some terms depending on the velocity of the solid phase were canceled out in the original derivation of the mixture formulation. Could this be the problem?

tetra-eder March 17, 2016 05:54

I solved the problem now. There were several reasons. The first was that I was not using the updated solver with u_heatsource. Since I have compiled Openfoam also in debug mode, I have two versions of my solver and only one has been updated during compilation and unfortunately this was not the one set in my environment :o. After this the velocity in the solid was constant as expected but for some reason the liquid fraction didn't move. And the problem was that I used a flag, which was set to false to de- or activate the convection of the liquid fraction in the heat equation. Now everything works as expected but I would like to discuss the boundary conditions with you. What I am using now is

Code:

FoamFile
{
    version    2.0;
    format      ascii;
    class      volScalarField;
    location    "0";
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 0;

boundaryField
{
  inFlow
    {
        type            fixedFluxPressure;
        gradient      uniform 0;
        value          uniform 0;
    }
    outFlow
    {
        type            fixedFluxPressure;
        gradient      uniform 0;
        value          uniform 0;
    }
    sides
    {
        type            fixedFluxPressure;
        gradient      uniform 0;
        value          uniform 0;
    }
    heatSource
    {
        type            fixedFluxPressure;
        gradient        uniform 0;
        value          uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}

Code:

FoamFile
{
    version    2.0;
    format      ascii;
    class      volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField  uniform (0 0.001 0);

boundaryField
{
    inFlow
    {
      type            fixedValue;
      value          uniform (0 0.001 0);
    }
    outFlow
    {
      type            zeroGradient;
    }
    sides
    {
      type            slip;
    }
    heatSource
    {
      type            fixedValue;
      value          uniform (0 0 0);
    }
    defaultFaces
    {
        type            empty;
    }
}

Code:

FoamFile
{
    version    2.0;
    format      ascii;
    class      volVectorField;
    location    "0";
    object      Us;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField  uniform (0 0.001 0);

boundaryField
{
    inFlow
    {
      type            fixedValue;
      value          uniform (0 0.001 0);
    }
    outFlow
    {
      type          fixedValue;
      value        uniform (0 0.001 0);
    }
    sides
    {
      type          fixedValue;
      value        uniform (0 0.001 0);
    }
    heatSource
    {
      type            fixedValue;
      value          uniform (0 0 0);
    }
    defaultFaces
    {
        type            empty;
    }
}

Do you have any improvement suggestions, especially for p_rgh?

han0459 July 14, 2016 06:25

Help
 
Dear everyone,

Since I am a newbie in the openFoam, I could not compile the meltFoam in my openFoam 4. Is anybody has successfully compile the meltFoam in the openFoam 4? If so, please help me how to debug this problem. And is there any other good ideas to simulate the solidification using openFoam? Because I'd like to develop the algorithm of solidification in the future.

Thanks,

Alex

fabian_roesler July 15, 2016 03:57

Hi

Actually you do not have to compile anything, as the whole solver for solidification and melting is there as a new fvOption module since version 3.

Check this out:
http://cpp.openfoam.org/v4/a02447.html#details

Cheers

Fabian

han0459 July 15, 2016 23:48

Hi Fabian,

Thank you so much for your information. I hope I can use this solver successfully. If there is a problem, I will ask for help in the forum again:p.

Yours,

Alex Han

han0459 July 18, 2016 22:00

Sorry, I still have a little problem. There is no tutorial for the solidification, right?

Alex

fabian_roesler July 20, 2016 10:04

The solver doesn't differ between melting and solidification.
When a boundary is colder than the melting point, solidification starts.

Cheers

Fabian

han0459 July 21, 2016 01:13

Hello Fabian,

I am sorry to ask a dull question. You said I do not have to compile for the module. But how can I use the melting/solidification source? Because I have searched for the tutorial or guide in the internet to use melting/solidification FvOption. But unfortunately I failed. Could you please give some suggestions how to employ that module?

Thanks a lot !

Alex Han

kera November 9, 2016 04:30

Hello Foamers,

I have been a silent spectator of this forum for a while and it seems that now is the time for me to ask for questions and help.

I am also thankful to Fabian, Shakil and others for their valuable suggestions which I found very helpful.

Quote:

Originally Posted by kuechenrole (Post 545645)
Here is a revised version of my 3 phase attempt. Following the hints of SSSS the main changes are:
  • reduced size, 2D domain inkl. symmetryPlane
  • rho and cp is equal for liquid and solid
  • rho becomes a volScalarField featuring the Boussinesq approach
  • added TEqn.relax()

The TEqn still doesn't work. Just cooling the cube results in higher temperatures than initiated and filling the isothermal cube results in floating Point error after raising and dropping T and liquid fraction.

Two ideas are still open:
  • Do anybody know an alternative approach to convMeltFoam, which doesn't use the convective term fvc::div( *** , fl )?
  • How can terms in the TEqn, which depend on the melt fraction (alpha1 = liquid + solid) be limited, other than just multiplying it with alpha1?

Thanks,

ole

Hello Ole,

Were you able to solve your problem? I am actually trying to do the same and using the gallium case to bench mark my solver.

If I am not considering the gas phase "alpha2" the solver works just fine. The moment I introduce the gas phase in my case, I face 2 problems:

1) Initially the TEqn takes too many iterations to converge
2) The melt fraction gamma has a very weird profile at the interface (i.e., (interface of alpha1 (solid+melt) and alpha2 (gas phase).

Would greatly appreciate if I could look into your code.

regards,
Ricky

kera November 21, 2016 09:11

Hello again,

I was able to solve the above mentioned problems up to an extent. But are there any test cases I could use to check my solver where all the 3 phases are present (i.e., phase1 air, phase2 (solid+melt)).

regards,
Ricky

akidess November 21, 2016 09:15

Tan, L. H., Leong, S. S., Leonardi, E., Barber, T. J., Jan. 2006. A numerical study of solid-liquid phase change with marangoni effects using a multiphase approach. Progress in Computational Fluid Dynamics, an International Journal 6 (6), 304-313.

kera November 22, 2016 05:30

Thank you, this is what I was looking for.

regards,
Ricky

pmdelgado2 November 22, 2016 17:21

@kera: Can you share what you did? Can you share your solver?

kera November 23, 2016 10:05

Hi Paul,

I am afraid, I am not allowed to share any thing related to this project for the time being, besides I have just tried to implement the algorithm presented by Kidess et al, in interFoam (which was not easy in my case as I am new to FOAM programming but thanks to all the help I could get from this thread) btw most of his presentations are availabe online. If you still need any help, I will try and help you out with my limited knowledge in this field.

Hello everyone.

Currently I am struggling to validate my solver, as the surface tension (\sigma) and its variation with temp (d\sigma/ dT) are merely used as user defined values in my code and I am not sure about the values which I have in hand, whether they are right or not and at the same time I am also stuck with initializing the intial conditions of the domain.

any help is really appreciated.

regards,
Ricky

kera December 5, 2016 14:32

Hello again!

I am not able to work on this project full time and hence not able to update the scenario more often.

I started working with the project from where I left off, I have implemented the velocity equation as below

\frac{\partial\left(\rho \mathbf{u}\right)}{\partial t}+\nabla\cdot\left( \rho \mathbf{u}\mathbf{u} \right)=-\nabla p + \nabla\cdot\left( \mu\nabla\mathbf{u}\right) + \rho(1-\beta(T-T_{ref}))\mathbf{g} + S_{Darcy} + F_{CSF} + F_{marangoni}

where F_{marangoni} = \frac{\partial \sigma}{\partial T} [\nabla T - \vec{n}(\vec{n} \cdot \nabla T)] and \vec{n} I have defined it in the same was as it is define in the curvature term.

I tried to validate the my solver with the above paper suggested by Anton Kidess but failed. I have attached an image of my result and as you can see its miles away from the results obtained by Tan et al and the interface diffusion is overwhelming .

Any help in this topic is really appreciated.

regards,
Ricky

PS: the liquid melt fraction is a linear function.

akidess December 12, 2016 05:51

In the paper you are trying to reproduce, there is metal on the bottom and gas at the top. Hence it's wrong that you initialized gamma to be zero at the right (the solid material).

kera December 12, 2016 10:09

Hello Anton,

Thank you for your reply.

Yes you are right. I didn't had to initialize gamma to be "zero" on the right.

I still think that there is a problem in my implementation of Marangoni force along with melting algorithm in interfoam.

I am currently following the work of Dr. Saldi, and when I tested my solver against his "Two-Phase Marangoni driven flow" (discussed in his thesis) the results were way off, i.e., the surface height at left and right walls were 0.194 and 0.204 respectively, which are still no where near to the experimental solution obtained by Sen and Davis (1982) .

I used the following boundary condition and schemes in my case:

Boundary Conditions:
Temp: Linear temperature distribution throughout the domain with left and right walls at T_{hot} and T_{cold} resp.

Velocity: fixedValue 0

Alpha: with 90deg contact angle at the walls and zeroGradient otherwise.

Schemes
Temp, Velocity: limitedLinear
alpha: vanLeer

I almost forgot to mention, I recalculated the fluid and thermal properties using the dimensionless numbers.

Thanks once agian.

Regards,
Ricky

akidess December 12, 2016 13:19

Quote:

Originally Posted by kera (Post 629292)
I am currently following the work of Dr. Saldi, and when I tested my solver against his "Two-Phase Marangoni driven flow" (discussed in his thesis) the results were way off, i.e., the surface height at left and right walls were 0.194 and 0.204 respectively, which are still no where near to the experimental solution obtained by Sen and Davis (1982) .

Sen and Davis do not present experimental results if I remember correctly. The results in that paper are based on a semi-analytic model. Based on what you put in you will obtain other results. What were the values you expected, for which dimensionless numbers?


All times are GMT -4. The time now is 13:35.