4.1. BL–plume interface

Before proceeding to a detailed analysis of the integral model (3.6), a rigorous definition of the BL–plume region interface, $\delta _i$, is required. Since the flow is driven by buoyancy, we use the heat transfer characteristics to identify the different mechanics of the BL-like and plume-like region of the boundary layer. For the vertical NCBL flow, heat is supplied by the isothermal wall and is conducted away from the wall through molecular diffusion $-\rho c_p \kappa \partial \bar \theta /\partial y$ and turbulent heat transport $\rho c_p \overline {v'\theta '}$. In our previous work (Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020), a constant heat flux layer, in which region the total transport of heat $\mathcal {Q}$ is approximately constant, such that

(4.1)\begin{equation} \mathcal{Q}\equiv \rho c_p \left(\overline{v'\theta'}-\kappa \frac{\partial \bar \theta}{\partial y}\right) =\mathcal{Q}_w, \end{equation}

has been identified for the temporally developing NCBL up to $\textit {Gr}_\delta =7.7\times 10^7$. The constant heat flux layer is consistent with the previous studies for turbulent spatially developing NCBLs (George & Capp Reference George and Capp1979; Hölling & Herwig Reference Hölling and Herwig2005; Wells & Worster Reference Wells and Worster2008). Here, the heat flux on the wall is given by

(4.2)\begin{equation} \mathcal{Q}_w ={-}\rho c_p \kappa \left.\frac{\partial \bar \theta}{\partial y}\right\rvert_{y=0}, \end{equation}

and $c_p$ is the specific heat of the fluid. In the present study, such a constant heat flux layer is found for $y/L_s<40$ (or $y^+<180$ in wall units) at $\textit {Gr}_\delta =1.21\times 10^8$ for DNS-C, as demonstrated by figure 6(a). Figure 6(b) also depicts a constant forcing region in the vicinity of the isothermal wall as introduced by Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020), so that the wall shear stress

(4.3)\begin{equation} {\tau}_w =\rho \nu \left.\frac{\partial \bar u}{\partial y}\right\rvert_{y=0}, \end{equation}

balances the local shear stress and the buoyancy force,

(4.4)\begin{equation} {\tau}_w = \rho\left( \nu\frac{\partial \bar u}{\partial y}-\overline{u'v'}\right)+\int_0^y \bar b \,\mathrm{d} y, \end{equation}

for $y/L_s<40$ ($y^+<180$) at $\textit {Gr}_\delta = 1.21\times 10^8$. The existence of the constant flux and forcing layers given by (4.1) and (4.4) indicates that there exists a slowly varying region in the close vicinity of the wall where the acceleration of the flow has negligible influence on the dynamics.

Figure 6.

Components of the (a) near-wall heat flux and (b) near-wall forcing in the near-wall region, showing the constant heat flux region and the constant momentum forcing region. Data obtained using DNS-C at $\textit {Gr}_\delta =1.21\times 10^8$.

Figure 7 depicts the mean streamwise velocity and temperature profiles for the temporally developing NCBL for both DNS-A and DNS-C. The profiles are normalised in wall units, where

(4.5a,b)\begin{equation} u^+\,{=}\, \frac{\bar{u}}{u_\tau}, \quad \theta^+\,{=}\, \frac{\theta_w-\bar{\theta}}{\theta_\tau}, \end{equation}

and $u_\tau$, $\theta _\tau$ are the friction velocity and temperature, respectively, defined by

(4.6a,b)\begin{equation} u_\tau = \sqrt{\frac{\tau_w}{\rho}},\quad \theta_\tau = \frac{\mathcal{Q}_w}{\rho c_p u_\tau}. \end{equation}

The mean streamwise velocity, as shown in figure 7(a), demonstrates a strong $\textit {Gr}_\delta$ dependency due to the additional buoyancy effect and follows a modified log-law as suggested by Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020), whereas the mean temperature, as shown in figure 7(b), appears independent of the $\textit {Gr}_\delta$ and features a temperature log-law (Tsuji & Nagano Reference Tsuji and Nagano1988; Ng et al. Reference Ng, Ooi, Lohse and Chung2017; Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020). In the present study, we identify the BL–plume interface $\delta _i$ to be the maximum mean streamwise velocity location for the results presented, i.e. the $y$ location at which $\bar {u}= \bar {u}_m$; this interface is found to be very close to the edge of the constant flux (and forcing) region. For instance, at $\textit {Gr}_\delta =1.2\times 10^8$ the BL–plume interface ($\bar {u}_m$ location) is located at $\delta _i^+ = 171.21$ in wall units; whereas the constant flux (and forcing) region extends up to $y^+\approx 180$. Nevertheless, in Appendix B, we revert the definition of $\delta _i$ to be the wall-normal location at which the molecular diffusion matches the turbulent heat flux (e.g. Wells & Worster Reference Wells and Worster2008), and the integral model is shown to be robust and insensitive to the exact definition of $\delta _i$, provided $\delta _i$ falls in the constant flux region.

Figure 7.

The mean profiles of the temporally developing NCBL at different $\textit {Gr}_\delta$ in wall units. (a) Mean streamwise velocity profile and (b) mean temperature profiles. The black dotted line in panel (b) indicates the logarithmic profile.

4.2. A self-similar model

Based on the definition of the BL–plume interface $\delta _i$, the buoyancy BC for the integral model is then obtained by

(4.7)\begin{equation} \left.{g\alpha } {q}_{\delta_i}\equiv \left(-\kappa\frac{\partial \bar b}{\partial y}-\overline{v'b'}\right)\right\rvert_{y=\delta_i}\approx{g\alpha } {q}_w, \end{equation}

where ${q} = \mathcal {Q}/{\rho c_p}$. The rate of change in the total buoyancy (3.6b) is therefore directly linked to the wall heat flux ${q}_w$,

(4.8)\begin{equation} \frac{\partial B}{\partial t} = \hat{u}\hat{b}(\beta_{\delta_i}-\mathcal{L}_B) ={g\alpha} {q}_{\delta_i} - \hat{u}\hat{b}\mathcal{L}_B={g\alpha } {q}_w- \hat{u}\hat{b}\mathcal{L}_B. \end{equation}

In the present study, the BCs of the integral model (3.7) are obtained by a linear interpolation of the DNS data onto the interface $y=\delta _i$. Figure 8 demonstrates that $\beta _{\delta _i}$, being the dimensionless buoyancy BC for the integral model, is approximately constant for $\textit {Gr}_\delta >4\times 10^6$. This suggests that the heat flux through the interface $q_{\delta _i}$ scales with the local characteristic scales,

(4.9)\begin{equation} q_{\delta_i} \propto \hat{u}\hat{b}. \end{equation}

In addition, (4.7) and (4.9) further indicate that the heat flux on the wall scales with the plume-like region scales,

(4.10)\begin{equation} q_w\propto\hat{u}\hat{b}. \end{equation}

The implication of (4.10) is that the heat transfer of the turbulent NCBL on the wall is driven by the characteristic velocity and temperature scales in the plume-like region. As depicted by figure 8, the normalised wall heat flux shows excellent agreement with the buoyancy BC $\beta _{\delta _i}$ and also appears constant for $\textit {Gr}_\delta >4\times 10^6$. Similarly, the shear and integral production BCs for the integral model, scaled by the top-hat characteristic quantities, are shown in figure 9. It is clear that both BCs scale well with the characteristic scales since they are constants in the range $\textit {Gr}_\delta >10^7$. In addition, the dimensionless Leibniz terms are also appear constant in this Grashof number range. Figure 10 shows the temporal evolution of the Leibniz term $\mathcal {L}_Q$ for the volume flux $Q$. Evidently, in the fully turbulent regime $\mathcal {L}_Q$ appears constant with some level of temporal oscillation. The Leibniz terms for buoyancy $\mathcal {L}_B$ and momentum $\mathcal {L}_M$ (although not shown here) also demonstrate similar temporal fluctuations. Nevertheless, taking a temporal average for $\textit {Gr}_\delta >10^7$ yields $\mathcal {L}_Q\approx 0.0029$, $\mathcal {L}_B\approx 0.0041$ and $\mathcal {L}_M\approx 0.0042$. The constant Leibniz term $\mathcal {L}_Q$ further implies that the spreading rate of the interface $\partial \delta _i/\partial t$ scales with the velocity scale $\hat {u}^2/\bar {u}_m$. In the present study, the maximum streamwise velocity $\bar {u}_m$ is found to scale with the characteristic velocity in the plume-like region (see, later in § 4.6), suggesting

(4.11)\begin{equation} \frac{\partial \delta_i}{\partial t} = E_1\hat{u} = E_2 \bar{u}_m, \end{equation}

where $E_1$ and $E_2$ are the proportionality constants. Equation (4.11) takes a similar form to the entrainment relation (Morton et al. Reference Morton, Taylor and Turner1956), and may interpreted as the entrainment from the turbulent plume-like region to the near-wall BL-like region.

Figure 8.

Temporal development of the dimensionless heat (buoyancy) BC $\beta _{\delta _i}=g\alpha q_{\delta _i}/\hat {u}\hat {b}$ and the normalised wall heat flux $g\alpha q_w/\hat {u}\hat {b}$.

Figure 9.

Normalised BCs for the integral model: (a) shear BC, $S_{\delta _i} = \tau _{\delta _i}/\hat {u}^2$; (b) integral production BC, $M_{\delta _i} = \varPi _{\delta _i}/\hat {u}^3$. The black dash–dotted lines indicates the values listed in table 2.

Figure 10.

The temporal development of the normalised Leibniz term $\mathcal {L}_Q$ in (3.8). The black dash–dotted line indicates $\mathcal {L}_Q=0.0029$.

Figure 11.

Profile coefficients for the integral model: (a) contribution of total shear production to the MKE, $\gamma _p$; (b) contribution of viscous diffusion to the MKE, $\gamma _f$; (c) contribution of viscous dissipation to the MKE, $\gamma _d$; (d) contribution of buoyancy to the MKE, $\gamma _b$. The black dash–dotted lines indicates the values listed in table 2.

Profile coefficients of the integral model (3.6), shown in figure 11, are numerically obtained by integrating the DNS data from $\delta _i$ to $L_y$. In figure 11, the profile coefficients $\gamma _p$ and $\gamma _b$ are clearly independent of $\textit {Gr}_\delta$ for $\textit {Gr}_\delta >10^7$, suggesting the plume-like region is self-similar; while the viscous terms $\gamma _d$ and $\gamma _f$ are of order $O(10^{-4})$ and $O(10^{-5})$, respectively, and gradually reducing in magnitude with increasing $\textit {Gr}_\delta$. The physical interpretation of figure 11 is straightforward: the bulk vigorous turbulent NCBL is accelerated and maintained by the buoyancy $\gamma _b$ while a small portion of the MKE is transferred to turbulent kinetic energy through shear production $\gamma _p$, and diffused $\gamma _f$ and dissipated $\gamma _d$ through viscosity. With increasing $\textit {Gr}_\delta$ the flow become more turbulent and the viscous effects diminish; while the driving buoyancy and shear remains significant at high $\textit {Gr}_\delta$ turbulent flow. It is noted that the buoyancy contribution $\gamma _b$ needs to be multiplied by $\textit {Ri}$ before it can be compared with the other three terms. In the present study, the Richardson number becomes a constant $\textit {Ri}\approx 0.048$ in the turbulent regime $\textit {Gr}_\delta >10^7$, as shown in figure 12, further supporting that the plume-like region of the turbulent NCBL is self-similar for $\textit {Gr}_\delta >10^7$. This indicates that the total buoyancy contribution to the MKE ($\gamma _b\textit {Ri}\approx 0.098$) is of leading order, while the viscous terms are negligible in the MKE balance as expected since the flow becomes more turbulent at high $\textit {Gr}_\delta$, and the plume-like region is away from the wall.

The close agreement of the DNS-A and DNS-C profile and BC coefficients shown in figures 8–12 in the range $5\times 10^7<\textit {Gr}_\delta <7\times 10^7$ shows our initialisation method for DNS-C is not disturbing the flow in any significant way, and the profile coefficients are domain size independent for the domains considered here. Results for the constant profile and BC coefficients obtained from DNS are summarised in table 2. Note the values listed in table 2 depend on the choice of $\delta _i$, however, we will show later in Appendix B that the self-similar integral model proposed here is robust and insensitive to the choice of the BL–plume interface $\delta _i$.

Figure 12.

The development of Richardson number $\textit {Ri}$ for the plume-like region. The black dash–dotted line indicates $\textit {Ri}=0.048$.

Table 2.

Average profile constants obtained from DNS for $10^7<\textit {Gr}_\delta <1.2\times 10^8$ based on the maximum streamwise velocity location.

4.3. Entrainment

The canonical entrainment correlation in the sense of Morton et al. (Reference Morton, Taylor and Turner1956) needs reformulation for our NCBL here since the continuity equation (2.1a) for a temporally evolving parallel flow vanishes. In the present study, we define the entrainment velocity for plume-like region as the growth (spreading) rate of the boundary layer with respect to the interface $y=\delta _i$,

(4.12)\begin{equation} v_e = \frac{\partial \hat{l}} {\partial t} = E \hat{u}, \end{equation}

where $v_e$ is the entrainment (expansion or spreading) velocity of the boundary layer and $E$ is the entrainment (spreading) coefficient as in Morton et al. (Reference Morton, Taylor and Turner1956) and van Reeuwijk, Vassilicos & Craske (Reference van Reeuwijk, Vassilicos and Craske2021). This definition, (4.12), is consistent with the unified description of the turbulent entrainment (van Reeuwijk et al. Reference van Reeuwijk, Vassilicos and Craske2021) and has been adopted by a number of temporally developing flows, for example Jonker et al. (Reference Jonker, van Reeuwijk, Sullivan and Patton2013) for a shear-driven flow with buoyancy stratification and Krug et al. (Reference Krug, Chung, Philip and Marusic2017) for a planar plume. With the definition of the top-hat characteristic velocity and length scales (3.4), the entrainment coefficient is obtained in terms of the integral quantities,

(4.13)\begin{equation} E = \frac{Q}{M}\frac{\partial}{\partial t}\left(\frac{Q^2}{M}\right)=\frac{Q}{M} \left(\frac{2Q}{M}\frac{\partial Q}{\partial t} -\frac{Q^2}{M^2}\frac{\partial M}{\partial t}\right). \end{equation}

By substituting (3.6a) and (3.6c) into (4.13), the entrainment coefficient now reads

(4.14)\begin{equation} E= \frac{2Q^2}{M^2}-2(S_{\delta_i}+\mathcal{L}_Q) - (M_{\delta_i}+\varSigma_M - \mathcal{L}_M) = E_b + E_s + E_p + E_{\nu}, \end{equation}

where $E_b$ is the net buoyancy contribution to the entrainment coefficient, $E_s$ the boundary effects on the entrainment coefficient, $E_p$ the shear production effect on the entrainment coefficient and $E_\nu$ the viscous contribution to the entrainment coefficient. These terms are expressed as

(4.15a)\begin{gather} E_{b} = Ri(2-\gamma_b), \end{gather}

(4.15b)\begin{gather}E_{s} = \mathcal{L}_M-(2S_{\delta_i}+2\mathcal{L}_Q+M_{\delta_i}), \end{gather}

(4.15c)\begin{gather}E_{p} ={-}\gamma_p, \end{gather}

(4.15d)\begin{gather}E_{\nu} ={-}(\gamma_f+\gamma_d). \end{gather}

The entrainment coefficient for the temporally developing NCBL obtained by (4.14) makes no assumption beyond Reynolds averaging and is in general consistent with the entrainment relation for spatially developing jets and plumes (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015; van Reeuwijk & Craske Reference van Reeuwijk and Craske2015) but with the additional BC $E_s$ and viscous effects $E_\nu$.

The total entrainment coefficient $E$ is then obtained using (4.14) with the profile coefficients given in table 2. Results are shown in figure 13(a), where $E$ calculated from (4.14) is compared with that from (4.12), i.e. $E=v_e/ \hat {u}$. Evidently, from figure 13(a), the calculated $E$ shows good agreement between the two methods. Figure 13(b) demonstrates the development of the individual contributions, (4.15), to the entrainment coefficient against $\textit {Gr}_\delta$. It is evident that all the components tend to be constants beyond $\textit {Gr}_\delta =10^7$. The shear production term $E_p$, among all the components, is of the leading order. In the present study, we have found $E_p\approx 0.041$, which is approximately half of the values reported by previous studies for free shear flows: in an energy-consistent analysis, van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015) reported that the shear production contribution is insensitive to the Richardson number, where $E_p \approx 0.08$, for the spatially developing axisymmetric jets ($\textit {Ri}=0$) and plumes ($0.1< Ri<0.2$). Krug et al. (Reference Krug, Chung, Philip and Marusic2017) numerically obtained $E_p\approx 0.08$ for a temporally developing planar plume at $\textit {Ri}=0.2$ and $E=E_p\approx 0.12$ for a temporally developing planar jet ($\textit {Ri}=0$) using DNS. This discrepancy is essentially due to the shear production profile difference between the free shear flows and the NCBL, as shown in figure 14. It is worth noting that the data presented in Krug et al. (Reference Krug, Chung, Philip and Marusic2017) was originally normalised by a centreline velocity of the plume $\bar {w}_{cl}$, where such centreline velocity is not applicable in the current NCBL flow. As such, we have assumed the mean profile of their plume data to be Gaussian, so that $\bar {w}_{cl}=\sqrt {2}\hat {u}$ is used to compare with our data in figure 14. As shown in figure 14(a), the NCBL develops two shear layers. The first peak of the shear production occurs in the near-wall region since the NCBL is bounded by the wall (see also in figure 6a, the maximum shear production occurs when $-\overline {u'v'} = \nu \partial \bar u /\partial y$). The near-wall peak, when scaled by the characteristic scales, increases with increase $\textit {Gr}_\delta$, whereas the second peak takes place in the plume-like region ($y>\delta _i$). In the plume-like region, our DNS data shows fairly good agreement between the relatively low $\textit {Gr}_\delta$ case (DNS-A) and high $\textit {Gr}_\delta$ case (DNS-C), convincingly showing that the statistics of the plume-like region are in the self-similar state. We stress that the near-wall shear production contributes minimally to the integral shear production term $E_p=-\gamma _p$ given the near-wall region is extremely thin, for example, the difference between integrating from $y=0$ and $y=\delta _i$ to $y=L_y$ at $\textit {Gr}_\delta =1.2\times 10^8$ is less than $5\,\%$. This is more clearly demonstrated in figure 14(b), where the shear production is premultiplied by the dimensionless wall-normal distance $y/\hat {l}$. Although different in shape, the contributions to the integral (area under the curve) for the present study roughly matches the temporal plume data reported by Krug et al. (Reference Krug, Chung, Philip and Marusic2017) in figure 14(b). Note in Krug et al. (Reference Krug, Chung, Philip and Marusic2017) the shear production contribution to the entrainment is defined by the integral in the cross-stream direction, i.e. $-\infty$ to $\infty$, whereas figure 14 only shows $0$ to $\infty$. As a result, the total shear production in the plume-like region is approximately half of the value reported by Krug et al. (Reference Krug, Chung, Philip and Marusic2017). More recently, Parker et al. (Reference Parker, Burridge, Partridge and Linden2020) compared the entrainment due to shear production between a spatially developing wall plume ($\textit {Ri}=0.09$) and a free planar plume ($\textit {Ri}=0.23$) and found the total shear production for the wall plume ($E_p=0.067$) is slightly higher than half of that for the free plume ($E_p=0.103$). The discrepancy between the NCBL flow reported here and the wall plume may be due to differences in the driving mechanisms and flow conditions. Nevertheless, both planar (temporal) and axisymmetric (spatial) flows suggest the entrainment process is dominated by the shear production (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015; van Reeuwijk & Craske Reference van Reeuwijk and Craske2015; Krug et al. Reference Krug, Chung, Philip and Marusic2017), consistent with the current study.

Figure 13.

Development of the entrainment coefficient. (a) The total entrainment coefficient $E$ and (b) buoyancy $E_b$, shear $E_p$, BCs $E_s$ and viscous $E_\nu$ contributions to the entrainment coefficient $E$.

Figure 14.

Dimensionless (a) shear production and (b) premultiplied shear production for the NCBL. Grey solid lines show the instantaneous shear production profiles for the turbulent NCBL from both DNS-A and DNS-C, in the range $10^7<\textit {Gr}_\delta <1.2\times 10^8$; solid lines also show the ensemble average of the dimensionless shear production in DNS-A (blue, over the range $10^7<\textit {Gr}_\delta <7.7\times 10^7$) and DNS-C (red, over the range $5\times 10^7<\textit {Gr}_\delta <1.2\times 10^8$); black dashed line indicates the shear production profile take from Krug et al. (Reference Krug, Chung, Philip and Marusic2017) for a temporally evolving planar plume; vertical line indicates the starting point of the plume-like region (maximum velocity location) $y=\delta _i$.

The effect of buoyancy on the entrainment, given in (4.15a), can be rewritten in the form of

(4.16)\begin{equation} E_b = Ri(2-\gamma_b) = \frac{2}{\hat{u}^3}\left(\frac{\displaystyle\int^{\infty}_{\delta_i}{\bar{u}}^2\, \mathrm{d} y \int^{\infty}_{\delta_i}\bar b \,\mathrm{d} y}{\displaystyle\int^{\infty}_{\delta_i}\bar u \,\mathrm{d} y} - \int^{\infty}_{\delta_i} \bar u \bar b \,\mathrm{d} y\right), \end{equation}

where the first term in the bracket is associated with the rate of change in volume flux $\partial Q/\partial t$ due to the total buoyancy and the second term in the bracket represents the rate of change in the momentum flux $\partial M / \partial t$ due to the mean buoyancy flux in the plume-like region. The whole term (4.16) can only become zero when the mean velocity and mean buoyancy (temperature) profiles share the same width and shape (see also Krug et al. Reference Krug, Chung, Philip and Marusic2017). This is consistent with the $\alpha$-MSG model derived by van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015), where the spreading rate ($\partial \hat {l}/\partial z$) is found independent of $E_b$, as the velocity and buoyancy are assumed to take Gaussian profiles of equal width in the model. Parker et al. (Reference Parker, Burridge, Partridge and Linden2020) obtained the buoyancy effect $E_b\approx 0.043\text {--}0.047$ for both spatially developing free and wall plumes. Their measurements are consistent with the experimental measurements reported by Paillat & Kaminski (Reference Paillat and Kaminski2014) $E_b\approx 0.04$. For a temporally developing planar plume, Krug et al. (Reference Krug, Chung, Philip and Marusic2017) suggests the buoyancy effect contributes approximately $15\,\%$ to the total entrainment coefficient, which is $E_b \approx 0.015$, smaller than that reported by Parker et al. (Reference Parker, Burridge, Partridge and Linden2020). In their DNS study, they point out that a smaller $E_b$ is attained due to the fundamental difference between the spatial and temporal flows – in spatially developing flows, the buoyancy additionally affects the streamwise transport of the momentum and energy, causing a net-buoyancy contribution to the entrainment, whereas such modification is absent in the temporal flows (Krug et al. Reference Krug, Chung, Philip and Marusic2017). In the present study, we have found $E_b\approx -0.0022$ negatively contributing to the spreading of the flow, and is one order of magnitude smaller than the values reported in the aforementioned studies. Such a result does not imply our results contradict previous studies, for the buoyancy contribution $E_b$ fully depends on the shapes of the profiles of velocity and buoyancy of the flow as given in (4.16). Figure 15 compares the velocity and buoyancy profiles for the plume and NCBL. According to figure 15(a), the streamwise mean velocity profiles for the plumes from different sources (Ramaprian & Chandrasekhara Reference Ramaprian and Chandrasekhara1989; Paillat & Kaminski Reference Paillat and Kaminski2014; Krug et al. Reference Krug, Chung, Philip and Marusic2017; Parker et al. Reference Parker, Burridge, Partridge and Linden2020) collapse into a single curve, while the wall plume DNS data (Parker et al. Reference Parker, Burridge, Partridge and Linden2020) demonstrates a slight difference in the near-wall region ($y^*/\hat {l}<0.1$) due to the wall friction. Here, the corrected coordinate is given by

(4.17)\begin{equation} y^*=\begin{cases} y-\delta_i, & \text{for NCBL},\\ y, & \text{for plumes}. \end{cases} \end{equation}

The DNS data for the NCBL, however, exhibits a different shape in comparison with the plume velocity profiles. Similar behaviour can also be observed in figure 15(b), where the mean buoyancy profile of the NCBL is shown along with those of the plumes. The wall plume case (Parker et al. Reference Parker, Burridge, Partridge and Linden2020) again diverges from the free shear plumes up to $y^*/\hat {l}=0.5$ due to the presence of the non-slip wall. It is clear that the plume buoyancy profiles given in figure 15(b) are quite different from the NCBL flow. A small buoyancy contribution (compared with $15\,\%$ in Krug et al. (Reference Krug, Chung, Philip and Marusic2017) and over $30\,\%$ in Parker et al. (Reference Parker, Burridge, Partridge and Linden2020)) is not surprising since the discrepancy in the profiles directly lead to a different buoyancy contribution to the entrainment. Assuming a self-similar mean flow, van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015) showed the buoyancy effect can both positively and negatively contribute to the entrainment coefficient, with the net-effect depending on the profile coefficients, for spatial plumes. In the present temporally developing NCBL, the profile coefficient $\gamma _b \approx 2.045$ is always greater than $2$. According to (4.15a), this determines a negative net-effect on the entrainment coefficient.

Figure 15.

Profile shape comparison of (a) mean streamwise velocity and (b) mean buoyancy. Filled markers represent experimental measurements for pure plumes, yellow circles, Paillat & Kaminski (Reference Paillat and Kaminski2014); blue triangles, Ramaprian & Chandrasekhara (Reference Ramaprian and Chandrasekhara1989). Lines represent DNS studies: red solid, planar plume reported by Krug et al. (Reference Krug, Chung, Philip and Marusic2017); blue solid and red dashed lines are the present study. Here $y^*$ is a corrected coordinate, given by $y^* = y-\delta _i$ for the NCBL and $y^*=y$ for the plumes.

The viscous and BC contributions to the entrainment coefficient are found to be $E_\nu \approx 2\times 10^{-4}$ and $E_s\approx -0.0037$, respectively, for the NCBL. With increasing $\textit {Gr}_\delta$, the viscous effect $E_\nu$ continues to lose its significance (see also in figure 11b,c); whereas the magnitude of $E_s$ appears independent of $\textit {Gr}_\delta$. By subtracting $E_b$ and $E_p$ from the spreading rate ($\partial \hat {l}/\partial z$), Parker et al. (Reference Parker, Burridge, Partridge and Linden2020) obtained a viscous component of $E_s^*\approx -0.03$ in their spatially developing wall plume. Note the viscous component $E_s^*$ in Parker et al. (Reference Parker, Burridge, Partridge and Linden2020) contains the combined effect of the BC $E_s$ and the viscous effect $E_\nu$ since it includes the wall shear and viscous dissipation. However, as noted by van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015), the entrainment and the spreading rate of the spatially developing plumes are usually different – except in the case of fully self-similar Gaussian plumes.

A brief comparison of the entrainment coefficient between the NCBL and previous studies is summarised in table 3. The current study suggests a constant entrainment coefficient of $E\approx 0.0353$ for the turbulent NCBL, giving good agreement with the experimental measurement ($E=0.04$) obtained by Bonnebaigt, Caulfield & Linden (Reference Bonnebaigt, Caulfield and Linden2018) for a plume with vertically distributed buoyancy source. It is evident that the production of turbulent kinetic energy $E_p$ is mainly responsible for the spreading of the plume-like region in the NCBL, whereas buoyancy $E_b$, viscous $E_\nu$ and BCs $E_s$ have negligible effects.

Table 3.

Entrainment (spreading) coefficient components for the temporally developing NCBL, plume and jets. Supercripts: $^{a}$ this is effectively $E_\nu +E_s$; $^{b}$ a plume with vertically distributed buoyancy source along the wall.

4.4. Solution to the self-similar system

As demonstrated in §§ 4.2 and 4.3, the outer plume-like region of the NCBL flow appears to be self-similar in the current $\textit {Gr}_\delta$ range and can be fully described by the integral model given in (3.6). It is convenient to rewrite (3.6) here in terms of the local (instantaneous) top-hat characteristic scales in order to obtain the solutions to the flow system:

(4.18a)\begin{gather} \frac{\partial }{\partial t} [\hat{u}\hat{l}] = \hat{b}\hat{l} - {\hat{u}^2}(S_{\delta_i}+\mathcal{L}_Q), \end{gather}

(4.18b)\begin{gather}\frac{\partial }{\partial t} [\hat{b}\hat{l}] = {\hat{u}\hat{b}}(\beta_{\delta_i}-\mathcal{L}_B), \end{gather}

(4.18c)\begin{gather}\frac{\partial }{\partial t} [\hat{u}^2\hat{l}]= \hat{u}^3(M_{\delta_i}+\gamma_p+\gamma_d+\gamma_f-\mathcal{L}_M)+ \hat{u}\hat{l}\hat{b}\gamma_b, \end{gather}

The dimensional consistency of (4.18) allows the following scaling relations for the characteristic scales in terms of total buoyancy $B$ and time $t$ (see Appendix A for detailed derivation):

(4.19a)\begin{align} \hat{u} = k_{\hat{u}}\sqrt{B}, \end{align}

(4.19b)\begin{align} \hat{b} = k_{\hat{b}}\sqrt{B}/t, \end{align}

(4.19c)\begin{align} \hat{l} = k_{\hat{l}}\sqrt{B}t, \end{align}

where $k_{\hat {u}}$, $k_{\hat {b}}$ and $k_{\hat {l}}$ are constants based on the profile coefficients,

(4.20a)\begin{gather} k_{\hat{u}} = \left[\frac{1-\gamma_b}{S_{\delta_i}+\mathcal{L}_Q+M_{\delta_i}+\gamma_p+\gamma_d+\gamma_f-(\beta_{\delta_i}-\mathcal{L}_B)/2}\right]^{1/2}, \end{gather}

(4.20b)\begin{gather}k_{\hat{b}} = \frac{1}{k_{\hat{l}}}, \end{gather}

(4.20c)\begin{gather}k_{\hat{l}} =\frac{1-k_{\hat{u}}^2(S_{\delta_i}+\mathcal{L}_Q+\beta_{\delta_i}-\mathcal{L}_B)} {k_{\hat{u}}}. \end{gather}

The scaling relations obtained here in (4.19) take a similar form to those suggested by Krug et al. (Reference Krug, Chung, Philip and Marusic2017) for the temporally developing planar plume. Note in Krug et al. (Reference Krug, Chung, Philip and Marusic2017), $B$ is a constant independent of time since the buoyancy is conserved; whereas in our temporal NCBL buoyancy is accumulating due to the existence of heat flux from the isothermal wall and fully depends on the local characteristic scales, i.e. $B(t)=\hat {b}\hat {l}$.

Since the profile coefficients are found to be constant for the current range of $\textit {Gr}_\delta$ (e.g. table 2), the set of ordinary differential equations (ODE) can be recast in the following form which encapsulates the entrainment assumption (4.12):

(4.21a)\begin{gather} \frac{\partial \hat{u}}{\partial t} = \hat{b}-(E+S_{\delta_i}+\mathcal{L}_Q)\frac{\hat{u}^2}{\hat{l}}, \end{gather}

(4.21b)\begin{gather}\frac{\partial \hat{b}}{\partial t} = (\beta_{\delta_i}-E-\mathcal{L}_B)\frac{\hat{u}\hat{b}}{\hat{l}}, \end{gather}

(4.21c)\begin{gather}\frac{\partial \hat{l}}{\partial t} = E \hat{u}. \end{gather}

System (4.21) yields a power-law solution (see the more detailed derivation in Appendix A),

(4.22a)\begin{gather} \hat{u} = A_{\hat{u}} t^{\xi} = k_{\hat{u}}A_{\hat{b}}^{1/2}A_{\hat{l}}^{1/2} t^{\xi}, \end{gather}

(4.22b)\begin{gather}\hat{b} = A_{\hat{b}} t^{\xi-1 }, \end{gather}

(4.22c)\begin{gather}\hat{l} = A_{\hat{l}} t^{ \xi+1 }, \end{gather}

where the exponent $\xi$ is given by

(4.23)\begin{equation} \xi = \frac{\beta_{\delta_i}-\mathcal{L}_B}{2E-\beta_{\delta_i}+\mathcal{L}_B}, \end{equation}

and $A_{\hat {b}}$ and $A_{\hat {l}}$ are constants from the integration which can be determined by the initial conditions. In the present study, the power-law solutions (4.22) are obtained using the constant profile coefficients in table 2 and the initial conditions $\hat {u}= 30.42$, $\hat {b}=0.1278$ and $\hat {l}=350.5$ from DNS-A at $\textit {Gr}_\delta =3.75\times 10^7$. The resultant integral constants and exponents for the power-law solution are tabulated in table 4. Note the initial condition is chosen arbitrarily from the DNS data, a different choice of the initial condition would result in a slightly different value for $A_b$ and $A_l$ due to the temporal nature of our DNS (for instance, at $\textit {Gr}_\delta =2.11\times 10^7$ gives $A_b=1.375$ and $A_l=0.0163$). Figure 16 compares the characteristic scales for the plume-like region predicted by the power-law solution (4.22) and the numerically generated data by DNS. The Grashof number is estimated by the characteristic length scale,

(4.24)\begin{equation} \textit{Gr}_\delta \equiv \frac{g\alpha \theta_w\delta^3}{\nu^2} = \frac{g\alpha \theta_w\hat{l}^3}{\nu^2} \frac{1}{\lambda_\delta^3}, \end{equation}

where $\lambda _\delta$ is the ratio of the characteristic length scale $\hat {l}$ to the velocity integral boundary layer thickness $\delta$. In the present study, it has been found that

(4.25)\begin{equation} \lambda_\delta \equiv \frac{\hat{l}}{\delta} \approx 1.33, \end{equation}

is a constant in the self-similar regime ($\textit {Gr}_\delta >10^7$). Shown by figure 16, the characteristic scales predicted by the ODE system exhibit excellent agreement with the DNS data, indicating the plume-like region can be well described by the self-similar model (4.21) using constant profile coefficients given in table 2.

Table 4.

Coefficients, integral constants and exponents for the scaling relation (4.20) and power-law solution (4.22) based on the profile coefficients given in table 2.

Figure 16.

Development of the top-hat characteristic scales for the plume-like region of the NCBL against $\textit {Gr}_\delta$: (a) the characteristic velocity scale $\hat {u}$; (b) the characteristic buoyancy scale $\hat {b}$; (c) the characteristic length scale $\hat {l}$.

4.5. Wall heat transfer characteristics

With the ODE system given in (4.21), it is also possible to predict the evolution of the wall heat flux with the growth of the boundary layer thickness for the turbulent NCBL. The wall heat transfer is typically characterised by the Nusselt number, which is defined as

(4.26)\begin{equation} \textit{Nu}_\delta \equiv \frac{\delta q_w}{\kappa \theta_w}. \end{equation}

Since the wall heat flux $q_w$ is directly linked to the BC for the integral model as given in (4.8), the Nusselt number in (4.26) can be written in terms of the local (instantaneous) top-hat characteristic scales,

(4.27)\begin{equation} \textit{Nu}_\delta = \frac{\beta_{\delta_i}}{ g\alpha \kappa \theta_w \lambda_\delta} \hat{u}\hat{b}\hat{l} \propto \hat{u}\hat{b}\hat{l}. \end{equation}

Figure 17(a) compares the $\textit {Nu}_\delta$ result of the present DNS study with the ODE prediction (4.27) and previous studies. For the spatially developing NCBL (Tsuji & Nagano Reference Tsuji and Nagano1988; Nakao et al. Reference Nakao, Hattori and Suto2017), the flow undergoes turbulent transition at $\textit {Gr}_\delta =10^6$, and the Nusselt number then has the empirical correlation,

(4.28)\begin{equation} \textit{Nu}_\delta = 0.107\textit{Gr}_\delta^{1/3}, \end{equation}

in the turbulent regime. Abedin et al. (Reference Abedin, Tsuji and Hattori2009) showed that their DNS data for a temporally developing NCBL closely follows the empirical correlation for the spatial flows. However, the temporally developing NCBL in the present study is not yet turbulent at $\textit {Gr}_\delta =10^6$. The discrepancy may be due to the fact that statistics in Abedin et al. (Reference Abedin, Tsuji and Hattori2009) are obtained by taking ensemble average of multiple simulations with different initialisation. At $\textit {Gr}_\delta =10^7$, our temporally developing NCBL becomes turbulent and the plume-like region appears self-similar. As shown in figure 17(a), our DNS data follow a similar trend to (4.28) in the turbulent regime. However, the power-law regression of our DNS data at $\textit {Gr}_\delta >10^7$ yields a slightly higher exponent than $1/3$, giving

(4.29a)\begin{gather} \textit{Nu}_\delta=0.0512\textit{Gr}_\delta^{0.369},\quad \text{for DNS-A}, \end{gather}

(4.29b)\begin{gather}\textit{Nu}_\delta=0.0399\textit{Gr}_\delta^{0.382},\quad \text{for DNS-C}, \end{gather}

with $95\,\%$ confidence bounds. The difference between the power-law scalings (4.28) and (4.29b) is rather subtle in the current DNS range, as shown in figure 17(b). In addition, the turbulence shows up in the Nusselt number due to the temporal nature of our NCBL, making the detection of an accurate power-law scaling difficult in the current $\textit {Gr}_\delta$ range. We, therefore, seek a different approach to accurately determine the scaling exponent for the Nusselt number in the current Grashof number range.

Figure 17.

Nusselt number development for the temporally developing NCBL: (a) comparing with the data reported by Tsuji & Nagano (Reference Tsuji and Nagano1988) and Nakao et al. (Reference Nakao, Hattori and Suto2017) for spatially developing NCBLs and Abedin et al. (Reference Abedin, Tsuji and Hattori2009) for a temporally developing NCBL; (b) compensated Nusselt number development for the temporally developing NCBL in the self-similar turbulent region.

The $\textit {Nu}_\delta$ predictions made from the ODE based, power-law solution (4.27), show excellent agreement with the data from both DNS-A and DNS-C, indicating the wall heat flux is indeed controlled by the characteristic scales from the plume-like region. Given by the power-law solution, (4.22), and making use of (4.24) and (4.27), the Nusselt number follows

(4.30)\begin{equation} \textit{Nu}_\delta \propto \textit{Gr}_\delta^{ 3\xi/(3\xi+3)}\approx 0.0408\textit{Gr}_\delta^{0.381}. \end{equation}

This $0.381$-power-law relation is consistent with the ultimate heat transfer scaling with a logarithmic correction reported by Ng et al. (Reference Ng, Ooi, Lohse and Chung2017), who obtained an empirical correlation

(4.31)\begin{equation} \textit{Nu}_\delta \propto \textit{Gr}_\delta^{0.37}, \end{equation}

after conditionally averaging the high wall shear events in their DNS study of an NCBL in a differentially heated vertical slot. Ng et al. (Reference Ng, Ooi, Lohse and Chung2017) suggest that the slightly higher $\textit {Nu}_\delta$–$\textit {Gr}_\delta$ power-law correlation in (4.31) is consistent with the ultimate turbulent regime power-law scaling with a logarithmic correlation. The exponent is also close to an experimental study by He et al. (Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012), who gives 0.38, for Rayleigh–Bénard convection at high Rayleigh number. A higher exponent has been reported by Eckert & Jackson (Reference Eckert and Jackson1950) – although the mean profiles were modelled by a $1/7$-power-law, they showed that the Nusselt number would follow $\textit {Nu}_\delta \propto \textit {Gr}_\delta ^{0.429}$ if the boundary layer is turbulent in the sense of Prandtl and von Kármán (i.e. momentum-dominated boundary layers).

4.6. Logarithmic correction for the ultimate turbulent regime

The Nusselt number correlation given by (4.30) suggests that the heat transfer of the temporally developing NCBL is driven by the bulk thermal convection in the plume-like region and limited by near-wall processes as seen in Rayleigh–Bénard convection. In the so-called ultimate turbulent regime, the boundary layer is dominated by the bulk thermal convection in the sense of Grossmann & Lohse (Reference Grossmann and Lohse2000,  Reference Grossmann and Lohse2011), and the $\textit {Nu}_\delta$ follows

(4.32)\begin{equation} \textit{Nu}_\delta \propto \textit{Re} \mathcal{L}(\textit{Re}), \end{equation}

as in homogeneous Rayleigh–Bénard convection, where $\textit {Re}$ is the bulk Reynolds number, and $\mathcal {L}(\textit {Re})$, being a function of Reynolds number, is a logarithmic-correction term for the turbulent kinetic boundary layer. The Reynolds number $\textit {Re}$ is then related to the Rayleigh number $\textit {Ra}$, giving $\textit {Re}\propto \textit {Ra}^{1/2}$ in the ultimate turbulent regime, for their Rayleigh–Bénard convection. It is noted in the vertical NCBL flow the gravity vector is perpendicular to the driving heat flux, whereas in the Rayleigh–Bénard convection the gravity is aligned with the driving heat flux. In our previous study (Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020), we have shown that the mean profiles follow a buoyancy-modified law of the wall for the NCBL (see also in figure 7). Following Grossmann & Lohse (Reference Grossmann and Lohse2011), the temperature log-law can be rewritten into

(4.33)\begin{equation} \frac{\theta_\tau}{\theta_w - \theta_{log}} = \frac{K_\theta}{\ln\left(\hat{\textit{Re}}\dfrac{u_\tau}{\hat{u}}e^{B_\theta K_\theta}\right)}, \end{equation}

where $\theta _{log}$ denotes the temperature at the end of the log-law region, $K_\theta$ is the thermal von Kármán constant, given by the ratio of von Kármán constant $K_v$ to the turbulent Prandtl number $\textit {Pr}_t$, $B_\theta$ is the constant from the temperature log-law, and the bulk Reynolds number $\hat {\textit {Re}}$ is given by the velocity scale in the plume-like region $\hat {u}$ and the length of the log-law $l_m$, as follows:

(4.34)\begin{equation} \hat{\textit{Re}} \equiv \frac{\hat{u}l_m}{\nu}. \end{equation}

In the present study, we assume the log-law extends up to the maximum velocity location (see also in Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020)), i.e. $l_m\approx \delta _i$. Making a substitution of (4.6) and (4.33) into the Nusselt number definition, (4.26), yields

(4.35)\begin{equation} \textit{Nu}_\delta=\frac{\hat{\textit{Re}}\textit{Pr}\dfrac{u_\tau}{\bar{u}_m}\dfrac{\bar{u}_m}{\hat{u}}K_\theta} {\ln\left(\hat{\textit{Re}}\dfrac{u_\tau}{\bar{u}_m}\dfrac{\bar{u}_m}{\hat{u}}\,\textrm{e}^{B_\theta K_\theta}\right)}\frac{\delta}{\delta_i} \frac{(\theta_w-\theta_{log})}{(\theta_w-\theta_\infty)}. \end{equation}

For asymptotically high $\textit {Gr}_\delta$ NCBL flows, the near-wall buoyancy effect would become negligible and the velocity is expected to follow a similar log-law to the momentum dominated flows (Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020), in the form of

(4.36)\begin{equation} \frac{u_\tau}{\bar{u}_m} = \frac{K_v}{\ln\left(\hat{\textit{Re}}\dfrac{u_\tau}{\bar{u}_m}\dfrac{\bar{u}_m}{\hat{u}}e^{B_u K_v}\right)}, \end{equation}

where $B_u$ is a constant for the velocity log-law. The Nusselt number scaling (4.35) can be further simplified into

(4.37)\begin{equation} \textit{Nu}_\delta=\hat{\textit{Re}}\textit{Pr} \mathcal{L}(\hat{\textit{Re}}), \end{equation}

which takes the same form as (4.32) given by Grossmann & Lohse (Reference Grossmann and Lohse2011), with the logarithmic correction term

(4.38)\begin{equation} \mathcal{L}(\hat{\textit{Re}}) = \frac{\mathcal{B}}{\ln\left(\hat{\textit{Re}}\dfrac{u_\tau}{\bar{u}_m}C_u\right) \ln\left(\hat{\textit{Re}}\dfrac{u_\tau}{\bar{u}_m}C_\theta\right)}, \end{equation}

and the coefficients $\mathcal {B}$, $C_u$ and $C_\theta$ are given by

(4.39a–c)\begin{equation} \mathcal{B} =K_v K_\theta \frac{\delta \bar{u}_m(\theta_w-\theta_{log})}{\delta_i \hat{u}(\theta_w-\theta_\infty)}, \quad C_u = \frac{\bar{u}_m}{\hat{u}} e^{K_v B_u}, \quad C_\theta = \frac{\bar{u}_m}{\hat{u}}e^{B_\theta K_\theta}. \end{equation}

Equation (4.37) indicates that $\textit {Nu}_\delta$ scales with the bulk Reynolds number $\hat {\textit {Re}}$ in the plume-like region with the logarithmic correction (4.38).

The logarithmic correction term (4.38), however, depends on the bulk Reynolds number $\hat {\textit {Re}}$, the velocity ratio ${u_\tau }/{\bar {u}_m}$ and the coefficients $\mathcal {B}$, $C_u$ and $C_\theta$. The velocity ratio ${u_\tau }/{\bar {u}_m}$ is implicitly given by the continuous fraction (4.36) and has analytical solution (Grossmann & Lohse Reference Grossmann and Lohse2011)

(4.40)\begin{equation} \frac{u_\tau}{\bar{u}_m} = \frac{K_v}{\omega(\hat{\textit{Re}}K_v C_u)}, \end{equation}

where $\omega (\phi )$ denotes the Wright omega function of $\phi$. In the present study, the bulk velocity and length scales are found to scale with the characteristic scales for the plume-like region, i.e. $\bar {u}_m/\hat {u}\approx 1.475$ and $\delta _i/ \delta \approx 0.1035$ are approximately constant in the current $\textit {Gr}_\delta$ range, as shown in figure 18(a,b). The $\delta _i/\delta$ value reported here is also close to the measurement of Tsuji & Nagano (Reference Tsuji and Nagano1988) and Abedin et al. (Reference Abedin, Tsuji and Hattori2009) at lower Grashof number range, which is $\delta _i/\delta = 0.07$. The temperature ratio $(\theta _w-\theta _{\textit {log}})/(\theta _w-\theta _\infty )$, on the other hand, is thought to be varying slowly since the temperature difference across the logarithmic region would increase slowly with $\textit {Gr}_{\delta }$ (we emphasise the very large $\textit {Gr}_{\delta }$ range here). As shown in figure 18(c), the temperature ratio $(\theta _w-\theta _{\textit {log}})/(\theta _w-\theta _\infty )$ demonstrates weak $\textit {Gr}_{\delta }$-dependency in the current $\textit {Gr}_{\delta }$ range, therefore, it appears reasonable to assume a constant temperature ratio in our $\textit {Gr}_{\delta }$ range. This suggests that $\mathcal {B}$, $C_\theta$ and $C_u$ are all approximately constant. With substitution of (4.40) into (4.38), the logarithmic correction term depends solely on the bulk Reynolds number $\hat {\textit {Re}}$. Unfortunately, the $\hat {\textit {Re}}$-dependency of (4.38) is difficult to obtain analytically. However, since $K_v$, $K_\theta$ and the bulk scale ratios ($\bar {u}_m/\hat {u}$ and $\delta /\delta _i$) are all approximately constants in the current Grashof number range Grashof number range (values tabulated in table 5), we are able to calculate (4.38) numerically. The local effective scaling exponent $\mathrm {d} \log \mathcal {L}(\hat {\textit {Re}})/\mathrm {d} \log \hat {\textit {Re}}$, as shown in figure 19(a), is found to be varying slowly with $\hat {\textit {Re}}$. The gradually reducing magnitude of the effective scaling exponent indicates the $\hat {\textit {Re}}$-dependency becomes weaker with increasing $\hat {\textit {Re}}$. We stress that (4.36) and (4.40) are approximations for asymptotically high $\textit {Gr}_{\delta }$ flow since the buoyancy effect (which will vary with $\textit {Gr}_{\delta }$ and $\hat {\textit {Re}}$) is not accounted for (Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020), and only at this asymptotic $\textit {Gr}_{\delta }$ range will these relations hold. Nevertheless, we present figure 19(a), based on the asymptotic relations (4.36) and (4.40), for illustrative purposes, as the general behaviour of the logarithmic correction term is the same.

Figure 18.

Temporal development of the scale ratios: (a) $\bar {u}_m/\hat {u}$; (b) $\delta _i/\delta$; (c) $(\theta _w-\theta _{\textit {log}})/(\theta _w-\theta _\infty )$. Horizontal black dash–dotted lines indicates the temporal average: $\bar {u}_m/\hat {u}=1.475$; $\delta _i/\delta =0.1035$; $(\theta _w-\theta _{\textit {log}})/(\theta _w-\theta _\infty )=0.735$.

Making use of (4.25) and $\delta _i\propto \delta$ in the current range of Grashof number, the bulk Reynolds number then scales with the ODE solution,

(4.41)\begin{equation} \hat{\textit{Re}} \propto \hat{u}\hat{l} \propto \textit{Gr}_{\delta}^{(2\xi+1)/(3\xi+3)}. \end{equation}

Based on the values tabulated in table 4, (4.41) is equivalent to a $\hat {\textit {Re}}\approx 0.39 \textit {Gr}_{\delta }^{0.46}$ scaling for the current Grashof number range, close to the $\hat {\textit {Re}}\propto \textit {Gr}_{\delta }^{1/2}$ as suggested by Grossmann & Lohse (Reference Grossmann and Lohse2011) for Rayleigh–Bénard flow. With (4.41), it is also possible to obtain the development of a local effective scaling exponent in terms of Grashof number ($\mathrm {d} \log \mathcal {L}(\hat {\textit {Re}})/\mathrm {d} \log \textit {Gr}_\delta$), as shown in figure 19(b). Again we note that the mean velocity profile in the NCBL at non-asymptotic Grashof number is adjusted by buoyancy and takes a more complicated log-law (cf. Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020), the extra buoyancy terms) than (4.36), resulting in a buoyancy-modified logarithmic correction term.

In the present study, the local effective scaling exponent for the logarithmic correction term is obtained empirically by

(4.42)\begin{equation} \mathcal{L}(\hat{\textit{Re}}) = \frac{\textit{Nu}_\delta}{\hat{\textit{Re}}\textit{Pr}}, \end{equation}

using DNS data. As shown in figure 20, for the current $\textit {Gr}_{\delta }$ range the scaling exponent is found to be $-0.08$. In addition, the dimensional consistency between (4.30) and (4.37) further suggests

(4.43)\begin{equation} \mathcal{L}(\hat{\textit{Re}}) \propto \frac{\textit{Gr}_{\delta}^{3\xi/(3\xi+3)}}{\textit{Gr}_{\delta}^{(2\xi+1)/(3\xi+3)}} = \textit{Gr}_{\delta}^{(\xi-1)/(3\xi+3)}. \end{equation}

For the current range of $\textit {Gr}_{\delta }$, $\xi =0.6156$ and the scaling exponent is found to be $-0.079$, consistent with the empirical results shown in figure 20. This indicates the Nusselt number is given by

(4.44)\begin{equation} \textit{Nu}_\delta\propto \hat{\textit{Re}}\textit{Pr} \mathcal{L}(\hat{\textit{Re}}) \propto \textit{Gr}_\delta^{0.46} \textit{Gr}_\delta^{{-}0.08} \approx \textit{Gr}_\delta^{0.38}. \end{equation}

This effective $\textit {Nu}_\delta$–$\textit {Gr}_\delta$ scaling exponent, higher than the empirical $1/3$-power-law, is consistent with our $0.381$-power-law from the ODE prediction and DNS regression results in the current Grashof number range.

Table 5.

Constants for the mean velocity and temperature log-law.

Figure 19.

The temporal development of the local effective scaling exponent for the logarithmic correction term; panel (b) shows same information as panel (a) but with Grashof number. The Grashof number $\textit {Gr}_\delta$ is obtained by (4.41).

Figure 20.

The compensated logarithmic correction scaling.

However, as depicted by figure 19, over a large range of $\textit {Gr}_\delta$ the logarithmic correction $\mathcal {L}(\hat {\textit {Re}})$ is expected to vary very slowly with $\textit {Gr}_{\delta }$ (or $\hat {\textit {Re}}$), and would eventually approach to some constant at asymptotically high $\textit {Gr}_{\delta }$. This further implies that the exponent $\xi$ and, therefore, the profile coefficients and entrainment coefficient $E$ for the integral model, are also slowly varying quantities with $\textit {Gr}_{\delta }$ – again we emphasise the very large $\textit {Gr}_{\delta }$ range here. At asymptotically high $\textit {Gr}_{\delta }$, $\xi$ approaches to unity since (4.43) would become independent of Grashof number, and the Nusselt number would follow

(4.45)\begin{equation} \textit{Nu}_\delta\propto \hat{\textit{Re}}\textit{Pr} \mathcal{L}(\hat{\textit{Re}}) = \textit{Gr}_{\delta}^{1/2} \textit{Gr}_{\delta}^0 = \textit{Gr}_{\delta}^{1/2}, \end{equation}

consistent with the $1/2$-power-law proposed by Wells & Worster (Reference Wells and Worster2008) and Grossmann & Lohse (Reference Grossmann and Lohse2011) for asymptotically high Grashof number. We stress that, similar to the logarithmic correction term, all the coefficients are very slowly varying, and can be approximately taken as constants for the current Grashof number range ($10^7<\textit {Gr}_\delta <1.2\times 10^8$). Nevertheless, making use of the values given in table 5, we are able to predict the Nusselt number based on (4.37). As demonstrated by figure 21, (4.37) still shows good agreement with the DNS data.

Figure 21.

Nusselt number prediction using the logarithmic correction scaling (4.37).

This idea is further justified by identifying the local heat transfer in the plume-like part of the flow separately from the logarithmic region. The heat transfer within the plume-like region can be expressed in terms of a local Nusselt number and a local Grashof number,

(4.46a,b)$$\begin{gather} \textit{Nu}_{\hat{b}} \equiv \frac{q_{\delta_i} \hat{l}}{\kappa\hat{b}}\approx\frac{q_{w} \hat{l}}{\kappa\hat{b}}\propto \frac{\hat{u}\hat{l}}{\kappa}, \quad \textit{Gr}_{\hat{b}} \equiv \frac{ \hat{b}\hat{l}^3}{\nu^2}. \end{gather}$$

With the ODE power-law solution (4.22), it immediately follows that

(4.47)\begin{equation} \textit{Nu}_{\hat{b}} \propto \textit{Gr}_{\hat{b}}^{(2\xi+1)/(4\xi+2)} =\textit{Gr}_{\hat{b}}^{1/2}. \end{equation}

This $1/2$-power-law scaling is confirmed by our DNS data, as shown in figure 22, where the heat transfer into the plume-like region clearly demonstrates the scaling given in (4.47) since the near-wall logarithmic correction is absent in this region, indicating the plume-like region is dominated by the turbulent mixing in the sense of Grossmann & Lohse (Reference Grossmann and Lohse2000,  Reference Grossmann and Lohse2001).

Figure 22.

The local Nusselt number development in the plume-like region.

In the limit of asymptotically high Grashof number, (4.45) suggests that the Nusselt number would follow $\textit {Nu}_\delta \propto \textit {Gr}_{\delta }^{1/2}$ since the logarithmic region is so large and the temperature across this logarithmic region increases at an asymptotically small rate. Similar results have been conjectured by Wells & Worster (Reference Wells and Worster2008), who used empirical correlations for the forced convection, suggesting at very high Grashof number the wall shear stress and heat flux follows

(4.48a)\begin{gather} \tau_w = K_f\hat{u}^2, \end{gather}

(4.48b)\begin{gather}q_w=\textit{St} \hat{u}\Delta T_{BL}, \end{gather}

where $\Delta T_{BL} = T_{w}-\hat {b}/g\alpha$ is the near-wall temperature difference between the isothermal wall and the plume-like region, $K_f$ is the skin friction coefficient, $\textit {St}$ is the Stanton number. Based on the experimental measurements of Tsuji & Nagano (Reference Tsuji and Nagano1988), Wells & Worster (Reference Wells and Worster2008) suggest $K_f=0.0046$ and $\textit {St} = 0.0035$ are both constants. The Reynolds analogy allows (4.48) to be rewritten into

(4.49)\begin{equation} \frac{q_w}{(T_{w}-\hat{b}/g\alpha)\hat{u}} \propto \frac{\textit{Nu}_\delta}{\textit{Re}_{\delta}\textit{Pr}} \propto 2K_f, \end{equation}

where the Reynolds number is given by

(4.50)\begin{equation} \textit{Re}_{\delta} \equiv \frac{\hat{u}\delta}{\nu}. \end{equation}

Reynolds analogy has also been adopted by Eckert & Jackson (Reference Eckert and Jackson1950), where they assumed the NCBL shares a similar near-wall shear and heat transfer characteristics to those of a forced convection boundary layer, resulting a $0.429$-power-law scaling for the Nusselt number. Rearranging (4.49) yields

(4.51)\begin{equation} \textit{Nu}_\delta \propto 2K_f \frac{\delta}{l_m} \hat{\textit{Re}}\textit{Pr}. \end{equation}

Equation (4.51) is, in general, consistent with the ultimate heat transfer correlation (4.37) suggested by Grossmann & Lohse (Reference Grossmann and Lohse2000,  Reference Grossmann and Lohse2001,  Reference Grossmann and Lohse2011): for a hydraulically smooth surface, $K_f$ is usually a function of $\textit {Re}_{\delta }$ (cf. Eckert & Jackson Reference Eckert and Jackson1950), i.e. $K_f(\textit {Re}_{\delta })$; while a constant $K_f$, as suggested by Wells & Worster (Reference Wells and Worster2008), would indicate a hydraulically rough surface.

4.7. Wall shear stress characteristics

Unlike the Nusselt number, the wall shear stress $\tau _w$ is not directly linked to the characteristic scales for the self-similar model. Given by (4.4), the wall shear stress is driven by the shear exerted by the plume-like region $\tau _{\delta _i}$ and the total buoyancy supplied by the isothermal wall in the near-wall region $B_i$ (see also in figure 6b),

(4.52)\begin{equation} \tau_w = \tau_{\delta_i} + B_i, \end{equation}

where the near-wall total buoyancy $B_i$ is given by

(4.53)\begin{equation} B_i \equiv \int^{\delta_i}_0 \bar{b}\,\mathrm{d} y. \end{equation}

In the present study, the shear stress on the BL–plume interface is found to scale with the characteristic scales in the plume-like region,

(4.54)\begin{equation} \tau_{\delta_i} = S_{\delta_i}\hat{u}^2, \end{equation}

with the profile coefficient $S_{\delta _i}$ being approximately a constant in the current Grashof number range. The near-wall buoyancy, however, can be expressed by

(4.55)\begin{equation} B_i = B_t - B, \end{equation}

where $B$ is the total buoyancy within the plume-like region, as given in (3.3), and $B_t$ is the total buoyancy for the entire flow, i.e.

(4.56)\begin{equation} B_t \equiv \int^{\infty}_0 \bar{b} \,\mathrm{d} y. \end{equation}

Although (4.55) and (4.56) are not directly linked to the characteristic scales, the total buoyancy of the flow accumulates with time due to the wall heat flux,

(4.57)\begin{equation} \frac{\partial B_t}{\partial t} = q_w. \end{equation}

Since the BL–plume interface lies in the constant heat flux region, with (4.7), (4.57) can be simplified into

(4.58)\begin{equation} \frac{\partial B_t}{\partial t} =q_w\approx q_{\delta_i} =\beta_{\delta_i}\hat{u}\hat{b}. \end{equation}

The temporal evolution of the near-wall buoyancy is, therefore, given by

(4.59)\begin{equation} \frac{\partial B_i}{\partial t}=\frac{\partial (B_t - B)}{\partial t}=\frac{\partial B_t}{\partial t}-\frac{\partial B}{\partial t} = \mathcal{L}_B \hat{u}\hat{b}. \end{equation}

Here the Leibniz term $\mathcal {L}_B$ also appears constant in the current Grashof number range (see table 2). As shown in (4.22), the dimensional consistency of the integral model also suggests that

(4.60a,b)\begin{equation} \hat{u} =A_{\hat{u}} t^\xi,\quad \hat{b}= A_{\hat{b}} t^{\xi-1}. \end{equation}

Equation (4.59) can be integrated into

(4.61)\begin{equation} B_i =\frac{\mathcal{L}_BA_{\hat{u}}A_{\hat{b}}}{2\xi} t^{2\xi} + C_B = K_B \hat{u}^2+C_B, \end{equation}

where $K_B = \mathcal {L}_BA_{\hat {b}}/(2\xi A_{\hat {u}})$ is a proportionality coefficient and $C_B$ is an integral constant.Figure 23 demonstrates the temporal development of the near-wall buoyancy $B_i$ using DNS data and a prediction made from (4.61) based on the constant coefficients given in tables 2 and 4. The integral constant is empirically found to be $C_B=5.5$ by fitting (4.61) to the DNS data.

Figure 23.

The temporal development of the total buoyancy in the near-wall BL-like region $B_i$.

Figure 24.

The compensated wall shear stress scaling (4.62).

Equations (4.52), (4.54) and (4.61) indicate that the wall shear stress of the NCBL follows

(4.62)\begin{equation} \tau_w - C_B = (S_{\delta_i}+K_B)\hat{u}^2, \end{equation}

as shown in figure 24. At asymptotically high Grashof number (or, Reynolds number), the wall shear stress becomes large and

(4.63)\begin{equation} \tau_w-C_B \approx \tau_w. \end{equation}

This suggests that the wall shear stress, at asymptotically high Grashof number, scales with the bulk velocity scale

(4.64)\begin{equation} \tau_w \approx (S_{\delta_i}+K_B)\hat{u}^2, \end{equation}

as in the momentum-dominated flows.