L I N E A R S T A B I L I T Y OF A B E R M A N F L O W I N A C H A N N E L PARTIALLY FILLED WITH A POROUS M E D I U M by CHUNTAODENG B.Eng., South China University of Technology, 1998 M.Sc, Asian Institute of Technology, 2000 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF CHEMICAL AND BIOLOGICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA November 2004 © Chuntao Deng 2004 ABSTRACT In this thesis, the coupled flow of a Newtonian fluid both above and through a porous medium is considered. This work was motivated from a pulp and paper application, namely twin-wire forming. In the fluid-only region, the two-dimensional flow field is governed by the Navier-Stokes equation. In the porous region, the flow field is governed by the Brinkman-extended Darcy law relationship. Inertial terms are retained in both regions and the interface conditions between the two domains are those as outlined by OchoaTapia and Whitaker (Int. J. Heat Mass Transfer 38, 2635 1995). The model equations were solved using two independent methods. In the first method we develop a similarity transform and reduce the governing equations to two, coupled, non-linear ordinary differential equations to form a three-point boundary value problem. This was solved numerically and validated analytically by examining asymptotic cases. Three characteristic solutions were identified and the stability of each was examined by the method of normal modes. In the second numerical approach, the governing equations were re-posed as a one-domain problem, using the procedure outlined by Basu and Khalili (Phys. Fluids 11, 1031 1999), so that the conditions at the interface need not be considered. The resulting equation was solved directly, in primitive variable form, using afinitevolume formulation. Finally, an experimental device was constructed to compare to the numerical predictions. Eight test cases were performed, using two different porous media, in which the velocity profile of the fluid was measured using pulsed ultrasound doppler anemometry (PUDA). Good agreement was found between the similarity numerical predications and the experimental measurements. ii Contents ABSTRACT ii Table of Contents iii List of Figures vi List of Tables xi Nomenclature xii Acknowledgment 1 xvii INTRODUCTION 1 2 LITERATURE REVIEW 6 2.1 Rheology of Fibre Suspensions 6 2.2 Papermaking 13 2.3 Objectives of the Research 15 3 M O D E L EQUATIONS 4 17 3.1 Introduction 17 3.2 Model Equations 21 SIMILARITY SOLUTIONS 23 4.1 Similarity Formulation 24 4.2 Steady-State Solutions 25 iii 4.2.1 Asymptotic Solutions 26 4.2.2 Series Solutions 30 4.2.3 Numerical Solutions 31 4.3 Stability of the Similarity Solutions 39 4.4 Summary 40 5 ONE DOMAIN PROBLEM 47 5.1 Introduction 47 5.2 One-Domain Formulation 47 5.3 Summary 49 6 EXPERIMENTAL WORK 54 6.1 Introduction 54 6.2 Experimental Apparatus 54 6.3 Pulsed Ultrasound Doppler Anemometry 58 6.3.1 Measuring Principle 58 6.3.2 Adjustment of Parameters . 59 6.4 Experimental Procedure 62 6.5 Data Evaluation 63 6.6 Results and Discussion 64 7 S U M M A R Y AND CONCLUSIONS 76 8 RECOMMENDATIONS 78 A DERIVATION OF T H E PRESSURE B O U N D A R Y CONDITION 91 A.l Normal Stress 91 A.l.l Fluid Region 92 A. 1.2 Porous Region 93 A.2 Pressure Drop 93 B T H E A S Y M P T O T I C SOLUTIONS iv 95 C C O D E VALIDATION 97 D C O D E FOR T H E SIMILARITY SOLUTION 99 D . l Shooting Method 99 D. 2 The Reverse Solver E T H E A L G O R I T H M AND ' C O D E FOR ONE DOMAIN P R O B L E M E. l Introduction 125 125 E.l.l Formulation of the Problem 125 E.l.2 Setup of the Boundary Conditions 130 E.1.3 Linearization and Solving the Linear Equations F 108 132 DETAILED DRAWINGS OF T H E A P P A R A T U S 154 F. l 154 Detailed Drawing of the Apparatus F.2 Description of the Pump 165 F.3 Seeding particles 165 v List of Figures 1.1 The upper figure is a sketch of the first industrial paper machine, the Fourdrinier former, and the lower figure is a schematic drawing of a modern twin-wire hybrid former [1] 3 1.2 A twin-wire former of a pilot paper machine (Paprican) 4 1.3 The channel flow considered 4 2.1 Yield stress versus crowding number from the data given by Bennington et al. [2]. The line shown is a regression using Eq. (2.6). Two different pulp types are shown: the SBK pulp is shown as a circle (o) and the T M P as a square ( • ) . Bennington et al. report experimental uncertainties as large as 100% 10 2.2 A schematic of the geometry in twin-wire forming [1] 13 2.3 A schematic of the growth of the web during one-sided dewatering [1]. . . . 13 3.1 The channel flow considered 18 4.1 (i): The x-component of the velocity profiles f . y and g (normalized by the y averaged velocity V—xa,t position x) estimated in the limit of Re —> 0. Four curves are shown: (a) Da/e —> 0, (b) Da/e = 5 x 1 0 , (c) Da/e = 1 x 1 0 , - 4 - 2 and (d) Da/e —>• oo. (ii): An estimate of the boundary layer thickness in the porous medium. The boundary layer thickness has been normalized to the thickness of the porous medium. These estimates were made with'/? = 0.2, X = 0.8 and e = 0.9 28 vi 4.2 The effect of Re for the case Da = 1 x 10~ on the the x component velocity 3 profiles, f y 4.3 The 6 th and g , here 0 = 0.2, x — 0-8, e = 0.9 29 y order series solution (-) compared to the similarity numerical solu- tion (•) with Da = 1 x 1(T , Re = 5, 0 = 0.2, x = 0.9 32 3 4.4 Example solutions: The plot represents the three different families of solutions obtained, here (a): Re = 10.93, (b): Re = 16.68 and (c): Re = 32.21. These simulations were conducted at Da = 0.05, 0 = 0.2, x — 0.8 and e = 0.9. The shaded area represents the porous region 4.5 36 Example solutions: The figure represents the shear stress at the interface —fyy(x) (here the circles indicate the series solutions given in §4.2.2). These simulations were conducted at Da = 0.05, 0 = 0.2, x = 0.8 and e = 0.9. The type of solution is indicated, where the line "|" shows the transitions from type-I to type-II and from type-II to type-Ill solutions. Here the solid-lines indicate stable solutions while the dotted region represents unstable solutions. 37 4.6 A 3-dimensional plot of the solutions at Da and Re considered, here 0 = 0.2, x 4.7 = 0.8, e = 0.9 38 The stability analysis for Da — 0.013: (i) is the steady solution, the detailed view in Region I is shown in (ii), whilst (iii) shows the real part of the corresponding eigenvalues in Region / . Here the dashed lines indicate unstable solutions 4.8 41 (a) Summary of solution behavior for simulations with 0 = 0.2, x — 0-8 and e = 0.9. Nine different solution regions are found and the detailed information about each region is tabulated in Table 4.1. (b) Region C is enlarged 4.9 42 Example solutions: the shear stress function —f (x) yy versus Re for (a) Da = 0.013, (b) Da = 0.026 and (c) Da = 0.0375. The type of the solutions are denoted, and the multiple solution regions are labelled A, B, etc. representing different type, stability and number of solutions. Here solid-line denotes stable solution and dotted-line, the unstable solution. . . vii 44 5.1 The flow fields as estimated by the two different numerical approaches, here AP X is evaluated at the interface y = 0.8 and AP is evaluated at x = 7.0. y The porous medium extends from y = 0.8 to y = 1, with Da = 1 x 10~ , 3 Re = 5 and e = 0.9. In the lower panel of the figure, the fully-developed profile using the one-domain approach at x = 7.0 (shown as a dashed line) is compared to the similarity solution (shown as a solid line) with J3 = 0.2. 5.2 50 The flow field determined as a function of (5 using the similarity approach. The normalized x-component of velocity (f for fluid region and g for porous y y region) is given in the figure on the left and the y-component (/ for fluid region and g for porous region) is given on the right. The porous medium (shown as the shaded area) extends from y = 0.8toy = l . This simulation was performed with Da = 1 x 10~ , Re — 5, and e = 0.9 3 51 5.3 The factors affecting the empirical constant f3 in the shear stress jump condition (the tested cases are for x — 0-9, e = 0.9) 52 6.1 the schematic of the experimental setup 56 6.2 the assembled test section in (a) exploded view and (b) assembled view . . 57 6.3 The working mechanism of pulsed ultrasound Doppler anemometry (image taken from Signal-Processing SA) 59 6.4 The interface of the pulsed ultrasound Doppler anemometer Dop2000 from Signal Processing SA 61 6.5 The velocity measured off the probe is the component velocity of particles moving in the ultrasound beam direction. Here v denotes the flow velocity and v denotes the velocity of flow in the ultrasound beam direction. It p shall be noted that 9 denotes the Doppler angle, where the refraction of the sound beam due to the sound speed difference in plexiglass and fluid has been accounted for 63 6.6 The schematic of the projections of the x component of the numerical velocity v , y component v on the ultrasound beam direction, the PUDA x y velocity v 65 p viii 6.7 The similarity numerical solutions shown as solid line (here (5 = 0.13) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 1 in Table 6.2 6.8 The similarity numerical solutions shown as solid line (here (3 — 0.3) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 2 in Table 6.2 6.9 The similarity numerical solutions shown as solid line (here (3 = 0.9) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 3 in Table 6.2 6.10 The similarity numerical solutions shown as solid line (here P = 1.5) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 4 in Table 6.2 6.11 The similarity numerical solutions shown as solid line (here P = 0.15) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 5 in Table 6.2 6.12 The similarity numerical solutions shown as solid line (here p = 0.58) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 6 in Table 6.2 ix 6.13 The similarity numerical solutions shown as solid line (here (3 = 1.1) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 7 in Table 6.2 74 6.14 The similarity numerical solutions shown as solid line (here /J = 1.1) for the . velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for the two steady state solutions for Experiment No. 8 in Table 6.2 75 B . l The shear stress function —f (x) + the interface for different Re, it coma yy pares the similarity numerical solution (o) with the asymptotic solution (-). The simulation is for the case Da = 1 x 10~ , x = 0.8, (3 = 0.2 and e = 0.9. 3 96 C . l The shear stress function —g (l) at the porous wall for different-Re, it yy compares the similarity numerical solution (•) for Da = 1000, x = 0.8, (3 — 0 and e = 1 with Berman solution (-) [3, 4] 98 F.l the 3-D view of the cover of the upstream-tank . 155 F.2 the dimensions of cover of the upstream-tank 156 F.3 the 3-D view of the lid of the test channel 157 F.4 the dimensions of the lid of the test channel 158 F.5 the 3D view of the screen plate (filter) for the suction 159 F.6 the dimensions of the screen plate for uniform suction 160 F.7 the 3-D view of the suction element (diffusor) of the device 161 F.8 the dimensions of the suction element (diffusor) 162 F.9 the 3D view of the upstream and downstream tanks connected with a channell63 F.10 the dimensions of the upstream and downstream tanks connected with a channel 164 List of Tables 4.1 The number, type and stability of the solutions for each solution region as shown in Fig. 4.8 43 6.1 The porous fabrics used in the experiment 55 6.2 65 The experimental conditions tested xi Nomenclature a yield stress constant b yield stress constant B binary constant to turn on/off the sink term c velocity of sound in a liquid Cf form drag constant C mass concentration of wood pulp suspensions C sediment concentration d(£) the interface between the paper web and the suspension D nominal diameter of the fibre D(<p) a function of the permeability and compressibility of the fiber network Da Darcy number, Da = K/H f the emitting frequency of the ultrasound pulse fd the frequency shift of the received echo from the particles in PUDA measurement f(y) similarity velocity function in fluid region F(y) the perturbed similarity function for the fluid region g(y) similarity velocity function in porous region g(£) fabric position g acceleration due to gravity G(y) the perturbed similarity function for steady case problem in the porous region h height for fluid region H total height for the channel m s e 2 xii J(£) the dewatering rate (flux) k(4>) a function to describe the fiber-fiber contact stress K permeability of the porous medium I fiber length Lj interaction length for crowding number m empirical constant of the power-law function M inter-phase momentum term (vector) n empirical constant of the power-law function n number of contacts between fibers N crowding number c Ni gel crowding number p unit vector that indicates fiber orientation ge p(li) P fiber length distribution pressure P pressure at position just above the interface P_ pressure at position just below the interface Pd the depth of the particle measured by ultrasound Pmax the maximum depth of the particle that measurable by PUDA PRF pulse repetition frequency q rescaled position of porous wall r aspect ratio of the fiber + d r the radius of curvature Re suction/injection Reynolds number, Re = U H/v S specific surface of the fiber defined as the ratio of the surface area to its volume t time td the time delay between an emitted burst and the echo issued from the particle tprj the time between two emissions % 0 w xiii T fabric tension u velocity vector in porous region u velocity component in porous region «i x velocity component in porous region Uy y velocity component in porous region U fabric moving velocity u uniform suction velocity at porous wall V mean axial velocity at x = 0 Vsed settling velocity of fiber V velocity vector in fluid region V velocity component in fluid region w V p velocity component in the ultrasound beam direction v x velocity component in fluid region Vy y velocity component in fluid region x ^max maximum velocity that can be measured using PUDA W velocity vector for one-domain approach Wk water retention ration X x coordinate y y coordinate z rescaled y coordinate Subscripts X derivative in x direction y derivative in y direction Superscripts * the corresponding dimensional properties, SI unit. the time dependent functions Greek letters j3 shear stress jump coefficient j3 the quadratic shear stress jump coefficient 2 P the artificial compressibility constant X the dimensionless interface position between the fluid and porous region 5 phase shift of the received echo for PUDA (5(£) web thickness Ap density difference between fiber and surrounding fluid e porosity of the porous media 4> solidity or volume concentration of the fibers (f) Gel concentration rj the n direction in a nonorthogonal coordinate A rescaled similarity velocity function for fluid region p dynamic viscosity of fluid [if apparent viscosity of the pulp suspension p effective viscosity v kinematic viscosity 0(0) function in Forchheimer model to describe M g e p a r fluid density the eigenvalues of the linear stability equations fiber-fiber r . contact stress tensor normal stress r yield stress of pulp suspension LO fiber coarseness Vi\ domain for flowing suspension region 0 domain for paper web region £ the £ direction in a nonorthogonal coordinate H(^) function in Forchheimer model to describe M xv ib rescaled similarity velocity function for porous region \I/ stream function <; intermediate constant, equal to •d a dimensionless constant for slip-velocity boundary condition at the interface y/e/Da xvi ACKNOWLEDGEMENT I wish to thank those who have offered their help and support to this work. In particular, I sincerely thank my supervisor Mark Martinez for his guidance, encouragement and help. My sincere thanks also goes to Dr. Bruce Bowen, Dr. James Olson, Dr. Richard Branion and Dr. Richard Kerekes, for their time and suggestions in this work. I would like to thank Peter Taylor, Doug Yuen, Tim Paterson and Alex Thng for helping the experiment setup; Ken Wong, Brenda Dutka and Lisa Hudson for their technical and administrative assistance. My friends Edmoncl, Jacky, Leila, Monica, Robin, Satya, Yaoguo, Yubing, and Zhaolin, deserve my great thanks for being so supportive. At last, this work could not be accomplished without the consistent support and love from Lili, my wife; and encouragement from my brother, sister and their families, and my parents. The financial support from Natural Sciences and Engineering Research Council of Canada (NSERC), University of BC Graduate Fellowship (UGF) and Advanced Papermaking Initiative (API) are acknowledged. xvn Chapter 1 INTRODUCTION The earliest recording of paper manufacture has been credited to the Chinese around 100 A.D., where a suspension of silk and mulberry was filtered in a bamboo mold. The first industrial paper machine was developed around 1800. Since this time, paper machines have been continually under development leading to the modern machines called twin-wire formers (see Fig. 1.1); the reader is referred to Norman's review for further details [5]. Here the term "wire" (or more recently termed "fabric") refers to a woven, usually plastic, permeable membrane through which the papermaking suspension is drained. During processing, the papermaking suspension is initially fluidized by turbulence created locally from a sudden change in geometry in an apparatus called a "headbox". In this case, the fibre network is broken down into smaller floes and single fibres with weakly correlated velocities [6]. Suspension fluidization and the associated breakage of floes is attained by inducing turbulence [7, 8, 9] and is aided by the addition of chemical deflocculants [10, 11, 12, 13, 14, 15, 16]. The fluidized fibre suspension is then formed into a plane liquid jet and distributed in the gap formed between the two moving wires (see Fig. 1.2). The velocity of the jet is between 10 — 35 m/s and it is 5 — 10 m wide and 1 cm thick. Often, however, this fluidized state is transient and the network re-flocculates within a few milliseconds upon the reduction of the applied shear [6, 17]. 1 1 W i t h non-traditional forming processes, namely high consistency (HC) forming, previous workers have attempted to freeze the structure of the sheet before it leaves the headbox. In the design put forth by Grundstrom et al. [18], the pulp suspension is first dispersed by turbulence, generated by numerous bends in the flow channel, and then allowed to be reflocculated in a small straight outlet channel before impinging 1 CHAPTER 1. INTRODUCTION 2 During roll-blade forming, water is filtered from the suspension in two distinct steps. In the first step water is forced through the wire by the momentum of the jet and by the tension and convergence of the fabrics passing over the roll. The fabrics are then passed over a series of blades, in the second step, where greater water removal is achieved. In both these steps the change in curvature of the fabrics generates a pressure field to balance the fabrics' tension. Early attempts to characterize the action of twin-wire formers were largely concerned with estimating this pressure. Hergert and Sandford [22] measured the pressure distribution on the surface of a roll during one-sided dewatering and measured gauge pressures as high as 5 kPa. Zhao and Kerekes [23] reported gauge pressures as high as 10 kPa in blade formers. In the late sixties and early seventies both Baines [24] and Meyer [25] made rigorous attempts to model the twin-wire dewatering process in order to better understand the changes in sheet structure. Although most of their assumptions were reasonable, they did not provide any solution to the resulting equations due to the complexity of the problem. In-subsequent work, a number of authors have simplified these analyses by considering the suspension as a Newtonian fluid and neglecting the variation of properties in the drainage direction [23, 26, 27, 28, 29]. These models have engineering value but give no insight into the development of paper properties. A full, theoretical analysis of the forming process would be extremely challenging. Furthermore, validation of such a model would be difficult, if not impossible. An alternative approach is to focus only on one aspect of the process, and neglect the other factors. In our case, the aspect addressed in this thesis is the interface between the suspension and the paper web. The geometry considered in this work is shown in Fig 1.3. The relevant literature will be reviewed in Chapter 2 of this thesis. At the end of this chapter the current state of knowledge is summarized and the objectives of the research on a moving wire. A similar concept was advanced by Nomura and his co-workers [19, 20, 21] but in this case turbulence was generated by a sudden change in geometry achieved in a saw-tooth shaped channel. Both of these concepts are similar in that turbulence intensity is strongly coupled to the mean velocity of the suspension. Hence fioc dispersion is achieved through an increase in jet velocity. Because of this distinctive forming method the structure of the H C paper differs from that of conventional paper. CHAPTER 1. INTRODUCTION headbox 3 forming table suction foils wet flat boxes forming board breast roll table rolls foils 1 suction flat boxes suction couch roll Y m i o T m T C n UUUUUUUU trays stretch guide wash forward drive roll return rolls Figure 1.1: The upper figure is a sketch of the first industrial paper machine, the Fourdrinier former, and the lower figure is a schematic drawing of a modern twin-wire hybrid former [1]. CHAPTER 1. INTRODUCTION Breast Roll Figure 1.2: A twin-wire former of a pilot paper machine (Paprican) ooooooooooooooooooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooooooooooooooooooo O O O O O P QP.O 0_Q O O P P O O P O P P P P P P Q P P P P P P P P O Q P P P P P P P P Q P 0_P_OQO_O Q Figure 1.3: The channel flow considered. CHAPTER 1. INTRODUCTION 5 work are stated. The model equations are presented in Chapter 3. In Chapter 4 the model equations are transformed to a three-point boundary value problem and solved using a similarity ansatz. The stability of the solutions is also considered in this chapter. In Chapter 5, the full set of equations is solved using a standard computational fluid dynamic (CFD) approach and validated experimentally in Chapter 6. In Chapter 7 we summarize the major findings in this work and the recommendations are given in Chapter 8. The appendices at the end provide additional information about the CFD algorithms, the numerical codes and the experimental apparatus. Chapter 2 LITERATURE REVIEW 2.1 Rheology of Fibre Suspensions Papermaking fibres are hollow, flexible, rod-like particles that have a wide distribution in both length and diameter depending upon species and growing conditions; when liberated from trees, they can curl or kink. North American fibres contain typically 40-50% cellulose, 20 to 35% hemicellulose and 15 to 35% lignin, with the remaining fraction consisting of resins, tannins, ash and miscellaneous compounds [30]. Papermaking fibre suspensions are a mixture of wood fibres, water, clay, as well as small quantities of inorganic salts and polymeric additives. This physical make-up leads to complex rheological behavior with the suspensions showing complex stress/strain-rate relationships. The major feature of this relationship is that it conforms to that of an inelastic non-Newtonian fluid. Under typical processing conditions, papermaking fibre suspensions aggregate into coherent networks which possess measurable mechanical strength [7, 8, 10, 11]. The reviews of Kerekes [31], Jokinen and Ebeling [12], and Powell et al. [32] and the references contained therein summarize our current knowledge of fibre suspension rheology. Although fibre networks form by surface charges, Mason [33] found that mechanical entanglement rather than colloidal forces was the principle source of fibre flocculation. Fibres collide because of the relative motion induced by the flow field. They entangle, bend, and remain networked from frictional forces transmitted by fibres that are locked into bent configurations [34]. Soszynski and Kerekes [35] confirmed the influence of bending stresses on 6 CHAPTER 2. LITERATURE REVIEW 7 floe strength through stress relaxation experiments. They defined the propensity for papermaking fibres to flocculate in terms of a dimensionless number, termed the crowding number, N, i.e., "-Hiff"^ (21) where C is the concentration of the solid material in the suspension (kg/m ); Lf is an 3 m interaction length based upon the length of a fibre (I); co is the fibre coarseness defined as the mass per unit length of fibre; D is the diameter of the fibre; and <fi is the solidity or volume concentration of the fibres. 1 It should be noted that, for polydisperse suspensions, Kerekes and Schell [37] recommend the use of the first moment of the fibre length distribution to estimate the interaction length, i.e., Lf = J>(Ji)Ji (2-2) This relationship however has yet to be justified. The physical significance of N for monodisperse suspensions is given by its definition; it reflects the number of fibres in a spherical volume of diameter equal to Lj. Flocculation does not generally occur when N < 1 as on average each fibre can rotate freely. The range of 1 < N < 60 represents a regime of fibre behavior described by Soszynski [38] as one of 'forced collisions' between fibres in simple shear flow. Finally, iV can be related to the average number of contacts n between fibres c in a network using [39] ( 2 - 3 ) Most papermaking operations in British Columbia occur in the range between 17 < N < 60 [40]. Understanding the motion of dilute (N <gC 1) fibre suspensions is difficult. Insight into this phenomenon can be gained by first examining the simplest case of the motion of a single isolated fibre settling under Stokes' flow conditions, i.e. [41], ^ sed 1 = ^^[( l n 2 ^ + 0.193 + O(ln2r )- )g + ( l n 2 r - 1.807 + O(ln2r,)- )(p-g)p1, (2.4) 1 d 1 d Soszynski [36] investigated the relationship between the mass and volume concentration of wood pulp suspensions, i.e. (f> « ( | + Wk)C , m where Wk is the water retention value. CHAPTER 2. LITERATURE REVIEW 8 where Ap is the difference in density between the fibre and the surrounding fluid; r is the d aspect ratio of the fibre defined by l/D; u. is the viscosity of the fluid; g is the acceleration due to gravity; and p is a unit vector that indicates fibre orientation. Unlike spheres, fibres can have significant motion perpendicular to gravity with a drift velocity strongly dependent upon its orientation. Ross and Klingenberg [42] considered this case for papermaking fibres by dropping one fibre at a time into a quiescent fluid to create a sediment. Using techniques similar to those employed by Doi and Chen [43] in polymer dynamics, these authors modelled papermaking fibres as a chain of prolate spheroids connected through elastic ball and socket joints in order to account for flexibility effects. With semi-dilute suspensions (1 < N < 60), this situation becomes more complex as each individual fibre moves under the influence of the long-range hydrodynamic disturbances of the other particles. As an example, during sedimentation, these disturbances lead to inhomogeneous settling rates [44] and local floe formation. With monodisperse suspensions, the latter has been observed by Kumar and Ramarao [45], and confirmed by Herzhaft et al. [46] (called clumps by these authors), who reported the existence of floes for monodisperse glass fibres settling at low crowding numbers. Kumar and Ramarao also noted that, as the fibre concentration increased, the number of floes increased and caused greater hindrance effects. With bidisperse suspensions, Davis [47] observed size segregation up to a critical concentration, above which particle interlocking prevented further segregation. It was found that the smaller the ratio of the large to small particle diameters, the smaller the critical concentration. Quite possibly due to these complexities, there are only a few mathematical models of the motion of interacting particles [41, 48, 49]. Koch and Shaqfeh [49] considered the sedimentation of stiff, axisymmetric, nonspherical particles in a dilute suspension and demonstrated that particle flocculation can occur. Moreover, from a stability analysis, they suggested that the sedimenting suspension should segregate into regions of high particle density and cause a relatively strong downward convection of particles. With N > 60, many authors consider the network to behave like a compressible porous medium [50, 51]. This approximation stems from the fact that at this level offibrecrowding, there are, on average, 3-4 interfibre contacts [34, 52, 53, 39] and according to Steenberg et CHAPTER 2. LITERATURE REVIEW 9 al. [54], this is the minimum number of contacts required to create a network. In this regime, it is widely known that pulp suspensions do not flow until a certain critical (or yield) shear stress is exceeded. Suspension yield stresses have been either measured directly using viscometers, or inferred from dynamic shear and tensile measurements [31]. For pulp suspensions, the yield stress r has been typically correlated with the suspension 0 mass concentration using equations of the form [55] r 0 = a(C (2.5) - C) b m s where a and b depend on the test procedure employed, the pulp type, and the degree of pulp treatment (e.g. beating) prior, to testing; and C is the sediment concentration. s Values of the constant b have been reported in the literature to range between 1.25 to 3.02. Few, if any of these works report the yield stress in terms of the crowding number. As an example, we do so, for the data given in [2], see Fig. 2.1, and correlate this data set with an equation of the form r = 5.5 x l O " ^ - N ) \ 4 2 0 where N geL gel R = 0.96 (2.6) 2 is the gel crowding number defined by Martinez et al. [56] to be N geL = 16. At this level of fibre crowding, the transport of mass and momentum can be expressed by the Eulerian volume-averaged two-phase flow equations [57, 58, 59]: 'V-(<££*) = 0 (2.7) V-[(l-c/>K] = 0 (2.8) - 0 V P * - V -f + pf(j>g- M (2.9) (pPftf-VvT (1 - <p)p w* • Vw* w = = + -(l-</>)VP* + p (l-<f>)g + V • [ p ( l - 0 ) ( V t Z f + W ) ] r w M (2.10) where u* and w* are the dimensional velocities (vector quantities, and * indicates a dimensional property) of the fibres and the water, respectively; P* is the fluid pressure; M is the inter-phase momentum term (vector quantity); pj and p w are the densities of the fibre and water, respectively; and f represents the fibre-fibre contact stress tensor as the CHAPTER 2. LITERATURE REVIEW 10 Crowding Number, N (dimensionless) Figure 2.1: Yield stress versus crowding number from the data given by Bennington et al. [2]. The line shown is a regression using Eq. (2.6). Two different pulp types are shown: the SBK pulp is shown as a circle (o) and the T M P as a square (•). Bennington et al. report experimental uncertainties as large as 100%. CHAPTER 2. LITERATURE REVIEW 11 fibre phase is conceptually envisioned as a continuum [60]. Surface tension effects between the phases have been neglected. In the one-dimensional filtration studies of Ingmanson et al. [50] and Meyer [51], where the relative fluid velocity was about 1 m/s, they described M with a Forchheimer-type 2 equation. These expressions have been extended to two dimensions by Zahrai [1], i.e., 6(0) = -(<[>) , __ { v v , ) 0~ / [l + 570 ] 3 1 2 L 3 J 0<O.3 (2.12) 0>O.3 = 0.1y/Q(#) (2.13) where the quantity S is the specific surface of the fibre defined as its surface area to volume a ratio. The coefficients given in the above expression are taken from [51], but the reader is referred to a number of references found in the literature [50, 51, 61, 62, 63, 64, 65, 66, 67]. Finally, the remaining term which needs to be defined is the particulate stress term, r. There is evidence in the literature of two different modes of behavior forfibroussuspensions. i) With regards to the first mode of behavior, a large body of literature is available which ascribes f as an implicit function of the local solidity of the network, i.e., T — k(<fi)I [50, 63, 68], where / is the identity tensor. A number of authors have represented the functional relationship between k and cb by power-law functions of the form [69, 70, 71]: fc(0) = m(0" - 0^) (2.14) k{(b) = m(0 - 0 ) (2.15) or n 9 where 0 is the "gel" point and is defined as the value below which k(<p) cannot g be distinguished from zero. Here, m and n are empirical constants determined from static compressibility tests and are defined uniquely for each relationship given above. It should be noted that this formulation for r neglects the shear stress that comes about from network deformation. It thereby assumes isotropic material behavior. 2 T h i s represents a Reynolds number of about 35, based upon the diameter of a fibre. CHAPTER 2. LITERATURE 12 REVIEW These formulations neglect plastic deformation and time-dependent behavior that have been shown to be important. ii) On the other hand, pulp suspensions fluidize under applied shear. Bennington and Kerekes [8] reported the apparent viscosity of a semi-bleached kraft suspension (SBK) can be approximately represented by a function of solidity according to 3.1 (2.16) As a result, the particle stress term for the fluid-like model may be estimated as [1] r = p 0(W* + W * ) T 7 (2.17) To help explain the differences between these expressions for r, it is helpful to reiterate the assumptions upon which the equations are based. Clearly, both imply that f varies with solidity. However, the fluid-like description assumes that r varies with the velocity gradient of the fibres. The fluid-like description assumes a complete dispersion of the fibre suspension; it is free of entanglements. This implies that fibre interaction occurs through the fluid phase. The power law model, on the other hand, tacitly assumes mechanical interaction between the fibres. In other words, the fibres in the suspension do not necessarily have to move at the same velocity as the fluid in the suspension; their motion can be hindered by entanglements. In summary, i) The mathematical description of the motion of dilute (N <C 1) and semi-dilute suspensions (1 < N < 60) is limited. In essence, the description of the long-range hydrodynamic force, fibre collisions and entanglement have yet to be tackled from a theoretical standpoint. ii) With concentrated suspensions (N > 60) the equations of motion are not a closed set of equations. The success of a mathematical model in this case is dependent on the accuracy of the closure equations for inter-phase momentum transport and the fibre-fibre contact stress. In the present state of the literature, the currently 2. LITERATURE CHAPTER REVIEW 13 R Fi gure 2.2: A schematic of the geometry in twin-wire forming [1]. 1 Figure 2.3: A schematic of the growth of the web during one-sided dewatering [1]. available correlations for these parameters assume isotropic material behavior and neglect plastic deformation and time-dependent behavior. 2.2 Papermaking In general the models used to describe the dewatering process during papermaking are remarkably similar and are summarized below. They began by considering the passage of two permeable wires over a curved surface with radius R as shown in Fig. 2.2. The fabrics move under tension T and at a constant velocity U. A fibrous suspensionflowsbetween the two fabrics. It is assumed that the fabric and suspension move at nearly the same velocity far upstream. The formation of the paper web is shown schematically in Fig. 2.3. The CHAPTER 2. LITERATURE REVIEW 14 coordinate system is denned with the origin at the point where the lower fabric touches the curved surface. The position of the fabric is denoted by g ( £ ) ; the interface between the web and the suspension is given by d(£). The web thickness is therefore 5(£) — g ( £ ) — d(£). In these modelling studies, the forming web is considered as a separate domain [25, 72]. For clarity, two domains, Q i and Q2, are usually defined to denote the regions describing the flowing suspension and the paper web respectively. In Q , the model equations, in l their most general form, are based upon the following assumptions: i) The flow of the suspension is considered to be Newtonian, incompressible and irrotational; ii) The wires are pre-stressed with a constant tension, are inertialess, inextensible, ideally flexible and of negligible thickness. The incompressibility condition and the assumption that the flow is irrotational can be expressed by a Laplace equation for the stream function \& V * =0 (2.18) 2 The suspension velocity u* is related to the stream function as: = ( 2 1 9 ) where the indices x and y represent differentiation with respect to x and y. Based upon a force balance, the curvature of the fabric is related to the pressure drop across the fabric (AP*) and the applied tension T. From classical mechanics this is given by T AP* = - (2.20) where r\ is the local radius of curvature of the fabric. The pressure drop across the fabric can be estimated from the Bernoulli equation as AP* = ^ - | ( ^ + * ) = ^ 2 where U is the upstream velocity and p is the density of the fluid. (2.21) CHAPTER 2. LITERATURE REVIEW 15 In 0 , many authors advance the argument that the drainage resistance of the fiber 2 mat varies in the machine direction. To help address this issue, Moch [73] assumed an empirical relation. Green and Kerekes [74] assumed that the web grows linearly with distance. Martinez [72], on the other hand, derived a diffusive-type differential equation for the solidity 4> in the web using the equations of motion as a starting point. This is given by (2.22) where D((p) is a function of the compressibility and permeability of the network; and </(£) is the dewatering rate. In this derivation, Martinez assumed that the fibres and water move at essentially the same velocity in the machine direction. By doing so, Martinez neglects the viscous stress in his formulation. Zahrai et al. [75] presented the only numerical solution of the coupled problem between domains Cl and f^.' x Finally, although these models have engineering value, and can predict the pressure acting along the length of the paper machine and the subsequent dewatering rate, they give no insight into the development of the physical properties of paper. Jansson [76], for example, demonstrated that fibres are randomly oriented on the surface and in the center of paper but become more preferentially oriented away from these three planes of the sheet. Similarly, the optical appearance of the uniformity, traditionally termed the formation, appears less grainy under certain drainage conditions. Both these phenomena can not be explained by the mathematical formulation given above as the conditions between the suspension and the fibre mat have not been adequately considered in the previous formulations. 2.3 Objectives of the Research The objectives of this work are: 1. To establish a model for the flow of fluid both through and above a porous medium. 2. To develop a robust algorithm to solve the governing equations. CHAPTER 2. LITERATURE REVIEW 16 3. To verify the model experimentally by measuring velocity profiles in a laboratory device. Chapter 3 MODEL EQUATIONS 3.1 Introduction The geometry considered is shown in Fig. 3.1. The channel consists of a fluid (without any fibers) bounded below by an isotropic non-deformable porous fiber mat, of permeability K and porosity e, and bounded above by a rigid impermeable plate. The total height of the channel is H and the thickness of the clear fluid region is h. At the bottom of the porous material (the porous wall), the fluid is withdrawn or injected at a uniform velocity U . w It is to be noted that the case considered here, although initiated from a pulp and paper application, also occurs in a wide range of fluid processes such as membrane separation, journal bearing lubrication, biological transport processes (blood flow), and solidification of metal alloy melts during processing. The underlying phenomena and the interactions between the clear fluid and the porous system are diverse and complex. Before presenting the model equations, it is instructive to examine the equations of motion for flow in a porous medium, see for instance the review of Nield and Bejan [77]. Based upon the observations of steady unidirectionalflowin uniform porous media, Darcy [78] suggested a linear relationship between velocity of discharge and the pressure gradient in the fluid. In vector notation, this relationship can be expressed as V P * = —^u* K (3.1) where u* is the so-called Darcy velocity or seepage velocity (the average of the fluid velocity over a volume element containing both the fluid and solid matrix); Darcy's law is limited 17 CHAPTER 3. MODEL EQUATIONS 18 0. - x h V H ogooooooooooooooooooooooooooooooooooooooooooooooooooo^ oooogoooooooooooooooooooooooooooooooooooooooooooooooo oogoooogoogoooooooooooooooooooooooooooooooooooooooooo gggggggooggoooooooooooooooooooooooooooooooooooooooooo OOO POOOPOOO OpoQQOPPPQOOOOPOPOOOOOOOOPOOPOOOOOPPPOQPIP Q l l l l l l j l l l l l i l l l l l i w y Figure 3.1: The channel flow considered. to flow where the Reynolds number, based upon the average grain or pore diameter, is of order unity or smaller. It is widely known that it is difficult to apply Darcy's law near the interface between a highly porous medium and a viscous fluid. The difficulty in this case is that in the clear fluid region, the flow is governed by the Navier-Stokes equations that are secondorder differential equations, whereas the Darcy equation is only first-order. An intuitive means of resolving this difficulty is to modify the Darcy's law relationship by including a second-order viscous term. Brinkman [79] was the first to propose this modification and the corresponding equation is known as the Brinkman-extended Darcy equation. Wooding [80] has extended this equation further by adding a convective acceleration term. He does this as the curvature of the streamlines in the interface area may be of the order of the pore or particle diameter and this may lead to non-zero values of the convective acceleration term (u* • \7)u*. In our notation, the equations of motion for steady incompressible flow in a porous medium are given by V-u* ^(tT-V)tT e = 0 (3.2) = - e V P * + e/x W e - e-^vT (3.3) K These equations closely resemble the form of the Navier-Stokes equations. The coefficient jj, is the "effective viscosity" of the porous medium and is not, in general, the same as e the fluid viscosity. Ochoa-Tapia and Whitaker [81, 82] suggested that fj, = p/e. For e CHAPTER 3. MODEL 19 EQUATIONS cases where high local velocities prevail in the porous medium, i.e., when the pore-scale Reynolds number is increased in the range 1 to 10, the linearity of the Darcy surface drag breaks down. Experimental results are available that show the transition to nonlinearity [83]. This transition has been characterized by adding a quadratic term in the right-hand side of the above equation. In vector notation, this term, is usually referred to as the Forchheimer term, and is expressed as pC'fK~ / \u*\u*, l 2 where Cj is a dimensionless form drag constant [77, 83, 84]. This has been studied by Lage [85] and Manole and Lage [86] and they advance the argument that the convective inertial term is normally small when compared to the Forchheimer term when Reynolds number is higher than 1. This term is not considered in this work. At the porous medium /fluid interface, it has been experimentally confirmed by Beavers and Joseph [87, 88] that when a viscous fluid flows past the surface of a porous medium, the effects of viscous shear in the clear fluid will penetrate beneath the permeable surface to form what is effectively a boundary layer region in the porous medium. The Darcy equation is clearly not compatible with the existence of a boundary layer region in the porous medium because no macroscopic shear term is associated with the equation [89]. The work of Beavers and Joseph [87, 88] was one of the first to address this issue. Beavers and Joseph circumvent this fundamental inadequacy of the Darcy equation near the interface region by proposing an empirical mathematical model in which the boundary layer could be approximated by a slip velocity, i.e., the velocity across the interface need not be continuous. The magnitude of this slip was postulated to be directly proportional to the prevailing shear stress, namely •d dy* (3.4) where v* and u* denote the x-component of velocity within the channel and porous medium, respectively; and i9 is a dimensionless constant experimentally determined to range between 0.1 < $ < 4. Some theoretical support for this condition is provided by Taylor [90] and Richardson [91] based upon an analogous model of a porous medium, and by the statistical treatments of Saffman [92]. Neale and Nader [89], as well as Vafai and Thiyagaraja [93] and Vafai and Kim [94], CHAPTER 3. MODEL EQUATIONS 20 have questioned the validity of the "slip velocity", i.e., there cannot physically .be any slip between the fluid in the channel and the fluid within the porous medium in an actual velocity profile. They proposed that both the velocity and shear stress must be continuous at the interface, i.e. u* = v* Out dv* x (3.5) x = < - > 3 6 To incorporate this extra boundary condition, these authors employed the Brinkmanextended Darcy model to describe the flow through the porous medium. Recently, OchoaTapia and Whitaker [81] utilized a sophisticated volume averaging technique to extend this formulation. They recognized that part of the viscous stress in the fluid region is taken by the porous matrix in porous region across the interface, i.e., the Darcy term related stress. They deduced the form of this Darcy stress (so called the stress jump) from a rigorous stress function analysis for the fluid and porous regions. In their study the jump in shear stress is inversely proportional to the square root of the permeability of the porous medium, < ldvl_dyl e dy* dy* = v* (3.7) x = *JK X where j3 is an empirical constant. Ochoa-Tapia and Whitaker derived that the value of (5 is of the order of the ratio of the thickness of the boundary region near the interface to the square root of the permeability, or its inverse, therefore should be 0(1), and it may be positive or negative as the difference in stress functions can be positive or negative. It follows that Eq. (3.6) is a special case of Eq. (3.8) when (5 = 0 and p = p/e. Ochoae Tapia and Whitaker's [81] proposed interface conditions have been adopted recently by Kuznetsov [95, 96, 97], and Alazami and Vafai [98]. Ochoa-Tapia and Whitaker [99] also includeed a quadratic jump condition /^pu , for flow where inertial effects are evident. This 2 term is not considered in this work. CHAPTER 3.2 3. MODEL EQUATIONS 21 Model Equations A Cartesian co-ordinate system (x*, y*) is used to describe the channel shown in Fig. 3.1. Here x* measures the distance along the length of the channel and y* 6(0, H). The assumptions made in this analysis are as follows: 1. The variation of quantities in the span-wise direction are neglected and the problem is considered two-dimensional; 2. The fluid is Newtonian; 3. The fluid is viscous and the flow laminar; 4. The fluid is incompressible and the body forces are negligible; 5. The porosity e is uniform; 6. The governing equation for the fluid region is the Navier-Stokes equation; the Brinkman extended Darcy's law with the inertial terms retained is used in the porous region [77, 100]; 7. The interface boundary conditions used are those outlined by Ochoa-Tapia and Whitaker [81]. When the independent variables are scaled with H, the flow components in each domain by U , and the pressure with pU , the equations of motion in dimensionless form reduce 2 w to V • v =0 f +(v-V)v t 2 V -u = 0 § + l(u.V)u = V(x,0<y<x) = - V P + ^ \7 v V ( x , x < y < 1) (3.9) - e V P + ^ - ^ where v and u represent the flow velocities in the clear fluid and porous regions respectively; x and y represent the channel coordinates (dimensionless); Da is the Darcy number defined by Da = K/H ; 2 Re is the Reynolds number (Re = U Hjv)\ w and x 1S the interface CHAPTER 3. MODEL EQUATIONS 22 position defined as x = h/H. The boundary conditions at the upper and lower walls are straightforward and given by v(x,0,t) = (0,0) 1 ' > u(x,l,t) = (0,l) J V (3.10) i.e. no slip. At y = x, we use the conditions outlined by Ochoa-Tapia and Whitaker [81]. In non-dimensional form, these are given by v = u (3.11) ldux dv (3 e dy dy ^Da x T+ • = r_ (3.13) where f3 is an empirical constant [81] which must be determined independently; and the subscripts + and — refer to the position above and below the interface, respectively. Equation (3.13) indicates an equilibrium of the normal stresses at the interface. Equations (3.9) through (3.13) represent the model equations. Chapter 4 SIMILARITY SOLUTIONS Before a solution to the model equations is given, it is instructive to briefly review other geometries which resemble the flow considered. This will aid in understanding the form of the similarity ansatz considered in this thesis. To begin, this literature is categorized into two groups. The major works in category (I) [81, 82, 87, 89] are fully-developed flows with U = 0. In general, these works focus on understanding the underlying phenomena and w the interactions between the free fluid and the porous system. This literature has been reviewed in Chapter 3. In the second category (II) are flows with h = H. This body of knowledge begins with the classic work of Berman [3] who investigated the problem of flow in a channel with both channel walls equally permeable. Berman assumed that the solutions were symmetrical about the centreline of the channel and reduced the two-dimensional NavierStokes equation to a fourth-order non-linear ordinary differential equation using a similarity transform. This equation depends upon a single non-dimensional number, the transverse Reynolds number Re defined in terms of the channel width and suction velocity. Similarity solutions were tested in early experiments, where evidence of good agreement between measurements and theoretical predictions are reported [101]. Many authors have extended Berman's symmetric series solution [102, 103, 104, 105, 106, 107, 108], and a very rich structure of solutions is evident in this literature. Recently, Cox and his co-workers [4, 109, 110, 111] considered in detail the asymmetric case of h = H with suction or injection at only one wall, while the other wall is 23 CHAPTER 4. SIMILARITY SOLUTIONS 24 impermeable. After scaling the Navier-Stokes equation by H/2 and U , they formulate w the steady-state problem using the stream function fy of a form originally suggested by Berman, i.e. fy(x,y) = —xf(y). After elimination of pressure from the Navier-Stokes equation, the separable component f(y) must satisfy the following equation /"" + Re(f'f" - ff") = 0 (4.1) subject to the boundary conditions / ( - ! ) = / ' ( - 1 ) = /'(1) = 0 /(1) = 1. (4.2) > Here y = — 1 represents the impermeable wall and Re = U H/2u. w Cox has shown that there is a unique steady solution for all values of Re, except in the range 7.05 < Re < 7.31 where there are two stable solutions. In this chapter, we build on the work of Cox and his co-workers [4, 109, 110, 111], and extend their work by considering coupled-flow both through and above a porous medium (i.e., h H) by extending the similarity ansatz. It is often found that similarity solutions lead to regions of no solutions or multiple solutions. Keeping in mind the limitations of similarity solutions as outlined by Brady and Acrivos [112], we present in §4.1 a general similarity formulation that depends on a number of independent dimensionless parameters. In §4.2, the steady-state solutions are considered in which both low-Reynolds number expansions and numerical solutions are presented. A linear stability analysis is conducted using the method of normal modes in §4.3. Concluding remarks are given in §4.4. 4.1 Similarity Formulation We reduce our model to something more tractable than Eqs. (3.9)-(3.13) by assuming the form of the solution. Following the procedure given by Berman [3], we postulate that the stream function fy(x,y,t) may be written as a product V(t,x,0<y< ) X V{t,x, <y X < 1) (4.3) CHAPTER 4. SIMILARITY 25 SOLUTIONS where V denotes the average (dimensionless) axial velocity at x = 0; and / and g are dimensionless time dependent functions which will be subsequently determined. It can be shown from an integral mass balance that V — x represents the average axial velocity in the channel at any downstream location x. The velocity components follow immediately 1 from the definition of \&, i.e. v= ( (V-x)f (y,t)J(y,t)y K u= y ((V - x)g (y, t), g(y, t)) (4.4) y where the subscript y or t denotes the partial derivative with respect to y or time t, respectively. After elimination of pressure from the equations of motion by cross differentiation, we find that the functions / and g must satisfy fyyt 9yyt fyyyy (^fyfyy R • (4-5) ^ ^ \9y9yy ~ 99yyy) ~ r ) J l ^ lf^ = ffyyy^j vvyy a l^-OJ y y e The boundary conditions at the upper and lower surfaces are /(0,i) = 0, f (0,t)=0, y g(l,t) = l, g (l,t) [fyy + , = 0. y (4.7) and at the interface y = Xy 9 = f, J_f Jvw Re 9y = fy, 9yy = £ _ f f 4. f _ f _ 9yyy JTvv + Jy h t 2 e R e 99yy 9 2 + 2 (-) 4 g y £ £ 1. 9yt- y D a R e 8 , , l^.Jj € Equations (4.5) through (4.9) represent a considerable reduction in the complexity of the model equations. It should be noted that the derivation of Eq. (4.9) is difficult and details of this are given in Appendix A. 4.2 Steady-State Solutions In this section we describe the steady solutions for f(y) and g(y). In §4.2.1 we present some asymptotic results for Re —> 0 using regular perturbation methods. By doing so, we 1 Assume that the the dimensional average axial velocity at position (x*,y*,t*) is v*, the mass balance show that V*H = v*H + U x*, hence v* = V* — U x*/H, which in dimensionless form, is v = V — x. w w x CHAPTER 4. SIMILARITY SOLUTIONS 26 provide results about the basic flow and show that the solution diverges when Da/e —» 0. In §4.2.2 we have also obtained solutions for small Re using a power series approximation. The prerequisite of this method, however, depends on an appropriate re-scaling of the original equations. In §4.2.3 we describe a numerical technique for solving the equations and show that, at large Re, the solution bifurcates. 4.2.1 Asymptotic Solutions To begin, we seek a steady solution of the similarity equations of the form oo f(y) = £ Re f (y) n n (4.10) g{y) = £ Re g (y). n n . n=0 With this, the leading terms in these series are readily found to be fo(y) = \C,y 2 + |CV, + Cy + g (y) =C 0 3 A C e ^ +C e~ ^ y 5 % (4.11) where C are the constants of integration. It should be noted that we have applied the l boundary conditions at the upper wall. The point of this exercise is that upon examination of the leading terms, we would anticipate difficulties in the limit of Da/e —> 0. Clearly, under this condition the form of the solution dictates that C must be zero; otherwise, g 5 would grow exponentially with y making it unmatchable with the boundary conditions at y = 1. In the limit of Da/e —»• 0, it can be shown that, after a modest number of algebraic manipulations, the flow field is given by /o(y) = 3 ( | ) - 2 ( | ) 2 When the magnitude of e^ ^ c Da 3 and g (y) = 1 0 (4.12) is not large, the constants of integration are deter- mined from the remaining boundary conditions. This leads to a linear system of algebraic CHAPTER 4. SIMILARITY SOLUTIONS 27 equations, i.e. X /2 X /6 -1 -x X X /2 0 -1 e ?X 0 0 1 0 0 0 0 1 X 0 0 0 1 2 3 2 0 c -5 B 1 e 2 e 2 0 3 0 c 4 0 c 5 1 c (4.13) 0 where intermediate constants in the above matrix are defined for brevity, i.e. <7 = y/e/Da, Bi = B = ( ^ 2 (V~e/3 - l)cV*, + i)c e-«. 2 This system of algebraic equations can be solved using standard techniques and the leading order estimates of the x component of the flow field, represented by f and g are given y y in Fig. 4.1 (i). Here, it is shown that the velocity in the clear fluid region joins that of flow in the porous medium across the interface region. The curvature of the profile at the interface decreases with increasing Da/e. Due to the momentum of the fluid in the clear fluid region, fluid penetrates into the porous medium and decays rapidly to zero within the porous wall for a low permeability medium. We term the depth of the fluid penetration as the thickness of a "boundary layer", i.e., the thickness from the interface position to where the x component velocity g (y) is 0.01. This thickness is normalized with the total y thickness of the porous medium and shown in Fig. 4.1(h). Clearly with Da/e > 10~ , the 3 porous medium offers very little resistance to the fluid's motion, and the fluid penetrates at least 90 percent into the porous region. When Da/e —> oo the flow in the channel behaves like a fully-developed Poiseuille flow. On the other hand, with Da/e —>• 0, there is no longitudinal velocity in the porous medium and the boundary layer diminishes to zero. It is interesting at this point to examine the solution at higher Re and we do so by computing one further term in the expansion. This was performed in Maple, a symbolic CHAPTER 4. SIMILARITY SOLUTIONS 10" 5 10" 28 10" 4 3 Da/e 10" 10' 10" 2 1 Figure 4.1: (i): The x-component of the velocity profiles f y and g (normalized by the y averaged velocity V - x at position x) estimated in the limit of Re —> 0. Four curves are shown: (a) Da/e -> 0, (b) Da/e = 5x 1(T , (c) Da/e = 1 x lfT , and (d) Da/e -> 0 0 . (ii): 4 2 An estimate of the boundary layer thickness in the porous medium. The boundary layer thickness has been normalized to the thickness of the porous medium. These estimates were made with ,6 = 0.2, x = 0.8 and e = 0.9. CHAPTER 4. SIMILARITY 0 29 SOLUTIONS 0.5 1 Figure 4.2: The effect of Re for the case Da = 1 x 10 profiles, f and g , here (3 = 0.2, x — 0.8, e = 0.9. y y 1.5 Velocity (dimensionless) 3 2 on the the x component velocity CHAPTER 4. SIMILARITY SOLUTIONS 30 solver, and the resulting algebraic expressions to the general case are too large to be reproduced here. To illustrate this point, the leading order and the first order approximations are given in Appendix B for one particular case, namely with Da = 1 x 10~ , x ~ 0-8, 3 e = 0.9 and P = 0.2. The effect of Re is also shown in Fig. 4.2. Here we see that, with increasing Re, the location of the maximum longitudinal velocity in the channel shifts towards the porous wall. As we proceed with higher order approximations, we found that the exponential form of the solution within the porous region makes the explicit solutions difficult. 4.2.2 Series Solutions Equations (4.5)-(4.9) are difficult to solve over wider range values of Re and Da. When Da is low, the equations become stiff and consequently difficult to be solved directly. The solution method considered here relies crucially on re-scaling f(y), g(y) and y so that the governing equations do not contain the Reynolds number and Darcy number. If V = z/q f(y)=q\{z)/Re g{y) = eqip{z) / Re (4.14) where q = \J'e/Da, the governing equations reduce to Xzzzz = XX Ipzzzz Iplpzzz-lpzlpzz ZZZ = (4-15) ~ XX Z ZZ (4.16) + lpzz subject to A(0) = 0, A (0) = 0, 2 \( ) xq = eip(x<l), K{X0) = ipzzz(XQ) = eip ( q), z Kzz m = ^, eq X ~ AA + \\ + tpip - ip + ip 2 2Z zz z z CHAPTER 4. SIMILARITY SOLUTIONS A(q) 31 = 0. (4.17) If we seek a solution of the form oo X(z) = ^ a z (4.18) n n n=0 oo ib(z) = Y. ^ h (- ) n 4 19 n=0 where z = z — xQ- When the boundary conditions at the impermeable wall are evaluated, the first few terms of these series are given by ib{ ) z = \c z = c + cz + where C ,C\,C ,C , 0 2 3 + \c z*-±& z* 2 o x + 0(z*) 0 (4.20) lcrj + ^-(c c -c c + c yz 3 + lc,z + 2 l 0 24 ~(C C -C C C + C C + C -C )Z 2 2 3 2 2 2 5 2 3 4 2 4 5 5 4 4 5 3 + 4 i 0(z ) 6 C and C are constants which need to be determined via the remain4 5 ing six boundary conditions using Newton's method, i.e., a method described in [113]. A representative solution is given in Fig. 4.3. We find that the series solution yielded sensible results in the region [Re, Da] = [—5, 5] x [10~ ,1]. 4 4.2.3 Numerical Solutions Both solutions presented thus far are useful when — 5 < Re < 5, where a negative Reynolds number indicates an injection case, due to the sign change of the velocity at the bottom wall. Outside this region, a higher order approximation is required. As a result we pursue a numerical solution for a wider range of Darcy and Reynolds numbers. In general, we recognize that this is a three-point boundary value problem. We developed two different numerical methods in an attempt to solve this problem. The first algorithm is based on a combination of shooting and Newton's methods [113]. This method is simple and does not require good initial guess of the solution. This method, however, was found useful only in a range slightly larger than the series solution presented in the previous section. As a result, we developed a second algorithm, based upon the work of-Terrill [104], which we call the "reverse solver". The algorithm was found to be robust for most values of Re and Da. Both algorithms will be discussed below. CHAPTER 4. SIMILARITY 0 Figure 4.3: The 6 th SOLUTIONS 0.5 1 1.5 Velocity (dimensionless) 2 order series solution (-) compared to the similarity numerical solut' (•) with Da = 1 x IO" , Re = 5,P = 0.2, 3 x = 0.9. CHAPTER 4. SIMILARITY SOLUTIONS 33 Method 1 - The Shooting Method The first algorithm considers the steady case only, where all the non-dimensional parameters (e, X) P, Da and Re) are specified. The numerical procedure is as follows: 1. It is useful to note that Eq. (4.5) may be integrated once in y to give (4.21) where Ai is the constant of integration (which is unknown at this point). In physical terms f (0) yyy = A . Further, we define a second constant f (0) x yy = A which will be 2 used subsequently. 2. If we assume trial values for A\ and A , we can treat Eq. (4.21) as an initial value 2 problem and integrate this relationship with a fourth-order embedded Runge-Kutta solver over the domain 0 < y < x [H4]3. At y = x the interface conditions are evaluated using Eqs. (4.8) and (4.9) for a prescribed value of (3. 2 4. Equation (4.6) is now integrated over the domain x < y < 1 using the same RungeKutta method. 5. The solution with the trial values of A and A are compared to the boundary x 2 conditions at y = 1. Newton's method, a root finding procedure [113], is used to update the values of A and A until the boundary conditions at y = 1 are satisfied x 2 to within a tolerance of 1 x 10~ . 6 Results from this technique yield very good agreement with the solutions obtained from the second numerical method for —60 < Re < 16. Method 2 - Reverse Solver The "reverse solver" considered here is designed to solve the re-scaled Eqs. (4.15) to (4.17). By "reverse", we mean that, in the traditional manner of solution, the parameters 2 A method to determine (3 will be considered in Chapter 5. CHAPTER 4. SIMILARITY SOLUTIONS 34 Re and Da are specified and the equations are integrated. The reverse approach involves a transformation of the governing equations and constructing plausible solutions without specifying Re and Da beforehand. Once a solution is constructed, the corresponding Re and Da are determined. The algorithm is summarized below: 1. Equation (4.15) may be integrated once in z to give (4.22) where A is the constant of integration. In physical terms X (0) 3 = A . At this zzz 3 point, we define a second constant A^(0) — A . 4 2. We begin the numerical procedure by choosing trial values of q, A and A . 3 3 4 3. With known values of A and A , Eq. (4.22) may be solved as an initial value problem 3 4 using a fourth-order Runge-Kutta method over the domain 0 < z < qx4. At z = qx the interface conditions are evaluated using Eq. (4.17) for a prescribed value of p. 5. Equation (4.16) may now be integrated over the domain qx < z < q as an initial value problem using the same method mentioned above. 6. At this point we compare the computed value of ip(q) and ip (q) with the boundary z conditions given in Eq. (4.17). We use Brent's method [115], a root finding procedure to update the final value of q until we find the point where ip (q) — 0. It should be z noted that under certain combinations of A and A , we can find two roots. Each 3 4 root represents a solution to the problem as discussed by Cox and King [111]. 7. The solutions with the trial values of A and A are now compared to the targeted 3 4 Darcy and Reynolds numbers and a similar iterative procedure is then invoked until convergence is achieved. • The initial guess was estimated from the results of the shooting method. A linear/nonlinear (e.g., 3 "spline" method) extrapolation was used to estimate A 3 and A . 4 The technique used is similar to a branch tracing method commonly used in solving nonlinear equations. CHAPTER 4. SIMILARITY SOLUTIONS 35 With this method we are able to find solutions for a wide range of Reynolds numbers and Darcy numbers. We validated our solution scheme with solutions from (i) King and Cox [4] for the case when Da is high (e.g., Da — 1000, see Fig. C.l), (ii) the asymptotic expansion in the region [Re, Da] = [—4.5, 4.5] x [10 ,1] (see e.g., Fig. B.l) and (iii) the series solution -4 when [Re, Da] = [—5, 5] x [10 ,1] (see Fig. 4.5, for example). Good agreement was found -4 for all cases. In fact, we were able to reproduce identically the bifurcation reported by King and Cox [4], A series of numerical experiments were conducted over the domain [Re, Da] = [0,170] x [0.013, 0.05] as this represents the region of interest for papermaking applications [116, 117]. In general, three types of solution behaviors were found. • type I - no backflow, • type II - backflow in the channel, no backflow in the porous region, • type III - backflow in both domains. Examples of these are shown in Fig. 4.4(a)-(c). The interesting feature here is that under certain conditions, the solutions found were not unique. In other words, we found multiple steady states for certain combinations of Re and Da. To help illustrate this behavior, a plot of the viscous stress ( — f (x)) yy 4 is given in Fig. 4.5. In this figure, the shear stress is shown as a function of Re at a constant Da. In regions where multiple steady states were found, there exists multiple values for the shear stress. As shown, multiple solutions were found over the entire domain with the upper branch displaying a turning point bifurcation in the region 16.541 < Re < 17.169 and at Re = 54.722. The number and the stability of the solutions depend upon both Re and Da. The stability of each branch will be discussed shortly. Included in this figure are the regions in which each type of solution exists. Part of the results for all simulations conducted is shown in Fig. 4.6 as a surface plot. To help simplify this figure only type-I and II solutions are included. We do so as type-Ill solutions are mostly unstable. 4 F o r a two-dimensional flow, the shear stress can be expressed as r* = M(§^ + g^")- When T* is scaled with pU^, the non-dimensional form of the shear stress becomes r = ^ ( ^ ^ + ^jr)- Substitution of Eq. (4.4) yields T(X) = T&(V — )fyy(x) x shear stress at the interface is —f (x)yy a t * n e interface. Hence, a reasonable representation of the CHAPTER 4. SIMILARITY SOLUTIONS 36 Type I x component ( f ) \ y component (f) \ \ y component ( g ) \ _ •V \ x component (g ) \ \ \ \ \ \ \ \ P\ \ 0.5 1 Velocity 1.5 Type II ( / v^x component (f ) . y component (f)\ \ -2 1 ^M\\\ \ \ y component (q)\ \ \ \ \ \ \.\ \ z i 0 \/ i 2 *± Velocity 0 Type III i.2 i.4 i.6 ^ x component (f ) y component (f)\ i.8 -"Vx^oTnponent (q ) \ \ \\ \ \ \ \ p 30 -25 0 Velocity 25 r<\\ 50 Figure 4.4: Example solutions: The plot represents the three different families of solutions obtained, here (a): Re = 10.93, (b): Re = 16.68 and (c): Re = 32.21. These simulations were conducted at Da = 0.05, (3 — 0.2, x = 0.8 and e = 0.9. The shaded area represents the porous region. CHAPTER 4. SIMILARITY SOLUTIONS 37 Figure 4.5: Example solutions: The figure represents the shear stress at the interface — fyy(x) (here the circles indicate the series solutions given in §4.2.2). These simulations were conducted at Da = 0.05, p = 0.2, x = 0.8 and e = 0.9. The type of solution is indicated, where the line "|" shows the transitions from type-I to type-II and from type-II to type-Ill solutions. Here the solid-lines indicate stable solutions while the dotted region represents unstable solutions. CHAPTER 4.3 4. SIMILARITY 39 SOLUTIONS Stability of the Similarity Solutions The temporal stability of the steady solutions shown in the previous section is analyzed by introducing small amplitude time-dependent perturbations of the same similarity form suggested in Eq. (4.4). The method we used is called the method of normal modes, which has also been used in the work of [107] and [109]. We consider the solutions are of the form: f(y,t) = f(y) + e° F(y) (4.23) 9(y,t) = g(y) + e° G(y) (4.24) t t where a = o + io~i is complex. Upon substitution of these relationships into Eqs. (4.5)r (4.9), after linearization, the following characteristic value problem for the perturbations F(y) and G(y) is obtained. Fyyyy e Gyyyy if^yyy fyyy^ fy^yy fyy-^y) oReFyy (4.25) 6 (d^yyy 9yyyG 9y^yy 9yy^*y) oReGyy (4.26) Re ~~Eo, Gyy with F{0) = F (0) = G(l) = G (l) = 0 v (4.27) y and F = G (4.28) F = G (4.29) y y + -^=F )j - oF = ^ - ^G ~ D ^ e Gyy = e^F yy (^ f-fFyy-fyyF) R + 2f F y ~e^ 9yyG y y G y (4.30) x 9 ~ l (4.31) yy 'i~ 1 G y + 2 Z at the interface y = X- These equations form an eigenvalue problem. For the suction problem, a given solution is stable if for every eigenvalue a < 0 is obtained. In the case r of injection, instead, all eigenvalues must have a positive real part for the solution to be CHAPTER 4. SIMILARITY 40 SOLUTIONS stable. This is due to the sign change in the reference time (h/U ), w since U < 0 in this w case. These equations were solved using a root finding procedure similar to that outlined by Durlofsky and Brady [106] and Zaturska et al. [107], that is, a shooting method combined with Newton's method. A representative example is shown in Fig. 4.7. The shear stress is shown in the top panel of the figure. There are two regions of interest where multiple solutions were found on each branch. These are labelled Region I and Region II. Region I is magnified in Fig. 4.7(h) and, for this region, a is shown in Fig. 4.7(iii). Clearly, for r the dotted region in the range 19.731 < Re < 19.833, a > 0 and, hence, this represents r an unstable solution. The analysis was conducted over the entire domain. Solutions which were found to be stable are presented as a solid line in Fig. 4.7(i). A summary of the solution behavior is shown in Fig. 4.8 for all Re and Da tested. Included with this is a supplemental table (Table 4.1) detailing the nature, number and stability of each solution region. The first observation that can be made from this figure is that multiple solutions exist for all Da and Re considered. Unique stable solutions were found in regions A, B, D, G and H. Regions with multiple stable solutions were generally type-II solutions. These were found in regions C, E and F. No solutions could be found in region I. To help illustrate the behavior as one crosses over a boundary in this diagram, the shear stress at three different Darcy numbers are shown in Figs. 4.9(a)-(c). In each figure, multiple solution regions labelled A, B, etc. have different type, stability and number of solutions. 4.4 Summary In this section we solved the equations of motion for flow in a channel partially filled with a porous medium using a similarity ansatz. Our numerical analysis has been fairly complete - we have applied regular perturbation techniques, power series approximations and a shooting method for small Re, together with a numerical algorithm for higher Reynolds number solutions. Although not stated in the text, we have attempted to approximate the solution using singular perturbation methods in the limit of Re —> —oo and Re —> + 0 0 but Figure 4.7: The stability analysis for Da = 0.013: (i) is the steady solution, the detailed view in Region I is shown in (ii), whilst (iii) shows the real part of the corresponding eigenvalues in Region I. Here the dashed lines indicate unstable solutions. CHAPTER, 4. SIMILARITY 42 SOLUTIONS (a) 0.05 —i 1—• 1 1 1 r 0.045 0.04 0.035 Q D I 0.03 0.025 0.02 E F 0.015 20 40 60 80 100 120 140 160 Re Figure 4.8: (a) Summary of solution behavior for simulations with (3 = 0.2, x — 0-8 e a n d = 0.9. Nine different solution regions are found and the detailed information about each region is tabulated in Table 4.1. (b) Region C is enlarged. CHAPTER 4. SIMILARITY 43 SOLUTIONS Table 4.1: The number, type and stability of the solutions for each solution region as shown in Fig. 4.8. Type of Solutions No. soln . a III II I Region No. stable 6 No. soln. No. stable No. soln. No. stable A 1 1 . 0 0 1 0 B 0 0 1 1 1 0 C 0 0 3 2 1 0 D 0 0 1 1 1 0 E 0 0 3 2 1 0 F 0 0 2 1 2 1 G 0 0 0 0 2 1 H 0 0 2 1 0 0 I 0 0 0 0 0 0 a Total number of solutions at given Da and Re. b Total number of stable solutions at given Da and Re. CHAPTER 4. SIMILARITY 44 SOLUTIONS (a) 20 II 10 i y° ~" ,--'\ -10 -20 ' in .c A H F ; D E 1 50 10 10 20 I !! Re ^30 150 100 40 50 60 Figure 4.9: Example solutions: the shear stress function —f (x) versus Re for (a) Da = yy 0.013, (b) Da = 0.026 and (c) Da = 0.0375. The type of the solutions are denoted, and the multiple solution regions are labelled A, B, etc. representing different type, stability and number of solutions. Here solid-line denotes stable solution and dotted-line, the unstable solution. CHAPTER 4. SIMILARITY SOLUTIONS 45 were unsuccessful. Our results at moderate Reynolds numbers are worth highlighting, as we observed bifurcation. The temporal stability of the steady solutions for the problem has been analyzed using the method of normal modes. In this analysis, we represent the time dependent solutions as a sum of a base steady solution and a small time disturbance of similarity form. After linearization, we show that the system of equations is an eigenvalue problem. A shooting method combined with Newton's root finding procedure was employed to solve the corresponding equations. Nine discrete solution regions were found. In each region, the number of the solutions, as well as the type and stability of these solutions were summarized. It should be noted that, the similarity variable we selected for the similarity approach can not be extended for models which include the Forchheimer term and the nonlinear shear stress jump term in the interface condition, because, these two terms are quadratic which would leave the transformed differential equations dependent on the x component velocity. Finally, we use the above results to estimate the shear stress at the interface of a twin-wire machine. As has been discussed, the shear stress at the interface at position (x, y) can be represented as a product of the channel average velocity at x and the second derivative of the function / , which depends only on the channel position y, i.e., by T*(X) = jfe(V — )f y{x)pUwx y Here V is the ratio of the inlet mean velocity of the suspension jet to the mean suction velocity. This expression shows that a high dewatering Re reduces the shear activity near the interface and hence a higher jet wire speed difference is required to create more shear. At the early stage of the dewatering, when the formed paper web is very permeable, there may exist a certain suction Re where no shear occurs at the interface, no matter what the jet/wire speed ratio is, as indicated in Fig. 4.5. On the other hand, when the web formed is less permeable, there seems to exist a wide region where high shear occurs at the interface region, as shown in Fig. 4.7(i). It is quite illustrative that, for any suction higher than Re — 16, there could be a backflow in the channel, which means flow recirculations that may contribute to fiocculation and deterioration of paper formation. Nevertheless, we roughly estimate from our available information that, the shear stress at the interface will be between 0 and 0.33 kPa for a twin-wire machine running at 1000 CHAPTER 4. SIMILARITY SOLUTIONS m/min with a jet/wire velocity ratio of 1.1 and a mean drainage velocity of 20 Chapter 5 ONE D O M A I N P R O B L E M 5.1 Introduction In Chapter 4, we divided the channel into two domains composed of an unobstructed fluid and a porous region, and solved the coupled problem by matching the constraints at the interface for an arbitrary given value of p. In this chapter, we consider that these two domains can be collapsed into one composite domain by re-posing the governing equations. By doing so, the interface conditions need not be considered. In this approach, a finite volume computational fluid dynamic (CFD) method was employed to directly solve the original steady transport equations for a finite length channel. To distinguish the two solution methods we employed, the method described in Chapter 4 is called the "similarity approach", and the one in Chapter 5 is termed the "one-domain approach". By comparing the output from both methods, we should be able to determine p. 5.2 One-Domain Formulation In the one-domain approach, we re-write the steady equations of motion for a two-dimensional channel as V-w 1 -(w-V)w e = 0 (5.1) = -eVP+-—V w-——w Re 1 Be 2 47 ReDa (5.2) CHAPTER 5. ONE DOMAIN 48 PROBLEM where w is the vector velocity in the channel, ( { 0, in fluid region 1, in porous region 1, in fluid region 0 < e < 1, in porous region We use the following boundary conditions: uniform inlet and fully developed flow at the outlet; no-slip condition on the upper impermeable wall boundary; and uniform suction on the lower porous boundary. The equations of motion were solved using a standard finitevolume formulation based on an artificial compressibility method (where a time derivative term is added in the formulation to facilitate the computing). We use the centered flux evaluations and implicit Euler time advance scheme for the discretization in space and in time respectively. Since the equations are. nonlinear, we apply an approximate factorization method to convert them into two sets of coupled linear equations in the x and y directions (where for the porous region the sink term is added to the right-hand side fluxes and its Jacobian component is added to the left-hand side Jacobian matrix). The transformed sets of equations are essentially a 3 x 3 block tridiagonal problem which can then be readily solved using a block Thomas algorithm. The convergence criterion is defined such that all changes in solution should be less than 1 x 1CT . We use a uniform grid with the number 4 of grid points being 500 in the flow direction and 80 in the normal direction, at which the grid independence is achieved. Details of the algorithm are presented in Appendix E. A numerical study was conducted in which the solutions from the one-domain and two-domain approaches were compared in order to determine the relationship (3(Re, Da). 1 To do so, f3 was determined by minimizing the sum of the squared differences of the component velocities from these two numerical solutions. An example of this is shown in the lower panel of Fig. 5.1. The comparison between these two methods is .based on the assumption that a similarity of the solution profiles is achieved after a certain length of flow. The pressure drop in the x direction AP calculated from the similarity solution 2 X 'The validation of the code is given in [118]. 2 Note that V needs to be set sufficiently large so that a fully-developed flow is achieved to be compared with the similarity solution. Here we set V = 21 for the suction conditions we considered. CHAPTER 5. ONE DOMAIN 49 PROBLEM assumes a near-linear shape when the suction flow is small; it is plotted along with AP 3 X from the one-domain solution. As shown, with x > 4 the pressure drop from each solution is virtually indistinguishable. For this work we compare the flow profiles only when x > 4. At a position x > 4, the pressure drop in y direction AP also shows good agreement, y as both the similarity and the one-domain solutions predict that the pressure drop in y direction occurs mainly in the porous region, see Fig. 5.1. The sensitivity of the solution to (5 is shown using the similarity solution method in Fig. 5.2. Clearly large differences in the x component of the velocity profile and the shear at the interface exist, when ft varies between — 1 < (3 < 1. The best fit (3 values from the comparison of the two solution methods are shown in Fig. 5.3. For the cases presented, (3 was found to increase with decreasing Da and increasing Re. In all cases tested (3 ~ 0(1) and is mostly positive. 5.3 Summary In this chapter we developed a one-domain approach to solve the steady governing equations in a primitive variable form. In this method, the clear fluid domain and the porous region are considered as a composite domain governed by a single transport equation, the Brinkman extended Darcy's law with an inertial term. This equation has the same form as the Navier-Stokes equation, but with a sink term that can be turned on and off by a binary constant B. Using this method, the interface boundary conditions need not be considered, i.e., there is no need to select the unknown empirical constant (3 a prion. The similarity solutions presented in the previous chapters were found to agree very well with the current one-domain solutions for low suction conditions, i.e. Re < 12, where only the type-I similarity solution was validated. We estimated (3 through a curve-fitting approach. The results indicate that P is a function of Re and Da and was found to be 0(1). There are limitations with this approach. It must be noted that, when 12 < Re < 16, it seems that the entrance effect, i.e., the developing length of the flow, increased and the flow did not quite develop to a similarity solution. With Re > 16, the algorithm became unstable in regions where we would expect type-II and type-Ill similarity solutions. 3 T h e derivation of the pressure drop AP , AP X y for the similarity approach is shown in Appendix §A.2. CHAPTER 5. ONE DOMAIN 50 PROBLEM Figure 5.1: The flow fields as estimated by the two different numerical approaches, here AP X is evaluated at the interface y = 0.8 and AP is evaluated at x = 7.0. The porous y medium extends from y = 0.8 to y = 1, with Da = 1 x 10~ , Re = 5 and e = 0.9. In 3 the lower panel of the figure, the fully-developed profile using the one-domain approach at x = 7.0 (shown as a dashed line) is compared to the similarity solution (shown as a solid line) with p = 0.2. CHAPTER 5. ONE DOMAIN 51 PROBLEM Velocity (dimensionless) Figure 5.2: The flow field determined as a function of j3 using the similarity approach. The normalized x-component of velocity (f for fluid region and g for porous region) is given y y in the figure on the left and the y-component (/ for fluid region and g for porous region) is given on the right. The porous medium (shown as the shaded area) extends from y = 0.8 to y = 1. This simulation was performed with Da = 1 x 10 , Re = 5, and e = 0.9. -3 CHAPTER 5. ONE DOMAIN PROBLEM 52 Figure 5.3: The factors affecting the empirical constant j3 in the shear stress jump condition (the tested cases are for x = 0-9, e = 0.9). CHAPTER 5. ONE DOMAIN PROBLEM 53 Finally, it should be noted that one may alternatively compute (5 solely using the one- • domain approach. The slope of the velocity profile near the interface within the fluid region and the porous region could be computed using the velocity data in the computing cells just above and below the interface. The empirical constant j3 could then be obtained from the difference of the computed slope according to the interface boundary condition Eq. (3.12). Chapter 6 EXPERIMENTAL WORK 6.1 Introduction The purpose of the experimental work is to confirm our numerical predictions. In particular we are interested in measuring regions with stable type-I, II or III solutions and regions with multiple steady states. This chapter begins with a description of the experimental apparatus, the methodology of the measurement and the experimental procedure. Typical test results are presented, and then compared to our predictions. 6.2 Experimental Apparatus The experimental channel consists of parallel plexiglass plates separated by a gap of 1.20 cm (see Fig. 6.1 and 6.2). The length and width of the channel are both 12.70 cm. The top plate is impermeable. A large number of holes have been drilled into the lower plate (see Fig. F.6). The open area of the lower plate has been measured to be 70%. Upstream of the channel is a sealed tank which is fed by an elevated reservoir, i.e. the head tank." The downstream tank is open to atmosphere and each tank has a baffle. Detailed drawings are given in Appendix F. In these experiments, loosely woven fabrics, i.e. the porous material, was secured to the bottom of the channel. Two different fabrics (AstenJohnson Inc. (http://www.astenjohnson.com)) were used in this study and their thicknesses and permeabilities are given in the Table 6.1 below. A glycerine-water solution 54 CHAPTER 6. EXPERIMENTAL WORK 55 Table 6.1: The porous fabrics used in the experiment Fabric Thickness (mm) Porosity* Permeability (cfm/ft )* Da A 1.92 0.38 500 2.66 x IO" B 1.01 0.41 305 1.90 x 10~ 2 3 3 ' porosity was determined by measuring the displaced volume of water in a graduated cylinder. * permeability was estimated by Darcy's law using the air permeability data in cubic feet per minute per square feet of mesh, at a 0.5 in. water pressure drop, and 23°C measured on a Frazier permeometer. with a kinematic viscosity of 0.37 cm /s was used in this apparatus. The use of glycerine in 2 the solution ensures a laminar flow in the channel. With this solution, the Reynolds number in the channel, i.e., based upon the channel velocity, is less than 2000 for all experimental conditions tested. The viscosity of the solution was measured using a H A A K E Rotary Viscometer. The suction flow rate was measured using a bucket and stopwatch method. The maximum suction flow velocity was 74 mm/s. CHAPTER 6. EXPERIMENTAL WORK • Head Tank R e s e r v o i r Figure 6.1: the schematic of the experimental setup S3 05 3 pa rft Figure 6.2: the assembled test section in (a) exploded view and (b) assembled view -a CHAPTER 6.3 6. EXPERIMENTAL 58 WORK Pulsed Ultrasound Doppler Anemometry Particle Image Velocimetry (PIV) and Laser Doppler Velocimetry (LDV) are the most common methods used to obtain experimental velocity and flow field measurements due to their ability to yield non-invasive real-time results. However, these two techniques are expensive, difficult to set up and not suitable for opaque systems. In comparison, a simple, non-intrusive measurement technique, an acoustic method, was found to be more advantageous for our situation. The commercialized unit we used is called the Pulsed Ultrasound Doppler Anemometer (PUDA) (product name: DOP2000) obtained from Signal Processing SA (Lausanne, Switzerland, http://www.signal-processing.com) and it has been applied previously in a number of studies [119, 120, 121, 122]. 6.3.1 Measuring Principle In pulsed Doppler ultrasound, a short ultrasonic burst is sent to the fluid periodically and a receiver continuously collects echoes issuing from targets (e.g., neutrally buoyant particles) that may be present in the path of the ultrasonic beam (see Fig. 6.3). From knowledge of the time delay td between an emitted burst and the echo issuing from the particle, the depth of the particle can be computed as Pd = ct^/2, where c is the sound velocity of the ultrasonic wave in the liquid. If the particle is moving at an angle 9 in relation to the axis of the ultrasonic beam, its velocity can be measured by computing the variation of its depth between two emissions separated in time by t j: pr {P - Pi) = viscose = ^{t 2 2 - h) (6.1) Since the time difference {t — t\) is always very short, i.e. less than a microsecond, it is 2 easier to measure it with the phase shift of the received echo 5, where t — t\ — 5/(2irf ) and 2 e f is the emitting frequency. With this information the velocity of the target is expressed e by: v = 2/ e = C O S 9t f pr ° fd 2/ e cos 9 where fd corresponds to the frequency shift of the received echo. (6.2) CHAPTER 6. EXPERIMENTAL WORK 59 Transducer Figure 6.3: The working mechanism of pulsed ultrasound Doppler anemometry (image taken from Signal-Processing SA) This technique requires particles present in the fluid to issue echoes. Thus for a fluid without sufficient particles, small particles need to be seeded into the fluid. The ideal particles should be small to minimize gravity effects. As well, since the information from the ultrasound burst is available only periodically, there are limitations of the velocity and the depth that can be measured. The maximum velocity that exists for each pulse repetition frequency (PRF) is 4/ e COS (6.3) Otprf and the maximum depth is P •*- max = (6.4) For this technique, there are a number of parameters that need to be adjusted before a reliable profile is obtained. 6.3.2 Adjustment of Parameters A n example readout is shown in Fig. 6.4. The profile on the right is the echo profile. The two strong echoes at about 4 and 16 millimeters indicate the top wall position and the interface between the porous fabric and the fluid in the channel. The profile on the left, 1 1 Note that the wall thickness adjacent to the ultrasound probe is about 4 millimeters and the velocity of flow in the porous region cannot be measured using P U D A . CHAPTER 6. EXPERIMENTAL WORK 60 on the other hand, shows the velocity profile of the fluid flow in the channel. It is clear that the two strong echoes correspond roughly to the zero velocities. The display information on the left of the velocity profile records the number of the measuring profile and the time this profile was recorded. When the moving average filter is on, the average of a total of 128 profiles is displayed. Typically, before a test, depending on the application, there are many parameters that need to be set up before reliable measurements can be obtained. Knowledge of the flow system helps to make this setup quicker. As stated, the limitations of the PUDA method require that for a shallow channel with relatively low velocities, a high emitting frequency is required. To make the result suited for high flow rate measurement, a low incident angle (or high Doppler angle) of the probe should be used. However too high a Doppler angle results in signal diffusion. Resolution is usually not an issue (0.203 mm) as the number of measuring points (gates) can be chosen to fit the depth of the channel. There are other parameters, e.g., the emitting power, the speed of sound in the fluid and the Doppler angle. The Doppler angle is a measure of the incident angle and the refractive index of the materials of the wall and the fluid. The Doppler angle of our application was determined to be 81° and the sound speed in the glycerine/water mixture is 1621 m/s. About 4 Status Profile time 11.4 187.8 9.G 20 + 15 4- 10 + ms Profile number 9 -Full- •300 Record cto +-100 -200 0 100 200 Exit Emission | Reception ] TGC PRF / profiles ft Processing conditions Compute Profiles to record ; Velocity scale 128 Equalize range _£j Scale factor And skip r i r 1500 300 mm/s H> t E 2000 modulus * Advanced compute | Filters ] Display | Preferences | Trigger | Multiplexer | Cinema j Module Scale <? 1 Sound speed 1S21 _ t j m/s c r r 4 Measure... Number of gates Doppler angle r Use it j ^ Auto — i : No more r j r \ than 9 5 Aut0 ~Z' ^ 8 Figure 6.4: The interface of the pulsed ultrasound Doppler anemometer Dop2000 from Signal Processing SA CT. CHAPTER 6.4 6. EXPERIMENTAL Experimental 62 WORK Procedure Before the start up of any experiment, the reservoir is filled with the water-glycerine solution containing a small amount of fine (~ 60 pm) polymeric seed particles (see Appendix §F.3 for a detailed description). Once the fluid is mixed thoroughly, the pump is started and the fluid is delivered to the head tank, measurements may commence. When the suction rate is high, tiny air bubbles accumulate on the top wall, which lowers the visibility and causes acoustic impedance. This can be readily resolved by suddenly opening the outflow whereby the tiny bubbles are flushed away. The system usually needs about five minutes to stabilize and can be run continuously for two to three hours before the temperature of the fluid (about 15°C) is heated up by two degrees. Once the system is stable, the echo profile should indicate two strong echoes at positions near the four millimeter and fifteen millimeter marks as these correspond to the upper wall position and the position of the interface between the fluid and porous fabrics. The strength of the 2 echo profiles are adjusted by controlling the time gain (in the PUDA software panel) such that there is no saturation at any channel depth and there is still enough echo everywhere in the channel. To avoid signal aliasing, the velocity profiles are monitored to be within the maximum velocity by choosing the right velocity mode and correct pulsed repetition frequency (PRF). The range of the depth of the channel and the velocity of the flow makes the 10 MHz emitting frequency an ideal one; t F PR is typically 64 to 100 /LAS, and the number of gates is less than 100. The 10 MHz ultrasound probe, attached on a specially made stand-alone protractor, is positioned on the top wall at a 30 degree angle of incidence to the top wall surface, where the region of contact between the probe and the wall surface is sealed with ultrasound gel from the surrounding air. From our one-domain numerical predictions, it has been shown that the similarity of the velocity should have been developed when x > 4. The position of the Doppler probe is positioned away from the entrance region, at about x = 8 (here x is dimensionless). Once a stable velocity profile is displayed, the profile is 2 T h e depth of the fluid region measured in the direction of the Doppler angle is slightly higher than the channel height by l . l / s i n ( 8 1 ° ) - 1.1 = 0.014 cm. CHAPTER 6. EXPERIMENTAL 63 WORK ooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooo Figure 6.5: The velocity measured off the probe is the component velocity of particles moving in the ultrasound beam direction. Here v denotes the flow velocity and u denotes p the velocity of flow in the ultrasound beam direction. It shall be noted that 9 denotes the Doppler angle, where the refraction of the sound beam due to the sound speed difference in plexiglass and fluid has been accounted for. recorded in binary and ASCII format files. Note that usually about 10 to 20 profiles are' recorded for a given flow condition. Thus, the average and the standard deviation of the velocity data are computed and subsequently nondimensionlized with the average velocity in the direction of the ultrasound beam in order to allow comparison with the similarity numerical predictions. We chose to compare the measurements to the similarity solution since the algorithm is robust over a wide range of Re. 6.5 Data Evaluation When trying to compare the results from the PUDA measurements with the numerical solutions, some complications occur. In the case when there is no suction, i.e., Re — 0, the flow field is one-dimensional and the x-direction velocity component v is given by v = x v /cos 9, p x where v is the measured velocity in the direction of the ultrasound beam (nonp dimensionlized by the suction velocity U ). This is shown in Fig. 6.5. On the other hand, w CHAPTER 6. EXPERIMENTAL 64 WORK when Re ^ 0, the flow field is two-dimensional and it is difficult to compare accurately the experimental result with the numerical solutions for pure x or y component velocities. A better alternative, is to use the numerical solutions to compute the velocity component in the direction of the ultrasound beam. With this, the depth that the ultrasound beam travels in the channel is apparently Hj sin^. Near the impermeable wall, the streamlines will be nearly parallel to the wall. Near the porous wall, however, the streamlines are curved. In this thesis, the numerical estimates will be projected onto the ultrasound beam direction over which the velocities are measured. The projection is given by v = v p x cos 9 + v sin 9, (6.5) y see Fig. 6.6, where 9 is the Doppler angle, v = (V — x)f and v = f. Upon examination y x y of this projection we see that at y = 0 the measured velocity should be zero as v = v = 0, x y while at y = 1 the velocity should be finite and given by v = v sin 9. Here V is taken p y to be the ratio of the average velocity V* and the mean suction velocity U at a position w where we define x = 0. The numerical velocity v is normalized with its mean velocity in p this direction, which is subsequently compared with the PUDA velocity data. 6.6 3 Results and Discussion A small number of experiments were conducted to verify two aspects of the similarity numerical predictions, namely (i) the transition from stable type-1 to type-II flow with 4 increasing Re and (ii) the existence of multiple steady solutions. In total 8 experiments were conducted and a summary of the experiments and predictions is shown in Table 6.2. The experimental results are shown in Figs. 6.7 - 6.14, the numerical result is shown as 5 solid line and the mean velocity data from PUDA are shown as dotted line. The standard 3 T h e P U D A velocity data are also normalized with its mean velocity. 4 Type-III solution indicates backflow in the porous region. However, since the P U D A method can not be used to measure the flow in the porous medium, the type-Ill flow solutions cannot be verified in this work. 5 I t should be noted that j3 used in each case of the numerical experiments is mostly interpolated or extrapolated from the best j3 fit in Fig. 5.3. For Experiment 4, the Re number is far from the region shown in this figure, therefore we assign (3 = 1.5 for this case. Ochoa-Tapia and Whitaker [81] suggested CHAPTER 6. EXPERIMENTAL 65 WORK / / / / / / / / 0>r J/////////////// J\y Si v ooooooooooooooooooooooooooooooooooooo ooooooooooooooooooooooooooooooooooooo OOOOOOOOOOOOOOOOOQQOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOO000ooooooooooooooo 0 0 OOOOOO POO OOP POO OOP OOP 0 0 OOP OOP POP O J 2 £ ) 1 I I I I 1I 11I I I I Figure 6.6: The schematic of the projections of the x component of the numerical velocity v , y component v on the ultrasound beam direction, the PUDA velocity v . x y p Table 6.2: The experimental conditions tested. X e Prediction 0.84 0.38 type-I 0.84 0.38 type-I 13.1 2.66 x IO" 0.84 0.38 type-I A 23.4 2.66 x IO" 5 B 0 1.90 x I O 6 B 9.5 1.90 x IO" 7 B 13.3 1.90 x IO" 8 B 15.7 1.90 x IO" Experiment Fabric 1 A 0 2.66 x IO" 2 A 6.35 2.66 x IO" 3 A 4 Re Da 3 3 3 3 0.84 0.38 type-II -3 0.92 0.41 type-I 0.92 0.41 type-I 3 0.92 0.41 type-I 3 0.92 0.41 type-II, multiple steady states 3 deviations of the experimental observations are also shown in these figures. Clearly, in general there is good agreement between the numerical predictions based upon the similarity method and the experimental measurements. In particular, excellent agreement was P should be of order 1 and they found — 1 < /? < 1.47 in a set of experiments conducted by Beavers and Joseph [87, 88]. CHAPTER 6. EXPERIMENTAL WORK 66 found when there is no suction or when suction Re is small (i.e., type-I solutions). Both results predict a shift of the maximum velocity from the mid-plane of the channel towards the porous medium as the suction Re increases. Included in thesefiguresare the predicted velocity profiles in the porous medium. Our experimental method is incapable of measuring these profiles. We included this as it is still interesting to see that increased suction changes the curvature of the velocity profiles near the interface. In one case, Experiment 4, we confirmed our predictions regarding backflow in the channel. Clearly the fluctuations in the flow, as given by the error bars in this figure, were found to be larger for type-II than type-I flows. Similarly in Experiment 8 we observed one stable solution and then after about 30 seconds the flow changed to another profile. Under this suction condition, the velocity profile underwent slow oscillations between two different backflow profiles, one with a small backflow and the other with a stronger backflow. Qualitatively, it seems that the two bifurcated numerical solutions agree with the PUDA results. In general the agreement is not as good with type-II flow solutions as with type-I solutions. Experimentally, stable backflows can be observed only when the fluid flow in the apparatus system is started up gently and kept stable under a combination of low inflow and low suction conditions, i.e., barV is small. With this, the disturbances from the surrounding environment, e.g., pump vibrations, upstream flow disturbances, etc., are minimized. Even when backflows are found, they seem to be very sensitive to disturbances. When the suction flow is high, the risk of introducing more disturbance is possible. This is evident in Experiment 4, for which the standard deviations are much higher than in any other cases. The relatively "unstable" backflow-type solutions can be understood from the theory of the instability of the boundary layer, e.g., the flat plate flow. In essence, backflow indicates that there is a point of inflection of the velocity profile in the flow, with an excessive adverse gradient, and the main flow separates from the wall. The backflow-type profile is vulnerable to transitions to instability or turbulence. Furthermore, it seems that the flow is more or less time-dependent; thus a numerical solution of the time-dependent flow might better approximate the experimental results. In this case, the solution is probably of a limit cycle form, i.e., a spatially non-uniform solution that is periodic in time, similar to that suggested by King and Cox [4]. CHAPTER 6. EXPERIMENTAL WORK 67 Regardless of the agreement, it should be noted that there are potential errors in the determination of Re values, since the channel height, the viscosity and the density of the fluid, as well as the mean suction velocity, all have inherent error. An error analysis showed that the accuracy of Re based on the uncertainty of viscosity, channel height, suction velocity and density of the fluid is approximately ±5 — 11% of Re. In addition to this, it is experimentally difficult to measure V at x = 0, and it is hard to estimate how uniform the suction flow U is, and how it affects the similarity of the flow. Moreover, w there are local variations of the channel heights h and H. We may have overestimated the porosity e with our measuring method, since at the edges of the fabrics, the pore openings are larger. As well, the permeability of the porous fabric was obtained from the product supplier, and may not be accurate. Measuring a Darcy permeability is possible, but may not fit our application, since the flow that pastes through the thin porous medium may well exceed a pore-scale Reynolds number of order 1 [77]. Finally, the wall effects in the PUDA measurement have sometimes made measuring flow velocity near the wall difficult, although by varying the angle of the probe with respect to the wall plane, it is possible to minimize this effect. Furthermore, because of the wall effects, and the number of the gates (or resolution) in measuring the velocity profile, the position of the wall or the surface of the porous fabric determined with the echo profiles generally correspond to the channel depth /i/sin#, but uncertainties exist and we estimated them to be within ±1.5 mm. Nevertheless, the measurements indicate the instability caused by the onset of the reverse-type flows. CHAPTER 6. EXPERIMENTAL WORK 68 Figure 6.7: The similarity numerical solutions shown as solid line (here (3 = 0.13) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 1 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 69 Figure 6.8: The similarity numerical solutions shown as solid line (here P = 0.3) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 2 in Table 6.2. CHAPTER 6. EXPERIMENTAL 0 c o r- 0.2 • 70 WORK 1 3 1 0.4 —i CL 0 § 0.6 - sz CO O 0.8 _i—• 0 0 5 Velocity 1 5 Figure 6.9: The similarity numerical solutions shown as solid line (here (5 = 0.9) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (clotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 3 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 0 71 VelocitV Figure 6.10: The similarity numerical solutions shown as solid line (here P = 1.5) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 4 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 72 Figure 6.11: The similarity numerical solutions shown as solid line (here (3 — 0.15) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 5 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 73 Figure 6.12: The similarity numerical solutions shown as solid line (here (3 = 0.58) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 6 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 74 Figure 6.13: The similarity numerical solutions shown as solid line (here ft = 1.1) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for Experiment No. 7 in Table 6.2. CHAPTER 6. EXPERIMENTAL WORK 75 Figure 6.14: The similarity numerical solutions shown as solid line (here (3 = 1.1) for the velocity profile in the direction of the ultrasound beam compared with the experimental results from PUDA (dotted line). The error bars denote one standard deviation from the mean velocity values. The simulation is for the two steady state solutions for Experiment No. 8 in Table 6.2. Chapter 7 S U M M A R Y A N D CONCLUSIONS Twin wire forming is a papermaking process in which a fibre suspension is drained between two permeable wires to form a continuous paper web. This process (somewhat) resembles cross-flow filtration. The purpose of this study is to estimate the shear at this interface. In this thesis we focused our efforts to understand the coupled two-dimensional flow of a Newtonian fluid both above and through a porous medium. In the fluid-only region, the two-dimensional flow field was governed by the Navier-Stokes equation. We consider the Brinkman-extended Darcy law relationship in the porous medium. Inertial terms are retained in the formulation and the interface conditions between the two domains are those given by Ochoa-Tapia and Whitaker [81]. It should be noted that these interface conditions are formulated with an empirical constant (3 that is unknown a priori. The model equations were solved using two independent methods. In the first method we pose a similarity variable and reduce the governing equations to two, coupled, non-linear ordinary differential equations. Depending on the suction Reynolds and Darcy numbers, three types of solutions were found, i.e. no backflow, backflow near the solid wall, and backflows near the solid wall and in the porous medium. Typically, backflow-type solutions exist at large Re. In addition, the solution bifurcates and multiple steady states were found. The characteristics of the non-unique solutions were studied via a linear stability analysis, i.e., the method of normal modes. In the next section of the work, the governing equations were re-posed as a one-domain problem, using the procedure outlined by Basu and Khalili [100] and solved in primitive 76 CHAPTER 7. SUMMARY AND CONCLUSIONS 77 variable form, using a finite volume formulation. Finally, experiments were conducted to measure the velocity profiles and generally good agreement with the numerical predictions was found. Chapter 8 RECOMMENDATIONS The numerical methods presented in Chapter 4 are limited to Re < 170. This limits the usefulness of the solution. Attempts were made to perform a singular perturbation expansion to analyze the case when Re —> +/ — oo. However, this problem proved to be too difficult and we were unsuccessful. In addition, we explored an extension of the similarity transform for the time dependent case. This posed a problem for the numerical algorithm used. It would be interesting to explore a spectral method to tackle this problem. Finally, to make this thesis more relevant to the papermaking industry, a particulate phase must be included. At this point we feel that it may be impossible to use the similarity ansatz to solve the flow field when the fluid is a fibre suspension. However, this method may be useful to understand the motion of single particles or even the particle deposition mechanism during papermaking. 78 Bibliography [1] S. Zahrai. On the Fluid Mechanics of Twin-Wire Formers. PhD thesis, Royal Institute of Technology, S-100 44 Stockholm, Sweden, October 1997. [2] C.P.J. Bennington, R.J. Kerekes, and J.R. Grace. The yield stress of pulp suspension. Can. J. Chem. Eng., 68:748-757, 1990. [3] A.S. Berman. Laminar flow in channels with porous walls. J. Appl. Phys., 24:12321235, 1953. [4] J.R. King and S.M. Cox. Asymptotic analysis of the steady-state and time-dependent Berman problem. J. Engng. Maths., 39:87-130, 2001. [5] B. Norman. On the mechanisms of dewatering in the twin-wire and press sections. Nordic Pulp Pap. Res. J., 6:39-46, 1987. [6] H. Karema, M . Kataja, M . Kellomaki, J. Salmela, and P. Selenius. Transient fluidization of fiber suspensions in straight channel flow. Tappi International Paper Physics Conference, San Diego, Ca, USA:369, 1998. [7] J. Gullichsen and E Harkonen. Medium consistency technology I. fundamental data. Tappi J., 64(6):69-72, 1981. [8] C.P.J. Bennington and R.J. Kerekes. Power requirements for pulp suspension fluidization. Tappi J., 79(2):253-258, 1996. 79 80 BIBLIOGRAPHY [9] J. Hietaniema and J. Gullichsen. Floe disruption in medium consistency fiber suspensions: Empirical results. International Paper Physics Niagra-on-the- Conference, Lake, Ontario, Canada:29-38, 1995. [10] A. Kuznetsov, A. Bagaev, and A. Dolgikh. Use of surfactants to improve the rheological properties of high consistency fiberboard pulp (engl. trans). Promst., Derevoobrab. 5:2-4, 1995. [11] C.T. Scott, S. Zauscher, and D.J. Klingenberg. Rheology and extrusion of low-grade paper and sludge. Tappi International Environmental Conference, Book 2:685-690, 1999. [12] 0. Jokinen and K. Ebeling. The flocculation tendency of papermaking fibres. Paperi ja Puu, 67(5):317-325, 1985. [13] C.T. Scott and S. Zauscher. Pulp extrusion at ultra-high consistencies. Tappi International Environmental Conference, Book 2:739-743, 1997. [14] A. Swerin, R.L. Powell, and L. Odberg. Linear and nonlinear dynamic viscoelasticity of pulp fibre suspensions. Nordic Pulp Pap. Res. J., 7(3):126-132, 1992. [15] A. Swerin. Flocculation Flocculated and Fibre Network by a Retention-Aid System. Strength in Papermaking Suspensions PhD thesis, Royal Institute of Technology, Stockholm, Sweden, 1995. [16] T. Wikstrom and A. Rasmuson. Yield stress of pulp suspensions. The influence of fibre properties and processing conditions. Nordic Pulp Pap. Res. J., 13(3):243-249, 1998. [17] F.E. Farrington. More fundamental approach to the problem of high consistency forming. In Tappi Engineering Conference, volume 2, pages 709-717, The Westin Atlanta, GA, 1986. Tappi Press. [18] K.-J. Grundstom, B.G. Norman, and D. Wahren. High consistency forming of paper. Tappi J., 56(7):81-84, 1973. 81 BIBLIOGRAPHY [19] T. Nomura, K. Wada, and T. Shimuzu. High consistency forming. Part 1: Research and development of headboxes. Tappi J., 72(1):115—122, 1989. [20] T. Nomura, K. Wada, and T. Shimuzu. High consistency forming. Part 2: Pilot plant tests. Tappi J., 72(4):171-176, 1989. [21] T. Nomura, K. Wada, and T. Shimuzu. High consistency forming. Part 3: Sheet quality and engineering data. Tappi J., 72(5):187-192, 1989. [22] R.E. Hergert and C L . Sandford. Pressure measurements in the forming zone of a twin-wire tissue machine. Proc. Int. Water Removal Symp., Vancouver:47-50, 1982. [23] R.H. Zhao and R.J. Kerekes. Pressure distribution between forming fabrics in blade gap formers. J. Pulp Pap. Sci., 21(3):97-103, 1995. [24] W.D. Baines. Flow in the formation zone of a twin-wire machine. Pulp Pap. Can., 68(10):T497-505, 1967. [25] H. Meyer. Hydrodynamics of the sheet forming process. Tappi J., 54(9):432-436, 1971. [26] J. Koskimies, J. Perkinen, H. Puolakka, E. Schultz, and B. Wahlstrom. A drainage model for the forming zone of a twin-wire former. Papper och Tra, 4:137-146, 1972. [27] D. Wahren, D. Dufva, and B. Wahlstrom. Mechanics of water removal in webstertype formers. Paper Technology and Industry, 16(2), 1975. [28] P.F. Turnbull, N.C. Perkins, W W . Schultz, and P.D. Beuther. One dimensional dynamic model of a paper forming process. Tappi J., 80(l):245-253, 1997. [29] H. Iwata, M. Hasuike, T. Adachi, T. Bando, and K. Sakamoto. Development of mitsubishi new former for papermaking machine. Mitsubishi Heavy Industries Technical Review, 30:1-6, 1993. [30] A.J. Panshin and C. de Zeeuw. Textbook of Wood Technology. edition, 1980. McGraw Hill, 4 h l 82 BIBLIOGRAPHY [31] R.J. Kerekes, R.M. Soszynski, and P.A. Tarn Doo. The flocculation of pulp fibres. Fund. Res. Symp., Oxford(l):265-310, 1985. [32] R.L. Powell, M . Weldon, S. Ramaswamy, and M.J. McCarthy. Characterization of pulp suspensions. Tappi Engineering Conference, pages 525-534, 1996. [33] S.G. Mason. Fibre motions and flocculation. Pulp Pap. Mag. Can., 55(13):96-102, 1954. [34] R. Meyer and D. Wahren. On the elastic properties of three-dimensional networks. Svensk Papperstid., 67(10):432-436, 1964. [35] R.M. Soszynski and R.J. Kerekes. Elastic interlocking of nylon fibres suspended in liquid. Part 2. process of interlocking. Nordic Pulp Pap. Res. J., 3(4): 180-184, 1988. [36] R.M. Soszynski. The apparent volumetric concentration of wood pulp fibres in suspension. Appita J., 42(5):362-363, 1989. [37] R.J. Kerekes and C.J Schell. Effects of fibre length and coarseness on pulp flocculation. Tappi J., 78(2): 133-139, 1995. [38] R.M. Soszynski. Fibre Flocculation in Shear Flow. PhD thesis, University of British Columbia, Vancouver, BC, Canada, 1988. [39] C.T.J. Dodson. Fibre crowding, fibre contacts, and fibre flocculation. Tappi J., 79(9):211-216, 1996. [40] D.M. Martinez, H. Kiiskinen, A-K. Ahlman, and R.J. Kerekes. On the mobility of flowing papermaking suspensions and its relationship to formation. J. Pulp Pap. Sci., 29(10), 2003. [41] M.B. Mackaplow and E.S.G Shaqfeh. A numerical study of the sedimentation of fibre suspensions. J. Fluid Mech., 376:149-182, 1998. [42] R.F Ross and D.J. Klingenberg. Simulation of flowing wood fibre suspensions. J. Pulp Pap. Sci., 24(12):388-392, 1998. 83 BIBLIOGRAPHY [43] M . Doi and D. Chen. Simulation of aggregating colloids in shear flow. J. Chem. Phys., 90:5271-5279, 1989. [44] M.A. Turney, M.K. Cheung, M.J. McCarthy, and Powell R.L. Hindered settling of rod-like particles measured with magnetic resonance imaging. AIChE J., 41:251-257, 1995. [45] P. Kumar and B.V. Ramarao. Enhancement of the sedimentation rates of fibrous suspensions. Chem. Engng. Comm., 108:381-401, 1991. [46] B. Herzhaft, E. Guazzelli, M.B. Mackaplow, and E.S.G Shaqfeh. An experimental investigation of the sedimentation of a dilute fibre suspension. Phys. Rew. Lett., 77(2):290-293, 1996. [47] R. Davis. The experimental study of the differential settling of particles in suspension at high concentration. Powder Technol, 2:43-51, 1968. [48] I.L. Claeys and J.F. Brady. Suspensions of prolate spheroids in stokes flow, part 3. hydrodynamic transport properties of crystalline dispersions. J. Fluid Mech., 251:479-500, 1993. [49] D.L Koch and E.S.G Shaqfeh. The instability of a dispersion of sedimenting spheroids. J. Fluid Mech, 209:521-542, 1989. [50] W.L. Ingmanson, B.D. Andrews, and R.C. Johnson. Internal pressure distributions in compressible mats under fluid stress. Tappi J., 42(10):840-849, 1959. [51] H. Meyer. A filtration theory for compressible fibrous beds formed from dilute suspensions. Tappi J., 45(4):296-310, 1962. [52] N. Pan. A modified analysis of the microstructural characteristics of general fiber assemblies. Textile Res. J., 63(6):336-345, 1993. [53] J. Ringner. The influence of fiber length distribution on the network strength of fiber suspensions. Master's thesis, Chalmers Univ., Goteborg, 1995. 84 BIBLIOGRAPHY [54] G. Steenberg, N. Thaln, and D. Wahren. Formation and properties of fiber networks. Consolidation of the Paper Web, Trans, of the Symposium, Cambridge, 1965. [55] N. Thalen and D Wahren. Shear modulus and ultimate shear strength of some paper pulp fiber networks. Svensk Papperstidn., 67(7):259-264, 1964. [56] D.M. Martinez, K. Buckley, S. Jivan, A. Lindstrom, R. Thurigaswamy, T.J. Ruth, and R.J. Kerekes. Characterizing the mobility of fibres during sedimentation. In The Science of Papermaking, Transactions of the 12th Fundamental Research Symposium, volume 1, pages 225-254, Oxford, September 2001. [57] F.M. Auzerais, R. Jackson, W.B. Russel, and W.F Murphy. The transient settling of stable and flocculated suspensions. J. Fluid Mech., 221:613-639, 1990. [58] J. Ystrom. On the Numerical Modeling of Concentrated Suspensions and of Viscoelastic Fluids. PhD thesis, Royal Institute of Technology (KTH), Stockholm, Sweden, 1996. [59] H. Enwald, E. Peirano, and A.E. Almstedt. Eulerian two-phase flow theory applied to fiuidization. Int. J. Multiphase Flow, 22:21-66, 1996. [60] I. Tosun, M.S. Willis, F. Desai, and G.G. Chase. Analysis of drag and particulate stress in porous media flows. Chem. Eng. Sci., 50(12): 1961—1969, 1995. [61] A.A. Robertson and S.G. Mason. Specific surface of cellulose fibres by the liquid permeability method. Pulp Pap. Mag. Can., 12:103-110, 1949. [62] P. Nilsson and K.O. Larsson. Paper web performance in a press nip. Pulp Pap. Mag. Can., 12:T438-T445, 1968. [63] S.T. Han. Compressibility and permeability of fibre mats. Pulp Pap. Mag. Can., 5:T134-T146, 1969. [64] E.R. Ellis. Compressibility Pulp and Its Application and Permeability to the Prediction of Never Dried Bleached Softwood Kraft of Wet Pressing University of Maine, Orono, Me, USA, 1981. Behaviour. PhD thesis, BIBLIOGRAPHY [65] G. Carlsson and Lindstrom. Permeability to water of compressed pulp fibre mats. Svensk Paper stidnmg, 86(12):128-134, 1983. [66] G. Ljungkvist. Pulp Characterization by Permeability Measurement. PhD thesis, Chalmers University of Technology, Gothenburg, Sweden, 1983. [67] H. Vomhoff. Dynamic Compressibility Local Stress Variations of Water Saturated Networks in Wet Pressing. and Influence of PhD thesis, Royal Institute of Technology, Stockholm, Sweden, 1998. [68] P.E. Wrist. The present state of our knowledge on the fundamentals of wet pressing. Pulp Pap. Mag. Can., 65(7):T284-T296, 1964. [69] H.D. Wilder. The compression and creep properties of wet pulp mats. Tappi J., 43(8):715-720, 1960. [70] W.L. Ingmanson and R.P. Whitney. The filtration resistance of pulp slurries. Tappi J., 37(ll):523-534, 1954. [71] K. A-S. Jonsson and B.T.L Jonsson. Fluid flow in compressible porous media I: Steady state conditions. AIChE J., 38(9): 1340-1348, 1992. [72] D.M. Martinez. Characterizing the dewatering rate in roll gap formers. J Pulp Pap Sci, 24(1):7-13, 1998. [73] N. Moch. On the Variable Permeability chines, During Dewatering Process in Paper Ma- undergraduate thesis, Department of mechanics, Royal Institute of Technol- ogy (KTH), Stockholm, Sweden, 1995. [74] S.I. Green and R.J. Kerekes. Numerical analysis of pressure pulses induced by blades in gap formers. Tappi J., 81(4), 1998. [75] S. Zahrai, D.M. Martinez, and A.A Dahlkild. Estimating the thickness of the web during twin-wire forming. J. Pulp Pap. Sci., 24(2):67-72, 1998. 86 BIBLIOGRAPHY [76] M. Jansson. Fiberriktningsanisotropi-variationer i z-led. undergraduate thesis, Stock- holm, Sweden, Department of Paper Technology, Royal Institute of Technology (KTH), 1999. [77] D.A. Nield and A. Bejan. Convection in Porous Media. Springer, 2 nd edition, 1995. [78] H.P.G. Darcy. Les Fontaines Publiques de la Ville de Dijon. Victor Dalmont, Paris, 1856. [79] H.C. Brinkman. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl. Sci. Res., A(l):27-34, 1947. [80] R.A. Wooding. Steady state free convection of liquid in a saturated porous medium. J. Fluid Mech., 2:273-285, 1957. [81] J. Ochoa-Tapia and S. Whitaker. Momentum transfer at the boundary between a porous medium and a homogeneous fluid I: Theoretical development. Int. J. Heat Mass Transfer, 38:2635-2646, 1995. [82] J. Ochoa-Tapia and S. Whitaker. Momentum transfer at the boundary between a porous medium and a homogeneous fluid II: Comparison with experiment. Int. J. Heat Mass Transfer, 38:2647-2655, 1995. [83] J.C. Ward. Turbulent flow in porous media. ASCE J. Hydraul. Dw., 90:1-12, 1964. [84] K. Vafai and C. Tien. Boundary and inertia effects on flow and heat transfer in porous media. Int. J. Heat Mass Transfer, 24:195-203, 1981. [85] J.L. Lage. Effect of the convective inertia term on benard convection in a porous medium. Numer. Heat Transfer A, 22:469-485, 1992. [86] D.M. Manole and J.L. Lage. The inertial effect on the natural convection flow within a fluid saturated porous medium. Int. J. Heat Fluid Flow, 14:376-384, 1993. [87] G.S. Beavers and D.D. Joseph. Boundary conditions at a naturally permeable wall. J. Fluid Mech, 30:197-207, 1967. 87 BIBLIOGRAPHY [88] G. Beavers, E. Sparrow, and R. Magnuson. Experiments on coupled parallel flows in a channel and a bounding porous medium. J. Basic Engineering, Series D, Trans. ASME, 92:843-848, 1970. [89] G. Neale and W. Nader. Practical significance of Brinkman's extension of Darcy's law: Coupled parallel flows within a channel and a bounding porous medium. Can. J. Chem. Eng., 52:475-478, 1974. [90] G.I. Taylor. A model for the boundary condition of porous material: Part 1. J. Fluid Mech., 49:319-326, 1971. [91] S. Richardson. A model for the boundary condition of porous material: Part 2. J. Fluid Mech., 49:327-336, 1971. [92] P.G. Saffman. On the boundary condition at the surface of a porous medium. Studies in Appl. Maths., 50:93-101, 1971. [93] K. Vafai and R. Thiyagaraja. Analysis of flow and heat transfer at the interface region of a porous medium. Int. J. Heat Mass Transfer, 30:1391-1405, 1987. [94] K. Vafai and S.J. Kim. Fluid mechanics of the interface region between a porous medium and a fluid layer - An exact solution. Int. J. Heat Fluid Flow, 11:254-256, 1990. [95] A.V. Kuznetsov. Influence-of the stress jump condition at the porous-medium/clearfluid interface on a flow at a porous wall. Int. Comm. Heat Mass Transfer, 24:401-410, 1997. [96] A.V. Kuznetsov. Analytical investigation of Couette flow in a composite channel partially filled with a porous medium and partially with a clear fluid. Int. J. Heat Mass Transfer, 41:2556-2560, 1998. [97] A.V. Kuznetsov. Fluid mechanics and heat transfer in the interface region between a porous medium and a fluid layer: A boundary layer solution. J. Porous 2(3):309-321, 1999. Media, 88 BIBLIOGRAPHY [98] B. Alazami and K. Vafai. Analysis of fluid flow and heat transfer interfacial conditions between a porous medium and a fluid layer. Int. J. Heat Mass Transfer, 44:17351749, 2001. [99] J. A. Ochoa-Tapia and S. Whitaker. Momentum jump condition at the boundary between a porous medium and a homogeneous fluid: Inertial effects. Porous Media, Journal of 1(3):201-217, 1998. [100] Basu A. J. and A. Khalili. Computation of flow through a fluid-sediment interface in a benthic chamber. Phys. Fluids, 11(6):1395-1405, 1999. [101] R.D. Bundy and H.L. Weissberg. Experimental study of fully developed laminar flow in a porous pipe with wall injection. Phys. Fluids, 13:2613-2615, 1970. [102] P. L. Donoughe. Analysis of laminar incompressible flow in a semiporous channels. National Advisory Committee for Aeronautics Technical Notes 3759, 1956. [103] E. Eckert, P. Donoughe, and B. Moore. Velocity and friction characteristics of laminar viscous boundary layer and channel flow over surfaces with ejection or suction. National Advisory Committee for Aeronautics, Technical Note 4102, 1957. [104] R.M. Terrill. Laminar flow in a unifomly porous channel. Aeronaut. Q., 15:299-310, 1964. [105] F.M. Skalak and C.-Y Wang. On the non-unique solutions of laminar flow through a porous tube or channel. SIAM J. Appl. Maths, 34:535-544, 1978. [106] L. Durlofsky and J.F. Brady. The spatial stability of a class of similarity solutions. Phys. Fluids, 27:1068-1076, 1984. [107] M.B. Zaturska, P.G. Drazin, and W.H.H. Banks. On the flow of a viscous fluid driven along a channel by suction at porous walls. Fluid Dyn. Res., 4:151-178, 1988. [108] S. Ferro and G. Gnavi. Spatial stability of similarity solutions for viscous flows in channels with porous walls. Phys. Fluids, 12(4):797-802, 2000. 89 BIBLIOGRAPHY [109] S.M. Cox. Two dimensional flow of a viscous fluid in a channel with porous walls. J. Fluid Mech, 227:1-33, 1991. [110] S.M. Cox. Analysis of steady flow in a channel with one porous wall, or with accelerating walls. SI AM J. Appl. Math, 51:429-438, 1991. [Ill] S.M. Cox and A.C. King. On the asymptotic solution of a high order non-linear ordinary differential equation. Proc. R. Soc. London A, 453:711-728, 1997. [112] J. F. Brady and A. Acrivos. Closed-cavity laminar flows at moderate reynolds numbers. J. Fluid Mech., 115:427-442, 1982. [113] G. Lindfield and J. Penny. Numerical Methods Using Matlab. Prentice Hall, 2 edition, 2000. [114] J. Cash and A. Karp. A variable order runge-kutta method for initial-value problems with rapidly varying right-hand sides. A CM Transactions on Mathematical Software, 16:201-222, 1990. [115] R. P. Brent. Algorithms for minimization without derivatives. Prentice-Hall, Engle- wood Cliffs, New Jersey,, 1973. [116] E.W.K. Young. The permeability and compressibility of semi-dilute pulp fibre suspensions: inversely solving the governing pde of sedimentation. Master's thesis, the University of British Columbia, 2003. [117] S. Zahrai and F. Bark. On the fluid mechanics of twin-wire blade forming in paper machines. Nordic Pulp Pap. Res. J., 4(10):245-252, 1995. [118] C. Deng and D.M. Martinez. Viscous flow in a channel partially filled with a porous medium and with wall suction, in press. Chem. Eng. Sci., 2004. [119] Y . Takeda, W.E. Fisher, J. Sakakibara, and K. Ohmura. Experimental observation of the quasiperiodic modes in a rotating couette system. Physical Review E, 47:41304134, 1993. BIBLIOGRAPHY 90 [120] P. Petitjeans, J.E. Wesfreid, and J.-Cl Attiach. Vortex stretching in a laminar boundary layer flow. Exp. fluids, 22:351-353, 1997. [121] D. Brito, H. Nataf, P. Cardin, J. Aubert, and J. Masson. Ultrasonic doppler velocimetry in liquid gallium. Exp. Fluids, 31:653-663, 2001. [122] S. Eckert and G. Gerbeth. Velocity measurements in liquid sodium by means of ultrasound doppler velocimetry. Exp. Fluids, 32:542-546, 2002. Appendix A D E R I V A T I O N OF T H E P R E S S U R E B O U N D A R Y CONDITION A.l N o r m a l Stress As described in the fluid mechanics textbooks, a dimensional Newtonian normal stress at the interface between the fluid r* and porous regions r£ are respectively r*_ = -P* dv* 2p-^ dy* du* + v ; = -p; + 2 ^ T (A.l) where v*, v*, u* and u* denote the y component velocities in fluid and porous regions respectively. When r* and P* are scaled with pU ,.we obtained the dimensionless conditions 2 as following: r + = -P+ + 2 du,, Re dy From Eq. (4.4), the above equations reduce to r.=-P- + | / . 91 (A.3) APPENDIX A. DERIVATION OF THE PRESSURE BOUNDARY CONDITION 92 In the following we will derive Eq. (4.9) from the normal stress matching condition r = r_. + To begin with, we observe that since g {x,t) = f {x,t) y (from Eq. (4.8)), Eq. r = r_ y + reduces to P = P_. + A . 1.1 Fluid Region The equations of motion for each velocity component of fluid region yield OP- dv dx dt dv x = _ OP- dt x fyy Re = By definition ~ dP f d P - Re dx dy 1 yy dy 2 +^Vy f) .O 2 Re dx K (A.4) Vy ' dy 2 TTT + v 7T7 + 757 (~~xl7 + 7£T7") v -ft + k-ffv J = J-dx= y rs 2 yyy V 2 T j-J dVy -fl " dx ^ V x d v, 2 dy y dVy+ yt 1 ,d v . x dx x (V-x)(-f dVy dy dv x dX f + ("/vt + ^ d P - J (A-5) J-dV dy - ffyy + ft) J(V- X)dx + J( fe- ffy ~ f ft)d V (A.6) Upon integration of dP, P- = i'fyt +^ - / 7 + fi)(Vx - ^x ) + A _ 2 W lp _ | / t d y + C„ (A.7) where Co is the constant of integration. Assume the pressure at x = 0 and interface y — x is Po, we obtained that Po = {-fyt + § f - ffyy + f )0 + ^ -\f -JMx)dy 2 2 y e + Co (A.8) where / is evaluated at interface y = x, thus Co = Po - ^ + \f {x) 2 + j ft(x)dy (A.9) Hence, P- - Po = {-fyt + ^ - ffyy + f )(Vx 2 -^(/ 2 - f {x)) 2 - \x ) + ^ 2 - [{ft- ft{x))dy (A.10) APPENDIX A. DERIVATION A.1.2 OF THE PRESSURE BOUNDARY CONDITION 93 Porous Region Similarly, for porous region, dP+ _ dx ~ , + dy __^x_}_ du^ ( dt = ( ^_ _ \ e dP = e dx 2[Ux l ) ( _I f a + ,9mL^ik eRe a du 1 x + U y dy ) ^ _ ^ dx 2 S_^. 1 2 x [ + du 2 eRe + £ du x dy + 2 DaR } Ux e (A.n) ) £_ DaRe e 2 (A 12) ' V ; Integration of dp gives + Hence, , + & -gy(x) g -9 (X) 2 ~7Pe _____ 2 ^ - D a - R e j 2 ( \ i zi \ H,*, *^ 1 { g x + l 9 t ) d y ( ' A M ) Equate Eq. (A. 10) and Eq. (A. 14) at y = x, i.e., P+ = P_ yields + Pe "™ + - ^Re~^~ + + ^~Da~Re ( A ' 1 5 ) Equation (A.15) is identical to Eq. (4.9). It is to be noted that this condition is in effect equivalent to the condition of pressure gradient matching in x direction. A.2 Pressure D r o p From the expressions derived in Eqs. (A. 10) and (A. 14), the steady case pressure drop in x and y directions, AP at the interface y — x -d AP at a position x are respectively: an X y • for fluid region AP, = P( y)-P(0,y) = A.{Vx-j) X) = (-Lj -ff yyy yy + f )(Vx-j) 2 (A.16) APPENDIX A. DERIVATION OF THE PRESSURE BOUNDARY CONDITION 94 and AP y = P(x, y) - P(x, 0) = ±-f y - ^ (A.17) • for porous region AP, 1 p(r\ n\ P(x, n,\ y) — - P(0, y) — = C{^-g ~ Ree = yyy yyyy ^2 e ,T-> X 2 1 1 o - j^BZ9y DaRe 1 + ' e>»» J*™"'" mi w.-, £ . 2 ~ 2 (A.18) i V x - j ) and A P = P ,y) B {x - P(x,0) = - ^ £ 9 - | j (A.19) where Ai and A are the previous-mentioned integration constants of the steady equations 2 of Eqs. (4.5) and (4.6) in the numerical procedure. It could be easily identified that Ai = A /e. It has to be noted that the pressure drop AP computed from the similarity 2 X solution is added a reference pressure in order to compare to the pressure drop from the one-domain solution. Similarly a reference pressure is added to the pressure drop AP of y the similarity solutions at position x. Appendix B T H E A S Y M P T O T I C SOLUTIONS The leading order and first order asymptotic approximation functions for f , j\ and g , 0 0 g_ with Da = 1 x lCT , £ = 0.8, e = 0.9 and (5 = 0.2 are: 3 fo(y) = -3.5171y + 4.3602y /i(y) = -0.1767y + 0.5112y - 0.6337 y + 0.5535y - 0.2700y g (y) = 0.9796 + 0.0211 y - 6.7389 x 10- e ^ - 1.7872 x 10 e- ^ (y) = 1.5513 x lO" e ' - 0 gi 3 (B.l) 2 7 6 5 17 7 30 J 3 30 8 - 5.8370 x 10~ - » + 30 32 3.6291 x 10" - - 3.4234 x 10" ye ° + 6.1728 x 10~ y - 9.1245 x 10" y e ^ - 9.8032 x 1 0 - e - + 8.1680 x W- e 3 0 ? y + 6 e 30 17 8 18 30y 30 (B.3) 9 0 e + 19 3 0 (B.2) 2 ° + 8.4463 x 10" 6 3 0 y + 3 e - 3.6965 x 10" 30y+3 2 6 30y+30 32 y 0 ? ; + 3 2 ° 3 0 y + 3 0 e ° y - 0.003321 y + 3.7519 ' l 0 ~ 2 x + 0.003203 3 0 e 3 0 y + 3 0 (B.4) The solutions are compared with the numerical solutions as shown in Fig. B . l , good agreement was found in —4.5 < Re < 4.5. 95 APPENDIX B. THE ASYMPTOTIC SOLUTIONS 96 Fi gure B . l : The shear stress function —f (x) t the interface for different Re, it compares a yy the similarity numerical solution (o) with the asymptotic solution (-). The simulation is for the case Da = 1 x I O , x = 0.8, (5 = 0.2 and e = 0.9. -3 Appendix C CODE VALIDATION We attempt to validate our numerical routine by examining conditions which should approximate the classical theory given by [3]. To do this, we set Da = 1000 and p = 0. When Da is large, the sink term in Eq. (4.6) vanishes and the problem is reduced to Berman's classical treatment. The validation is shown in Fig. C l , good agreement was achieved. The discrepancy is higher when Re increases, in this case the sink term becomes increasingly important. 97 APPENDIX C. CODE 98 VALIDATION 140r 120- Figure C . l : The shear stress function —g (l) yy at the porous wall for different Re, it compares the similarity numerical solution (•) for Da = 1000, with Berman solution (-) [3, 4]. = X 0-8, ft = 0 and e = 1 Appendix D C O D E FOR T H E SIMILARITY SOLUTION The code written in Matlab 6.5 for the similarity solution is included here. Two methods are presented. The first one is based on a regular shooting method in conjunction with Newton's method and the second is based on our modification of Terill's method [104], which we called "the reverse solver". D.l Shooting Method The main function files, are listed in the following. • similarity_shooting.m • eqfluid.m • eqporous m : • odesolve.m • rkck.m • gaussj2.m 99 APPENDIX D. CODE FOR THE SIMILARITY 0/ 0/ 0/ 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / SOLUTION 0/0/ 0/ 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0/0/ o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / ID h o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o /o / °/ Similarity °/ Last 0 0 Solution - updated - o / o / o / o / o / o / o / o / o / h/o/o/o/o Shooting Method 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / /a /.la h o / o / o / o / o / o / o / o / o /o / o / °/ x 0 May 2 9 2 0 0 4 % 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / o/ o/ o/ o/ o/ o/ o/ o/ o/ o/o/ o/ o/ % 10 lo lo lo lo lo lo lo lo lo lo lo lo lo to lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo lo 'Similarity_shooting.nl' is 7, c a l c u l a t e the similarity the main s c r i p t solutions 7o P a r a m e t e r s ( a l l dimensionless) : 7, Re: suction Reynolds 7, D a : D a r c y number using file to shooting of the porous medium the porosity of the fluid, by d e f a u l t 7o p o r o s i t y m : the porosity of the porous medium, beta: 7o height: e m p i r i c a l constant the 7o I n t e r m e d i a t e height the functions and Newton's to method. number 7o p o r o s i t y f : 7, call of the of the fluid shear is '1' <1 stress jump conditon region variables: 7» x x l : the 7o A l : The integration constant for fluid 7o A 2 : The integration constant for porou initial guess xx=[Al f_yy(0)] equation, Al=f_yyy(0) equation. 7o Soltot: Final solution [y f f_yy f_yyy] for fluid 7o Soltot2: Final solution [y g g_y g _ y y g _ y y y ] for porous f_y region. region. clc clear 7o c l e a r global old workspace Re D a x x l A 2 S o l t o t 7. T h e r e q u i r e d solutioins Repa=[13.2]; at porosityf=1; porositym=.9; iter=0; Soltot2 various Re n u m b e r are stored in 'Repa' vector Da=2.66e-3; beta=0.2; valfac=l; nn=length(Repa);Alpa=zeros(l,nn); height=.85; 7. h e i g h t = . 9 9 9 9 9 9 9 9 9 9 ; 7o I n i t i a l 7o a t for guess various Re of 7. f o r validation A l = - 1 a n d gammaa=l are found to be good f o r 'ODEsolve' numbers, ii=l:nn 7. g r i d Soltot=zeros(50,1);Soltot2=zeros(50,1); Re=Repa(ii) ; Al=-1; 7, A l l o w gammaa=l; N e w t o n ' s maximum i t e r a t i o n maxit=1000; store 7o u p d a t e d following i n the '/, r e g i o n times to be is in fluid/porous. 50 1000 tol=le-9; 7o L e t x x l v e c t o r xxl=[Al number gammaa]; initial guess inner maxdx=lelO; A l and gammaa=gl' ' (0) 'while' l o o p by Newton's and it method, is APPENDIX D. while CODE iter<maxit iter=iter+l; 7, S t a r t of computing FOR THE SIMILARITY SOLUTION & maxdx>tol maxdx=0; Jacobin matrix 'a' numerically [gO,nfun,hstart]=odesolve('eqfluid_shoot',0.height,[0 0 xxl(2)] le-6.height,le-5*l.height); z_l=porosityf*(xxl(1)*Re/porosityf-... Re/porosityf~2*(g0(2)~2-g0(l)*g0(3))); gp0_2=porositym*(gO(3)+l/sqrt(Da)*g0(2)*beta); gpO_3=porositym*Re*(z_l/Re+gO(2)/Da/Re*valfac-gO(l)*gO(3)+... g0(2)~2*(l-l/porositym"2)+g0(l)*gp0_2/porositym~2); A2=gp0_3/Re+(g0(2)~2-g0(l)*gp0_2)/porositym-... porositym/Da/Re*gO(2)*valfac; [ggO,nfun.hstart]=odesolve('eqporous_shoot'.height,1,... [gO(D g0(2) gp0_2],le-6,l,le-5*l,l-height); xxl(l)=xxl(l)+le-6*xxl(l); [gl,nfun.hstart]=odesolve('eqfluid_shoot0.height,[0 0 xxl(2)] le-6.height,le-5*l.height); z_10=porosityf*(xxl(1)*Re/porosityf-Re/porosityf~2*(gl(2)"2-... gl(l)*gl(3))); gpl_2=porositym*(gl(3)+l/sqrt(Da)*gl(2)*beta); gpl_3=porositym*Re*(z_10/Re+gl(2)/Da/Re*valfac-gl(1)*gl(3) +.. . gl(2)"2*(l-l/porositym~2)+gl(1)*gpl_2/porositym~2); A2=gpl_3/Re+(gl(2)~2-gl(l)*gpl_2)/porositym... -porositym/Da/Re*gl(2)*valfac; [ggl,nfun.hstart]=odesolve('eqporous_shoot'.height,1,... fgl(l) gl(2) gpl_2],le-6,l,le-5*l,l-height); pdeAl=((ggl(l)-l)-(ggO(l)-l))/(le-6*xxl(D) ; pdeA2=((ggl(2)-0)-(ggO(2)-0))/(le-6*xxl(D) ; xxl(l)=xxl(l)-le-6*xxl(l); xxl(2)=xxl(2)+le-6*xxl(2); [g2.nfun.hstart]=odesolve('eqfluid_shoot',0,height, [0 0 xxl(2)] le-6.height,le-5*l.height); z_12=porosityf*(xxl(1)*Re/porosityf-Re/porosityf"2*(g2(2)"2-... g 2(l)*g2(3))); gp2_2=porositym*(g2(3)+l/sqrt(Da)*g2(2)*beta); gp2_3=porositym*Re*(z_12/Re+g2(2)/Da/Re*valfac-... g2(l)*g2(3)+g2(2)"2*(l-l/porositym"2)+... g2(l)*gp2_2/porositym~2); A2=gp2_3/Re+(g2(2)"2-g2(l)*gp2_2)/porositym-... porositym/Da/Re*g2(2)*valfac; [gg2,nfun.hstart]=odesolve('eqporous_shoot'.height,1,... [g2(l) g2(2) gp2_2],le-6,l,le-5*l,l-height); APPENDIX D. CODE FOR THE SIMILARITY SOLUTION pdegaml=((gg2(l)-l)-(gg0(l)-l))/(le-6*xxl(2)); pdegam2=((gg2(2)-0)-(ggO(2)-0))/(le-6*xxl(2)); xxl(2)=xxl(2)-le-6*xxl(2); a=[pdeAl '/, E n d o f pdeA2; pdegaml computing Jacobin matrix 7, U s e G a u s s - J o r d o n m e t h o d % 'dx= [ D e l t a A l ; 7, m e e t t h e of to solve (increment it d i m e n s i o n agreement between dx=gaussj2(a); xxl=xxl+dx; for such form that sufficiently small, -gg0(2)+0]'; 'a' 'gaussj2' Deltagamma]' % The s o l u t i o n dx i s '/, d x i s pdegam2;-gg0(l)+l i.e., solution need to xx vectors, less than the update) be t r a n s p o s e d xx i s allowed updated to until tolerance. dx=dx'; maxdx=max(abs(dx)); end if maxdx<=tol I iter>=maxit Alpa(l)=xxl(l); else error ('Newtons method d i d not converge when Re i s 7 4.0f', 0 Re) end '/, B y n o w , '/, t o the right solve the 7. S o l v e t h e a=0; initial initial b=height; nm=n-l; for values value value tol=le-6; found, Below is a standard procedure problem: in fluid region: n=50; dtol=tol/nm; hmin=le-3*dxx; y0(l,:)=[0 xxl is problem- dxx=(b-a)/nm; hstart=dxx; x0(l)=a; initial hmax=dxx; 0 xxl(2)]; nfuntot=0; i=l:nm zO(l)=porosityf*(xxl(1)*Re/porosityf-Re/porosityf~2*(y0(l,2)"2-... y0(l,l)*y0(l,3))); ip=i+l; x0(ip)=x0(i)+dxx; [yO(ip,:),nfun.hstart]=odesolve('eqfluid_shoot',xO(i),x0(ip),y0(i,:) dtol.hstart,hmin,hmax); zO(ip)=porosityf*(xxl(1)*Re/porosityf-Re/porosityf2*(y0(ip,2)"2-... y0(ip,l)*y0(ip,3))); '/. gl"'(eta) nfuntot=nfuntot+nfun; end Sol=[x0' yO z 0 ' ] ; Soltot= [Soltot S o l ] ; % store a l l solutions 7. n o t e t h e first into column i s 'Soltot', a 0 vector. APPENDIX D. % Solve the CODE initial a=height; nm=n-l; value b=l; x(l)=a; region: n=50; dtol=tol/nm; hmin=le-3*dxx; y(l,:)=[gO(l) 103 SOLUTION problem- i n porous tol=le-6; dxx=(b-a)/nm; hstart=dxx; for FOR THE SIMILARITY hmax=dxx; g0(2) gp0_2]; nfuntot=0; i=l:nm z(l)=porositym*(A2*Re/porositym-Re/porositym~2*(y(1,2)"2-... y(l,l)*y(l,3))+l/Da*y(l,2)*valfac); ip=i+l; x(ip)=x(i)+dxx; [ydp,:),nfun.hstart]=odesolve('eqporous_shoot',x(i),x(ip),y(i,:),.. dtol,hstart,hmin,hmax); z(ip)=porositym*(A2*Re/porositym-Re/porositym~2*(ydp,2)"2-... y d p , l ) * y ( i p , 3 ) ) + l / D a * y ( i p , 2 ) * v a l f ac) ; nfuntot=nfuntot+nfun; end Sol2=[x' y z>]; Soltot2=[Soltot2 Sol2]; 7, s t o r e °/ n o t e the 0 °/, a l l solutions first column i s 'Soltot2', a 0 fder2=Soltot(l,5)*Re*Da~1.5/porositym~1.5 7 calculate fder3=Soltot(l,6)*Re*Da~2/porositynT2 % values 7, 7, into calculate the boundary layer thickness i n the vector. the transformed porous initial i n terms of... parameters region: alpha=sqrt(porositym/Da) xi=interpl(Soltot2(:,4),Soltot2(:,2),0.01) end •/.•/.•/.•/.•/.'/./.7./././././././.7.'/.'/. , , , , , , , , End of Main Program 7;/;/:/:/:/:/;/;/.n7,n7.nn%7.ny:/, 7o e q f l u i d _ s h o o t .m function z=eqfluid_shoot(x,g) 7, ' e q f l u i d _ s h o o t . m ' 7, f o r the fluid 7, x x l ( l ) = A l , global is region Al is the a function m-file (the integrated integration that defines the ODE e q u a t i o n form) constant. Re x x l D a z(l)=g(2); z(2)=g(3); z(3)=xxl(l)*Re-Re*(g(2)"2-g(l)*g(3)); 7;/.nn7,7,7,7.7.nn%7.% End of e q f i u i d _ s h o o t .m 7.7o7.r/.7;/x/.7.7.7.7.7.7x/:/;/.nn7.7. . APPENDIX D. CODE FOR THE SIMILARITY SOLUTION 7. e q p o r o u s _ s h o o t . m function z=eqporous_shoot(x,g) 7, ' e q p o r o u s _ s h o o t . m ' i s a f u n c t i o n m - f i l e t h a t d e f i n e s t h e % f o r t h e p o r o u s r e g i o n ( t h e i n t e g r a t e d form) '/, A2 i s t h e i n t e g r a t i o n c o n s t a n t . ODE e q u a t i o n g l o b a l Re A2 Da porosity=.9; valfac=l; valfac2=l; z(l)=g(2); z(2)=g(3); z(3)=A2*Re-valfac2*Re/porosity*(g(2)~2-g(l)*g(3))+l/Da*g(2)*porosity*valf 7X/X/.7J.7.7.7.mnra End o f e q p o r o u s _ s h o o t .m 7or/.7o%7oyo%7.yo/o7o7./c7./.%r/./o%%7o/o°/o 0 0 0 0 0 7o o d e s o l v e . m function [yb,nfun.hstart]=odesolve(f,a,b,ya,tol,hstart,hmin,hmax) 7, ' o d e s o l v e . m ' i s a f u n c t i o n m - f i l e t h a t u s e s t h e "embedded" R u n g e - K u t t a 7, method t o s o l v e IVP o f s i m u l t a n e o u s f i r s t o r d e r d i f f e r e n t i a l 7. of t h e equations form: */. 7. d y l / d x = f l ( x , y l , y 2 , . . .ym); dy2/dx=f2(x,yl,y2,...ym); ' 7. 7. dym/dx=fm(x,yl , y 2 , . . .ym); 7, where I V c o n d i t i o n s a r e : 7. y l ( a ) = a l p h a l , y 2 ( a ) = a l p h a 2 , . . . , ym(a)=alpha(m) ; % a < x < b. 7, The u s e r - p r e s c r i b e d a c c u r a c y of t h e s o l u t i o n i s c h e c k e d by a d a p t i v e l y 7. a d j u s t i n g t h e s t e p - s i z e of the 7o g l o b a l t o l e r a n c e 7, I n p u t f is s o l u t i o n s u b i n t e r v a l to ensure t h a t met. arguments: =dummy name of r h s functions = s t a r t o f s o l u t i o n i n t e r v a l x=a a =end of s o l u t i o n i n t e r v a l x=b b y tol = i n i t i a l v a l u e s of v e c t o r y [ l : m ] a t hstart hmin =initial =minimum step-size step-size hmax =maximum step-size a 7o Output = g l o b a l a c c u r a c y r e q u i r e d at arguments: x=b x=a a APPENDIX D. CODE FOR THE SIMILARITY 7, °L yb nfun =values of =number o f 7c hstart =last v e c t o r y [ l : m ] a t x=b rhs function evaluations successful value 105 SOLUTION needed o f hnew b e f o r e return m=length(ya); 7, I n i t i a l i z e v a r i a b l e s , hold=hstart; xold=a; 7, L o o p u n t i l yb i s while where yoldfya; xold is a scalar, nfun=0; determined to yold is a vector success=0; required tolerance success==0 7. S i n c e y a i s 7o w o u l d b e a vector 6 times with l e n g t h m, n u m b e r m everytime call to of f u n c t i o n of evaluations 'rkck.m' nfun=nfun+6*m; [ynew.ydiff]=rkck(f,xold,yold,hold); gamma=(tol*hold/(ydiff*(b-a)))"0.25; hnew=0.8*gamma*hold; 7, I f if gamma<l, adjust step size and l o o p t o recalculate ynew gamma<l if hnew<hold/10 hnew=hold/10; end if hnew<hmin fprintf('\n odesolve failed at x=7,8.3e XnXn^old',xold) error end hold=hnew; 7. I f gamma>=l, accept ynew and add h o l d to xold else if hnew>5*hold hnew=5*hold; end if hnew>hmax hnew=hmax; end if xold+hold<b xold=xold+hold; 7o I f hold=hnew; xold+hold>=b, yold=ynew; calculate yb and return APPENDIX D. CODE FOR THE SIMILARITY 106 SOLUTION else success=l; hstart=hnew; hold=b-xold; nfun=nfun+6*m; [yb.ydiff]=rkck(f ,xold,yold',hold); end end end mr/.my.r/:/x/.nn End of odesoive.m yx/x/x/x/x/x/x/x/x/,yx/x///, % rkck.m function [y2, y d i f f ] = r k c k ( f , x l , y l , h ) % 'rkck.m' % method, is a function m-file which */, u s i n g t h e is due 5th-order 5th-order methods to that Cash and Karp formula, is evaluates (1990). while the determined to the "embedded" The difference help update Runge-Kutta solution is advanced between 4th the the step-size and in '/. ' o d e s o i v e . m ' . °/ I n p u t 0 argument s : °/ f = dummy name °/, xl = value 0 of % yl = values */, h = solution % Output of the scalar of rhs x at vector functions start y at interval of start the equation solution of interval of solution sets to be solved. interval step-size arguments: 7. y2 °/, ydiff= = value rhs of y at function % k l , k 2 , . . . k 6 are end of solution evaluations vectors with interval needed length the same a s y l and y2. kl=h*feval(f,xl,yl); k2=h*feval(f,xl+h/5,yl+kl/5); k3=h*feval(f,xl+3*h/10,yl+3*kl/40+9*k2/40); k4=h*feval(f,xl+3*h/5,yl+3*kl/10-9*k2/10+6*k3/5); k5=h*feval(f,xl+h,yl-ll*kl/54+5*k2/2-70*k3/27+35*k4/27); k6=h*feval(f,xl+7*h/8,yl+1631*kl/55296+175*k2/512+575*k3/13824+... 44275*k4/110592+253*k5/4096); y2=37*kl/378+250*k3/621+125*k4/594+512*k6/1771+yl; ydiff0=abs((37/378-2825/27648)*kl+(250/621-18575/48384)*k3+... (125/594-13525/55296)*k4-277*k5/14336+(512/1771-1/4)*k6); */, T h e v a l u e % between of each scalar element 'ydiff' of is vectors taken as y2(5th the maximum o f order) the and y 2 ( 4 t h differences order). It is APPENDIX D. CODE FOR THE SIMILARITY 7. i m p o r t a n t to give ydiffO '/, y d i f f = y 2 ( 5 t h ) - y 2 ( 4 t h ) , 7o a t 'odesolve at line a agebraic the 48', form, computer which i s SOLUTION will otherwise, give assumed to if just put d i v i d e d by zero warning be errors. the roundoff ydiff=max(ydiff0); End of /o/o/o/t/o/o/o/ti/ti/ti/o/o/o/o/o/o/ti/o/ti/o/o/o/o/o/o/o/fi/ti/o rkck.m IttLlnlnLltiliilfiLLItiLLIitLLLLLtttt'iltiLloltiLltiLL 7o g a u s s j 2 . m function 7 0 7. x=gaussj2(A) 'gaussj2' uses the Gauss-Jodan method to solve simutaneous linear equations. 7, T h e e l i m i n a t i o n l o o p '/, i f 7. t h e i~=k is conducted makes e l i m i n a t i o n o n l y diagonal. Note t h a t 7o b e e n n o r m a l i z e d t o 1, for n times conducted during the therefore n-1 and the o n e l e m e n t s who a r e not on process the final the instead of A(k,k) result is has x=A(i,n+l). n=size(A,1); '/, p a r t i a l pivot selection 7, i n f i l e \ g a u s s j .m i n t h e for with scaling; submited the assignment. k=l:n pivot=k; for maxaij=zeros(1,n); spivot=0; i=k:n maxaij(i)=abs(A(i,k)); for j=k:n if abs(A(i,j))>=maxaij(i) maxaij(i)=abs(A(i,j)); end end abssi=abs(A(i,k)/maxaij(i)); if spivot <abssi spivot=abssi; • pivot=i; end end if pivot~=k for j=k:n+l temp=A(pivot,j); A(pivot,j)=A(k,j) ; A(k,j)=temp; end end same a s described APPENDIX D. CODE FOR THE SIMILARITY 7, GJ e l i m i n a t i o n SOLUTION 108 w i t h n o r m a l i z a t i o n of A(k,k) . akk=A(k,k); for j=l:n+l A(k,j)=A(k,j)/akk; end for i=l:n i f i~=k aik=A(i,k); for j=l:n+l A(i,j)=A(i,j)-aik*A(k,j); end end end end if A(n,n)==0 errorCzero p i v o t c o e f f i c i e n t e n c o u n t e r e d i n l a s t row!') end 7. F i n a l result x=A(l:n, n+1); y»7o7oyoyoyo7o7o%7o7a7o7.7o7.7o%°/o7.°/o°/o°/o°/o°/»7.%7o°/o D.2 End o f g a u s s j 2. m yX7X7XVXVXV/X7XW/XW/XyXVh°IX°h The Reverse Solver The code is presented as the main script file following by the various supporting function files. • stability_shootMain_auto.m • solveatconstDaMain.m • coeff.m • Boundeqn.m • solveDa2.m • coxsolvefun.m • brent, m • eqfmid.m APPENDIX D. CODE FOR THE SIMILARITY SOLUTION • eqporous.m • odesolvel.m • odesolve2.m • rkckl.m • rkck2.m • solvetoanyposition.m • solvetoanyposition_p.m • eqstability.m • eqstability_p.m APPENDIX D. CODE FOR THE SIMILARITY SOLUTION V °F"/ °/ °/ V V °/ V V °/ V °/ V V °/ V °/ °/ V °/ V V °/ V V °/ °/ V °/ V V V °/ V V °/ V V V V V V V V °/ V V V °/ V °/ V / / / / / / / / / / / / / / / / / / / / / / 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 / Q / O / O / O / O / O / O / O / O / Q / O / D / O / O / O / O / O / O / ^ % Similarity S o l u t i o n - A Reverse Last updated - July Solver & Stability analysis % 06 2 0 0 4 % 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o % Main program f i l e '/, a n d S t a b i l i t y of the Reverse Analysis. Solver.to solve for Note t h a t the eigenvalues % Parameters ( a l l dimensionless): '/, B l _ r : the i n t e g r a t i o n constant '/, lambda_r: 7, 0 /, eigenvalue R e l : the /,Totf l u i d _ c a l : 0 the Re n u m b e r the targetDa: °/ the targeted actual real in fluid roots region. solved solutions r e s c a l i n g constant D a l : the 0 s t a b i l i t y eqn solutions are $\sigma$. intermediate b p o s i t i o n l : the of the steady obtained Da at 'q', which equals sqrt (epsilon/Da) . solution is solved. s o l v e d Da, an i n t e r m e d i a t e value. % f2der: the initial guess, equals to lambda_yy(0). '/. f3der_ini: the initial guess, equals to psi_yyy(0) '/, beta: the e m p i r i c a l shear '/, inter: the dimensionless height '/, Totfluid: final solution for fluid '/, Totfluid2: final solution for porous stress jump of constant fluid region 'h'. region. region. clear clc global B l _ r lambda_r R e l T o t f l u i d _ c a l targetDa Dal f2der Tot_Totfluidcal porositym beta bpositionl... inter... record_xx Tot_Totfluidcal_p record_xx_p show_tot= [ ] ; load ini_guess.dat % to l o a d the % to ini_guess % guess '/, t h e contains integration program is obtain the previous guess of constant designed to s o l u t i o n as solution for lambda_zz(0), B l _ r and e i g e n v a l u e be able to the a slightly initial different lambda_zzz(0), guess Re. and the lambda_r. self-iterate to another $Re$. size=length(ini_guess(:,1)); tol=le-6; maxit=300; beta=0.2; porositym=0.9; inter=0.8; targetDa=.0026; guess=ini_guess(1,4:5); for i=l:size Tot_Totfluidcal=[]; record_xx=[];Tot_Totfluidcal_p=[]; f3der_ini=ini_guess(i,2:3); [Totfluid, Totfluid2, resultv_l]=solveatconstDaMain... (targetDa, °/o c o m p u t e iter=0; maxdx=lel0; record_xx_p=[]; f2der=ini_guess(i,1); f2der, s o l u t i o n at f3der_ini); targeted Da. APPENDIX D. CODE FOR THE SIMILARITY Bl_r=guess(l) ; lambda_r=guess(2); x=[lambda_r while '/, s t a r t Bl_r] ; iter<maxit 11 SOLUTION compute °/ e q u a t i o n eigenvalue of using Newton's 0 Stability method & maxdx>tol iter=iter+l; maxdx=0; a=coeff('Boundeqn',x); dx=gaussj2(a) x=x+dx'; maxdx=max(abs(dx)) ; Tot_Totfluidcal=[] ; record_xx=[]; Tot_Totfluidcal_p=[] ; record_xx_p=[] ; end if maxdx<=tol fprintf('\n Successful s o l u t i o n i n °/ 3.0f 0 iterations \n\n', iter) else error('Newtons method d i d not converge') end show=[resultv_l' x(2) show_tot=[show_tot; x(l) guess]'; guess=[x(2) x(l)]; % result output. show]; dlmwrite('show_tot.xls',show_tot,'t'); f p r i n t f C7„4.13f\n' f printf ,show_tot) C\n') end y//.y//.yx/////.y.y.y.7.y.y.y.y.y:/.y.yx/././././././.7. , , , , , , Mam end of Program 7.7.7.7.y.y.y.y.y.y/.y.y/.y.y.y.y.y.y.y.y.y.y. 7o s o l v e a t c o n s t D a M a i n . m function [Totfluid, f2der, global Totfluid2, resultv_l]=solveatconstDaMain(targetDa,... f3der_ini) Dal Rel porositym Totfluid 7. ' ' s o l v e a t c o n s t D a M a i n . m ' ' i s 7, t h i s method 7o d e r i v a t i v e , assumes and grid_nfluid=60; at the constant the Totfluid2 driver second same t i m e grid_nporous=40; file get the f3rdder = 7 the 0 above solveDa2.m at solutions 7, g r i d _ n s h o u l d b e z=0, at set but to find target Da at same t i m e the f i n d i n g methods numbers. 34.27405585,1,le-3,le-6); fzero('solveDa2',[f3der_ini(1),f3der_ini(2)],options); two r o o t third solveDa2. X[f3rdder,irr]=brent('solveDa2',34.23405585, = optimset('Tolx',le-13); for derivative °/, f i l e options beta c o u l d be u s e d alternatively. in APPENDIX % fzero D. CODE function FOR THE SIMILARITY i s a matlab built SOLUTION i n function to find stress_inter_l=-Totfluid(grid_nfluid,5); % the stresswall_l=-Totfluid2(grid_nporous,5); % -the s h e a r resultv_l=[Dal f2der f3rdder 7.7.7.7.7.7o7.7o/o7o/.y./.%°/o7o7o%7o%7. 0 0 0 7o c o e f f .m i s t h e function global Jacobian root '-f_yy(chi) stress'-f_yy(l)' stresswall_l]'; end of solveatconstDaMain.m augmented real shear % output y„7:/:/.7:/.n%7:/.7./7.n7:/.7.%7.7././. 0 0 0 0 matrix a=coeff(f,x) B l _ r lambda_r n=length(x); for Rel stress_inter_l a single interface np=n+l; delx=le-4*x; i=l:n Bl_r=x(2); lambda_r=x(1); a(i,np)=-feval(f,i,x); for j=l:n xtemp=x(j); x(j)=x(j)+delx(j); Bl_r=x(2); lambda_r=x(l); ftemp=feval(f,i,x); x(j)=xtemp; a(i,j)=(ftemp+a(i,np))/delx(j); end end 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / j n e n d OI /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/oA 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / CuSII .m /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 7o B o u n d e q n . m 7. e q u a t i o n s function t o be e v a l u a t e d by coeff.m and s o l v e d by Newton's method: y=Boundeqn(i,x) g l o b a l ' B l _ r lambda_r porositym 7o i n t e g r a t i o n R e l D a l T o t f l u i d _ c a l T o t f l u i d _ c a l m e d B2 b e t a ... inter f r o m y=0 t o y=interface [W_f,nfun.hstart]=odesolve1('eqstability',0,inter,[0 0 1],le-6,inter,... le-6*inter,inter); W_f_3der=Bl_r+Rel*lambda_r*W_f(2)+Rel*(W_f(3)*Totfluid_cal(2)-... 2*Totfluid_cal(3)*W_f(2)+Totfluid.cal(4)*W_f(1)); 7o E v a l u a t e the boundary conditions at the interface. W_p(l)=W_f(1) ; W_p(2)=W_f(2); W_p(3)=porositym*(W_f(3)+beta/sqrt(Dal)*W_f(2)); W_p_3der=porositym*W_f_3der+(Rel*Totfluid_cal(2)-... Totfluid_cal(2)*porositym*Rel)*W_f(3)+(2*porositym*Rel*Totfluid_cal(3) Rel*Totfluid_cal(2)*beta/sqrt(Dal)+porositym/Dal-... APPENDIX D. CODE FOR THE SIMILARITY SOLUTION 2*Totfluid_cal(3)*Rel/porositym)*W_f(2)... +(-porositym*Rel*Totfluid_cal(4)+Rel*(Totfluid.cal(4)+... beta/sqrt(Dal)*Totfluid_cal(3)))*W_f(1)-(porositym-1)*Rel*W_f(2)*lambda_ B2=W_p_3der-porositym/Dal*W_p(2)-Rel/porositym*(Totfluid_cal(2)*W_p(3)... -2*Totfluid_cal(3)*W_p(2) +(Totfluid_cal(4)+beta/sqrt(DaD*. .. Totfluid_cal(3))*porositym*W_p(l))-Rel*lambda_r*Totfluid_cal(3); % Integrate from the interface to the porous wall. [W_p_atl,nfun,hstart]=odesolve2('eqstability_p',inter,1,W_p,le-6,(1-inter),. le-6*(1-inter),1-inter); if i==l y=(W_p_atl(D); else y=(W_p_atl(2)); end y.y;/.y.7.y.y.y.7:/:/.7.7.7.7:/.y.n%7.n7.77.7.7.7.7.7.7. 0 7. ' ' s o l v e D a 2 . m ' ' 7a b y b r e n t ' s 7. f o r m u l a root function global at of too. 7, g r i d _ n p o r o u s be to the Note be file to mmmy.mmnm /. /.y. , Boundeqn.m find the 'bpositionl' target solutions. to get larger Da upon position c a n be In t h i s detailed (41 is success good of routine enough), the help the set the which predicted to solution, at Rel Totfluid interface Totfluid2 f2der solutions Nosolnf=zeros(grid_nfluid,6); A l A2 I n d e x O p h i bpositionl... Soltot2=zeros(grid_nporous,1); Nosolnp=zeros(grid_nporous,6); Totfluid2=zeros(grid_nporous,0); Totfluidm=zeros(grid_nfluid,0); Totfluidm2=zeros(grid_nporous,0); Al=gamma; optimset('Tolx',le-6); y.bpositionl 7, s o l v e = f z e r o (' c o x s o l v e f u n ' , [7.746,7.74592]) position ; 'q' [bpositionl,ir]=brent('coxsolvefun',18.6052,150,1,.0000002,le-6) if obtained targetDa... Totfluid=zeros(grid_nfluid,0); = Index0==l the solutions iterations. grid_nporous=40; Soltot=zeros(grid_nfluid,1); the grid_nfluid, IndexO=0; grid_nfluid=60; ^options by determine detail , f_y(0)=0 zzz=solveDa2(gamma) porositym phi_l Dal driver f i n d i n g method, guess computed 7, w i l l the of ''bpositionl=sqrt(porositym/target_Da)'' 7o e f f i c i e n t 7. i s is end APPENDIX D. CODE FOR THE SIMILARITY 114 SOLUTION IndexO=0; fprintfCno root found Totfluid=[Totfluid \ n \ n f2der=7.8.4f f3der=°/.8.4f \ n ' , fsder, Al) Nosolnf]; Totfluid2=[Totfluid2 Nosolnp]; Totfluidm=[Totfluidm Nosolnf]; Totfluidm2=[Totfluidm2 Nosolnp]; else fsder=f2der; f3der=Al; 7, s o l u t i o n contains D a a n d Re ( i n t e r m e d i a t e values) : Dal=porositym/bpositionl~2; Rel=porositym~l.5/sqrt(Dal)*phi_l(1) ; 7, c o m p u t e a=0; detail nm=n-l; lambda(z) tol=le-6; dxx=(b-a)/nm; hstart=dxx; xO(l)=a; for solution b=interface; n=grid_nfluid; dtol=tol/nm; hmin=le-6*dxx; hmax=dxx; y O ( l , : ) = [0 0 f s d e r ] ; nfuntotO; i=l:nm zO(l)=Al-(yO(l,2)~2-yO(l,l)*yO(l,3)); ip=i+l; xO(ip)=xO(i)+dxx; [yO(ip,:),nfun.hstart]=odesolve('eqfluid',xO(i),xO(ip),yO(i,:),... dtol.hstart,hmin,hmax); z O ( i p ) = A l - ( y O ( i p , 2 ) - 2 - y O ( i p l ) * y O ( i p , 3 ) ) ; 7. ) gl"'(eta) nfuntot=nfuntot+nfun; end 7o r e s c a l e lambda(z) back to the desired form: f (y) xO=xO/sqrt(porositym/Dal); y0(:,l)=yO(:,1)*sqrt(porositym/Dal)/Rel; y0(:,2)=y0(:,2)*porositym/Dal/Rel; y0(:,3)=y0(:,3)*porositym~1.5/Rel/Dal"1.5; zO=zO*porositym~2/Rel/Dal~2; Sol=[xO' yO z O ' ] ; Soltot= [Soltot S o l ] ; Totfluid=[Totfluid Soltot]; Soltot=zeros(grid_nfluid,1); '/, c o m p u t e detail a=interface; nm=n-l; dxx=(b-a)/nm; hstart=dxx; x(l)=a; for solution psi(z) b=bpositionl; n=grid_nporous; dtol=tol/nm; hmin=le-6*dxx; y(l,:)=phi; tol=le-6; hmax=dxx; nfuntot=0; i=l:nm z(l)=(A2-(y(l 2)-2-y(l,l)*y(l,3))+y(l,2)); ) APPENDIX D. CODE FOR THE SIMILARITY 115 SOLUTION ip=i+l; x(ip)=x(i)+dxx; [y(ip,:), nfun.hstart]=odesolve('eqporous',x(i),x(ip),y(i,:),... dtol,hstart,hmin,hmax); z(ip)=(A2-(y(ip,2)~2-y(ip,l)*y(ip,3))+y(ip,2)); nfuntot=nfuntot+nfun; end °/ r e s c a l e 0 psi(z) back to the desired form: g(y) x=x/sqrt(porositym/Dal); y(:,l)=y(:,1)*sqrt(porositym/Dal)*porositym/Rel; y(:,2)=y(:,2)*porositym~2/Dal/Rel; y(:,3)=y(:,3)*porositym~2.5/Rel/Dal~l.5; z=z*porositym~3/Rel/Dal~2; Sol2=[x' y z>] ; Soltot2=CSoltot2 Sol2]; Totfluid2=[Totfluid2 Soltot2]; Soltot2=zeros(grid_nporous,1); end zzz=targetDa-Dal % the difference 7/oy.yoyc%7o7o7/oyoyoyoyoyoyoy.yoy«y.y.yo7o7.y.y/.o/o7o°/. c between the end of intermediate s o i v e D a 2 .m and t a r g e t e d da. y,7,7.7.7.7.7o7.7o7o7.7.7.7c7o7o7o/oy.%yoyoy.y.7«yoyoy.7. o 7, c o x s o l v e f u n .m function z=coxsolvefun(omega) 7, R o o t f i n d i n g global for b p o s i t i o n l where A l A2 p o r o s i t y m p h i _ l f_y(l)=0 interface f2der or phi psi_z(q)=0. IndexO b e t a inter interface=inter*omega; [lambda,nfun,hstart]=odesolve('eqfluid',0,interface,[0 0 f2der],le-6,... interface,le-6*interface,interface); 7. E v a l u a t e interface conditions lambda3der=Al-lambda(2)"2+lambda(l)*lambda(3); phi(l)=lambda(l)/porositym; phi(2)=lambda(2)/porositym; phi(3)=lambda(3)+l/sqrt(porositym)*beta*lambda(2); phi3der=lambda3der+lambda(2)"2-lambda(l)*lambda(3)+phi(l)*phi(3)-... phi(2)~2+phi(2); A2=phi3der+phi(2)~2-phi(2)-phi(1)*phi(3); % Solve to porous wall y=l [phi_l,nfun.hstart]=odesolve('eqporous',interface, omega, phi, le-6,... omega-interface,le-6*(omega-interface),omega-interface); Datem=porositym/omega~2; Retem=porositym~l.5/sqrt(Datem)*phi_l(1); APPENDIX D. CODE FOR THE SIMILARITY 116 SOLUTION z=phi_l(2)*porositym~2/Datem/Retem i f a b s ( z ) > = 10000000 Index0=l; end myx/.y.%my//x/x///.y.y.y.%y.7.%yx/.7. end of c o x s o l v e f u n . m %7,7x/:/.n7,7.7;/;/.7xm7.7x/.n% 7. b r e n t .m function [R,ir]=brent(f,xi,xf,nr,dx,tol) 7, ' B r e n t ' f u n c t i o n u s e s B r e n t ' s method t o f i n d m u l t i p l e r e a l r o o t s f o r an 7o 7. 7. 7. a l g e b r a i c e q u a t i o n . The i n c r e m e n t a l s e a r c h p r o c e d u r e i s e m p l o y e d t o f i n d b r a c k e t o f a r e a l r o o t o r a s i n g u l a r i t y . The b i s e c t i o n method i s u s e d t o n a r r o w t h e i n t e r v a l as w e l l t o f i n d w h e t h e r s i n g u l a r i t y i s p r e s e n t . I f no s i n g u l a r i t y f o u n d , B r e n t ' s method w i l l be u s e d t o " p o l i s h " t h e r o o t v a l u e . 7. I n p u t a r g u m e n t s : 7, f - dummy name of a l g e b r a i c e x p r e s s i o n whose r o o t s 7. x i - start of s e a r c h are sought interval 7, x f - end o f s e a r c h 7, 7, n r - number of r e a l r o o t s s o u g h t dx - i n i t i a l s e a r c h i n c r e m e n t u s e d t o l o c a t e b r a c k e t f o r a s i n g l e 7, interval t o l - absolute root tolerance 7. O u t p u t a r g u m e n t s : 7. R - A r r y of r o o t s 7, found i r - number of r o o t s f o u n d i n i n t e r v a l 7o I n i t i a l i z e t h e root. or accuracy [ x i xf] variable g l o b a l IndexO xl=xi; y l = f e v a l ( f , x l ) ; ir=0; 7, F i n d i f t h e f i r s t if end yl==0 ir=ir+l; r o o t found i f R(ir)=xl; f(xi)=0 xl=xl+dx; y l = f e v a l ( f , x l ) ; 7. Loop t o f i n d r e m a i n i n g r o o t s when i r < = n r and x<=xf while ir<nr 7o I n c r e m e n t a l s e a r c h t o l o c a t e i n t e r v a l 7. s i n g u l a r i t y interval_found=0; while interval_found==0 [ x l , x3] h a v i n g s i n g l e r e a l r o o t or APPENDIX D. CODE FOR THE SIMILARITY SOLUTION x3=xl+dx; if (dx>0 & x 3 > x f ) fprintf C\n I (dx<0 & x 3 < x f ) Only °/ 5 . Of r o o t s 0 found. \n',ir) return end y3=feval(f,x3); if yl*y3>0 xl=x3; if yl=y3; 7, f o r (abs(IndexO-l)<le-6) solvepositionbauto ir=l; 7. f o r solvepositionbauto only R(ir)=0; 7. f o r solvepositionbauto only return 7. f o r solvepositionbauto only 7. f o r end solvepositionbauto only only else interval_found=l; end end 7, B i s e c t i o n to narrow interval for three times 7» s i n g u l a r i t y root_found=0; singularity=0; dyold=abs(y3-yl); for i=l:3 x2=(xl+x3)/2; if y2=feval(f,x2); yl*y2>0 xl=x2; yl=y2; x3=x2; y3=y2; else end end x2=(xl+x3)/2; y2=feval(f,x2); dy=abs(y3-yl); if dy>dyold singularity=l; else root_found=0; end 7c B r e n t ' s while method to polish root_found==0 r=y2/y3; & s=y2/yl; the root singularity==0 t=yl/y3; q=(r-l)*(s-l)*(t-l); p=s*(t*(r-t)*(x3-x2)-(1-r)*(x2-xl)); x4=x2+p/q; y4=feval(f,x4); if abs((p/q)/x4)<=tol ir=ir+l; break end R(ir)=x4; root_found=l; and determine existence APPENDIX D. if CODE FOR THE SIMILARITY I (dx>0 & x l < x 4 & x4<x2) x3=x2; elseif y3=y2; x2=x4; (dx<0 & x l > x 4 & x4>x2) y2=y4; I (dx>0 & x2<x4 & x4<x3) xl=x2; yl=y2; x2=x4; 1 SOLUTION (dx<0 & x2>x4 & x4>x3) y2=y4; end end xl=x3; yl=y3; end 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 6 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / „„J /JJJJJJJJJJJJJJJJJJJJJJJJJJJt/JJt end function '/, 0 °/. of „ 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / 0 / IJJJJJJJJJJJJJJJJJJJJJJJJJJJi b r e n t .m z=eqfluid(x,g) 'eqflu-id.m' i s 7 the T ^ ^ „ „ j _ J2 n original a function m-file transformed from the that defines the equations set for 3 r d o r d e r ODE xxl(l)=Al global A l z(l)=g(2); z(2)=g(3); z(3)=Al-(g(2)~2-g(l)*g(3)); y>vhv/xv/xv/xv/xv/xv/°/°/x function °/ 0 global eqfluid.m xx/iyoy.y.y.y.xy.yo/cyoyoy.xyoXxxxxxxxy.yoyoXy.ycyo'/oyoyo z=eqporous(x,g) 'eqporous.m' % the end of original is a function m-file transformed from the that defines the equations set for 3 r d o r d e r ODE A2 valfac=l; valfac2=l; z(l)=g(2); z(2)=g(3); z(3)=A2-valfac2*(g(2)~2-g(1)*g(3))+g(2)*valfac; eqporous.m function % [yb,nfun.hstart]=odesolvel(f,a,b,ya,tol,hstart,hmin,hmax) 'odesolvel.m' % here it h lo h h h /a h L h h h h h h lo h h lo la h lo lo L h h lo /a h calls is a slight $rkckl$ m o d i f i c a t i o n of routine instead of 'odesolve.m' $rkck$ m=length(ya); hold=hstart; xold=a; yold=ya; nfun=0; success=0; APPENDIX while D. CODE FOR THE SIMILARITY SOLUTION 119 success==0 nfun=nfun+6*m; [ynew,ydiff]=rkckl(f,xold,yold,hold); gamma=(tol*hold/(ydiff*(b-a)))"0.25; hnew=0.8*gamma*hold; if gamma<l if hnew<hold/10 hnew=hold/10; end if hnew<hmin fprintf('\n odesolve failed at x=%8.3e \ n \ n , x o l d x o l d ) error end hold=hnew; else if hnew>5*hold hnew=5*hold; end if hnew>hraax hnew=hmax; end if xold+hold<b xold=xold+hold; hold=hnew; yold=ynew; else success=l; hstart=hnew; hold=b-xold; n f u n = n f un+6 *m; [yb,ydiff]=rkckl(f,xold,yold,hold); end end end %7o%7.y././.%%/.7./o/.7.7.%7.%/.7o7o/o/. oo o function o o o o 0 end of odesoivei.m t [yb,nfun,hstart]=odesolve2(f,a,b,ya,tol,hstart,hmin,hmax) 7, ' o d e s o l v e 2 . m ' is 7o h e r e $rkck2$ routine it y:/;/;/;/.n7.y.y.y.y.7:/:/:/:/.n7.y.7.7.y.7;/:/.7.n7:/:/.y;/ calls a slight m o d i f i c a t i o n of instead of 'odesolve.m' $rkck$ m=length(ya); hold=hstart; xold=a; yold=ya; nfun=0; success=0; APPENDIX while D. CODE FOR THE SIMILARITY SOLUTION 120 success==0 nfun=nfun+6*m; [ynew.ydiff]=rkck2(f,xold,yold,hold); gamma=(tol*hold/(ydiff*(b-a)))"0.25; hnew=0.8*gamma*hold; if gamma<l if hnew<hold/10 hnew=hold/10; end if hnew<hmin fprintf C\n odesoive failed at x=°/ 8.3e \ n \ n , x o l d ' , x o l d ) 0 error end hold=hnew; else if hnew>5*hold hnew=5*hold; end if hnew>hmax hnew=hmax; end if xold+hold<b xold=xold+hold; hold=hnew; yold=ynew; else success=l; hstart=hnew; hold=b-xold; nfun=nfun+6*m; [yb,ydiff]=rkck2(f,xold,yold,hold); end end end o /a /o7o7o%7 °t/o77 .o7o7o7«7o7o7°o ./°o /o /7o°o / 0000 00 function global 7. ' r k c k l . m ' 7o t h e end of odesolve2.m 7o7o7o7o7o77o7o7»7o7.o7o7o7o7o7o7o77o77o7o7o7o777o77o7o777o7o [y2, y d i f f ] = r k c k l ( f , x l , y l , h ) Totfluid_cal make u s e d perturbed of solutions the base F(y) solutions in fluid f(y) region. Totfluid_cal=solvetoanyposition(xl); kl=h*feval(f,xl,yl); Totfluid_cal=solvetoanyposition(xl+h/5); in calculating APPENDIX D. CODE FOR THE SIMILARITY SOLUTION k2=h*f e v a K f , x l + h / 5 , y l + k l / 5 ) ; Totfluid_cal=solvetoanyposition(xl+3*h/10) ; k3=h*feval(f,xl+3*h/10,yl+3*kl/40+9*k2/40) ; Totfluid_cal=solvetoanyposition(xl+3*h/5); k4=h*fevaKf,xl+3*h/5,yl+3*kl/10-9*k2/10+6*k3/5); Totfluid_cal=solvetoanyposition(xl+h); k5=h*feval(f,xl+h,yl-ll*kl/54+5*k2/2-70*k3/27+35*k4/27); Totfluid_cal=solvetoanyposition(xl+7*h/8); k6=h*f evaKf,xl+7*h/8 yl+1631*kl/55296+175*k2/512+575*k3/13824+... ) 44275*k4/110592+253*k5/4096); y2=37*kl/378+250*k3/621+125*k4/594+512*k6/1771+yl; ydiff0=abs((37/378-2825/27648)*kl+(250/621-18575/48384)*k3+... (125/594-13525/55296)*k4-277*k5/14336+(512/1771-l/4)*k6); ydiff=max(ydiff0); LLLLLLLLLLLLLLLLLLLLLLL function global of rkckl.m LLhhhhLLLLLLLhhhhLhLLhhLhhhLlnLLLhhh ydiff]=rkck2(f,xl,yl,h) Totfluid_calmed '/, ' r k c k 2 . m ' °/, t h e [y2, end make u s e d perturbed of solutions the base G(y) solutions i n medium f(y) and g(y) in calculating region. Totfluid_calmed=solvetoanyposition_p(xl); kl=h*feval(f,xl,yl); Totfluid_calmed=solvetoanyposition_p(xl+h/5); k2=h*feval(f,xl+h/5,yl+kl/5); Totfluid_calmed=solvetoanyposition_p(xl+3*h/10); k3=h*feval(f,xl+3*h/10,yl+3*kl/40+9*k2/40); Totfluid_calmed=solvetoanyposition_p(xl+3*h/5); k4=h*fevaKf,xl+3*h/5,yl+3*kl/10-9*k2/10+6*k3/5); Totfluid_calmed=solvetoanyposition_p(xl+h); k5=h*fevaKf,xl+h,yl-ll*kl/54+5*k2/2-70*k3/27+35*k4/27); Totfluid_calmed=solvetoanyposition_p(xl+7*h/8); k6=h*fevaKf,xl+7*h/8,yl+1631*kl/55296+175*k2/512+575*k3/13824+... 44275*k4/110592+253*k5/4096); y2=37*kl/378+250*k3/621+125*k4/594+512*k6/1771+yl; APPENDIX D. CODE FOR THE SIMILARITY 122 SOLUTION ydiff0=abs((37/378-2825/27648)*kl+(250/621-18575/48384)*k3+... (125/594-13525/55296)*k4-277*k5/14336+(512/177l-l/4)*k6); % The v a l u e 7. b e t w e e n of scalar each 7. i m p o r t a n t 'ydiff element to of give ydiffO , the 7, a t 48', at line taken as y2(5th a agebraic 7c y d i f f = y 2 ( 5 t h ) - y 2 ( 4 t h ) 'odesoive is vectors computer which i s the maximum o f order) form, will otherwise, give assumed to the and y 2 ( 4 t h if differences order). just It is put d i v i d e d by zero warning be errors. the roundoff ydiff=max(ydiffO); %7o7o7o7o7o7»7o7»7.77.7o7.7o7o7o7o7.7o7.7o7o 7o end of r k c k 2 .m 7O7O7O7»7.7O7O7O7O7O7O%7.7»77»7O7O7.7O7O7O7D7O7O7O7O7O7O7O7O7.7O7»7O solvetoanyposition.m function 7. t o Totfluid_cal=solvetoanyposition(xx) get base 7o b a s i c a l l y global if f2der s o l u t i o n at it is an any d e s i r e d initial value position i n the fluid region, solver. Al bpositionl Rel inter Tot_Totfluidcal record_xx length(Tot_Totfluidcal)==0 hit=0; else hit=length(Tot_Totfluidcal(:,1)); end if hit>=10 for i=l:hit if xx==record_xx(i) Totfluid_cal=Tot_Totfluidcal(i,:); return end end end if xx==0 x0=0; y0=[0 0 f2der]; zO=Al-(yO(2)-2-yO(l)*yO(3)); else [yO,nfun.hstart]=odesolve('eqfluid',0,bpositionl*xx,[0 0 f2der],le-6,... bpositionl*inter,le-6*bpositionl*inter,bpositionl*xx); x0=xx; zO=Al-(yO(2)~2-yO(l)*yO(3)); end yO(1)=y0(1)*bpositionl/Rel; yO(2)=yO(2)*bpositionl~2/Rel; yO(3)=yO(3)*bpositionl~3/Rel; zO=zO*bpositionl~4/Rel; Totfluid_cal=[xO yO z O ] ; APPENDIX D. CODE FOR THE SIMILARITY Tot_Totfluidcal= [Tot_Totfluidcal; record_xx=[record_xx SOLUTION Totfluid_cal]; xx]; y.y.y.y.yx/x/.yx/x/.yx/.yx/.°/x/. end of s o i v e t o a n y P y,y.y.yx/.yx/x/x/.y.yx/.yx/.y.y. osition.m solvetoanyposition_p.m function /. t o 0 get Totfluid_calmed=solvetoanyposition_p(xx) base 7, b a s i c a l l y global solution it is at A2 b p o s i t i o n l value position i n the porous region, solver. Rel Dal porositym Tot_Totfluidcal_p if any d e s i r e d an i n i t i a l phi interface inter... record_xx_p length(Tot_Totfluidcal_p)==0 hit=0; else hit=length(Tot_Totfluidcal_p(:,1)); end if hit>=10 for i=l:hit if xx==record_xx_p(i) Totfluid_calmed=Tot_Totfluidcal_p(i,:); return end end end if xx==inter xO=xx; yO=phi; zO=A2-(yO(2)~2-y0(1)*y0(3))+y0(2); else [yO,nfun.hstart]=odesolve('eqporous',interface, xx*bpositionl,phi,le-6, xx*bpositionl-interface,le-6*(bpositionl*xx-interface) bpositionl*xx-interface); xO=xx; zO=A2-(yO(2)~2-yO(l)*yO(3))+yO(2); end yO(l)=yO(1)*sqrt(porositym/Dal)*porositym/Rel; yO(2)=yO(2)*porositym"2/Dal/Rel; yO(3)=yO(3)*porositym~2.5/Rel/Dal~l.5; zO=zO*porositym~3/Rel/Dal~2; T o t f l u i d _ c a l m e d = [xO yO z O ] ; Tot_Totfluidcal_p=[Tot_Totfluidcal_p; record_xx_p=[record_xx_p y.y.yx/x/.yx/x/,y.yx/x/.°/x/. e n d o f s o i v e t o a n y, e q s t a b i l i t y .m Totfluid_calmed]; xx]; y P osition_p. m y.y„y.yx/.y.y.y.y.y.y.y.y.°/x/x/.y,y. APPENDIX function D. CODE FOR THE SIMILARITY 124 SOLUTION z=eqstability(x,g) 7, ' e q s t a b i l i t y . m ' d e f i n e s t h e s t a b i l i t y e q u a t i o n i n f l u i d g l o b a l B l _ r lambda_r R e l T o t f l u i d _ c a l z(l)=g(2); z(2)=g(3); region, z(3)=Bl_r+Rel*(g(3)*Totfluid_cal(2)-2*Totfluid_cal(3)*g(2)+Totfluid_cal(4)*... g(l))+Rel*(lambda_r)*g(2); %7.y.%y.y.yo%y.%y.7oy////o7.7./o/./oyo/o7o/.%r/.7./„°/o% 0 oo o o o end o f e q s t a b i i i t y . m y//////.y//.y;///;///;///////x///././.7./.7./ , , , , i 7. e q s t a b i l i t y _ p .m function z=eqstability_p(x,g) 7. ' e q s t a b i l i t y . m ' d e f i n e s the s t a b i l i t y equation i n f l u i d g l o b a l lambda_r R e l D a l T o t f l u i d _ c a l m e d z(l)=g(2); z(2)=g(3); region. B2 p o r o s i t y m z(3)=B2+Rel/porositym*(g(3)*Totfluid_calmed(2)-2*Totfluid_calmed(3)*g(2)+... Totfluid_calmed(4)*g(l))+Rel*(lambda_r)*g(2)+porositym/Dal*g(2); 7//////////X///X/.7//////.7//X/X///.7.7////X/////.7. end o f e q s t a b i l i t y . p . m 7X/X/.7X/X/X/X/X/X/X/X/. Appendix E CODE FOR T H E ONE DOMAIN P R O B L E M E.l Introduction Here we present the detailed computer algorithm and code written in FORTRAN 90 for the one-domain approach outlined in Chapter 5. E.l.l Formulation of the Problem Equation The equations of dimensionless form is as follows: dw 1 dw x 1 dw x dP 1 dw dw Be dt e dx e dy dy Re dx dy DaRe x 2 2 x 2 x 2 where Da = K/H , Re — pU H/ p, and w and w denote the x and y velocity components 2 w x y respectively. The reference parameters are: the reference length is the total channel height "H", reference velocity: "U ", reference time: t = H/U and pressure P is nondimensionw lized with "p/U ". 2 w Porosity e = 1 in fluid region and e < 1 in porous medium region. B is a binary constant equals 0 for fluid region and equals 1 for porous region. 125 APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 126 To solve the equations for steady, incompressible flow using artificial compressibility method, the pressure term in continuity equation is added an artificial compressibility (3, so that the equation becomes: dP dw dw. ~dt ' pedx x dw 1 dw 1 ~dt e dx e dw.y 1 d(W Wy) 2 x dt x ' e dx = + (3edy (E.4) 0 dP d(w w ) x y dy 1 dw _ 2 y = ' e dy —e dw Be dy DaRe Be 2 x dx Re dx dP 1 ,d w, d w. dy Re dx dy' 2 2 2 2 l 2 ~Da~Re £ Wy (E.5) (E.6) This equation set can be written into the form of dU dF dG dt dx dy (E.7) = S where ( P \ U = v J \ w \ ( i£a ef3 F p_ i ^ WxVJy e ( t " J _ dw dx x Re r 1 dw v Re dx ^JL 1 dw e Re dy >K) £ Re dy eP G = x _ S = Da Re i W _ x \ Da% y ) W Equation (E.7) in control volume finite form is: Fi-l/2,j Gi j+i/2 — Gij-i/2= 5, Ay t Ax (E. APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 127 A fully discrete form using implicit Euler time advance scheme yields rrn+l i,j jrn XTT + n ~ i,j U pn+l l i,J U 0U ~T±T = At r nn+l ^ i j + 1/2 __ pn+l i+l/2,j i-l/2j r = _ /~<n+l ^i.j-1/2 Ax~~ /-p x Q Ay + Q - 5 ( E - It is reasonable to assume that F _\\_ . is a function of UJ ? and UJ\^}_, as well for G } n is a function of L / ^ 1 and 17™+^ and +1 1 9 ) n+ l n is a function of U^ . 1 The quantities at time step n + 1 can be approximated with quantities at older time step n by Taylor series expansion, from which we know that for example, r pn+l _ i+l/2,j ~~ p/jjn+l \ r = U Tjn+l i,j ' F{U^ U \ i+l,j> + 5 U ^ \ U ? + h 3 + 5 U ^ dF n 8F n Hence, h J _ ~ST - 0 U / 1 - ^ P { n , F ^ U r i+l/2,j S - m — + 6 T T n + , 1 ^ u U 8F l T? n U i 8G i 8G . , U n n + 1 1 T 7 T n + i _ ^ } i-l/2,_ , n + l , 8G n n ac . n n + ^^ ,7'--^ X N i-l/2,j l i + l / 2 j 8 F N i _L r ^ + 1 / 0 1 / 9 ; (E.11) v Therefore At Ax { i dUij8 F Ax dUij T n Ax ° U i + 1 * i PlC U K j Ay c T T n+l M _ r 1 h:J n ^ U + W n i+l/2,j K1 8U.J K n F i,j-l/2 Ay 8U , 1 8C 8U 8 J^ Ax dUi-.j n 1 Ay n r i-l/2,j Ay 8 U W ° C ^i+l/2,3 Ax U ^ + 1 C n n ^i,j-l/2 Ay l J v Equation (E.12) could be written as (7 + AtB + AtB )5U^ 1 x y +ALW™.^ +AtC^C/^- + A t ^ l ^ + = Al(FI)l AtC 5ir?£_ v (E.13) APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 128 where A — x B 9F7_ Ax dUi-ij 1/2ij 1 = x 1 C = x Ax dU j 1 Ax ^/2» dUi+ij Ay = B.„ = 4 dR i-l/2,3 ldS? d 2 da Ax dUij it 9F 1 Ay ^ - i / 1 *%+l/2 2 dUij-x Ay i at/ij 5G: Ay 1 ^ 2 at/, <9[/„ Ay a^J+a C — v r FI = ^+l/2,i G. i-l/2j fin U Ax (E.14) t J - l / 2 Ay Since / u ' x i + l / 2 , j x i + Pi+l/2,i _1_ l/2,j >3 >x J+1/2J V y 2 1 J+1/2J w dwx Re <9x ^V ,.? 2 dlVy Re 9x i(PjJ + P- j + l ^ - X U j ; 1 Re Wx i,j+U>x i + l j Wy 2 Ax 1 y i+i,j- v i,j Re Ax w i j + W y 2e w 2e/3 G M + l/2 i j + W i i,.7 + 1 \ *u » , j + ^ " ; y y y 2 7,7+1 A 1 2 J 2 s,. +l 2e 1 Re Wx i . j + l Ay e(Pi,,-+Pij +i ) I £" r -Be 2 Re Ay APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 129 2ej3 GV 7 - 1 / 2 1 - 2e -__\ 2 1 , i ( / ' • , - P,.,. e w x W j j - W y x 7,j — 1 Ay Re :) 1 2 w y i , j ~ y w —l Ay Re J+T-Ux z-1 2e(3 1/2 1 -VJ Re Ax i-1 X / 2 j 2e Re Ax V Using the method of approximate factorization, Solving Eq. (E.13) equals to solve the equation shown as following: {I + AtA E„ x + AtB lfi + AtC E )(I x x + AtA E^ lfi + y AtB n+l + y AtC E )5U% y 0il (E.15) where Eu^j is subscript operator, This leads to 0 A<r — A Ax ReAx i f e y + ^ - i j ) 0 1 Ax 0 0 0 0 1 Ax ReAx (W - x W_ i+l,j) - 2DalTe Ax 0 4e ReAx 0 i_ A 0 2e/3 0 2<r £ 1 Ay 2E 1 1 4 A, ReAx \y w ij "f" W y i + 1 J7 4E j) 0 0 0 2 j7 4e 0 (w x i,j W 1 — x ReAx 1 2e/3 (Wy i j - l + W y ij) + 1 ReAj/ 1 4E W 2e y <w x —l W *>i+i "b y w X I J ; ij-i/ ReAj/ 2DaRe Ax APPENDIX E. THE ALGORITHM o 1 Ay By = AND CODE FOR ONE DOMAIN PROBLEM 130 o o - ReAy Be 4e ( V >J~ W i V 1 W i J + ) 1 0 0 1 + W 4£ 2DaRe x 2 2e ( W y M - 1 ReAy Be 2DaRe Ay i 0 0 2e/3 Ay' § E.1.2 o 2e ( V M "F ! / W W + 1 ReAy Setup of the Boundary Conditions The boundary conditions include implicit and explicit conditions, where the explicit ones are used to determine the ghost cell values and the implicit ones are for the update of the cell values. The boundary conditions include upper impermeable wall and lower porous wall; determination of inflow and outflow conditions. • Upper Wall The upper wall is a solid impermeable wall, therefore, the boundary conditions are: W x = W y — 0 and ^ = 0. The values in the ghost cells for our coordinate could be set as: (E.16) P(U) P(i,0) = w (i,0) = -w (i,l) (E.17) w (i,Q) = -w (i,l) (E.18) x y x y The corresponding implicit boundary condition is: ( 1 \ I 5PP+ 1 \ ( 1 - i 1 V 1 \ V -1 This condition is used to construct the coefficient matrix for ghost cells at the upper wall. APPENDIX • Bottom E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 131 Wall The bottom wall is a porous wall, which implies a uniform suction in y direction and zero velocity at x direction. The explicit pressure boundary condition is derived from the y momentum equation, i.e., 1 d(w w ) 1 dw. e e dy x y dx 1 ,d w , dP e 2 dy d w. Be 2 1 + Re^ dx 2 (E.19) DaRe Wy dy Since 4- = 0 because w = 0 throughout x direction, from continuity we thus have x also = 0. Therefore, the above equation reduces to: dH where w.y wall (E.20) DaRe dy Wy = 1. This means that the explicit conditions become: P(i, grialy + 1) = +P(i, gridy) — + l) = -w (i, w (i, gridy + 1) = -w (i. w (i,gridy x y x y Ay (E.21) ReDa (E.22) gridy) gridy) + 2w yw a (E.23) U The corresponding implicit conditions are: \ / <5P , n+1 i,gndy \ \ ( J \ yi,gridy 0UJ n+1 i,gridy+l - 1 x i,gridy V SP x i,gridy+l - 1 ) UUJ y i,gridy+l / ULU which is the same with the implicit upper wall condition. • Inflow In this case, we specify the inflow velocity w velocity w xi n x i n and w y i n , where w y i n = 0. The can be specified to different values, it is reasonable to assume uniform inflow of velocity U = w i /U . in x n suc A pressure gradient can be assumed at the inflow position, which finally gives: P(0,j) = +P(l,j)-Ax dP dx , (E.24) ir w {0,j) x = -w (l,j) x + 2Ui, (E.25) (E.26) APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 132 The corresponding implicit boundary conditions are: / 1 \ r pn+l \ SKt ( 1 / >n+l 1 n+l -1 1 V • ,n+l -1 V \ / Outflow At channel outlet we assume the flow will be fully-developed, i.e., the pressure, velocity components are not the functions of position x anymore; Following the above, we have: We could assume the pressure P , out and both the velocity gradients along the x and y directions should be zeros. These give: P(gridx = w (gridx + 1, j) = w (gridx,j) (E.28) w (gridx + l,j) = w {gridx,j) (E.29) x y -P{gmdx,j) + (E.27) + l,j) 2P out x y This corresponds to the implicit conditions as following: \ gridx,j \ x gridx,j \ ( Xpn+1 gridx + l,j +1 x gridx+l,j ULU V E.1.3 \ ^ y gridx,j ) W + +1 \ ^ y gridx+l,j W + ) Linearization and Solving the Linear Equations The approximate factorization method is used to reduce the problem to solve a system of a block tri-diagonal matrix for constant' "i" and "j" values respectively. For example, for constant " i " , we have: APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 133 -1 1 1 1 0 1 1 Bj 0 RHSi Cj RHS A, 2 -1 1 0 0 1 1 1 1 where the right hand side RHSi is the sum of source term and flux integral at time step n. In this solution, to dampen the fluctuation of the pressure terms, an artificial total pressure gradient is also added on the right hand side .expression. This treatment is found necessary to have converged solution. The above block tri-diagonal matrix is solved using standard Thomas algorithm. The control of the source term and the porosity values in different flow domains are operated through turning on/off the source term values of cell number falling into the region. The code for the One-domain problem is is written in Fortran90. We presented it here as onedomaincode. f 90. The main program, is listed first followed by the subroutines. APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN One-Domain S o l u t i o n Fluid flow i n a channel a porous This program partially The codes filled algorithm preceding program filled Deng J u l y 6, 2004 Navier-Stokes Equation for w i t h a porous and the descriptions i n the fluid medium by a o n e - d o m a i n parameters with in this code flow in a channel approach. c o u l d be found in the chapter. onedomaincode implicit integer none i,j,k,m,gridx,gridy,maxiter,iterno parameter double the updated for partially medium by Chuntao last PROBLEM (gridx=480,gridy=60,maxiter=100) precision p(0:gridx+1,0:gridy+1), & u(0:gridx+l,0:gridy+l),v(0:gridx+1,0:gridy+1), & x(0:gridx+l),y(0:gridy+1),divp(3,gridx,gridy), & errors(3), deltaU(3,gridx,gridy), & Iden(3,3), flcom(3,gridx,gridy), & RHSxx(3,0:gridy+1),RHSyy(3,0:gridx+l), & LHSxx(3,3,3,0:gridy+l),LHSyy(3,3,3,0:gridx+l), & xdis, & ydis.deltax,deltay,deltat,conv_cri,errmax, ORfactor, LHSyA pO,uO,vO,pi,LHSx(3,3,3,2),LHSy (3,3,3,2), (3,3,gridx,gridy),LHSyB (3,3,gridx,gridy), LHSyC(3,3,gridx,gridy), LHSxA(3,3,gridx,gridy), LHSxB(3,3,gridx,gridy),LHSxC(3,3,gridx,gridy), sij(3,gridx, gridy), Da, lambda,Re parameter(Re=5.dO,Da=l.e-3,lambda=l.dO) parameter(p0=l.d0, u0=l.dO,v0=l.dO, parameter(xdis=8.OdO, pi=3.1415926535897932d0) ydis=l.OdO,deltat=0.02d0,conv_cri=l.d-4) parameter(ORfactor=l.2d0) gridx : grid number in x direction number in y direction gridy : grid maxiter : a l l o w e d maximum i t e r a t i o n p : pressure u,v : v e l o c i t y i n x and y d i r e c t i o n times respectively & & & & 134 APPENDIX E. THE ALGORITHM f lcom flux divp source term deltat length of x.y x deltax sizes AND CODE integral and y FOR ONE DOMAIN PROBLEM 135 computed added time to suppress pressure oscillation step coordinates of cells i n x and y d i r e c t i o n respectively deltay LHSxx,LHSyy: lefthand LHSxx.LHSyy: righthand deltaU changes Iden 3 by 3 i d e n t i t y side matrix side in of vector Block of over-relaxation factor L2 norm of component solution tolerrance xdis,ydis domain ! define identity Tridiagonal Algebric matrix errors conv_cri Block solution ORfactor errmax Tridiagonal Algebric each update length w and h e i g h t 1 matrix Iden(l:3,l:3)=0.0d0 Iden(l,l) =1.0d0 Iden(2,2) =1.0d0 Iden(3,3) =1.0d0 ! meshing the call ! domain meshing(xdis,ydis,gridx,gridy,x,y,deltax,deltay) I n i t i a l i z e the do i=0, do quantity i n each cell (p,u,v) with Sec.1.1 gridx+1 j=0,gridy+l p(i,j)=pO*dcos(pi*x(i))*dcos(pi*y(j)) u(i,j)=u0*dsin(pi*x(i))*dsin(2.d0*pi*y(j)) v(i,j)=v0*dsin(2.dO*pi*x(i))*dsin(pi*y(j)) enddo enddo ! set call values i n ghost cells by boundary conditions BC_cavity(p,u,v,gridx,gridy,deltax, deltay) open(unit=ll,file='TestO.dat',status='unknown') open(unit=10,file='Error.dat',status='unknown') errmax=l.OdO iterno=0 do w h i l e ((errmax.ge.conv_cri) .and. (iterno.It.maxiter)) data Equation Equation APPENDIX ! ! ! E. compute THE ALGORITHM Jacobian CODE FOR ONE DOMAIN coeff_matrixx(LHSxA,LHSxB,LHSxC,gridx,gridy,u,v,deltax) call coeff_matrixy(LHSyA,LHSyB,LHSyC,gridx,gridy,u,v,deltay) integral call flux_comput(p,u,v,gridx, call source(u,v,gridx, compute deltay,flcom) addedterm(p,gridx,gridy,deltax,deltay,divp) constant do gridy,deltax, gridy,sij) A*deltax~2*divP call ! Along PROBLEM coefficient call compute f l u x ! At AND i , solve for deltaU~ i=l,gridx boundary call boundary_LHSx(LHSx) LHSxx(1:3,1:3,1:3,0)=LHSx(1:3,1:3,1:3,1) LHSxx(1:3,1:3,1:3,gridy+1)=LHSx(1:3,1:3,1:3,2) RHSxx(l:3,0)=0.d0 RHSxx(l:3,gridy+l)=0.d0 ! At interior do j=l,gridy do k=l,3 RHSxx(k,j)=deltat*(flcom(k,i,j)+divp(k,i,j)+sij(k,i,j)) do !Ax=-LHSyA/deltay; m=l,3 Bx=-LHSyB/deltay; Cx=LHsyC/deltay LHSxx(k,m,1,j)=-deltat*LHSyA(k,m,i,j)/deltay LHSxx(k,m,2,j)=-deltat*LHSyB(k,m,i,j)/deltay+Iden(k,m) LHSxx(k,m,3,j)= end end end call ! prepare do deltat*LHSyC(k,m,i,j)/deltay do do do SolveBlockTri(LHSxx, to solve deltaU j=l,gridy do k=l,3 RHSxx, gridy+2, gridy+2) 136 APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM deltaU(k,i,j)=RHSxx(k,j) end end end ! along do do do constant do j j=l,gridy call boundary_LHSy(LHSy) LHSyy(1:3,1:3,1:3,0)=LHSy(1:3,1:3,1:3,1) LHSyy(1:3,1:3,1:3,gridx+1)=LHSy(1:3,1:3,1:3,2) RHSyy(l:3,0)=0.dO RHSyy(1:3,gridx+1)=0.d0 do i=l,gridx do k=l,3 RHSyy(k,i)=deltaU(k,i,j) do m=l,3 LHSyy(k,m,1,i)=-deltat*LHSxA(k,m,i,j)/deltax LHSyy(k,m,2,i)=-deltat*LHSxB(k,m,i,j)/deltax+Iden(k,m) LHSyy(k,m,3,i)= end end end deltat*LHSxC(k,m,i,j)/deltax do do do call do SolveBlockTri(LHSyy, RHSyy, gridx+2, gridx+2) i=l,gridx p(i,j)=p(i,j)+RHSyy(1,i)*0Rfactor u(i,j)=u(i,j)+RHSyy(2,i)*0Rfactor v(i,j)=v(i,j)+RHSyy(3,i)*0Rfactor do k=l,3 deltaU(k,i,j)=RHSyy(k,i) end end end do do do ! boundary call fit BC_cavity(p,u,v,gridx,gridy,deltax,deltay) errmax=0. call errcomput(deltaU,gridx,gridy,errors, errmax) iterno=iterno+l ! output L2 norm v e r s u s write(11,105) iteration iterno, time steps errors(1),errors(2),errors(3) enddo ! call Usymmetry (u, gridx, gridy,Usym) 137 APPENDIX ! E. THE ALGORITHM FOR ONE DOMAIN PROBLEM output ! ! open(unit=10,file='ErrorO.dat',status='unknown') do i=l,gridx do j=l,gridy write(10,100) end end ! AND CODE p(i,j),u(i,j),v(i,j) do do write(10,120)(x(i),i=l,gridx) ! do j=l,gridy ! write(10,120)(p(i,j),i=l,gridx) ! end do ! write(10,120)(y(j),j=l,gridy) ! write(10,120)(Usym(j),j=l,gridy) close(10) close(ll) 105 format(i6,3el8.8) 100 format(3el8.8) 120 format(lell.4) stop end I ! mesh Main Program ydisO, gridxO.gridyO, generation subroutine implicit integer double END o f meshing(xdis0, i , j , gridxO,gridyO precision xO(0:gridx0+l),y0(0:gridy0+l),xdis0,ydis0,deltax0,deltay0 deltax0=xdis0/dfloat(gridxO) deltayO=ydisO/dfloat(gridyO) do i=0,gridxO+l xOd) = (1.0d0*i-0.5d0)*deltax0 enddo do j=0,gridyO+l yO(j) enddo x0,y0,deltax0,deltay0) none = (1.0d0*j-0.5d0)*deltay0 138 APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 139 return end I END o f ! Compute coefficients subroutine implicit integer double for subroutine meshing Ay.By.Cy coeff_matrixy(LHSyA,LHSyB,LHSyC, gridx,gridy,u,v,deltay) none i,j,gridx, gridy precision beta, LHSyA (3,3,gridx,gridy),u(0:gridx+1,0:gridy+1), v(0:gridx+1,0:gridy+1), LHSyC parameter (beta=l.d0, Re, (3,3,gridx,gridy), Re=5.d0, Da, & deltay,LHSyB(3,3,gridx,gridy), epsl, lambda, h, tern, B Da=l.e-3) tem=dfloat(gridy) h=int(0.9d0*tem) do i do = 1, gridx j=l, if gridy (j>=h) then epsl=.9d0 lambda=l.dO/epsl B=l.d0 else epsl=l.d0 l a m b d a = l .dO B=0.d0 endif LHSyA(l,l,i,j) = O.dO LHSyA(l,2,i,j) = O.dO LHSyA(1,3,i,j) = 1.d0/2.dO/beta/epsl LHSyA(2,l,i,j) = O.dO LHSyA(2,2,i,j) = 1.dO/4.dO*(v(i,j-l)+v(i,j))/epsl LHSyA(2,3,i,j) = LHSyA(3,l,i,j) = LHSyA(3,2,i,j) = O.dO LHSyA(3,3,i,j) = 1.dO/2.dO*(v(i,j)+v(i,j-1))/epsl & +1.dO/Re/deltay*epsl*lambda 1.dO/4.dO*(u(i,j)+u(i,j-1))/epsl 1.dO/2.d0*epsl & +1.dO/Re/deltay*epsl*lambda LHSyB(l,l,i,j) = O.dO LHSyB(1,2,i,j) = O.dO LHSyB(1,3,i,j) = O.dO LHSyB(2,l,i,j) = O.dO LHSyB(2,2,i,j) = 1.dO/4.dO*(v(i,j-1)-v(i,j+1))/epsl -2.dO/Re/deltay*epsl*lambda -1.dO/2.dO*(B*(epsl/Da/Re))*deltay LHSyB(2,3,i,j) = 1.dO/4.dO*(u(i,j-1)-u(i,j+1))/epsl & & & APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN LHSyB(3,1,i,j) = LHSyB(3,2,i,j) = O.dO LHSyB(3,3,i,j) = 1. d O / 2 . d O * ( v ( i , j - 1 ) - v ( i , j + l ) ) / e p s l PROBLEM 0.dO & -2.dO/Re/deltay*epsl*lambda -1.dO/2.dO*(B*(epsl/Da/Re))*deltay LHSyC(l,l,i,j) = O.dO LHSyC(l,2,i,j) = O.dO LHSyC(l,3,i,j) = 1.dO/2.dO/beta/epsl LHSyC(2,l,i,j) = O.dO LHSyC(2,2,i,j) = LHSyC(2,3,i,j) = LHSyC(3,l,i,j) = 1.dO/4.dO*(v(i,j)+v(i,j+l))/epsl -1.dO/Re/deltay*epsl/lambda 1.dO/4.dO*(u(i,j)+u(i,j+l))/epsl 1.dO/2.dO*epsl LHSyC(3,2,i,j) = O.dO LHSyC(3,3,i,j) = 1.dO/2.dO*(v(i,j)+v(i,j+l))/epsl & -1.dO/Re/deltay*epsl*lambda end end do do return end ! END o f ! Compute coefficients subroutine implicit integer double for subroutine coeff_matrixy Ax,Bx,Cx coeff_matrixx(LHSxA,LHSxB,LHSxC,gridx,gridy,u,v,deltax) none i,j,gridx, precision gridy beta, u(0:gridx+1,0:gridy+1), v(0:gridx+l,0:gridy+l), LHSxA LHSxC parameter (beta=l.dO, tem=dfloat(gridy) h=int(0.9d0*tem) do i = 1, gridx do j = l , if - gridy (j>=h) then epsl=.9d0 B=l.dO lambda=l.dO/epsl else B=0.dO epsl=l.dO lambda=l.dO end i f Re, & deltax,B, (3,3,gridx,gridy),LHSxB (3,3,gridx,gridy), (3,3,gridx,gridy), Re=5.d0, Da=l.e-3) Da, epsl, lambda,h,tern & & 140 APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN LHSxA(l,l,i,j) O.dO LHSxA(l,2,i,j) 1 LHSxA(l,3,i,j) 0 . dO LHSxA(2,l,i,j) 1. d 0 / 2 . d 0 * e p s l LHSxA(2,2,i,j) I.d0/2.d0*(u(i,j)+u(i-l,j))/epsl PROBLEM .d0/2.d0/beta/epsl +1.dO/Re/deltax*epsl*lambda LHSxA(2,3,i,j) O.dO LHSxA(3,l,i, j) 0 . dO LHSxA(3,2,i, j) 1. LHSxA(3,3,i,j) I.d0/4.d0*(u(i-l,j)+u(i,j))/epsl d0/4.d0*(v(i-l,j)+v(i,j))/epsl +1.dO/Re/deltax*epsl*lambda LHSxB(l,l,i,j) O.dO LHSxB(l,2,i,j) O.dO LHSxB(l,3,i,j) O.dO LHSxB(2,l,i,j) 0 . dO LHSxB(2,2,i,j) 1. dO/2.dO*(u(i-1,j)-u(i+l,j))/epsl -2.dO/Re/deltax*epsl*lambda -1.dO/2.dO*(B*(epsl/Da/Re))*deltax LHSxB(2,3,i,j) O.dO LHSxB(3,l,i,j) 0 . dO LHSxB(3,2,i,j) 1. d O / 4 . d O * ( v ( i - 1 , j ) - v ( i + 1 , j ) ) / e p s l LHSxB(3,3,i,j) 1.dO/4.dO*(u(i-l,j)-u(i+l,j))/epsl -2.dO/Re/deltax*epsl*lambda -1.dO/2.dO*(B*(epsl/Da/Re))*deltax 0 . dO LHSxC(l,l,i,j) LHSxC(l,2,i,j) 1. LHSxC(l,3,i,j) 0 . dO LHSxC(2,l,i,j) 1. d 0 / 2 . d 0 * e p s l LHSxC(2,2,i,j) d0/2.d0/beta/epsl I.d0/2.d0*(u(i,j)+u(i+l,j))/epsl -1.dO/Re/deltax*epsl*lambda LHSxC(2,3,i,j) O.dO LHSxC(3,l,i,j) 0 . dO LHSxC(3,2,i,j) 1. LHSxC(3,3,i, j) 1.dO/4.dO*(u(i,j)+u(i+1,j))/epsl d0/4.d0*(v(i,j)+v(i+1,j))/epsl -1.dO/Re/deltax*epsl*lambda enddo end do return end END o f ! compute flux subroutine subroutine coeff_matrixx integral flux_comput (pp,uu,vv,gridxO, gridyO,deltax, deltay,flcomO) 141 APPENDIX E. implicit integer double THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM 142 none k, i , j , gridxO, gridyO, ip, im, jp, jm precision flcom0(3,gridxO, gridyO),pp(0:gridxO+1,0:gridyO+1), uu(0:gridxO+l,0:gridyO+l),vv(0:gridxO+1,0:gridyO+1), Re, parameter beta, (Re=5.d0, deltax, beta=l.dO, deltay.Da, epsl, lambda, & & h,tem Da=l.e-3) tem=dfloat(gridyO) h=int(0.9d0*tem) do k=l,3 do i=l, gridxO ip=i+l; im=i-l do j = l , gridyO jp=j+l; if jm=j-l (j>=h) then epsl=.9d0 lambda=l.dO/epsl else epsl=l.dO lambda=l.dO endif if (k==l) then flcomO(k,i,j)=-(uu(ip,j)-uu(im,j))/2.dO/beta/deltax/epsl -(vv(i,j elseif & p)-vv(i,jm))/beta/deltay/2.dO/epsl (k==2) then flcomO(k,i,j)=-((uu(i+l,j)+uu(i,j))**2/4.dO/epsl & -(uu(i,j)+uu(i-l,j))**2/4.d0/epsl & +(pp(i+1,j)-pp(i-1,j))/2.dO*epsl & -(uu(i+l,j)-2.0d0*uu(i,j)+uu(i-l,j)) & /(Re*deltax)*epsl*lambda)/deltax & -((uu(i,j+l)+uu(i,j))*(vv(i,j+l) & +vv(i,j))/4.d0/epsl & -(uu(i,j-l)+uu(i,j))*(vv(i,j-1) & +vv(i,j))/4.d0/epsl & -(uu(i,j+l)-2.d0*uu(i,j)+uu(i,j-1)) & /(Re*deltay)*epsl*lambda)/deltay else fIcomO(k,i,j)=-((uu(ip,j)+uu(i,j))/2.dO*(vv(ip,j) & +vv(i,j))/2.d0/epsl-(uu(im,j)+uu(i,j))/2.d0*(vv(im,j) & +vv(i,j))/2.d0/epsl-(vv(ip,j)-2.dO*vv(i,j)+vv(im,j)) & /(Re*deltax)*epsl*lambda)/deltax & -((vv(i,jp)+vv(i,j))**2/4.d0/epsl & -(vv(i,j)+vv(i,jm))**2/4.d0/epsl & +(pp(i,jp)-pp(i,jm))/2.d0*epsl & APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN PROBLEM -(vv(i,jp)-2.dO*vv(i,j)+vv(i,jm)) 143 & /(Re*deltay)*epsl*lambda)/deltay endif enddo enddo enddo return end END o f subroutine implicit integer double subroutine flux_comput source(uu,vv,gridxO, gridyO,sij) none k, i , j , gridxO, precision gridyO,h sij(3,gridxO, gridyO), & uu(0:gridxO+l,0:gridyO+l),vv(0:gridxO+1,0:gridyO+1), Re, parameter B, epsl, (Re=5.d0, Da.tem Da=l.e-3) tem=dfloat(gridyO) h=int(0.9d0*tem) do k=l,3 do i = l , gridxO do j = l , if gridyO (j>=h) then B=l.dO epsl=.9d0 else B=0.dO epsl=l.dO endif if (k==l) then sij(k,i,j)=0.dO elseif (k==2) then sij(k,i,j)=-B*epsl/Da/Re*uu(i,j) else sij(k,i,j)=-B*epsl/Da/Re*vv(i,j) endif enddo enddo enddo return end END o f subroutine source & APPENDIX E. THE ALGORITHM subroutine implicit integer do FOR ONE DOMAIN PROBLEM 144 errcomput(deltaU,gridx,gridy,errors,errmax) none m,i,j, double AND CODE gridx, gridy precision deltaU(3,gridx,gridy), errors(3), errmax m=l,3 errors(m)=0. do i=l,gridx do j=l,gridy errors(m) end end = errors(m)+deltaU(m,i,j)**2.dO do do errors(m)=dsqrt(errors(m)/dfloat(gridx)/dfloat(gridy)) end do errmax=errors(1) ! errmax=errors(2) do m=2,3 if (errors(m) .gt. errmax) then errmax=errors(m) end end if do return end I ! END o f set ghost cell subroutine implicit integer subroutine errcomput values BC_cavity(p,u,v,gridx,gridy,deltax, deltay) none i , j , gridx,gridy double p r e c i s i o n p ( 0 : g r i d x + l , 0 : g r i d y + l ) , u(0:gridx+1,0:gridy+1), v(0:gridx+l,0:gridy+l), Vw, parameter(Vw=l.dO, Rat, k i , Re, Re=5.d0, deltax, Da, tem, deltay, & Uin, & h Da=l.e-3) Uin=21.d0 Rat=Uin/Vw ki=100.d0 tem=dfloat(gridy) h=int(0.9d0*tem) do i=l,gridx p(i,gridy+l)= ! p(i,0)= p(i,gridy)-deltay/Re/Da p(i,l) u(i,gridy+1)=-u(i,gridy) v(i,gridy+1)=-v(i,gridy)+2.dO ! at bottom w a l l , suction APPENDIX E. THE ALGORITHM p(i,0)= p(i,l) AND CODE ! at upper FOR ONE DOMAIN PROBLEM 145 wall u(i,0)=-u(i,l) v(i,0)=-v(i,l) end do do j=l,gridy p(0,j)=p(l,j)-deltax*ki ! p(0,j) ! i f = (j.eq.l) ! if left pd.j) p(0,j)=-p(1,j) (j>=h) ! ! at Ifixed pressure at one point then u(0,j)=-u(l,j) ! else u(0,j)=-u(l,j)+Rat*2.dO ! endif ! u(0,j)=-u(l,j) v(0,j)=-v(l,j) p(gridx+1,j)=-p(gridx,j)+.08d0*2.d0 ! ! at right wall, fully developed p(gridx+1,j)=-p(gridx,j) ! u(gridx+l,j)= u(gridx,j) v(gridx+1,j)= v(gridx,j) p(gridx+1,j)= p(gridx,j) ! u(gridx+l,j)=-u(gridx,j) ! v(gridx+1,j)=-v(gridx,j) end ! at right wall, fully developed do return end I END o f ! Suppresion for subroutine implicit integer pressure subroutine ossilation, term routine addedterm(p,gridx,gridy,deltax,deltay,divp) none precision divp(3,gridx,gridy),p(0:gridx+1, 0:gridy+l), deltax,deltay, A=0.5d0 divp(1:3,1:gridx,1:gridy)=0. i=l,gridx do add a source i,j,gridx,gridy double do BC_cavity j=l,gridy A & APPENDIX E. THE ALGORITHM divp(l,i,j)= AND CODE FOR ONE DOMAIN end & -2.d0*p(i,j)+p(i-l,j))/deltax/deltax & +I.d0/2.d0*(p(i,j+l) & do do return end I END o f ! Sanity check subroutine implicit for subroutine addedterm symmetry Usymmetry (u, gridx, gridy, u_sym) none integer j , m, g r i d x , gridy double p r e c i s i o n u(0:gridx+l,0:gridy+1), u_sym(gridy) m=gridx/2 do j=l,gridy u_sym(j)=(u(m,j)+u(m+1,j))*0.5d0 end do return end I __]_ ! Ax,Bx,Cx and A y , B y , C y subroutine implicit double 146 A*deltax*deltax*(l.d0/2.d0*(p(i+l,j) -2.d0*p(i,j)+p(i,j-1))/deltay/deltay) end PROBLEM at 0 f subroutine boundary boundary_LHSy(LHSy) none p r e c i s i o n LHSy (3,3,3,2) LHSy(l,1,1,1) = 0 dO LHSy(l,2,l,l) LHSy(1,3,1,1) = 0 dO = 0 dO LHSy(2,1,1,1) = 0 dO LHSy(2,2,l,l) LHSy(2,3,l,l) = 0 dO = 0 dO LHSy(3,1,1,1) = 0 dO LHSy(3,2,1,1) = 0 dO LHSy(3,3,l,l) = 0 dO LHSy(l,1,2,1) = 1 dO LHSy(l,2,2,l) = 0 dO LHSy(l,3,2,l) = 0 .dO LHSy(2,1,2,1) = 0 .dO LHSy(2,2,2,l) = 1 .dO Usymmetry APPENDIX ALGORITHM L H S y ( 2 , 3 , 2 , 1) = 0 . dO L H S y ( 3 , 1, 2 , 1) = 0 . dO L H S y ( 3 , 2 , 2 , 1) = 0 . dO L H S y ( 3 , 3 , 2 , 1) = 1. dO L H S y ( l , 1, 3 , 1) =- -1 dO L H S y C l , 2 , 3 , 1) = 0 dO L H S y ( 1 , 3 3 , 1) = 0 dO L H S y ( 2 , 1, 3 , 1) = 0 dO L H S y ( 2 , 2 3 , 1) = L H S y ( 2 , 3 3 , 1) = 0 dO L H S y ( 3 , 1 3 , 1) = 0 dO L H S y ( 3 , 2 3 , 1) = 0 dO L H S y ( 3 , 3 3 , 1) = 1 dO L H S y C l , 1 1, 2 ) = 1 dO 1 dO LHSyCl 2 1, 2 ) = 0 dO LHSyCl 3 1, 2 ) = 0 dO LHSy(2 1 1 2) = 0 dO LHSy(2 2 1 2) =- -1 LHSy(2 3 1 2) = 0 dO LHSy(3 1 1 2) = 0 dO dO LHSy(3 2 1 2) = 0 dO LHSy(3 3 1 2) =- -1 dO LHSyCl 1 2 2) = 1 dO LHSyCl 2 2 2) = 0 dO LHSyCl 3 2 2) = 0 dO LHSyC2 1 2 2) = 0 dO LHSy(2 2 2 2) = LHSy(2 3 2 2) = 0 dO LHSyC3 1 ,2 2) = 0 dO L H S y ( 3 2 ,2 2) = 0 dO LHSyC3 3 ,2 2) = L H S y ( 1 ,1 , 3 2) = 0 .dO L H S y C l ,2 ,3 ,2) = 0 .dO L H S y C l ,3 ,3 ,2) = 0 .dO L H S y C 2 ,1 ,3 , 2 ) = 0 .dO L H S y ( 2 ,2 ,3 ,2) = 0 .dO LHSy(2 ,3 ,3 ,2) = 0 .dO L H S y ( 3 ,1 ,3 ,2) = 0 .dO L H S y ( 3 ,2 ,3 ,2) = 0 .dO L H S y ( 3 ,3 ,3 ,2) = 0 .dO return end E. THE 1 dO 1 dO AND CODE FOR ONE DOMAIN PROBLEM 147 APPENDIX E. THE ALGORITHM END o f subroutine implicit none LHSy (3,3,3,2) L H S y ( 1, 3 , 1, 1) = 0 dO = 0 dO = 0 dO L H S y ( 2 , 1, 1, 1) = 0 dO L H S y ( 2 , 2 , 1, 1) L H S y ( 3 1, 1, 1 ) = 0 dO = 0 dO = 0 dO L H S y ( 3 2 , 1, 1) = 0 dO L H S y ( 3 3 , 1, 1) 0 dO L H S y ( 1 1, 2 , 1) = 1 dO L H S y ( 1 2 , 2 , 1) = 0 dO L H S y ( 1 3 , 2 , 1) = 0 dO LHSy ( 2 L H S y ( 1, 1, 1, 1) L H S y ( 1, 2 , 1, 1) L H S y ( 2 , 3 , 1, 1) LHSy ( 2 2 2 1) LHSy ( 2 3 2 1) = 0 dO = 1 dO = 0 dO L H S y (,3 1 2 1) = 0 dO L H S y ;3 2 2 1) 0 dO L H S y :3 3 2 1) = 1 dO LHSy : i 1 3 1) LHSy : i 2 3 1) =- -1 dO = 0 dO LHSy : i 1, 2 , 1) 3 3 1) = 0 dO L H S y [2 1 3 1 ) = 0 dO 2 3 1) = 1 dO 'a 3 3 1) = 0 dO L H S y :3 1 3 1) L H S y ^3 2 3 1) = 0 dO = 0 dO L H S y :3 3 3 1) = 1 dO LHSy LHSy LHSy LHSy subroutine boundary_LHSx(LHSy) double p r e c i s i o n LHSy AND CODE a a a 1 1 2) =--1 dO 2 1 2) 3 1 2) = 0 dO = 0 dO L H S y <2 1 1 2 ) L H S y (.2 , 2 1 2) L H S y (2 , 3 1 2) = 0 dO = 1 .dO = 0 .dO L H S y (3 , 1 , 1 , 2 ) = 0 .dO L H S y (3 , 2 , 1 , 2 ) = 0 .dO LHSy (3 , 3 ,1 , 2 ) 1 .dO L H S y (1 , 1 , 2 , 2 ) 1 .dO L H S y ( 1 ,2 , 2 , 2 ) = 0 .dO FOR ONE DOMAIN boundary_LHSy PROBLEM 148 APPENDIX E. THE ALGORITHM LHSy(2,1,2,2) = 0 dO = 0 dO LHSy(2,2,2,2) = 1 dO LHSy(2,3,2,2) = 0 dO LHSy(3,1,2,2) = 0 dO LHSy(3,2,2,2) = 0 dO LHSy(3,3,2,2) 1 dO LHSy(1,1,3,2) 0 dO = 0 dO LHSy(l,3,2,2) LHSy(l,2,3,2) LHSy(2,2,3,2) = 0 dO = 0 dO = 0 dO LHSy(2,3,3,2) = 0 dO LHSy(3,1,3,2) LHSy(3,2,3,2) = 0 dO = 0 dO LHSy(3,3,3,2) = 0 dO LHSy(l,3,3,2) LHSy(2,1,3,2) AND CODE FOR ONE DOMAIN PROBLEM return end ! END o f ! following subroutines subroutine integer NRows, integer double LHS(3, ! j ! ! = 1, Compute t h e Scale 3, 3, algorithm MaxSize) MaxSize), Inv(3,3), RHS(3, TempMat(3,3), TempMat2(3,3), MaxSize) inverse of the main block Invert3x3(LHS(l,l,2,j), the right-most block Mult3x3(Inv, call Copy3x3(TempMat, the right-hand call CopyVec(TempVec, Left-multiply ! subtract from ! term the ! First the by the RHS(l.j), the jth the j+lst RHS o f diagonal inverse. TempMat) inverse. TempVec) RHS(l,j)) row by t h e the by the LHS(1,1,3,j)) side MultVec(Inv, diagonal. Inv) LHS(1,1,3,j), call and TempVec(3) TVec2(3,3) NRows-1 call Scale Tomas B l o c k t r i a n g l e j precision call ! boundary_LHSx MaxSize double p r e c i s i o n 10 for S o l v e B l o c k T r i ( L H S , RHS, NRows, double p r e c i s i o n do are subroutine row. jth LHS m a n i p u l a t i o n sub-diagonal This row. involves on t h e the j+lst row super-diagonal and 149 APPENDIX E. call Add3x3(LHS(1,1,2,j+l), call Copy3x3(TempMat2, RHS call AddVec(RHS(l,j+l), call TempMat2) CopyVec(TVec2, RHS(l,j), TempVec, -l.dO, TempVec) TVec2) RHS(1,j+l)) do with forward elimination inverse of the loop last main block diagonal, = NRows call Invert3x3(LHS(l,l,2,j), Scale the right-hand call MultVec(Inv, call CopyVec(TempVec, Now d o do ! -l.dO, TempMat) manipulation MultVec(LHS(1,1,1,j+l), Compute t h e ! TempMat, PROBLEM LHS(1,1,2,j+l)) call Done ! CODE FOR ONE DOMAIN Mult3x3(LHS(l,l,l,j+l), LHS(1,1,3,j), end j AND call Now t h e 10 THE ALGORITHM 20 j the by the RHS(l,j), inverse. TempVec) RHS(l,j)) back-substitution. = NRows - Matrix-vector RHS(l,j) side Inv) 1, 1, multiply = RHS(l,j) -1 and - subtract. (LHS(l,l,3,j)*RHS(l,j+l) LHS(1,2,3,j)*RHS(2,j+l) + + - & & LHS(1,3,3,j)*RHS(3,j+l)) RHS(2,j) = RHS(2,j) - (LHS(2,1,3,j)*RHS(1,j+l) LHS ( 2 , 2 , 3 , j ) * R H S ( 2 , j + l ) + +' & & LHS(2,3,3,j)*RHS(3,j+l)) RHS(3,j) = RHS(3,j) - (LHS(3,1,3,j)*RHS(1,j+l) LHS(3,2,3,j)*RHS(2,j+l) LHS(3,3,3, j)*RHS(3, 20 end + do end double SpewMatrix(Source) precision Source(3,3) write(6,*) Source(1,1), Source(2,l), Source(3,l) write(6,*) Source(l,2), Source(2,2), Source(3,2) write(6,*) Source(l,3), Source(2,3), Source(3,3) return end subroutine SpewVector(Source) & & j+D) return subroutine + 150 APPENDIX double E. THE ALGORITHM precision write(6,*) AND CODE FOR ONE DOMAIN Source(3) Source(l), Source(2), Source(3) return end subroutine double CopyVec(Source, Target) precision Source(3), Target(l) = Source(l) Target(2) = Source(2) Target(3) = Source(3) Target(3) return end subroutine double Copy3x3(Source, Target) precision Source(3,3), Target(l.l) = Source(l.l) Target(1,2) = Source(l,2) Target(1,3) = Source(1,3) Target(2,1) = Source(2,1) Target(2,2) = Source(2,2) Target(2,3) = Source(2,3) Target(3,1) = Source(3,1) Target(3,2) = Source(3,2) Target(3,3) = Source(3,3) Target(3,3) return end subroutine double MultVec(A, Vec, ResultO) precision A(3,3), Vec(3), Result0(3) ResultO(l) = A(l,l)*Vec(l) + A(l,2)*Vec(2) Result0(2) = A(2,l)*Vec(l) + A(2,2)*Vec(2) + A(2,3)*Vec(3) + A(l,3)*Vec(3) Result0(3) = A(3,l)*Vec(l) + A(3,2)*Vec(2) + A(3,3)*Vec(3) return end subroutine Mult3x3(A, B , C) double precision A(3,3), B(3,3), C(3,3) C(l,l) = A(1,1)*B(1,1) + A(1,2)*B(2,1) + A(1,3)*B(3,1) C(l,2) = A(1,1)*B(1,2) + A(1,2)*B(2,2) + C(l,3) = A(1,1)*B(1,3) + A(1,2)*B(2,3) + A(1,3)*B(3,3) A(l,3)*B(3,2) PROBLEM 151 APPENDIX E. THE ALGORITHM AND CODE FOR ONE DOMAIN C(2,l) = A(2,1)*B(1,1) + A(2,2)*B(2,1) + C(2,2) = A(2,1)*B(1,2) + A(2,2)*B(2,2) + A ( 2 , 3 ) *B (3 , 2) C(2,3) = A(2,1)*B(1,3) + A(2,2)*B(2,3) + A(2,3)*B(3,3) C(3,l) = A(3,1)*B(1,1) + A(3,2)*B(2,1) + A(3,3)*B(3,1) PROBLEM A(2,3)*B(3,1) C(3,2) = A(3,1)*B(1,2) + A(3,2)*B(2,2) + A(3,3)*B(3,2) C(3,3) = A(3,1)*B(1,3) + A(3,2)*B(2,3) + A(3,3)*B(3,3) return end subroutine Add3x3(A, B, Factor, double precision A(3,3), C(l,l) = A ( l , l ) + Factor C(l,2) = A(l,2) + Factor C(l,3) = A(l,3) + Factor C(2,l) = A(2,l) + Factor C(2,2) = A(2,2) + Factor C(2,3) = A(2,3) + Factor C(3,l) = A(3,l) + Factor C(3,2) = A(3,2) + Factor C(3,3) = A(3,3) + Factor C) B(3,3), Factor, * * * B(l,l) * * * B(2,l) * * * B(3,l) C(3,3) B(l,2) B(l,3) B(2,2) B(2,3) B(3,2) B(3,3) return end subroutine AddVec(A, double p r e c i s i o n B, Factor, A(3), B(3), C(l) = A(l) + Factor * B(l) C(2) = A(2) + Factor * B(2) C(3) = A(3) + Factor * B(3) C) Factor, C(3) return end subroutine double Invert3x3(Block, precision Block(3,3), Inverse) Inverse(3,3) double p r e c i s i o n Detlnv Detlnv Block(l,l)*Block(2,2)*Block(3,3) = 1. / (+ & + Block(l,2)*Block(2,3)*Block(3,l) & + Block(l,3)*Block(2,l)*Block(3,2) & - Block(l,3)*Block(2,2)*Block(3,l) & - Block(l,2)*Block(2,l)*Block(3,3) & - Block(l,l)*Block(2,3)*Block(3,2)) 152 APPENDIX E. Expand THE ALGORITHM by m i n o r s Inverse(l,l) to AND compute the = + Detlnv * CODE FOR ONE DOMAIN PROBLEM inverse (Block(2,2)*Block(3,3) - & (Block(2,l)*Block(3,3) - & (Block(2,l)*Block(3,2) - & (Block(l,2)*Block(3,3) - & (Block(l,l)*Block(3,3) - & (BlockQ, l)*Block(3,2) - & (BlockQ,2)*Block(2,3) - & (BlockQ , 1)*Block(2,3) - & (BlockQ , 1)*Block(2,2) - & Block(3,2)*Block(2,3)) Inverse(2,l) = - Detlnv * Block(3,l)*Block(2,3)) Inverse(3,l) = + Detlnv * Block(3,l)*Block(2,2)) Inverse(l,2) = - Detlnv * Block(3,2)*Block(l,3)) Inverse(2,2) = + Detlnv * Block(3,l)*Block(l,3)) Inverse(3,2) = - Detlnv * Block(3,l)*Block(l,2)) Inverse(l,3) = + Detlnv * Block(2,2)*Block(l,3)) Inverse(2,3) = - Detlnv * Block(2,l)*Block(l,3)) Inverse(3,3) = + Detlnv * Block(2,l)*Block(l,2)) return end I END o f subroutine SolveBlockTri 153 Appendix F D E T A I L E D D R A W I N G S OF T H E A P P A R A T U S F.l D e t a i l e d D r a w i n g of t h e A p p a r a t u s 154 APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS Figure F . l : the 3-D view of the cover of the upstream-tank 155 APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS 156 APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS Figure F.3: the 3-D view of the lid of the test channel APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS Figure F.5: the 3D view of the screen plate (filter) for the suction APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS 160 OOOOOOOOOOOOOOOOOOOOOOOOO o o o o o o o o o o o o o o o o o o o o o o o o o OOOOOOOOOOOOOOOOOOOOOOOOO o o o o o o o o o o oo o o o o o o o o o o o o o OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO O O OO OO OO OOOOOOOOOO OOO O OO O OOOOOOOOOOOOOOOOOOOOOOOOO O O O O O OO O O O O O OO OO O O O O O O O O o OOOOOOOOOOOOOOOOOOOOOOOOO OO OO O O O O O O O O O 0 O O O O O O O O O O O OOOOOOOOOOOOOOOOOOOOOOOOO O OO OOO OO OO OOOOOOOOOO OO OO O OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO O O O O O OO O O O O O O O OO O O O O O O O O O O OOOOOOOOOOOOOOOOOOOOOOO o OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOOOOOOOOOOOO OOP OOOOOOOOOOOOOOOOOOOOOOOOO O O O C ' O O O O O O O O O O O O O O O O O O O OrO OOOOOOOOOOOOOOOOOOOOOOCK^) CO buO APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS Figure F.7: the 3-D view of the suction element (diffusor) of the device APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS 162 APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS 163 Figure F.9: the 3D view of the upstream and downstream tanks connected with a channel APPENDIX F. DETAILED DRAWINGS OF THE APPARATUS 164 APPENDIX F.2 F. DETAILED DRAWINGS OF THE APPARATUS 165 D e s c r i p t i o n of the P u m p The centrifugal pump we used is an all bronze pump purchased from Pricepump (http: //www.pricepump.com), the following is the brief description: • Model CD150-AB-362-6A111-100-36-1T6 • 1 horse power, motor speed 3600 R P M • the range of flow rate is 25 to 160 G P M • the optimum working region is: 75-120 G P M at the pressure head of 30-50 feet F.3 Seeding particles The seeding particles are the polymeric particles produced by EMS-CHEMIE AG in Switzerland (http://www.ems-group.com/ en/indexl.html). The product is named Griltex 2A PI, whose average particle size is around 60 microns with a density a little bit above that of water (1.04) but lower than our fluid (1.085). It has to be noted that before adding the particles into the fluid, a little bit of soap was added to provide reasonable wetting of the particles, together with particles into a beaker of water and to thoroughly mix before adding to the reservoir. The amount of the particles to be provided has to be sufficient to get enough echo at medium emitting power.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Linear stability of a Berman flow in a channel partially...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Linear stability of a Berman flow in a channel partially filled with a porous medium Deng, Chuntao 2004
pdf
Page Metadata
Item Metadata
Title | Linear stability of a Berman flow in a channel partially filled with a porous medium |
Creator |
Deng, Chuntao |
Date Issued | 2004 |
Description | In this thesis, the coupled flow of a Newtonian fluid both above and through a porous medium is considered. This work was motivated from a pulp and paper application, namely twin-wire forming. In the fluid-only region, the two-dimensional flow field is governed by the Navier-Stokes equation. In the porous region, the flow field is governed by the Brinkman-extended Darcy law relationship. Inertial terms are retained in both regions and the interface conditions between the two domains are those as outlined by Ochoa- Tapia and Whitaker (Int. J. Heat Mass Transfer 38, 2635 1995). The model equations were solved using two independent methods. In the first method we develop a similarity transform and reduce the governing equations to two, coupled, non-linear ordinary differential equations to form a three-point boundary value problem. This was solved numerically and validated analytically by examining asymptotic cases. Three characteristic solutions were identified and the stability of each was examined by the method of normal modes. In the second numerical approach, the governing equations were re-posed as a one-domain problem, using the procedure outlined by Basu and Khalili (Phys. Fluids 11, 1031 1999), so that the conditions at the interface need not be considered. The resulting equation was solved directly, in primitive variable form, using a finite volume formulation. Finally, an experimental device was constructed to compare to the numerical predictions. Eight test cases were performed, using two different porous media, in which the velocity profile of the fluid was measured using pulsed ultrasound doppler anemometry (PUDA). Good agreement was found between the similarity numerical predications and the experimental measurements. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0058953 |
URI | http://hdl.handle.net/2429/17303 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2005-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2005-994651.pdf [ 9.48MB ]
- Metadata
- JSON: 831-1.0058953.json
- JSON-LD: 831-1.0058953-ld.json
- RDF/XML (Pretty): 831-1.0058953-rdf.xml
- RDF/JSON: 831-1.0058953-rdf.json
- Turtle: 831-1.0058953-turtle.txt
- N-Triples: 831-1.0058953-rdf-ntriples.txt
- Original Record: 831-1.0058953-source.json
- Full Text
- 831-1.0058953-fulltext.txt
- Citation
- 831-1.0058953.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0058953/manifest