[RavingRaven]

Hi everyone,


I'm a new member here - hope to meet a few of you!



I am busy with an interesting project in my free time - modelling the acoustic excitation (sonification) of bubbles undergoing cavitation. In a nutshell, soundwaves cause pressure fluctuations which cause bubble growth and Rayleigh collapse.

The image below illustrates the concept well:



You have a sinusoidal driving pressure (such as from a piston), at some frequency or amplitude, causing pressure fluctuations in the liquid medium. The bubble reacts by growing, then collapsing, then growing, then collapsing, before ultimately collapsing in a shockwave.

The mesh and geometry look like this:

domain.jpg
Attachment 98650
domain3.jpg

- blockMesh is used
- 2D and axisymmetric
- Only a single bubble is modelled
- The centre of the bubble sits at the origin, and quarter-symmetry is assumed (i.e one quarter of the bubble's area is modelled)
- At the centre is a fine radial mesh, while closer to the opposite boundaries it is more rectangular
- The bubble initial radius is R0 = 3e-6m. The radial/circular mesh extends up to 20 * R0, while the entire domain is (200 * R0 by 200 * R0).

The specific boundary and initial conditions I am trying to model are as follows:

- Two phases (air and water)
- The gas phase is ideal (I model it as perfectGas) and the liquid phase follows the Tait EOS (I use adiabaticPerfectFluid)
- Initial pressures (both p and p_rgh) of 148KPa and 100 KPa for the air and water, respectively
- Uniform temperature of 300K throughout the domain

The boundary conditions are as follows:

boundaryconditions.jpg

Where the wall that provides the driving pressure, or acoustic excitation, is uses the following code:

{
type codedMixed;
refValue uniform 0;
refGradient uniform 0;
valueFraction uniform 1;
source uniform 0;
value uniform 0;
name codedPatchBC_movingWall;
code #{
const scalar t = this->db().time().value();
const scalar pi = constant::mathematical:i;

// U = A * s * sin(2 * pi * f * t) * t + l
// f -> frequency 25000;
// A -> amplitude 75000;
// s -> scale 1;
// l -> level 100000;

this->refValue() = 100000 - 75000 * 1 * sin(2 * pi * 25000 * t);
this->refGrad() = 0;
this->valueFraction() = 1;
#};
}

I am currently using the compressibleInterFoam solver (I also tried compressibleInterIsoFoam).

The field is initialised using setFieldsDict, where I use a sphere to set the pressure (p and p_rgh) and alpha.water fraction:

bubble_alpha.jpg

I have been struggling with solver divergence and crashing. The most common crash is from a floating point exception, where there is a negative absolute temperature. Bizarre pressure and temperature discontinuities and singularities appear at the bubble edges:

problem1.png

After adjusting the timestep and Courant number and not getting any results, I felt that perhaps the problem was the setFieldsDict entry, where maybe the alpha.water boundary between the water and air had jagged nooks and 'corners' where the singularities formed. And so I used a setExprFieldsDict instead and used an arctan function to initialise the field with a slightly smoother transition between the water and air.

This had some success, and the singularities no longer appear en masse like that. However, I still have pressure and temperature divergence and I just cannot seem to create a stable solution.

Does anyone see anything obvious that I am doing wrong? Would you recommend modelling the driving pressure using a coded BC? Does the approach I am taking make sense? Is the quarter-symmetry and axisymmetry assumption a sensible one?

Any suggestions or advice would be greatly appreciated.

[Tobermory]

That sounds like a really interesting problem Ryan. It looks like the pressure solver is struggling - do you get acceptable convergence on the p equation before the solution starts blowing up, or is it hitting the maxIter limits? Would additional pressure correctors or PIMPLE iterations help?

I don't have any experience with compressibleInterFoam but here are a few other suggestions: try make the p solver work harder and see if that keeps the solution under control; add a limiter on the T field with the limitT function object in fvOptions to prevent unphysical temperatures; add in some non-orthogonal correctors, if you have not already done so.

Good luck and let us know how you get on!

[RavingRaven]

Originally Posted by Tobermory

That sounds like a really interesting problem Ryan. It looks like the pressure solver is struggling - do you get acceptable convergence on the p equation before the solution starts blowing up, or is it hitting the maxIter limits? Would additional pressure correctors or PIMPLE iterations help?

I don't have any experience with compressibleInterFoam but here are a few other suggestions: try make the p solver work harder and see if that keeps the solution under control; add a limiter on the T field with the limitT function object in fvOptions to prevent unphysical temperatures; add in some non-orthogonal correctors, if you have not already done so.

Good luck and let us know how you get on!

Hi Tobermory,


Thank you for the response! In line with your advice, I went to the fvSolution file and increased the number of correctors for the PIMPLE solver:

PIMPLE
{
momentumPredictor no;
transonic no;
nOuterCorrectors 3;
nCorrectors 2;
nNonOrthogonalCorrectors 2;
}

And in doing so increased the number of iterations.


While there was a minor increase in stability, unfortunately, the problem persists:

--> FOAM FATAL ERROR: (openfoam-2306)
Negative initial temperature T0: -2534.49

From Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::adiabaticPerfectFluid<Foa m::specie> >; Type = Foam::sensibleInternalEnergy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::adi abaticPerfectFluid<Foam::specie> >, Foam::sensibleInternalEnergy>]
in file ./src/thermophysicalModels/specie/lnInclude/thermoI.H at line 57.

FOAM aborting

What is odd is that the p_rgh residuals are not high. In fact, the number of iterations goes to zero, with tiny residuals on the order of e-16. However, this may be because the timestep also drastically plummets down to the e-12 range and lower (at this scale, "normal" is in the region of e-8 or so).

Again, strange pressure singularities appear:

pressure_problem.jpg


I even hacked the fvSolution file to turn off the energy equation implicitly (by simply setting the solver for T, h and e to have maxIter 0), but this made little difference.


All in all, it just seems to be a persistent instability issue, where pressure singularities form either at the bubble-liquid interface or near the mesh boundaries. I doubt it is the mesh, because checkMesh only returns "OKs".



So, I am afraid I am still at a loss.

[Severus]

Hi,
I have the following suggestions,
1) You could try to increase nAlphaSubCycles in fvSolution.

2) I see that you increased the outercorr to 3, you could try a case with a larger number of outercorr.

3) We sometime get a similar crash in multiphaseEulerFoam (Euler-Euler solvers) when the outlet is quite close to the free surface of the bubble column, so the outlet needs to be kept at a larger distance from the free surface

4) Your initial radius seems really small. I am not sure if that is the reason for your solver crashing, but I would suspect it. So you initially patch a bubble radius of 3e-6 using setFields? I would suggest you to check if the solver doesn't crash if you keep increasing the bubble size, then we can see if that is the reason.



Thanks