CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Bubble Radius in Schnerr Sauer cavitation model (mistake?) (https://www.cfd-online.com/Forums/openfoam-solving/128068-bubble-radius-schnerr-sauer-cavitation-model-mistake.html)

zarox January 2, 2014 05:28

Bubble Radius in Schnerr Sauer cavitation model (mistake?)
 
Hi every one,

Perhaps, I am mistaken or I have misunderstood something but I don't agree with the formulation of the reciprocal of the bubble radius.

Code:


  return pow
    (
        ((4*constant::mathematical::pi*n_)/3)
      *limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1),
        1.0/3.0
    )

Thanks to the Radius defintion : Rb = (alpha/(n*4/3*pi*(1-alpha)))^(1/3), The correct reciprocal radius should be :

Code:


  return pow
    (
        ((4*constant::mathematical::pi*n_)/3)
      *(1.0 + alphaNuc() - limitedAlpha1)/limitedAlpha1,
        1.0/3.0
    )

Perhaps, I forget something ...

What do you think?

zarox January 3, 2014 10:06

Definition is the key!
 
Hi everyone,

Ok, the answer is : it depends on the definition of alpha. If alpha is the volume fraction of liquid so the equation is correct... The formula given previously came from Simulating Cavitating Flows With LES in openfoam ECCOMAS CFD 2010 and alpha was the vapor volume fraction (as usual no ?), so I missunderstood the formula in the context. For OpenFoam, alpha1 definition depends on the user parameter definition, but the examples files usually have a definition of alpha1 with greater density than alpha2, so it stand for liquid.

All right!

shipman February 21, 2014 04:32

DEar Zarox,

Actually nowadays i am checking the code of interphasechangefoam with the schneerSauer trasport model. I am a bit confused about the some terms in the code. could you tell me please

1. when i look at the original paper of schner, he says: R=(3*(1-alpha1))/(4pi*n*alpha)^(1/3) it seems we dont need the alphaNuc() in the below equation? it is defined in the .h file as nucleation site volume fraction? what does it mean exactly? it means Vnuc/(1+Vnuc)??? and why we use here?

return pow
(
((4*constant::mathematical::pi*n_)/3)
*limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1),
1.0/3.0
)

2. In the below equation; in the original of the paper "p - pSat()" parameter should be inside of first sqrt like:sqrt(2*(p - pSat())/(3*rho1()))? then why it is put and multiplied with the mixture density rho? is it true?

return
(3*rho1()*rho2())*sqrt(2/(3*rho1()))
*rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
}

3. also in the same eqation what does it mean mag(p - pSat()) mag means magnitude? and why we add 0.01*pSat())) this term?really interesting?

4. again why we need to use alphaNuc() in the following equations. bcs i cant find in the original paper of schneer this term.

Cv_*(1.0 + alphaNuc() - limitedAlpha1)*pCoeff*min(p - pSat(), p0_)
);
(-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff

5. final question is related with the definitian of mDotP() and mDotAlphal() ? what do they correspond to exactly?

thanks in advance

Peter Müller February 21, 2014 11:41

Hi Zarox

I guess the alphaNuc() is just to be able to start cavitation. When everything is water at the beginning you will always have a multiplication by zero (from alpha_vapor) even when p < pSat. So you add a small portion of vapor in the water to not have this problem.

The reason for the sqrt(1/(p - pSat)) instead sqrt(p-pSat) is because you multiply by p - pSat in the solver so you obtain the same expression as Schnerr-Sauer. I wen't through all the equation and I obtained the formulation of Schnerr-Sauer, so to me it seams correct.

The reason for this 0.01*pSat I've no clue. I remember that I've seen it somewhere in a paper but I can not remember which it was.

Your final question goes hand in hand with sqrt(1/(p - pSat)). In order of generality the equations are formed in the most general way so you can switch the cavitation model and you always pas the same stuff to the pressure and alphaEquation. If you open these equations in interPhaseChangeFoam solver you will see what I mean.

Hope I could help you little bit.
Nice Weekend

shipman February 22, 2014 01:01

Hi Peter,

Thanks so much for your kind reply. yeah you are right The reason for the sqrt(1/(p - pSat)) instead sqrt(p-pSat) is because it is multiplied by p - pSat in the solver.

About one point i am still confused. In this definition: sqrt(mag(p - pSat()) + 0.01*pSat())); what does it mean "mag" exactly? it makes result pozitif afer this (p-PSat)?


And, in my cavitation simulation i use water whose Psat=2300kPa. According to this equation, if the local pressure p is negatif in the recirculation zone, this term sqrt(mag(p - pSat()) + 0.01*pSat())) will give always floating point error due to negative inside of sqrt, isnt it?

if you can let me know will be appreciated.

Thanks in advance.


Baris

Peter Müller February 23, 2014 11:26

Hi Boris

The function "mag" takes the calculates the magnitude of a given field, in this case a field of a scalar value, the pressure. It is like using "abs" which you probably know better.
It is just returning the size of the value and is therefore always positive.

To your second question: I think your pressure should never be lower than zero since for cavitation you are solving for the absolute value of the pressure and not as for incompressible calculations only the difference in pressure is relevant.
The strong density gradient at the interface region may lead to your negative pressure but physically this is wrong to my opinion.

Nice evening

shipman February 23, 2014 22:04

Hi Peter,

thank so much for your reply. I got it point now. After your answer, i checked the multiphase tutorials which i used before like (cavitatingFoam[compressible], interfoam, bubbleFoam, and interphasecahngeFoam[incompressible] all of them are using absolute pressure and it doesnt drop below zero even some of them are incompressible.

Could you tell me pls what you mean exactly with your this comment: ...and not as for incompressible calculations only the difference in pressure is relevant.

Last thing; do you think that my previous post http://www.cfd-online.com/Forums/ope...foam_help.html about interphasefoam which shows an error about the disappering of cavitation even though mass flow is increasing, is related with this pressure problem?

Thanks in advance.

Baris.

Peter Müller February 24, 2014 02:04

Hi Baris

When you run incompressible calculations the only effect on the flow due to pressure is based on the local pressure difference. So for example when you have a simple pipe case with one in and one outlet it doesn't matter whether you impose 0Pa or 10000Pa at the outlet if you have a zeroGradient for p at the inlet. The difference in pressure from in to outlet should be the same in the converged solution and this is what drives your flow.

To you second question: Having a look at your boundary specification I guess it is over constraint. I do not know what inletOutlet does, must switch from zeroGradient to fixedValue dependent on phi I guess. Imposing the pressure at in an outlet to me is a bit hard. I would impose the pressure only at the outlet. You should then modify your inlet velocity to obtain the desired pressure drop. I do not think that you can say it's an error if cavitation is disappeared when massflow increases since it's still only an effect of the pressure field and you should have a look at this. Perhaps at the beginning the flow was not converged and pressure was in cavitating conditions while at the end the local pressure was not low enough.

Nice Week

shipman February 24, 2014 03:16

hi
 
Dear Peter,

Thanks so much for your reply. So even in multiphase flow (for interphasechangefoam solver) which is also incompressible solver, i can set the outlet pressure 0 or 100000Pa (while inlet pressure is zerogradient) which will not change the results, isnt it?

So i set the my new BCs as follows, try it again. will let you know the results.

==> 0/U
Code:

internalField  uniform (0 3.104 0);
boundaryField
{
    inlet
    {
        type            fixedValue;
        value          $internalField;
    }
      outlet
    {
        type            zeroGradient;
    }
    wall
    {
        type            fixedValue;
        value            uniform (0 0 0);
    }

==> 0/P
Code:

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

Last thing is abut your this comment: as you said
HTML Code:

Perhaps at the beginning the flow was not converged and pressure was in cavitating conditions while at the end the local pressure was not low enough.
if the flow is not converged at the begining what i should do? Could you give to me so advices?

Thanks in advance

Here(in japan) already monday evening, so for you i think have a nice day Pete:)

Baris

Peter Müller February 24, 2014 04:42

Hei Baris

It's not completely correct for interPhaseChangeFoam. Clearly for momentum equation it should not make a difference but since you want to resolve cavitation it is important that you obtain the correct absolute value of pressure since this is needed to make the model work correctly. So even if interPhaseChangeFoam is basically incompressible just with two densities it's important to use absolute values for pressure and it should never be negative.

interPhaseChangeFoam is a transient solver which should converge in one inner iteration. So at the beginning if the solver is not able to converge in one time step you could lower the time step or you don't care for the moment since the transient phenomena you are interested are probably not the once during start-up of the simulation. You could also initialize the flow field perhaps by not solving for cavitation at the beginning until you think the fields are looking as you expect them:

nAlphaCorr 0;
nAlphaSubCycles 0;

Good luck

shipman February 26, 2014 22:26

Hi Peter,

Thanks so much for your advice. I am trying now and when i get the result will share with you.

I am just confused about one thing which you said: "So even if interPhaseChangeFoam is basically incompressible just with two densities it's important to use absolute values for pressure and it should never be negative"

Some days ago i have discussed this topic 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. Then, indeed if the absolute pressure never becomes (-);

1. so it means that i cant use interphasechangefoam solver for the liquids who have - pSat. isnt it? (bcs i already tried with negative pSat, it gave floating point error)

2. Also i read checked the book of Leighton, T., The acoustic bubble. Academic press, p. 67-129, 1994. He said indetical thing with you : "In reality, it is impossible to have a negative liquid pressure around an explosively expanding cavity for even a fraction of a second."

3. However, in the book of Brennen "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...


About this point if you can give be some insight i would like to be appreciated.

Thanks in advance.

Baris

isabel September 18, 2014 12:25

Dear everybody,

In the file SchnerrSauer.C, I think that this line extracts the value of alpha1 from the main code:

const volScalarField& alpha1 = alpha1_.db().lookupObject<volScalarField>("alpha1" ) ;

I do not understand if this other line is right:

const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");

I think the line above is wrong, the correct line would be this:

const volScalarField& p = p_.db().lookupObject<volScalarField>("p");

Am I right?

zarox September 24, 2014 08:05

object Registry
 
Hello Isabel,

Coming back from holidays, I showed your post. I don't have all the details explanation in mind. But roughly, it can be shown that alpha1_ is an inheritance from phaseChangeTwoPhaseMixture--> incompressibleTwoPhaseMixture--> two PhaseMixture. So this object is already available by the sub-class SchnerrSauer. The object alpha1_ is an volscalfield IOobject so contains object registery from mesh.
Code:

 
  alpha1_
    (
        IOobject
        (
            dict.found("phases") ? word("alpha" + phase1Name_) : alpha1Name,
            mesh.time().timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    ),

So, writing

Code:

const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
You call the register from alpha1_ [db], usefully you have already the alpha1_ object access from sub-class. And you look for the field "p".
Now, I hope it is helpful to understand why call the register from "p" can't be done : the p object is not available at first so we use the available data containning register for mesh data :: alpha1_.

I should say, at first look at these lines working to this solvers, I had the same reflex than you. But now the logic is quite understanble.

I hope it should help. Perphaps, I have once again miss something, but I am quite confident about the explanation. If you can think about that and look from your side I will be interresting in what you puzzle out.

Have a nice day.

Emeline


All times are GMT -4. The time now is 19:16.