"Science, Faculty of"@en . "Earth, Ocean and Atmospheric Sciences, Department of"@en . "DSpace"@en . "UBCV"@en . "DeGiuli, Eric"@en . "2009-09-16T13:18:54Z"@en . "2009"@en . "Master of Science - MSc"@en . "University of British Columbia"@en . "The problem of turbulent \now in a rough pipe of arbitrary shape is considered.\nThe classical Izakson-Millikan argument for a logarithmic velocity\nprofile is presented, and matched asymptotic expansions are introduced.\nScaled, dimensionless equations are produced and simplified. A simple mixing\nlength turbulence model is presented, which closes the problem. To calibrate\nthe model, the mechanical problem is solved in the case of a circular\npipe. Excellent agreement with engineering relations is obtained. The mechanical\nproblem for a non-circular pipe is posed, and the boundary layer\nproblem is solved. This leaves unknown the wall stress, which is sought\nthrough approximate methods of solution in the outer region. These are\npresented and the approximate solutions thus obtained are compared to full\nnumerical solutions and data for a square, elliptical, and semi-elliptical pipe.\nThe approximations are vindicated, but agreement between the numerical\nsolutions and data is only moderate. Discrepancies are explained in terms\nof the neglected secondary \now.\nThe thermal problem is posed, with scalings taken for intended application\nin glaciology. The problem is solved for a circular pipe. Heat transfer\nresults are presented and compared with empirical relations. The general\nproblem for a non-circular pipe is posed, and approximate methods of solution\nare motivated, in analogy to those used for the mechanical problem.\nThese are used to obtain approximate solutions, which are compared with\nnumerical solutions, to good agreement. Possible applications of these solutions\nare discussed."@en . "https://circle.library.ubc.ca/rest/handle/2429/12802?expand=metadata"@en . "1964228 bytes"@en . "application/pdf"@en . "Turbulent Flow in Geophysical Channels by Eric DeGiuli B.Sc., University of Toronto, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Geophysics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) September, 2009 c\u00C2\u00A9 Eric DeGiuli 2009 Abstract The problem of turbulent flow in a rough pipe of arbitrary shape is con- sidered. The classical Izakson-Millikan argument for a logarithmic velocity profile is presented, and matched asymptotic expansions are introduced. Scaled, dimensionless equations are produced and simplified. A simple mix- ing length turbulence model is presented, which closes the problem. To cal- ibrate the model, the mechanical problem is solved in the case of a circular pipe. Excellent agreement with engineering relations is obtained. The me- chanical problem for a non-circular pipe is posed, and the boundary layer problem is solved. This leaves unknown the wall stress, which is sought through approximate methods of solution in the outer region. These are presented and the approximate solutions thus obtained are compared to full numerical solutions and data for a square, elliptical, and semi-elliptical pipe. The approximations are vindicated, but agreement between the numerical solutions and data is only moderate. Discrepancies are explained in terms of the neglected secondary flow. The thermal problem is posed, with scalings taken for intended applica- tion in glaciology. The problem is solved for a circular pipe. Heat transfer results are presented and compared with empirical relations. The general problem for a non-circular pipe is posed, and approximate methods of so- lution are motivated, in analogy to those used for the mechanical problem. These are used to obtain approximate solutions, which are compared with numerical solutions, to good agreement. Possible applications of these solu- tions are discussed. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . 9 1.1 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Turbulence Closure . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Circular Pipe Solution . . . . . . . . . . . . . . . . . . . . . . . 16 4 An Approximation Scheme for Non-Circular Pipes . . . . 22 4.1 Constant-stress Approximation . . . . . . . . . . . . . . . . . 29 4.1.1 Square Pipe . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Zero-stress Approximation . . . . . . . . . . . . . . . . . . . 34 4.2.1 Elliptical Pipe . . . . . . . . . . . . . . . . . . . . . . 35 4.2.2 Semi-Elliptical Pipe . . . . . . . . . . . . . . . . . . . 38 5 The Thermal Problem . . . . . . . . . . . . . . . . . . . . . . . 39 5.1 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2 Circular Pipe . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.3 Non-Circular Pipe . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3.1 Constant Heat Flux Approximation . . . . . . . . . . 50 iii Table of Contents 5.3.2 Zero Heat Flux Approximation . . . . . . . . . . . . . 51 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Appendices A Exact Solution for the Velocity Profile in a Circular Pipe 60 B Geometrical Transformations . . . . . . . . . . . . . . . . . . 61 B.1 Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 B.2 Wall Coordinates . . . . . . . . . . . . . . . . . . . . . . . . 63 B.3 Coordinate Derivation . . . . . . . . . . . . . . . . . . . . . . 64 B.4 Example: Anisotropic Superellipses . . . . . . . . . . . . . . 65 C Velocity Profile in an Infinite Plane Channel . . . . . . . . 68 D Viscous Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 iv List of Tables 6.1 Glossary of Symbols . . . . . . . . . . . . . . . . . . . . . . . 56 v List of Figures 1 Hydraulic radius analogy applied to Blasius\u00E2\u0080\u0099 mean velocity correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Sketch of secondary flow in pipes of triangular and rectangular cross-section . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 Sketch of idealized wall roughness possibilities. . . . . . . . . 5 1.1 Sketch of geometry. . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1 Mixing-length l(n). . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Exact, inner, outer, composite, and law-of-the-wall solutions for the mean velocity profile in a fully rough, circular pipe, in semi-log coordinates. . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 Exact, inner, outer, composite, and law-of-the-wall solutions for the mean velocity profile in a fully rough, circular pipe. . 21 4.1 Sketch of (n, \u00CE\u00B8) coordinates. . . . . . . . . . . . . . . . . . . . 22 4.2 Examples of anisotropic superellipses. . . . . . . . . . . . . . 28 4.3 Wall stress for a square pipe; secondary flow velocity. . . . . . 32 4.4 Streamwise velocity in a square pipe. . . . . . . . . . . . . . . 33 4.5 Wall stress for an elliptical pipe; secondary flow velocity. . . . 36 4.6 Streamwise velocity in an elliptical pipe. . . . . . . . . . . . . 37 4.7 Wall stress for a semi-elliptical pipe . . . . . . . . . . . . . . . 38 5.1 Stanton number correlation . . . . . . . . . . . . . . . . . . . 48 5.2 Heat flux in a square pipe. . . . . . . . . . . . . . . . . . . . . 52 5.3 Heat flux in a semi-elliptical pipe. . . . . . . . . . . . . . . . 53 B.1 Curvature sketch . . . . . . . . . . . . . . . . . . . . . . . . . 62 B.2 Derivation of G(n, \u00CE\u00B8) . . . . . . . . . . . . . . . . . . . . . . . 63 B.3 Preimage of \u00E2\u0084\u00A6 in (n, \u00CE\u00B8) plane . . . . . . . . . . . . . . . . . . 66 D.1 Wall roughness boundary layers . . . . . . . . . . . . . . . . . 69 D.2 Transition from smooth wall to rough wall . . . . . . . . . . . 71 vi List of Figures D.3 Boundary layer thickness in the transitional regime. . . . . . 73 vii Acknowledgements It is trite, but true, that one\u00E2\u0080\u0099s first experience in graduate school is a be- wildering and difficult, but ultimately satisfying, formative experience. My past thirty months in Vancouver can only be described as such. To all whose paths I\u00E2\u0080\u0099ve crossed in this journey, I owe thanks. It is difficult to separate pro- fessional from personal acknowledgments, as I have been fortunate enough to see both sides of my supervisors and colleagues. For this, I am indebted to Garry Clarke, who has forged a uniquely social glaciology group at UBC. To all its members, past and present\u00E2\u0080\u0093 Christian Schoof, Alex Jarosch, Faron Anslow, Christian Reuten, Tim Creyts, Andrew Schaeffer, Tanya Stickford\u00E2\u0080\u0093 I offer thanks for stimulating conversation and, often, fine cuisine. To Garry I also owe immense gratitude for helping me navigate the once-unknown worlds of graduate school and glaciology. Somehow, Garry was never too busy to answer my questions. As time progressed, my waxing interest in taking a mathematical ap- proach to geophysical problems led to discussions with, and eventually for- mal supervision by, Christian Schoof, who helped me understand the utility and importance of asymptotic methods, which form the basis for this thesis. As a new professor and head of a nascent field program, Christian was (and is) immensely busy, but he always found time in his schedule to answer my incessant queries and cryptic remarks. Still, it would have been impossible to reconcile abstract mathematical constructs with physical reality without the guiding insight of Neil Balmforth, whose group meetings provided ample opportunity for the presentation of my research-in-progress. Thanks, Neil\u00E2\u0080\u0093 I hope the next two years will be just as inspiring. Neil\u00E2\u0080\u0099s meetings were also an opportunity to meet fellow students working on similar problems in Mathematics and Mechanical Engineering. Thanks to Jeff Carpenter, Ian Chan, Neville Dubasch, and Anja Slim for making these enjoyable learning experiences. Finally, modern science inevitably means dealing with com- puter problems\u00E2\u0080\u0093 thanks to Cornel Pop and Marek Labecki for helping me resolve those I\u00E2\u0080\u0099ve had. viii To those who helped me understand: When one is young, one venerates and despises without that art of nuances which constitutes the best gain of life, and it is only fair that one has to pay dearly for having assaulted men and things in this manner with Yes and No. Everything is arranged so that the worst of tastes, the taste for the unconditional, should be cruelly fooled and abused until a man learns to put a little art into his feelings and rather to risk trying even what is artificial \u00E2\u0080\u0093 as the real artists of life do. The wrathful and reverent attitudes characteristic of youth do not seem to permit themselves any rest until they have forged men and things in such a way that these attitudes may be vented on them \u00E2\u0080\u0093 after all, youth in itself has something of forgery and deception. Later, when the young soul, tortured by all kinds of disappointments, finally turns suspiciously against itself, still hot and wild, even in its suspicion and pangs of conscience \u00E2\u0080\u0093 how wroth it is with itself now! how it tears itself to pieces, impatiently! how it takes revenge for its long self-delusion, just as if it had been a deliberate blindness! In this transition one punishes oneself with mistrust against one\u00E2\u0080\u0099s own feelings; one tortures one\u00E2\u0080\u0099s own enthusiasm with doubts; indeed, one experi- ences even a good conscience as a danger, as if it were a way of wrapping oneself in veils and the exhaustion of subtler honesty \u00E2\u0080\u0093 and above all one takes sides, takes sides on principle, against \u00E2\u0080\u009Cyouth.\u00E2\u0080\u009D \u00E2\u0080\u0093 Ten years later one comprehends that all this, too\u00E2\u0080\u0093 was still youth. \u00E2\u0080\u0094Friedrich Nietzsche[19] ix Introduction From river meandering to glacial outburst floods to lava tube formation, myriad geophysical phenomena result from the interaction of a turbulent flow with boundary processes, whether these be erosion, melting, or so- lidification. To theoretically study the resulting evolution of the channel shape, it is necessary to construct models which do not impose the plane channel or cylindrical pipe assumptions that are widely used in engineering applications. Because the phenomena of interest involve physical processes in addition to those of turbulent flow, it is desirable to have simple mod- els that capture the most important geometrical properties of a non-trivial channel shape. It is the goal of this thesis to provide such models. Turbulence is the spatially and temporally chaotic motion of a fluid that spontaneously results when inertia dominates viscosity. Predicting turbulent phenomena directly from the Navier-Stokes equations is widely considered the greatest problem of classical physics. For geophysical and engineering applications of turbulent flow, usually it is mean, rather than instantaneous, quantities that are of interest. In this context, the \u00E2\u0080\u009Cproblem of turbulence\u00E2\u0080\u009D is to satisfactorily deal with the loss of information inherent in averaging. To illustrate, consider the Reynolds averaging procedure, whereby flow quanti- ties are considered as sums of mean and fluctuating components, viz., ui = u\u00CC\u0084i + u\u00E2\u0080\u00B2i, P = P + P \u00E2\u0080\u00B2, for velocity and pressure, respectively. When these expressions are substi- tuted into the Navier-Stokes (NS) equations and time averaged, one obtains the Reynolds-averaged Navier-Stokes equations1[23]: \u00E2\u0088\u0082u\u00CC\u0084i \u00E2\u0088\u0082xi = 0 u\u00CC\u0084j \u00E2\u0088\u0082u\u00CC\u0084i \u00E2\u0088\u0082xj = \u00E2\u0088\u00921 \u00CF\u0081 \u00E2\u0088\u0082P\u00CC\u0084 \u00E2\u0088\u0082xi + \u00CE\u00BDM \u00E2\u0088\u00822u\u00CC\u0084i \u00E2\u0088\u0082xj\u00E2\u0088\u0082xj \u00E2\u0088\u0092 \u00E2\u0088\u0082 \u00E2\u0088\u0082xj ( u\u00E2\u0080\u00B2iu \u00E2\u0080\u00B2 j ) . (1) 1To perform the time averaging, one assumes that fluctuating components have zero- mean, and that mean quantities do not vary on the averaging time scale. 1 Introduction These take the form of the original NS equations, except for the additional terms involving \u00CF\u0081u\u00E2\u0080\u00B2iu \u00E2\u0080\u00B2 j , known as Reynolds stresses. The equations are un- closed, since there are ten unknowns but only four equations. Therefore, additional hypotheses on the Reynolds stresses are needed. Finding them is referred to as the closure problem. In this thesis, we will consider partic- ularly simple forms of closure, known as eddy viscosity models. These are discussed in chapter 2. Because of the closure problem, a large body of experimental work has attempted to quantify mean properties of turbulent flow using empirical relations. Since non-circular pipes and non-planar channels arise frequently in engineering applications, empirical relations have been extended to these cases, often using the hydraulic radius analogy. This states that empirical correlations for the mean velocity of turbulent flow in a circular pipe may be used for any pipe shape when the radius is replaced by the hydraulic radius R \u00E2\u0089\u00A1 2A Pw , where A is the cross-section\u00E2\u0080\u0099s area, and Pw its perimeter. For example, Figure 1 shows the validity of Blasius\u00E2\u0080\u0099 mean velocity correlation for pipes of triangular, rectangular, and annular cross-section. A related analogy is also valid for open-channel flows, but the hydraulic radius is defined without the factor of two, and the perimeter taken is the \u00E2\u0080\u009Cwetted perimeter,\u00E2\u0080\u009D the total less that of the free surface. The hydraulic radius analogy is very successful in computing mean ve- locities and temperatures, but it says nothing about the azimuthal variation of these quantities, which is essential for the aforementioned geophysical problems. For example, the evolution of a subglacial meltwater channel de- pends on the variation of heat flux around the channel. To date, very little theoretical work has probed the structure of azimuthal variation. A major difficulty is the ubiquitous presence of secondary flow. Prandtl observed that in pipes of any non-circular cross-section, there always exists streamwise vorticity, which he termed \u00E2\u0080\u009Csecondary flow of the second kind\u00E2\u0080\u009D[24].2 In regular geometries, the vorticity forms a pattern of steady, stable rolls which are divided by the lines of symmetry of the cross- sectional shape. (see Figure (2)) As shown by Launder and Ying [12], this motion cannot be predicted by the simplest eddy viscosity turbulence closures. Therefore, with few exceptions, research on secondary flow has been experimental and numerical, 2The \u00E2\u0080\u009Cfirst kind\u00E2\u0080\u009D is streamwise vorticity caused by streamwise curvature of the pipe. 2 Introduction Figure 1: Hydraulic radius analogy applied to Blasius\u00E2\u0080\u0099 mean velocity correla- tion; from Schlichting [26]. The plots of different shapes have been staggered for clarity. Note that the analogy fails for laminar flow. Figure 2: Sketch of secondary flow in pipes of triangular and rectangular cross-section; from Schlichting [26]. using advanced Reynolds stress closure models and, more recently, direct numerical simulation (DNS). This is somewhat discouraging, since the Reynolds stress models are largely empirical, especially in the treatment of the near-wall region.[23] Furthermore, for those models that accurately predict the secondary flow, when the employed constitutive relations are tested a priori with DNS data, 3 Introduction they can fail profoundly, indicating that the right answer is often obtained for the wrong reasons.[17] For many applications, the secondary flow is only of interest to the extent to which it affects the azimuthal variation of wall quantities.3 In this case, one may attempt to avoid entirely the complexity of modelling the secondary flow. Since the secondary flow velocities amount to only a few percent of the streamwise velocity, apparently independent of Reynolds number [8], one may consider the secondary flow as a perturbation on the main stream- wise flow, and pose a series solution for velocity in a suitable dimensionless parameter. To pursue this approach, in the next section we use simple similarity con- siderations to construct a model based on local equilibrium, after Townsend [31]. Using this \u00E2\u0080\u009Clocal model\u00E2\u0080\u009D we infer azimuthal variation of velocity and temperature, which can be used to compute wall stress and wall heat flux. If desired, the secondary flow can be computed as higher-order terms in this expansion. The geophysical channels under consideration are typically rough. That is, on the finest scales, from centimetres down to microns, small asperities disrupt the viscous shear layer at the surface. To avoid the difficulty of explicitly modelling the flow around the asperities, we imagine an ideal- ized situation in which the given cross-sectional shape has the roughness smoothed out, and we attempt to parameterize the roughness by one or several length scales [10, 25]. To account for the effects of asperities on the flow, the roughness parameters must be incorporated into the constitutive model of turbulence. In this thesis, we will consider only a single rough- ness length scale.4 By calibration with the data of Nikuradse [21], this can be identified with an effective sand diameter, in reference to Nikuradse\u00E2\u0080\u0099s classical experiments on pipes with sand glued to their surface. Even with a single roughness length scale, there are many ways of defin- ing the boundary of the idealized smoothed pipe. For a circular pipe, choos- ing the boundary is a one degree of freedom problem; engineers usually set this \u00E2\u0080\u009Cwall offset\u00E2\u0080\u009D by calibrating the velocity profile to a logarithmic form [28], as discussed below. However, it is preferable for the geometry to be defined independently of flow data. Herwig [6] argued that a natural way of placing the boundary is to ensure that the volume of the idealized pipe is equal to the volume of the true, rough pipe, as illustrated in Figure 3. Since this may easily be extended to non-circular pipes, we will adopt this 3An exception is river meandering, where the secondary flow advects sediment. 4The effect of roughness type on the roughness parameters is reviewed in ([6, 10, 25]). 4 Introduction method. Figure 3: Sketch of idealized wall roughness possibilities parameterized by a single length scale h. The volume-preserving equivalent smooth wall (esw) is shown as a dashed line. Strictly, the placement of this line depends on the large-scale curvature of the boundary; shown here is the plane-channel case. However, in applications where the roughness is small compared to the radius of curvature, the line tends to its plane-channel position. From Herwig et al [6]. The problem of wall-bounded turbulence is formidable. Even without the novelties discussed above, the motion of a turbulent fluid through a smooth circular pipe is a subject of intense research interest.[13, 23] One simplifi- cation that can be used to great advantage in theoretical work exploits the boundary layer nature of the flow. At high Reynolds number, nearly all of the shear in the mean velocity is concentrated in a thin boundary layer near the wall of the channel, while the bulk of the flow moves nearly uniformly. This is mathematically expressed as a separation of the domain into dis- tinct regions where different terms in the governing equations are dominant. In particular, viscous stresses are comparable to Reynolds stresses in the boundary layer, but negligible in the bulk of the flow. This observation, due to Prandtl, led to the development of boundary layer theory, a major advance in the study of wall-bounded turbulence [26]. Millikan [16] and Izakson [9] used the boundary layer observation, along with dimensional arguments, to estimate the form of the velocity profile close to the wall. Because their simple argument leads naturally to the more sophisticated technique of matched asymptotic expansions, which features heavily in this thesis, we reproduce it here. For simplicity, we will consider a smooth, circular pipe with an imposed pressure gradient \u00E2\u0088\u0086P/\u00E2\u0088\u0086x. It is 5 Introduction conventional to use wall units u+ = u\u00CC\u0084 U\u00CF\u0084 , y+ = (R\u00E2\u0088\u0092 r) \u00CE\u00B4v , with \u00CE\u00B4v = \u00CE\u00BDM/U\u00CF\u0084 a viscous length scale. The velocity scale U\u00CF\u0084 is the friction velocity U\u00CF\u0084 \u00E2\u0089\u00A1 \u00E2\u0088\u009A \u00CF\u0084w/\u00CF\u0081 = \u00E2\u0088\u009A R 2\u00CF\u0081 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0086P\u00E2\u0088\u0086x \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3, where \u00CF\u0084w is the mean wall stress, and the second equality follows from global momentum balance. By the Buckingham Pi theorem, the wall-normal velocity derivative may be written du+ dy+ = 1 y+ \u00CF\u0086 ( y+, y R ) . (2) Since the viscous stress is unimportant in the outer layer, one may suppose that du/dy is independent of \u00CE\u00BDM there: \u00CF\u0086 ( y+, y R ) \u00E2\u0086\u0092 \u00CF\u0086O ( y R ) y/R \u00E2\u0088\u00BC 1. (3) Conversely, within the thin viscous boundary layer, one may suppose that the geometry of the pipe should be unimportant; that is, R should not appear in the dimensional form of du/dy. Then5 \u00CF\u0086 ( y+, y R ) \u00E2\u0086\u0092 \u00CF\u0086I ( y+ ) y+ \u00E2\u0088\u00BC 1. (4) Izakson[9] and Millikan[16] proposed, independently, that there should be an overlap region in which both forms of \u00CF\u0086 are valid. This is only possible if \u00CF\u0086 ( y+, y R ) \u00E2\u0086\u0092 1 \u00CE\u00BA , (5) for some constant \u00CE\u00BA, in the overlap layer. Then (2) can be integrated to yield u+ \u00E2\u0088\u00BC 1 \u00CE\u00BA log(y+) +B (6) in the overlap region, the celebrated logarithmic profile [23], known as the law of the wall. It is well supported by data in the range y+ & 100, y . 0.15R with the constants \u00CE\u00BA \u00E2\u0089\u0088 0.421 (the Von Ka\u00CC\u0080rma\u00CC\u0080n constant) and B \u00E2\u0089\u0088 5This requires that U\u00CF\u0084 be considered as \u00E2\u0088\u009A \u00CF\u0084w/\u00CF\u0081; that is, the velocity scale is set by the wall stress. 6 Introduction 5.6 [14]. For many applications, (6) is sufficient as a solution to the fluid problem. This is primarily because as \u00CE\u00B4v/R \u00E2\u0086\u0092 0, asymptotically all of the shear lies within the overlap region in which (6) is valid. One consequence of this is that integrating (6) from the edge of the boundary layer to the centre of the pipe gives an expression for mean velocity as a function of boundary layer thickness (or Reynolds number) in excellent agreement with experimental data [23]. This is true despite the fact that the arguments leading to (6) are invalid near the centre of the pipe. It is clear from the derivation that the assumption of a circular pipe was inessential. In fact, the universal form (6) is also observed in plane channels and slowly-growing flat plate turbulent boundary layers [23]. For non-circular pipes, sufficiently far from corners, a logarithmic profile for the wall-normal velocity is also observed, but the constants \u00CE\u00BA and B must be replaced with functions of a coordinate running parallel to the wall, and the mean wall stress should be replaced with its local value [1, 5, 20, 32]. Townsend [31] extended the Izakson-Millikan argument to the case of rough walls, by proposing that the outer flow should not depend on the roughness length scale ks. The derivation is then entirely analogous, except that B must be replaced by a function of the roughness Reynolds number k+s = ksU\u00CF\u0084/\u00CE\u00BDM . 6 The major fault of the above analysis is that no indication is given as to the accuracy of (6), or indeed why there should even be an overlap layer in which \u00CF\u0086 tends to a constant. To clarify these issues with the heuristic boundary-layer analysis initiated by Prandtl, in the 1950\u00E2\u0080\u0099s applied mathe- maticians developed the technique of matched asymptotic expansions [7]. In this, scaled, dimensionless equations are produced in the inner and outer regions, and a suitable small parameter is identified, usually a dimensionless boundary layer thickness, \u00CE\u00B4. Since the regions are characterized by different physical balances, these equations will necessarily differ in their dependence on \u00CE\u00B4. The existence of a boundary layer implies that \u00CE\u00B4 \u001C 1, so that solutions may be sought as asymptotic expansions in \u00CE\u00B4. In general, these solutions will contain undetermined constants or functions, which are sought through the process of matching, which ensures that the solution varies smoothly from the inner to the outer domain. In the language of matched asymptotic expansions, the law of the wall is the velocity profile in the matching re- gion; the Izakson-Millikan argument produces the leading order behaviour by ensuring that it is possible to match the inviscid outer solution with the 6\u00CE\u00BA is observed to be universal, i.e., independent of k+s . 7 Introduction geometry-independent inner solution. Several matching techniques exist. Introducing n, the dimensionless dis- tance from the wall, the inner region is characterized by n . \u00CE\u00B4, while the outer region is n \u00E2\u0088\u00BC 1. The most general form introduces an intermediate coordinate \u00CE\u00B6 such that \u00CE\u00B6/n\u00E2\u0086\u0092 0 and \u00CE\u00B4/\u00CE\u00B6 \u00E2\u0086\u0092 0 as \u00CE\u00B4 \u00E2\u0086\u0092 0. Then, one imposes the condition that, in \u00CE\u00B6 coordinates, the inner and outer solutions are equal as \u00CE\u00B4 \u00E2\u0086\u0092 0. A simpler way of matching is Van Dyke\u00E2\u0080\u0099s rule, which states, loosely, that in the limit \u00CE\u00B4 \u00E2\u0086\u0092 0, the inner expansion of the outer solution should equal the outer expansion of the inner solution. The matched inner and outer solutions may be used to obtain a com- posite expansion that is asymptotically valid in the entire flow domain. A simple prescription is to add the inner and outer solutions, and subtract their common part. As emphasized by Panton [22], the composite solution is the best representation of the solution at relatively large \u00CE\u00B4. This is par- ticularly important when applied to a turbulent boundary layer, since the oft-used law of the wall is merely the common part of inner and outer solu- tions, and can be quite inaccurate as a velocity profile when \u00CE\u00B4 is relatively large. Matched asymptotic expansions were first applied to the problem of a turbulent boundary layer by Yajnik [33] and Mellor [15]. These authors did not impose a turbulence closure, so that their solution is necessarily incom- plete. However, by assuming that the Reynolds stress admits an asymptotic expansion of a particular form, they were able to reproduce several empirical results, including the law of the wall. This had the important effect of clar- ifying (6) as an asymptotic result, valid only in the limit \u00CE\u00B4 \u00E2\u0086\u0092 0, and only in a limited range of n. When precise tests are made of the form of the veloc- ity profile, as in the Princeton Superpipe experiment [14], the distinction is crucial. For example, deviations from the logarithmic profile have been inter- preted by some researchers as evidence for additional \u00E2\u0080\u009Cmeso\u00E2\u0080\u009Dlayers, which, physically, should represent different physical balances. However, these cor- rections can be more economically interpreted as higher order terms in a standard asymptotic expansion [34]. Not only does this simplify the phys- ical interpretation, but it also reduces the number of degrees of freedom which need to be constrained by data, since the corrections can be obtained self-consistently from the governing equations. These early works did not obtain a complete solution to the velocity profile because no closure was imposed. However, they show that matched asymptotic expansions are the natural tool with which to solve boundary layer problems. In chapter 4, we will employ the technique to partially solve the difficult problem of a turbulent boundary layer in a non-circular pipe. 8 Chapter 1 Governing Equations As shown in Figure 1.1, we will consider an infinitely long pipe (or channel), aligned with x\u00CC\u0082, of general cross-section \u00E2\u0084\u00A6, with boundary \u00E2\u0088\u0082\u00E2\u0084\u00A6 of arbitrary shape. (See Table 1 for a glossary of symbols.) Figure 1.1: Sketch of geometry. The Reynolds-averaged Navier-Stokes equations for an incompressible fluid [23] are \u00E2\u0088\u0082ui \u00E2\u0088\u0082xi = 0, uj \u00E2\u0088\u0082ui \u00E2\u0088\u0082xj = \u00E2\u0088\u00921 \u00CF\u0081 \u00E2\u0088\u0082P \u00E2\u0088\u0082xi + \u00CE\u00BDM \u00E2\u0088\u00822ui \u00E2\u0088\u0082xj\u00E2\u0088\u0082xj \u00E2\u0088\u0092 \u00E2\u0088\u0082 \u00E2\u0088\u0082xj ( u\u00E2\u0080\u00B2iu \u00E2\u0080\u00B2 j ) . (1.1) 9 Chapter 1. Governing Equations In the above, and throughout, the Einstein summation convention is used, with indices running from 1 to 3 unless specified otherwise. Also, it will be convenient to use the dual notation u = (u1, u2, u3) = (u, v, w) and x = (x1, x2, x3) = (x, y, z). As boundary conditions we will assume a prescribed axial pressure gradient \u00E2\u0088\u0092 |\u00E2\u0088\u0086P/\u00E2\u0088\u0086x| and the no-slip condition at the wall. In the case of an open channel, we will assume that the free surface is flat, with normal z\u00CC\u0082. In this case the vertical component of gravity is assumed already absorbed into the pressure; the axial component gives rise to dP/dx. To simplify the discussion, we will only consider explicitly the case of a pipe. There is no loss of generality, since for any given channel, we can create an equivalent pipe by reflecting the channel shape about the free surface. The stress-free condition then reduces to regularity of the solution along the axis of symmetry. An important quantity not directly appearing in (1.1) is the wall stress \u00CF\u0084w, which must be defined as (\u00CF\u0081\u00CE\u00BDM\u00E2\u0088\u0082u/\u00E2\u0088\u0082y \u00E2\u0088\u0092 \u00CF\u0081u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2) \u00C2\u00B7 n\u00CC\u0082, in which y = (y, z), u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 \u00E2\u0089\u00A1 u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 y\u00CC\u0082+u\u00E2\u0080\u00B2w\u00E2\u0080\u00B2 z\u00CC\u0082, and n\u00CC\u0082 is a dimensionless unit vector normal to the wall. The inclusion of Reynolds stresses here is necessary because, unlike the case of a smooth wall, for a rough wall they need not vanish on the surface \u00E2\u0088\u0082\u00E2\u0084\u00A6. Near the wall, the velocity scale is set by \u00CF\u0084w. Its mean value \u00CF\u0084w is known a priori since by integrating (1.1) over \u00E2\u0084\u00A6, we see that |\u00E2\u0088\u0082\u00E2\u0084\u00A6| \u00CF\u0084w = \u00E2\u0088\u0092|\u00E2\u0084\u00A6|\u00E2\u0088\u0086P/\u00E2\u0088\u0086x. The appearance here of the hydraulic radius R = 2|\u00E2\u0084\u00A6|/|\u00E2\u0088\u0082\u00E2\u0084\u00A6| suggests that it is the relevant length scale for the outer region. In fact, there is another important outer region length scale, more nat- urally used to scale the tranverse coordinates: the maximum distance from the wall, Nmax. As scales for nondimensionalizing the equations, the differ- ence between R and Nmax is merely cosmetic, but later we will propose a constitutive relation for a turbulent mixing length. This should be, as much as possible, independent of geometry. This can only be so in units made dimensionless by Nmax.7 We will consider fully developed flow, in which the velocity is independent of x [23]. Initially, we will seek a solution with non-zero velocity only in x\u00CC\u0082, so that (u, v, w) = (u(y, z), 0, 0). This neglects the secondary flow, to which we will return in Chapter 4. Under these assumptions, one can show that momentum conservation in the transverse directions requires only that the pressure is constant through the cross-section. Since the continuity equation is satisfied identically, the only equation we need consider is conservation of 7For many shapes, such as squares, circles, and isosceles triangles, the two length scales are identical. The important exception is the infinite plane channel, for which the hydraulic radius is twice the maximum distance from the wall. 10 1.1. Scaling streamwise momentum. 1.1 Scaling Since the basic pipe flow results from the balance of the imposed pressure gradient with the Reynolds shear stresses, a natural choice for the velocity scale of turbulent fluctuations is UT = \u00E2\u0088\u009A R 2\u00CF\u0081 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0086P\u00E2\u0088\u0086x \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 = \u00E2\u0088\u009A\u00CF\u0084w/\u00CF\u0081 \u00E2\u0089\u00A1 U\u00CF\u0084 , the friction velocity. As a mean velocity scale, we will set U \u00E2\u0089\u00A1 U\u00CF\u0084/\u00CE\u00BA, where \u00CE\u00BA . 1 is treated as an empirical constant. For agreement with data, we will find that \u00CE\u00BA = 0.421, the Von Ka\u00CC\u0080rma\u00CC\u0080n constant.8 An appropriate scaling for the flow variables in the outer region is then u = Uu\u00E2\u0088\u0097, (u\u00E2\u0080\u00B2, v\u00E2\u0080\u00B2, w\u00E2\u0080\u00B2) = UT (u\u00E2\u0080\u00B2\u00E2\u0088\u0097, v \u00E2\u0080\u00B2 \u00E2\u0088\u0097, w \u00E2\u0080\u00B2 \u00E2\u0088\u0097), x = Lx\u00E2\u0088\u0097, (y, z) = Nmax(y\u00E2\u0088\u0097, z\u00E2\u0088\u0097), P = 1 2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0086P\u00E2\u0088\u0086x \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3LP\u00E2\u0088\u0097, \u00CF\u0084w = \u00CF\u0081U2T \u00CF\u0084w\u00E2\u0088\u0097, with L a streamwise length scale. Note that the fluctuating velocities have been scaled isotropically, in accordance with their assumed random nature. These scalings imply that the dimensionless mean wall stress is unity. Immediately dropping the asterisk decoration, the dimensionless stream- wise momentum equation is 0 = 2\u00CE\u00BE +\u00E2\u0088\u0087 \u00C2\u00B7 ( \u00E2\u0088\u0092u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 + 1 \u00CE\u00BARe\u00CF\u0084 \u00E2\u0088\u0087u ) , (1.2) where we have introduced the notation \u00E2\u0088\u0087 = ( \u00E2\u0088\u0082\u00E2\u0088\u0082y , \u00E2\u0088\u0082\u00E2\u0088\u0082z ) and u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 \u00E2\u0089\u00A1 u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 y\u00CC\u0082 + u\u00E2\u0080\u00B2w\u00E2\u0080\u00B2 z\u00CC\u0082. The dimensionless constants are Re\u00CF\u0084 = U\u00CF\u0084Nmax \u00CE\u00BDM , \u00CE\u00BA = UT U , \u00CE\u00BE = Nmax R . For all turbulent flows, Re\u00CF\u0084 \u001D 1/\u00CE\u00BA, and the viscous term may safely be ignored, in the outer region. However, as mentioned in the introduction, the 8Values are quoted in the range 0.37 to 0.43 [14]. It is clear from boundary layer con- siderations that we should use a value obtained at the largest possible Reynolds number. We will use the value determined in the Princeton Superpipe experiment at Re > 106, \u00CE\u00BA = 0.421.[14] 11 1.1. Scaling shear is concentrated in a thin boundary layer near the wall, within which variables need to be rescaled to reflect different physics. Since rough walled turbulent flows are affected by viscosity in addition to roughness, to choose the correct boundary layer rescaling we must first be precise about exactly which boundary layer is important. For any wall-bounded turbulent flow, there is a viscous sublayer, adja- cent to the wall, in which viscous stresses dominate the Reynolds stresses. However, for a rough wall, this layer may be irrelevant to determining mean flow properties if the roughness length scale ks is much bigger than the vis- cous length scale \u00CE\u00BDM/U\u00CF\u0084 . In this case, drag on the flow is caused by pressure differences between the fore and aft sides of roughness elements; viscous drag is negligible. We will consider this hydraulically rough (or fully rough) limit. Then we can ignore the viscous stress entirely, even in the boundary layer. In this case, (1.2) reduces to \u00E2\u0088\u0092 2\u00CE\u00BE = \u00E2\u0088\u0087 \u00C2\u00B7 (\u00E2\u0088\u0092u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2). (1.3) The modifications due to viscosity are considered in Appendix D. Evidently, the turbulence model will depend on the roughness length scale ks. We will use the simple scaling ks = Nmaxks\u00E2\u0088\u0097. In the fully rough limit, the dimensionless boundary layer thickness \u00CE\u00B4 will be proportional to ks: \u00CE\u00B4 = cks. It is convenient to choose9 c \u00E2\u0089\u0088 0.0279 so that, for a circular pipe, the law of the wall takes the simple form10 u \u00E2\u0088\u00BC log(n/\u00CE\u00B4). (1.4) We have derived (1.3) and justified neglect of the viscous term. To solve (1.3), we must now specify a turbulence closure. 90.0279 \u00E2\u0089\u0088 exp(\u00E2\u0088\u00928.5\u00CE\u00BA) 10This also requires that we choose \u00CE\u00BA \u00E2\u0089\u0088 0.421 12 Chapter 2 Turbulence Closure The simplest turbulence closure, due to Boussinesq [3], is the eddy viscosity model \u00E2\u0088\u0092 (u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2) = \u00CE\u00BDT\u00E2\u0088\u0087u, (2.1) where now the eddy viscosity \u00CE\u00BDT must be modelled. For a general flow, this constitutive assumption is not justified, and leads to erroneous results [23]. However, for convectively steady problems, such as we are considering here, with appropriate modelling of \u00CE\u00BDT , (2.1) is very successful [23]. Dimensionally, \u00CE\u00BDT is the product of a velocity scale with a length scale. Adopting the eddy phenomenology of Kolmogorov [11], and making an anal- ogy with the molecular viscosity of a gas, Prandtl [24] proposed that the length scale could be interpreted as a mixing length of dominant eddies, and that the velocity scale could be obtained as the product of the local velocity gradient with the mixing length, viz.,11 \u00CE\u00BDT = l2 |\u00E2\u0088\u0087u| , (2.2) where l is the mixing length. This implies the closure \u00E2\u0088\u0092 (u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2) = l2 |\u00E2\u0088\u0087u| \u00E2\u0088\u0087u. (2.3) To specify l, Prandtl invoked a dimensional argument in analogy to Izakson- Millikan. There are five length scales in the problem: R, Nmax, n, \u00CE\u00BDM/U\u00CF\u0084 , and ks. Prandtl proposed that there should be an overlap layer in which the only available length scale is the distance to the wall, so that l \u00E2\u0088\u009D n. It is easily seen that this prescription implies the law of the wall, so that the over- lap layer can be identified with the matching region of the boundary layer structure discussed in the introduction. Then, as in the Izakson-Millikan argument, in the fully turbulent outer layer, one would generally expect l to be geometry-dependent. Conversely, in the boundary layer, l could depend on viscous or roughness length scales. 11We have generalized from the plane channel considered by Prandtl to an arbitrary shape by replacing |du/dy| with |\u00E2\u0088\u0087u| 13 Chapter 2. Turbulence Closure In the symmetric plane channel case considered by Prandtl, l is a func- tion of distance to the wall, since this is the only available coordinate. For a non-circular pipe, there is the possibility that l could depend on two trans- verse coordinates. However, the mixing length is interpreted as the size of characteristic eddies \u00E2\u0080\u0093 rotating, circular vortices. Since the distance from the wall, n(x), is the radius of the largest circle, centred at x, that is con- tained within the pipe, the identification of n with eddy size is natural and entirely independent of the cross-sectional shape. We may therefore still expect l = l(n). Although heuristic, the mixing length approach has been very successful in relating experimental results to eddy phenomenology [31]. For example, using a sensible parametrization of l derived from DNS data, L\u00E2\u0080\u0099vov and others [13] have shown that the shear stress and mean velocity profiles of circular pipe and plane channel flow can be reproduced extremely accurately with few adjustable parameters. However, the closure (2.3) suffers from an important physical defect: as the centre of the pipe is approached, the Reynolds shear stress tends to zero linearly, while the velocity gradient also tends to zero linearly [22]. To satisfy (2.3), this forces the mixing length to diverge [13]. If l is truly a physically relevant quantity, this cannot be permitted. This can be alleviated by using an alternative velocity scale in the defi- nition of \u00CE\u00BDT . Kolmogorov [11] proposed that the appropriate velocity scale is provided by the turbulent kinetic energy, kT = u\u00E2\u0080\u00B22, so that \u00CE\u00BDT = \u00CE\u00BAkl \u00E2\u0088\u009A kT , with \u00CE\u00BAk a constant. To model kT , the turbulent energy equation is con- sidered under the approximation of local balance [13]. This implies that production of turbulent energy, \u00E2\u0088\u0092u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 \u00C2\u00B7\u00E2\u0088\u0087u, balances dissipation of turbulent energy, \u000FT . The latter is estimated as \u000FT = \u00CE\u00BA\u000Fk 3/2 T /l, with \u00CE\u00BA\u000F a constant; this may be taken as the definition of l. Defined in this way, l is well behaved throughout the entire pipe. Then \u00CE\u00BDT = \u00CE\u00BAkl \u00E2\u0088\u009A kT = \u00CE\u00BAkl ( l\u000FT \u00CE\u00BA\u000F )1/3 = \u00CE\u00BA\u00CE\u00BD(l4\u000FT )1/3, (2.4) with \u00CE\u00BA\u00CE\u00BD = \u00CE\u00BAk/\u00CE\u00BA 1/3 \u000F . Local balance of turbulent energy implies \u000FT = \u00E2\u0088\u0092u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 \u00C2\u00B7 \u00E2\u0088\u0087u = \u00CE\u00BDT |\u00E2\u0088\u0087u|2. (2.5) Solving (2.4) and (2.5), we find \u000FT = \u00CE\u00BA3/2\u00CE\u00BD l 2|\u00E2\u0088\u0087u|3, \u00CE\u00BDT = \u00CE\u00BA3/2\u00CE\u00BD l2|\u00E2\u0088\u0087u|. (2.6) 14 Chapter 2. Turbulence Closure To match with known behaviour in the logarithmic region [23], we must take \u00CE\u00BAk \u00E2\u0089\u0088 0.55, \u00CE\u00BA\u000F \u00E2\u0089\u0088 0.553, so that \u00CE\u00BA\u00CE\u00BD = 1. Therefore, the prescription for \u00CE\u00BDT is identical to the original proposal, (2.3). However, here we used local balance of turbulent energy, which is only strictly valid in the logarithmic region. Towards the centre of the pipe, the energy flux terms become important, so that our algebraic model for energy balance breaks down. The failure of the algebraic approximation leads to the divergence of l in (2.3). Fortunately, this has little impact on the mean velocity, which is nearly constant there. Hitherto, all parameterizations of l have been intended for smooth-walled flows, for which the no-slip condition requires that l \u00E2\u0086\u0092 0 as n \u00E2\u0086\u0092 0.[23] For a rough flow, we propose instead that l tends to a constant at the wall. Phenomenologically, eddies are prevented from decreasing in size indefinitely by the roughness asperities at the wall. Using the refined definition of l as \u00CE\u00BA\u000Fk 3/2 T /\u000FT , the mixing length is ob- served in DNS data to saturate to a constant in the outer region.12 We will model this behaviour. A simple form with the desired properties is l(n) \u00E2\u0089\u00A1 n+ c3ks 1 + c4n , (2.7) where ks is now dimensionless c4 & O(1) for the desired behaviour. With these constitutive relations, the streamwise momentum equation takes the form \u00E2\u0088\u0092 2\u00CE\u00BE = \u00E2\u0088\u0087 \u00C2\u00B7 (( n+ c3ks 1 + c4n )2 |\u00E2\u0088\u0087u| \u00E2\u0088\u0087u ) . (2.8) 12In direct numerical simulation, the Reynolds stress can be computed directly from the fluctuating velocity fields u\u00E2\u0080\u00B2 and v\u00E2\u0080\u00B2, so that the implied profile of mixing length can be found by rearranging (2.4). 15 Chapter 3 Circular Pipe Solution To calibrate the constants c3 and c4 with well-verified empirical relations, we now solve (1.3) in the case of a fully rough circular pipe, using matched asymptotic expansions. To maintain continuity with earlier arguments and with later analysis, as a coordinate we use the (dimensionless) distance from the wall, n, which is related to the usual (dimensionless) radial coordinate r by n = 1\u00E2\u0088\u0092 r. In this case, (1.3) with (2.3) and (2.7) becomes \u00E2\u0088\u0092 2 = 1 1\u00E2\u0088\u0092 n d dn ( (1\u00E2\u0088\u0092 n) ( n+ c3/c \u00CE\u00B4 1 + c4n )2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dudn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dudn ) . (3.1) Since variables have been scaled for the outer (fully turbulent) region, we na\u00CC\u0088\u00C4\u00B1vely pose an outer solution u = u0 + \u00CE\u00B4u1 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 . Inserting this into (3.1) and integrating once, we find the leading order equation \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3du0dn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 du0dn = (1\u00E2\u0088\u0092 2n+ n2 + C0)(1 + c4n)2n2(1\u00E2\u0088\u0092 n) . (3.2) At the centre of the pipe, n = 1 and du/dn = 0, so that the integration constant C0 = 0. Assuming du0/dn > 0 and integrating again, the outer region velocity is u0 = log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009A1\u00E2\u0088\u0092 n 1 + \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n ) + 2 \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n\u00E2\u0088\u0092 2 3 c4(1\u00E2\u0088\u0092 n)3/2 + C1, (3.3) where the constant C1 needs to be determined by matching with the bound- ary layer. To determine the inner (boundary-layer) equations, we need to rescale the variables n and u. Near the wall, the velocity scale is set by the wall stress, so that the appropriate scale is the friction velocity U\u00CF\u0084 . Therefore, an appropriate rescaling in the boundary layer is N = 1 \u00CE\u00B4 n, U = u, (3.4) 16 Chapter 3. Circular Pipe Solution leading to \u00E2\u0088\u0092 2\u00CE\u00B4 = 1 1\u00E2\u0088\u0092 \u00CE\u00B4N d dN ( (1\u00E2\u0088\u0092 \u00CE\u00B4N) ( N + c3/c 1 + c4\u00CE\u00B4N )2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dUdN \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dUdN ) . (3.5) Letting U = U0 + \u00CE\u00B4U1 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 as above, we find C2 = (N + c3/c) 2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dU0dN \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dU0dN . (3.6) At N = 0, the wall stress is unity in outer variables, so that C2 = 1. Assuming dU0/dN > 0 and integrating, we find the inner solution U0 = log ( cN c3 + 1 ) , (3.7) where, by the no-slip condition, there is no integration constant. Solutions (3.3) and (3.7) need to be matched. We will apply Van Dyke\u00E2\u0080\u0099s rule, which states that the outer expansion of the inner solution must equal the inner expansion of the outer solution. The former is u0m = log ( cn c3\u00CE\u00B4 ) +O(\u00CE\u00B4) (3.8) Expressed in outer variables, the inner expansion of the outer solution is u0m = log (n 4 ) + 2\u00E2\u0088\u0092 2 3 c4 + C1 +O(\u00CE\u00B4 log \u00CE\u00B4). (3.9) To obtain agreement to O(\u00CE\u00B4 log \u00CE\u00B4), we should take C1 = log(4c/c3\u00CE\u00B4)\u00E2\u0088\u00922+ 23c4. As mentioned in the introduction, the leading order velocity profile in the matching region is the law of the wall. Therefore, comparing (3.8) with (1.4), we see that c3 = c \u00E2\u0089\u0088 0.0279. This calibrates the wall-offset constant c3. The appearance of log(\u00CE\u00B4) in C1 indicates that our original expansion for u was not quite correct: written as an explicit asymptotic series in \u00CE\u00B4, we have, in the outer region, u0 = \u00E2\u0088\u0092 log(\u00CE\u00B4) + [ log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009A1\u00E2\u0088\u0092 n 1 + \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n ) + 2( \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n\u00E2\u0088\u0092 1) \u00E2\u0088\u00922 3 c4((1\u00E2\u0088\u0092 n)3/2 \u00E2\u0088\u0092 1) + log(4) ] . (3.10) 17 Chapter 3. Circular Pipe Solution The leading order term, \u00E2\u0088\u0092 log(\u00CE\u00B4), is asymptotically different from 1 as \u00CE\u00B4 \u00E2\u0086\u0092 0. Since we have only made asymptotic expansions in one dependent variable, there is no need to revise the expansion and match again. Note that the leading order contribution is produced entirely in the boundary layer. To calibrate the constant c4, which determines the saturation of l(n) for large n, we need to use a quantity sensitive to the velocity in the outer region. A suitable measure is the friction factor, a dimensionless measure of mean velocity[23]: f \u00E2\u0089\u00A1 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0086P\u00E2\u0088\u0086x \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 4R\u00CF\u0081U2\u00E3\u0080\u0088u0\u00E3\u0080\u00892\u00E2\u0084\u00A6 (3.11) = 8\u00CE\u00BA2 \u00E3\u0080\u0088u0\u00E3\u0080\u00892\u00E2\u0084\u00A6 . (3.12) Computation of \u00E3\u0080\u0088u0\u00E3\u0080\u0089\u00E2\u0084\u00A6 requires integrating the solution over the entire flow domain. For this purpose, a composite solution is invaluable. The additive composite solution resulting from (3.3) and (3.7) is u0c = log ( 1 + \u00CE\u00B4 n ) + log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009A1\u00E2\u0088\u0092 n 1 + \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n ) + 2 \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n\u00E2\u0088\u0092 2 3 c4(1\u00E2\u0088\u0092 n)3/2 + C1. (3.13) Using u0c in f , we find f = 8\u00CE\u00BA2 [ 1 pi \u00E2\u0088\u00AB 2pi 0 \u00E2\u0088\u00AB 1 0 (1\u00E2\u0088\u0092 n) u0c(n) dn d\u00CE\u00B8 ]\u00E2\u0088\u00922 (3.14) = 8\u00CE\u00BA2 [ \u00E2\u0088\u009216 15 \u00E2\u0088\u0092 8c4 21 + log ( 4 \u00CE\u00B4 ) \u00E2\u0088\u0092 2 + 2c4 3 +O(\u00CE\u00B4 log \u00CE\u00B4) ]\u00E2\u0088\u00922 (3.15) = [ 1 2 \u00E2\u0088\u009A 2\u00CE\u00BA ( log ( 4 cks ) \u00E2\u0088\u0092 46 15 + 6c4 21 +O(\u00CE\u00B4 log \u00CE\u00B4) )]\u00E2\u0088\u00922 . (3.16) The experimental result, due to Nikuradse [21] and later improved by Schlicht- ing [26], is f = [ 2.00 log10 ( 1 ks ) + 1.74 ]\u00E2\u0088\u00922 . (3.17) We can obtain agreement to O(\u00CE\u00B4 log \u00CE\u00B4) if we take13 c4 = 0.6. This completes calibration of the turbulence model. The resulting prescription for the mix- ing length is shown in Figure 3.1. Analogous results and computations for a plane channel are given in Appendix C. 13Note that log(10)/(2 \u00E2\u0088\u009A 2\u00CE\u00BA) \u00E2\u0089\u0088 1.93 for \u00CE\u00BA = 0.421. Better agreement can be obtained by using \u00CE\u00BA = 0.406, but it seems preferable to obtain agreement with the newest data; that is, with the Superpipe experiment that yields \u00CE\u00BA = 0.421 [14, 29]. 18 Chapter 3. Circular Pipe Solution 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 n l(n ) Figure 3.1: Mixing-length l(n) (solid) from (2.7) with c3 = 0.0279, c4 = 0.6, ks = 10\u00E2\u0088\u00923. For comparison, also shown (dashed) is n+ \u00CE\u00B4. Incidentally, with the simple mixing length model we have proposed, the exact solution to 3.1 can be obtained without the use of matched asymptotic expansions. The derivation of this solution is presented in Appendix A. The exact, inner, outer, composite, and law-of-the-wall solutions are compared in Figures 3.2 and 3.3, for \u00CE\u00B4 = 0.01 (top) and \u00CE\u00B4 = 0.001 (bottom). For \u00CE\u00B4 = 0.01, although the law of the wall is the indeed the outer asymptote of the inner solution and vice versa, it is a poor approximation to the velocity profile compared with the composite solution, which is within 0.1% of the exact solution everywhere. For \u00CE\u00B4 = 0.001, the law of the wall is a good approximation within a limited region, but fails for large and small values of n. The composite solution is within 0.001% of the exact solution everywhere. 19 Chapter 3. Circular Pipe Solution 10\u00E2\u0088\u00923 10\u00E2\u0088\u00922 10\u00E2\u0088\u00921 100 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 n u exact inner outer composite law of the wall 10\u00E2\u0088\u00920.834 10\u00E2\u0088\u00920.825 2.76 2.77 2.78 10\u00E2\u0088\u00923 10\u00E2\u0088\u00922 10\u00E2\u0088\u00921 100 0 1 2 3 4 5 6 7 n u exact inner outer composite law of the wall 10\u00E2\u0088\u00921.3381 10\u00E2\u0088\u00921.3372 3.853 3.8535 3.854 Figure 3.2: Exact, inner, outer, composite, and law-of-the-wall solutions for the mean velocity profile in a fully rough, circular pipe for \u00CE\u00B4 = 0.01 (top), and \u00CE\u00B4 = 0.001 (bottom), in semi-log coordinates. The exact solution is barely distinguishable from the composite solution (see inset). 20 Chapter 3. Circular Pipe Solution 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 n u exact inner outer composite law of the wall 0.594 0.596 0.598 4.068 4.07 4.072 4.074 4.076 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 n u exact inner outer composite law of the wall 0.518 0.52 6.235 6.24 6.245 Figure 3.3: Exact, inner, outer, composite, and law-of-the-wall solutions for the mean velocity profile in a fully rough, circular pipe for \u00CE\u00B4 = 0.01 (top), and \u00CE\u00B4 = 0.001 (bottom). The exact solution is barely distinguishable from the composite solution (see inset). 21 Chapter 4 An Approximation Scheme for Non-Circular Pipes Figure 4.1: Sketch of (n, \u00CE\u00B8) coordinates. For a non-circular pipe, the streamwise momentum equation \u00E2\u0088\u0092 2\u00CE\u00BE = \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087u) (4.1) is a partial differential equation (PDE), for which two independent coordi- nates are needed in the transverse direction (see Figure 4.1). To distinguish between wall-parallel and wall-normal directions, it is convenient to use n, the distance from the wall, supplemented with an azimuthal coordinate \u00CE\u00B8, perpendicular to n.14 The formulae needed to use these coordinates are developed in Appendix B. In these coordinates, (4.1) is \u00E2\u0088\u0092 2\u00CE\u00BE = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082n (\u00E2\u0088\u009A G \u00CE\u00BDT \u00E2\u0088\u0082u \u00E2\u0088\u0082n ) + 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B8 ( \u00CE\u00BDT\u00E2\u0088\u009A G \u00E2\u0088\u0082u \u00E2\u0088\u0082\u00CE\u00B8 ) , (4.2) in which G \u00E2\u0089\u00A1 |\u00E2\u0088\u0082v/\u00E2\u0088\u0082\u00CE\u00B8|2 characterizes the local geometry. For any \u00E2\u0084\u00A6, it may written as G(n, \u00CE\u00B8) = Gw(\u00CE\u00B8) (1\u00E2\u0088\u0092 n kw(\u00CE\u00B8))2 (4.3) 14\u00CE\u00B8 is normalized to the interval (0, 2pi). Although (n, \u00CE\u00B8) coordinates are valid in all of \u00E2\u0084\u00A6, the domain in the (n, \u00CE\u00B8) plane may be non-smooth, in which case extra regularity conditions are needed to obtain physical solutions. We will return to this point in section 4.1. 22 Chapter 4. An Approximation Scheme for Non-Circular Pipes in terms of its value at the wall, Gw, and the wall curvature, kw. The transformations that lead to equations (4.2) and (4.3) are given in Appendix B. For a circular pipe, kw = 1, while for a plane channel, kw = 0. Equation (4.2) is a difficult nonlinear partial differential equation. Even using matched asymptotic expansions, in general the outer problem must be solved numerically. We have done so with COMSOL, a finite element package. These numerical solutions are used to validate the approxima- tions developed below. Instead of relying entirely on numerical solutions, by examining the boundary layer structure of the problem, we can motivate approximate methods of solution. We will return to this after considering the boundary layer problem. To rescale for the inner region, we supplement (3.4) with V = \u00CE\u00BDT /\u00CE\u00B4. Then, (4.2) becomes \u00E2\u0088\u0092 2\u00CE\u00BE\u00CE\u00B4 = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082N (\u00E2\u0088\u009A GV \u00E2\u0088\u0082U \u00E2\u0088\u0082N ) + \u00CE\u00B42\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B8 ( V\u00E2\u0088\u009A G \u00E2\u0088\u0082U \u00E2\u0088\u0082\u00CE\u00B8 ) . (4.4) From this we see that, to O(\u00CE\u00B4), only the wall-normal turbulent flux is im- portant. Posing a solution U = U0 + \u00CE\u00B4U1 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , we find \u00CF\u0084w(\u00CE\u00B8) = V \u00E2\u0088\u0082U0 \u00E2\u0088\u0082N (4.5) with V = (N + 1)2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0082U0\u00E2\u0088\u0082N \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3+O(\u00CE\u00B4). The equation \u00CF\u0084w(\u00CE\u00B8) = (N + 1)2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0082U0\u00E2\u0088\u0082N \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00E2\u0088\u0082U0\u00E2\u0088\u0082N (4.6) is then the boundary layer equation analogous to (3.6).15 The important difference is that \u00CF\u0084w is now an unknown function of \u00CE\u00B8. Equation (4.6) can, of course, be integrated to U0 = \u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8) log(N + 1). (4.7) Recall that U was nondimensionalized with the friction velocity U\u00CF\u0084 = \u00E2\u0088\u009A \u00CF\u0084\u00CC\u0084w/\u00CF\u0081. From (4.7) we see that the appropriate velocity scale is actually the local fric- tion velocity \u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8)/\u00CF\u0081. Therefore, to leading order, azimuthal dependence 15Recall that C2 = 1, c3 = c. 23 Chapter 4. An Approximation Scheme for Non-Circular Pipes in the boundary layer is restricted to the variation of wall stress, which sets the appropriate velocity scale. The critical simplification that has allowed us to find a solution to the boundary layer problem is neglect of the term involving \u00E2\u0088\u0082U/\u00E2\u0088\u0082\u00CE\u00B8, which rep- resents the Reynolds shear stress along \u00CE\u00B8\u00CC\u0082. This simplification turns (4.2) from a partial differential equation into an ordinary differential equation, with parametric \u00CE\u00B8 dependence. This is justified by the small thickness of the boundary layer in comparison to relevant transverse length scales. In the outer layer, there is no such separation of length scales. However, there are several heuristic reasons to expect that the most important variation of u is along rays of constant \u00CE\u00B8, that is, \u00E2\u0088\u0087u \u00E2\u0089\u0088 \u00E2\u0088\u0082u \u00E2\u0088\u0082n \u00E2\u0088\u0087n, (4.8) a \u00E2\u0080\u009Cradial\u00E2\u0080\u009D approximation. First, the streamwise velocity is empirically known to increase mono- tonically from the wall to the centre of the pipe, where the velocity is, of course, independent of direction. This geometrical constraint suggests that azimuthal variation is maximal near the wall, where (4.8) is asymptotically valid. Second, we have already motivated a mixing length turbulence model in which the mixing length l depends only on n, the distance from the wall. Physically, this corresponds to the identification of l as the size of dominant eddies. In this model, azimuthal dependence is expected to derive from geometrical quantities that vary between rays of constant \u00CE\u00B8; namely, wall curvature, kw(\u00CE\u00B8), and the maximum distance from the wall, nmax(\u00CE\u00B8). For cross-sectional shapes with these quantities weakly varying, (4.8) may be a good approximation. It is important to recall that we have neglected the secondary flow. The secondary flow adds important \u00E2\u0080\u009Cnon-radial\u00E2\u0080\u009D physics that raises serious doubt about the validity of (4.8). We will comment further on this below. 24 Chapter 4. An Approximation Scheme for Non-Circular Pipes Applying (4.8) to (4.1), and using (B.14) from Appendix B, we have \u00E2\u0088\u00922\u00CE\u00BE = \u00E2\u0088\u0087 \u00C2\u00B7 ( \u00CE\u00BDT \u00E2\u0088\u0082u \u00E2\u0088\u0082n \u00E2\u0088\u0087n ) = \u00E2\u0088\u0087 \u00C2\u00B7 ( \u00CE\u00BDT \u00E2\u0088\u0082u \u00E2\u0088\u0082n e\u00CC\u0082n ) = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082n (\u00E2\u0088\u009A G\u00CE\u00BDT \u00E2\u0088\u0082u \u00E2\u0088\u0082n ) = 1 1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n \u00E2\u0088\u0082 \u00E2\u0088\u0082n ( (1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n)l(n)2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0082u\u00E2\u0088\u0082n \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00E2\u0088\u0082u\u00E2\u0088\u0082n ) . (4.9) It is easy to see that we could have obtained this equation directly from (4.11) by neglecting the turbulent shear stress along e\u00CC\u0082\u00CE\u00B8. Integrating this equation, \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0082u\u00E2\u0088\u0082n \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00E2\u0088\u0082u\u00E2\u0088\u0082n = \u00CF\u0084(n, \u00CE\u00B8)1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n ( 1 + c4n n )2 , (4.10) where we have introduced the reduced stress \u00CF\u0084(n, \u00CE\u00B8) = \u00CF\u0084w(\u00CE\u00B8)\u00E2\u0088\u0092 2\u00CE\u00BEn+ \u00CE\u00BEkw(\u00CE\u00B8)n2. In exceptional cases, we may have \u00E2\u0088\u0082u/\u00E2\u0088\u0082n < 0. However, the approximations we will later make require \u00E2\u0088\u0082u/\u00E2\u0088\u0082n \u00E2\u0089\u00A5 0. Therefore, u(n, \u00CE\u00B8) = C + \u00E2\u0088\u00AB n 0 \u00E2\u0088\u009A \u00CF\u0084(n\u00E2\u0080\u00B2, \u00CE\u00B8)\u00E2\u0088\u009A 1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n\u00E2\u0080\u00B2 ( 1 + c4n\u00E2\u0080\u00B2 n\u00E2\u0080\u00B2 ) dn\u00E2\u0080\u00B2. (4.11) The constant C is determined by matching to the boundary layer solution (4.7). Assuming the validity of this functional form for u, the remaining \u00E2\u0080\u009Cbound- ary condition\u00E2\u0080\u009D that determines the unknown function \u00CF\u0084w(\u00CE\u00B8) is regularity of the solution in the outer region. In particular, u and \u00E2\u0088\u0087u should be smooth throughout the outer region. This is guaranteed everywhere except along the curve n = nmax(\u00CE\u00B8), where the coordinates are discontinuous. (see Figure 4.2) Along this curve, we must therefore prescribe the jump conditions [u]+\u00E2\u0088\u0092 = 0, [ \u00E2\u0088\u0082u \u00E2\u0088\u0082y ]+ \u00E2\u0088\u0092 = 0, [ \u00E2\u0088\u0082u \u00E2\u0088\u0082z ]+ \u00E2\u0088\u0092 = 0. (4.12) The three conditions are not independent, since [u]+\u00E2\u0088\u0092 = 0 implies that the tangential component of \u00E2\u0088\u0087u is continuous along the curve. Therefore, we 25 Chapter 4. An Approximation Scheme for Non-Circular Pipes only need to impose differentiability of u normal to the curve. In fact, the system formed by equations (4.11) and (4.12) form an integro-differential equation that is not amenable to solution without further approximation. To make progress, we will consider two alternative methods of approximation. The first is to approximate the form of u, and apply the exact boundary condition (4.14), while the second is to directly approximate an expression for the reduced stress \u00CF\u0084 . Because these are uncontrolled approximations, they must be duly tested with any given geometry. To focus the discussion, it is helpful to restrict attention to a family of geometries collectively known as anisotropic superellipses. These are shapes parameterized by\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3y0 a \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3my + \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3z0 b \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3mz = 1 \u00E2\u0087\u0094 { y0 = a | cos(\u00CE\u00B8)|2/my sgn(cos(\u00CE\u00B8)) z0 = b | sin(\u00CE\u00B8)|2/mz sgn(sin(\u00CE\u00B8)) . (4.13) We will let 0 < my,mz \u00E2\u0089\u00A4 2. As special cases, this family includes ellipses (my = mz = 2) and diamonds (my = mz = 1). A plane channel is re- produced in the singular limit mz = 2, b/a \u00E2\u0086\u0092 0. We will let the shape either be \u00E2\u0080\u009Cfull\u00E2\u0080\u009D, parameterized by \u00CE\u00B8 \u00E2\u0088\u0088 (0, 2pi), or \u00E2\u0080\u009Csemi\u00E2\u0080\u009D, where \u00CE\u00B8 \u00E2\u0088\u0088 (0, pi), and the bottom of the shape is the line z = 0. (see Figure 4.2) The former covers a wide variety of pipe shapes, and, through the symmetry argument mentioned earlier, a variety of realistic open-channel shapes. The \u00E2\u0080\u009Csemi\u00E2\u0080\u009D shapes are reasonable pipe shapes for many geophysical applications [18].16 The formulae for y(n, \u00CE\u00B8), z(n, \u00CE\u00B8), Gw(\u00CE\u00B8), kw(\u00CE\u00B8), and nmax(\u00CE\u00B8) can all be analytically derived for anisotropic superellipses. These results are presented in Appendix B. We will apply the first approximation, which uses (4.12), only to shapes that are symmetric about the curve n = nmax(\u00CE\u00B8). All \u00E2\u0080\u009Cfull\u00E2\u0080\u009D anisotropic superellipses are of this type. In this case, the jump conditions (4.12) sim- plify. Continuity of u is guaranteed, and instead of requiring no jump in the u derivative normal to the curve n = nmax(\u00CE\u00B8), we have instead that the derivative vanishes. Therefore, the jump conditions become local boundary conditions. Restricting our attention to one side of the curve, we may con- sider n = nmax(\u00CE\u00B8) as the level set of n\u00E2\u0088\u0092nmax(\u00CE\u00B8). Then, setting the derivative 16We have normalized by the maximum distance from the wall. Therefore, for a \u00E2\u0080\u009Cfull\u00E2\u0080\u009D shape, 1 = min(a, b). 26 Chapter 4. An Approximation Scheme for Non-Circular Pipes normal to the curve to zero, we have17 0 = \u00E2\u0088\u0087(n\u00E2\u0088\u0092 nmax(\u00CE\u00B8)) \u00C2\u00B7 \u00E2\u0088\u0087u = \u00E2\u0088\u0082u \u00E2\u0088\u0082n \u00E2\u0088\u0092 1 G dnmax(\u00CE\u00B8) d\u00CE\u00B8 \u00E2\u0088\u0082u \u00E2\u0088\u0082\u00CE\u00B8 . (4.14) 17Here we use the fact that |\u00E2\u0088\u0087n| = 1, |\u00E2\u0088\u0087\u00CE\u00B8| = 1/\u00E2\u0088\u009AG; see Appendix B. 27 Chapter 4. An Approximation Scheme for Non-Circular Pipes y z \u00E2\u0088\u00924 \u00E2\u0088\u00923 \u00E2\u0088\u00922 \u00E2\u0088\u00921 0 1 2 3 4 0 0.5 1 1.5 2 y z \u00E2\u0088\u00924 \u00E2\u0088\u00923 \u00E2\u0088\u00922 \u00E2\u0088\u00921 0 1 2 3 4 0 0.5 1 1.5 2 y z \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 1.5 y z \u00E2\u0088\u00921.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 1.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.8 \u00E2\u0088\u00920.6 \u00E2\u0088\u00920.4 \u00E2\u0088\u00920.2 0 0.2 0.4 0.6 0.8 1 Figure 4.2: Examples of anisotropic superellipses. Contours of n (black) and \u00CE\u00B8 (red) are shown, along with the surface n = nmax(\u00CE\u00B8) (magenta). (top) half superellipse with my = 1,mz = 2. (middle) half superellipse with my = mz = 1. (bottom left) square with my = mz = 1, a/b = 1. (bottom right) ellipse with my = mz = 2, a/b = 1.5. 28 4.1. Constant-stress Approximation 4.1 Constant-stress Approximation For circular pipes at small boundary layer thicknesses, the logarithmic pro- file is a reasonable approximation to the velocity profile throughout the outer region, although the arguments leading to it are invalid there (see Figure 3.2). Indeed, it is used in many engineering applications when more sophis- ticated analysis is unwarranted. Therefore, a simple ansatz for the velocity profile in a non-circular pipe is to extend the logarithmic profile into the outer region, viz., u0 = \u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8) log (n \u00CE\u00B4 ) . (4.15) It is simple to check that this results from (4.11) under neglect of wall curvature, mixing length saturation, and pressure gradient. Whereas in a circular pipe the logarithmic profile is associated with constant stress [31], in a non-circular pipe, (4.15) is associated with stress constant along each ray of constant \u00CE\u00B8. Substituting (4.15) into the boundary condition (4.14), after some re-arranging, leads to d d\u00CE\u00B8 log(\u00CF\u0084w) = 2Gw(1\u00E2\u0088\u0092 kwnmax)2 ( nmax dnmax d\u00CE\u00B8 log (nmax \u00CE\u00B4 ))\u00E2\u0088\u00921 , (4.16) and therefore determination of \u00CF\u0084w is reduced to quadrature. The constant of integration is set by the requirement of unit mean wall stress. A well-known shortcoming of the logarithmic profile is that it has a non- zero slope at the centre of the pipe, which is forbidden by regularity.18 It is easy to see from (4.14) that this may lead to problems in a non-circular pipe. In particular, if dnmax/d\u00CE\u00B8 = 0 somewhere, then 1/G du0/d\u00CE\u00B8 is required to diverge. This severely restricts the applicability of the logarithmic profile approximation. However, for shapes where dnmax/d\u00CE\u00B8 6= 0 everywhere, (4.16) is often a good approximation. To illustrate, we will derive an explicit expression for \u00CF\u0084w for a square pipe, and compare this with the DNS data of Huser and Biringen.[8] 18For circular pipes, this is one reason to prefer a solution of the type presented in Chapter 3. 29 4.1. Constant-stress Approximation 4.1.1 Square Pipe By symmetry, we can restrict our discussion to the region \u00CE\u00B8 \u00E2\u0088\u0088 (0, pi/2). For a square pipe19, kw = 0 and Gw(\u00CE\u00B8) = 4 sin2(2\u00CE\u00B8). Then for \u00CE\u00B8 < pi/4, n = nmax(\u00CE\u00B8)\u00E2\u0087\u0094 z = 0 \u00E2\u0087\u0094 nmax(\u00CE\u00B8) = 2 sin2(\u00CE\u00B8), (4.17) using (B.22) from Appendix B. Similarly, for \u00CE\u00B8 \u00E2\u0088\u0088 (pi/4, pi/2), we have nmax(\u00CE\u00B8) = 2 cos2(\u00CE\u00B8). Substituting these expressions into (4.16) and integrating, we find \u00CF\u0084w(\u00CE\u00B8) = \u00CF\u00840 [ log ( nmax(\u00CE\u00B8) \u00CE\u00B4 )]2 , (4.18) where the constant \u00CF\u00840 is set by the requirement that \u00E3\u0080\u0088\u00CF\u0084w\u00E3\u0080\u0089\u00E2\u0088\u0082\u00E2\u0084\u00A6 = 1. There is a minor problem with this expression: for nmax < \u00CE\u00B4 the result is an unphysical divergence of the wall stress. This indicates that there is a corner boundary layer. Since the approximations made do not warrant a rigorous treatment, a simple fix, inspired by (4.7), is to replace nmax/\u00CE\u00B4 by nmax/\u00CE\u00B4 + 1. The approximate solution for u0 is then u0(n, \u00CE\u00B8) = \u00E2\u0088\u009A \u00CF\u00840 log ( nmax(\u00CE\u00B8) \u00CE\u00B4 + 1 ) log (n \u00CE\u00B4 + 1 ) . (4.19) In Figure 4.3, the expression (4.18) is compared to the solution of the full partial differential equation (1.3), the direct numerical simulation (DNS) data of Huser and Biringen[8]20, and the zero-stress approximation (next section). Only one eighth of the perimeter is shown, from a corner (left) to the centre-line (right). The units are arbitrary. It is apparent that (4.15) is an excellent approximation to the full PDE, equation (2.8), justifying the constant stress approximation. However, outside of the corner region, s . 0.05, agreement with data is only marginal. This indicates that either the turbulence model is inaccurate, or that our neglect of the secondary flow was unjustified. Since simple mixing-length turbulence models have been shown to reproduce data extremely well, as in Chapter 3, we seek an explanation in terms of the secondary flow. In particular, we may suppose 19This derivation can easily be generalized to an arbitrary diamond shape (a/b > 1 in 4.13). 20The DNS of Huser and Biringen was computed for a smooth wall at Re\u00CF\u0084 = 300, corresponding to \u00CE\u00B4 = 3 10\u00E2\u0088\u00924. Since the only impact of the wall on the outer flow should be through \u00CE\u00B4, we can compare the rough wall theory to data by considering a fully rough pipe at the same value of \u00CE\u00B4. 30 4.1. Constant-stress Approximation that the velocity (u, 0, 0) is only the leading order term in a more complete solution that includes the secondary flow. The latter is induced because wall stress variations are unstable to isotropizing perturbations [31]. Therefore, the secondary flow is in such a direction as to minimize wall stress variation. This is corroborated by the data presented in Figures 4.4 and 4.3. Since the wall stress must be zero at the corner and have unit mean (shown as a dashed line), one expects that as \u00CE\u00B4 \u00E2\u0086\u0092 0, \u00CF\u0084w will increase its slope at the corner. This is predicted by equation (4.18). Also, consider the intersection of the DNS data with the theoretical curves. To the left of this point, the secondary flow has increased the wall stress, so u should increase and the flow should be converging toward the wall. This is supported by the DNS data. (see Figure 4.4) Note, in particular, that the intersection point is aligned with the centre of the lower secondary flow vortex. 31 4.1. Constant-stress Approximation 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 s \u00CF\u0084 w (s) zero stress const. stress DNS PDE Figure 4.3: (top) \u00CF\u0084w for a square pipe, in terms of arclength s; (bottom) Secondary flow velocity field from DNS of Huser and Biringen.[8] 32 4.1. Constant-stress Approximation Figure 4.4: Streamwise velocity in a square pipe. (left) Constant-stress approximate solution, (right) DNS of Huser and Biringen.[8] 33 4.2. Zero-stress Approximation 4.2 Zero-stress Approximation In a circular pipe and plane channel, the shear stress decreases linearly from the wall to the centre-line, where it vanishes. A simple approximation for a general \u00E2\u0084\u00A6, assuming (4.8), is therefore to let the analogous quantity, the reduced stress \u00CF\u0084(n, \u00CE\u00B8), vanish on n = nmax(\u00CE\u00B8). It may be seen from (4.10) that this implies \u00E2\u0088\u0082u/\u00E2\u0088\u0082n = 0 on n = nmax(\u00CE\u00B8), which may or may not be a good approximation to (4.14), depending on \u00E2\u0084\u00A6. This gives a simple expression for the wall stress, \u00CF\u0084w(\u00CE\u00B8) = \u00CE\u00BE nmax(\u00CE\u00B8) (2\u00E2\u0088\u0092 kw(\u00CE\u00B8)nmax(\u00CE\u00B8)) , (4.20) which should be contrasted with (4.18). In particular, note that, here, \u00CF\u0084w has no \u00CE\u00B4 dependence. Also, there are no adjustable constants. Fortunately, |\u00E2\u0088\u0082\u00E2\u0084\u00A6|\u00E3\u0080\u0088\u00CF\u0084w(\u00CE\u00B8)\u00E3\u0080\u0089\u00E2\u0088\u0082\u00E2\u0084\u00A6 = \u00E2\u0088\u00AB 2pi 0 ( 2\u00CE\u00BE nmax(\u00CE\u00B8)\u00E2\u0088\u0092 \u00CE\u00BEkw(\u00CE\u00B8)n2max(\u00CE\u00B8) ) ds = 2\u00CE\u00BE \u00E2\u0088\u00AB 2pi 0 ( nmax(\u00CE\u00B8)\u00E2\u0088\u0092 12kw(\u00CE\u00B8)n 2 max(\u00CE\u00B8) ) \u00E2\u0088\u009A Gw(\u00CE\u00B8) d\u00CE\u00B8 = 2\u00CE\u00BE \u00E2\u0088\u00AB 2pi 0 \u00E2\u0088\u00AB nmax(\u00CE\u00B8) 0 (1\u00E2\u0088\u0092 nkw(\u00CE\u00B8)) \u00E2\u0088\u009A Gw(\u00CE\u00B8) dn d\u00CE\u00B8 = 2\u00CE\u00BE \u00E2\u0088\u00AB 2pi 0 \u00E2\u0088\u00AB nmax(\u00CE\u00B8) 0 \u00E2\u0088\u009A G(n, \u00CE\u00B8) dn d\u00CE\u00B8 = 2\u00CE\u00BE|\u00E2\u0084\u00A6|. In dimensionless terms, \u00CE\u00BE = |\u00E2\u0088\u0082\u00E2\u0084\u00A6|/(2|\u00E2\u0084\u00A6|), so that \u00E3\u0080\u0088\u00CF\u0084w(\u00CE\u00B8)\u00E3\u0080\u0089\u00E2\u0088\u0082\u00E2\u0084\u00A6 = 1, as desired. This identity is an indication that (4.20) is a sensible approximation. To find the velocity profile, we will ignore the minor effect of wall cur- vature, noting that the same velocity profile results for a circular pipe and plane channel (see Appendix C). Then the integral (4.11) can be computed analytically: u(n, \u00CE\u00B8) = C + \u00E2\u0088\u00AB n 0 \u00E2\u0088\u009A \u00CF\u0084(n\u00E2\u0080\u00B2, \u00CE\u00B8)(1 + c4n\u00E2\u0080\u00B2) n\u00E2\u0080\u00B2 dn\u00E2\u0080\u00B2 = C + 2 \u00E2\u0088\u009A \u00CF\u0084(n, \u00CE\u00B8) ( 1\u00E2\u0088\u0092 c4\u00CF\u0084(n, \u00CE\u00B8) 3\u00CE\u00BE ) \u00E2\u0088\u0092\u00E2\u0088\u009A\u00CF\u0084w log (\u00E2\u0088\u009A \u00CF\u0084w + \u00E2\u0088\u009A \u00CF\u0084(n, \u00CE\u00B8) \u00E2\u0088\u009A \u00CF\u0084w \u00E2\u0088\u0092 \u00E2\u0088\u009A \u00CF\u0084(n, \u00CE\u00B8) ) . (4.21) This matches to the inner solution when C = \u00E2\u0088\u00922\u00E2\u0088\u009A\u00CF\u0084w ( 1\u00E2\u0088\u0092 c4\u00CF\u0084w 3\u00CE\u00BE ) \u00E2\u0088\u0092\u00E2\u0088\u009A\u00CF\u0084w log ( 2\u00CF\u0084w \u00CE\u00B4 ) . 34 4.2. Zero-stress Approximation To illustrate, we consider an elliptical pipe. 4.2.1 Elliptical Pipe For an elliptical pipe with a > b = 1, kw(\u00CE\u00B8) = 8a \u00CF\u0086(\u00CE\u00B8)3 , Gw(\u00CE\u00B8) = \u00CF\u0086(\u00CE\u00B8)2 4 , nmax(\u00CE\u00B8) = \u00CF\u0086(\u00CE\u00B8) 2a , (4.22) in terms of \u00CF\u0086(\u00CE\u00B8) = \u00E2\u0088\u009A 4a2 sin2(\u00CE\u00B8) + 4 cos2(\u00CE\u00B8). Then, after some manipulation, \u00CF\u0084w(\u00CE\u00B8) = \u00CE\u00BE a (1 + 2(a2 \u00E2\u0088\u0092 1) sin2(\u00CE\u00B8))\u00E2\u0088\u009A a2 sin2(\u00CE\u00B8) + cos2(\u00CE\u00B8) . (4.23) In Figure 4.2.1, this approximation is compared to the solution of the full PDE (2.8), and to the DNS of Nikitin and Yakhot.[20]21 Although the qualitative behaviour of the approximation is correct, it overestimates the anisotropy of \u00CF\u0084w by more than a factor of two compared to the PDE. The PDE is a reasonable approximation to the DNS, but it too overestimates the anisotropy of \u00CF\u0084w. As with the square pipe, this can be explained in terms of the secondary flow, which is directed towards the elongated end of the ellipse. Owing to the absence of corners, here the secondary flow is less pro- nounced. The velocity distributions from the PDE and DNS are compared in Figure 4.6. 21The DNS was performed for a smooth wall at Re\u00CF\u0084 = 157, corresponding to \u00CE\u00B4 = 6\u00C3\u0097 10\u00E2\u0088\u00924. Again, we consider the equivalent fully rough wall. 35 4.2. Zero-stress Approximation 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 1.2 s \u00CF\u0084 w (s) zero stress DNS PDE Figure 4.5: (top) \u00CF\u0084w for an elliptical pipe of aspect ratio 2, in terms of arclength; (bottom) Contours of secondary flow velocity field from DNS of Nikitin and Yakhot.[20] Arrows indicate the direction of secondary flow. 36 4.2. Zero-stress Approximation Figure 4.6: Streamwise velocity in an elliptical pipe. (top) Solution of full PDE, (bottom) DNS of Nikitin and Yakhot.[20] 37 4.2. Zero-stress Approximation 4.2.2 Semi-Elliptical Pipe Semi-superellipses with my = 2 (see top of Figure 4.2) have dnmax/d\u00CE\u00B8 = 0 at the centre of the pipe, and therefore will not admit the constant-stress approximation. Fortunately, for these shapes, the zero-stress approximation is very accurate. For general mz and b = 2, we have kw(\u00CE\u00B8) = 4am2z \u00CF\u0086(\u00CE\u00B8)3 | sin(\u00CE\u00B8)|2\u00E2\u0088\u00924/mz ( 2\u00E2\u0088\u0092 2 mz cos2(\u00CE\u00B8)\u00E2\u0088\u0092 sin2(\u00CE\u00B8) ) , nmax(\u00CE\u00B8) = 2| sin(\u00CE\u00B8)|2/mz 1 + amz \u00CF\u0086(\u00CE\u00B8) | sin(\u00CE\u00B8)|2\u00E2\u0088\u00922/mz (4.24) in terms of \u00CF\u0086(\u00CE\u00B8) = \u00E2\u0088\u009A a2m2z| sin(\u00CE\u00B8)|4\u00E2\u0088\u00924/mz + 16 cos2(\u00CE\u00B8). The implied formula for \u00CF\u0084w = \u00CE\u00BEnmax(2\u00E2\u0088\u0092 kwnmax) does not simplify very much. For mz = 1, it is compared in Figure 4.7 to the solution of the PDE (2.8). Agreement is quite good. There are no DNS results for this shape, but we can imagine that the secondary flow would reduce the anisotropy in \u00CF\u0084w, in a similar manner as in the square pipe. 2 4 6 8 10 12 14 16 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 s \u00CF\u0084 w (s) zero stress PDE Figure 4.7: \u00CF\u0084w for a semi-elliptical pipe of aspect ratio 2, with linear cor- ners, in terms of arclength. The left half corresponds to the top half of the boundary. 38 Chapter 5 The Thermal Problem For many geophysical problems involving turbulent flow, in addition to the velocity distribution it is necessary to obtain the temperature distribution in the flow and heat fluxes into the wall. As an application of the model, so- lutions and methods of the previous chapters, we will pose the heat transfer problem in a non-circular pipe and develop an approximation scheme, with an eye toward applications in glaciology. In particular, floods under, through, or adjacent to glaciers involve the flow of turbulent water through a conduit or channel of ice [4]. Owing to melting at the conduit wall, the conduit shape may evolve down the flow path. If the length scale along which the conduit evolves is much longer than its hydraulic radius, then the infinite pipe assumption considered in this thesis may be appropriate. In this case, one would like to know the heat flux at the wall, so that the melt rate can be determined, and, therefore, the evolution of the conduit shape. The situation may be novel from a physical point of view because, in the case of a subglacial flood, the water is only infinitesimally warmer than the surrounding ice. This affects the relevant physics, as explained below. In many ways, the thermal problem is analogous to the mechanical prob- lem. Molecular diffusion of heat, which may be neglected in the core of the turbulent flow, is important at the wall, leading to a thermal boundary layer. Since the thermal problem depends on the velocity field, this is com- plicated by the presence of the viscous boundary layer. For a smooth wall, the boundary layer thickness depends on the thermal diffusivity of the fluid; when the Prandtl number Pr = \u00CE\u00B1/\u00CE\u00BDM 22 is greater than unity, the thermal boundary layer is located within the viscous boundary layer, and universal behaviour is observed [27]. For a rough wall, there are three boundary layers. In the fully rough regime, this actually simplifies the problem, because the most important properties of the velocity and temperature distributions are determined by the same roughness boundary layer\u00E2\u0080\u0093 at least, in the context of the simplest 22\u00CE\u00B1 is the thermal diffusivity of the fluid. 39 5.1. Scaling turbulent heat flux models. We will consider this situation. As before, we decompose the temperature into mean and fluctuating components, viz., T = T\u00CC\u0084 + T \u00E2\u0080\u00B2. (5.1) We perform the Reynolds-averaging procedure and drop the bar decora- tion immediately. The resulting Reynolds-averaged energy equation for an incompressible fluid is cpuj \u00E2\u0088\u0082T \u00E2\u0088\u0082xj = \u000F+ \u000FT + cp\u00CE\u00B1 \u00E2\u0088\u00822T \u00E2\u0088\u0082xj\u00E2\u0088\u0082xj \u00E2\u0088\u0092 cp \u00E2\u0088\u0082 \u00E2\u0088\u0082xj ( u\u00E2\u0080\u00B2jT \u00E2\u0080\u00B2 ) . (5.2) In these equations, \u000F and \u000FT are the viscous energy and turbulent energy dissipations, given by \u000F = 1 2 \u00CE\u00BDM ( \u00E2\u0088\u0082ui \u00E2\u0088\u0082xj + \u00E2\u0088\u0082uj \u00E2\u0088\u0082xi )2 , \u000FT = 1 2 \u00CE\u00BDM ( \u00E2\u0088\u0082u\u00E2\u0080\u00B2i \u00E2\u0088\u0082xj + \u00E2\u0088\u0082u\u00E2\u0080\u00B2j \u00E2\u0088\u0082xi )2 . (5.3) The dissipation terms are often neglected, but they are important in glacio- logical applications; we will, therefore, retain them. As a model for \u000FT , we will use the simple model23 discussed in Chapter 2, \u000FT = l2|\u00E2\u0088\u0087u|3. (5.4) We are interested in the heat flux at the wall, which must be defined as qw \u00E2\u0089\u00A1 \u00E2\u0088\u0092\u00CF\u0081cw\u00CE\u00B1\u00E2\u0088\u0082T \u00E2\u0088\u0082n + \u00CF\u0081cwu\u00E2\u0080\u00B2nT \u00E2\u0080\u00B2, (5.5) in terms of the inward facing normal n. 5.1 Scaling We will introduce separate scales for the mean flow and turbulent fluctua- tion temperatures. The former is measured from a datum T0, which we will assume is the temperature at the pipe boundary. To determine the scale 23It should be noted that this model for \u000FT can also be motivated by considering the Reynolds-averaged Navier Stokes equations as modelling a complex fluid. Then, using methods of continuum mechanics, one may derive the implied energy equation. The dissipation term that results is \u000FT . 40 5.1. Scaling for turbulent temperature fluctuations, \u00E2\u0088\u0086TT , we will assume that the tur- bulent energy flux balances the turbulent dissipation. The former scales as cpUT\u00E2\u0088\u0086TT /Nmax, while the latter scales as24 (\u00CE\u00BANmax)2(U/Nmax)3. Thus \u00E2\u0088\u0086TT = Nmax cpUT ( (\u00CE\u00BANmax)2 ( U Nmax )3) = \u00CE\u00BAU2 cp , (5.6) where we have used \u00CE\u00BA = UT /U . This is an unusual situation, not often considered in the literature. It is relevant for the glaciological application because the water temperature in subglacial conduits does not deviate much from 0oC. For example, consider a subglacial flood with water running at 10m/s. Then the above scaling predicts \u00E2\u0088\u0086TT \u00E2\u0088\u00BC 0.01K. The mean flow temperature scale, \u00E2\u0088\u0086T , is assumed proportional to the turbulent scale, viz., \u00E2\u0088\u0086T = (PrT /\u00CE\u00BA)\u00E2\u0088\u0086TT . The turbulent Prandtl number, PrT , is an empirical constant, usually taken in the range 0.7-0.9 [27]. Scaling the new variables T = T0 + \u00E2\u0088\u0086T T\u00E2\u0088\u0097, T \u00E2\u0080\u00B2 = \u00E2\u0088\u0086TT T \u00E2\u0080\u00B2\u00E2\u0088\u0097, qw = \u00CF\u0081cpUT\u00E2\u0088\u0086TT qw\u00E2\u0088\u0097 \u000F = \u00CE\u00BDM U2 N2max \u000F\u00E2\u0088\u0097, \u000FT = UU2T Nmax \u000FT\u00E2\u0088\u0097, (5.7) and dropping asterisks leads to the dimensionless energy equation \u00CE\u00B3PrT \u00CE\u00BA2 u \u00E2\u0088\u0082T \u00E2\u0088\u0082x = 1 \u00CE\u00BA2Re\u00CF\u0084 \u000F+ \u000FT \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 v\u00E2\u0080\u00B2T \u00E2\u0080\u00B2 \u00E2\u0088\u0092 \u00CE\u00B3 \u00E2\u0088\u0082 \u00E2\u0088\u0082x u\u00E2\u0080\u00B2T \u00E2\u0080\u00B2 + PrT PrRe\u00CF\u0084\u00CE\u00BA2 ( \u00E2\u0088\u00872T + \u00CE\u00B32\u00E2\u0088\u0082 2T \u00E2\u0088\u0082x2 ) . (5.8) The new dimensionless numbers are \u00CE\u00B3 = Nmax L , PrT = \u00CE\u00BA \u00E2\u0088\u0086T \u00E2\u0088\u0086TT , P r = \u00CE\u00B1 \u00CE\u00BDM . (5.9) We will assume fully rough conditions, so that the terms involving Re\u00CF\u0084 can be dropped. Then, to complete the model, a turbulence closure for the tur- bulent heat fluxes must be specified. An expedient model often considered in the literature is the natural extension of the eddy viscosity hypothesis, v\u00E2\u0080\u00B2T \u00E2\u0080\u00B2 = \u00E2\u0088\u0092\u00CE\u00BDT\u00E2\u0088\u0087T, u\u00E2\u0080\u00B2T \u00E2\u0080\u00B2 = \u00E2\u0088\u0092\u00CE\u00B3\u00CE\u00BDT \u00E2\u0088\u0082T \u00E2\u0088\u0082x . (5.10) 24Note that l \u00E2\u0088\u00BC \u00CE\u00BANmax. 41 5.1. Scaling Physically, the assumption is that the same turbulent motion that advects velocity fluctuations also advects temperature fluctuations; therefore, their correlations should take an analogous form, with the same eddy viscosity. They may differ by a constant of proportionality; this has been incorporated in PrT . With these assumptions, (5.8) becomes \u00CE\u00B3PrT \u00CE\u00BA2 u \u00E2\u0088\u0082T \u00E2\u0088\u0082x = \u000FT +\u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087T ) + \u00CE\u00B32\u00CE\u00BDT \u00E2\u0088\u0082 2T \u00E2\u0088\u0082x2 (5.11) This is a linear, inhomogeneous problem for the temperature field T . Fur- thermore, we know that u, \u00CE\u00BDT and \u000FT are independent of x. We therefore seek a solution in the form T = T\u00CC\u0084 (y, z) + e\u00CE\u00BBx T\u00CC\u0082 (y, z), (5.12) which is consistent with the homogeneous boundary conditions on T on \u00E2\u0088\u0082\u00E2\u0084\u00A6.25 The first component, T\u00CC\u0084 , contains the effect of dissipative heating on the temperature field, and satisfies \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087T\u00CC\u0084 ) = \u00E2\u0088\u0092l2|\u00E2\u0088\u0087u|3. (5.13) Evidently, this operates independently of the streamwise advection of heat. By integrating this equation over \u00E2\u0084\u00A6, we see that the mean wall heat flux due to dissipation is \u00E3\u0080\u0088q\u00CC\u0084w\u00E3\u0080\u0089\u00E2\u0084\u00A6 = \u00E2\u0088\u00921|\u00E2\u0088\u0082\u00E2\u0084\u00A6| \u00E2\u0088\u00AB \u00E2\u0084\u00A6 \u00CE\u00BDT\u00E2\u0088\u0087u \u00C2\u00B7 \u00E2\u0088\u0087u dS = \u00E2\u0088\u00921 |\u00E2\u0088\u0082\u00E2\u0084\u00A6| \u00E2\u0088\u00AB \u00E2\u0088\u0082\u00E2\u0084\u00A6 u \u00CE\u00BDT\u00E2\u0088\u0087u \u00C2\u00B7 n\u00CC\u0082 dl + 1|\u00E2\u0088\u0082\u00E2\u0084\u00A6| \u00E2\u0088\u00AB \u00E2\u0084\u00A6 u \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087u) dS = \u00E2\u0088\u00922\u00CE\u00BE |\u00E2\u0088\u0082\u00E2\u0084\u00A6| \u00E2\u0088\u00AB \u00E2\u0084\u00A6 u dS = \u00E2\u0088\u00922\u00E2\u0088\u009A2\u00CE\u00BA f , (5.14) where f = 8\u00CE\u00BA2/\u00E3\u0080\u0088u\u00E3\u0080\u0089\u00E2\u0084\u00A6 is the friction factor. This identity is a consequence of conservation of energy among the three levels of description: mean flow, turbulent fluctuations, and heat. In particular, the imposed pressure gradi- ent supplies energy to the mean flow; this is converted to turbulent kinetic energy by straining of the Reynolds stress along velocity gradients; turbu- lent energy is dissipated into heat; all the heat must escape at the wall as heat flux. 25To solve a given physical problem, in general a superposition of solutions is necessary. That is, e\u00CE\u00BBx T\u00CC\u0082 must be replaced by \u00E2\u0088\u00AB A\u00CE\u00BBe \u00CE\u00BBx T\u00CC\u0082\u00CE\u00BB d\u00CE\u00BB. 42 5.1. Scaling The temperature field due to advection, T\u00CC\u0082 , satisfies \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087T\u00CC\u0082 ) = \u00CE\u00BB\u00CE\u00B3PrT \u00CE\u00BA2 uT\u00CC\u0082 + \u00CE\u00B32\u00CE\u00BB2\u00CE\u00BDT T\u00CC\u0082 . (5.15) In the treatment of the velocity problem, we have already implicitly assumed that \u00CE\u00B3 \u001C 1. Since we also have \u00CE\u00BA \u00E2\u0089\u00A4 1, the streamwise diffusion term is much smaller than the advection term and may be discarded. If \u00CE\u00B3 \u001C \u00CE\u00BA2, as we expect, then both terms involving \u00E2\u0088\u0082T/\u00E2\u0088\u0082x may be discarded, and T\u00CC\u0083 is unnecessary. However, to facilitate comparison with the dissipation-induced temperature field, we note that T\u00CC\u0082 (y, z) satisfies \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CE\u00BDT\u00E2\u0088\u0087T\u00CC\u0082 ) = c\u00CE\u00BBuT\u00CC\u0082 , (5.16) with c\u00CE\u00BB = \u00CE\u00BB\u00CE\u00B3PrT /\u00CE\u00BA2. The equations governing \u00CE\u00BDT\u00E2\u0088\u0087T\u00CC\u0082 and \u00CE\u00BDT\u00E2\u0088\u0087T\u00CC\u0084 should be contrasted. The former will be largest where u is \u00E2\u0080\u0093 namely, near the centre of the pipe \u00E2\u0080\u0093 whereas the latter will be largest where gradients of u are large, namely, in the boundary layer. As with the momentum problem, we will first seek solutions for a circular pipe. 43 5.2. Circular Pipe 5.2 Circular Pipe For a circular pipe, (5.13) becomes d dn ( (1\u00E2\u0088\u0092 n)l(n)2dT\u00CC\u0084 dn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dudn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3) = \u00E2\u0088\u0092(1\u00E2\u0088\u0092 n)l(n)2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dudn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A33 . (5.17) We will use the technique of matched asymptotic expansions. Posing a naive expansion T\u00CC\u0084 = T\u00CC\u00840 + \u00CE\u00B4T\u00CC\u00841 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , (5.18) we find the leading order outer problem d dn ( (1\u00E2\u0088\u0092 n)3/2n (1 + c4n) dT\u00CC\u00840 dn ) = \u00E2\u0088\u0092(1\u00E2\u0088\u0092 n) 5/2(1 + c4n) n , (5.19) which can be integrated to dT\u00CC\u00840 dn = \u00E2\u0088\u0092(1 + c4n) r3/2n [ C1 + log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009Ar 1 + \u00E2\u0088\u009A r ) + 2 \u00E2\u0088\u009A r + 2 3 r3/2 + 2 5 r5/2 \u00E2\u0088\u0092 2 7 c4r 7/2 ] . (5.20) Here, we have used r = 1\u00E2\u0088\u0092n to shorten some of the expressions. Finite heat flux at n = 1 (r = 0) requires that C1 = 0. This can again be integrated, to T\u00CC\u00840 =C2 \u00E2\u0088\u0092 12 [ log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009Ar 1 + \u00E2\u0088\u009A r )]2 \u00E2\u0088\u0092 2(1 + c4)\u00E2\u0088\u009A r log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009Ar 1 + \u00E2\u0088\u009A r ) \u00E2\u0088\u0092 ( 16 15 \u00E2\u0088\u0092 16c4 7 ) log(n)\u00E2\u0088\u0092 ( 2 5 \u00E2\u0088\u0092 20c4 21 ) r + 12 35 c4r 2 \u00E2\u0088\u0092 2 21 c24r 3. (5.21) The constant C2 needs to be determined by matching to the boundary layer solution. For a fully rough pipe, the thickness of the thermal boundary layer is proportional to ks, just like the momentum boundary layer. Therefore, we may use the same boundary layer rescaling n = \u00CE\u00B4N . Denoting the temperature in the boundary layer by TI , we will expand TI as T\u00CC\u0084I = T\u00CC\u0084I0 + \u00CE\u00B4T\u00CC\u0084I1 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , (5.22) giving the leading order equation d dN ( (N + 1) dT\u00CC\u0084I0 dN ) = \u00E2\u0088\u0092 1 N + 1 . (5.23) 44 5.2. Circular Pipe Hence dT\u00CC\u0084I0 dN = \u00E2\u0088\u0092q\u00CC\u0084w \u00E2\u0088\u0092 log(N + 1) N + 1 , (5.24) which can be integrated to T\u00CC\u0084I0 = \u00E2\u0088\u0092q\u00CC\u0084w log(N + 1)\u00E2\u0088\u0092 12 [log(N + 1)] 2 . (5.25) Since T = 0 on n = 0, there is no integration constant. Matching with Van Dyke\u00E2\u0080\u0099s rule, as N \u00E2\u0086\u0092\u00E2\u0088\u009E, this becomes T\u00CC\u0084I0 \u00E2\u0088\u00BC \u00E2\u0088\u0092q\u00CC\u0084w log (n \u00CE\u00B4 ) \u00E2\u0088\u0092 1 2 [ log (n \u00CE\u00B4 )]2 +O(\u00CE\u00B4 log(\u00CE\u00B4)2) (5.26) This must be equated to the inner expansion of the outer solution, T\u00CC\u00840 \u00E2\u0088\u00BC C2 \u00E2\u0088\u0092 12 [ log (n 4 )]2 \u00E2\u0088\u0092 2(1 + c4) log (n4)\u00E2\u0088\u0092 ( 16 15 \u00E2\u0088\u0092 16c4 7 ) log(n) \u00E2\u0088\u0092 ( 2 5 \u00E2\u0088\u0092 20c4 21 ) + 12 35 c4 \u00E2\u0088\u0092 221c 2 4 +O(\u00CE\u00B4). (5.27) This implies q\u00CC\u0084w = log(\u00CE\u00B4) + c6, C2 = 1 2 log(\u00CE\u00B4)2 + c6 log(\u00CE\u00B4) + c7, (5.28) where c6 and c7 are \u00CE\u00B4-independent constants.26 One may easily check that the identity (5.14) is satisfied. To find these solutions, we neglected advection. In reality, the dissipated heat produced in the outer core of the flow is advected downstream, so that the temperature at the centre of the pipe continually increases. For our neglect of advection to be sound, we would hope that the bulk of dissipated heat is produced in the boundary layer and diffused out at the boundary. It is therefore worth considering the distribution of dissipation. We will consider an artificial separation of \u00E2\u0084\u00A6 into n < \u00CE\u00B4\u00CE\u00B7, and n > \u00CE\u00B4\u00CE\u00B7, with 0 < \u00CE\u00B7 < 1, so that the separation is in the matching region. Denoting the contribution of dissipation in the boundary layer and outer region by QBL, and QO, 26c6 = \u00E2\u0088\u0092 log(4)\u00E2\u0088\u0092 46/15\u00E2\u0088\u0092 42c4/7, c7 = 1/2 log(4)2 \u00E2\u0088\u0092 2(1 + c4) log(4) + 2/5\u00E2\u0088\u0092 20c4/21\u00E2\u0088\u0092 12c4/35 + 2c 2 4/21. 45 5.2. Circular Pipe respectively, we have, using (5.4), QBL = \u00E2\u0088\u00AB \u00CE\u00B4\u00CE\u00B7 0 \u00E2\u0088\u00AB 2pi 0 (1\u00E2\u0088\u0092 n)\u000FT d\u00CE\u00B8 dn \u00E2\u0089\u0088 \u00E2\u0088\u00AB \u00CE\u00B4\u00CE\u00B7 0 \u00E2\u0088\u00AB 2pi 0 (n+ \u00CE\u00B4)2 ( 1 n+ \u00CE\u00B4 )3 d\u00CE\u00B8 dn = 2pi \u00E2\u0088\u00AB \u00CE\u00B4\u00CE\u00B7 0 dn n+ \u00CE\u00B4 = 2pi log(\u00CE\u00B4\u00CE\u00B7\u00E2\u0088\u00921) +O(1). (5.29) The contribution from the outer region is QO = \u00E2\u0088\u00AB 1 \u00CE\u00B4\u00CE\u00B7 \u00E2\u0088\u00AB 2pi 0 (1\u00E2\u0088\u0092 n)\u000FT d\u00CE\u00B8 dn \u00E2\u0089\u0088 \u00E2\u0088\u00AB 1 \u00CE\u00B4\u00CE\u00B7 \u00E2\u0088\u00AB 2pi 0 (1\u00E2\u0088\u0092 n) ( n 1 + c4n )2(\u00E2\u0088\u009A1\u00E2\u0088\u0092 n(1 + c4n) n )3 d\u00CE\u00B8 dn = 2pi \u00E2\u0088\u00AB 1 \u00CE\u00B4\u00CE\u00B7 (1\u00E2\u0088\u0092 n)5/2 1 + c4n n dn = 2pi [ \u00E2\u0088\u00922c4 7 (1\u00E2\u0088\u0092 n)7/2 + 2 5 (1\u00E2\u0088\u0092 n)5/2 + 2 3 (1\u00E2\u0088\u0092 n)3/2 +2 \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n+ log ( 1\u00E2\u0088\u0092\u00E2\u0088\u009A1\u00E2\u0088\u0092 n 1 + \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n )]1 \u00CE\u00B4\u00CE\u00B7 = 2pi [ + 2c4 7 \u00E2\u0088\u0092 2 5 \u00E2\u0088\u0092 2 3 \u00E2\u0088\u0092 2\u00E2\u0088\u0092 log ( \u00CE\u00B4\u00CE\u00B7 4 ) +O(\u00CE\u00B4\u00CE\u00B7 log \u00CE\u00B4\u00CE\u00B7) ] = 2pi log(\u00CE\u00B4\u00E2\u0088\u0092\u00CE\u00B7) +O(1). (5.30) Which contribution is greater depends on the value of \u00CE\u00B7. For 0 < \u00CE\u00B7 < 1/2, QBL > QO, while for \u00CE\u00B7 > 1/2, QBL < QO. This sensitive dependence on \u00CE\u00B7 indicates that most of the dissipation occurs in the matching region. As \u00CE\u00B4 \u00E2\u0086\u0092 0, this region tends toward the wall, so that, indeed, most of the dissipation is produced near the wall. Heat transfer results are usually presented as a heat transfer coefficient, or Stanton number correlation.27 Abusing notation, this is, in dimensional terms, St = \u00E2\u0088\u0092qw/(cp\u00E3\u0080\u0088u\u00E3\u0080\u0089(T (n = 1) \u00E2\u0088\u0092 T (n = 0))). We will use an overbar 27Note that St = Nux/(RexPr), where Nux and Rex are the Nusselt and Reynolds numbers determined by an arbitrary length scale x, and Re is determined with the mean velocity. 46 5.2. Circular Pipe to indicate that this is for the dissipation-induced temperature field. In our notation, we have S\u00CC\u0084t \u00E2\u0089\u00A1 \u00E2\u0088\u0092\u00CF\u0081cpUT\u00E2\u0088\u0086TT q\u00CC\u0084w \u00CF\u0081cpU\u00E3\u0080\u0088u\u00E3\u0080\u0089\u00E2\u0084\u00A6\u00E2\u0088\u0086T T (n = 1) = \u00E2\u0088\u0092\u00CE\u00BA \u00E2\u0088\u009A f/8 PrT q\u00CC\u0084w T (n = 1) = \u00CE\u00BA \u00E2\u0088\u009A f/8 PrT \u00E2\u0088\u0092 log(\u00CE\u00B4)\u00E2\u0088\u0092 c6 1 2 log(\u00CE\u00B4) 2 + c6 log(\u00CE\u00B4) + c7 + 4(1 + c4) = 2\u00CE\u00BA \u00E2\u0088\u009A f/8 PrT ( log(1/\u00CE\u00B4)\u00E2\u0088\u0092 c6 + c 2 6 \u00E2\u0088\u0092 2c7 \u00E2\u0088\u0092 8(1 + c4) log(\u00CE\u00B4) + c6 )\u00E2\u0088\u00921 = \u00E2\u0088\u009A f3/16 PrT ( 1\u00E2\u0088\u0092 c 2 6 \u00E2\u0088\u0092 2c7 \u00E2\u0088\u0092 8(1 + c4) (log(\u00CE\u00B4) + c6)2 )\u00E2\u0088\u00921 . (5.31) As \u00CE\u00B4 \u00E2\u0086\u0092 0, the term in parentheses tends to unity, giving a simple form for S\u00CC\u0084t. Although there are accepted engineering relations for St, most are in- tended for smooth wall flows, where the Prandtl number plays a critical role in determining the relative size of thermal and viscous boundary layers. More importantly, accepted heat transfer relations are for advection domi- nated temperature fields. Still, it is useful to compare the relations. We will denote these advection-induced Stanton numbers by S\u00CC\u0082t. A simple relation for high Re that follows from the Izakson-Millikan argument, applied to the thermal problem, is [27] S\u00CC\u0082t = f/(8PrT ) (5.32) Engineers favour the Dittus-Boelter equation28 S\u00CC\u0082t = 0.023 ( 2 \u00CE\u00BE \u00E2\u0088\u009A 8 f Re\u00CF\u0084Pr 3 )\u00E2\u0088\u00921/5 . (5.33) Since the Prandtl number appears here as the ratio of viscous and thermal boundary layer thicknesses for a smooth wall, the most straightforward way to extend this result to rough walls is to set Pr = 1 and replace Re\u00CF\u0084 with exp(\u00E2\u0088\u0092\u00CE\u00BAB)/\u00CE\u00B4. 28This is usually presented as Nu = 0.023Re0.8Pr0.4, where the Nusselt and Reynolds numbers are based on the hydraulic diameter, and the Reynolds number is based on the mean velocity. 47 5.2. Circular Pipe 10\u00E2\u0088\u00926 10\u00E2\u0088\u00925 10\u00E2\u0088\u00924 10\u00E2\u0088\u00923 10\u00E2\u0088\u00922 10\u00E2\u0088\u00921 10\u00E2\u0088\u00924 10\u00E2\u0088\u00923 10\u00E2\u0088\u00922 10\u00E2\u0088\u00921 100 101 \u00CE\u00B4 St Present model Izakson\u00E2\u0088\u0092Millikan Dittus\u00E2\u0088\u0092Boelter Figure 5.1: Comparison of Stanton number correlations as a function of boundary layer thickness \u00CE\u00B4. The three relations are compared in Figure 5.1 with PrT = 0.8, \u00CE\u00BE = 1, and Pr = 1. All three show similar behaviour, but at small \u00CE\u00B4, the present result is more than five times smaller than the engineering relations. This indicates that, for a given temperature difference between the centre of the pipe and the wall, the dissipation-induced heat flux at the wall is much smaller than that caused by advection. For a general problem in which advection and dissipation are both im- portant, we have seen that the temperature fields are additive, viz., T = 48 5.3. Non-Circular Pipe T\u00CC\u0084 + e\u00CE\u00BBxT\u00CC\u0082 .29 Therefore, the Stanton number is St = \u00E2\u0088\u0092\u00CE\u00BA \u00E2\u0088\u009A f/8 PrT q\u00CC\u0084w + e\u00CE\u00BBxq\u00CC\u0082w T\u00CC\u0084 (n = 1) + e\u00CE\u00BBxT\u00CC\u0082 (n = 1) = \u00E2\u0088\u0092\u00CE\u00BA \u00E2\u0088\u009A f/8 PrT [ q\u00CC\u0084w T\u00CC\u0084 (n = 1) ( T\u00CC\u0084 (n = 1) T\u00CC\u0084 (n = 1) + e\u00CE\u00BBxT\u00CC\u0082 (n = 1) ) + q\u00CC\u0082w T\u00CC\u0082 (n = 1) ( e\u00CE\u00BBxT\u00CC\u0082 (n = 1) T\u00CC\u0084 (n = 1) + e\u00CE\u00BBxT\u00CC\u0082 (n = 1) )] = S\u00CC\u0084t ( T\u00CC\u0084 (n = 1) T\u00CC\u0084 (n = 1) + e\u00CE\u00BBxT\u00CC\u0082 (n = 1) ) + S\u00CC\u0082t ( e\u00CE\u00BBxT\u00CC\u0082 (n = 1) T\u00CC\u0084 (n = 1) + e\u00CE\u00BBxT\u00CC\u0082 (n = 1) ) (5.34) This indicates that the Stanton numbers for the dissipation and advection induced temperature fields are additive, provided they are weighted with the temperature difference carried by each field. As indicated by the factor e\u00CE\u00BBx, this will depend on the downstream distance. In particular, sufficiently far downstream, if \u00CE\u00BB > 0, the advection effects should dominate, while if \u00CE\u00BB < 0, dissipation will dominate. 5.3 Non-Circular Pipe The difficulties that appear in solving the thermal problem in a non-circular pipe are analogous to those earlier dealt with for the velocity problem. Ac- cordingly, we will follow the same procedure in determining approximate solutions, and omit some details. We will neglect streamwise advection of heat. In (n, \u00CE\u00B8) coordinates, (5.13) is 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082n (\u00E2\u0088\u009A G \u00CE\u00BDT \u00E2\u0088\u0082T \u00E2\u0088\u0082n ) + 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B8 ( \u00CE\u00BDT\u00E2\u0088\u009A G \u00E2\u0088\u0082T \u00E2\u0088\u0082\u00CE\u00B8 ) = \u00E2\u0088\u0092l2|\u00E2\u0088\u0087u|3. (5.35) Again, we attempt solution with matched asymptotic expansions. We will first consider the boundary layer problem. Writing TI = TI0 + \u00CE\u00B4TI1 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 for the solution in the inner region, the leading order problem is \u00E2\u0088\u0082 \u00E2\u0088\u0082N ( (N + 1) \u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8) \u00E2\u0088\u0082TI0 \u00E2\u0088\u0082N ) = \u00E2\u0088\u0092\u00CF\u0084w(\u00CE\u00B8)3/2 N + 1 . (5.36) 29As earlier, for a general problem e\u00CE\u00BBx must be replaced by \u00E2\u0088\u00AB A\u00CE\u00BBe \u00CE\u00BBxd\u00CE\u00BB. 49 5.3. Non-Circular Pipe Hence \u00E2\u0088\u0082TI0 \u00E2\u0088\u0082N = \u00E2\u0088\u0092qw(\u00CE\u00B8)\u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8)(N + 1) \u00E2\u0088\u0092 \u00CF\u0084w(\u00CE\u00B8) N + 1 log(N + 1), (5.37) which may be integrated to TI0 = \u00E2\u0088\u0092qw(\u00CE\u00B8)\u00E2\u0088\u009A \u00CF\u0084w(\u00CE\u00B8) log(N + 1)\u00E2\u0088\u0092 \u00CF\u0084w(\u00CE\u00B8) 2 [log(N + 1)]2. (5.38) An immediate consequence of this expression is that the wall heat flux must vanish whenever the wall stress does. Like the wall stress in the mechanical problem, the wall heat flux must be determined by matching to the outer problem. The general outer problem is intractable. As before, we have found nu- merical solutions with COMSOL, a finite element package. These solutions are compared to analytical results in what follows. To make analytical progress, we will employ similar approximations as in chapter 4. Proceeding as before, we will discard the explicit \u00CE\u00B8 derivatives. Then, (5.35) becomes an ordinary differential equation in n with parametric \u00CE\u00B8 dependence. Integrating once, q \u00E2\u0089\u00A1 \u00E2\u0088\u0092\u00CE\u00BDT \u00E2\u0088\u0082T \u00E2\u0088\u0082n = \u00E2\u0088\u0092C(\u00CE\u00B8)\u00E2\u0088\u009A G + 1\u00E2\u0088\u009A G \u00E2\u0088\u00AB n 0 \u00E2\u0088\u009A G l2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u0082u\u00E2\u0088\u0082n \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A33 dn\u00E2\u0080\u00B2. (5.39) Here, q is a reduced heat flux, analogous to the reduced stress of the mo- mentum problem. The unknown function C(\u00CE\u00B8) must be related to qw by matching, and then determined by regularity of the solution in the outer region. As before, we will only consider this regularity condition explicitly for symmetric shapes. In this case, on n = nmax(\u00CE\u00B8), T satisfies 0 = \u00E2\u0088\u0082T \u00E2\u0088\u0082n \u00E2\u0088\u0092 1 G dnmax(\u00CE\u00B8) d\u00CE\u00B8 \u00E2\u0088\u0082T \u00E2\u0088\u0082\u00CE\u00B8 . (5.40) The system of equations (5.39) and (5.40) is still too difficult to solve an- alytically. Therefore, as in the momentum problem, we will first consider approximating the velocity and temperature fields, and applying the exact boundary condition. 5.3.1 Constant Heat Flux Approximation For the same reasons as in chapter 4, a reasonable and simple ansatz for the form of the temperature field is to assume that the inner solution is valid in the outer region, that is, to use equation (5.35). Physically, this 50 5.3. Non-Circular Pipe corresponds to assuming a constant heat flux in the outer region. If making this assumption, it is reasonable to also assume that the velocity takes the analogous form, (4.15). Then, note that T0 = \u00E2\u0088\u0092qw\u00E2\u0088\u009A \u00CF\u0084w log (n \u00CE\u00B4 + 1 ) \u00E2\u0088\u0092 1 2 u20. (5.41) Under these assumptions, (5.40) becomes 0 = \u00E2\u0088\u0092qw\u00E2\u0088\u009A \u00CF\u0084w 1 nmax \u00E2\u0088\u0092 1 G dnmax d\u00CE\u00B8 log (n \u00CE\u00B4 + 1 ) \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B8 (\u00E2\u0088\u0092qw\u00E2\u0088\u009A \u00CF\u0084w ) \u00E2\u0088\u0092 u0 ( \u00E2\u0088\u0082u0 \u00E2\u0088\u0082n \u00E2\u0088\u0092 1 G dnmax d\u00CE\u00B8 \u00E2\u0088\u0082u0 \u00E2\u0088\u0082\u00CE\u00B8 ) . (5.42) This is a first order, linear equation in qw/ \u00E2\u0088\u009A \u00CF\u0084w, whose solution can im- mediately be written down. However, if we assume that the wall stress is determined by (4.14), as is consistent with the assumed form of u0, then the source term vanishes identically, and we see that qw/ \u00E2\u0088\u009A \u00CF\u0084w satisfies the same equation as \u00E2\u0088\u009A \u00CF\u0084w. Since the equation admits the trivial scaling symmetry, we have, immediately, \u00E2\u0088\u0092qw\u00E2\u0088\u009A \u00CF\u0084w \u00E2\u0088\u009D \u00E2\u0088\u009A\u00CF\u0084w. Setting the constant of proportionality so that the correct mean wall heat flux is attained, we have qw = \u00E2\u0088\u00922\u00E2\u0088\u009A2\u00CE\u00BA f \u00CF\u0084w, (5.43) a surprising result, considering that the temperature and velocity fields dif- fer. This result allows the solutions for \u00CF\u0084w from Section (4.1) to be im- mediately applied. As an example, in Figure 5.2, qw for a square pipe, as predicted by (5.43), is compared to the solution of the full PDE (5.13). The agreement is very good. 5.3.2 Zero Heat Flux Approximation An alternative approximation, analogous to letting the reduced stress vanish on n = nmax, is to let the reduced turbulent heat flux q vanish there. To obtain an expression for q, we need to compute the dissipation integral in (5.39). A simple approximation that makes this possible is again to ignore the minor effect of wall curvature. Then q = \u00E2\u0088\u0092C(\u00CE\u00B8)\u00E2\u0088\u009A Gw \u00E2\u0088\u0092 c4 5\u00CE\u00BE \u00CF\u0084\u00CC\u00835/2 + 2 3 \u00CF\u0084\u00CC\u00833/2 + 2\u00CF\u0084w \u00E2\u0088\u009A \u00CF\u0084\u00CC\u0083 \u00E2\u0088\u0092 \u00CF\u00843/2w log (\u00E2\u0088\u009A \u00CF\u0084w + \u00E2\u0088\u009A \u00CF\u0084\u00CC\u0083 \u00E2\u0088\u009A \u00CF\u0084w \u00E2\u0088\u0092 \u00E2\u0088\u009A \u00CF\u0084\u00CC\u0083 ) , (5.44) 51 5.3. Non-Circular Pipe 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.2 0.4 0.6 0.8 1 1.2 s q w (s) / < q w > PDE Constant Heat Flux Figure 5.2: Heat flux in a square pipe, normalized by its mean. One eighth of the pipe boundary is shown. The corner is at right. with \u00CF\u0084\u00CC\u0083 = \u00CF\u0084w \u00E2\u0088\u0092 2\u00CE\u00BEn the approximate reduced wall stress. This must be matched to the outer asymptote of the inner solution (5.37), qI0 \u00E2\u0088\u00BC qw + \u00CF\u00843/2w log (n \u00CE\u00B4 ) . (5.45) Hence C\u00E2\u0088\u009A Gw = \u00E2\u0088\u0092qw \u00E2\u0088\u0092 c45\u00CE\u00BE \u00CF\u0084 5/2 w + 2 3 \u00CF\u00843/2w + 2\u00CF\u0084 3/2 w \u00E2\u0088\u0092 \u00CF\u00843/2w log ( 2\u00CF\u0084w \u00CE\u00BE\u00CE\u00B4 ) . (5.46) Letting q = 0 when \u00CF\u0084\u00CC\u0083 = 0, we find30 C = 0 and therefore qw(\u00CE\u00B8) = \u00E2\u0088\u0092\u00CF\u00843/2w ( log ( 2\u00CF\u0084w \u00CE\u00BE\u00CE\u00B4 ) + c4 5\u00CE\u00BE \u00CF\u0084w \u00E2\u0088\u0092 83 ) . (5.47) 30This is slightly inconsistent, in that we earlier assumed that \u00CF\u0084 = 0 on n = nmax. We cannot make that assumption here, because (5.44) requires \u00CF\u0084\u00CC\u0083 \u00E2\u0089\u00A5 0. The best way to proceed would be to retain \u00CF\u0084 in the dissipation integral. This can be done, also somewhat inconsistently, by ignoring only the effect of wall curvature in the denominator of \u00E2\u0088\u0082u/\u00E2\u0088\u0082n. This leads to a much longer formula for q, without much improving the approximation. 52 5.3. Non-Circular Pipe In Figure 5.3, this approximation is tested for a semi-superellipse against the full PDE (5.13). Agreement is only slightly worse than the accuracy of the \u00CF\u0084w approximation (see Figure 4.7). 2 4 6 8 10 12 14 16 0 0.5 1 1.5 2 s q w (s) / < q w > PDE Zero Heat Flux Figure 5.3: Heat flux in a semi-elliptical pipe, normalized by its mean. 53 Chapter 6 Conclusion In this thesis, we have considered turbulent flow in a rough pipe of arbitrary shape. Although many engineering problems fall into this category, these are usually of a few standard shapes, such as the square pipe, which have been extensively studied. Instead, we focus on the universal aspects of the problem, so that solutions may be obtained in evolving geometries, as occurs in geophysical applications. One such universal aspect is the boundary layer structure of the solu- tion. First, using the heuristic Izakson-Millikan argument, and later, using matched asymptotic expansions, we have exploited this structure to solve both the mechanical and thermal problems in the boundary layer. It was found that, for a general geometry, this leads to functions of \u00CE\u00B8, which must be determined by matching to the outer problem. To make progress with the outer problem, we motivated approximate methods of solution, based on simple physical hypotheses. For the mechan- ical problem, these were, first, assuming that the shear stress is constant in the outer region, and second, assuming that the shear stress goes to zero at the centre of the pipe. Although neither hypothesis strictly holds, they were shown to furnish good approximations to the distribution of shear stress at the wall. Unfortunately, neither approximation is suitable for all geome- tries. The constant-stress approximation was found to be successful for a square pipe, while the zero-stress approximation worked for an ellipse and a semi-superellipse. For the thermal problem, the analogous approximations are, first, as- suming a constant heat flux in the outer region, and second, assuming that the heat flux goes to zero at the centre of the pipe. The first approximation was shown to work for a square pipe, and the second for a semi-superellipse. Discrepancy between the solution of the full partial differential equations and data was explained in terms of the neglected secondary flow. No attempt at modelling this flow was attempted, because the simple eddy viscosity approach used here is insufficient [23]. However, even without the secondary flow, the solutions may be useful for geophysical applications. One way to apply these results is to pose a physical problem in which 54 Chapter 6. Conclusion the pipe shape is restricted to a certain family of shapes; for example, semi- ellipses of varying corner exponent, or ellipses of varying aspect ratio. Then, one can couple the analytical results for wall quantities with the physics that determines the pipe shape evolution. For example, in the glaciological application, the melt rate is proportional to the heat flux at the wall. Once this is known, the time step can be advanced, and the new shape fit to the family of shapes. Also, the results for the simple circular pipe are novel in several respects. First, they have been analytically integrated, from the boundary layer all the way to the centre of the pipe. Usually, results are quoted as quadratures [27], or are only accurate in limited regions. Because they are so simple, they form a model problem for matched asymptotic expansions. Second, the results are presented for fully rough pipes, which appear much more infrequently in the literature than smooth pipes. Third, the thermal problem was solved for the novel case in which the turbulent dissipation is important. As discussed above, this led to some differences from the normal advection dominated thermal problem. The results may be relevant in glaciology. From a theoretical point of view, one might hope to glean some insight from the approximate solutions presented for non-circular pipes. The success of the \u00E2\u0080\u009Cradial\u00E2\u0080\u009D approximation suggests that, indeed, the most important coordinate is the distance from the nearest wall. This is reflected in the presented solutions, which heavily feature nmax(\u00CE\u00B8). Within the context of mixing-length turbulence models, this is easy to interpret: velocity and temperature solutions depend strongly on the size of the largest nearby eddies. 55 Chapter 6. Conclusion Table 6.1: Glossary of Symbols \u00E2\u0084\u00A6 Pipe cross-section \u00E2\u0088\u0082\u00E2\u0084\u00A6 Boundary of \u00E2\u0084\u00A6 u Streamwise mean flow velocity v = (v, w) Transverse mean flow velocity u\u00E2\u0080\u00B2 Streamwise fluctuating velocity (v\u00E2\u0080\u00B2, w\u00E2\u0080\u00B2) Transverse fluctuating velocity P Pressure T Mean flow temperature T \u00E2\u0080\u00B2 Fluctuating temperature \u00CE\u00BDT Turbulent viscosity \u00CE\u00BDM Molecular viscosity \u00CE\u00B1 Thermal diffusivity cp Heat capacity at constant pressure \u00CF\u0081 Fluid density 56 Bibliography [1] A. Aly, A. Trupp, and A. Gerrard. Measurements and prediction of fully developed turbulent flow in an equilateral triangular duct. J. Fluid Mech., 85:57\u00E2\u0080\u009383, 1978. [2] G. B. Arfken and H. J. Weber. Mathematical Methods for Physicists. Academic Press, Cambridge, U.K., fourth edition, 1994. [3] J. Boussinesq. Theorie de lecoulement tourbillant. Mem. Pres. Acad. Sci., XXIII:46, 1877. [4] G. Clarke. Hydraulics of subglacial outburst floods: new insights from the Spring-Hutter formulation. J. Glac., 49(165), 2003. [5] S. Gavrilikis. Numerical simulation of low-reynolds-number turbulent flow through a straight square duct. J. Fluid Mech., 244:101\u00E2\u0080\u0093129, 1992. [6] H. Herwig, D. Gloss, and T. Wenterodt. A new approach to understand- ing and modelling the influence of wall roughness on friction factors for pipe and channel flows. J. Fluid Mech., 613:35\u00E2\u0080\u009353, 2008. [7] E. J. Hinch. Matched Asymptotic Expansions. Cambridge, Cambridge, U.K., first edition, 1991. [8] A. Huser and S. Biringen. Direct numerical simulation of turbulent flow in a square duct. J. Fluid Mech., 257:65\u00E2\u0080\u009395, 1993. [9] A. Izakson. On the formula for the velocity distribution near walls. Tech. Phys. USSR IV, 2:155\u00E2\u0080\u0093162, 1937. [10] J. Jime\u00CC\u0081nez. Turbulent flows over rough walls. Annu. Rev. Fluid Mech., 36:173\u00E2\u0080\u0093196, 2004. [11] A. Kolmogorov. The equations of turbulent motion in an incompressible fluid. Izv. Acad. Sci. USSR; Phys., 6:56\u00E2\u0080\u009358, 1942. [12] B. Launder and W. Ying. Prediction of flow and heat transfer in ducts of square cross-section. Proc. Inst. Mech. Engrs., 187:455\u00E2\u0080\u0093461, 1973. 57 Chapter 6. Bibliography [13] V. S. L\u00E2\u0080\u0099vov, I. Procaccia, and O. Rudenko. Universal model of finite Reynolds number turbulent flow in channels and pipes. Phys. Rev. Lett., 100(054504), 2008. [14] B. J. McKeon, J. Li, W. Jiang, J. F. Morrison, and A. J. Smits. Further observations on the mean velocity distribution in fully developed pipe flow. J. Fluid Mech., 501:135\u00E2\u0080\u0093147, 2004. [15] G. L. Mellor. The large reynolds number asymptotic theory of turbulent boundary layers. Int. J. Eng. Sci, 10:851, 1972. [16] C. Millikan. A critical discussion of turbulent flows in channel and circular tubes. In Proc. Fifth Intl Congr. of Applied Mechanics, pages 386\u00E2\u0080\u0093392. Wiley, 1939. [17] G. Mompean, S. Gavrilakis, L. Machiels, and M. O. Deville. On pre- dicting the turbulence-induced secondary flows using nonlinear k \u00E2\u0088\u0092 \u000F models. Phys. Fluids, 8:1856\u00E2\u0080\u00931868, 1996. [18] F. S. L. Ng. Mathematical modelling of subglacial drainage and erosion. PhD thesis, Oxford Univ., Oxford, U.K., 1998. [19] F. Nietzsche. Beyond good and evil. On- line: http://records.viu.ca/~johnstoi/Nietzsche/ beyondgoodandevil_tofc.htm, 2009. Trans. Ian Johnston. [20] N. Nikitin and A. Yakhot. Direct numerical simulation of turbulent flow in elliptical ducts. J. Fluid Mech., 532:141\u00E2\u0080\u0093164, 2005. [21] J. Nikuradse. Gesetzma\u00CC\u0088ssigkeit der turbulenten stro\u00CC\u0088mung in glatten rohren. Forsch. Arb. Ing.-Wes., (356), 1932. [22] R. Panton. Composite asymptotic expansions and scaling wall turbu- lence. Phil. Trans. R. Soc. A, 365:733\u00E2\u0080\u0093754, 2007. [23] S. B. Pope. Turbulent Flows. Cambridge Univ. Press, Cambridge, U.K., first edition, 2000. [24] L. Prandtl. Essentials of Fluid Dynamics. Blackie & Son, 1952. [25] M. R. Raupach, R. Antonia, and S. Rajagopalan. Rough-wall turbulent boundary layers. Appl. Mech. Rev., 44:1\u00E2\u0080\u009325, 1991. [26] H. Schlichting. Boundary Layer Theory. Springer, seventh edition, 1979. 58 [27] H. Schlichting and K. Gersten. Boundary Layer Theory. Springer, eighth edition, 2000. [28] M. P. Schultz and K. A. Flack. The rough-wall turbulent boundary layer from the hydraulically smooth to the fully rough regime. J. Fluid Mech., 580:381\u00E2\u0080\u0093405, 2007. [29] M. A. Shockling, J. J. Allen, and A. Smits. Roughness effects in tur- bulent pipe flow. J. Fluid Mech., 564:267\u00E2\u0080\u0093285, 2006. [30] C. Speziale, R. Abid, and E. C. Anderson. Critical evaluation of two- equation models for near-wall turbulence. A.I.A.A., 30(2):324\u00E2\u0080\u0093331, 1992. [31] A. A. Townsend. The Structure of Turbulent Shear Flow. Cambridge Univ. Press, Cambridge, U.K., second edition, 1976. [32] H. Xu. Direct numerical simulation of turbulence in a square annular duct. J. Fluid Mech., 621:23\u00E2\u0080\u009357, 2009. [33] K. S. Yajnik. Asymptotic theory of turbulent shear flows. J. Fluid Mech., 42:411\u00E2\u0080\u0093427, 1970. [34] K. S. Yajnik. Analysis of turbulent pipe and channel flows at moderately large Reynolds number. J. Fluid Mech., 61:23\u00E2\u0080\u009331, 1973. 59 Appendix A Exact Solution for the Velocity Profile in a Circular Pipe For a fully rough, circular pipe, the streamwise momentum equation is \u00E2\u0088\u0092 2 = 1 1\u00E2\u0088\u0092 n d dn ( (1\u00E2\u0088\u0092 n) ( n+ \u00CE\u00B4 1 + c4n )2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dudn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dudn ) , (A.1) where we have taken c3 = c as in Chapter 3. Integrating this once yields\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3dudn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 dudn = (1\u00E2\u0088\u0092 2n+ n2 + C0)(1 + c4n)2(n+ \u00CE\u00B4)2(1\u00E2\u0088\u0092 n) . (A.2) At the centre of the pipe, n = 1 and du/dn = 0, so that the integration constant C0 = 0. Assuming du/dn > 0 and integrating again, the velocity is u = \u00E2\u0088\u009A 1 + \u00CE\u00B4 (1\u00E2\u0088\u0092 c4\u00CE\u00B4) log (\u00E2\u0088\u009A 1 + \u00CE\u00B4 \u00E2\u0088\u0092\u00E2\u0088\u009A1\u00E2\u0088\u0092 n\u00E2\u0088\u009A 1 + \u00CE\u00B4 + \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n ) + 2(1\u00E2\u0088\u0092 c4\u00CE\u00B4) \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n\u00E2\u0088\u0092 2c4 3 (1\u00E2\u0088\u0092 n)3/2 + C1. (A.3) Applying the no-slip condition at n = 0, we find C1 = \u00E2\u0088\u0092 \u00E2\u0088\u009A 1 + \u00CE\u00B4 (1\u00E2\u0088\u0092 c4\u00CE\u00B4) log (\u00E2\u0088\u009A 1 + \u00CE\u00B4 \u00E2\u0088\u0092 1\u00E2\u0088\u009A 1 + \u00CE\u00B4 + 1 ) \u00E2\u0088\u0092 2(1\u00E2\u0088\u0092 c4\u00CE\u00B4) + 2c43 . (A.4) By expanding in \u00CE\u00B4, the exact solution (A.3) can be used to obtain the outer solution of the matched asymptotic expansion, (3.3). Similarly, by rescaling for the boundary layer and expanding in \u00CE\u00B4, the inner solution can also be obtained. 60 Appendix B Geometrical Transformations Herein we present the transformations necessary to use the (n, \u00CE\u00B8) coordinate system. A general orthogonal coordinate system (\u00C2\u00B5, \u00CE\u00B7) defined in the Cartesian y = (y, z) plane is described by a line element ds2 = Ed\u00CE\u00B72 +Gd\u00C2\u00B52, (B.1) in which G = |\u00E2\u0088\u0082y/\u00E2\u0088\u0082\u00C2\u00B5|2 and E = |\u00E2\u0088\u0082y/\u00E2\u0088\u0082\u00CE\u00B7|2. The unit vectors in the direction of increasing \u00C2\u00B5 and \u00CE\u00B7 are e\u00CC\u0082\u00C2\u00B5 = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082y \u00E2\u0088\u0082\u00C2\u00B5 , e\u00CC\u0082\u00CE\u00B7 = 1\u00E2\u0088\u009A E \u00E2\u0088\u0082y \u00E2\u0088\u0082\u00CE\u00B7 . (B.2) The volume element dS = \u00E2\u0088\u009A EGd\u00C2\u00B5d\u00CE\u00B7 is the area of the rectangle spanned by \u00E2\u0088\u009A Gd\u00C2\u00B5e\u00CC\u0082\u00C2\u00B5 and \u00E2\u0088\u009A Ed\u00CE\u00B7e\u00CC\u0082\u00CE\u00B7. These relations can be used to define the familiar gradient, divergence, and curl operators. In particular, for a scalar function \u00CF\u0086(\u00C2\u00B5, \u00CE\u00B7), the gradient \u00E2\u0088\u0087\u00CF\u0086 is defined by [2] \u00E2\u0088\u0087\u00CF\u0086 = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082\u00CF\u0086 \u00E2\u0088\u0082\u00C2\u00B5 e\u00CC\u0082\u00C2\u00B5 + 1\u00E2\u0088\u009A E \u00E2\u0088\u0082\u00CF\u0086 \u00E2\u0088\u0082\u00CE\u00B7 e\u00CC\u0082\u00CE\u00B7. (B.3) For a vector function \u00CE\u00A6(\u00C2\u00B5, \u00CE\u00B7) = (\u00CE\u00A6\u00C2\u00B5,\u00CE\u00A6\u00CE\u00B7), the divergence \u00E2\u0088\u0087 \u00C2\u00B7\u00CE\u00A6 is defined by [2] \u00E2\u0088\u0087 \u00C2\u00B7 \u00CE\u00A6 = 1\u00E2\u0088\u009A EG \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00C2\u00B5 (\u00E2\u0088\u009A E\u00CE\u00A6\u00C2\u00B5 ) + 1\u00E2\u0088\u009A EG \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 (\u00E2\u0088\u009A G\u00CE\u00A6\u00CE\u00B7 ) . (B.4) We will not need the curl operator. B.1 Curvature It will be helpful to have expressions for the curvature of the contours of \u00C2\u00B5 and \u00CE\u00B7. To derive these, consider a contour \u00CE\u00B7 =constant. e\u00CC\u0082\u00C2\u00B5 is a unit tangent to this curve. (see Figure B.1). Since |e\u00CC\u0082\u00C2\u00B5|2 = 1, we have e\u00CC\u0082\u00C2\u00B5 \u00C2\u00B7 \u00E2\u0088\u0082e\u00CC\u0082\u00C2\u00B5/\u00E2\u0088\u0082\u00C2\u00B5 = 0, 61 B.1. Curvature Figure B.1: Geometrical setup for curvature derivation. and therefore \u00E2\u0088\u0082e\u00CC\u0082\u00C2\u00B5/\u00E2\u0088\u0082\u00C2\u00B5 \u00E2\u0088\u009D e\u00CC\u0082\u00E2\u008A\u00A5\u00C2\u00B5 = e\u00CC\u0082\u00CE\u00B7. The curve has curvature k\u00CE\u00B7, defined by \u00E2\u0088\u0082e\u00CC\u0082\u00C2\u00B5 \u00E2\u0088\u0082\u00C2\u00B5 = k\u00CE\u00B7 \u00E2\u0088\u009A Ge\u00CC\u0082\u00CE\u00B7. (B.5) Since the contour is part of a coordinate system, there are analogous contours of \u00C2\u00B5 at right angles to the contours of \u00CE\u00B7. This fact can be exploited by expanding the identity \u00E2\u0088\u00822y \u00E2\u0088\u0082\u00CE\u00B7\u00E2\u0088\u0082\u00C2\u00B5 = \u00E2\u0088\u00822y \u00E2\u0088\u0082\u00C2\u00B5\u00E2\u0088\u0082\u00CE\u00B7 \u00E2\u0087\u0094 \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 ( \u00E2\u0088\u009A Ge\u00CC\u0082\u00C2\u00B5) = \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00C2\u00B5 ( \u00E2\u0088\u009A Ee\u00CC\u0082\u00CE\u00B7) \u00E2\u0087\u0094 \u00E2\u0088\u0082 \u00E2\u0088\u009A G \u00E2\u0088\u0082\u00CE\u00B7 e\u00CC\u0082\u00C2\u00B5 + \u00E2\u0088\u009A G \u00E2\u0088\u0082e\u00CC\u0082\u00C2\u00B5 \u00E2\u0088\u0082\u00CE\u00B7 = \u00E2\u0088\u0082 \u00E2\u0088\u009A E \u00E2\u0088\u0082\u00C2\u00B5 e\u00CC\u0082\u00CE\u00B7 + \u00E2\u0088\u009A E \u00E2\u0088\u0082e\u00CC\u0082\u00CE\u00B7 \u00E2\u0088\u0082\u00C2\u00B5 (B.6) Taking the inner product of this expression with e\u00CC\u0082\u00C2\u00B5, \u00E2\u0088\u0082 \u00E2\u0088\u009A G \u00E2\u0088\u0082\u00CE\u00B7 = \u00E2\u0088\u009A E e\u00CC\u0082\u00C2\u00B5 \u00C2\u00B7 \u00E2\u0088\u0082e\u00CC\u0082\u00CE\u00B7 \u00E2\u0088\u0082\u00C2\u00B5 = \u00E2\u0088\u009A E ( \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00C2\u00B5 (e\u00CC\u0082\u00C2\u00B5 \u00C2\u00B7 e\u00CC\u0082\u00CE\u00B7)\u00E2\u0088\u0092 e\u00CC\u0082\u00CE\u00B7 \u00C2\u00B7 \u00E2\u0088\u0082e\u00CC\u0082\u00C2\u00B5 \u00E2\u0088\u0082\u00C2\u00B5 ) = \u00E2\u0088\u0092 \u00E2\u0088\u009A EG k\u00CE\u00B7. (B.7) By symmetry, we may obtain the analogous expression for k\u00C2\u00B5, defined by \u00E2\u0088\u0082e\u00CC\u0082\u00CE\u00B7/\u00E2\u0088\u0082\u00CE\u00B7 = \u00E2\u0088\u009A Gk\u00C2\u00B5e\u00CC\u0082 \u00E2\u008A\u00A5 \u00CE\u00B7 . Then 31 k\u00CE\u00B7 = \u00E2\u0088\u00921\u00E2\u0088\u009A EG \u00E2\u0088\u0082 \u00E2\u0088\u009A G \u00E2\u0088\u0082\u00CE\u00B7 , k\u00C2\u00B5 = 1\u00E2\u0088\u009A EG \u00E2\u0088\u0082 \u00E2\u0088\u009A E \u00E2\u0088\u0082\u00C2\u00B5 . (B.8) 31Note that the handedness of the coordinates changes under the symmetry transfor- mation. This leads to the sign change. 62 B.2. Wall Coordinates B.2 Wall Coordinates In the particular coordinate system in which \u00CE\u00B7 = n is distance from the wall, and \u00C2\u00B5 = \u00CE\u00B8 is a coordinate perpendicular to n, the above expressions simplify considerably. It is well known that n satisfies the Eikonal equation |\u00E2\u0088\u0087n| = 1. Using (B.3) we see that this implies E \u00E2\u0089\u00A1 1. Immediately, k\u00CE\u00B8 = 0, so that contours of \u00CE\u00B8 are straight lines, as expected. This fact may be used to obtain the functional form of G. Let the pipe boundary \u00E2\u0088\u0082\u00E2\u0084\u00A6 be parameterized by Figure B.2: Geometrical setup for G derivation. yw(\u00CE\u00B8) and consider a closed loop \u00CE\u0093 along contours of n and \u00CE\u00B8, as shown in Figure B.2. Since the contours of \u00CE\u00B8 are straight lines, e\u00CC\u0082n is independent of n. Then 0 = \u00E2\u0088\u00AB \u00CE\u0093 d` = \u00E2\u0088\u0092 \u00E2\u0088\u00AB \u00CE\u00B80+\u00CE\u00B4\u00CE\u00B8 \u00CE\u00B80 \u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8\u00E2\u0080\u00B2 d\u00CE\u00B8\u00E2\u0080\u00B2 + \u00E2\u0088\u00AB n0 0 e\u00CC\u0082n(\u00CE\u00B80)dn\u00E2\u0080\u00B2 + \u00E2\u0088\u00AB \u00CE\u00B80+\u00CE\u00B4\u00CE\u00B8 \u00CE\u00B80 e\u00CC\u0082\u00CE\u00B8 \u00E2\u0088\u009A Gd\u00CE\u00B8\u00E2\u0080\u00B2 \u00E2\u0088\u0092 \u00E2\u0088\u00AB n0 0 e\u00CC\u0082n(\u00CE\u00B80 + \u00CE\u00B4\u00CE\u00B8)dn\u00E2\u0080\u00B2 = ( \u00E2\u0088\u0092\u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8 (\u00CE\u00B80)\u00E2\u0088\u0092 n0\u00E2\u0088\u0082e\u00CC\u0082n \u00E2\u0088\u0082\u00CE\u00B8 (\u00CE\u00B80) + e\u00CC\u0082\u00CE\u00B8(n0, \u00CE\u00B80) \u00E2\u0088\u009A G(n0, \u00CE\u00B80) ) \u00CE\u00B4\u00CE\u00B8 +O((\u00CE\u00B4\u00CE\u00B8)2) (B.9) 63 B.3. Coordinate Derivation Taking the limit \u00CE\u00B4\u00CE\u00B8 \u00E2\u0086\u0092 0 and taking the inner product with e\u00CC\u0082\u00CE\u00B8, we find \u00E2\u0088\u009A G = e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8 + n0e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082e\u00CC\u0082n \u00E2\u0088\u0082\u00CE\u00B8 = e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8 \u00E2\u0088\u0092 n0e\u00CC\u0082n \u00C2\u00B7 \u00E2\u0088\u0082e\u00CC\u0082\u00CE\u00B8 \u00E2\u0088\u0082\u00CE\u00B8 = e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8 + n0 \u00E2\u0088\u009A Gkn = e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082yw \u00E2\u0088\u0082\u00CE\u00B8 + n0 \u00E2\u0088\u0082 \u00E2\u0088\u009A G \u00E2\u0088\u0082n (B.10) The first term on the right hand side is independent of \u00CE\u00B80. Since this equation is true for all n0 and \u00CE\u00B80, we may drop the subscripts and consider it as an ordinary differential equation in n. Evaluating this equation at n0 = 0, we see that e\u00CC\u0082\u00CE\u00B8 \u00C2\u00B7 \u00E2\u0088\u0082yw/\u00E2\u0088\u0082\u00CE\u00B8 = \u00E2\u0088\u009A Gw(\u00CE\u00B8), the value of \u00E2\u0088\u009A G at the wall. The equation can be integrated to\u00E2\u0088\u009A G(n, \u00CE\u00B8) = \u00E2\u0088\u009A Gw(\u00CE\u00B8)(1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n). (B.11) kw is the constant of integration. However, the curvature of contours of n is then kn(n, \u00CE\u00B8) = \u00E2\u0088\u00921\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u009A G \u00E2\u0088\u0082n = kw(\u00CE\u00B8) 1\u00E2\u0088\u0092 kw(\u00CE\u00B8)n. (B.12) At the wall, kn = kw, justifying the notation; i.e., kw is indeed the wall curvature. The gradient \u00E2\u0088\u0087\u00CF\u0086 of a scalar function \u00CF\u0086(n, \u00CE\u00B8) is \u00E2\u0088\u0087\u00CF\u0086 = \u00E2\u0088\u0082\u00CF\u0086 \u00E2\u0088\u0082n e\u00CC\u0082n + 1\u00E2\u0088\u009A G \u00E2\u0088\u0082\u00CF\u0086 \u00E2\u0088\u0082\u00CE\u00B8 e\u00CC\u0082\u00CE\u00B8. (B.13) For a vector function \u00CE\u00A6(n, \u00CE\u00B8) = (\u00CE\u00A6n,\u00CE\u00A6\u00CE\u00B8), the divergence \u00E2\u0088\u0087 \u00C2\u00B7 \u00CE\u00A6 is \u00E2\u0088\u0087 \u00C2\u00B7 \u00CE\u00A6 = 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082n (\u00E2\u0088\u009A G\u00CE\u00A6n ) + 1\u00E2\u0088\u009A G \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B8 (\u00CE\u00A6\u00CE\u00B8) . (B.14) B.3 Coordinate Derivation Given a shape \u00E2\u0084\u00A6, it is a simple task to derive explicit expressions for y(n, \u00CE\u00B8), and hence Gw and kw. Let \u00E2\u0088\u0082\u00E2\u0084\u00A6 be parameterized by a general curve y0(n, \u00CE\u00B8) and consider an interior point (y, z) \u00E2\u0088\u0088 \u00E2\u0084\u00A6. By definition of n, n2 = (y \u00E2\u0088\u0092 y0(\u00CE\u00B8))2 + (z \u00E2\u0088\u0092 z0(\u00CE\u00B8))2, (B.15) 64 B.4. Example: Anisotropic Superellipses for some unknown value of \u00CE\u00B8. In fact, we may consider this as defined for all \u00CE\u00B8, with the added restriction that the distance be minimal: \u00E2\u0088\u0082n/\u00E2\u0088\u0082\u00CE\u00B8 = 0. Then 0 = (y \u00E2\u0088\u0092 y0)y\u00E2\u0080\u00B20 + (z \u00E2\u0088\u0092 z0)z\u00E2\u0080\u00B20, (B.16) with \u00E2\u0080\u00B2 = d/d\u00CE\u00B8. This may be used to define a winding angle tan(\u00CE\u00B1(\u00CE\u00B8)) = \u00E2\u0088\u0092z\u00E2\u0080\u00B20/y\u00E2\u0080\u00B20. We find, simply, y = y0(\u00CE\u00B8) + n sin(\u00CE\u00B1) sgn(tan(\u00CE\u00B1)), z = z0(\u00CE\u00B8) + n cos(\u00CE\u00B1) sgn(tan(\u00CE\u00B1)). (B.17) Using these relations, it is matter of algebra to compute E = 1, G = (y\u00E2\u0080\u00B220 + z \u00E2\u0080\u00B22 0 ) ( 1 + n y\u00E2\u0080\u00B220 (y\u00E2\u0080\u00B220 + z\u00E2\u0080\u00B220 )3/2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3z\u00E2\u0080\u00B20y\u00E2\u0080\u00B20 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0080\u00B2)2 , (B.18) giving explicit expressions for Gw and kw Gw(\u00CE\u00B8) = y\u00E2\u0080\u00B220 + z \u00E2\u0080\u00B22 0 , kw(\u00CE\u00B8) = \u00E2\u0088\u0092y\u00E2\u0080\u00B220 (y\u00E2\u0080\u00B220 + z\u00E2\u0080\u00B220 )3/2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3z\u00E2\u0080\u00B20y\u00E2\u0080\u00B20 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0080\u00B2 . (B.19) It is not yet clear that these coordinates are orthogonal. However, a com- putation shows F \u00E2\u0089\u00A1 \u00E2\u0088\u0082y/\u00E2\u0088\u0082n \u00C2\u00B7 \u00E2\u0088\u0082y/\u00E2\u0088\u0082\u00CE\u00B8 = 0 identically. It should be noted that the preimage of \u00E2\u0084\u00A6 in the (n, \u00CE\u00B8) plane is not, in general, rectangular. This is the price paid for orthogonality. Generically, the domain is the area bounded by \u00CE\u00B8 = 0, \u00CE\u00B8 = 2pi,n = 0, and n = nmax(\u00CE\u00B8), a continuous but perhaps not differentiable function. For example, the preim- age of the half-superellipse appearing at the top of Figure 4.2 is shown in Figure B.3. In some cases, the natural transformation that produces \u00E2\u0084\u00A6 will have folds\u00E2\u0080\u0093 regions in which y(n, \u00CE\u00B8) is not one-to-one. The graph (nmax(\u00CE\u00B8), \u00CE\u00B8) is the locus of points for which, given y1, y1 = y(n, \u00CE\u00B8) has multiple solutions. It is easiest obtained not by algebraic meth- ods, however, but with geometrical intuition. We will illustrate this below, for anisotropic superellipses. B.4 Example: Anisotropic Superellipses In this thesis, we consider the family of anisotropic superellipses (see Figure 4.2). These are shapes parameterized by\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3y0 a \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3my + \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3z0 b \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3mz = 1 \u00E2\u0087\u0094 { y0 = a | cos(\u00CE\u00B8)|2/my sgn(cos(\u00CE\u00B8)) z0 = b | sin(\u00CE\u00B8)|2/mz sgn(sin(\u00CE\u00B8)) . (B.20) 65 B.4. Example: Anisotropic Superellipses Applying the above formulae, we find, after some algebra, y(n, \u00CE\u00B8) = y0(\u00CE\u00B8)\u00E2\u0088\u0092 nbmy \u00CF\u0086(\u00CE\u00B8) | cos(\u00CE\u00B8)|2\u00E2\u0088\u00922/my sgn(cos(\u00CE\u00B8)) z(n, \u00CE\u00B8) = z0(\u00CE\u00B8)\u00E2\u0088\u0092 namz \u00CF\u0086(\u00CE\u00B8) | sin(\u00CE\u00B8)|2\u00E2\u0088\u00922/mz sgn(sin(\u00CE\u00B8)) Gw(\u00CE\u00B8) = ( 2 mymz | sin(\u00CE\u00B8)|2/mz\u00E2\u0088\u00921| cos(\u00CE\u00B8)|2/my\u00E2\u0088\u00921\u00CF\u0086(\u00CE\u00B8) )2 , (B.21) kw(\u00CE\u00B8) = abm2ym 2 z \u00CF\u0086(\u00CE\u00B8)3 | sin(\u00CE\u00B8)|2\u00E2\u0088\u00924/mz | cos(\u00CE\u00B8)|2\u00E2\u0088\u00924/my ( 1\u00E2\u0088\u0092 1 mz cos2(\u00CE\u00B8)\u00E2\u0088\u0092 1 my sin2(\u00CE\u00B8) ) , in which \u00CF\u0086(\u00CE\u00B8)2 = a2m2z| sin(\u00CE\u00B8)|4\u00E2\u0088\u00924/mz + b2m2y| cos(\u00CE\u00B8)|4\u00E2\u0088\u00924/my . 0 1 2 3 4 5 6 7 0 0.2 0.4 0.6 0.8 1 \u00CE\u00B8 n Figure B.3: The preimage of \u00E2\u0084\u00A6 in the (n, \u00CE\u00B8) plane is shown for the case of a semi-superellipse with my = 2,mz = 1. Equivalently, this is the graph of n = nmax(\u00CE\u00B8). The form of nmax(\u00CE\u00B8) depends on my and mz. We will assume a > b, and consider only \u00CE\u00B8 \u00E2\u0088\u0088 (0, pi). Then, for \u00E2\u0080\u009Cfull\u00E2\u0080\u009D anisotropic superellipses, when my = 2 (see bottom right of Figure 4.2), n = nmax \u00E2\u0087\u0094 z = 0, which may easily be solved to find nmax(\u00CE\u00B8) = b\u00CF\u0086(\u00CE\u00B8) amz | sin(\u00CE\u00B8)|4/mz\u00E2\u0088\u00922, |\u00CE\u00B8 \u00E2\u0088\u0092 pi/2| \u00E2\u0089\u00A5 \u00CE\u00B8c (B.22) When my < 2 (see bottom left of Figure 4.2), this is only valid for |\u00CE\u00B8 \u00E2\u0088\u0092 pi/2| >= \u00CE\u00B8c, where the ray at pi/2\u00E2\u0088\u0092 \u00CE\u00B8c goes through the centre of the pipe. For example, for a square pipe, \u00CE\u00B8c = pi/4. In the region |\u00CE\u00B8 \u00E2\u0088\u0092 pi/2| <= \u00CE\u00B8c, n = nmax \u00E2\u0087\u0094 y = 0, which may also easily be solved: nmax(\u00CE\u00B8) = a\u00CF\u0086(\u00CE\u00B8) bmy | cos(\u00CE\u00B8)|4/my\u00E2\u0088\u00922, |\u00CE\u00B8 \u00E2\u0088\u0092 pi/2| \u00E2\u0089\u00A4 \u00CE\u00B8c (B.23) 66 B.4. Example: Anisotropic Superellipses The expression for \u00CE\u00B8c can be obtained from equating (B.22) and (B.23). We will not consider the case my,mz > 2, for which the curve is more complicated. For \u00E2\u0080\u009Csemi\u00E2\u0080\u009D superellipses, when my = 2 (see top of Figure 4.2), clearly n = nmax \u00E2\u0087\u0094 z = nmax. In this case, nmax(\u00CE\u00B8) = b| sin(\u00CE\u00B8)|2/mz 1 + amz \u00CF\u0086(\u00CE\u00B8) | sin(\u00CE\u00B8)|2\u00E2\u0088\u00922/mz . (B.24) 67 Appendix C Velocity Profile in an Infinite Plane Channel In the turbulence literature, there are two canonical examples of pressure- driven flow: the circular pipe, and the infinite plane channel. It is often claimed that, while the boundary layer physics is universal, differences in the velocity profile between the two cases may be expected in the outer region. For a plane channel, kw = 0 and \u00CE\u00BE = 1/2, so that equation (1.3) takes the form \u00E2\u0088\u0092 1 = d dn ( l(n)2 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3du0dn \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 du0dn ) (C.1) Integrating this once yields du0 dn = \u00E2\u0088\u009A 1\u00E2\u0088\u0092 n l(n) , (C.2) which is identical to equation (3.2) for the circular pipe. Since they are also subject to the same boundary conditions, their solutions will coincide. In comparison with data, however, the best fits are produced with slightly different values of the constants \u00CE\u00BA and c4.32 This suggests, first, that to leading order, wall curvature has no effect on the velocity profile, and second, that a refined mixing length model should depend explicitly on \u00CE\u00BE and/or kw. 32See, for example, L\u00E2\u0080\u0099vov et al [13]. 68 Appendix D Viscous Effects Figure D.1: Schematic diagram of three roughness regimes. The wall is at the bottom; the dark line is a nominal \u00E2\u0080\u009Cedge\u00E2\u0080\u009D of the viscous boundary layer, say 5\u00CE\u00BDM/U\u00CF\u0084 . (left) Fully rough regime; (centre) Transitional regime; (right) Hydraulically smooth regime. At high Reynolds number, viscosity affects Reynolds-averaged quanti- ties in two ways, distinct from a modelling point of view. Directly, there are the viscous stresses, negligible in the outer layer, but significant in the boundary layer if the flow is hydraulically smooth or transitional. Indirectly, viscosity affects the Reynolds stresses; in our model, the mixing length must be damped near the wall. To see why, consider a Taylor expansion of the fluctuating velocities at the wall, n = 0. We will abuse notation and let v\u00E2\u0080\u00B2 be the velocity in the wall-normal direction, and w\u00E2\u0080\u00B2 the velocity in the wall-parallel direction. Then u\u00E2\u0080\u00B2 = A0 +A1n+A2n2 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , v\u00E2\u0080\u00B2 = B0 +B1n+B2n2 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , (D.1) w\u00E2\u0080\u00B2 = C0 + C1n+ C2n2 + \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 , where each of the {Ai, Bi, Ci} are functions of time and \u00CE\u00B8, if the flow is fully- developed. By the no-slip condition, A0 = B0 = C0 = 0. The continuity equation is 0 = \u00E2\u0088\u0082 \u00E2\u0088\u0082n (\u00E2\u0088\u009A Gv\u00E2\u0080\u00B2 ) + \u00E2\u0088\u0082w\u00E2\u0080\u00B2 \u00E2\u0088\u0082\u00CE\u00B8 = \u00E2\u0088\u0092 \u00E2\u0088\u009A Gwkwv \u00E2\u0080\u00B2 + \u00E2\u0088\u009A G \u00E2\u0088\u0082v\u00E2\u0080\u00B2 \u00E2\u0088\u0082n + \u00E2\u0088\u0082w\u00E2\u0080\u00B2 \u00E2\u0088\u0082\u00CE\u00B8 . (D.2) 69 Appendix D. Viscous Effects Since v\u00E2\u0080\u00B2 = w\u00E2\u0080\u00B2 = 0 along the wall, which is a contour of n, we also have \u00E2\u0088\u0082w\u00E2\u0080\u00B2/\u00E2\u0088\u0082\u00CE\u00B8|n=0 = 0, and hence \u00E2\u0088\u0082v\u00E2\u0080\u00B2/\u00E2\u0088\u0082n|n=0 = B1 = 0. We may now form the velocity product u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 = A1B2n3 +O(n4), showing the fast decay of the wall- normal shear stress.33 The turbulence model should include this behaviour. This is usually accomplished by empirical wall damping functions [23]. The degree to which the mixing length l should be damped depends on which definition is taken for l. Since \u00E2\u0088\u0082u/\u00E2\u0088\u0082n and \u000FT both tend to a finite constant at the wall, we have \u00CE\u00BDT \u00E2\u0088\u00BC u\u00E2\u0080\u00B2v\u00E2\u0080\u00B2 \u00E2\u0088\u00BC n3 near the wall. The naive definition of l, equation (2.2), then predicts l(n) \u00E2\u0088\u00BC n3/2. Unfortunately, this is inconsistent with the more refined definition (2.4), which predicts l(n) \u00E2\u0088\u00BC n9/4. Although this is not a problem in practice, it is evidence of unphysical behaviour in mixing length phenomenology.34 Since the universal behaviour u \u00E2\u0088\u00BC log(n/\u00CE\u00B4) is observed for all wall- bounded turbulent flows, the transition from a smooth wall to a rough wall can be quantified by the thickness of the boundary layer, \u00CE\u00B4 as a function of the ratio of roughness to viscous length scales, k+s = ksRe\u00CF\u0084 . Empirically, this function is usually presented in a slightly different form, by writing u \u00E2\u0088\u00BC log(n/\u00CE\u00B4v) + \u00CE\u00BA5.6 + \u00CE\u00BA\u00E2\u0088\u0086U+(k+s ). Then \u00E2\u0088\u0086U+ \u00E2\u0086\u0092 0 for a smooth wall and \u00E2\u0088\u0086U+ \u00E2\u0088\u00BC log(k+s ) for a rough wall. The interesting behaviour is in between; this is shown in Figure D.2. To date, there have been no theoretical predic- tions of the function \u00CE\u00B4(k+s ) or, equivalently, \u00CE\u00B4U +(k+s ). We will produce one below by modelling the damping of the mixing length. When the viscous term is included and the eddy viscosity closure adopted, the streamwise momentum equation 1.2 is \u00E2\u0088\u0092 2\u00CE\u00BE = \u00E2\u0088\u0087 \u00C2\u00B7 (( l2|\u00E2\u0088\u0087u|+ 1 \u00CE\u00BARe\u00CF\u0084 ) \u00E2\u0088\u0087u ) . (D.3) After rescaling for the boundary layer, expanding u in \u00CE\u00B4 and taking the leading order equation, this leads to the quadratic equation \u00E2\u0088\u0092 1 = l(n)2 ( du dn )2 + \u00CE\u00B4v du dn , (D.4) 33One may, of course, form the other velocity products, and derived quantities, such as the turbulent kinetic energy and turbulent dissipation. See Pope [23]. 34At the very least, one expects two mixing lengths, as pursued by L\u00E2\u0080\u0099vov et al [13]. In fact, this behaviour persists in the popular k \u00E2\u0088\u0092 \u000F model, which forms a turbulence model with transport equations for the turbulent kinetic energy and turbulent dissipation. As shown by Speziale et al, the known boundary layer behaviour of the turbulent kinetic energy cannot be reproduced without ad-hoc wall functions [30]. 70 Appendix D. Viscous Effects Figure D.2: Law-of-the-wall offset \u00E2\u0088\u0086U+ as a function of k+s , for various experimental data (points) and empiricisms (lines) (from Jimenez [10]). with \u00CE\u00B4v = 1/(\u00CE\u00BARe\u00CF\u0084 ). Solving the quadratic, du dn = 2 ( \u00CE\u00B4v + \u00E2\u0088\u009A \u00CE\u00B42v + 4l(n)2 )\u00E2\u0088\u00921 . (D.5) We will consider a damped mixing length model by directly modelling the right hand side of this equation. A simple form with the correct asymptotic behaviour is 1 2 ( \u00CE\u00B4v + \u00E2\u0088\u009A \u00CE\u00B42v + 4l(n)2 ) = n3 n2 +K2\u00CE\u00B42v(n/\u00CE\u00B42 + 1) + \u00CE\u00B42, (D.6) where \u00CE\u00B42 = \u00CE\u00B4v + c3ks and K is a constant that determines the decay of l. In particular, consider a smooth wall (ks = 0). For n\u001C K\u00CE\u00B4v, we have 1 2 ( \u00CE\u00B4v + \u00E2\u0088\u009A \u00CE\u00B42v + 4l(n)2 ) \u00E2\u0088\u00BC n3 + \u00CE\u00B4v, (D.7) 71 Appendix D. Viscous Effects giving l(n) \u00E2\u0088\u00BC n3/2 as required by (2.2). Conversely, for a fully rough flow with ks/\u00CE\u00B4v \u00E2\u0086\u0092\u00E2\u0088\u009E, this predicts l(n) \u00E2\u0088\u00BC \u00EF\u00A3\u00B1\u00EF\u00A3\u00B2\u00EF\u00A3\u00B3 n+ c3ks n\u001D \u00CE\u00B4v n3 n2 +K2\u00CE\u00B42v + c3ks n \u00E2\u0088\u00BC \u00CE\u00B4v . (D.8) This reduces to the inviscid model at large n, and includes a small damping factor near n = 0. Inserting this mixing length model into the expression for du/dn, we find u(n) = \u00CE\u00B422 \u00CE\u00B423 log ( 1 + n \u00CE\u00B42 ) + K2\u00CE\u00B42v 2\u00CE\u00B423 log ( 1 + n2 K2\u00CE\u00B42v ) + K3\u00CE\u00B43v \u00CE\u00B42\u00CE\u00B423 arctan ( n K\u00CE\u00B4v ) , (D.9) with \u00CE\u00B423 = \u00CE\u00B4 2 2 +K 2\u00CE\u00B42v . As n/\u00CE\u00B43 \u00E2\u0086\u0092\u00E2\u0088\u009E, this implies u(n) \u00E2\u0088\u00BC log(n/\u00CE\u00B4) (D.10) provided log ( \u00CE\u00B4 \u00CE\u00BA\u00CE\u00B4v ) = \u00CE\u00B422 \u00CE\u00B423 log ( \u00CE\u00B42 \u00CE\u00BA\u00CE\u00B4v ) + K2\u00CE\u00B42v 2\u00CE\u00B423 log ( K \u00CE\u00BA ) \u00E2\u0088\u0092 K 3\u00CE\u00B43vpi 2\u00CE\u00B42\u00CE\u00B423 . (D.11) This is a theoretical prediction for the thickness of the boundary layer. To obtain agreement with the smooth wall result, we find, numerically, K \u00E2\u0089\u0088 3. This function is compared to recent data in Figure D.3. Although the model predicts the correct asymptotic behaviour, the fit in the transitional regime is poor. In fact, the naive expression \u00CE\u00B4 = max(c2\u00CE\u00B4v, c3ks), which assumes an infinitely sharp transition from viscous to rough behaviour, fits much better. This suggests that the physics of the near-wall region is more nonlinear than is implied by the simple mixing length model above. 72 Appendix D. Viscous Effects 100 101 102 10\u00E2\u0088\u00921 100 k s + \u00CE\u00B4 Figure D.3: Theoretical prediction of \u00CE\u00B4(k+s ) (solid, magenta) as compared to the naive fit (dashed, black) and the data of Shockling [29] (circles, dark) and Schultz and Flack [28] (circles, light). 73"@en . "Thesis/Dissertation"@en . "2009-11"@en . "10.14288/1.0052882"@en . "eng"@en . "Geophysics"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-NoDerivatives 4.0 International"@en . "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en . "Graduate"@en . "Turbulent flow in geophysical channels"@en . "Text"@en . "http://hdl.handle.net/2429/12802"@en .