    // Solve the Energy equation

    fvScalarMatrix TEqn
    (
            fvm::ddt(T)
          + fvm::div(mdotf, T)
          - fvm::laplacian(K, T)
          + fvm::SuSp(-fvc::div(mdotf),T)
   );

     TEqn.relax();
     TEqn.solve();



