Large inconsistency between versions in two phase flows with significant viscosity differences between phases

[DanielFractal]

When running simulations on incompressibleVoF when there are phases with large differences in viscosities, the simulations present instabilities and different results from v12.

I was able to pinpoint that this is due to the new way of computing the divergence of the stress tensor (in particular the transpose of the velocity gradient part). As an example with the simplest model (linear viscous stress), it used to be the case that the divergence would interpolate the volume field alpha*mu*grad u^T (as a multiplication of volume fields) :

const fvVectorMatrix divDevTauCorr
    (
        this->divDevTauCorr
        (
            -(this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))),
            U
        )
    );

Now, in v13, alpha*mu and the grad u^T fields are both interpolated to faces separately and then multiplied as surface fields:

const surfaceScalarField alphaRhoNuEff
    (
        fvc::interpolate(this->alpha_*rho*this->nuEff())
    );

    const fvVectorMatrix divDevTauCorr
    (
        this->divDevTauCorr
        (
           -alphaRhoNuEff
           *fvc::dotInterpolate(this->mesh().Sf(), dev2(T(fvc::grad(U)))),
            U
        )
    );

This idea of interpolating then multiplying (rather than multiplying then interpolating) causes a large difference when computing this term, proportional to viscosity and/or velocity gradient jumps on the interface. I am also pretty sure this is the cause as I have built OF with and without this interpolation of the mu scalar field.

I have attached a serial working example where the difference is visible (the case I am working on) of squeeze flow with a very large viscosity that does not present any problem in v12. The differences are also noticeable when the viscosity is not as large as well, but I use nu=2.0 for demonstrating purposes. The mesh is orthogonal and reasonably refined, and the implicit momentum predictor step does not even properly run in the new version due to this difference in results. When using either CICSAM or PLIC variants, the differences from the previous version are very large, as opposed to convective schemes like upwind, which I suspect is due to the large diffusion of interface that happens in the latter case.

simpleSqueeze2D_MWE.zip

Still, even after reading the given justification for commit https://github.com/OpenFOAM/OpenFOAM-13/commit/6592798ca0704aa6eaaaaf74ea3cc6a9553317dc#diff-49d145f386e7e961f9251be71e7b40baf74dfad71e1f9618f101af55c97a0570, I wonder if this change makes sense. Returning a surface field is not a problem, but the order of interpolation should matter.