Adding Crank-Nicolson off-centred flux correction to multiphaseInterFoam

Functionality to add/problem to solve

An explicit volumetric flux correction is added to multiphaseMixture [1] (used by multiphaseInterFoam). Currently, if time schemes with order > 1 (e.g. CrankNicolson, backward) are used with multiphaseInterFoam, the execution silently proceeds without the correction which brings about discrepancies in the results shown by the test case below.

In the proposed patch, a Crank-Nicolson off-centered volumetric flux correction is calculated and applied in the multiphaseMixture.C [1], equivalent to the implementation in alphaEqn.H [2] used by interFoam (and derivatives).

Analogous to interFoam, the patch implements the correction for CrankNicolson and disables the backward scheme by raising a FatalErrorInFunction.

Target audience

Users of multiphaseInterFoam interested in time schemes of order > 1, i.e. CrankNicolson.

Proposal

git patch multiphaseInterFoam-CN-Corr.patch

Changes

Main changes

• the main implementation consists of setting the off-centering coefficient according to the ddt scheme and explicitly correcting the fluxes for CrankNicolson. An error is raised for a time scheme other than Euler or CrankNicolson.

• it was necessary to make alphaPhiCorrs_ a member to allow use of alphaPhiCorr.oldTime(), which is needed for CN correction calculation.

Additional changes consistent with interFoam

• interface compression at non-coupled boundary faces is omitted
• restart handling as in interFoam (alphaPhiCorr fields are read if present - analogous to alphaPhi0 in interFoam)
• made naming consistent (alphaPhi -> alphaPhiCorr) in the final loop over the phases in method solveAlphas(const scalar cAlpha) in multiphaseMixture.C [1]

build

• run openfoam/applications/solvers/multiphase/multiphaseInterFoam/Allwmake

What does success look like, and how can we measure that?

Droplet test case

• droplet falling at fixed velocity, periodic boundary conditions (top/bottom).
• case files: test-multiphaseInterFoam-CN-correction.tar.gz
• 'before' result calculated on master (commit 97c8e62f)

Results

Notes

• Paraview state used (in case root) postProc.pvsm
• the difference can grow in time, especially for 'real-world' cases.

Links / references

[1] openfoam/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

[2] openfoam/applications/solvers/multiphase/VoF/alphaEqn.H