UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A viscous-inviscid interaction procedure Stropky, Dave 1988

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

Item Metadata


831-UBC_1988_A7 S77.pdf [ 6.05MB ]
JSON: 831-1.0097925.json
JSON-LD: 831-1.0097925-ld.json
RDF/XML (Pretty): 831-1.0097925-rdf.xml
RDF/JSON: 831-1.0097925-rdf.json
Turtle: 831-1.0097925-turtle.txt
N-Triples: 831-1.0097925-rdf-ntriples.txt
Original Record: 831-1.0097925-source.json
Full Text

Full Text

A VISCOUS-INVISCID INTERACTION PROCEDURE By  Dave Stropky B.A.Sc,  A  T H E S I S  T H E  University of British Columbia, 1983  S U B M I T T E D  I N  P A R T I A L  R E Q U I R E M E N T S  M A S T E R  O  F  F  O  R  T  A P P L I E D  H  F U L F I L L M E N T  E  D E G R E E  O  O  F  S C I E N C E  in T H E  F A C U L T Y  O  F  G R A D U A T E  S T U D I E S  Department of Mechanical Engineering  We accept this thesis as conforming t o the required s t a n d a r d  T H E  U  N  I  V  E  R  S  I  T  Y  O  F  B R I T I S H  April 1988 © Dave Stropky, 1988  C  O  L  U  M  B  I  A  F  In  presenting  this thesis in partial fulfilment  of the  requirements for an advanced  degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department  or  by  his  or  her  representatives.  It  is  understood  that  copying  or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department of  /VJ&T'/M/CAI  E7V6W  The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3  DE-6G/81)  A  bstract  A new viscous-inviscid semi-inverse (VISI) interaction method has been developed for predicting the flow field arising from a combination of inviscid potential flow and viscous flow. The technique involves matching the bounding velocities for each region by iteratively solving for the displacement thickness, 6*(x). The formula used to update S*(x) after each iteration is generated by linearly perturbing the governing equations. Application of the VISI procedure to predict the unseparated flow past a flat plate gives excellent results, producing numerical solutions essentially indistinguishable from the appropriate analytical solution in less than 0.5 seconds of CPU time on an Amdahl 5850 computer. Application of the technique to external flow over a backward facing step (BFS) indicates that the region of strong interaction between the viscous and inviscid flows is much larger than reported for internal flow. Calculated reattachment lengths, LR, are clearly influenced by the thickness of the boundary layer upstream of the step, thicker boundary layers producing longer reattachment lengths. Good accuracy is achieved for a relatively coarse distribution of control points, and rapid convergence (< 2 seconds on the Amdahl 5850) usually occurs. Finite-difference predictions using an elliptic code (TEACH-T), modified at the outer boundary to simulate external flow, have also been made for the BFS, largely as a basis of comparison for the VISI results. Comparison of results for the two models (VISI and TEACH) gives similar trends in LR as a function of Rh and x, (a measure of the displacement thickness at the step). The values of LR obtained with the VISI method, however, are 15-80% longer than those from T E A C H . Direct comparison with experiments is difficult because the experimental data does not clearly identify the effects of x, in the resulting values of LR. Trends appear to be the same for all computed and observed cases however. Disagreement between the VISI and T E A C H results is thought to be due to a combination of neglecting velocities in the recirculation region in the VISI model, and numerical error and inaccurate boundary conditions in the T E A C H code. ii  C ontents  ii  Abstract  viii  Nomenclature  xi  Acknowledgements 1  2  3  INTRODUCTION  1  1.1  Problem Description  4  1.2  Literature Review  6  1.2.1  Theory  6  1.2.2  Experiments  8  A NEW  VISCOUS-INVISCIDP R O C E D U R E  10  2.1  Current Interaction Schemes  10  2.2  A New VISI Technique  12  2.2.1  13  A Simple Example of the VISI Technique  2.3  The General VII Method  18  2.4  The Governing Equations  29  2.4.1  The Inviscid Region  29  2.4.2  The Viscous Region  31  A P P L I C A T I O N TO 3.1  THE BLASIUS PROBLEM  37  Problem Geometry  38  3.1.1  Boundary Conditions  39  3.1.2  Asymptotic Limits  41  iii  CONTENTS  iv  3.2  Problem Formulation  43  3.2.1  Velocity Matching Locations  43  3.2.2  Iteration Method  45  3.2.3  Velocity Calculations  46  3.2.4  Perturbation Calculations  48  3.3  4  3.3.1  Radius of Convergence and Stability  56  3.3.2  Solution Accuracy, Grid Refinement, and Execution Speed . .  59 62  4.1  Problem Geometry  63  4.1.1  64  4.3  The Displacement Line  Modifications to the Governing Equations  65  4.2.1  The Inviscid Equation  65  4.2.2  The Viscous Equations  66  Parameter Calculations  74  4.3.1  Velocity Calculations  74  4.3.2  Velocity Perturbation Calculations  75  TEACH  78  5.1  Computational Procedure  78  5.1.1  Control Volume Formulation  79  5.1.2  Hybrid Differencing and Linearization  80  5.1.3  Solution Procedure  81  5.2  6  56  APPLICATION TO A BACKWARD FACING STEP  4.2  5  Results and Analysis  Problem Geometry  81  5.2.1  Special Cells  83  5.2.2  Previous Work and Present Refinements  83  RESULTS A N D ANALYSIS  85  6.1  VII Results  85  6.1.1  Grid Refinement  86  6.1.2  Parametric Variation  90  6.2  T E A C H Results  95  v  CONTENTS  6.3  7  6.2.1  Grid Refinement  95  6.2.2  Velocity Profiles  96  6.2.3  Parametric Variation  97  Comparison of Results  101  6.3.1  Reattachment Lengths  102  6.3.2  Other Characteristics  104  6.3.3  Execution Times  113  CONCLUSIONS A N D RECOMMENDATIONS  115  7.1  Conclusions  115  7.2  Recommendations  120  References  122  APPENDICES  124  A  124 A.l  Asymptotic Solution to the Boundary Layer Equations Integrated along A* = a x — x 2 V RK C  B B.l  Calculation of Asymptotic Value for UPOT{X)  B. 2 Calculation of U T{X) for A*(x) = Ky/x PO  C  124 127 127 130 131  Cl  Upori^k) Calculations  C. 2  AL^por( *)p x  Calculations  131 133  List  of  Figures  1.1  VI Flow  3  1.2  Backward Facing Step  5  2.1  VI Procedures  11  2.2  Ellipse and Square  13  2.3  Spliced Cubics  20  2.4  Endpoint Slopes  21  2.5  Spliced Parabolae  24  2.6  Modified Perturbation Curve  25  3.1  Blasius Problem Geometry  38  3.2  Matching Point Locations  44  3.3  Perturbation Curve for n = N  48  3.4  p ^ N Subcase Locations  53  3.5  p = N Subcase Locations  54  3.6  Iteration History  58  3.7  Blasius Grid Refinement  60  4.1  BFS - Problem Geometry  63  4.2  The Displacement Line  65  4.3  Viscous Regions  66  4.4  Separation Profiles  68  4.5  Reattachment Singularity  73  5.1  Scalar Cell  79  vi  LIST OF FIGURES  vii  5.2  Boundary Conditions  82  5.3  Grid Layout  83  6.1  V I I Grid Refinement - Effect of x  6.2  VII Grid Refinement - Effects of x  6.3  VII Grid Refinement - Pressure Distribution  89  6.4  Pressure Coefficient - Effect of x,  91  6.5  Pressure Coefficient - Effect of T  92  6.6  Pressure Coefficient - Effect of R  94  6.7  Velocity Vectors near the Step  96  6.8  T E A C H - Pressure Coefficient - Effects of x,  97  6.9  T E A C H - Pressure Coefficient - Effects of R  99  6.10  T E A C H - Pressure Coefficient - Case 2  100  6.11  Reattachment Lengths  102  6.12  Wall  6.13  Refined W a l l Pressure Comparison  106  6.14  Displacement Line Comparison  107  6.15  Shape Factor Comparison  109  6.16  Separation Streamline and Zero Velocity Line Comparison . . . .  Ill  6.17  86  N  88  0  h  h  Pressure Comparison  Shear Stress Coefficient Comparison  104  112  Nomenclature  Symbols a  ellipse area constant  a , a, a, a n  s  e  w  coefficients of finite difference equation  a, b, c, d  coefficients of cubic in ^-space  A, B, C, D  coefficients of cubic in x-space  A,  areas of ellipse and square  As  E  A  velocity profile coefficient  C, D  flux coefficients  Cf  shear stress coefficient  C i , . . . , C-j  integral-momentum coefficients  IU-'-IIN  V'-space cubics  n  n  T\,-.-,7  x-space cubics  7o, ?E, 9B  zero pressure gradient curves  h  step height  N  hi,..., h  N  cubic endpoint values  tU.  viscous perturbation coefficient  I, J  bounding velocity functions  LR  reattachment length  p  pressure  P  Peclet number (=  e  puAx/T)  Q  cubic expansion ratio  r  ellipse and square area parameter  r*  solution value for r  R  Reynolds number (= U^h/v)  h  S,  source terms  p  u, v  finite-difference  velocities  UE  bounding velocity  U  undisturbed velocity  EO  U'  E  x, x, y  velocity disturbance cartesian coordinates  viii  NOMENCLATURE  ix  x,  step location  x  reattachment coordinate  R  x, x  interaction region endpoint locations  x  equation switch location  N  0  e  x'  virtual origin  a  Pohlhausen curve coefficient  T  Diffusion coefficient for 4>  6, A  boundary layer thickness  A/  recirculation bubble height  A  viscous boundary coordinate  0  61, A i  displacement line  6*, A*  displacement thickness  AX  cell width  AXI, ..., AX  N  cubic lengths  r?  velocity profile variable (= y/8)  0, 0  momentum thickness  A  pressure gradient parameter  fi, v  viscosity  IT  3.1415...  p  fluid density  a  constant  r  shear stress  <f>  general transport variable  ip  normalizing transformation variable  Jl  perturbation coefficient  Subscripts A, B, C, D  equation switch values  BL  boundary layer  E  ellipse  crit  critical value  IN  inviscid  NOMENCLATURE  x  k  velocity matching location index  n  cubic function index  N  number of cubics  n, s, e, w  value at cell face  N, 5, E, W  neighbouring cell value  p  perturbation location index  POT  potential  S  square  T  velocity profile order  VIS  viscous  oo  free stream condition  Prescripts  K  influence coefficient  A  perturbation  Abbreviations  BFS  backward facing step  b.l.  boundary layer  VI  viscous-inviscid  VII  viscous-inviscid interaction  VISI  viscous-inviscid semi-inverse  A cknow ledgements  I would like to express my sincere thanks to both of my supervisors, Professor I. S. Gartshore and Professor M . Salcudean. Professor Gartshore, with his infinite patience, allowed me to go my merry way with this thesis, and always provided knowledgable guidance when I strayed off course. Professor Salcudean provided a keenly balanced source of knowledge, stimulation and encouragement during some of the most critical stages of this development. I would also like to thank Dr. Djilali for the references he provided and for his invaluable time-saving tips. I thank Paula for her love and support (and for doing more than her share of the dishes) and last but certainly not least, I thank God that this work is completed.  xi  Chapter 1  I N T R O D U C T I O N  Flow separation and reattachment are important considerations in the design of a wide variety of engineering components.  Consider the design of airfoils, tur-  bine blades, automobile shells and diffusers. In these cases flow separation is to be avoided under normal operating conditions as it is disruptive to the operating efficiency of the device. In some situations separation is desired, such as during aircraft landing, where high drag is desired. If the performance of these devices is to be well understood for any operating range, an understanding of separated flow is essential. Separated flows are among the most complex of flow situations encountered. Mathematical modelling of these flows has been studied for many decades, but solutions to the more complex cases were elusive until the advent of large digital computers. The computer has enabled numerical solutions of the non-linear Navier Stokes equations. Most numerical techniques involve solving the governing differential equations in discretized form, a method which is very powerful. Numerical solutions to several separated flow problems have been determined by this means. Although numerical methods of solving the Navier-Stokes equations have many 1  CHAPTER  1.  INTRODUCTION  2  advantages they also have some disadvantages. The techniques involve manipulation of large matrices and can get relatively expensive for complex problems. In addition, if the matrices are very large the method cannot be used on anything but a large computer.  Use of discrete methods as a tool to optimize design parameters is  more cumbersome than closed form analytic methods. The results of numerical parametric variation are difficult to predict and so optimization is performed by combining physical intuition with trial and error. Closed form analytic solutions to separated flow problems would be desirable, but finding these solutions seems unlikely owing to the complexity of the equations. In the last twenty years methods have been developed to reduce the order of computation needed to solve flows containing regions of separation. 'Viscous-inviscid Interaction' (VII) methods are probably the most widely used of these methods. The idea behind VII methods is illustrated in Figure 1.1. Most numerical methods use one governing set of equations for the entire flow field regardless of the nature of various regions of the flow. In many practical situations the flow can be divided into separate regions where, (i) viscous effects are important, and (ii) viscous effects are negligible. The separate regions are solved using equations which govern the flow in each corresponding region, with the restriction that these solutions must match at the interface between regions. These equation sets are generally reduced in order from the original equation set thus making their solution less complex. The benefits of this type of analysis are a reduction in both execution time and the amount of computer memory required. Purely inviscid and viscous boundary layer (b.l.) flows have been analysed in the last hundred years and much is known about their dynamics. Highly efficient analytical and numerical techniques have been developed to solve these two quite  CHAPTER  1.  INTRODUCTION  3  INVISCID  Figure 1.1  VI Flow  different problems. Many practical flow situations are ideally suited to the application of VII techniques. External boundary layer flows are a good example, as the large computational domain with its two distinct regions makes the use of global equation solvers very inefficient. In this situation VII techniques are ideal, and for example are the basis for many highly efficient and effective unseparated airfoil design and analysis codes. It seems possible, in principle at least, to extend the use of VII methods to situations of separated flow. In fact several internal separated b.l. flows have been solved usingfinite-differencetechniques for the viscous regions and panel methods for the inviscid core. It has also been shown, for fully developed flow, that the viscous regions could be solved by integral b.l. methods. Integral methods are simple, extremely efficient and yield accurate results in many cases. By the process of making the separated flow model as simple as possible (while still retaining the essential physics), new insight into the problem can be achieved. It is usually beneficial to solve complex problems in steps, starting with the simplest models and working towards the most complex. With the availability of large com-  CHAPTER  1.  INTRODUCTION  4  puters and many complex Navier-Stokes solvers, the initial steps of this 'buildingblock' scheme have been bypassed. Few attempts to make a simple mathematical model to solve separated flow problems have been made. The aim of this study is to develop such a simple model to gain further insight into the nature of separated flows.  1.1  Problem Description  The VII technique developed is quite general, but only the simplest types offlowsare considered here. The current model is developed for 2-D, steady, incompressible, laminar flow. Turbulent flow was not considered as turbulence modelling would add an extra dimension of uncertainty which would be difficult to separate from the inherent problems of model development. Turbulent conditions can be added after the model is working properly. Comparatively little work has been done using VII methods for external flow and even less for incompressible separated flow. As far as the author is aware, no work has been done on what forms the basis of this thesis: solving external incompressible separated flow problems using VII methods which employ integral techniques and linearized inviscid theory. The basic VII model will be developed and first tested on a simple case of unseparated flow. This basic model will then be modified to solve a specific case of separated flow. Many flow geometries could be used for the separated flow case, however the 'Backward Facing Step' (BFS) is singled out for the following reasons: • the geometry is simple. • the separation point is fixed.  CHAPTER  1.  INTRODUCTION  5  • the step size is a useful and convenient length scale. • there are no curved solid surfaces. • experimental data is available. • the BFS is a common test case for numerical methods. In addition the BFS contains a discontinuity in geometry which occurs in many situations and requires special consideration when using integral methods. Figure 1.2 shows several key features of external flow over a BFS. The inviscid region is  3  INVISCID  13  V/SCOUS  RECION RECION  Figure 1.2  A B C D  BOUNDARY LAYER EDCE DISPLACEMENT THICKNESS SEPARATION STREAMLINE ZERO VELOCITY LINE  Backward Facing Step  everywhere outside the b.l. whose definition is determined by the particular integral model chosen. In the viscous region, a boundary layer develops on approach to the step where it separates as a free shear layer. This shear layer reattaches some distance downstream and a new boundary layer develops on the lower wall. In a  CHAPTER  1.  INTRODUCTION  6  region surrounding the step, the b.l. interacts strongly with the inviscid flow. This makes normal b.l. theory inadequate and generates the need for VII methods. The flow patterns are similar for both internal and external flow but comparison between the two is not simple and will be addressed later in the thesis.  1.2 1.2.1  Literature Review Theory  Briley and M Donald [l] have conducted an excellent survey of VII theory. They give c  reference to the validation of VII solutions for several flow situations, including the B F S . In view of the number of validated results, it can be asserted that VII methods can and do work in certain situations involving separated flow. The reasoning behind interactive methods and the situations to which they apply are given in a paper by J . C . Williams [2]. The above two papers also summarize reoccurring problems and solution techniques associated with VII methods.  Three problems  which have caused much difficulty in the development of interactive methods are, (i) the Goldstein singularity, (ii) the velocity profile critical point, and (iii) development of a stable method for interacting the viscous and inviscid regions. Problem (i) involves a singularity encountered at separation which occurs when the boundary layer equations (used for the viscous region) are integrated with the pressure distribution prescribed (the 'direct' method). S. Goldstein [3] found an analytic solution, valid near separation, and singular at separation so that numerical methods could be employed through separation. Williams [2] describes how solutions free of this singularity were obtained by prescribing the displacement thickness  CHAPTER  1.  INTRODUCTION  7  rather than the pressure distribution (the 'inverse' method). Although a few authors such as Crimi and Reeves [4] have successfully used the direct method, the majority of solutions to separated flow problems are found through the use of inverse boundary layer methods. Problem (ii) is associated with integral-momentum methods.  Many integral  methods employ velocity profiles which are a function of more than one independent parameter. In these cases a singularity termed the 'velocity profile critical point' can arise during integration of boundary layer equations.  Shamroth [5] gives a  summary of the problems associated with integral methods, including the velocity profile critical point. Essentially the critical point arises when, at a certain point in the flow, the independent parameters in the profile function become linearly dependent. It occurs at different locations for different profile functions and can be avoided in the solution domain with a judicious choice of the profile function. The singularity is the result of the profile function chosen and is not a feature of integral methods in general. An example of the complicating effect of this singularity is given in a paper by Green [6]. In the present study velocity profiles are used which are a function of a single parameter and so this problem is not encountered. Having eliminated singularities from the VI solution, the problem (iii) of devel1  oping an interaction method or formula still remains. There are several currently available procedures, but these are generally problem specific. These methods can be placed into two major categories, (a) Those that solve the b.l. equations in differential form, (b) Integral b.l. solvers. Each of these can be further subdivided into internal and external flow. While no work has been done using VII methods for a BFS in external flow, Kwon and Pletcher [7] have used VII methods to suc1  Numerical solutions to the Navier Stokes equations do not contain these singularities  CHAPTER  1.  INTRODUCTION  8  cessfully solve this problem for internal flow. B.R. Williams [8] accurately predicts the post stall behavior of an airfoil, an example of the application of VII methods to separated external flow. Several other authors [4,6,9,10] have also employed VII methods, but their interaction techniques are rather cumbersome. Techniques such as graphical solutions, variable relaxation factors throughout the field, and even computer graphic interaction were required for solution in these cases. After a review of the literature it became apparent that the interaction procedures developed for VII methods needed improvement.  1.2.2  Experiments  Relatively little work has been done for a BFS in laminar flow. One of the first investigations was done by Moore [11] who found a reattachment length of 22 step heights for a Reynolds number of 338 based on displacement thickness at the step. He claimed that if this Reynolds number is less than 400 the flow will be laminar up to reattachment. R.J. Goldstein et al [12] did a careful study of laminarflowover a BFS and found two criteria must be met for the flow to remain laminar through reattachment. First the ratio of displacement thickness to step height must be greater than about 0.4, and second the Reynolds number based on step height must be less than about 520. Goldstein's results on reattachment length are unclear as to whether reattachment length is a function of both Reynolds number based on step height and b.l. thickness at the step. Armaly et al [13] claim that for internal flow the reattachment length is likely a function of three or more variables. Armaly's results, as with those of Denham and Patrick [14], are for fully devel-  CHAPTER  1.  INTRODUCTION  9  oped flow at the step and for area expansion ratios of 1:1.94 and 1:1.33 respectively. Conversely, Moore and Goldstein have thin boundary layers at the step compared to the channel height. Their area ratios were 1:1.02 and 1:1.06 respectively. The scatter of these results for reattachment length plotted against Reynolds number based on step height suggests that reattachment length is not a function of a single variable as many authors assume.  Chapter 2 A  N  E  W  V I S C O U S - I N V I S C I D P R O C E D U R E  2.1  Current Interaction Schemes  Viscous-inviscid methods offer an attractive way of dealing with separated flow, but these methods have been applied to this flow situation with some difficulty. Several iterative algorithms have been developed to compute the interaction between the viscous and inviscid layers. Each method presents its own particular set of difficulties and none are completely general. Almost all of the methods can be placed into three categories: direct, inverse and semi-inverse. A short dicussion of the meaning, use and limitations of each method follows. Figure 2.1 shows a schematic of the direct, inverse, and semi-inverse techniques. Most aerodynamicists are familiar with the direct VII procedure when applied to unseparated flow. In the direct scheme the displacement thickness, 6*(x), is applied to the inviscid equations which then yield the boundary layer edge velocity distribution, UE (X). IN  10  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  11  INVISCID  DIRECT  SOLVER INVISCID SOLVER INVERSE  CORRECTION  6*  FORMULA FOR  6"  <5*  INVERSE  VISCOUS viscous  SOLVER  SOLVER DIRECT SEMI-INVERSE  Figure 2.1  VI Procedures  This velocity distribution is then applied to the viscous equations to yield a new S*(x) and the process is repeated until convergence is obtained. The inverse method switches the roles of 6*(x) and UE{X). In this method 6*(x) is applied to the viscous equations which yield UE [X). VIS  This velocity distribution  is then applied to the inviscid equations to yield a new estimate of 6*(x) and so forth. The semi-inverse method applies an initial estimate of 6*(x) to both sets of equations which yield two estimates UE[X). These estimates are combined to generate a correction formula (see figure 2.1) which is used to update 6*(x) and the process repeats to convergence. The literature review shows that the viscous (b.l.) equations encounter a singularity at separation when solved in a direct fashion. Although this is not generally  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  12  true for integral methods, the direct method is limited as a possible model choice because of lack of generality. The inverse method presents difficulties when dealing with the inviscid equations in external flow, the large flow domain causing convergence problems. B.R. Williams [8] used Fourier analysis to study the stability of the three schemes. He concludes, "... the direct scheme is stable for attached flows but in general is unstable for separated flows. The fully inverse method is stable for separated flows but the relaxation factor is inversely proportional to the size of the computational domain which results in slow convergence for external aerodynamics where a large but finite computing domain is used. . . . the semi-inverse method of a direct inviscid calculation coupled with an inverse calculation of the boundary layer is the suitable choice for a mixture of attached and separated flows in external aerodynamics." For these reasons the semi-inverse method is chosen here as a basis of model construction. A new problem not encountered in the direct or inverse methods arises in that a correction procedure (see figure 2.1) must be developed to update 6*(x). A major part of the time spent on this thesis was used in developing such a procedure. The new viscous inviscid semi-inverse (VISI) procedure is quite general although it is better suited to integral methods.  2.2  A New VISI Technique  In deriving the proposed model, a few analytic and semi-empirical correction formulae suggested by other authors were tried. It was found that these iterative formulae worked marginally or not at all. In one case it was found that a permutation of the  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  13  formula worked but the unmodified formula did not. Iteration formulae have not been developed aiming at generality so that they can not be applied to different problems with confidence. It was undertaken to derive an iteration formula which would be general enough to use in many types of problems. The new technique is based on linearized theory and has all the limitations of linearized theory. The method is straightforward and can be applied to any set of equations governing the viscous and inviscid regions. The idea is to perturb the governing equations so that they themselves generate the correction formula. The method will be illustrated with a simple example in the following section.  2.2.1  A Simple Example of the VISI Technique  r r+a  ELLIPSE  Figure 2.2  I  r  SQUARE  Ellipse and Square  The technique will be applied to a non-fluid dynamic problem so that the equations can be kept simple. Figure 2.2 shows an ellipse of major semi-axis r + a, minor semi-axis r, and a square of side 2r. The area of the ellipse, A , and the area of the E  square, As, will be equated through variation of r using the proposed technique.  CHAPTER  2.  A NEW  VISCOUS-INVISCID  14  PROCEDURE  From analytic geometry the corresponding areas are, AE  =  irr(r + a) — irr + ixra  (2.1)  A  =  4r  (2.2)  S  2  2  The two area functions are then perturbed to generate their respective influence coefficients. The perturbations are done by giving a small displacement A r to the independent variable r and then invoking linearized theory with the assumption that A r <C r. Perturbing equations 2.1 and 2.2 for the areas of the ellipse and square and denoting the changes in these areas as AAE and AAS, AE + AAE  A +AA S  S  =  7r(r + A r )  + ixa(r + A r )  =  [7rr(r + a)] + [7r(2r+  =  4(r + A r )  =  [4r ] + [8rAr +  2  a)Ar +  7r(Ar) ] 2  2  2  4(Ar) ] 2  Neglecting powers of A r greater than unity and subtracting the areas from both sides, AAE  «  KEAT  A As  «  KsAr  where, K  E  =  7r(2r + a)  (2.3)  K  =  8r  (2.4)  s  KE and Ks are the linear influence coefficients with respect to r for the areas of the  CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE  15  ellipse and square respectively. The following steps outline the iterative procedure in which the solution, As = A , is obtained using the derived perturbation parameters. s  1. Guess r 2. Compute AE, AS, K , KS from equations 2.1, 2.2, 2.3, 2.4 E  3. If AE — As go to step 7. 4. Since AE 7^ As a perturbation to r is required which will make the equality hold, ie: A  E  + AA = A + AA s  E  S  or AE  + KE&T = As + KsAr  Rearranging and substituting for AE, AS, KE, _ A  5. Add A r to r, ie: r  neu)  f  AE-AS  _ ?rr(r + a) - 4r  " l T 5 - t f j s ~  = r  old  and Ks, 2  8 r - j r ( 2 r + a)  1  J  +A r  6. Go to step 2. 7. The solution is AE = As at r = r*. Numerical Results and Analysis of the Sample Problem In this section various numerical results will be shown for two different initial guesses (the parameter a is arbitrarily set to 1). The extreme simplicity of this problem allows calculation of an exact analytical solution: r* = 7r/(4 — tr) « 3.65979. It is  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  16  also noted that this problem has a degenerate solution, r * = 0, a solution that is to be avoided. For the first run, r = 40 will be chosen as an initial guess to show some characteristics of the iterative procedure. The following table summarizes the results of the computation procedure outlined above. ITERATION  r  A  E  A  s  A r  0  40.0000  6400.0000  5152.2120 -19.0412  1  20.9588  1757.0871  1445.8570  -9.4769  2  11.4819  527.3342  450.2387  -4.6525  3  6.8294  186.5601  167.9790  -2.1648  4  4.6645  87.0307  83.0077  -0.8267  5  3.8379  58.9164  58.3298  -0.1702  6  3.6677  53.8077  53.7829  -0.0079  7  3.6598  53.5768  53.5768  -0.0002  Convergence is somewhat slow at first, but increases rapidly near the solution. This is because at first A r = r so that the approximations made in the linearized theory are no longer accurate. It can be shown that the solution converges for any initial guess r which is greater than the exact solution value. What happens if the initial guess is too low? The following table illustrates the case of initial guess r = 1.  CHAPTER  2. ANEW  VISCOUS-INVISCID  ITERATION  r  PROCEDURE  A  A  A r  s  E  17  0  1.0000 6.2832 4.0000  -1.6025  1  -.6025  -.7524  1.4520  +0.5279  2  -.0746  -.2169  0.0223  +0.0732  3  -.0015  -.0046  ~ 0  +0.0015  4  ~ 0  ~ 0  ~ 0  ~ 0  Here something happened which was to be avoided, the solution converged to the degenerate solution of zero area. The convergence is not monotonic and during the iteration process physically unrealistic negative values of r and A From equation 2.5,  A r  —> oo as 8r—7r(2r+a) —+0, or  r  =  r  c  n  =  t  E  are obtained.  \ x [7ra/(4 — 7r)]  -  | r * .  This is the critical value of r below which the solution converges to the degenerate solution,  r* =  0. Also, if the initial guess is very close to  r  c r  ,t  the  A r  computed will  be highly inaccurate. If the problem had been a linear one, the solution would be found in one iteration and there would be no critical points. As the problem gets more complex and nonlinear there can exist multiple solutions and multiple critical points. If there are critical points too near the desired solution this method is ineffective as the intial guess will require great care. Fortunately, many complex non-linear problems have single solutions and critical points far from the desired solution (ie. large radius of convergence). For these problems this method becomes highly effective because of its speed and relative simplicity.  CHAPTER  2.3  2. ANEW  VISCOUS-INVISCID  PROCEDURE  18  The General VII Method  How does the example of the previous section relate to a real VII problem? First, replace the independent variable r with the displacement thickness function 6*(x). Next replace the area functions A  E  and As with the viscous and inviscid boundary  layer edge velocity distributions, U (x) Evis  and U (x). ElN  In the simple example a  single variable r is perturbed by a small amount A r to generate the linear influence coefficients necessary for the VISI solution method. If 6*(x) is a function of JV independent variables h  n  (the value of 6*(x) at n different locations say), then  these h can be perturbed by A/I„ in an analogous fashion to yield the linearized n  perturbation function, A6*. This perturbation function can be subdivided into linear functions of the A/I„ whose coefficients are the linear influence coefficients similar to K  E  and Ks in the simple example. The boundary layer edge velocities can then  be matched using the new VISI technique, in a manner analogous to matching the areas of the ellipse and square. Solving for 6* (x) exactly makes the solution of a VII problem extremely complex. A common technique used in reducing the order of complexity is to satisfy the equations at a finite number of points and use interpolation between these points. This method is adopted for use in the general VII procedure. The general procedure which follows is quite straightforward but the algebra neccessary for its derivation can be a distraction. The following list is an outline which can be referred to while following the derivations. [Step 1] Devise a function 6*(x) whose shape can be controlled by a discrete number of independent parameters, h . n  [Step 2] Calculate the perturbation function A6*(X) in terms of the h and Ah . n  n  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  19  [Step 3] Linearize the perturbation function with respect to the Ah in order to n  calculate the influence coefficients necessary for solution.  Devising ^*(x) If 6*(x) is a simple function like 8*(x) = h\x + h , hi and hi being the perturba2  tion parameters (analogous to r), the perturbation function and the corresponding influence coefficients are quite simple to calculate, ie:  8* + A8*  =  (hi + Ahi) x + (h + Ah )  =  (h x + h ) + (A/IX x + A/I )  2  x  2  2  2  Subtracting 8* from both sides, the influence function, A8*(X) = XAhi + A/12, is obtained. This function is already linear in Ahi and Ah so by inspection the linear 2  influence coefficients for the parameters hi and h  2  are x and 1 respectively.  A  function like hix + h is very simple to use but unfortunately does not have enough 2  degrees of freedom to solve a typical VII problem. A method of 'spliced cubic polynomial functions' (also known as 'cubic splines') is devised to give a variable degree of freedom to 8*(x) while keeping necessary algebra to a minimum. Figure 2.3 shows 8*(x) constructed of spliced cubics in the solution domain. From figure 2.3, 8*(x)  =  J (x)  (x -i<x<x )  n  n  (2.6)  n  where, /„ (x) = A + B x + C x  2  n  n  n  + Dx n  z  (2.7)  CHAPTER  2. A NEW VISCOUS-INVISCID  Figure 2.3  20  h  hn  h n-t  h  PROCEDURE  Spliced Cubics  The cubics are constructed so as to have continuous slope and value at their endpoints. This gives four boundary conditions for each cubic, ie:  >  The coefficients A ...D n  n  (2.8)  of equation 2.7 are completely determined if the slope  and value at each end of the curve 7 (x) are supplied. The values at each end of the n  curve (see figure 2.3) are given by h and h -\. n  n  If the slope at each end is expressed  entirely as a function of the /i's then B*(x) will be completely determined by these /l's. Figure 2.4 illustrates the method used in calculating the slope at the endpoints using endpoint values. The slope at any endpoint is taken as the average slope of the straight lines drawn from that endpoint to the two neighbouring endpoints, ie:  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  21  4  Figure 2.4  Endpoint Slopes  similarly, d  \  7(  1  AX„_i  (2.10)  AX„  where, AXy =  X,- — X , _ i  j'=n-l,n,n+l  Now 7 (x) (and thus 6*(x)) can be expressed entirely as a function of / i n  h , and / i n  n  +i.  Expressing 6*(x) in this manner makes the h  n  n  fe -i>  _ ,  n  2  the independent  parameters (analogous to r) rather than the A ... D of equation 2.7. This is done n  n  to reduce necessary algebra and give more intuitive 'feel' to the problem. To reduce subsequent algebra further, a linear normalizing transformation is introduced, x — X -i n  Cn-1  < x< x  n  (0 <  rp <  (2.11)  1)  AX„  This transformation scales all the cubics to the same unit length. The entire VII problem will be solved in ^-space instead of the physically real x-space. The linear transformation of equation 2.11 transforms the cubic polynomials 7 [ ) into x  n  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  22  cubic polynomials, f (tp) , where n  fnW =  a + b ip + c ^ + d n  n  n  (2.12)  t/j  3  n  with the relations to J„(x), 7n(x)  (2.13)  From equations 2.6 and 2.13,  X -1  < X < X  n  The next step is to write f (if>) as a function of  /i„_  n  left endpoint is  fe„_i  h.  and the right is  (0 <  n  2  ...  %p <  fen+i-  1)  The value at the  The slopes at each end are given by  n  equations 2.9 and 2.10. These give the following four boundary conditions (in ipspace) sufficient to determine  [1]  fn{0)  [2]  fn(l) =  [3]  />)  =  =  f {4>) n  K_ 1 = a« K = a + b + c + d n  =  n  n  A Z „ — J(x -i)  Qn n  n  ~ fe„- ) + (^ -  (fen-l  Ax -^7(x ) n  n  = b  n  • i l  W  completely.  fen-l)  2  = b + 2c + 3d n  n  n  1 (fen+1  2  fen)  +  Q n  =  (fen  where and  A X  n  _ i  fen-l)  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  23  Solving the above for a ... d , n  On  =  b  =  n  C  «  =  d  n  n  ^ n - l  §[(A»-/l„-l) +  | ( ^ n  ~  fen-l)  -  Qn(/l»-l-*n-2)]  2QI I(^"+ +  1  ~  ^n) ~ Qn[^n-1 ~ h - )  (2.14)  n 2  = -(h - h -i) + t^{hn+l ~ hn) + ^(hn-i - h - ) n  n  n  2  Step [1] of the development outline on page 18 is now complete in that S*(x) is expressed entirely as a function of the h . Equation 2.14 needs to be modified n  for n = 1,2, N. These modifications depend on the boundary conditions for the problem being analysed. Examples are shown in Chapters 3 and 4.  The Since  perturbation function A6*(X) / (t/') n  represents 6*(x) and is a function of the h , perturbing the h in n  n  equation 2.14 gives the perturbation function t\f (ip) and thus A6*(X) . This method n  produces what are termed local perturbation curves. This means that if a single h  n  is perturbed by A/I„ then only a small portion of 6* (x) surrounding the perturbation is affected. This was the reason for selecting spliced cubics as they are the lowest order polymomial which give continuous slope and value as well as local perturbation curves.  Spliced parabolae were tried as they have continuous slope and value,  but they produce global perturbation curves. Figure 2.5 shows this global effect when 6*(x) is constructed of spliced parabolae. The smooth line is S*(x) before the perturbation, and the wavy line is 6*(x) after one of the parabolae endpoints is perturbed by an amount Ah. This waviness results because the unperturbed  CHAPTER  2. A NEW VISCOUS-INVISCID  Figure 2.5  PROCEDURE  24  Spliced Parabolae  endpoints are fixed and continuity of value and slope is always enforced. Note that although only one parabola is affected upstream of the perturbation, all parabolae downstream are influenced. Note also the large change in slope that occurs at the end of the last parabola, where some boundary condition must be applied. Together these two effects produced slow, oscillatory-type convergence.  The convergence  rate was substantially improved by using spliced cubics, and in almost every case convergence was monotonic. Perturbing equation 2.14 leads to perturbation curves which affect two cubics on each side of the perturbation. This is good but there is a way to get the perturbation even more localized. The perturbation curve can be constructed such that only one cubic on each side of the perturbation is affected. This reduces necessary algebra by a factor of two and increases convergence speeds over perturbation curves produced by equation 2.14. Figure 2.6 shows how the two-cubic perturbation curve is implemented. The index p refers to the particular cubic endpoint height that  CHAPTER  2. A NEW VISCOUS-INVISCID  Figure 2.6  PROCEDURE  25  Modified Perturbation Curve  is being perturbed. The perturbation functions, A / and of as the incremental changes to the cubics f  and /  p  p + 1  can be thought  A/ +I,  p  P  when h is given a disp  placement A/I . The perturbation curves are formulated as cubic polynomials with p  coefficients which are the linear increments to the coefficients of the original cubics. The perturbation curves are given as,  A  A/  /p(V0 P + 1  (0)  AO + A6 V > + A C 0 + A d V' 2  =  p  P  =  Aa  p + 1  8  P  + A6  p + 1  p  rp + A C  rp + A d  (2.15)  2  P + 1  p + 1  V  s  If <5*(x) is to remain unchanged outside the perturbation region (x -i < x < x i ) p  p+  and continuous everywhere, the slope and value must not change at x _i and x i . p  These conditions give four boundary conditions,  A/ (0)  =  0  A/>)  =  0  =  0  p  A/  p+1  (1)  p+  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  26  Continuity of A<5*(X) in the perturbation region at x gives the additional three p  constraints,  A/ (l) P  A/p (0) +1  A/I(l)  These seven constraint equations are insufficient to determine the eight coefficients of A / and p  A/  P +  I  in terms of  A/I , p  thus one extra independent equation is required.  It is clear that the perturbation to A6*(X) must be small so that linear theory is applicable. Minimizing the arc length of the perturbation curves is perhaps the best method, but this leads to an elliptic integral equation and so is discarded for reasons of complexity. Another possibility is minimizing the r.m.s. area under the perturbation curves. This leads to the following equation with the coefficient A C P being chosen as the independent parameter although any of the eight coefficients in equations 2.15 would have sufficed.  [^(A/p)  2  dAC.  /Wp)  2  dxP + Q  o  p+1  +1  W  = 0  The factor Q +\ is different from 1 if the two perturbation cubics are of different p  lengths. Using the seven constraint equations to write  AOP  ... A d  p +  i  in terms of  ACP,  substituting into the equation above and solving, Q ACE  2 + 1  (ll + 6 Q 2(1 + Q  3 P  +1  p + 1  )  )-5  x  A/I„  (2.16)  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  27  The parameter c is now a linear function of Ah . Next the other seven coefficients p  p  are written as linear functions of &h using equation 2.16 and the original seven p  constraint equations.  Ad  p  Ab Ad  p  p  Adp+i  The functions Af (ip) p  -  0  =  0  =  Ah  =  Ahp  =  Q p + l  =  -2Ah  =  Ah, (2 + 3 Q  — AC  p  p  (2.17) (3a/lp — A C p ) p  (1 + 2 Q p +  p +  i ) + 2Q  P + 1  Q iAc  i) -  p+  AC  P  p  and A/ (T/>) and hence (because the perturbation is local) P+1  A6*(X) are now completely determined by A/I alone. The preceding formulae will p  be modified at the end of the interaction region (p=N) depending on the nature of the problem being investigated. Examples of these modifications are given in Chapter 3. Step [2] of the outline on page 18 is now complete. Influence Coefficients The influence coefficients for Ah are easily obtainable. For example, the variable p  c has a linear influence coefficient, Kc , due to Ah and is given by the equation p  p  p  ACp = Kc Ah . From equation 2.16 by inspection, p  p  Kc  Q  2 p+1  (ll +  r  2(1 +  6Q  P  +  QUi)  1  )-5  (2.18)  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  Similarly the influence coefficients for the variables a ...d p  28 are obtained from  p+1  equations 2.17. Further, the influence function of / (V>), corresponding to the perP  turbation Ah , can be found using the equation Af (tp) p  p  [ Kf (ip) ] Ah , A/p(0)  —  p  p  being the linear influence function. From equations 2.15 and 2.17 , /P(V>) = AC p 02 + [Ah - AC ) ip  3  A  p  or, since AC = Kc Ah p  P  p  p  Af (i>) = [l + Kc W -^)}  X Ah  2  p  p  p  and therefore, Kf {1>) = 1 +  Kc ^ -^)  (2.19)  2  p  p  Similarly, = 1 - 302 + 20s + Q p + i[ 3(0 - 202 + 0s) - Kc (^ + 202 _ ^ ) ]  Kf ty) p+1  p  (2.20) Knowing that Ah only affects the cubics / p  p  and /  p + 1  , it is clear that all influence  functions Kf {ip) corresponding to the perturbation Ah are zero unless n = p, p+1. n  p  Having described <5*(x) and its influence functions entirely as functions of the /i's, the method is now ready for application to a real VII problem. Up to now the type of VII problem to be solved has not influenced any of the derivations. In this light the method is completely general. In order to best illustrate the new VISI method, specific flow conditions will be imposed. As discussed on page 4 the conditions will be 2-D, steady, laminar, incompressible, external flow, with the added constraint that the boundary layer equations be valid in the viscous region.  CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE  2.4  29  The Governing Equations  The proposed semi-inverse procedure requires that the interface velocity distributions UE  VIS  [X) and UE, (X) be matched at a discrete number of points in the solution N  domain. Only problems with thin viscous regions and small separation regions will be examined, so that boundary layer theory can be assumed to apply. A discussion of this assumption pertaining to regions of separated flow is given in Chapter 4. A modified von Karman integral method will be used to obtain UE (X). In the invisVIS  cid region it is assumed that the external flow UE (X) will not be greatly disturbed JN  by the presence of a small separation region so that linearized potential theory is thought to be applicable.  T h e Inviscid Region  2.4.1  From linearized potential theory (ref 1, pg.147), (2.21)  where UE,O {X) is the undisturbed velocity distribution and Uj (x) is the perturIN  N  bation to this distribution. In the present study only problems with U p ( ) — x  E  IN  Uoo , the free stream condition, are considered. The disturbance velocity, Uj (x) , is due to the presence of the b.l. and is given N  in incompressible linearized theory by, (2.22)  Now U (x) = U + U' ( ) therefore, X  ElN  00  IN  4-{U {x) ElN  6*(X)} = (Ulx + U^)  d6*_  dx  +  dU' 'IN dx  CHAPTER  2. ANEW  VISCOUS-INVISCID  PROCEDURE  30  d6*_  dx Substituting into equation 2.22,  -I  di  TT J-c  Substituting into equation 2.21,  E,A ) u„  U  X  =  i +  i  / - ° 5e t +0  0  X—  7T ./-oo  Introducing the non-dimensional variables, x  x h  =  —  UPOT{*) (2.23) U  P O T  (x) 6*{x)  A'(x)  where h is a convenient fixed length scale of the problem, equation 2.21 is nondimensionalized to give,  UPOTW  =  1 +  1  r +  co^{A*(Q}  IT J-  C/  Por  (2.24)  (x)  Equation 2.24 is in the required form for the semi-inverse procedure; it will be used to match the boundary layer edge velocity distribution given by the viscous equations.  CHAPTER  2.4.2  2. A NEW  VISCOUS-INVISCID  PROCEDURE  31  The Viscous Region  The governing equations for the viscous region are the classical boundary layer equations for steady flow found in Schlichting [15]. du x-momentum  u——V  dx continuity  du dU v— = UE — dy dx du dv — + — = 0 ax ay  1 dr 1— — p ay  E  ,  . (2.25)  . (2.26) '  The pressure term on the R.H.S of equation 2.25 has been replaced by utilizing Bernoulli's equation for steady flow, _\&V_ _ p dx  JJ E  dU dx  E  where UE is synonymous with the boundary layer edge velocity distribution, U  Evis  (x).  Equations 2.25 and 2.26 can be integrated across the boundary layer (in unseparated flow) thus, rr , du du JO ox s s  du du ay  dUn dU ax  T  E  W  p  and fv du , v = - / — dy Jo dx These are combined to give the well known momentum-integral equation for 2-D, incompressible boundary layers (ref. 15, pp. 158-160),  ^ i > +^ where,  roo  f = 7 ii  - = /„ p - ^ *  (2 27)  '  "  6  (2.28)  Jo UE  U  E  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  32  Several approximate methods have been developed to solve equation 2.27. Among the simplest of these is the method due to Th. von Kaxman and K. Pohlhausen (ref. 15, pp.206-222). Following Schlichting, equation 2.27 is solved by assuming a suitable form of the velocity profile, u(y). This profile must be capable of handling various pressure gradients, satisfying kinematic boundary conditions (including separation profiles) and should be relatively simple. Polynomial profiles of the form,  £  -  E  A (f)'  (0 < | < 1)  (2.29)  satisfy the above requirements if T is large enough. It is also assumed that for y > S, u = U . The boundary conditions for equation 2.29 are: E  at y = 6  u = U du  E  m  =  dy  at y — 0  0  m=l,2...T-3  u = 0 du  dU dx  2  dy  =  E  2  E  du z  dy->  =  °  (  T  >  4  )  Two values of T were chosen in the model development, T = 4,5. Using the transform variable, -  U.  Js^.  6  du _ ^ du dr) dy  and solving for the coefficients At of equation 2.29 from the boundary condtions, the following equation for ^ is obtained.  = UE  F {v) + A G {r,) T  T  (2.30)  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  33  where, 6 dU v dx 2  =  (2.31)  E  Here T = 4 for the fourth order profile and T = 5 for the fifth order profile. The functions FT and Gx are given as, = 2 (n - n + \n )  F,{n)  s  A  (2.32) G*{l)  = |(r?-2r  2 ?  + 2r -r ) 4  ?  5  ?  Solving for 6* and 5 from equations 2.28 using equations 2.30 and 2.32,  or, 6* — _ c  -\- C r  2  (2.33)  A  similarly, 0  TGS +  6  A + TCS A  2  where the coefficients TC\. . .jC^ are given as,  Ci  r  rC C  T  C  T  2  = + A l - F r f a ) ] dr? *o  = - f G (n) drj  Jo  T  3  = + [ F {r})[l-F {ri)]  4  =  l  Jo  T  T  + / G ( r ; ) [ l - 2 F r ( r ? ) ] dr; 1  r  JO  r  C  5  =  dr)  - f  1  Jo  G (n) dn 2  n  (2.34)  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  The following table gives the values of TC\ .. . C T  5  34  for the fourth and fifth order  profiles. Fourth Order Profile  Fifth Order Profile  (r=4)  (T = 5)  3/10  1/3  -1/120  -1/60  Cz  37/315  775/6237  TC\  -1/945  -25/16632  TCB  -1/9072  -47/110880  TC\ T  T  C  2  Next, the shear stress parameter, r /p, will be calculated with the aid of the equaw  tion, r  w  =  H  du dy  H du r?=0  Using equation 2.30 this can be written, v U P  E  dFr . dGr + A—— drj dt]  v U  E  [TOQ +  A]  (2.35)  n=0  The coefficients y C and rCV are tabulated below. 6  (r=4)  (T = 5)  CQ  2  5/3  Cn  1/6  1/4  T  T  Equation 2.27 can now be solved using the equations developed in this section. It is convienient to cast the equations into non-dimensional form. The following is a list of variables used to tranform the viscous equations into non-dimensional form.  CHAPTER  2. A NEW VISCOUS-INVISCID  PROCEDURE  35  Again h is the same convenient fixed length scale for the problem.  UE_  h 0  0  h  A  h x h  X  Equation 2.27 is to be solved on the computer, using a Runge-Kutta algorithm. The unknown variables U L, A , 0, and A are solved by the following system of equations B  derived from equations 2.27, 2.31, 2.33, 2.34 and the non-dimensionalizations given above. The coefficients iC\..  Cj have been replaced by C\... C7 and it is implied  that either fourth or fifth order profile coefficients can be used.  de dx dU dx BL  [C + ( C - C - 2 C ) A - ( C + 2C ) A - 2 C A ] RHUBLA 2  6  7  1  3  2  4  3  5  A RHA* (2.36)  dA dx  where, R  -  °°  kU  V  CHAPTER  2.  A NEW  VISCOUS-INVISCID  PROCEDURE  36  The algebraic relations which are implicit in equations 2.36 are given in equations 2.37 below. A* _  £  =  d  =  C  + C  s  + C  2  4  (2.37)  A A + C  5  A  2  The shear stress coefficient is found from equation 2.35.  (2.38)  Equations 2.36 can now be solved using the proposed semi-inverse method. These equations are parabolic, so that only initial values of the four unknowns are required. The values of UE  VIS  to values of UE  IN  obtained are matched at corresponding locations  (calculated from equation 2.24) and a solution is obtained.  The equations derived in this section are used for cases of unseparated flow. The extension to separated flow requires some modifications. Chapter 4 illustrates these modifications for a particular case of separated flow. The following chapter is an application using the equations in the form derived in this chapter.  Chapter 3  APPLICATION TO T H E BLASIUS  P R O B L E M  In the development of any complex model, testing is done in a stepwise fashion, from simple to complex cases. Thus the new VII theory is applied first to a simple unseparated flow problem: uniform flow past a semi-infinite flat plate in zero pressure gradient. H. Blasius first studied this flow at the turn of the century using Prandtl's then new boundary layer theory. The Blasius solution is an exact solution of Prandtl's boundary layer equations, but these equations are an asymptotic approximation of the Navier-Stokes equations and so the Blasius solution is in reality an approximate solution. In fact higher order b.l. theory (ref. 15, pp.145-148) shows that Prandtl's equations are first order approximations to the Navier-Stokes equations with e = l/y/R^ as a perturbation parameter. The proposed VII method is equivalent to the second order approximation with the restriction of zero wall curvature. The aim of this chapter is to test the VII method by finding the second order solution to the Blasius problem.  37  CHAPTER  3.1  3. APPLICATION  TO THE BLASIUS PROBLEM  38  Problem Geometry  Figure 3.1 shows the geometry of the Blasius problem with the displacement thickness curve, A*(x), divided into sections used by the procedure. The inviscid equation  Figure 3.1  Blasius Problem Geometry  (equation 2.24) requires that A* be defined from —oo < x < +00. With origin of the flat plate located at x = y = 0, the following definition of A* is devised,  A*(x) = -  0  x<0  J.(x)  0 < x < x„ X -1 < X < X n  Xjv < X  The region x„ < x < x  N  <  (3.1) n  OO  is termed the interaction region and is the region in  which the VII method is used. The functions /„(V>) are the spliced cubics derived in Chapter 2. The curves 7 (x) and J^(x) are fitted smoothly into each end of 0  the interaction region for calculation of the inviscid equation. The shape of these  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  39  curves is derived from the governing equations for the viscous region. A review of the literature shows that very few authors have meaningful solutions outside the interaction region, especially on the downstream side. Much attention was paid to the development of the curves £,(x) andjjs(x) to ensure a meaningful solution outside of the interaction region.  3.1.1  Boundary Conditions  The viscous equations are examined first to determine the initial values (at x = x ) 0  of the independent variables. A zero pressure gradient will be assumed at the start of the interaction region. From Bernoulli's equation this gives  = 0. With this  condition the viscous equations (equations 2.36) degenerate to, dQ  dx A^ A  Ci  0 A  C  >  3  A = 0 Taking derivatives of the second and third of equations 3.2,  A  dQ  dx  0  dA dx  Combining these with equations 3.2 to eliminate A and 0,  (3.2)  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  40  where, 2 CQ CJ C  2  3  The solution to this simple differential equation is, A*(x)  =  a  x-C  (3.4)  UBL Rh  where C is the constant of integration. Equation 3.4 is the non-interactive Pohlhausen approximation of A*(x) for the Blasius problem. Equation 3.4 will be used to generate the curves J (x) and 7 (x), as the second order solution should not be much 0  E  different from the Blasius solution. Additional mathematical reasons for doing this will be shown in the next section. With A*(x) = 0 at x = 0, equation 3.4 gives 7 (x) 0  as,  The external velocity, UE (X), will be assumed to be very nearly VIS  over the  entire length of T (x), so little error is introduced by setting 1*BL(X) = 1. 7 (x) then 0  0  becomes, (3.5) This equation gives the necessary information for calculating the initial values of the independent variables in terms of A*(x). From equations 3.1 and 3.5,  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  41  Substituting into equations 3.2, the initial values of the independent variables are, 1 A(x )  0  0  (3.6)  A(x ) 0  The next step is to calculate the curve ^(x).  Again it is assumed that out-  side the interaction region the solution may be approximated by the first order Pohlhausen solution. In the case of  7 (>Cj the origin was known (x = 0). 0  For complete  generality the origin of 7E(X) will be given a non-zero value, x', to be determined as part of the solution. This 'virtual' origin is useful for describing flows in which the b.l. separates, reattaches, then redevelops, such as the backward facing step problem of Chapter 4. In addition the boundary layer edge velocity will be assumed to be equal to the free stream value, Uoo at some large distance downstream of the interaction region. From equation 3.4 with C = x' and UBL = 1, (3.7) It should be pointed out that the initial and final spliced cubics must fit smoothly into the curves 7 (x) and 3*E(X) respectively. This constrains A* somewhat, but it 0  is necessary for a proper solution.  3.1.2  Asymptotic Limits  The value of UBL[*) at the end of the interaction region will in general be different from the freestream value of 1. One of the assumptions just made was that U (x) BL  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  42  returned to the freestream value at some large distance downstream. The arguments for the derivation of equation 3.7 are not proof that £^,(x) —• 1 as x—• oo. Appendix A contains a proof of this fact. This proof is valid for all problems in which the zero pressure gradient freestream is given a small disturbance. The interaction region may contain separated or unseparated flow; this does not affect the proof as the disturbance is assumed to be localized. All numerical codes which solve elliptical problems must have downstream conditions specified on boundaries at finite distances, but all exact solutions to the b.l. equations have asymptotic behavior for large x. In this respect the current method is more robust than these discrete methods. Next the elliptic inviscid equations will be examined to check asymptotic limits. The inviscid equation is given by  where A* is defined by equations 3.1.  The evaluation of UPOT{*) as x —> co is  somewhat involved so is moved to Appendix B, section B l . The result shows that both the viscous and inviscid solutions follow the trend of exact solutions to the b.l. equations. The solutions upstream of the interaction region are given for the viscous equations (ie: free stream conditions) but not for the inviscid equation. Both the viscous and inviscid equations are invalid very near the origin so analysis there would be meaningless. The beginning of the interaction region will be placed sufficiently far upstream from any disturbance so that solutions there are close to the freestream conditions (this is not a concern for the Blasius problem). Upstream of the interaction zone A*(x) = 7 (x). If the function 7 is chosen correctly, the conditions 0  0  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  43  upstream of the interaction zone will remain near to their freestream values. The form of 7 chosen, namely <*\Jj[^, has the desired properties. This particular func0  tion has another quite remarkable characteristic. If A*(x) = £,(x) for 0 < x < oo then UPOT{*) = ^BL(X) = 1 for all x ! In other words J is an exact solution to the 0  Blasius problem for the VII method. A proof of this claim follows:  Equation 3.4  is derived by setting A = A' = 0 in the b.l. equations. If the initial conditions are UBL = 1 and A = 0 then these conditions will not change when integrated along J because it was derived from the fact that A = A' = 0. The solution to UPOT{*) 0  involves some algebra and so is moved to Appendix B, section B2. Schlichting (ref. 1, pg. 147) shows that the solution of second order for the flat plate using l/y/R^ as a perturbation parameter gives the same result. This result provides a check on the accuracy of using spliced cubics in the interaction region to approximate the exact solution. With A*(x) given by equation 3.1, initial conditions of the b.l. equations given by equation 3.6, and asymptotic limits checked, the coefficients and functions necessary for the VII method can now be calculated.  3.2  Problem Formulation  3.2.1  Velocity Matching Locations  In section 2.3, A* was derived as a function of hi,h ,...,h , 2  N  these being the values  of A* at the endpoints of the spliced cubic polynomials. In order to solve for the h , U (x) and Up (x) must be matched at N different values of x. These values n  BL  OT  of x (termed x*) could be distributed anywhere along A*(x). The problem is set up so that strong coupling between the viscous and inviscid equations occurs only in  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  44  the interaction region, and so the x* will be confined to that region. It seems logical to locate one x per cubic and in the same relative location within each cubic, the k  former giving equal 'weighting' to each cubic and the latter minimizing necessary algebra. Grid refinement can be localized by choosing individual cubics of different lengths, any induced algebra being eliminated by use of the unity transformation of equation 2.11. The locations of the x are chosen to correspond with the endpoints k  of each cubic, for reasons which follow. Figure 3.2 shows the portion of A*(x) where the last cubic, f [ip), splices into 7E{*)- Two possible locations of x* within fN are N  X  Figure 3.2  _  *  I  Matching Point Locations  shown. Locating x at © ensures that the external velocities match at x = x , while k  N  this is not true for location © . For location © , the fact that the b.l. equations are integrated along a section of cubic after the last x , causes a 'kink' in the solution k  for A*(x). No kink was observed when x was moved to ® . Results of numerical k  experiments showed that the kink became smaller if (i) the last cubic was shortened, or (ii) the last x was moved closer to the end of the last cubic. When the last x is k  k  located at © , two constraints are imposed. The external velocities must match and  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  45  the last cubic must splice smoothly into ^E(X). Relaxing of the second constraint at ® allows the velocities to match at some value other than the desired one by locally distorting 6*(x). The remainder of this chapter is devoted to applying the linearized method developed in Chapter 2 to the Blasius problem.  3.2.2  Iteration Method  Prior to any derivations, an algorithm analogous to the one on page 15 is given below. This gives the reader a better view of what needs to be formulated, rather than trying to follow the actual derivations. 1) Guess  hi,...,h  N  2) Calculate a ,..., d from equations 2.14 n  n  3) Calculate Upori^k) and U i(x ) B  k  from equations 3.9, 3.10 and 3.11 (to be de-  rived). 4) If U (x ) P0T  = U {x ) for k = l,N  k  BL  Stop  k  5) Calculate the linearized perturbation functions AUPOT[XIC) and /\U i,(x ) and B  thus obtain the influence functions KUpor{^k)  P  and K U L{x ) B  k  k  given by equa-  p  tions 3.22, 3.23 and 3.25 ,and 3.26 (to be derived). 6) Calculate the perturbation parameters Ah from the (N) equations. p  N  (3.8)  = 52[{KUpoT{xt)p-KU {x ) } Ah ]  U {x )-U T{xk) BL  k  7) Update the h , ie: h p  8) Goto step 2)  BL  PO  Pnew  =h  PoU  + &h  p  k  p  p  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  46  The equation necessary for step 2) has been derived in Chapter 2 so that only steps 3) - 6) need development.  In the following section, equations for the boundary  layer edge velocities at the matching point locations are developed.  3.2.3  Velocity Calculations  Calculation of Upoxjxk) From equation 2.24,  l r  +co  UpoT^k) = 1 +  dA' IT  - / - ^ d TT J-oo Xk — Z  t  Using equations 3.1, 3.5, 3.7, 2.12, and the transform equation 2.11, the above equation is written as the sum of its constituent parts.  +  UPOT(U) = 1 + - f 7T JO  +  ~I  1  N  IT ^  1  r  AX„  1  b  n  Jo  + 2c„  e + 3d  n  ?  ipk-t  di  dt  x TT Jx N  (xjfc - $)  y/t-X  1  where, and  0jb =  x* —  X  n  - x  AX„  Integration of Upori^k) must be done with care. When £ = tpk a singularity arises in the integral. This only occurs in the interaction region (all x* are located there) so that the Cauchy principal part of the integral surrounding the singularity must be evaluated. Since the x* are located at the intersection between two cubics (except in the last case), the principal part becomes slightly more difficult to evaluate. The integration is done in Appendix C, section C l . The resulting calculations are: k 7^ n  CHAPTER  3. APPLICATION  TO THE BLASIUS  U OT{X>C) =  J {X-k) + JEM  1 +  P  PROBLEM + Ik,k+1  0  (3.9)  N  Jfc-1  n=l  47  n=*+2  k= n N-l IN,N+I  (3.10)  n=l  All the functions  (J ,  J  0  N  ,J  E  , I , +iy k  k  and  IN,N+I)  a  r  e  given in Appendix C, section  CI.  U (x )  Calculation of  BL  k  In the interaction region A* consists of the spliced cubics f {ip)- Thus it is convenient n  to solve for UBL in i/>-space instead of x-space. The boundary layer equations 2.36 can be transformed into V'-space by use of the equation, d dx  (x„_! < x < x ) n  Ax dtp n  where it is understood that all variables in the transformed set of equations are functions of ip. From the second of the transformed equations is an expression for the derivative of U L>  s  o  B  that, £ J r dU 1  BL  dip where U \(tp ) = U (x ). BL  k  BL  Substituting for  k  at the cubic endpoints (rp = l), k  J  (  ,  rt* dU  BL  Jo  dtp  and noting all matching points are  CHAPTER  3. APPLICATION  TO THE BLASIUS  Equation 3.11 is symbolic as U  BL  PROBLEM  48  is solved as part of the system of b.l. equations.  Step 3) of the outline on page 45 is now completed.  3.2.4  Perturbation Calculations  Perturbations to 5* (A6*) yield perturbations to U (&U ) for each region. Before E  E  any calculations of &U can be made, the perturbation curve for the case p = N B  must be developed, as discussed on page 27. Figure 3.3 shows the perturbation curves when h  N  is given a small displacement Ah . N  Figure 3.3  A/ (V0 W  =  A  a  N  + A&N V* +  A  C  The boundary conditions for  Perturbation Curve for n = N  N V" + A d tp are, 2  3  N  A/*(0)  0 0  A/„(l)  CHAPTER Solving for  3. APPLICATION  TO THE BLASIUS  PROBLEM  49  . . . Ad ,  ACL  N  N  Aa  N  Ab  N  AC  N  Ad  N  -  0  =  0  =  3Ah  =  (3.12) -  N  AX  ^ A J  N  -2Ah + A X N  £  ( X ) |  X W  ^A7 (X) \XN £  N  The curve downstream of the interaction region is altered when h  N  is perturbed  This altered curve (denoted ^ B ( X ) in figure 3.3) must be of the form £ B ( X ) = —— if the asymptotic conditions are to hold after the perturbation. Letting C =  x  a1  g {x)  =  E  a  (3.13)  R  h  When the curve 7E{X)  1  S  perturbed, its linear perturbation curve  alent to the linearized difference between £ . E ( X ) and rameter in 7E{ ) is ' X  x  s  o  JE(X).  ATE{*)  The only variable pa-  that it becomes the perturbation parameter, ie:  7 (x) + A7 {X)  =  e  E  a  (x' +  JX -  AX')  Rh X —  AX'  X'  a  x — x'  RH  Using the binomial theorem, /  AX  AX'  '1  x  :  —  X'  »  2(x - x')  This gives the linearized perturbation function as, a AJ  £  (X)  AX  =  2^]R  h  (X-X')  is equiv-  CHAPTER  3. APPLICATION  TO THE BLASIUS  The term A I ' can be written in terms of Ah  N  Ah ,  A ^ X , , ) =  N  PROBLEM  50  by using the boundary condition  ie:  (3.14)  Because  AJE(X)  is the linearized difference between §E(*) and 7E(X), a linear esti-  mate of x' comes from the equation, G  J (x ) + A/I„ = £  N  9E(* ) N  or, (3.15)  In step 7) of the outline on page 45, x^ is the update to x', to be used in the subsequent iteration. Finally, using equation 3.14 to solve equations 3.13, Aa A6  N  n  =  0  =  0 (3.16)  AC„  = (3-- Cl) Ah  Ad  =  N  N  (n  - 2) Ah  N  where, AX,  2  (3.17)  ( x „ - x ' )  Following along similar lines to the influence coefficient calculations of Chapter 2, the perturbation coefficients and functions corresponding to Ah are, N  Kc Kd  N  =  3- n  N  = 0 —2  (3.18)  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  51  with Kf„{1>)  =  Kc ^  + Kd ip  2  N  s  N  and (3.19) Having now derived the case p = N, calculations for steps 5) and 6) can be completed.  Calculation of the Inviscid Velocity Perturbation &UPQT[XIC) To obtain AUporix-k), equation 2.24 for UPOT{*IC) is perturbed, ie:  TT  f \±  1 -u  TT ( \  UpoTKXk) + &U (x ) POT  k  1  r°° ^ [ ^ ( 0 + A A ' ( Q ]  = 1 + - / 7T  -s Xfc —  J-oo  Q  Separating out the perturbation terms and noting that the perturbation curve AA* is zero if x < x , 0  AUPOTM  =  - / ~ * J*o  Xfc -  — di i  Again AA* can be divided into its constituent parts. Also remembering that that the perturbation function due to a single A/I is comprised only of the two cubics surrounding the perturbation,  AUPOTMP  =  " ^ 7  -22—1 7T  n  =  p  AX„  JO  1pk-  JZ  dC + ~ / 7T Jx„  T — d£ Xfc - C  3  3.20  where A£/por(2;fc)p is the perturbation velocity at x = x* due to perturbation of a single h, namely h at x = x . For clarity, the list below describes the meaning of p  §  This term is zero if p ^ n  p  CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM  52  each index. n  Index of cubic function Index of velocity matching point location Index of perturbation location  The perturbation function AUPOT i k)i due to all perturbed /i's is given by the sum x  of the individual perturbations, ie: N  N  P=I  p=i  (3.21)  Analogous to the calculation of Upori^k), singularities are carefully avoided by considering the principal value of the sum of the two integrals surrounding the singularity. Again there are two main cases to consider: [l] p -fi N, [2] p = N. These are different from the cases for Upori^k) as there the case was dependent on the value of k (the velocity matching location) while here the case depends on the value of p (the perturbation location). For case [l] there are four subcases to consider, each depending on the relative location of X* to x . In each case it is p  important to note that the only contribution to the integral for At/por(xjt) comes from the two cubic perturbation curves surrounding the h . The perturbation cubics p  were carefully designed so that all other perturbation cubics were nonexistent and so do not contribute to A£/por(x*)Figure 3.4 shows the locations of the four subcases labelled ® , © , ® , and ® . Calculation of KU T(^k)p involves some algebra and so is moved to Appendix C, PO  section C2. The appendix shows fairly detailed calculations of the Cauchy principal part of the particular integral under consideration.  The resulting formulae are  CHAPTER  3. APPLICATION  TO THE BLASIUS  Figure 3.4  PROBLEM  53  p ^ JV Subcase Locations  summarized in equation 3.22.  i [ A J ( 0 ) + AJ p  t  p+1  (0 )]  case (T)  t  case (?)  KUpoT{xk) — P  AX,  AX  Vf+i  m +  1  _  \Kd +i+2Kc r  AX  P  +  AXp+i  v  case (3) case (T)  1  (3.22) The functions  A J  P  and A  J  R+X  are given in Appendix C, section C2.  For case [2] where p = JV there are three subcases to consider, as shown in figure 3.5. Again the calculations for KUporixk)* are fairly complex and so are moved to Appendix C section C2. It is noted that the complexity is due to the fact that the matching locations are located at the end of every cubic, making the Cauchy principal part of the singular integrals more difficult. The resulting calculations of are given in equation 3.23 below.  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  54  1  1  1  ©  O *  *  *  1  ^  X„  Figure 3.5  p — N Subcase Locations p= N  *«MV»*) 1 J[-  + &J {*k)]  AJ (x J s  t  B  * —  case © case  (2)  (3.23)  case ©  J[A J ^ X * ) ]  The functions J (ip ), AJ (x ), and A J ^ ( x ) are given in Appendix C, section C2. N  k  B  k  t  This leaves one final perturbation calculation, A L ^ i ( x i ) . Calculation of the Viscous Velocity Perturbation, &UBL,{*k) From the previous viscous velocity calculation, dU  BL  drp  A R A A  X  N  2  X „ - i < x < X„  h  The R.H.S of this expression needs to be written in terms of A* as it is the function whose perturbation is known. Substituting for A from the first of equations 2.37,  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  55  combining these two equations,  B  _  =  dxjj  C R A 2  1  „  [A*  AX„  dU L  2  h  This equation is then perturbed giving, d{U  AU )  +  BL  A*  AX  R  BL  di>  AA*  +  .A  AA  +  C R (A + A A ) 2  2  h  Because the perturbation A A is small compared with A, the expression can be written, A* d{U  + AU )  BL  A  BL  AA*  +  »  X  dtp  C R 2  A  h  2  By this method of linearization, the unknown function A A is eliminated. Subtracting U L from both sides gives the perturbation, B  d  AX  1J * {  Ubl)  =  AA*  N  X n - l < X < X„  cTihA-*  (3.24)  For a single perturbation A/I„ at x = x„,  A* =  A  Xp-i  {  KfP+IM A/ip x  where Kf [xp) and Kf \(ijj) p  p+  < x<  x < x< x p  X„  p+1  are calculated in section 2.3. Because the boundary  layer equations are parabolic, any perturbations downstream of x* will not influence U L (xjt) so from equation 3.24 and the above equation, B  0  k<p-l AXp  KU (x ) BL  k  P  f l  (3.25)  C R Jo  = '  2  g k  h  AXp  n=p C R Jo 2  v h  A  3  k>p  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  56  Analogous to the inviscid example the perturbation function AUp (x;t) is the sum or  of all the individual perturbations, or N  (3.26) p=l  Having completed the formulae necessary for steps 5) and 6) of the guide on page 45, equation 3.8 can now be solved to yield the correction parameters A/I . P  Although the algebra necessary for the derivations is somewhat lengthy, the results are fairly general, as will be shown in the next chapter where an application to separated flow will be examined. Results of this analysis are shown in the following section. A discussion of the limitations and stability of the method is included. Initial conditions for A*, grid refinement, and run times are also analysed.  3.3  Results and Analysis  From the analysis in section 3.1.2 it is shown that U L{*) = ^por(x) = 1 is an exact B  solution using the present VISI method for the case A*(x) = ? (x) for 0 < x < oo. 0  This offers a convenient way of checking the accuracy of the present method. In essence the method approximates 7 {x) using spliced cubics. 0  3.3.1  R a d i u s of Convergence and Stability  A complex mathematical analysis is required to analytically determine the radius of convergence and the stability of an iterative scheme. An analysis of this nature is beyond the scope of this study, however numerical experiments were performed to simulate the results. The radius of convergence determined from these experiments  CHAPTER  3.  APPLICATION  TO THE BLASIUS  57  PROBLEM  was fairly large. Any reasonable guess for A* (x) converged smoothly and rapidly to the solution for all cases considered in the Blasius application. Critical points in the viscous equations are found in the third of equations 2.36. Utilizing the relations 2.37, dA dx  when  A  For T = 4 this corresponds to A = 12,60 and for T = 5, A = f, *5Q. it i only s  possible to have an initial guess outside the larger of the critical values if A*/A is negative.  The smaller critical value represents the limit of the profile shape for  steady flow (U/UE becomes > 1 for A greater than this value). It was determined through numerous test cases that a very unusual guess was required for the solution to diverge. In many of these cases divergence occurred for one or two iterations followed by full convergence. Figure 3.6 shows the iteration history for a poor initial guess. The symbols mark the locations of the cubic endpoints. In this case the domain is divided into 10 cubics of equal length. The x-axis has the dimensions of length. For example with A* measured in mm and the fluid being air at 20° C and a freestream speed of 10 m/s, the domain is from x = .68 m to x — 1.36 m. The solid line through the endpoints for iteration zero is the initial shape of A*(x) produced by the spliced cubics. The smooth line is the exact solution discussed in section 3.3. Convergence is rapid and monotonic even though the intial guess is wavy.  CHAPTER  3. APPLICATION  TO THE BLASIUS  Figure 3.6  PROBLEM  Iteration History  CHAPTER  3.3.2  3.  APPLICATION  TO THE BLASIUS  59  PROBLEM  Solution Accuracy, G r i d Refinement, and Execution Speed  An exact solution provides a precise measure of the accuracy of the present method. The solutions obtained were so near to the exact solution that special plotting schemes are required to distinguish the order of accuracy for a given solution. The exact solution yields values of x' = 0 and UBL — 1- The deviation of x and U~BL from these exact values is used to measure the accuracy of a particular solution. Figure 3.7 shows how these values are used to assess grid refinement. Grid refinement was performed by successively dividing the interaction region into 2 equal length N  cubics. The accuracy of the method is quite good. The value of x' is very sensitive to the value of A* at x . Despite this, and using only 4 cubics, the value found for N  x' is less than 1. This is surprising considering that x' is found by extrapolation from x at x  N  N  (1000 in this case). For example, a 1% reduction in the solution of A*  causes x' to increase by a factor of twenty. Increasing the number of cubics  increases the accuracy of solution as expected, since the shorter cubics can simulate the exact solution curve with greater accuracy. This trend is also reflected in the parameter (l — UBL)MAX, which measures the maximum deviation of UBL from the exact solution. This parameter shows the velocity is accurate to nearly 1 part in 1000 with only two cubics. The amount of CPU time required for solution is dependent on the number of cubics, the number of iterations required to attain the desired accuracy, and step size used in the integration routine for the inviscid equations. The example shown in figure 3.6 (N = 10) converged to [&h/h]  MAX  < 10  -3  in 5 iterations and 0.293  seconds on the UBC Computings Center's Amdahl 5850 mainframe computer.  CHAPTER  3. APPLICATION  TO THE BLASIUS  (x =iq 0  10"  x =iooo N  PROBLEM  R =ioo) h  —i  io-H  io"  H  io'  H  io" H  +  (i-uBL)t io"  -H  io-  H  +  + +  16  N  Figure 3.7  Blasius Grid Refinement  32  CHAPTER  3. APPLICATION  TO THE BLASIUS  PROBLEM  61  It is noted that the program has not been optimized for speed and could conceivably become significantly faster. For this particular case the new VISI method worked very well, but a more useful application is given in the following chapter.  Chapter 4  APPLICATION B A C K W A R D  TO A  FACING  STEP  In the previous chapter the VII method was applied to a simple situation of unseparated flow. This was done to ensure that the basic components of the model were functioning correctly. The use of the VII model, however, as a solution procedure for the Blasius problem is unnecessary. Many effective non-interactive techniques, such as the classical Pohlhausen method, give equally accurate results with very much less computational effort. A more useful application of the VII method is one containing a region of separated flow. The reasons for selecting the Backward Facing Step (also known as a 'Downstream Facing Step') as a model geometry are listed on page 4. The viscous equations developed in Chapter 2 are modified for the region of separated flow. A simple modification of Pohlhausen's method is developed for this purpose. As discussed in Chapter 1, the VII method is applied to the Backward Facing Step (BFS) in steady, 2-D, laminar, incompressible flow.  62  CHAPTER  4.1  4. APPLICATION  TO A BACKWARD  63  FACING STEP  Problem Geometry  I  1 A L  r  L.  1 :m i i  L„  -j  x  n  x„  X  Figure 4.1  B F S - Problem Geometry  Figure 4.1 shows the geometry used in modelling the B F S using the VII method. Dimensions in the y-direction have been greatly exaggerated to make clear the different regions and the parameter definitions.  Region I is the inviscid region.  This region is comprised of uniform flow disturbed by the presence of a boundary layer growing over the step. Regions II and III are the viscous regions, with III being the region of recirculating flow. Region III is bounded by the 'separation streamline', a line through which no fluid passes. The parameters used to describe the geometry are listed below. A l l quantities are dimensionless, having been scaled with respect to the step height, h. x  Distance from origin of b.l. to beginning of interaction region.  0  x  N  x  s  =>  Distance from origin of b.l. to end of interaction region.  =>  Distance from origin of b.l. to step.  CHAPTER x L  4. APPLICATION  =>•  R  TO A BACKWARD  FACING  STEP  64  Distance from origin of b.l. to reattachment. Reattachment length.  R  A,(x)  = r * Height of recirculation region.  A (x)  =>•  c  Distance from x-axis to outer edge of b.l.  A*(x)  Displacement thickness.  A(x)  Boundary layer thickness.  Ai(x)  Displacement line.  =>  A = A — A, 0  A = A , + A* x  The effect of the size of the interaction region is not known a priori, but can be determined through variation of the parameters x and x . 0  4.1.1  N  T h e Displacement Line  In the Blasius problem A* is the measure of streamline displacement due to the viscous region. In the BFS problem it is A i (not A*) which gives the streamline displacement. Figure 4.2 shows Ai(x) subdivided into components, analogous to the Blasius problem. Ai(x) is described by the same functions as in equation 3.1, except for 7 (x). Because the upper wall is displaced a unit distance from the x-axis, 0  the function 7 (x) must also be displaced a unit distance, ie: 0  1  x < 0 = 1+Q  A (x) = x  \/iE  0 < x< x  0  = o„ + b ip + c„V + d„0 2  n  = a  Ix - x'  V *  3  x  n-l ^  X  N  x  < n x  < X < OO  (4.1)  CHAPTER  4. APPLICATION  TO A BACKWARD  —  FACING  STEP  65  • — ^ —  0  X  X  1  2 ' / / / / X  Figure 4.2  n  X  N  The Displacement Line  All asymptotic properties calculated in Chapter 3 hold for the BFS problem, as shown in Appendices A and B.  4.2 4.2.1  Modifications to the Governing Equations The Inviscid Equation  The perturbation U (x) of equation 2.21 can be thought of as the perturbation to N  uniform flow caused by the presence of a boundary layer flowing over a step. With A i being the amount the streamlines are displaced for the BFS problem, the same development which led to equation 2.24 is valid here. From equation 2.24 with A i  CHAPTER  4. APPLICATION  TO A BACKWARD  FACING  STEP  66  replacing A*, U {x) POT  =  1+  1 f+oo^AxU)}  ix J-  (4.2)  U (x) P0T  4.2.2  T h e Viscous Equations  Figure 4.3 shows the boundary layer divided along x into three separate regions. Regions V I and V3 consist entirely of unseparated flow and so are governed by  Figure 4.3  Viscous Regions  the same set of equations 2.36 as in the Blasius problem. In region V3, A* in equations 2.36 is replaced by A i and in region V I , A* is replaced by A i — 1. In region V2, A * cannot be replaced by A j — A , because equations 2.36 were derived by taking derivatives and A , is not a constant. The equations for region V2 will be rederived to account for the variablity of A , . When the integration passes from one region into the next, the equations sets must be switched. A discussion of how the switch is implemented is given after modification of the equations in the region  CHAPTER 4. APPLICATION  TO A BACKWARD FACING STEP  67  of separated flow. In region V2 the ' F L A R E ' approximation derived by Reyhner et al [16] is implemented. This is done to permit forward marching of the boundary layer solution in a region which contains backflow. Using this approximation, the z-momentum equation is written as, dU  1 dr  dx  p dy  E  where C = 1 for u > 0 and C = 0 for u < 0. In the viscous model used here the ' F L A R E ' approximation is taken one step further. It makes no less sense to neglect all convective terms of the momentum equation in the reverse flow region, for the reason that if u is zero everywhere from continuity v is also zero. The integral of mass flow across the recirculation region is zero and so with the small velocities in that region, the integral of momentum should also be quite small. With the region of reverse flow being replaced with essentially 'dead' air, the zero velocity line becomes identical with the separation streamline. The separated shear layer is then treated as a boundary layer flowing over this region of dead air, and a modified form of the Pohlhausen integral method used in Chapter 3 is used to solve the equations of motion. Figure 4.4 shows a typical profile in the separation region with the dead air approximations incorporated. The velocity profile shown is obviously in error but since the velocities in the recirculation region for the BFS have been shown to be small, the results of the approximation should not have significant effects. The assumptions will be checked by comparison to other more sophisticated and validated methods. In his book on incompressible flow, Curie [17] points out that the Pohlhausen method is not very accurate in regions where the pressure gradient is adverse. Using fourth order profiles, errors in the separation location on a fiat plate in linearly  CHAPTER  4.  APPLICATION  TO A BACKWARD  Figure 4A  FACING  STEP  68  Separation Profiles  retarded flow of up to 30% compared with exact solutions are mentioned. This error decreases to 20% when fifth order profiles are used. In favorable pressure gradients, however, the fourth order profiles are more accurate. No single Pohlhausen profile is accurate over a wide range of pressure gradients, and so two profiles are chosen to determine the effects. If the differences in solution are significant, then this is sufficient evidence to justify the effort of using a more sophisticated b.l. solution technique.  Modified Pohlhausen Method for the Separation Region For region V2 (see figure 4.3) both fourth and fifth order velocity profiles will again be used for the separated shear layer, with n = y/6 repaced by r) = [y — 6,)l\6 — 6,} 0  (see section 2.4.2). This change to rj does not affect the succeeding derivations in section 2.4.2 so that the momentum-integral equation can be written (in non-  CHAPTER  4. APPLICATION  TO A BACKWARD  FACING  STEP  69  dimensional form, in 0-space), A X  with,  n  [ C  6  +  C A] 7  (4.3)  A* A  0  ^  =  C3 + C 4 A+  RH A  A =  2  dU  C5A  (4.4)  2  BL  dtp  and where the coefficients Cy ... C^ retain their previous values. In the 'dead air' (reverse flow) region, with u = v = 0, the integrated xmomentum equation becomes, Si  Jo  dr  p Jo Jo dy  dx  dy  Integrating, TT  D  U  E  X  ax  1  I  \  (4.5)  p  where T is the shear stress at the wall and T> is the shear stress at y = 6,. As a W  further (but consistent) assumption, the shear stress at the wall will be neglected. The assumption that u is small in the reverse flow region does not imply that is small there as well, although  is zero at reattachment. If this rather drastic  assumption causes significant inaccuracy then some model for r will have to be w  developed. Writing the laminar form for r , 7  du dy  T = p. {  ix du 6 dn  fx U  E  n=o  [C  6  +  C A] 7  CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP  70  Substituting into equation 4.5, non-dimensionalizing, and transforming into ipspace,  du _ dip ' BL  Ai  _ ^ Rh A  =  [  C  6 + C  t  A  _  ]  (4  6)  with, A, -  Ax-A*  (4.7)  Combining equations 4.3, 4.4, 4.6, and 4.7, the following system of equations suited to the use of the Runge-Kutta solver is obtained.  dQ dip dU dip BL  dA dip  AX„  [C + (Cf - C - 2 C ) A —(C + 2C ) A - 2 C A ] 2  6  :  3  2  4  3  5  RHUBL' A X  n  A  Rh A  2  dAi dip A c  +  / A A d0 Ve) dip  Ce  ^  (4.8)  (C + 2C A) 4  5  A ~ 2  dA dip  A 0  dQ  .._  - - A ( C  4  +  ^„ ,. dA  2C A)5  and where equations 4.4 give the algebraic relations between A*, A , 0, and A, and equation 4.7 gives A , . Switching Equation Sets There are two locations where equations sets are switched: [l] At the step, where equations for region VI are switched to equations for region V2. [2] At reattachment, where equations for region V2 are switched to equations for region V3. For the first switching location a condition on both sets of equations is that  CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP  71  A , = 1 (see figure 4.1). From equations 4.6 and the second of equations 4.8, A, = - A  rc6  + c,  (4.9)  so that for A , = 1, A =  -  (4.10)  . C + C AJ 6  7  Now, because A i is prescribed before any calculations begin in order to obtain the influence coefficients, it must not be altered by the equation set switch. With A* = Ai — 1 for region VI and A* = A^ — A , for region V2, it becomes clear that A* also cannot change at the switch point as A, = 1 here. Subscripting the variables at the end of region V I with the letter A, and those at the beginning of region V2 with the letter B, A ; = A ; = A a f d + CaAa] Substituting for A from equation 4.10, B  -A,  A :  C\ + C A 2  . C + Cy A 6  t  £  Solving for A in terms of (known) A * , B  \C C A* ) 2 Co 1+  Here the  7  A  4C,  (4.11)  2  has been chosen because from Bernoulli's equation it is seen that  the pressure gradient is positive across the step and therefore A is negative. Equation 4.10 is then used to calculate A using the now known A . B  B  CHAPTER  4. APPLICATION  TO A BACKWARD  FACING  STEP  72  Also, from the second of equations 4.4, 0  B  = A  [C + C A + C A ] 2  B  s  4  fl  5  (4.13)  Equations 4.11-4.13 represent the fact that A, A , and 0 are all discontinuous across the step. Experiments and other sophisticated codes show that A is nearly discontinuous at the step. No measurements of A or 0 for the free shear layer were found in the literature so comparisons for these two variables is difficult. The T E A C H code in Chapter 5 will be used to estimate the free shear layer for comparison purposes. The second location for switching of equations occurs at the reattachment point (x = x ). Analysis of equations 4.8 shows that A' becomes infinite at this point. R  This is caused by the dead air velocity profile and is not a general property of the integrated boundary layer equations. The singularity is of minor consequence as it can be bypassed with a slight loss of accuracy. Figure 4.5 shows diagramatically how this is accomplished. The equations of region V2 are integrated very close to reattachment, where A has a sufficiently small value, e. At this point the viscous 7  equations are switched to those of region V3. Remembering that A i cannot be altered, and using subscript C to denote the variables of region V2 just before the switch, and D for the variables of region V3 just after the switch,  or  (4.14)  If the switching point is sufficiently close to reattachment, negligible error will arise from setting A to its reattachment value. With reattachment located by the point D  CHAPTER 4. APPLICATION  TO A BACKWARD  FACING STEP  73  A,  ^  ^  ^  ^  A,  tJ/)))}/}/)}}))  Figure 4.5  X  t  x „  Reattachment Singularity  of zero shear stress, equation 2.38 yields, A  -  C  e  (4.15)  Substituting equations 4.14 and 4.15 into the first of equations 4.4, (4.16)  From equations 4.15 and 4.16, and the second of equations 4.4,  (4.17)  At both switching locations some of the variables contain discontinuities which are a cause of concern. The discontinuities at the second switching location can be  CHAPTER 4. APPLICATION  TO A BACKWARD FACING STEP  74  made negligible by careful integration. No such procedure will help at the first switching location. This is not too distressing however, when it is realized that the underlying assumptions for boundary layer theory are not valid at either switching location. The effect of the discontinuities on the solution can only be compared to a model which does not use boundary layer theory as a basis. This is done in Chapter 6, where the present model is compared to the T E A C H code used in Chapter 5. Now that the governing equations have been formulated, the linearized perturbation method is applied to produce the necessary functions and coefficients.  4.3 4.3.1  Parameter Calculations Velocity Calculations  Calculation of UpoT{*k) The inviscid calculations do not change as equation 4.2 is identical to equation 2.24 (with A i replacing A*) so that all prior derivations hold. For the backstep problem Upor(xk)  is given by equations 3.9 or 3.10 which were developed for the Blasius  problem. It is worth noting that Upori^k) is given by these same equations for any problem in which linearized potential theory is applicable and the far field conditions are those of uniform flow.  Calculation of UaLJXk) Analogous to equation 3.11, the symbolic equation for U i is written, B  CHAPTER  4. APPLICATION  TO A BACKWARD  FACING  STEP  75  The viscous equations are integrated parabolically to x*, with the equation set switching procedure developed in section 4.2.2 invoked if integration crosses into regions V2 and/or V3.  4.3.2  Velocity Perturbation Calculations  Calculation of AUpoT&k) Again the derivations are identical to those of Chapter 3 and so AUpori^k)  is  calculated from equation 3.21 and 3.22 or 3.23.  Calculation of A UBL (**) The perturbation function is calculated for each of the three regions V1-V3 in figure 4.3.  For region V3 the derivation is identical to the Blasius case and so  AL Bi (xfc) is given by equations 3.25 and 3.26. For region V I , A* is replaced by 7  /  A — 1 so that because 1 is a constant the perturbation does not change, and again x  equations 3.25 and 3.26 can be used. The perturbation function in region V2 is not the same because A* = Ay — A, here and A, is not a constant. Some extra algebra is required to get the perturbation function AUBL in terms of A A i only. From the second of equations 4.8, dU _ d%l> ~ BL  perturbing this equation,  AX A RA n  2  h  CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP  76  Extracting the perturbation, I  AX Rh A  TT \  N  dip  2  (4.18)  AA  The next step is to find AA in terms of A A J . Combining equations 4.6 and 4.7 with the first of equations 4.8, Ai"  C A + C i — Cj — 2  A-C« = 0  2  Perturbing, A  C ( A + A A ) + C\ — Cj — 2  2  X  + AAX  ( A + AA) = 0  A + AA  Neglecting higher order terms and solving for AA in terms of A A A  AA =  1 5  AAI  Ax  2C A + C\ — Cj 2  but, Ax  C2 + —= A  2C? A + C\ — Cn  A  A J 2  Substituting the last two equations into equation 4.18,  d  AX  V 7 ( A ^ )  =  R  (4.19)  AAX  c + 2  AJ 2  Finally, equation 3.25 from the Blasius problem can be modified to give,  k<p-l AX  KU {x ) BL  k  f  P  C Rh 2  p+i  A  k=  !  X(A) Kf (iP) n  n=p 0 Rh JO 2  dip  P  k>p  (4.20)  CHAPTER  4. APPLICATION  TO A BACKWARD  FACING  STEP  77  where 1 1+  x <x<x 0  C  1  s  or  x <x<x R  N  n -1  &  C A J 2  x <x<x s  R  2  After computing the influence coefficients, KUsLi^k^pi from equation 4.20, equation 3.26 is used to solve for A U B ^ X * ) . The guide on page 45 is again used along with the parameters derived in this chapter to solve the backstep problem. Parametric variation, initial configurations for Ai(x), and grid refinement are analyzed in Chapter 6. A discussion of the assumptions and limitations of the model are also found in Chapter 6, along with comparisons to experiment and the T E A C H solution from Chapter 5.  Chapter 5  TEACH  VII methods are used to calculate approximate solutions to the Navier-Stokes equations. More exact solutions are found through numerical integration of these equations. In this chapter a numerical code (TEACH-T) is applied to the case of a BFS in steady, incompressible, laminar flow. A short discussion of the methods used by the code is included. A detailed description of the methods is found in a text by Patankar [18]. It should be noted that the T E A C H - T code was developed for internal flow and is a fully elliptic method. More efficient parabolic-elliptic external flow solvers exist, but they are not considered in this study.  5.1  Computational Procedure  For steady, 2-D, laminar flow the continuity equation and the components of the Navier-Stokes momentum equations can be written in the following general form.  78  CHAPTER 5. TEACH  79  T is a general diffusion coefficient and £^ is a general source term. The values of <f>  t  T and  for the three governing equations are given in the following table.  CONSERVATION EQUATION r Mass  l  x-momentum  u  y-momentum  V  s*  0  0 dp dx  ,d u dx> 2  +  fl{  dv dxdy 2  +  du dxdy 2  5.1.1  dy  Control Volume Formulation  it-  Figure 5.1  Scalar Cell  The principle of the T E A C H - T code is to divide the computational domain into a discrete number of 'cells' or control volumes. The control volumes for the scalar variables (p, p) are staggered from those of u and v so that the velocities are  CHAPTER 5.  TEACH  80  evaluated on the boundaries of the scalar cells. The scalar cells are located such that any in contact with a physical boundary will have its face(s) coincident with the boundary. A typical control volume is shown in figure 5.1 with its four adjacent neighbors (N,S,E,W), its four faces (n,s,e,w) and its associated u and v cells. Equation 5.1 can be integrated over the control volume and then transformed using Gauss's divergence theorem to yield, / / Sfdxdy J Jcv  5.1.2  H y b r i d Differencing and Linearization  Following Patankar, equation 5.2 is discretized using hybrid differencing, a combination of central and upwind differencing. The central method assumes a linear variation of (f> between cells while the upwind scheme assumes a step distribution of <j> whose value depends on the sign of the convective velocity. This hybrid scheme is necessary for higher Reynolds number convection-dominated flows. Central differencing alone leads to numerical instabilities when the Peclet number, P = pu AX/T T  is greater than 2. The source term is linearized using central differencing to yield, (5.3) Using equation 5.3 and the hybrid differencing scheme, equation 5.2 can be written,  (a — S ) (j) = a 4> + a <f> + a (f) + a 4> + S p  p  p  N  N  s  s  B  where, a  p  = d  N  + a + a + a s  B  w  E  w  w  u  (5.4)  CHAPTER  5.  TEACH  81  and for example, = max(|C /2|, D ) n  n  C /2 n  where, dx  D.-J  Xy,  Jx Central differencing is always used for the diffusive terms (ie. D ) and upwind n  differencing is used in the convective terms (ie. C ) if \P \ > 2. n  5.1.3  e  Solution Procedure  The system of equations 5.4 is solved by a method which calculates the velocities from the momentum equation based on the pressure field from the previous iteration. The pressure and velocity fields are then corrected using an approximate pressure correction equation. The corrected pressure is used to calculate the velocities again and the process repeats until continuity of mass is satisfied. The system of equations 5.4 is solved by an efficient tri-diagonal matrix solver after each pressure correction.  5.2  Problem Geometry  Figure 5.2 shows an unsealed schematic of the computational domain and its associated boundary conditions. The conditions are imposed so as to simulate a boundary layerflowingover a small step in external flow. The boundary conditions on the lower walls and step side wall are those of a non-permeable no-slip surface. The inlet conditions are the same as those for the VII problem at the start of its computational domain (ie. at x = x ). The velocity profile at the inlet is that of 0  CHAPTER  5.  TEACH  Figure 5.2  82  Boundary Conditions  a zero pressure gradient Pohlhausen profile x /h stepheights from the origin of the 0  b.l. The top surface is assumed to be far enough away from the small disturbance at the step that the conditions have returned to those of the free stream. Abdullah's [19] zero pressure gradient condition is implemented here. This condition causes the flow at the top wall not to accelerate as it moves downsteam, thus simulating the free stream. The conditions far downstream of the step should be those of a boundary layer in zero pressure gradient. This condition is difficult to simulate and so is replaced by the condition of zero velocity gradient. Use of a large solution domain means the boundary conditions are more accurate but unfortunately the number of cells needed becomes prohibitively large. For the majority of results presented in this study, the size of the computational domain and the corresponding scalar cell arrangement are shown in figure 5.3. The grid has 37 x 27 active cells with an outer layer of fictitious cells (not shown) used to implement the boundary conditions. The size of the domain was chosen in reference  CHAPTER  5. TEACH  83  Figure 5.3  Grid Layout  to the work of Kwon and Pletcher [7] and tested in the study discussed section 5.2.2.  5.2.1  Special Cells  Staggering of the u and v cells with respect to the scalar cells causes some difficulty at the corner of the step.  Following Djilali [20], the u and v cells which have  half-faces in contact with a wall require special treatment. Calculation of fluxes corresponding to these half-faces must use half the face area and a normal velocity equal to that at the outer edge of the half-face.  5.2.2  Previous W o r k and Present Refinements  Prior to using the geometry and boundary conditions discussed in section 5.2, a detailed study was conducted with a slightly different geometry and a different  CHAPTER  5.  TEACH  84  boundary condition on the top surface. The inlet height was 15h and the top surface was an impermeable free-slip wall. The inlet velocity condition was that of slug flow. Two different lengths (20h and 40h) were used to generate two different boundary layer thicknesses at the step. The outlet lengths and boundary conditions are the same as those in section 5.2. Grid refinement was done on one case (R = 100) from n  325 to 1200 cells. The present case consists of 999 cells. This is refined for a single case to 3120 for comparison. The outlet length is then increased from 35h to 47h for the same case with a total of 3216 cells.  Chapter 6  RESULTS  A N D ANALYSIS  Discussions in this chapter pertain only to the BFS as results for the Blasius case are summarized at the end of Chapter 3. The results will be presented as follows: (i) results of the new VISI method will be analyzed, (ii) results of the T E A C H code will be analyzed, (iii) results from (i) and (ii) will be compared. Some experimental results will also be shown.  6.1  VII Results  Results presented in this section include grid refinement and parametric variation of x and R^. Results of the grid refinement study will be analysed in terms of the e  displacement line, as it is (graphically) the most detailed. Results of the parametric study will be analysed in terms of the pressure coefficient (with some comments about reattachment) as it yields the most information. Additional results such as reattachment lengths, shape factors, and shear stress distributions will be discussed in section 6.3.  85  CHAPTER  6.1.1  6. RESULTS  AND  ANALYSIS  86  G r i d Refinement Xs=60  Rh=100  100  (x-x  0  Figure 6.1  )/h  VII Grid Refinement - Effect of x  N  Figure 6.1 illustrates the effects of grid spacing (cubic length) and the position of x  N  on the solution in terms of Ai(x). In each case the values of x„, x,,and Rh  are the same. The same grid expansion factor was used in each case, each cubic being 10% longer in the x-direction than the previous one. Reattachment occurs at [x — x )/h w 31. The results show that grid independence is realized for quite 0  coarse mesh sizes. The cases N = 9, 16, 20 had an initial cubic size (AXX) of 5h. The curves N= 22 and JV = 11 had initial cubic sizes of 2h and 8h respectively. The  CHAPTER  6. RESULTS  AND  ANALYSIS  87  cases with AXy — bh show a remarkable independence of Ax on x . This verifies w  the hypothesis that the boundary layer returns to a Blasius profile (ie: A becomes x  ^(x)) at some distance downstream of the step. The disturbance region is seen to extend ~ 36/i downstream of reattachment, which is much larger than the 14/i reported in the literature [7] for internal flow. No analysis was performed to check if the disturbance region varied with Rh or x . a  Figure 6.1 shows a peculiar 'dip' in the solution near the start of the interaction region. This was due to the fact that the beginning of the interaction region was inside the disturbance region. Figure 6.2 shows the effect of moving the beginning of the interaction region upstream. The dip is caused by the solution having to compensate for the fact that J (x) is not close enough to the real solution near the 0  step. The dip does not appear to affect the solution downstream of the step. For example, the change in L for the two curves shown is less than 1%. The effects R  of the dip on the pressure distribution are shown in figure 6.3. The hollow squares coincident with the zero C line are not calculated by the method, they come from p  the fact that upstream of the interaction region the conditions are set to those of the free stream. It it interesting to note that the pressure variation for this case (x = 100) immediately compensates to reach the same pressure minimum at the 0  step as the x = 50 case. Again as in figure 6.2 the effects are negligible downstream 0  of the step. This analysis shows that the disturbance zone extends more than 30h upstream of the step which is again much larger than reported in the literature [7] for internal flow. The upstream interaction is weak, however, and small changes to the upstream conditions do not significantly affect the flow downstream of the step.  CHAPTER  6. RESULTS  AND  ANALYSIS  6. RESULTS AND  CHAPTER  ANALYSIS  ©  w 6 +  a o o. <  On  " «•  X  © I  •  0  *S  100  050  -BO  Figure 6.3  -40  I  40  R  h  T  130 200  i  130 200 4  80  VII Grid Refinement - Pressure Distribution  120  CHAPTER  6. RESULTS  AND  ANALYSIS  90  Parametric Variation  6.1.2  The VII model for the BFS has two independent parameters, x and R . e  h  For an  internal BFS there is an extra parameter, the expansion ratio, which is a measure of the channel height relative to the step size. Some authors use the displacement thickness at the step with no step present, 6*, instead of x , but they are related s  through Rh and thus are interchangable. Some authors (eg. [7,12,14]) ignore or neglect the effect of x, on the flow downstream of the step. Conversely, Armaly et al[l3] conclude that the reattachment length, L , is not a function of Rh alone, and is indeed a function of the parameters r  listed above. Figure 6.4 shows clearly that the effect of x, is not negligible in the VII model. For the x, = 40 case A* « 0.78 and L a  L  R  R  « 12.6 while for x = 350, A* « 2.32 and s  « 18.2. The longer reattachment length is associated with smaller pressure rises  and a pressure peak which occurs further downstream. This observation can be explained with the aid of the inviscid velocity equation 4.2. As the boundary layer gets thicker (ie: larger x„) the ratio of h to Si becomes smaller. Thus any changes that the step induces in A become smaller as x, increases, and in consequence the x  derivative of A i is reduced. From equation 4.2 this causes Upor to shift towards 1, and since C« — 1 — ^POT pressure is lower. The lower pressure allows the slow moving fluid near the bottom of the detached shear layer to penetrate further downstream before it is forced to reattach. Figure 6.5 shows the effects of both x and T (the order of the velocity profile s  polynomial) on the pressure distribution. For each case (T = 4 or 5) the results are the same as those in figure 6.4. An interesting point to note, which is not easily  CHAPTER  C  6. RESULTS AND ANALYSIS  °-  X  D  I -20  X  S  R  h  T  300  350  200 4  020  040  200 4  I  20  40  (x-x  s  Figure 6.4  0  60  )/h  Pressure Coefficient - Effect of x,  60  CHAPTER  6. RESULTS  AND  ANALYSIS  6•  +  O  w5e  •  a •  o  d  A  4,  AO ft  I  -20  I  20  40  (x-x  s  Figure 6.5  T  Xo  *S  R  080  110  100  5  •  040  060  100  5  •  h  020  040  100  5  A 080  110  100  4  Q 040  060  100  4  O 020  040  100  4  60  )/h  Pressure Coefficient - Effect of T  60  CHAPTER  6. RESULTS  AND  ANALYSIS  93  ascertained from the figure, is that reattachment always occurs near the point of maximum pressure gradient. This observation is supported through computation of the parameter  near reattachment (not shown). This parameter is observed  to reach a minimum near reattachment in each case, so that from the second of equations 4.8,  is also at a minimum there. Bernoulli's equation shows that  this corresponds to a C gradient maximum. Figure 6.5 also shows the pressure p  gradient returns to zero in an asymptotic manner. This is consistent with the downstream boundary condition developed in Chapter 3. The order of velocity profile has a small but consistent effect on the solution. The pressure gradients for the case T = 4 are steeper than those for T = 5. This causes the reattachment lengths to be typically 5% longer in the latter case. This consistent difference warrents further investigation using Timman's [21] single parameter velocity profiles. These profiles were discovered by the author after the completion of this work. This was somewhat disappointing as Timman's profiles were claimed to (probably) be as accurate as a single parameter profile can be. Figure 6.6 shows the variation of C with Rh, with x held fixed. Upstream of p  4  the step the symbols show a definite independence of Reynolds number. This is quite remarkable considering that given the same values for x , the values of A i at 4  corresponding x-locations will be greater for lower Reynolds numbers. This may be caused by the weak interaction upstream of the step. In the Blasius case, where there is no interaction, the pressure remains fixed at zero, independent of R^. The parameter x', describing the virtual origin of the redeveloping boundary layer, was introduced so to allow for greater freedom in the solution, x' showed little variation with Reynolds number, but definitely increased with increasing x . 3  Values of x' ranged from 7 to 50 for x ranging from 40 to 130. Since x' is positive in s  CHAPTER  6. RESULTS  AND  ANALYSIS  (x-x  s  Figure 6.6  )/h  Pressure Coefficient - Effect of Rh  CHAPTER  6. RESULTS AND ANALYSIS  95  all cases, the boundary layer is effectively thinned by the presence of the step. This shows that the amount of fluid flow out of the boundary layer in the interaction region is greater than that of a zero pressure gradient boundary layer. The amount of fluid lost seems to be in some way proportional to R because x' remains nearly h  constant as Rh is varied. Further discussion of the VII results for the BFS is postponed to section 6.3. From the results shown so far the model seems to be yielding logical results. Some insight into the nature of the flow separation has been achieved because of the relative simplicity of the model.  6.2  T E A C H Results  Results presented in this section include a typical velocity field plot and parametric variation analogous to that in section 6.1.2. Some comparisons with VII results will be made, but the majority of comparisons are postponed to section 6.3. Grid refinement is discussed, but results are again postponed to section 6.3.  6.2.1  Grid R e f i n e m e n t  Grid refinement was performed using the conditions described in section 5.2.2. The resulting reattachment length varied from 4.2 (315 cells) to 6.1 (1200 cells). The variation between 1000 to 1200 cells was 3% so it was felt that 1000 cells would be adequate for this study. Hackman's [22] work on grid refinement shows that in order to have an accurate solution, a higher order differencing scheme must be used in conjuction with grid refinement. It must be realized therefore that the present T E A C H solutions may suffer from numerical error due to what Hackman refers to  CHAPTER  6.  RESULTS  AND ANALYSIS  96  as 'False Diffusion' (see also Djilali, ref. 20).  6.2.2  Velocity Profiles Xo=20  Figure 6.7  Rh=200  Y-EXPANSION 5:1  Velocity Vectors near the Step  Figure 6.7 shows a typical velocity vector distribution in a region near the step. The y-dimensions have been expanded 5 x for clarity. The center of the arrowhead coincides with the location where the vector is calculated (center of the scalar cell) and the length of the tail is proportional to the the magnitude of the velocity. A l l velocities are normalized with respect to the free stream.  The region of recircu-  lating fluid is shown clearly as is a point very near to reattachment. Assumptions concerning the use of the unmodified F L A R E approximation are verified as most of the velocities in the recirculation region are less than 10% of the free stream value.  CHAPTER  6. RESULTS  AND  ANALYSIS  97  The assumption of neglecting T (ie: |^) cannot be verified as the grid is too coarse W  to discern the gradient very close to the wall.  6.2.3  Parametric Variation  « =400 R« = 4 0 0 =100 =100 R  0 •  *  i~\ 0  i  i  i  i  i  10  20  30  40  50  (x-x  0  Figure 6.8  =100 =100  Xo= Xo= Xo= Xo= Xo= Xo=  110 20 100 40 20 5  60  )/h  T E A C H - Pressure Coefficient - Effects of x.  Figure 6.8 shows results analogous to figure 6.5. In all cases of T E A C H results the step is 20h downstream of x , so that varying x is equivalent to varying x,. In 0  0  figure 6.8 there are two cases of Rh, 100 and 400. The results are consistent with those of the V I I method in that increasing x, decreases the pressure gradient and  CHAPTER  6. RESULTS  increases L . R  AND  ANALYSIS  98  The discontinuity in pressure gradient inherent in the VII model is  also shown in the T E A C H results. Two descrepancies are apparent in the figure and require explanation. The first is the initial portion of the pressure curve for the case x = 5. This is due to the relatively coarse grid and thin boundary layer, 0  the combination of which causes the discretized form of the initial velocity profile to inadequately represent the true shape. The second discrepancy is the outlet pressure values. The pressures are obviously not returning to their freestream values as they should. The problem could be caused by incorrect outlet conditions, outlet conditions applied too far upstream, false diffusion, or a combination of these. In section 6.3 the grid is refined and the outlet conditions are moved downstream to evaluate some of these effects. Figure 6.9 is a plot of C vs R with x fixed, analogous to figure 6.6. This plot p  h  0  shows even more clearly the insensitivity of the pressure distribution upstream of the step on R^. In this case, however, the difference between the pressure maxima is significant. This could easily be influenced by the fact the the exit pressures have significantly different values. Results from the study discussed in section 5.2.2 are shown in figure 6.10. The conditions are similar to those of figure 6.9 except the top wall is a free slip impermeable boundary instead of a line of zero pressure. This shows clearly (note the difference in scale) that the two conditions have quite different effects even though they are imposed far enough from the step that the difference in their effects at the step could be thought to be negligible. In figure 6.10 the upstream independence of C on R is gone, and the adverse gradients are less than those of figure 6.9. This p  n  results in a reduction of 8% in L for the R = 100 case. It was found that for the R  h  case of the free slip top wall that the u-velocities near the top wall were accelerated  CHAPTER  6. RESULTS AND ANALYSIS  •  10  r 20  t *  i * ** T •  »  30  »  I 40  •  R„=400  Xo= 20  •  R„=200  Xo= 20  *  R,=100  Xo= 20  K= 50  Xo= 20  60  SO  (x-x )/h 0  Figure 6.9  T E A C H - Pressure Coefficient - Effects of R  h  CHAPTER  6. RESULTS  AND  ANALYSIS  100  O  o  °  c O  9  9  9  J  8  f  _ . .  K  e  *  4  * ' '  0>  o o  Si OH  4  KE*100  *  R E " SOC  O  RE" SOC  X  10  20  30  X  Figure 6.10  40  /  50  60  h  T E A C H - Pressure Coefficient - Case 2  70  CHAPTER  6. RESULTS  AND  ANALYSIS  101  10% from inlet to outlet. In the zero pressure case, acceleration was less than | %. The reduced pressure gradients in the free slip case are a reflection of the induced acceleration, this acceleration being caused by the imposed boundary conditions and not by the presence of the step.  6.3  Comparison of Results  The ultimate test for any model is a direct comparison with the physical process it is meant to simulate. If the physical process cannot be described or measured accurately then evaluation of the model's performance becomes more difficult. In addition to measurements from physical experiments, other models may be used for comparison. The T E A C H model used in Chapter 5 is a standard model for computing certain classes of flow, the BFS being one such class.  Although the  T E A C H code was developed for internal flow, it was felt that it could be modified to simulate external flow conditions. In the following section the T E A C H code is used to compare with results of the present VII model. Aside from the reattachment length plot (figure 6.11), all comparisons were made at R^ = 200, x, = 40, these being in the mid-range of the parametric variation. Experiments with high (> 10) channel height to step size ratios and large inviscid cores are also used for comparison. These conditions were the closest to simulating external flow. No open air experiments were found in the literature.  CHAPTER  6.3.1  6. RESULTS  AND  ANALYSIS  102  Reattachment Lengths  The most significant length scale in laminar flow over a BFS is the reattachment length and thus it usually receives the most attention. Figure 6.11 shows a plot of L  R  vs Rh for several values of x,. Results of some experimental observations are  Lr  500  Figure 6.11  Reattachment Lengths  also shown. The wide scatter of results shows that flow over a BFS is not solely a function of Reynolds number based on step size. The experimental results for Leal et al and O'Leary et al were taken from the paper by Denham and Patrick [14]. The equivalent values of x for these experiments are not given as their upstream t  CHAPTER  6. RESULTS  AND  ANALYSIS  103  conditions were different. For Goldstein et al [12] the range of x was from 4 to 30. e  Results from the VII model follow those of T E A C H although the reattachment lengths are considerably longer (15-80%). Agreement is better at higher Reynolds numbers. This makes some sense as the VII model has governing equations which are more accurate at higher Reynolds numbers. These differences could be caused by the deficiency of the VII model (because of the approximations used) and/or because of numerical errors in the T E A C H solution. Hackman [22] shows for a specific case (Rh « 300), L  R  is underpredicted by about 10% using hybrid differencing with a  relatively fine mesh. The effects of x are the same for both models although these effects are more a  pronounced in the VII model. It is difficult to give exact reasons for these differences, but they are probably the same reasons stated above. Although Goldstein et al [12] vary x„, they do not show explicitly the dependency of LR on x,, but only indicate that increasing x, increases L . Their data is R  presented as a least squares average with the range of results shown by the scatter bar.  The range of x used in the experiment lies below that used in the models. e  This tends to make the agreement look worse than if comparable values were used. The present initial guess for A (x) employed in the VII model is inadequate when x  high Rh and low x are required simultaneously, and in consequence no results are s  obtained for this case. The method of initial guess needs to be improved so that direct comparison with Goldstein's results can be made. The results of the other experiments shown infigure6.11 give fair agreement with the VII model. These results show clearly that a more precise definition of flow over a BFS is required if results are to be compared. Experiments should be run in an open air tunnel (if laminar flow can be maintained) over a wide range of x, values to provide  CHAPTER  6. RESULTS AND ANALYSIS  104  a fair comparison with the results of the VII model.  6.3.2  Other Characteristics  In this section comparisons between the models are given in terms of the pressure coefficient on the lower wall and some of the boundary layer characteristics. As previously mentioned, a mid-range case (R^ = 200, x, = 40) is used in all comparisons. Figure 6.12 shows the pressure distribution on the lower walls. The reattaches  n 6  =40  Xo =20  Rh =200  I  Figure 6.12  Wall Pressure Comparison  CHAPTER  6. RESULTS  AND  ANALYSIS  105  merit lengths for the VII and T E A C H models are 12.6 h and 9.5 h respectively. An examination of the T E A C H pressure field (not shown) indicates that the pressure in the viscous region varies only slightly with y, even in the recirculation region. This reinforces the VII assumption that as far as the pressure variation is concerned, the separated shear layer behaves like a boundary layer. The comparison between the pressure distributions is quite good considering some of the assumptions made in the VII model. Note that at the step both curves have approximately the same value and slope. The VII model then proceeds to increase pressure less rapidly than the T E A C H , thus causing the longer reattachment lengths. It may be that the VII model is unable to increase pressure fast enough due to the restrictions of the Pohlhausen profiles. The difference may also be due to the assumptions made in the recirculation region. Some of these assumptions are discussed later. The discrepancy in pressure tendencies far downstream was investigated through grid refinement in the T E A C H model. The grid density downstream of the step was doubled in both directions and the resulting pressure curve was almost identical to that of figure 6.12. A 3% increase in L also resulted. Next the outlet was moved R  downstream by 12 h to check the effects of proximity in applying the downstream boundary condition. The results, plotted in figure 6.13, show a better agreement downstream of the pressure peak. This enforces the earlier hypothesis that the interaction zone for the BFS is longer in external flow. The following two figures give more evidence that the present T E A C H model is not behaving as expected downstream of reattachment. Figure 6.14 shows the displacement line, A i , from 20 h upstream to 40 h downstream of the step, where the boundary layer should be redeveloped. Both methods compare favorably until just slightly upstream of reattachment. The VII curve continues smoothly downstream  CHAPTER  6. RESULTS AND ANALYSIS  Figure 6.13  Refined W a l l Pressure Comparison  106  CHAPTER  6. RESULTS  AND  ANALYSIS  CHAPTER  6. RESULTS  AND  ANALYSIS  108  with the boundary layer appearing to redevelop about 30 h downstream of the step or about 18 h downstream of reattachment for this case. The T E A C H curve turns abruptly just before reattachment and remains nearly constant until again about 18 h downstream of reattachment it begins to rise. Another point to note about these curves is that their derivatives are all quite small. Thus linearized inviscid theory is indeed applicable. A better indication of the state of the boundary layer is given by the shape factor, A*/©, shown in figure 6.15. In the recirculation region the shape factor is changed to A i / 0 so that no argument arises over the thickness of the shear layer for the T E A C H model. The horizontal line at A*/0 w 2.55 corresponds to a T = 4 Pohlhausen profile in zero pressure gradient. The agreement is quite good upstream of the step, which is expected as Pohlhausen profiles work well in favorable pressure gradients for attached flow. Agreement in the recirculation region is also fairly good which suggests that the neglect of velocities in this region is not as drastic as it may seem. The agreement becomes progressively worse downstream of reattachment where the T E A C H shape factor continues to decrease. The VII model produces a curve that would be expected for a boundary layer which redevelops after encountering a small disturbance. The outlet shape factor produced by the T E A C H model (~ 1.6) could never be realized by the Pohlhausen profile , as it has a minimum at ~ 2.25. The low value corresponds to a proportionately high favourable pressure gradient. Perhaps introduction of the pressure condition at the top wall in conjunction with the zero velocity gradient condition at the exit in the T E A C H code is the cause. In the study discussed in section 5.2.2, the condition of an impermeable no-slip top wall caused the inviscid fluid near that wall to accelerate. Perhaps the zero pressure condition at the top wall, in compensation for not accelerating the fluid near that  CHAPTER  6. RESULTS  AND  ANALYSIS  A' =20  A" =40  0  s  •  109  R -Z00 h  •  vii  •  TEACH  • •  OS  o  •  E—  o  w < as m •  V- • • •  20  T"  40  60  X  Figure 6.15  •  60  /1\  Shape Factor Comparison  100  CHAPTER  6. RESULTS AND ANALYSIS  110  wall, has accelerated the fluid far beneath it causing the boundary layer to distort. It is apparent that a detailed study of the effects of different boundary conditions is required for the T E A C H code. Most of the simplifying assumptions used in construction of the VII model are for the recirculation region. The modified F L A R E approximation (see section 4.2.2) was used to permit forward marching of the viscous solution in this region. A consequence of this assumption is that the zero velocity line becomes coincident with the separation steamline. Figure 6.16 shows these lines computed from the T E A C H and VII models. The longer reattachment length from the VII method is reflected in the graph. Notice that neither curve from the T E A C H extrapolates to the top of the step. This is most likely caused by the coarseness of the grid. If the x-axis is rescaled as (x — x )/L a  R  the VII curve lies about midway between the  T E A C H curves. Perhaps the most drastic of the assumptions used in the VII is neglecting the shear stress at the walls in the recirculation region. Figure 6.17 shows a plot of the shear stress coefficient on the lower wall. If the T E A C H model has fair accuracy, the assumption of zero shear stress appears to be somewhat inaccurate. Equation 4.5 gives some insight into the effects of the assumption. Examination of T, (the shear stress at the bottom of the free shear layer) shows that its magnitude is about equal to r at the point where T is a minimum. From equation 4.5, with T m —r, w  W  w  the right of the equation increases by a factor of two. Of the variables on the left side of equation 4.5, the only one which would change by a large amount would be ^f.  Increasing the magnitude of ^f- would bring the VII curve in figure 6.12  closer to the T E A C H curve in the recirculation region. The above effect is only present near the center of the recirculation region as the shear stress goes to zero at  CHAPTER 6. RESULTS AND ANALYSIS  X =40 S  X =20 O  111  R =200 H  O  TEACH SEPARATION STREAMLINE  A  TEACH ZERO VELOCITY UNE  oo o  Y/h A  O A A  T"  10  (x-x  s  Figure 6.16  15  20  )/h  Separation Streamline and Zero Velocity Line Comparison  CHAPTER  6. RESULTS AND ANALYSIS  (x-x  s  Figure 6.17  112  )/h  Shear Stress Coefficient Comparison  CHAPTER  6. RESULTS  AND  ANALYSIS  113  either end. Some semi-empirical model for T could be included in the VII model W  to improve the results.  6.3.3  Execution Times  Computation of the parabolic boundary layer equations in the VII model involves the use of library integration routine. This routine works well when applied to the Blasius problem, but causes much difficulty with the BFS application. The routine is essentially a 'black box' and does not allow integration to cease before a prespecified location. The method of equation switching (section 4.2.2) allows the singularity at separation to be bypassed. The location of separation is not known a priori and thus the routine will not stop before reaching the singularity. A method of using very small integration sections was devised to overcome this, but it added considerably to the CPU time. A single iteration was performed using the above method and the resulting reattachment length was used to recompute the same case in a normal manner. The result was a reduction of 40% in CPU time. Even with the present method employed, the integration routine often reached the singularity because the step size was too large. It is difficult to equate the degree of convergence for the T E A C H and VII solutions. The convergence criteria were freely selected as [&h/h]  < 10~ for the s  MAX  VII model and, for the T E A C H , a residual error of < 10~ for the maximum of the s  mass, u, and v momentums when normalized with respect to the inlet momentum. For comparison, reducing each iteration criterion by a factor of 10 caused a 0.04% change in L  R  for the VII method and a 0.37% change for the T E A C H . For cases  which did not encounter the separation singularity, the VII method converged in  CHAPTER  6. RESULTS  AND  ANALYSIS  114  0.4 to 1.9 seconds for grid sizes ranging from N = 9—22. In comparison the T E A C H method converged in 22-41 seconds for the 37 x 27 mesh, and 150-180 seconds for the 65 x 48 mesh. All calculations were made on the UBC Computing Center's Amdahl 5850 mainframe computer. It should be reiterated that the T E A C H code was not designed for external viscous-inviscid flows and thus is not optimized to solve these types of problems. It is unlikely, however, that the convergence speeds could be increased by an order of magnitude by such an optimization. In light of the vast improvements in CPU time alone, further development of the VII model merits some attention. An integration routine should be found or written which will allow integration of the boundary layer equations up to a dynamically specified point.  Chapter 7  CONCLUSIONS A N D RECOMMENDATIONS  7.1  Conclusions  A new viscous-inviscid interaction (VII) method has been developed for predicting the flow field arising from a combination of inviscid potential flow and viscous flow. The technique involves matching the bounding velocities for the inviscid and viscous regions by iteratively solving for the displacement thickness, 6*(x). The formula used to update S*(x) after each iteration is generated by linearly perturbing the governing equations. When the bounding velocities match at a selected number of control points to within a specified accuracy, convergence is achieved, and the flow field prediction is complete. The new method is of a type known as 'semi-inverse', and in comparison to other methods is quite general and quite robust. In the technique as developed here, linearized potential theory is used for the inviscid flow and a boundary layer momentum integral model is used for the viscous  115  CHAPTER  7. CONCLUSIONS  AND  RECOMMENDATIONS  116  flow region. Neither of these assumptions is essential to the method, but both are useful here provided: i) the potential flow region is only slightly disturbed by the presence of the viscous region ii) the regions of separation, if they exist, are thin and of finite length (ie: strong interaction between the regions is of finite length) Further and somewhat drastic assumptions have been made in the development of the momentum integral model to allow the prediction of separated regions in a simple manner. Again the assumptions are not necessary for the general technique or the integral model and could be improved if greater complexity was clearly necessary. In principle however, the present assumptions provide a simple way of rapidly obtaining approximate predictions of complex flow regions. In addition to the above requirements, a method is derived for ensuring that the flow returns to its undisturbed state far from the region of strong interaction. Again this method is presently developed forflowswith the conditions (i) and (ii) given above, but these are not required in general. Application of the VII procedure to predict the unseparated flow past a flat plate gives excellent results, producing numerical solutions essentially indistinguishable from the appropriate analytical solution in less than 0.5 seconds of CPU time on an Amdahl 5850 computer. Application of the technique to external flow over a backward facing step leads to encouraging trends. The following results have been observed in the solution: 1. The region of interaction between the viscous and inviscidflowsis found to be very large (> 30h upstream of the step and ~ 40h downstream), much larger  CHAPTER  7. CONCLUSIONS  AND  RECOMMENDATIONS  117  than has been reported for internal flow [7]. 2. Downstream of the step, the technique successfully predicts a redeveloping boundary layer in the interaction region downstream of reattachment, as expected. 3. In all cases studied, values calculated for the virtual origin, x', of the redeveloping boundary layer show that the viscous region is effectively thinned by the presence of the step. 4. The calculated reattachment lengths, LR, are clearly influenced by the thickness of the boundary layer upstream of the step, thicker boundary layers producing longer reattachment lengths. This is also reflected in the pressure distribution as thicker boundary layers have shallower pressure gradients. 5. The reattachment length (for a fixed b.l. origin) becomes nearly proportional to the Reynolds number based on step height (Rh) for Rh > 150. Below this value, LR decreases more rapidly with decreasing Rh. 6. For a given fixed b.l. origin, the pressure distribution upstream of the step and the pressure peaks are nearly independent of Rh. 7. The interaction upstream of the step has a weak influence on the region of fluid downstream. 8. Minor alterations of the velocity profiles used in the momentum integral model have a small but consistent effect on the results. For the observed results, good accuracy is achieved for a relatively coarse distribution of control points  (AX  m t n  w  6/i),  and rapid convergence (less than 2 seconds  on the Amdahl 5850 computer) usually occurs.  CHAPTER  7. CONCLUSIONS  AND  RECOMMENDATIONS  118  Largely as a basis for comparison with the VII method, finite-difference predictions using an elliptic computer code (TEACH-T) have also been made for the backward facing step. The code has been modified to simulate external flow by imposing a zero pressure gradient on the upper boundary. This condition is more accurate than a free slip or no slip surface for external flow. A velocity profile representing laminar Blasius flow is imposed at the inlet and a zero velocity gradient condition is used at the outlet. The latter is somewhat inaccurate, but results can be improved by imposing this exit condition further downstream. Grid refinement from 999 to 3120 cells produced a slight (3%) change in LR. The numerical results are likely to suffer from false diffusion [20,22] and underestimates of 10% in LR could exist because of this deficiency. The VII andfinite-differencemethods, when applied to the same BFS problem, give similar trends in LR as a function of Rh and x„ (distance from origin of b.l. to step - a measure of the b.l. thickness). The values of LR obtained with the VII method, however, are 15-80% longer than those from T E A C H , closer agreement occurring at higher Rh. Direct comparison with experimental data is confused by the fact that the latter data does not clearly identify the effects of x, in the resulting values of LR. Trends appear to be the same for all computed and observed cases however. Disagreement between the VII, T E A C H , and experimental results are thought to be due to a combination of the following effects: 1. Neglecting the velocities in the region of flow reversal in the VII model. This was done to permit forward marching of the solution for the viscous region and is a modified form of the F L A R E [16] approximation.  CHAPTER  7.  CONCLUSIONS  AND  RECOMMENDATIONS  119  2. Neglecting the wall shear stress (r ) in the recirculation region in the VII w  model. Again this was done to construct as simple a model as possible, however the results of the T E A C H code show that T is not negligible in the W  central region of the recirculation bubble. Simple analysis of the VII model shows this assumption leads to longer reattachment lengths. 3. Inadequate velocity profile representation in the VII model. It has been shown [17,21] that the Pohlhausen profiles used in construction of the momentum integral model often give inaccurate results in regions of adverse pressure gradient. 4. False diffusion effects in the T E A C H results. This numerical error arises when low order differencing schemes are used in regions where the flow in inclined to the grid and strong transverse gradients exist. 5. Inaccurate exit boundary conditions in the T E A C H code. Use of a zero velocity gradient at the exit is inconsistent with the conditions of a redeveloping boundary layer in external flow. The condition of zero velocity gradient is equivalent to accelerated flow conditions in which the boundary layer does not grow. 6. Incomplete experimental control and/or specification of the upstream b.l. thickness effects on measured LR. In some cases the upstream conditions were different than those used by the VII and T E A C H models, making comparisons more difficult.  CHAPTER  7.2  7. CONCLUSIONS  AND  RECOMMENDATIONS  120  Recommendations  1. A theoretical or semi-empirical model of the wall shear stress in the recirculation region should be incorporated into the VII model. 2. Timman's [21] velocity profiles should be tried in the VII momentum integral model in place of the present Pohlhausen profiles. 3. A family of velocity profiles should be developed that are capable of representing the reverse flow region in the separation bubble. This would automatically give a model for the wall shear stress. These profiles will most likely be a function of more than one parameter, so that care must be taken to ensure that velocity profile critical points reported in the literature [5] do not occur in the computational domain. 4. Improve the method of making the initial guess for Ay(x). The present method employs a single cubic as an initial guess, but this curve becomes too distorted (ie: contacts the wall) for certain parameter ranges. A slightly more constrained curve should be implemented. 5. Apply the VISI method to external turbulent flow over a backward facing step. Following the idea behind the present model, use the simplest turbulence models first, and more complex ones if needed. 6. Apply the VISI method to a different class of flows, such as uniform flow past a blunt rectangular section. This would require modification of the inviscid equations and possibly equations of higher order than the boundary layer equations for the viscous region. 7. Modify the T E A C H - T model to be parabolic in the viscous region (outside  CHAPTER  7. CONCLUSIONS  AND  RECOMMENDATIONS  121  the recirculation region), so as to make the code more efficient for solving viscous-inviscid flows. Implement Abdullah's [19] pressure condition at the outlet as well as at the top wall. Use a higher order differencing scheme to reduce numerical error. 8. Conduct experiments which closely simulate the conditions of external laminar flow over a backward facing step. Extend for example the range of x„, to get a clear indication of its influence on the reattachment length.  References [l]  Briley, R. and M Donald, H. A Survey of Recent Work on Interacted Boundary Layer Theory for Flow with Separation. Num. Phys. Aspects of Aero. Flows II pp.141-162 Springer-Verlag  [2]  Williams, J.C. Incompressible Boundary Layer Separation. Ann. Rev. Fluid Mech. Vol. 9, pp.141-162 (1977)  [3]  Goldstein, S. Q.J. Mech. Appl. Math Vol. 43 (1948)  [4]  Crimi, P. and Reeves, B.L. Analysis of Leading Edge Separation Bubbles on Airfoils. AIAA Journ. Vol. 13, No.ll (1976)  [5]  Shamroth, S.J. On Integral Methods for Predicting Shear Layer Behavior. J. Appl. Mech. Vol. 36, pp.643-681 (1969)  [6]  Green, J.E. Two-Dimensional Turbulent Reattachment as a Boundary Layer Problem. Agard ch 4. (1966)  [7]  Kwon, O.K. and Pletcher, R.H. A Viscous-inviscid Interaction Procedure Part 1: Method for Computing Two-Dimensional Incompressible Separated Channel Flows. J. of Fluids Eng. Vol.108, pp.64-70 (1986)  [8]  Williams, B.R. The Prediction of Separated Flow Using a Viscous-inviscid Interaction Method. Aero. Journ. pp.185-197 (May 1986)  [9]  Klineberg, J.M. and Steger, J.L. Calculation of Separated Flows of Subsonic and Transonic Speeds. Third Int. Conf. Numer. Methods in Fluid Mech., Proceed. Vol.11, pp.161-168 (1973), Springer-Verlag  c  [10] Acrivos, A. and Schrader, M.L. Steady Flow In a Sudden Expansion at High Reynolds Numbers. Phys. Fluids Vol.25, No.6 (June 1982)  122  REFERENCES  123  [11] Moore,T.F. Some Experiments on the Reattachment of a Laminar Boundary Layer Separating From a Rearward Facing Step on a Flat Plate Aerofoil. J. Royal Aero. Society Vol.64, pp.668-672 (Nov. 1960) [12] Goldstein, R.J., Eriksen, V.L., Olson, R.M., and Eckert, E.R.G. Laminar Separation, Reattachment and Transition of the Flow over a Downstream-Facing Step. J. Basic Eng. pp.732-741 (Dec. 1970) [13] Armaly, B.F., Durst, F., Pereira, J.C.F., and Schonung, B. Experimental and Theoretical Investigation of Backward-Facing Step Flow. J. Fluid Mech. Vol. 127, pp.473-496 (1983) [14] Denham, M.K. and Patrick, M.A. Laminar Flow over a Downstream-Facing Step in a Two-Dimensional Flow Channel. Trans. Inst. Chem. Engrs. Vol. 52, pp.361-367 (1974) [15] Schlichting, H. Boundary Layer Theory. Seventh Edition, pg. 131 , McGraw Hill. [16] Reyhner, T . A . and Flugge-Lotz, I. The Interatction of a Shock Wave with a Laminar Boundary Layer. Int. J. of Nonlinear Mech. Vol. 3, No. 2, pp.173-199 (1968) [17] Curie, N. and Davies, H.J. Modern Fluid Dynamics, Volume 1: Incompressible Flow. pp. 209-210 , D. Van Nostrand Co. Ltd. [18] Patankar, S.V. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Co., New York (1980) [19] Abdullah, Z. A personal communication. (Jan. 1988) [20] Djilali, N. An Investigation of Two-Dimensional Flow Separation with Reattachment. Ph.D. Thesis, pp. 86-87, U.B.C. (July 1987) [21] Rosenhead, L. (Editor) Fluid Motion Memoirs - Laminar Boundary Layers, pp. 292-332, Oxford University Press [22] Hackman, L.P. A Numerical Study of the Turbulent Flow over a Backward Facing Step using a Two-Equation Turbulence Model. Ph.D. Thesis, University of Waterloo (1982)  Appendix A  A.l  Asymptotic Solution to the Boundary Layer Equations Integrated along A* = a  x-x' UR c  h  The solution to the b.l. equations as x —* oo will be examined. The variable U in c  the equation for A* will be shown to equal the constant free stream value of U  BL  at  oo. From the second of equations 2.36: dU dx.  A A R  BL  (A.l)  2  h  From the integration path equation for A* it is seen that A* oc \fx as x —• co . In order for the velocity profiles of the Pohlhausen method to make physical sense, the value of A must be bounded. From the first of equations 2.36, A* —  = ci + c  2  A  Because A is bounded, it follows that A oc y/x as x —> co.  (A.2) Substituting into  equation A . l ,  dV~BL  —— dx  A  a x  124  .  (A.3  APPENDIX  125  A.  Again, because A is bounded, it can be written as a constant plus terms which become negligible near oo. The value of A U  B  L  0  corresponds to the asymptotic value of  .  A=  A  +  0  0 (l)  as  x —• oo  Substituting into equation A.3, dU A . 0(1) —.— oc 1 — ax x x BL  Integrating, U Now, because U L is bounded, A B  OC yl lnx + 0(0)  BL  0  0  = 0 at x = oo and therefore,  A —• 0  as  x —* oo  Differentiating equation A.2, A A - A' A* = c A' A" 2  2  (A.4)  Now since A*, A are both oc -y/x then A*' and A' are both oc x~a as x — > do (because A' is also bounded). In the limit as x — > oo equation A.4 becomes: A' oc x , -1  or A' —• 0  as  x —* co  Substituting A = A' = 0 into the boundary layer equations (see section 3.1), lim  X—oo  A  A  *  x — X =  a  UL  Rh  B  This equation is identical to the asymptotic integration path except for the variable U. Because the above relation must be true, U L must assume the constant value B  APPENDIX U  C  A.  as the integration path relation must also hold, ie:  UL B  So by choosing U  C  —>  U  E  as  x — • co  = 1 it is necessarily true that asymptotically  constant and A = A' = 0.  Q.E.D.  UBL{*)  Appendix B  B.l  Calculation of Asymptotic Value for c7p0r(x)  In this section it will be shown that Upor( ) ~~+ 1 as x—•oo. The inviscid velocity x  distribution is given as,  Ppor(x) = 1 + - / J-oo  7T  -*  — d£  X- f  The displacement thickness function, A*(x), can be divided up into its constituent parts (see figure 3.1. In doing so the integral for UPOT[X) is subdivided into integrals corresponding to the constituent curves 7 (x), J„(x), and 7E(X). 0  Upoi{x]  =  u*-mf  1+  TT  X  JO  £  di +  ^hr-mm^ ^  X-£  TT JXn-!  In  i r °° 7T 7x I*N +  w  d  X x— - £  ^  '  IE  It is easily seen from the above equation that if the sum of all integrals vanish as x —»• oo then  U OT P  will  D  e  equal to 1. The following is a proof that each separate  integral vanishes as x —> oo and consequently 'Upor ~~*• !• Substituting for the 127  APPENDIX B.  128  function 7 (see section 3.1.1) into the integral I , 0  0  j  a  f°  _  di  x  a  In Jo  2n\fR~ y/x  \/x +  In  h  From the above equation in the limit, lim T X-»oo i  Next, substitute for J  in  n  a  =  0  Inl  = 0  I. n  1 /*- (j3 + 2 C £ + 3 0 £ ) di f J  I  B  n  n  TT Jxn -i  The limit of this equation will be taken while it is still in integral form, ie: f  m  xl oo/n  =  {B + 2C i +  Xn  — 7T X  n  3D i )di 2  n  n  JX„-i  The integral is simply the value of A*(x) at the right endpoint of the n  t h  cubic  minus the value at the left endpoint and so is finite. Taking the limit, lim X-oo  l  1  T  n  TT  [ A*(x ) - A*(x _ ) ] = 0 n  n  x  CO  Finally it will be proven that in the limit the integral I  E  is zero. Care must be  taken to evaluate the Cauchy principal part of the integral as there is a singularity at x = i. I I  E  =  E  is written as the sum of two integrals thus, a  f  x  -  di  £  lTT\/R~h J*N  h J** (x -  i) Vt-*' + n  _7r\/x —;  ln  y/i^-y/x^x Li  1  a  di  r  c  1T\/R~h Jx-  2*VRh  Jx+t  X-£  +  In TT\f) x —X'  (x -  i) vi-x  1  y/j^+Vx-X' y/i^-Vx-X'  X+e  APPENDIX B.  129  The lower limit of L\ is zero because Xoo » x and so the In argument tends to 1. N  The upper limit of L is zero because as f —» oo the argument of the In function 2  tends to 1 (here too ^ *oo ^ * ) and the term preceding the In function tends to N  ^ . This leaves the limits surrounding the singularity, or IE  1  =  7Ty/x — X  (y/x-x'-e + y/x-x') (y/x-x' + e -  In  y/x-x')  (y/x-x'-e - y/x-x') (y/x-x' + e + y/x-x')  1  Multiplying out the terms and letting 7 = x — x', 1 7Ty/x — X  In 1  y/7 - 6 - y/7 + ^ ( y / T + f ~ y / T ^ ) 2  2  2  v V - -  V7 2  y/7  ( y/7+e  -  y/l~ ) €  Using the binomial theorem and the fact 7 ^> e, 1 e  2  ~e  2  \/7  2  2  «  V/T " ^ -  2  7  and  Substituting the above in equation B . l , IE  1  =  1 -  In  1  £  2 t  •1 _ l i  7Ty/*  2 7  X  so that finally, lim l»m , Tf—»oo e—»0 IE  —  U  Putting it all together, lim T T 1 \ X — 0 0 Up or (Xj  , +  lim lim X—oo i  =  , 1  =  1 + 0+0+0  =  1  T  0  . +  TV lim lim r i X—oo 2^ i n + n=l  Q.E.D.  lim r X—00 i l  i  m  (B.l)  APPENDIX  B.2  B.  130  Calculation of  Kyfi  for A*(x) =  U T(X) PO  In this section it will be shown that UPOT(*) = 1 for 0 < x < oo with K = const, The equation for the potential velocity distribution gives,  7T  di  X —£  J-oo  di  K_ r+°  1 +  2n Jo  Taking the Cauchy principal part of the integral,  £W(x)  =  X-e  K  1 +  27T /X  In  +  V  K  .  27r  \/x  ln  V*~+y/i X+c  The limits 0 and oo both give zero function values as they both cause the argument of the ln function to tend to 1. Substituting the two remaining limits and combining,  UPOT{*) = *™o  (y/x + \fx-e) 1  + T^7=  l n  {y/x -  y/x+e)  [y/x - y/x^e) [y/x + y/x + l)  Using the binomial theorem and neglecting higher orders of e,  UPOT(X)  — «™o  1 +  = In 2xy/x  ( V^ " ife) 2  Cancelling terms and letting e —• 0,  ^POT(X) =  1 +  K 2-Kyfx  1+0  Q.E.D  ln  2y/x~ 2yjl  Appendix C  Cl  Calculations  UpoT{*k)  UpoT[*k) =  dt  1 + - /  7T JO  r°°  a  + 7T '*« Jx  0 s/t +  (x* -  1 "  1  - t — f 7T ~ AX„  rb 1  n  + 2c t + n  3d t  2  n  dt  rpk-t  Jo  dt  (x* - t)  y/i-x'  Integrating, and taking care to evaluate the principal part of the integral under the summation sign and also noting the expression for Jf;(x ) is valid for x'<Xjfc  }  fc  Uporixk) =  1 +  In  1" \Jx~k \JXk~—\/Xo  7T  y/x - X' k  Jo{x )  TT ~  N  k  y/x  X'  - x'-y/Xk - X'  N  £  (6 + 2c e+3d„e ) In 2  n  y/x - X' + y/x J (Xjt)  k  + -T —  In  n  3d(0jk + -) - 2 c „ n  AX„  (Cl) This equation is valid for any value of ip except 0 or 1. The matching points were k  selected to coincide with the cubic endpoints, thus giving tpk = 0 when k — n—1. In If x' > x  4  then JEM  2o  = * yV  —  Xi  - arctan \  131  APPENDIX  C.  132  other words a singularity is created in equation C l at x (between the two cubics f k  k  and fk+i). This singularity can be removed if the principal sum of the two integrals surrounding the singularity is evaluated. There are two distinct cases to consider: 1) The singularity is between two cubics 2) The singularity is between the last cubic and ^JS(X).  For case l) let Ik,k+i be the principal sum. 1  T  lim  J-k,k+i =  {b +2c {l-E ) + 3d {l-E ) }\n  -3d {%-e )-2c  2  «i-»0  k  1  k  k  1  k  1  k  AX 7T t  1  lim e -»0 2  A X  i  +  1  {b  k+1  + 2c £ k+1  £2  + 3d £l}\n  2  k+1  l-e  7 T  3d {±  + e ) - 2c  k+1  2  k+l  2  From continuity of slope and the scale factor between the cubics /* and fk+i, b i  _ b + 2c + 3d  k+  k  k  AXfc i  AXfc  +  Substituting into  k  £ i _  A x  t  E  AX  f c  2  +  1  and taking limits,  (C.2)  For case 2) where the singularity (matching point) is at x , N  1 AX 7T  lim e,-0  {b + 2c {l-e ) + N  N  1  3d (l-e ) }ln 2  N  1-ei  1  - 3d {* - ) - 2c N  N  , 1  H —  7T  lim « ->0 2  y/x ~X N  f  + E  In  2  \ / x - X ' + \/x -X'  +E  y/x -X'  +E  N  N  - y/x -x'  N  N  2  2  From continuity of slope and the scaling factor between f b + 2c + 3d N  N  AX,  N  and  E  X  — E  2  N  and 7E, 1  A X  W  £l  N  APPENDIX C.  133  From the binomial theorem,  y/x -X' + S N  Vt-N-X'  ~  2  +  £2 2y/x -X< N  Substituting into I , N  N+1  and taking limits,  7T y/x —  X'  N  In  4 (x - x') N  AX,  (2c„ + §<f„) AX.  (C.3)  The various principal sums of equation C . l have been calculated. Chapter 3 contains the final form of the equation for Upor{xk)-  C.2  A[/por(xfc)p Calculations  Calculations for At/por(xit)p are divided into two main cases: l) p 7^ N 2) p = N. Each of these cases is further subdivided depending of the relative location of the velocity matching point, X * , to the perturbation position, x . p  p^N For this case the perturbation curve is entirely within the interaction region because of the local nature of the perturbation curve. Thus the equation for AUpor{xk)p can be written as,  AUpoT{Xk)p  —  TT ~  AX„  J  (CA)  ^k-i  KUporixkii where, from Chapter 2,  Kf {rP) = Ka + Kb i> + Kc i>l + Kd ^ n  n  n  n  n  z k  APPENDIX  C.  134  There are four subcases to consider, labelled © . . . © . Subcase © : x ^ x^  , x, x  k  p  p + 1  In this subcase equation C.4 contains no singularities and so can be integrated in a straightforward manner to give,  KUp {*k)p = H  AJp(^Jt)  OT  AJ  +  p+1  (0 ) ] fc  where,  {Kb + 2Kc ib n  n  + 3Kd ib )  ln  2  k  n  k  AX„  4k  1 -fb  - 3Kd {rp + i ) - 2Kc n  k  n  k  Subcase © : x = x _, k  p  In this subcase the singularity at x* must be avoided. To take the principal part of equation C.4 the lower limit of the first integral must be replaced with a small value, e, which is then shrunk to zero , ie: AZ7p (x ) 0r  fc  p  - M  r>  i  = «™o  7T  Jt  AX„  j mo  di  Ah„  4k-i  KUpQT ( X j t ) j :  + Substituting for Kf (i) p  7T A X  and Kf (i), p+i  1 KUpoT{*-k)p = 7T  Subcase ® : x = x k  Dp  ++  ,  JJ O o  1  Ahr,  di  i  1pk~  integrating and taking limits, \Kd AJ {4 ) p+1  k  -  2  p  -  + 2Kc  p  V  P  A X  -  C  p  For this subcase limits must be taken on both integrals of equation C.4, as the  APPENDIX  C.  135  singularity is between f and / +i, ie: p  p  lim  e—0  A^por(xjfc),  7T  AX  d£  Jo  D  KUpoT{xk)i  lim  +  e->0  7T A X  P+1  de  ./«  A/I„  Integrating and taking limits using the fact that there is continuity of slope between  f and A / p  p + 1  at  x = Xjt,  Kb  Subcase ® : x* = x  ^  p+1  KUpoT{*-k)p =  7T  +  \Kd + 2Kc ^ \Kd p  p  AX  AX P + I  D  +  p+1  + 2ifc  p+1  AX , P+  p  This subcase is analogous to Subcase © except that the upper limit of the second integral in equation C.4 is replaced by 1 — e. The resulting expression is, \Kd  + 2ift  v+1  KUpoT{*-k)p =  •K  A J  P  ( ^ * )  p + 1  AX  p=N In this case the perturbation curve extends from x _ to oo. The equation for N  AUp (xk)p OT  t  is, d  AUpoT{Xk)  p  —  7T AX N JO  zrt (t\ Ipk- <  KUpoT{*k)p  For this case there are three subcases to consider. Subcase © : x < x _ t  N  t  .„ d  Ah  N  J*N  (C.5)  APPENDIX C.  136  For this subcase x  is completely outside the perturbation region so that equa-  k  tion C.5 contains no singularities and can be integrated in a straightforward manner to give, KUpo {Xk) T  = -[AJ (ip )  N  N  k  +  AJ (x )\ E  k  where A J E ^ ) is the second term in equation C.5. Using standard integral tables,  x* - x' I AJjs(Xjfc)  =  2  ln  A+ l A-l  x* > x'  -)  -  I  1-A' X ' - X i  7T  — — arctan A' 2 ]}  x* < x'  where, A=  and  V xjb-x'  IX  A'  M  x  V x'-x*  Subcase © : x^ = x _, N  For this subcase the principal part of equation C.5 is calculated by replacing the lower limit of the first integral with e:, integrating, and then taking the limit of the resulting expression as e —• 0, ie: KUporix-k)*! = Subcase ® : x^ = x  Kd + 2Kc N  N  AX,  7T  N  For this subcase the principal part of equation C.5 is calculated by replacing the upper limit of the first integral with l — £\ and replacing the lower limit of the second integral with x + e - The limit of the resulting expression is calculated as ei and N  2  e tend to zero. The calculations below are done in detail to provide an example of 2  evaluating the Cauchy principal part of an expression. From equation C.5, 1-ei  KUpoT(X-k)p —  «i.«a-»o  —  7T  AX„  JO  Kf U)  K7 (0  N  B  JXN+e2  X  f c  —  APPENDIX C.  137  Using integral tables, inserting the limits, and recalling that Kb = 0, N  KUpoT{x ) k  1  lim «i-o  +  -  N  =  (2ifc {l-e } + 3ia {l-e } ) 1  N  1  AX,  7T  V +x ^ v + ^ /l  lim  «2-»0  3/a {f-e }-2iirc  2  N  2TT{X -X'}  In  +  /  N  In  w  1  K  AX,  1  r  1  Using the binomial theorem,  '1 +  ^2  X  W  1 +  X  £2  2 {x„ - x'}  therefore,  KUPOT{X)C)N  = 7T  AX  {2Kc + 3Kd ) ln N  - \Kd - 2Kc  N  N  N  N  + 2TT{X„-X'}  ln  4 {x„ - x'}  -2  ^2  From continuity of slope between A/ and A J E at x = x , and from the scaling N  N  equation for xjj, 2Kc + 3Kd N  AX,  e 1 — = e AX,  -1  N  2 K-x'}  x  and  2  Substituing into -K^or KUpoT[xk) where,  N  =  -[AJ (*N)] E  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items