Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Linear stability of a Berman flow in a channel partially filled with a porous medium Deng, Chuntao 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items