CFD Online Discussion Forums - Modifying pEqn: "cannot scale a matrix containing a faceFluxCorrection"

This post is related to the two-phase flow solver isoAdvector in OpenFoam v1706.

I am working to implement a new method to capture surface tension and hence, I am currently playing around with pEqn.H. Here is the pressure matrix I am trying to create.

Code:

fvScalarMatrix p_rghEqn

(   

    isLiquid * fvm::laplacian(rAULf, p_rgh)

    + (1 - isLiquid) * fvm::laplacian(rAUVf, p_rgh)

    ==  

     - isLiquid * fvc::laplacian(rAULf * DPField, isInterfaceOwnerLiquid)

     - (1 - isLiquid) * fvc::laplacian(rAUVf * DPField, isInterfaceOwnerLiquid)

     + isLiquid * fvc::laplacian(rAULf * DPField, isInterfaceOwnerVapor)

     + (1 - isLiquid) * fvc::laplacian(rAUVf * DPField, isInterfaceOwnerVapor)

     + fvc::div(phiHbyA) 

);

"isLiquid" is a volScalarField

Code:

volScalarField isLiquid

(

    IOobject

    (   

        "isLiquid",

        runTime.timeName(),

        mesh,

        IOobject::NO_READ,

        IOobject::AUTO_WRITE

    ),  

    mesh,

    dimensionedScalar("zero",dimless,0)

);

On executing this code, it throws the following error. I have attached a copy of the output log as well.

Code:

--> FOAM FATAL ERROR: 

cannot scale a matrix containing a faceFluxCorrection



    From function void Foam::fvMatrix<Type>::operator*=(const Internal&) [with Type = double; Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::Internal = Foam::DimensionedField<double, Foam::volMesh>]

    in file /usr/local/openfoam/v1706/OpenFOAM-v1706/src/finiteVolume/lnInclude/fvMatrix.C at line 1242.



FOAM aborting



#0  Foam::error::printStack(Foam::Ostream&) in "/usr/local/openfoam/v1706/OpenFOAM-v1706/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so"

#1  Foam::error::abort() in "/usr/local/openfoam/v1706/OpenFOAM-v1706/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so"

#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at /usr/local/openfoam/v1706/OpenFOAM-v1706/src/OpenFOAM/lnInclude/errorManip.H:85

#3  Foam::fvMatrix<double>::operator*=(Foam::DimensionedField<double, Foam::volMesh> const&) at /usr/local/openfoam/v1706/OpenFOAM-v1706/src/finiteVolume/lnInclude/fvMatrix.C:1242

#4  ? at /usr/local/openfoam/v1706/OpenFOAM-v1706/src/finiteVolume/lnInclude/fvMatrix.C:2247

#5  ? at ~/OpenFOAM/rbardia-v1706/applications/solvers/temperatureFoam/pEqn.H:50

#6  __libc_start_main in "/lib64/libc.so.6"

#7  Foam::UOPstream::write(float) in "/home/rbardia/OpenFOAM/rbardia-v1706/platforms/linux64GccDPInt32Opt/bin/temperatureFoam"

Aborted

The error is pointing to the factor "isLiquid", which I am multiplying with the "fvm::laplacian" operation in the formation of the pEqn matrix. This multiplication apparently is not allowed for matrices that use face flux correction often implemented for non-orthogonal matrices by using the corrected laplacian scheme as shown here. The fvSchemes file I am using is also attached.

Code:

laplacianSchemes

{

    default         Gauss linear corrected;

}

However, my question stems from the fact that I am writing a similar scalar matrix for another variable "T", which is not posing this issue. TEqn.H is attached and the matrix is shown below:

Code:

fvScalarMatrix TEqn

(

    fvm::ddt(rhoCp ,T) 

    - isLiquid * fvm::laplacian(kLEff, T)

    - (1 - isLiquid) * fvm::laplacian(kVEff, T)

    ==  

    isLiquid * fvc::laplacian(kLEff * DTField, isInterfaceOwnerLiquid)

    + (1 - isLiquid) * fvc::laplacian(kVEff * DTField, isInterfaceOwnerLiquid)

    - isLiquid * fvc::laplacian(kLEff * DTField, isInterfaceOwnerVapor)

    - (1 - isLiquid) * fvc::laplacian(kVEff * DTField, isInterfaceOwnerVapor)

);

Can someone help me with the reason why the faceFluxCorrection is not thrown for the TEqn but it happens for pEqn? Is there a way this multiplication can be done for the pEqn matrix as well.

Help is much needed and greatly appreciated.

Thanks.