[WaterHammer1985]

Hello,

I'm trying to validate a weir coefficient using interIsoFoam but have run into an issue (see attached image). If I use the flow depth interface (alpha=0.5), the weir coefficient is much lower than anticipated. If I use the pressure to calculate the head, the weir coefficient is in the expected range. So I am unsure why the pressure matches what would be anticipated but the water surface depth does not.

Questions:
1. What would cause the pressure at the interface to be negative? There are ~7 mesh cells between 0Pa and the alpha interface so it isn't negative in just 1 cell (see the line plot in attached image)

2. Is there a way to have the interface correspond to 0Pa? If I change the alpha phase to 0.99 the interface is still negative.

0/U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{

    inlet_water
    {
        type                      flowRateInletVelocity;
        volumetricFlowRate        constant 30;
        value                     uniform (0 0 0);
    }

    outlet_air
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
    outlet_water
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

    "(inlet_wall|conc)"
    {
        type            noSlip;
    }

    "(wall1|wall2)"
    {
        type            symmetry;
    }



    "(bed1|bed2)"
    {
        type            noSlip;
    }


    top
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

    defaultFaces
    {
        type            empty;

0/alpha.water

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet_water
    {
        type            fixedValue;
        value           uniform 1;

    outlet_air
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }

    outlet_water
    {
        type            inletOutlet;
        inletValue      uniform 1;
        value           uniform 1;
    }


    "(inlet_wall|conc)"
    {
        type            zeroGradient;
    }

    "(wall1|wall2)"
    {
        type            symmetry;
    }


    "(bed1|bed2)"
    {
        type            zeroGradient;
    }


    top
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }

    defaultFaces
    {
        type            empty;
    }
}

0/p_rgh

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform $pressure;

boundaryField
{
    inlet_water
    {
        type            fixedFluxPressure;
    }

    outlet_water
    {
        type            totalPressure;
        rho             rho;
        psi             none;
        gamma           1;
        p0              uniform 4286970;
        value           uniform 4286970;
    }

    outlet_air
    {
        type            totalPressure;
        rho             rho;
        psi             none;
        gamma           1;
        p0              uniform 0;
        value           uniform 0;
    }


    "(inlet_wall|conc)"
    {
        type            fixedFluxPressure;
    }

    "(wall1|wall2)"
    {
        type            symmetry;
    }


    "(bed1|bed2)"
    {
            type            fixedFluxPressure;
    }


     top
    {
        type            totalPressure;
        p0              uniform 0;
    }


    defaultFaces
    {
        type            empty;
    }
}

fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(rhoPhi,U)  Gauss upwind; //Gauss limitedLinearV 1;
    div(phi,alpha)  Gauss upwind; //Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
    div(phi,omega) Gauss upwind;
    div(phi,p_rgh)  Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

fvSolution

solvers
{
    "alpha.water.*"
    {
        isoFaceTol      1e-6;
        surfCellTol     1e-6;
        nAlphaBounds    3;
        snapTol         1e-12;
        clip            true;

        nAlphaSubCycles 1;
        cAlpha          1; // Note: cAlpha is not used by isoAdvector but must
                           // be specified because interfacePropertes object
                           // reads it during construction.
    }

    "pcorr.*"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-10;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        smoother        DICGaussSeidel;
        tolerance       1e-07;
        relTol          0.05;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-07;
        relTol          0;
    }



    "(U|k|omega|epsilon).*"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }



    "(U|k|epsilon)Final"
    {
        $U;
        relTol          0;
    }


}

PIMPLE
{
    momentumPredictor no; 
    nCorrectors     3;
    nOuterCorrectors 1;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


relaxationFactors
{
    fields
    {
    }
    equations
    {
        ".*" 0.1;
    }
}

[WaterHammer1985]

Additional information:

I'm using snappyHexMesh for the domain and mesh quality is good. The problem persists using interFoam and also doesn't appear to change based on fvSchemes order (1st vs 2nd). And none of the tutorial cases (weirOverflow, waterChannel, damBreak) have the same issue.

Any suggestions on where else to look to track this down would be greatly appreciated!