// Solve the Momentum equation

fvVectorMatrix UEqn
(
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  + fvc::div
    (
        (Alpha/(scalar(1.001) - Alpha))*(sqr(rhoc)/rho)*Vdj*Vdj,
        "div(phiVdj,Vdj)"
    )
  - fvm::laplacian(mu, U, "laplacian(muEff,U)")
);

//solve(UEqn == -fvc::grad(p) - fvc::grad(rho)*gh);
