Near wall Prandtl number effects on velocity gradient invariants and flow topologies in turbulent Rayleigh–Bénard convection

The statistical behaviours of the invariants of the velocity gradient tensor and flow topologies for Rayleigh–Bénard convection of Newtonian fluids in cubic enclosures have been analysed using Direct Numerical Simulations (DNS) for a range of different values of Rayleigh (i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Ra=10^7-10^9$$\end{document}Ra=107-109) and Prandtl (i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pr=1$$\end{document}Pr=1 and 320) numbers. The behaviours of second and third invariants of the velocity gradient tensor suggest that the bulk region of the flow at the core of the domain is vorticity-dominated whereas the regions in the vicinity of cold and hot walls, in particular in the boundary layers, are found to be strain rate-dominated and this behaviour has been found to be independent of the choice of Ra and Pr values within the range considered here. Accordingly, it has been found that the focal topologies S1 and S4 remain predominant in the bulk region of the flow and the volume fraction of nodal topologies increases in the vicinity of the active hot and cold walls for all cases considered here. However, remarkable differences in the behaviours of the joint probability density functions (PDFs) between second and third invariants of the velocity gradient tensor (i.e. Q and R) have been found in response to the variations of Pr. The classical teardrop shape of the joint PDF between Q and R has been observed away from active walls for all values of Pr, but this behavior changes close to the heated and cooled walls for high values of Pr (e.g. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pr=320$$\end{document}Pr=320) where the joint PDF exhibits a shape mirrored at the vertical Q-axis. It has been demonstrated that the junctions at the edges of convection cells are responsible for this behaviour for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pr=320$$\end{document}Pr=320, which also increases the probability of finding S3 topologies with large negative magnitudes of Q and R. By contrast, this behaviour is not observed in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pr=1$$\end{document}Pr=1 case and these differences between flow topology distributions in Rayleigh–Bénard convection in response to Pr suggest that the modelling strategy for turbulent natural convection of gaseous fluids may not be equally well suited for simulations of turbulent natural convection of liquids with high values of Pr.

www.nature.com/scientificreports/ The relative thicknesses of hydrodynamic and thermal boundary layers is dependent on Pr, which is known to affect the scalar spectrum and it is possible to obtain roll-off of the scalar spectrum in the inertial range and an inertial-diffusive range is obtained for Pr < 1 4 . By contrast, the roll-off of the scalar spectrum is obtained for the length scales smaller than the Kolmogorov length scale and a viscous-diffusive range is observed for the scalar spectrum for Pr ≫ 1 4 . As the temperature distribution in turbulent flows is affected by Pr, it can be expected that the velocity distribution in natural convection will also be affected by Pr because the flow is induced by the temperature difference.
The Prandtl number has indeed been shown to affect the turbulent kinetic energy spectrum in a recent analysis by the present authors 5 . However, the Prandtl number of the fluid does not only affect the distribution of turbulent kinetic energy in Rayleigh-Bénard convection but also has the potential to alter the distribution of flow topologies, as they can be categorised in terms of the invariants of the velocity gradient ∂u i /∂x j tensor (i.e. P, Q and R) where u i is the ith component of the velocity vector 6,7 . Depending on the values of the invariants of the velocity gradient ∂ u i /∂ x j tensor, 8 different topologies (i.e. S1-S8 topologies) can be identified in the threedimensional P, Q, R phase space. The velocity-gradient tensor can be split into symmetric and skew-symmetric parts: A ij = ∂u i /∂x j = S ij + W ij , where S ij = 0.5(A ij + A ji ) and W ij = 0.5 A ij − A ji are the symmetric and skew-symmetric components, which are referred to as strain and rotation rates, respectively. Three eigenvalues, 1 , 2 and 3 , of A ij can be obtained from solutions of the characteristic equation 3 + P 2 + Q + R = 0 where P, Q, R are the invariants of A ij 6,7 : The discriminant, D = 27R 2 + 4P 3 − 18PQ R + 4Q 3 − P 2 Q 2 /108 , of the characteristic equation divides the P − Q − R phase-space into two regions depending on the sign of the discriminant. For D > 0 (D < 0) , a focal (nodal) topology is obtained 6,7 and the velocity gradient tensor exhibits one real eigenvalue and two complex conjugate eigenvalues for focal topologies, whereas three real eigenvalues are obtained for nodal topologies. The solutions of D = 0 are given by two surfaces in the P − Q − R phase space 6,7 : r 1a = P Q − 2P 2 /9 /3 − 2 −3Q + P 2 3/2 /27 and r 1b = P Q − 2P 2 /9 /3 + 2 −3Q + P 2 3/2 /27 . For a positive discriminant (i.e. D > 0 ), the A ij tensor has purely imaginary eigenvalues on the surface r 2 , which is given by R = PQ . The surfaces r 1a , r 1b and r 2 , divide the P − Q − R phase space into eight flow topologies. The first invariant P = −∂u i /∂x i of the velocity gradient tensor vanishes for incompressible fluids, and therefore only topologies S1-S4 are observed for P = 0 , as shown in Fig. 1. Therefore, in Rayleigh-Bénard convection of incompressible fluids the flow topologies are determined by the behaviours of the second and third invariants (i.e. Q and R) of the velocity gradient tensor 6,7 and only S1-S4 topologies can be seen.
The flow structures associated with S1-S4 topologies are schematically shown in Fig. 1c. One aspect of this work focuses on flow topologies close to the active walls and therefore only the upper part (above the horizontal plane crossing the origin of the coordinate system) of the velocity field is shown in Fig. 1, representative of a situation close to the lower wall, but in principle it can be mirrored at the horizontal plane (see e.g. Ref. 8 ). It has been demonstrated by Perry and Chong 6 and Soria et al. 9 that S4 topologies are obtained predominantly for positive values of Q, whereas Blackburn et al. 10 demonstrated that the topologies S2 and S4 are predominantly obtained in the regions away from the wall in boundary layer flows. The 'teardrop' structure in the joint probability density function (PDF) between Q and R has been demonstrated by Chong et al. 7 and Chacin and Cantwell 11 . The analysis by Ooi et al. 12 and experimental evidences 11,13 suggested that the same qualitative behaviour is observed in a range of different incompressible turbulent flows indicating some degree of universality in the joint PDFs between Q and R. The theoretical justifications of the 'teardrop' shape of the Q − R joint PDF for incompressible flows have been provided by Elsinga and Marusic 14 and the loss of 'teardrop' structure was shown to be a mark of intermittency in some previous analyses 15 . Tsinober 16 postulated that the enstrophy production is large in S4 topology whereas the strain rate production is concentrated in regions of S1 topology. The flow topology distributions in Rayleigh-Bénard convection, where temperature and velocity fields are intrinsically coupled, are yet to be analysed in detail [17][18][19][20] in comparison to the vast body of literature (e.g. Refs. [5][6][7][8][9][10][11][12][13][14][15][16] ) on other wall-bounded flows.
The analyses by Dabbagh et al. [17][18][19] revealed the existence of the teardrop shape in the bulk region away from the walls in Rayleigh-Bénard convection but the small scale structures in the vicinity of the hot and cold walls have not been discussed there in terms of Q and R. Xi et al. 20 reported a transition of flow topologies from a quadruple structure to a dipole structure based on Rayleigh number in turbulent Rayleigh-Bénard convection, which has implications on the Nusselt number (or heat transfer rate). A recent analysis revealed that large-scale circulation in Rayleigh-Bénard convection is affected by Prandtl number 21 . However, the effects of Prandtl number on the flow topology are yet to be analysed and the present work addresses this gap in the existing literature. In this respect, the main objectives of the present analysis are: (a) to demonstrate and explain the effects of Prandtl number on the statistical behaviours of Q and R and their joint PDFs and (b) to indicate the implications of the above findings on flow topology distribution for Rayleigh-Bénard convection of Newtonian fluids. According to Buckingham's pi theorem 22 , the Nusselt number Nu = h L/k , where h represents the convective heat transfer coefficient, can be taken to be a function of Ra and Pr (i.e. Nu = f (Ra, Pr) ) for Rayleigh-Bénard convection in a cubic enclosure. Therefore, three-dimensional Direct Numerical Simulations (DNS) of Rayleigh-Bénard convection in a cubic enclosure for different values of Ra (i.e. Ra = 10 7 − 10 9 ) and Prandtl number ( Pr = 1 and 320) have been carried out in order to meet these objectives. It is worth noting that Pr = 320 corresponds for example to silicone oil at 20 o C which exhibits Newtonian rheological behaviour 23 .
The conservation equations for mass, momentum and energy for incompressible Newtonian fluids under transient conditions take the following form: (1) Scientific RepoRtS | (2020) 10:14887 | https://doi.org/10.1038/s41598-020-71665-9 www.nature.com/scientificreports/ The last term on the right-hand side of Eq. (3) originates due to Boussinesq's approximation and the temperature difference between horizontal walls is considered to be small enough so that this approximation remains valid. Also, in Eq. (3), ν is the kinematic viscosity and the Kronecker delta ( δ i2 ) indicates that buoyancy forces affect the flow only in the vertical direction (i.e. x 2 direction). The reference temperature ( T ref ) is taken to be the cold wall temperature (i.e. T ref = T C ). In Eq. (4), α=k/(ρ cp) is the thermal diffusivity of the fluid. Equations (2-4) are solved in a coupled manner in conjunction with the following boundary conditions. The simulation configuration is schematically shown in Fig. 1a which demonstrates that the differentially heated horizontal walls are subjected to constant wall temperature boundary conditions (i.e. T = T H at x 2 = 0 and T = T C at x 2 = L where T H > T C ). All the other walls are considered to be adiabatic (i.e. ∂T/∂x 1,3 = 0 at x 1,3 = 0, L ). Finally, no-slip and impermeability conditions are specified for all walls (i.e. u 1,2,3 = 0 at x 1,2,3 = 0, L).  Figure 2 indicates that distributions of temperature are affected by Pr, which is also reflected in a moderate influence of Pr on the mean Nusselt number (i.e. Nu ), which can be substantiated from Table 1.  where subscripts 'wf ' refers to the condition of the fluid in contact with the wall, T wall is the wall temperature and T ref is the appropriate reference temperature, which can be taken to be the temperature of the other wall. In Eq. (5), the mean value of the Nusselt number on the walls (i.e. Nu ) is obtained by averaging convective heat transfer coefficients ( h (A,t) ) over horizontal active walls and over time here.
The distributions of Q and R normalised by their standard deviations Q + = Q/σ (Q) and R + = R/σ (R) in the central x 1 − x 2 plane for Ra = 10 7 and 10 9 are shown for both Pr = 1 and 320 in Fig. 3. For incompressible fluids Q is given by: 0.5 −S ij S ij + W ij W ij and thus a positive (negative) value of Q is indicative of the vorticitydominated (strain rate-dominated) regions. It can be seen from Fig. 3 (top row) that the vorticity-dominated regions (i.e. Q + > 0 ) are predominantly obtained in the central core of the enclosure away from the wall. By contrast, the strain rate dominated (i.e. Q + < 0 ) regions are predominantly concentrated in the vicinity of the wall within the boundary layer.
The differences in the Q + distributions between the near-wall region and away from the wall regions can be explained by the large scale or plume dominated mean flow as it is shown in Fig. 2. The flow in the bulk region is driven by the buoyant force, whereas large velocity gradients can be found predominantly in the boundary layer region (i.e. strain rate-dominated areas). Therefore, this difference between bulk and boundary layer region changes depending on the dominant physical mechanisms. The dissipation rate of kinetic energy E = 2νS ij S ij can be expressed as: E = ν(−4Q S ) with Q S = −S ij S ij /2 which suggests that Q S = 0.25(−E)/ν assumes large negative values in the regions where the dissipation rate of kinetic energy remains large within the boundary layer 25 . The dissipation rate of kinetic energy weakens with increasing wall normal distance, and thus the magnitude of negative Q S values decreases with increasing wall normal distance, which in turn increases the propensity of the positive semi-definite values of Q W = W ij W ij /2 (i.e. Q W ≥ 0 ) to overcome negative values of Q S to yield positive values of Q = Q S + Q W away from the wall (i.e. in the bulk flow region in the middle of the enclosure).
For P = 0 , the third invariant R takes the form: where ω i is the i th component of vorticity. It is important to note that ( −S ij S jk S ki ) contributes to dissipation rate E = (2µ /ρ)S ij S ij = τ ij (∂u i /∂x j )/ρ generation (with τ ij being the component of the viscous stress tensor), whereas ( ω i ω j S ij ) contributes to the production rate of enstrophy (i.e. � = ω i ω i /2 ). Therefore, the sign of the R indicates the competition and relative strengths of the enstrophy production rate and the dissipation rate generation 16 . A comparison between Q + and R + fields in Fig. 3 reveals that in particular large negative values of Q + are mostly associated with large positive values of R + in the bulk region of the enclosure for both Pr = 1 and 320 and this is particularly prominent for Pr = 1 . However, the near-wall behaviour is different for different Prandtl numbers. In order to demonstrate this behaviour, the contours of joint probability density functions of Q + and R + in the bulk region (defined as  Fig. 4 in the case of both Pr = 1 and 320. Although Fig. 4 shows the expected and well-known teardrop shape suggesting predominance of S4 and S2 (and to some extent S1) topologies both in the bulk region away from the walls and at the hot and cold boundaries for Pr = 1 , the conventional teardrop shape is obtained only in the bulk region for Pr = 320 and the lower tail of the joint PDF flips from an unstable node-saddle-saddle S2 topology towards www.nature.com/scientificreports/ a stable node-saddle-saddle topology S3 at the hot and cold walls. These differences in joint PDFs between Q and R are expected to have implications on the distribution of the flow topologies within the enclosure. The variations of the volume fraction VF of flow topologies averaged in x 1 − x 3 planes in the vertical direction (i.e. x 2 -direction) are shown in Fig. 5 for Ra = 10 7 and 10 9 in the case of both Pr = 1 and 320. It can be seen from Fig. 5 that Pr = 320 (and to a lesser extent for Pr = 1 ), the volume fraction of S1 and S4 topologies is about 40% in the bulk region, whereas the volume fraction of S2 and S3 is about 10% each within the bulk region of the domain regardless of Ra values. However, it is worth mentioning that the sum of both unstable and both stable topologies always seems to be very close to 50% for all values of Pr and Ra. This behaviour is consistent with the theoretical estimates by Hasslberger et al. 26 assuming a symmetric population in Q − R space. Close to the boundaries the volume fractions of nodal topologies (S2, S3) increase, whereas the volume fractions of focal topologies (S1,S4) decrease.
Although results are shown only for Ra = 10 8 in Fig. 4, the distributions of Q − R joint PDFs are qualitatively similar for the range of Rayleigh numbers considered in this study ( 10 7 ≤ Ra ≤ 10 9 ). This can be substantiated from the distributions of flow topologies in the central x 1 − x 2 plane for Ra = 10 7 and 10 9 at Pr = 1 and 320, which are shown in Fig. 6. It can be seen from Fig. 6 that the focal topologies S1 and S4 are predominantly obtained in the bulk region at the interior of the domain, whereas the nodal topologies S2 and S3 are dominant in the vicinity of the hot and cold walls. This behaviour does not change with the variation in Rayleigh number within the Ra = 10 7 − 10 9 range considered here. However, the topology distribution in Fig. 6 suggests that the small-scale structures become more frequent for larger values of Rayleigh number and smaller values of Prandtl number which implies an increasing Grashof number.
In order to analyse the origin of the flipping of the tail of the joint PDF between Q and R towards the S3 quadrant, the region close to the bottom heated wall at x 2 = 0 is investigated further. The iso-surfaces of nondimensional temperature θ = 0.65 (for better visibility of the structures θ = 0.65 is used in Fig. 7a) coloured by non-dimensional vertical velocity V for Ra = 10 7 , 10 8 and Pr = 1, 320 are shown in Fig. 7 (first row) together with the corresponding distributions of Q + and R + and flow topologies on the θ = 0.85 isosurface (second to fourth row). The ridge like structures for Pr = 320 correspond to plume regions with large wall normal velocities directed away from the wall. The peaks where two ridges meet each other will subsequently be called junction points and a comparison with the flow topologies reveals that the junctions can be associated with topology S3. However, it can be seen that the S3 topology can also be observed in the valleys in between the ridges. In order to identify the origin of the reverse tail of the joint PDF between Q + and R + it is instructive to identify the regions with large negative values of Q + and R + and these quantities are mapped onto the iso-surfaces of non-dimensional  www.nature.com/scientificreports/ temperature θ = 0.85 as well in Fig. 7. It becomes obvious from Fig. 7 that for Pr = 320 there is a one-to-one relation between the junctions and the locations responsible for exhibiting S3 topology which is observed for the lower-left tail of the joint PDF between Q + and R + (i.e. for large negative Q + and R + ). Figure 2 reveals that isolated plumes drive the convection process in the Pr = 320 case, whereas frequent roll ups in the Pr = 1 case are indicative of a large-scale circulation. This is further illustrated in Fig. 8 showing the variation of non-dimensional temperature iso-surfaces θ = 0.4 and θ = 0.6 together with path lines coloured by non-dimensional vertical velocity magnitude for different Pr, exemplarily at Ra = 10 7 . This behaviour is consistent with previous findings by van der Poel et al. 21 and Verzicco and Camussi 27 . Figure 7 further reveals that it is rare to obtain simultaneous occurrences of large negative values of Q + and R + in the near wall region for reveals that the plumes in the Pr = 320 case are representative of the S3 topology and thus the occurrence of S3 with large negative values of Q + and R + are more likely in this value of Pr, which is also reflected in the flipping of the tail of the joint PDF between Q + and R + to the quadrant of S3 topology. As the convection is driven by large-scale circulation in the Pr = 1 case, the occurrences of S3 topology with large negative values of Q + and R + are either rare or absent.

Summary
The statistical behaviours of the invariants of the velocity gradient tensor and flow topologies in turbulent Rayleigh-Bénard convection of Newtonian fluids in cubic enclosures have been investigated using three-dimensional Direct Numerical Simulations (DNS) for a range of different Rayleigh (i.e. Ra = 10 7 − 10 9 ) and Prandtl numbers (i.e. Pr = 1 and 320). It has been found that the convection in the case of large values of Pr is plume-dominated, whereas a large-scale circulation in the enclosure has been found for Pr = 1 . The focal topologies S1 and S4 have been found to be the two dominant topologies in the bulk region, whereas the probabilities of obtaining nodal topologies S2 and S3 increase in the vicinity of the active hot and cold walls. The proportion of different flow topologies in the bulk region has been found to be consistent with a previous analytical study 27 assuming a symmetric population in the two-dimensional phase space based on second and third invariants of the velocity gradient tensor (i.e. Q and R). The classical teardrop shape of the joint PDF between Q and R has been observed away from active walls, but this behavior changes close to the heated and cooled walls where the joint PDF depicts a shape mirrored at the vertical Q-axis for large values of Prandtl number (e.g. Pr = 320 ). It has been demonstrated that junctions, at the edges of convection cells arising from plume-dominated convection, are responsible for this behavior for Pr = 320 . This behaviour is not present in the Pr = 1 case due to the convection driven by a large-scale circulation. The differences between flow topology distributions in Rayleigh-Bénard convection in response to Pr suggest that turbulence modelling of natural convection possibly should explicitly account for Prandtl number effects. Therefore, the modelling methodologies for simulations of turbulent natural convection of gaseous fluids may not be equally well-suited for natural convection of liquids with large values of Prandtl number. www.nature.com/scientificreports/

Methods
The governing equations of mass, momentum and energy conservation equations have been solved in a finitevolume framework using an open-source CFD package OpenFOAM. For these computations pressure-velocity coupling has been addressed by the use of the PIMPLE algorithm. The convective and diffusive fluxes are evaluated using second-order central difference schemes. The temporal advancement has been carried out using the second-order Crank-Nicolson scheme in conjunction with adaptive time-stepping for the sake of computational economy. It has been ensured that the Courant number is always below unity so that the simulations have enough temporal resolution. Ideally, the grid size x for DNS of Rayleigh-Bénard convection should satisfy �x < min( η K , η B ) where η K is the Kolmogorov length scale and η B is the Batchelor length scale. These length scales are defined as: where �ǫ� = �2ν S ′ ij S ′ ij � is the mean kinetic energy dissipation rate with the angled bracket indicating a volume-time averaging technique. Here, S ′ ij is the strain rate component based on the fluctuating velocity field and ν is the kinematic viscosity. The mean Kolmogorov scale η K is the smallest length scale of turbulence when ν ≤ α (i.e. Pr ≤ 1 ). In case of ν > α (i.e. Pr > 1 ), the smallest length scale is determined by the dissipation rate �ǫ T � = �α (∇T ′ ) 2 � of scalar variance, and is represented by the Batchelor scale η B 4 . Criteria proposed by Grötzbach 28 and Shishkina et al. 29 have been used for determining the initial non-dimensional mesh size ( �x/L ), along with the corresponding Nu correlations suggested by Refs. 23,30 .
Additionally, Grötzbach 28 suggested that at least three grid points should be placed in the viscous boundary layers to estimate the realistic Nusselt number. Verzicco and Camussi 31 and Stevens et al. 32 found that the number of nodes that should be placed in the thermal and hydrodynamic boundary layers increases with increasing Gr. Based on this, the Pr = 1 cases are more critical in terms of grid resolution than the Pr = 320 cases, because the values of Gr are much higher for Pr = 1 cases than the one in the Pr = 320 cases for the same set of Ra (i.e. Gr = Ra/Pr ). Lam et al. 33 proposed a correlation for the normalised viscous layer thickness as δ/L = 0.65 Gr −0.16 Pr 0.08 = 0.65Ra −0.16 Pr 0.24 based on experimental analysis of turbulent Rayleigh-Bènard convection of Newtonian fluids in a single convection cell of unity aspect ratio for 10 6 ≤ Ra ≤ 10 11 and 6 ≤ Pr ≤ 10 3 . Using this normalised viscous layer thickness correlation 33 more than 15 grid points have been placed in the Pr = 1 cases for all the values of Ra ( 10 7 − 10 9 ) considered in the current analyses.
Finally, it has been ensured a-posteriori that the chosen grid resolution is sufficient by ensuring y + and y + √ Pr (where y + = ρu τ y/µ with u τ = √ τ w /ρ , τ w and y being the friction velocity, wall shear stress magnitude and wall normal distance of the wall adjacent grid point, respectively) remains smaller than unity. For the present analysis, Cartesian grids of 250 3 , 460 3 and 700 3 ( 150 3 , 230 3 and 490 3 ) have been used for Pr = 1 ( Pr = 320 ) simulations of Ra = 10 7 , Ra = 10 8 and Ra = 10 9 , respectively.