A V I S C O U S - I N V I S C I D I N T E R A C T I O N P R O C E D U R E By Dave Stropky B.A.Sc, University of British Columbia, 1983 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D 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 to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A April 1988 © Dave Stropky, 1988 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 b s t r a c t 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 indistinguish-able 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 o n t e n t s A b s t r a c t ii N o m e n c l a t u r e vi i i A c k n o w l e d g e m e n t s xi 1 I N T R O D U C T I O N 1 1.1 Problem Description 4 1.2 Literature Review 6 1.2.1 Theory 6 1.2.2 Experiments 8 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 1 0 2.1 Current Interaction Schemes 10 2.2 A New VISI Technique 12 2.2.1 A Simple Example of the VISI Technique 13 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 3 A P P L I C A T I O N T O T H E B L A S I U S P R O B L E M 37 3.1 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 Results and Analysis 56 3.3.1 Radius of Convergence and Stability 56 3.3.2 Solution Accuracy, Grid Refinement, and Execution Speed . . 59 4 A P P L I C A T I O N T O A B A C K W A R D F A C I N G S T E P 62 4.1 Problem Geometry 63 4.1.1 The Displacement Line 64 4.2 Modifications to the Governing Equations 65 4.2.1 The Inviscid Equation 65 4.2.2 The Viscous Equations 66 4.3 Parameter Calculations 74 4.3.1 Velocity Calculations 74 4.3.2 Velocity Perturbation Calculations 75 5 T E A C H 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 Problem Geometry 81 5.2.1 Special Cells 83 5.2.2 Previous Work and Present Refinements 83 6 R E S U L T S A N D A N A L Y S I S 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 CONTENTS v 6.2.1 Grid Refinement 95 6.2.2 Velocity Profiles 96 6.2.3 Parametric Variation 97 6.3 Comparison of Results 101 6.3.1 Reattachment Lengths 102 6.3.2 Other Characteristics 104 6.3.3 Execution Times 113 7 C O N C L U S I O N S A N D R E C O M M E N D A T I O N S 115 7.1 Conclusions 115 7.2 Recommendations 120 References 122 A P P E N D I C E S 124 A 124 A . l Asymptotic Solution to the Boundary Layer Equations Integrated along A* = a x — x VCRK 2 124 B 127 B. l Calculation of Asymptotic Value for UPOT{X) 127 B. 2 Calculation of UPOT{X) for A*(x) = Ky/x 130 C 131 C l Upori^k) Calculations 131 C. 2 AL^por( x*)p Calculations 133 L i s t o f F i g u r e s 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 VII Grid Refinement - Effect of xN 86 6.2 VII Grid Refinement - Effects of x 0 88 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 Rh 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 Rh 99 6.10 T E A C H - Pressure Coefficient - Case 2 100 6.11 Reattachment Lengths 102 6.12 W a l l Pressure Comparison 104 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 . . . . I l l 6.17 Shear Stress Coefficient Comparison 112 N o m e n c l a t u r e Symbols a ellipse area constant an, as, ae, aw coefficients of finite difference equation a, b, c, d coefficients of cubic in ^-space A, B, C, D coefficients of cubic in x-space AE, As areas of ellipse and square A velocity profile coefficient Cn, Dn flux coefficients Cf shear stress coefficient C i , . . . , C-j integral-momentum coefficients IU-'-IIN V'-space cubics T\,-.-,7N x-space cubics 7o, ?E, 9B zero pressure gradient curves h step height hi,..., hN cubic endpoint values tU. viscous perturbation coefficient I, J bounding velocity functions LR reattachment length p pressure P e Peclet number (= puAx/T) Q cubic expansion ratio r ellipse and square area parameter r* solution value for r Rh Reynolds number (= U^h/v) Sp, source terms u, v finite-difference velocities UE bounding velocity UEO undisturbed velocity U'E velocity disturbance x, x, y cartesian coordinates viii NOMENCLATURE ix x, step location xR reattachment coordinate x 0, xN interaction region endpoint locations x e equation switch location x' virtual origin a Pohlhausen curve coefficient T Diffusion coefficient for 4> 6, A boundary layer thickness A / recirculation bubble height A 0 viscous boundary coordinate 61, A i displacement line 6*, A* displacement thickness AX cell width AXI, ..., AXN 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 c k n o w l e d g e m e n t s 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 so-lutions 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 differ-ential 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 com-putation 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 tech-niques. External boundary layer flows are a good example, as the large computa-tional 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 using finite-difference techniques for the viscous regions and panel methods for the inviscid core. It has also been shown, for fully developed flow, that the vis-cous 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 'building-block' 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 of flows are 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 RECION 13 V/SCOUS RECION A BOUNDARY LAYER EDCE B DISPLACEMENT THICKNESS C SEPARATION STREAMLINE D ZERO VELOCITY LINE Figure 1.2 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 Literature Review 1.2.1 Theory Briley and M c Donald [l] have conducted an excellent survey of VII theory. They give reference to the validation of VII solutions for several flow situations, including the BFS. 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 an-alytic solution, valid near separation, and singular at separation so that numerical methods could be employed through separation. Williams [2] describes how solu-tions 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 au-thors 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 singularities1 from the VI solution, the problem (iii) of devel-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 suc-1 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 re-view of the literature it became apparent that the interaction procedures developed for VII methods needed improvement. 1.2.2 E x p e r i m e n t s 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 laminar flow over 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, UEIN(X). 10 CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 11 <5* DIRECT I N V I S C I D S O L V E R INVERSE INVERSE v i s c o u s S O L V E R DIRECT 6* I N V I S C I D S O L V E R V I S C O U S S O L V E R CORRECTION FORMULA FOR 6" 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 UEVIS[X). 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 equa-tions 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 singu-larity 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 conver-gence 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 unsta-ble 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 exter-nal 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 formu-lae 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 I r E L L I P S E S Q U A R E Figure 2.2 Ellipse and Square The technique will be applied to a non-fluid dynamic problem so that the equa-tions 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, AE, and the area of the square, As, will be equated through variation of r using the proposed technique. CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 14 From analytic geometry the corresponding areas are, AE = irr(r + a) — irr2 + ixra (2.1) AS = 4r2 (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 = 7r(r + A r ) 2 + ixa(r + A r ) = [7rr(r + a)] + [7r(2r+ a ) A r + 7 r ( A r ) 2 ] AS+AAS = 4 ( r + A r ) 2 = [4r2] + [8rAr + 4 ( A r ) 2 ] Neglecting powers of A r greater than unity and subtracting the areas from both sides, AAE « KEAT A As « KsAr where, KE = 7r(2r + a) (2.3) Ks = 8r (2.4) 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 = As, is obtained using the derived perturbation parameters. 1. Guess r 2. Compute AE, AS, KE, KS from equations 2.1, 2.2, 2.3, 2.4 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: AE + A A E = As + A A S or AE + KE&T = As + KsAr Rearranging and substituting for AE, AS, KE, and Ks, _ AE-AS _ ?rr(r + a) - 4r2 A f " l T 5 - t f j s ~ 8r-jr (2r + a) 1 J 5. Add A r to r, ie: r n e u ) = rold + 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 AE As 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 PROCEDURE 17 ITERATION r AE As A r 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 AE are obtained. From equation 2.5, A r —> oo as 8r—7r(2r+a) —+0, or r = r c n t = \ 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 non-linear 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. ANEW VISCOUS-INVISCID PROCEDURE 18 2.3 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 AE and As with the viscous and inviscid boundary layer edge velocity distributions, UEvis(x) and UElN(x). 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 hn (the value of 6*(x) at n different locations say), then these hn can be perturbed by A/I„ in an analogous fashion to yield the linearized 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 KE 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 num-ber of independent parameters, hn. [Step 2] Calculate the perturbation function A6*(X) in terms of the hn and Ahn. CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 19 [Step 3] Linearize the perturbation function with respect to the Ahn in order to calculate the influence coefficients necessary for solution. Devising *^(x) If 6*(x) is a simple function like 8*(x) = h\x + h2, hi and hi being the perturba-tion parameters (analogous to r), the perturbation function and the corresponding influence coefficients are quite simple to calculate, ie: 8* + A8* = (hi + Ahi) x + (h2 + Ah2) = (hxx + h2) + (A/IX x + A/I2) Subtracting 8* from both sides, the influence function, A8*(X) = XAhi + A/12, is obtained. This function is already linear in Ahi and Ah2 so by inspection the linear influence coefficients for the parameters hi and h2 are x and 1 respectively. A function like hix + h2 is very simple to use but unfortunately does not have enough 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) = Jn(x) (xn-i<x<xn) (2.6) where, /„ (x) = An + Bnx + Cnx2 + Dnxz (2.7) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 20 h h h n h n-t Figure 2.3 Spliced Cubics The cubics are constructed so as to have continuous slope and value at their end-points. This gives four boundary conditions for each cubic, ie: The coefficients An ...Dn of equation 2.7 are completely determined if the slope and value at each end of the curve 7n(x) are supplied. The values at each end of the curve (see figure 2.3) are given by hn and hn-\. 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: > (2.8) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 21 4 similarly, where, Figure 2.4 Endpoint Slopes d 7( \ 1 A X „ _ i AXy = X,- — X , _ i AX„ j'=n-l,n,n+l (2.10) Now 7n(x) (and thus 6*(x)) can be expressed entirely as a function of / i n _ 2 , fen-i> hn, and / i n + i . Expressing 6*(x) in this manner makes the hn the independent parameters (analogous to r) rather than the An... Dn of equation 2.7. This is done to reduce necessary algebra and give more intuitive 'feel' to the problem. To reduce subsequent algebra further, a linear normalizing transformation is introduced, x — Xn-i A X „ Cn-1 < x< xn (0 < rp < 1) (2.11) 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 7n[x) into CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 22 cubic polynomials, fn(tp) , where fnW = an + bn ip + cn ^ + dn t/j3 (2.12) with the relations to J„(x), 7n(x) (2.13) From equations 2.6 and 2.13, Xn-1 < X < Xn (0 < %p < 1) The next step is to write fn(if>) as a function of / i „ _ 2 ... fen+i- The value at the left endpoint is fe„_i and the right is hn. The slopes at each end are given by equations 2.9 and 2.10. These give the following four boundary conditions (in ip-space) sufficient to determine fn{4>) completely. where [1] fn{0) = K_ 1 = a« [2] fn(l) = K = an + bn + cn + dn [3] / > ) = A Z „ — J(xn-i) = bn • i l Qn ( f e n - l ~ fe„-2) + (^ - fen-l) W = Axn-^7(xn) = bn + 2cn + 3dn 1 2 (fen+1 fen) + (fen fen-l) and Q n = A X n _ i CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 23 Solving the above for an ... dn , O n = ^ n - l bn = § [ ( A» - / l „ - l ) + Q n ( / l » - l - * n - 2 ) ] C « = | ( ^ n ~ fen-l) - 2 Q I + I ( ^ " + 1 ~ ^n) ~ Qn[^n-1 ~ hn-2) dn = -(hn - hn-i) + t^{hn+l ~ hn) + ^ (hn-i - hn-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 hn. Equation 2.14 needs to be modified 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 perturbation function A6*(X) Since / n (t/ ') represents 6*(x) and is a function of the hn, perturbing the hn in equation 2.14 gives the perturbation function t\fn(ip) and thus A6*(X) . This method produces what are termed local perturbation curves. This means that if a single hn 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 (2.14) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 24 Figure 2.5 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 pertur-bation 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 PROCEDURE 25 Figure 2.6 Modified Perturbation Curve is being perturbed. The perturbation functions, A / p and A / P + I , can be thought of as the incremental changes to the cubics fp and / p + 1 when hp is given a dis-placement A/I p. The perturbation curves are formulated as cubic polynomials with coefficients which are the linear increments to the coefficients of the original cubics. The perturbation curves are given as, A/p(V0 = AO P + A6p V> + AC P 02 + Ad p V'8 A / P + 1 (0) = A a p + 1 + A6 p + 1 rp + A C P + 1 rp2 + A d p + 1 V s If <5*(x) is to remain unchanged outside the perturbation region (xp-i < x < x p +i) and continuous everywhere, the slope and value must not change at x p_i and x p + i . These conditions give four boundary conditions, A/p(0) = 0 A / > ) = 0 A / p + 1 (1) = 0 (2.15) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 26 Continuity of A<5*(X) in the perturbation region at xp gives the additional three constraints, A / P ( l ) A/p + 1 (0) A /I(l) These seven constraint equations are insufficient to determine the eight coefficients of A / p and 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 CP being chosen as the independent parameter although any of the eight coefficients in equations 2.15 would have sufficed. dAC. [^(A/p)2 dxP + Qp+1 /oWp+1)2 W = 0 The factor Qp+\ is different from 1 if the two perturbation cubics are of different lengths. Using the seven constraint equations to write A OP ... A d p + i in terms of A C P , substituting into the equation above and solving, A CE Q 2 + 1 ( l l + 6 Q p + 1 ) - 5 2(1 + Q P 3 + 1) x A / I „ (2.16) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 27 The parameter cp is now a linear function of Ahp. Next the other seven coefficients are written as linear functions of &hp using equation 2.16 and the original seven constraint equations. Adp - 0 Abp = 0 Adp = Ahp — ACp Adp+i = Ahp = Q p + l (3a/lp — ACp) = -2Ahp (1 + 2 Q p + i ) + 2 Q P + 1AC P = Ah, (2 + 3 Q p + i ) - Qp+iAcp (2.17) The functions Afp(ip) and A/P+1(T/>) and hence (because the perturbation is local) A6*(X) are now completely determined by A/Ip alone. The preceding formulae will 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 Ahp are easily obtainable. For example, the variable cp has a linear influence coefficient, Kcp, due to Ahp and is given by the equation ACp = Kcp Ahp. From equation 2.16 by inspection, Kcr Q2p+1 ( l l + 6 Q P + 1 ) - 5 2(1 + QUi) (2.18) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 28 Similarly the influence coefficients for the variables ap...dp+1 are obtained from equations 2.17. Further, the influence function of /P(V>), corresponding to the per-turbation Ahp , can be found using the equation Afp(tp) — [ Kfp(ip) ] Ahp , A/p(0) being the linear influence function. From equations 2.15 and 2.17 , A/P(V>) = ACp 02 + [Ahp - ACp) ip3 or, since ACP = Kcp Ahp Afp(i>) = [l + KcpW2-^)} X Ahp and therefore, Kfp{1>) = 1 + Kcp^2-^) Similarly, Kfp+1ty) = 1 - 302 + 20s + Qp+i[ 3(0 - 202 + 0s) - Kcp(^ + 202 _ ^ ) ] (2.20) Knowing that Ahp only affects the cubics / p and / p + 1 , it is clear that all influence functions Kfn{ip) corresponding to the perturbation Ahp are zero unless n = p, p+1. 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. (2.19) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 29 2.4 The Governing Equations The proposed semi-inverse procedure requires that the interface velocity distribu-tions UEVIS [X) and UE,N (X) be matched at a discrete number of points in the solution 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 UEVIS(X). In the invis-cid region it is assumed that the external flow UEJN (X) will not be greatly disturbed by the presence of a small separation region so that linearized potential theory is thought to be applicable. 2.4.1 T h e I n v i s c i d R e g i o n From linearized potential theory (ref 1, pg.147), where UE,OIN{X) is the undisturbed velocity distribution and UjN(x) is the pertur-(2.21) bation to this distribution. In the present study only problems with UEpIN(x) — Uoo , the free stream condition, are considered. The disturbance velocity, UjN(x) , is due to the presence of the b.l . and is given in incompressible linearized theory by, (2.22) Now UElN(x) = U00 + U'IN(X) therefore, 4-{UElN{x) 6*(X)} = (Ulx + U^) d6*_ dx dU' 'IN + dx CHAPTER 2. ANEW VISCOUS-INVISCID PROCEDURE 30 d6*_ dx Substituting into equation 2.22, Substituting into equation 2.21, UE,AX) -I TT J-c di u„ = i + i / - + 0 ° 5e t 0 7T ./-oo X — Introducing the non-dimensional variables, x = — UPOT{*) U P O T (x) A'(x) x h 6*{x) (2.23) where h is a convenient fixed length scale of the problem, equation 2.21 is non-dimensionalized to give, UPOTW = 1 + 1 r + c o ^ { A * ( Q } IT J-C/ P o r(x) (2.24) 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. A NEW VISCOUS-INVISCID PROCEDURE 31 2.4.2 The Viscous Region The governing equations for the viscous region are the classical boundary layer equations for steady flow found in Schlichting [15]. du du dUE 1 dr , . x-momentum u——V v— = UE — 1— — (2.25) dx dy dx p ay du dv . continuity — + — = 0 (2.26) ax ay ' The pressure term on the R.H.S of equation 2.25 has been replaced by utilizing Bernoulli's equation for steady flow, _\&V_ _ JJ dUE p dx E dx where UE is synonymous with the boundary layer edge velocity distribution, UEvis (x). Equations 2.25 and 2.26 can be integrated across the boundary layer (in unsep-arated flow) thus, rs , du du dUn rs du du dUE TW JO ox ay ax 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 > + ^ f = 7 (2'27) where, roo ii " 6- = /„ p - ^ * Jo UE UE (2.28) 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 sep-aration 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 = UE. The boundary conditions for equation 2.29 are: at y = 6 u = UE dmu = 0 m=l,2...T-3 dy at y — 0 u = 0 d2u = dUE dy2 E dx dzu dy-> = ° ( T > 4 ) Two values of T were chosen in the model development, T = 4,5. Using the transform variable, - U. Js^ . du _ ^ du 6 dr) dy and solving for the coefficients At of equation 2.29 from the boundary condtions, the following equation for ^ is obtained. = FT{v) + A GT{r,) (2.30) UE CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 33 where, = 62 dUE v dx (2.31) Here T = 4 for the fourth order profile and T = 5 for the fifth order profile. The functions FT and Gx are given as, F,{n) = 2 (n - ns + \nA) G*{l) = | ( r? -2r ? 2 + 2r ? 4 - r ? 5 ) Solving for 6* and 5 from equations 2.28 using equations 2.30 and 2.32, or, similarly, 6* — _ -\- r C 2 A c 0 6 TGS + A + TCS A2 where the coefficients TC\. . .jC^ are given as, r C i = + Al - F r f a ) ] dr? * o r C 2 = - f GT(n) drj Jo TC3 = + [l FT{r})[l-FT{ri)] dr) Jo TC4 = + / 1 G r (r;)[ l-2Fr(r?)] dr; J O r C 5 = - f1 G2n(n) dn Jo (2.32) (2.33) (2.34) CHAPTER 2. A NEW VISCOUS-INVISCID PROCEDURE 34 The following table gives the values of TC\ .. .TC5 for the fourth and fifth order profiles. Fourth Order Profile (r=4) Fifth Order Profile (T = 5) TC\ 3/10 1/3 TC2 -1/120 -1/60 TCz 37/315 775/6237 TC\ -1/945 -25/16632 TCB -1/9072 -47/110880 Next, the shear stress parameter, rw/p, will be calculated with the aid of the equa-tion, rw = H du dy H du r?=0 Using equation 2.30 this can be written, v UE P dFr . dGr + A — — drj dt] v UE [TOQ + A] (2.35) n=0 The coefficients yC 6 and rCV are tabulated below. (r=4) (T = 5) TCQ TCn 2 1/6 5/3 1/4 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 A 0 0 h h x X h Equation 2.27 is to be solved on the computer, using a Runge-Kutta algorithm. The unknown variables UBL, A, 0, and A are solved by the following system of equations 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 dUBL [C6 + ( C 7 - C 1 - 2 C 3 ) A - ( C 2 + 2C4) A 2 - 2 C 5 A 3 ] RHUBLA A dx RHA* dA dx (2.36) 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 equa-tions 2.37 below. A* _ = d + C 2 A £ = C s + C 4 A + C 5 A 2 The shear stress coefficient is found from equation 2.35. 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 UEVIS obtained are matched at corresponding locations to values of UEIN (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. (2.37) (2.38) Chapter 3 A P P L I C A T I O N T O T H E B L A S I U S 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 sim-ple 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 or-der 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 approxima-tion 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. APPLICATION TO THE BLASIUS PROBLEM 38 3.1 Problem Geometry Figure 3.1 shows the geometry of the Blasius problem with the displacement thick-ness 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, 0 x < 0 A*(x) = -J.(x) 0 < x < x„ X n -1 < X < X n Xjv < X < O O (3.1) The region x„ < x < xN 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 70(x) and J (^x) are fitted smoothly into each end of 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 B o u n d a r y Cond i t i ons The viscous equations are examined first to determine the initial values (at x = x0) 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 C i > (3.2) 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, CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 40 where, 2 CQ CJ 2 C 3 The solution to this simple differential equation is, A*(x) = a x-C UBL Rh (3.4) 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 gen-erate the curves J0(x) and 7E(x), as the second order solution should not be much 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 70(x) as, The external velocity, UEVIS(X), will be assumed to be very nearly over the entire length of T0(x), so little error is introduced by setting 1*BL(X) = 1. 70(x) then becomes, 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, (3.5) CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 41 Substituting into equations 3.2, the initial values of the independent variables are, A(x0) 0 1 A(x0) (3.6) 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 70(>Cj the origin was known (x = 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, It should be pointed out that the initial and final spliced cubics must fit smoothly into the curves 70(x) and 3*E(X) respectively. This constrains A* somewhat, but it is necessary for a proper solution. 3.1.2 A s y m p t o t i c L i m i t s 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 UBL(x) (3.7) 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 con-ditions 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 equa-tions (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 inter-action zone A*(x) = 70(x). If the function 70 is chosen correctly, the conditions CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 43 upstream of the interaction zone will remain near to their freestream values. The form of 70 chosen, namely <*\Jj[^, has the desired properties. This particular func-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 J0 is an exact solution to the 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 J0 because it was derived from the fact that A = A' = 0. The solution to UPOT{*) 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,h2,...,hN, these being the values of A* at the endpoints of the spliced cubic polynomials. In order to solve for the hn, UBL(x) and UpOT(x) must be matched at N different values of x. These values 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 xk per cubic and in the same relative location within each cubic, the 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 xk are chosen to correspond with the endpoints of each cubic, for reasons which follow. Figure 3.2 shows the portion of A*(x) where the last cubic, fN[ip), splices into 7E{*)- Two possible locations of x* within fN are X _ * I Figure 3.2 Matching Point Locations shown. Locating xk at © ensures that the external velocities match at x = xN, while 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 xk, causes a 'kink' in the solution for A*(x). No kink was observed when xk was moved to ® . Results of numerical experiments showed that the kink became smaller if (i) the last cubic was shortened, or (ii) the last xk was moved closer to the end of the last cubic. When the last xk is 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,...,hN 2) Calculate an,..., dn from equations 2.14 3) Calculate Upori^k) and UBi(xk) from equations 3.9, 3.10 and 3.11 (to be de-rived). 4) If UP0T(xk) = UBL{xk) for k = l,N Stop 5) Calculate the linearized perturbation functions AUPOT[XIC) and /\UBi,(xk) and thus obtain the influence functions KUpor{^k)P and K UBL{xk)p given by equa-tions 3.22, 3.23 and 3.25 ,and 3.26 (to be derived). 6) Calculate the perturbation parameters Ahp from the (N) equations. N UBL{xk)-UPOT{xk) = 52[{KUpoT{xt)p-KUBL{xk)p} Ahp] (3.8) 7) Update the hp, ie: hPnew = hPoU + &hp 8) Goto step 2) 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 V e l o c i t y Ca lcu l a t i ons Calculation of Upoxjxk) From equation 2.24, dA' UpoT^k) = 1 + l r+co IT - / - ^ d t TT J-oo Xk — Z 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 TT Jx where, dt + 1 N 1 r IT ^ A X „ Jo 1 bn + 2c„ e + 3dn ? ipk-t di xN (xjfc - $) y/t-X1 and 0jb = x* — X n - x A X „ 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 PROBLEM 47 UPOT{X>C) = 1 + J0{X-k) + JEM + Ik,k+1 Jfc - 1 N n=l n=*+2 (3.9) k = n N-l IN,N+I n=l (3.10) All the functions (J0, J N , J E , Ik,k+iy and IN,N+I) a r e given in Appendix C, section CI. C a l c u l a t i o n o f UBL(xk) In the interaction region A* consists of the spliced cubics fn{ip)- Thus it is convenient 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 (x„_! < x < x n) dx Axn dtp 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 UBL> s o that, £ J r1 dUBL J ( , rt* dUBL dip Jo dtp where UBL\(tpk) = UBL(xk). Substituting for and noting all matching points are at the cubic endpoints (rpk = l), CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 48 Equation 3.11 is symbolic as UBL 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 P e r t u r b a t i o n Ca lcu la t ions Perturbations to 5* (A6*) yield perturbations to UE (&UE) for each region. Before any calculations of &UB can be made, the perturbation curve for the case p = N must be developed, as discussed on page 27. Figure 3.3 shows the perturbation curves when hN is given a small displacement AhN. The boundary conditions for Figure 3.3 Perturbation Curve for n = N A/ W(V0 = A a N + A&N V* + A C N V"2 + A d N tp3 are, A/*(0) 0 0 A/„(l) CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 49 Solving for ACLN . . . AdN, AaN - 0 AbN = 0 ACN = 3 A h N - AXN ^ A J £ ( X ) | X W -2AhN + A X N ^ A 7 £ ( X ) \ X N AdN = (3.12) The curve downstream of the interaction region is altered when hN 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 = xa1 gE{x) = a Rh (3.13) When the curve 7E{X) 1 S perturbed, its linear perturbation curve ATE{*) is equiv-alent to the linearized difference between £ .E(X) and J E ( X ) . The only variable pa-rameter in 7E{X) is x ' s o that it becomes the perturbation parameter, ie: 7E(x) + A7e{X) = a JX - (x' + AX') Rh a Using the binomial theorem, / A X ' '1 : » X — X' RH A X ' x — x' A X x — X' 2(x - x') This gives the linearized perturbation function as, A J£( X ) = a AX 2^]Rh (X-X') CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 50 The term A I ' can be written in terms of AhN by using the boundary condition A ^ X , , ) = AhN, ie: (3.14) Because A J E ( X ) is the linearized difference between §E(*) and 7E(X), a linear esti-mate of x' G comes from the equation, J £ ( x N ) + A / I „ = 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, A a N = 0 A 6 n = 0 A C „ = (3-- Cl) AhN AdN = (n - 2) AhN (3.16) where, A X , 2 ( x „ - x ' ) (3.17) Following along similar lines to the influence coefficient calculations of Chapter 2, the perturbation coefficients and functions corresponding to AhN are, KcN = 3 - n KdN = 0 — 2 (3.18) CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 51 with Kf„{1>) = KcN^2 + KdN ips and (3.19) Having now derived the case p = N, calculations for steps 5) and 6) can be com-pleted. Calculation of the Inviscid Velocity Perturbation &UPQT[XIC) To obtain AUporix-k), equation 2.24 for UPOT{*IC) is perturbed, ie: TT f \ ± TT ( \ 1 -u 1 r°° ^ [ ^ ( 0+AA ' ( Q ] UpoTKXk) + &UPOT(xk) = 1 + - / -s 7T J-oo Xfc — Q Separating out the perturbation terms and noting that the perturbation curve AA* is zero if x < x0, AUPOTM = - / ~ — di * J*o Xfc - 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 = - 2 2 — 1 " ^ 7 J- dC + ~ / T — d£ 3 3.20 7T n = p A X „ JO 1pk- Z 7T Jx„ Xfc - C where A£/por(2;fc)p is the perturbation velocity at x = x* due to perturbation of a single h, namely hp at x = xp. For clarity, the list below describes the meaning of § This term is zero if p ^ n 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 ixk)i due to all perturbed /i's is given by the sum 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 p. In each case it is important to note that the only contribution to the integral for At/por(xjt) comes from the two cubic perturbation curves surrounding the hp. The perturbation cubics 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 KUPOT(^k)p involves some algebra and so is moved to Appendix C, 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 PROBLEM 53 Figure 3.4 p ^ JV Subcase Locations summarized in equation 3.22. KUpoT{xk)P — i [ A J p ( 0 t ) + AJ p + 1 (0 t ) ] A X , + 1 m V f + i A X A X p + i _ \Kdr+i+2Kcv A X P + 1 case (T) case (?) case (3) case (T) (3.22) The functions A JP and A JR+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 1 J[-J[ *«MV»*) + &JB{*k)] case © AJ s (x t J * — case (2) AJ^X* ) ] case © (3.23) The functions JN(ipk), AJB(xk), and A J ^ ( x t ) are given in Appendix C, section C2. 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, dUBL A X N A X „ - i < x < X„ drp Rh A 2 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, dUBL = _ dxjj C2 Rh A 2 This equation is then perturbed giving, [A* „ 1 AX„ d{UBL + AUBL) AXR A* + AA* . A + A A di> C2 Rh (A + A A ) 2 Because the perturbation A A is small compared with A, the expression can be written, d{UBL + AUBL) A X » A* + AA* dtp C2 Rh A 2 By this method of linearization, the unknown function A A is eliminated. Subtract-ing UBL from both sides gives the perturbation, d A XN AA* 1J{* Ubl) = cTihA-* For a single perturbation A/I„ at x = x„, X n - l < X < X„ (3.24) A A* = { X p - i < x < X„ KfP+IM x A/ip xp < x < x p + 1 where Kfp[xp) and Kfp+\(ijj) are calculated in section 2.3. Because the boundary layer equations are parabolic, any perturbations downstream of x* will not influence UBL (xjt) so from equation 3.24 and the above equation, 0 k<p-l KUBL(xk)P = ' AXp f l C2 Rh Jo g AXp k n=p C2Rh v Jo A 3 k>p (3.25) CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 56 Analogous to the inviscid example the perturbation function AUpor(x;t) is the sum of all the individual perturbations, or N p=l (3.26) 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/IP. 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 UBL{*) = ^por(x) = 1 is an exact solution using the present VISI method for the case A*(x) = ?0(x) for 0 < x < oo. This offers a convenient way of checking the accuracy of the present method. In essence the method approximates 70{x) using spliced cubics. 3.3.1 R a d i u s of Convergence and S t a b i l i t y 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 PROBLEM 57 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, For T = 4 this corresponds to A = 12,60 and for T = 5, A = f, *5Q. it i s only 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. dx dA when A CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM Figure 3.6 Iteration History CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 59 3.3.2 S o l u t i o n A c c u r a c y , G r i d Ref inement , and E x e c u t i o n 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. Fig-ure 3.7 shows how these values are used to assess grid refinement. Grid refinement was performed by successively dividing the interaction region into 2N equal length cubics. The accuracy of the method is quite good. The value of x' is very sensitive to the value of A* at xN. Despite this, and using only 4 cubics, the value found for x' is less than 1. This is surprising considering that x' is found by extrapolation from xN (1000 in this case). For example, a 1% reduction in the solution of A* at xN 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 PROBLEM (x 0=iq xN=iooo R h=ioo) ( i-uB L)t 1 0 " —i i o - H io" H io' H io" H i o " -H io- H + + + + 16 32 N Figure 3.7 Blasius Grid Refinement CHAPTER 3. APPLICATION TO THE BLASIUS PROBLEM 61 It is noted that the program has not been optimized for speed and could con-ceivably 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 A P P L I C A T I O N T O A B A C K W A R D F A C I N G S T E P In the previous chapter the VII method was applied to a simple situation of unsep-arated 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 de-veloped 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. APPLICATION TO A BACKWARD FACING STEP 63 4.1 Problem Geometry 1 A L I r 1 : L . i i m X L„ -j x n x„ Figure 4.1 BFS - Problem Geometry Figure 4.1 shows the geometry used in modelling the BFS 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. All quantities are dimensionless, having been scaled with respect to the step height, h. x0 Distance from origin of b.l. to beginning of interaction region. xN => Distance from origin of b.l. to end of interaction region. x s => Distance from origin of b.l. to step. CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 64 x R =>• Distance from origin of b.l. to reattachment. LR Reattachment length. A,(x) =r* Height of recirculation region. A c(x) =>• Distance from x-axis to outer edge of b.l. A*(x) Displacement thickness. A(x) Boundary layer thickness. A = A 0 — A , Ai(x) => Displacement line. A x = A , + A* The effect of the size of the interaction region is not known a priori, but can be determined through variation of the parameters x 0 and xN. 4.1.1 T h e Disp lacement L i n e 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 70(x). Because the upper wall is displaced a unit distance from the x-axis, the function 70(x) must also be displaced a unit distance, ie: 1 x < 0 A x(x) = -= 1 + Q\/iE 0 < x < x0 = o„ + bnip + c„V2 + d„03 x n- l ^ x < x n Ix - x' = a V * XN < X < OO (4.1) CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 65 — • — ^ — X0 X 1 2 ' / / / / X n XN Figure 4.2 The Displacement Line All asymptotic properties calculated in Chapter 3 hold for the BFS problem, as shown in Appendices A and B. 4.2 Modifications to the Governing Equations 4.2.1 The Inviscid Equation The perturbation UN(x) of equation 2.21 can be thought of as the perturbation to 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*, UPOT{x) = 1 + 1 f + o o ^ A x U ) } ix J-UP0T(x) (4.2) 4.2.2 T h e V i s c o u s Equa t ions Figure 4.3 shows the boundary layer divided along x into three separate regions. Regions VI 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 VI, 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 'FLARE' 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, where C = 1 for u > 0 and C = 0 for u < 0. In the viscous model used here the 'FLARE' 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 dUE 1 dr dx p dy CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 68 Figure 4A 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\60 — 6,} (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), with, A* A 0 A X n [ C 6 + C 7 A ] ^ = C3 + C 4 A + C 5 A 2 A = RH A 2 dUBL dtp (4.3) (4.4) and where the coefficients Cy ... C^ retain their previous values. In the 'dead air' (reverse flow) region, with u = v = 0, the integrated x-momentum equation becomes, Jo dx p Jo Si dr dy dy Integrating, TT D U E X 1 I \ ax p (4.5) where TW is the shear stress at the wall and T> is the shear stress at y = 6,. As a 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 rw will have to be developed. Writing the laminar form for r 7, T{ = p. du dy ix du 6 dn fx UE [ C 6 + C 7 A ] n=o CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 70 Substituting into equation 4.5, non-dimensionalizing, and transforming into ip-space, duBL_Ai = _ ^ [ C 6 + C t A ] (4_6) with, dip ' Rh A A , - A x - 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 A X „ [C 6+ (Cf - C : - 2 C 3 ) A — (C 2 + 2C4) A 2 - 2 C 5 A 3 ] dip RHUBL' dUBL A X n A dip Rh A 2 dAi / A A d0 dA dip Ve) dip dip A c +Ce A 2 ~ ^ (C 4 + 2C5A) dA dip A 0 dQ . . _ ^„ ,. dA - - A ( C 4 + 2 C 5 A ) -(4.8) 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 reattach-ment, 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 r c 6 + c, (4.9) so that for A , = 1, A = - (4.10) . C 6 + C 7 AJ 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 VI 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 B from equation 4.10, A : - A , C\ + C2 At . C 6 + Cy A £ Solving for A B in terms of (known) A*, \C1 + C7A*A) 2 Co 4C, 2 (4.11) Here the 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. Equa-tion 4.10 is then used to calculate A B using the now known A B . CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 72 Also, from the second of equations 4.4, 0 B = A B [Cs + C 4 A f l + C 5 A 2] (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 discon-tinuous 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 = xR). Analysis of equations 4.8 shows that A' becomes infinite at this point. 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 7 has a sufficiently small value, e. At this point the viscous 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 D to its reattachment value. With reattachment located by the point CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 73 A, ^ ^ ^ ^ A, tJ/)))}/}/)}})) X t x „ Figure 4.5 Reattachment Singularity of zero shear stress, equation 2.38 yields, A - C e Substituting equations 4.14 and 4.15 into the first of equations 4.4, From equations 4.15 and 4.16, and the second of equations 4.4, (4.15) (4.16) (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 Parameter Calculations 4.3.1 V e l o c i t y Ca lcu l a t i ons 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 UBi is written, 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 V e l o c i t y P e r t u r b a t i o n Ca lcu la t ions 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 AL7Bi/(xfc) is given by equations 3.25 and 3.26. For region VI, A* is replaced by A x — 1 so that because 1 is a constant the perturbation does not change, and again 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 AAi only. From the second of equations 4.8, dUBL _ A X n A d%l> ~ RhA2 perturbing this equation, CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 76 Extracting the perturbation, I TT \ AX N dip Rh A 2 AA (4.18) The next step is to find AA in terms of AAJ . Combining equations 4.6 and 4.7 with the first of equations 4.8, A i " C 2 A 2 + C i — Cj — A - C « = 0 Perturbing, C 2 ( A + A A ) 2 + C\ — Cj — A X + AAX (A + AA) = 0 A + AA Neglecting higher order terms and solving for AA in terms of AA 1 5 A AA = but, 2C 2 A + C\ — Cj A x 2C? A + C\ — Cn — = A A A x C2 + AAI A 2 J Substituting the last two equations into equation 4.18, d V 7 ( A ^ ) = AXR c 2 + A 2 J AAX Finally, equation 3.25 from the Blasius problem can be modified to give, KUBL{xk)f AXP C2 Rh p+i n=p 02Rh JO A ! X(A) Kfn(iP) k<p-l k=P dip k>p (4.19) (4.20) CHAPTER 4. APPLICATION TO A BACKWARD FACING STEP 77 where 1 x0<x<xs or xR<x<xN C& 1 n -1 1 + C 2 A 2 J xs<x<xR After computing the influence coefficients, KUsLi^k^pi from equation 4.20, equation 3.26 is used to solve for AUB^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 T E A C H VII methods are used to calculate approximate solutions to the Navier-Stokes equa-tions. More exact solutions are found through numerical integration of these equa-tions. 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 TEACH-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 s* Mass l 0 0 x-momentum u dp dx ,d2u + fl{dx> + d2v dxdy y-momentum V dy d2u dxdy 5.1.1 C o n t r o l V o l u m e F o r m u l a t i o n i t -Figure 5.1 Scalar Cell The principle of the TEACH-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, 5.1.2 H y b r i d Differencing and L i n e a r i z a t i o n Following Patankar, equation 5.2 is discretized using hybrid differencing, a combi-nation 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 differ-encing alone leads to numerical instabilities when the Peclet number, PT = pu AX/T is greater than 2. The source term is linearized using central differencing to yield, Using equation 5.3 and the hybrid differencing scheme, equation 5.2 can be written, J Jcv / / Sfdxdy (5.3) (ap — Sp) (j)p = aN4>N + as <f>s + aB(f)E + aw 4>w + Su (5.4) where, ap = dN + as + aB + aw CHAPTER 5. TEACH 81 and for example, = max(|C n/2|, Dn) - Cn/2 where, D.-J Jx Xy, dx Central differencing is always used for the diffusive terms (ie. Dn) and upwind differencing is used in the convective terms (ie. Cn) if \Pe\ > 2. 5.1.3 S o l u t i o n P r o c e d u r e 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 iter-ation. The pressure and velocity fields are then corrected using an approximate pressure correction equation. The corrected pressure is used to calculate the veloci-ties 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. 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 layer flowing over 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 = x0). The velocity profile at the inlet is that of 5.2 Problem Geometry CHAPTER 5. TEACH 82 Figure 5.2 Boundary Conditions a zero pressure gradient Pohlhausen profile x0/h stepheights from the origin of the 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 S p e c i a l Ce l l s 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 P r e v i o u s 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 (Rn = 100) from 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 R E S U L T S A N D A N A L Y S I S 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 xe and R^. Results of the grid refinement study will be analysed in terms of the 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. RESULTS AND ANALYSIS 6.1.1 G r i d Refinement 86 Xs= 6 0 Rh= 1 0 0 100 (x-x0 )/h Figure 6.1 VII Grid Refinement - Effect of xN Figure 6.1 illustrates the effects of grid spacing (cubic length) and the position of xN 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 — x0)/h w 31. The results show that grid independence is realized for quite 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 xw. This verifies the hypothesis that the boundary layer returns to a Blasius profile (ie: A x becomes ^(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 xa. 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 J0(x) is not close enough to the real solution near the step. The dip does not appear to affect the solution downstream of the step. For example, the change in LR for the two curves shown is less than 1%. The effects of the dip on the pressure distribution are shown in figure 6.3. The hollow squares coincident with the zero Cp line are not calculated by the method, they come from 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 (x0 = 100) immediately compensates to reach the same pressure minimum at the step as the x0 = 50 case. Again as in figure 6.2 the effects are negligible downstream 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 CHAPTER 6. RESULTS AND ANALYSIS © w 6 + o o . © I a < O n " « • -X 0 *S R h T • 100 130 200 i 050 130 200 4 -BO -40 I 40 80 120 Figure 6.3 VII Grid Refinement - Pressure Distribution CHAPTER 6. RESULTS AND ANALYSIS 90 6.1.2 P a r a m e t r i c V a r i a t i o n The VII model for the BFS has two independent parameters, xe and Rh. 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 xs, but they are related through Rh and thus are interchangable. Some authors (eg. [7,12,14]) ignore or neglect the effect of x, on the flow down-stream of the step. Conversely, Armaly et al[l3] conclude that the reattachment length, Lr, is not a function of Rh alone, and is indeed a function of the parameters 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*a « 0.78 and LR « 12.6 while for xs = 350, A* « 2.32 and LR « 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 x become smaller as x, increases, and in consequence the 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 xs and T (the order of the velocity profile 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 6. RESULTS AND ANALYSIS C °--20 I 20 X 0 X S R h T D 300 350 200 4 020 040 200 4 4 0 I 60 60 (x-xs )/h Figure 6.4 Pressure Coefficient - Effect of x, CHAPTER 6. RESULTS AND ANALYSIS 6 • + o d O • w5e a • A 4, AO ft Xo * S R h T 080 110 100 5 • 040 060 100 5 • 020 040 100 5 A 080 110 100 4 Q 040 060 100 4 O 020 040 100 4 -20 I 20 I 40 (x-xs )/h 60 60 Figure 6.5 Pressure Coefficient - Effect of T 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 Cp gradient maximum. Figure 6.5 also shows the pressure 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 pa-rameter 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 Cp with Rh, with x4 held fixed. Upstream of the step the symbols show a definite independence of Reynolds number. This is quite remarkable considering that given the same values for x4, the values of A i at 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 x3. Values of x' ranged from 7 to 50 for xs ranging from 40 to 130. Since x' is positive in CHAPTER 6. RESULTS AND ANALYSIS (x-xs )/h Figure 6.6 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 Rh because x' remains nearly 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 TEACH 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 Ref inement 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 Rh=200 Y-EXPANSION 5:1 Figure 6.7 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 TW (ie: |^) cannot be verified as the grid is too coarse to discern the gradient very close to the wall. 6.2.3 P a r a m e t r i c V a r i a t i o n R « =400 Xo= 110 0 R« =400 Xo= 20 • =100 Xo= 100 =100 Xo= 40 * =100 Xo= 20 =100 Xo= 5 i~\ i i i i i 0 10 20 30 40 50 60 (x-x0 ) / h Figure 6.8 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 x0, so that varying x 0 is equivalent to varying x,. In figure 6.8 there are two cases of Rh, 100 and 400. The results are consistent with those of the VII method in that increasing x, decreases the pressure gradient and CHAPTER 6. RESULTS AND ANALYSIS 98 increases LR. 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 x0 = 5. This is due to the relatively coarse grid and thin boundary layer, 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 Cp vs Rh with x 0 fixed, analogous to figure 6.6. This plot 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 imper-meable 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 Cp on Rn is gone, and the adverse gradients are less than those of figure 6.9. This results in a reduction of 8% in LR for the Rh = 100 case. It was found that for the case of the free slip top wall that the u-velocities near the top wall were accelerated CHAPTER 6. RESULTS AND ANALYSIS i * * * t T • » » • * • R„=400 Xo= 20 • R„=200 Xo= 20 * R,=100 Xo= 20 K= 50 Xo= 20 10 r 20 30 I 40 SO 60 (x-x0)/h Figure 6.9 T E A C H - Pressure Coefficient - Effects of Rh CHAPTER 6. RESULTS AND ANALYSIS 100 O ° o c O 0> o o 9 9 9 J 8 f K _ . . e * 4 * ' ' Si OH 4 KE*100 * RE" SOC O RE" SOC X 10 20 30 40 50 60 70 X / h Figure 6.10 T E A C H - Pressure Coefficient - Case 2 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. RESULTS AND ANALYSIS 102 6.3.1 Rea t t achmen t 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 LR 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 xt for these experiments are not given as their upstream CHAPTER 6. RESULTS AND ANALYSIS 103 conditions were different. For Goldstein et al [12] the range of xe was from 4 to 30. 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 TEACH solution. Hackman [22] shows for a specific case (Rh « 300), LR is underpredicted by about 10% using hybrid differencing with a relatively fine mesh. The effects of xa are the same for both models although these effects are more 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 depen-dency of LR on x,, but only indicate that increasing x, increases LR. Their data is presented as a least squares average with the range of results shown by the scatter bar. The range of xe used in the experiment lies below that used in the models. This tends to make the agreement look worse than if comparable values were used. The present initial guess for A x(x) employed in the VII model is inadequate when high Rh and low xs are required simultaneously, and in consequence no results are 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 in figure 6.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 O t h e r Charac te r i s t i c s 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 pre-viously 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 reattach-es =40 Xo=20 Rh=200 n 6 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 LR also resulted. Next the outlet was moved 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 106 Figure 6.13 Refined W a l l Pressure Comparison 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 TEACH 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 109 OS o E— o w < as m 20 A"s=40 A' 0 =20 Rh-Z00 • • • • • • vii • TEACH V- • • • • 40 T" 60 60 100 X /1\ Figure 6.15 Shape Factor Comparison 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 FLARE 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 — xa)/LR 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 TEACH 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 rw at the point where TW is a minimum. From equation 4.5, with Tw m —r, 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 111 XS=40 XO=20 RH=200 Y/h o o o A O A A O TEACH SEPARATION STREAMLINE A TEACH ZERO VELOCITY UNE 10 T" 15 (x-xs )/h 20 Figure 6.16 Separation Streamline and Zero Velocity Line Comparison CHAPTER 6. RESULTS AND ANALYSIS 112 (x-xs )/h Figure 6.17 Shear Stress Coefficient Comparison CHAPTER 6. RESULTS AND ANALYSIS 113 either end. Some semi-empirical model for TW could be included in the VII model to improve the results. 6.3.3 E x e c u t i o n T imes 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 so-lutions. The convergence criteria were freely selected as [&h/h]MAX < 10~s for the VII model and, for the T E A C H , a residual error of < 10~s for the maximum of the 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 LR 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 TEACH 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 C O N C L U S I O N S A N D R E C O M M E N D A T I O N S 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 sim-ple manner. Again the assumptions are not necessary for the general technique or the integral model and could be improved if greater complexity was clearly neces-sary. 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 for flows with 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 inviscid flows is 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 ex-pected. 3. In all cases studied, values calculated for the virtual origin, x', of the redevel-oping 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 thick-ness 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 distri-bution of control points ( A X 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 pre-dictions 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 rep-resenting 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 and finite-difference methods, 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 TEACH, 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 FLARE [16] approximation. CHAPTER 7. CONCLUSIONS AND RECOMMENDATIONS 119 2. Neglecting the wall shear stress (rw) in the recirculation region in the VII model. Again this was done to construct as simple a model as possible, how-ever the results of the T E A C H code show that TW is not negligible in the 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 ve-locity 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 com-parisons more difficult. CHAPTER 7. CONCLUSIONS AND RECOMMENDATIONS 120 7.2 Recommendations 1. A theoretical or semi-empirical model of the wall shear stress in the recircu-lation 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 represent-ing 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 func-tion 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 dis-torted (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 TEACH-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. R e f e r e n c e s [l] Briley, R. and McDonald, 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., Pro-ceed. Vol.11, pp.161-168 (1973), Springer-Verlag [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 Sep-aration, 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 Publish-ing Co., New York (1980) [19] Abdullah, Z. A personal communication. (Jan. 1988) [20] Djilali, N. An Investigation of Two-Dimensional Flow Separation with Reat-tachment. 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 ' UcRh The solution to the b.l. equations as x —* oo will be examined. The variable Uc in the equation for A* will be shown to equal the constant free stream value of UBL at oo. From the second of equations 2.36: dUBL A dx. A 2 Rh (A.l) 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 + c2 A (A.2) Because A is bounded, it follows that A oc y/x as x —> co. Substituting into equation A . l , dV~BL A . — — a - (A.3 dx x 124 APPENDIX A. 125 Again, because A is bounded, it can be written as a constant plus terms which become negligible near oo. The value of A 0 corresponds to the asymptotic value of U B L . A = A0 + 0 (l) as x —• oo Substituting into equation A.3, Integrating, dUBL A . 0(1) —.— oc 1 — ax x x UBL OC yl0lnx + 0(0) Now, because U B L is bounded, A 0 = 0 at x = oo and therefore, A —• 0 as x —* oo Differentiating equation A.2, A A - A' A* A2" = c 2 A' (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), l i m A * X — o o A = a x — X UBL Rh This equation is identical to the asymptotic integration path except for the variable U. Because the above relation must be true, U B L must assume the constant value APPENDIX A. UC as the integration path relation must also hold, ie: UBL —> UE as x — • co So by choosing UC = 1 it is necessarily true that asymptotically UBL{*) constant and A = A' = 0. Q . E . D . Appendix B B.l Calculation of Asymptotic Value for c7p0r(x) In this section it will be shown that Upor(x) ~~+ 1 as x—•oo. The inviscid velocity distribution is given as, Ppor(x) = 1 + - / -* — d£ 7T J-oo 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 70(x), J„(x), and 7E(X). Upoi{x] = 1 + u*-mfdi + ^hr-mm^ TT JO X £ ^ TT JXn-! X-£ In i r+°° d^ 7T 7xw X — £ I*N x - ' IE It is easily seen from the above equation that if the sum of all integrals vanish as x —»• oo then UPOT 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 70 (see section 3.1.1) into the integral I0, j a fx° di _ a In Jo In \/x + 2n\fR~h y/x From the above equation in the limit, lim T X - » o o i0 = a Inl = 0 Next, substitute for Jn in In. 1 /*- (j3B + 2 C n £ + 3 0 n £ J ) I f TT Jxn-i di The limit of this equation will be taken while it is still in integral form, ie: x lmo o /n = — fXn {Bn + 2Cni + 3Dni2)di 7T X 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, 1 lim T X - o o ln TT C O [ A*(xn) - A*(xn_x) ] = 0 Finally it will be proven that in the limit the integral IE is zero. Care must be taken to evaluate the Cauchy principal part of the integral as there is a singularity at x = i. IE is written as the sum of two integrals thus, IE = a f x - £ lTT\/R~h J*N di + h J** (x - i) Vt-*' 2*VRh Jx+t (x - i) vi-x1 a rc 1T\/R~h Jx-di _7r\/x —; ln y/i^-y/x^x1 Li n X - £ + TT\f) In x —X' y/j^+Vx-X' y/i^-Vx-X' X+e APPENDIX B. 129 The lower limit of L\ is zero because Xoo » xN and so the In argument tends to 1. The upper limit of L2 is zero because as f —» oo the argument of the In function tends to 1 (here too ^ *oo ^ *N) and the term preceding the In function tends to ^ . This leaves the limits surrounding the singularity, or 1 IE = In 7Ty/x — X1 Multiplying out the terms and letting 7 = x — x', 1 (y/x-x'-e + y/x-x') (y/x-x' + e - y/x-x') (y/x-x'-e - y/x-x') (y/x-x' + e + y/x-x') In y/72 - 6 2 - y/7 2 + ^(y/T+f ~ y / T ^ ) v V - - V7 2 - y/7 ( y/7+e - y / l ~ € ) 7Ty/x — X 1 Using the binomial theorem and the fact 7 >^ e, \ / 7 2 ~ e2 « V/T 2" 1 e2 ^ - 2 7 and Substituting the above in equation B. l , 1 IE = 7 T y / * In 1 - 1 £ 2 t • 1 _ l i X 2 7 so that finally, Putting it all together, lim T T 1 \ , , lim T . lim r i l i m r X — 0 0 Up or (Xj = 1 + X—oo i 0 + X—oo 2^ i n + X — 0 0 i l i m l»m , Tf—»oo e—»0 IE — U l i m l i m TV l i m n=l = 1 + 0 + 0 + 0 = 1 (B.l) Q.E.D. APPENDIX B. 130 B.2 Calculation of UPOT(X) for A*(x) = Kyfi 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 J-oo X — £ di 1 + K_ r+° 2n Jo di Taking the Cauchy principal part of the integral, £ W ( x ) = 1 + K 27T V/X In X-e + K . 2 7 r\/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 1 + T ^ 7 = l n (y/x + \fx-e) {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 ( 2V^ " ife) Cancelling terms and letting e —• 0, ^POT(X) = 1 + K 2-Kyfx 1 + 0 ln 2y/x~ 2yjl Q.E.D Appendix C C l UpoT{*k) Calculations UpoT[*k) = 1 + - / 7T JO dt + (x* - 0 s/t a r°° dt 7T Jx + 1 " 1 r1bn + 2cnt + 3dnt2 - t — f 7T ~ AX„ Jo rpk-t 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;(xfc) is valid for x '<Xjfc } Uporixk) = 1 + 1" \Jx~k In \JXk~—\/Xo 7T y/xk - X' In y/xN - X' + y/xk - X ' y/xN - x'-y/Xk - X' Jo{xk) + -T — TT ~ A X „ (6n + 2cne+3d„e2) In J £ ( X j t ) 3dn(0jk + -) - 2 c „ (Cl) This equation is valid for any value of ipk except 0 or 1. The matching points were selected to coincide with the cubic endpoints, thus giving tpk = 0 when k — n—1. In 131 If x' > x 4 then JEM = 2o * y V — Xi - arctan \ APPENDIX C. 132 other words a singularity is created in equation C l at xk (between the two cubics fk 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. T 1 lim J-k,k+i = «i-»0 A X t 7 T 1 lim e 2-»0 A X i + 1 7 T {bk+2ck{l-E1) + 3dk{l-E1)2}\n {bk+1 + 2ck+1£2 + 3dk+1£l}\n £ 2 l-e2 -3dk{%-e1)-2ck 3dk+1{± + e2) - 2ck+l From continuity of slope and the scale factor between the cubics /* and fk+i, bk+i _ bk + 2ck + 3dk £ i _ A x t + 1 A X f c + i AXfc Substituting into and taking limits, E2 A X f c (C.2) For case 2) where the singularity (matching point) is at xN, 1-ei 1 lim e,-0 A X N 7 T {bN + 2cN{l-e1) + 3dN(l-e1)2}ln - 3dN{* - £l) - 2cN , 1 lim H — « 2->0 7T y/xN~Xf + E2 In \ / x N - X ' + \/xN-X' + E2 y/xN-X' - y/xN-x' + E2 From continuity of slope and the scaling factor between fN and 7E, 1 bN + 2cN + 3dN EX and — A X , E2 A X W APPENDIX C. 133 From the binomial theorem, y / x N - X ' + S2 ~ Vt-N-X' + Substituting into IN,N+1 and taking limits, £2 2y/xN-X< 7T y/xN— X' In 4 (xN - x') A X , (2c„ + §<f„) A X . (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, xp. 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, TT ~ A X „ J AUpoT{Xk)p — ^k-i (CA) KUporixkii where, from Chapter 2, Kfn{rP) = Kan + Kbni> + Kcni>l + Kdn^kz APPENDIX C. 134 There are four subcases to consider, labelled © . . . © . Subcase © : xk ^ x^ , x p, x p + 1 In this subcase equation C.4 contains no singularities and so can be integrated in a straightforward manner to give, KUpOT{*k)p = H A J p ( ^ J t ) + A J p + 1 ( 0 f c ) ] where, A X „ {Kbn + 2Kcnibk + 3Kdnibk2) ln 4k 1 -fbk - 3Kdn {rpk + i ) - 2Kcn Subcase © : xk = xp_, 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 0 r(x f c) p = «™o i r> j mo 4 k - i - M 7T A X „ Jt di Ah„ KUpQT ( X j t ) j : + 7T A X D + , J O di Ahr, p + 1 o 1pk~ i Substituting for Kfp(i) and Kfp+i(i), integrating and taking limits, 1 KUpoT{*-k)p = 7T \Kdp + 2Kcp AJp+1{4k) - 2 V- P-A X C Subcase ® : xk = xp For this subcase limits must be taken on both integrals of equation C.4, as the APPENDIX C. 135 singularity is between fp and /p+i, ie: A^por(xjfc), lim e—0 7T AXD Jo d£ KUpoT{xk)i lim + e->0 7T AX P + 1 ./« de A/I„ Integrating and taking limits using the fact that there is continuity of slope between fp and A / p + 1 at x = Xjt, KUpoT{*-k)p = Subcase ® : x* = x p 7T Kbp+1 ^ + \Kdp + 2Kcp ^ \Kdp+1 + 2if c p + 1 AX P+I AXD + AXP+, 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, KUpoT{*-k)p = •K A J P ( ^ * ) \Kdv+1 + 2 i f t p + 1 AX p = N In this case the perturbation curve extends from xN_t to oo. The equation for AUpOT(xk)p is, d zrt (t\ .„ d AUpoT{Xk)p — 7T AX N JO Ipk- < J*N AhN (C.5) KUpoT{*k)p For this case there are three subcases to consider. Subcase © : x t < xN_t APPENDIX C. 136 For this subcase xk is completely outside the perturbation region so that equa-tion C.5 contains no singularities and can be integrated in a straightforward manner to give, KUpoT{Xk)N = -[AJN(ipk) + AJE(xk)\ where A J E ^ ) is the second term in equation C.5. Using standard integral tables, AJjs(Xjfc) = -x* - x' I 2 ln A + l I X ' - X i 1 - A ' A - l 7T - ) — — arctan A' 2 ]} x* > x' x* < x' where, A = and A' I X M x V xjb-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)*! = 7T KdN + 2KcN AX, Subcase ® : x^ = x 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 xN + e2- The limit of the resulting expression is calculated as ei and e2 tend to zero. The calculations below are done in detail to provide an example of evaluating the Cauchy principal part of an expression. From equation C.5, KUpoT(X-k)p — «i.«a-»o — 7T A X „ JO 1 - e i KfNU) JXN+e2 X f c — K7B(0 APPENDIX C. 137 Using integral tables, inserting the limits, and recalling that KbN = 0, KUpoT{xk)N = l im 1 «i - o -7T l im + «2 -»0 (2ifcN{l-e 1} + 3iaN{l-e1}2) In AX, 3/aw{f-e1}-2iircK AX, In V/ l + x ^ + 1 2TT{XN-X'} Using the binomial theorem, v / r + ^ - 1 '1 + ^ 2 X W X 1 + £ 2 2 {x„ - x'} therefore, KUPOT{X)C)N = 7T AX N + {2KcN + 3KdN) ln ln 2TT{X„-X'} 4 {x„ - x'} ^ 2 - \KdN - 2KcN - 2 From continuity of slope between A/N and AJE at x = xN, and from the scaling equation for xjj, 2KcN + 3KdN AX, - 1 ex 1 and — = 2 K-x'} Substituing into -K^or KUpoT[xk)N = -[AJE(*N)] e2 AX, where,
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A viscous-inviscid interaction procedure
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A viscous-inviscid interaction procedure Stropky, Dave 1988
pdf
Page Metadata
Item Metadata
Title | A viscous-inviscid interaction procedure |
Creator |
Stropky, Dave |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | 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, δ*(x). The formula used to update δ*(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 TEACH. 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 TEACH 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 TEACH code. |
Subject |
Fluid dynamics -- Approximation methods |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0097925 |
URI | http://hdl.handle.net/2429/28521 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1988_A7 S77.pdf [ 6.05MB ]
- Metadata
- 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
- 831-1.0097925-fulltext.txt
- Citation
- 831-1.0097925.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0097925/manifest