CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

melting problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree74Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 20, 2013, 10:25
Default solid/liquid phase change solver convMeltFoam
  #81
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
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
Attached Files
File Type: gz convFinMeltFoam.tar.gz (8.8 KB, 122 views)
File Type: gz convMeltFoam.tar.gz (7.5 KB, 167 views)
File Type: gz lduMatrix.tar.gz (62.6 KB, 116 views)
fabian_roesler is offline   Reply With Quote

Old   December 20, 2013, 11:41
Default
  #82
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
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?
ahmmedshakil is offline   Reply With Quote

Old   December 21, 2013, 09:46
Default
  #83
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
Hallo Fabian,

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

Viele Grüße
Florian
riesotto is offline   Reply With Quote

Old   December 21, 2013, 15:40
Default
  #84
Member
 
Anja Miehe
Join Date: Dec 2009
Location: Freiberg / Germany
Posts: 48
Rep Power: 7
AnjaMiehe is on a distinguished road
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.
ahmmedshakil likes this.
AnjaMiehe is offline   Reply With Quote

Old   December 24, 2013, 09:04
Default
  #85
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
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
riesotto is offline   Reply With Quote

Old   December 24, 2013, 13:01
Default
  #86
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi riesotto,
Please see the post #75, and I guess you could have idea about that.
ahmmedshakil is offline   Reply With Quote

Old   December 25, 2013, 07:47
Default
  #87
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
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
riesotto is offline   Reply With Quote

Old   December 25, 2013, 09:55
Default
  #88
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
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
Attached Images
File Type: jpg TEqn_1.jpg (26.6 KB, 83 views)
File Type: jpg TEqn_2.jpg (27.3 KB, 75 views)
riesotto is offline   Reply With Quote

Old   December 26, 2013, 05:17
Default Different derivation of enthalpy
  #89
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
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
fabian_roesler is offline   Reply With Quote

Old   December 26, 2013, 05:26
Default
  #90
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
Quote:
Originally Posted by AnjaMiehe View Post
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.
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
fabian_roesler is offline   Reply With Quote

Old   December 28, 2013, 19:30
Default temperature dependent properties
  #91
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
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:
Originally Posted by AnjaMiehe View Post
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.
ahmmedshakil is offline   Reply With Quote

Old   January 2, 2014, 10:15
Default
  #92
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
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
riesotto is offline   Reply With Quote

Old   January 3, 2014, 13:41
Default
  #93
New Member
 
Florian
Join Date: Jan 2011
Location: Mannheim, Germany
Posts: 24
Rep Power: 6
riesotto is on a distinguished road
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
Attached Files
File Type: zip PCMSolidPisoFoam.zip (5.0 KB, 59 views)
File Type: zip PCMSolidPisoFoamCase.zip (59.9 KB, 38 views)
File Type: zip rhoPCMSolidPisoFoam.zip (9.2 KB, 49 views)
File Type: zip rhoPCMSolidPisoFoamCase.zip (60.6 KB, 34 views)
shuisheng and farzinkhan like this.
riesotto is offline   Reply With Quote

Old   January 3, 2014, 15:11
Default
  #94
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
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
fabian_roesler is offline   Reply With Quote

Old   January 9, 2014, 04:43
Default
  #95
New Member
 
Join Date: Mar 2012
Posts: 25
Rep Power: 5
callahance is on a distinguished road
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.
callahance is offline   Reply With Quote

Old   January 10, 2014, 09:38
Default
  #96
Member
 
Thomas Vossel
Join Date: Aug 2013
Location: Germany
Posts: 45
Rep Power: 4
ThomasV is on a distinguished road
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...
Attached Images
File Type: jpg Graph_Temp_alpha.jpg (30.5 KB, 50 views)
File Type: jpg Graph_weightFracSi.jpg (28.6 KB, 37 views)
ThomasV is offline   Reply With Quote

Old   January 11, 2014, 05:42
Default melting with third phase
  #97
Senior Member
 
fabian_roesler's Avatar
 
Fabian Roesler
Join Date: Mar 2009
Location: Bad Friedrichshall, Germany
Posts: 164
Rep Power: 8
fabian_roesler is on a distinguished road
Quote:
Originally Posted by callahance View Post
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 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
callahance likes this.
fabian_roesler is offline   Reply With Quote

Old   January 13, 2014, 05:42
Default
  #98
New Member
 
Join Date: Mar 2012
Posts: 25
Rep Power: 5
callahance is on a distinguished road
Alright then... thanks for the reply.. I'll try to figure something out
callahance is offline   Reply With Quote

Old   January 16, 2014, 07:02
Default adding radiation?
  #99
dzi
Member
 
Join Date: Nov 2011
Location: Berlin
Posts: 31
Rep Power: 5
dzi is on a distinguished road
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:
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) 
);
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.
Code:
fvScalarMatrix hEqn
    (
        fvm::ddt(cp, T)
      + fvm::div(phi*fvc::interpolate(cp), T)
      + hs*exp(-pow((T-Tmelt)/Tdim,2))/Foam::sqrt(constant::mathematical::pi)/Tdim*fvm::ddt(T)
      + hs*exp(-pow((T-Tmelt)/Tdim,2))/Foam::sqrt(constant::mathematical::pi)/Tdim*(U & fvc::grad(T))
      - fvm::laplacian(lambda/rho, T)
      + radiation->Sh(thermo) 
    );
which leads at runtime to
HTML Code:
--> FOAM FATAL ERROR: 
incompatible fields for operation 
    [T] + [h]
Hm...looks like a dimension problem.

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)
Output:
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)
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 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
dzi is offline   Reply With Quote

Old   January 16, 2014, 09:02
Default Phase change with temperature dependent thermo properties
  #100
Senior Member
 
Mohammad Shakil Ahmmed
Join Date: Oct 2012
Location: AU
Posts: 123
Rep Power: 5
ahmmedshakil is on a distinguished road
Hi,
I am still stuck with this problem. In my cases, with temperature dependent properties the solution oscillates.
\frac{d (\rho Cp\, T)}{dt} +  div(\rho Cp\, U\, T) = laplacian(K\, T) + ST

In the above equation, I've experienced the problem is basically coming from the transient part .i.e \frac{d(\rho Cp, T)}{dt}. 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:
\frac{d \rho h}{dt} + \nabla .(\rho V h) = \nabla (\alpha \nabla h) - \frac{d (\rho \Delta H)}{dt}
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 h^{n+1}, have to calculate \Delta H^{n+1} as usual way (mentioned by Voller updating scheme)
(4) Then from h^{n+1} I have to calculate T^{n+1}
(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.
ahmmedshakil is offline   Reply With Quote

Reply

Tags
melting openfoam

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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 Se-Hee 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


All times are GMT -4. The time now is 21:06.