.. K v. tau-wiv : 33.. r t. ... 5 - {Di- . ¢n . Wu. 2 . . . .5. an JAN}... .5555.“ 16:3: .5 .. «Judah. zuofion... 35.... . . nunfiufl 5..ny . 2- lio- . . it: hm; . I. v :— Q 713...! .firuvlv-H‘ banana). .. . 2. 4...? . . “3&- ..i) .t 4 ‘I 1.11 5.. 6;! tr. :5. 3.3... 1.5.... a 5:0. .25.“:- .515: . a) K 1.16112... yvttted.’ .. "than“... .51 it’ll-"9. usual I. . £28.03“ : xi: .i...I ngzxnuxeiialtufiz .13.! .15.»..5 3:. a}? . .. 130.113.. 3. it I)... . air? 2.3.1.5Eg , x .2... J .1}. .1 AF”) l... up... 5.55...» .. “55«MW#V9531§M.W5.%§M¥%2 .2 _ 5. ..L. 2, r .5; . I, .1. . . . .t. 55.5%.... .. , waggg . . 3m... x . . THEBlS ZECQ LIBRARY Michigan State University This is to certify that the dissertation entitled WAVELETS AND THE NUMERICAL SOLUTION OF HEAT TRANSFER AND NEWTONIANINON-NEWTONIAN FLUID FLOW PROBLEMS presented by MSW has been accepted towards fulfillment of the requirements for Doctor of Philosophy . Mechanical Engineering degree m _, Dr. Andre Benard '1" M1 {Q/A- 1 m. " V Ma'pr professor Date May 31, 2000 MSU is an Affirmative Action/Equal Opportunity Institution 0-12771 PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 11100 W.“ WAVELETS AND THE NUMERICAL SOLUTION OF HEAT TRANSFER AND NEWTONIAN/NON-NEWTONIAN FLUID FLOW PROBLEMS By Ahmed S. Sowayan A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2000 ABSTRACT WAVELETS AND THE NUMERICAL SOLUTION OF HEAT TRANSFER AND NEWTONIAN/NON-NEWTONIAN FLUID FLOW PROBLEMS By Ahmed S. Sowayan Wavelet-based methods have demonstrated great potential for solving partial dif- ferential equations of various types. In this work, the capabilities of wavelet-based methods are explored by solving various heat transfer and fluid flow problems. In- herent to these techniques are difficulties encountered in implementing the boundary conditions. Two methods are therefore studied to implement Dirichlet boundary con- ditions. These methods are the Fictitious Boundary Approach and the Fictitious Domain / Penalty Formulation. These methods are evaluated by solving a number of one-and two—dimensional problems. A detailed description of the procedure used is provided for each method. It is found that both methods provide an efficient and good approximation to the solution of differential equations. The Fictitious Domain / Penalty method is, however, found to be superior because the resulting system of equations is block circulant and can be inverted easily. The Fictitious Domain / Penalty method also exhibits less errors in most of the considered examples. In the Fictitious Domain/ Penalty Formulation, the resulting large system of equation is solved iteratively via the Conjugate Gradient Method and the Preconditioned Conjugate Gradient Method. The problems modified with the Fictitious Boundary Approach are solved using simple Gaussian elimination. The Fictitious Domain/ Penalty formulation along with the Preconditioned Con- jugate Gradient Method are then used to solve heat conduction problems, as well as Newtonian and non-Newtonian fluid flow problems. The fluid flow problems in the present study are formulated in such a manner that the solution of the continuity and momentum equations are turned to solutions of a series of Poisson equations. This is achieved by using methods known as the Conjugate Gradient and Steepest Descent Methods and a segregation procedure of the dependent variables. In addi- tion, standard methods can be used to solve the nonlinear problems associated with non-Newtonian fluid. The Picard iterative method was used here. This formulation provide accurate results rapidly using personal computers, regardless of the geome- try of the problem. The fluid flow problems studied consist of the lid-driven cavity box and rotating concentric cylinders. These problems are solved for both Newto— nian and non-Newtonian fluids. It was therefore shown that wavelets can be used to solve nonlinear non-Newtonian fluid flow problems using Galerkin / F ictitious Domain formulation combined with more established techniques. I dedicate this work to my parents, my wife, and my children. iv ACKNOWLEDGMENTS It would not have been possible to complete this work without the help of many people. I am so grateful to my advisor Dr. André Bénard for introducing me to the exciting world of wavelets and PDEs. Also, I am thankful for his guidance, support, and understanding throughout my research. His stimulating comments and arguments have been a constant source of inspiration. I would like to express my appreciation and thanks to my committee members: Drs. Alejandro Diaz, Michael Frazier, Craig Somerton, and John McGrath. I also truly thank the “General Organization for Technical Education and Voca- tional Training (GOTEFT)” in The Kingdom of Saudi Arabia for supporting me and sending me to this exciting country to pursue my higher education. I also would like to thank the “Saudi Arabian Cultural Mission to the USA.” in Washington DC. for their endless help in communicating with GOTEFT. The Saudi Student House (SSH) in Lansing, Michigan has given me a wonderful time to be able to deal with the frustration encountered in the Ph.D. research. I am also grateful to the stafl of the Islamic Center of the Grater Lansing Area for providing a warm and a comfortable atmosphere for the worship services, which kept me on track with my research. I also would like to thank my familial members who are too many to list, but I shall attempt anyway. First and foremost to thank are my parents, Salih (dad) and Hussah (mom). My brothers, Mohammad, Abdallah, Abdulaziz, Suliaman, and Ali are also to be thanked. My sweet sisters, Norah and Muznah, are also thanked for their true trust in me to finish this work. I am also truly thankful for my beloved wife Galya Sulim, who tirelessly helped me, understood me, and encouraged me during the course of the last ten years. I am deeply grateful for you Galya. I also thank my sweet daughters Sarah and LeeAnn whose wonderful smiles and noise kept me going toward completion of this work. Last but not least, all thanks be to Allah for his guidance. vi TABLE OF CONTENTS LIST OF TABLES x LIST OF FIGURES xi 1 Introduction 1 1.1 Problem Statement and Objectives ..................... 1 1.2 Motivations for Using Wavelet Functions .................. 3 1.3 Wavelets and the Solution of Partial Differential Equations ........ 5 1.4 Dissertation Outline ............................. 7 2 Review of Wavelet Functions and the Solution of Differential Equa- tions 9 2.1 Introduction .................................. 9 2.2 Daubechies Wavelets ............................. 10 2.3 Application of Wavelet Functions ...................... 11 2.4 Properties of Wavelets ............................ 12 2.5 Examples of Orthogonal Wavelet Functions ................ 14 2.6 Wavelet-Based Methods for Solving Differential Equations ........ 17 2.7 The Wavelet-Galerkin Method .................. , ...... 19 2.7.1 Connection Coefficients .......................... 21 3 Solving Differential Equations Using the Wavelet Galerkin Method 24 3.1 Overview of the Wavelet-Galerkin Method ................. 24 3.2 Brief Review of the Fictitious Boundary Method in one-dimension . . . . 25 3.2.1 Example I: One-dimensional Poisson equation .............. 27 3.2.2 Simple Harmonic Oscillator ........................ 31 3.3 Application Of the Fictitious Boundary Method in Two Dimensions . . . 35 3.3.1 Two Dimensional Steady State Heat Conduction Problem ....... 40 3.4 The Penalty Formulation ........................... 46 3.5 Fictitious Domain/ Penalty Formulation ................... 47 3.6 Solution of problems with the Fictitious Domain Penalty method . . . . 51 3.6.1 One-Dimensional Poisson Equation .................... 51 3.6.2 Simple Harmonic Motion Equation .................... 56 3.6.3 Steady State Heat Conduction Equation ................. 57 3.7 Summary ................................... 65 vii 4 Solving Newtonian Fluid Flow Problems Using a Wavelet Galerkin Method 67 4.1 Introduction .................................. 67 4.2 Governing Equations ............................. 68 4.3 Fictitious Domain Method for Fluid Flow Problems ............ 71 4.4 Penalty Formulation of the Governing Equations ............. 72 4.5 Wavelet-Galerkin Method and Discretization of the Governing Equations 76 4.6 Solution of the Resulting System of Equations ............... 79 4.7 Solution Methods ............................... 80 4.7.1 A Conjugate Gradient Method (CGM) .................. 80 4.7.2 The Steepest Descent Method (SDM) ................... 83 4.8 Examples ................................... 84 4.8.1 Lid-Driven Cavity Problem ........................ 85 4.8.2 Flow between concentric cylinders .................... 87 4.9 Discussion and Summary ........................... 96 5 Solving Non-Newtonian Fluid Flow Problems Using Wavelet Galerkin Method 97 5.1 Introduction .................................. 97 5.2 Governing Equations ............................. 99 5.2.1 Constitutive equations ........................... 100 5.2.2 Governing equations in scalar form .................... 101 5.3 Fictitious Domain/ Penalty Formulation of the Governing Equations . . . 102 5.4 Wavelet-Galerkin Method and Discretization of the Momentum Equations 107 5.5 Solution of the Resulting System of Equations ............... 110 5.5.1 Picard Iteration Method .......................... 111 5.5.2 Calculation of the Viscosity Function 17’" ................. 112 5.6 Examples ................................... 114 5.6.1 Flow Between Rotating Concentric Cylinders .............. 115 5.6.2 Lid-Driven Cavity Problem ........................ 118 5.7 Discussion and Summary ........................... 119 6 Conclusions 122 6.1 Summary and Contributions ......................... 122 6.2 Findings of the Study ............................ 123 6.3 Future Research Topics ............................ 124 APPENDICES 126 A Construction of Generalized Wavelet Systems 127 A.1 Construction of The Filter Coeflicients ................... 128 A2 Construction of the Scaling Function .................... 131 A.3 Wavelet Multiscale Representation of Functions .............. 133 B Incorporation of the boundary measure function Haw“ 139 viii C Algorithms 142 BIBLIOGRAPHY 144 ix LIST or TABLES 3.1 Performance of PCG method and CG method ............... 65 1.1 1.2 2.1 2.2 2.3 2.4 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 LIST on FIGURES Schematic representation of the theoretical and computational progress that lead to solution of viscoelastic flow problems ............ Schematic of the reasoning behind the use of wavelet functions to solve PDE’S .................................... Daubechies D6 scaling (left) and wavelet functions (right). ........ Illustration the numerous areas where wavelets have been applied ..... Meyer scaling and wavelet functions (a,b) and Coiflet scaling and wavelet functions (c,d). .............................. D8 scaling and wavelet functions (a,b) and D12 scaling and wavelet func- tions (c,d) .................................. Solution of the one-dimensional Poisson equation using the Fictitious Boundary Method at two resolution levels. ............... Error (u — uexad) in the wavelet solution for difl'erent levels of dilation using the Fictitious Boundary Method. ................. Solution of the differential equation describing a simple harmonic oscillator using the Fictitious Boundary Method. ................. Error (:1: — mama) in the wavelet solution for different levels of dilation using the Fictitious Boundary Method. ................. In the fictitious boundary method, a domain of arbitrary shape w is en- larged in all directions to form a new domain {2 of simple geometry suitable for computations with an insulated boundaries ......... Geometry and boundary conditions used for solving the two-dimensional heat conduction example. ........................ Nonzero entries of matrix A (Level=5) .................... Temperature contours at different levels of dilation using the Fictitious Boundary Method. ............................ Error (u—uemct) contours for different levels of dilation using the Fictitious Boundary Method. ............................ 3.10 In the Fictitious Domain method, a domain w of arbitrary geometry is embedded into an auxiliary, simple domain (2 with periodic boundary conditions .................................. 3.11 Nonzero entries of matrix A (Level=4) .................... 3.12 Fictitious domain/ penalty solution of the 1D Poisson Equation. ..... 3.13 Error (u — new“) in the wavelet solution of the one—dimensional Poisson equation at different levels of dilation ................... xi 11 12 15 16 32 33 36 37 39 41 42 44 45 48 53 54 55 3.14 3.15 3.16 3.17 3.18 3.19 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 5.1 5.2 5.3 5.4 5.5 Solution to the Simple Harmonic motion equation using the Fictitious Domain / Penalty Method shown at various levels. ........... Error (:1: — rem“) in the Wavelet Galerkin Fictitious Domain/ Penalty Method solution of the Simple Harmonic Oscillator equation ...... Illustration of the nonzero entries of matrix A (Level=4) .......... Temperature contours at different levels of dilation using the Fictitious Domain/ Penalty Method. ........................ Error (u—uexact) contours for diflerent levels of dilation using the Fictitious Domain/Penalty Method (5 = 10‘“). .................. Convergence of the solution using the Conjugate Gradient method. In the Fictitious Domain method, a domain w of arbitrary geometry is embedded into an auxiliary, simple domain (2 with periodic boundary conditions .................................. Domain and boundary conditions for the lid-driven cavity problem. . . . Comparison of the velocity vectors between SDM and CGM for the driven cavity problem (level=4) .......................... Position of the zero velocity point in both SDM and CGM (level-=6). . . Comparison in the residual history of both SDM and CGM for level=5. . Iterations required for the PCG solver along the SDM iteration for the u component of the velocity. ........................ Iterations required for the PCG solver along the SDM iteration for the 1) component of the velocity. ........................ Streamlines obtained by the SDM for the driven cavity problem (level=6,7). Residual history of the SDM at various levels for the lid-driven box. . . . Decomposition of the tangential velocity into Cartesian coordinate system. Velocity vectors obtained by the SDM for the flow between two concentric cylinder problem (level=5,6) ........................ Streamlines obtained by the SDM for the flow between two concentric cylinder problem (level=5,6) ........................ Streamlines obtained by the SDM for the flow between two concentric cylinder problem (level=7,8) ........................ Comparison of the velocity magnitudes of the wavelet and exact solutions at the crossection y = 0, (0 = 0, 7r in polar sense). ........... Comparison of the velocity magnitudes of the wavelet and exact solutions at the crossection y = O, (0 = 0, 7r in polar sense). ........... Velocity vectors for the flow between two concentric cylinder problem (n = 0.8). .................................... Comparison between the wavelet and exact solutions at the crossection y = 0, (0 = 0, 7r in polar sense)(n = 0.8). ................ Velocity vectors for the lid-driven cavity problem (level=5). ....... Streamlines for the lid-driven cavity problem (level=6), left Newtonian fluid, and right non-Newtonian fluid. .................. Residual history of the Picard iteration method ............... xii 57 58 60 62 63 71 86 86 87 88 89 90 91 91 92 9'3 93 94 94 95 116 A.1 Haar Scaling and Wavelet Functions .................... 136 A2 Multilevel representation of the function f (1:) = sin(27r:z:) ......... 137 xiii Chapter 1 Introduction 1.1 Problem Statement and Objectives The objectiveof this study is to develop and adapt a methodology based on wavelets to solve non-Newtonian fluid flow problems. Wavelet functions are considered be- cause of their potential ability to capture solutions with sharp gradients such as those encountered in non-Newtonian viscoelastic fluid flows. It is also demonstrated that wavelet-based methods have the capabilities of solving large problems, such as en- countered in non-Newtonian flow problems, with modest resources. The amount of work required to compute a solution for the most elementary vis- coelastic flow is considerable. This is due to the large number of dependent variables, coupled governing equations, and the presence of steep stress gradients. For example, a two-dimensional problem require six variables: two velocity components, pressure, and three stress components, and in some formulations three additional rate of strain components. Viscoelastic flow Problems? Non-Newtonian Fluid Flow Problems Newtonian Fluid Flow Problems Heat Conduction Problems Fictitious Boundary Method Fictitious Domain Method Penalty Solution Wavelet Analysis & Theoretical Background Figure 1.1: Schematic representation of the theoretical and computational progress that lead to solution of viscoelastic flow problems. A first step toward solving viscoelastic fluid flow problems, therefore, consists in developing a solution procedure for generalized non-Newtonian fluid flow problems within the context of wavelets and multiresolution analysis. This will further one’s understanding of the capabilities of wavelet-based methods and multilevel techniques. The steps required to achieve this goal are shown in Figure (1.1), where the various building blocks required are put in a pyramided form, with the most important ones at the bottom. The following overall steps are required to solve the problems of interest in this work: 0 To study, understand, and develop wavelet-based methods for solving partial diflerential equations in the context of fixed scale analysis 0 To study various methods for the imposition of boundary conditions (espe- cially the Dirichlet boundary conditions) in the context of the wavelet Galerkin method 0 To apply a wavelet-based method to common engineering problems such as heat conduction problems, Newtonian, and non-Newtonian fluid flow problems (nonlinear) The very recent discovery of wavelet functions, and their use for solving nonlinear engineering problems makes this study very novel. It should also be noted that the techniques presented in this dissertation have the following advantages: There is no need to generate a complex computational grid in the case of irreg- ular domain The numerical implementation of the methodology is relatively simple, and independent of the geometry of the problem These techniques can be applied to any number of dimensions 0 Their extension to nonlinear problems is straightforward, as shown below 1.2 Motivations for Using Wavelet Functions Problems with sharp gradients are typically difficult to analyze, and a fine grid is required to capture the gradient. The discretization of these problems leads to a very large number of degrees of freedom. Classical methods such as finite element or finite difference method sufler from large memory requirements, especially when three-dimensional problems are considered. In view of this large memory requirement, new techniques are needed for tackling such large problems. Wavelet-based methods can provide an alternative since they were recently demonstrated to provide low in-core memory, bounded condition number of the resulting system of equations, efficient preconditioning, and to provide an accurate solution. For example, ordinary differential equations problems, when solved with a wavelet—Galerkin technique, give rise to a circulant matrix [1, 2, 22], e.g. a1 (12 03 a4 a5 a5 0.1 (12 a3 a4 0 = (14 a5 a1 02 a3 (1'1) (13 a4 a5 (11 (12 02 a3 a4 a5 01 ' 'NxN which requires only N entries for storage rather than N xN and which can be inverted easily and rapidly with techniques such as a Fast Fourier Transform or a Fast Wavelet Transform. Similar benefits appear when wavelets are used for solution partial dif- ferential equations. The reasoning behind the proposed approach is depicted in the flow chart shown in Figure (1.2). This flow chart can be interpreted as follows: the current formulation significantly reduces the computational process required and allows to solve relatively large problems using desktop computer capabilities. Solving large problems requires, however, the use of an iterative solver, e.g. Conjugate Gradient methods. Also, using Start Mesh Refinement or Small Grid Large Problems With Need Solution Via Iterative Such As Conjugate Gradient Sharp Gradeints Methods Methods ACMIC and Fm F351 8‘ EffiClCfll Solution Solution Need Come Periodic Functions ‘7'— Good Preconditioners om Such As Wavelet Functions Figure 1.2: Schematic of the reasoning behind the use of wavelet functions to solve PDE’S. periodized wavelet functions results in the efficient and rapid solution of the resulting system of equations. The efficiency and rapidity found in solving this system of equations come from the well-behaved stiffness matrix, i.e. it has a bounded condition number, and it also allows good preconditioning of the resulting system of equations. 1.3 Wavelets and the Solution of Partial Differen- tial Equations Practical uses of wavelets were initially developed in the field of signal and image processing. It was soon realized also that ordinary diflerential equations could be solved using wavelets. The application of wavelet-based methods to the numerical solution of partial differential equations (PDE’s) has now been studied both from the theoretical and computational point views in a number of publications. The literature on wavelet-based methods applied on the solution of PDE’s can be divided into four major trends [16]: 1. Adaptive techniques in which wavelets are used to detect where the grid has to be refined or coarsened to optimally represent the solutions. Two of these methods are the o Wavelet Element Method (WEM) [39, 40] o Wavelet Finite Difference Optimized Method (WFDOM) [23] 2. Solutions to PDEs are also achieved via the construction of some basis function by shift-invariant refinable spaces, which is often referred to as multiresolution analysis [8, 25]. 3. Wavelet collocation methods which are based on the use of the autocorrelation function of some orthogonal compactly supported wavelet functions. In this trend the compression properties of wavelet are used [5]. 4. Wavelet-Galerkin methods, in which scaling function bases may be used as alternative bases in Galerkin methods. In these methods different versions of the fictitious boundary/domain techniques are presented [2, 15, 29, 43]. This last category has exhibited potential for solving large problems very effi- ciently. The work presented in this dissertation is based on the fourth trend, which is the wavelet-Galerkin Method. It is worthwhile mentioning that the implementation of the boundary conditions of solving PDE’S in this trend can be subdivided into: 0 Appending boundary conditions via the so-called Capacitance Matrix Method [1,2,31,42] o Appending boundary conditions using Fictitious Domain/ Lagrange Multipliers [15,16,25] 0 Implementing boundary conditions via Fictitious Domain / penalty method [22, 20,37,43] A wavelet Galerkin formulation for solving non-Newtonian fluid flow problems has not been presented to date in the literature. However, wavelets have been used to solve various partial differential equations. It is therefore worthwhile to review in the next chapter wavelets and some of the related papers to this study as they provide the background for this work. 1.4 Dissertation Outline This dissertation is divided into six chapters and two appendices. In chapter two, the wavelet functions in terms of their construction, properties, and use are presented briefly along with their use to solve differential equations. In chapters three and four, various techniques and implementations of the Dirichlet boundary conditions in the context of wavelet-Galerkin method are presented. These techniques are used to solve engineering problems such as heat conduction problems, Newtonian and non- Newtonian fluid flow problems in chapters five and six. Conclusions of the study are provided in chapter seven. Chapter 2 Review of Wavelet Functions and the Solution of Differential Equations 2. 1 Introduction Wavelet functions have generated significant interest from both theoretical and ap- plied research over the last few years. The name wavelet comes from the requirement that they should integrate to zero, waving above and below the x-axis. A number of ideas in the field of wavelet analysis come from work done in subband coding, group theory, and coherent states in physics [44]. The unifying concepts for understanding wavelets were provided recently by Meyer, Mallat, Daubechies, Battle, and many others. Since then, the number of applications where wavelets have been used has exploded. Many different types of wavelet functions have been presented over the past few years. In this research, the Daubechies family of wavelets will be considered due to their useful properties. These properties are presented in the next section. 2.2 Daubechies Wavelets Daubechies wavelets are compactly supported functions. This means that they have non zero values within a finite interval and have a zero value everywhere else. This makes them quite useful for representing solutions of partial differential equations. Their orthogonality, compact support, and their capability of providing an.exact representation of polynomials of a fixed degree allow the proper calculation of regions of strong or sharp gradient. The Daubechies scaling function is given by the following relation: (15(9)) = ak¢(2x — k) (2.1) where N is an even positive integer number that determine the number of non-zero co- efficients, k is an integer that translates the scaling function ¢(a:), and the coefficients ak are called the filter coeflicients. A wavelet function is defined by N-l m) = Z(—1)"b~-1_k¢(2x — k) (2.2) k=0 The relation between the coefficients ak and bk is given by 10 W) WI) Figure 2.1: Daubechies D6 scaling (left) and wavelet functions (right). bk = (—1)k+1aN_1_k (2.3) Additional information on the construction of these type of wavelets is presented in Appendix (A). Figure (2.1) shows the scaling function and the associated Daubechies wavelet of six filter coefficients D6, which has support on [0,5], i.e. the genus for this wavelet is 6. Within the Daubechies family of wavelets, there are subclasses of wavelets distinguished by the number of filter coefficients. 2.3 Application of Wavelet Functions Wavelet functions are used in a variety of areas in science and technology. Figure (2.2) indicates the areas where wavelets have been used. In the context of this dissertation, wavelet functions are used to represent the solution of partial differential equations. Wavelet functions have a very important property, which is the localization in 11 Time frequency analysis Others Surfaces Differential equation / Multign'd Image compression Subband filtering Coherent states Splines Transient Analsis Multiresolution Analysis Adaptive gn'dding Figure 2.2: Illustration the numerous areas where wavelets have been applied. both frequency (or scale) via dilation and in time (or space) via translation. This localization is often advantageous in many cases. Like Sines and cosines in Fourier analysis, wavelets are used as basis functions in representing other functions. The representation of functions with peaks can be done more efliciently with a localized function having a compact support. Many classes of functions can then be represented by wavelets in a more compact way. For example, functions with discontinuities and functions with sharp spikes usually take substantially fewer wavelet basis functions than sine-cosine functions to achieve a comparable approximation. 2.4 PrOperties of Wavelets As seen above, there are many properties of wavelet functions that make them unique. These properties are presented below and commented upon briefly. 12 o Orthogonality The wavelet transform in orthogonal wavelets is a unitary transformation (i.e. its adjoint is its inverse). Consequently its condition number is equal to 1. Also the orthogonality property links the L2 norm of a function to the norm of its wavelet coeflicients. 0 Compact support This implies that they have non zero values within a finite interval and zero everywhere else. 0 Rational coefficients For computer implementations, it is of use if the filter coefficients ak and bk are rational or dyadic rational. On a computer, multiplication by a power of 2 corresponds to shifting bits, which is a very fast operation. 0 Symmetry This prOperty is important in signal processing applications and the solution of differential equations. 0 Smoothness The smoothness of wavelets plays an important role in compression applica- tion. Smooth basis function are desired in numerical analysis applications where derivatives are involved and often must be continuous up to a given order. 0 Number of vanishing moments of the wavelet functions 13 This property determines the convergence rate of wavelet approximations of smooth function. 0 Analytical expressions Explicit analytical expression for wavelet or scaling functions are not always available. 2.5 Examples of Orthogonal Wavelet Functions The most simple examples of orthogonal functions are the Haar scaling and wavelet functions shown in Figure (A2). A more interesting example is the Meyer wavelet and scaling function shown in Figure (2.3). These functions belong to C°° and have faster than polynomial decay. These wavelets are not compactly supported, but their effective support is [-8,8], i.e. their Fourier transform is compactly supported. Their effective support means that they do not have a zero value outside the interval [-8,8] but have a very small value. The scaling function and wavelet are symmetric around 0 and -1/2 respectively, and the wavelet has an infinite number of vanishing moments. Coiflets wavelet shown in Figure (2.3) are compactly supported with order N (N = 1, 2, .., 5). They are nearly symmetric with compact support of (6N — 1). The most frequently used wavelets are the original Daubechies wavelets. They are a family of orthogonal wavelets indexed by N E Z, where N is the number of vanishing wavelet moments. They are supported in an interval of length (2N — 1). A disadvantage of i this class of wavelets is that they are not symmetric or antisymmetric except when N = 1. Figures (2.1, 2.4) show different Daubechies wavelet and scaling functions. 14 (a) 1 .5 f (b) 1.5 : NX) —o.s 4 _, . -10 -5 o 5 10 -1o —5 o 5 10 x x (d) 15 , .......................................... 1’ W ............................... o.5»--~ 3? ‘5 ....................... W) ............................... ............................. ...... 15 20 10 15 Figure 2.3: Meyer scaling and wavelet functions (a,b) and Coiflet scaling and wavelet functions (c,d). 15 (a) 1.5 . W) W) (d) 1.5 3 : .' W) W) -1.5 1 'L o Figure 2.4: D8 scaling and wavelet functions (a,b) and D12 scaling and wavelet functions (c,d). 16 2.6 Wavelet-Based Methods for Solving Differen- tial Equations Among the first researchers to use wavelets to solve differential equations were Amaratunga et al. [1, 2], who solved the one-dimensional Helmholtz equation ( uu(x) + au(T) :2 f (x) ), and the two-dimensional Poisson equation ( outlay) + uyy(:z:, y) = f (:13, y) ) using the wavelet Galerkin method. Daubechies compactly sup- ported wavelets were used in these papers, and this allowed refining the solution in regions of high gradient, e.g. stress concentrations without having to regenerate the mesh for the entire problem. The solution to the one-dimensional Helmholtz equation, via the wavelet Galerkin method with periodic and non-periodic boundary conditions, was compared to a sim- ple finite difference solution. It was found in reference [2] that the rate of convergence of the wavelet solutions compared extremely favorably to the finite difference solu- tions. Qain et al. [31] presented a numerical method for the solution of partial dif- ferential equations using the so-called capacitance matrix method. It was found in reference [31] that the numerical solutions has a convergence related to the genus of the Daubechies wavelet basis function. The rate of convergence was also found to be independent of the geometry. The method in reference [31] was tested for two- problems. First, the unsteady, two-dimensional Navier-Stokes equations which were written in the stream function vorticity formulation. In this case a rectangular do- main was chosen. The second problem was the periodic Helmholtz equation (2.4), 17 which has the form (A + a)U = F (2.4) where a is constant and A is the Laplacian differential operator. In this case an L-shape domain was used to illustrate the idea of using the method of wavelet capac- itance matrix method. Weiss [42] solved numerically the flow of an inviscid, incompressible, two- dimensional fluid based on compactly supported wavelets (wavelet Galerkin method). Weiss applied this method to study the long-time evolution of non periodic Navier Stokes flow in an L-shape region. Restrepo et al. [34] solved the wave equation by the wavelet Galerkin approach. A comparison of the wavelet Galerkin method to the standard finite difl'erence and the Fourier pseudo-spectral methods was made. The comparison was based on the computational efficiency, which is the reciprocal product of the wall clock time and the storage requirement. It was found that the wavelet Galerkin method is comparable to the finite difference method, requiring less storage but more time than the finite difference. The Fourier-pseudo spectral method is found to be the most efficient method for solving the wave equation. Glowinski et a1. [22, 20, 43] presented a new formulation for the solution of Dirich- let boundary value problems. This formulation is called the fictitious domain / penalty method. This method eliminates the difficulties associated with generating a com- plex mesh in the case of an irregular domain. This is based on the fact that one 18 can expand the boundary measure under the chosen basis which lead to a fast, ap- proximate calculation of the boundary integral. Also in the papers by Glowinski et al.[22, 20], a multigrid preconditioner is presented for the conjugate gradient method, which provides an efficient solver for the linear system arising from a wavelet-Galerkin discretization of a Dirichlet boundary value problem. The preconditioner was chosen to be a wavelet-based multigrid method for solving an elliptic equation, however over the fictitious domain and with periodic boundary conditions. The Neumann problem is solved in reference [20] with a slightly diflerent formulation than with the Dirichlet problem. A formulation using Lagrange multipliers to enforce the boundary conditions was presented in [15, 16] to solve large topology optimization problems. This study was tested by solving problems with Simple geometry. Approximate solutions for the resulting system of equations in [15] were achieved via a preconditioned conjugate gradient method. It was shown in [15, 16] that the rate of convergence of the con- jugate gradient method is independent of the problem size by using an appropriate preconditioning matrix. 2.7 The Wavelet-Galerkin Method As discussed above, a number of methodologies are available for solving differential equations. The wavelet-Galerkin method is briefly discussed here since it will be used later. 19 For a partial differential equation of the form G(V,VI, V1,“. . .) = 0 (2.5) An approximate solution is defined as t V = Z G. ¢(:z: — k) (2.6) kz—s where s and t are integers that are the translates of d>(x — k). The solution in this case is projected onto the subspace spanned by (s,t) = {¢(x — k) : k = —s, . . .,t} (2.7) The coefficients Ck in the expansion equation (2.6) are found by substituting the expression l7 into the partial differential (2.5) and again projecting the resulting expression into the subspace @(s, t). The coefficients Ck, are found by evaluating the following relation: /°° an: — k)G(V, v,, v“, . . . ) d2: = o (2.8) In this study, we used a class of compactly supported scaling functions that are introduced by Daubechies [13]. The field variables are projected into the subspace of the trial function. When a test function from the same space is used, a system of 20 algebraic equations for the coefficients of the field variable results. This is achieved by evaluating the integral (2.8), and making use of the orthogonality property for the scaling function. To evaluate expression (2.8), we must evaluate integrals of the form /¢(:z:) ¢x(z — k1) . . . ¢(:c —- k2) . . .dx (2.9) and the results of these computations (2.9) are called the connection coefficients. 2.7.1 Connection Coefficients Using the wavelet Galerkin method to solve the governing equation requires the com- putation of integrals of the forms 03:2 = / was) «5339) dx (2.10) 033.13.“: = / am 7'29) if: dz (2.11) where the superscript d,- refers to the number of differentiations with respect to a: of the scaling function q5(:c). These formulas are called two-term and three-term connection coefficients. The name of these inner products, connection coefficients, was coined by Latto et al. [26]. Latto et al.[26] developed an exact computational method for evaluating these inner products [35]. This method replaces a quadrature problem by a linear algebra problem, namely solving a simultaneous linear system of 21 equations instead of calculating integrals. Following the convention in [26], one refers to the inner products as connection coefficients as 02,2): ($2.232) =1: 22(2: (2.12) C°3°.=(¢.,¢.¢>i) =/_: «52(2: )aw )¢.(z)dx (2.13) C2301: ((1524),,- ¢z) =/: ¢k($ 9W )¢z($) d1: (2.14) It should be noted that the translates of (9(2) is conventionally defined as 2(1) = (WE " 3) Using integration by parts along with some change of variables, the following rela- tionships for the connection coeflicients are found: 0222 =22. (2.15) and 03:12 ___ (_1)d1 Com,d2+d1 (216) 22 where the superscript d,- refers to the number of differentiations with respect to a: of the scaling function ¢(:r), m and l are translation integers. The connection coefficients depends on the choice of the basis function. They can be calculated prior to the solution and then used as input data. Software and algorithms are available in the literature for evaluating the connection coefficients [10, 26, 35]. 23 Chapter 3 Solving Differential Equations Using the Wavelet Galerkin Method 3.1 Overview of the Wavelet-Galerkin Method The use of wavelets as bases in a Galerkin-type method to solve differential equations requires a computational domain of simple shape e.g. a square in two-dimensions. The two methods discussed below use the fairly old idea of extending the (typically complicated) domain into a larger auxiliary domain of simple shape. Both methods differ in the implementation of the Dirichlet boundary conditions and the methodol- ogy used for extending the domain. The first method is called the “Fictitious Boundary Approach” and was intro- duced in [27, 28], while the second method is named the “Fictitious Domain Penalty 24 Formulation” and was developed in [22, 20, 43]. In this chapter both techniques are briefly presented followed by solved examples. 3.2 Brief Review of the Fictitious Boundary Method in one-dimension A differential equation of the form L(u(a:)) = 0, in w, x 6 [0,8] 3 E Z (3.1) subject to u(0) = a and u(s) = fl is considered, where L is a second order linear differential operator. The so-called “Fictitious Boundary” method consists of extending the domain (1) of the solution u(x) ([0, 3]) to Q with an interval wider than u), for example [—t, s + t]. The weak form of the new problem is then /L(ua($))v(:r)d:c = 0 in Q where v E V (3.2) where V is H 1 for example. It is subjected to the boundary conditions 32a 02a 62: lx=—t= “51‘:- lx=s+t= 0 SJ 6 Z (3.3) 25 where ua(:r) is the approximate solution to (3.1). The left end point where a: = —t, corresponds to the first point on the first test function, conversely :1: = s+t corresponds to the last point on the last test function. A wavelet expansion for u is defined as 21 . . ua($) = 201222232: — k), k, j e z (3.4) k=l where ()6 is the Daubechies scaling function [13] with six filter coeflicients (N = 6), and the unknown dk’s are the wavelet expansion coefficients. The domain to of the solution of u(:z:) is [0, s] and is widened to a new domain (2. The test functions 21’s are chosen to be 22(2) 2 (M232 — l), where l is an integer called the translate of ’U(.’L‘). The solution to (3.1) in Q has assumed fictitious boundary conditions of the form (3.3). The boundary conditions given by equations (3.3) should have a minimal impact on the solution within the domain [0, s] for this methodology to be accurate. Substituting the values of u and v in (3.2) results in simultaneous algebraic equa- tions in dk’s. The solution of ii in the original domain to is achieved by imposing the original boundary conditions in (3.1), i.e. appending them to the resulting system of equations meal.-.) = a'=Zd.2%¢(2ij.=o—k) (3.5) I: 21.22:.) = a = 2 (1.2222122, — k) (3.6) I: It should also be noted that in evaluating equation (3.2), certain integrals must 26 be evaluated. These integrals are referenced to in the literature as the connection coefficients C21,”? [26, 35] and defined as Cm = [(221 (21¢ - k)¢“2(2j2: - l)d2: (3.7) where d1 and d2 refers to the number of differentiation of the scaling function 45(3) with respect to x. Equation (3.7) is called the two-term connection coefficient. The connection coefficients depends on the choice of the basis function and not on the problem. As disscused in the previous chapter, they can be calculated prior to the solution and then used as an input data. 3.2.1 Example I: One-dimensional Poisson equation To illustrate the “Fictitious Boundary Approach” method, the following simple dif- ferential equation is considered: u”(a:) = f(:z:) :1: 6 [0,1] (3.8) subject to u(0) = a and u(1) = 6. The weak form of equation (3.8) in Q is /_ u’(:c)v'(:c)d:z: = — f(a:)v(x)d2: \7’ v(:z:) E H1 (3.9) t -—t where 21(2) is the test function. The boundary conditions in the new domain (2, which is a: E [—t, 1 + t], are given in equation (3.3). 27 The expansion of u and f (x) are chosen to be 21' . = 2 (1.222213; — k) (3.10) k=1 2' _ = 2: 2222212: — k) (3.11) 1:21 where ek=/f() (2jx—k)dx which results from the orthogonality condition of (M232: — k) [13]. The test functions v’s are chosen to be 'u(:r) = d)(2jrr - l). Substitution of the wavelet expansions (3.10) and (3.11) into (3.9) yields 21' Z (1.03,): = — Z 202;," (3.12) i l where C? kand C10,? are defined as Czk =/¢( (2jx—i)¢ (2’x—k)d$ (3.13) Clok 0=/¢( 2jx—l)¢ ¢(2jx—k)dx (3.14) 28 Due to the orthogonality of the scaling function and its translates, equation (3.14) gives . . 1 ifl = k 03;? = / N211: — z)¢(2’:v - k)dx = 51k = 0 otherwise Also a change of variables C3,}: = Cfifi is used to simplify equation (3.12) into Edict—Ii = 6]; (3.15) where k can vary from 1 to 21. The system of equations given in (3.15) for different values of k is equivalent to the following system: Anannxl = bnxl (3.16) where n = 2j. The entries of A, X ,and b are AiJc = 02:1; bk = 9k Xi : d1 The matrix A is sparse and symmetric. To obtain the correct solution in the domain [0,1], the boundary condition equations must be included in (3.16). One approach is to force the solution to satisfy ua(:l:|I=0) = a = ng’éaamxzo —— 2') and u.(x|.=1) = a = Z d.2%'¢(2"xl.=1 - 2'). The system (3.16) then becomes Amannxl = bmxl , m = n + 2 (3.17) and Amxn is now a rectangular matrix, i.e. an overdetermined system of equations. A possible approach to solve the rectangular system (3.17) is with a least square method [36]. Another approach is to replace the equations from the system (3.16) that are related to the boundary points, for example, equation (3.5), by the boundary conditions equations. The first approach is used below. In the least square method, the system given by (3.17) is multiplied by the trans- pose of A. This yields a new system of equations that can be solved using Gaussian elimination, i.e. Anannxl : bnxl 30 where Anxn : AnxmAmxn and bnxl = Anmemxl The solution to problem (3.8) is shown in Figure (3.1) at different values of j, the dilation level factor. The following parameters were used in this example oz 2 0, fl = 1, and f (:13) = 1. The error in the wavelet solution is shown also for different levels of dilation in Figure (3.2). It can be noted from Figure (3.2) that the solution error diminishes as the dilation level j increases. The solution obtained with a level equal to 6 exhibits a very small 81' TOT. 3.2.2 Simple Harmonic Oscillator A simple harmonic oscillator can be described by the following differential equation: if(t)+w,2,:r(t) = 0 tE[0,1] (3.18) with 33(0) = 0.117(1) = 5 31 — wavelet solution exact solution Level=5 Level=6 1 - f . - 1.4 . - . - , 1.2 » 0.8 > 4 . 1 _ x X ‘5 ‘5 0.4 ~ 0.6 ~ 4 0.4 > < 0.2 > 0.2 ~ 0 * 0 L A 0 0 2 0 4 0 6 0 8 1 0 0 2 0 4 0 6 0 8 1 x x Figure 3.1: Solution of the one-dimensional Poisson equation using the Fictitious Boundary Method at two resolution levels. 32 0.05 r r 1 I I 1 T r l ' - I ' , - -. ........ Level=5 0.045- __ Lake - 0.04- _ 0.035 — a 0.03 — - 3 0.025 - J IJJ 0.02 - — 0.015 - ~ 0.01 - ~ 0.005 \/ _ 0 1 i i 1 1 l 4 L I' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.2: Error (u — uezact) in the wavelet solution for different levels of dilation using the Fictitious Boundary Method. 33 where can is the natural frequency of the system motion. Following the procedure outlined above, one can obtain the discretized form of equation (3.18) in wavelet space as Zd-C‘" +de —0 319 . 1 k—i a? 10“ (° ) where k is chosen to vary from 1 to 2i. The system of equations in (3.19) for different values of k is equivalent to the following system: Anannxl : 0 (320) where n is equal to 2i. The entries of A and X are 1,1 A1,): = Ck—i To obtain the solution in the domain t E [0, 1], the boundary conditions in equation (3.18) are included in system (3.20) by forcing the discretized solution to satisfy xa(t|¢=0) = a = Zd,2%¢(2it|.=o — 1) (3.21) 34 20(1),:1) = 5 = Z d,~2%q>(21'1|,_._I — 2‘) (3.22) The system of equation given by(3.20) then becomes Amannxl : bmxl (323) where m = n + 2, and A is rectangular. The system can be converted to a square matrix so that it can be solved. To do this, Lu et a1. [27, 28, 29] suggested that some of these equations in (3.23) could be eliminated since they are not required in order to obtain a solution in the domain of interest. Figure (3.3) shows the solution to this problem (3.18) at different values of j the dilation level factor. This example is tested for the following parameters: a = 1, fl = O,and 02,, = 7.57r. The eliminated equations are associated with k = 1 and k = 21', i.e. the equations at the edge of the extended domain. The error in the wavelet solution is shown in Figure (3.4). 3.3 Application Of the Fictitious Boundary Method in Two Dimensions The formulation of this method in two-dimensional problems is similar to the ap- proach used for the one-dimensional problems discussed previously. For a differential 35 Level=6 Level=5 x(t) x(t) 0 0.5 1 0 0.5 1 t 1 Level=7 Level=8 1.5 fl 1 .5 - x(t) X(t) -1.5 ‘ -1 .5 0 0.5 1 0 0.5 1 Figure 3.3: Solution of the differential equation describing a simple harmonic oscillator using the Fictitious Boundary Method. 36 ‘5 I I I I l I I I l ...... bvebs ° 19119136 1 h - - - level=7 ‘ - - lovel=8 0.5 " 5: ‘1. '1 '1. - I ’_\ .5", _'_' I I}: “x ’ 0“ II; ’ ‘\\*~. ,. ’7’, \‘.\ ’71, *§*. I’l 2.; : . . . - - . r m .2 . . . a o -0.5 ~ -1 '- 4 -1.5 - ~ _2 l 1 J l l 1 1 L 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.4: Error (:2: — mama) in the wavelet solution for different levels of dilation using the Fictitious Boundary Method. 37 equation of the form L(u(a:,y)) :2 0 inw (3.24) andu = g in 80) where L is a linear second order differential operator. The domain 02 is a domain of arbitrary shape that is enlarged into 0, a new domain of simple shape e.g. a square. This procedure is illustrated in Figure (3.5). In the following example, 02 is a domain that fits in a square of dimensions [0, 3] x [0, s] that is enlarged to 9 which is a square with boundaries defined by [-—t, s + t] x [—t, s + t] where s and t are integers. The problem given in (3.24) can be written in its weak form and consists in finding 11 E V(Q) such that [2L(u(a:, y))v(z,y)dxdy = O Vv(:c, y) E V(§2) (3.25) subject to the boundary conditions 611(2', y) an Ian: 0 (3.26) where V is H1 for example, and 21(23, y) is the test function which was chosen to be WI, 11) = ¢(2j$ - p)¢(2jy — q). 38 Figure 3.5: In the fictitious boundary method, a domain of arbitrary shape (.0 is enlarged in all directions to form a new domain 0 of simple geometry suitable for computations with an insulated boundaries. 39 A wavelet expansion for u(:r, y) can be defined as 2.1 00(2, y) = Z d,,k2j¢(2ix — i)(b(2jy — k) 2:, y e [0, s]2 (3.27) i,k=l Substituting the expansion of it given by equation (3.27) and the test functions v’s into (3.25) yields a discretized equation in wavelet space in terms of the connection coefficients Cfl’d’. The solution in the original domain 02 is achieved by imposing the original boundary condition in wavelet expansion, i.e. forcing the solution to satisfy ua(xl6w1y|6w) = Z di,k2j¢(2jxlaw — Z)¢(2jy|0w _ k) (328) 1']: 3.3.1 Two Dimensional Steady State Heat Conduction Prob- lem In this example, the steady state heat conduction equation with constant properties is solved in the square domain shown in Figure (3.6). The problem is stated as , V2u(:c, y) : 0 in 0) subject to u(x, 1) = 1 ) 1' (3.29) 11(1), 0) = 0 u(0, y) = O HO. 9) = 0 J 40 u(:z:,l)=1 (1’1) ”(01y)=0 u(11y)=0 V’u(I. y) = f(Ly) ‘ '7 u(:r,0) = 0 I Figure 3.6: Geometry and boundary conditions used for solving the two-dimensional heat conduction example. It should be pointed out that u(:r, 1) = 1 is enforced on the top corners shown in Figure (3.6). The test function is chosen to be v(x, y) = (15(230: — p)¢(2jy — q). Substituting the expansion of it given by equation (3.27) and 1) into the integral equation (3.25) yields the following discretized heat equation: 21' 21' Z dk,,,C;:},, + 2: (4,0,2, = 0 (3.30) k=1 (:1 This equation can be written in a matrix form as Anannxl : O (331) 41 20 2.5%.... ‘35....“ “3...... 5.5%.... 3.5.... .1 _ “a... .55.... #5....- °.-....... ‘-.....‘. $5”...- '.-._.'..I mu...” _, 120 . 13‘ 1 140 - “'3‘ “5%.. 5.“... "°.. ”‘1. ..“-. $.53 . I 160 - \ - '-._.. '15... 01.... 's... ‘-.._ '5... "a. "-. "-. \‘a '3'. ..". a... x... ‘5. 180 - «a l «'5'...‘ '.-..:...L .555... ..‘u._'...\ 0 50 1 00 1 50 nz = 2772 Figure 3.7: Nonzero entries of matrix A (Level=5). where n = (2j)2, and A = A1 + A2. The entries of the matrices Aland A2 are A1 0;}, and A2,, =Cgf, (3-32) k,p _ To illustrate the structure of the system of equations, the nonzero entries of the matrix A are shown in Figure (3.7) by adding non-symmetric off-diagonal entries. The final system of equations that will provide the solution in the domain of interest is obtained by replacing the equations (3.30) associated with the points (13),, yb) 42 which are in 0w by the following equations: u($bayb) : 9 Where (Eb : E and yb = l 2] 2] (3.33) and g = 2,. d.,k2i¢(2ixb — z')¢(2jyb — k) Appending these boundary conditions (3.33) to the original system has, however, the detrimental effect of destroying the nice properties of the matrix shown in Figure (3.7). In this example, 9 has the values shown in Figure (3.6). Figure (3.8) shows the temperature contours for different levels of dilation. The errors contours are shown in Figure(3.9). One might notice that the solution is not perfectly symmetric, and this is simply due to the asymmetric of the scaling function gb(2j:r — i)¢(2jy — k), as seen in Figure (2.1). 43 Figure 3.8: Temperature contours at different levels of dilation using the Fictitious Boundary Method. 44 Level=5 Level=6 0.4- - o.4~ 0.2» - 0.2. 4 o . - o . ‘ - o 02 04 06 03 1 o 02 04 06 03 1 X X Figure 3.9: Error (u—uemct) contours for different levels of dilation using the Fictitious Boundary Method. 45 3.4 The Penalty Formulation The penalty formulation method arose from the theory of constrained minimization. In this method, a functional I is to be minimized on a space V subject to the constraint Bu = g, (3.34) where B differential operator taking V into another space Q. For a Dirichlet problem, the operator B is equal to 1. This problem can be viewed as one of finding the minimum of I in some constraint set K = {'0 EVle = g}. The idea behind penalty methods for constrained mini- mization problems is to append to the functional being minimized a penalty term, which gets larger in magnitude the more severely the constraint Bu = g is violated. In other words, if a function 22 is tested as a candidate for a minimizer of I, the further 7) is away from satisfying the constraint, the greater the penalty must be paid. A new functional IE is therefore introduced as 1 15(1)) = 1(1)) + E1901), 2) E V (3.35) where P(v) is a penalty functional. P must satisfy certain conditions in order to be considered as penalty functional. These conditions are listed in [4]. It was shown in [4] that us will be a solution to the following variational problem: 1 < 61(u5),v >v +E < 6P(u5),v >v= 0, V7) 6 V (3.36) 46 where 6 is the variation. Implementation of the boundary condition by the penalty formulation has the advantage of introducing no new variables to the system of equation. On the other hand, other methods to implement the boundary conditions such as with Lagrange multipliers require additional unknowns to the system of equations. 3.5 Fictitious Domain / Penalty Formulation The Fictitious Domain Penalty Method to the boundary value problem (3.37) is used in light of work done in [22, 20, 43]. This method consists in first embedding the domain w into a simple geometry, named the domain 9, as shown in Figure (3.10). The Galerkin method is then used to solve the governing equations in the large domain 9 with periodic boundary conditions. A modified functional is used with a penalty term added to the weak form of (3.37) to ensure the implementation of the original boundary conditions. This requires the transformation of some integrals of the form fwu(:r, y) (13 to f“ u(x,y)dQ. The problem to solve can be expressed with L(u(:1:,y)) = 0 in w E R (3.37) Bulaw = 9 where L is a linear second order differential operator, B is an Operator associated with the boundary conditions. For a Dirichlet problem, B is equal to 1. 47 0w u(:1:, u) Ian periodic 652 Figure 3.10: In the Fictitious Domain method, a domain w of arbitrary geometry is embedded into an auxiliary, simple domain 0 with periodic boundary conditions. 48 The differential equation given by (3.37) and its boundary conditions are repre- sented in terms of wavelet functions defined in the bigger domain Q. The solution to original problem in w is achieved by adding a penalty term in w to the weak formula- tion of (3.37) in Q. The penalty term involves a penalty parameter 5 which is chosen to be a positive but small number i.e. 8 << 1 in order for the solution to converge to the desired solution in w. The form of (3.37) that is minimized in the fictitious domain penalty formulation is given by I(u,v) = LL(ua(x)) 11(33) d1: + ‘2}; a (Bun — g)2v(:r)dx (3.38) Taking the variation of [(u, v) (61(u, v) = 0), yields / L(ua(x)) v(:r) d3: + :/ (Bua — g) v(a:)d:r = O for v 6 H1362) (3.39) 9 8w where H1162) is the space of H1(Hilbert space) functions in O which are periodic in :r and v is the test function which was chosen to be the Daubechies scaling function with six filter coefficients. This minimization procedure will be discussed in more detail in chapter four. It was shown in [43] that in the limit when 8 ——> 0, ua converges to the original solution u, the desired solution, which satisfies the boundary conditions in 0w. The second integral term of equation (3.39) is evaluated in the larger domain by using the 49 following formula [43]: f ds=/ f Haw” d!) (3.40) 6w R2 where Haw” is called the boundary measure function [43] (see Appendix B for an expression of Haw“). In one-dimensional problems, the fictitious domain is a line of length 23', where j is a fixed integer that specifies the level of dilation in the approximation. It will be referred to each unit in this line by point k. In the two-dimensional case, it was assumed that the fictitious domain 9 is a square of dimensions 2j x 23'. This is done without loss of generality because any domain can be mapped with a simple change of variables. Equation (3.39) can be written in terms of the boundary measure function “aw“ ofw [22,43] as l 1 f9 L(ua(:r)) m) dx + E Luau. — g) v(:c)||6w|| dz = 0 W e 11,,(9) (3.41) an and g are expanded in terms of scaling function. The test function 1) is chosen to be ¢(2j:c — k) for one-dimensional examples and is chosen to be ¢(2j:z: — p) (M2321: — q) for two—dimensional problems. The expressions of ac, g, and v are then substituted into the weak/fictitious penalty formulation equation (3.41). This will yield a discretized form of (3.37) in terms of the connection coefficients 0:25;". 50 3.6 Solution of problems with the Fictitious Do- main Penalty method The same examples as in the previous sections are used for the purpose of comparing the Fictitious Boundary method and the Fictitious Domain Penalty method. 3.6.1 One-Dimensional Poisson Equation For the one-dimensional Poisson equation (3.8) the weak penalty fictitious domain formulation is , , 1 ' fnu (:1:)v (2:)dx + E f()(uh‘) — g)v(z)||3w||d:r = fnf(x)v(x)dx (3.42) Taking 22(3) = (25(2ja: — k), g(:r) = Z, g,2%¢(2jx — 2'), and substituting the wavelet expansions (310,311) for u and f (3:) yield the following discretized Poisson equation: . 2' 3d,. ifk c 8w . ggk ifk c 8w 21 2 30,12, + = (2.2% + (3.43) i=1 0 elsewhere 0 elsewhere where 7 is chosen to be equal to 1. Since the solution is periodic in Q, the connection coefficients Cfi’?’ have non zero values for Daubechies scaling function with three vanishing moments when —4_<_k—z'g4 (3-44) 51 OI' k—4gz‘gk+4 (3-45) Also, the summation in equation (3.43) can be simplified as _ k+4 1d,. k c 6w . 19;. k C 6w 21 Z d,C,1’_1,.+ 5 = (2sz + ‘ (3.46) i=k‘4 O elsewhere 0 elsewhere or in matrix form 1 1 . . Am", an1+ EPan an1 = rnx1+ Esnxl, where n = 2’ (3.47) where r, and s are the forcing vectors, P is diagonal matrix of entries of either one or zero, and X contains the unknown coefficients d’s. The matrix A has the structure shown in Figure (3.11). A is a circulant matrix with entries of the connection coefficients 0:13,. The entries of the vector 3 are simply a, 5, or zeroes. This problem was solved with the same values of a, B, and f (as) as in example 1 of the first method with penalty factor 8 equal to 10’5. Figure (3.12) shows the solution of this problem at different level of dilation. The error of the approximate solution is shown in Figure (3.13) at different levels of approximation. 52 ooooo coco 2-000000 000 0000000 00 4-00000000 o 10» 000000000 12» o o o o o o o o o o o o o o o o o o 14- o o o o o o o o o o o o o o o o o o 16- o o o o o o o o o O E j 6 3 1‘0 112 1‘4 {6 nz=144 Figure 3.11: Nonzero entries of matrix A (Level=4). 53 —— wavelet solution exact solution Level=5 Level=6 1 . . r f 1 . - . a 0.8) 0.8 - 0.6 ' 0.6 > 3? 3? ‘5’ ‘5 0.4 - 0.4 . 0.2 - 0.2 ’ O A A A A O A A A 0 0.2 0.4 0.6 0.8 1 O 0.2 0.4 0.6 x 0.8 Figure 3.12: Fictitious domain/penalty solution of the 1D Poisson Equation. 54 0.025 T , , , , , r q 1 ................. LeveIé ' — LeveI=6 0.02~ _ 0.015 — ‘ § 0: 0.01 - - Figure 3.13: Error (u — uemct) in the wavelet solution of the one-dimensional Poisson equation at different levels of dilation. 55 3.6.2 Simple Harmonic Motion Equation The Galerkin weak/ penalty formulation of equation (3.18) is given by the following expression: f(u(:r:1:)u+wnu((z)v))d$+-/(u)v:1:( )||6w||d:z:=0 (3.43) Taking v(:1:)- - (25(2jx — k), (at): Z g,22¢( (25(2333 — 2'). Substituting the wavelet expansions (3.10) into equation (3.48) and taking the advantage of the periodicity nature of the solution in the fictitious domain (2, allows equation (3.48) to be written as k+4 11 w" id).c k C 602 fig;c k C (902 Z al.-Ck;- + ”27"" + = (3-49) i:k-4 O elsewhere 0 elsewhere Equation (3.49) can also be written in matrix form as 1 1 j Anxn anl + Eann anl : Esnxl Tl : 2 where s is the forcing vector associated with the boundary conditions, X contains the unknown coefficients (1’3, and P is a diagonal matrix of entries of either one or zero. The matrix A has a structure similar to the one shown in Figure (3.11). It is a circulant matrix with entries made of the connection coefficients CE). The entries of the vector 3 are simply a, B, or zeroes. Figure (3.14) shows the solution to this example with penalty factor 5 equal to 10'7. The error in the solution is presented 56 Level=5 Level=6 1 .5 . . . . 1.5 . a . X0) x(t) 0 0.2 0.4 0.6 0.8 1 t Level=7 Leval=8 1 .5 . r . . 1 .5 . . . x(t) x(t) 0.2 0.4 0.6 0.8 1 '0 0.2 0.4 0.6 0.8 1 1 t Figure 3.14: Solution to the Simple Harmonic motion equation using the Fictitious Domain/ Penalty Method shown at various levels. in Figure (3.15). 3.6.3 Steady State Heat Conduction Equation Equation (3.29) can be written in a weak/ penalty form as Bu 8'0 Bu 61) 1 [(5555: + 07/537 ) dxdy + E f(u(:r,y) — g)||8w|| datdy —- O (3.50) 57 0.6 I r 7 f l l r 1 I ........ Ieve'é ° - level=6 0.4 - — — — level=7 ‘ —-—— Ieve|=8 Error —O.4 r- ~ -0.6 " " “0.8 l l l l l l l l l 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.15: Error (:1: - mama) in the Wavelet Galerkin Fictitious Domain/ Penalty Method solution of the Simple Harmonic Oscillator equation. 58 u and 9 can be written in wavelet space as 21 240:, y) 2 2: 4.4214034 — 04423.4 — k) (3.51) i,k=l and 9(1', 11) ’-‘-‘ Z 94,42’¢(ZJ$ - z')¢>(2’y — ’6) (3-52) i,k=l Substituting expression (3.51) and (3.52) into equation (3.50) and taking v(:r, y) = (15(2303 - p)¢(2j y -— q), yield the following discretized heat equation: q+4 P+4 551d” if (p,q) C (90) 19” if (p, q) C 30) Z dk,pC;'_lk + Z dq,(C;:_ll + = E k=q-4 l=P-4 0 elsewhere 0 elsewhere (3.53) Equation (3.53) can be written in matrix form as ARXanx1+§Pananl = is“), where n = 22f (3.54) 01' AX = b (3.55) where A = A + P, and b = is. The matrix A is block circulant and has the structure shown in Figure (3.16). 59 100 150 250 R zoolha 50 $23.... 3,55%“... 3.35%... «Ea-u... '-......_ 1 00 150 200 250 nz = 4352 Figure 3.16: Illustration of the nonzero entries of matrix A (Level=4). 60 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 3.17: Temperature contours at different levels of dilation using the Fictitious Domain/Penalty Method. 61 Level=5 Level=6 7% 1M fl O'BU 08> w 0.6’ 1 0.6[ >. > 0.4' ‘ 0.4’ 0.2’ ' 02* 0‘ 4 0 A 0 02 O4 06 08 1 0 02 04 06 08 1 X X Figure 3.18: Error (11 — new“) contours for different levels of dilation using the Fic- titious Domain/ Penalty Method (5 = 10‘s). 62 X 10‘3 L6V6l=5 X 10'3 L9V9|=6 1.2 ~ ~ 1.2 L J 1 - 1 1 ‘ 0 8 - 0 8 > ~ E To 3 D 9 .12 8 8 m 0.6 ' '1 m 0.6 l '1 0.4 d 0.4)1 - 0.2 ~ - 0.2 ~ - 0 l A L“ l 1 0 4L ‘ 1 0 50 100 150 200 250 0 50 100 150 200 250 Iterations Iterations Figure 3.19: Convergence of the solution using the Conjugate Gradient method. 63 The diagonal blocks of A have a structure similar to the one shown in Figure (3.11). The penalty matrix P is a diagonal matrix of entries zero or 1. The vector 3 has entries of zero or gm which could be found from equation (3.52). For this example the same boundary conditions values of example (3) of the first method are used, therefore the entries of the vector 3 is either zero or 5,17. The system of equation (3.54) was solved iteratively using the standard conjugate gradient method [36] (see Appendix (C) for algorithm). Iterations were stopped when the residual, define as Mfg—bl, was found to be 10‘7. Figure (3.17) shows the temperature contours at different levels of dilation. The error contours are shown in Figure(3.18). The residual history is shown in Figure (3.19) for different levels of approximation. It can be noted from Figure (3.19) that the solution converges rapidly with the first fifty iterations. The system of equations given by (3.54) is also solved with the preconditioned conjugate gradient method. Different preconditioners were tried to implement the preconditioned conjugate gradient scheme. The selection of the preconditioner ma- trices is based on two aspects. First, they must be positive definite. Second, they must not require significant effort to invert. Some of the preconditioners tried are: . M=A+c [1 1 1---1]T><[1 1 1...” o M = A + diagQ) o M = diag(A) where M is the preconditioner matrix, c is a constant chosen in such a way to make M positive definite, A and A are defined previously. All these three preconditioner 64 are positive definite and easy to invert. Of these, the only preconditioner that gives a reasonable reduction in the calculations is when M = diag(A). Table (3.1) shows the performance of this preconditioner in preconditioned conjugate gradient (PCG) method against the standard conjugate gradient (CG) method without a precondi— tioner (see Appendix (C) for algorithm). Table 3.1: Performance of PCG method and CG method I Level CG Iterations PCG Iterations I: 42 27 E 110 [ 62 ] [ 6 242 l 127 | f 7 F102 T256 ] f 8 [ 776 | 370 ] 3.7 Summary Two methods for solving differential equations with Dirichlet boundary conditions using a wavelet—Galerkin method are reviewed and discussed, namely the so—called “Fictitious Boundary ” and “Fictitious Domain/Penalty” methods. Evaluation of the two methods is performed by solving several one- and two-dimensional examples using fixed scale expansion of the unknowns. For each problem the error is computed so that the accuracy of the solution can be evaluated. It is found that the Fictitious Domain/Penalty method shows better agreement with the exact solutions than the Fictitious Boundary method as it introduces less computational errors due to the methodology used to implement the boundary conditions in the extended domain. 65 Also, the Fictitious Domain / Penalty method is superior since the resulting system of equations is block circulant. The Fictitious Domain/ Penalty method is amenable to its use within the framework of multiscale techniques. 66 Chapter 4 Solving Newtonian Fluid Flow Problems Using a Wavelet Galerkin Method 4.1 Introduction In this chapter the formulation and discretization of the equations governing the behavior of Newtonian fluids using a wavelet-Galerkin are presented. The fictitious domain/penalty method, which was presented in chapter three, is used to formulate the modified governing equations. Formulation of these problems using this method allows efficient and accurate solutions since the fictitious domain'penalty method introduces less computational errors than the fictitious boundary method. Also, this method allows solution to large problems using modest resources such as a personal computer. In this chapter, we consider the solution of a creeping flow problems, 67 i.e., flow with low Reynolds number which is defined as Re 2 E}, where L and v are characteristic length and velocity, )1 is the viscosity. The methodology is evaluated by solving various examples. Those examples are the lid-driven cavity fluid flow problem and flow between rotating concentric cylinders. 4.2 Governing Equations The description of the motion of a continuous medium is based on the conservation of mass, conservation of momentum, and the conservation of energy. The associated equations of state and constitutive relations are also essential for quantifying the motion of a continuous media. To reduce the complexity of the problem, certain assumptions are considered. These assumptions are as follows: 0 Steady state fluid flow 0 Incompressible flow Isothermal flow 0 Two-dimensional flow Creeping flow (Re 3 O). The equations of interest for a flow under the previous assumptions are provided below [6]. The first equation is given by conservation of mass V-sz (4.1) We then have the conservation of momentum pv-Vv+V-T+Vp=0, (4.2) where p is fluid density, v is the fluid velocity, 7' is the shear stress tensor, and p is the pressure. The momentum equation can be written in dimensionless form as Rev-Vv+V-T+Vp=0 (4.3) OI‘ V-T+Vp=0 ifRezo (4.4) where Re is the Reynolds number and is defined as Re = %L, where L and v are characteristic length and velocity, and p is the dynamic viscosity. For Newtonian fluid flow the shear stress tensor 7’ is given by 7' = -u’7 (4.5) where is ’7 is called the rate of shear strain tensor. For two-dimensional problems ’7 is 1: (W) + (Vv)T (4.6) 69 or in Cartesian coordinates . 7n ’ny 1: (4.7) ”7.1/1 Vyy and 244 04 + 494 ‘ 8 8 8 7: I I y (4.8) Q2 @ .31 By + 61: 28y where u and v are the :1: and y-components of the velocity vector v. The substitution of equations (4.5,4.8) into equation (4.4) results in the following equation of motion in scalar form. The continuity equation is given by Bu 81) and the x-momentum and the y-momentum are written as 6211 8211 0p I“? 5372') - '5; - (4-10) (922) 621) 8p 14$ 5y?) — 8—3; - (4-11) 70 Periodic boundary conditions on 89 BQ Problems 0) and subject to the boundary conditions v=gin8w 71 Figure 4.1: In the Fictitious Domain method, a domain 02 of arbitrary geometry is embedded into an auxiliary, simple domain Q with periodic boundary conditions. 4.3 Fictitious Domain Method for Fluid Flow We consider a fluid governed by equations (4.9,4.10,4.11) in a closed domain called (4.12) This problem is solved in a larger periodic domain (2 subject to periodic boundary conditions in 89 as shown in Figure (4.1). The original boundary conditions are applied via a penalty formulation. This method was discussed in chapter three. The essence of this method is illustrated in Figure (4.1). In the implementation process, the fictitious domain Q consists of .2" points in both a: and y directions, where j is the level of dilation. The dots in the bottom diagram of Figure (4.1) represent the points where the original boundary conditions in 60) are applied using a penalty formulation. 4.4 Penalty Formulation of the Governing Equa- tions The weak form of equations (4.9,4.10,4.11) in the fictitious domain are as follows: Bu 00 [9(5; + 5;) 1141311) d0 = 0 (4.13) 8211 6221 0 [1045“; + a; -' 55—)w1($,y)d9 = 0 (4.14) 82v 62v 6 [004555 + 6?) — 5;) 101(33, y) (19 = 0 (4.15) 72 where w(:r:, y), w1(:r,y), and w2(:r,y) are the test functions and 10,101,102 6 Vp. A typical choice for VP is H1(Q) and the space 140(9) is defined by V,,(f2) = {w(:r, y) E H1(Q) : w(:r, y) is periodic on 60} (4.16) To reduce the order of differentiation on the variables (11, v, p), integration by parts for equations (4.14,4.15) is needed, that is :0 (due to periodicity) =0_(due to periodicity) ‘ Bu Bu p—wl [if —fy f pig—Eldxdy+ /p —'w1 3’: d2: fy 8:1: 6 a ”By ) (4.17) =0 (due to periodicity) — 1.132186121122114 (,pr 1:: 44 +f f pig-""61“?! = 0 . Therefore, the weak form of the x-momentum equation is Bu 8w] Bu Owl “0101 —— dQ= —dQ . #1582: 62: +_ 8y 0y —) 0:1: (418) Similarly, the weak form of the y-momentum equation is 012 3102 80 ng pawg —— d9: — dQ ’ 4. [9((912 8:1: +— 8y 8y ) pay ( 19) In the weak forms equations (4.18,4.19), we have the following linear and bilinear 73 forms: Bu(u,w1) = l(w1) Bv(v, 1122) = [(102) where 1(w1)= 9.05:,ng [(2122): 917981133619 Bu(u,w1) = “[5253ng + g15§agy1 ) d9 and _ 61)an 61) 8202 Bv(v,w2) ”‘wa 6:1: + By By )dQ (4.20) (4.21) (4.22) (4.23) (4.24) (4.25) This problem can be reformulated as an unconstrained problem by modifying the bilinear equations (420,421) in a penalty formulation. The quadratic functionals of 74 equations (420,421) are given by references [32] and [33] Of -2 222 222 _/,, 22 1(“)‘2/,,(ax) +(6y)d§2 a pa do 811 2 — =-gffl(g:—t)2 g—S) (19— 0128de The minimization of the following functionals is now considered: 0(4) = 1(a) + 51; [<4 — 40211841140 1,,(3): Iran-215- [(4 — g)? 114421140 (4.26) (4.27) (4.28) (4.29) (4.30) where e is small positive number, [[802]] is the boundary measure [43] which is discussed in Appendix B. The necessary conditions for the minimum of Ip(u) and 1,,(0) are (4.31) __ Bu 0(611) Bu 8(624) 8(611) 1 __ arr/”(116$ 6:17 “‘03, 33, p 82: +E(u g)6u|[6w||)dQ—0 (4.32) It should be noted that 611 = 101 and 62) = 102, therefore the weak form of x-momentum equation, with a penalty formulation, is given by 81181111 624 61111 [a Owl p/fl(ax—— Ba: +6—y— By )d§2+— :/fl(u—g)w1H3wlldQ= 196de (4.33) Similarly, the weak form of the y-momentum written in penalty formulation is 81) 8__w2 82) (9—102 fay 6102 ,u/ 0:1: 31: +33, 37, )dQ+ i/nw— g=)w2||0w||df2 p dQ (4.34) 4.5 Wavelet-Galerkin Method and Discretization of the Governing Equations The first step in applying the wavelet-Galerkin method is to expand the dependent variables 11, v, p in terms of scaling functions as follows: N—l N—l ”(33 11):: Z Uik¢(X -k) (4.35) i=0 k: 0 N-l N—l ”(33, y) = Vik¢(X — 0450/ — k) (435) j=0 k=0 76 2 l. 2 .1. 19(1‘, 31) = EMMX - 0950’ - k) (437) 0 H o a. H i where X = 21:12, Y = 213), N = 21', U,,. = 210,-,“ V,,. = 210,-. , P,,. = 21pm, and j is the level of discretization. In the Galerkin formulation the test functions w’s are chosen from the same space as u, v, p, i.e. w = w] = 1122 = (15(X -p)q§(Y — q) where p and q are integer translates of the scaling functions 05(X) and ¢(Y) respectively. The boundary conditions in equation (4.12) is also expanded in wavelet space as 90.4): 2:131 123‘ emu—040%) (4.33) Gm : ngpq : ff9($ay)¢(x ‘ P)¢(Y "' q)d:1:dy The coefficients U’s, V’s, P’s are found by substituting equations (435,436,437) into equations (413,433,434). The discretized continuity equation is given by p+4 q+4 Z Uiqufi, + Z Vkafik = (4.39) i=p—4 kzq—4 where the unknowns are the coefficients Uiqandek. 77 The discretization of the x-momentum equation is as follows: :C“_. 25kg #251323:_010../3(x-)(xmax/3(2— (040-044 :6, =03; + 4422; .0: U0. ¢(X -i)¢(-X p)dx ¢(Y- k)¢(Y- aldy / / 3% + 221201 k: 0 Uikllaw” ¢(X"1)¢ p)d:r (MY —k)¢(Y—-q)dy / _ / 1,5 — 212:3250104 ”awn/4m -:)¢(X-P)d$/¢(Y-k)¢(Y-'Q)dy :C‘IOi =6kq — 21:12:53.)“ —4)4(X— p)dx/¢(Y— —q)dy= o (4.40) The delta function 6 will cancel one of the summation in each term of this equation (4.40). For example, 6kg will eliminate the summation associated with the index k. Also taking advantage of the compact support of the connection coefficients 0,21,- and the periodicity of the solution in 9, equation (4.40) can be written in a simpler form as p+4 q+4 p+4 8 8 #( Z Uiani + 2 Up kCllk) +____ [____l :HUP :2 P14101391 + ”—9qu (441) i=p—4 kzq- -4 i=p— —4 Similarly, the y-momentum can be written as pH 11 (1+4 11k +L—la:”Vp 0+4 10 [lawllG “(Z 14ch + Z V. WC :2 1).-,0 _,.+———G ,, (4.42) i=P'4 (C: q- -4 1:: —-q -4 The indices p, q are chosen to vary from 0 to 21 — 1, also the coefficients qu are 78 normally found from the boundary conditions of the problem. 4.6 Solution of the Resulting System of Equations Equations (439,441,442) form a complete set of simultaneous linear equations in terms of the unknown coefficients U ’s, V’s, P’s. This set of equations can be written in matrix form as A B‘ v G = (4.43) B 0 p 0 This system of equation has been studied in the literature extensively. For exam- ple, the work in the references [3, 7, 25, 36] solved this system direct. On the other hand, there are number of studies in the literature that solved this system iteratively by omitting the continuity equations and iterating over the pressure. Some of these methods are Variant Methods, Conjugate Gradient methods, Penalty methods, Decomposition methods, and Augmented Lagrangian methods [17, 21,18,19,41] In this study, two methods are considered. These methods are: the conjugated gradient method, and the steepest descent method or sometimes called the Gradient method with fixed step size. These methods are based on decomposing the system of equations (4.43) into a series of Poisson equations in the velocity components 11, and This approach is considered in part because of the deveIOpment of an efficient and 79 fast solver for the Poisson problems was achieved in chapter three via the wavelet- Galerkin fictitious domain penalty method with a Preconditioned Conjugate Gradient (PCG) solver. The system of equations (4.43) therefore requires a preconditioner for an efficient solution. This preconditioner to deal with the presence of the penalty entries in the diagonal of A. These entries cause problems during the inversion of A due to their large difference with other non-penalized entries. 4.7 Solution Methods In this section a conjugate gradient and the steepest descent methods are presented and discussed. A comparison of the two methods is available with solving examples with both methods. 4.7.1 A Conjugate Gradient Method (CGM) The algorithm for solving the Poisson equations (410,411) via a conjugate gradient method is based on the work done in [17, 18]. Such algorithm involves the following steps: 0 Step 0: Initialization p0 E L2(Q) arbitrarily given (4.44) 80 then solve uAvoz Vpo onw,v=gin8w (4.45) by the wavelet-Galerkin fictitious domain penalty method in Q, i.e. solve pH 0 11 (1+4 11 +_[__lawllUp pH 0 10 +||8w|| 423 4,3,. .- 2U 22> 4 =2 9.2+ . _..~.. =P-4 k=q- -4 i=P-4 (4.46) p+4 q+4 q+4 80) 80) ..( Z 1430,11 + Z 143,011+ +” “14,:2 123,033+ —” 5 “GP, i=p—4 k=q—4 k=q-4 (4.47) and set 20 = V - v0 (4.48) 0° = 2° (4.49) Then for n > 0, compute p"+1,z"+1, y"+1 from p", z", y" as follows 0 Step 1: Descent Zn 2 ll ”132(0) (450) — (a y", y")L2(n) 81 "+1 = p" — any" (4.51) check for convergence 0 Step 2: New descent direction z"+1 = z" — anay" (4.52) H Zn“ ”12(9) 1" z n znni’ (4'53) 142(9) yn+1 = Zfl+l + 7713/” (4.54) n = n + 1 go to equation (4.50). To implement this algorithm, it is necessary to know ay". From Theorem 5.10 in reference [17], one can evaluate ay" by solving the following Poisson equations in x subject to homogeneous boundary conditions, that is qu": Vy" in w, x = 0 in 6w (4.55) ay" V - x" (4.56) Thus each iteration requires solving two uncoupled Poisson problems. 82 4.7.2 The Steepest Descent Method (SDM) This method is easier to implement than the conjugate gradient method presented previously. The algorithm for this method is as follows: p0 E L2(§2) arbitrarily given (4.57) for n 2 0, compute 13"“ by solving ,uAv“= Vp" in Q (4.58) 10"“ = p" - puV - V“ (4.59) where N is the number of dimensions of the problem, e.g. for two dimensional problems N = 2, p is a damping parameter which can be chosen form the range 0 < p < —,%[18]. For the examples which were solved p was chosen to be equal to 0.6. In this method, at each iteration, one needs to solve two uncoupled Poisson problems. In the implementation stage, it should be noted that these algorithms are written in terms of scaling function expansions. For example, the steepest descent algorithm can be written in scaling function coefficients as follows: ng E R arbitrarily given (4.60) 83 for n 2 0, compute PM.+1 by solving equations (446,447) via the wavelet-Galerkin fictitious domain penalty method. Then update the pressure coefficients PM by p+4 q+4 n 1_ n n 10 n 10 Pp: —— qu —- pp ( Uiq Cp_,- + E pk 0,14,) (4.61) i=p—4 k=q—4 4.8 Examples The examples that are considered in this dissertation are the lid-driven cavity problem and flow between two concentric cylinders with the outer cylinder rotating. The first example, the lid-driven cavity problem, was solved by both SDM and CGM methods for low levels. It was then solved for high levels by the SDM method. The second example, i.e., flow between two concentric cylinders, was solved via the SDM method. 84 It“ 4.8.1 .Lid-Driven Cavity Problem This problem was solved with the boundary conditions shown in Figure (4.2). With- out loss of generality, the dimensionless viscosity p is set to equal to 1. The results are provided by velocity vectors for low levels of dilation, and streamlines for high levels of dilation. The streamlines are calculated by a simple linear differencing scheme. For example, from point 11 to point n + 1 the stream function has the form 1 1 wn-H — wn‘ ‘- 2(yn+1 — yn)(un + “n+1) _ 5(xn-H " $n)(vn + vn+l) (4-62) where w is the stream function, u is the x-component of the velocity, and v is y- component of the velocity. The residual criterion for termination of the computation is defined as Residual = Z (122+1— Pg.) (4.63) i,j=0 The residual history is provided for the SDM method. A comparison of the residual history for both SDM and CGM is provided for level=5. 85 v=(1 0)‘ v = (0 0)‘ V: (0 O)‘ Va v = (O 0)‘ Figure 4.2: Domain and boundary conditions for the lid-driven cavity problem. Conjugate Gradient Method Steepest Descent Method ' . v ' 1m\v ”W 1\ l ~ I I ”—.-.~”“‘ \ ‘ ' a’»»-.-._.—.—.“—.‘ I l I I I I l ’ ‘ ‘ ‘ ~ 1 I l I \ I I a .. .., .. ~ ~ \ ' u l ‘ l t t 0-8: . 1 , . 0-81 '411111 , . . ~ | | l \ o I l ‘ \ \ \ s - I I / I I I ' c Q \ \ a I o | \ \ \ \ \ . .., r I I I l I . 0H6 \ \ \ s . . , , , , 06. . \ \ \ \ x - .. I I l ’ I l H s \ s s - , ' I I o \ \ \ \ \ ~ - I I I I I ’ 4 s s .. - , , , ‘ ‘ ‘ ‘ "“ " " v' I I l I a 0.4' ~ ~ - - . , 0.4’ ‘ ‘ ‘ ‘ "' ' I I I I 0.2' 0.2. . . . . . . - - . . 0 ‘ ‘ * ‘ “ o . . . . . O 0.2 0.4 0.6 0.8 1 o 0.2 0.4 0.6 0.8 1 Figure 4.3: Comparison of the velocity vectors between SDM and CGM for the driven cavity problem (level=4). 86 Steepest Descent Method Conjugate Gradient Method 1 .2 l fir I I 1 .4 I l I r 1 - l 1.2 - - 0.8 "' :> - 1 - . .:’ (0.5009,o.7572) 2 0.6 - ‘ 0.8 - _-;’ ' : (05000.0.7570) g 0.4 ~ - 0.6 3% - a: = 2 0.2 - E - 0.4 ~ :1: 2 :2. 1L: : - - 200 - - 200 - - m E. 150 ~ - .§ 150 - - :2 2 2 2 o o 8 0 a. 100 - ~ 100 » J 50 ~ - so - - 0 J_ 1 1 0 J 1 l o 200 400 600 o 200 400 600 SDM Iterations SDM Iterations Figure 4.6: Iterations required for the PCG solver along the SDM iteration for the u component of the velocity. 89 Leve|=6 Levei=7 300 fi fi . 300 I . . 290 - j 290 - - 280 - - 280 ~ 4 m 0) .5 270 - - .5 270 ~ ‘é E 2 2 8 8 n- 260 " d o- 260 ' " 250 ~ ~ 250 - 2 240 - - 240 - - 230 l l 1 230 l l l 0 200 400 600 O 200 400 600 SDM Iterations SDM Iterations Figure 4.7: Iterations required for the PCG solver along the SDM iteration for the 2) component of the velocity. 90 SDM Level=7 SDM Leve|=6 1F ,2 , ,. --' - -' \ ‘ ‘ ‘ 0 0 0.2 0.4 0.6 0.8 1 O 0.2 0.4 Figure 4.8: Streamlines obtained by the SDM for the driven cavity problem (level=6,7). - level=4 5 - level=6 — level=5 — level=7 Residual Residual 0.5 0 200 400 600 o 200 400 600 800 Iterations Iterations Figure 4.9: Residual history of the SDM at various levels for the lid-driven box. 91 Va -u H V9 1 Figure 4.10: Decomposition of the tangential velocity into Cartesian coordinate sys- tern. ‘ where V9,, is the outer tangential velocity, V9,. is the inner tangential velocity, To is the outer radius, and r, is the inner radius. The results for this example are in terms of velocity vectors and streamlines. The wavelet solution for this example shows some radial flow in some regions of the domain especially when lower levels considered, this might be due to the representation of the boundary is not accurate. The represen- tation for the boundary in this example is achieved via the Haar wavelet system. A better representation is to use a smoother wavelet function, e.g., Daubechies wavelets. 92 Level=5 Leve|=6 1 0.5 ~ 0.5 0 - 0 .05 —0.5 “11 —o:5 0 0:5 1 '11 —o:s 6 0:5 1 Figure 4.11: Velocity vectors obtained by the SDM for the flow between two concentric cylinder problem (level=5,6). Level=5 Leve|=6 1 1 . 0.8 0.8 ) 0.6 o. 0.4 > o 4 0.2 o 2 -o.2 -0 2 4,4 —O.4 —0.6 -0- —o.a -°-3 ' -o.5 o 0.5 1 -o.5 o 0.5 1 X X Figure 4.12: Streamlines obtained by the SDM for the flow between two concentric cylinder problem (level=5,6). 93 Leve|=8 Level=7 1 ~ 1 0.81 0.8: 0.6 r 0.6 ' 0.4 r 0.4 r 0.2 ' 0.2 r >. O 2 >. 0 l -02 . -O.2 ’ -O.4 - -0.4 * -O.6 » -0.6 r —0.8 - -0.8 ) —O:5 0 04.5 1 -O:5 Figure 4.13: Streamlines obtained by the SDM for the flow between two concentric cylinder problem (level=7 ,8). Level=5 --- Exact Solution — Wavelet Solution .0 .o .o b O) on Velocity Magnitude .0 M Leve|=6 1 . . 1 0.8 ,‘ m . O .é' 1 g 0.6 r 1 (B 2 .é‘ 8 0.4 .I B I > . _I 0.2 - ., 1 O A A 1 -1 0 0.5 1 x Figure 4.14: Comparison of the velocity magnitudes of the wavelet and exact solutions at the crossection y = 0, (0 = 0, 7r in polar sense). 94 -—- Exact Solution -— Wavelet Solution Level=7 Level=8 1.4 1.4 . . 4 1.2 . 1 2 § ‘ ’ § ‘ ’ 310.8 310.8 1 2 2 $0.6 E 0.6» § § 0 d) > 0.4 > 0.4 0.2 1 0.2 l- l 0 ‘ o . -1 0 1 -1 0 1 X x Figure 4.15: Comparison of the velocity magnitudes of the wavelet and exact solutions at the crossection y = 0, (0 = 0, 7r in polar sense). 95 4.9 Discussion and Summary In this chapter, the solution of viscous fluid flow problems is achieved via the wavelet Galerkin method. The modeling of this problem is simple, since it does not require geometric decomposition of the computational domain. The solution methods consid- ered, for examples SDM, gives good results compared to the exact solution as shown in Figure (4.14,4.15). For the lid-driven cavity problem, the solution agrees well with other solutions found in the literature, e.g., reference [9]. Figures (4.6,4.7) show the iterations required, at each step of the SDM, to solve the resulting system of equations by the Preconditioned Conjugate Gradient (PCG) method which was used in chapter three. For the 1) component of the velocity in Figure (4.7) the number of iterations for the PCG slightly decreases as the SDM iterations continue. On the other hand, for the u component of the velocity shown in Figure (4.6) there are nearly four regions in the plots (4.6). Those regions are: in the beginning, there is a constant number of iterations, then there is a dramatic decrease. Then this is followed by a slightly gradual increase in the number of iterations. The last region is a constant number of iterations. The residual behavior of the SDM is very similar for the examples which were solved. Both solution methods, the SDM and CGM, give very similar results for the problems considered. However, the residual of the SDM is gradually decreasing as we iterate over the pressure. Conversely, the residual of the CGM decreases dramatically the first few iteration but then keep decreasing very slowly. It seems that the CGM is slightly costlier to implement. 96 Chapter 5 Solving N on-N ewtonian Fluid Flow Problems Using Wavelet Galerkin Method 5.1 Introduction Fluids that are not described by the Newtonian relations are commonly encountered in a wide variety of industrial processes. Some examples of non-Newtonian fluids are motor oils, blood, high molecular weight liquids, pastes, and polymers. The process- ing, operations, and transportations of such fluids are vital and essential problems in the food, plastics, chemical, petroleum, and polymer industries. Non-Newtonian behavior is identified in a number of different ways. Most of these fluids have a viscosity that is dependent on the shear rate. A decreasing viscosity with an increase in the shear rate, which is known as “shear thinning”, is the most common 97 behavior of such fluids. Elasticity and memory of the fluid (recoil for example) are observed in many situation for these fluids. In these fluids, differences in the normal stress components are also encountered in many flows. This lead to well-known effects such as rod climbing, known as Weissenberg effect, and the curvature of the free surface in an Open channel flow. A comprehensive list and discussion of these and other non-Newtonian effects are explained in the book by Bird el al.[6]. Non-Newtonian fluids are divided into two distinct categories: 0 inelastic fluids or fluids without memory 0 viscoelastic fluids, in which memory effects are significant The inelastic fluids can be viewed as a generalization, in some sense, of the New- tonian fluids. The viscosity in this category depends on the rate of deformation of the fluid and, thus, allows for the shear thinning effects to be modeled. Viscoelastic fluids, on the other hand, represent a significant departure from the Newtonian limit in terms of both the physical behavior and computational complexity. The primary difficulty in solving these flow problems is to account for the “memory” of the fluid. In other word, the motion of the fluid particles depends not only on the present stress rate but also on the deformation history of the fluid particles. In this study only the first kind of the non-Newtonian fluid, that is the inelastic fluid, will be considered. 98 5.2 Governing Equations The governing equations for a non-Newtonian fluid differ from a Newtonian fluid in the constitutive relationship between the stress tensor and the rate of strain tensor. The quantification of a non-Newtonian fluid flow is based on the conservation of mass, conservation of momentum, and the conservation of energy. In this chapter, a number of assumptions will be considered to reduce the complexity of the flow conditions. These assumptions are: c Steady state fluid flow 0 Incompressible flow Isothermal flow Two-dimensional flow Creeping flow The equations of interest for a flow under the previous assumptions in Cartesian coordinates and vectorial notations are as follows [6]: o conservation of mass V-v=0 (5.1) 0 conservation of momentum pv-Vv+V-T+Vp=0 (5.2) 99 where p is fluid density, v is the fluid velocity, T is the shear stress tensor, and p is the pressure. The momentum equation can be written in dimensionless form as Rev-Vv+V-T+Vp=0 (5.3) 01' V-T+Vp=0 ifRezO (5.4) where Re is the Reynolds number and is defined as Re = 1%, where L and U are characteristic length and velocity, 770 is the zero shear rate viscosity. 5.2.1 Constitutive equations There is no mathematical model that can actually cover all non-Newtonian fluids under all flow situations, unlike the Navier Stokes equations, which cover all New- tonian fluids. However, among the numerous constitutive models and for different types of non-Newtonian fluids, some are more popular than others. As mentioned previously, the viscosity in the non-Newtonian fluids have a dependency on the shear rate. This dependency leads to complex constitutive equations. The simplest and the most familiar non-Newtonian viscosity model is the power-low model which has the form 17 = m(2I‘)"’l (5.5) 100 r = [% (1:1)] % (5.6) where I‘ is the rate of strain tensor ’7 which is defined as Bu 6v 6n - 72:2: 73:31 2— '7; + _ 1: = 6x a a” (5.7) 611 6v 6v 71/1: 7.111; 733 + 5} 25:1; and m is a parameter called the consistency index with units of Pa - s" and n is dimensionless exponent. When n = 1 and m = p, the Newtonian fluid is recovered. If n < 1, the fluid is said to be shear thinning and if n > 1, the fluid is said to be shear thickening. The shear stress is given by the following relation: 7' = - WI (5.8) where 17 is the non-Newtonian viscosity which is a function of the magnitude of the shear strain tensor[14]. 5.2.2 Governing equations in scalar form The continuity equation is given by the following equation: 611 011 a; + "a; — O (5.9) 101 Substitute equation (5.8) into equation (5.3) gives Rev-Vv-V-n'1+Vp=O (5.10) where n is the non-dimensional non-Newtonian viscosity 1) = 1136' Therefore, the x and y-momentum components in dimensionless form Bu 011 a Bu 6‘ Bu 6v 8]) _ Re(u5;+v5;)— [55(2n—)+— "(7+— )] +— —0 (5.11) + — = o (5.12) 221222-2211+1221}? where Re, m, 770 are defined previously. 5.3 Fictitious Domain / Penalty Formulation of the Governing Equations The governing equations given by (5.11,5.12) are formulated using the fictitious do- main method which was discussed in chapters three and four. A non-Newtonian fluid flow problem described by equations (5.11,5.12) is considered in a closed domain 62 102 subject to the boundary conditions v = g in 06) (5.14) This problem is solved in a larger, but simpler, periodic domain 9 subject to periodic boundary conditions in 812 as it is shown in Figure (4.1). The original boundary conditions are applied via a penalty formulation as discussed before in chapters three and four. The continuity equation (5.9) is treated as an additional constrained and will be formulated, to be related to the pressure p, also in a penalty fashion as -ep+V-v=0 (5.15) where e is small positive number i.e. 5 << 1. The weak formulation of the governing equations (5.11,5.12) for low Reynolds number Re z 0) is derived as follows: —/9 [5152394256133 3912622243: —w.(scy)=de 0 (5.16) _/[63 zed—3+ 3%))4. gang—v )] 1122(2: y )dQ+/g: —-,w2(:r y)dQ=0 (5.17) where w1($,y) and w2(:r,y) are the test functions and w(:c,y) 6 VP. Typical choice 103 for Vp is H162) and the space V,,(Q) is defined by VP(Q) = {w(:r, y) E H1(Q) : w(x,y) is periodic on 69} (5.18) To distribute the differentiation equally among the variables (11,12) and the test functions 1121 and 102, integration by parts for equations (5.16,5.17) is needed. For example, the x-momentum equation (5.16) is integrated as follows: :0 (due to periodicity) 611 — / 2175—1111(x, y)dy + L, 2172—:6‘” acid!) an 17 :0 (due to periodicity) N ’ 6a 611 ,,, _ [9171;(55+5;)w1(56,y)dx + fan(g—:+%§)%—yldfl =0 (due to periodicity) fl + P1111017, y)dy — fflp—id9= 0 an Therefore, the weak form of the x-momentum equation is —] dQ= ”129312112 / l: 611 6w__1_ 611 612 61121 277— 6:1: — +— 7762: 6:1: + ”(63/ 6:1: 6y Similarly, the weak form of the y-momentum equation is 611+ 6_:6w__2 ”6111/6102 p6w2 (5.19) (5.20) (5.21) From the weak form equations (520,521), the following bilinear and linear forms 104 are found: B(U, 'U, 1121) = [(1111) B(u, v, 1122) = [(102) where 61.01 [(1111) — ”19-8—2:— d9 81112 I w — —d0 (2) ”Pay _ 611 6w1 611 and 611 611 61122 B(uivaw2):[2 [ME/+5; E+2 (5.22) (5.23) (5.24) 623 (5.26) (5.27) To implementing the boundary conditions, the procedure discussed in chapter four should be followed. Therefore, the unconstrained momentum equations can be 105 written as ‘ f, [2.73—3.22 +72%; + $2,232] an +§ fn(u — g) w1|l6w|| d9 — (917%‘19 = 0 v (5.25) (5.29) —fnp‘%";1d§2 = 0 where “662“ is the boundary measure function which is discussed in appendix (B). Substituting equation (5.15) into equations (528,5.29) yields the final momentum equations written in a penalty formulation as 611611) 6n 6v 6w fa [27757—521 + "(5; + a—J—L] d9 1 5' MIL - glwlllawll d9 61) 39 l 6' 6w 331d!) ) fat‘S-l‘; + 106 (5.30) and the y-momentum is 5.4 8v 8:: 6w 3 )79‘22 +2776: 9—5092] an (5.31) Wavelet-Galerkin Method and Discretization of the Momentum Equations In the wavelet-Galerkin method, the dependent variables u, v are expanded in terms of scaling functions of Daubechies type as follows: ”(56, y) = 22(1), y) = where X = 2%, Y :2 2jy, N = N—l N—l Z Z Umx — i)¢(Y — k) (5.32) i=0 k=0 N—l N—l mm — i)¢(Y — k) (533) ll O a. ll 0 H. 23', Ujk = 23m,“ ij = 23.12,)c , and j is the level of discretization. In the Galerkin formulation the test functions 10's are chosen from the same space as u, v, i.e. ml 2 wz = ¢(X — p)¢(Y — q) where p and q are integer that are the translates of the scaling functions ¢(X) and ¢(Y) respectively. The boundary 107 conditions in equation (5.14) is also expanded in wavelet space as N— lN-l 9(a: 3): Z Zen-tam ¢(-Y k) (5.34) i=0 1:: 0 where Gm = 2j9pq : //g(x,y)¢(X — P)¢(Y " q)d:1:dy (5.35) The coefficients U’ s and V’s are found by substituting equations (532,533,534) into equations (530,531). This will result in having the momentum equations in scaling function expansions. For example, the x-momentum equation is discretized 108 as follows: :01: —5kq 2nz..U.1/¢'(X— )'—¢(X p)dxf¢(Y— k)¢(Y— q)dy =‘fip =Cfilk +77 23,1. ch / ¢(X - z')<(>(X -p)d$/ ¢’(Y - k)¢'(Y - aldy :Cvloi 0:211 mam/('(X —z') X:)dx/¢(Y- —q)dy +- 22.1%. Haw“ / ¢(X «was p)dx / ¢(Y —k)¢(Y — q)dy = :5,” —— :2..sz Haw” / ¢(X — i)¢(X p)dx/ ¢(Y )(Y— q)dy :01: i —% E... U... f ('(X — z')¢'(X —- pm] ¢(Y — k)¢(Y — my =(j'10p =C30k J fl -%§:.—,mk/¢(X—z‘)¢'(X —p)dxf¢(—Y k)¢(Y m) =0 > l (5.36) Taking advantage of the delta function, equation (5.36) is written in a simpler manner as [a II H H pH 4):":4 V'kczigicligq + I: qu "' a: Gm l i=P— k— q-4 ’ _l n+4 11 1 p+4 0+4 10 10 _ ' UiQpC — ‘ 4Zk:q—4 ikCi— qu—k —0 5 z=p—-4 —p_ J 109 (5.37) Similarly the discretized y-momentum equation is 4 4 4 ”Zip—4 22:1,- 4 UikC.-1°p01°k+2nZZ:,_4 Vkaqlk +77 ZpH—V Viquli + “6%)“le ’ @11qu ( (5'38) “221251—41 2123-114 Uikczio £1230 q_El-'ZZ::— 4 Vkauk— — 5.5 Solution of the Resulting System of Equations Equations (537,538) form a nonlinear coupled system of equations in terms of the coefficients U’ s and V’ s. This system of equation is written in matrix form as AU 0 U Fu = (5.39) 0 By V Fv This system can be solved iteratively using Newton’s method, or more easily by using the Picard iteration method. The system of equation (5.39) can be converted to an easier one to solve, by segregating the dependent variables U and V from equations (537,538) [33] as p+4 q+4 llawllU (277 — -) Z UiqC“ + 7] Z Ukaflk + —'€—'—U pq — '"uF (5.40) i=P— —4 k=q— —4 110 where p+4 9+4 p+4 q+4 “at—Ulla = -7, Z Z V-kCIOC _, +— :Z Z V; 0.10 1,010,, + ——€—G ,, (5.41) i=p—4k:q i=p— 4k:q— —4 The dependent variables U and V in the y-momentum equation (5.38) can be segre- gated as p+4 q+4 ”aw €l___l 1? Z 14,011+ (247— -) Z V, 011—14,, = F, (5.42) i=p—4 k=q-4 where p+4 (1+4 p+4 q+4 ”aw” _ ——-n E E U 010 ,Cl°,+: Z 2 U C,§1,-C,1‘1,+ —G,, (5.43) i=p— 4k=q— 4 Ei=p— 4k=q- -4 8 This segregation scheme produces equations (540,542) which are series of non- linear Poisson equations and can be solved iteratively by Picard Iteration method. 5.5.1 Picard Iteration Method The iteration procedure of this method is simple. For example, equation (5.40) can be written in matrix form as AWU=RW) W“) B (77W = F507) (5-45) 111 It should be noted that n is a function of U and V. For the system in (5.44) the Picard algorithm is given by 44(77"‘)U""+1 = F1107“) (5-46) where the superscript indicates the iteration number. This method is known to hav a reasonably large radius of convergence. An improvement in the convergence rate can sometimes be obtained by use of a relaxation formula [33], A(n’")U* = F501“) (5-47) where U'"+1 = aU’" + (1 — a)U* (5.48) The matrices A(17m) and B(77"‘) are not circulant, nearly symmetric. 5.5.2 Calculation of the Viscosity Function 17’” The viscosity function 71’” is calculated using a scaling function expansion. Equation (5.13) can be written as "1-22 2(91‘1)2+2(Qfl)2+(3_u’1 2,. 21212:)..(232) 7;" 77 —770 82: 8y By By 826 8:13 (5.49) 112 The following substitution can be used to simplify the evaluation of 17: , ni"=a§;" 775":ng né"=%3 772"=65’:) (5.50) The following wavelet expansions are used for 771", 17;", 17g",and 17;” at each iteration m: ‘ 77in = Eli/filo 61 "L‘MX - 0450/ “ k) 77$" = 2.17.2052, (X — z‘)¢(Y— k) 115" = 2.”, ‘0 ea",,.¢(X — z‘)¢(Y — k) ”in = Zia-=0 e4 mk¢(X — i)¢(Y _ k) J (5.51) The coefficients e’s are found by substituting equations (551,532,533) into equations (5.50) and taking the inner product of both sides by ¢(X — p)¢(Y — the evaluation of 31.4 is as follows: 6,, 5kg 2.2,. [3(x_.)¢(x_ p)dz/¢(Y— kw —4)dy= ___Cl0i 6kg A L 24,1530 .-’,’i./¢’(X-i)¢ (-X p)d$/¢(Y- ”MY-(1)61; 113 q). For example, (5.52) Therefore Similarly and N—l 6;" = E UimC'10 Pq vq p_i i=0 N—l m _ m 10 62139 — Up,qu—k k=0 N—l m _ m 10 83W — ‘fp,qu_k k=O N—l 8;" :2 :VmCIO m m p—i i=0 (5.53) (5.54) (5.55) (5.56) With the knowledge of the coefficients e’s at each iteration, we substitute those e’s into equations (5.51) to evaluate nf‘,n§",n§",and 17;" at each iteration step m. The nfn’s are then substituted into equation (5.49) for the calculation of 17'" at each point (p, 9)- 5.6 Examples The examples that were considered in Chapter Four will be considered again here, namely the flow between concentric cylinders and lid-driven cavity problems. Those 114 examples will be computed for different values of the dimensionless exponent n. These examples are as follows. 5.6.1 Flow Between Rotating Concentric Cylinders In this example, the configuration of this problem in the fictitious domain method with the boundary conditions is shown in Figures (4.1,4.10). This problem has an exact solution with the following assumptions: o Steady state fluid flow 0 Incompressible flow 0 Fully developed flow. The exact solution for this problem is as follows: W =3 011.3 + 027‘ (5.57) n — 3 s — n +1 (5.58) V ,-— V. o 01 = ”J ”‘T (5.59) " (5.60) Leve|=5 Leve|=6 1 r 1 . ""'""ZZ\\\ //::::::::~~\\\\\ {II/Irv ''''' ‘\\\\\\\ /////”' ..... ~~\\\\\ \ {III/II! ........ \\\\\:\\ 0.5 , {III/I” .......... “‘\\\\ ‘ 0.5 [Illi' \\\ ‘\‘ [551/11 \\\ ‘ IIIII’ u\\:‘ I’ll' o‘\“ I ’1‘. \\‘ ‘ ’5'. .o:' {‘0' .u j T 0’ "" iii. ‘ >~ o t::: o.’;$ .. DU, \ \‘:“ III’II §§\\' oII’,’ \§\\\\“' oli”,$;, \\\\\\\Q.' ..aol”;,, -0.5 \\\\\\\s ~~~~~~~~~~ ”/;’/ -0'5 \\\\\\‘5 ....... '1’), \\\\\\\s nnnnnn aI’IIII \\\\\\‘ ...... a’l/z/ \\\\\s-.—.—o—~a’/’// \‘WM” ‘~“-J’ '1 A C ‘ —1 . . m X X Figure 5.1: Velocity vectors for the flow between two concentric cylinder problem (71. = 0.8). where V290 is the outer tangential velocity, V3, is the inner tangential velocity, 7‘, is the outer radius, and r,- is the inner radius. Figure (5.1) shows the velocity vectors of the flow field. Figure (5.2) shows a comparison between the wavelet solution and the exact solution for this example for different levels. 116 --- Exact Solution — Wavelet Solution Level=5 Leve|=6 1 . 1 . 0.8 0.8 ° :1) 3 g "g 0.6 * g 0.6 '5" an 5 2 a? 3? .3 0.4 § 0.4 0 a > > 0.2 0.2 n O 1 -1 0 _1 0 Figure 5.2: Comparison between the wavelet and exact solutions at the crossection y = 0, ((9 = 0, 7r in polar sense)(n = 0.8). 117 1.2 1 1 I l I 1 1 - W — ””fla—i—z—Wx‘ I f f / / l I ..o .—-' _. _. .., ..,, ‘ x \ \ \ \ ‘ ‘ r f f I 1' I I ; .. _. _. .. .. \ \ \ \ \ l ‘ ‘ 0.8 .., \ 1 1 ‘ 1 i t o ‘ . ‘ ‘ I I l I l .. \ \ \ \ \ \ \ \ s .. - , , ’ I I I I I ’ I \ \ \ \ \ \ \ \ x .. .. , I I I / l I I I , \ \ \ \ \ \ \ \ \ .. .. , I I / I I I I ’ , 0.6 _ 5 \ \ \ \ \ \ \ \ .. .. ,. I I I / I I I , . - \ \ \ \ \ \ \ c. .. .., .— r I I I I I I , 0 \ \ \ \ \ x s .. ._ .. I I I I I , l a . \ \ \ \ x \ s .. .. .. ,. , , , , , , 0 O4 '- ‘ ‘ ‘ ‘ ‘ ‘ * - '- 0’ v r I I I I a '4 0.2 — - - - - - - — - — - O .- — _0.2 I 1 M 1 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Figure 5.3: Velocity vectors for the lid-driven cavity problem (level=5). 5.6.2 Lid-Driven Cavity Problem This problem is solved using the boundary conditions shown in Figure (4.2). Re- sults are provided for this problems by velocity vectors and stream lines for different values of the dimensionless exponent 72. Figure (5.3) shows the velocity vectors, Fig- ures (5.4) shows the stream lines for this problem. It should be pointed out that when the dimensionless exponent n equals 1 the Newtonian flow is recovered, see Figure (5.4). The iteration history of this problem is shown in Figure (5.5). It should be pointed out the iteration was stopped when the residual satisfied the required 118 0.2 0.4 0.6 0.8 X Figure 5.4: Streamlines for the lid-driven cavity problem (level=6), left Newtonian fluid, and right non-Newtonian fluid. tolerance. The residual is defined as Residual = (/Z(vn+1 — w)? (5.61) 5.7 Discussion and Summary In this chapter, a wavelet-Galerkin method was presented to solve some non- Newtonian fluid flow problems. This method is developed with a penalty formulation of the relevant functional. A segregated solver was used to de~couple the governing equations. This lead to the conversion of complicated non-linear momentum equa- tions into a series of uncoupled non-linear Poisson equations. Those Poisson equations were solved iteratively by the Picard iteration method. The method presented in this 119 Residual 100 90 80 70 60 50 40 30 20 10 l l l l l 20 40 60 80 100 Iterations Figure 5.5: Residual history of the Picard iteration method. 120 120 chapter requires the solution of two Poisson equations at each iteration in the Picard method. This procedure was tested by various examples. Results of these examples suggest that the proposed solution method is accurate and rapid. 121 Chapter 6 Conclusions 6.1 Summary and Contributions Wavelet-based techniques were used in this work to solve heat transfer, Newtonian and non-Newtonian fluid flow problems. Inherent to these techniques are the diffi- culties encountered in implementing the boundary conditions and two methods to enforce them were discussed. Clarifications were made on the techniques used to implement the boundary conditions. Also, the solution of Newtonian and nonlinear non-Newtonian fluid mechanics problems with non-periodic boundary conditions is achieved efficiently and a variety of iterative methods are used to solve these prob- lems. A segregated formulation based on the solution of a Poisson equation allows to solve these complex problems efficiently. All this is done using a modest desktop computer. 122 6.2 Findings of the Study It is found in this study that both the Fictitious Boundary and the Fictitious Do- main / Penalty methods provide good approximations to the solution of the differen- tial equations, but that the Fictitious Domain/Penalty Formulation method exhibits slightly less errors than the Fictitious Boundary Approach. Also, the errors found in the examples, which were solved by the Fictitious Domain/Penalty Formulation method, can be further reduced by using a smaller penalty parameter 5. Furthermore, the resulting stiffness matrix with the Fictitious Domain/ Penalty Formulation has a very nice properties that allow its rapid inversion. In the Fictitious Domain/ Penalty formulation, the resulting system of equation is solved iteratively using the Conjugate Gradient Method and the Preconditioned Con- jugate Gradient Method. A simple diagonal preconditioner matrix is suggested in this study. Using the Preconditioned Conjugate Gradient Method reduced the necessary numbers of iterations to almost half of the standard Conjugate Gradient Method. The Fictitious Domain/ Penalty formulation can be applied within the framework of multiresolution analysis in solving partial differential equations. The use of the Fictitious Domain/ Penalty formulation and Preconditioned Conju- gate Gradient method, allows achieving an efficient and quick solutions to Newtonian and non-Newtonian fluid flow problem using the capability of a modest personal computers (Pentium Pro 200 MHz with 64 Mb of RAM). In addition, in the non-Newtonian flow problems, a penalty method was used to implement both the boundary conditions and enforce continuity. A segregated 123 scheme was used in the non-Newtonian fluid flow problem and this turned compli- cated, non-linear momentum equations into a series of Poisson equation which are solved iteratively and efficiently using the wavelet-Galerkin method. The wavelet formulation of the problems in this work allow modeling and solving large problems without the need of generating a complex geometry. This was demonstrated with a few examples. 6.3 Future Research Topics The results of this dissertation are only a small step in the evolution of wavelets and their use to solve solution of partial differential equation related to viscoelastic fluid flow problems. To continue on this work, the following topics should be investigated. 0 Efficient Preconditioning Finding more efficient preconditioner than the simple diagonal preconditioner used her is an essential step toward the solution of very large problems. 0 Application of wavelet-based methods to the solution of viscoelastic fluid flow problems Applying the current method to solve standard viscoelastic flow problems and study the performance of wavelets as compared to other techniques. 0 Extension of the formulation of fluid flow problems to problem with high Reynolds Number. 0 Use of divergence free wavelets for solving incompressible fluid flow problems 124 This can be exploited by choosing the proper bases functions for velocity and pressure fields to satisfy the Ladysenskaja-Babuska-Brezzi (LBB) condition [8, 11]. o Multiresolution wavelet method Applying the Mallat algorithm [30] to solve some heat transfer/ fluid flow prob- lems using the multiscale analysis. 0 Implementation of other types of boundary conditions, e.g., Neumann and Robin by using the appropriate modified functionals. 125 APPENDICES 126 Appendix A Construction of Generalized Wavelet Systems The Haar scaling function discussed earlier is a special case of a more general class of functions. In general a scaling function is a solution to a dilation equation of the form 4(4) = Z 414(24 — 4) (A.1) k A wavelet function is defined by 4(4) = Z(—1)*a~-1_.¢(2z - k) (A2) I: where N is an even positive integer, and the coefficients a), are called the filter coef- ficients. Normally, there is a finite number of these coefficients which are non-zero. 127 AA rl-‘l-xll'l The construction of a wavelet function is based on the following steps: 1. finding the filter coefficients, 2. finding the scaling function. In order to derive the filter coefficients certain conditions on the scaling and the wavelet functions have to be imposed: 1: r-T‘ 1. the scaling function and its translates should be orthonormal, : 2. the wavelet function should be orthogonal to the scaling function, 3. the scaling function is required to be able to exactly represent polynomial of a certain order, say 12. A.1 Construction of The Filter Coefficients Consider the normalization of the scaling function (6, such that +00 ¢(x)d:r = 1 (A.3) —00 this yields the following condition for the filter coefficients ak: Zak = 2 (A.4) k 128 For the scaling function to be orthogonal to its integer translates, the wavelet coeffi- cients must satisfy the additional requirement that +00 (:0) (15(1: — l)d:r = 501 (A.5) —CD this yields the condition 2 a), ak+21 = 2601 (A.6) k where 6 is the Kronecker delta function. The scaling function and its wavelet function have compact support when a finite number of wavelet coefficients are non zero. Equations (A.4,A.6) will produce, in an N coefficients system, 5.1,,- +1 equations. This is not enough to determine a unique set of filter coefficients. Daubechies [12, 13] enforced the requirement that the scaling functions are able to exactly represent polynomials of order up to, but not greater than, p. For a polynomial of the form f(fL‘) = Co +Cl$+02$2 + +Cp_1IL'p-l (A.7) can be exactly represented by an expansion of the form 1(4) = Z 4.40: _' k) (A8) 129 Taking the inner product of equation (A8) with the wavelet function 111(23) gives (due to the orthogonality condition) (f($).10($))= 20144505 - 1413(3)): 0 (A-9) Thus, from equation (A.7), co/¢(2:)dx+61/1/2(Lr)xdx+02/7,b(x)x2dx+---+cp_1/1/1(:c)xp-l deO (A.10) This is valid for all c,-(z' = O, 1, ..... , p — 1)[12, 13, 38]. By choosing c, = 1 and all other c,- = 0, gives the moment equation +00 (1(3) 45' d1: = 0, l: 0,1,2,...,p — 1 (A.11) —00 It was shown in [13] that for N coefficients system, p = %. The final conditions to get a unique set of coefficients {ah I: = 1, 2, ..., N — 1} comes form the substitution of equation (A.2) into equation (A.11) N E :(-1)kak k’ = 0, 1: 0, 1,2, —2— — 1 (A.12) k 130 "..— A.2 Construction of the Scaling Function The scaling functions, in general, do not have closed form solutions. They are gener- ated recursively from the dilation equation (A.1). Without going into a complicated mathematical analysis, we can write equation (A.1) as 06(22) = ao¢(2$) + a1¢(2:c — 1) + a2¢(2a: — 2) + - . - + aN_1¢(2x -— (N — 1)) (A.13) Writing this equation for all integer values of 25,-(i = 0, 1,2, ....,N — 1), produce the following system of equations: 0005(0) M” = 00¢(2) + 01¢(1) + 0245(0) S A O V II 49(2) = 0095(4) + 0145(3) + 02¢(2) + 0345(1) '1' 04¢(0) (A.14) ¢(N — 2) = aN_3¢>(N — 1) + aN-2¢(N — 2) + GN_1¢(N — 3) (MN-1) = aN—1¢(N— 1) This can be written in a matrix form as (M — I) = 0 (A15) 131 where I is the identity matrix. The vector (1) and the matrix M are as follows: (15(0) 45(1) (1): 9(2) ¢(N-2) _¢(N-1>. F . a0 0 0 0 0 0 a2 a1 a0 0 O 0 a3 a2 a1 O O 0 M : 0 0 0 (IN—3 aN—4 (IN—5 0 0 0 951—1 931—2 @1141 L 0 O 0 0 0 (110-1 . The solution to equation (A.15) is not unique, and so a normalization condition is required in order to uniquely determine the vector . That is 2‘1“") =1. k E Z (A.16) ’6 Hence, the values of the scaling functions at the integers are given by the solution of 132 equation (A.15) normalized by equation (A.16). Now the values of 05(3) are known at the integer points of x, the values at the half integers can be found from equation (A.1) 44;) = Z 414(24 — (4) (AI?) 1: This process is repeated as many times as necessary to find the values of 05(10) at all dyadic points {2,L; i,j E Z}. A.3 Wavelet Multiscale Representation of Func- tions Multiresolution analysis consists of deve10ping a representation of a function f (2:) at different level of resolution. This is achieved by expanding the given function in terms of a function 05(25) called the basis function. The basis function can be scaled to give multiple representation of the original function [30, 38, 44]. This require the use of a scaling function 051(2) which is given by 4(4) = 2441(4) (A.12)) 1: 41(4) = 4(21'4 — 4) (1)19) 133 where k is an integer called the translate of ¢(:c), j is called the level of repre— sentation, and c; are called the filter coefficients. For example, this scaling function can be defined as in [24]: set all c}, on level 3' equal to zero except for 0,7 equal to one. Using the interpolating subdivision scheme results in ¢{(:1:). Also, using linear superposition and the subdivision scheme at level j, yield that [38] ((4) = 24.4034 — k) = 2441(4) (MO) I: I: It can be seen from equation (A.19) how all scaling functions are simply translates and dilates of one fixed function. In order to develop a multiresolution representation of a function in L2(R), the square of integrable functions, one seeks a sequence of embedded subspaces Vj. Define the space V], such that V, =span { 4;;(4) |k e Z } (A.21) This implies that 3.115;: P1f(ar) = f (x). where P,- is the L2 orthogonal projection from L2(R) to V,- [31, 38]. For example, if (0(5):) is the Haar scaling function, the V} is the space of functions I: 1 which are piecewise constant on the interval [55, —'—2t-) with k E Z. The different VJ- spaces satisfy the following prOperties which make them a Multiresolution analysis 134 [38], i.e., wavelets space satisfy the following conditions: 1. Nestedness: V,- C V,'+1. 2. Translation: ifg(2j:c) E V,- then g(2j:1: + k) E V}. 3. Dilation: the subspaces are related by scaling law: g(2j:r) E V,- then g(2j+l:2:) E 14..., 4. Completeness: every function g(:1:) in L2(R) can be approximated with arbitrary precision with a function from V,- for suitably high j. Sometimes this property is also expressed as: UjezV} is dense in L2(R). The translation and dilation properties follow from the definition of V,- equa- tion (A.21).The nestedness property follows immediately from the dilation equation (A.18). If 05(93) can be written as a linear combination of the ¢(2x — k), then 051(3) also can be written as a linear combination of W; +1(gr) and thus V,- C V,-+1. In the literature , this often refers to a scaling function W (:17) E V,- such that its integer translates [ W(a: -— k), k 6 Z} is from a Riesz basis for the space V,. Using the Haar scaling function Figure (A.1), the demonstration of the multires- olution of functions is presented next using the Haar scaling function. The Haar scaling function ¢(.2:)is defined as 1 0 _<_ :1: g 1 05(4) = (A.22) 0 otherwise 135 1 4 a 1 5 1.2> 1 i 1 0.5 10.8 < ._ I X v v O» ° 0.6 < 5 .05) 0.4» 1 0.2- 4 '1 o 1 -15 - - o 0.5 1 1.5 o 0.5 1 1.5 X K Figure A.1: Haar Scaling [and Wavelet Phnctions f (:r) E L2(R) can be approximated by its projection onto the space V0 such that 1%) = Z 424(2°4 — 4) (A23) 1: or in general a function may be approximated by its projection onto the space V,- as f’(:v) = Z 414(214 — k) (4.24) I: Figure (A.2) shows different approximations to a function at different levels. The wavelet function, on the other hand, may also be used to define orthogonal subspaces [38] such that let W,- = Span{¢i(z‘) : k E Z}, then Ujez W,- = L2(R) (A.25) where 136 level=5 Original Function 1(x) level=6 '0‘) f(x) X Figure A.2: Multilevel representation of the function f (11:) = sin(27r:z:) v... = Y.- 4 W.. V.- 1 W.- (426) where 69 represent a direct sum. This implies the Mallat transformation [44] V0 C V1 C V2 ------ C V141 (A27) and V}+1 = V0 613 W0 69 W1 """ G9 Wj (A428) 137 From the relations (A26) and (A.27) above, one can notice that wavelets provide bases for W,- to describe the difference between to successive subspaces V,- and V,-+1. We then introduce a wavelet function (0(22) such that (Ma: — It) forms a Riesz basis for the subspace W0: 4),: = 244(22): — k) (A29) is a Riesz basis for W,-, which leads to the relation (A.25). It should be pointed out that since W0 is contained in the space V1, we can write the wavelet function in terms of the scaling function at the next higher level, 40(9) = Zbk¢(21$ — 14) (A.30) k The Haar wavelet functions shown in Figure (A.1) satisfy equation (A30) with coefficients ()0 == 1, and bl = —1. So if ¢(:1:) is the Haar scaling function, the corre- sponding wavelet function 111(33) is 10(33) = 4(255) - 45(217 - 1) (A-31) 138 Appendix B Incorporation of the boundary measure function [[802] | The boundary measure function ”80)” is defined in [43] as llawll = -Vx.. - n (8.1) 1 if (2:, y) on 00) 0 otherwise where n is the outward normal to the boundary, and wa is the gradient of a char- acteristic function xw. Typically, [[602]] has zero values in all points of the domain {2 except for the points near to the boundary (902. The boundary of (12, which is an open 139 bounded set in R2, is defined as aw = {(x, y) e R2 : F(:I:.y) = 7} for some constant 7 and for some Lipschitz function F defined in R2 of 6w, with (12 having a finite perimeter. The boundary measure function [[601]] can be approximated by wavelet expansion Haw“ = -Vx.. - n = pr,q¢(2jx - 1045(ij - q) (B3) =Zq:XP,($-q¢ My“ (1) (13.4) F; F144 9(12 -(p)¢ 9— 9) (3.5) where XM = éxwfifl 2,- , 2,- Sii), and FM = -21—,4F(321—;3,%9), with c = fx¢(z)dx. The expansions for the partial derivative of X4.) and F are as follows: “69;; y) = 2%)... 4(4 — (4)40) — 4). (8-6) P14 140 W = 2(9’5).,.4(4 — 44(4) - 4), (8.7) By M 33/ 811;: y) = ;(:—:)p.9 (“1' ‘ PM“?! _ (1)) (BB) ——8F;:’ y) = §(g—§),,q 9(27 - p)¢(y - 9), (39) where, for example, the coefficients gxxp q can be calculated as 0X (a—xhw : Eclianmfl- (13.10) Therefore 8 61" 5 BF ( .5)...(-..;)... + (55)...(5,;).,.) ((fi(%§)...)2 + ((35)...)2) For implementation purposes, the approach suggested in [15] in dealing with (3.11) (“10.9 = “ ”800“ is used. In [15] [[602]] is set to be equal to constant 7 in the points near to the boundary ('90) and zero elsewhere. The value of 7 is chosen to make the integral f5, ”60)” d9 equal to perimeter of 0). For the two-dimensional problem Haw“ is 7 if [n,n+1]x[m,m+1]flaw#0 ”840”,... c: (8.12) 0 elsewhere 141 Appendix C Algorithms o Conjugate Gradient Compute r0 = b — Axe, p0 = 20 For j = 0,1, ..... , until convergence Do: a. = (Tz',1‘z') J (A pjvpj) MH=%+%M 71+): 7'2“ — 02' P1 ’8. = (Tl-+11 rid-1) J (’1': '1') Pj+1 = r,-+1 + 153' Pj EndDo o Preconditioned Conjugate Gradient Compute To = b — A130, 20 = M-lTo, p0 = 20 For j = 0,1, ..... , until convergence Do: 142 a_ : (T,‘,Z,‘) ‘7 (A pjapj) au=a+wm 73+) = 7‘1 - 09' P1 _ -l Zj+1 — M 7‘j+1 16' : (Tz'+l: 21H) J (71.2)) Pj+1 = Zj+1 + 5]" Pj EndDo where A is the stiffness matrix, M is the preconditioner matrix, :1: is the vector of the unkowns, r is the residual vector, p and z are intermediate vectors defined in the algorithms. The coefficients a and 3 are scalar quantities. 143 BIBLIOGRAPHY 144 Bibliography [1] K. Amaratunga and J. William. Wavelet based green’s function approach to 2d pde’s. Engineering Computations, 10:349—367, 1993. [2] K. Amaratunga, J. Williams, S. Qian, and J. Weiss. Wavelet-galerkin solutions for one-dimensional partial differential equations- International journal for nu- merical methods in engineering, 37:2703—2716, 1994. [3] O. Axelson. Iterative Solution Methods. Cambridge University Press, Cam- bridge, 1996. [4] E. B. Becker, F. G. Carey, and J. T. Oden. Finite Elements. Englewood Cliffs, NJ : Prentice-Hall, 1981. [5] S. Bertoluzza and G. Naldi. A wavelet collocation method for the numerical solution of partial differential equations. Applied and Computational Harmonic Analysis, (3):1—9, 1996. [6] R. B. Bird, R. C. Armstrong, and O. Hassager. Dynamics 0f Polymeric Liquids. John Wiley and Sons, Inc., New York, 1987. [7] J. H. Bramble and J. E. Pasciak. A preconditioning technique for indefinite systems resulting from mixed approximations for elliptic problems. Math. Comp, (50):1—17, 1988. [8] S. Dahlke and I. Weinreich. Wavelet-galerkin-method: an adapted biorthogonal wavelet basis. Technical report, No. A-91-25, Freie University, Berlin, 1991. [9] W. Dahmen, A. Kurdila, and P. Oswald. Multiscale Wavelet Methods For Partial Difierential Equations. Academic Press, San Diego, 1997. [10] W. Dahmen and C. A. Micchelli. Using the refinement equation for evaluating integrals of wavelets. SIAM J. Numer. Anal, 30:507—537, 1993. [11] W. Dahmen and R. Schneider. Composite wavelet bases for Operator equation. Preprint TU Chemnitz, 1996. [12] I. Daubechies. Orthonormal bases of compactly supported wavelets. Commun. Pure Appl. Math., (41):909—996, 1988. 145 [13] I. Daubechies. Ten Lecture on Wavelets. CBMS-NSF Regional Conference Series, SIAM, Philadelphia, 1992. [14] William M. Deen. Analysis of Transport Phenomena. Oxford University Press, Inc., New York, 1998. [15] Alejandro Diaz. A wavelet-galerkin scheme for analysis of large scale problems on simple domains. International journal for numerical methods in engineering, 44(11):1599—1616, 1999. [16] DeRose Jr. Giuseppe C. A. Solving Topology Optimization Problems Using Wavelet-Galerkin Techniques. PhD thesis, Michigan State University, 1998. [17] R. Glowinski. Numerical Methods for Nonlinear Variational Problems. Springer- Verlag, Inc., New York, 1984. [18] R. Glowinski, T. W. Pan, and J. Periaux. A fictitious domain method for dirich- let problem and applications. Computer Methods in Applied Mechanics and Engineering, (111)2283-303, 1994. [19] R. Glowinski, T. W. Pan, and J. Periaux. A fictitious domain method for external incompressible viscous flow modeled by navier-stokes equations. Comp. Meth. in Appl. Mech. and Eng, (112)1133-148, 1994. [20] R. Glowinski, T. W. Pan, R. O. Wells, and X Zhou. Wavelet and finite element solutions for the neumann problem using fictions domains. Technical Report, Computational Mathematics Laboratory, Rice University, 1994. [21] R. Glowinski and O. Pironneau. Energy Methods in Finite Element. Wiley- Interscience Pub., 1979. [22] R. Glowinski, A. Rieder, R. O. Wells, and X Zhou. A wavelet multigrid pre- conditioner for dirichlet boundary-value problems in general domains. Technical Report, Computational Mathematics Laboratory, Rice University, 1994. [23] L. Jameson. On the wavelet optimizes finite difference method. Technical report, Institute for Computing Applications in Science and Engineering, 1994. [24] B. Jawerth and W. Sweldens. An overview of wavelet based multiresolution analysis. Preprint, World Wide Web, 1996. [25] A. Kunoth. Multilevel preconditioning—appending boundary conditions by la- grange multipliers. Advances in Computational Mathematics, Special Issue in Multiscale Methods, 4:145—170, 1995. [26] A. Latto, H. L. Resnikoff, and E. Tenenbaum. The evaluation of connection coefficients of compactly supported wavelets. Aware Technical Report AD910708 (July 1991) and the Proceedings of the French- USA Workshop on Wavelets and Turbulence, Princeton University, June 1991, Editor, Y. Maday, Springer- Verlag, 1992. 146 ' [27] D. Lu, T. Ohyoshi, and K. Miura. Treatment of boundary conditions in one- dimensional wavelet-galerkin method. JSME International Journal, Ser. A, 40(1):382-388, 1997. [28] D. Lu, T. Ohyoshi, and L. Zhu. Treatment of boundary conditions in the ap- plication of wavelet-galerkin method to sh wave problem. Int. J. of The Soc. of Mat. Eng. for Resources, 5(1):15—25, 1997. [29] Dianfeng Lu. Study of Boundary Value Problems in Application of Wavelet- Calerkin Method to Wave Motion Analysis. PhD thesis, Akita University, Japan, 1998. [30] S. G. Mallat. Multiresolution approximations and wavelet orthonormal bases of L2(R). Trans. Amer. Math. Soc., 315(1):69—87, 1989. [31] S. Qian and J. Weiss. Wavelets and the numerical solution of partial differential equations. J. of Computational Physics, (106):155—175, 1993. [32] J. N. Reddy. An Introduction To The Finite Element Method. McGraw-Hill, Inc., New York, 1993. [33] J. N. Reddy and D. K. Gratling. The Finite Element Method in Heat Transfer and Fluid Dynamics. CRC Press, Inc., Florida, 1994. [34] J. M. Restrepo and G. K. Leaf. Wavelets-galerkin discretization of hyperbolic equations. Journal of Computational Physics, (122):118-128, 1995. [35] J. M. Restrepo and G. K. Leaf. Inner product computations using periodized daubechies wavelets, 1997. [36] Yousef Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, 1996. [37] A. Sowayan, A. Bénard, and A. Diaz. Wavelets and the numerical solution of heat transfer problems: A discussion of two methods. To be appear. [38] W. Sweldens. The Construction and Application of Wavelets in Numerical Anal- ysis. PhD thesis, University of South Carolina, 1995. [39] A. Tabacco, C. Canuto, and K. Urban. The wavelet element method, part i: construction and analysis. Preprint No. 1038, Instituto del Analisi Numerica del C.N.R Pavia, to appear in ACHA, 1997. [40] A. Tabacco, C. Canuto, and K. Urban. The wavelet element method, part ii: realization and additional features in 2d and 3d. Preprint N o. 1052, Instituto del Analisi Numerica del C.N.R Pavia, to appear in ACHA, 1997. [41] F. Thomasset. Implementation of Finite Element to Navier-Stokes Equations. Springer-Verlag, Inc., New York, 1981. 147 I ll [42] J. Weiss. Wavelet and the study of two dimensional turbulence. Aware Technical Report AD910’708 ( July 1991) and the Proceedings of the French- USA Workshop on Wavelets and Turbulence, Princeton University, June 1991, Editor, Y. Ma- day, Springer- Verlag, 1992. [43] R. O. Wells and X Zhou. Wavelet solution for the dirichlet problem. Technical Report 90-92, Computational Mathematics Laboratory, Rice University, 1992. [44] J .R. Williams and K. Amaratunga. Introduction to wavelet in engineering. In- ternational journal for numerical methods in engineering, 37:2365—2388, 1994. 148