Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution

Results on the Prandtl-Blasius type kinetic and thermal boundary layer thicknesses in turbulent Rayleigh-B\'enard convection in a broad range of Prandtl numbers are presented. By solving the laminar Prandtl-Blasius boundary layer equations, we calculate the ratio of the thermal and kinetic boundary layer thicknesses, which depends on the Prandtl number Pr only. It is approximated as $0.588Pr^{-1/2}$ for $Pr\ll Pr^*$ and as $0.982 Pr^{-1/3}$ for $Pr^*\ll\Pr$, with $Pr^*= 0.046$. Comparison of the Prandtl--Blasius velocity boundary layer thickness with that evaluated in the direct numerical simulations by Stevens, Verzicco, and Lohse (J. Fluid Mech. 643, 495 (2010)) gives very good agreement. Based on the Prandtl--Blasius type considerations, we derive a lower-bound estimate for the minimum number of the computational mesh nodes, required to conduct accurate numerical simulations of moderately high (boundary layer dominated) turbulent Rayleigh-B\'enard convection, in the thermal and kinetic boundary layers close to bottom and top plates. It is shown that the number of required nodes within each boundary layer depends on Nu and Pr and grows with the Rayleigh number Ra not slower than $\sim\Ra^{0.15}$. This estimate agrees excellently with empirical results, which were based on the convergence of the Nusselt number in numerical simulations.

not considerably change with increasing grid resolution [16,17,18,19,20,21] or by guaranteeing (e.g. in ref. [22,21]) that the Nusselt numbers calculated from the global energy dissipation rate or thermal dissipation rate well agree with that one calculated from the temperature gradient at the plates or the ones obtained from the overall heat flux. The knowledge that the profiles are of Prandtl-Blasius type offers the opportunity to a priori determine the number of required grid points in the BLs for given Rayleigh number and Prandtl number, valid in the boundary layer dominated ranges of moderately high Ra numbers.
In section 2 we will first revisit the Prandtl-Blasius BL theory -see refs. [9,10,11,12,13,7] or for more recent discussions in the context of RB refs. [6,23] -and derive the ratio between the thermal boundary layer thickness δ θ and the velocity boundary layer thickness δ u as functions of the Prandtl number Pr extending previous work (section 3). We will also discuss the limiting cases for large and small Pr, respectively. The transitional Prandtl number between the two limiting regimes turns out to be surprisingly small, namely Pr * = 0.046. The crossover range is found to be rather broad, roughly four orders of magnitude in Pr. In section 4 we note that the Prandtl-Blasius velocity BL thickness is different from the velocity BL thickness based on the position of the maximum r.m.s. velocity fluctuations (widely used in the literature), but well agrees with a BL thickness based on the position of the maximum of an energy dissipation derivate that was recently introduced in ref. [24,21]. We then derive the estimate for the minimum number of grid points that should be placed in the boundary layers close the top and bottom plates, in order to guarantee proper grid resolution. Remarkably, the number of grid points that must have a distance smaller than δ u from the wall increases with increasing Ra, roughly as ∼ Ra 0. 15 . This estimate is compared with a posteriori results for the required grid resolution obtained in various DNSs of the last three decades, finding good agreement. Section 5 is left to conclusions.

Prandtl boundary layer equations
The Prandtl-Blasius boundary layer equations for the velocity field u(x, z) (assumed to be two-dimensional and stationary) over a semi-infinite horizontal plate [9,10,11,12,13,7] read with the boundary conditions u x (x, 0) = 0, u z (x, 0) = 0, and u x (x, ∞) = U . Here u x (x, z) is the horizontal component of the velocity (in the direction x of the largescale circulation), u z (x, z) the vertical component of the velocity (in the direction z perpendicular to the plate), and U the horizontal velocity outside the kinetic boundary layer (wind of turbulence). Correspondingly, the equation determining the (stationary) temperature field T (x, z) reads with the boundary conditions T (x, 0) = T plate and T (x, ∞) = T bulk , which under Oberbeck-Boussinesq conditions is the arithmetic mean of the upper and lower plate temperature. Applying these equations to RB flow implies that we assume the temperature field to be passive. The dimensionless similarity variable ξ for the vertical distance z from the plate measured at the distance x from the plate's edge is Since the flow in Prandtl theory is two-dimensional, a streamfunctionΨ can be introduced, which represents the velocity field. The streamfunction is nondimensionalized as Ψ =Ψ/ √ xνU , and the temperature is measured in terms of ∆/2, giving the non-dimensional temperature field Θ. Rewriting eqs. (1) and (2) in terms of Ψ and Θ one obtains Here the boundary conditions are The temperature and velocity profiles obtained from numerically solving equations (4)-(7) (for particular Prandtl numbers) are already shown in textbooks [12,7,13] and in the context of RB convection in refs. [23,25]: From the momentum equation (6) with above boundary conditions one immediately obtains the horizontal velocity dΨ/dξ. The dimensionless kinetic boundary layer thicknessδ u can be defined as that distance from the plate at which the tangent to the function dΨ/dξ at the plate (ξ = 0) intersects the straight line dΨ/dξ = 1 (see figure 1 a). As equation (4) and the boundary conditions (6) contain no parameter whatsoever, the dimensionless thicknessδ u of the kinetic boundary layer with respect to the similarity variable ξ is universal, i.e., independent of Pr and U or Re,δ Solving numerically equation (5) with the boundary conditions (7) for any fixed Prandtl number, one obtains the temperature profile with respect to the similarity variable ξ (see figure 1 b). Note that in contrast to the longitudinal velocity dΨ/dξ, the temperature profile Θ depends not only on ξ but also on the Prandtl number, since Pr appears in equation (5) as the (only) parameter. The distance from the plate at which the tangent to the Θ profile intersects the straight line Θ = 1 defines the dimensionless thickness of the thermal boundary layer, where C(Pr) is a certain function of Prandtl number. E.g., one numerically finds C ≈ 3.417, 1.814, and 1.596 for Pr = 0.7, 4.38, and 6.4, respectively (see figure 1 b).
From (8) and (9) one obtains the ratio between the (dimensional) thermal boundary layer thickness δ θ and the (dimensional) kinetic boundary layer thickness δ u : As discussed above, the constant A and the function C = C(Pr) are found from the solutions of equations (4)-(7) for different Pr. A and C(Pr) reflect the slopes of the respective profiles, With (3) the physical thicknesses are δ u =δ u / U xν and δ θ =δ θ / U xν , generally depending on U and the position x along the plate. The physical thermal BL thickness then is Thus, explicitly it depends neither on U nor on the position x along the plate. Reminding the definition of the thermal current J = u z T −κ∂ z T , we get ∂Θ , on x-average we have δ θ is the so-called slope thickness, see Sect. 2.4 of reference [23]. In contrast to the thermal BL thickness δ θ the physical velocity BL thickness δ u = A −1 / U xν depends explicitly both on the position x and on the wind amplitude U . In a Rayleigh-Bénard cell we choose for x a representative value x =ãL =ãΓH. Then the famous Prandtl formula [9] results Here a = ãΓ The constant a has been obtained empirically [5], based on the experimental measurements by [26] performed in a cylindrical cell of aspect ratio one, filled with water. The result was [5] a ≈ 0.482.
We note that this value probably depends on the aspect ratio, on the shape of the RB container, and can also be different for numerical 2D Rayleigh-Bénard convection [27,28,29]. It will also be different for the slope thickness as considered here or other definitions as e. g. the 99% -thickness. It seems worthwhile to note that similarly to the case of δ θ also δ u can be expressed by a profile slope at the plate. Analogously to the temperature case one calculates for the kinetic thickness δ u = U/ ∂ux ∂z (0). Here U appears explicitly and the derivative may depend on x. The denominator is the local stress tensor component, whichafter averaging -describes the momentum transport, just as the temperature profile derivative at the plate characterises the heat transport. In combination with eq. (14) it says that the kinetic stress behaves as ∂ux ∂z (0) ∼ U √ Re/(aH). From eqs. (10) and (14) we also find the useful (and known) relation for Prandtl-Blasius boundary layers From solving equations (4)-(7) together with relations (11) one obtains that the BL thickness ratio (10) has two limiting cases, namely δ θ /δ u ∼ Pr −1/2 for very small Pr 1 and δ θ /δ u ∼ Pr −1/3 for very large Pr 1. We thus present the ratio of the thermal and kinetic boundary layer thicknesses normalised by Pr −1/3 in figure 2 for different Pr from Pr = 10 −6 to 10 6 . The figure confirms that the scaling of the ratio between the thermal and kinetic boundary layer thicknesses in the low and high Prandtl number regimes is Pr −1/2 and Pr −1/3 , respectively. Between these two limiting regimes there is a transition region, whose width is about 4 orders of magnitude in Pr. In the next section we will derive analytic expressions for the ratio δ θ /δ u in the respective regimes, which will be used in the remainder of the paper to analyse the resolution properties for DNS in the BLs of the Rayleigh-Bénard system. In the Prandtl-Blasius theory the asymptotic velocity amplitude U is a given parameter; the resulting heat current Nu is a performance of the boundary layers only. In contrast, in Rayleigh-Bénard convection the heat transport is determined by the BLs together with the bulk flow. Therefore in RB convection the wind amplitude U no longer is a passive parameter, but U and Nu are actively coupled properties of the full thermal convection process.
The Reynolds number Re is defined as the dimensionless wind amplitude, From the law for the kinetic BL thickness (14), the thermal BL thickness δ θ (13), and the BL thickness ratio (10) one obtains This Re ∼ Nu 2 law is in perfect agreement with the GL theory [3,4,5,6]. In that theory several sub-regimes in the (Ra, Pr) parameter space are introduced, depending on the dominance of the BL or bulk contributions. In regimes I and II the BL of the temperature field dominates, while in III and VI it is the thermal bulk. Regimes I and II differ in the velocity field contributions: It either is the u-BL (I) or the u-bulk (II) which dominates; analogously the pair III and IV is characterized. The labels (for lower Pr) and u (for upper Pr) distinguish the cases in which the thermal BL is thicker or smaller than the kinetic one. All ranges in the GL theory, which are thermal boundary layer dominated, show the Re ∼ Nu 2 behaviour, namely I l , I u , II l , II u . In the thermal bulk dominated ranges of RB convection the relation between Re and Nu is different.
In III u we have Re ∼ Nu 4/3 , in IV l it is Re ∼ Nu, and in IV u also Re ∼ Nu 4/3 holds; but here the Prandtl-Blasius result (18) is not applicable, since the heat transport mainly depends on the heat transport properties of the bulk. In the range I ∞ , although boundary layer dominated, also a different relation (Re ∼ Nu 3 ) holds; here the upper and the lower kinetic BLs fill the whole volume and therefore there is no free flow outside the BLs, in contrast to the Prandtl-Blasius assumption of an asymptotic velocity with the LSC amplitude U .
3. Approximations for the ratio δ θ /δ u of the temperature and velocity boundary layer thicknesses In this section we will derive analytical approximations for the ratio δ θ /δ u for the three regimes identified in the previous section, cf. figure 2. We start by discussing the low (Pr < 3 × 10 −4 ) and the high (3 < Pr) Prandtl number regimes, before we discuss the transition region 3 × 10 −4 ≤ Pr ≤ 3.

Approximation of
In the case of very small Prandtl number, Pr 1, the thickness of the velocity boundary layer is negligible compared with the thickness of the temperature boundary layer, i.e., δ θ δ u . Hence, in most of the thermal boundary layer it is u x ≈ U . Introducing the similarity variable as in ref. [13] one obtains the following equation for the temperature as a function of η: The solution of this boundary value problem is the Gaussian error function According to (3) and (19), the similarity variable ξ used in the Prandtl equations and the similarity variable η used in the approximation for Pr 1 are related as follows Applying now the formulae (20), (21) and (11) we obtain the following equalities: This leads to the approximation for the function C(Pr) = √ πPr −1/2 for very small Pr.
In figure 2 one can see that for very small Prandtl numbers, Pr < 3 × 10 −4 , the approximation (22) is as expected indistinguishable from the numerically obtained δ θ /δ u .
3.2. Approximation of δ θ /δ u for Pr > 3 Meksyn [12], based on the work by Pohlhausen [11], derived that the solution of the temperature equation (5), together with relation (7) equals The constant D can be found as usual from the boundary condition at infinity and was approximated in [11,12] for Pr > 1 as follows
In figure 2 the approximation (24) is presented as a blue dash-dotted curve. For Pr > 3 the function (δ θ /δ u )Pr 1/3 almost coincides with the constant E.
3.3. Approximation of δ θ /δ u in the crossover range 3 × 10 −4 ≤ Pr ≤ 3 As one can see in figure 2, the approximation (22) well represents δ θ /δ u in the region Pr < 3×10 −4 , while (25) is a good approximation of δ θ /δ u for Pr > 3. An approximation of the ratio of the thermal and kinetic boundary layer thicknesses in the transition region 3 × 10 −4 ≤ Pr ≤ 3 is obtained by applying a least square fit to the numerical solutions of the Prandtl-Blasius equations (4)-(7). One finds: As seen in figure 2, this relation is a good fit of the full solution in the transition regime.

Summary
For the ratio δ θ /δ u of the thicknesses of the thermal and kinetic boundary layers, which depends strongly (and only) on Pr, we find according to (22), (25), and (26) The crossover Prandtl number Pr * between the asymptotic behaviours, cf. first and last line of (27), is defined as the intersection point Pr * = 0.046 of the asymptotic approximations. Note that this crossover between the small-Pr behaviour δ θ /δ u ∝ Pr −1/2 and the large-Pr behaviour δ θ /δ u ∝ Pr −1/3 does not happen at a Prandtl number of order 1, but at the more than 20 times smaller value Pr * = 0.046. In this sense most experiments are conducted in the large Pr regime. However, also note that other definitions of the BL thicknesses lead to other crossover Prandtl numbers. Finally, we also give the thickness of the kinetic BL in the three regimes, as obtained from (27) and (13), namely We compare this Prandtl-Blasius result (28) for the kinetic boundary layer thickness in terms of Nu and Pr (thus valid if the heat transport is BL dominated) with the estimate given in reference [21], where the kinetic boundary layer thickness in a cylindrical cell is identified as two times that height at which the averaged quantity " u := u · ∇ 2 u t,φ,r has a maximum, because it was empirically found that the maximum of " u is approximately in the middle of the velocity boundary layer. Here u is the velocity field and the averaging is over time t, the azimuthal direction φ, and over the radial direction 0.1R < r < 0.9R, with R the radius of the cylindrical convective cell. The restricted range for the radial direction has been used in order to exclude the singularity region close to the cylinder axis and the region close to the sidewall, where the definition misrepresents the kinetic boundary layer thickness. Figure 3 shows that there is a very good agreement between the theoretical Prandtl-Blasius slope boundary layer thickness and that obtained using (29). The figure also shows that the position of the maximum r.m.s. velocity fluctuations is not a good indicator for the velocity boundary layer edge; it rather seems to identify the position where the LSC is the strongest.

Resolution requirements within the boundary layers in DNS
We now come to the main point of the paper: What can we learn from the Prandtl-Blasius theory for the required mesh resolution in the BLs of DNS of turbulent RB convection? Obviously, a "proper" mesh resolution should be used in order to obtain accurate results. In a perfect DNS the local mesh size should be smaller than the local Kolmogorov η K (x, t) and Batchelor η B (x, t) scales (see e.g. ref. [30]), and the resolution in the boundary layers should be also sufficient, see e.g. [31,16,25,32,21]. It indeed has been well established that the Nusselt number is very sensitive to the grid resolution used in the boundary layers; when DNS is underresolved, the measured Nusselt number is too high [31,33,16,34,35,36,21]. Hitherto, the standard way to empirically check whether the mesh resolution is sufficient is to try a finer mesh and to make sure that the Nusselt number is not too different. In this way the minimal number of grid points that is needed in the boundary layer is obtained by trial and error: Grötzbach [31] varied the number of grid points in the boundary layer between 1 and 5 in simulations up to Ra = 3 × 10 5 with Pr = 0.71 and found that 3 grid points in the boundary layers should be sufficient. Verzicco and Camussi [33] tested this at Ra = 2 × 10 7 and Pr = 0.7 and stated that at least 5 points should be placed in the boundary layers. Stevens et al. [21] tested the grid resolution for Ra = 2 × 10 6 to 2 × 10 11 and Pr = 0.7. They found that for Ra = 2 × 10 9 the minimum number of nodes in the boundary layers should be around 10 and that this number increases for increasing Ra. Together with the earlier series of papers the data clearly suggest that indeed there is an increase of required grid points in the BL with increasing Rayleigh number.
However, one must be careful. The empirical determination of the required number of grid points in the BL is not only intensive in computational cost, but also difficult.
The Nusselt number obtained in the simulations not only depends on the grid resolution in the BLs at the top and bottom plates, but also on the grid resolution in the bulk and at the side walls where the thermal plumes pass along [21]. So obviously a general theory-based criterion for the required grid resolution in the thermal and kinematic boundary layers will be helpful for performing future simulations. In this section we will derive such a universal criterion, harvesting above results from the Prandtl-Blasius boundary layer theory.
We first define the (local) kinetic energy dissipation rates per mass, Its time and space average for incompressible flow with zero velocity b.c. is u t,V = It is connected with the Nusselt number through the exact relation This follows directly from the momentum equation for Rayleigh-Bénard convection in Boussinesq approximation [37]. Here, · t,V denotes averaging over the whole volume of the convective cell and over time and (later) · t,A denotes averaging over any horizontal plane and time.
We start with the well established criterion that in a perfect DNS simulation the (local) mesh size must not be larger than the (local) Kolmogorov scale [38] η K (x, t), which is locally defined with the energy dissipation rate of the velocity, η K is the length scale at which the inertial term ∼ u 2 r /r and the viscous term ∼ νu r /r 2 of the Navier-Stokes equation balance, where u r ∼ ( u r) 1/3 has been assumed for the velocity difference at scale r. A corresponding length scale η T follows from the balance of the advection term ∼ u r T r /r and the thermal diffusion term κT r /r 2 in the advection equation; it is However, for large Pr the velocity field is smooth at those scales at which the temperature field is still fluctuating. Then the velocity difference u r ∼ u /νr and advection term and thermal diffusion term balance at the so-called Batchelor scale [39] η B , which is defined as For small Pr < 1 obviously η T > η B > η K and for comparison with the grid resolution, the Kolmogorov scale η K seems to be the most restrictive (i.e., smallest) length scale. In contrast, for large Pr > 1 it η T < η B < η K and one may argue that η T is the most restrictive length scale. This indeed may be the case in the Prandtl number regime in which the velocity field can still be described through Kolmogorov scaling u r ∼ ( u r) 1/3 , but for even larger Pr the velocity field becomes smooth u r ∼ u /νr and then the grid resolution should be compared to the Batchelor scale η B as smallest relevant length scale. In below analysis, for Pr > 1 we will restrict ourselves to this limiting case.
We now define global Kolmogorov and Batchelor length scales η global , respectively, (and also the global length scale η global ). Using the exact relation (31), one can find how the global Kolmogorov length η global K depends on Ra, Pr, and Nu, namely The admissible global mesh size h global should clearly be smaller than both η global K and η global B , which implies that one is on the safe side provided that H for Pr ≤ 1 (36) or with the relation (34) between the Kolmogorov and Batchelor length A similar way to estimate mesh requirements in the bulk was suggested for the first time by Grötzbach [31]. Note that with these estimates for the required bulk resolution for most times and locations one is on the safe side, as equation (31) is an estimate for the volume averaged energy dissipation rate, which is localized in the boundary layers. However, not only the background field but also plumes detaching from the boundary layers do require an adequate resolution.
To estimate the number of nodes that should be placed in the boundary layers, we will first estimate the area averaged energy dissipation rate in a horizontal plane in the velocity BL, u t,A∈BL . Employing eqs. (17), (14) and (30), one can find a lower bound for this quantity, namely From eqs. (31), (38), (18) and (27) For the Kolmogorov length η BL K in the velocity BL one can therefore write The mesh size h BL in the BL must be smaller than η BL K and η BL B , i.e., one is on the safe side if according to (39), (40), (36) and (37). From the relations (41), (27) and (13) one can estimate the minimum number of nodes of the computational mesh, which must be placed in each thermal and kinetic boundary layer close the plates. We find that this minimum number of nodes in the thermal boundary layers is while the minimum number of nodes in the kinetic boundary layers is The number of nodes in the thermal boundary layer looks very restrictive for very low Pr; however, one should realise that for very low Pr the thermal boundary layer also becomes much thicker than the velocity boundary layer. Hence, the criterion for the number of nodes in the thermal boundary layers determines the ideal distribution of nodes above the viscous boundary layer. For very high Pr the kinetic boundary layer becomes much thicker than the thermal boundary layer, and hence the restriction for the velocity boundary layer determines the ideal distribution of nodes above the thermal boundary boundary layer. Note that for large Pr equation (42) suggests that the number of grid points in the thermal boundary layer becomes independent of Pr (for fixed Nu): Indeed, as the velocity field is smooth anyhow, with increasing Pr no extra grid points are necessary in the thermal BL.
Note that the numerical pre-factors in these estimates significantly depend on the Prandtl number and on the empirically determined (ref. [5]) value of a, cf. eq. (15). The minimum node numbers for other values of Pr can be calculated directly using the relations (42)- (43). Apparently the scaling exponent depends much less on Pr. -All these estimates only give lower bounds on the required number of nodes in the boundary layers.
As discussed at the beginning of this section, previous studies by Grötzbach [31], Verzicco and Camussi [33], and Stevens et al. [21] found an increasing number of nodes that should be placed in the thermal and kinetic boundary layers. The theoretical results thus confirm all above studies, because the increasing number of nodes was due to the increasing Ra number at which the tests were performed. To be more specific: according to the estimates (44) and (45) for Pr = 0.7 the minimum number of nodes that should be placed in the thermal and kinetic boundary layers is N ≈ 2.3 for Ra = 3 × 10 5 , N ≈ 4.4 for Ra = 2 × 10 7 , and N ≈ 8.7 for Ra = 2 × 10 9 . The empirically found values at the respective Ra with Pr ≈ 0.7 are 3 for Ra = 3 × 10 5 , 5 for Ra = 2 × 10 7 , and 10 for Ra = 2 × 10 9 . Thus there is very good agreement between the theoretical results and the empirically obtained values, especially if one considers the difficulties involved in determining these values empirically, and the empirical value for the constant a (15) that is used in the theoretical estimates. We want to emphasize that not only the boundary layers close to the plates, but also the kinetic boundary layers close to the vertical walls must be well resolved.
To sum up, the mesh resolution should be analysed a priori using the resolution requirements in the bulk (36), (37) and in the boundary layers (42), (43). Having conducted the DNS, the Kolmogorov and Batchelor scale should be checked a posteriori, to make sure that the mesh size was indeed small enough (as it has been done, for example, in refs. [19,20]).

Conclusion
In summary, we used laminar Prandtl-Blasius boundary layer theory to determine the relative thicknesses of the thermal and kinetic boundary layers as functions of Pr (27).
We found that neither the position of the maximum r.m.s. velocity fluctuations nor the position of the horizontal velocity maximum reflect the slope velocity boundary layer thickness, although many studies use these as criteria to determine the boundary layer thickness. In contrast to them, the algorithm by Stevens et al. [21] agrees very well with the theoretical estimate of the kinetic slope boundary layer thickness.
We used the results obtained from the Prandtl-Blasius boundary layer theory to derive a lower bound on the minimum number of nodes that should be placed in the thermal and kinetic boundary layers close to the plates. We found that this minimum number of nodes increases not slower than ∼ Ra 0.15 with increasing Ra. This result is in excellent agreement with results from several numerical studies over the last decades, in which this minimum number of nodes was determined empirically. Hence, the derived estimates can be used as guideline for future direct numerical simulations.