CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Cavitation and Diverging Solution

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

Like Tree1Likes
  • 1 Post By Tobermory

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 27, 2024, 14:09
Default Cavitation and Diverging Solution
  #1
New Member
 
Ryan
Join Date: Feb 2024
Location: South Africa
Posts: 5
Rep Power: 2
RavingRaven is on a distinguished road
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:

Quote:
{
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.
RavingRaven is offline   Reply With Quote

Old   February 28, 2024, 06:04
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 670
Rep Power: 14
Tobermory will become famous soon enough
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 likes this.
Tobermory is offline   Reply With Quote

Old   February 28, 2024, 12:39
Default
  #3
New Member
 
Ryan
Join Date: Feb 2024
Location: South Africa
Posts: 5
Rep Power: 2
RavingRaven is on a distinguished road
Quote:
Originally Posted by Tobermory View Post
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:


Quote:
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:


Quote:
--> 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.
RavingRaven is offline   Reply With Quote

Old   February 29, 2024, 14:44
Default
  #4
Member
 
Shravan
Join Date: Mar 2017
Posts: 63
Rep Power: 9
Severus is on a distinguished road
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
Severus is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
mixture model for cavitation mech5190 FLUENT 12 December 7, 2016 03:24
Checking Size of Element in CFX or CFD Post dreamchaser CFX 26 November 6, 2014 01:02
Cavitation with CFX Eric CFX 0 January 10, 2006 09:53


All times are GMT -4. The time now is 12:30.