
[Sponsors] 
February 26, 2014, 06:14 

#21 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Dear Nakagawa,
Thanks so much for your reply. Could you tell me the key point that how can i understand the necessity of absolute pressure? In which condition i need to use it? On the other hand you said that "pSat is expressed in the absolute pressure". I checked the solver it is defined also as constant 2300Pa. So how can we say that in interphasechangeFoam pSat is defined as absolute pressure? And one more thing: Above post at #18, Bruno said that Pabsolute=Preference*rho+Patm, So could you tell me the exact equation of absolute pressure in interphasechangeFoam. And you said that in interfoam no neeed to use absolute pressure, but i checked the dimension of p in interfoam which is dimensions [1 1 2 0 0 0 0]; so again acc. to Bruno's post above at #18, absolute pressure is used in interfoam too? or Do i misunderstand something? Please could you explain clearly. Thanks in advance Baris 

February 26, 2014, 09:11 

#22  
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 285
Rep Power: 7 
Quote:
1. pressure in interPhaseChangeFoam is defined in pEqn as: Code:
p_rgh = p  rho*gh; The reason for this is that in phaseChangeTwoMixture we want to use pSat as defined in Pa, which is a property of the fluidvapour mixture that you are working with. This is also important as our densities will change depending on the volume fraction value in the cell (I suppose it would be possible to convert all of this into relative units but I don't think it's worth the extra effort). Hope that clarifies things a bit Peace A 

February 26, 2014, 09:11 

#23 
Member
shinji nakagawa
Join Date: Mar 2009
Location: Japan
Posts: 32
Rep Power: 9 
hi shipman,
Let's think physical meaning of the saturationvapourpressure. Its value is not relative. We can not choose arbitrary value. We have to pick the value from a steam table or Phase diagram. http://en.wikipedia.org/wiki/Water_(...#Phase_diagram if you want to simulate phase change, your solver should have something like phase diagram or table (or equations corresponding to them). So, you have to use the absolute pressure. Negative value is unphysical. The dimensions [1 1 2 0 0 0 0] means pressure itself. The dimensions [0 2 2 0 0 0 0] means that pressure is divided by density. If density is constant at everywhere in your simulation model, you can use [0 2 2 0 0 0 0] . This means the NS equation is divided by constant density. Please check Ueqn of the solver. If density is NOT constant at everywhere in your simulation model, you will use [1 1 2 0 0 0 0]. Because you can not simplify the NS eauation by dividing with density. In interFoam solver, density depends on alpha. So the dimensions of pressure is [1 1 2 0 0 0 0]. This is not related to absolute/relative pressure we are talking. So, you can use the relative pressure value with interFoam. Tutorials for interFoam use relative value, 0, isn't it? Compressible solver is always use [1 1 2 0 0 0 0] because density changes. In previous discussion, topics of p_rgh == p rho*gh, abusolute/relative (Prelative = Pabs  Patm; gauge pressure in engineering), and p/rho (two sets of dimensions) are mixed up, and we got confused. Maybe, I should say OpenFOAM is using the absolute pressure in every solver. We can use relative value when no properties depend on pressure. I appreciate any correction or modification of this comment. 

February 27, 2014, 05:01 
Hi

#24 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Hi Arthur and Nakagawa,
First, thanks so much for the elucidative informations. Now, i got the real meaning of absolute value which is used in openfoam. ==> Nakagawa, yeah i think that you are right. there was a misunderstading between reference pressure and definition of p/rho. i got it now. i will have 2 more question. please let me know your idea: 1. Some days ago i have discussed about absolute pressure with my prof. he said that inside of the nozzle in the recirculation zone sometimes vel. reaches ~100m/s and result in negative absolute pressure. and also He said that there are some liquids which have  negative pSat. But when i tried negative pSat inside of interphasechangefoam it gives floating point error. then i have checked 2 books: A. HTML Code:
. Leighton, T., The acoustic bubble. Academic press, p. 67129, 1994. He said indetical thing with you : B. HTML Code:
"Brennen, C. E., Cavitation and bubble dynamics, Oxford University Press on Demand, p. 4890, 1995." ==> So please let me know your ideas about absolute pressure whether it can be negative or not? 2. And, i have already posted about the my main problem using interphasechangefoam. http://www.cfdonline.com/Forums/ope...foam_help.html At the begining i see the cavitation formation at the recirlation zone, however, although mass flow is increasing (means vel. inside nozzle also increasing) cavitation region disappears soon. Then i saw that although velocity increases at the recirculation zone, i found out that pressure also increases... Any help to solve this problem will be very appreciated. thanks in advance... Best regards. Baris. 

February 27, 2014, 05:36 

#25 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 285
Rep Power: 7 
Hi,
I have had the same problem as you with the floating point error using SchnerrSauer model (Kunz was fine). I traced it down to the following line in SchnerrSauer.C: Code:
/( rho*sqrt( mag(p  pSat()) + 0.01*pSat() ) ); Kunz uses a slightly different formulation so you can set pSat to any value and as long as the relative difference will be the same you should be OK in terms of floating point errors (the behaviour might become unrealistic though, you'd need to look at the source terms in the transfer model's code to verify that). As far as the other part of your question goes, my understanding of the problem is quite biased towards Prof. Leighton's version but that might be because I am more familiar with it (I took lectures from him). I'm afraid I'm not qualified enough to make some bald statements on this topic. Still, from the mass transfer modelling point of view in OpenFOAM I noticed the following lines in both Kunz and SchnerrSauer models: Code:
SchnerrSauer: Cc_*limitedAlpha1*pCoeff*max(p  pSat(), p0_), Cv_*(1.0 + alphaNuc()  limitedAlpha1)*pCoeff*min(p  pSat(), p0_) Kunz: mcCoeff_*sqr(limitedAlpha1)*max(p  pSat(), p0_)/max(p  pSat(), 0.01*pSat()), mvCoeff_*min(p  pSat(), p0_) Code:
0_("0", pSat().dimensions(), 0.0) Peace, A 

February 27, 2014, 07:04 

#26 
Member
shinji nakagawa
Join Date: Mar 2009
Location: Japan
Posts: 32
Rep Power: 9 
hi,
About 1st question from shipman, I also vote to Prof. Leighton. I'm not saying Prof.Brennen or your supervisor is wrong. I have never read the book you sited. We should be careful about the context. Using gauge pressure is quite common in engineering. Absolute pressure is 0 at perfect vacuum. This is the lowest pressure. http://en.wikipedia.org/wiki/Vacuum#Measurement 

February 27, 2014, 07:24 

#27 
Member
shinji nakagawa
Join Date: Mar 2009
Location: Japan
Posts: 32
Rep Power: 9 
Hi Artur,
thank you for sharing information. Could you explain why you think the term " mag(p  pSat()) + 0.01*pSat() " could be negative? What I think is; The term "mag(ppSat())" is positive or zero. pSat is positive. (unless you give negative value in dictionary) Therefore "mag(p  pSat()) + 0.01*pSat()" should be positive. pSat is given from dictionary, not calculated in SchnerrSauer. Is this correct? I have never used phasechange solver itself and am not familiar with the models. But I'm interested in the phasechange solver. 

February 27, 2014, 07:35 

#28  
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 285
Rep Power: 7 
Quote:
EDIT: I had a closer look at the behaviour of the source terms in the SchnerrSauer model depending on the relative values of p and pSat and confirmed what I've written above: their values will not change as long as the relative difference between the p and pSat remain the same and the FP error will occur when (ppSat) tends to 0. Have a look at the attached Python script if you are interested. Last edited by Artur; February 27, 2014 at 10:38. 

February 27, 2014, 22:16 
Hi

#29 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Hi to all,
Thanks so much for your replies. 1. For just last confirmation: The term "mag(ppSat())" is always positive or zero. pSat is positive. (unless it is given negative value in dictionary) and if the relative difference ppsat tends to 0 while pSat is negative then we can get Floating point error. We all agreed about this point right? 2. ==>to Arthur: I have checked your written code, and following term will not give any floating point error bcs of your absolute setting for pSat, even though it is set to negative or relative difference ppsat tends to 0. Code:
(rho*np.sqrt(np.abs(ppSat_) + 0.01*np.abs(pSat_))) # changed here Code:
Cc_a = Cc_ * alpha1 * pCoff * max([ppSat_,p0_]) Cv_a = Cv_ * (1.+alphaNuc()alpha1) * pCoff*min([ppSat_,p0_]) 3. Finally, about the my main problem, when i set the intiial conditions as attached directly i obtain cavitation inside nozzle then after maybe 250 iteration later it is disappeared. Just i am confusing about one point. When i checked the velocity inside nozzle which is increasing, but at the same point pressure also increasing. I dont know why? Also i used kwSST turbulence model. if we need further data please let me know to solve my problem. thanks Baris. 

February 28, 2014, 05:11 

#30 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 285
Rep Power: 7 
Just change the control parameters at the top of the file to something like:
Code:
# =========== # CHANGE THESE p = 1001. pSat_ = 1000 # =========== I'll have a look at your case later today when I find some time 

February 28, 2014, 12:27 

#31 
Senior Member
Artur
Join Date: May 2013
Location: Southampton, UK
Posts: 285
Rep Power: 7 
OK, about your case setup:
1. I don't think you need to use the following with zeroGradient on alpha1 outlet: Code:
inletValue $internalField; 2. I noticed you don't have the nut BC's in your 0 folder  from my experience it is useful to have just for the sake of specifying a nutWallFunction the way you want it to be (again,, this probably won't change much in terms of your problems) 3. I noticed you have fixedValue BC on pressure and inletOutlet on U. From what I know usually one should use pressureInletOutletVelocity when specifying a fixedVAlue pressure otulet (see the cavitatingBullet tutorial, for instance: Code:
U: outlet { type pressureInletOutletVelocity; phi phi; value $internalField; } p_rgh: outlet { type fixedValue; value $internalField; } 4. Usually when I use multiphase solvers I specify fixedFluxPressure or buoyantPressure BCs on noslip walls. The reason for doing so is to keep the flux consistent with the U BC (looking at multiphaseFixedFluxPressureFvPatchScalarField.H: "This boundary condition adjusts the pressure gradient such that the flux on the boundary is that specified by the velocity boundary condition") Again, see the cavitating bullet tutorial for an example: Code:
bullet { type fixedFluxPressure; } Peace, A 

March 1, 2014, 06:35 
hi

#32 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Hi Arthur,
First thanks so much for your reply. As a final, i set the my condition as follow; 1. 0/alhpa1 similar with tutorial Code:
internalField uniform 1; boundaryField { inlet { type fixedValue; value $internalField; } outlet { type inletOutlet; inletValue $internalField; } wall { type zeroGradient; } Code:
internalField uniform 0; boundaryField { inlet { type calculated; value uniform 0; } outlet { type calculated; value uniform 0; } wall { type nutkWallFunction; value uniform 0; } Code:
internalField uniform (0 3.104 0); boundaryField { inlet { type fixedValue; value $internalField; } outlet { type pressureInletOutletVelocity; phi phi; value $internalField; } wall { type fixedValue; value uniform (0 0 0); } Code:
internalField uniform 1e5; boundaryField { inlet { type zeroGradient; } outlet { type fixedValue; value $internalField; } wall { type fixedFluxPressure; } Finally, based on your last comment, i never tried initialising by using like simplefoam etc...Could you give me some advice how can i do it? You mean that run my case by setting according to simplefoam, and use final data as initial date for interphasefoam?? if like this how can i understand the final stable result in simpleFoam? Sorry never tried before so no idea. thanks so much baris 

March 2, 2014, 00:00 
hi

#33 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Hi Arthur,
According above setting i tried my case, unfortunately, it gives same results. At the begining of the calculation, the cavitation takes place at the entrance of nozzle, then extend until same degree then disappeared. This is my main problem which i tried to solve it for almost 2 months. When i sent you my file bcs of the over size i couldnt send the my geometry. now i am sending it too. When you have time could you look at please. https://www.dropbox.com/s/b7e8s5xx5n...onstant.tar.gz Thanks in advance Baris. 

March 2, 2014, 16:28 

#34  
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,651
Blog Entries: 39
Rep Power: 99 
Greetings to all!
@Baris: Quote:
The description of the solver interPhaseChangeFoam is as follows: Quote:
Best regards, Bruno
__________________


March 3, 2014, 01:48 

#35 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Hi Bruno,
Yeah my dimensions are correct as you said. My nozzle thickness is 1.94mm, inlet width is ~8mm and outlet is 1.94mm. Actually, i think that i mistyped something about words of my Prof. He just said that absolute pressure can be negative value in the case of real type nozzle throat, velocity reaches the order of 100 m/s. However, for my case (nozzle) throat mean velocity reaches max 25 m/s. And, according to bernouilli i checked, never became negative. I have another question. Maybe, problem is my case is complex and need to initiliase the my flow field. I never tried this but i read some post which is said that use potential flw to initiliaze the velosity field. Do u think that after that is it better to use simpleFoam to initilize the pressure, k and omega ? Another question is that after initializing to convert the pressure case should i multiple by rho all the fields value?? and in the incompressible case pressure can be calculated negative. during converting these value for interphasechangefoam should i use directly as negativee?? If you can help will be appreciated. Thanks. Baris 

March 4, 2014, 11:05 

#36 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,651
Blog Entries: 39
Rep Power: 99 
Hi Baris,
My advice is to initialize the case with the same solver, but with an inlet velocity 100 to 10000 times smaller than the one you currently have. Then use changeDictionary after the case has somewhat converged, to change the inlet speed and then continue the simulation. The tutorial "multiphase/interFoam/ras/damBreakPorousBaffle" seems to briefly demonstrate how to use changeDictionary. You can find more examples by running: Code:
find $FOAM_TUTORIALS $FOAM_UTILITIES name changeDictionaryDict Best regards, Bruno
__________________


March 7, 2014, 06:26 
hi

#37 
Senior Member
Baris (Heewa)
Join Date: Jan 2013
Location: Japan
Posts: 130
Rep Power: 5 
Dear Bruno,
Indeed, thanks so much for your advice. I am still trying and i think it is going fine until now. After completion the result, i will post here. n addition, i have 2 questions: 1. According to your above advice: what do you mean exactly by saying HTML Code:
you should do a study of results vs inlet speed, to assess at which point the solver is no longer able to properly simulate the case 2. Next question is to change the kinematic eddy viscosity equation (nut) in kepsilon turbulence model which is defined: Code:
nut_ = Cmu_*sqr(k_)/epsilon_; Code:
nut_ = f(rho)*Cmu_*sqr(k_)/epsilon_; f(rho) = rho1+ (rho2rho)/[(rho2rho1)]^5 Could you help to me please. Thanks in advance. 

March 23, 2014, 13:43 

#38  
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,651
Blog Entries: 39
Rep Power: 99 
Hi Baris,
Quote:
You could also try with a ramp inlet: ramp inlet velocity initial condition using timeVaryingMappedFixedValue Quote:
And I have the strange feeling that this question has either being asked in the past or that OpenFOAM already has this kind of calculation somewhere. Good luck! Best regards, Bruno
__________________


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
GroovyBC the dynamic cousin of funkySetFields that lives on the suburb of the mesh  gschaider  OpenFOAM  300  October 29, 2014 19:00 
c++ libraries and solver compiling  vaina74  OpenFOAM Installation  13  February 3, 2012 18:43 
Saving ParaFoam views and case  sail  OpenFOAM Paraview & paraFoam  9  November 25, 2011 16:46 
DecomposePar links against liblamso0 with OpenMPI  jens_klostermann  OpenFOAM Bugs  11  June 28, 2007 17:51 
user defined function  cfduser  CFX  0  April 29, 2006 10:58 