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 Eric DeGiuli 2009 Abstract The problem of turbulent flow in a rough pipe of arbitrary shape is considered. 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 mixing length turbulence model is presented, which closes the problem. To calibrate the model, the mechanical problem is solved in the case of a circular pipe. Excellent agreement with engineering relations is obtained. The mechanical 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 application 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 solution 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 solutions 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 . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 11 2 Turbulence Closure . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Circular Pipe Solution . . . . . . . . . . . . . . . . . . . . . . . 16 4 An Approximation Scheme for Non-Circular 4.1 Constant-stress Approximation . . . . . . . . 4.1.1 Square Pipe . . . . . . . . . . . . . . 4.2 Zero-stress Approximation . . . . . . . . . . 4.2.1 Elliptical Pipe . . . . . . . . . . . . . 4.2.2 Semi-Elliptical Pipe . . . . . . . . . . Pipes . . . . . . . . . . . . . . . . . . . . 5 The 5.1 5.2 5.3 . . . . . Thermal Problem . . . . . . . . . . . . . Scaling . . . . . . . . . . . . . . . . . . . . Circular Pipe . . . . . . . . . . . . . . . . . Non-Circular Pipe . . . . . . . . . . . . . . 5.3.1 Constant Heat Flux Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 29 30 34 35 38 . . . . . . . . . . . . . . . . . . . . . . . . . 39 40 44 49 50 iii Table of Contents 5.3.2 Zero Heat Flux Approximation . . . . . . . . . . . . . 6 Conclusion 51 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Appendices A Exact Solution for the Velocity Profile in a Circular Pipe B Geometrical Transformations . . . B.1 Curvature . . . . . . . . . . . . . B.2 Wall Coordinates . . . . . . . . . B.3 Coordinate Derivation . . . . . . . B.4 Example: Anisotropic Superellipses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 61 63 64 65 . . . . . . . . 68 . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 C Velocity Profile in an Infinite Plane Channel D Viscous Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 iv List of Tables 6.1 Glossary of Symbols . . . . . . . . . . . . . . . . . . . . . . . 56 v List of Figures 1 3 Hydraulic radius analogy applied to Blasius’ mean velocity correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sketch of secondary flow in pipes of triangular and rectangular cross-section . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sketch of idealized wall roughness possibilities. . . . . . . . . 3 5 1.1 Sketch of geometry. . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1 3.2 Mixing-length l(n). . . . . . . . . . . . . . . . . . . . . . . . . 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. . . . . . . . . . . . . . . . . . . . . . . . Exact, inner, outer, composite, and law-of-the-wall solutions for the mean velocity profile in a fully rough, circular pipe. . 19 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Sketch of (n, θ) coordinates. . . . . . . . . . . . . . . . . Examples of anisotropic superellipses. . . . . . . . . . . Wall stress for a square pipe; secondary flow velocity. . . Streamwise velocity in a square pipe. . . . . . . . . . . . Wall stress for an elliptical pipe; secondary flow velocity. Streamwise velocity in an elliptical pipe. . . . . . . . . . Wall stress for a semi-elliptical pipe . . . . . . . . . . . . . . . . . . . 22 28 32 33 36 37 38 5.1 5.2 5.3 Stanton number correlation . . . . . . . . . . . . . . . . . . . Heat flux in a square pipe. . . . . . . . . . . . . . . . . . . . . Heat flux in a semi-elliptical pipe. . . . . . . . . . . . . . . . 48 52 53 B.1 Curvature sketch . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Derivation of G(n, θ) . . . . . . . . . . . . . . . . . . . . . . . B.3 Preimage of Ω in (n, θ) plane . . . . . . . . . . . . . . . . . . 62 63 66 D.1 Wall roughness boundary layers . . . . . . . . . . . . . . . . . D.2 Transition from smooth wall to rough wall . . . . . . . . . . . 69 71 2 3.3 . . . . . . . . . . . . . . 3 20 21 vi List of Figures D.3 Boundary layer thickness in the transitional regime. . . . . . 73 vii Acknowledgements It is trite, but true, that one’s first experience in graduate school is a bewildering and difficult, but ultimately satisfying, formative experience. My past thirty months in Vancouver can only be described as such. To all whose paths I’ve crossed in this journey, I owe thanks. It is difficult to separate professional 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– Christian Schoof, Alex Jarosch, Faron Anslow, Christian Reuten, Tim Creyts, Andrew Schaeffer, Tanya Stickford– 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 approach to geophysical problems led to discussions with, and eventually formal 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– I hope the next two years will be just as inspiring. Neil’s 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 computer problems– thanks to Cornel Pop and Marek Labecki for helping me resolve those I’ve 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 – 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 – 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 – 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’s own feelings; one tortures one’s own enthusiasm with doubts; indeed, one experiences even a good conscience as a danger, as if it were a way of wrapping oneself in veils and the exhaustion of subtler honesty – and above all one takes sides, takes sides on principle, against “youth.” – Ten years later one comprehends that all this, too– was still youth. —Friedrich 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 solidification. 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 models 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 “problem of turbulence” is to satisfactorily deal with the loss of information inherent in averaging. To illustrate, consider the Reynolds averaging procedure, whereby flow quantities are considered as sums of mean and fluctuating components, viz., ui = ūi + u0i , P = P + P 0, for velocity and pressure, respectively. When these expressions are substituted into the Navier-Stokes (NS) equations and time averaged, one obtains the Reynolds-averaged Navier-Stokes equations1 [23]: ∂ ūi =0 ∂xi ∂ ūi 1 ∂ P̄ ∂ 2 ūi ∂ 0 0 =− + νM − ūj ui uj . ∂xj ρ ∂xi ∂xj ∂xj ∂xj (1) 1 To perform the time averaging, one assumes that fluctuating components have zeromean, 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 ρu0i u0j , known as Reynolds stresses. The equations are unclosed, 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 particularly 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≡ 2A , Pw where A is the cross-section’s area, and Pw its perimeter. For example, Figure 1 shows the validity of Blasius’ 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 “wetted perimeter,” the total less that of the free surface. The hydraulic radius analogy is very successful in computing mean velocities 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 depends 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 “secondary flow of the second kind”[24].2 In regular geometries, the vorticity forms a pattern of steady, stable rolls which are divided by the lines of symmetry of the crosssectional 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, 2 The “first kind” is streamwise vorticity caused by streamwise curvature of the pipe. 2 Introduction Figure 1: Hydraulic radius analogy applied to Blasius’ mean velocity correlation; 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 streamwise 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 considerations to construct a model based on local equilibrium, after Townsend [31]. Using this “local model” 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 idealized 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 roughness length scale.4 By calibration with the data of Nikuradse [21], this can be identified with an effective sand diameter, in reference to Nikuradse’s classical experiments on pipes with sand glued to their surface. Even with a single roughness length scale, there are many ways of defining the boundary of the idealized smoothed pipe. For a circular pipe, choosing the boundary is a one degree of freedom problem; engineers usually set this “wall offset” 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 3 4 An exception is river meandering, where the secondary flow advects sediment. The 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 simplification 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 distinct 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 ∆P/∆x. It is 5 Introduction conventional to use wall units u+ = ū (R − r) , y+ = , Uτ δv with δv = νM /Uτ a viscous length scale. The velocity scale Uτ is the friction velocity s p R ∆P Uτ ≡ τw /ρ = , 2ρ ∆x where τw 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+ 1 + y = φ y , . (2) dy + y+ R Since the viscous stress is unimportant in the outer layer, one may suppose that du/dy is independent of νM there: y + y φ y , → φO y/R ∼ 1. (3) R R 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 y φ y+, → φI y + y + ∼ 1. (4) R Izakson[9] and Millikan[16] proposed, independently, that there should be an overlap region in which both forms of φ are valid. This is only possible if 1 y → , (5) φ y+, R κ for some constant κ, in the overlap layer. Then (2) can be integrated to yield 1 u+ ∼ 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 κ ≈ 0.421 (the Von Kàrmàn constant) and B ≈ 5 This requires that Uτ be considered as wall stress. p τw /ρ; that is, the velocity scale is set by the 6 Introduction 5.6 [14]. For many applications, (6) is sufficient as a solution to the fluid problem. This is primarily because as δv /R → 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 κ 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 ks+ = ks Uτ /νM .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 φ tends to a constant. To clarify these issues with the heuristic boundary-layer analysis initiated by Prandtl, in the 1950’s applied mathematicians 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, δ. Since the regions are characterized by different physical balances, these equations will necessarily differ in their dependence on δ. The existence of a boundary layer implies that δ 1, so that solutions may be sought as asymptotic expansions in δ. 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 region; the Izakson-Millikan argument produces the leading order behaviour by ensuring that it is possible to match the inviscid outer solution with the 6 κ is observed to be universal, i.e., independent of ks+ . 7 Introduction geometry-independent inner solution. Several matching techniques exist. Introducing n, the dimensionless distance from the wall, the inner region is characterized by n . δ, while the outer region is n ∼ 1. The most general form introduces an intermediate coordinate ζ such that ζ/n → 0 and δ/ζ → 0 as δ → 0. Then, one imposes the condition that, in ζ coordinates, the inner and outer solutions are equal as δ → 0. A simpler way of matching is Van Dyke’s rule, which states, loosely, that in the limit δ → 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 composite 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 δ. This is particularly 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 solutions, and can be quite inaccurate as a velocity profile when δ 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 incomplete. 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 clarifying (6) as an asymptotic result, valid only in the limit δ → 0, and only in a limited range of n. When precise tests are made of the form of the velocity profile, as in the Princeton Superpipe experiment [14], the distinction is crucial. For example, deviations from the logarithmic profile have been interpreted by some researchers as evidence for additional “meso”layers, which, physically, should represent different physical balances. However, these corrections can be more economically interpreted as higher order terms in a standard asymptotic expansion [34]. Not only does this simplify the physical 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̂, of general cross-section Ω, with boundary ∂Ω 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 ∂ui = 0, ∂xi 1 ∂P ∂ 2 ui ∂ 0 0 ∂ui uj =− + νM − ui uj . ∂xj ρ ∂xi ∂xj ∂xj ∂xj (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 − |∆P/∆x| 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 ẑ. 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 τw , which must be defined as (ρνM ∂u/∂y − ρu0 v 0 ) · n̂, in which y = (y, z), u0 v 0 ≡ u0 v 0 ŷ+u0 w0 ẑ, and n̂ 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 ∂Ω. Near the wall, the velocity scale is set by τw . Its mean value τw is known a priori since by integrating (1.1) over Ω, we see that |∂Ω| τw = −|Ω|∆P/∆x. The appearance here of the hydraulic radius R = 2|Ω|/|∂Ω| suggests that it is the relevant length scale for the outer region. In fact, there is another important outer region length scale, more naturally used to scale the tranverse coordinates: the maximum distance from the wall, Nmax . As scales for nondimensionalizing the equations, the difference 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̂, 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 7 For 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 s p R ∆P UT = = τw /ρ ≡ Uτ , 2ρ ∆x the friction velocity. As a mean velocity scale, we will set U ≡ Uτ /κ, where κ . 1 is treated as an empirical constant. For agreement with data, we will find that κ = 0.421, the Von Kàrmàn constant.8 An appropriate scaling for the flow variables in the outer region is then u = U u∗ , (u0 , v 0 , w0 ) = UT (u0∗ , v∗0 , w∗0 ), (y, z) = Nmax (y∗ , z∗ ), P = x = Lx∗ , 1 ∆P LP∗ , 2 ∆x τw = ρUT2 τw∗ , 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 streamwise momentum equation is 1 0 0 0 = 2ξ + ∇ · −u v + ∇u , (1.2) κReτ ∂ ∂ where we have introduced the notation ∇ = ( ∂y , ∂z ) and u0 v 0 ≡ u0 v 0 ŷ + u0 w0 ẑ. The dimensionless constants are Reτ = Uτ Nmax , νM κ= UT , U ξ= Nmax . R For all turbulent flows, Reτ 1/κ, and the viscous term may safely be ignored, in the outer region. However, as mentioned in the introduction, the 8 Values are quoted in the range 0.37 to 0.43 [14]. It is clear from boundary layer considerations 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 , κ = 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, adjacent 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 viscous length scale νM /Uτ . 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 − 2ξ = ∇ · (−u0 v 0 ). (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 = Nmax ks∗ . In the fully rough limit, the dimensionless boundary layer thickness δ will be proportional to ks : δ = cks . It is convenient to choose9 c ≈ 0.0279 so that, for a circular pipe, the law of the wall takes the simple form10 u ∼ log(n/δ). (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. 9 10 0.0279 ≈ exp(−8.5κ) This also requires that we choose κ ≈ 0.421 12 Chapter 2 Turbulence Closure The simplest turbulence closure, due to Boussinesq [3], is the eddy viscosity model − (u0 v 0 ) = νT ∇u, (2.1) where now the eddy viscosity νT 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 νT , (2.1) is very successful [23]. Dimensionally, νT is the product of a velocity scale with a length scale. Adopting the eddy phenomenology of Kolmogorov [11], and making an analogy 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 νT = l2 |∇u| , (2.2) where l is the mixing length. This implies the closure − (u0 v 0 ) = l2 |∇u| ∇u. (2.3) To specify l, Prandtl invoked a dimensional argument in analogy to IzaksonMillikan. There are five length scales in the problem: R, Nmax , n, νM /Uτ , 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 ∝ n. It is easily seen that this prescription implies the law of the wall, so that the overlap 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. 11 We have generalized from the plane channel considered by Prandtl to an arbitrary shape by replacing |du/dy| with |∇u| 13 Chapter 2. Turbulence Closure In the symmetric plane channel case considered by Prandtl, l is a function 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 transverse coordinates. However, the mixing length is interpreted as the size of characteristic eddies – rotating, circular vortices. Since the distance from the wall, n(x), is the radius of the largest circle, centred at x, that is contained 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’vov 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 definition of νT . Kolmogorov [11] proposed that the appropriate velocity √ scale 02 is provided by the turbulent kinetic energy, kT = u , so that νT = κk l kT , with κk a constant. To model kT , the turbulent energy equation is considered under the approximation of local balance [13]. This implies that production of turbulent energy, −u0 v 0 · ∇u, balances dissipation of turbulent 3/2 energy, T . The latter is estimated as T = κ kT /l, with κ a constant; this may be taken as the definition of l. Defined in this way, l is well behaved throughout the entire pipe. Then p lT 1/3 νT = κk l kT = κk l = κν (l4 T )1/3 , κ (2.4) 1/3 with κν = κk /κ . Local balance of turbulent energy implies T = −u0 v 0 · ∇u = νT |∇u|2 . (2.5) Solving (2.4) and (2.5), we find T = κ3/2 l2 |∇u|3 , ν νT = κν3/2 l2 |∇u|. (2.6) 14 Chapter 2. Turbulence Closure To match with known behaviour in the logarithmic region [23], we must take κk ≈ 0.55, κ ≈ 0.553 , so that κν = 1. Therefore, the prescription for νT 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 → 0 as n → 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. 3/2 Using the refined definition of l as κ kT /T , the mixing length is observed 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) ≡ n + c3 ks , 1 + c4 n (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 ! n + c3 ks 2 − 2ξ = ∇ · (2.8) |∇u| ∇u . 1 + c4 n 12 In direct numerical simulation, the Reynolds stress can be computed directly from the fluctuating velocity fields u0 and v 0 , 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 − r. In this case, (1.3) with (2.3) and (2.7) becomes ! 1 d n + c3 /c δ 2 du du −2= (1 − n) . (3.1) 1 − n dn 1 + c4 n dn dn Since variables have been scaled for the outer (fully turbulent) region, we naı̈vely pose an outer solution u = u0 + δu1 + · · · . Inserting this into (3.1) and integrating once, we find the leading order equation du0 du0 (1 − 2n + n2 + C0 )(1 + c4 n)2 = . (3.2) dn dn n2 (1 − n) 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 √ √ 1− 1−n 2 √ u0 = log + 2 1 − n − c4 (1 − n)3/2 + C1 , (3.3) 3 1+ 1−n where the constant C1 needs to be determined by matching with the boundary 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τ . Therefore, an appropriate rescaling in the boundary layer is 1 N = n, δ U = u, (3.4) 16 Chapter 3. Circular Pipe Solution leading to d 1 − 2δ = 1 − δN dN (1 − δN ) N + c3 /c 1 + c4 δN 2 dU dU dN dN ! . (3.5) Letting U = U0 + δU1 + · · · as above, we find C2 = (N + c3 /c)2 dU0 dU0 . dN dN (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 cN +1 , (3.7) U0 = log c3 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’s rule, which states that the outer expansion of the inner solution must equal the inner expansion of the outer solution. The former is cn u0m = log + O(δ) (3.8) c3 δ Expressed in outer variables, the inner expansion of the outer solution is n 2 u0m = log + 2 − c4 + C1 + O(δ log δ). (3.9) 4 3 To obtain agreement to O(δ log δ), we should take C1 = log(4c/c3 δ)−2+ 23 c4 . 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 ≈ 0.0279. This calibrates the wall-offset constant c3 . The appearance of log(δ) in C1 indicates that our original expansion for u was not quite correct: written as an explicit asymptotic series in δ, we have, in the outer region, √ √ 1− 1−n √ u0 = − log(δ) + log + 2( 1 − n − 1) 1+ 1−n 2 3/2 − c4 ((1 − n) − 1) + log(4) . (3.10) 3 17 Chapter 3. Circular Pipe Solution The leading order term, − log(δ), is asymptotically different from 1 as δ → 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≡ = ∆P 4R ∆x ρU 2 hu0 i2Ω 8κ2 . hu0 i2Ω (3.11) (3.12) Computation of hu0 iΩ 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 √ √ δ 1− 1−n 2 √ + 2 1 − n − c4 (1 − n)3/2 + C1 . u0c = log 1 + + log n 3 1+ 1−n (3.13) Using u0c in f , we find Z 2π Z 1 −2 2 1 f = 8κ (1 − n) u0c (n) dn dθ (3.14) π 0 0 −2 16 8c4 4 2c4 2 = 8κ − − + log −2+ + O(δ log δ) (3.15) 15 21 δ 3 −2 1 4 46 6c4 = √ log + + O(δ log δ) . (3.16) − cks 15 21 2 2κ The experimental result, due to Nikuradse [21] and later improved by Schlichting [26], is −2 1 f = 2.00 log10 + 1.74 . (3.17) ks We can obtain agreement to O(δ log δ) if we take13 c4 = 0.6. This completes calibration of the turbulence model. The resulting prescription for the mixing length is shown in Figure 3.1. Analogous results and computations for a plane channel are given in Appendix C. √ Note that log(10)/(2 2κ) ≈ 1.93 for κ = 0.421. Better agreement can be obtained by using κ = 0.406, but it seems preferable to obtain agreement with the newest data; that is, with the Superpipe experiment that yields κ = 0.421 [14, 29]. 13 18 Chapter 3. Circular Pipe Solution 0.8 0.7 0.6 l(n) 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 n Figure 3.1: Mixing-length l(n) (solid) from (2.7) with c3 = 0.0279, c4 = 0.6, ks = 10−3 . For comparison, also shown (dashed) is n + δ. 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 δ = 0.01 (top) and δ = 0.001 (bottom). For δ = 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 δ = 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 4.5 4 2.78 2.77 3.5 2.76 u 3 −0.834 −0.825 10 10 2.5 2 1.5 exact inner outer composite law of the wall 1 0.5 0 −3 10 −2 −1 10 0 10 10 n 7 6 3.854 3.8535 5 3.853 4 10 −1.3381 −1.3372 u 10 3 exact inner outer composite law of the wall 2 1 0 −3 10 −2 10 −1 10 0 10 n 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 δ = 0.01 (top), and δ = 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 4.5 4 3.5 4.076 4.074 3 u 4.072 2.5 4.07 4.068 2 0.594 1.5 0.598 exact inner outer composite law of the wall 1 0.5 0 0.596 0 0.2 0.4 0.6 0.8 1 n 7 6 5 6.245 6.24 u 4 6.235 0.518 3 0.52 exact inner outer composite law of the wall 2 1 0 0 0.2 0.4 0.6 0.8 1 n 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 δ = 0.01 (top), and δ = 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, θ) coordinates. For a non-circular pipe, the streamwise momentum equation − 2ξ = ∇ · (νT ∇u) (4.1) is a partial differential equation (PDE), for which two independent coordinates 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 θ, perpendicular to n.14 The formulae needed to use these coordinates are developed in Appendix B. In these coordinates, (4.1) is 1 ∂ νT ∂u 1 ∂ √ ∂u √ − 2ξ = √ G νT +√ , (4.2) ∂n G ∂n G ∂θ G ∂θ in which G ≡ |∂v/∂θ|2 characterizes the local geometry. For any Ω, it may written as G(n, θ) = Gw (θ) (1 − n kw (θ))2 (4.3) 14 θ is normalized to the interval (0, 2π). Although (n, θ) coordinates are valid in all of Ω, the domain in the (n, θ) 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 approximations 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 = νT /δ. Then, (4.2) becomes √ 1 ∂ δ2 ∂ V ∂U ∂U √ − 2ξδ = √ +√ . (4.4) GV ∂N G ∂N G ∂θ G ∂θ From this we see that, to O(δ), only the wall-normal turbulent flux is important. Posing a solution U = U0 + δU1 + · · · , we find τw (θ) = V with V = (N + 1)2 ∂U0 ∂N ∂U0 ∂N (4.5) + O(δ). The equation τw (θ) = (N + 1)2 ∂U0 ∂U0 ∂N ∂N (4.6) is then the boundary layer equation analogous to (3.6).15 The important difference is that τw is now an unknown function of θ. Equation (4.6) can, of course, be integrated to p U0 = τw (θ) log(N + 1). (4.7) p Recall that U was nondimensionalized with the friction velocity Uτ = τ̄w /ρ. From (4.7) wep see that the appropriate velocity scale is actually the local friction velocity τw (θ)/ρ. Therefore, to leading order, azimuthal dependence 15 Recall 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 ∂U/∂θ, which represents the Reynolds shear stress along θ̂. This simplification turns (4.2) from a partial differential equation into an ordinary differential equation, with parametric θ 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 θ, that is, ∇u ≈ ∂u ∇n, ∂n (4.8) a “radial” approximation. First, the streamwise velocity is empirically known to increase monotonically 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 θ; namely, wall curvature, kw (θ), and the maximum distance from the wall, nmax (θ). 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 “non-radial” 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 ∂u −2ξ = ∇ · νT ∇n ∂n ∂u ên = ∇ · νT ∂n 1 ∂ √ ∂u =√ GνT ∂n G ∂n 1 ∂ ∂u ∂u = (1 − kw (θ)n)l(n)2 . (4.9) 1 − kw (θ)n ∂n ∂n ∂n It is easy to see that we could have obtained this equation directly from (4.11) by neglecting the turbulent shear stress along êθ . Integrating this equation, ∂u ∂u τ (n, θ) 1 + c4 n 2 = , (4.10) ∂n ∂n 1 − kw (θ)n n where we have introduced the reduced stress τ (n, θ) = τw (θ) − 2ξn + ξkw (θ)n2 . In exceptional cases, we may have ∂u/∂n < 0. However, the approximations we will later make require ∂u/∂n ≥ 0. Therefore, Z n p τ (n0 , θ) 1 + c4 n0 p u(n, θ) = C + dn0 . (4.11) n0 1 − kw (θ)n0 0 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 “boundary condition” that determines the unknown function τw (θ) is regularity of the solution in the outer region. In particular, u and ∇u should be smooth throughout the outer region. This is guaranteed everywhere except along the curve n = nmax (θ), where the coordinates are discontinuous. (see Figure 4.2) Along this curve, we must therefore prescribe the jump conditions [u]+ − = 0, ∂u ∂y + = 0, − ∂u ∂z + = 0. (4.12) − The three conditions are not independent, since [u]+ − = 0 implies that the tangential component of ∇u 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 τ . 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 y0 m y z0 m z y0 = a | cos(θ)|2/my sgn(cos(θ)) . (4.13) + =1 ⇔ z0 = b | sin(θ)|2/mz sgn(sin(θ)) a b We will let 0 < my , mz ≤ 2. As special cases, this family includes ellipses (my = mz = 2) and diamonds (my = mz = 1). A plane channel is reproduced in the singular limit mz = 2, b/a → 0. We will let the shape either be “full”, parameterized by θ ∈ (0, 2π), or “semi”, where θ ∈ (0, π), 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 “semi” shapes are reasonable pipe shapes for many geophysical applications [18].16 The formulae for y(n, θ), z(n, θ), Gw (θ), kw (θ), and nmax (θ) 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 (θ). All “full” anisotropic superellipses are of this type. In this case, the jump conditions (4.12) simplify. Continuity of u is guaranteed, and instead of requiring no jump in the u derivative normal to the curve n = nmax (θ), 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 consider n = nmax (θ) as the level set of n−nmax (θ). Then, setting the derivative 16 We have normalized by the maximum distance from the wall. Therefore, for a “full” shape, 1 = min(a, b). 26 Chapter 4. An Approximation Scheme for Non-Circular Pipes normal to the curve to zero, we have17 0 = ∇(n − nmax (θ)) · ∇u 1 dnmax (θ) ∂u ∂u − . = ∂n G dθ ∂θ 17 (4.14) √ Here we use the fact that |∇n| = 1, |∇θ| = 1/ G; see Appendix B. 27 Chapter 4. An Approximation Scheme for Non-Circular Pipes 2 1.5 z 1 0.5 0 −4 −3 −2 −1 0 y 1 2 3 4 2 1.5 z 1 0.5 0 −4 −3 −2 −1 0 y 1 2 3 4 1 1.5 0.8 1 0.6 0.4 0.5 z z 0.2 0 0 −0.2 −0.4 −0.5 −0.6 −1 −1.5 −0.8 −1 −0.5 0 y 0.5 1 −1 −1.5 −1 −0.5 0 y 0.5 1 1.5 Figure 4.2: Examples of anisotropic superellipses. Contours of n (black) and θ (red) are shown, along with the surface n = nmax (θ) (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 profile 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 sophisticated 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., n p u0 = τw (θ) log . (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 θ. Substituting (4.15) into the boundary condition (4.14), after some re-arranging, leads to n −1 d dnmax max 2 log(τw ) = 2Gw (1 − kw nmax ) nmax log , dθ dθ δ (4.16) and therefore determination of τw 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 nonzero 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θ = 0 somewhere, then 1/G du0 /dθ is required to diverge. This severely restricts the applicability of the logarithmic profile approximation. However, for shapes where dnmax /dθ 6= 0 everywhere, (4.16) is often a good approximation. To illustrate, we will derive an explicit expression for τw for a square pipe, and compare this with the DNS data of Huser and Biringen.[8] 18 For 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 θ ∈ (0, π/2). For a square pipe19 , kw = 0 and Gw (θ) = 4 sin2 (2θ). Then for θ < π/4, n = nmax (θ) ⇔ z = 0 ⇔ nmax (θ) = 2 sin2 (θ), (4.17) using (B.22) from Appendix B. Similarly, for θ ∈ (π/4, π/2), we have nmax (θ) = 2 cos2 (θ). Substituting these expressions into (4.16) and integrating, we find nmax (θ) 2 , τw (θ) = τ0 log δ (4.18) where the constant τ0 is set by the requirement that hτw i∂Ω = 1. There is a minor problem with this expression: for nmax < δ 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 /δ by nmax /δ + 1. The approximate solution for u0 is then n √ nmax (θ) u0 (n, θ) = τ0 log + 1 log +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 19 This derivation can easily be generalized to an arbitrary diamond shape (a/b > 1 in 4.13). 20 The DNS of Huser and Biringen was computed for a smooth wall at Reτ = 300, corresponding to δ = 3 10−4 . Since the only impact of the wall on the outer flow should be through δ, we can compare the rough wall theory to data by considering a fully rough pipe at the same value of δ. 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 δ → 0, τw 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 1.8 1.6 1.4 τw(s) 1.2 1 0.8 0.6 zero stress const. stress DNS PDE 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 s 32 Figure 4.3: (top) τw for a square pipe, in terms of arclength s; (bottom) Secondary flow velocity field from DNS of Huser and Biringen.[8] 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 Ω, assuming (4.8), is therefore to let the analogous quantity, the reduced stress τ (n, θ), vanish on n = nmax (θ). It may be seen from (4.10) that this implies ∂u/∂n = 0 on n = nmax (θ), which may or may not be a good approximation to (4.14), depending on Ω. This gives a simple expression for the wall stress, τw (θ) = ξ nmax (θ) (2 − kw (θ)nmax (θ)) , (4.20) which should be contrasted with (4.18). In particular, note that, here, τw has no δ dependence. Also, there are no adjustable constants. Fortunately, Z 2π |∂Ω|hτw (θ)i∂Ω = 2ξ nmax (θ) − ξkw (θ)n2max (θ) ds 0 Z 2π p 1 2 = 2ξ nmax (θ) − kw (θ)nmax (θ) Gw (θ) dθ 2 0 Z 2π Z nmax (θ) p = 2ξ (1 − nkw (θ)) Gw (θ) dn dθ 0 Z 0 2π Z = 2ξ 0 nmax (θ) p G(n, θ) dn dθ 0 = 2ξ|Ω|. In dimensionless terms, ξ = |∂Ω|/(2|Ω|), so that hτw (θ)i∂Ω = 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 curvature, 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: Z np τ (n0 , θ)(1 + c4 n0 ) u(n, θ) = C + dn0 n0 0 ! p √ p τw + τ (n, θ) √ c4 τ (n, θ) p = C + 2 τ (n, θ) 1 − − τw log √ . 3ξ τw − τ (n, θ) (4.21) This matches to the inner solution when √ √ c4 τw 2τw C = −2 τw 1 − − τw log . 3ξ δ 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 (θ) = in terms of φ(θ) = 8a φ(θ)2 φ(θ) , G (θ) = , nmax (θ) = , w 3 φ(θ) 4 2a (4.22) p 4a2 sin2 (θ) + 4 cos2 (θ). Then, after some manipulation, τw (θ) = ξ (1 + 2(a2 − 1) sin2 (θ)) p . a a2 sin2 (θ) + cos2 (θ) (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 τw 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 τw . 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 pronounced. The velocity distributions from the PDE and DNS are compared in Figure 4.6. 21 The DNS was performed for a smooth wall at Reτ = 157, corresponding to δ = 6 × 10−4 . Again, we consider the equivalent fully rough wall. 35 4.2. Zero-stress Approximation 1.2 1 τw(s) 0.8 0.6 0.4 zero stress DNS PDE 0.2 0 2.5 3 3.5 4 4.5 s Figure 4.5: (top) τw 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θ = 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 2 4am2z 2−4/mz 2 2 2− | sin(θ)| cos (θ) − sin (θ) , kw (θ) = φ(θ)3 mz 2| sin(θ)|2/mz (4.24) amz 1+ | sin(θ)|2−2/mz φ(θ) p in terms of φ(θ) = a2 m2z | sin(θ)|4−4/mz + 16 cos2 (θ). The implied formula for τw = ξnmax (2 − kw nmax ) 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 τw , in a similar manner as in the square pipe. nmax (θ) = 1.6 1.4 1.2 τw(s) 1 0.8 0.6 0.4 zero stress PDE 0.2 0 2 4 6 8 10 12 14 16 s Figure 4.7: τw for a semi-elliptical pipe of aspect ratio 2, with linear corners, 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, solutions 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 problem. 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 complicated 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 P r = α/νM 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– at least, in the context of the simplest 22 α 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̄ + T 0 . (5.1) We perform the Reynolds-averaging procedure and drop the bar decoration immediately. The resulting Reynolds-averaged energy equation for an incompressible fluid is cp uj ∂T ∂ 0 0 ∂2T = + T + cp α − cp uj T . ∂xj ∂xj ∂xj ∂xj (5.2) In these equations, and T are the viscous energy and turbulent energy dissipations, given by 1 = νM 2 ∂uj ∂ui + ∂xj ∂xi 2 , 0 ∂u0j 2 ∂ui 1 + T = νM . 2 ∂xj ∂xi (5.3) The dissipation terms are often neglected, but they are important in glaciological applications; we will, therefore, retain them. As a model for T , we will use the simple model23 discussed in Chapter 2, T = l2 |∇u|3 . (5.4) We are interested in the heat flux at the wall, which must be defined as qw ≡ −ρcw α ∂T + ρcw u0n T 0 , ∂n (5.5) in terms of the inward facing normal n. 5.1 Scaling We will introduce separate scales for the mean flow and turbulent fluctuation temperatures. The former is measured from a datum T0 , which we will assume is the temperature at the pipe boundary. To determine the scale 23 It should be noted that this model for T 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 T . 40 5.1. Scaling for turbulent temperature fluctuations, ∆TT , we will assume that the turbulent energy flux balances the turbulent dissipation. The former scales as cp UT ∆TT /Nmax , while the latter scales as24 (κNmax )2 (U/Nmax )3 . Thus 3 ! Nmax U κU 2 2 ∆TT = (κNmax ) = , (5.6) cp UT Nmax cp where we have used κ = 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 0o C. For example, consider a subglacial flood with water running at 10m/s. Then the above scaling predicts ∆TT ∼ 0.01K. The mean flow temperature scale, ∆T , is assumed proportional to the turbulent scale, viz., ∆T = (P rT /κ)∆TT . The turbulent Prandtl number, P rT , is an empirical constant, usually taken in the range 0.7-0.9 [27]. Scaling the new variables T 0 = ∆TT T∗0 , T = T0 + ∆T T∗ , = νM U2 2 Nmax ∗ , T = U UT2 Nmax qw = ρcp UT ∆TT qw∗ T ∗ , (5.7) and dropping asterisks leads to the dimensionless energy equation γP rT ∂T 1 ∂ u = 2 + T − ∇ · v 0 T 0 − γ u0 T 0 2 κ ∂x κ Reτ ∂x 2 P rT 2 2∂ T + ∇ T + γ . P rReτ κ2 ∂x2 (5.8) The new dimensionless numbers are γ= Nmax , L P rT = κ ∆T , ∆TT Pr = α . νM (5.9) We will assume fully rough conditions, so that the terms involving Reτ can be dropped. Then, to complete the model, a turbulence closure for the turbulent heat fluxes must be specified. An expedient model often considered in the literature is the natural extension of the eddy viscosity hypothesis, v 0 T 0 = −νT ∇T, 24 u0 T 0 = −γνT ∂T . ∂x (5.10) Note that l ∼ κNmax . 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 P rT . With these assumptions, (5.8) becomes γP rT ∂T ∂2T 2 u = + ∇ · (ν ∇T ) + γ ν T T T κ2 ∂x ∂x2 (5.11) This is a linear, inhomogeneous problem for the temperature field T . Furthermore, we know that u, νT and T are independent of x. We therefore seek a solution in the form T = T̄ (y, z) + eλx T̂ (y, z), (5.12) which is consistent with the homogeneous boundary conditions on T on ∂Ω.25 The first component, T̄ , contains the effect of dissipative heating on the temperature field, and satisfies ∇ · (νT ∇T̄ ) = −l2 |∇u|3 . (5.13) Evidently, this operates independently of the streamwise advection of heat. By integrating this equation over Ω, we see that the mean wall heat flux due to dissipation is Z −1 hq̄w iΩ = νT ∇u · ∇u dS |∂Ω| Ω Z Z −1 1 = u νT ∇u · n̂ dl + u ∇ · (νT ∇u) dS |∂Ω| ∂Ω |∂Ω| Ω Z −2ξ = u dS |∂Ω| Ω √ −2 2κ = , (5.14) f where f = 8κ2 /huiΩ 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 gradient supplies energy to the mean flow; this is converted to turbulent kinetic energy by straining of the Reynolds stress along velocity gradients; turbulent energy is dissipated into heat; all the heat must escape at the wall as heat flux. 25 To solve a given physical problem, R in general a superposition of solutions is necessary. That is, eλx T̂ must be replaced by Aλ eλx Tˆλ dλ. 42 5.1. Scaling The temperature field due to advection, T̂ , satisfies ∇ · (νT ∇T̂ ) = λγP rT uT̂ + γ 2 λ2 νT T̂ . κ2 (5.15) In the treatment of the velocity problem, we have already implicitly assumed that γ 1. Since we also have κ ≤ 1, the streamwise diffusion term is much smaller than the advection term and may be discarded. If γ κ2 , as we expect, then both terms involving ∂T /∂x may be discarded, and T̃ is unnecessary. However, to facilitate comparison with the dissipation-induced temperature field, we note that T̂ (y, z) satisfies ∇ · (νT ∇T̂ ) = cλ uT̂ , (5.16) with cλ = λγP rT /κ2 . The equations governing νT ∇T̂ and νT ∇T̄ should be contrasted. The former will be largest where u is – namely, near the centre of the pipe – 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 − n)l(n)2 dT̄ du dn dn = −(1 − n)l(n)2 du dn 3 . (5.17) We will use the technique of matched asymptotic expansions. Posing a naive expansion T̄ = T̄0 + δ T̄1 + · · · , (5.18) we find the leading order outer problem ! d (1 − n)3/2 n dT̄0 (1 − n)5/2 (1 + c4 n) =− , dn (1 + c4 n) dn n (5.19) which can be integrated to √ √ 1− r dT̄0 (1 + c4 n) 2 3/2 2 5/2 2 7/2 √ C1 + log = − 3/2 + 2 r + r + r − c4 r . dn 3 5 7 1+ r r n (5.20) Here, we have used r = 1 − n to shorten some of the expressions. Finite heat flux at n = 1 (r = 0) requires that C1 = 0. This can again be integrated, to √ 2 √ 1 1− r 2(1 + c4 ) 1− r √ √ √ T̄0 =C2 − log − log 2 1+ r r 1+ r 16 16c4 2 20c4 12 2 − − log(n) − − r + c4 r2 − c24 r3 . (5.21) 15 7 5 21 35 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 = δN . Denoting the temperature in the boundary layer by TI , we will expand TI as T̄I = T̄I0 + δ T̄I1 + · · · , giving the leading order equation d dT̄I0 1 (N + 1) =− . dN dN N +1 (5.22) (5.23) 44 5.2. Circular Pipe Hence dT̄I0 −q̄w − log(N + 1) = , dN N +1 which can be integrated to T̄I0 = −q̄w log(N + 1) − 1 [log(N + 1)]2 . 2 (5.24) (5.25) Since T = 0 on n = 0, there is no integration constant. Matching with Van Dyke’s rule, as N → ∞, this becomes T̄I0 ∼ −q̄w log n δ − 1 h n i2 + O(δ log(δ)2 ) log 2 δ (5.26) This must be equated to the inner expansion of the outer solution, n 16 16c 1 h n i2 4 T̄0 ∼ C2 − log − − log(n) − 2(1 + c4 ) log 2 4 4 15 7 2 20c4 12 2 − − + c4 − c24 + O(δ). (5.27) 5 21 35 21 This implies q̄w = log(δ) + c6 , C2 = 1 log(δ)2 + c6 log(δ) + c7 , 2 (5.28) where c6 and c7 are δ-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 Ω into n < δ η , and n > δ η , with 0 < η < 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 , 26 c6 = − log(4) − 46/15 − 42c4 /7, c7 = 1/2 log(4)2 − 2(1 + c4 ) log(4) + 2/5 − 20c4 /21 − 12c4 /35 + 2c24 /21. 45 5.2. Circular Pipe respectively, we have, using (5.4), Z δη Z 2π (1 − n)T dθ dn 3 Z δη Z 2π 1 2 (n + δ) ≈ dθ dn n+δ 0 0 Z δη dn = 2π n+δ 0 = 2π log(δ η−1 ) + O(1). QBL = 0 0 (5.29) The contribution from the outer region is Z 1 Z 2π QO = (1 − n)T dθ dn δη 0 1 Z 2π 2 √ 3 n 1 − n(1 + c4 n) ≈ (1 − n) dθ dn 1 + c4 n n δη 0 Z 1 1 + c4 n = 2π (1 − n)5/2 dn n η δ 2 2 2c4 (1 − n)7/2 + (1 − n)5/2 + (1 − n)3/2 = 2π − 7 5 3 √ 1 √ 1− 1−n √ +2 1 − n + log 1 + 1 − n δη η 2c4 2 2 δ η η = 2π + − − − 2 − log + O(δ log δ ) 7 5 3 4 = 2π log(δ −η ) + O(1). (5.30) Z Which contribution is greater depends on the value of η. For 0 < η < 1/2, QBL > QO , while for η > 1/2, QBL < QO . This sensitive dependence on η indicates that most of the dissipation occurs in the matching region. As δ → 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 = −qw /(cp hui(T (n = 1) − T (n = 0))). We will use an overbar 27 Note that St = N ux /(Rex P r), where N ux 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 −ρcp UT ∆TT q̄w ρcp U huiΩ ∆T T (n = 1) p κ f /8 q̄w =− P rT T (n = 1) p κ f /8 − log(δ) − c6 = 1 2 P rT 2 log(δ) + c6 log(δ) + c7 + 4(1 + c4 ) p −1 2κ f /8 c2 − 2c7 − 8(1 + c4 ) log(1/δ) − c6 + 6 = P rT log(δ) + c6 p −1 f 3 /16 c26 − 2c7 − 8(1 + c4 ) = 1− . P rT (log(δ) + c6 )2 ¯ ≡ St (5.31) As δ → 0, the term in parentheses tends to unity, giving a simple form for ¯ St. Although there are accepted engineering relations for St, most are intended 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 dominated temperature fields. Still, it is useful to compare the relations. We will ˆ denote these advection-induced Stanton numbers by St. A simple relation for high Re that follows from the Izakson-Millikan argument, applied to the thermal problem, is [27] ˆ = f /(8P rT ) St (5.32) Engineers favour the Dittus-Boelter equation28 r −1/5 ˆ = 0.023 2 8 Reτ P r3 St . ξ f (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 P r = 1 and replace Reτ with exp(−κB)/δ. 28 This is usually presented as N u = 0.023Re0.8 P r0.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 1 10 0 10 −1 St 10 −2 10 −3 10 Present model Izakson−Millikan Dittus−Boelter −4 10 −6 10 −5 10 −4 10 −3 δ 10 −2 10 −1 10 Figure 5.1: Comparison of Stanton number correlations as a function of boundary layer thickness δ. The three relations are compared in Figure 5.1 with P rT = 0.8, ξ = 1, and P r = 1. All three show similar behaviour, but at small δ, 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 important, we have seen that the temperature fields are additive, viz., T = 48 5.3. Non-Circular Pipe T̄ + eλx T̂ .29 Therefore, the Stanton number is p κ f /8 q̄w + eλx q̂w St = − P rT T̄ (n = 1) + eλx T̂ (n = 1) " ! p κ f /8 q̄w T̄ (n = 1) =− P rT T̄ (n = 1) T̄ (n = 1) + eλx T̂ (n = 1) !# q̂w eλx T̂ (n = 1) + T̂ (n = 1) T̄ (n = 1) + eλx T̂ (n = 1) ! ! λx T̂ (n = 1) T̄ (n = 1) e ¯ ˆ = St + St T̄ (n = 1) + eλx T̂ (n = 1) T̄ (n = 1) + eλx T̂ (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λx , this will depend on the downstream distance. In particular, sufficiently far downstream, if λ > 0, the advection effects should dominate, while if λ < 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. Accordingly, we will follow the same procedure in determining approximate solutions, and omit some details. We will neglect streamwise advection of heat. In (n, θ) coordinates, (5.13) is 1 ∂ √ ∂T 1 ∂ νT ∂T √ √ G νT +√ = −l2 |∇u|3 . (5.35) ∂n G ∂n G ∂θ G ∂θ Again, we attempt solution with matched asymptotic expansions. We will first consider the boundary layer problem. Writing TI = TI0 + δTI1 + · · · for the solution in the inner region, the leading order problem is p ∂ ∂TI0 −τw (θ)3/2 (N + 1) τw (θ) = . (5.36) ∂N ∂N N +1 29 As earlier, for a general problem eλx must be replaced by R Aλ eλx dλ. 49 5.3. Non-Circular Pipe Hence ∂TI0 τw (θ) −qw (θ) − =p log(N + 1), ∂N τw (θ)(N + 1) N + 1 (5.37) which may be integrated to τw (θ) −qw (θ) log(N + 1) − [log(N + 1)]2 . TI0 = p 2 τw (θ) (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 numerical 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 θ derivatives. Then, (5.35) becomes an ordinary differential equation in n with parametric θ dependence. Integrating once, ∂T C(θ) 1 q ≡ −νT =−√ +√ ∂n G G Z n√ G l2 0 ∂u ∂n 3 dn0 . (5.39) Here, q is a reduced heat flux, analogous to the reduced stress of the momentum problem. The unknown function C(θ) 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 (θ), T satisfies 0= ∂T 1 dnmax (θ) ∂T − . ∂n G dθ ∂θ (5.40) The system of equations (5.39) and (5.40) is still too difficult to solve analytically. 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 n 1 −qw T0 = √ log + 1 − u20 . (5.41) τw δ 2 Under these assumptions, (5.40) becomes n ∂ −q 1 dnmax −qw 1 w − log +1 0 =√ √ τw nmax G dθ δ ∂θ τw ∂u0 1 dnmax ∂u0 − u0 . (5.42) − ∂n G dθ ∂θ √ This is a first order, linear equation in qw / τw , whose solution can immediately 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 / τw satisfies the same √ equation as τw . Since the equation admits the trivial scaling symmetry, we have, immediately, √ −qw √ ∝ τw . τw Setting the constant of proportionality so that the correct mean wall heat flux is attained, we have √ −2 2κ qw = τw , (5.43) f a surprising result, considering that the temperature and velocity fields differ. This result allows the solutions for τw from Section (4.1) to be immediately 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 √ ! √ √ τw + τ̃ C(θ) c4 5/2 2 3/2 3/2 √ q = −√ − τ̃ + τ̃ + 2τw τ̃ − τw log √ , (5.44) 5ξ 3 Gw τw − τ̃ 51 5.3. Non-Circular Pipe 1.2 qw(s) / <qw> 1 0.8 0.6 0.4 0.2 PDE Constant Heat Flux 0 0 0.1 0.2 0.3 0.4 0.5 s 0.6 0.7 0.8 0.9 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 τ̃ = τw − 2ξn the approximate reduced wall stress. This must be matched to the outer asymptote of the inner solution (5.37), n qI0 ∼ qw + τw3/2 log . (5.45) δ Hence C c4 5/2 2 3/2 2τw 3/2 3/2 √ = −qw − τw + τw + 2τw − τw log . (5.46) 5ξ 3 ξδ Gw Letting q = 0 when τ̃ = 0, we find30 C = 0 and therefore c4 8 2τw 3/2 qw (θ) = −τw log + τw − . ξδ 5ξ 3 (5.47) 30 This is slightly inconsistent, in that we earlier assumed that τ = 0 on n = nmax . We cannot make that assumption here, because (5.44) requires τ̃ ≥ 0. The best way to proceed would be to retain τ in the dissipation integral. This can be done, also somewhat inconsistently, by ignoring only the effect of wall curvature in the denominator of ∂u/∂n. 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 τw approximation (see Figure 4.7). PDE Zero Heat Flux 2 qw(s) / <qw> 1.5 1 0.5 0 2 4 6 8 10 12 14 16 s 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 solution. 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 θ, 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 mechanical 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 geometries. 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, assuming 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, semiellipses 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 “radial” 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 (θ). 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 Ω ∂Ω u v = (v, w) u0 (v 0 , w0 ) P T T0 νT νM α cp ρ Table 6.1: Glossary of Symbols Pipe cross-section Boundary of Ω Streamwise mean flow velocity Transverse mean flow velocity Streamwise fluctuating velocity Transverse fluctuating velocity Pressure Mean flow temperature Fluctuating temperature Turbulent viscosity Molecular viscosity Thermal diffusivity Heat capacity at constant pressure 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–83, 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–129, 1992. [6] H. Herwig, D. Gloss, and T. Wenterodt. A new approach to understanding and modelling the influence of wall roughness on friction factors for pipe and channel flows. J. Fluid Mech., 613:35–53, 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–95, 1993. [9] A. Izakson. On the formula for the velocity distribution near walls. Tech. Phys. USSR IV, 2:155–162, 1937. [10] J. Jiménez. Turbulent flows over rough walls. Annu. Rev. Fluid Mech., 36:173–196, 2004. [11] A. Kolmogorov. The equations of turbulent motion in an incompressible fluid. Izv. Acad. Sci. USSR; Phys., 6:56–58, 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–461, 1973. 57 Chapter 6. Bibliography [13] V. S. L’vov, 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–147, 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–392. Wiley, 1939. [17] G. Mompean, S. Gavrilakis, L. Machiels, and M. O. Deville. On predicting the turbulence-induced secondary flows using nonlinear k − models. Phys. Fluids, 8:1856–1868, 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. Online: 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–164, 2005. [21] J. Nikuradse. Gesetzmässigkeit der turbulenten strömung in glatten rohren. Forsch. Arb. Ing.-Wes., (356), 1932. [22] R. Panton. Composite asymptotic expansions and scaling wall turbulence. Phil. Trans. R. Soc. A, 365:733–754, 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–25, 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–405, 2007. [29] M. A. Shockling, J. J. Allen, and A. Smits. Roughness effects in turbulent pipe flow. J. Fluid Mech., 564:267–285, 2006. [30] C. Speziale, R. Abid, and E. C. Anderson. Critical evaluation of twoequation models for near-wall turbulence. A.I.A.A., 30(2):324–331, 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–57, 2009. [33] K. S. Yajnik. Asymptotic theory of turbulent shear flows. J. Fluid Mech., 42:411–427, 1970. [34] K. S. Yajnik. Analysis of turbulent pipe and channel flows at moderately large Reynolds number. J. Fluid Mech., 61:23–31, 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 ! n + δ 2 du du 1 d (1 − n) , (A.1) −2= 1 − n dn 1 + c4 n dn dn where we have taken c3 = c as in Chapter 3. Integrating this once yields du du (1 − 2n + n2 + C0 )(1 + c4 n)2 = . dn dn (n + δ)2 (1 − 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 √ √ √ 1+δ− 1−n √ u = 1 + δ (1 − c4 δ) log √ 1+δ+ 1−n √ 2c4 + 2(1 − c4 δ) 1 − n − (1 − n)3/2 + C1 . (A.3) 3 Applying the no-slip condition at n = 0, we find √ √ 1+δ−1 2c4 C1 = − 1 + δ (1 − c4 δ) log √ − 2(1 − c4 δ) + . 3 1+δ+1 (A.4) By expanding in δ, 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 δ, the inner solution can also be obtained. 60 Appendix B Geometrical Transformations Herein we present the transformations necessary to use the (n, θ) coordinate system. A general orthogonal coordinate system (µ, η) defined in the Cartesian y = (y, z) plane is described by a line element ds2 = Edη 2 + Gdµ2 , (B.1) in which G = |∂y/∂µ|2 and E = |∂y/∂η|2 . The unit vectors in the direction of increasing µ and η are 1 ∂y 1 ∂y êµ = √ , êη = √ . (B.2) G ∂µ E ∂η √ The√volume element √ dS = EGdµdη is the area of the rectangle spanned by Gdµêµ and Edηêη . These relations can be used to define the familiar gradient, divergence, and curl operators. In particular, for a scalar function φ(µ, η), the gradient ∇φ is defined by [2] 1 ∂φ 1 ∂φ ∇φ = √ êµ + √ êη . G ∂µ E ∂η (B.3) For a vector function Φ(µ, η) = (Φµ , Φη ), the divergence ∇ · Φ is defined by [2] 1 ∂ √ 1 ∂ √ √ √ ∇·Φ= EΦµ + GΦη . (B.4) EG ∂µ EG ∂η We will not need the curl operator. B.1 Curvature It will be helpful to have expressions for the curvature of the contours of µ and η. To derive these, consider a contour η =constant. êµ is a unit tangent to this curve. (see Figure B.1). Since |êµ |2 = 1, we have êµ · ∂êµ /∂µ = 0, 61 B.1. Curvature Figure B.1: Geometrical setup for curvature derivation. and therefore ∂êµ /∂µ ∝ ê⊥ µ = êη . The curve has curvature kη , defined by √ ∂êµ = kη Gêη . ∂µ (B.5) Since the contour is part of a coordinate system, there are analogous contours of µ at right angles to the contours of η. This fact can be exploited by expanding the identity ∂2y ∂2y = ∂η∂µ ∂µ∂η √ ∂ ∂ √ ⇔ ( Gêµ ) = ( Eêη ) ∂η ∂µ √ √ √ ∂êµ √ ∂êη ∂ G ∂ E ⇔ êµ + G = êη + E ∂η ∂η ∂µ ∂µ Taking the inner product of this expression with êµ , √ ∂êη ∂ G √ = E êµ · ∂η ∂µ √ ∂êµ ∂ = E (êµ · êη ) − êη · ∂µ ∂µ √ = − EG kη . (B.6) (B.7) By symmetry, the analogous expression for kµ , defined by √ we may obtain 31 ∂êη /∂η = Gkµ ê⊥ . Then η √ √ 1 ∂ E −1 ∂ G kη = √ , kµ = √ . (B.8) EG ∂η EG ∂µ 31 Note that the handedness of the coordinates changes under the symmetry transformation. This leads to the sign change. 62 B.2. Wall Coordinates B.2 Wall Coordinates In the particular coordinate system in which η = n is distance from the wall, and µ = θ is a coordinate perpendicular to n, the above expressions simplify considerably. It is well known that n satisfies the Eikonal equation |∇n| = 1. Using (B.3) we see that this implies E ≡ 1. Immediately, kθ = 0, so that contours of θ are straight lines, as expected. This fact may be used to obtain the functional form of G. Let the pipe boundary ∂Ω be parameterized by Figure B.2: Geometrical setup for G derivation. yw (θ) and consider a closed loop Γ along contours of n and θ, as shown in Figure B.2. Since the contours of θ are straight lines, ên is independent of n. Then Z 0= d` Γ Z θ0 +δθ Z n0 Z θ0 +δθ √ Z n0 ∂yw 0 0 0 =− dθ + ên (θ0 )dn + êθ Gdθ − ên (θ0 + δθ)dn0 0 ∂θ θ 0 θ0 0 0 p ∂yw ∂ên (θ0 ) − n0 (θ0 ) + êθ (n0 , θ0 ) G(n0 , θ0 ) δθ + O((δθ)2 ) = − ∂θ ∂θ (B.9) 63 B.3. Coordinate Derivation Taking the limit δθ → 0 and taking the inner product with êθ , we find √ ∂yw ∂θ ∂yw = êθ · ∂θ ∂yw = êθ · ∂θ ∂yw = êθ · ∂θ G = êθ · ∂ên ∂θ ∂êθ − n0 ên · ∂θ √ + n0 Gkn √ ∂ G + n0 ∂n + n0 êθ · (B.10) The first term on the right hand side is independent of θ0 . Since this equation is true for all n0 and θ0 , we may drop the subscripts and consider it as an ordinary differential equation in n. Evaluating p √ this equation at n0 = 0, we see that êθ · ∂yw /∂θ = Gw (θ), the value of G at the wall. The equation can be integrated to p p G(n, θ) = Gw (θ)(1 − kw (θ)n). (B.11) kw is the constant of integration. However, the curvature of contours of n is then √ −1 ∂ G kw (θ) kn (n, θ) = √ = . (B.12) ∂n 1 − kw (θ)n G At the wall, kn = kw , justifying the notation; i.e., kw is indeed the wall curvature. The gradient ∇φ of a scalar function φ(n, θ) is ∇φ = ∂φ 1 ∂φ ên + √ êθ . ∂n G ∂θ (B.13) For a vector function Φ(n, θ) = (Φn , Φθ ), the divergence ∇ · Φ is 1 ∂ √ 1 ∂ ∇·Φ= √ GΦn + √ (Φθ ) . G ∂n G ∂θ B.3 (B.14) Coordinate Derivation Given a shape Ω, it is a simple task to derive explicit expressions for y(n, θ), and hence Gw and kw . Let ∂Ω be parameterized by a general curve y0 (n, θ) and consider an interior point (y, z) ∈ Ω. By definition of n, n2 = (y − y0 (θ))2 + (z − z0 (θ))2 , (B.15) 64 B.4. Example: Anisotropic Superellipses for some unknown value of θ. In fact, we may consider this as defined for all θ, with the added restriction that the distance be minimal: ∂n/∂θ = 0. Then 0 = (y − y0 )y00 + (z − z0 )z00 , (B.16) with 0 = d/dθ. This may be used to define a winding angle tan(α(θ)) = −z00 /y00 . We find, simply, y = y0 (θ) + n sin(α) sgn(tan(α)), z = z0 (θ) + n cos(α) sgn(tan(α)). (B.17) Using these relations, it is matter of algebra to compute E = 1, G = (y002 + z002 ) 1 + n y002 z00 02 02 (y0 + z0 )3/2 y00 0 2 , (B.18) giving explicit expressions for Gw and kw Gw (θ) = y002 + z002 , kw (θ) = −y002 z00 0 . (y002 + z002 )3/2 y00 (B.19) It is not yet clear that these coordinates are orthogonal. However, a computation shows F ≡ ∂y/∂n · ∂y/∂θ = 0 identically. It should be noted that the preimage of Ω in the (n, θ) plane is not, in general, rectangular. This is the price paid for orthogonality. Generically, the domain is the area bounded by θ = 0, θ = 2π,n = 0, and n = nmax (θ), a continuous but perhaps not differentiable function. For example, the preimage 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 Ω will have folds– regions in which y(n, θ) is not one-to-one. The graph (nmax (θ), θ) is the locus of points for which, given y1 , y1 = y(n, θ) has multiple solutions. It is easiest obtained not by algebraic methods, 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 y0 my z0 m z y0 = a | cos(θ)|2/my sgn(cos(θ)) + =1 ⇔ . (B.20) z0 = b | sin(θ)|2/mz sgn(sin(θ)) a b 65 B.4. Example: Anisotropic Superellipses Applying the above formulae, we find, after some algebra, nbmy | cos(θ)|2−2/my sgn(cos(θ)) φ(θ) namz | sin(θ)|2−2/mz sgn(sin(θ)) z(n, θ) = z0 (θ) − φ(θ) 2 2 2/mz −1 2/my −1 | sin(θ)| | cos(θ)| φ(θ) , Gw (θ) = (B.21) my mz abm2y m2z 1 1 2−4/mz 2−4/my 2 2 1− | sin(θ)| | cos(θ)| cos (θ) − sin (θ) , kw (θ) = φ(θ)3 mz my y(n, θ) = y0 (θ) − in which φ(θ)2 = a2 m2z | sin(θ)|4−4/mz + b2 m2y | cos(θ)|4−4/my . 1 0.8 n 0.6 0.4 0.2 0 0 1 2 3 θ 4 5 6 7 Figure B.3: The preimage of Ω in the (n, θ) plane is shown for the case of a semi-superellipse with my = 2,mz = 1. Equivalently, this is the graph of n = nmax (θ). The form of nmax (θ) depends on my and mz . We will assume a > b, and consider only θ ∈ (0, π). Then, for “full” anisotropic superellipses, when my = 2 (see bottom right of Figure 4.2), n = nmax ⇔ z = 0, which may easily be solved to find nmax (θ) = bφ(θ) | sin(θ)|4/mz −2 , amz |θ − π/2| ≥ θc (B.22) When my < 2 (see bottom left of Figure 4.2), this is only valid for |θ − π/2| >= θc , where the ray at π/2 − θc goes through the centre of the pipe. For example, for a square pipe, θc = π/4. In the region |θ − π/2| <= θc , n = nmax ⇔ y = 0, which may also easily be solved: nmax (θ) = aφ(θ) | cos(θ)|4/my −2 , bmy |θ − π/2| ≤ θc (B.23) 66 B.4. Example: Anisotropic Superellipses The expression for θc 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 “semi” superellipses, when my = 2 (see top of Figure 4.2), clearly n = nmax ⇔ z = nmax . In this case, nmax (θ) = b| sin(θ)|2/mz . amz 1+ | sin(θ)|2−2/mz φ(θ) (B.24) 67 Appendix C Velocity Profile in an Infinite Plane Channel In the turbulence literature, there are two canonical examples of pressuredriven 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 ξ = 1/2, so that equation (1.3) takes the form d 2 du0 du0 −1= l(n) (C.1) dn dn dn Integrating this once yields du0 = dn √ 1−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 κ 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 ξ and/or kw . 32 See, for example, L’vov 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 “edge” of the viscous boundary layer, say 5νM /Uτ . (left) Fully rough regime; (centre) Transitional regime; (right) Hydraulically smooth regime. At high Reynolds number, viscosity affects Reynolds-averaged quantities 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 0 be the velocity in the wall-normal direction, and w0 the velocity in the wall-parallel direction. Then u0 = A0 + A1 n + A2 n2 + · · · , v 0 = B0 + B1 n + B2 n2 + · · · , (D.1) w 0 = C0 + C1 n + C2 n 2 + · · · , where each of the {Ai , Bi , Ci } are functions of time and θ, if the flow is fullydeveloped. By the no-slip condition, A0 = B0 = C0 = 0. The continuity equation is 0= p √ ∂v 0 ∂w0 ∂ √ 0 ∂w0 Gv + = − Gw kw v 0 + G + . ∂n ∂θ ∂n ∂θ (D.2) 69 Appendix D. Viscous Effects Since v 0 = w0 = 0 along the wall, which is a contour of n, we also have ∂w0 /∂θ|n=0 = 0, and hence ∂v 0 /∂n|n=0 = B1 = 0. We may now form the velocity product u0 v 0 = A1 B2 n3 + O(n4 ), showing the fast decay of the wallnormal 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 ∂u/∂n and T both tend to a finite constant at the wall, we have νT ∼ u0 v 0 ∼ n3 near the wall. The naive definition of l, equation (2.2), then predicts l(n) ∼ n3/2 . Unfortunately, this is inconsistent with the more refined definition (2.4), which predicts l(n) ∼ 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 ∼ log(n/δ) is observed for all wallbounded turbulent flows, the transition from a smooth wall to a rough wall can be quantified by the thickness of the boundary layer, δ as a function of the ratio of roughness to viscous length scales, ks+ = ks Reτ . Empirically, this function is usually presented in a slightly different form, by writing u ∼ log(n/δv ) + κ5.6 + κ∆U + (ks+ ). Then ∆U + → 0 for a smooth wall and ∆U + ∼ log(ks+ ) for a rough wall. The interesting behaviour is in between; this is shown in Figure D.2. To date, there have been no theoretical predictions of the function δ(ks+ ) or, equivalently, δU + (ks+ ). 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 1 2 − 2ξ = ∇ · l |∇u| + ∇u . (D.3) κReτ After rescaling for the boundary layer, expanding u in δ and taking the leading order equation, this leads to the quadratic equation − 1 = l(n)2 du dn 2 + δv du , dn (D.4) 33 One may, of course, form the other velocity products, and derived quantities, such as the turbulent kinetic energy and turbulent dissipation. See Pope [23]. 34 At the very least, one expects two mixing lengths, as pursued by L’vov et al [13]. In fact, this behaviour persists in the popular k − 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 ∆U + as a function of ks+ , for various experimental data (points) and empiricisms (lines) (from Jimenez [10]). with δv = 1/(κReτ ). Solving the quadratic, −1 p du = 2 δv + δv2 + 4l(n)2 . (D.5) dn 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 p 1 n3 δv + δv2 + 4l(n)2 = 2 + δ2 , (D.6) 2 n + K 2 δv2 (n/δ2 + 1) where δ2 = δv + c3 ks and K is a constant that determines the decay of l. In particular, consider a smooth wall (ks = 0). For n Kδv , we have p 1 δv + δv2 + 4l(n)2 ∼ n3 + δv , (D.7) 2 71 Appendix D. Viscous Effects giving l(n) ∼ n3/2 as required by (2.2). Conversely, for a fully rough flow with ks /δv → ∞, this predicts n δv n + c3 ks n3 l(n) ∼ . (D.8) 2 + c3 ks n ∼ δv 2 2 n + K δv 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 δ22 n K 2 δv2 n2 K 3 δv3 n u(n) = 2 log 1 + + log 1 + 2 2 + arctan , δ2 K δv Kδv δ3 2δ32 δ2 δ32 (D.9) with δ32 = δ22 + K 2 δv2 . As n/δ3 → ∞, this implies u(n) ∼ log(n/δ) (D.10) provided log δ κδv δ2 = 22 log δ3 δ2 κδv K 2 δv2 + log 2δ32 K κ − K 3 δv3 π . 2δ2 δ32 (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 ≈ 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 δ = max(c2 δv , c3 ks ), 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 0 δ 10 −1 10 0 1 10 10 2 10 k+s Figure D.3: Theoretical prediction of δ(ks+ ) (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
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Turbulent flow in geophysical channels
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Turbulent flow in geophysical channels DeGiuli, Eric 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Turbulent flow in geophysical channels |
Creator |
DeGiuli, Eric |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | The problem of turbulent ow in a rough pipe of arbitrary shape is considered. 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 mixing length turbulence model is presented, which closes the problem. To calibrate the model, the mechanical problem is solved in the case of a circular pipe. Excellent agreement with engineering relations is obtained. The mechanical 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 ow. The thermal problem is posed, with scalings taken for intended application 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 solution 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 solutions are discussed. |
Extent | 1964228 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-09-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052882 |
URI | http://hdl.handle.net/2429/12802 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_degiuli_eric.pdf [ 1.87MB ]
- Metadata
- JSON: 24-1.0052882.json
- JSON-LD: 24-1.0052882-ld.json
- RDF/XML (Pretty): 24-1.0052882-rdf.xml
- RDF/JSON: 24-1.0052882-rdf.json
- Turtle: 24-1.0052882-turtle.txt
- N-Triples: 24-1.0052882-rdf-ntriples.txt
- Original Record: 24-1.0052882-source.json
- Full Text
- 24-1.0052882-fulltext.txt
- Citation
- 24-1.0052882.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0052882/manifest