THEORETICAL INVESTIGATION OE THE PARTICLE BOUNDARY LAYER IN PARTICLE-DRIVEN GRAVITY CURRENTS by Elaine Yu-Ling Lin BSc., University of Northern British Columbia, 2002 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in MATHEMATICAL, COMPUTER AND PHYSICAL SCIENCES (MATHEMATICS) THE UNIVERSITY OF NORTHERN BRITISH COLUMBIA August 2004 © Elaine Yu-Ling Lin, 2004 1^1 Library and Archives Canada Bibliothèque et Archives Canada Published Heritage Branch Direction du Patrimoine de l'édition 395 W ellington Street Ottawa ON K 1A 0N 4 Canada 395, rue W ellington Ottawa ON K 1A 0N 4 Canada Your file Votre référence ISBN: 0-494-04698-8 Our file Notre référence ISBN: 0-494-04698-8 NOTICE: The author has granted a non­ exclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distribute and sell theses worldwide, for commercial or non­ commercial purposes, in microform, paper, electronic and/or any other formats. AVIS: L'auteur a accordé une licence non exclusive permettant à la Bibliothèque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par télécommunication ou par l'Internet, prêter, distribuer et vendre des thèses partout dans le monde, à des fins commerciales ou autres, sur support microforme, papier, électronique et/ou autres formats. The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author's permission. L'auteur conserve la propriété du droit d'auteur et des droits moraux qui protège cette thèse. Ni la thèse ni des extraits substantiels de celle-ci ne doivent être imprimés ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conformément à la loi canadienne sur la protection de la vie privée, quelques formulaires secondaires ont été enlevés de cette thèse. While these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n'y aura aucun contenu manquant. Canada Abstract Particle-driven gravity currents are commonly modelled by assuming a ver­ tically well-mixed particle distribution ([4], [9], [10], [11], [17], [21]), but theo­ retical investigation has indicated that the volume fraction does change with depth, particularly for larger sizes of particles, [18]. Experimental concentra­ tion measurements have been made ([1], [7], [8]), and a highly stratified density profile of a low-concentration flow has been observed, especially near the surface of solid bodies such as a river bed. The particle boundary layer (PEL), a new terminology, is defined as the region in which the concentration profile of par­ ticles has a noticeable variation with depth. A closer look at the experimental results showed in [13] and [15], it is theorized that a PEL exists. To fully under­ stand the dynamics of the flow within the PEL, a new mathematical model is constructed to investigate how the particle distribution and scalings in the gov­ erning equations are affected by relaxing the unrealistic vertical independence assumption. A typical thickness of the PEL is estimated algebraically from the particle transport equation, and is then compared with that of a steady state viscous boundary layer. The ratio of the thickness of the two layers justifies the removal of the viscous terms from the Navier-Stokes equations, resulting in a system of first-order partial differential equations. With specified upper layer volume fraction and velocities, theoretical solutions of the system are ob­ tained and discussed in several cases using perturbation theory and method of characteristics. 11 C ontents A b stract ii C ontents iii List o f Figures v List o f Tables v A cknow ledgem ents vi 1 1 2 3 4 Introd u ction 1.1 B a c k g ro u n d .................................................................................................... 1 1.2 Literature S u rv e y .......................................................................................... 8 1.3 Outline of T h e s i s .......................................................................................... 16 Form ulation o f th e M od el E quations 18 2.1 Governing Equations W ithin The P E L ................................................... 19 2.2 Flow within the upper l a y e r ...................................................................... 27 2.3 Dimensional A n a ly sis................................................................................... 28 A nalysis O f T h e P article Transport E quation 38 3.1 Effects of the Particle Transport E q u a ti o n ............................................ 39 3.2 Approximate Thickness of a P E L ............................................................. 40 3.3 Inviscid System Equations for the PEL ................................................ 43 3.4 Analytical Solution of the PEL E q u a tio n ................................................ 44 A nalysis O f T h e N avier-Stokes and M ass C onservation E quations 55 4.1 Solutions of the Velocity, Pressure, and Volume Fraction Equations iii . 56 5 4.2 E x a m p le ............................................................................................................ 63 D iscu ssion 70 5.1 Summary of the Model Assumptions ........................................................ 72 5.2 Limitations of the M o d e l............................................................................... 73 5.3 Identified Problems for Future Research 74 ................................................... A p p en d ix A 76 N o ta tio n 77 R eferences 78 IV List of Figures 1 Visualization of a one-layer fluid with physical quantities defined . . . 19 2 Visualization of a fluid in a flxed re g io n ..................................................... 21 3 Visualization of Pressure Field within Each L a y e r ................................. 30 4 Changes of the PEL Thickness (solid line) and Upper Layer Volume Fraction (dotted line) Using S i ( t ) .............................................................. 5 Changes of the PEL Thickness (solid) and Upper Layer Volume Frac­ tion (dashed) Using % ( ( ) .............................................................................. 6 48 Comparison of the Two Calculated Values for 6 { t ) Using S 3 (dotted) 50 (solid) and ..................................................................................................... 51 7 Ô (solid) and 4>up (dotted) profiles at x = 0 . 5 ........................................... 52 8 Volume Fraction Profiles for small time scales at x = 0 . 5 .................... 64 9 Volume Fraction Profiles for large time scales at x = 0.5 65 10 Volume Fraction Profiles of Silicon Carbide Particles of Sizes 37 iim and 53 /im at t = 1 0 .................................................................................... 66 11 Pressure P r o f ile ............................................................................................... 67 12 Pressure Profiles of po + epi with e = 0.02 68 .............................................. List of Tables 1 The Thickness Ratio of a Viscous Eoundary Layer to a PEL of Silicon Carbide P a r tic le s ........................................................................................... 42 A cknow ledgem ents I would like to express my deep felt gratitude to my supervisor, Dr. Patrick Mont­ gomery, for his great support, encouragement, supervising, and assistance throughout this thesis. His professional comments give me an insight into this research problem. I would also like to thank him for refining the algorithm and improving my writing skills. W ithout his tremendous patience and help, this thesis would not have made it this far. I also want to thank all of my committee members. Dr. Charles Grant Brown, Dr. Lee Keener, and Dr. Margot Mandy, their comments make this thesis to be presented more properly and completely. Last but not least, I especially thank Dr. Peter Jackson for attending my oral defence as the external examiner on short notice. This research was supported by the National Sciences and Engineering Research Council in the form of a Discovery Grant for Dr. Patrick Montgomery and the financial assistances of the University of Northern British Columbia Graduate Stud­ ies (Teaching Assistantships, Scholarship, and graduate student travel grant) and CMS/CAIMS travel awards. VI 1 Introduction 1.1 Background Gravity currents [2] occur when fluid flows primarily laterally into another of different density. The importance of the motion of gravity currents did not attract the attention of many scientists until 1940, when von Karman wrote the first article to analyze the wind speed in response to the American military concerns of released nerve gas back onto friendly troops during War World II [9]. This type of gravity current is called compositionally driven flow, where the driving forces are from the compositional density differences between ambient and intruding fluids, for instance, the release of salty water into fresh water. The excess density providing the driving forces may also be due to differences in temperature, salinity contrast, or suspended sediments. A particle-driven flow, sometimes called a turbidity current [21], is a grav­ ity current in which the density difference is induced principally by solid particles suspended in the current via turbulent diffusion. Particle sedimentation may create a sedimentary layer at the bottom of the flow. Such currents can occur in lakes, seas, as well as oceans and may be initiated by earthquakes or other large mass movements such as the flow of sand, ash, silt, or thick lava, resulting in potential environmental hazards, for example submarine landslides. Cases with homogeneous fluids may be modelled by applying the ideas of similarity solutions to the governing equations ([2],[9],[14]). Although there are few theoretical solutions for gravity currents, and numerical calculations are prevalent for all but several simple cases, one can still borrow some ideas and results from homogeneous flows while constructing mathematical models for particle-driven flow. Similar to particle-driven gravity currents, a buoyantly unstable flow called a particle-bearing gravity current [21 ] is formed in one of the two ways: when the suspensions of parti­ cles are initially or sufficiently dilute through a settling process, or when an interstitial fluid of a released particle suspension has a density less than or equal to th at of the ambient fluid. For example, the intrusion of fresh water rivers containing mud into oceans is a common example seen around the world. Initially, the driving forces of particle-driven flows are solely due to the presence of suspended particles and the current flows underlying the ambient fluid, but, as the majority of particles settle out, the current may stop, or in some cases lift off like a buoyant plume and ascend through the ambient layer. This scenario is termed buoyancy reversal and recently has been studied by Sparks et a /.[29] and Hogg et a /.[10], where the creation of a surface gravity current is considered in which the suspended materials play little or no role in the fluid dynamics. Since the reversing buoyancy currents are indistin­ guishable from compositional flows, one can apply the shallow-water framework to model the dynamics of these flows. In this thesis, the focus is on particle-driven gravity currents where, due to fluidized particles, the buoyancy forces th at drive the flow are continually decreasing as the particles settle out and a transport equation for the particle volume fraction is therefore of key importance in describing the interaction between the flow of the current and the transportation and sedimentation of particles. The local density of these currents strongly depends on the distribution of particles, as do the buoyancy forces within the momentum equations. The concentration of particles also may vary with time and position due to sedimentation and possible re-entrainment of parti­ cles. The dynamics of the current, therefore, are strongly related to both particle transportation and distribution. This thesis aims to investigate this concept through the consideration of a particle boundary layer (PBL). A thorough understanding of turbidity currents is of great relevance as it can be applied in a number of different laboratory, industrial, and natural situations, for instance, pollutant dispersal, vol­ canic hazards, and sewage treatment. To begin to discuss the idea of a particle boundary layer, some terminology is introduced. Fluids are said to be incompressible if their densities are unaffected by changes in pressure within the flows. A homogeneous (constant density) steady flow (all fluid properties do not change in time) is a sufficient condition for incompress­ ibility, but an incompressible flow does not imply a constant density as long as the divergence of the velocity field is zero. Entrainment is a phenomenon in which settled particles re-enter the currents due to the shear-induced turbulent diffusion effect lead­ ing to a self-accelerating current^, considered by Pantin [26] and Parker [27]. Several useful parameters also help the description of fluid flow, and two im portant ones are the Fronde Number and Reynolds Number. The dimensionless Froude number (Pr), a flow property parameter, reflects the ratio of inertial to gravitational forces acting on a fluid flow. For an inviscid homogeneous current moving through a relatively deep ambient fluid, Benjamin found th at Fr % \/2, which differs from the experimental value, Fr = 1.19, since at the nose of the currents, viscous forces are not entirely negligible. The nondimensional Reynolds number, defined as Re = where U is the velocity scale of the flow for a length L and u is the kinematic viscosity, reflects the ratio of inertial to viscous forces acting on a fluid flow and is used in momentum, heat, and mass transfer equations to classify flows th at are dynamically similar. ^The erosion of a sedimentary bed can also cause the same problems. In the initial stages of a gravity current, the dynamics are balanced between in­ ertial and buoyant forces, but as the velocity decreases, viscous forces become more im portant than inertial forces, and eventually they will dominate the dynamics of gravity currents ([4]). For high Reynolds number flows (Re ~ 10^), v is very small and hence viscous forces may be neglected with the consequence th at buoyancy forces may be assumed to be balanced by inertia forces. This measured discarding of viscous effects is generally found to be true far from the boundaries of flow fields and an ideal fluid is therefore generally considered to have viscosity u = 0 (or, equivalently Re = oo). On the other hand, the dominant balance in the flow is between the pressure and viscous forces for a fluid with low Reynolds number (such as the spreading of lava domes). As a result, the flow field of a turbidity current can be divided into an outer region where the flow is assumed to be inviscid and irrotational, and a viscous boundary layer where viscous effects are significant. Recent studies ([9], [10], [18], and [20 ]) have shown th at by ignoring the existence of the thin boundary layer, the flow in the outer region can be well predicted, and the results are generally in fair agreement with published experimental data where these exist. Once the flow in the inviscid (outer) region has been determined, the governing equations describing viscous flows within the boundary layer can then be solved and coupled with the solutions in the outer layer. To predict the motion of a gravity current and its deposited sediment, several models providing either numerical or analytical solutions have been constructed under reasonable assumptions. One of the most common models assumes th at steady fluid flows are inviscid around bodies of various shapes. However, analytical solutions do not always agree with experiments. In 1905, Ludwig Prandtl found th at near bodies of various shapes, a thin layer called boundary layer, exists when the Reynolds number is much larger than 1 so th at the velocity varies rapidly enough for the viscous forces to be im portant and the no-slip condition, i.e. the velocity of the flow is equal to zero at the surface of the solid body, has to be satisfied in this region [14]. In addition, the thickness of the boundary layer changes along the surface of solid body, fluid properties such as velocity components and particle flux are continuous across the layers, and the bulk velocity within the layer reaches its maximum at the edge of the layer. To simplify and solve the governing equations, two models are generally used to describe two dimensional problems: shallow-water and simple box models ([9], [10]). In shallow-water models, the verification is typically compared to experiments in which currents are produced by suddenly releasing a fixed volume of fluid contain­ ing suspensions of particles in water at rest. After a sufficient time, or distance of propagation, the height of the current is seen to be much smaller than its length, and vertical accelerations in the motion are negligible compared with horizontal accelera­ tions. This observation leads to the assumption of a hydrostatic pressure distribution: the pressure at any point in the current is equal to the static pressure due to the surface elevation. A consequence of the hydrostatic assumption is th at for a con­ stant density inviscid fluid the horizontal velocity field is independent of depth. The system of shallow-water equations is formulated by integrating the mass balance equa­ tion over depth z and rewriting the horizontal Navier-Stokes momentum equation as a function of height (h), horizontal velocity (u), and particle volume fraction {

2 %. Temperature may be important in forming gravity currents. dynamics of gravity currents in the upper layer have been investigated extensively, and the theoretical expectations and experimental observations are in fair agreement. 1.2 Literature Survey This subsection contains a survey of several studies in fluid dynamics pertinent to this thesis, such as the validity of assumptions and mathematical techniques used in solving the governing equations. 1) B oundary layer on a flat p late (Fluid Mechanics, by P. J. Kundu (1990) [14]) In chapter 10 of [14] Blasius’ theory of a boundary layer on a flat plate is dis­ cussed. For flow near a boundary, a thin (viscous) boundary layer occurs in which viscous forces are dominant, th at is, when the bulk flow has a low Reynolds num­ ber. To investigate the motion within this region, a term describing the effects of the forces should be taken into account in both momentum equations. For a homo­ geneous fluid with low Reynolds number flowing over a flat plate whose thickness is ignored, Blasius showed th at a similarity solution could be obtained by solving the equations generated from the conservation laws coupled with boundary conditions, for example the no-slip condition, and the continuity of fluid properties at the edge of the boundary layer. Both velocity components are found to increase with height within the boundary layer since the drag forces, due to the viscosity between the flow and the solid body, slow down the velocity near the surface of the plate. This result yields the observation th a t if velocity depends on depth, the distribution of particles in particle driven flow may also be depth dependent. For turbidity currents, this idea of separating the flow into regions is applicable, but instead of dividing it by the strength of the viscous forces, a PBL is a region at which particles are not dis­ 8 tributed equally in depth due to the deposit diffusion. Therefore, in the PBL region, a transport equation describing the advection of particles should be included in the system of the governing equations. Also, to examine the significance of the viscous effects within the PBL, the height of the viscous boundary layer defined in this thesis is compared with the thickness of the PBL. 2) P a rtic le -d riv e n g ra v ity c u rre n ts (Fluid Mech, by R.T. Bonnecaze, H.E. Huppert, and J.R. Lister (1993) [4]) In this study, a two-dimensional, fixed volume gravity current is created through the intrusion of a heavier fluid with particles in suspension into an ambient layer, where the density of the interstitial fluid is equal to th at of the ambient. The monodis­ persed particles are assumed to be dilute, non-cohesive, and with equal settling ve­ locities where they are precipitated only through a viscous sublayer, and vertically well-mixed by strong turbulence induced by the process of sudden release of the intruding fluid. In addition, the flow is assumed to be sufficiently slow such that re-entrainment may be neglected, and the Reynolds number sufficiently large such th a t viscous effects are ignored. After a sufficient, but short, propagation time, a low aspect ratio (height to length) is also assumed and vertical acceleration is neglected when compared to horizontal acceleration such th at a hydrostatic distribution is ob­ tained and the horizontal velocity has no vertical dependence. The above assumptions lead to the shallow-water equations by rewriting the conservation of mass and NavierStokes equations which then couple with appropriate boundary conditions and are solved numerically using the two-step Lax-Wendroff method, and finally, the height, velocity, and density profiles are produced for both one-layer and two-layer flows. Experimental d ata are reported in order to determine lengths and densities of sed­ iments of two-layer gravity currents for different sizes of particles. The experimental results are in good agreement with numerical expectations for small particles, but for larger particles, it is only true initially. Once viscous forces are no longer negligible, which happens at later stages, the results diverge because viscous effects are not in­ cluded in this model. The paper suggests several improvements and modifications. First, this model is also valid for poly disperse particles (particles are adhere to each other). The total deposit density can be redefined as in [9]. Second, when the intersti­ tial fluid is less dense than the ambient, a reversing buoyancy current occurs, and the shallow water equations should be modified in this case. The detailed modification can be found in Appendix B of [10]. The effects of slope and entrainment are also critical to the dynamics and deposition of gravity currents. The main difference of this model to the one constructed in this thesis is th at the volume fraction of particles will have vertical structure; it is not assumed to be vertically well-mixed within the PBL. As such, the standard shallow-water equations are suitable to describe the up­ per layer portion of the flow in the case of a dilute suspension, but the motion within the PBL must be described by including additional terms from the original NavierStokes equations. The experimental setup, such as the initial height of the current and type and sizes of particles, can be used to estimate the relative thickness ratio of the two boundary layers. All of the suggested modifications are also applicable for the new model. 3) A m ath em atical fram ework for th e analysis of particle-driven gravity currents (Proc. R. Soc. Land, by T.C. Harris, A.J. Hogg, and H.E. Huppert (2001) [10]) 10 In this paper, the technique of perturbation series expansions is used to solve one-layer equations for both shallow-water and box models. This technique permits the derivation of analytical solutions to the shallow-water model which are typically integrated numerically. The results indicate the validity of models and explain how gravity currents are different from compositional currents in which similarity solu­ tions exist because the density difference is preserved for a homogeneous fluid, but it does not occur for a turbidity current due to the continuous sedimentation in which the driving buoyancy forces progressively decrease. Similar to the second pa­ per summarized above, the evolution of particle density and the rate of propagation are employed followed by non-dimensionalizing variables with suitable scalings and initial conditions. Once these two equations are reduced to first-order equations, accurate approximations of current length and fluid velocity can be obtained by in­ troducing incomplete gamma functions. Using similar ideas and scalings, the original shallow-water equations are expressed in terms of new spatial and time variables so th a t the approximate solutions from Taylor series expansions (up to three terms) can be gained accurately for the shallow-water model. The solutions are compared with numerical results discussed in the first paper [14] to show an excellent agreement between the two methods when the value of the new time variable is less than the value at which the shock forms. For larger values of the time variable, more terms of Taylor series can be included. The analytical result of the length of currents for both models are very similar by applying this technique, which is why the simplified box model is widely used. W ith slight modiflcation, this technique may be used to analyze other different kind of gravity currents such as modelling fluids in a rotating region. 11 The initial assumptions for the flow are similar to what we considered in the model, except the vertically well-mixed property. The approximate solutions from the box model can be used as boundary conditions at the interface when solving the thickness of the PBL as an example. The ideas of expanding nondimensional equations and dependent variables to calculate the leading-order terms are also applicable for our model. 4) H ydraulic th eory and particle-driven gravity currents (Stud. Appl. Math, by T.B. Moodie (2000) [19]) To study the motions of gravity currents, some scientists apply shallow-water the­ ory by coupling the depth-independent horizontal velocity field with the hydrostatic pressure distribution. In this paper, it is explained th at these two assumptions are not always valid when dealing with fluids consisting of particles in suspension. In this case, the models based on shallow-water theory used to examine the behavior of such flow fields may contain errors questioning results from the past few decades. Due to this impact, the essential details are discussed. For a compositional current, according to the hydrostatic distribution, the pressure can be expressed as: p = Pig{H + r] - h) + p 2 g{h - z) (1) such th at where pi and are the densities of the ambient fluid and the intruding fluid, re­ spectively. Also note th at H is defined as the mean total depth, 77 represents the 12 dispacement of the top free surface from its undisturbed state, and h is the thickness of the gravity current. If the horizontal pressure gradient has no vertical structure, the horizontal velocity is independent from the depth of the current. W ith this con­ clusion, shallow-water theory is applicable for a deep ambient layer together with a zero initial horizontal velocity field in the upper layer. The contradiction arises when particles are present, as in this case p 2 in equation ( 1) should be replaced by: {Pp — P2)4> + P2 - (3) If we take the partial derivative of the equation (3) with respect to x, it is obvious th at the horizontal pressure gradient in the horizontal momentum equation produces # the vertical structure; hence, it is inconsistent to have depth independent horizontal velocity which compliments the research question th at it is impossible to formulate the shallow water equations for such currents. In the second part of this paper, the author justified th at the z-dependent term in the equation of horizontal pressure has a greater effect than the z-independent term by both providing numerical data and explaining the result in mathematical terms. In addition, after a su&cient time, as more and more particles settle down to the bot­ tom of the current, the volume fraction of particles decreases, as does the density of the current. Another influence may be a result of turbulent effects due to the uneven riverbed, for example. When currents hit barriers, they start to climb up and try to pass through the barriers, but some particles may be trapped in the rough surface of the barriers, and the concentration of suspended particles after passing through the barriers is less than before which again reduces the bulk density of the lower layer 13 mixture. As a consequence of these effects, the ratio of the ^-dependent term to the 2-independent term becomes larger demonstrating the fact th a t 2-dependent term dominates the equation of horizontal acceleration of the fluid. The characteristic of a depth dependent horizontal velocity field provides a major paradigm shift from using the traditional model governed by shallow water theory to a new stage as dealing with particle-driven gravity currents. Accordingly, to adopt shallow water theory when modelling the motion within the upper layer, a small ini­ tial volume fraction of particles is essential to be assumed. However, not all scientists working on fluid dynamics accept this new investigation because the research results they obtained by applying standard shallow water equations are reasonable, accept­ able, and applicable to describing and explaining phenomena in the laboratory, and recent research ([18], [21]) has been published based on this new model. 5) Sedim ent tran sp ort and d ep osition from a tw o-layer fluid m odel of gravity currents on sloping b ottom s (Stud. Appl. Math, by T.B Moodie, J.P. Pascal, and G.E. Swaters (1998) [21]) For a homogeneous fluid, a low aspect ratio of height to length guarantees the validity of shallow-water theory such th at the horizontal velocity is independent of depth. When particles are present, it is sufficient to include another key assumption, the smallness of the initial volume fraction, th at is, suspensions are dilute. After appropriate scalings, the authors found th at when the density of the interstitial fluid is greater than th at of the ambient, shallow-water theory is applicable since the mo­ mentum equations and the equation of the particle conservation are decoupled by adopting the two key assumptions. In contrast, when the fluids have the same den­ 14 sity, no m atter how small the aspect ratio or the volume fraction is, shallow-water theory cannot be used due to the presence of the depth-dependent particle concentra­ tion term in the vertical momentum equation. One cannot conclude, therefore, th at the horizontal velocity field is depth independent, which leads to the shallow-water theory. In the numerical investigation, a relaxation scheme is used to solve the sys­ tems of conservation laws. Furthermore, a brief discussion of how the formation of a rear internal bore or hydraulic jump is affected by the variable slope bottom and initial depth of the intruding fluid is also included in this paper. However, this model is based on a two-layer fluid, which is not the case in this thesis. In the upper layer, after defining a different reduced gravity, once the initial volume fraction is small, the vertical derivative of the nonhydrostatic portion of the pressure may be seen to be negligible, leading directly to shallow-water theory regardless of the relative densities of the ambient and interstitial fluids. 6) M od elin g sedim ent d ep osition p attern s arising from suddenly re- lecised fixed-volum e turbulent suspensions ( Stud. Appl. Math, by T.B. Moodie, J.P. Pascal, and J.C. Bowman (2000) [18]) Beginning from the conservation laws and the assumptions leading to shallowwater theory, one of the most im portant conclusions obtained from this paper is th at the volume fraction has vertical structure in contrast to the assumption accepted by most papers th a t the released fluid is vertically well-mixed. The turbulence is not strong enough to keep particles being well mixed through out the entire current. The strength of this dependence is proportional to the particle size. For fine particles, the variation is not very pronounced so th at it is acceptable to assume the particles are vertically well-mixed as in previous studies. For larger particles (with faster settling 15 velocity), a dram atic distribution varying in depth is obtained. This is the starting point of this thesis in which it is assumed th at the concentration profile starts to change in depth within a specific region close to the bottom since the turbulence due to the sudden release of the current has a greater impact near the surface of the flow. The authors took into account the variations of the particle settling velocities near the endwall because the initial turbulent energy may inhibit particle sinking. This is not considered in our case, the settling velocity is assumed to be constant. Although they also included the z —dependent term in the volume fraction as we do, but to work in the realm of shallow water theory, they specified a particular function for the particle volume fraction such th at the concentration can only increase with depth while our model allows it to decrease. This paper also helps in formulating equations by introducing very similar scalings to simplify and nondimensionalize the governing equations. 1.3 O utline of Thesis The outline for the remainder of this thesis is as follows. In the second chapter, the general governing equations are formulated according to the ideas of conservation laws, non-dimensionalized with proper scalings, and an equation for particle trans­ port is derived and used to calculate the height of the PBL. The systems of equations for the upper layer and PBL are summarized subject to several boundary conditions at the end of the chapter, and it is noted th at the second order PDEs may be further simplified to first order if viscosity is small. To examine the influences of viscous effects, in Chapter 3 the relative thickness of the viscous boundary layer to the PBL is considered by approximating the thickness of the PBL from the transport equation and then taking the ratio of the two boundary layers. The resulting governing equa­ 16 tions are rewritten and the transport equation is solved as an example with specified upper layer volume fraction and velocities. Once the thickness of the PBL is ob­ tained, the remaining equations are solved in Chapter 4 using a perturbation method and the method of characteristics. But, to successfully adopt the volume fraction and velocities of the upper layer from the existing papers ([10],[11]), we must assume th a t the fluid is dilute such th at shallow-water theory is consistent when modelling the upper layer flow. Finally, solutions are plotted using Maple for several cases and examples. 17 2 Formulation of th e M odel Equations In this chapter, a mathematical model is constructed from the basic physical principles of mass conservation, momentum, and particle transportation. A timedependent system of partial differential equations describing the motion within the particle boundary layer (PBL) for one-layer particle-driven gravity currents in two spatial dimensions is derived. The general cross-section of a gravity current, ac­ cording to the definition of the PBL, is divided into two regions: the outer layer and the PBL. Several models ([4],[9],[10],[11],[18],[21]) have been proposed to pre­ dict the dynamics within the outer region, in which the volume fraction of particles is z-independent by adopting shallow-water theory, and the mathematical analyses compare fairly well with experimental observation ([9],[10],[11]). By using these avail­ able results ([10],[11]) for the outer layer problem as boundary conditions for the PBL, this chapter will concentrate on establishing the model equations for the PBL. To reduce the complexity of the model, several standard assumptions are made initially based on physical considerations and established practice, for example ne­ glecting surface tension effects and the Coriolis force. The resulting equations are expressed as a system of five equations with five dependent variables. The physical variables are defined in the first section of this chapter followed by a derivation of the conservation of particle transport in the PBL. The thickness of the PBL is to be calculated based on this equation. The second section summarizes the dimensionless governing equations for the upper layer, which are quoted from [17]. Finally, the model equations formulated in the first section for the PBL are nondimensionalized in the last part of this chapter. Dimensional analysis reveals less significant terms 18 from equations, and by removing those terms, a simplified system is obtained and summarized. 2.1 G overning E quations W ith in T he PEL Consider a turbidity current produced by instantaneously releasing a fixed-volume suspension of bulk density p with relatively heavy particles of density pp into a deep ambient fluid of density pa- Let Uup, Wupi 4>i ^ be the volume fraction of parti­ cles, horizontal, and vertical velocities for the upper region and the PEL, respectively. All of the PEL variables are functions of x, z, t only, with y removed by a symmetry consideration. Also, denote the height of the current as h{x,t), and the thickness of PEL as S(x,t) . The physical configuration of the variables is depicted in Figure 1. Deep Ambient Pa Turbidity Current Upper Layer PEL i|)(x,z,t) Figure 1; Visualization of a one-layer fluid with physical quantities defined In general, the total rate of change in mass for a fluid with density p and velocity u = {u, V, w) is described by the mass balance equation [14] dt + V • {pu) = 0. 19 (4) A simplifying assumption typically made [14] is to assume th at the fluid is incom­ pressible, or V • = 0. Using this, the above equation can be simplifled to In the upper layer and the PEL, p = Pp(j) + pi{l — 4) where pi is the density of the background fluid carrying the particles. In the ambient layer, p = Pa- Two more equations representing conservation of momentum for an incom­ pressible fluid in the horizontal and vertical directions are du du du 1 dp ,d‘^u d“ ^u , . and dw dw dw Idp ,d'^w d'^w. where p is the pressure and u is the kinematic viscosity [14]. Equations (5) — (7) and fluid incompressibility represent four equations for, in essence, m, w , p, and (j). Another equation is needed to close the system. This is now derived, as the equation is not a standard result available in most texts concerning fluid dynamics. To formulate the transport equation, we assert th at the source of particles is only from the upper layer due to the sedimentation and these particles are not leav­ ing the PEL: they accumulate at the bottom of the PEL. Consider a fluid layer of height 5{x,t) flowing through a region fixed between Xi and x^ (see Figure 2). At a 20 Input Source z = 5(x,t) Inflow Mass Outflow Mass Figure 2: Visualization of a fluid in a flxed region given time ti, the mass inside the region Xi < x < nr X 2 rS d(x,ti) / / Jxi Jo and 0 < z < 6 is given by p{x, z ,t i) dzdx. ( 8) Then, the total mass between t\ < t <1^ within this region is ÇX2 pS{x,t2) / / p(x, z , t 2 ) dzdx — / / Jx\ Jo Jxi Jo PX2 P pO S {( xX , t l ) p(x, z , t i ) dzdx. (9) The change of mass due to the inflow mass flux across the vertical boundary at Xi is Pt2 ps{xi,t) / / p { x i ,z ,t ) u { x i ,z ,t ) dzdt, Jti Jo ( 10) and th at of the outflow mass flux at X2 is Pt2 - p S { x 2 ,t) / Ju Jo p{x 2 ,Z,t)u{x 2 ,Z,t) dzdt. 21 ( 11) Since particles settle out of the upper layer and enter into the PEL, the input source is defined as pt ft 2 rX2 / / Pp4>upVs dxdt, Jh ' t l Jx\ Jxi where Vg = ^ aipp-Pa) jg (12) Stokes free-fall velocity [9] (which depends on particle size, particle and fluid densities, and fluid viscosity), a is a representative length scale of the particle, and g = 9.8 m /s^ is the gravitational acceleration. Using the principle of conservation of mass, the total mass inside the region must equal the sum of the inflow mass, outflow mass, and input source. This is given as P X 22 p5{x,t2) r0{x,t2) / rX / Jxi Jo rS(xi,t) / / J t\ J0 J pt2 pX j ‘X 22 p { x , z , t 2 ) dzdx — pS{x,tl) rO (x,tl) p { x , z , ti) dzdx = / Jxi Jo rt2 p { x i , z , t ) u { x i , z , t ) dzdt - r5{x2,t) / p{x 2 , z , t ) u { x 2 , z , t ) dzdt + J t \ J (J pX2 / (13) P p (j)u p V sd xd t. Jx\ Combining terms, the previous equation gives X2 / / ro[x, t 2 ) {j rt2 / / Jti rS(xi,t) I / \ Jo + ro{x, t i) p{x, z , t 2 ) dz - J \ p {x,z,ti ) dz ] dx = rS { x 2 ,t) p { x i , z , t ) u { x i , z , t ) dz - J p{x 2 , z , t) u { x 2 , z , t ) dz Jo Pp(f>upVsdx'^dt. (14) 'Xl Here, we must assume th at all of the functions are differentiable. The first integrand 22 in the first line of equation (14) can be written according to the following steps: ro(x,t rS{x,t2) 2) pro(x,ti) 5 {x,ti) p{x,zM )dz - / Jo Jo r5{x,h) = p{x,z,ti)dz / rS{x,t 2 ) p{x,z,h)dz + / / Jo r&{x,t\) p{x,z,t2 )d z - J 8{x,t\) Jo r&(x,ti) = rSixM) {p{x,z,t 2 ) - p{x, z, ti)) d z + / p{x,z,ti)dz / p{x,z,t2 ) d z Jo J 5{x,ti) ç 5{xM) = Jo p5{x,tl) = y / 0 { —p {x,z,t)d t]d z+ / \ Jti Ot Js(x,ti) ft2 \ rSixM) J p{x,zM )dz y f) ^ ’ - dtdz + p{x,5*,t2){5{x,t2) - 5{x,ti)), y where 5* = 5{x, t*) and ti < t * < t 2 , by the Mean Value Theorem for Integrals [30]. The previous calculation may be continued as f (i -Q^dz + p{x,ô*,t2)-Ql ) dt. (15) P art of the right hand side of equation (14) may be rearranged similarly. This follows: rô{xi,t) ro{xi,t) / Jo fS /•0{x2,t) {X2 ,t) p { x i , Z , t ) u { x i , Z , t ) d z / Jo - p { x 2 , Z , t ) u { x 2 , Z , t ) d z f‘S{x2,t) = / Jo l>S{x2,t) { p { x i , z , t ) u { x i , z , t ) rS{x2,t) = -J J - p { x 2 , Z , t ) u { x 2 , Z , t ) ) d z - I-X2 Q rX2 — {pu)dxdz - p{xi,5**,t)u{xi,5'^,t) j / Js(xi,t) p { x i , Z , t ) u { x i , Z , t ) d z dx ■f(i (16) 23 where S** = 6{x**,t) and Xi < x** < X2 . Substituting the results (15) and (16) into the equation (14) gives IL [I rt2 PX2 *2 / / rX2 J rà(x,ti) g / g g \ % rS{x2,t) \ J = g d ô \ J g^ipu) dz + p{xi,5**,t)u{xi,5**,t)— j d x d t + J Pp(f)upVs dxdt. (17) These terms may all be combined as one double integral which is zero for an arbitrary region [xi,X 2 ] x [ti,t 2 \. Therefore, since the integrand is continuous, it follows th at the integrand vanishes as well. In the limit as tg, ti t, X2 , Xi —»x, and 5*, 6** S(x, t), dz + /)(a;, Ô, t)'u(a;, Ô, t) — — pp4>upVs. eqn (17) gives J at of J — This may be simplified to J lâ f ~d^ J p(z,(^,f)4z,<^,()— = Pp(f>upVs(18) Equation (18) holds in the PEL, but it is general and may also be used for the upper layer. The same concepts can be applied to derive an equation for the conser­ vation of particles mentioned in [17]. T hat is, if Ô and p are replaced by h and Pp4> respectively, and coupled with the following equality g — J rh{x,t) Pp4)udz = J Ph{x,t) g gi^ ^{pp(j)u) dz -h pp(f)u{x,h,t)— , 24 (19) which gives I - d T - &X '“’*)âï- Equation (18) becomes ^ J Pp4> d z — Pp^~Q^ + ^ y ppÇu d z — Pp(f)u{x,h,t) — ,dh 1 / u + P p ^ ~ g ^ + P p 9 'fJ‘{ x , h , t j — , — — ppÇupVs Simplify the terms to get nyx,t) rh{x,t) Q / Q Pp(j) rn (x,t) rh(x,t) dz + — J Pp(j)U d z = -p p (t}u p V s- W ith the ^^-independent integrands, this equation can be further simplified to Q Q rn (x ,t) rh{x,t) — (Pp0h.) + ^ y Pp4>u d z = — P p 4>up'^s- The negative sign of Pp4>upVs indicates the direction of particle sedimentation, and the above equation may be observed to have a similar form as the equation (3.5) on p.68 of [17], validating eqn (18). Returning to the PEL, substitution of mass conservation (4) into (18) yields J Q d5 - — {pw)dz + p { x ,5 ,t)— -f 25 85 p { x , 5 , t ) u { x , 5 , t ) — = Pp(f>upVs, or p ( x , 0 , t ) w ( x , 0 , t ) - p ( x , 6 , t ) w ( x , S , t ) + p ( x , ô , t ) ^ + p ( x , â , t ) u ( x , 6 , t ) ^ = Pp<^>upXn ( 20 ) Therefore, at the interface, since w(x, 0,t) = 0 is a physical boundary condition. + o t O X ( 21 ) = p and p = Pup, u = Uup, w = Wup due to the continuity of these variables at the interface. In summary, the dimensional equations for the PEL are: du dw Ï • -E • -S du du Idp dw dw Idp - « d'^u , d ‘^w d^u d‘^w. and f + u{x,S ,t)^ = at dz Pup Note th at eqn (51) is the fluid incompressibility. 26 + w{x,S,t). (26) 2.2 Flow w ithin th e upper layer W ithin the upper layer, particles are assumed to be vertically well-mixed and drift down into the PEL during the flow and no re-entrainment occurs at the interface. The non-dimensional governing equations for the upper layer are a special case of (22) - (26). They are typically stated [17] as 1+ 9pup . dx d(j)np dz J rh (4>uph) + (^up Uupd^ + 4>up-0 — 0, (30) + U up{x,h,t)^, (31) and Wup{x,h,t) = ^ where g' — All the variables describing the upper layer flow, denoted with the subscript ‘up’, are functions of x and t only. Note th at since the PEL requires th a t the upper layer lies in the range of ^ < z < h, eqn (30) is not precisely valid. However since we are focusing on the flow within the PEL, a correct restatement of eqn (30) is not needed. An approximation of the typical thickness of the PEL to be discussed in Chapter 3 shows th a t the upper layer lies far above the viscous boundary layer, and hence viscosity is excluded from eqns (27) — (31). In addition, if the current is dilute in suspension, eqns (27) — (31) can be reduced to the standard shallow-water equations 27 as stated in [4] and [11]. 2.3 D im ensional A nalysis Nondimensionalizing equations is a process used to eliminate variable units for two main purposes. Experiments are easier to perform by perm itting the design of a smaller model than the actual size of the problem, and a standard solution domain aids numerical calculation. In addition, nondimensionalization helps to find the rel­ ative significance of terms in the governing equations. It is often the case th at small parameters may be identified, and the equations may be simplified by neglecting these small terms. First, the nondimensional variables with star notation are introduced as follows. The length variables are defined as X = x*L , 0 < X* < 1, z = z 'jy , 0 < z* < 1. and Similarly, the time, volume fraction, horizontal velocity, vertical velocity, particle density, Reynolds number, and Stokes free fall velocity are defined as t = t * T , t *>0 (j) — is the general equation for the horizontal momentum. The choice of pressure scale P = is made following the convention taken by Moodie [19]. As i?e —> 00 , —> 0, consequently, the horizontal momentum equation is simplified to g J \d t dx dz J Furthermore, if fiuid is dilute, U'^ dx \ g J dz^ (33) < < 0 (1) is negligible. For example, mud or clay has ^ ~ 1.6, [6], and if e = 2%, then du dx du du = 0.032. In this case eqn (33) is reduced to dp g'H d(j> d'^u Note th a t e, U, H are still variables. V ertical M om entum The pressure field has been defined as p = Pa~ pgz + p, then in general dp dp 32 dp , , Substituting this result into equation (24) gives dp / dw dw dw' P I -ITT + u —— h wdt dx dp'' / d^w d^w' Dividing both sides of the above equation by pi, write I dp 1 — — I — dz Pi dp dz + 1+ in a manner similar to the horizontal momentum equation. Introducing the nondi­ mensional variables gives Wdw* 9' + C/W .dw* 1 -u dx* H L .dw* ' ----dz* H dz' p, H dz* + Now, dividing the above equation by multiplying it by ^ , and dropping * notation to get the general nondimensional momentum equation as 1 9' f ^ dw dw \ L2 dp + g'H d(t) 1 d^w Re dx"^ + Under the assumption of a small aspect ratio, 1 d^w Re dz^ . (35) and large Reynolds number, equation (35) is simplified to (36) 33 M ass C onservation Substituting p = {Pp — Pi)(j) + p% into equation (22) gives d(f) , 0 0 ^0 Introducing the nondimensional variables into the above equation followed by dividing it by ^ and dropping the * notation, dt ox (37) dz is obtained. P article Transport At the interface of the upper layer and the PBL, Uup = u and w^p = w. T hat is, u* = and w* = when expressed in nondimensional variables, where Uup and Wup are velocity scales of the upper layer flow. Let 0„p = ^up4>up- Equation (26) gives T at* + u* L dt* Fix a = Pp^up4*up^^s (pp - Pi)€up0* + A L az* -u + dx* Pp^up4*up^^s {Pp — A*)Gup0up + p, drop * notation, divide the above equation by ^ ^ dt ^ U ^^^dx P p 4 >upVs_____________ 1^ up _ {pp - P i ) e u p 4 >up + P i is formulated. 34 W so th at w.up (38) If dividing the numerator and denominator of the first term on the right hand side of eqn (38) by pj, ^eup4>up can be eliminated assuming a dilute suspension, and the particle transport equation is given as 95 Uup dô Pp4^up'^a I^up m + ■ £ r “ ”>’â ï + I V ” "'" (39) Note th at W depends on [/, L, and H and is not fixed independently. Boundary conditions The boundary condition of pressure at the interface of the upper layer and the ambient layer is Pupis^i h, t) Replacing p ^ „ by U ^ ^ p t p l ^ , by Pupph. + P i ( l - e.pÿ;,) and k by H ^ f h ' gives Pupip^i h, t) — h (40) T ^ after dropping the * notation, where Hup is the height scale for the upper layer. Similarly, at the interface of the upper layer and the PBL, p{x, 5, t) = Pup9 {h - 5) + pgS. Therefore, the pressure condition at this specific height is U2 SJ) = 1^1 + ^e„p0 upj 35 Hu (41) since pup = p ai z = ô. Finally, two more nondimensional boundary conditions for the vertical velocity are the kinematic conditions of no net flow at z = 0, w{x,0,t) = 0, (42) and continuity at the interface, — w{x, S, t) = S, t). (43) W ith the initial assumptions of low Reynolds number and small aspect ratio, the nondimensional governing equations (32), (33), (36), (37), and (38) of the PBL coupled with boundary conditions (40) - (43) describe the motion of a particle-driven gravity current flowing under a deep ambient fluid. Equations (27) — (31) have been well examined in [18] and, therefore, are adopted to describe the dynamics in the upper layer. These upper layer solutions are used to match with th at of the PBL due to the continuity of fluid properties. Ideally, the upper layer and the PBL problems should be solved simultaneously to match the conditions at the interface. But since the upper layer motions have been well investigated, to focus on the dynamics in the PBL, all upper layer fluid variables are assumed to be known and used as boundary conditions when seeking PBL solutions. Note th at two height scales are used for different layers: H = for the PBL but Hup — ^ the heights are extremely different, H « for the upper layer since the magnitudes of Hup due to the presence of small viscosity, I 0 ~^rn?/s for water at 20°C. Later in Chapter 3, examining the heights of the two boundary layers shows th a t the viscous forces are in fact negligible which redefines the H in the same manner as Hup. To avoid confusion, the same height scale is used 36 for both layers, resulting in velocity scales U and W . Moreover, further simplified version of horizontal momentum and particle transport equations are stated when the current has small initial volume fraction. W ith this assumption, the governing equations of the upper layer flow are simplified to standard shallow-water equations, and results from published papers ([10], [11]) are borrowed to satisfy the continuity property at the interface when solving the PBL problem (to be discussed in Chapter 3 and 4). 37 3 A nalysis O f The Particle Transport Equation The focus of this chapter is the impact of the particle transport equation (21), on turbidity current models. The relative thickness of the PBL to the viscous boundary layer is of importance in determining not only the validity of the model used in the upper layer, but also to permit further simplification of the model equations. For example, if the viscous boundary layer is thick compared to the PBL, the shallowwater model used in the upper layer may not fully describe its dynamics, and viscous forces should be taken into account. On the other hand, if the ratio is relatively small, viscous effects may be neglected and one can apply the results obtained from the shallow-water model as boundary conditions required to solve the PBL problem. A short first section contains some physically sensible effects of the particle trans­ port equation. The second part of this chapter develops a basic algebraic approxi­ mation to the typical thickness of a particle boundary layer. Adopting the average height of the viscous boundary layer quoted from [14], the ratio of the two boundary layers is estimated. To calculate this ratio, scales of height, velocity, and length of a flow are specified using the experimental set up data from [4]. A table illustrates the results for silicon carbide in various sizes of particles. The ratio reveals the insignif­ icance of viscous effects played in Navier-Stokes equations compared with the effect of the particles and permits simplification of the second-order horizontal momentum equation (33) to a first-order PDF. As a result, the dimensional governing equations and boundary conditions are rescaled using the same height and velocity scalings as the upper layer. 38 The thickness of the PBL, relying only on the fluid variables of the upper layer and some physical properties such as the density and settling velocity of particles, is solved using the method of characteristics when the upper layer velocity (uup) and volume fraction are specified. Solutions are discussed for two cases: Uup = constant and Uup = Uup{x, t) in § 3.4. Two examples are given to illustrate the variation of 5 for the box model [10] and shallow-water model [11]. 3.1 Effects o f th e Particle Transport Equation This section investigates physically sensible effects of the particle transport equa­ tion (21). As this equation is not usually used in this situation, and hasn’t been derived in any papers cited in the references, it deserves a brief investigation. A direct observation of the equation = (44) tells th at the dynamics inside the PBL do not affect its height due to its derivation by vertical integration. However, the upper layer velocity and the upper layer volume fraction have a great impact on the variation of S. As an illustrative example, first consider a spatially independent (horizontal) pro­ file, 5 = 6 {t) only. Then, equation (44) integrates as + Physically, if m; I d r. ~ — +w > 0, particles enter the upper layer and result in the growth 39 of the boundary layer, and 0{t) > 5(0). Conversely if ^ +w < 0, the layer height shrinks, which happens, for example, when re-entrainment occurs. These ob­ servations are sensible physically, and help to verify the applicability of eqn (44). -t-w = 0 and u = constant. Here, eqn (44) Another simple example is i f becomes which has travelling wave solutions S{x,t) = So{x — ut) representing a profile 5q moving downstream with fiuid velocity u. This case represents the case in which the particles settling into the PBL is offset with a vertical velocity w = — 3.2 A pproxim ate Thickness of a PBL The aim of this section is to estimate how thick the PBL is relative to th at of the viscous boundary layer. Equations (33) and (36) were derived for a viscous particledriven gravity current. If the viscous layer is relatively thin, the current could be modelled as an inviscid flow, neglecting the viscous terms from the momentum equa­ tions so th a t the first-order Navier-Stokes equations are generated, and the shallow water model is valid under certain assumptions used to describe the dynamics in the upper layer. On the other hand, if the viscous layer is thick or similar in size to the PBL, the entire fiuid should be modelled as a fully viscous flow and shallow-water theory would not be suitable in this case. To approximate a typical value of 5, the dimensional equation (26) is considered. Let e, T, L, U, a , W , H, and Ôpbl represent respectively typical magnitudes of: 40 volume fraction, time, length, horizontal velocity, Stokes’ settling velocity, vertical velocity, height, and thickness of the PBL of a flow. Then equation (26) can be used to provide the dimensional scaling estimate Substituting W = ^ ^PBL U SPBL _ T L p and T = ^ into the above equation gives 1 / L pp€a To compare with the viscous boundary layer consider an accepted viscous boundary layer thickness [14], let ôyis = where Uoo is the free-stream velocity far away from the boundary. The ratio of Spbl ^ S^s is then given as ( J J l ( — Assuming th at Uoo > U, this simplifies to SpBL ^ 1 ( I U f Ppca \ I U Equivalently, assuming both sides are positive. yÿ ( T ) + Due to the small kinematic viscosity (e.g. water at 20°C has i/ = 10 ® m ^/s 41 àvis SpBL Particle Sizes estimate 53 n m 23 yum 37 ixm e = 1 % 0.03439 0.03177 0.02590 E = 1.5 % 0.03436 0.03170 0.02583 E = 2.0 % 0.03433 0.03165 0.02576 Table 1: The Thickness Ratio of a Viscous Boundary Layer to a PBL of Silicon Carbide Particles < 0.04. Table 1 displays the relative ratios, [14]), it is seen th at calcu­ lated by eqn (45) for silicon carbide, which is commonly used in experiments with density 3217 kgfm^, for a variety of particle sizes with the stokes settling veloc­ ity given in parentheses: 23 g,m (0.00063917 m /s), 37 ixm (0.00165411 m /s), and 53 yum (0.00339402 m /s) . For the calculation, different initial volume fractions are used as e = 1%, 1.5%, and 2.0% such th at (p « 1 , U = \/ g'H m /s, g' = 0.229 m /s^, H = 0.3 m, and L = 7, 6, 4 m for 23, 37, 53 yum of particles, in view of the ex­ perimental parameters in [4]. The values in Table 1 show th at in all cases the visons boundary layer height is less than 4 % of the PBL height. As a typical example, clay, one of the most common materials found in rivers, has a density of 2600 kg/mP and particle size < 0.002 mm in diameter [28]. Choosing a particle size of 37 yum, a = 0.00119377 m /s, and e = 1.5 % gives the relative ratio of 0.03178. Therefore, is inversely proportional to the size and density of particles. Due to the small ratio size, it is assumed th at the current can be modelled as an inviscid flow, especially for larger particles. Experimentally, particle sizes are chosen in a range of 20 yum 200 yum in diameter, and small ratios are implicitly assumed in this model application. 42 3.3 Inviscid System Equations for th e PBL The conclusion obtained in § 3.2 validates the assumption th at the viscous terms can be neglected in the dimensional momentum equations (6) and (7). The dimensionless equations may then be obtained from eqns (6) and (7) directly, and are stated as 'du du du\ dp dx d(j) dx (46) and (47) In (46) and (47), there is a redefined velocity scale U = where H is the total height including the ambient layer. This scale compresses z to be in a smaller range of 0 < z < < 1. If the height scale of the upper layer is also defined in the same manner. Hup = H = ^ , then equation (38) is equivalent to Pp4^up'^s dt dx (P p up + w. P i ) ^ u p 4 ^up T Pi (48) The resulting governing equations of the model, therefore, are (32), (37) (46),(47), and (48) with boundary conditions. (49) and ^ p (x , S,t) = h( ^ l + j 6 up(f>up I , (50) w(x,0,t) =Wup(x,0,t)- (51) 43 Conditions (49), (50), and (51) are obtained from (40), (41), and (43) , respectively. by employing the equalities Uup = U = \/ g'H and Wup = W. In the case of small e, eqn (46) and eqn (48) can be simplified to du du du dp d(j) and respectively. 3.4 A nalytical Solution o f th e PB L Equation The particle transport equation (48) is independent from the others, to the extent th a t it may be solved based on the fluid variables of the upper layer and the known constants pi, Pp, and once the functions of Uup and (pup are specified. To adopt the required functions from existing research ([10], [11]), a small initial volume fraction must be assumed to permit the upper layer description via a shallow-water model. In this model, Uup = Uup{x,t), then Wup = — and the equality (53) holds at z = S. The particle transport equation can then be expressed in the general form of a first order linear nonhomogeneous partial differential equation, where Uup{x,t), 5 ^ ^ = -Wup, and ^ ( p u p { x , t ) are (assumed) known continuous functions for all (x,t) G 2) of points in the (æ, t ) —plane. 44 To solve the equation (54), the method of characteristics [25] is an appropriate approach in order to find the solution, 5. As such, the variables x, t are transformed to the characteristic variables Ç and rj such th at and the determinant of the Jacobian (J) of the transformation does not vanish in D. The original equation (54) is rewritten using the new variables, tjj, ^ and resulting in a simpler PDE, which can then be integrated. Details are shown for the following two cases: constant Uup, and when is an arbitrary differentiable function of x and t. C ase I: Uup = constant Consider the case th at there is a source term keeping the upper layer velocity, Uup, at a positive constant value. Characteristic curves may be described in the form ^ —X Uiipt, and T] — t such th at d et(J) ^ 0. Applying the chain rule to equation (54) and expressing term by F{x ,t ) and ^^0 „p (a;,t) term by f { x , t ) give ^ ^), ;?))V' = /(a;($, 77), (((, ;?))- This integrates as V-K,, ) = where /i = gJo f ((+"«pT,T) dr jg + (56) integrating factor and G is any function of After + 45 performing the integration in (55), the solution of the original variable is given by 5{x,t) = ip{x - Uupt,t), where G(Ç) is replaced by 5(a:,0). The thickness of the PBL is calculated once the initial value of Ô is given. Note th at if w^p = 0, then [ f { ^ + UupT,T)dr + G{^). Jo (56) C ase II: Uy,p = Uup{ x, t ) In most cases, the upper layer velocity likely varies in x and t . In this case, it is difficult to derive a general formula for the thickness of the PBL since the characteristic curves, ^ = Uup{x,t), cannot be generally determined without specifying the u„p function to permit the integration as in the previous case. As such, numerical methods (such as the finite difference approach) may be effective to solve (54). However, it is useful to write out the solution of â when Uup = Uup{t) which is an assumption used in box models. The solution is obtained in a similar manner to Case I, but with a different variable ^ defined as ^ — X - J U u p {t)d t. Using the box model introduced in [10] as an example, Uup is determined by the variable S{t) defined in terms of a current length, S{t) = [l{t)]^, and expressed as an 46 incomplete gamma function, \ 5 /3 This integral equation is solved using the method of iterated kernels with So{t) = 0. Assuming Wup = 0 and Uup{t) = ^ = g-S'(t) (1 - S { t ) ) , the right hand side of eqn (54) ^^(f)up{x,t) = 0.83(1 - S( t) Ÿ Pi by choosing silicon carbide particles in a size of 37 fim, where (/)up(t) = (1 —S{t))‘^ is approximated [10]. As in [10] an approximation of ' q \ 5/3 is chosen, which allows 5{t) 5{t) =J {l (^r\ dr to be calculated as stated in (56) to get \ 5/3' 2 0.83 - .5 J 1 with initial condition (5(0) = 0. Figure 4, generated from Maple, shows the relations of 5 (solid line) and volume fraction of the upper layer (dots) in nondimensional time (t). Note th at the constant term stretches or compresses the height of the PBL. For example, 53 ixra of sil47 2.5 0.5 2.5 0.5 t Figure 4: Changes of the PBL Thickness (solid line) and Upper Layer Volume Fraction (dotted line) Using Si{t) icon carbide has a height 1.5 times^ higher than the one shown in Figure 4, which means th a t the larger the particle size the more distinguishable the PBL. Also, an infinite increase of the volume fraction and 5 when t > 1.7 illustrates the inaccuracy of this approximation to S{t) at larger time scales. This also may be caused by the occurrence of re-entrainment, but without source terms at the bottom boundary this could not occur in the model, which is contrary to a modelling assumption. To increase the validity of the solution in time, the method of iterated kernels and Taylor series expansion is used to obtain the third iterative solution (more accurate than % (()) of S %(() = ^The stokes’ settling velocity is given as a dimensional value in § 3.2, therefore, to calculate one must nondimensionalize % by taking the ratio the dimensional settling velocity over the scales of^. 48 where r is a nonlinear time variable given by r = ( |i ) Note th at r(|) and r(|, r) are complete and incomplete gamma functions, which are defined as e dt and Q Q _22. T5+" 5 ^ M!(|+n)' respectively. The third iteration of S{t) is slightly different than the one given in [10] where part of the incomplete gamma functions are missing. This mistake can be is found to be [ | (r(|) —r(|,r))] ® easily checked by plotting S{t) and (p{t). from a direct calculation. Figure 5 depicts the variations of S and (j)up using % ((). The graph shows th at the thickness of the PBL grows initially with the onset of particle sedimentation, but eventually Ô reaches a maximum value when all of the particles settle down to the PBL. A comparison of the PBL thickness from Figure 4 and 5 is replotted in Figure 6. The two 6 functions calculated using (line) and S 3 (dots) are in fair agreement with each other at initial nondimensional time t < 1 , but diverge for later time. This shows th at including S 3 is crucial for a long time solution. The above example is valid when Wup = 0 is assumed. In the box model, Upp = Uup(t) and Wup needs to be calculated in a manner consistent with the box model. T hat is, Wup{t) = ^ (box models are independent of x), and since hi = A where l{t) is the current length and A is the constant volume per unit width (without loss of generality, assume A = 1), then dh dt dt^l{t)^ 49 dt' 0.4- 0.2 t Figure 5: Changes of the PBL Thickness (solid) and Upper Layer Volume Fraction (dashed) Using %(() Due to the large negative values of Wup at small time scale, 5, solved from eqn (54), is negative at any time when adopting S{t) = Sz{t) from [10], which doesn’t make sense physically. This example shows the inconsistency for a box model applied in the upper layer. To remove the limitation on w^pi another upper layer velocity Uup with both temporal and spatial variations must be chosen. The asymptotic expansion developed in [11] for Uup and is used, and is stated as _ ^ - 3t ' and / + 50 27 2 0 7 /„(ÿ )y ’ 0. 2 - 0.5 t Figure 6: Comparison of the Two Calculated Values for 5{t) Using Si (solid) and S 3 (dotted) where X y , T = PK /? = 6 X 10 ^ (37 fj,m of silicon carbide), and rr/ X 1 1 1 2 - 4 + 4%/ ' 3/ 5 This expansion is valid for t < < % 40.3. W ith u solved from y - Ü = 0, 51 Uup(,x,t), Wup can be 5 4 3 2 --------------------0 Figure 7; 5 (solid) and (dotted) profiles at a; = 0.5 ÔVb w^p{x,z,t) = w{x,h,t ) + — {h - z) dh dh . ,duup du up dx since + -§^{uh) = 0, [11]. Eqn (54) becomes d5 2x dS 26 Therefore, X 4(%,() = V » (^ ,z ), where /J' /a (0.83(^up((T^/^ r )) dr + G ((, 0) 52 (57) Figure 7 gives the profiles of 6 and the upper layer volume fraction at x = 0.5 with the initial condition S(x, 0) = 0 assumed. Note th at after time t = 20 the volume fraction becomes negative, and clearly the solution is not valid. Estimation of the PBL thickness shows this layer is much thicker compared to the viscous boundary layer, implying shallow-water theory is valid for the upper layer with the initial assumption of a dilute suspension, and the viscous effects are negligible in the PBL. This results in a simplification of the model equations from second-order PDEs to five first-order equations. In general, the PBL and the upper layer are to be modelled simultaneously, but to focus on the PBL problem, the upper layer velocities and volume fraction dX z = 5 are assumed to be known such th at the thickness of the PBL can be solved independently from the other four model equations since it only depends on the fluid properties of the upper layer. Strategies of finding 5 are discussed in cases of constant and non-constant upper layer horizontal velocities using the char­ acteristic method. The analytic solution may not be able to be explicitly stated in general, and thus it is suggested th at a numerical approach may be of interest. Two special examples are provided for which 5 is solved analytically for silicon carbide with particle size of 37 iim. The first example with the upper layer velocity expressed as an incomplete gamma function from a box model introduced in [10] illustrates the be­ haviour of the thickness and the upper layer volume fraction. The result, from Figure 5, is th at without particles settling down to the PBL, 5 = 0 implies the non-existence of the PBL, but with the onset of sedimentation, the layer thickness increases until all particles are removed from the upper layer, giving a maximum for Ô. However, 53 the undetermined w^p in box model reveals the difficulty of getting the exact solution of 5 when the effects of Wup are taken into account. Choosing an (x, f)—dependent function of u^p as defined in [11] for the second example, it is seen th at the height of the PBL increases initially, but near the end it starts to decrease. The magnitude of 5 in Figure 7 is larger than the one observed from Figures 5 due to the influence of Wup at the interface and the slower decay of 4>up- The upper layer volume fraction in the box model (Figures 5) decreases rapidly and reaches 0 near t — 3 which stops the growth of Ô, and Wup = 0 is assumed. All the results displayed in this chapter are for the case of small e. Once 5 is specified, one can proceed to solve the remaining four equations which are to be discussed in next chapter. 54 4 A nalysis O f The N avier-Stokes and M ass Con­ servation Equations The aim of this chapter is to find solutions for velocity , pressure, and volume fraction within the PBL. Although, the equations are more difficult to manipulate than the particle transport equation considered in Chapter 3, the small parameter € included in the system, permits an application of direct perturbation theory [12]. Equations and solutions are expanded in a power series in e to generate sets of systems for each order of e” , n = 0,1,2,3... Each solution can be determined iteratively subject to appropriate boundary conditions by using the method of characteristics in two or three dimensions. The general solution form is F = where F denotes any fluid property. General approaches of getting solutions are discussed in cases of 0 = (f>{z,t) and 0 = (f){x,z,t). In general, one cannot obtain analytical solutions due to the complexity of the system of equations. By employing the upper layer velocity and its volume fraction from [11] (details of how to obtain Wup are shown in § 3.4) as an example, the 0 (1) solutions are obtained in the case of = (j){x, z, t). 55 4.1 Solutions o f th e Velocity, Pressure, and V olum e Fraction Equations The equations (52), (47), (37), and (32) are stated as the system /u\ 0 0 0 0 d dt 0 0 0 1 / ■u 0 1 —ez 0 0 0 0 0 0 0 u 1 0 0 0 \ + 0 0 0 0 u \ + dx ( \ u w 0 0 0 0 0 -1 ez w 0 0 w p 0 1 0 0 0 \^/ \« / (58) Multiplying the entire system by the inverse matrix of the third 4 x 4 matrix gives /. 0 0 0 0 0 0 0 w \ 0 0 0 €£ w ,0 0 0 Ï) / \ u d dt w fu w + 0 w w 1 0 0 0 P 0 ( 0) 1° 0 ezu 0 0 -w / 1 0 w ( \ u d dx w ( \ u \ 1 0 0 0 + p U ) 0 1 0 0 0 0 1 0 0 0 V d dz 0 w 0 — 0 p kV /y (59) when w ^ 0. Note th a t to solve the system as in example in § 4.2, eqn (52) is used as the horizontal momentum equation instead of eqn (46). The singularity of the first m atrix of (59) reveals the difficulty involved when calculating exact solutions. (If all matrices in (59) are diagonalizable, then variables may be decoupled and can be solved individually.) As such, the perturbation method is a good alternative approach when a small parameter e is involved. It is noted th at if 4>is known, p can be solved from eqn (47) and consequently u and w remain to be solved simultaneously. Hence, it is worth examining the system for several cases of (f), and following this logical progression of calculating solutions. C ase I: 4> = C onstant 56 A constant value of volume fraction 0 reduces the four equations down to three as du dx du dw dz du du dp and - î - The above equations are actually oversimplified and do not describe the PBL. In fact the PBL does not exist in this case since = 0. C ase II: 0 = (j){z) In this case, = 0 given from equation (37) leads to the non-existence of the PBL region as in Case I. C ase III: 0 = (j){z, t ) Note th at at time t , the variation of 0 in x is much less than in z within a small fixed region so th a t 0 = 0(z, t ) is a reasonable assumption. In this case, a perturbation series is applied as follows 00 0(e, = 0n(z, f)e", (60a) n=0 GO tt(e, X , z , t ) = Y ^ Un( x, z , f)e", (60b) n—0 oo w { e , X, z , t ) = Wn{ x, z , f)e", (60c) ^ P n { x , z , t)e” . (60d) n=0 oo p(e, X, z , t ) = n=0 A detailed but routine substitution of the expansions (60) into system (58) yields sets 57 of problems isolated, according to the powers of e. The zeroth-order problem, 0(1 ), is stated as S + Ç dpQ dt - . « = 0, dz (61c) = (61d) Combining equations (61b) and (61c) simplifies the problem to Since (62b) implies uq = U o { x , t ) if it is initially so, then essentially shallow water theory is valid to leading order. Thus u o = u o { x , t ) = U u p and po = = Pup ( “up” represents the fluid properties of the upper layer). W ith the (assumed) known function of Uup{x,t), Wq and po are solved from (62a) and (62b) to give dUiip Wo{x,Z,t) = - z - ^ , and p„(x, *) = - £ ( 5 ï ^ + « ,(( . respectively. 58 + p„(0, t), Note th at po{0,t) can be solved from the boundary condition (50) since po doesn’t depend on z. Although po = Pup, but in most cases, Pup is not specified. Once Wo has been found, 0o is obtained from (62c) using the characteristic method dx by integrating dz duup(x,t) For r) — t, the general solution of (f>o is then given by 4t(zwO = (63) where the function G can be determined from the boundary condition at the interface (at z = ^, ^o(<^,() = ^wp(()) The first order terms in the perturbation solution are solved from the 0(e) problem once all of the zeroth order terms are known. The 0(e) system may be calculated as 0Z Note th at uq = Uup(x,t), so ^ #Z = 0 in (64b). The first step is to integrate eqn (64c) to get Pi{x,z,t) = - J S - ^ d S + Pi{x,s,t) 59 At the interface, Pup{x,5,t) = height, Pq = ~ Po+ 0{e). Since Pq does not depend on due to continuity. That is, Pn{x, 5,t) = 0 for n > 1. Equation (64b) is rewritten as dui dt dui dx duo{x,t)dui_ ^ dx dz dpi dx duo{x,t) dx The characteristic equations in three dimensions are dt & = dx S = “ “’ dz duo{x,t) du\ _ dpi duo{x,t) Given a curve F, the solution is a hypersurface in (x, y, z, Ui)-space passing through F which is parameterized as z = x{q, r), t = t{q, r), z = z{q, r), and ui = ui{q, r). The family of characteristic curves (solutions of the characteristic equations) is given as X — x{k, q, r), t = t{k, q, r), z = {k, q, r), and ui = ui{k, q,r) a t k — 0 and x = x{q, r), t = t{q,r), z = z{q,r), and Ui = Ui{q,r). The solution of the initial value problem satisfying the parametrized initial data is solved once q and r are solved in terms of x , z , t [31]. As an example, consider the case of a three-dimensional PDE with the initial value of U i { x , 0 , t ) = G{x,t). The family of characteristic curves are given as t = k + X, X — J Uo{ k -h X)dk + t , 60 z = k with T and A as the parameters. The solution ui{ x ,z ,t) , therefore, is ui = G{ t , A) = G{x — zuo{t),t - z). Once Ui is known, wi is integrated directly from (64a) and then can be solved from (64d) using the characteristic method in two dimensions. This example will be discussed in more detail later in § 4.2. The general procedure used for the 0(e) solution may be applied to solve the O(e^) problem. du2 IT du2 , dui , “0 ^ + ^ ^ = 0, OX dz duo , du2 , dui , sT + + “ ‘ a? + = ^ + . . ^ + % OZ (65a) duo dp2 0. (6 5 c ) . , ^ + . # = 0. GZ OZ (65d) As before, solve p 2 from ej)i in (65c) , U2 from p 2 in eqn (65b), wg from eqn (65a), and finally, (j) 2 from W2 in eqn (65d). C ase IV : (j) = (j){x, z, t) The last case considered is when |^ , and ^ are completely general. The 0(1) 61 problem is duQ at duQ dwo d x + dz duo duo (66a) “ dpo dx ’ (66b) (66c) = 0. dt (66d) The system (66) is very similar to (61), except for an additional term in eqn (66d). Going through the same argument as for case 111, uq, Wq and po are solved precisely as before, and the three dimensional characteristic method permits 0o to be solved from (66d) once Wo and uq are known. W ith the zeroth order solution determined, one can seek the next order solutions to the 0(e) problem: " ^ ' + # = 0. dz dz (67c) As before, direct integration of (67c) gives p\. The three dimensional characteristic method is then used to find Ui from (67b), Wi from (67a), and finally from (67d). One can observe th at each 0 (e ”) problem in case IV is almost the same as the 0 (e ”) problem in case III but with additional ^ terms. Using the same strategies as in 0 (e), the order solutions can all be solved iteratively but with more and more 62 terms involved in the system equations. 4.2 Exam ple In this section , a sample calculation is completed using the ideas from § 4.1 to obtain an approximate solution. Using the same upper layer properties [11] as in § 3.4, perturbation expansions to the system (66) are sought satisfying boundary condition (50) rewritten as 1, p{x, 5,t) = h (68) \P p — P i/ where pi = 1000 kg/rn^ (water), pp = 3217 kg/m^ (silicon carbide). Note th at 5 should be solved primarily as discussed in § 3.4 as stated in eqn (57). As solved in the previous section, uq = tt„p, po = Pup, and Wq is found (66a) as where the boundary condition w { x , 0 ,t) = 0 and wo{x, 6 ,t) = Wup{x,5,t) have been used. Once Uq and Wq are found, 0o can be solved from (66d) as ^2/3 \ 3/2' / ' I y See Appendix A for the derivation of 0o. Figure 8 shows the volume fraction profiles for small time scales t = 0.5, 1, 5 63 0.98 O y 1 e 0# O I <ü 0.94 > 0.92 0.9-. 0.5 Z Figure 8: Volume Fraction Profiles for small time scales at a: = 0.5 (from left to right by viewing at the end points) at re = 0.5. One can observe th at the amount of particles decrease with depth, and at smaller times the horizontal profile has very little changes in z. This model gives results consistent with the assumed vertically well-mixed particles in the upper layer because it takes time to stratify the distribution of particles. Note th at the end point value of z on each curve represents the thickness of the PBL at the corresponding time. Figure 9 gives the profiles at larger time scales t = 5, 10, 15, 20. At a fixed height, the volume fraction increases in time, but decreases at later stage. This may happen when the rate of particles flowing in is less than the sedimentary rate at a specific height or when the thickness of the PBL starts to shrink. Figure 10 gives the particle distributions of the two different sizes at t = 10. The larger size of particles has more observable variation with depth, this is in accord with the result reported in [18]. 64 t= 10 0.8 'o 1 t = 15 t = 20 0.2 - z Figure 9: Volume Fraction Profiles for large time scales at x = 0.5 Combining (66b) and (66c) gives duQ dt duo + Uq- ft- = dx dpo dx and hence, po is given (since Pup is unknown) as X Po{x,t) — +po{Q,t), where po{0,t) is calculated from (68). The pressure distributions of silicon carbide are given in Figure 11 for x = 0.1,0.5,0.7. Three curves overlap and are difficult to distinguish, but it ’s obvious to observe th at the pressure decreases in time. Physi­ cally, this can be understood by noting th at the gradual sedimentation of particles reduces the bulk density, and hence the hydrostatic pressure, within the current. The above results don’t take into account the influences of e z ^ and e z ^ from 65 0.9 2 37 wm 0.8 g > 0.5 0.4 0.34 z Figure 10: Volume Fraction Profiles of Silicon Carbide Particles of Sizes 37 n m and 53 fj,m at t = 10 eqn (52) and eqn (47), respectively, which are im portant in the PBL problem. To include the effects of these terms, we solve the 0(e) terms by considering the case when (f) = and (pup = 0«p(O, t) is chosen. Then, 5 = 6 {t) if its initial value is assumed to be zero. By using eqns (62c) po is given by with the properties th a t po = Pup aod = 0 for n > 1. The eqn (64c) gives where Pi{S,t) satisfies the boundary condition (68). Therefore, eqn (64b) can be written as dui dt + 2x dui 3t dx 2 2z dui 3t 3t d z + —Ml — 66 = 0, 0,5 P,0 0.4 0.3 0.2 30 Figure 11: Pressure Profile which gives Wi = 0 satisfying the boundary condition U i { x , 6 , t ) — 0. Consequently, = W i { x , t ) = 0 because w { x , 0 , t ) = W o { x , 0 , t ) + W i { x , 0 , t ) e and w o { x , 0 , t ) = 0. Moreover, solved from d (j)i 2 z d (j)i dt 3t dz = 0 has solution 4>i = 0. The disadvantages of this example are: the e z ^ term is still missing in the horizontal momentum equation, and all 0(e) solutions are equal to 0 except Pi due to the boundary conditions z = 5. Therefore, only the solution for pressure is improved in this example (Figure 12), while the other terms remain unchanged. The three shorter, horizontal behaviour like, curves in Figure 12 are the pressure at t — 0.5,1,5 from top to bottom , and the other three represent the pressure at t = 20,15,10 from top to bottom when viewing at z = 0. The end point value of z on each curve corresponds to the thickness of the PBL at each specified time scale. Pressure decreases initially in time, but starts to increase after a certain stage due to the Pi term. The size of e also affects the magnitude of the pressure, th at is, the 67 0,51 0.4- 0.2 t=20 t=5 t=10 z Figure 12: Pressure Profiles of po + epi with e = 0.02 larger the e, the bigger the pressure. Solving the remaining four variables requires an application of the perturbation method by expanding each solution and equation in a small parameter e. The equa­ tions are then solved iteratively using the method of characteristics to create a series solution for each variable coupled with pre-calculated S and appropriate boundary conditions at z = Ô using an asymptotic series expansion to express the upper layer volume fraction and velocity [11]. This expansion is based on the shallow-water model, and a small initial volume fraction must be assumed, or the shallow-water theory cannot be used to model the upper layer flow due to the presence of particles. A remarkable zeroth-order volume fraction profile varying in depth is observed in Figure 9 for the case of (f) = 4>{x,z,t), and the pressure distribution is displayed in Figure 12 neglecting the variation of ^ 'va x. Figure 10 shows th at when dealing with larger sizes of particles, the changes of particle distribution in z has greater impact on 68 the fluid dynamics within the PBL than small particles. More accurate results may be obtained by incorporating more terms in the expansion, or by taking a numerical approach to solving the original equations. 69 5 D iscussion Both theoretical investigation in turbidity currents [18] and separate experimental results ([13],[15]) have noted a vertical variation in particle concentration, especially near the base of the currents. In this paper, such a region was introduced and defined as a particle boundary layer (PBL). This study was undertaken to construct a model used to predict the dynamics within the PBL and estimate the thickness of the layer for monodisperse, non-cohesive particles suspended in gravity currents produced by the sudden release of fixed-volume suspensions. The current was divided into two layers: the upper layer and the PBL. Governing equations for both layers consisted of mass conservation and the Navier-Stokes equations using appropriate boundary conditions and neglecting re-entrainment, surface tension, and Coriolis forces. An additional equation to capture the changes of particle volume fraction in depth was highlighted for the PBL problem, and a particle transport equation was formulated based on the law of conservation of mass over a fixed spatial region. This equation is not standard in other reported research modelling particle-driven gravity currents. The relative thickness ratio of the viscous boundary layer to the PBL (approximated from the dimensional particle transport equation) provided evidence to remove the viscous terms from the Navier-Stokes equations, resulting in a simplified system of first-order PDEs. The system of nondimensional equations with small aspect ratio and dilute suspension are summarized as Upper Layer: d'^up ~S7 I dwup ~dT ~ 70 ^ r, 0 I dpup d ^ ’u p - w^,p[x, h,t) = — + Uup[x, h, t ) — , d 9 duup I dupp ^ - ^ 4>uph + f f \ ( 0 wp I Ug Uupdz j + (pup^j d u p p __ dppp ^ 9(ppp PBL: du dw â ; + & ■ “’ • - I - I du du du dp d(j) Particle Transport Equation at z = 5: + Upp— — h Wpp . Interface: 9 —Pup{ x, h, t ) = h, — p ( x , 6, t ) = h, Bottom: w{x, 0, t ) = 0. Note th at the reduced gravity g' was defined as g' = for the particle transport equations. 71 and u = Upp, w = Wpp For specified upper layer velocity and volume fraction, the flow behaviour within the PBL was solved (since the required upper layer functions are mostly stated in cases of shallow-water and box models) as an example using perturbation and characteristic methods with the small parameter e. Unfortunately, the box model approach in the upper layer is unsuitable since Wup is neither specified nor calculable and thus 5 cannot be obtained analytically. But, when the shallow-water model was applied by assuming a dilute suspension of particles, the solutions were obtained and depicted in Figures 7 — 12. Note th at the magnitude of 6 is physically unreasonable, it is expected to have a value less than 1 since Uup and Wup were calculated based on the height from 0 to h in the upper layer [11], but it is in a range of 5 to h in this model. Moreover, these results didn’t take into account the influence of the e z ^ term shown in the horizontal momentum equation, and only the zeroth-order problem was solved, except pressure, due to the complexity of equations in higher orders of e. If more terms were adopted from the asymptotic expansions [11] when solving the system, results of greater accuracy could be anticipated. 5.1 Sum m ary o f th e M odel A ssum ptions The model introduced in this thesis was constructed based on the following assump­ tions. 1) A two-dimensional, low-aspect ratio gravity current is generated by the sudden release of a fixed volume of fluid with small particles in dilute suspension into a deep ambient fluid (one-layer flow is considered in this case) flowing over horizontal bottom topography and passes through a horizontal surface. 72 2) Particles are monodisperse, non-soluble, non-cohesive, and the settling velocities are small and assumed to be constant. 3) Particles drift down from the upper layer into the PBL, and accumulate on the bottom of the PBL. 4) Re-entrainment, surface tension, and Coriolis forces are negligible. 5.2 L im itations o f th e M odel When particles are highly concentrated, interactions among particles are significant and must be taken into account when modelling the dynamics of the fluid. Such in­ teractions are excluded in this thesis as only dilute suspensions have been considered. If particles are polydisperse, the settling velocities and the initial volume fractions vary with particle size, and hence more complicated expressions for these two fluid parameters (ug and 0) should be defined. In the case of a shallow ambient layer ^ , the model should be rebuilt based on a two-layer flow capturing the influences due to the motion of the ambient fluid. In addition, the initial turbulent energy may inhibit particle sinking such th at the variation of the particle settling velocity at the beginning of the flow is another factor to be considered [18]. Moreover, modifications of the model are required when modelling currents flowing over a sloping bottom or when re-entrainment is of interest. Direct measurement of the thickness of the PBL is difficult to perform in labora­ tories. Comparing with the observations reported in [8], however, the concentration profile of particles behaves relatively similar as the one obtained from our model ®The initial height of the gravity current is more than 70% of the depth of the ambient fluid [4] 73 (both profiles increase in depth) although the model in [8] is for a steady flow passing through a sloping bottom with re-entrainment and polydisperse particles. Addition­ ally, the influence of the particle size on the volume fraction over depth is in accord with the result shown in [18]. 5.3 Identified Problem s for Future Research The solutions obtained in the examples contain many problems and limitations as de­ scribed above. As such, computing numerical solutions are of interest as an alternate solution method. By setting up appropriate computer code using a numerical method such as finite differences, the systems of equations for both layers can be solved at the same time such th at z and h lie in specific ranges (z G [0,5] and 5 < h « 1) giving sufficiently physical meanings for all fluid variables. Note th at only the hor­ izontal velocities and volume fractions can be solved using finite difference schemes since the temporal variations of pressure and vertical velocities are not specified in the governing equations, but, direct integrations of vertical momentum and incom­ pressible flow equations give the iterative solutions. By removing limitations made for the upper layer equations, systems with both small and large values for e can all be simulated. The system of equations for large e is similar to the one for small e, but with different horizontal momentum equations, particle transport equation, and the boundary conditions for pressure. Upper Layer: 9 1+ I \ f 9uy,p M ^ dUup duup\ dpup j = - i f + 74 dup I ) — p Ç x , 8 f t ^ — h ^1 ~ e .iip < f u p j ) Examining the influences of e will answer the question of how the initial concentration of particles will affect the distribution of volume fraction in the PBL. Direct observa­ tion from the particle transport equation shows th at Ô is also affected by the settling velocity of particles, th at is, if different sizes of particles are chosen, the volume frac­ tion profiles which depend on Ô can then be compared to examine if they are in close accord with the results reported in [18]. 75 A ppendix A 4>q{x , z, t) is solved from eqn (66d) with boundary condition 5, t) = 4>up{x, t) using the method of characteristics in three dimensions. Let A: be a running parameter along a curve F which is parametrized hy x — X, z = 5, and t = t . The characteristic equations are found and their solutions are stated as follows dt , , , — = 1 =» t = k + T dk 3(fc + T) ^ z Rearrange the last two equations to get the solutions of A and r such th at (^o(T, Z, t) = <^up(A, T) - ■#'.p I T ’ f “ T 76 j N otation The following symbols are used in this paper: Symbol Ft Re u h u w Pp Pa P X z t p S Vs 9 L H U W T e P 9 ' J Definition Froude number Reynolds number kinematic viscosity current height horizontal velocity in the PBL vertical velocity in the PBL volume fraction of particles particle density density of ambient layer bulk density longitudinal axis vertical axis, normal to x time pressure thickness of the PBL particle settling velocity gravitational acceleration length scale height scale horizontal velocity scale vertical velocity scale time scale initial volume fraction of particles pressure scale reduced gravity Jacobian Page described on 3 3 3 5, 18 5, 18 18 5 ,1 8 18 18 18 18 18 18 19 19 20 21 26 26 26 26 26 26 28 28 40 Note: the upper layer fluid properties are denoted with a subscript “up” . 77 R eferences [1] M.S. Altinakar, W.H. Graf, and E.J. Hopfinger, “Flow structure iu turbidity currents” , J. Hydraul. Res 34 (1996) 713-718. [2] O.K. Batchelor, H.K. Moffatt, and M.G. Worster, Perspectives in Fluid Dynamics, (Cambridge U.R, Cambridge, England, 2000). [3] T.B. Benjamin, “Gravity currents and related phenomena” , J. Fluid Mech. 31 (1968) 2CKb248. [4] R.T. Bonnecaze, H.E. Huppert, and J.R. Lister, “Particle-driven gravity cur­ rents” , J. Fluid Mech. 250 (1993) 339-369. [5] W.B. Dade and H.E. Huppert, “A box model for non-entraining suspensiondriven gravity surges on horizontal surface” , Sedimentology 42 (1995) 453-471. [6] Density of water, PageWise,Inc. Online, Internet, 20 Jan. 2004, Available http: / / allsands.com/Science/densityofwater-anq-gn.htm [7] M.H. Garcia, “Hydraulic jumps in sediment-driven bottom currents” , J. Hydraul. Eng 119 (1993) 1094-1117. [8] M.H. Garcia, “Depositional turbidity currents laden with poorly sorted sedi­ ment”, J. Hydraul. Eng 120 (1994) 1240-1263. [9] M.A. Hallworth, A.J. Hogg, and H.E. Huppert, “Effects of external flow on compositional and particle gravity currents” , J. Fluid Mech. 359 (1998) 109142. [10] T.C. Harris, A.J. Hogg, and H.E. Huppert, “A mathematical framework for the analysis of particle-driven gravity currents”, Proc. R. Soc. Land. A457 (2001) 1241-1272. [11] A.J. Hogg, M. Ungarish, and H.E. Huppert, “Particle-driven gravity currents: asymptotic and box model solutions” , Eur. J. Mech. B - Fluids 19 (2000) 139165. [12] D.W. Jordan and P. Simth, Nonlinear Ordinary D i f ferential Equations : A n Introduction to Dynamical Sy stems (Oxford University Press, 1999). [13] B. Kneller and C. Buckee, “The structure and fluid mechanics of turbidity cur­ rents: a review of some recent studies and their geological implications” , Sedi­ mentology 47 (Suppl.l) (2000) 62-94. [14] P.K. Kundu, Fluid Mechanics (Academic P. Inc., San Diego, Galifornia, 1990). [15] G.V. Middleton, “Sediment deposition from turbidity currents” , Annu. Rev. Earth Planet. Sci. 21 (1993) 89-114. 78 [16] P.J. Montgomery and T.B. Moodie, “Analytical and numerical results for flow and shock formation in two-layer gravity currents” , J. Austral. Math. Soc. Ser. B40 (1998) 35-58. [17] T.B. Moodie, “Gravity Currents”, J. Computational and Appl. Math. 144 (2002) 49-83. [18] T.B. Moodie, J.P. Pascal, and J.C. Bowman, “Modeling sediment deposition pat­ terns arising from suddenly released flxed-volume turbulent suspensions” , Stud. Appl. M ath 105 (2000) 333-359. [19] T.B. Moodie, “Hydraulic theory and particle-driven currents” , Stud. Appl. M ath 105 (2000) 115-120. [20] T.B. Moodie and J.P. Pascal, “Nonhydraulic effects in particle-driven gravity currents in deep surroundings”. Stud. Appl. Math. 107 (2001) 217-251. [21] T.B. Moodie, J.P. Pascal, and C.E. Swaters, “Sediment transport and deposition from a two-layer fluid model of gravity currents on sloping bottom s” , Stud. Appl. Math. 100 (1998) 215-244. [22] Y. Nino, F. Lopez, and M. Carcia, “Threshold for particle entrainment into suspension”, Sedimentology, 50 (2003). [23] A. Nir and A. Acrivos, “Sedimentation and sediment flow on inclined surfaces” , J. Fluid Mech. 212 (1990) 139-153. [24] R.J. LeVeque, Numerical Methods f o r Conservation Laws (BirkhauserVerlag, Basel, 1992). [25] P. V. O ’Neil, Beginning Partial D i f ferential Equations, (Wiley and Sons, 1999). [26] H.M. Pantin, “Interaction between velocity and effective density in turbidity flow; phase-plane analysis, with criteria for autosuspension” , Mar. Ceol 31 (1979) 5999. [27] C. Parker, “Conditions for the ignition of catastrophically erosive turbidity cur­ rents” , Mar. Ceol 46 (1982) 307-327. [28] Soil sizes, NASA’s Coddard Space Flight Center, Online, Internet, 22 August 2001, Available http://ltpwww.gsfc.nasa.gov/globe/activ-98/soilsizs.htm [29] R.S.J. Sparks, R.T. Bonnecaze, H.E. Huppert, J.R. Lister, M.A. Hallworth, H. Mader, J. Phillips, “Sediment-laden gravity currents with reversing buoyancy” , Earth Planet. Sci. Lett. 114 (1993) 243-257. 79 [30] James Stewart, Single Variable Calculus, 3rd ed. (Brooks/Cole Publishing Com­ pany, 1995) [31] E. Zauderer, Partical Differential Equations of Applied M athematics (Wiley and Sons, Inc, 1998). 80