# Diverging result for Temperature field in interFoam

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

 October 12, 2009, 10:14 Diverging result for Temperature field in interFoam #1 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Hi all.. I added an energy equation to interFoam, applied to a falling film flow on a cylindrical surface. The result shows a temperature field growing unbounded. I have no idea what could be wrong...maybe the implementation itself or my choice of scheme. Has anyone tried to add an energy equatioin to interFoam? Are there any peculiar problems I should be aware of? Please help... Thanks

 October 13, 2009, 06:19 #2 Member   Edin Berberovic Join Date: Mar 2009 Posts: 31 Rep Power: 16 How do you calculate the face fluxes in the temperature equation? They should be consistent with the fluxes calculated for the momentum equation.

 October 13, 2009, 11:53 #3 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Hi thanks for the response. The energy equation I added is: fvm::ddt(rho*cp,T) + fvm::div(rhoPhi*cpf,T) - fvm::laplacian(kappaf,T) rhoPhi comes from the momentum equation, cpf and kappaf are the cell face values of cp (specific heat) and kappa (thermal conductivity) respectively. First I calculated face values of gamma (vol fraction of phase 1) as gammaf : surfaceScalarField gammaf = min(max(fvc::interpolate(gamma), scalar(0)), scalar(1)); then I used this value to compute face values for cp and kappa which I called cpf and kappaf as follows: surfaceScalarFiled cpf = gammaf*cp1 + (scalar(1) - gammaf)*cp2 surfaceScalarFiled kappaf = gammaf*kappa1 + (scalar(1) - gammaf)*kappa2. Something tells me that my implementation is suspect though. I am not just quite sure... Thanks. Kummi likes this.

 October 14, 2009, 03:54 #4 Member   Edin Berberovic Join Date: Mar 2009 Posts: 31 Rep Power: 16 It depends on how you calculate your rhoPhi. It is updated after the solution of the gammaEqn. You should update your energy flux (rhof*Cpf*phi) to be consistent with this volume flux, i.e. not doing a simple linear interpolation of Cp. Another hint: when you are dealing with different materials (having different diffusion coefficients for energy), the you must apply harmonic face interpolation)

 October 14, 2009, 10:58 #5 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Thanks for the reply Edin. Regarding rhoPhi, I use the values directly from the updates given by the gammaEqn with no modifications whatsoever. Now two questions, 1.) How do I update the values for Cp without doing a simple linear interpolation? Should I employ a hihger order scheme for this? 2.) Please could you explain how harmonic face interpolation is implemented? Is there a dedicated scheme for this in OpenFOAM? Thanks.

 October 14, 2009, 20:03 #6 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Hi Edin, I have reviewed the results closely and its apparent that the instability results from the sharp change in thermal properties at the interface. I noticed that the temperature changes sharply at the region where gamma changes from its liquid value to its value for the gas phase. And as the computation marches in time, the instability propagates within the computation domain. So it may come down to a handling the change in thermal properties at the interphase. Does anyone have a clue? Thanks.

 October 14, 2009, 21:27 #8 Member   Eelco Gehring Join Date: Mar 2009 Posts: 70 Rep Power: 16 Ovie, I think you need to work in Cp in rhoPhi, it needs to be consistent with the volume flux calculate in gammaEqn.H. By that I mean, Cp needs to be taken into account at that point in the calculation. Chiven, why do you incorporate "turbulence->nut()/Pr"? I wouldn't expect to see this in a solver like interFoam. Why don't you define a kappa for both fluids?

 October 14, 2009, 21:51 #9 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Thanks Feijooos. However, I am not quite sure I understand what you meant by working in Cp in rhoPhi. Is there some other way to compute Cp for the "mixture" other than interpolating using the values of gamma and the constant cp values for the individual phases? Or are you suggesting that cp somehow has to be incoporated in the gammaEqn when solving for gamma? Please be kind to clarify. Chiven, thanks for the insight. I am trying to review the suggested code. It looks reasonable though. Thanks all..

 October 14, 2009, 21:52 #10 Senior Member   J. Cai Join Date: Apr 2009 Posts: 180 Rep Power: 16 Dear feijooos, thank you for the comments. Maybe you are right, and i shall accept your suggestions, and do some calculations. best regards, Chiven

October 14, 2009, 21:57
#11
Senior Member

J. Cai
Join Date: Apr 2009
Posts: 180
Rep Power: 16
Quote:
 Originally Posted by ovie Chiven, thanks for the insight. I am trying to review the suggested code. It looks reasonable though. Thanks all..

Dear Ovie, my revision of the codes have NOT been validated. I am doing some calculations, but it is calculating very slowly. I am not sure whether it is OK. If I have a new progress, I shall post it. Good luck.
Chiven

 October 15, 2009, 03:22 #12 Member   Edin Berberovic Join Date: Mar 2009 Posts: 31 Rep Power: 16 @ovie 1) do not update Cp alone, you don' need it. simply update the whole energy flux as well as rho*Cp instead, using the same way as mass flux and rho are updated. 2) you may simply use Gauss harmonic for the diffusion term, or implement it on your own like 1.0/fvc::interpolate(1/lambda).

 October 15, 2009, 09:59 #13 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 @ Edin, Thanks for the insight. I would review the suggested approach and see how it impacts the convergence of my result.

 October 17, 2009, 08:57 #15 Member   Edin Berberovic Join Date: Mar 2009 Posts: 31 Rep Power: 16 ovie, this is correct, congratulations. you are still doing linear interpolation of kappa, but I guess its influence is very small. enjoy it.

 October 17, 2009, 10:59 #16 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Hi Edin, Thanks for the review. But a quick one: which of these two expressions do you suppose would be the more appropriate expression for calculating kappaf in the function kappaf() in twoPhaseMixture.C? alpha1f*rho1_*cp1_*(1/fvc::interpolate(1/((1/Pr1_)*(nuModel1_->nu())))) + (scalar(1) - alpha1f)*rho2_*cp2_*(1/fvc::interpolate(1/((1/Pr2_)*(nuModel2_->nu())))) or alpha1f*(1/fvc::interpolate(1/((rho1_*cp1_*(1/Pr1_)*(nuModel1_->nu()))))) + (scalar(1) - alpha1f)*(1/fvc::interpolate(1/((rho2_*cp2_*(1/Pr2_)*(nuModel2_->nu()))))). Or is either inappropriate? Thanks.

 October 18, 2009, 09:05 #17 Member   Edin Berberovic Join Date: Mar 2009 Posts: 31 Rep Power: 16 hi ovie. maybe it would be a little complicated to write it all in one expression, because face values must be obtained from cell center values of kappa which should already be some kind of averages between phases. one way would be to update kappa() at cell centers also as a weighted average, just like mu() is updated, and then evaluate the conductivity at faces simply as 1.0/fvc::interpolate(1/kappa()). the other way would be to update kappa() at cell centers as a harmonic mean, but then interpolate it linearly at faces. i think this is something you may play with for a while, you can test it on purely diffusion problems. however, i am thinking whether the viscosity should be also evaluated using the same procedure. updating kappa() as weighted average is consistent with the diffusion coefficient for the UEqn, i.e. mu(), but the value at cell faces should definitely be evaluated using harmonic interpolation. this was shown e.g. in patankar(1980). it is clear that density is a weighted average, but i am not happy with the same definition for viscosity (and generally diffusion coefficients, like conductivity). nevertheless, as i see that all people are doing it exactly in the same way like it is in OpenFOAM, so i think the eventual errors are really small. regards. Elham and amolrajan like this.

 October 18, 2009, 12:47 #18 Member   Eelco Gehring Join Date: Mar 2009 Posts: 70 Rep Power: 16 Ovie, I was reading over this thread again and I saw that I missed your question on how to incorporate cp into the face fluxes. I guess you already found that out I have a question for you though, why do you use face values for kappa? Why not just use cellcentered values? Thanks, Eelco

 October 18, 2009, 13:41 #19 Member   Ovie Doro Join Date: Jul 2009 Posts: 99 Rep Power: 16 Hi Eelco, Those were my thoughts precisely: to use cell-centered values for kappa in the heat equation. I dont have any logical explaination for my choice of face centered values except that in the momentum equation, I noticed that face centered values were used for mu. And I simply followed suit. I dont know how this affects the accurary and maybe convergence of the results though.

 October 18, 2009, 13:49 #20 Member   Eelco Gehring Join Date: Mar 2009 Posts: 70 Rep Power: 16 hmm interesting. I'm wondering if you can just use: kappa == gamma*kappa1 + (scalar(1) - gamma)*kappa2; Don't see how this would affect the solution, but you need to specify the different conductivities for both fluids.

 Tags interfoam energy, temperature field

 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post azman OpenFOAM Running, Solving & CFD 3 June 7, 2022 05:21 Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 07:20 yapalparvi OpenFOAM Running, Solving & CFD 8 October 14, 2009 21:18 matteo_gautero OpenFOAM Running, Solving & CFD 0 February 28, 2008 07:51 cricke OpenFOAM Running, Solving & CFD 0 December 10, 2007 08:17

All times are GMT -4. The time now is 18:51.