Hi Fabian,
Thank you for the effort which you have put into the meltFoam solver. I am interested in using your solver in the development of a thermal storage model. In your paper you indicate the choice of constants C and b in the porosity function for the meltFoam test case based on the work of Gau and Viskanta. In your own experiment, you follow the work of Shmueli et al. in guide the selection of the large constant C. I intend to model the location of the phase front in a melting Silicon body. Could you please point me in the right direction regarding the selection of the constants C and b? Brent and Voller et at. (1988) state "Further work needs to be carried out in the area of mushy region phase change to establish guiding principles in assigning appropriate values for both the constants C and b". ...I am hoping there has been some development in this area since 1988! |
Hi Liam
Unfortunately, there are not much more references I can recommend to you. Shmueli et al. performed a parametric study on the large constant C. However, they modeled a VOF gas phase above the melt. This may lead to wrong results. In my present study I found out, that the large constant is not the only parameter to vary. Next to the permeability (which is the large constant C), the viscosity of the melt in the transition region has to be taken into account. Moreover, by using the liquid fraction as porosity function, porosity is changing linear. This is not necessarily the case for all PCM. As you use a pure PCM (silicon) you should choose a large C for the Darcy source term. You could also use the switch-of technique instead of the Darcy term. Soon I will publish a paper on all this stuff concerning the Darcy type source term. Also keep in mind that my solver is for non-isothermal phase change of mixtures. Regards Fabian |
Hi Fabian,
Thanks for getting back to me on this. Your advice has been helpful so far - I'll be sure to share my findings with you. Just to keep you in the loop, I also emailed Vaughan Voller with a similar question and this was his response: With a pure material where one would physically expect a sharp interface between the solid and liquid the choice of the parameters in the porosity model are not critical--provided they are sufficient to "switch of" the liquid velocity in the computational cell undergoing the phase change So in a first case I would use the Gallium melting values. But there would no harm in seeing what happens when you adjust these values--I would not expect adjustment up to 1/2 and order to have much impact--but worth and easy to check. Thanks again for your help. Liam Montgomery |
Parametric study for Darcy cosntant
Hi Liam
Thanks that you share your knew findings with us. Vaughan Voller is for sure wright when he says, that the C has no big impact when it comes to isothermal melting. A higher C does not change much. Is C to low, you will see your solid PCM sinking downwards. I would suggest a parametric study varying C from a low value lime 100 to say 10e8. This will give you a feeling. Apropos, when you use different densities for solid and liquid, C will become more and more important. Regards and thanks again for the discussion Fabian |
Thanks for all your contributions!
I just downloaded the erfMelt piece and the meltFOAM piece. Will try to use this to model lava!
anyone has any kind of documentation for these? no worries if not. Just figured it will save me some time if someone already wrote something. |
I just downloaded the erfMelt piece and the meltFOAM piece. Will try to use this to model lava!
anyone has any kind of documentation for these? no worries if not. Just figured it will save me some time if someone already wrote something. |
No documentation jet
Hi
I did not write any documentation for the solver, since it is not 100 percent finished. The erfMeltSolver reduces the non-linearity effects in the energy conservation equation. However, it does not entirely get rid of it. For simple problems, the results converge with a very small error for the energy conservation. With increasing convective transport, the error increases and some iteration to account for the non-linearity have to be performed. F. Rösler, D. Brüggemann (2011): Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments. Heat and Mass Transfer, Vol. 47 Issue 8 , 1027-1033, DOI: 10.1007/s00231-011-0866-9 http://www.springerlink.com/content/b1tp01k2u7q8j432/ In my previous work, I use a linear liquid fraction function and the linear source method proposed by Voller. Regards Fabian |
Quote:
if I have time I ll rebuild the solver for OpenFOAM-2.2.x At the Moment I rebuild the flamelet solver and build a transient solver for combustion based on the flamelet thermodynamic model. After that (and please write me a message if I dont Keep you updated ) I try to rebuild the solver for 2.2.x Tobi |
problem with liquid fraction
Hi,
I'm working with the meltIcoFoam solver which is basically cloned from the meltFoam. In my cases, I have seen the liquid fraction at the boundary gives some unusual value. Is there anyone who has worked with the same code of meltFoam? It will be great if anyone can upload validated case for me for the meltFoam. Cheers! Quote:
|
Hello Mohammad,
could you post your solver and test case? Without more information, helping you is a difficult task. Regards, Anja |
melting problem
1 Attachment(s)
The solver I'm using is exactly the meltIcoFoam "http://www.cfd-online.com/Forums/openfoam/93620-melting-problem.html", and the test case is as attached.
|
Hello Mohammad,
I do not see anything in your set up for the test case that might force the solver to crash. As I do not know, how you changed the solver to work with a newer version of OpenFOAM than 1.4 (which I guess is the one posted in #2 of this thread) I can only guess.
Regards, Anja |
boundary condition of liquid fraction
Hi
I did not see any boundary condition for the liquid fraction in your test case. It should be zeroGradient everywhere except for inlets. Check the erfMeltFoamSolver in the thread. You do not necessarily need the boundary condition as the liquid fraction (alpha) is a function of temperature but then it has to be handled a little bit different in the solver code. All fields depending on the liquid fraction, so basically all thermo physical properties, should use the boundary condition of the liquid fraction. Code:
volScalarField xyz Fabian |
Hello Fabian,
thanks for the reply. Again, I learned something new. Regards, Anja |
Temperature and liquid fraction update with Voller method
Hi,
Has anyone already implemented the linear source method proposed by Voller ? Will it be possible to post code for that (Just the iteration loop between temperature and liquid fraction)? Cheers shakil |
fl-T update cycle
Hello Shakil,
here is the frame of the way I do it. For the hEqn.h file, it is: Code:
int iter=0; The additional solving parameters iterMax and convergence are declared as follows in a file like readAdditional.h. Code:
dictionary flUpdate = mesh.solutionDict().subDict("flUpdate"); Code:
flUpdate Have a nice day, Regards, Anja |
Hi Fabian,
I'm back onto using your solver to model a melting problem. Essentially I am trying to model the melting of a converging prism, with a heat flux applied to the top (eg: http://imgur.com/RPsaWHa) The melting occurs as expected in the x-z plane with convection visibly effecting the melting front as shown here: http://imgur.com/iEl6kfu However, since your solver only appears to work in 2D (one of the assumptions in shell-and-tube paper), the melting is not correctly modelled in the x-y plane: http://imgur.com/9J6bAbM How much effort would be involved in extending the solver to 3D if at all possible? Would I be better off using 2D results and trying to extrapolate to 3D? Thanks again for your help. Kind regards, Liam Montgomery |
3D melting problem
Hi
the solver is fully capable of solving 3D cases :). The assumption in the paper is only stated to reduce cell number and so the calculation time. However, I have to state again, that the energy conservation in my error function solver is only valid for small time steps as one would actually need an iterative correction procedure like Voller is using it. Regards Fabian |
Ok, thanks - that's great news. Do you have any suggestions as to why convection is not appearing the same in both the x-z and x-y planes as shown in the images?
|
Hi everyone!
I'm going to enter the fray of writing a melting/solidification solver too... :) That's why I'd like to ask for some heads-up on the solver for OF2.1 which Fabian kindly shared with us. I understand that it's based upon the enthalpy-porosity technique by Voller et al and the fluid flow is modelled using a PIMPLE solution. I'd like to ask though what general idea the simulation of the flow is based upon. Is it a VOF approach, an Euler-Euler approach or even something different? |
History excurse
Hi
Well it is the enthalpy method. You use one eulerian grid for both, fluid and solid PCM. Oleinik [1] proposed a weak formulation of the Stefan problem in the 1960s. Based on this approach, Rose [2] and Solomon [3] developed a method based on enthalpy. This is the mother of all enthalpy methods developed since. [1] O. A. Oleinik: A method of solution of the general Stefan problem, Soviet Mathematics, vol. 1, p. 1350–1354, 1960. [2] M. E. Rose: A Method for Calculating Solutions of Parabolic Equations with a Free Boundary, Mathematics of Computation, vol. 14, p. 249–256, 1960. [3] A. D. Solomon: Some Remarks on the Stefan Problem, Mathematics of Computation, vol. 20, p. 347–360, 1966. Hope this is what you were searching for. :) Regards Fabian |
Thanks for your answer... :)
I guess I was a bit confused with the model and things like the buoyancy term and thus thought that there was some kind of multiphase-flow going on. As the model simplifies things into a one-phase flow there is nothing like that though... ;) |
Hi everybody,
I would like to update the Rosler's solver to the version 2.2.0 of O.F. Do you know what kind of modifications I have to made for that? Thanking you in advance. Eric |
Hi,
I have an stupid question about post processing. I can plot the alpha (phase fraction) contour in paraview. But I'm wondering how to plot the graph (and export) for comparing as it is not possible in Paraview. |
Hello Mohammad,
well, first of all, there are no stupid questions :-). A way I use is only to enable alpha to be shown > in Volume Fields. Then I do a "contour" like you with only one value, 0.5 for example. Split the view and click "Spreadsheet View" for the new one. For the spreadsheet view, only the eye in front of "contour" should be enabled. With the spreadsheet view still active, do "Save Data" and choose "csv" as file type, for example. You will obtain a column file with x,y,z,alpha (not necessarily in that order) that you can plot with gnuplot or other tools to compare the results. If your cases do not vary in geometry, you could compare even a distorted contour plane in 3D using paraview. Only open both cases, create the contour for each case, colour them differently and let them be shown (active eye") in a single 3D view. I hope, this does help you Regards, Anja |
Thanks AnjaMiehe.
|
Voller Source based enthalpy method
1 Attachment(s)
Hi,
I'm trying to update the solid-liquid fraction as mentioned by Voller source based enthalpy method. First, for checking the case I have used just heat conduction. The way I'm trying to solve the energy equation is: volScalarField gamma0 = gamma; //-storing the prev time value do { gamma.storePrevIter(); volScalarField dFdT = 1/(Tl-Ts); volScalarField Finv = gamma*(Tl-Ts)+Ts; volScalarField Sp = -((rho*L)/runTime.deltaTValue())*dFdT; volScalarField Sc = -Sp*Finv + ((rho*L)/runTime.deltaTValue())*(gamma0-gamma); fvScalarMatrix TEqn ( fvm::ddt(rhoCp, T) -fvm::laplacian(K,T) - fvm::Sp(Sp, T) = Sc ); TEqn.solve(); gamma = gamma + dFdT*(T-Finv); gamma = gamma = min(max(gamma,scalar(0)),scalar(1)); Cp = gamma*Cpl+(1-gamma)*Cps; K = gamma*Kl +(1-gamma)*Ks; //- residual calculation of g }while (condition for convergence of g) This the main code for energy equation, and implemented exactly the same case of Gallium error melting just for Test. I'm having the PROBLEM as, PROBLEM (1) Temperature is not updating with time advancing (2) I have switch off (Sp and Sc )to check from where the problem is coming, I saw it's coming from Sp and Sc. (3) I'm wondering whether the loop has any problem ??? I have attached the solver..... PLEASE ADVICE ME TO COME OUT OF THE PROBLEM |
Hello Mohammad,
I don't have time right now to look into it, i will try it tomorrow. However, there are some things that make me think about it already now: - Could you provide the equations you based your code on here. I don't have the paper at hand an Voller wrote lots of them ;-) - for dFdt I suppose a linear slope of the fraction of solid, right? But even thought, this is not 1/(Tl-Ts) everywhere, only within the mushy zone. Thus: Code:
dFdt=pos(Tl-T)*pos(T-Ts)/(Tl-Ts) - The old time value is used gamma.oldTime(), no need to store it in an additional variable So much for now, Regards, Anja |
1 Attachment(s)
Hi AnjaMiehe,
I have attached the equation that I am solving. The paper that I am working on is: http://www.tandfonline.com/doi/abs/1...2#.Up63WkAW1hg Attachment 27158 Quote:
|
2 Attachment(s)
Hello Mohammad,
to begin with, I am not quite sure if you will like my solution. I have always understood OpenFOAM to code in using grad, div, laplacian ... and not trying to put the discrete equation together with a_p ... Using the Voller paper and your equations I did the following: with Knowing that the source term S can be put together in the following way The equation to be coded is with The equation is coded in the following way Code:
fvScalarMatrix TEqn Code:
gamma = (T-Ts)/(Tl-Ts); Sorry, that I cannot assist you with the exact implementation of the Voller paper. If you want to move on to more complex correlations of gamma to temperature, you could try to read in an interpolation table. See this link for more information: http://www.cfd-online.com/Forums/ope...ies-false.html Regards, Anja |
Well I'm currently writing my masters thesis about this topic and I got something running that works (although has to be revamped at certain places to include additional features). My problem is about the flow and solidification of an alloy though but as I'm right now changing the temperature calculations in order to correctly take care of the eutectic solidification I can post how I'm doing things right now. My solver's design also is based on Voller's enthalpy-porosity-model but I'm not really copying a certain paper 1:1...
But oh well here is the code of my TEqn.H file which calculates the new temperature field: Code:
// initialize and solve energy equation Code:
forAll(fracSol, celli) |
Hi,
I have implemented the method that I mentioned. It's working now and giving the correct value. Now I have some things to dig out: (1)Anyone please suggest me, how can I find out the interface velocity ? (I want to calculate the interface velocity in the solver ) (2) Have anyone worked with the non equillibrium solidification considering the under cooling ? cheers |
Quote:
|
Hi,
I'm little bit of puzzle of a term. Please help me to correct the definition: In the energy equation: d/dt(rho*Cp T) + div(rho U Cp T) = laplacian(K T) + ST, I am a bit of concern about the source term, whether it be... (1) ST = rho*L*d/dt(gamma), ONLY THE LATENT HEAT PRODUCT OR (2) ST = dH*d/dt(gamma), where dH = rho*(Cl-Cs)*(T-Tref)+rho* L. Where, Cl = liquid specific heat and Cs is the solid specific heat. |
Question about Fabian Test Case (Gallium)
Hi,
I have a query about the Fabian Gallium test case. I have seen that, in the fvSolution file, the solver has been chosen for p_rgh and T , but not for U. Can any one please explain the reason behind it ? in fvSolution file: Code:
solvers |
Hello Mohammad,
do please take a look at posts No #31,#32 of this thread. Fabian is giving an idea about it there. Regards, Anja |
melting/solidification with temperature dependent thermo-properties
2 Attachment(s)
Hi Anja Miehe,
Sorry it's me again. In my solver, I want to account for temperature dependent thermo-physical properties during melting/solidification. The governing equation that I'm solving is similar to Fabian and as follows, and the source term is updated as iterative procedure as stated by Voller, like Sp and Sc. And at each iteration I'm updating the liquid fraction. Until now it works fine, and I have also validated the case with Gallium melting. The result is good. But the problem I'm facing now is to work with the temperature dependent thermo-physical properties. The code diverges after 2/3 time steps. I'm sure, the way I am working with the temperature dependent properties are not correct. Please give me some advice to deal with the temperature dependent thermo-physical properties in the case of melting/solidification problem. **My solver and case file is attached. Thanks in advance |
Hello Mohammad,
no need to be sorry, I am glad when I can help. Good to see that you found a way to code the Voller scheme like you intended. Well done. I did not try to compile your code as I am already in my Christmas Holidays and have only my little home laptop with me. Nonetheless, there are some things striking me.
Code:
relaxationFactors I hope, this will do it already. Regards and happy Xmas, Anja |
solid/liquid phase change solver convMeltFoam
1 Attachment(s)
Hi all
First of all, sorry for backing out from the very nice discussion in this thread in the last few months. I was packing together my thesis. But now it is officially submitted and two months of vacation are over. I appreciate that the solid/liquid phase change topic is still strongly represented in this forum and I am happy that Anja and Mohammad shared their solvers with us. I went through all of your posts and there is the topics of the iterative correction of the liquid fraction that should be discussed. The method used in Anjas Solver (Post 70) is iterating the energy conservation equation but is not correcting the liquid fraction. I recommend this early Voller paper to understand the necessity of the iterative update: Quote:
(1) is discretized as follows: (2) As gamma is calculated from temperature, this algebraic equation is non linear. Voller linearized it using an iterative method: (3) This equation can now be solved with any code (like OpenFOAM). However, as the source term was linearized, temperature,liquid fraction and enthalpy don not necessarily have to correlate. For this purpose, within the iterative procedure, a correction of the liquid fraction has to be performed. Have a look at the attached picture. When enthalpy is conserved, the following is correct: or (4) Substituting the enthalpy temperature correlation (5) into (4) leads to the equation to calculate the corrected liquid fraction: or (6) The steeper the slope of the enthalpy-temperature during phase change, the more unstable gets this method. Thus, equation (6) has to be under-relaxed leading to (7) with omega in the range of 0.6 to 0.8. However, one can consider the whole term in front of the brackets as under-relaxation factor which is much smaller. During iteration, the bracket will turn to zero as (8) So this is the iterative correction of the liquid fraction. However, this method will lead to nonsense when the energy conservation equation is not solved at least once per iteration step. In standard OpenFOAM, one can not give a minimum number of iterations in the fvSolution dictionary. So there are two possibilities to still achieve at least one iteration:
Never the less I would like to share the basic solver on which the solver from my thesis is founded as well as the necessary changes in the ldu matrix solvers to use the code in an efficient way. Extract the data inside the convMeltFoam.tar.gz and run the install script. This will provide the changes inside the ldu matrix solvers, compile OpenFOAM src to make the changes and install the solver. However, be careful as this script will only run with OF 2.2.x. Cheers :) |
All times are GMT -4. The time now is 19:13. |