[fabian_roesler]

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

[ahmmedshakil]

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?

[riesotto]

Hallo Fabian,

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

Viele Grüße
Florian

[AnjaMiehe]

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.

[riesotto]

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

[ahmmedshakil]

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

[riesotto]

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]

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

[fabian_roesler]

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]

Originally Posted by AnjaMiehe

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

[ahmmedshakil]

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 ?

Originally Posted by AnjaMiehe

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.

[riesotto]

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]

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

[fabian_roesler]

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

[callahance]

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.

[ThomasV]

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...

[fabian_roesler]

Originally Posted by callahance

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]

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

[dzi]

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:

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.

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

--> 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:

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:

#################   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 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

[ahmmedshakil]

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.