[andrew-c]

I am trying to simulate hydrostatic boiling in a 3D domain. It involves two tall, connected columns of liquid water, which boil above a certain height due to a reduction in pressure. At the bottom of one column is a pressure outlet, and a pressure inlet at the bottom of the other (see attached image of mesh). Gauge pressure is specified at both (lower than 1 atm, to simulate the columns being above ground level). I'm using a symmetry condition at the midplane to reduce cost.

I am using Eulerian multiphase (implicit), realizable k-e, and species model. Mass transfer mechanism is species mass transfer, with a piecewise linear saturation pressure. The domain is 3D with gravity enabled, and is fairly large - about 2 m in z, and about 0.3 m in x and y. Mesh is fairly fine (I think) - largest size is 1e-3 m, smallest size is 2e-5 m.

I am trying to run in transient until the water in the columns boils down to a certain level - should take about 20-30 s (see attached image for the farthest I've gotten). The water level will correspond to the saturation pressure at a given temperature (right now it's all at 300 K). For now, I'm running with the energy equation off just to get flow convergence first, and plan to turn it on afterwards. However, I can't seem to achieve convergence.

The residuals are high (e-1 to e+1) - especially for continuity and vf-water. In all cases, these two residuals don't drop below e-3 for every time step, which seems like a red flag for convergence. Most cases eventually diverge, even with energy off. The few that don't diverge with energy off will diverge very quickly after enabling energy.

The best result I've gotten so far has converged for the first ~10-100 timesteps, but then the residuals start to steadily increase afterwards - screenshot is attached below. This eventually leads to divergence.

I believe the species model is the root of my issues. When I run without it (or even just with a different mass transfer mechanism, such as evaporation-condensation), it works much better. The trouble is, I need the species model, because I am eventually hoping to simulate diffusion due to a gradient in saturation pressure between the two tubes.

I've listed everything I've tried so far below. I would appreciate any suggestions to help achieve convergence. I'm still pretty new to Fluent, so anything helps. Thanks so much!!

Troubleshooting steps I've tried:

Mesh: I've refined the mesh down to the current level (max sizing 1 mm). My system can't handle much finer without it taking forever each iteration. Min orthogonal quality is 0.26, max skewness is 0.8, average is 0.23. I've tried both tetrahedral and polyhedral. I've also tried "improve mesh" in Fluent, improving ~1-10% of the worst cells (default settings).

Solver: The flow itself is minimal and slow, so I've tried laminar and turbulent. Laminar diverges very quickly when energy is turned on - not sure if this has to do with species model/thermal diffusion. For turbulence, I've tried k-e and k-omega, with just about all potential options (realizable, RNG, standard, enhanced wall functions, etc.). I've tried implicit and explicit Eulerian. I've also tried mixture model instead of Eulerian with no luck.

Solution methods: I've tried first order, second order, SIMPLE and coupled. I've tried some with energy off at first, and others with energy on from the start. I've tried timesteps from 0.1 s all the way down to 1e-8 s. I've increased max iterations/timestep to promote convergence (up to 50-100 in most cases, but usually I still don't see convergence each timestep). I've tried reducing URFs (both slightly and dramatically), but it only seems to stave off divergence to a later point.

BCs: All walls are at 0 heat flux (I've tried const. temperatures too). Inlet and outlet are both pressure BCs (0 gauge pressures).

Multiphase options: I've tried all mass transfer mechanisms under Raoult's law (ranz-marshall, hughmark, etc.). I'm using interphase heat transfer in the two-resistance model (other models don't seem to help). Omitting interphase HT leads to very rapid divergence when energy is turned on.

Other options: I've tried dynamic mesh and adaptive timestepping, both with no success. Both standard and hybrid initialization produce the same results. FMG initialization isn't compatible with multiphase. I don't have any URFs, and material properties are all constant (besides saturation pressure).