## Mathematical Model  
### Turbulence Model \(k\)-\(\Omega\)-\(k_{\theta}\)-\(\Omega_{\theta}\)  
The Reynolds averaged Navier-Stokes equations for incompressible flows can be written as [8, 22]  
\[\mathbf{\
abla}\cdot\mathbf{u}=0\,, \tag{1}\] \[\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\mathbf{\
abla} \mathbf{u}=-\frac{1}{\rho}\mathbf{\
abla}P+\mathbf{\
abla}\cdot\left[(\
u+\
u_{t}) \left(\mathbf{\
abla}\mathbf{u}+\mathbf{\
abla}\mathbf{u}^{T}\right)\right]+\mathbf{g}\,,\] (2) \[\frac{\partial T}{\partial t}+\mathbf{u}\cdot\mathbf{\
abla}T=\mathbf{ \
abla}\cdot\left[(\alpha+\alpha_{t})\,\mathbf{\
abla}T\right]+\frac{Q}{\rho C_{p }}\,. \tag{3}\]  
where \(\mathbf{u}\) is the average fluid velocity, \(P\) the pressure, \(T\) is the average temperature and \(Q\) represents the volume heat source or sink. In Equations (2) and (3) Reynolds stresses \(\mathbf{\tau}^{R}\) and the turbulent heat flux \(\mathbf{q}^{R}\) are calculated using the mean velocity \(\mathbf{u}\), the mean temperature \(T\) and two new variables, the eddy viscosity \(\
u_{t}\) and the eddy heat diffusivity \(\alpha_{t}\) by setting [25, 26]  
\[\mathbf{\tau}^{R}=\rho\overline{\mathbf{u}^{\prime}\mathbf{u}^{\prime }}=-\rho\
u_{t}\left(\mathbf{\
abla}\mathbf{u}+\mathbf{\
abla}\mathbf{u}^{T}\right)+ \rho\frac{2k}{3}\,\mathbf{I}\,, \tag{4}\] \[\mathbf{q}^{R}=\rho C_{p}\overline{\mathbf{u}^{\prime}T^{\prime}} =-\rho C_{p}\alpha_{t}\,\mathbf{\
abla}T\,, \tag{5}\]  
where \(\mathbf{u}^{\prime}\) and \(T^{\prime}\) are the values of fluctuating velocity and temperature field, respectively. The hypothesis in Equation (4) is known as the Boussinesq hypothesis, wherein \(k\) represents the turbulent kinetic energy. The hypothesis in Equation (5) is known as diffusion gradient hypothesis. When \(\alpha_{t}\) is a constant we have simple eddy diffusivity model (SED), but it is possible to consider \(\alpha_{t}\) as a function of the velocity field, leading to the generalized diffusion gradient model. In this paper, in order to extend the range of the turbulence model validity with respect to the Prandtl number, we employ a four-parameter logarithmic turbulence model \(k\)-\(\Omega\)-\(k_{\theta}\)-\(\Omega_{\theta}\) (KLW model). It is obtained through a change of state variables in the \(k\)-\(\varepsilon\)-\(k_{\theta}\)-\(\varepsilon_{\theta}\) turbulence model [8, 16, 17, 18]. Let \(k\) be the turbulent kinetic energy, \(\varepsilon\) be its dissipation rate, \(k_{\theta}\) be the mean-squared temperature fluctuations value and \(\varepsilon_{\theta}\) be its dissipation rate defined as  
\[k:=\frac{1}{2}\overline{\mathbf{u}^{\prime 2}}\,,\qquad\varepsilon:=v\| \overline{\mathbf{\
abla}\mathbf{u}^{\prime}}\|^{2}\,,\qquad k_{\theta}:=\frac{1 }{2}\overline{T^{\prime 2}}\,,\qquad\varepsilon_{\theta}:=\alpha\|\overline{\mathbf{\
abla}T^{ \prime}}\|^{2}\,. \tag{6}\]  
To compute the time scales of turbulence we propose to employ logarithmic values of \(k\) and \(k_{\theta}\) specific dissipation rate \(\omega\) and \(\omega_{\theta}\),  
\[\Omega=\ln\left(\frac{\varepsilon}{C_{\mu}k}\right)\,,\qquad\Omega_{\theta}= \ln\left(\frac{\varepsilon_{\theta}}{C_{\mu}k_{\theta}}\right)\,, \tag{7}\]  
where \(C_{\mu}\) is a model constant. The motivation to use logarithmic variables is well known [8, 23]. During numerical computations the state values of turbulence equations, such that \(k_{\theta}\), \(\varepsilon_{\theta}\), \(k\) and \(\varepsilon\), may take negative values near the boundary. Laplacian operator with negative coefficients looses its stability and convergence of the Navier-Stokes system may not be attained. This is one of the main problems with \(k\)-\(\varepsilon\) models and numerical limiters are very common [23]. Furthermore logarithmic variables take smoother variations at the wall where turbulence state variable profiles take larger changes [23].