Spectral analysis of boundary layers in Rayleigh-Bénard convection
A combined experimental and numerical study of the boundary layer in a 4:1 aspect-ratio Rayleigh-Bénard cell over a four-decade range of Rayleigh numbers has been undertaken aimed at gaining a better insight into the character of the boundary layers. The experiments involved the simultaneous Laser Doppler Anemometry (LDA) measurements of fluid velocity at two locations, i.e. in the boundary layer and far away from it in the bulk, for Rayleigh numbers varying between and . In parallel, direct numerical simulations (DNS) have been performed for the same configuration for Rayleigh numbers between and . The temperature and velocity probability density functions and the power spectra of the horizontal velocity fluctuations measured in the boundary layer and in the bulk flow are found to be practically identical. Except for the smallest Rayleigh numbers, the spectra in the boundary layer and in the bulk central region are continuous and have a wide range of active scales. This indicates that both the bulk and the boundary layers are turbulent in the Ra number range considered. However, molecular effects can still be observed and the boundary layer does not behave like a classical shear-driven turbulent boundary layer.
pacs:47.20.Bp, 47.27.nb, 44.25.+f
Rayleigh-Bénard (RB) convection in a fluid trapped between two horizontal plates of unequal temperature with the bottom wall being warmer than the top wall, has long served as a paradigm of thermal convection. Despite numerous studies (see e.g. overviews by Siggia (1994); Kadanoff (2001)), it continues to pose challenges because of a number of still unresolved and controversial issues. We recall that RB convection is characterised by the Rayleigh number , and the Prandtl number , where is the thermal expansion coefficient, is the temperature difference between the cold and the hot plate, is the vertical distance between the plates, is the kinematic viscosity, and is the thermal diffusivity. In practical RB set-ups the lateral size, , of the plates is finite. This is taken into account by a third dimensionless number known as the aspect ratio .
A topic of intense debate is the scaling law for the dimensionless heat transfer, , especially for the range (Niemela et al., 2000). Here Nu is the Nusselt number , stands for the convective heat transfer coefficient, is the thermal conductivity, and is the averagedtemperature. Heslot et al. (1987); Castaing et al. (1989); Wu et al. (1990) proposed scaling laws for the heat transfer and temperature statistics in RB convection, which differed from the classical theories, such as that of Howard (1963). In the search for general scaling behaviour, Grossmann and Lohse (2000) presented a unifying theory dividing the plane into several regions, each with different scaling. Crucial assumptions in this theory are the existence of a large-scale circulation, also known as the wind, and a Blasius type laminar boundary layer. The assumption that the boundary layers of Rayleigh-Bénard convection behave as a laminar Blasius boundary layer is supported by several observations: (i) the scaling of the friction factor is compatible with that of a laminar boundary layer past a flat plate at moderate Ra (Chavanne et al., 1997, 2001; Amati et al., 2005); (ii) the typical Reynolds number Re of the flow is low at moderate Ra: for water at and , which is generally accepted to be too low to sustain turbulence. Grossmann and Lohse (2000) argued that at large enough Ra this laminar boundary layer breaks down and becomes turbulent, presuming that the transition to a turbulent boundary layer is shear triggered. The critical Ra number for this transition was estimated to be in the order of for .
Despite this apparently convincing evidence, it is not clear how a boundary layer which is highly unsteady due to continuous plume impingement and detachment can behave quantitatively as a laminar boundary layer. This would indicate that the plumes only cause unsteady disturbances and a consequent time-dependence in the boundary layer but are passive otherwise.
The aim of this paper is to enhance the understanding of the dynamics of the boundary layers by comparing velocity statistics in the boundary layer with those determined in the bulk. It has widely been regarded that the main prerequisite for unsteady velocity fluctuations to be qualified as turbulence, in contrast to other uncorrelated, chaotic and random, fluctuations, is a power spectrum that is continuous over a wide range of scales. For this reason we focused on the analysis of spectral density functions (sdf), complemented with the probability density functions (pdf), aimed at clarifying the nature of the fluctuations in the boundary layer. For this purpose we use in parallel complementary numerical and experimental techniques: direct numerical simulations (DNS) are performed for the lower Ra values, Laser Doppler Anemometry (LDA) measurements are conducted for the medium and higher Ra values. All experiments and simulations combined cover a Ra range between and . In the midrange both techniques overlap over about one decade.
The experiments were conducted in a () cell filled with water (described in more detail in Verdoold et al. (2006)). In brief, two copper plates at the top and bottom are kept at a constant temperature by passing water through the plates’ internal channels, drawn from two basins containing constant temperature water (inaccuracy less than K). The plates thus impose a controlled temperature difference on the working fluid. Throughout this paper we will use a Cartesian coordinate system with its origin at the centre of the bottom wall, Fig. 1. The - and - coordinates are measured along the bottom wall, parallel to the side walls of the RB cell. The -coordinate measures the wall normal distance.
Velocities were measured by using two one-component laser Doppler anemometers manufactured by Dantec. The green ( nm) and blue ( nm) colours of a W Argon-ion laser are used to simultaneously measure the horizontal component in the -direction in the bulk and in the boundary layer. The measured velocity component is in the direction of the large-scale circulation. To enable the detection of instantaneous flow reversals one beam of each colour was frequency pre-shifted by a Bragg cell. The light scattered by m polyamid seed particles was collected in on-axis forward direction, which, in comparison to back-scatter, results in a superior signal-to-noise ratio of the Doppler signals. The photomultiplier output signals were first downmixed electronically to an effective pre-shift of 30 kHz and then fed to two Dantec BSA signals processors to determine the instantaneous velocities. The LDA system measured the instantaneous velocities of the seed particles with an inaccuracy of %. The length and diameter of both measurement volumes are mm and mm, respectively, which is sufficiently small to resolve the smallest length scales of the flow. For each Ra number, the measurement time was at least hours to have sufficient data for an accurate statistical description of the flow.
The mean data rate varied between Hz and Hz, depending on Ra and distance from the wall. For moderate Ra numbers, the relatively low characteristic velocities result in a relatively low data rate. Measurement locations close to the wall have an inherently low data rate due to the low fluid volume flux.
Table 1 shows the experimental parameters as well as the relevant non-dimensional numbers. The Rayleigh number is varied by changing the temperature difference between the bottom and top plates from nearly K at to K at . For the highest Ra numbers, the relatively large temperature differences cause larger density differences, but these remain small enough to neglect non-Boussinesq effects. The Prandtl number will change as well, but the ratio between the smallest and largest Pr number is only 1.7.
The direct numerical simulations cover the lower to moderate Ra number range between and . The numerics are based on staggered finite differences, with a central second order scheme for advection, and a second order Adams-Bashforth time integration scheme. To effectively resolve the boundary layers, grid stretching in the wall-normal direction has been implemented, and special care has been taken to preserve the skew-symmetry of the advective operator (Verstappen and Veldman, 2003). The code, as reported in van Reeuwijk (2007), has been extended to support sidewalls. In the situation without sidewalls, the Poisson equation is Fourier-transformed in the two homogeneous directions. As a result, the Poisson equation decouples for these directions, and the pressure amplitudes can be solved per component by solving a tridiagonal system in the wall-normal direction. To incorporate sidewalls, the Fourier-transform was replaced by a cosine transform, which automatically satisfies a Neumann boundary condition for pressure. The sidewall boundary conditions for the velocity and temperature could be enforced directly due to the explicit time-integration.
Dirichlet boundary conditions are applied at the top and bottom plate, i.e. no-slip and isothermal boundary conditions for velocity and temperature, respectively. On the sidewalls, Neumann boundary conditions are applied for both temperature and velocity, i.e. free-slip and insulating walls for velocity and temperature. The free-slip (stress free) velocity boundary conditions were used to avoid having to resolve the boundary layer at the sidewalls. Although this changes the flow physics near the sidewalls, its global effect is limited for large aspect ratio domains, as preliminary simulations showed. Neumann boundary conditions were applied for pressure at all boundaries.
In the simulations, the properties of water are used except that is kept fixed for all Ra. The grid resolution (Table 2) is sufficient to resolve the smallest turbulent scales, i.e. the Kolmogorov scale and the Corrsin scale . At , the Corrsin scale dominates the grid requirements, and the resolution roughly corresponds to . To resolve the boundary layers at the top and bottom wall properly, grid stretching has been applied in the wall-normal direction, such that at least eight points are within the thermal boundary layer. It should be noted that the simulation at the highest Ra number () may be slightly underresolved in the horizontal directions. High quality simulation as reported in Hanjalić (2005) at and required points per horizontal direction, but reduction to points was the only way to simulate over a sufficiently long time interval to obtain reliable one-point statistics required for the comparison with the experiments. The additional information shown in Table 2 are the dimensionless timestep , the number of simulated turnovers and the average Nusselt number Nu. The typical convective turnover time with the typical velocity based on the free fall scaling. The average Nusselt number is calculated as , with averaged over the entire bottom surface and over many turnovers.
ii.3 Measurement locations
Simultaneous measurements have been performed in the boundary layer and in the bulk of the cell for various Ra. The measurement location in the bulk is at , with a fixed wall distance , see Fig. 1. The measurement location in the boundary layer has the same - and -coordinates, but the wall distance is different, i.e., , where is the kinematic boundary layer thickness determined as in Grossmann and Lohse (2000):
Table 3 lists the physical measurement locations, where the DNS values have been rescaled to . Figure 2 shows that the wall distance for the boundary layer measurements slightly deviates from linear behaviour (on a double-log plot) due to minor changes in Pr for the experimental points. In separate experiments, LDA was used to determine mean-velocity and rms-velocity profiles above the wall (at ). The results of these experiments confirmed that the measurement locations at were all located in the boundary layer for the full Ra number range considered in this study. The DNS values for , denoted by solid circles in Fig. 2, support this observation. Clearly, (1) overestimates the scaling exponent of , so that the prefactor in the BL measurement points is a decreasing function of Ra. Hence, as Ra increases, the measurements are performed deeper into the kinematic boundary layer. The disparity between the theoretical prediction (1) and the measurements for will be discussed in more detail at the end of section III.2.
The DNS time series are short in respect to the LDA series. To enable a comparison between LDA and DNS results, the DNS data are extracted from four points in the domain. These points are symmetrically placed around the origin with horizontal coordinates , , as shown in Fig. 1. The four obtained time series are regarded as nearly statistically independent, and are then averaged to effectively improve the statistics of the DNS time series.
iii.1 Probability density functions
Figure 3 shows the measured and computed variance of the horizontal velocity in the boundary layer and in the bulk. The graph shows that for increasing Ra there is a gradual transition of the variances between DNS and LDA results giving confidence in the procedure followed.
Temperature pdfs obtained from DNS for the boundary layer and the bulk are depicted in Figure 4(b) and show an exponential distribution even for . The straight lines (on logarithmic scale) are generally associated with the hard turbulent regime, which sets in at for aspect-ratio-1 cells (e.g. Heslot et al. (1987); Castaing et al. (1989); Qiu and Tong (2001)). Recall that the transition between soft and hard turbulence takes place at relatively low Ra number in the aspect-ratio-4 cell, see Kerr (1996). The boundary layer pdf is skewed towards positive temperatures. This is generally attributed to the interaction of thermal plumes with the flow. Plumes leaving the boundary layer result in positive temperature excursions.
Figure 5 depicts the probability density functions (pdfs) of the computed and measured component. As expected, the pdfs of both horizontal velocity components were found to have the same shape in the simulations. For brevity only the component is therefore depicted in Figure 5. For a better comparison, the pdfs are scaled with the rms value . The results show that the pdfs in both the boundary layer and in the bulk exhibit a nearly Gaussian-like shape for the complete Ra number range. Only the pdfs at deviate largely from this trend, mainly due to the presence of a strong oscillatory component in the velocity time series, as will be shown later when discussing the power spectra.
Vincent and Meneguzzi (1991); Mouri et al. (2002) make a distinction between sub-Gaussian, Gaussian and hyper-Gaussian stages in the development of turbulence. They argue that the nearly Gaussian shape of the pdfs in the present study indicates that turbulence is in a fully developed stage. The pdfs for the highest Ra values have non-zero mean values. This is caused by the large-scale circulation or “wind” within the RB cell. For the mean value is approximately . The wind is rather weak in the sense that its mean value is small compared to the magnitude of the rms velocity fluctuations, in accordance with other findings reported in van Reeuwijk et al. (2005); Verdoold et al. (2006); Niemela and Sreenivasan (2006); van Reeuwijk et al. (2007a). It is also observed that the wind is not always in the same direction. In some experiments (Fig. 5LABEL:) it appears to be clockwise, while in others (Figs. 5LABEL:, LABEL:) it is counterclockwise.
iii.2 Spectral density functions
The obtained velocity time series are investigated further by computing the spectral density function (sdf). The sdf is defined as the Fourier transform of the autocorrelation function , where () is the fluctuating component of the turbulent velocity, is time, is the time lag, and the overline denotes ensemble averaging. All sdfs are calculated using the same algorithm, to make a fair comparison for both DNS and LDA data. This algorithm, described in Tummers and Passchier (2001), can account for the randomness of the sampling times in LDA measurements.
Figure 6 shows the sdfs of the horizontal component in the -direction determined in the boundary layer and in the bulk for the full range of Ra considered in this study. The sdfs determined from DNS data are accurate at high frequencies, but less reliable in the low frequency range due to the relatively small measurement times. In contrast, the low frequency region in the LDA sdfs is accurate, while the high frequency range suffers from a relatively high statistical scatter caused by the random sampling process in LDA. The level of the statistical scatter (which strongly depends on the mean data rate) sets a lower limit to the spectral density that can be accurately determined from LDA data.
When comparing the various sdfs in Figure 6 it is clear that the sdf for is different from all others. At this low Ra numbers, the flow is not turbulent but in a state of spatio-temporal chaos. As mentioned earlier, the velocity time series is characterised by strong oscillations. As a consequence, the sdf in Fig. 6(a) (and to a lesser extent the sdf in Fig. 6(b)) is dominated by a limited number of discrete peaks and the sdf rolls off rapidly at higher frequencies.
Peaks in the sdf also occur in the very low frequency range in Figs. 6(i) and 6(j). The peak frequencies correspond to periods of s and s at and , respectively. These oscillations are unrelated to the ones observed at low Ra, and have their origin in the alternating growth and decay of rolls that form the large scale circulation as discussed in Verdoold et al. (2006).
For all sdfs measured in the bulk have a wide range of active scales, and this range increases with Ra. These sdfs are continuous, and are not dominated by a small number of frequencies as one might expect in spatio-temporal chaos. Clearly these sdfs indicate that the bulk in the RB cell is turbulent for . Although an inertial subrange is absent for the measurements at lower Ra, the sdf at (Fig. 6(j)) shows an inertial subrange spanning roughly one decade. Note that the presence of an inertial subrange is not a necessary condition for turbulence. An inertial subrange can only form when there is a range of scales which is neither influenced by macro- or micro-scales, i.e. at high Re when the production and dissipation scales are sufficiently separated. Therefore, an inertial subrange is absent despite the turbulence at lower Re (see also Wu et al., 1990; Camussi and Verzicco, 1998).
The most striking observation is that the sdfs in the bulk and in the boundary layer are practically identical for . As the sdfs for the bulk indicate that the flow is turbulent, one cannot escape from the conclusion that the flow in the boundary layers must also be characterised as turbulent. It is difficult to reconcile this observation with the conceptual image of a laminar, Blasius-type boundary layer along the walls of the RB cell. For a Blasius boundary layer, the laminar-turbulent transition occurs at a relatively high Reynolds number, i.e., depending on the intensity and nature of disturbances in the flow (Schlichting and Gersten, 1997). Here, is a Reynolds number based on plate length and outer velocity. However, even at the highest Ra considered in this study, the value of based on the cell width and the wind velocity is only , which is far below the critical value for laminar-turbulent transition in a flat plate boundary layer. Thus, the turbulence in the boundary layer cannot be due to a shear-triggered transition (occuring inside the boundary layer), but more likely due to the impingement and detachment of plumes (travelling into and out of the boundary layer), which is corroborated by the great similarity between the bulk and the boundary-layer spectra.
Naturally, one can always envision that the boundary layer only quantitatively behaves as a laminar Blasius boundary layer: the plumes introduce a time-dependence on time-scales that prevent pure laminarity. The underlying assumption is that on the average, the thermals and plumes are passive when it comes to scaling of various important integral parameters of the kinematic boundary layer such as the friction factor and the boundary layer thickness . Indeed, scales as , as one would expect for a laminar boundary layer past a flat plate at moderate Ra (Chavanne et al., 1997, 2001; Amati et al., 2005; van Reeuwijk et al., 2007b).
The scaling of the boundary layer thickness is slightly more complicated. Assuming a laminar boundary layer, scales as (Schlichting and Gersten, 1997), which is related to Ra as upon assuming . Experiments show that a scaling exponent holds well for the boundary layer thickness of the sidewalls (Qiu and Xia, 1998a), but not for the scaling of the bottom and top boundary layers (Qiu and Xia, 1998b; Kerr, 1996; van Reeuwijk et al., 2007b) which scales with exponents from to . The simulations in this paper confirm the weak Ra number dependence of (Fig. 2). This disparity indicates that plumes impinging on and emerging from the boundary layer are not entirely passive, but actively influence the boundary layer thickness. Given the important role of in the theory of Grossmann and Lohse (2000), it seems to us that a detailed study of how the interplay of plumes, the wind and the boundary layer influences is desirable. A first step in this direction is taken in a forthcoming study van Reeuwijk et al. (2007a, b), of which the results indicate that the scaling of cannot be derived from laminar boundary layer theory and scales as .
Although the sdfs point towards the presence of turbulence in the boundary layer, the characteristics of the boundary layer are of a rather different nature than those of classical (i.e. shear-driven) turbulent boundary layers. Indeed, viscous forces are still significant in the boundary layer, as can be deduced from the domain height in viscous units , where is the friction velocity Schlichting and Gersten (1997). Assuming that , the ratio which is the shear-Reynolds number , can be approximated by
The simulations show that at , and the viscous thickness of the boundary layer is only a fraction of this already very low value. Therefore, viscous effects are significant in the boundary layer.
It is interesting to note that the scaling of the boundary layer thickness for flat plate flow changes from (laminar) to (turbulent, see White (1991)), where is the Reynolds number based on the distance to the leading edge of the plate, and the thickness is the location where the velocity is 99% of the outer velocity. The decreased exponent resembles larger entrainment and thus a quicker growth of the boundary layer as a function of . Upon assuming that the typical development length is and that , the anticipated scaling for a turbulent boundary layer is . If is the scaling exponent from the simulations, then , so the exponent may resemble an intermediate state.
However, it is not clear whether the entrainment mechanism of a developing boundary layer past a flat plate is transferable to RB convection; the active role of plumes which transfer fluid (and heat) into and out of the boundary layer may significantly affect the entrainment characteristics. It is our suspicion that the relatively weak forcing in the horizontal direction due to the wind (which is responsible for the friction factor ) is only weakly coupled to the relatively strong forcing in the wall-normal direction due to buoyancy. This may be the reason why the boundary layer has an sdf which is typical of developed turbulence but has many laminar features at the same time. Due to the weak coupling, the conceptual image of a laminar boundary layer may be appropriate for but not for in the boundary layer in RB convection. Therefore, the straightforward application of well-known characteristics of the Blasius-type boundary layer (such as the criteria for laminar-turbulent transition and the scaling of with Re) to the boundary layers in RB convection is questionable without additional justification.
A combined numerical and experimental investigation was carried out to study the boundary layers in RB convection in a 4:1 aspect-ratio cell filled with water. The study covers a wide range of Ra numbers varying between and . The results for the lower Ra numbers were obtained from direct numerical simulations, while those for the higher Ra numbers followed from laser Doppler anemometry. In the midrange both techniques have an overlap of about one decade.
The probability density functions of the horizontal velocity components have a Gaussian-like shape for the complete Ra number range, except at the lowest Ra where the flow is in a state of spatio-temporal chaos. Probability density functions that are scaled with the velocity rms value overlap for the entire Ra number range, with the exception of the largest Ra numbers where the large-scale circulation, or “wind”, is predominantly visible in the histograms.
The spectral density functions of the horizontal velocity component that were determined in the bulk flow and in the boundary layer are surprisingly similar. Spectra in both regions of the flow are continuous and have a wide range of active scales for . This indicates that both the bulk and the boundary layers are turbulent in the Ra number range considered, while (i) molecular effects are still noticeable in the boundary layer and (ii) the boundary layer does not behave like a classical shear-driven turbulent boundary layer.
- Siggia (1994) E. D. Siggia, Annu. Rev. Fluid Mech. 26, 137 (1994).
- Kadanoff (2001) L. P. Kadanoff, Phys. Today 54, 34 (2001).
- Niemela et al. (2000) J. J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly, Nature 404, 837 (2000).
- Heslot et al. (1987) F. Heslot, B. Castaing, and A. Libchaber, Phys. Rev. A 36, 5870 (1987).
- Castaing et al. (1989) B. Castaing, G. Gunaratne, F. Heslot, L. Kadanoff, A. Libchaber, S. Thomae, X.-Z. Wu, S. Zaleski, and G. Zanetti, J. Fluid Mech. 204, 1 (1989).
- Wu et al. (1990) X.-Z. Wu, L. P. Kadanoff, A. Libchaber, and M. Sano, Phys. Rev. Lett. 64, 2140 (1990).
- Howard (1963) L. N. Howard, J. Fluid Mech. 17, 405 (1963).
- Grossmann and Lohse (2000) S. Grossmann and D. Lohse, J. Fluid Mech. 407, 27 (2000).
- Chavanne et al. (1997) X. Chavanne, F. Chillà, B. Castaing, B. Hébral, B. Chabaud, and J. Chaussy, Phys. Rev. Lett. 79, 3648 (1997).
- Chavanne et al. (2001) X. Chavanne, F. Chillà, B. Chabaud, B. Castaing, and B. Hébral, Phys. Fluids 13, 1300 (2001).
- Amati et al. (2005) G. Amati, K. Koal, F. Massaioli, K. R. Sreenivasan, and R. Verzicco, Phys. Fluids 17, 121701 (2005).
- Verdoold et al. (2006) J. Verdoold, M. J. Tummers, and K. Hanjalić, Phys. Rev. E 73, 056304 (2006).
- Verstappen and Veldman (2003) R. W. C. P. Verstappen and A. E. P. Veldman, J. Comput. Phys. 187, 343 (2003).
- van Reeuwijk (2007) M. van Reeuwijk, Ph.D. thesis, Delft University of Technology (2007), URL http://repository.tudelft.nl/file/525273/372306.
- Hanjalić (2005) K. Hanjalić, Int. J. Heat Fluid Fl. 26, 828 (2005).
- Grossmann and Lohse (2001) S. Grossmann and D. Lohse, Phys. Rev. Lett. 86, 3316 (2001).
- Qiu and Tong (2001) X.-L. Qiu and P. Tong, Phys. Rev. E 64, 036304 (2001).
- Kerr (1996) R. M. Kerr, J. Fluid Mech. 310, 139 (1996).
- Vincent and Meneguzzi (1991) A. Vincent and M. Meneguzzi, J. Fluid Mech. 225, 1 (1991).
- Mouri et al. (2002) H. Mouri, M. Takaoka, A. Hori, and Y. Kawashima, Phys. Rev. E 65, 056304 (2002).
- van Reeuwijk et al. (2005) M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalić, Phys. Fluids 17, 051704 (2005).
- Niemela and Sreenivasan (2006) J. J. Niemela and K. R. Sreenivasan, J. Fluid Mech. 557, 411 (2006).
- van Reeuwijk et al. (2007a) M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalić, Submitted to Phys. Rev. E (2007a), URL http://arxiv.org/abs/0709.0304.
- Tummers and Passchier (2001) M. J. Tummers and D. M. Passchier, Meas. Sci. Technol. 12, 1641 (2001).
- Camussi and Verzicco (1998) R. Camussi and R. Verzicco, Phys. Fluids 10, 516 (1998).
- Schlichting and Gersten (1997) H. Schlichting and K. Gersten, Grenzschicht-Theorie (Springer-Verlag Berlin Heidelberg, 1997), 9th ed.
- van Reeuwijk et al. (2007b) M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalić, Submitted to Phys. Rev. E (2007b), URL http://arxiv.org/abs/0709.1891.
- Qiu and Xia (1998a) X.-L. Qiu and K.-Q. Xia, Phys. Rev. E 58, 486 (1998a).
- Qiu and Xia (1998b) X.-L. Qiu and K.-Q. Xia, Phys. Rev. E 58, 5816 (1998b).
- White (1991) F. M. White, Viscous fluid flow (McGraw-Hill, 1991).