[shipman]

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

[Artur]

Originally Posted by shipman

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?

So could you tell me the exact equation of absolute pressure in interphasechangeFoam.

Hi, going bottom-to-top on your list:

1. pressure in interPhaseChangeFoam is defined in pEqn as:

p_rgh = p - rho*gh;

And has the units of [1 -1 -2 0 0] or [kg m-1 s-2] or Pa in other words (hence we know its absolute and not relative to the density as in isimpleFoam, for instance ([0 2 -2 0 0 0 0]). Here p is the total pressure that you compare with your pSat values (from Bernoulli you could say that H = 0.5 rho U^2 + p).

The reason for this is that in phaseChangeTwoMixture we want to use pSat as defined in Pa, which is a property of the fluid-vapour 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

[snak]

hi shipman,

Let's think physical meaning of the saturation-vapour-pressure. 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 un-physical.



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.

[shipman]

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.

. Leighton, T., The acoustic bubble. Academic press, p. 67-129, 1994. He  said indetical thing with you :

" He stated that In reality, it is impossible to have a negative liquid pressure around an explosively expanding cavity for even a fraction of a second." ==> So, this judgement is same with most of people said in forum.

B.

"Brennen, C. E., Cavitation and bubble dynamics, Oxford University Press on Demand, p. 48-90, 1995."

he said that it should be noted that pressure inside the bubble is always positive even though outside oressure around of bubble could be negative for a fraction of a second. ==> this is totally opposite of leighton.

==> 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.cfd-online.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.

[Artur]

Hi,

I have had the same problem as you with the floating point error using Schnerr-Sauer model (Kunz was fine). I traced it down to the following line in SchnerrSauer.C:

/( rho*sqrt( mag(p - pSat()) + 0.01*pSat() ) );

If your (p-pSat) term tends to zero the value inside the square root becomes -ve and you get the error.

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 Schnerr-Sauer models:

Schnerr-Sauer:
        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_)

Apparently, depending on the direction of the mass transfer (or whether we're talking about Cc or Cv coefficient) the driving term (p-pSat) may be either +ve or -ve (since p0_ is defined as 0:

0_("0", pSat().dimensions(), 0.0)

I'd be very happy to hear someone more knowledgeable address the second issue you have raised too. Still, I think you will need to provide a bit more details about your case.

Peace,

A

[snak]

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

[snak]

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(p-pSat())" 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 phase-change solver itself and am not familiar with the models. But I'm interested in the phase-change solver.

[Artur]

Originally Posted by snak

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(p-pSat())" 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 phase-change solver itself and am not familiar with the models. But I'm interested in the phase-change solver.

Yes, this is exactly as you say. But shipman mentioned that the FP error occurred when he set pSat to a negative value in the dictionary specifically. Hence this is the only place where the equations may fall over because of doing so (I did a quick test and confirmed this by using mag(0.001*pSat) and setting it as < 0 -> it will no longer give a FP error).

EDIT: I had a closer look at the behaviour of the source terms in the Schnerr-Sauer 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 (p-pSat) tends to 0. Have a look at the attached Python script if you are interested.

[shipman]

Hi to all,

Thanks so much for your replies.

1. For just last confirmation:
The term "mag(p-pSat())" is always positive or zero.
pSat is positive. (unless it is given negative value in dictionary) and if the relative difference p-psat 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 p-psat tends to 0.

(rho*np.sqrt(np.abs(p-pSat_) + 0.01*np.abs(pSat_))) # changed here

I think that if we set the pSat to negative then another problem will arise related with the absolute pressure whether it can be negative or not. Because when we look at the final source terms equations (taken from your code): lets say if pSat is set to -68000Pa

Cc_a = Cc_ * alpha1 * pCoff * max([p-pSat_,p0_])
Cv_a = Cv_ * (1.+alphaNuc()-alpha1) * pCoff*min([p-pSat_,p0_])

In order to start the vaporization Pabsolute should drop to below of the psat (-68000Pa). However i have never obtained neg. absolute pressure value in interphasechangefoam. What you think about this point?

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.

[Artur]

Just change the control parameters at the top of the file to something like:

# ===========
 # CHANGE THESE
 p = -1001.
 pSat_ = -1000
 # ===========

And you will see NaN square root terms (as I said, the difference between the two must tend to zero).

I'll have a look at your case later today when I find some time

[Artur]

OK, about your case setup:

1. I don't think you need to use the following with zeroGradient on alpha1 outlet:

inletValue      $internalField;

but that probably won't change anything since the keyword is not used (I think)

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:

U:
    outlet
    {
        type            pressureInletOutletVelocity;
        phi             phi;
        value           $internalField;
    }

p_rgh:
    outlet
    {
        type            fixedValue;
        value           $internalField;
    }

See a nice discussion in this thread: http://www.cfd-online.com/Forums/ope...tvelocity.html

4. Usually when I use multiphase solvers I specify fixedFluxPressure or buoyantPressure BCs on no-slip 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:

    bullet
    {
        type            fixedFluxPressure;
    }

None of the above seem radically wrong to me, they are just different to what I use. Perhaps give it a go with the above mentioned changes and see if that helps. Also, in my cases it helps a lot when I initialise a case from a steady-state solution with no cavitation (just simpleFoam), the stability is much better. Do you ever do that too?

Peace,

A

[shipman]

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

internalField   uniform 1;
boundaryField
{
      inlet
    {
       type            fixedValue;
        value           $internalField;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
    }
    wall
    {
        type            zeroGradient;
    }

2. Actually nut file already existed in the 0 folder but during compress it gives oversize file error so i deleted it. it is set : 0/nut:

internalField   uniform 0;
boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
    wall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

3. Acc. to your advice set the 0/U as follows:

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);
    }

4. for p_rgh:

internalField   uniform 1e5;
boundaryField
{
   inlet
    {
         type           zeroGradient;
    }
    outlet
    {
        type            fixedValue;
         value           $internalField;
    }
    wall
    {
        type           fixedFluxPressure;
    }

I already started the case after setting like this, and will let you know result.

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

[shipman]

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.

[wyldckat]

Greetings to all!

@Baris:

Originally Posted by shipman

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.

I've taken a quick look into your geometry and boundary conditions... are you certain the dimensions of the mesh are correct? I ask this because:

1. The inlet is 0.00794 m by 0.00194 m.
2. The outlet is 0.00194 m by 0.00194 m.
3. The intake speed is 3.1 m/s.

Taking into account the comment by your professor, I would say that you're trying to use a solver that is not meant for simulations near the sonic region!


The description of the solver interPhaseChangeFoam is as follows:

Solver for 2 incompressible, isothermal immiscible fluids with phase-change (e.g. cavitation). Uses a VOF (volume of fluid) phase-fraction based interface capturing approach.


The momentum and other fluid properties are of the "mixture" and a single momentum equation is solved.
The set of phase-change models provided are designed to simulate cavitation but other mechanisms of phase-change are supported within this solver framework.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

I believe that you need a solver that implements formulation that is meant for (near) sonic flow. The only thing that comes to mind is to attempt to merge the solvers cavitatingFoam and sonicFoam, but that's not an easy task

Best regards,
Bruno

[shipman]

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

[wyldckat]

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:

find $FOAM_TUTORIALS $FOAM_UTILITIES -name changeDictionaryDict

Honestly, 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. I say this because 25 m/s for such a small geometry is still a pretty high flow speed.

Best regards,
Bruno

[shipman]

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

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

actually before using changeDictionary case, first i set the intlet velocity 1000 times smaller than real case(Vreal=3.1m/s ---- so i set it 0.0031) then i checked the velocity field until it gives proper distribution, then change the inlet velocity.

2. Next question is to change the kinematic eddy viscosity equation (nut) in kepsilon turbulence model which is defined:

nut_ = Cmu_*sqr(k_)/epsilon_;

i want to change as:

 nut_ = f(rho)*Cmu_*sqr(k_)/epsilon_;

f(rho) = rho1+ (rho2-rho)/[(rho2-rho1)]^5

i want to change it by multiply a function f based on the mixture density (rho) which is used inside of interphasechangefoam. rho1 and rho2 are constant so there is no problem. However, i just dont know how can i define the function f(rho) which can use the calculated mixture density at each time step in the turbulence equation.

Could you help to me please.

Thanks in advance.

[wyldckat]

Hi Baris,

Originally Posted by shipman

1. According to your above advice: what do you mean exactly by saying

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

actually before using changeDictionary case, first i set the intlet velocity 1000 times smaller than real case(Vreal=3.1m/s ---- so i set it 0.0031) then i checked the velocity field until it gives proper distribution, then change the inlet velocity.

OK, that's the idea I was trying to indicate. Start with slow flow speed at the inlet and then gradually try with greater flow speeds, one change per run.
You could also try with a ramp inlet: http://www.cfd-online.com/Forums/ope...ixedvalue.html

Originally Posted by shipman

2. Next question is to change the kinematic eddy viscosity equation (nut) in kepsilon turbulence model which is defined:

nut_ = Cmu_*sqr(k_)/epsilon_;

i want to change as:

 nut_ = f(rho)*Cmu_*sqr(k_)/epsilon_;

f(rho) = rho1+ (rho2-rho)/[(rho2-rho1)]^5

i want to change it by multiply a function f based on the mixture density (rho) which is used inside of interphasechangefoam. rho1 and rho2 are constant so there is no problem. However, i just dont know how can i define the function f(rho) which can use the calculated mixture density at each time step in the turbulence equation.

That's a programming question and should be asked in the respective forum Namely: http://www.cfd-online.com/Forums/ope...g-development/

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