Horizontal convection: Effect of aspect ratio on Rayleigh number scaling and stability

Applied Mathematical Modelling 35 (2011) 1647–1655

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Horizontal convection: Effect of aspect ratio on Rayleigh number scaling and stability Gregory J. Sheard a,⇑, Martin P. King b a b

Fluids Laboratory for Aeronautical and Industrial Research (FLAIR), Department of Mechanical and Aerospace Engineering, Monash University, VIC 3800, Australia School of Engineering, Sunway Campus, Monash University, Jalan Lagoon Selatan, Bandar Sunway, 46150 Selangor Darul Ehsan, Malaysia

a r t i c l e

i n f o

Article history: Received 27 August 2010 Accepted 17 September 2010 Available online 29 September 2010 Keywords: Horizontal convection Differentially heated box Nusselt number scaling Overturning circulation Boussinesq Spectral-element method

a b s t r a c t Horizontal convection in a rectangular enclosure driven by a linear temperature profile along the bottom boundary is investigated numerically using a spectral-element discretization for velocity and temperature fields. A Boussinesq approximation is employed to model buoyancy. The emphasis of this study is on the scaling of mean Nusselt number and boundary layer quantities with aspect ratio and Rayleigh number. At low Rayleigh number, Nusselt number and boundary layer thickness are found to be independent of Rayleigh number but dependent on aspect ratio. At higher Rayleigh numbers, convective flow dominates, and Nusselt number, boundary layer thickness and peak boundary layer velocity become independent of aspect ratio. In this regime, the Rayleigh number scaling of these quantities agrees well with exponents predicted by theory, with respective values of 1/5, 1/5 and 2/5. Unsteady flow develops at a critical Rayleigh number independent of aspect ratio, and the development of unsteady flow is found to lead to an increase in the Nusselt number scaling exponent from 0.2 to approximately 0.3, which is closer to the theoretical upper bound than has yet been reported in the study of horizontal convection flows. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction Horizontal convection refers to the heat and fluid flows established in an enclosure due to differential heating along just one horizontal boundary [1]. This is in contrast with other forms of convection which are often driven by a temperature differential imposed between two opposite boundaries (see, for instance [2]). Whether horizontal convection is achieved by an applied horizontal temperature gradient or by an applied heat flux, unstable convective flow is forced in one side of the enclosure while the rest of fluid is convectively stable. Therefore, in contrast to the extensively studied Rayleigh–Bénard convection, whereby the cooling and heating both promote convective overturning, the strength of overturning in horizontal convection is ultimately limited by heat diffusion. Motivation for the study of horizontal convection comes from geophysical and geological flows. For example, despite heavy simplifications, studies into horizontal convection are providing understanding and insight into meridional (North–South) overturning circulation in the oceans, where they are heated and cooled along a thin horizontal layer. Interest in horizontal convection is also emerging amongst researchers including engineers, applied mathematicians and oceanographers. For reviews which discuss recent advances and outstanding questions in this subject, the readers are referred to Hughes and Griffiths [1], Wunsch and Ferrari [3], as well as references therein.

⇑ Corresponding author. E-mail addresses: [email protected] (G.J. Sheard), [email protected] (M.P. King). 0307-904X/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2010.09.041

1648

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

A brief description of horizontal convection follows. For further details readers are directed to Hughes and Griffiths [1] describe horizontal convection in detail, but a brief overview is included here. At low Rayleigh numbers, flow driven by horizontal convection is dominated by diffusion, and is stable in time. It comprises a nearly symmetrical overturning circulation of fluid driven by a buoyancy destabilization on the heated boundary with flow moving along the bottom boundary from the cold end to the hot end (or the hot to cold end if the top boundary is the heated boundary). Buoyant fluid then rises (or descends) in a narrow vertical plume adjacent to a side wall, before returning to complete the circulation in a diffusive horizontal return flow. As the Rayleigh number is increased, studies have shown that convective effects begin to dictate the fluid and heat transfer behaviour, with thermal and velocity boundary layers developing along the heated boundary. Beyond some critical Rayleigh number, the flow eventually becomes unsteady, which is particularly visible in the vicinity of the vertical plume [4]. A scaling analysis by Paparella and Young [5], in which dissipation was shown to vanish as kinematic viscosity and thermal diffusion go to zero, was used to present an argument that horizontal convection was inherently non-turbulent. Regardless of whether the flow is formally classified as turbulent, it does feature the convection of irregular small-scale buoyant flow structures from the heated wall boundary layer into the vertical plume beyond some critical Rayleigh number. Experiments by Mullarney et al. [4] with water in an enclosure with height-to-width aspect ratio of 0.16 showed that beyond the diffusion-dominated regime, the Nusselt number scaled with approximately the 1/5th power of Rayleigh number. However, Siggers et al. [6] used a variational analysis to determine that an upper bound on Nusselt number scaling was a 1/3rd power of Rayleigh number. To the authors’ knowledge, scaling exponents above 1/5th have not previously been reported for horizontal convection flows. Chiu-Webster et al. [7] studied horizontal convection in the infinite-Prandtl number limit relevant to very viscous fluids, at a range of aspect ratios and Rayleigh numbers. That study also found the Nusselt number to scale with the 1/5th power of Rayleigh number, and presented evidence of an aspect ratio independence beyond Rayleigh numbers of approximately 107. Despite these past investigations, the aspect ratio-dependence of features such as the transition to the convection-dominated regime, and the onset of unsteady flow remain poorly understood. In this paper we investigate the effect of the aspect ratio of the enclosure on the scaling relationships for heat transfer, boundary layer thicknesses, boundary layer velocities, overturning, and the transition to unsteady flow. 2. Model description 2.1. Problem definition The problem considered in this paper is the two-dimensional horizontal convection of fluid in a rectangular enclosure of width L and height D. The flow is driven by a linear temperature profile applied along the bottom wall of the enclosure, as illustrated in Fig. 1. The side and top walls are insulated (a zero temperature gradient is imposed normal to the walls), and a no-slip condition is imposed on the velocity field on all walls. A Boussinesq approximation of the fluid buoyancy is employed, whereby density differences in the fluid are neglected with the exception of the gravity contribution. A scalar field representing the fluid temperature (which relates linearly to the density via a thermal expansion coefficient, a) is evolved via an advection–diffusion operation in conjunction with the velocity field. 2.2. Governing equations and parameters The equations governing a Boussinesq fluid may be written as

@u ^ T; ¼ ðu  rÞu  rp þ Prr2 u  Pr Ra g @t r  u ¼ 0; @T ¼ ðu  rÞT þ r2 T; @t

ð1Þ ð2Þ ð3Þ

Fig. 1. A schematic diagram of the system. The origin of the coordinate system is positioned at the bottom-left corner, gravity acts vertically downward, and a temperature difference of dT is imposed along the bottom wall.

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

1649

^a where u is the velocity vector, p the kinematic static pressure, t is time, Ra is the Rayleigh number, Pr the Prandtl number, g unit vector in the direction of gravity, and T is the temperature. In Eqs. (1)–(3), lengths have been scaled by the enclosure width L, velocities by L/jT (where jT is the thermal diffusivity of the fluid), time by jT/L2, and temperature by dT (the imposed temperature difference imposed across the bottom wall). The (horizontal) Rayleigh number is defined as

g adTL3

Ra ¼

mjT

;

where g is the acceleration due to gravity and m is the kinematic viscosity of the fluid. If a heat flux is defined as

F T ¼ jT q0 cp

@T ; @y0

where q0 is a reference fluid density, cp the specific heat capacity of the fluid, @T=@y is the mean temperature gradient along the bottom wall over 0 6 x 6 L/2, then a flux Rayleigh number may be defined as

RaF ¼

g aF T L4 : q0 cp j2T m

In the experimental study by Mullarney et al. [4], horizontal convection was driven by applying a heat flux over half of the bottom wall, and a constant temperature along the other half. Those flows were conveniently scaled using the flux Rayleigh number, whereas the horizontal Rayleigh number (Ra) is more appropriate to use for horizontal convection driven using a linear temperature profile as applied in the present study. The relationship between fluid viscosity and thermal diffusivity is parameterized by the Prandtl number

Pr ¼

m : jT

Throughout this study the Prandtl number is maintained at Pr = 6.14, consistent with water at room temperature. Finally, the mean Nusselt number represents the ratio of convective to conductive heat transfer, and may be defined as

Nu ¼

FT L

q0 cp jT dT

:

2.3. Numerical approach The Boussinesq flow described by Eqs. (1)–(3) is computed on a two-dimensional domain using a high-order in-house solver, which employs a spectral-element method for spatial discretization and a third-order time integration scheme based on backwards-differencing. The fluid component of this solver was employed in recent studies into wake dynamics around arresting bodies [8], and as the basis for stability analysis computations of bluff-body wake flows in [9,10]. To model a Boussinesq flow, two modifications to the standard fluid flow solver are required: Firstly, the transport of a scalar field representing temperature is computed concurrently with the velocity field (Eq. (3)); secondly, the momentum equation (Eq. (1)) is modified by the addition of a gravity term describing the contribution of buoyancy. Meshes were constructed for enclosures with aspect ratios D/L = 2, 1, 0.625, 0.333 and 0.16. Care was taken to ensure that the flow was resolved in the vicinity of each of the walls, and in particular on the heated boundary, with coarser meshing employed in the interior. A total of between 412 and 1229 spectral-elements were used in the meshes, and for all simulations elements with a polynomial degree of 8 were used. In the finite-volume numerical simulations of [4], 12-element-wide boundary layer meshes with a thickness of 0.0372L were used. By comparison, in the meshes used in this study, the nearest 12 collocation points to the heated boundary were distributed over distances from the boundary ranging from 0.037L down to 0.0030L. The present study spans a very wide range of Rayleigh number, which for higher Rayleigh number placed considerable limitations on permissible time steps. Beyond Ra  O(105), the maximum allowable time step in the computations scaled approximately with 1/Ra. Following the boundary- and thermal-layer scaling analysis [4], higher spatial resolution was also required as fluid scales reduced with increasing Ra. 3. Results With an increase in Rayleigh number, the horizontal convection flow passes from a diffusion-dominated regime to a steady-state convection-dominated regime, before subsequently developing unsteady flow, which is concentrated in the vicinity of the vertical plume rising from the hot end of the heated boundary. Fig. 2 plots temperature contours and streamlines for horizontal convection in an enclosure with D/L = 0.625 in each of these regimes.

1650

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

Fig. 2. Contour plots of temperature overlaid with velocity streamlines for horizontal convection in an enclosure with D/L = 0.625 at Rayleigh numbers Ra = (a) 1.43  103, (b) 1.43  108 and (c) 1.43  109 Light and dark contours show arbitrary levels of warm and cool fluid, respectively.

The same progression through these regimes was found for all enclosure aspect ratios in the range investigated, 0.16 6 D/ L 6 2, though the Rayleigh numbers marking the transition between neighboring regimes do exhibit a dependence on D/L, which will be explored subsequently. Flows are computed in this study over a wide range of Rayleigh numbers 4.36  104 6 Ra 6 8.52  1011. The quantities of most interest in horizontal convection are the Nusselt number, the thermal and velocity boundary layer thicknesses, and the peak velocity within the boundary layer adjacent to the heated boundary. Scalings for these quantities proposed by [4] with flux Rayleigh number (RaF) are given as

Nu / Ra1=6 F ; UL=jT / RaF1=3 ; h=L / RaF1=6 : In the present study, convection is controlled not by the flux Rayleigh number, but by the horizontal Rayleigh number Ra, and thus the expected scaling relationships are recast by RaF = Nu Ra as

Nu / Ra1=5 ; UL=jT / Ra2=5 ; h=L / Ra1=5 : The mean Nusselt number is calculated in the present configuration by computing the average temperature gradient normal to the wall along the cooler half of the heated boundary, as defined in Section 2.2. The calculated Nusselt numbers are plotted against Rayleigh number in Fig. 3(a). The plot shows that at low Rayleigh number, the Nusselt number is independent

(a)

(b)

1

2

d(log10 Nu)/d(log10 Ra)

0.8

log10 Nu

1.5 1 1 5

0.5 0

0.6

0.4

0.2

0

-0.5 2

4

6

log10 Ra

8

10

12

2

4

6

8

10

12

log10 Ra

Fig. 3. (a) A plot of log10Nu against log10Ra, for D/L = 2 (h), 1 (M), 0.625 (}), 0.333 (r), and 0.16 (s). Akima splines are fitted to the data for guidance. A dotted line shows the empirical trend proposed by Mullarney et al. [4]. (b) A plot of the gradient of the curves in (a), calculated using finite differences. A dotted line illustrates the theoretical gradient of 1/5 from [4], and a dash-dot line shows the theoretical upper bound proposed by Siggers et al. [6].

1651

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

of Rayleigh number, but is a function of D/L, with Nu decreasing with decreasing D/L. As Ra is increased, the data for each aspect ratio collapse onto a single trend, which is nearly linear on a log–log plot, with a gradient of 1/5. This gradient is consistent with theoretical scaling [4]. The trend in Fig. 3(a) exhibits a slight increase in gradient beyond approximately log10Ra  9.5. It is noted that this corresponds to the development of unsteady flow in the enclosure. To further elucidate the scaling of the data shown in Fig. 3(a), gradients were computed by fitting Akima splines to the data (Akima splines are less susceptible to the wiggle artifacts which affect other curve-fitting functions such as polynomial interpolation or cubic splines [11]), and calculating the gradients using finite differences. This data is plotted in Fig. 3(b). In enclosures with D/L J 1, the gradient increases from zero to 1/5 at log10Ra  4. With decreasing D/L, the collapse to the gradient of 1/5 occurs at higher Ra; i.e. when D/L = 0.16, log10 Ra  8. Fig. 3(b) confirms that beyond log10Ra  9.5, the calculated gradients increase from approximately 0.2 to values in the range 0.25 to 0.30. This is a significant observation, as such a gradient was not detected in the measurements of Mullarney et al. [4], yet Siggers et al. [6] performed an analysis which proposed an exponent of 1/3 as the upper bound on Nu–Ra scaling for horizontal convection. The elevated gradients detected in these simulations may signify a transition to a previously undetected regime of horizontal convection. The simulations performed in this study employ a higher spatial resolution than that employed in the numerical simulations conducted by Mullarney et al. [4], which may explain why that study did not report scaling exponents beyond 1/5. The increase in gradient detected here occurs as unsteady flow develops in the enclosure shows that heat transfer is enhanced by the development of unsteady flow in this configuration. When interpreting these results it should be noted that they are obtained from a two-dimensional model, and therefore do not resolve the threedimensional effects which could emerge at higher Rayleigh numbers. Thermal and velocity boundary layer thicknesses are calculated at x = L/2. These quantities are plotted against Rayleigh number in Fig. 4. The thermal boundary layer thickness is taken to be the point at which the temperature is 5% less than the temperature at the top wall, and the velocity boundary layer thickness is taken to be at the point of maximum velocity in the boundary layer. These definitions are consistent with Mullarney et al. [4] with the exception that in that study the thermal boundary layer thickness was taken with reference to the temperature at mid-height instead of the top wall. In the low Rayleigh number diffusion regime, the boundary layer thicknesses are independent of Rayleigh number. In this regime, enclosures with larger D/L have a larger boundary layer thickness. With increasing Ra, the hthermal/L data overshoots the empirical trend measured by Mullarney et al. [4], before collapsing onto a single trend with a gradient of approximately 1/5, consistent with theory. Similar Rayleigh number dependence and collapse behaviour is observed for hvelocity/L in Fig. 4(b). The peak boundary layer velocity on the centreline of the enclosure displays two regimes of linear behaviour with Rayleigh number on a log–log plot as shown in Fig. 5. For all D/L, a unit gradient is found in the low-Ra regime, which persists to approximately log10Ra  3.5 to 6.5 for D/L  2 down to 0.16. At higher Rayleigh numbers, the data again collapses to a linear trend, this time with a gradient of 2/5, which is consistent with theory. Furthermore, the data in this convective regime agrees very well with the empirical trend reported by Mullarney et al. [4]. To characterize the enclosure aspect ratio dependence on transition from the diffusion-dominated regime to the convective regime, a criterion was established whereby deviation of more than 5% from the Ra-independent values of Nu and hvelocity/L identified the critical Rayleigh number. Fig. 6 plots the critical Rayleigh numbers as a function of D/L for the data shown in Fig. 3(a) and Fig. 4(b). It is found that for Nu, the critical Rayleigh number denoting transition to convective flow is insen-

(a)

(b)

0

-0.5

log10 hvelocity/L

log10 hthermal/L

-1

-1

-1.5

1 -2

1 -1.5

5

5 0

2

4

6

log10 Ra

8

10

0

2

4

6

8

10

log10 Ra

Fig. 4. Plots of (a) log10hthermal/L and (b) against log10Ra for various D/L. Symbols are as per Fig. 3. Dotted lines shows the empirical trends proposed by Mullarney et al. [4], and a gradient of 1/5 is provided for comparison with theory.

1652

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

4

2

2

log10 umaxL/κT

5 0

-2 1 -4

1

-6 -4

-2

0

2

4

6

8

10

12

log10 Ra Fig. 5. A plot of log10umax/jT) against log10Ra for various D/L. Symbols are as per Fig. 3. The dotted line shows the empirical trend proposed in [4], and a gradient of 2/5 is provided for comparison with theory.

7

log10 Racrit

6

5

4

3

0

0.5

1

1.5

2

D/L Fig. 6. A plot of the logarithm of critical Rayleigh number against D/L, for Nu (h) and hvelocity/L (s) data. Curves are spline fits to the data for guidance.

sitive to enclosure aspect ratio above D/L  1. At lower aspect ratios, the critical Rayleigh number increases hyperbolically as D/L ? 0. This trend is shared by the boundary layer thickness, though at consistently higher critical Rayleigh numbers, with the exception that above D/L  1, the critical Rayleigh number is seen to decrease appreciably with increasing D/L. Due to the sparsity of the original data (hvelocity/L was computed at power-of-10 Ra intervals), it is unclear whether the Racrit data point at D/L = 0.333 undershoots an otherwise hyperbolic trend, or whether the low-D/L values of 0.16 and 0.333 form a separate branch to the data at D/L P 0.625. It is known that beyond some Rayleigh number in the convective regime, the horizontal convection flow develops unsteady flow. An example of this can be seen in Fig. 2(c), where the vertical plume in the bottom-right corner of the enclosure exhibits a time-dependent pulsing. Similar observations can be made from the visualizations in Mullarney et al. [4]. Analysis of time histories of heat flux through the bottom wall permitted the temporal characteristics of the saturated flows computed in this study to be determined. It was found that somewhere in the range 3.5  108 [ Ra [ 8.5  108, the flow transitioned from a steady to a time-dependent state. This appears to be consistent with the transition from the steady regime (regime I) to the entrainment regime (regime III) from the regime diagram in Hughes and Griffiths [1]. Fig. 7 plots the log10Ra-D/L parameter space computed in this study, identifying steady and time-dependent cases. It is found that there exists little or no aspect ratio dependence on the transition Rayleigh number for unsteady flow. Given that the Nu and boundary layer data at various D/L have been shown to collapse onto a single curve, implying independence on D/L, it is interesting to

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

1653

Unsteady flow

log10 Ra

10 Convectiondominated flow

5

Diffusiondominated flow

0

0

0.5

1

1.5

2

D/L Fig. 7. A plot of the log10Ra–D/L parameter space computed in this study. Symbols denote points at which data was acquired, with open and filled symbols showing steady and unsteady solutions, respectively. Solid and dash-dot lines reproduce the critical Rayleigh number curves given in Fig. 6 for comparison, and a dashed line marks the approximate boundary between steady and unsteady regimes.

observe that a similar independence on enclosure aspect ratio is found for the transition to unsteady flow. This result strongly suggests that the mechanism leading to the transition to unsteady flow is closely tied to the thermal and velocity boundary layers along the heated boundary. In a two-dimensional flow, the difference between stream function values on any two streamlines returns the volumetric flux between those streamlines. Thus in an enclosed flow, a measure of the strength of recirculating turnover of fluid within the enclosure is given by the maximum absolute value of the stream function within the enclosure (Wmax, where the stream function is defined as zero at the walls of the enclosure). Fig. 8 plots the logarithm of Wmax against the logarithm of Ra for each of the enclosure aspect ratios considered in this study. At low Rayleigh numbers, Wmax is dependent on D/L, and exhibits a scaling exponent of 1 (i.e. Wmax / Ra). With increasing Rayleigh number, the different D/L trends each collapse onto a single trend which exhibits a scaling exponent of 1/5 (i.e. Wmax / Ra1/5), which is apparently due to boundary layer transport [1]. At the same Rayleigh number that unsteady flow was first detected within the enclosures, the stream function displays a transition to a scaling behaviour with an exponent of approximately 2/5 beyond log10 Ra  8.3, approximately in unison. This implies that the transition is independent of D/L, which in turn suggests that the different scaling behaviour is caused by an alteration to the flow within the boundary layer on the heated boundary. It is proposed that the development of transient

3 5 2

2 5

log10 Ψmax

1

1

0

-1 1

-2 1 -3

2

4

6

8

10

12

log10 Ra Fig. 8. A plot of log10Wmax against log10Ra, for D/L = 2, 1, 0.625, 0.333, and 0.16, where symbols and guidance lines are as per Fig. 3. For guidance of scaling in the diffusion-dominated regime and the steady and unsteady convective regimes, gradients of 1, 1/5 and 2/5, respectively, are supplied.

1654

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

(a)

(b)

0 -0.2

-0.4

-0.4

-0.6

log10 Δy/L

log10 Δx/L

-0.6 -0.8 1

-1

5

-1.2

5 -0.8

1

-1

-1.4

-1.2

-1.6 -1.8

-0.2

2

4

6

log10 Ra

8

10

-1.4

2

4

6

8

10

log10 Ra

Fig. 9. Plots of the logarithms of (a) the horizontal and (b) the vertical distances of the point of Wmax to the side wall adjacent to the vertical plume and heated boundary, respectively, against the logarithm of Rayleigh number. Data from each enclosure is represented by symbols and guidance lines as per Fig. 3. For guidance, gradients of 1/5 are provided in (a) and (b), respectively.

upwelling plumes from the heated boundary is responsible for this change in behaviour. A recent study [12] proposed a model for horizontal convection featuring a vertical turbulent line plume. When recast in terms of the conventional Rayleigh number, that model predicted Wmax to scale with Ra3/10, which is somewhat below the scaling behaviour suggested in Fig. 8. One interpretation of what is observed here beyond log10Ra  8.3 is that the increased gradient could mark part of a migration to a subsequent scaling regime at higher Rayleigh numbers, and that the scaling exponents associated with that regime may not necessarily remain consistent with the initial gradients detected here. The location of the point of maximum streamfunction denotes the point about which the flow in the enclosure is circulating. For all enclosures studied here, this point was located at the centre of the enclosure at low Rayleigh numbers (the diffusion-dominated regime), and once the convective regime was reached, this point was found to migrate towards the corner where the buoyant fluid exited the heated boundary layer through the vertical plume. Presumably, the rate at which this point of maximum streamfunction approaches the bottom and side walls will be dictated by the scaling of the horizontal boundary layer thickness and the vertical plume width, respectively. Defining the horizontal and vertical distances of the point of maximum streamfunction to the horizontal and vertical boundaries extending from this corner as Dx and Dy, respectively, the logarithm of these quantities (normalized by L) are plotted in Fig. 9. The plots in Fig. 9 confirm that throughout the low Rayleigh number regime, the point of maximum streamfunction remains stationary, very near the centre of the enclosure. This is reflected in the streamlines shown in Fig. 2(a). In Fig. 9(a) the curves are coincident at low Rayleigh numbers as a result of the normalisation, as the position was consistently L/2 from the side wall for all enclosures. In Fig. 9(b) the initial vertical distances are seen to be dependent on D/L, as the initial position was consistently D/2 from the heated horizontal boundary. However, in both plots the curves are observed to collapse onto an approximately linear trend once the convective regime is reached. These linear trends both adopt gradients consistent with a Rayleigh number scaling of Ra1/5. This is consistent with the expected scaling for the boundary layer thickness on the heated boundary reported earlier and elsewhere [4]. The adoption by the plume width of the same scaling behaviour suggests that the plume scaling is governed by transport within the boundary layer on the heated wall. As noted earlier, a marked aspect ratio dependence is observed for the onset of the convective regime, with lower D/L enclosures transitioning to the steady convective regime at higher Rayleigh numbers. This is so pronounced that the enclosure with D/L = 0.16 almost completely bypasses the steady convective regime on the way to the unsteady convective regime. Finally, it is noted that beyond log10Ra  8 the data in Fig. 9 is disrupted by the development of unsteady flow at equilibrium, which acts to vary the location of the point of maximum streamfunction with time. This data was extracted from instantaneous snapshots at equilibrium rather than from time-averaged sequences, and thus should not be treated as a definitive representation of these quantities beyond log10Ra  8. 4. Conclusion Horizontal convection has been computed at high spatial resolution using a spectral-element method, over a wide range of Rayleigh numbers and enclosure aspect ratios, at a Prandtl number of 6.14, which is consistent with water at room temperature.

G.J. Sheard, M.P. King / Applied Mathematical Modelling 35 (2011) 1647–1655

1655

At low Rayleigh number, the mean Nusselt number and boundary layer thickness demonstrate Rayleigh number independence, though they do vary with aspect ratio. Above some critical Rayleigh number, mean Nusselt number, boundary layer thickness, and boundary layer velocity each collapse to single curves independent of aspect ratio, and in agreement with theory, i.e. Nu / Ra1/5, h/L / Ra1/5, and UL/jT / Ra2/5, respectively (see Figs. 3(a), 4 and 5). At higher Rayleigh numbers there is evidence for transition towards the theoretical upper bound of 1/3 for the mean Nusselt number scaling exponent (see Fig. 3(b)). This increase in the exponent for mean Nusselt number scaling with Rayleigh number from 0.2 towards 1/3 occurs with the development of unsteady flow in the enclosure. Unsteady flow develops above a critical Rayleigh number in the range 3.5  108 [ Ra [ 8.5  108, which appears occur independent of aspect ratio (see Fig. 7). This transition is associated with an increased scaling of maximum streamfunction from Wmax / Ra1/5 to Ra2/5. In the convective regime, the distances of the location of Wmax from the heated boundary and the side wall are found to scale with Ra1/5, suggesting that the plume scales with the boundary layer. In enclosures with D/L > 1 (i.e. tall, slender enclosures), the critical Rayleigh number for transition from diffusion dominated to convection dominated flow is independent of aspect ratio. This likely occurs as a result of the top boundary being sufficiently far from the bottom boundary to wield no influence on the horizontal convection quantities monitored in this study. Acknowledgements The authors are grateful to Dr. Graham Hughes of the Research School of Earth Sciences, ANU, Canberra, Australia, for his insight on streamfunction scaling in horizontal convection. This research was undertaken in part using the NCI National Facility in Canberra, Australia, thanks to a Merit Allocation Scheme grant. NCI is supported by the Australian Commonwealth Government. G.J.S. received financial support through a Faculty of Engineering Small Grant. References [1] [2] [3] [4] [5] [6] [7] [8] [9]

G.O. Hughes, R.W. Griffiths, Horizontal convection, Annu. Rev. Fluid Mech. 40 (2008) 185–208. J.J. Niemela, K.R. Sreenivasan, Confined turbulent convection, J. Fluid Mech. 481 (2003) 355–384. C. Wunsch, R. Ferrari, Vertical mixing, energy, and the general circulation of the oceans, Annu. Rev. Fluid Mech. 36 (2004) 281–314. J.C. Mullarney, R.W. Griffiths, G.O. Hughes, Convection driven by differential heating at a horizontal boundary, J. Fluid Mech. 516 (2004) 181–209. F. Paparella, W.R. Young, Horizontal convection is non-turbulent, J. Fluid Mech. 466 (2002) 205–214. J.H. Siggers, R.R. Kerswell, N.J. Balmforth, Bounds on horizontal convection, J. Fluid Mech. 517 (2004) 55–70. S. Chiu-Webster, E.J. Hinch, J.R. Lister, Very viscous horizontal convection, J. Fluid Mech. 611 (2008) 395–426. G.J. Sheard, T. Leweke, M.C. Thompson, K. Hourigan, Flow around an impulsively arrested circular cylinder, Phys. Fluids 19 (2007) 083601. G.J. Sheard, M.J. Fitzgerald, K. Ryan, Cylinders with square cross section: wake instabilities with incidence angle variation, J. Fluid Mech. 630 (2009) 43– 69. [10] H.M. Blackburn, G.J. Sheard, On quasi-periodic and subharmonic Floquet wake instabilities, Phys. Fluids 22 (2010) 031701. [11] H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, J. Assoc. Comput. Mach. 17 (1970) 589–602. [12] G.O. Hughes, R.W. Griffiths, J.C. Mullarney, W.H. Peterson, A theoretical model for horizontal convection at high Rayleigh number, J. Fluid Mech. 581 (2007) 251–276.