
[Sponsors] 
December 20, 2013, 10:25 
solid/liquid phase change solver convMeltFoam

#81 
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8 
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 

December 20, 2013, 11:41 

#82 
Senior Member
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5 
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 conductionconvection phase change problem and for some complex problems. Have you ever tried with the temperature dependent properties? 

December 21, 2013, 09:46 

#83 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
Hallo Fabian,
vielen Dank fürs uploaden. Zieh mir den Solver direkt rein. Viele Grüße Florian 

December 21, 2013, 15:40 

#84 
Member
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7 
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. 

December 24, 2013, 09:04 

#85 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
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+(1beta)*cs]*T+beta*(cscl)*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+(1beta)cs)): d(cp*T)/dt + div(v*cp*T) + L*d(beta)/dt + L*div(v*beta)  Tm*(cscl)*d(beta)/dt + Tm*(cscl)*div(v*beta)  div(lambda/rho*grad(T)) = 0 (Eq3) This equation is very similar to yours except [Tm*(cscl)*d(beta)/dt ] and [Tm*(cscl)*div(v*beta)] Now my question Why is term 5 in your solver: Tm*(cscl)*d(beta)/dt = Tm*d(cp)/dt and term 6 in your solver: Tm*(cscl)*div(v*beta) = Tm*div(v*cp) ??? is (cscl)*beta=cp I have the same Problem with h (TEqn.H Z.40). Here: h=cp*(TTmelt)+alpha1*L why not: h=cp*T + alpha1*(cscl)*Tmelt + alpha1*L Thanks again for your solver and Frohe Weihnachten Florian 

December 24, 2013, 13:01 

#86 
Senior Member
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5 
Hi riesotto,
Please see the post #75, and I guess you could have idea about that. 

December 25, 2013, 07:47 

#87 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
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 

December 25, 2013, 09:55 

#88 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
Hi,
until now i don't understand why: Tm*(cscl)*d(beta)/dt =  Tm*d(cp)/dt and Tm*(cscl)*div(v*beta) =  Tm*div(v*cp) in the energyequation. 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*(cpscpl)*fvc::ddt(alpha1) + Tmelt*(cpscpl)*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 

December 26, 2013, 05:17 
Different derivation of enthalpy

#89 
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8 
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*(cpscpl)*fvc::ddt(alpha1) and + Tmelt*(cpscpl)*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 

December 26, 2013, 05:26 

#90  
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8 
Quote:
Fabian 

December 28, 2013, 19:30 
temperature dependent properties

#91  
Senior Member
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5 
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 ? Quote:


January 2, 2014, 10:15 

#92 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
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 speciestransport. 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 temperaturediffusion in the fin layer. kind regards Florian 

January 3, 2014, 13:41 

#93 
New Member
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6 
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_0rho_0*beta*(TTl) 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 

January 3, 2014, 15:11 

#94 
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8 
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 

January 9, 2014, 04:43 

#95 
New Member
Join Date: Mar 2012
Posts: 25
Rep Power: 5 
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 (waterair) so I have the three phases icewaterair. 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. 

January 10, 2014, 09:38 

#96 
Member
Thomas Vossel
Join Date: Aug 2013
Location: Germany
Posts: 45
Rep Power: 4 
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... 

January 11, 2014, 05:42 
melting with third phase

#97  
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8 
Quote:
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 

January 13, 2014, 05:42 

#98 
New Member
Join Date: Mar 2012
Posts: 25
Rep Power: 5 
Alright then... thanks for the reply.. I'll try to figure something out


January 16, 2014, 07:02 
adding radiation?

#99 
Member
Join Date: Nov 2011
Location: Berlin
Posts: 31
Rep Power: 5 
Dear meltingfoamers,
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: Code:
fvScalarMatrix hEqn ( fvm::div(phi, h)  fvm::Sp(fvc::div(phi), h)  fvm::laplacian(turbulence>alphaEff(), h) ==  fvc::div(phi, 0.5*magSqr(U), "div(phi,K)") + radiation>Sh(thermo) ); 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. Code:
fvScalarMatrix hEqn ( fvm::ddt(cp, T) + fvm::div(phi*fvc::interpolate(cp), T) + hs*exp(pow((TTmelt)/Tdim,2))/Foam::sqrt(constant::mathematical::pi)/Tdim*fvm::ddt(T) + hs*exp(pow((TTmelt)/Tdim,2))/Foam::sqrt(constant::mathematical::pi)/Tdim*(U & fvc::grad(T))  fvm::laplacian(lambda/rho, T) + radiation>Sh(thermo) ); HTML Code:
> FOAM FATAL ERROR: incompatible fields for operation [T] + [h] Next try: Understanding the dimensions. For this I put the terms into variables and print out dimensions: Code:
Foam::fvScalarMatrix fsm_ddt_cpT = fvm::ddt(cp, T); (1) Foam::volScalarField vsf_ddt_cpT = fvc::ddt(cp, T); (2) Foam::fvScalarMatrix fsm_radiationSource = rad.Sh(thermo); (3) Info<< "################# dimensions: [kg m s K kgmol A cd]" << nl << endl; Info<< "################# fsm_ddt_cpT.dimensions() " << fsm_ddt_cpT.dimensions() << nl << endl; (4) Info<< "################# vsf_ddt_cpT.dimensions() " << vsf_ddt_cpT.dimensions() << nl << endl; (5) Info<< "################# fsm_radiationSource.dimensions() " << fsm_radiationSource.dimensions() << nl << endl;(6) Code:
################# dimensions: [kg m s K kgmol A cd] ################# fsm_ddt_cpT.dimensions() [0 5 3 0 0 0 0] // = W/kg/m^3 (7) ################# vsf_ddt_cpT.dimensions() [0 2 3 0 0 0 0] // = N*m/kg*s = W/kg (8) ################# fsm_radiationSource.dimensions() [1 2 3 0 0 0 0] // = N*m/s = W (9) Additionaly I think i have to use the thermopysical model library in parallel (polynomialTransport) or build a new one. (eg like in How creating new thermo physical model) 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: How to add Radiation to existing buoyantBoussinesqSimpleFoam solver and writing radiation source term in fireFoam ) Thank you for help, dirk Last edited by dzi; January 16, 2014 at 07:11. Reason: links dont work 

January 16, 2014, 09:02 
Phase change with temperature dependent thermo properties

#100 
Senior Member
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5 
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. Last edited by ahmmedshakil; January 17, 2014 at 09:22. 

Tags 
melting openfoam 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
conduction problem  venkataramana  OpenFOAM  3  December 1, 2013 08:30 
natural convection problem for a CHT problem  SeHee  CFX  2  June 10, 2007 06:29 
Adiabatic and Rotating wall (Convection problem)  ParodDav  CFX  5  April 29, 2007 19:13 
Melting Problem  M  FLUENT  0  April 29, 2007 16:07 
Is this problem well posed?  Thomas P. Abraham  Main CFD Forum  5  September 8, 1999 14:52 