// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

volTensorField gradU(fvc::grad(U));

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + MRF.DDt(rho, U)
  + (
     - fvc::div(diff.mu()*Foam::dev2(Foam::T(gradU)))
     - fvm::laplacian(diff.mu(), U)
    )
 ==
    rho * thermo.composition().Qc() * E
  + fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
    K = 0.5*magSqr(U);
}