
[Sponsors] 
October 18, 2009, 13:05 

#21 
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
That was my initial approach but of course the result didnt converge. Now I dont know if the lack of convergence resulted from the incorrect approximation of the divergence term or the unsteady term or the diffusion term in the heat equation. I would have to do one at a time to see which actually caused the problem.
Again, reading through Edin's comments, he seems to argue that its not correct to apply the same linear interpolation to kappa as is done to density. So there are two issues here, one is the use of face centered values as against the cell centered values and the method of computing the "mixture" average: linear or harmonic interpolation. 

October 18, 2009, 15:54 

#22 
Member
Eelco Gehring
Join Date: Mar 2009
Posts: 70
Rep Power: 9 
I don't see why it's not possible to do a simple linear interpolation, since the calculation is based on gamma, which has already been calculated and corrected.
Edin, would you mind explaining why it is not possible to do the same interpolation for kappa as is done with rhoCp? Thanks, Eelco 

October 19, 2009, 03:51 

#23 
Senior Member
J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 9 
Hi, I meet a problem. That is I want to update the rho1 according to the liquid temperature, thus rho1 should be changed to volScalarField from dimensionedScalar. However, my C++ knowledge is limited, and now I have no idea about this. Would you help me? thanks.
Chiven part of the files of twoPhaseMixture.C and twoPhaseMixture.H are attached. for the file of twoPhaseMixture.C rho1_(nuModel1_>viscosityProperties().lookup("rho")), rho2_(nuModel2_>viscosityProperties().lookup("rho")), bool twoPhaseMixture::read() { if (transportModel::read()) { if ( nuModel1_().read(subDict(phase1Name_)) && nuModel2_().read(subDict(phase2Name_)) ) { nuModel1_>viscosityProperties().lookup("rho") >> rho1_; nuModel2_>viscosityProperties().lookup("rho") >> rho2_; return true; } else { return false; } } else { return false; } } for the file of twoPhaseMixture.C dimensionedScalar rho1_; dimensionedScalar rho2_; // Return constaccess to phase1 density const dimensionedScalar& rho1() const { return rho1_; } // Return constaccess to phase2 density const dimensionedScalar& rho2() const { return rho2_; }; 

February 8, 2010, 14:39 

#24 
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Hi Chiven,
Have you been able to implement the temperature dependent density in your code? Please let me know if you still need any help with that. Regards 

February 8, 2010, 19:36 

#25 
Senior Member
J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 9 
Hi Ovie, thank you, this time I am struggling for the interPhaseChangeFoam with condensation (vapour to water), I have finished the code implementation, but the calculated temperature fields are not right. I am thinking how to improve the codes.
How about your progress? Best regards, Chiven 

February 8, 2010, 20:29 

#26 
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Hi Chiven,
I am also working on phase change but evaporation. In my case the properties of my fluid are all dependent on temperature and another scalar variable (dry solids mass fraction) so I have another scalar transport equation in addition to the energy transport equation. Right now, I am having two challenges: 1.) I have implemented the temperature dependence of the various fluid transport properties but when I tried to validate the code on a single phase 2D channel flow the results dont converge. The temperature field looks reasonable but not "smooth" and same goes for velocity. I have included a thread on this at : Poor convergence with temperature dependent density on modified pisoFOAM I am still not quite sure how to fix this. 2.) The second issue is how to implement the production and annihilation source terms for the alpha1 equation for a temperature dependent phase change phenomena. I have also included a thread on this at: Source terms for alpha1 Equation for temeprature based interPhaseChangeFoam These are the tow challenges I am facing. I am sure though that you have gotten around implementing your viscosity models. If you have any ideas on how to proceed with my two questions, I would very much appreciate it. Regards and thanks... ovie 

May 5, 2010, 11:04 

#27 
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
Hej,
I am working on the same problem. I have a model for the phase change from vapor to water, already implemented and at that point i am struggling with the implementation for the temperature equation. The first idea was to implement the energy equation only on one side, since I am only interested on the water side at the moment. I am using saturated steam on subcooled water. Have zou guys come further with implementation for the energy equation in interPhaseChangeFoam?
__________________
~roman 

May 5, 2010, 12:08 

#28 
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Hey Romant,
If you are using the single field representation which the interFoam family of solvers is based on, then the energy equation must apply to the entire flow field which is both the liquid and vapour phases of the flow. So I am not quite sure what you meant by implementing the energy equation for only the liquid side as you indicated in your comment. The implementation of the energy equation is pretty straight forward just like any scalar associated with the flow. However, it must be expressed in conservative form. i.e. d(rho*cp*T)/dt + div(cp*T*rho*V) = div(kappa*grad T) + source terms The first term (energy stored) on the LHS of the equation can be rewritten as: d(rho*cp*T)/dt = d(rhoCp*T) where rhoCp = "rho*cp". rhoCp is for the mixture. To evaluate rhoCp, you would go to the alphaEqnSubCycle.H code. At the end, you find the line where rho is updated for the solver. Include a line that computes rhoCp by doing: rhoCp = alpha1*rho1*cp1 + (scalar(1)  alpha1)*rho2*cp2 The second term (energy flux) on the LHS is rewritten (based on OF terminology): div(cp*T*rho*V) = div(cp*T*rhophi). rhophi is the face mass flux. obtained from the alphaEqn.H code. Again, this can be written as: div(cp*T*rhophi) = div(rhophiCp*T) rhophi is updated for the solver in the alphaEqn.H code. To update rhophiCp, go the alphaEqn.H code and just after the line that computes rhophi, include the following line: rhoPhiCp = phiAlpha*(rho1*cp1  rho2*cp2) + phi*rho2*cp2; For the energy diffusion term i.e. first term on the RHS, you have to include a line that computes kappa in your code or read it off the transportProperties dictionary. The last term on the RHS are the source terms. These are the models which should describe energy consumption when liquid is converted to vapour (depending on the temperature and saturation pressure) and energy released when vapour reverts to liquid. So you would have something like source terms = (mass_of_vapour_to_liquid  mass_of_liquid_to_vapour) * latent heat of vapourization. This is how I have implemented my solver so far. However, I am yet to come up with appropriate models for the evaporation and condensation processes. Have you been able to get past that? I hope I have been able to help somewhat. Last edited by ovie; May 6, 2010 at 01:35. 

May 6, 2010, 03:59 

#29 
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 12 
Hej,
I think i am gonna try that as soon as possible as you described. at the moment I have considered saturated steam only and direct contact condensation, so there is no evaporation. the model that i am using is something like m = h * A (Tsat  Tbulk) / latentHeat where A is modelled based on the alpha per cell with Amax = sqrt(3) * Vcell^(2/3) and then the interface becomes Al = Amax * (2*alpha) for 0<alpha<0.5 and Al = Amax * 2 * (1alpha) for 0.5<alpha<1 for the heat transfer coefficient, h, I used a model based on steam falling onto a flat plate for the moment. of course this model is exchangeable, as well as the interface area. evaporation is set to 0. this is also the reason I also wanted to model the temperature on only the water side, because the other is supposed to stay constant.
__________________
~roman 

June 21, 2010, 07:04 

#30 
New Member
Join Date: May 2010
Posts: 27
Rep Power: 8 
Hey guys,
I am new to openfoam and c++ programming and I am trying to implement Ovie's temperature field code. As I am trying to compile it I keep receiving the following error message "metro@ubuntu:~/OpenFOAM/OpenFOAM1.6.x/applications/solvers/multiphase/interTempFoam$ wmake SOURCE=interTempFoam.C ; g++ m64 Dlinux64 DWM_DP Wall Wnostrictaliasing Wextra Wnounusedparameter Woldstylecast Wnonvirtualdtor O3 DNoRepository ftemplatedepth40 I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/transportModels I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/transportModels/incompressibleTemp/lnInclude I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/transportModels/interfaceProperties/lnInclude I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/turbulenceModels/incompressible/turbulenceModel I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude IlnInclude I. I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/OpenFOAM/lnInclude I/home/metro/OpenFOAM/OpenFOAM1.6.x/src/OSspecific/POSIX/lnInclude fPIC c $SOURCE o Make/linux64GccDPOpt/interTempFoam.o In file included from interTempFoam.C:58: createFields.H: In function ‘int main(int, char**)’: createFields.H:26: error: ‘alpha1’ was not declared in this scope createFields.H:26: error: ‘rho1’ was not declared in this scope createFields.H:26: error: ‘cp1’ was not declared in this scope createFields.H:26: error: ‘rho2’ was not declared in this scope createFields.H:26: error: ‘cp2’ was not declared in this scope createFields.H:42: error: ‘rhoPhi’ was not declared in this scope In file included from interTempFoam.C:85: TEqn.H:5: error: no matching function for call to ‘ddt(Foam::volScalarField&, <unresolved overloaded function type>)’ TEqn.H:6: error: no matching function for call to ‘div(Foam::surfaceScalarField&, <unresolved overloaded function type>)’ TEqn.H:7: error: no matching function for call to ‘laplacian(Foam::surfaceScalarField&, <unresolved overloaded function type>)’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:11: warning: unused variable ‘transonic’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:14: warning: unused variable ‘nOuterCorr’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:3: warning: unused variable ‘nCorr’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:8: warning: unused variable ‘momentumPredictor’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:11: warning: unused variable ‘transonic’ /home/metro/OpenFOAM/OpenFOAM1.6.x/src/finiteVolume/lnInclude/readPISOControls.H:14: warning: unused variable ‘nOuterCorr’ make: *** [Make/linux64GccDPOpt/interTempFoam.o] Error 1" Can anyone help me? Thanks in advance Regards Metro 

June 21, 2010, 12:31 

#31  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Quote:
You still need to make additional modifications to the codes for the interTempFoam solver like creating variables in the createFields.H. I am attaching a copy of my codes for your review to see if this would help. ovie Last edited by ovie; June 21, 2010 at 13:15. 

June 21, 2010, 13:17 

#32  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Quote:
In addition, rather than making hard coded changes to the twoPhaseMixture class, I implemented a derived class for the thermal properties which inherits the functions in the incompressibleTwoPhaseMixture class. This I called incompressibleTwoPhaseThermalMixture class. This directory should be in the $FOAM_SRC/transportModels/incompressible directory Dont forget to recompile this as well. ovie 

June 21, 2010, 16:17 

#33 
New Member
Join Date: May 2010
Posts: 27
Rep Power: 8 
Ovie,
This is great thank you for your response and sharing your files. Very much appreciated. I am just in the progress of reviewing it and finding what is wrong with mine. I noticed in the previous responses that you were considering adding mass transport due to evaporation and condensation to your interFoam solver. Did you have any success? 

June 21, 2010, 16:33 

#34  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Quote:
metro, I have not been able to implement that yet. I am still trying to benchmark the results for my falling film flow with experimental results. Are you working on something similar? ovie 

June 21, 2010, 18:03 

#35 
New Member
Join Date: May 2010
Posts: 27
Rep Power: 8 
Ovie,
Thank you again for attaching the additional files. As I said before I just need to review and compare it to my code. I am currently looking at condensation/evaporation in two phase horizontal pipes so in affect one can roughly look at it at a radial falling film problem. How accurate is our model compared to its benchmarks? Regards Metro 

June 22, 2010, 02:53 

#36  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
metro
No worries. We all need the help and the forum has truly been invaluable. Quote:
The part relating to phase change however is definitely similar so long as temperature is the driving force for evaporation/condensation. Quote:
1. D. Gao et al :Numerical simulation of wavy falling film flow using VOF method (Journal of Computational Physics) 2. Z. F Xu et al: Mass transfer across the falling film: Simulations and Experiments (Chemical Engineering Science) The results are pretty encouraging so far. 

June 28, 2010, 03:07 

#37 
New Member
Join Date: May 2010
Posts: 27
Rep Power: 8 
Hey Ovie,
I have implemented your code and I am now having stability issues. When I run my case, the temperature diverges and ends up reaching inf. I tested the case/boundary conditions with interFoam and no issues. When it came to implementing the div(rho*phi*cp,T) in fvScheme, I used a standard linear Gaussian discretisation. Is this correct? Or did you utilise another methodology? Can you think of any other reason why this is happening? Regards Metro 

June 28, 2010, 14:33 

#38  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Quote:
Quote:
Goodluck! 

June 28, 2010, 16:01 

#39 
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Metro,
please note that this solver is specifically for laminar flows. You have to make some modifications to extend to turbulent flows. Just needed to clarify this so u are using it for the right flow problem. rgds.. 

June 29, 2010, 02:57 

#40  
Member
Ovie Doro
Join Date: Jul 2009
Posts: 99
Rep Power: 8 
Quote:
Like I said in an earlier reply, if the instability you speak of is at the interfacial region then its almost a consequence of the sharp jump in material properties at the interface. In other words, if the difference in specific heat (cp) between the two fluids is significant then grid resolution might not solve the problem. Besides, since cp is such a huge value compared to the other fluid properties, it only serves to blow up the convective term in the discretized governing equation relative to the other terms. So the instability is persistent. However, you may modify the equations if certain conditions apply: if the specific heat of the fluids is not dependent on temperature or still weakly dependent or if the temperature differences within the flow (eg between the heated or cooled wall and the incoming fluid) is not "large" (note that you would have to evaluate the fluid properties at the mean temperature), you can move cp from the convective term to the diffusive term so that the governing equation to solve is d(rho*T)/dt + div(rho*V*T) = div((kappa/cp)*grad T) + source terms/cp This is a more stable problem. For your TEqn.H, you can do: volScalarField kappa = twoPhaseProperties.kappa(); volScalarField Cp = alpha1*cp1 + (scalar(1)  alpha1)*cp2; fvScalarMatrix TEqn ( fvm::ddt(rho, T) + fvm::div(rhoPhi, T)  fvm::laplacian((kappa/Cp), T) ); TEqn.solve(); rho and rhoPhi already come from alphaEqnSubCycle.H and alphaEqn.H respectively. In this case, you would no longer require rhoCp and rhoPhiCp. Hope this helps. Last edited by ovie; June 29, 2010 at 03:30. 

Tags 
interfoam energy, temperature field 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Moving mesh  Niklas Wikstrom (Wikstrom)  OpenFOAM Running, Solving & CFD  122  June 15, 2014 06:20 
Getting a concentration field around a bubble in InterFoam  azman  OpenFOAM Running, Solving & CFD  2  January 7, 2013 14:47 
Adding temperature field to InterFoam  yapalparvi  OpenFOAM Running, Solving & CFD  8  October 14, 2009 20:18 
Problem with rhoSimpleFoam  matteo_gautero  OpenFOAM Running, Solving & CFD  0  February 28, 2008 07:51 
Problems calculating field gh with interFoam  cricke  OpenFOAM Running, Solving & CFD  0  December 10, 2007 08:17 