[Kellis]

Good afternoon all,

As a precursor to a project involving phase change simulations I have been playing around with the cavitation models included in the compressibleInterFoam solver within OpenFOAM 10. My problem of interest involves predicting phase changes within a closed vessel as a response to changes in ambient temperature, and tracking the distribution of mass within the domain is important.

To establish a baseline before modifying the solver, I set up a case with a simple 10 cm cube which is initially filled half way with liquid water, and the remainder with water vapor. I set the initial pressure within the vessel below the saturation pressure to trigger vaporization of the liquid phase and then let the simulation run until it reaches a quasi-steady state. I am using the phaseChangeProperties dictionary and some solution/solver parameters from the ballValve tutorial case, which uses the Merkle cavitation model.

I am tracking the total mass within the domain using a coded function object of the following form:

dimensionedScalar mtot = gSum(alphaLiq*rhoLiq*v + (1.0 - alphaLiq)*rhoGas*v);

With my initial solution parameters, over 20% of the initial mass in the domain disappeared within the first 5s of simulation time, going from an initial total mass of 0.497 kg to a total mass of 0.391 kg (see attached figure). I have been trying to find a set of solution parameters which result in acceptable mass conservation but have not been able to achieve this despite a number of attempts. Starting from the baseline case, I've been modifying things one or two at a time to see what affects the solution accuracy and stability, with the following results:

1. Baseline case. 78% of initial mass remains after 5s simulation time.

The following two cases were copied from the baseline with modifications as follows:

2. Double cell count in all 3 directions. 83% of mass remains.
3. Increase from 3 to 10 outer correctors, and 0 to 2 non-orthogonal correctors. 93% of mass remains.

Each subequent case uses 10 outer correctors and 2 non-orthogonal correctors unless otherwise specified.

4. Increase from 2 to 5 non-orthogonal correctors. 93% of mass remains (no change).
5. Increase from 10 to 20 outer correctors. 93% of mass remains (no change).
6. Increase alpha correctors from 2 to 4 and alpha subcycles from 1 to 2. Unstable, negative temperature.
7. Decrease max Co and alpha Co from 0.25 to 0.1. Unstable, negative temperature.
8. Increase max Co and alpha Co from 0.25 to 0.9. 93% of mass remains (no change).
9. Change div(phi,p) from Gauss upwind to Gauss limitedLinear 1. 93% of mass remains (no change)
10. Change pcorr solver tolerance from 0.02 to 1e-4, increase alpha solver max iterations from 5 to 10, turn on momentum predictor. 93% of mass remains (no change).
11. Change initial timestep from 1e-4 to 1e-6s. Unstable, floating point exception.
12. Double cell count in all 3 directions. Unstable, negative temperature.

In summary, making the mesh finer produced a small benefit, and increasing the outer correctors to 10 gave a significant benefit, though increasing further did not yield better results. Trying to employ both a fine mesh and increased correctors produced an unstable simulation.

I am assuming that the cavitation model was validated before release to ensure mass conservation, but I've gone through what I feel like are the most obvious settings and can't seem to get the darn thing to work right. I've attached my baseline case if anyone wants to download things and have a look around - I would really appreciate any insight on this problem.

Cheers,
Kellis

[Kellis]

As a sanity check and to isolate the cavitation model, I removed the phaseChangeProperties file to prevent mass transfer. I then reran the case with the liquid initialized in the top half of the domain so it has to fall through the gas phase and settle on the bottom (initializing with the liquid on the bottom as before would give a trivial solution). With this setup and otherwise identical solution parameters, I get essentially zero mass change (as expected), so the issue definitely lies in the cavitation model.

Does anyone have a solution for this issue?