solid/liquid phase change solver convMeltFoam
 

Hi

Attached you will find the announced solver for solid/liquid phase change as well as an extended version which takes heat transfer through a fin in to account. Implementation is very simple using a third phase fraction for the fin material.

Cheers

Fabian

Hi Fabian Roesler,
Thanks for sharing the solver. I'm also following Voller Source Based Enthalpy method. Up to now its performing well. I have also made the thermophysical properties temperature dependent (Cp, K). Until now, the conduction case is performing well and solution converging with just two iterations. But, I'm wondering about the conduction-convection phase change problem and for some complex problems. Have you ever tried with the temperature dependent properties?

Hallo Fabian,

vielen Dank fürs uploaden. Zieh mir den Solver direkt rein.

Viele Grüße
Florian

Hello Fabian,

good to see you are back. I only want to correct a tiny thing: the solver on #70 is not mine, but the adapted one of Mohammad. I am about to write my PhD down as well, and will publish my coding as soon as possible. It is astonishing different to yours and similar as well. I do update the fraction, of course ;-) I am using Scheil model and temperature dependent properties, without updating and relaxing that would blow up immediately.

Hi Fabian,

i checked your solver and I have some problems with the TEqn.
In the book of V.R. Voller (http://sabotin.ung.si/~sarler/VOLLER/adheat.pdf) I found the following equation for enthalpie:

H=[beta*cl+(1-beta)*cs]*T+beta*(cs-cl)*Tm+beta*L (Eq1)

here beta is the liquid volume fraction, cl is the heat capacity of the liquid, cs is the heat capacity of the solid,Tm is the melting temperature and L is the latent heat.
Then we have the energy equation (here without Source) as follow:

dH/dt + div(v*H) - div(lambda/rho*gradT)=0 (Eq2)

If I put Eq1 into Eq2 I get the following (with mean heat capacitiy cp=(beta*cl+(1-beta)cs)):

d(cp*T)/dt + div(v*cp*T) + L*d(beta)/dt + L*div(v*beta) - Tm*(cs-cl)*d(beta)/dt + Tm*(cs-cl)*div(v*beta) - div(lambda/rho*grad(T)) = 0
(Eq3)

This equation is very similar to yours except [Tm*(cs-cl)*d(beta)/dt ] and [Tm*(cs-cl)*div(v*beta)]

Now my question :) Why is term 5 in your solver:

Tm*(cs-cl)*d(beta)/dt = Tm*d(cp)/dt

and term 6 in your solver:

Tm*(cs-cl)*div(v*beta) = Tm*div(v*cp)

???
is (cs-cl)*beta=cp

I have the same Problem with h (TEqn.H Z.40). Here:

h=cp*(T-Tmelt)+alpha1*L

why not:

h=cp*T + alpha1*(cs-cl)*Tmelt + alpha1*L

Thanks again for your solver and Frohe Weihnachten

Florian

Hi riesotto,
Please see the post #75, and I guess you could have idea about that.

Hi ahmmedshakil,

thx for your fast reply.
I'm a little bit confused about post #75, but I will give it a shot again.

I know from the paper of Fabian, that there is a need to handle the liquid fraction with a "continuous liquid fraction"-function.

I will think about all these things and I will read Voller again.

kind regards
Florian

Hi,

until now i don't understand why:

Tm*(cs-cl)*d(beta)/dt = - Tm*d(cp)/dt

and

Tm*(cs-cl)*div(v*beta) = - Tm*div(v*cp)

in the energy-equation. Even with post #75 and Voller. I've done two simulations with the testcase of Fabian. The first one with the TEqn:

fvScalarMatrix TEqn
(
fvm::ddt(cp, T)
+ fvm::div(phiCp, T)
+ L*fvc::ddt(alpha1)
+ L*fvc::div(phi, alpha1)
- Tmelt*fvc::ddt(cp)
- Tmelt*fvc::div(phiCp)
- fvm::laplacian(lambda/rho, T)
);

And the second one with the TEqn derived from Voller:

fvScalarMatrix TEqn
(
fvm::ddt(cp, T)
+ fvm::div(phiCp, T)
+ L*fvc::ddt(alpha1)
+ L*fvc::div(phi, alpha1)
+ Tmelt*(cps-cpl)*fvc::ddt(alpha1)
+ Tmelt*(cps-cpl)*fvc::div(phi, alpha1)
- fvm::laplacian(lambda/rho, T)
);

I changed the equation for h in the second case as well .

The results show significant differences. Which formulation is correct?? Attached you can find the pictures of the liquid fraction after 1080s for the first TEqn and the second TEqn. Speed of the solvers are nearly the same.

kind regards
Florian

Different derivation of enthalpy
 

Hi Florian,

if you compare the enthalpy temperature curve, let’s say in excel, you will see, that the curves are identical except the point of origin. So Vollers equations should give the same results as my equations. However, in Vollers equations, the two Terms

+ Tmelt*(cps-cpl)*fvc::ddt(alpha1)
and
+ Tmelt*(cps-cpl)*fvc::div(phi, alpha1)

depend on alpha1, the value we are trying to converge. This explains the different solutions. Witch equation is better for a faster solution depends on the case.

Regards

Fabian

That’s good. I already thought that this was not your work but sometimes it takes some time to understand why iteration alone will not do the trick. The linearization decouples liquid fraction and temperature. Thus, a consistent solution is no longer achieved. Good luck with writing your thesis and I look forward to see your code.

Fabian

temperature dependent properties
 

Hi AnjaMiehe,
I've followed the process that you had suggested me in the post #79 for temperature dependent properties. What I have seen so far in my problem is that the solution is oscillating and diverging. I've tried with different underrelax values but have no luck so far. Is it possible to share your ideas how you treat the phase change problem with the temperature dependent properties ?

Hi Fabian,

thx for your reply. At the moment I'm testing a little bit with your solver convFinMeltFoam. I used the equation of Voller and a piso algorithm instead of pimple. The result are not to bad at all.

I have a question about the diffusion term in TEqn.H in your solver. In your equation the diffusion of T is:
laplacian(lambda/rho,T)
so you divide the whole equation by rho (rho of the PCM). In my opinion there's a problem if the diffusion term is solved for the fin area. Here rho (for example for Aluminium) is much higher than the rho of PCM. Because of this, the diffusion of the temperature in the fin layer is much higher than in reality.

Currently I'm testing another version of your solver, but without the boussinq appx. I solve all the equations with variable rho (rho=f(T)), like aktive species-transport. I will upload all the solvers when I have finished first tests.
It would be nice if you can give me some feedback about the temperature-diffusion in the fin layer.

kind regards
Florian

Hi All,

attached you can find the solver (PCMSolidPisoFoam) with the equation of Voller and boussineq appx. It is the solver of Fabian Rösler (convFinMeltFoam), but with piso instead of pimple and with a little bit different TEqn.H.

The results looks not to bad. The testcase is the testcase of Fabian with the gallium.

Based on this solver I tried to programm a solver with rho coupled in the UEqn.H, TEqn.H and pEqn.H. without boussineq apprx. For the first shot I used
rho=f(T)=rho_0-rho_0*beta*(T-Tl)
that is the boussineq apprx. :) but it could be every rho = f(T). Rho is solved in all equation. Later I would like to implement cp=cp(T), nu = nu(T) and rho = rho(T). So far so good, the Problem is, that I get bad results with this solver. The problem occurs in the lower wall. Here I used instead of type buoyantPressure a zeroGradient boundary condition. Does anyone know why the solver makes this s***??? The results gets much better with very very fine mesh.

I attached the solver (rhoPCMSolidPisoFoam) and the testcase.

kind regards
Florian

Hi Florian

I am happy that you like he solver.
Concerning your first post you are 100 percent correct when you say, that density differs a lot between PCM and the fin material. This solver is just for a simple check, in which way, a fin or wall influences the melting process. What I did so far is changing the thermal conductivity in a way that the overall thermal diffusivity matches the correct value again.
For the solver in your second post: I once changed my solver the way you plan it. The necessary changes are not that tricky art all. You even could implement the new thermo physical properties class. What could lead to an error is a continuity problem. When the PCM changes phase, it expends in most cases. This increase in volume has to be taken into account by using an inlet/outlet boundary condition. May be this does the trick. Good luck.

Regards

Fabian

Hi,

I read all the posts in this threads and I was wondering if the code can be extended so that a third phase can be simulated. Basically i want to simulate a simple melting problem (say ice -->water) plus the water has a free surface (water-air) so I have the three phases ice-water-air. This means i should somehow combine meltFoam with interFoam. Unfortunately i have a small programming knowledge so any help on the procedure would be appreciated.

Hi!

I'm also still into this and currently am hunting a crash error which probably occurs in a reheating scenario due to recalescence and I'm not totally satisfied with the initial primary solidification and the final concentration of Si (as I'm looking into the solidification of an AlSi alloy)...

But nevertheless - here are some teasers... ;)

melting with third phase
 

Hi callahance

I did exactly what you propose in my thesis. However, I am still in the PhD process and want the solver, results and thesis be public not before my defense of the doctor's thesis. You are right; I combined the compressibleInterFoam solver with my own melting solver to allow a free surface between liquid and gas.

Regards

Fabian

Alright then... thanks for the reply.. I'll try to figure something out :)

adding radiation?
 

Dear melting-foamers,
I tried to add heat transfer from radiation into the solver.
There is a paper http://www.tfd.chalmers.se/~hani/kur...Foam_final.pdf
which explains the basic steps for bouyantSimpleFoam.

What we have in bouyantSimpleFoam, is a equation for entalpy:
This solver works in the buoyantSimpleRadiationFoam case from tutorial.

The melting solver uses follwing hEqn (its not the latest version of it), where I added a radiation term.
It compiles, but at runtime it stops.
which leads at runtime to
Hm...looks like a dimension problem.

Next try:
Understanding the dimensions. For this I put the terms into variables and print out dimensions:

Output:Now I'm confused because of different dimensions between (1) and (2).
Additionaly I think i have to use the thermopysical model library in parallel (polynomialTransport) or build a new one.
(eg like in http://www.cfd-online.com/Forums/ope...cal-model.html)

A lot of text, but simple question :):
Does somebody know the right way to add radiation to the solver?
(there are simlar parallel threads, which did not help me yet:
http://www.cfd-online.com/Forums/ope...am-solver.html
and
http://www.cfd-online.com/Forums/ope...-firefoam.html
)

Thank you for help,
dirk

Phase change with temperature dependent thermo properties
 

Hi,
I am still stuck with this problem. In my cases, with temperature dependent properties the solution oscillates.


In the above equation, I've experienced the problem is basically coming from the transient part .i.e . If I placed the rho*Cp term outside of the ddt term then it converges. So, I'm wondering may be the energy equation have to be written in different way. If anyone have idea ?
I have another queries about following form of energy equation:

Please correct me if I am wrong for solving this:
(1) At old time (say n) I have to calculate h = CpT,
(2) Solving the above equation, h will be found at n+1 time
(3) From , have to calculate as usual way (mentioned by Voller updating scheme)
(4) Then from I have to calculate
(5) Iteration of 2 to 4 until convergence....
Now, my question is that:
(a) Isn't the temperature is one step lagging?
(b) In step (4) what is the best way to calculate temperature from h ?

Thanks in advance.

Hi,

Concerning your dimension problem, if you take a look at for ex. EulerDdtScheme.C, definition for fvc::ddt is

in your case rho is cp, while for fvm::ddt is

Guess it's where dimensions difference comes from.

About radiation term addition:
Can't you just write it like:

as you enthalpy equation is basically temperature equation times specific heat?

thank you alexeym, for your hints!

-> that is also what I thought.

To test the way via another variable, in buoyantSimpleRadiationFoam i defined

and in the solver hEqn.H
which works.

but when I use radiation in melt solver in similar way:
it runs into this (dimension? field?) error

thanks again, i think i will try something else...

Hello everyone,
I am trying to include species eqn to the enthalpy porosity solver and I want to add a loop in case of 0 < alpha < 1 (by using if statement?).
Since I am new to openfoam, I don't know how to do it.
Can somebody help me? Thanks in advance :)

Robert

Would please write down the species equation and condition you want to impose?

Well, actually I am using solver created by Mr. Fabian Roesler (convMeltFoam).
After solving the energy equation, I want to calculate the equilibrium concentration at the interface (0<alpha1<1) by using phase diagram.
Would you give me some advice please?

Regards
Robert

Hi!

I'd like to ask for advice for a solver problem of mine. I'm looking into an alloy phase change problem and ran into some problems with the energy equation. What works so far is this:

This gives me the graph I want as can be seen here with nice primary and eutectic solidification:
Attachment 28356

The problem now is that the above simulation was done for a single cell. When doing simulations with multiple cells and convection this won't work and give strange results. So one has to use a source term just like the first to lines of code with a time derivation and a divergence. I'll leave convection aside for now and do this:

This unfortunately gives me this result I also had presented in an early post of mine:
Attachment 28357

Primary and eutectic solidification sort of exist but the entire recalescence aspect is gone...

I now would like to know if you have any suggestions as to what I should do or give advice as to why the second result is so different. By design both should be the same equation as one might say:



So the source term expression for both approaches should be pretty much the same or am I mistaken here? Another difference is the sign of the source term one has to add or subtract depending on which variant is used (EDIT: This is ok though as I simply defined my temperatureChange field the other way round.)...

Does anyone have an idea what is wrong with the approach with OpenFOAM's fvm::ddt(Temp) and why the temperature won't rise but stay the same?

melting/solidification with temperature dependent thermo-properties
 

Hi,
Has anyone working with melting/solidification modelling considering temperature dependent thermo physical properties (i.e. Cp = Cp(T), K=K(T)) ?
Please give some advice how to deal with this type of problem ? I have used the source based enthalpy method as Fabian, but the problem happens when I make thermophysical properties as variable

Thanks in advance

Hello Shakil,

I am so sorry, I almost forgot about you. At present, I am writing my things together for my PhD and my contract is running out, so I am bit blinkered.

For the implementation of temperature dependent thermophysical properties in general, please take a look here:
http://www.cfd-online.com/Forums/ope...ies-false.html

Then, going on from post #56, use the interpolation procedure within the do-loop. When you implement it, start with the thermophysical post and get that one to run. After that, if you go for solidification / melting with dependent properties, start with a weak dependence and build up from there.

I am going to publish my code as soon as I have my PhD ready, but that might take some months.

Best Regards, Anja

Anybody can help me?when I compile meltFoam with wmake command,i face a problem.
 

Anybody can help me?when I compile meltFoam with wmake command,i face a problem about compiling meltFoam:

fuguang@fuguang-Veriton-D630:~/Downloads/meltFoam_solver$ wmake
SOURCE=meltFoam.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-100 -I/opt/openfoam230/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/openfoam230/src/OpenFOAM/lnInclude -I/opt/openfoam230/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/meltFoam.o
In file included from meltFoam.C:83:0:
pEqn.H: In function ‘int main(int, char**)’:
pEqn.H:8:11: error: ‘ddtPhiCorr’ is not a member of ‘Foam::fvc’
/opt/openfoam230/src/finiteVolume/lnInclude/readTimeControls.H:38:8: warning: unused variable ‘maxDeltaT’ [-Wunused-variable]
make: *** [Make/linux64GccDPOpt/meltFoam.o] Error 1

Version mismatch
 

Hi wumin

This sounds to me like a version mismatch. Which OpenFOAM version do you use? I bet you use OpenFOAM 2.3.x since the ddtPhiCorr in pEqn was changed to a more general ddtCorr. Have a look into other pEqn of standard solvers in OF and check for differences

pEqn.H:8:11: error: ‘ddtPhiCorr’ is not a member of ‘Foam::fvc’

Regards
Fabian

Hi Fabian,
thanks for your reply.You are right! The version I now use is OpenFoam-2.3.0.
My work is mainly to simulate the melting and solidification of selective laser melting（3D printing）.
I want to compile the meltFoam,but always fail.Can you tell me how to deal with the problem?
the system I install is ubuntu 12.04.

Hi Wumin,
Have tried with the meltFoam version provided by Fabian in post#19 (http://www.cfd-online.com/Forums/ope...g-problem.html) ?

Enthalpy linearization
 

Hi

I have been working on the melting problem for one month and more. I am using the solver from Fabian as a test bed to implement the enthalpy linearization. while in this method, it is a bit complex to implement it in OpenFOAM. Somebody has gone through this method in the past?

Best Regards
Rohith

Hi RaghavendraRohith,
It sounds cool that someone also trying with the enthalpy linearization method. I'm using linearization method based on Voller's source based method i.e. linearization of the term. I have also validated my case with the gallium melting and the results are cool! But now stuck with different types of problem...

#shakil

Hi Shakil

Actually I am not asking for source based method

Please go through it, let me know your idea if so.


Regards
Rohith

a wmake problem
 

Hello! Fabian
I have some problem to wmake meltFoam on OF2.2.0. Can you help me to solve this problem. I'm a new user and I want to know if I can use this Foam to simulate a solid particle dissolution in liquid.
Thanks very much!


Making dependency list for source file meltFoam.C
SOURCE=meltFoam.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-100 -I/opt/openfoam220/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/openfoam220/src/OpenFOAM/lnInclude -I/opt/openfoam220/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/meltFoam.o
/opt/openfoam220/src/finiteVolume/lnInclude/readTimeControls.H: In function ‘int main(int, char**)’:
/opt/openfoam220/src/finiteVolume/lnInclude/readTimeControls.H:38:8: warning: unused variable ‘maxDeltaT’ [-Wunused-variable]
g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-100 -I/opt/openfoam220/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/openfoam220/src/OpenFOAM/lnInclude -I/opt/openfoam220/src/OSspecific/POSIX/lnInclude -fPIC -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPOpt/meltFoam.o -L/opt/openfoam220/platforms/linux64GccDPOpt/lib \
-lfiniteVolume -lOpenFOAM -ldl -lm -o /home/hwy/OpenFOAM/hwy-2.2.0/platforms/linux64GccDPOpt/bin/meltFoam

hello! Have your problem been solved?

Hello Houwy

This is not fabian, but i can help you in this aspect

I do not see what actually is your problem, the solver is compiled now and may be you can simulate your case, yes meltFoam is a solver to solve melting of a solid particle or a solid body.

Go forward with your simulation once

Regards
Rohith

Solver Update Fabian Rösler OF2.3.0
 

Hello everyone,

here is Fabians convMeltFoam solver of post #81 for you compiling in OpenFOAM 2.3.0 with the test case, as one boundary condition for the pressure changed its name. Have fun using it.

Regards, Anja

help
 

Hello,
Fabian.Thank you for your solver about melting.But I have some questions with your solver. In the UEqn.H , you use a symbol "DC" calculated form DCl,DCs and alpha. Can you explain what are their meanings ? And can you give me your mathematical model?