Generally, incompressible viscous turbulent heat transfer energy equation is
[latex]\dfrac{\partial T}{\partial t} + \nabla\cdot(\boldsymbol{U} T) =
\alpha \nabla^2 T + \dfrac{\nu}{c_\mathrm{p}} [\nabla \boldsymbol{U} +(\nabla \boldsymbol{U})^\mathrm{T}]:\nabla \boldsymbol{U}
[/latex]
ddt(T) + div(U T) = alpha laplacian(T) + ((nu/cp (grad(U) + (grad(U)).T)) && grad(U))
where [latex]\alpha = \dfrac{k}{\rho c_\mathrm{p}}[/latex]
When turbulent stress considered, the equation becomes
[latex]\dfrac{\partial T}{\partial t} + \nabla\cdot(\boldsymbol{U} T) =
\alpha_\mathrm{eff}\nabla^2 T + \dfrac{\nu_\mathrm{eff}}{c_\mathrm{p}} [\nabla \boldsymbol{U} +(\nabla \boldsymbol{U})^\mathrm{T}]:\nabla \boldsymbol{U}
[/latex]
where [latex]\alpha_\mathrm{eff}= \nu/Pr + \nu_t/Pr_t[/latex], $latex Pr$ is Prandl number, $latex \nu_t$ is turbulent viscosity and $latex Pr_t$ is turbulent Prandl number.
With OpenFOAM, the incompressible viscous turbulence energy equation is written
volScalarField alphaEff
(
"alphaEff",
turbulence->nu()/Pr + turbulence->nut()/Prt
);
volTensorField gradU = fvc::grad(U);
volTensorField tauEff = turbulence->nuEff() / cp * (gradU + gradU.T());
{
// Add energy equation
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(alphaEff, T)
- (tauEff && gradU)
);
TEqn.solve();
}