Hi

Soon, there will be my PhD thesis available. I am in the final state of printing and it will also be available as eBook. I will post the link here. Moreover, my former colleges from the chair of engineering thermodynamics and transport processes (LTTT) at the University of Bayreuth an I prepare a new article.

Cheers Fabian

Hi Fabian

I did a quick test today and the short conclusion is that the convMeltFoam, which uses inner iterations to update alpha, gives different result to erfConvectionMeltFoam,which uses continuous alpha function. The plots given in my earlier thread were due to the former solver and the later solver gives result sort of similar to what you presented in your paper.

@fabian

Sure, would be happy to read your thesis. Kindly post it here when you are done with it.

@YS

Can I have access to the paper you are referring to?

Cheers

Hi Fabian,

I have attached two cases below, galliumErfMelting and LHTES, and they are exactly prepared according to your 2011 paper. I am using erfConvectiveMeltingPimpleFoam on OpenFOAM 2.1.1. The validation on galliumErfMelting is perfect but the result of the other case, LHTES, is somehow weird. After 10h of simulation time (36000 sec), the wax is hardly melting. Please help to have a quick look to see if there's anything I missed. Thanks a lot for the help!

Attachment 32856

Attachment 32857

Attachment 32858

I tried to replace the wax by gallium in the LHTES case and the material is melting normally.

convMeltFoam Problem
 

Dear Folks,

for my Master-Thesis i want to simulate the partial melting of a moving plate (quadratic) under a stationary heat source (laser-welding).

I successfully downloaded and compiled the convMeltFoam - Solver for OF 2.3. For a non-moving solid region everything works just perfect.
But when i added the modifications to the solver which has been posted to simulate a moving solid region the Temperature-distribution doesn't change and the pressure rises rapidly. Also the simulation crashes after a few timesteps.

Modifications in createFields.H
Modifications in UEqn.H

My boundarys are:
U
Us:

p_rgh

I also tried different variations with Us and p_rgh, without any success.

Thanks in advance for any help.

Best
Phil

@Anja Miehe: Thanks for the solver. I have OF2.3.0 and while compiling the solver I get following warnings:

There is no other warnings or errors. The tutorial attached is also running without any errors.

Do you think these warnings are serious? If so how can I modify the solver code to avoid them?

warning: unused variable
 

Hi,

those warnings appear because of the curly brackets around the TEqn. The variables are defined outside these brackets and are only used inside. Thus, the compiler does not recognized their usage. So everything is fine.

@ Phil: I need some time to study your problem

Cheers

Fabian

@fabian_roesler: Thank you for the prompt reply.

Hy Fabian,

i did proceed with my problem. The key was to deactivate the Momentum-Predictor and to use uncorrected div-schemes. Furthermore it is necessary to define the inlet and outlet velocity boundaries in the U-File. Only the internal field values from Us are needed.

Now the temperature distribution looks way better and the velocity magnitude (flow of the molten metal around a keyhole (capillary) in the middle of the plate ) looks more realistic than before. (See attached picture)
Nevertheless the pressure distribution looks quite unrealistic. I would expect a pressure rise in front of the capillary to force the fluid around the capillary.


P.S. The metal is completely fluent in the capillary region, thus a collision with the soild is not the problem.
P.S.2 If i change the inlet velocity boundary to fixedValue, everything explodes.

Best Phil

Hullo Fabian

I have conducted some test runs and compared with the experimental data available as a part of my master thesis. It was interesting to note that once the solid melted, a lot of noise in the temperature can be observed. I am not really sure if its numerical or physical. I checked for an other probe at a different location where the solid has not melted yet. Here no oscillations can be observed. Can you have a look at the plots I attached and kindly comment on it?

Cheers

Hi Akash,
Have you checked with the fined mesh ?

cheers
shakil

Hi

Earlier in this thread I posted that the erfConvectiveMeltingFoam is a method to decrease the nonlinear behavior of the energy conservation equation. However, you will not get rid of it entirely. Thus, when the influence of the convective term increases, the error of the energy conservation increases too. In short, cases with high conductivity of the PCM work well with this solver. Materials like wax or plastics with low thermal conductivity and proportionally high convective heat transfer give erroneous results. This is the reason I switched to the iterative method, proposed by Voller, implemented in the convMeltFoam. This method is quite slow but conserves energy when used correctly. I hope this clarifies your situation a bit.


This could be due to the Darcy source term. This source term forces velocity to zero in the solid and transition region. However, as you use PISO or PIMPLE with no momentum predictor, momentum is not solved for directly. You rather calculate the velocity from your pressure field. Thus, the pressure in solid cells with flow against them in the near fluid cells will increase to 'resist' the flow. In contrary, the pressure in cells with no or negative flow against them will be lower. Always imagine the solid as a kind of high viscous fluid.

Cheers

Fabian

Hi Ahmed

The mesh is fine but not finest. I did a grid study for my geometry and I took mesh levels of 0.5, 1 and 2 million. There is a noticeable difference from half to one million. But one to two million has little change. Considering the computational resources available to me, I opted for 1 million cells which looked like an optimum solution. Maybe taking 2 million cells will dampen the oscillitations but I am not sure if they will eliminate them entirely. What do you think?

Cheers
Akash

Hi Akash

From the plots attached, after the melting point of the solid. The solver tends to be chaotic due to the modeling of convection. Which algorithm are you using the one with error function or an iterative approach?

Best Regards
Rohith

Hi Rohith

The solver I am using is convMeltFoam from Fabian. Could you be more specific about what you mean by the algorithm employed? You mean the pressure corrector or something else?

Cheers
Akash

Hi Akash

The algorithm, i mentioned means the approach used to model the temperature and Enthalpy coupling. However, i see that you have used ConvMeltFoam which means the temperature follows source based iterative approach by Voller et al., i may suggest you to reduce the time step size as well with different mesh refinements to optimize the courant number.

The oscillations you mentioned above could be possible due to the iterations. But it seems a bit strange to me that the slope and curvature vary very large from the experimental results.

My call would be try to make the approximation more stiff in modeling temperature.

Regards
Rohith

Hi Rohith

Thank you for the suggestions. I would try to implement them and see how it goes.

As you indicated, the slopes are quite different but after having a word with the testing department we have come to conclusion that the experimental values are also error prone. Especially probe 2(the first picture) as it is located between the heaters (my geometry has two parallel heaters in a tank). It was quite difficult to measure with reasonable accuracy.

But probe 17(second picture) is far away from the heaters and shows good relation with experimental values. Nonetheless, it would be interesting to see what would happen with probe 17 once melting front reaches the location. I have not yet run that far.

Cheers
Akash

Hi guys,

I´m having some doubts about how to implement the source terms in Openfoam. Currently there seems to be two approaches:

1) Fabian Rösler's one, based on writing the enthalpy equation in terms of the liquid fraction writing the terms which depends on the liquid fraction as explicit ( fvc:: )

2) Voller's approach, based on linearized source terms as follows:







Whith both of them correcting the liquid fraction after each iteration

So which one would the best one to implement? I mean, Voller didn´t take account of convection, so the convection source term doesn´t appear in his equations. But Voller´s method has been validated in numerous testcases. Also I've validated Rösler approach in easy-geometry testcases, but I didn´t use it in complex-geometries.

Thank you all for the cooperation

Hi anonymous,
First of all, both of the methods are good. But for the implementation it depends on what type of problem you are solving. Linearised Source Based Method (LSM) is robust, and it just takes (at least) 3-4 iteration for convergence. Second method, which called Fictitious Source Method (FSM), takes a little bit more iteration to converge.
Now for the use, it depends on how the thermal gradient is changing. If thermal gradient changes strongly with times, then LSM may give some false results for the melt front and you may observe some false diffusion. In this case, my personal observation, it is better to go with FSM. Another problem for the LSM is to get the convergence of the results.
Also, none of the method is good for isothermal phase change.
If thermo-physical properties are function of temperature, then LSM will cause problem unless careful steps are not taken.

cheers,
#shakil

PhD Thesis on modeling and simulation of phase change processes in macro encapsula
 

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

Very nice! I might be able to use this at a later time so I'll keep it as reference.
Do you think it will be usable for highly transient phase changes, including evaporation?

Dear Ahmmed,

Thank you for your answer, I really appreciate it. So it's prefered yo use FSM implementation just as Rösler did in his convMeltFoam solver and correct the liquid fraction as explained in http://www.cfd-online.com/Forums/ope...tml#post467200 #80 Am I right?

Thank you

Rösler Thank you for all the work you did I will have a look in your thesis. Is there any way of having the index of your work so we can see the topics which are discussed?

method
 

Hi anonymous,
As I said earlier, it depends on the problem you are solving. And now I can reckon FSM will be good to use.

Cheers
#shakil

I am happy that you like my work. The solve is designed for solid/liquid phase change. All the Darcy term and stuff would not be necessary for evaporation. However, its implementation should not be that hard.

@ssss: I will post the index the next days.

Cheers

Fabian

I think that I have been a bit unclear.
I am interested in simulating a highly transient solid->liquid->gas phase transition of a polymer, not just liquid->gas. I'm generally interested in simulating ablation flows (i.e. determine amount of ablated mass and temperature at a wall) and corresponding wall temperatures of a polymer surface in contact with a plasma with T on the order of some 10^4 K.

I haven't looked into the theory any further than common heat conduction equations yet, as this is only a small aspect of my work, so I can't judge myself right now if your method is suitable for my cause. I would expect that the Darcy term is still required in this case because there will still be a solid phase which is not moving?

I have also seen something similar under the term "2-phase-Stefan-problem". This is for tracking the phase-changing interfaces in a 1-dimensional case and might also be an option for me, as I expect my problem to be more or less 1D (since energy flow is largely 1D in the direction normal to the wall).

Hi chriss85

I'd say that solid>liquid>gas transition should be possible. "The only thing" you have to do is implementing the liquid to gas transition into my solver. I put that into quotes because this is no easy task. You have to account for the interface between liquid and gas by using an appropriate method like VOF etc. and extend it for phase change between the two phases. Have a look into these threads:

http://www.cfd-online.com/Forums/ope...e-boiling.html
http://www.cfd-online.com/Forums/ope...mentation.html
http://www.cfd-online.com/Forums/ope...nge-vof-2.html

Bitan Shu programmed a solver for liquid>gas transitions using interFoam as starting solver. However, there could be another solver, closer to what you search. Dig a little bit into the forum and you might find what you need.

To make it short: You can not solve your problem by just using my solver. You will have to invent something on your own or at least extend some existing solver including mine. This will be no easy task. But keep foaming

Cheers

Fabian

Hi FOAMers',
Have anyone worked with the chtMultiregionFOAM for melting simulation ?( OR Have anyone add the melting solver with chtMultiRegionFoam ? ) If so, it will be great if you could share the ideas how you plugged the code....

Thanks in advance
#shakil

What's your interest for using chtMultiregionFOAM? Grid Moving so that liquid and solid could be treated alone? Adding some solid-fins to modify the heat transfer?

chtMultiregion+melting
 

Hi SSS,
I am solving a problem for multilayer systems. In my case I have two layer, and both of the layer has different properties and heat absorption. I have already model the case before melting. Now, I want to simulate the scenario when the top layer starts to melt; actually melting depth for the top layer. In my case, I am pretty much sure only the top layer will melt. Therefore, I have to plug in the melting code for top layer equation. If anyone can advice me how to do that, it will be great !.

Thanks in advance.
#shakil

Maybe you should get the chtMultiRegion code, get rid of the turbulence parts and write your own TEqn to be solved in each of the fluid Regions.

I'm not too experienced with chtMultiRegion, but you should be able to implement it, although you will need to modify the createFields.H and the readTransportProperties, in order to get the need values for the solidification TEqn

I would like to ask also a question about Voller Linearized Source Term implementation:

Whenever I include TEqn.H inside the PIMPLE loop (in the outerloop or in the pressure correction loop) my solution advances in time a lot, let's say it should take 1000s to have a liquid fraction of 95%, but puting TEqn.H inside loop takes the time down to 50s, although the cualitative results are fabolous, it goes much quicker.

If I put the TEqn.H outside the PIMPLE loops (both of the loops) then I obtain the results I expect to have.

This may be a small source of errors, because I'm not changing the temperature during both of the PIMPLE loops, and the temperature and the boussinesq density should change whenever the velocity field changes.

Anyone knows why is this happening?

Thank you very much.

Hello,

what are DCs and DCl in TEqn in meltFoam?

Mithil.

They form the constants in the porosity function:



The "C" in the equation is the DCl constant. The "b" in the equation is the "DCs" constant. The names might be confusing here as DCl means something like "Darcy constant large" and DCs means "Darcy constant small" - so nothing in terms of liquid or solid... ;)

Both constants are defined in the transport properties and DCs is just there to prevent a division by zero (for that purpose a tiny number gets added to the denominator) and DCl is a more or less random "big number" meant to lead to a high number for the porosity function which then leads to decreased velocities in the momentum equation. While its number is not really random as it influences how fast the velocities will decrease, it's not really a fixed number you can look up but something you adjust to your case so the results match your expectations for the respective material you're simulating...

Hi,

readTransportProperties.H
hEqn.H
DCl should be large (~10^8), DCs is used to avoid division by zero and should be small (~10^-15 for example).

species
 

Hi Yahoo
I am working on the segregation with openfoam. Have you been successful in the species calculation? If so, would you pleased to give me some advice, or show me your solver? Thank you !

Hi all,

I have a problem with modelling the melting/solidification of pcm in sphere by using the enthalpy method. I have to do this modelling using matlab and by using finite difference method.

Hello all,
I have went through all your posts and you are really doing a good job. Might my post is in the wrong place, but I decided to give it a try, since I have posted it in a different section.

I have read many papers and they mentioned I have to change the enthalpy governor equation to the following, the enthalpy term to a temperature or the temperature term to enthalpy to model it easily in matlab.

I have changed the temperature into enthalpy and now, all the terms in the implicit method equation are enthalpy. BUT, what about the boundary conditions, Am I need them, if so, I have them interms of temerature, how to change them into enthalpy.

any help will be greatly appreciated


Best

Hey Fabian,

I'm working on a problem you might have already solved in your PhD and therefore I would like to ask you for a hint.

My Master thesis is about the simulation of casting processes including the filling, temperature transport and solidification and thus exactly the reversal of your subject.
In order to add the gas phase I'm trying to combine interFoam and convMeltFoam, as you suggested earlier.

Using the model of your PhD, everything seems to work out fine except the solution of the TEqn, which exploits whenever I include the calculation of the liquid fraction fl.

The only thing I changed so far is the update equation of the liquid fraction. I replaced the heat capacity cp of the hole cell with a cp calculated from just the liquid and solid phase.

I suggest that the new possibility of changing the fl by filling a cell instead of just melting solid is causing the trouble. By filling a cell I can produce huge differences of fl between two time steps, which the iteration algorithm for T and fl can't handle anymore.

Do you come to that problem?

My supervisor mentioned a smoothing function to restrict the time step if huge fl changes happen could be a solution. I haven't tried it yet, but will post again, if something works.

I have to excuse, if the solution is described in your PhD. We already ordered a book through inter-library loan, but it seems to be loaned. Do you know about a library online version? If I haven't got access by next week I will buy it. (Until then I'm dealing with the googleBooks prev)

Thanks,

ole

Dear kuechenrule,

The TEqn is known to create troubles in interFoam, even without working melting problems. Are you working with the fields multiplied by the density? Something like rhoCp, rhoPhiCp, etc. Which version of interFoam are you using?

Are you ensuring, that the TEqn and the liquidFraction converge in each of the PIMPLE/PISO iterations? Did you try if your code works without gas phase in it?

If you want to limit the timeStep, you could just change the maxDeltaT in the controlDict to the maximun value you want to achieve.

Anyway I'm pretty sure you have taken care of looking at this small problems. I will be glad if I could help you trying to debug the errors that the solver gives

Dear SSSS,

thanks for your help! Even if I used OpenFoam a few times, I'm new in the field of programming and manipulation. Thus please mention all small problems or mistakes that you suggest I could have done. :)

I'm working with rhoCp, rhoCpPhi, etc, since I consider big density differences. between gas and melt (liquid+solid).
I'm using OF 2.3.1 (are there specific versions for each solver?)
I do ensure convergence of the TEqn and fl by a fl residuum of e-6 in and a maximum iteration loop number of 150. Teqn is solved inside the PIMPLE loop.
The code seems to work without a gas phase (testcase: filled open cube with cold walls; just T transport no movement at all).

I'll will play with the time step a bit and try to find out, if it works. Thanks for the advice.

An other question arises today for me. The PIMPLE loop features the following order: alphaEqn -> UEqn -> TEqn -> pEqn. At which point do I update rho and rhoCp etc, which are dependent on the melt fraction (alpha1) and the liquid fraction (fl)? If I update rho outside the TEqn loop, the TEqn solves a different rho as the UEqn. If I update rho inside the TEqn, UEqn solves a different rho as pEqn. I think you got the point. Or does it not matter since the PIMPLE loop runs until all rho's are nearly equal?

I attached my TEqn. Please don't get me wrong, I don't want you to do my work, but maybe I did some obvious mistakes you can see.

alpha1 - melt fraction;
fl - liquid fraction of melt;
fs - solid fraction of melt;
rho1 - density of pure liquid;
rho2 - density of pure gas;
rhos - density of pure solid;

{
int iter = 0;
scalar residual = 1;
scalar meanResidual = 1;

do
{
iter++;
fl.storePrevIter();

rhoCp = alpha1*(fl*rho1*cp1 + fs*rhos*cps) + alpha2*rho2*cp2;
lambda = alpha1*(fl*lambda1 + fs*lambdas) + alpha2*lambda2;
rhoAlpha1 = alpha1*(alpha1*(fl*rho1 + fs*rhos) + alpha2*rho2);

rhoCpPhi = tphiAlpha()*fvc::interpolate(rho2*cp2 - (fl*rho1*cp1 + fs*rhos*cps)) + phi*fvc::interpolate(fl*rho1*cp1 + fs*rhos*cps);
rhoAlpha1Phi = fvc::interpolate(alpha1)*(tphiAlpha()*fvc::interpo late(rho2 - (fl*rho1 + fs*rhos)) + phi*fvc::interpolate(fl*rho1 + fs*rhos));

fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T)
+ fvm::div(rhoCpPhi, T)
+ L*fvc::ddt(rhoAlpha1, fl)
+ L*fvc::div(rhoAlpha1Phi, fl)
- Tmelt*fvc::ddt(rhoCp)
- Tmelt*fvc::div(rhoCpPhi)
- fvm::laplacian(lambda, T)
);

TEqn.solve();

Tcorr = (Tl-Ts)*fl+Ts;
fl = max(min(fl+alphaRel*(fl*cp1 + fs*cps)/L*(T-Tcorr),scalar(1)),scalar(0));
fs = scalar(1) - fl;

residual = max(mag(fl.internalField()-fl.prevIter().internalField()));
meanResidual = sum(mag(fl.internalField()-fl.prevIter().internalField())*mesh.V())/sum(mesh.V()).value();

} while ((iter < minTCorr || residual > alphaTol ) && iter <= maxTCorr);

cp = alpha1*(fl*cp1 + fs*cps) + alpha2*cp2;
h = cp*(T-Tmelt)+alpha1*fl*L;
}

Thanks,

ole

Dear kuechenrule,

About the PIMPLE methods and the rho, this is the easiest way of dealing with it. With the PIMPLE or PISO iterative method you try to converge the equations, eventhough the densities may vary por each equation. Using small time-steps you could reduce the problem, because the smaller the time step the smaller the variations in density.

Moreover, I would now try to eliminate the density variation with the liquid and solid, just try to use the same density for them, as it may induce some problems. Furthermore, I would also use the same cp for both liquid and solid, so that you get rid of problems that can arise from using the interpolation methods.

Next, I would try to use another approach to the convMeltFoam. As you can see convMeltFoam taked into account the convective transport term for the liquidFraction, or "fl". Although this might work when there is no gas, and thus, when you are not using interFoam, once you add this term to interFoam, you are dealing with high gradients over the surfaces and thus, poblems whith the fvc::div term.

At last I would try to relax the TEqn, just adding TEqn.relax() after the TEqn.solve() code.

If you don't mind, I could have a look at your complete interFoam code, and try to fix it.

Thank you.