{
    volScalarField divPhi
    (
        "divPhi",
        fvc::div(phi)
    );

    // Momentum equation
    fvVectorMatrix UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

    rAU = 1.0/UEqn.A();

    // Insert the additional components. Note this will destroy the H and A

    UEqn += fvm::SuSp(-divPhi, U) + divPhi*U;
    UEqn.relax();

    UpEqn.insertEquation(0, UEqn);
}