A P P L I C A T I O N S O F A H Y B R I D A S Y M P T O T I C - N U M E R I C A L M E T H O D F O R C E R T A I N T W O - D I M E N S I O N A L S I N G U L A R P E R T U R B A T I O N P R O B L E M S by M I C H E L E S U S A N N E T I T C O M B E B.Math. (Mathematics) University of Waterloo, 1988 M.A.Sc . (Aerospace Engineering) University of Toronto, 1993 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A Apri l 1999 © Michele Susanne Titcombe, 1999 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract This thesis examines a hybrid asymptotic-numerical method for treating two-dimensional sin-gular perturbation problems whose asymptotic solution involves reciprocal logarithms of the small perturbation parameter, e, in the form The purpose of this hybrid asymptotic-numerical method is to treat the slow convergence prob-lems of asymptotic expansions of this form. For the applications that we consider in this thesis, we believe there is sufficient evidence of convergence of these expansions for small enough e. The hybrid method uses the method of matched asymptotic expansions to exploit the asymp-totic structure to reduce the problem to one that is asymptotically related to the original. In general, one must solve this related problem numerically. The hybrid related problem contains the entire infinite logarithmic expansion in its solution, thus removing the necessity of obtaining each coefficient in successive terms individually, as one would have to do using only the method of matched asymptotic expansions. The hybrid solution essentially sums the infinite expansion of reciprocal logarithms and in so doing, improves the accuracy of the solution since the error of the approximation is smaller than any power of ( —1/loge). A n important feature of the hybrid related problem is that it is non-stiff. Thus, it does not suffer from the difficulty in solving the original problem numerically of resolving the rapidly varying scale structure. Another advantage of the hybrid method solution is that the parameter dependence of the problem is reduced from that of the original. The reduction in parameter dependence means that the hybrid method solution is less computationally intensive than a full numerical solution. We show that singular perturbation problems containing infinite logarithmic expansions arise in a wide variety of contexts. Four chapters of this thesis are dedicated to the detailed application of the hybrid method to such singular perturbation problems occurring in fluid flow in a straight pipe with a core, skeletal tissue oxygenation from capillary systems, heat transfer convected from small cylindrical objects, and low Reynolds number fluid flow past a cylinder that is asymmetric to the uniform free-stream. Following the detailed analysis of these four problems, we remark on possible extensions to the general framework of applicable problems. For example, we discuss applications in black body radiation, multi-body low Reynolds number fluid flow, vibration of thin plates with small holes or concentrated masses, localized non-linear reactions on catalytic surfaces, low frequency scattering of light, diffusion of a chemical species out of an almost impervious container, and steady-state current flow from microelectrodes. 11 Table of Contents Abstract ii Table of Contents iii Acknowledgement v Chapter 1. Introduction and Background 1 Chapter 2. The Hybrid Asymptotic-Numerical Method 7 2.1 Conditions on Applicable Singular Perturbation Problems 7 2.2 Framework for Applicable Second-Order Steady Problems 8 2.3 Outline of the Hybrid Method on a Second-order Linear Problem 10 2.4 Effect on the Asymptotic Solution Structure of Relaxing Certain Conditions 17 2.5 Subdomain Shape-Dependent Parameter d 20 2.6 Advantages of the Hybrid Method 22 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core 24 3.1 Hybrid Method Solution for a Pipe and Core of Arbitrary Shape 25 3.2 Comparison of Hybrid Method Results to Exact and Numerical Solutions 27 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue 34 4.1 Oxygen Transport Model and Assumptions 36 4.2 One Circular Capillary in a Circular Tissue Domain 40 4.3 Solution using the Method of Matched Asymptotic Expansions : . . . . 43 4.4 Solution using the Hybrid Method 48 4.5. Time Estimate to Reach Steady State 50 4.6 Finding the Modified Green's Function 53 4.7 Computed Results and Discussion 56 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies 63 5.1 Singular Nature of the Problem 64 5.2 Matched Asymptotics for an Arbitrary Flow Field 66 5.3 A Higher-Order Term in the Expansion 69 5.4 Two Specific Flow Fields 72 5.4.1 Array of Cylindrical Bodies in Uniform Flow 72 5.4.2 Array of Cylindrical Bodies in Shear Flow 73 5.5 Results and Discussion 75 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder 81 6.1 The Stokes Region 83 6.2 The Oseen Region 84 6.3 The Hybrid Formulation 85 6.4 Numerical Solution of the Hybrid Related Problem in the Oseen Region 86 iii Table of Contents 6.5 Calculating the Lift Coefficient CL 89 6.6 Determining M for an Ellipse 91 6.7 The Linearized Oseen Solution, $ o c , and its Derivatives 93 6.8 Results and Discussion 96 Chapter 7. Other Applications of the Hybrid Method 101 7.1 A Non-linear Model Problem 101 7.2 Low Reynolds Number Flow Past an Array of Cylindrical Bodies 102 7.3 A Biharmonic Eigenvalue Problem 104 7.4 Extensions to the General Framework 110 Chapter 8. Conclusions 113 Bibliography 115 iv Acknowledgement I wish to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for its support in the form of a two-year postgraduate scholarship. A huge thank you to my supervisor, Michael Ward, for his generous financial support and for imparting some of his mathematical wizardry to me. Working with Michael was a truly rewarding experience. Thank you to Brian Seymour, my co-supervisor, for his monetary contribution to my cause as well as his participation in discussions regarding my thesis. Thank you to Uri Ascher (Computer Science), Leah Edelstein-Keshet, Anthony Peirce and Steven Rogak (Mechanical Engineering), the other members of my supervisory committee, especially for their suggestions at the time of my thesis proposal. I also thank Judy Maxwell for giving me boundless support, both as my friend and as administrative assistant of the Institute of Applied Mathematics. To my parents, Archie and Ruth, and to my siblings, Leslie, Carole, Richard, Geoffrey and Lianna: thank you for setting such a high standard for what it means to be a Titcombe. Thank you to Julie, for teaching me what friendship is all about. To Benoit, thank you for your infinite patience in waiting for me to finish my thesis. You are the best! There are so many other people that I wish to acknowledge for their support during my doctoral program, but to name them all individually, I would run the risk of having an acknowledgements section that is longer than my dissertation! To all of my friends: thank you for listening to my woes, for keeping my life in balance and most of all, for making me laugh. Finally, I thank my friend, Fred Weinberg, for giving me inspiration, confidence and a much-needed diversion from my thesis in our weekly flute duets. Fred passed away in March of this year before I had the chance to thank him for sharing music with me. I dedicate this thesis to his memory. v Chapter 1 Introduction and Background This thesis explores applications of a method combining techniques of asymptotic analysis and numerical analysis, which we call the hybrid asymptotic-numerical method, or simply, the hybrid method, to improve the accuracy of approximate solutions to certain singular perturbation problems in two dimensions. We construct mathematical problems that model certain physical processes, such as fluid flow or heat conduction, in an attempt to further our understanding of the world around us. Often these mathematical problems are extremely difficult (or even impossible) to solve exactly and in these circumstances, approximate solutions are necessary. In certain mathematical problems, the governing equations contain dimensionless parameters, such as the Reynolds number or the Peclet number, that can be either large or small. In these limits of extreme values of the dimensionless parameters, one can obtain an approximate solution through the use of perturbation methods. The basic idea of these methods is to begin with the solution to a simpler problem (as a first approximation) and then to obtain systematically better and better approximations. In this way, one finds an asymptotic solution to the problem as an expansion of these successive approximations. A specific type of mathematical problem that one can solve using these approximation tech-niques is a singular perturbation problem. This is a problem for which no single expansion is valid in the entire solution domain: it is necessary to construct more than one perturbation expansion, each valid in a certain subregion of the whole solution region. The singular perturbation problems that we consider in this thesis are two-dimensional prob-lems with strongly localized perturbations whose solutions involve infinite expansions in powers of ( — 1/loge), where s is the small perturbation parameter. (The negative sign is present merely for future notational convenience.) The removal of a small subdomain from a given solution domain is an example of a strongly localized perturbation. These are unusual singular perturbation problems in the sense that the small parameter does not multiply the highest derivative in the governing equation. Thus, the non-uniformity is not due to a reduction in order of the differential operator as in the more familiar singular perturbation problems. The non-uniformity of the singular perturbation problems considered in this thesis is due to a near singularity resulting from the removal of a small subdomain or from a region degenerating to a point in one of the solution domains. In their text on multiple scale and singular perturbation methods, Kevorkian & Cole [26] refer to this type of problem as a "singular boundary problem". A n analytic technique available to treat singular perturbation problems is the method of 1 Chapter 1. Introduction and Background matched asymptotic expansions: the asymptotic expansions in the various solution domains are required to match in some overlap domain of validity, thereby creating an approximate solution that is uniformly valid over the entire solution region. In his text on perturbation methods in fluid mechanics, Van Dyke [56] provides a brief historical account of the method of matched asymptotic expansions. He describes the evolution of the technique from its be-ginning as a generalization of the boundary layer theory of Prandtl in 1905 [43], through to its widespread use in the 1950's in studying viscous flows. In particular, Kaplun [24] used the method of matched asymptotic expansions in resolving the Stokes paradox in low Reynolds number fluid flow in two dimensions. We will revisit this classic problem in a subsequent chap-ter. The method of matched asymptotic expansions is limited in practicality due to the often increasing level of difficulty of obtaining the unknown coefficients at each subsequent order in an infinite recursive set of problems. Another approach to finding an approximate solution has been by solving these mathematical problems numerically. Specifically, numerical analysis of singular perturbation problems began receiving significant attention in the 1970's, with the first conference on the subject in 1978 [16]. More recently, Roos et al. [47] produced a text of numerical methods for singularly perturbed differential equations in convection-diffusion and flow applications. However, neither that text nor the papers in the conference proceedings consider any "singular boundary problems", which is indicative of the lack of a systematic numerical study of problems of this type. Although there has been significant development in computer codes for solving a wide variety of singularly perturbed partial differential equations, using finite difference methods or finite element meth-ods, for example, one should never use any computational package without first being aware of the potential difficulties of the equations to solve. In particular, the full numerical approach has limitations since it is difficult, without extensive computation, to isolate the parameter dependence and to resolve the rapidly varying scale structure that is inherent in the problems that we consider in this thesis. One of many methods that combines techniques from both of the major approaches to find-ing approximate solutions is the hybrid asymptotic-numerical method of Ward, Henshaw &: Keller [57], which they used to treat certain classes of eigenvalue problems in a bounded, two-dimensional domain with small perforations and with certain boundary conditions on the resulting holes. These types of problems are examples of strongly localized perturbation prob-lems. Using the method of matched asymptotic expansions, they demonstrated that in some cases the eigenvalue expansion for the singularly perturbed problems begins with an infinite logarithmic series in powers of (-1/ log(£c/)) . Here, e represents the size of the domain perfo-ration, and the constant d depends on the shape of the hole and on the boundary condition of the hole. Their hybrid method used asymptotic analysis to formulate a related problem, whose solution contained the entire infinite logarithmic series. The related problem was non-stiff and straightforward to solve numerically. By summing all the logarithmic terms in the series, the approximation gave accurate results since the error was smaller than any power of ( - l / log(ed)) . In 1995, Kropinski, Ward & Keller [29] appropriately modified the hybrid method to treat steady, two-dimensional, incompressible fluid flow at low Reynolds number past a symmet-ric, cylindrical body. The asymptotic expansions for the drag coefficient and for the velocity field in the limit of low Reynolds number also commence with infinite logarithmic series. They 2 Chapter 1. Introduction and Background implemented a straightforward finite difference scheme to solve the asymptotically related prob-lem, which was independent of the cross-sectional shape of the cylindrical body. This exploited Kaplun's equivalence principle [24], which established an asymptotic equivalence between cylin-ders of different cross-sectional shape, based on an "effective" radius of the cylinder. The crux of applying the hybrid method to strongly localized perturbation problems is to use the method of matched asymptotic expansions to formulate an asymptotically related problem, in which a singularity structure replaces the localized perturbation. The resulting related problem turns out to be quite easy to solve numerically. In solving this related problem, it is then possible to determine the asymptotic solution of the original problem that is correct to all logarithmic terms. The hybrid method has been successful in improving the approximate solution in the area of strongly localized perturbation problems that contain an infinite expansion S(e) of logarithmic terms for the quantity of interest of the form The traditional difficulty of asymptotic solutions containing infinite logarithmic expansions is that these solutions converge very slowly, which means that one must retain a great number of terms of the expansion in order to obtain sufficient accuracy at moderately small values of the perturbation parameter. The hybrid method applies techniques of asymptotic analysis and simple numerical analysis to circumvent the slow convergence difficulty and thus, to improve the accuracy of the approximate solution. In general, asymptotic expansions need not converge. By definition, the error made in truncat-ing the asymptotic solution is of the order of the first neglected term. Even when the expansions do converge, if they are of the form S(e), the slow convergence throws the usefulness of the asymptotic solution into question. The infinite series appearing in the solutions to the problems that we consider in this thesis are not only asymptotic but we believe also converge for small e. We will demonstrate the convergence problem of the reciprocal logarithmic asymptotic series by contrasting it to the more familiar type in powers of e. We define P(s) to be the power series in s of the form £ 0. (1.1) CO P(e) ~ ^ T a j e ' - 1 , £ -»• 0. We will take cij = 2 1 J in both S(e) and P(e). For s less than approximately 0.6, S(e) with aj = 2 1 _ J is a convergent infinite geometric series, which has the exact sum of Similarly, the exact sum of P(s), for the same cij, is P(£) = 2 2 - e 3 Chapter 1. Introduction and Background In the top graph of Figure 1.1, we plot the five-term partial sum, S5(e), and the exact sum S(e) versus e for the reciprocal logarithmic series. In the lower graph, we show the corresponding plots for the power series in e. One can see that the truncated reciprocal logarithmic series only compares well to the exact sum for small values of e, whereas the truncated power series is virtually indistinguishable from the corresponding exact sum curve. co 3 -2 -0 1 1 1 1 1 1 1 1 i i i I 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0 91 1 1 i i i i i i i i I 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 e Figure 1.1: Five-term partial sum of a convergent asymptotic series and the exact sum of the series versus e. Top graph: reciprocal logarithmic series S(e). Bottom graph: power series P(e). One may wonder if the hybrid method could be extended to treat three-dimensional singular perturbation problems. The answer is no. This is not due to a shortcoming of the method, but rather due to an important difference between two-dimensional and three-dimensional problems. The hybrid method was designed to overcome the slow convergence difficulties of asymptotic solutions containing infinite reciprocal logarithmic series in the small perturbation parame-ter. For example, the singularity structure of certain problems that we consider calls for a free-space Green's function representation of the Laplacian operator in the solution. In the two-dimensional case, the free-space Green's function for the Laplacian operator has the form log r, where r is the distance measured from the singularity, which allows for the reciprocal logarithmic form of the asymptotic expansion of the solution. In contrast, the corresponding three-dimensional free-space Green's function is of the form 1/r, which will not generate a re-ciprocal logarithmic expansion with its inherent convergence difficulties and hence, the hybrid method is not necessary. It also may seem that singular perturbation problems whose solutions involve infinite expansions in powers of (—1/loge) are uncommon, but we will see that such problems arise in a variety of contexts. The types of strongly localized perturbation problems involving infinite reciprocal 4 Chapter 1. Introduction and Background logarithmic expansions that comprise the bulk of this thesis fall into four categories: low Peclet number convective heat transfer, low Reynolds number fluid flow, oxygen diffusion in muscle tissue from capillaries, and fully developed laminar fluid flow in a straight pipe. This thesis continues the development of the hybrid method by applying it to singular per-turbation problems that arise in a wide range of disciplines. In Chapter 2, we elaborate on the class of applicable singular perturbations problems and describe the steps of the hybrid method using an illustrative example. In Chapters 3 and 4, we apply the hybrid method to certain problems on bounded domains. In Chapters 5 and 6, we apply the hybrid method to certain problems on unbounded domains. In Chapter 7, we outline extensions of the hybrid method applications to problems in scattering theory, plate vibration and others, and finish with conclusions in Chapter 8. Chapter 2 provides information on the hybrid asymptotic-numerical method. In particular, we state conditions that a singular perturbation problem must satisfy in order to produce the asymptotic structure of the solution that the hybrid method requires. We offer a common framework for the governing equations and boundary conditions for a second-order, steady, singular perturbation problem to which one may apply the hybrid method. We outline the steps of the hybrid method on a linear version of this framework on a bounded domain, with small detours to explain how to modify the method for non-linear problems. We explore the effect on the asymptotic solution structure of relaxing the necessary conditions on the singular perturbation problem. The hybrid method treats strongly localized singular perturbation prob-lems whose asymptotic solution contains an infinite series in terms of reciprocal logarithms of the form in (1.1). The subdomains that are removed from the solution domain are of arbitrary cross-sectional shape, with a unique shape-dependent parameter d for a particular shape. In the chapter describing the hybrid method, we introduce this parameter d more formally, and show how to determine its value. We distinguish the asymptotic and the numerical portions of the method, and describe its advantages over using only the method of matched asymptotic expansions or only a direct numerical computation. Chapter 3 entails the application of the hybrid method to fully developed laminar flow in a straight pipe containing a thin core. The small parameter e in this problem is a measure of the size of the core cross-section. We derive the pipe flow equations from the Navier-Stokes equations and apply the hybrid method to obtain an approximate solution for the axial velocity of the fluid in a pipe with a core that are both of constant but arbitrary cross-sectional shape. Since this application has a direct link to the linear, second-order singular perturbation problem on a bounded domain that we use in Chapter 2 to outline the hybrid method, we only repeat certain key details of the application. For the special pipe-core geometries of a concentric annulus and an eccentric annulus, we compare the hybrid solution results to those of an exact or fully numerical solution in terms of the mean flow velocity or the friction coefficient. Chapter 4 describes in detail the application of oxygen transport from multiple capillaries to skeletal muscle tissue. A mathematical model of the transport of oxygen from capillaries to skeletal muscle tissue is a diffusion problem in a two-dimensional, bounded domain with Neumann and mixed boundary conditions. We consider N capillaries of small but arbitrary cross-sectional shape and demonstrate that, for N > 1, this is a singular perturbation prob-lem that involves an infinite expansion of logarithmic terms of the small parameter £, which characterizes the size of the capillary cross-sections. We apply the hybrid method to solve ap-5 Chapter 1. Introduction and Background proximately for the steady-state oxygen partial pressure in the tissue. In general, our oxygen transport model incorporates the effects of tissue heterogeneities such as mitochondria, variable permeability of the capillary walls and the facilitation of oxygen transport by the presence of myoglobin. We demonstrate the asymptotic results with some specific examples to illustrate these and other effects, and in certain cases, compare with the exact or numerical solution. Chapter 5 shows the hybrid method applied to a convective heat transfer problem past multiple, arbitrarily shaped bodies. Convective heat transfer from an array of small, cylindrical bodies of arbitrary shape in an unbounded, two-dimensional domain is a singular perturbation problem involving an infinite logarithmic expansion in the small parameter e, representing the order of magnitude of the size of the bodies. Using the hybrid method, we formulate a related problem to solve for an approximate solution for the dimensionless, steady-state temperature. We assume that the velocity field of the fluid surrounding the bodies is arbitrary but known. From our asymptotic solution for an arbitrary velocity field, we present the results for two special cases: a uniform flow field and a simple shear flow field. We demonstrate the asymptotic results of the hybrid method through a number of examples and, in a particular case, we compare to an exact analytical solution. Chapter 6 describes the hybrid method applied to the non-linear problem of low Reynolds number fluid flow past a body asymmetric to the uniform free-stream. Low Reynolds number fluid flow past a small, cylindrical body of arbitrary shape in an unbounded, two-dimensional domain is a singular perturbation problem involving an infinite logarithmic expansion in the small parameter e, representing the Reynolds number. Using the hybrid method, we construct a related problem that we solve to obtain an asymptotic solution for the dimensionless, steady-state velocity field and for the coefficients of lift and drag. This application extends the work of Kropinski et al. [29] of low Reynolds number fluid flow past a symmetric, but otherwise ar-bitrarily shaped, cylinder. We modify their finite difference code to incorporate the asymmetry in the flow field, which we use to compute an asymptotic solution for the coefficient of lift that is correct to all logarithmic terms. Chapter 7 suggests possible extensions of the hybrid method to other applications that include unsteady problems, non-linear problems and fourth-order eigenvalue problems. One application is a non-linear problem that extends the convective heat transfer analysis of Chapter 5 and is also a rough model for steady viscous incompressible fluid flow. We describe an extension of the low Reynolds number problem of Chapter 6 to fluid flow past an array of cylindrical bodies that are symmetric to the free-stream. We provide some details of the analysis on a biharmonic eigenvalue problem that is a linear, fourth-order problem on a bounded domain which models the vibration of thin plates with small holes. To close the chapter on possible extended applications, we touch on problems that model such physical processes as non-linear reactions on catalytic surfaces and low frequency scattering of light. In the final chapter, Chapter 8, we reflect on the hybrid asymptotic-numerical method as a powerful tool in treating singular perturbation problems involving infinite reciprocal logarithmic expansions. We comment on the advantages of applying the hybrid method to these problems over using only the method of matched asymptotic expansions or only a full numerical approach. Finally, we summarize the wide variety of contexts in which singular perturbation problems of this type occur. 6 Chapter 2 The Hybr id Asymptotic-Numerical Method In this chapter, we answer certain questions about the hybrid asymptotic-numerical method. What conditions on singular perturbation problems are essential for the hybrid method to be of use? What happens to the asymptotic structure of the solution if we relax any of these conditions? What are the basic steps of the hybrid method? To aid in answering these questions, we provide a framework for the governing equation and boundary conditions of applicable singular perturbation problems. This framework includes steady, second-order, two-dimensional problems that are linear or non-linear, and that are on bounded or unbounded domains. We touch on some generalizations to the framework in this chapter but we leave the detailed discussion for later in Chapter 7. We state the necessary conditions on a strongly localized singular perturbation problem to produce certain features in the asymptotic structure of the solution. By enumerating the basic steps of the hybrid method on a second-order linear problem, we illustrate how the hybrid method exploits this asymptotic solution structure. Throughout the illustration, we discuss any modifications to the steps of the hybrid method for non-linear problems. After, we explore what happens when we relax any of the conditions on applicable singular perturbation problems. 2.1 Conditions on Applicable Singular Perturbation Problems We describe necessary conditions on singular perturbation problems to produce the asymptotic solution structure that makes the hybrid method useful. Firstly, applicable singular pertur-bation problems are two-dimensional whose governing equation involves an operator with a logarithmic fundamental solution (eg. the two-dimensional Laplacian operator). The problems must also be strongly localized singular perturbation problems, such as those on a solution domain with the removal of a small subdomain. A second condition is that the boundary condition on the subdomain border must contain a Dirichlet component. Finally, the unper-turbed solution of the global expansion in the singular perturbation problem must satisfy a non-degeneracy condition. These conditions are necessary to produce certain features in the asymptotic solution for which the hybrid method is applicable. One such feature is that the asymptotic solution must involve reciprocal logarithmic gauge functions of the form ( —1/loge), where e is the perturbation parameter. The other feature is that the local solutions, valid in a region close to the removed subdomain, must be constant multiples of a canonical local solution. In the next section, we 7 Chapter 2. The Hybrid Asymptotic-Numerical Method introduce a framework for applicable second-order steady problems. Then, we explain how the necessary conditions give rise to the essential features of the asymptotic solution structure on a linear version of the general framework. 2.2 Framework for Applicable Second-Order Steady Problems We present a second-order, steady singular perturbation problem for which we want to obtain an approximate solution for the unknown $(x; e), where s is the small perturbation parameter. A framework for the governing equation and boundary conditions of $(x;e) — , x-i\s), on a bounded or unbounded domain Q, is x G Q\Q,0 C TZ2, (2.1a) xe<9ft 0, & # 0 , (2.1b) x e dil, (2.1c) |x| —> oo. (2.Id) First, we describe the structure of the governing equation. In (2.1a), 7V(x,$, V $ ) is a scalar linear or non-linear function and c is of one sign. Of course, one would impose certain conditions (on A ( x , $ , V $ ) , for example) to guarantee the existence of a solution. To illustrate the hybrid method, we will simply assume that a solution exists to (2.1). The domain Q is a two-dimensional domain that is bounded or unbounded, from which we remove f i 0 , a finite collection of k subdomains whose distance apart is 0(1). If ft is a bounded domain, then the distance from the boundary of to each of the k subdomains is also 0(1). Unless we specify otherwise, boldface variables (eg. the spatial variable x = (x\,X2)) represent vectors in two dimensions. This framework of applicable singular perturbation problems is by no means fully general. For instance, we can modify this form to include eigenvalue problems (for example, the original application in Ward et al. [57] in 1993), and also extend it to include unsteady problems. We discuss these and other extensions to the general framework of applicable problems, including a fourth-order eigenvalue problem, in Chapter 7. Now, we move to the form of the boundary conditions on the subdomain border, <90o, of £ir> A requirement of the boundary condition in (2.1b) is that the Dirichlet component must be present, i.e. 6 ^ 0 . To describe a purely Dirichlet boundary condition, we would set b = oo. Also, in (2.1b), the function 3>o may depend on the spatial variable x and d/dn represents the outward normal derivative to the domain. In Section 2.4, we will demonstrate that the asymptotic solution structure cannot contain reciprocal logarithmic gauge functions if we relax the requirement of the Dirichlet component in the boundary condition on <9$7o- We will also explore the ramifications of relaxing the other conditions on applicable singular perturbation problems. If the solution domain fi is bounded, the problem to solve is (2.1a) with (2.1b) and the boundary condition on its outer boundary d£l in (2.1c). Again, d/dn is the outward normal derivative to the domain. The outer boundary condition could be purely Neumann with B = oo, or purely Dirichlet with B = 0, or a combination of both. V - [ c V $ ] + A ( x , $ , V $ ) = 0, — e— + fc($ - $ 0 ) = 0, on ft Bounded - r - + 5 $ = 0, an or, Q. Unbounded $ ~ $co, 8 Chapter 2. The Hybrid Asymptotic-Numerical Method If Q, is unbounded, we would solve (2.1a) with (2.1b) and we impose a far-field condition of the form in (2.Id), in which may depend on x. In the next four chapters, we delve into applications of the hybrid method on specific strongly localized singular perturbation problems. In Chapter 3, we solve for the axial velocity, w, of fully developed laminar flow in a straight pipe with a core. For this special case, the governing equation is linear and the solution domain is bounded. The problem to solve is (2.1a)-(2.1c), in which Also, for the pipe flow application, c = 1 in (2.1a) and $ 0 = 0 in (2.1b). Chapter 4 contains our detailed examination of a second application on a bounded domain that is also linear. In this application, we examine the oxygen partial pressure in a transverse section of skeletal muscle tissue that is supplied by multiple capillaries. We solve for the oxygen partial pressure, p, in the tissue and the corresponding variables of (2.1a)-(2.1c) for this problem are In (2.1b) for the oxygen transport application, and $ 0 = Pck are specified constants. Our first application on an unbounded domain is in Chapter 5, which is a linear convective heat transfer problem about an array of cylindrical bodies. We solve for the temperature, u(x;e), in (2.1a), (2.1b) and (2.Id) with For the convective heat transfer application, c = 1 in (2.1a). As well, in (2.1b), b = Kk and $ 0 = a f c are specified constants, and in (2.Id), = 1 represents the ambient temperature in the unbounded solution domain. The three problems that we have just mentioned are linear and second-order. In Chapter 6, we study low Reynolds number fluid flow past an asymmetric cylindrical body. This problem $(x; e) = w(x; e) (axial velocity component) N — (3 (positive constant) = D (bounded domain of the pipe) i7o = De (k — 1 small subdomain, the pipe core) b = oo (Dirichlet condition on subdomain border) B = oo (Dirichlet outer boundary condition). $(x; e) — p(x; e) (oxygen partial pressure) c = 'P(x) (spatially dependent) iV(x, V$) = -Ai(x) (spatially dependent) f i = D (bounded domain of tissue) Oo = U/Jfe (k small subdomains, the capillaries) b = K f c / ^ ( x ) (mixed boundary conditions on capillary walls) B — 0 (Neumann outer boundary condition). <fr(x; e) = M ( X ; e) JV(x,$, V$) = -v • V« Q = n2 fio = y/j>£ (temperature) (linear function with v known) (unbounded domain) (k small subdomains, the cylindrical bodies) (mixed boundary condition on cylindrical bodies). b = Kk 9 Chapter 2. The Hybrid Asymptotic-Numerical Method is non-linear and on an unbounded domain. We present the low Reynolds number problem as a fourth-order non-linear problem to solve for the stream function, although it is possible to formulate it as a second-order non-linear system in terms of velocity and pressure. In this way, we extend the general framework of applicable singular perturbation problems to include the fourth-order biharmonic operator. Later in this chapter, we elaborate on how to apply the hybrid method to non-linear problems. Much later in this thesis, in Chapter 7, we outline a biharmonic eigenvalue problem that models the vibration of a thin plate with small cutouts or concentrated masses. So far, we have provided a general framework for possible applicable singular perturbation problems. We mentioned that the problems must be two-dimensional and must contain a Dirichlet component in the boundary condition on the collection of subdomain boundaries. Also, we require that the unperturbed problem satisfy a non-degeneracy condition. We illustrate the necessity of these conditions on the asymptotic solution structure of the perturbation problems by stepping through the hybrid method on a linear version of (2.1) on a bounded domain. 2.3 Outline of the Hybrid Method on a Second-order Linear Problem Since three of the four major applications of the hybrid method in this thesis are second-order linear problems, we outline the steps of the method on a general strongly localized singular perturbation problem of a similar form. We illustrate the steps of the hybrid method by applying it to a second-order, linear problem involving the Laplacian operator on a bounded domain. In so doing, we clarify the role of the conditions on the singular perturbation problems that the hybrid method can treat. Also, we comment on necessary modifications for non-linear problems. We can basically summarize the hybrid method in five steps. The first step of the hybrid method is to apply the method of matched asymptotic expansions to the singular perturbation problem to obtain the asymptotic structure of the solution. In particular, this step determines if the asymptotic solution involves reciprocal logarithmic gauge functions. It also determines if the local solutions, valid close to the removed subdomain, are constant multiples of a canonical solution. These two features of the asymptotic solution are essential for applying the remaining steps of the hybrid method. We consider a linear version of (2.1a), where N is now a linear function (which we will henceforth refer to as L) and O is a bounded domain D. Also in (2.1a), we set c = 1, which means that the first term involves the two-dimensional Laplacian. Although in general fio represents a collection of k subdomains that are removed from Q, we will take k = 1 for ease of illustration, and call the small subdomain Also, we take 3>o m (2.1b) to be a constant C. Thus, the second-order, linear, singularly perturbed problem for our illustration of the hybrid 10 Chapter 2. The Hybrid Asymptotic-Numerical Method method is 0, x e D\D£ C V,2, (2.2a) 0, xedD£, b^O, (2.2b) 0, x e dD. (2.2c) Here and throughout, unless otherwise specified, A represents the two-dimensional Laplacian operator, A = d2 jdx\ + d2 jdx\. In (2.2a), L(x,u,v) is linear in u and v. We describe the five steps of the hybrid method in determining an asymptotic solution to (2.2). For the asymptotic solution, we define two solution regions. One is the global (outer) region, that is valid away from the removed subdomain, D£, located at £, where |x - £| = 0(1). The second is the local (inner) region, that is valid close to the removed subdomain, where |x — £| = 0(e). For an arbitrarily shaped subdomain, £ is the location of its centroid. At times, we will refer to the removed subdomain D£ as the "hole" in the domain D. Step 1. Apply the method of matched asymptotic expansions to ensure two features of the asymptotic solution structure. One feature is that the local problems are the same, so that each local solution is a multiple of a canonical local solution. The second feature is that the gauge functions of the asymptotic expansion are reciprocal logarithms of the form ( —1/loge), where £ is the perturbation parameter. This asymptotic solution structure allows us to define an effective radius of the subdomain D£ and to establish an asymptotic equivalence between solutions on subdomains of arbitrary and circular cross-sectional shape. We refer to this as Kaplun's equivalence principle, since he first remarked on it in his study of steady viscous flow past a stationary body [24]. In the global region, away from the removed subdomain D£, we expand the solution in a general asymptotic series as CO $(x; e) = $0(x) + ^ ( ^ ( x ) + £ ^ ( ^ ( x ) + • • • . (2.3) Here, <J?o(x) is the solution to the unperturbed problem, that is, the problem on the domain D without the removal of the subdomain D£ (and, of course, without the boundary condition on dD£). By definition, the gauge functions Fj(e) in an asymptotic expansion like (2.3) are such that Fj+1{e) = oitye)], or lim = 0. (2.4) In the local region, close to the hole at x = £, we define the local variables y = ^ 7 ^ > <t>(y;s) = <S>(ey + C,£). (2.5) We substitute these variables into the governing equation (2.2a) and the boundary condition on the hole (2.2b) and obtain 0, y £ A > , (2.6a) 0, y G dD0, b ± 0. (2.6b) A $ + Z ( x , $ , V $ ) = + 6 ( $ - C) = dr + 5 $ = A ^ + £ 2 i ( £ y + C , 0 , £ _ 1 V ^ ) on 11 Chapter 2. The Hybrid Asymptotic-Numerical Method Here, the subscript y indicates that differentiation is with respect to the local variables. In (2.6a), we assume that we obtain the region D£ by shrinking the distance from every point on the boundary of a fixed domain Do to the centre by the factor e, and so we write D£ = EDQ. Also, dD0 is the boundary of D0. We expand the local solution in a general asymptotic series as CO <Ky; e) = My) + h{e)My) + £ f^)My) + • • • , • (2.7) where we take </>o(y) = C and where the gauge functions fj(s) also satisfy (2.4). We substitute this expansion into (2.6). Provided that the gauge functions fj(e) also satisfy fj(e) > £2L(ey + £, </>, £ -1Vj,c/>), then the second term on the left-hand side of (2.6a) is a higher-order correction. For now, we assume that fj(s), for j > 1, satisfy this condition. Later, once we choose the form of the gauge functions, we will return to verify that this condition holds. Thus, the local functions 4>j(y) in (2.7) satisfy Ay<f>j = 0, y^D0, (2.8a) ^ + 6 ^ = 0, y e ^ o , b^O. (2.8b) on In addition, a matching condition replaces the missing boundary condition. Looking at (2.8), we see that the local problems for 4>j(y), where j > 1, are the same. That is, we can write each local solution as a constant multiple of some canonical local solution My) = a3My)- (2.9) We determine the constants, a,j, for j > 1, through the matching procedure. In (2.9), c/>c(y) is the canonical local solution satisfying (2.10a) b^O, (2.10b) (2.10c) Since <f>c satisfies Laplace's equation in two dimensions, the far-field behaviour of <f>c is that it grows at worst logarithmically plus some constant. In the far-field behaviour in (2.10c), we have written this constant in the form of (— log cZ) for notational convenience later on. Here, d is a constant that depends on b and on the shape of the hole D0. The solution to (2.10) uniquely determines the constant d for a given shape of Do. For general subdomain shapes, we must solve (2.10) numerically to determine d. In Section 2.5, we describe the numerical procedure to solve (2.10), and for certain subdomain shapes, such as a circle or an ellipse, we show how to determine d analytically. Matching the global and local solutions requires that the expression for $(x;e) a s x - + ( agree with the expression for <j>(y;e) as |y | —> oo in some region of overlap. For the global solution as x —• £, we use a Taylor expansion in (2.3) and obtain $ ~ MO + V $ o ( x ) | x = < : • (x - C) + • • • + i 7 , i (£)$ i (x -> 0 + • • • • (2.11) Ay(j>c = 0, y 0 D 0 , ^p- + b<f>c = 0, y e 8D0, on 4>c ~ log |y | - logci, |y | -* oo. 12 Chapter 2. The Hybrid Asymptotic-Numerical Method Here, the notation 3>i(x —• C) refers to the behaviour of the solution $ i ( x ) as x (. We mentioned at the beginning of the chapter that the unperturbed solution must satisfy a non-degeneracy condition. We will see that this condition on the unperturbed solution is MO ± C. (2.12) For the local solution behaviour as |y | —»• oo, in terms of global variables, we use (2.9) and (2.10c) in (2.7) and obtain <t> ~ C + / i (£ )a i [ log |y | - l o g d ] + f2(e)ci2[log |y | - log cl] + ••• ~ C + /i(£)ai[log |x - C| - log(ed)] + / 2 a 2 [ log |x - C| - log(ed)] + • • • . ( ' " ' Here, the (— log d) form of the constant in the far-field behaviour of 4>c from (2.10c) conveniently allows us to form the (log(£ci)) term. The matching proceeds without difficulty if we choose the gauge functions as Wi th this choice, at the 0(1) level, the matching requirement is M0 = C + a1. (2.15) Thus, the solution to the unperturbed problem determines the constant a\. Then, with ci\ known, we obtain ci2 from the solution to the problem for $ i ( x ) which is of the form A $ i + X ( x , $ 1 , V $ i ) = 0, x £ D \ { ( } , (2.16a) ^ + 5 $ i = 0, x G c l D , (2.16b) on $ x ~ a x l o g | x - C| + a 2 , x C- (2.16c) At this stage, we see the necessity for the form of the non-degeneracy condition on the unper-turbed solution in (2.12): without it, (2.15) would require that a\ = 0, which would eliminate the logarithmic behaviour of $ i as x ^ C- The matching procedure continues in this manner and we obtain an infinite recursive set of global problems to solve to determine the unknown coefficients cij in (2.9). At the 0{v3) level in the matching, the behaviour of the solution $ j (x) , for j > 1, as x —> C> is $j ~ a j l o g | x - C | + aj+i, x - * C - (2-17) Although the choice of gauge functions is not unique, it is only the reciprocal logarithmic form of gauge functions that would permit us to determine all the free constants in the asymptotic expansions. Also, this choice of local gauge functions, fj(e), satisfy the condition fj(s) ~> e2L that allowed us to neglect the second term on the left-hand side of (2.6a). Thus, we have ensured that the gauge functions are reciprocal logarithms of the small perturbation parameter £. Since we will see this reciprocal logarithmic form often in the applications of the hybrid method, we define log(ed) 13 Chapter 2. The Hybrid Asymptotic-Numerical Method In Chapter 1, we described the apparent difficulty of the method of matched asymptotic expan-sions in that the asymptotic solution containing such infinite logarithmic series converges very slowly unless e is very small. Thus, truncating the series is in general not very accurate even for moderately small e. As well, to determine the coefficients cij in (2.9), the method of matched asymptotic expansions requires that we solve an infinite set of problems for $ j (x) , for j > 1. In linear systems, with one hole in the domain, these calculations may not be so difficult. How-ever, for linear problems with multiple holes or for non-linear problems, there is an increasing degree of complexity in calculating each successive coefficient in the infinite expansion. For this reason, using the hybrid method to circumvent these recursive calculations is of great benefit. The subsequent steps reveal how the hybrid method exploits the asymptotic structure of the original problem to reduce it to an asymptotically related problem that is non-stiff and easier to solve than the original. We will show how the related problem also avoids having to calculate the coefficients Oj individually. In general, one must solve this related problem numerically. Step 2. Write the local solution expansion (2.7) from Step 1 in terms of an asymptotic expression A(ed) and the canonical local solution c/>c(y). The asymptotic expression A(ed) is A M ~ | > , ^ — J , e-+0. (2.19) We remark that v 4 ~ 0 ( l ) a s £ — > 0. A(ed) is an asymptotic expression for the infinite expansion of reciprocal logarithms of the small perturbation parameter e, and involves the subdomain shape-dependent parameter d. As we saw in Step 1, c/>c(y) is the canonical local solution to (2.10), which uniquely determines the parameter d for a given b in (2.10b) and for a specific shape of the subdomain DQ. We require the far-field behaviour of the local solution to determine the singularity structure of the hybrid method related problem. We will formulate this related problem in Step 3. To begin, we extend the local solution expansion in (2.7) to <*(y; e) = ¥>o(y; v) + evi(y; v) + eVa(y; v) + • • • • (2.20) Here, cpo(y; v) incorporates all of the logarithmic terms of the local solution, and the remaining terms are beyond all orders of reciprocal logarithms. In (2.20), <po(y; v) is ¥>o(y; v) = C + v{ed)ai<f>c(y) + [v(ed)]2a2<f>c(y) + • • • , (2.21) which we obtained by substituting (2.9) and (2.14) into (2.7), using the definition for v(ed) in (2.18). Again, </>c(y) is the canonical local solution to (2.10). Then, using the definition for A(ed) in (2.19), we write the local solution expansion of Step 1 as 0(y;£) = (Po(y;v) + ••• = C + <^c(y)z>(£(i){a1 + a2v(ed) + a3[u(ed)}2 + •••} (2.22) = C + A{ed)v{ed)4>c{y) + ••• . Now, we determine the far-field behaviour of the local solution in terms of global variables. We 14 Chapter 2. The Hybrid Asymptotic-Numerical Method substitute (2.10c) into (2.22), and use (2.5) and (2.18), to obtain <f> ~ C + A(ed)i/(ed)[log\y\- log d] H ~ t 7 + A(ed)u(sd)[log |x - C| - log(ed)} + ••• (2.23) ~ C + A(ed)z/(£d)[log |x - C| + f - 1 ( £ d ) ] + • • • , |y | oo. In the next step, we use the far-field behaviour of the local solution to determine the singularity structure of the asymptotically related problem that we formulate with the hybrid method. Step 3. Write the global expansion for <5(x; e) in terms of a function $o ( x ! u ) t n a t incorporates all of the logarithmic terms of the asymptotic solution and formulate a related problem for this function. For Step 3, we write the global expansion as $(x; £) = $5(x; v) + £$*(x; v) + £2<3>*(x; »)+••• . (2.24) Here, $o(x> v ) incorporates all of the logarithmic terms of the asymptotic solution and the remaining terms are beyond all orders of reciprocal logarithms. We will formulate a problem to solve for <I>Q that is asymptotically related to the original problem (2.2) for $(x; e). This related problem will contain the asymptotic expression, A(ed) from (2.19), in its solution. The solution of the.related problem for $Q then determines the value of A(sd), which essentially sums the entire infinite logarithmic series without having to determine each aj in (2.9) individually. To formulate a related problem for $ Q , valid in the global region, we substitute (2.24) into (2.2a) and (2.2c), and require that, as x —> £, the solution agree with the far-field behaviour of the local solution in (2.23). We obtain that the global related problem for <J>Q is A $ $ + Z ( x , $ $ , V $ $ ) = 0, x e D \ { C } , (2.25a) f9d>* ^ + £ $ * = 0, x G dD, (2.25b) $ 5 ~ C + A ( e d ) - r A ( e d ) i / ( £ d ) l o g | x - C | , x ->• C (2.25c) The governing equation and outer boundary condition for $ o ( x i v ) a r e °f ^ n e s a m e form as the original problem, except that the subdomain D£ has been shrunk to the point {£}. Also, the localized perturbation in (2.2b) has been replaced by the singularity structure (2.25c), which we determined from matching to the far-field behaviour of the local solution. The condition in (2.25c) specifies the form of both the singular part and the regular part of the singularity structure. In the final steps of the hybrid method, we use this over-determination in the singularity condition to find A(ed). Step 4. Decompose the global related problem solution $5(x> v) m * ° a regular part (satisfying the unperturbed problem) and a singular part (involving a Green's function). The behaviour of this decomposition provides a second expression for the singularity structure of the global solution. The first expression for the singularity structure of $Q is in (2.25c). For Step 4, we decompose the first term on the right-hand side of the global related problem solution in (2.24) as $*(x; v) = $5 f f (x) + 27rA(£ci>(£ci)G*(x; Q- (2.26) 15 Chapter 2. The Hybrid Asymptotic-Numerical Method Here, ( x ) 1 S a regular function that satisfies the unperturbed problem and G'(x; £) is a Green's function with the singularity structure G ~ i - l o g | x - C | + GH(C), x ^ C ( 2 - 2 7 ) Here, GR(C) is the regular part of the Green's function. For the oxygen transport application on a bounded domain that we discuss in Chapter 4, G is a modified Green's function. Wi th the decomposition in (2.26), we can write another expression for the singularity structure of the global related solution. We substitute (2.27) into (2.26), and let x —• C , to obtain $o ~ nH(0 + 2TTAUGr(C) + Av log |x - C | , x -> C- (2-28) Again, we require that the unperturbed solution satisfy the non-degeneracy condition, that is, nH(0 # c. In the corresponding Step 4 for a non-linear problem, where a non-linear function N(x,u,v) would replace £ ( x , u, v), we seek a singular solution to the global related problem in the form $o ~ R + 5 1 o g | x - C I - (2.29) Here, R is the regular part and a S is the coefficient of the singular part. The singular solution defines a function R = R(S). Comparing (2.25c) with (2.29), we find that in this case, R = C + A , S = Av. (2.30) In the next, and final, step of the hybrid method, we obtain an expression for A(ed) by com-paring the two forms of singularity structure for the global related problem solution. Step 5. Compare the two expressions for the singularity structure of the global related problem solution, from Steps 3 and 4 to determine an expression for A(ed). Once we obtain A(ed), then we have determined the global solution $(x;e) in (2.24) correct to all logarithmic terms. This final matching step requires that the two expressions for the singularity structure of the global solution in (2.25c) and (2.28) agree. The log |x — C | terms automatically agree and comparison of the remaining terms gives an expression for A(ed), Here, this procedure is valid provided that <&QH(C) 7^ C and that the denominator is non-zero. To ensure that the denominator does not vanish, we insist that the denominator be of one sign, which then gives a range of validity for the small parameter e. For the expression in (2.31) to be valid, we require that 0 < e < ecr = f r V 2 7 r G « ( c ) , (2.32) where d is the subdomain shape-dependent parameter in (2.10). As well, (ed) ^ 1 so that v(ed) in (2.18) is defined. In this linear case with one hole in the solution domain, the expression for A(ed) is in the form of a geometric series in powers of u(ed). However, this would not be true for non-linear problems or for problems with multiple holes in the domain. We examine two 16 Chapter 2. The Hybrid Asymptotic-Numerical Method multiple hole problems in the oxygen transport application of Chapter 4 and the convective heat transfer application of Chapter 5. For the non-linear version for Step 5, if a non-linear function N(x,u,v) replaces L(x.,u,v), we must assume that a solution exists to x£D\{C}, (2.33a) x e dD, (2.33b) x - > C - (2.33c) We can subtract off the singular part of QH and write the regular part R as R(S) = lim [$„(x ; S) - Slog |x - C|] • (2.34) From (2.30), we have S = Av = v(R-C). Thus, v = S/(R(S) - C) and A = Sjv. We compute solutions to (2.33a)-(2.33b) and then from (2.34), we obtain R(S) for various values of .5', which gives us v(S) and A(S) also. To obtain A(ed), we use v = — l/log(eci). In determining A(ed) through this process, we have essentially summed the infinite logarithmic expansion in (2.19). We employ a technique of this form in Chapter 6 when we examine the non-linear application of low Reynolds number fluid flow past an asymmetric cylindrical body. This completes the outline of the hybrid method on a linear, steady, second-order singular perturbation problem. The hybrid method exploits the asymptotic structure of the original problem to reduce it to an asymptotically related problem. This related problem is non-stiff and contains the infinite logarithmic expansion A(ed) from (2.19) in its solution. The specific asymptotic structure of the solution that the hybrid method requires is that it involves reciprocal logarithmic gauge functions and that the local problems are all the same. We stated necessary conditions on applicable singular perturbation problems to produce this asymptotic solution structure. In the next section, we investigate the effect on the asymptotic solution of relaxing these conditions. A $ „ + J V ( x , $ H , V $ „ ) = 0, - T ! L + B$H = 0, on $ H ~ S l o g l x - C|, 2.4 Effect on the Asymptotic Solution Structure of Relaxing Certain Conditions Now, we explore the effect on the asymptotic solution structure of a singular perturbation problem of relaxing the conditions that we claimed were essential for the hybrid method to be useful. The main purpose of the hybrid method is to improve the accuracy of approximation in the asymptotic solution of singular perturbation problems that involve reciprocal logarithms of e. We stated that an applicable singular perturbation problem must be two-dimensional, with a Dirichlet component in the boundary condition on the removed subdomain, and whose unperturbed solution satisfies a non-degeneracy condition. Without these conditions, we assert that the singular perturbation problem will not possess an asymptotic solution structure with reciprocal logarithms. First, we consider the effect of relaxing the Dirichlet condition on the subdomain boundary, i.e. setting b = 0 in (2.1b). For ease of illustration, we take c = 1 in (2.1a) on a bounded 17 Chapter 2. The Hybrid Asymptotic-Numerical Method domain with one removed subdomain D£. After reseating to the local variables in (2.5), the local problems for </>j(y), j > 1, in (2.7) are Av<f>j = °> y i A ) , (2.35a) dn 9 ^ J - 0, y G dD0. (2.35b) Here, D0 is the scaled subdomain, with D0 = e 1D£. If we impose a logarithmic growth in the far-field behaviour of the local solution of the form <My) ~ CJ l o § |y | + • • • , |yI oo, (2.36) then we show that Cj must be zero. By the Divergence theorem, we have A ^ d x = tim [ ^i-dS+ . [ ^dS. (2.37) R-*co J d\y\ J dn v ; y£Do \y\=R yedDo Here, R is the radius of a large circle surrounding Do and d/dn is the outward normal derivative to the solution domain. Using (2.35a) and (2.35b) in this expression, only the first term on the right-hand side remains. We use the form of (2.36) to obtain tim / ?pdS= tim %2%R = 0. (2.38) ft^co J d\y\ R-+co R y ' |y|=# Thus, the constant Cj must be zero. Without a logarithmic form to the local solutions, the gauge functions cannot be reciprocal logarithms of e. We see that the effect of relaxing the Dirichlet condition the subdomain boundary is that we cannot generate reciprocal logarithms in the asymptotic solution structure. Next, we show that it is impossible to have reciprocal logarithmic gauge functions in the asymp-totic solution of a singular perturbation problem that is not two-dimensional. We use the linear, second-order problem of Section 2.3 in three dimensions for this illustration. The matching con-dition is that the behaviour of the global solution as x £ in (2.11) agree with the behaviour of the local solution in (2.7) as |y| —> oo. The local problem for (f>o(y) satisfies A„c6o = 0, ygDo, (2.39a) ^ + ^ o - C ) = 0, yedDo. (2.39b) From the matching, <f>o(y) has the far-field behaviour <t>o - $o(C), |y| -* oo. (2.40) We write the solution to (2.39) in the form MY) = MO + (c - MO)Mv), ( 2-4i) 18 Chapter 2. The Hybrid Asymptotic-Numerical Method where vo(y) satisfies A„w 0 = 0, y £ D0, (2.42a) ~ + b(v0 - 1) = 0, y e dD0. (2.42b) on y Since Do is three-dimensional in this illustration, the solution t>o(y) has the asymptotic form t>o~(C-*o(0)j^r, | y H o o , (2.43) for some constant ao that we determine from the shape of D0 and the value of b. Since this solution is of order e in the global region, we set the gauge function Fi(e) = e in the global expansion and require that $ i ~ (C-$0 (0 )7 -371 , x ^ C (2.44) l x s I We see that since the local solution does not have a logarithmic form, the gauge functions cannot be reciprocal logarithms of e. We emphasize this point by considering the special case of a small subdomain of radius e. In two dimensions, if D£ = sDo is a circle of radius e, then the solution for (f>(y; s) is of the form c/> = log |y | + C + i . (2.45) Here, we remark that the Dirichlet component of the boundary condition must be present, i.e. 6 / 0. Since the fundamental solution of the Laplacian operator in two dimensions is logarithmic in form, the gauge functions in the expansion (2.7) will be reciprocal logarithms in terms of e. In contrast, if D£ = eDo is a sphere of radius e in three dimensions, their the solution for <6(y; e) is of the form <f>= -rK + C (2.46) |y | b The fundamental solution for the Laplacian operator in three dimensions is not logarithmic in form. Hence, the local solution will not grow logarithmically in the far-field and the gauge functions of the asymptotic solution cannot be reciprocal logarithms. We have shown here that the singular perturbation problem must be two-dimensional to produce the asymptotic solution structure required by the hybrid method. Finally, we consider the effect of relaxing the non-degeneracy condition in (2.12) on the un-perturbed solution on the asymptotic solution. We use the second-order, linear problem of Section 2.3 and now allow *o(C) = C. (2.47) The matching condition is that the behaviour of the global solution as x —> £ in (2.11) agree with the behaviour of the local solution in (2.7), with <f>o = C, as |y | —> 0 0 . Writing the global solution in terms of the local variable, |x — £| = ey, then we set f\(e) = e and match & ~ v$o(C)U=c • y , |y | -> oo- (2-48) 19 Chapter 2. The Hybrid Asymptotic-Numerical Method Tims, the local problem for 4>i(y) satisfies (2.8) and has the far-field behaviour of (2.48). Thus, the local solutions do not have a logarithmic form if we relax the non-degeneracy condition in (2.12), and so the asymptotic solution will not contain reciprocal logarithms. In fact, if we look at the hybrid method solution to the second-order linear problem in Section 2.3, we see this effect right away. If the unperturbed solution $ Q h equals the constant C at the subdomain location £, then in (2.31), A(ed) would be identically zero. This would mean that the infinite sum of the product of coefficients and reciprocal logarithmic gauge functions cancels out completely, which means that the effect of the domain perturbation is extremely weak, being smaller than 0((—1/log(ed)) J) for any j. Hence, it is a necessary condition on the singular perturbation problem that its unperturbed solution satisfy the non-degeneracy condition $o K (C) / C. In the next section, we describe how to determine the subdomain shape-dependent parameter d from the solution to (2.10). 2.5 Subdomain Shape-Dependent Parameter d We determine the solution to (2.10) when b = oo for a specific subdomain profile Do in terms of polar coordinates (p,0). A one-to-one mapping z = / ( £ ) , where z = pel°, transforms the profile Do in the z-plane to the unit circle in the £-plane where the two exterior planes correspond. Under the transformation, <f>c satisfies AoSc = 0, |£| > 1, <Pc = 0, |£| = 1. This problem has the solution <f>c — log|£ | . For |z| large (and thus, |£| large), z ~ PL |£| oo. Since the mapping is unique for each profile Do, this fixes the value of f3. Writing the solution as (f>c = log(|z|/|/3|) = log |z| - log \/3\ and matching to the far-field structure of the canonical local solution in (2.10c), we obtain that d=\P\. When b ^ oo, we cannot determine d using / (£ ) and instead we must compute it numerically. Ward el al. [57] and Kropinski et al. [29] describe a numerical procedure for determining d = d(b), in which they map the exterior of Do to the interior of a unit disc, and then use a finite difference scheme on a uniform polar grid to solve for the stream function, which is modified to remove its singular behaviour. They obtain d by matching the behaviour of the modified stream function near the origin. We can analytically derive d for a cylinder of elliptic shape and circular shape. First, for one elliptic cylinder with definition © 2 + © 2 ^ . ••>'•• 20 Chapter 2. The Hybrid Asymptotic-Numerical Method the solution to (2.10a)-(2.10b) with b = oo is <f>c(y) = i]i - n®. Here, (771,772) are elliptic coordinates defined by 7/1 = c cosh 771 cos 772, 7/2 = c sinh 771 sin 772, c— Va*2 — b*2. Also, 771 = 77° = tanh - 1(6*/a*) is the level line of the elliptic boundary. Matching this solution to the far-field structure of the inner solution in (2.10c), requires that that 77° = log(2d/c). Solving for d, we obtain that , a* + b* d — , 2 ' for an elliptic cross-section with major semi-axis a* and minor semi-axis b*. Next, for one circular cylinder with radius po£, the solution to (2.10a)-(2.10b) with b > 0 is 4>c{y) — log IyI + a. To satisfy (2.10b) at p = \y\ = p0, the constant must be a = - l o g / 9 0 + l/(bp0). Matching this to the far-field structure of the inner solution in (2.10c) requires that a = — log d, giving d = poe- 1^ 0, for a circular cross-section of radius p0e. For the particular case of b = 0 0 , Ransford [45] provides a table of the shape-dependent param-eter d, which we have reproduced in Table 2.1 for certain cross-sectional shapes. D0 d circle, radius r r ellipse, semi-axes a and b a + b 2 equilateral triangle, side h a 0.422/1 871"^ isosceles right triangle, short side h 33 / 4 r m 2 / i square, side h ; 4 L « 0.5902* Table 2.1: Shape-dependent parameter d (capacity) for cross-sectional shape Do — e 1D£, for b = 00 in (2.10b), (adapted from Ransford [45]). The subdomain shape-dependent parameter d is part of the canonical local solution for the second-order singular perturbation problems in Chapters 3, 4 and 5. In the low Reynolds number fluid flow application in Chapter 6, the canonical local problem involves a body shape-dependent matrix M for flow past an asymmetric body. In the case of low Reynolds number fluid flow past a symmetric body, the fourth-order canonical local problem involves an analogous parameter d. 21 Chapter 2. The Hybrid Asymptotic-Numerical Method 2.6 Advantages of the Hybrid Method As its name suggests, the hybrid asymptotic-numerical method combines techniques of asymp-totic analysis and numerical analysis. We describe the advantages of the hybrid method for treating certain, two-dimensional, strongly localized singular perturbation problems over using the two traditional techniques for approximate solutions: the method of matched asymptotic expansions and full numerics. For the asymptotic part, the hybrid method uses the method of matched asymptotic expansions as a means to reduce the original strongly localized singular perturbation problem to an asymptotically related problem. The related problem has a speci-fied singularity structure instead of the localized perturbation on the subdomain boundary that occurs in the original problem. The hybrid method exploits the asymptotic solution structure of the original problem, which must involve reciprocal logarithmic gauge functions of the per-turbation parameter e and must have local problems that are the same, so that their solutions are multiples of a canonical solution. From the solution to this canonical solution, we determine a parameter d that depends on the shape of the subdomain and on the constant b in the subdo-main boundary condition. In some problems, due to the conditions on the asymptotic solution structure, the gauge functions involve d and e only in terms of their product, ed. This feature allows for an asymptotic equivalence between cylinders of different cross-sectional shape, based on an "effective" radius of the cylinder, which is known as Kaplun's equivalence principle. To determine the unknown coefficients in the asymptotic solution, the method of matched asymptotic expansions requires the solution of an infinite recursive set of problems. Often, the problems increase in analytical complexity at each order in the expansion. For non-linear problems, it may be too difficult to determine analytically more than a few of these coefficients. We have seen in Chapter 1 that infinite reciprocal logarithmic expansions, if they converge, converge very slowly and so any truncation of the infinite series may sacrifice the desired accuracy. We also mentioned in Chapter 1 that we believe that these series do converge for small e. The fact that we are able to compute values of A(ed), which is asymptotic to the infinite reciprocal logarithmic series, provides evidence of such convergence. The hybrid method circumvents this slow convergence difficulty by formulating a related problem whose solution actually contains the infinite logarithmic expansion. Thus, by solving this related problem which is easier than for the original, we obtain the asymptotic solution correct to all logarithmic terms. The numerical part of the hybrid method involves solving the related problem, which is non-stiff and has a smaller parameter space than the original problem. For linear problems, we wrote the related problem solution in terms of the unperturbed solution, $Q h (x) , and a Green's function, G(x; £). Computing the unperturbed solution and the Green's function involves numerical techniques for solving partial differential equations such as the finite element method or the finite difference method. For a given domain and locations of the removed subdomains, we only need to compute the Green's function, G'(x;£), and its regular part, GR(Q, once. Then, for a given subdomain shape, Do, we compute a single value d from the solution to the canonical local problem in (2.10). The remaining part of the numerics for a linear problem is to compute the expression for A(ed) and hence, obtain the solution correct to all logarithmic terms, as a function of e. We note that for linear problems with N holes in the domain, we must numerically solve a linear system for Ai(ed), i = 1,... , N, which is the multiple-hole version of the asymptotic infinite sum in (2.19). For non-linear problems, we compute solutions to the related problem depending on a parameter that represents the strength of the singularity. 22 Chapter 2. The Hybrid Asymptotic-Numerical Method Again, in general, this would involve finite element or finite differences techniques. For both linear and non-linear problems, changing the shape of the removed subdomain or the measure of its size is very easy in the hybrid method solution. However, in the full numerical computation of the original problem, such a change would require a new definition of the geometry and of the solution grid. As well, full numerics often have difficulty resolving the rapid change in scale of singular perturbation problems with small holes in the solution domain. Using the hybrid method, it is possible to obtain asymptotic accuracy with less restrictive storage and processor time requirements than numerical solvers. We will demonstrate the advantages of the hybrid method through its detailed application to four problems. The problems that we have chosen as applications for the hybrid method are models in: fully developed steady laminar flow in a straight pipe, oxygen transport from capillaries to skeletal muscle tissue, convective heat transfer, and low Reynolds number fluid flow. The first two models are applications on bounded domains, and the last two on unbounded domains. In the following chapters, we discuss each application in more detail, all the while underlining the features that fink them together. 23 Chapter 3 Viscous F lu id Flow in a Straight Pipe wi th a T h i n Core We consider steady, incompressible, laminar flow in a straight pipe containing a thin core. Both the pipe and thin core have a constant cross-section of arbitrary shape. This problem is the first of two applications of the hybrid method on a bounded domain. The second application, which we describe in the next chapter, is that of oxygen transport from multiple capillaries to skeletal muscle tissue. Using the hybrid method, we determine an asymptotic solution for the flow velocity between the walls of the pipe and the pipe core. Due to our assumptions about the geometry and about the flow, the problem possesses the essential feature for the hybrid method of being two-dimensional. With an approximate solution for the flow velocity, we can compute such qualitative measures as the mean flow rate and the friction coefficient. Some special cases of the straight pipe-core geometry are a concentric annulus and an eccentric annulus. For these special cases, we compare our hybrid results to those contained in the account of fully developed flow in a straight pipe of constant cross-section by Ward-Smith [58]. To derive the equations for this pipe flow problem, we begin with the Navier-Stokes equations for steady, incompressible flow in three dimensions We show how the assumptions of the problem reduce the problem to its two-dimensional frame-work. In (3.1), u = (u, v, w) is the velocity vector and p is the pressure of the fluid, where both u and p are functions of the spatial variable x = (x,y, z). In (3.1a), p is the density and p is the dynamic viscosity of the fluid. The boundary conditions are that u = 0 on the pipe wall and on the core wall. We orient the reference frame so that the positive ^-direction is along the axis of the pipe in the direction of the flow. With our assumptions, the pipe flow is unidirectional and hence, the velocity vector has only the component in the axial direction, u = (0,0, if) . Then, the continuity equation (3.1b) gives that dw/dz = 0 and so we write w = w(x, y). From the x-and i/-momentum equations, the first and second equations of the vector equation (3.1a), we obtain that p = p(z). The ^-momentum equation of the Navier-Stokes equations their reduces p(u • V)u + V p = pAu, V • u = 0. (3.1a) (3.1b) 24 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core to a two-dimensional equation 1 dp d2w 32w udz dx2 dy2 (3.2) Since the left-hand side of this equation is a function of z only and the right-hand side is a function of (x,y) only, we can equate both sides to a constant. Thus, we define the positive constant P = --^, (3-3) / i dz involving the dynamic viscosity \x and dp/dz, the constant negative pressure gradient in the axial direction. The equation for w in (3.2) and the conditions that the axial velocity w vanish at the pipe and core walls comprise a two-dimensional singular perturbation problem on a bounded domain that we will solve approximately using the hybrid method in the next section. 3.1 Hybrid Method Solution for a Pipe and Core of Arbitrary Shape We define D to be the pipe cross-section, from which is removed a small subdomain, D£, representing the cross-section of the thin core. The small parameter £ is a measure of the size of the removed subdomain D£. The problem to solve for the axial component of the velocity, w = w(x; e), is Aw = - /3, x = (x1,x2) G D\D£, (3.4a) w = 0, x G 3D, (3.4b) w = 0, x G c9.De, (3.4c) where A is the two-dimensional Laplacian operator, and dD and dD£ are the boundaries of D and D£, respectively. In Chapter 2, we linked this pipe flow problem to the framework in (2.1) of applicable second-order, steady, singular perturbation problems. For the pipe flow problem, the governing equa-tion is linear and the solution domain is bounded. In (2.1a), we set c = 1 and the function N = j3, a positive constant. As well, both boundary conditions of the pipe flow problem are Dirichlet, for which we set b = oo and $ 0 = 0 in (2.1b) and B = oo in (2.1c). If we were to make a naive attempt at a regular perturbation expansion of (3.4), we would have difficulties at order e. We will briefly outline how a regular perturbation fails for the special case of an annular pipe-core geometry with /? = 1, and then return to the general pipe flow problem to proceed with the singular perturbation approach of the hybrid method. For the special case of an annular pipe-core geometry, with e < r < 1 and j3 = 1, suppose that we expand the solution to (3.4) in a regular perturbation expansion of the form w{r\e) - w0(r) + ewi(r) + e2w2(r) + ••• . (3.5) 25 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core Here, r = |x | . The function wo(r) satisfies the unperturbed problem, that is, the pipe flow problem on a pipe without a core. The unperturbed problem for this special case, whose solution is a regular function as r —> 0, is A w 0 = - 1 , r < 1, (3.6a) w0 = 0, r = 1, (3.6b) which has the solution Mr) = Y- Z^ -- (3.7) At the next order of the problem, which is order e, w\(r) must satisfy A w i = 0 with w\(e) — w\(l) = 0. The only solution to this problem is the trivial solution, Wi(r) = 0. Hence, at order e, we see that we cannot satisfy the condition w(e) = 0 using a regular perturbation approach. We now proceed with the singular perturbation approach for the asymptotic solution to the general pipe flow problem in (3.4) using the hybrid method. In Section 2.3 of the previous chapter, we outlined the hybrid method on a second-order, linear problem on a bounded domain, which is precisely the form of this pipe flow problem. To avoid repetition, we reiterate only certain details of the hybrid method, formulation. In the local (inner) region, close to the subdomain D£ located at x = £, we use the local variables in (2.5), with v(y;e) as the local axial velocity, to write the solution expansion as v(y, e) = V0(y; u) + eV1(y;v) + e2V2(y; »)+•••, (3.8) where the leading-order local solution is Vo(y;v) = v(ed)A(ed)vc(y), following the form in (2.22) with C - 0. The definitions for v(ed) and A(ed) are in (2.18) and (2.19), and. vc(y) is a canonical local solution that satisfies (2.10) with b = oo. Substituting (2.10c) into (3.8), and using (2.5) and (2.18), we obtain that the far-field behaviour of the leading-order local solution, in terms of global variables, is F 0 ( y ; ^ ) - ^ A [ l o g | x - C I + ^" 1] + - - - , |y | ->oo. (3.9) This equation has the same form as (2.23) with C = 0. Following (2.24), we write the expansion in the global region in the form w(x; e) = W0(x; u) + eWi(x; v) + s2W2(x; v) + • • - . (3.10) The function Wo(x; v) will incorporate all of the logarithmic terms through A(ed), which is the asymptotic infinite sum that we defined in (2.19). To formulate a related problem, for Wo, we substitute (3.10) into (3.4a) and (3.4b), and require that it match via (3.9). Thus, we obtain that the related problem for Wo is AW0 = -P, x G D\{C}, (3.11a) Wo = 0, x G 3D, (3.11b) W 0 ~ A ( e d ) + A(ed) i / (ed) log |x-C | , x -»• C (3.11c) 26 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core We decompose the solution for W 0 , which is correct to all logarithmic terms, as W0(x; v) = W O H ( X ) + 2irA(ed)v(ed)G(x; C). (3.12) This follows the form of (2.26), where W0H(x) is a regular function that satisfies the unperturbed problem, AW0H = -P, x € D, (3.13a) W0H = 0, x e dD, (3.13b) and where G(x; C) is the Green's function satisfying A G = 0, x e ^ j C } , (3.14a) G = 0, x G c9D, (3.14b) G ' ~ ^ l o g | x - CI + R(Q, x - > C - (3.14c) In (3.14c), R(C) is the regular part of the Green's function. For the final matching step of the hybrid method, we substitute (3.14c) into (3.12) and compare with (3.11c)'. This provides an expression for A(ed), which is A ( £ d ) ~ l-2xu(ed)R(C)' (3'i5) This expression for A(sd) is valid provided that neither the numerator nor the denominator vanish. Since (5 is a positive constant, we are guaranteed that WQH(£) / 0 by a maximum principle. Requiring that the denominator is non-zero gives us the range of validity for £ in (2.32). Wi th (3.15) and the solutions to (3.13) and (3.14) in (3.12), we have the asymptotic solution for the axial flow velocity, w(x;e), correct to all logarithmic terms. For a given, cross-sectional shape of the pipe and location C of the thin core, we obtain the Green's function (S'(x; C) &ird its regular part R(C) from (3.14); and with a specified constant /?, we find the unperturbed solution WoH(x) from (3.13). Obtaining the Green's function is independent of the shape of the thin core and of the measure e of its size. Then, for a given cross-sectional shape of the thin core, we determine the value of the subdomain shape-dependent parameter d. In general, we determine d from the solution to (2.10), but for b = oo, Table 2.1 in the previous chapter fists values of d for certain regular cross-sectional shapes of the small subdomain. Thus, for a given constant (3 and pipe-core geometry, we easily compute the asymptotic solution as a function of £, where e lies within the range of validity in (2.32). In the next section, we demonstrate the hybrid method results on certain special cases of the pipe flow problem in (3.4) for which we can obtain an exact solution. For cases where obtaining an exact solution is impossible, or too difficult, we compare the hybrid solution with the results from a direct numerical solution. 3.2 Comparison of Hybrid Method Results to Exact and Nu-merical Solutions We compare the results of the hybrid method to those special cases where it is possible to find an exact solution. Otherwise, we compare the hybrid method results with those of a direct 27 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core numerical solution of the original problem. We have used Matlab and its Partial Differential Equations Toolbox [36] for the numerics of the hybrid method and for the direct numerical computations. We indicate the direct numerical results as discrete points in the plots. We draw our comparisons in terms of the mean flow velocity, W, and the friction coefficient, / . The mean flow velocity is W = ^ - I wdx. (3.16) D\D£ Here, A D is the area of the cross-section of the pipe-core geometry, which we also call the duct cross-section. For laminar flow in ducts of non-circular cross-section, with or without cores, one can express the friction coefficient / as f=^C\ = fcsC\, (3.17) where Re = WLp/p is the Reynolds number. Laminar flow, for which there is no significant mixing between adjacent layers of the fluid, occurs for Reynolds numbers in the approximate range 0 < Re < 2000. In (3.17), fcs = 16/Re is the friction coefficient for a pipe of circular cross-section without a core, and C\ is a constant that depends on the pipe-core geometry. A collection of values of C\ for various non-circular duct geometries appears in Table C5 of Ward-Smith [58], a portion of which we reproduce in Table 3.1. In the definition for the Reynolds number, L is a characteristic diameter defined as L = 4AD/PD, where PD is the wetted perimeter of the pipe and core. One can also express the friction coefficient, / , in the form L (-%• / = A z ^ - . (3-18) 2pW Substituting this expression into (3.17), with the definition of (3 in (3.3), we can express the constant C\ as a function of the mean flow velocity, W, and the characteristic diameter, L, as C i = ±U=. (3.19) Z2W V ' Thus, with the mean flow velocity that we obtain using the hybrid method, we will compute this constant C\ and compare to tabulated values in [58] for a special case of the pipe-core geometry. Example: Concentric annulus. We consider flow between the walls of a circular pipe of radius r 0 , ro ^ 1, and a concentric thin core of radius e. Since the geometry is radially symmetric, we solve for w(r; e) in Aiv — -/3, e < r = |x| < r0, (3.20a) w = 0, r = £, (3.20b) w = 0, r = r0. (3.20c) 28 Figure 3.1: Mean flow velocity W versus pipe core radius e for a concentric annulus £ < r < 2. This has an exact solution of P 4 77, - T - TV, logv>0/r) log(ro/e) (3.21) Now, we consider the asymptotic solution for w(r;e) from the hybrid method. From. (3.14), with a circular pipe of radius TTJ and the location of the core at £ = 0, the Green's function is G = - ^ ( log r - logrr j ) . (3.22) Thus, its regular part is R = — (l/27r) log?'o, a constant. For a circular pipe of radius 7"o, the solution to the unperturbed problem in (3.13) is W0H(r) = t{rl - r 2 ) (3.23) Thus, we have that WOH(0) = Prl/4. In the case of a circular core, D£, of radius e, the solution to the canonical local problem in (2.10) with b = oo gives that the shape-dependent parameter is d = 1. Substituting all of this information into (3.15) gives that From (3.12), we get that the hybrid solution, valid for |x| >• 0(e), and correct to all logarithmic terms, is U>(X;E) = - rn-r - r, , log(r 0 / r) log (T-o/e) (3.24) 29 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core Hybrid d Tabulated C\ 0.0001 1.1216 1.1216 0.001 1.1669 1.1669 0.01 1.2518 1.2518 0.05 1.3496 1.3480 0.1 1.4065 1.3964 0.15 1.4554 1.4244 Table 3.1: Constant C\ in (3.19), from the hybrid method solution, and tabulated values of C\ from Ward-Smith [58] for laminar flow through ducts of concentric annular section. Comparing the exact solution in (3.21) with the hybrid solution in (3.24), we see that the hybrid method solution agrees well with the exact solution. We note that the hybrid solution does not capture the 0(e2) term. We compare the hybrid method results to those of the exact solution in terms of the mean flow velocity W in (3.16). Using (3.21) in (3.16), the exact mean flow velocity is l o g ( V £ ) (3.25) Again, A D in (3.16) is the area of the pipe-core cross-section, which in this special case is AD = Tr(r2-£2). In Figure 3.1, for a concentric annulus with ?'o = 2 and (3 = 1, we have plotted the mean flow velocity W from (3.16) versus the core radius e using the exact solution in (3.21) and using the hybrid solution in (3.24). The range of e in our plot lies well within the range of validity in (2.32). The plot indicates that the hybrid method results agree very well with the exact solution, with an error that is 0(e2). Figure 3.2: Geometry of the eccentric annular section of pipe and core. Using the mean flow velocity that we obtained from the hybrid method, we calculate the constant C\ in (3.19) and compare to tabulated values of this constant for laminar flow through ducts of concentric annular section in Ward-Smith [58]. Table 3.1 contains the hybrid method 30 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core values and tabulated values of C\ for several ratios of radii, S/TQ. We see that for ratios of the radii of the pipe and core up to 0.01, the values of C\ from the hybrid method are exactly the four-digit accurate tabulated values. Example: Eccentric Annulus. We consider flow between the walls of a circular pipe D of radius 7 - 0 = 2, containing a circular core D£ of radius s located at £ = (—1,0) and with (3 — 1. The offset of the centre of the subdomain D £ from the centre of D is the eccentricity e (see Figure 3.2). In our computation, e = 1.0. 0.45 0.44 : 0.42 r 0.41 Figure 3.3: Mean flow velocity W versus pipe core radius e for an eccentric annulus, with pipe radius TQ — 2 and eccentricity e = 1. From Ward-Smith [58] and from Piercy et al. [40], using our notation, the mean flow velocity, W, for an eccentric annulus in terms of the volume flow rate, Q, is Q=WAD = f o 4 e 2 M 2 b — a 1 n - £ -where M = (N2 - r2Y< 1, N + M a = — log . 2 N - M' Again, e is the eccentricity. - 8 e 2 M 2 V n e M ~ n ( b + a)) ^ sinh n(b — a) n=l v ' . N TV2, - £ 2 -(- e 2 2e : 1, N + M - e b = — log . 2 N - M - e (3.26) (3.27a) (3.27b) Using the hybrid method, we obtain the asymptotic solution to (3.4) and compare these results to the exact results of the mean flow velocity. The unperturbed solution is in (3.23), since it is 31 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core the same as in the concentric annulus example. For the solution to (3.14), we let G(x;C) = - ^ l o g | x - C | + fl(x;C)- (3-28) Thus, the function i2(x;£) is regular as x —> £ and satisfies AR = 0, xeD, (3.29a) R = - ^ - l o g | x - C | , xedD. (3.29b) We compute -ft(x; Q from (3.29) on a circle of radius ?-o = 2, with £ = ( —1, 0), and the value of R at x = £, using the Partial Differential Equations Toolbox [36] on a mesh of 4064 elements. We use this information, and the unperturbed solution from (3.23) evaluated at £, to compute the value of A(sd) = A(e) from (3.15). The value of £ corresponds to an eccentricity of e = 1. Since the core cross-section is circular, the shape-dependent parameter is d = 1. Thus, with all of this substituted into (3.12), we have the hybrid method solution for the axial velocity in an eccentric annular pipe section correct to all logarithmic terms, which we then use in (3.16) to compute an approximate value for the mean flow velocity. e Figure 3.4: Hybrid method and direct numerical results for the mean flow velocity W versus measure of size e of cross-sectional shape for a circular pipe of radius ro = 2 with a concentric core of three different cross-sectional shapes. In Figure 3.3, for an eccentric annulus with pipe radius r 0 = 2, (3 = 1 and eccentricity e — 1, we have plotted the mean flow velocity W versus the core radius £ from the exact solution and from the hybrid solution. For this example, the plot indicates that the hybrid method results compares reasonably well with the exact mean flow velocity. 32 Chapter 3. Viscous Fluid Flow in a Straight Pipe with a Thin Core Example: Annulus with Various Core Cross-sectional Shapes. We consider flow be-tween the walls of a circular pipe D of radius ?"o = 2 containing a concentric core Dc of various cross-sectional shapes and with f5 = 1. We use Table 2.1 for the shape-dependent parameter d for a square core, an elliptic core, and an equilateral triangular core. Using the notation in the table, we set the major and minor semi-axes of the ellipse as a — 2 and b = 1, and both the side of the square and the equilateral triangle as h = 1. To compute the hybrid method solution, we use the Green's function in (3.22) and unperturbed solution in (3.23) from the concentric annulus example since these solutions are independent of the shape and size of the core. We use this information in (3.15) to obtain A(ed). For a circular pipe of radius TQ — 2 containing a concentric core, Figure 3.4 contains curves of mean flow velocity, W, versus e, a measure of the size of the core, for three different cross-sectional shapes of the core. In the hybrid method, the change in shape and size of the core requires only that we vary the product ed, which allows us to compute the entire e curve very easily. In contrast, for each change of shape and size of the core in the direct numerical solution, we had to recreate the solution geometry and remesh the solution grid. For a core of elliptic cross-section, the figure shows that the hybrid method results agree very well with those of the direct numerical solution. The slight discrepancy in comparing the results for the other two core cross-sectional shapes, the square and equilateral triangle, could be due to the inability of the numerical method to resolve the corners of the core. In the next chapter, we explore a second application of the hybrid method on a bounded domain but with an array of removed subdomains. The application of Chapter 4 models oxygen transport from multiple capillaries to skeletal muscle tissue. 33 Chapter 4 Oxygen Transport from Capillaries to Skeletal Muscle Tissue We apply the hybrid asymptotic-numerical method to an oxygen transport problem in a bounded, two-dimensional domain representing a transverse section of skeletal muscle tissue that receives oxygen from an array of capillaries of small but arbitrary cross-sectional shape (Figure 4.1). Outside of the capillaries, we obtain an asymptotic solution for the oxygen partial pressure (i.e. pressure due to unbound oxygen molecules) in the tissue domain. This is the second of two major applications in this thesis of the hybrid method on a bounded domain. In the oxygen distribution process of the micro-circulation, oxygen binds to its carrier, haemo-globin, in red blood cells, which transports it through the arterioles, branching to the capillary networks, to the collecting venules. In the capillaries, the oxygen is released from its carrier and diffuses into the surrounding tissue. The oxygen transport model that we present in this chapter is an idealization of this full transport process, in which we retain its overall features while producing a mathematically tractable problem. The analytical study of tissue oxygenation from capillaries has been the focus of considerable research since the original work of August Krogh [28] in 1919. The reviews of Popel [42] and Fletcher [10] and references therein provide substantial information on the approaches and advancements of the theoretical work in this area. Hoofd [17] reviews extensions to the Krogh cylinder model of a single circular capillary surrounded by a concentric circular tissue domain (a) Actual skeletal muscle tissue (b) Mathematical idealization Figure 4.1: Capillary blood supply in skeletal muscle tissue. 3 4 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue (Figure 4.2). Following much of this groundwork, we adopt the approach that oxygen transport from capillaries to tissue is a passive process, driven by diffusion rather than by the consumption of energy. Assuming Fick's law, J = —Wc, relating the oxygen flux J to the gradient of oxygen concentration c, and Henry's law, c = ap, and that there is a balance of mass in the tissue, the equation governing the oxygen partial pressure p is a~ - V • [aWp] - M. (4.1) Here, a is the oxygen solubility coefficient, V is the oxygen diffusion coefficient and M is the oxygen consumption rate in the tissue. By balance of mass in the tissue, we mean that the time rate of change of the amount of oxygen per unit volume equals the net diffusion flux through the tissue boundaries plus the rate of chemical reaction within the volume minus the rate of consumption of oxygen. In addition to satisfying the governing equation, the oxygen partial pressure p must satisfy appropriate conditions at the capillary walls and on the outer tissue boundary. We demonstrate that, for N > 1 capillaries in the tissue domain, this is a singular perturba-tion problem whose solution contains an infinite expansion of logarithmic terms of the small parameter e, which characterizes the size of the capillary cross-sections. Figure 4.2: The Krogh cylinder: a capillary of radius, Rc, surrounded by tissue of radius, R. In Section 4.1, we present our mathematical model of oxygen transport from capillaries to skeletal muscle tissue and elaborate on the assumptions behind our model. In particular, we discuss how to include in the model the effects of tissue heterogeneities, such as oxygen-consuming mitochondria, and of enhancement or facilitation of oxygen transport to the tissue by the presence of myoglobin, an iron-protein compound that can reversibly bind up to one oxygen molecule. In Section 4.2, we examine a simple version of the full model problem, of one circular capillary contained in a circular tissue domain as in the original Krogh model cross-section (Figure 4.2), to reveal the form of the asymptotic solution. In Section 4.3, for the special case of capillaries with circular cross-sectional shape, we construct an asymptotic solution for the oxygen partial pressure in the tissue using the method of matched asymptotic expansions. In Section 4.4, for arbitrarily shaped capillary cross-sections, we apply the hybrid method to 35 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue construct an asymptotic solution for the oxygen partial pressure that is correct to all logarithmic terms. In Section 4.5, we compute a time estimate to reach steady state to justify our time-independent model and in Section 4.6, we show how to compute the modified Green's function that occurs in the hybrid solution. In Section 4.7, we demonstrate the asymptotic results with some specific examples that illustrate important physiological effects such as capillary interaction; the cross-sectional shape of the capillaries; variable permeability of the capillary wall; tissue heterogeneities and myoglobin facilitation. 4.1 Oxygen Transport Model and Assumptions We are interested in the steady-state solution to (4.1) in a bounded, two-dimensional domain. One reason that we can view the oxygenation process from capillaries in a two-dimensional domain is the regular longitudinal geometry of skeletal muscle tissue: the arrangement of the capillaries is, for the most part, parallel to the surrounding muscle fibres (Figure 4.3). In this type of muscle tissue, one can orient the £3-axis in the axial direction, running along the capillaries, and the .Ti.x-2-plane in a transverse cut perpendicular to the direction of the capillaries, which we have shown in Figure 4.1(b). We describe another reason for the two-dimensional framework shortly after we state our general mathematical model for the oxygen transport problem corresponding to Figure 4.4. j 1 . 41 »;: § f II ^ mi | 1 ; j] pjj*| H a »: "I ;! Figure 4.3: Photomicrograph of a longitudinally sectioned skeletal muscle, displaying the par-allel arrangement of the capillaries (arrows) between the fibres (reproduced from "Histology: A Text and Atlas" [48]). For our model, we assume that the process has reached steady state. In Section 4.5, we de-scribe the details of determining a lower bound for the time necessary to reach steady-state conditions. This bound depends on the size and shape of the capillary cross-section, the area of the surrounding tissue domain and the diffusivity, V. 3 6 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Figure 4.4: Our general model for i = 1,... , N capillaries of arbitrary shape of magnitude of order e in a tissue domain D with boundary dD. Now, we state the general, steady-state, dimensionless model of oxygenation of a two-dimensional skeletal muscle tissue domain D from N small capillaries Df of arbitrary cross-sectional shape (see Figure 4.4). The mathematical model to solve for the steady-state, dimensionless oxygen partial pressure p(x; e), where x = (xi,x2), including the boundary conditions at each capillary wall dDf and at the outer tissue boundary dD, is x G D\ U Df, (4.2a) j=i x G dDf, i = 1,.. . ,N, (4.2b) x G dD. (4.2c) Here, we rendered the quantities in our model dimensionless with respect to a characteristic length scale of the tissue domain, L*; a characteristic oxygen partial pressure, p*; and the transverse diffusivity, V*. With these scalings, the dimensionless oxygen consumption rate is M = L*2M* I (V£p*). The small, dimensionless parameter, e, represents the order of magnitude of the capillary cross-sections, which we assume to be independent of the axial position £ 3 . Our model problem can be viewed as the leading-order problem of a perturbation analysis of (4.3) in terms of the small perturbation parameter, fi, in (4.4). Following researchers like Hoofd [18], we have decoupled the capillary-tissue diffusion process of oxygen from the transport of oxygen within the capillaries. As we will describe shortly, the generality of our model allows us to analyze more fully the diffusive process in two dimensions before potentially embedding it in a more complex coupled process. In addition to the regular geometry of skeletal muscle tissue, another reason for the two-dimensional framework is that the axial diffusion term is small relative to the other terms in the governing equation. For simplicity, we illustrate the balance of terms for the special case of a tissue cylinder of radius 1 with constant diffusivity enclosing a single concentric capillary whose cross-sectional shape depends only on the axial variable, £ 3 = z. We can write the gov-erning equation for the oxygen partial pressure p(r, z; fi,e) in the tissue volume, which includes V • [VVp] = M, dp Z^TT + Ki(P - Pa) = 0, on ^ = 0 dn ' 37 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue the axial diffusion term, as dz2 ef(z) <r <1, 0 < z <1. (4.3) Here, n is a small parameter defined by ) 2 (4.4) Typically, in skeletal muscle tissue, the intercapillary separation, i * , is a few microns whereas the capillary length, L*, is on the order of a thousand microns. As well, diffusivity in the axial direction, V*, is small relative to that in a transverse cut, V*. The diffusivity, V = oD, is the product of the oxygen solubility and diffusion coefficients. The capillary has a radius of sf(z), where e is a measure of the size of the capillary cross-section. For the analysis that we show in Sections 4.3 and 4.4 to be valid, we require that u2 <C ( — l / l o g £ ) a , for any a. The boundary condition (4.2b) models the capillary wall as a finitely-permeable membrane, where Ki is the permeability coefficient of the ith capillary and pci is the oxygen partial pres-sure within the ith capillary (assumed constant). The limit —s- oo represents an infinitely-permeable capillary wall for the ith capillary, or equivalently, that the oxygen partial pressure p at the boundary of the ith capillary is equal to the constant capillary pressure, pci, of that capillary. In contrast, the limit K{ —> 0 leads to the case of a perfectly-insulating capillary. This limit is physically-unreasonable since such capillaries would not contribute to the transport of oxygen in the tissue. In (4.2b) and (4.2c), d/dn is the directional derivative along the outward normal to the tissue domain. We incorporate skeletal muscle tissue heterogeneities, such as oxygen-consuming mitochondria, through the oxygen diffusivity V in (4.2a) and (4.2b) and the consumption rate of oxygen A4 in (4.2a) in the tissue. In homogeneous tissue models, both V and A4 are assumed constant, which one can interpret as meaning that the mitochondria are regularly distributed throughout the tissue. The assumption that the rate of irreversible oxygen consumption in tissue is constant, i.e. M. = Mo, is known as zero-order kinetics. First-order kinetics involve a linear dependence of the consumption rate on the oxygen partial pressure, M. = M.{p) = cp for some constant c > 0. A piecewise-linear expression for such an oxygen consumption rate, known as mixed zero- and first-order kinetics, is for some positive constant c. The mixed-order kinetics is a piecewise-linear approximation to Michaelis-Menten kinetics, in which the oxygen consumption rate has the form Here, pM is the value at which the consumption rate M is half its maximum. For small values of oxygen partial pressure p, the Michaelis-Menten kinetics approach first-order kinetics, whereas for large p, they approach zero-order kinetics. (4.5) M = Mop (4.6) P + PM 38 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue We can also model the local oxygen consumption by mitochondria with a spatially-dependent M. For example, a particular model for the oxygen consumption rate in a cell with m mito-chondria is m / I /* I 2 \ M = M0 + Y , M < E X P [ - J )• (4-7) Here, we model each mitochondrion by a Gaussian distribution superimposed on a constant background consumption rate Mo- In this way, there is freedom for varying the form of the heterogeneous tissue through the parameters, which are: the location, the amplitude, M i (where M i > Mo); and the variance, cr; (for o~{ reasonably small), of the i th mitochondrion for i = 1,.. . , 7 n . One could also combine this Gaussian distribution form of the mitochondria with Michaelis-Menten kinetics, as in K A Mop ^ Mjp ( |x - d M = h > exp 5— P + PM ' f-f P + PM ' V cr'f The oxygen diffusivity V within the skeletal muscle tissue can vary with location, whether in cells or in the extracellular phase of the tissue, and within the cells, whether at discrete oxygen-consuming mitochondria. Certain theoretical models consider that the tissue is a two-phase medium (eg. [7] [52]), having a constant specified diffusivity within each phase. We incorporate this approach in defining 'P(x) to have a similar Gaussian distribution form as in where VQ is the diffusivity outside of the m mitochondria and V{ is the diffusivity within the zth oxygen-consuming mitochondrial region. The fact that myoglobin is able to bind reversibly to one oxygen molecule can enhance, or facilitate, the diffusion of oxygen into the tissue. To incorporate myoglobin facilitation into our model, we follow Fletcher [11] and define the myoglobin-facilitated pressure p* as p* = p + pF[s(pc)-s(?)}. (4.8) Here, pc is an average of the capillary oxygen partial pressures, and pF is the facilitation pressure, which is constant for constant diffusivity V and is zero if no myoglobin is present in the tissue. When oxygen and myoglobin are in equilibrium, the myoglobin saturation s(p) has the form - - P (4.9) P0.5 + P where po.s is the myoglobin half-saturation pressure. To obtain the myoglobin-facilitated oxygen partial pressure, p*(x; e), we solve for p in the absence of myoglobin from (4.2) and substitute into (4.8). There have been recent efforts in extending the single capillary Krogh model to a multiple capillary system (eg. [41], [5], and [18]), and of those, Clark et al. [5] and Hoofd [18] include myoglobin facilitation in their models. With greater freedom in the physical parameters, our 39 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue multiple capillary model aims to eliminate certain restrictions of previous multiple capillary models, such as circular cross-sectional shapes, periodicity in the capillary bed, and homoge-neous tissue properties. We will find an approximate solution to our model using the hybrid method, which is based on a systematic asymptotic analysis that provides a measure of the error in the approximation. As well, the hybrid method is computationally faster and has a smaller physical parameter space than a direct numerical method. In the next section, we consider a simple version of (4.2), having one concentric capillary of radius e in a circular tissue domain of radius 1, with a constant diffusivity V — 1. We note that the annular capillary-tissue geometry is that of the original Krogh model in Figure 4.2, with R — 1 and Rc = e. The solution to the one-capillary problem contains a finite number of terms, and not an infinite logarithmic expansion in £. Nevertheless, the basic form of the gauge functions in the asymptotic solution is the same, which gives us insight into the structure of the solution of the problem with N > 1 capillaries, which does involve an infinite logarithmic expansion. 4.2 One Circular Capillary in a Circular Tissue Domain We reveal the structure of the asymptotic solution to (4.2) by considering the simple model problem A p = M, £ < r < l , (4,10a) e^--K(p-Pc) = 0, r = e, (4.10b) dp „ = 0, r = 1. (4.10c) or Here, A4 > 0, K > 0 and pc > 0 are constants. We remark that the capillary boundary condition gives a negative radial pressure gradient so that the capillary acts as an oxygen source. In this special case, we compare the asymptotic approximation to (4.10) with the exact solution, pE{r)i which is M PE{r) = PC + Y r2 - e2 e2 - 1 + l o g ( -This solution tends to negative infinity as K —> 0 in (4.11). This is to be expected from (4.10b) since in this limit, the capillary wall becomes perfectly insulating. From (4.11), we can see that the dominant term as £ —> 0 is (9(log£). The dominant term in the exact solution suggests that, in the global (outer) region away from the capillary where r £, we expand the asymptotic solution to (4.10) in the form p(r;e) = (log e)po(»0 + P i (>) + ••• • (4.12) Substituting (4.12) into (4.10a) and (4.10c), we find that the global region problems to solve 40 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue are A p 0 = 0, £ < r < 1, (4.13a) T^T = 0, r = l, (4.13b) and A p i = M, e < r < 1, (4.14a) 7 J 7 = 0, r = 1, (4.14b) such that pi(r) is singular and po(r) is regular as r 0. From (4.13) and requiring that po(r) is bounded as r —> 0, we have that Po(r) = po = constant. (4.15) In the local (inner) region near the capillary, we introduce the local variables p = r/e and q(p;e) = p(ep;e). These are the radial version of the local variables in (2.5). Substituting the local variables into (4.10a) and (4.10b), to leading order, we get that q satisfies A p q = 0, p>l, (4.16a) da ^ - « ( ? - P c ) = 0, p=l. (4.16b) Here, A p indicates that the Laplacian operator is with respect to the local radial variable p. To match between the local and global region solutions, we require also that q grow logarithmically as p —T oo. We expand the local region solution as w here qo(p) satisfies The solution to (4.18) is q{p;e)=pc + q0(p) + o(l), (4.17) Apq0 = 0, p>l, (4.18a) ^ - K 9 O = 0, p = l, (4.18b) q0~a0logp, p-^oo. (4.18c) qo(p) = aologp + — • (4.19) Using (4.18c) in (4.17) as p —> oo, and comparing to (4.12) as r —>• 0, we find that po = —ao and that Pi ~ a 0 l o g r + pc + —, as r -> 0. (4.20) 41 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue We need to solve (4.14) together with (4.20) for pi. Or, equivalently, we can solve Ap1 = M + 2wa06(r), \r\ < 1, 97 = °' r = L (4.21a) (4,21b) In (4.21a), S(r) is the Dirac delta function. We determine the value of a0 by integrating (4.21a) over the region D, which is a circle of radius 1, and applying the Divergence theorem with the use of (4.21b). We get that MTT = M J dA= -2i:a0. D (4.22) Thus, a 0 = - M / 2 and so from (4.20), M, M Pi ~ - - y log r + pc - — , as r 0. The solution to (4.14) together with (4.23) is Mr2 M , M Pi = — + P c - T l o g r - - . To recapitulate, the asymptotic solution is , , M p(r;e) = pc + — r2 1 , /£ (4.23) (4.24) (4.25) Comparing (4.25) to (4.11), we see that the asymptotic solution agrees with the exact solution up to the 0(1) terms and is missing the 0(e2) correction. From the solution, we can see the role of the parameter K, the permeability coefficient of the cap-illary wall. As K —*• oo, we obtain the asymptotic solution for the Dirichlet boundary condition at the capillary, p — pc. As K —>• 0, the physically-unreasonable case of a perfectly-insulating capillary wall, we notice that the solution becomes increasingly negative. In Section 4.5, we determine an estimate for the minimum time necessary tcr for the diffusion process to reach steady state. In this time scale estimate, tcr —> oo as K —• 0, meaning that the steady-state condition would take forever to occur. The asymptotics of this simple model problem provide us with the approach for attacking the more complicated general problem in (4.2). In the next two sections, we consider multiple cap-illaries in an arbitrarily shaped tissue domain. First, we use the method of matched asymptotic expansions to show that for N > 1 capillaries, the asymptotic solution contains an infinite logarithmic expansion in e, the measure of the magnitude of the capillary cross-sections. Sec-ond, we use the hybrid method on the multiple capillary problem to sum essentially the entire logarithmic series in the asymptotic solution. 42 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue 4.3 Solution using the Method of Matched Asymptotic Expan-sions In this section, we obtain the asymptotic solution to (4.2) for the special case of N capillaries with equal, circular cross-sections of radius e. In the next section, we will apply the hybrid method to (4.2) with arbitrarily shaped capillary cross-sections. As in the simple model problem for one capillary of the previous section, we expand the global (outer) region solution in the form CO p(x;£) = ^- 1(e)p ( 0 )(x) + p(1)(x) + £ ^ - 1 ( £ ) p W ( x ) + -- - , (4.26) where we define v(e) = — 1/loge. In the local (inner) region, near the ith capillary located at x = we define the local variables x - £ y.- = — ^ - L , » (y ; ; e) = Key,- + £ ) - ( 4 - 2 7 ) The local variables in (4.27) are the extension of those in (2.5) for multiple subdomains located at , for i — 1,.. . , N. We expand the solution in the local region near the ith capillary as ft(y<; e) = Pd + ql1\yl) + "(e)q?\yi) + •••. (4.28) Here, pci is a specified constant that represents the oxygen partial pressure within the ith capillary, for i = 1,... ,N. Substituting the local variables into (4.2a)-(4.2b) and using V(eyi + £,•) ~ *P(£i) + 0(e), we Yi\ > 1, (4.29a) y i | = 1. (4.29b) To allow for matching to the global region solution, it is necessary for each qi to grow logarith-mically as |y;| —• oo, for i = 1,... , N. Substituting (4.28) into (4.29), we find that all the local region problems for qf^ are the same, so (?) (?) (c) (?) we can write q\ — a\ 'q\ for j > 1. Here, a\ are unknown constants that we will determine through the matching process, and q\C\yi) is the canonical local solution satisfying = 0» M > ! . (4-30a) da (c) i \ m ) d n ~ + Kiq< = °' | y i 1 = X' ( 4 ' 3 0 b ) 9 J c ) ~ l o g | y i | + ^ , | y 4 | - o o . (4.30c) find that, to leading order in e, qi(yf,s) satisfies Ayqi = 0, 'P(ti)-^1 + Ki(qi - Pci) = 0, 43 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue This is the multiple hole version of the canonical local problem that we first introduced in (2.10). In the global region, far from the capillaries, we substitute (4.26) into (4.2a) and (4.2c) to obtain the problems for p^\x.) for j > 0. The problem to solve for p(° '(x) is V • [PVpW] = 0, x e D\ U {&}, (4.31a) dpW -j— = 0, x G 3D. (4.31b) As well, we require that p( 0 ' is regular as x —> for i — 1,... , N. The solution to this problem is that p(°) is a constant. We will show later how to determine this constant. The remaining global region problems for p( J)(x), j > 1, are V-[PVpW] = ejlM(x), x e £ \ U { & } , (4.32a) i=i ^ — = 0, x e 3D. (4.32b) on In (4.32a), eik is the Kronecker delta function defined as elk = { l ~ k . (4.33) k \ 0, otherwise. v ' As well as satisfying (4.32), p^(x) must be singular as x —> for i = 1,.. . , N. In anticipation of matching the local and global region expansions, we write the global region expansion a s x - > £ i ? for i = 1,... , N, and the local region expansion as |y;| —> oo. As x — £ j , for i = 1,.. . , N, the global region expansion has the form CO p(x - fc; e) ~ + P ( 1 ) (x -» &) + v]-\e)v (3){* - Ci) + • • • • (4.34) Using (4.30c), we find that as |y,-| —» oo, the ith local region expansion, in terms of global region variables, is of the form 9i(y,-; e) ~Pct + af '{ log |x - &| + v-l{e) + P(C,-)/«i} + K £ ) a ! 2 ) { l o g |x - £.| 4- „ - ! (£ ) + 7>(&)/K,-} + • • • . (4.35) Once again, u(e) = —1/ loge; the constants p c;, V(£j) and K J , for i = 1,.. . , A , are known; and aa- are unknowns that we will determine through matching to the global region solution. Comparing the 0(v~v) terms in (4.34) and (4.35), we see that matching requires a ( i ) = 4 D = ... = o ( J ) = p ( 0 ) . (4.36) As well, the matching procedure provides the behaviour of the global region solution p^\x) for j > 1, which is ~ (log |x - £,| + ^ ) + a p 1 ' + e j l P c i + •••, x -> i = 1,.. . , N. (4.37) 44 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Using (4.32) and (4.37), we can write the global region problems for p(J')(x), j > 1, as V - [ P V p W ] = e i T A4(x) , x e / J \ . U { 6 } , (4.38a) dpti) = 0, x G 3D, (4,38b) ~ aj j )(log |x - fc| + V^i)/Ki) + a J J + 1 ) + e j l P c i , x - i = 1 , . . . , N. (4.38c) Combining (4.38a) with (4.38c), we can write N V-[PVp<'>] = e , i ^ W + 2 7 r ^ a t a V ( ^ ) ^ ( x - ^ ) , x E D . (4.39) Here, <5(x — £ J is the Dirac delta function. The solution to (4.39) together with (4.38b) is not unique. To address this point, we consider the corresponding eigenvalue problem V • [PVfa] = -Xk<f>k, x£D, (4.40a) d<j>k On = 0, x G 3D. (4.40b) The first eigenpair for this problem is (0 O ,A O ) = ( 1 , 0 ) , and the remaining eigenvalues, A&, are such that 0 < Ai < A2 < . . . < Afc < oo. We expand the solution p^ ' (x) , for j > 1, as a sum of the eigenfunctions CO . , ^ A/t < 4>k, 9k > Here, b3 = enM(x) + 2ir J^Li ^ ' ^ ( ^ ( x - & ) i s t h e right-hand side of (4.39) and < f,g > is the inner product of the functions / and g defined by < / , g >= JD f(x)g(x) dx. In (4.41), since A 0 = 0, we require that < bj,(f>0 > = < bj,l >= 0 for a solution to exist, and. so the solvability condition for pW(x), for J > 1, is x—> m , •, f — — f A f ( x ) d x , 7 = 1 Eo?)p(w= 27ri^ (4-42) i = i [ 0, j > 2. If (4.42) is satisfied, then we can write the solution p^\x), j > 1, as CO for any constant C j . We need to fix the constant Cj in order to obtain a unique solution p(-»')(x). Requiring that < p^\(pQ >= 0, where <J>Q = 1, is equivalent to < b;,<f>k > jp^\x)dx=Q. (4.44) 45 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Using < pti\<t>o >= 0 in (4.43) and applying the orthogonality condition of the eigenfunctions, we get that Cj = 0 for j > 1. Thus we can decompose the solution p^\x), for j > 1, as pW(x) =p(f)(x) + f4 j ). (4.45) In this decomposition, a^' for j > 1 are global constants, and pl J ' (x) is the unique solution to V • [PVpW] = ejlM(K), x £ j ) \ U {£•}, (4.46a) Q7p) = 0, x 6 3 D , (4.46b) C m 2#> - aJ j ) ( log |x - £ | + V(ti)/Ki) + ap ' + 1 ) + ejlPci, x - i = 1,... , N, (4.46c) j p (£ )(x)dx = 0. (4.46d) In (4.45) and (4.46c), the unknowns are a^' and of+1^ for i = 1,... , JV, which we will determine at each level j for j > 1. We decompose the unique solution p{P(x) as N j#>(x) = 27T ] T a p V ( ^ ) C ( x ; £ ) + e ^ M - (4.47) i=i In (4.47), G(x; £) is the modified Green's function satisfying 1 V . p > V G ] = - - L x G D \ { C } , (4.48a) dc - 7 — = 0, x G 3 D , (4.48b) an G ~ 2^(0 1 O S | X~ C I + jR (C )' X " C ' (4'48C) y G(x; 0 rfx = 0. (4.48d) D In (4.48a), A D is the area of the region D and in (4.48c), R(C) is the regular part of the modified Green's function. We determine R(C) uniquely from the solution to (4.48). Also in (4.47), Pft\x) is a regular function, as x —> z = 1,... , JV, that satisfies V •[VVp (k )] = M{x)-^- [ M(x)dx, x G D , (4.49a) A-D J D (1) = 0, x G 3 D , (4.49b) Or Jp (k\x)dx= 0. (4.49c) D 46 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue In (4.49a), we used (4.42) in the right-hand side. Substituting (4.47) into (4.45) and using (4.48c), we obtain an expression for the behaviour of p(J')(x) as x —* £j which is ~ a j j ) log |x - & | + 27raS JV(fc)ii(&) + 2 7 r £ 4 ^ ( 6 ^ ; + w k 1 ^ ) + 4 J ) + o(l). (4.50) k=l kjti Comparing (4.38c) and (4.50), we find that the log |x — £J terms automatically agree. From the remaining terms, we see that N 2*a (i j)V(ti)R(ti) + 2TT £ "PntkMti; tk) + ejlP^\d) + = k=i + eHPd + oiJ)^ (^-)/«.-. (4-51) for i = 1,. . . , JV and j > 1. To determine the constant p(°\ and hence the constants a j 1 ' for i = 1,. . . , JV, we use (4.42) with j = 1 and (4.36) to obtain / A*(x)dx p(°) = - - o . (4.52) 2 * E £ i P ( f c ) Thus, the constants aj1^ for i = 1,... ,7V are known, and (4.51) provides JV equations for the (JV + 1) unknowns: a^) and a - J + 1 ' for i — 1,... , JV and j > 1. The solvability condition in (4.42) for flj-J+1^ provides the last equation for these unknowns. In summary, with the known, we solve (4.51) and (4.42) recursively to find a\ J\ for i = 1,.. . , JV, and for each j > 2. The parameters to specify in the problem are £ order of magnitude of the capillary cross-sections JV number of capillaries pci oxygen partial pressure within the i th capillary K i permeability coefficient of the ith capillary wall location of the ith capillary ^ ' "P(x) diffusivity in the tissue A4(x) consumption rate of oxygen in the tissue D geometry of the tissue domain. For a given number and location of the capillaries in the tissue domain, we use (4.48) once only to compute the modified Green's function G(x ;£) and its regular part R(Q- We use (4.52) to compute the constant p( 0 ' and (4.49) to compute the regular function p^ ' (x ) . Then, we evaluate these expressions at the capillary locations In general, we must perform these computations numerically. For certain special cases, such as a circular tissue domain D with one concentric capillary, it is possible to solve for the modified Green's function G and the function p Q analytically (see Section 4.6). 47 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Thus, we have an infinite set of recursive problems for the global region function coefficients a t-7'. To obtain the global solution for p(x;e) in (4.26) for a given e, the order of magnitude of the capillary cross-sections, we would use the expressions in (4.45), (4.47) and (4.52), which require these coefficients for j > 1. In the next section, we examine the solution to the problem via the hybrid method which avoids the daunting task of determining each set of coefficients at each j'-level. 4.4 Solution using the Hybrid Method We now apply the hybrid method to solving (4.2) with N capillaries of arbitrary cross-sectional shape. Using (4.28) together with q\^ = q\c\ we expand the solution in the local region, near the zth capillary, as ft(y*;e) = Pet + Q\0)(yuVi,... ,uN) + £ O j 1 }(y ; ; ' A , . . . , vN) + • • • . (4.54) Here, we define V{ = — 1/log(£c(;) where di depends on the shape of the cross-section of the ith capillary. We take Q\0^ — A,g- C ' (y) , where A; = Ai(v\,... , 7>>/v), for i — 1,... ,JV, are to be determined. Also, q\°\y) is the canonical local solution satisfying y, g* A 0 , (4.55a) y, G 8D°, (4.55b) |y;| -> oo. (4.55c) Here, di is the shape-dependent parameter that is also known as the logarithmic capacity (see Garabedian [14]), and is the scaled ith capillary, such that Df = eD®, and 8D® is its boundary. This canonical local problem is a multiple subdomain version of (2.10). For the case of capillaries with circular cross-sectional shape of radius e from the previous section, A{ is asymptotic to the infinite logarithmic expansion (4.56) This expression is the multiple subdomain version of A{ed) in (2.19). We note that each A{ ~ 0(1) as £ —> 0 for i = 1,... , N. We recall that in the previous section the constants di are all the same and are unknown coefficients that are determined through matching for each j > 1. For capillaries of arbitrary cross-sectional shape, the d{ are not all the same but each A{ = Ai(v\,... , VN), for i = 1,... , A , sums a similar infinite logarithmic series. The hybrid method will exploit the asymptotic structure of the solution to sum essentially the entire logarithmic series. For a specific cross-sectional shape of the i th capillary, we compute the shape-dependent pa-rameter di uniquely from the solution to (4.55). In some cases, like for a circular or elliptic Ayq\c) = 0, (lf] ~log|y,- | - logd , - - - - , 48 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue cross-sectional shape, it is possible to determine d{ analytically (see Section 2.5). For a circular capillary of radius e, di = e x p ( - - P ( £ J ) / K ; ) . For the particular case of K,- = oo, we find di from the table of capacities from Ransford [45], which we reproduced in Table 2.1 in Chapter 2. Using QJ 0 ' = Aitif* and (4.55c) in (4.54), we have that the far-field behaviour of the i th local region solution, in terms of global region variables, is qiiYi^) ~ Pd + A;[log |x - & | + v-^edi)] + • • • , |y;| oo. (4.57) In the global region, we expand the solution as p (x ; £ ) = 7 G^ + P ( ° ) ( x ; , ; / 7 V ) + £ P ( 1 ) ( x ; ^ , . . . ,zv w )+ . . - . (4.58) Here, pG is a global constant that we will determine along with A ; , for i = 1,.. . , N, and again, Vi = -1/log(edj). We will see that pG ~ O(logs) as £ —> 0. Substituting (4.58) into (4.2a) and (4.2c), and requiring that the global region solution match to the local region solution, we get that P(° ) (x ; v\,... , VN) is the unique solution to V • [PVPW] = Af (x), x G D\ U (4.59a) i=l —z— = 0, x e dD, (4.59b) On P ( ° ) ~ Ailog\x-tl\ + A t v ^ + P c i - p G , x ^ ^ , i = l,...,N, (4.59c) y P ( 0 ) ( x ) r / x = 0. (4.59d) D There is a solution to (4.59) provided that N E A^') = - ^ / M(x)dx. (4.60) i = 1 D This condition is analogous to the solvability condition (4.42) of the previous section. We write the solution P(°) (x ; v\,... , VN) in the form N P(°) (x ; vt,... , vN) = 2TT Y M(yx, • • • , uN)V^)G^ fc) + P £ 0 ) ( X ) . (4.61) »=i Here, G' is the modified Green's function that satisfies (4.48). As well, P^° ' (x ) is a regular function, as x -> for i = 1,... , A , that satisfies (4.49). Using (4.48c) in (4.61), we obtain that the near-field behaviour of P ^ 0 ' is N P ( ° ) ~ Ailog |x - & | + 2 7 r A 1 P ( & ) P ( & ) + 2 7 r AkV(tk)G(ti-, £k) + P £ 0 ) ( £ ) , k^i x - & , i = l , . . . , A . (4.62) 49 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Comparing (4.59c) and (4.62), we see that the log |x - £J terms automatically agree. From the remaining terms, we obtain a set of N equations N 2wAlV(^)R(L) + 2nJ2 AkP^k)G(^; tk) + PR°\^) = A^f1 + pcl - pG. (4.63) k = l k^i Thus, (4.60) together with (4.63) provide (N + 1) conditions for the unknowns: pG and Ai, for i = l,...,N. We let e — 0 in (4.63) and find that, for a non-trivial solution for A,- in (4.63) to exist, we require that pG balance with A{V~X. This means that pG ~ O(log£) as e —» 0 since A,- ~ 0(1) as £ —> 0. If we neglect the off-diagonal terms in (4.63), that represent capillary interaction, then we obtain an expression for pG from the diagonal entries, 1 N (Pa)diag = (pd + Aiur1 - 2KAlV{Zi)R{Zi) - PR°XZi)) • (4.64) i=i In summary, (4.60) and (4.63) provide (JV + 1) equations for the unknowns pG and A i , for i = 1,.. . ,N. As in Section 4.3, we must specify the parameters in (4.53), as well as the cross-sectional shapes of the A r capillaries. For a given tissue geometry for which we specify the locations and the number of capillaries, we compute the modified Green's function G'(x; C) and its regular part R(C) once only from (4.48), and compute the function PJ^(x) from (4.49); and evaluate these expressions at the capillary locations For each specified cross-sectional shape and permeability coefficient K{ of the ith capillary, we determine the unique shape-dependent parameter di from (4.55), for i = 1,... , N. 4.5 Time Estimate to Reach Steady State Our goal here is to obtain an asymptotic estimate when £ <C 1 for the length of time necessary for the oxygen diffusion process to attain steady state. For simplicity, we consider only one capillary in the tissue domain. We consider the unsteady version of (4.2) with N = 1 for p(x,rj;e), in which we replace (4.2a) with ^ = y.[Wp]-M, xeD\D£. (4.65) We also include an initial condition, p(x, 0;e) = Po(x). We assume that the diffusivity V depends only on the spatial variable x, but we allow the oxygen consumption rate M. to be a function of the oxygen partial pressure, p(x,£;£) , such that dM/dp > 0. We will compare the steady-state time estimates for the case with M = A4(x), independent of p, versus the case with M = M{p). We substitute p(x, t; e) = p a t (x ; e) + e~Xtv(x; e) (4,66) 50 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue into (4.65), where v <C 1 and pst(x;e) satisfies the steady-state problem (4.2) with N = 1. Linearizing (4.65), with boundary conditions (4.2b) and (4.2c), we obtain that v(x; e) satisfies the eigenvalue problem V-[PVv]-M'(pst)v = -\v, x£D\D£, (4.67a) dn S V ^ - + KV = 0, xedD£, (4,67b) dv 0^ = 0, xedD, (4.67c) together with the normalization condition v2dx = 1. (4.68) D\D£ There are infinitely many eigenvalues A of (4.67). We are interested in an estimate for the smallest such eigenvalue, \Q{S) as e -> 0, which will provide us with the time scale necessary to reach steady-state conditions. First, let us consider the case with M = A4(x), independent of p. In this case, (4,67a) reduces to V • [Wv] = -Xv, x e D\D£. (4.69) We will determine the first eigenpair, v^°\x;e) and A(°'(£), of (4.67) with (4.67a) replaced by (4.69), as £ —> 0. Since the first eigenvalue for the unperturbed problem (i.e., with no capillary) is AQ°' = 0, we expect that the first eigenvalue for the corresponding perturbed problem will tend to zero as £ —> 0. Based on this, we expand v^°\x;s) and the smallest eigenvalue \(°\e) as v(°\x;e) ~ v (0°\x) + v(ed)v[°\x) + ••• , (4.70a) A<°)(£) ~ A 0 0 ) + v(ed)\ {° ) + . . - , (4,70b) as £ —> 0. Here, u(ed) is some unknown gauge function. We will construct global region and local region expansions for the eigenfunction v^°\x; e), corresponding to regions away from and near to the capillary respectively. Through the procedure of matching the solutions in the global and local regions, we will show that the gauge function, u(ed), has the reciprocal logarithmic form of (2.18). Our interest in constructing the matched asymptotic solution lies in obtaining A^ 0 ' , which we will use in (4.66) to calculate the steady-state time estimate. In the global region, away from the capillary, we substitute (4.70) into (4.69), (4.67c) and (4.68) and obtain that the global region functions v^°\ for j > 0, satisfy V • [PVvf] = -J2 A S ° - W ° \ x G D\{£0}, (4,71a) «=o 5 r ( o ) dn 0, x G dD, (4.71b) 4 0 M 0 ) dx=\1, * = J = 0 ' (4.71c) 0, otherwise. D 51 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue In the local region, near the capillary located at x = £, we define the variables y = £ _ 1 ( x - C) and w(°\y;e) = v(°\ey + C;e). We expand the local region solution w(°\y;e) as w ( 0 ) (y; e) = a0is(ed)wc(y) + • • •, (4.72) where a0 is a constant that we will determine through the matching procedure. We substitute (4.72) and the local variables into (4.69) and (4.67b) and find that wc(y) satisfies the canonical local problem in (4.55), to leading order in e. For the matching procedure, we compare the global region expansion of v(°) in (4.70a) as x —> £ with the local region expansion of in (4.72) as |y | —*• oo, using the far-field behaviour from (4.55c) in terms of the global region variables. The matching procedure requires that the gauge function is v{ed) = — l/log(ed) as in (2.18), and that v^0) = a 0 , (4.73a) i ; i 0 ) ( x ) ~ ^ 0 ) l o g | x - C | + - - - , x - C - (4.73b) Combining (4.71a) for j = 1 with (4.73b), and given that A{, 0 ) = 0, we can write V • [TS7v[0)] = -A< 0 ) t ;J 0 ) + 2TTV(0)V{C)S(X - C), X £ D. (4.74) Here, 6(x — Q is the Dirac delta function. Using (4.71b) and Green's identity, we obtain < Lv[°\v^0) >=< Lv(°\v{0) >, where L is the operator L = V • [PV]. Since Lv^0) = 0, and with Lv[°^ given above, we obtain A f = ^ . (4-75) Here, A D is the area of the domain D. We determine the value of the constant (and hence ao) from the normalization condition in (4.71c) with i = j = 0, which gives us that 40) = ( ^ ) - 1 / 2 - (4-76) Now, we construct the estimate for the time necessary for the process to reach steady state. The oxygen partial pressure is approximately the steady-state value, p(x,t) « p st(x) for t > (Af 0 ) )" 1 . Using v(ed) = -l/\og(ed) and (4.75) in (4.70b), with A^ 0 ) = 0, we can express the minimum time required to reach steady state as log(ed)AD 1 > - ~ t v k f { " T ) From this expression, we see that the time estimate to reach steady state is O(loge). This gives a slow decay to steady state, yet one that is more rapid than an 0{e~1) decay. Next, let us consider the case with M = M(p). The first eigenvalue, Ap°^, of the unperturbed problem of (4.67) is dx D 52 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue For M independent of p, the first eigenvalue A(°) is 0 ( —1/log(£eQ), which is smaller than Ap 0 ' , the first eigenvalue for the case with M — M(p), whose leading order is in (4.78). In (4.66), we see that for a larger first eigenvalue A ' 0 ) , there is a faster decay to steady state. Thus, the process in the case with M = M(p) decays to steady state more rapidly than in the case with M = .A4(x), independent of p. Example: Krogh Cylinder Geometry. We consider the special case of a circular cross-section of tissue of radius 1 containing a concentric circular capillary of radius £, with uniform diffusivity V = 1 in the tissue. For a circular capillary with constant diffusivity, A D = w and the solution to (4.55) gives that the shape-dependent parameter is (1(K) = exp( —1/K ) . Substituting these expressions into (4.77), we obtain that i » i ^ - l o g e + ^ , £ < 1 . (4.79) In the next section, we provide details of the numerical procedure for finding the modified Green's function solution to (4.48). 4.6 Finding the Modified Green's Function For certain special cases, it is possible to solve analytically (4.48) for the modified Green's function. One such special case is that of a circular cross-section of tissue having radius a, with constant diffusivity V = Vo in the tissue, and with a concentric capillary (singularity) located at the origin. In this case, the solution to (4.48) is ^ ) = 2 ^ 1 o g r + I ^ Q - 2 1 o g a - ^ ) . (4.80) The regular part evaluated at the capillary centre is ^ = * k i l - 2 ] c s a ) - (4'81) When it is not possible to obtain an analytic solution, we discuss how to numerically solve (4.48) for the modified Green's function, G (x ;£ ) , and its regular part, R(C)- We used the Partial Differential Equation Toolbox [36], a finite element code for solving for w(x) in elliptic partial differential equations of the form - V • [cVu] + au = f, x <E D, (4.82a) On c— + qu = g, x e dD. (4.82b) on We employ a regularization procedure to compute G in (4.48) by introducing a regularization parameter 6. To implement the guarantee of uniqueness in (4.48d), we solve (4.48) for G$(x; £, 6) with (4.48b) replaced by ^ + SGs = 0, xedD, (4.83) on 53 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue which has a unique solution for non-zero 8, where 8 < 1. We impose the singularity structure from (4.48c) by seeking a solution of the form Gs(x;C,S) = 2itV{Q l og |x -C | + i?5(x;C,<5). (4.84) Substituting (4.84) into (4.48), with (4.48b) replaced by (4.83), we obtain that Rs(x; £,8) is a regular function as x —• C that satisfies V-[VVRs} = -~, xeD, dRs m + SR5 = 1 d log |x - C| •logjx- C| L 2TTV(() dn "° > ] 2wV(C) J Rs(x; C) dx = - 2 7 r ^ ( f l y lo§ l x - CI dx. x e dD, (4.85a) (4.85b) (4.85c) We expand the solutions to (4.48), with (4.48b) replaced by (4.83), and (4.85) as 8 0. To determine the form of these expansions, we consider the corresponding eigenvalue problem for the special case of a circular domain D of radius a, with constant diffusivity V = 1 and with the singularity (capillary) located at the origin Afa + \s4>s = 0, 0 < r < a, (4.86a) ^ + 6<f>s = 0, r = a. (4.86b) Since A 0 = 0 and (f>0 — C, a constant, form a solution to (4.86), we expand the principal eigenvalue and the corresponding eigenfunction as Xos(8) ~ 8Xt + 62X2 + ••• , (4.87a) 4>os(r; 6)~C + 6<h(r) + •••, 6^0. (4.87b) Substituting (4.87) into (4.86), and imposing a solvability condition, we find that </3i(r) satisfies A4>i = -CXi, 0 < r < a, (4.88a) ^ 1 = -C, r = a. (4.88b) or The solution to (4.88) is ^ ( r ) = -Cr2/(2a), providing that X1 = 2/a. Thus, from (4.87a), we can write 28 Xos(S) ~ — + • • • • (4-89) We express the solution, R$(r), as a sum of eigenfunctions of the homogeneous operator in (4.85a). Then, using the orthogonality of eigenfunctions and Green's second identity, we get f^0<<l>js,<i>js> f^0 x]5<4>38,(f>38> 54 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Here, < u, v > is the integral inner product, and / is the right-hand side of (4.85a) and g is the right-hand side of (4.85b). Since A 0s ~ 0(6), and the leading-order eigenfunction is a constant C, the solution for Rs(r) would be 0(6~J) providing that the numerator does not vanish as 5 —> 0. For a circular domain of radius a, with constant diffusivity, the numerator in (4.90) is S(-Caloga). (4.91) Thus, we can see that for a circular domain of radius a, the numerator does vanish as 6 —> 0. Since the numerator vanishes, we expand the solutions for 6 —>• 0 as G 5(x; C, S) ~ G 0(x; C) + SG,(x; Q + • • • , (4,92a) Rs(x; C, 6) ~ R0(x; C) + W i ( x ; C) + • • • • (4.92b) We compute the regular solution to (4.85a)-(4.85b) for two values of 6 and we call the resulting regular solutions (Rs^c and (Rs2)c- The size of 6 we choose is bounded below by the tolerance 6 for our computed solution: e.g. 5 > 6cr = 0.001. We expand the computed solutions as (R6l)c~ R0c + 61Rlc+-- - , (4.93a) (Rs2)c ~ Roc + hRic + • • • • (4,93b) From an. extrapolation of (4.93), we obtain the solution ^ ^ ^ J c - W c , ( 4 . 9 4 ) h - oi This is the leading-order solution in the expansion for i2,$(x; £, 6) in (4.92b), although it is not unique. We still need to impose (4.85c) to determine the arbitrary additive constant. Substituting Rs = Roc + k into (4.85c), we obtain the constant k and thus can express the solution as Rs(x;C,5) = R0c- -y-° Example: Krogh cylinder geometry. We test this computation technique on a special case for which we know the exact solution: a circular tissue domain of radius 1, with diffusivity V = 1 and containing a single concentric capillary. For this special case, the exact solution to (4.48) for the modified Green's function is the purely-radial solution in (4.80) with a = 1 and VQ — 1, which simplifies to < ? M = i l o g r + i Q ^ ) . (4.96) We compare the regular part of the exact solution in (4.96) as r —• 0 to what we compute from (4.48), with (4.48b) replaced with (4.83). We note that, since the singularity (capillary) is located at the origin and since the outer boundary is at r = 1, the log |x — £| term vanishes 2irV(C) log |x — £| dx • Roc dx (4.95) D n 55 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Maximum absolute error L2 of error Computed R(0) {0.025,0.05} 8.42e-04 0.0017 0.1195 {0.01,0.025}, 8.26e-04 0.0016 0.1194 {0.01,0.05}, 8.30e-04 0.0016 0.1194 Table 4.1: Test of (5-dependence on R$(0) using three pairs of 8 values: maximum absolute error and L2 norm of the error between Rs(0) in (4.95) and the exact R(0) in (4.97). in (4.85b) and so the exact solution to the problem for the regular part, R$, is independent of 8. The regular part of the modified Green's function evaluated at the capillary location is 3 R(0) = — « 0.1194. (4.97) We test the dependence on 8 of our computed solution (Rs)c for three pairs of 8\ and 82. Table 4.1 contains the maximum absolute error, the Z 2 -norm of the error and the computed value of R(0) for each of the <5-pairs, on a computational mesh containing 4,088 elements. The values in the table indicate that that the computed solution is essentially independent of 8. In the final section of this chapter, we display the results of the hybrid solution for the oxygen transport application. 4.7 Computed Results and Discussion We present our results demonstrating the effect of homogeneous versus heterogeneous tissue with respect to oxygen diffusivity and consumption rate, and the effects of capillary inter-action, the cross-sectional shape of the capillaries, variable capillary wall permeability, and myoglobin-facilitated oxygen transport. Although we can consider tissue domains and capillar-ies of arbitrary shape in our model, we choose simple geometries in the examples to highlight certain effects that are unrelated to the tissue geometry. In certain cases, where possible, we compare to the exact solution or a direct numerical solution of (4.2). We have used Matlab, in particular the Partial Differential Equations Toolbox [36], for the simple numerics of the hybrid method and for the direct numerical computations. Unless otherwise specified, the discrete points in the plots below are the direct numerical results. In Table 4.2, we list ranges of dimensional parameter values for use in our computations that we obtained from a number of references. The dimensionless diffusivity V and dimensionless capillary pressure pc are measured relative to their dimensional counterparts (indicated with *). In the examples below, we choose values for our dimensionless parameters, such as the dimen-sionless oxygen consumption rate M. = L*2 M*/(p* T3*), that correspond to typical dimensional quantities. Example: Homogeneous versus heterogeneous tissue. For the special case of one circular capillary of radius e that is concentric with a circular tissue domain D of radius 1 with diffusivity V = 1, the exact solution of (4.2) for the oxygen partial pressure in the tissue is (4.11). In 56 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Dimensional Parameter Range of Values Units O2 consumption rate M* 3.3 X l f r 5 -> 3.1 x 10" 2 mL-02 /cm 3 -s Tissue diffusivity -p* 3.9 x 1CT 1 0 -> 4.8 x 1 (T 1 0 mL-02/cm-s-mm Hg Tissue O2 partial pressure P* 26 -* 30 mm Hg Tissue domain length scale L* 1.5 x 1 0 - 3 -> 1.8 x 1 0 - 2 cm Table 4.2: Range of dimensional parameter values; compiling data from Clark et al. [5], Ellsworth et al. [9], Hoofd [17], Hsu et al. [21] and Popel [42]. Exact Hybrid Numerical 4 Mitochondria 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Figure 4.5: Minimum oxygen partial pressure pmin versus e, the capillary radius, for a circular tissue domain of radius 1 with one concentric circular capillary. Parameter values: pc = 5, M = 0.5, K = 0 0 , V — 1. Mitochondria parameter values: Ci = (±0 .5 , ±0.5) , M o — 0.5, Mi = 50, V0 =l,Vi = 4, and cr; = 0.05, for i = 1,... , 4. Figure 4.5, we have plotted the minimum oxygen partial pressure pmin versus the capillary radius e to compare the hybrid method solution with the exact solution, and for certain values of e, with the direct numerical solution of (4.2). The (dimensionless) parameter values for this example are: capillary pressure pc = 5, oxygen consumption rate M = 0.5, and capillary permeability coefficient K = 00 (an infinitely-permeable capillary wall). In homogeneous tissue models, the oxygen diffusivity and oxygen consumption rate are constant. Also in Figure 4.5, we have included the corresponding pmin(e) curve for a heterogeneous tissue domain of the same geometry but containing 4 oxygen-consuming mitochondria. We incorporate the mitochondria using the Gaussian distribution form in (4.7) for both M(x.) and 7 3(x) with parameter values: mitochondria locations (^ = (±0 .5 , ±0.5); variance CT,- = 0.05; background oxygen consumption rate and diffusivity, M o = 0.5 and Vo = 1 (the same as for the homogeneous tissue); and oxygen consumption rate and diffusivity within the ith 57 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue mitochondria, Mi = 50 and Vi = 4, for i = 1,.. . ,4. The pmin{£) curves in Figure 4.5 are increasing functions of e, indicating that the larger the capillary radius, the more oxygen the tissue receives. The heterogeneous tissue p m ; n ( £ ) curve lies below that of the homogeneous tissue, given the same background oxygen consumption rate and diffusivity, showing that the presence of the mitochondria lowers the oxygenation in the tissue. The figure also reveals that the hybrid solution agrees well with the exact solution for values of £ up to approximately 0.2. Table 4.3 contains the time to compute a single point of the homogeneous tissue pmin(e) curve of Figure 4.5. The table shows that the direct numerical solution took approximately 13 times as long to compute as the hybrid method solution using meshes of comparable refinement. Pmm(0.05) Mesh size (elements) Time/point (seconds) Exact 4.3758 Hybrid 4.3760 7680 2.99 Numerical 4.3765 7424 40.49 4.3759 29696 187.88 Table 4,3: C P U time to compute one point, pmin(e) — p m i„ (0 .05) , for the homogeneous tissue case of Figure 4.5: hybrid method solution versus the direct numerical solution. Example: Capillary interaction. We consider JV = 4 capillaries of circular cross-section with radius e in a circular tissue domain D of radius 1, and vary the intercapillary spacing to display the effect of interaction on oxygenation of the tissue. Figure 4.6(a) shows the locations at £ j{ = ( ± J / 4 C O S ( T T / 4 ) , ±j74sin(7r/4)), for j = 1,2,3 of the i = 1,.. . ,4 capillaries. The parameter values for this example are: diffusivity V = 1 and oxygen consumption rate M = 0.3; capillary permeability coefficients K , - = oo; and intracapillary oxygen partial, pressure p c ; = 5, for i = 1,.. . ,4. With circular cross-sectional shapes of radius £ and with «;,• = oo, the shape-dependent parameters are cZ; = 1, for i — 1,... ,4. In the top graph of Figure 4.6(b), we have plotted the minimum oxygen partial pressure pmin versus s, the radius of the capillary cross-sections. For certain values of e, we have included the results from a direct numerical solution of (4.2). Of the -p m ; n (£ ) curves that include interaction effects, the j = 2 curve lies uppermost, showing that the tissue receives more oxygen when the capillaries are in this position, than in the other capillary spacings we considered. The hybrid method results compare well with those of the direct numerical computation. Again, the direct numerical solution was significantly more time-consuming since it involved redefining the geometry and the grid for each e and for each set of capillary locations. There are two ways in which we can test the no-interaction limit of the hybrid method solu-tion. One way is to include only the first term on the right-hand side in (4.58), that neglects the capillary interaction term in (4.61). The other way is to neglect the off-diagonal terms, representing interaction, in (4.63), which results in (4.64). For the J = 2 case, we included the corresponding pmin(£) curves of pG, the first term on the right-hand side of (4.58), and of (Pa)diag from (4.64). We see that the pmin(£) curve corresponding to pG lies above the other curves, indicating that the global effect of capillary interaction is to lower the oxygenation in the tissue. In their multiple capillary oxygen transport model, Clark et al. [5] found a reduction 58 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue Figure 4.6: (a) Locations of 4 circular capillaries at ( ± j / 4 cos(7r/4), ± j / 4 sin(7r/4)), for j = 1,2,3, within a circular tissue domain of radius 1. (b) Minimum oxygen partial pressure pm%n versus capillary radius £ for the tissue geometry shown in (a). Parameter values: pc = 5, M = 0.3, K = oo, V = 1. Top graph: j = 1,2,3 cases from hybrid solution (curves) and from numerical solution (discrete points). Bottom graph: for j = 2 case, hybrid solution including interaction; and hybrid solution excluding interaction by using only the first term in (4.58) and by using the diagonal term approximation in (4.64). in oxygen partial pressure levels in the tissue due to capillary interaction. On the other hand, the Pmin(s) curve corresponding to (pG)diag lies below all of the other curves, which is what one would expect for the effect of capillary interaction. Figure 4.6(b) is similar in spirit to that of Ward et al. [57] for the interaction of four circular perfectly-absorbing obstacles in an absorption time distribution problem. For the particular case of four capillaries of radius e = 0.04, Figure 4.7 displays the minimum oxygen partial pressure pmin for specific values of r, the radial spacing of the capillaries, that we computed using the hybrid method. Figure 4.7 shows that maximum oxygenation occurs when the radial spacing of the capillaries is approximately 0.6 in the tissue domain of radius 1. Example: Capillary cross-sectional shape. We demonstrate the effect of capillary cross-sectional shape for one capillary that is concentric with the rectangular tissue domain, - 2 < xx < 2, -1.5 <x2< 1.5. We assume that K — oo so that the capillary is infinitely permeable. Then, from Rans-ford [45], we have d; for five cross-sectional shapes: (i) circle, (ii) ellipse, (iii) equilateral triangle, (iv) isosceles triangle, and (v) square. With e as the radius of the circular cross-section, we 59 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue 4.95 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) (b) Figure 4.7: (a) Radial spacing r of the 4 capillaries, (b) Minimum oxygen partial pressure pmin versus r of 4 circular capillaries of radius e = 0.04 in a circular tissue domain of radius 1, using the hybrid method. Parameter values: pc = 5, M = 0.3, K = oo, V — 1. choose the side length scales for the remaining four shapes so that each of the five bound-aries enclose the same area. As in the previous example, we consider homogeneous tissue with diffusivity V = 1 and oxygen consumption rate M = 0.5, and the capillary pressure pc = 5. In Figure 4.8, we have plotted the minimum oxygen partial pressure pmin versus £ and see that the pmin(£) curve for the circular cross-sectional shape lies below all the other curves. This illustrates the isoperimetric inequality (see Garabedian [14]) that for various cross-sectional shapes enclosing the same area, the minimum d; occurs for a circular cross-section. Actual capillary cross-sectional shapes range from circular to somewhat triangular with rounded edges. The figure suggests that a slight variation from the circular cross-section will increase the oxygenation of the tissue. Example: Variable capillary wall permeability. To illustrate the effect of variable cap-illary wall permeability, we consider one circular capillary concentric with an elliptic tissue domain with semi-axes a = 2 and 6 = 1 . For a circular capillary cross-section, the shape-dependent parameter is di = exp(—V(£i)/ni), where £ t is the location and K 8- is the capillary wall permeability of the i th capillary. The parameters are: constant diffusivity V — 1, oxygen consumption rate A4 = 0.3, and. capillary pressure pc = 5. Figure 4.9 shows the minimum oxygen partial pressure pmin as a function of n for the capillary cross-section of radius £ = 0.1. There is a dramatic variation in minimum pressure for K up to approximately 10, above which, 60 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.1 Figure 4.8: Minimum oxygen partial pressure pmin versus e, a measure of the capillary cross-sectional size, for various capillary shapes that all enclose an area of we2, of one capillary concentric with a rectangular tissue domain with —2 < x\ < 2 and —1.5 < x2 < 1.5. Parameter values: p c = 5, M. = 0.5, K — oo, V = 1. Figure 4.9: Minimum oxygen partial pressure pm%n versus K for an elliptic tissue domain with semi-axes a = 2 and 6 = 1 , containing one concentric, circular capillary of radius e = 0.1. Parameter values: pc = 5, M = 0.3, V = 1. the capillary wall is essentially infinitely permeable. We note that as K - > 0, this solution approaches the physically unreasonable case of a fully insulated capillary wall. 61 Chapter 4. Oxygen Transport from Capillaries to Skeletal Muscle Tissue X Figure 4.10: Effect of myoglobin (Mb) facilitation on the oxygen partial pressure p(x, 0) versus radial position x in a circular tissue domain of radius 1 containing a concentric circular capillary of radius £ = 0.05. Parameter values: A4 = 0.3, pc = 5, pF = 0.5, po.5 = 1, K = oo, V = 1. Example: Myoglobin facilitation. To demonstrate the effect of myoglobin, whose presence in the tissue can enhance or facilitate oxygen transport from the capillaries, we have plotted in Figure 4.10 the oxygen partial pressure p(x,0) along the positive a;-axis radius of the circular tissue domain of radius 1 containing one concentric circular capillary of radius e = 0.05 in the presence or absence of myoglobin in the tissue. When myoglobin is present, we use (4.8) and (4.9) with facilitation pressure pF — 0.5 and myoglobin half-saturation pressure p0,5 = 1. For both cases, we use the parameters: oxygen consumption rate M — 0.3, diffusivity V = 1, capillary pressure pc — 5, and capillary permeability coefficient K = oo. The plot indicates that there is more oxygen in the tissue with myoglobin present than without. This oxygen transport problem was the second of two applications of the hybrid method to prob-lems on bounded domains. In the next two chapters, we consider two problems on unbounded domains. 62 Chapter 5 Convective Cyl indr ical Heat Transfer Past Small Bodies We present asymptotic solutions to a convective heat transfer problem in an unbounded, two-dimensional domain. We consider an array of small, cylindrical bodies of arbitrary shape in a fluid with a known velocity field and constant ambient temperature. Outside of these bodies, we solve for the temperature field. This is the first of two applications of the hybrid method on an unbounded domain. The second application on an unbounded domain involves low Reynolds number fluid flow past an asymmetric cylindrical body, which we discuss in Chapter 6. We assume that the surrounding fluid is incompressible and that we can neglect any effects of natural convection and viscous dissipation. With these assumptions, the equation governing the temperature T is where we have treated the thermal conductivity k, the density p, and the heat capacity at constant pressure cv as constants. We assume that we already have the solution to the Navier-Stokes equations for the velocity field u, completely independent of the temperature. We illustrate that this is a singular perturbation problem involving an infinite logarithmic expansion in the small parameter £, the order of magnitude of the size of the cylindrical bodies. The convective heat transfer problem is analogous to the problem of low Reynolds number flow past an array of cylindrical bodies. Both problems have the infinite logarithmic expansion structure in their asymptotic solutions, but the heat transfer problem is easier to analyze than the low Reynolds number problem for multiple bodies. In Chapter 7, we discuss an extended application of the hybrid method on low Reynolds number fluid flow past an array of symmetric cylindrical bodies. In Section 5.1, for a single cylindrical body, we illustrate the singular nature of the convective heat transfer problem and outline the hybrid asymptotic-numerical method that we use to con-struct an asymptotic solution for the steady-state, dimensionless temperature. In Section 5.2, for an array of cylindrical bodies, we generalize the convective heat transfer problem and con-struct an asymptotic solution to the temperature for an arbitrary velocity field. In Section 5.3, we derive a further term in the asymptotic solution expansion that reveals the asymmetry of the temperature field near each cylindrical body. In Section 5.4, we derive asymptotic solutions (5.1) 63 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies for two specific cases of the velocity field: a uniform flow field and a simple shear flow field. Finally, in Section 5.5, we present examples to illustrate the asymptotic results. 5.1 Singular Nature of the Problem We seek the steady-state solution to (5.1) that we non-dimensionalize with respect to the magnitude of the free-stream velocity, the ambient temperature and a characteristic length scale of the cylindrical bodies. The dimensionless and steady-state form of (5.1) is Aw = Pe(v • Vw), (5.2) where iv and v are the dimensionless temperature and velocity field respectively, and Pe = UoolpCp/k is the Peclet number. Here, is the magnitude of the free-stream velocity and / is the length scale of the bodies. We also note that the differentiation is now with respect to dimensionless variables. We will demonstrate the singular nature of the problem by attempting a regular perturbation expansion in the asymptotic limit Pe —> 0. For simplicity, we consider only one body of circular cross-section and radius 1 with a constant temperature a on its boundary. The problem to solve for w, with appropriate boundary conditions, is Aw = e(v-Vw), |y| > 1, (5.3a) w = a, |y| = 1, (5.3b) w ~ 1, |y| —»• oo, (5.3c) where e = Pe <C 1 and |y| is 0(1) in the Stokes region near the body. Leal [32] and Romero [46] examined the three-dimensional analogue of this problem in their studies of flow past a sphere. We assume a regular perturbation expansion for w(y; e) of the form w(y; e) = W^(y) + eW^(y) + •••. (5.4) To leading order in e, the problem to solve is AWW = 0, iy(°) ~ 1, y| > 1, (5.5a) y| = 1, (5.5b) y| —>• oo. (5.5c) The general solution to (5.5a) is = a + 6 log |y | , where we apply the boundary conditions to specify the constants a and b. However, in setting a = a to satisfy (5.5b), we find that it is not possible to satisfy (5.5c). Hence, a regular perturbation expansion fails. This difficulty is analogous to the Stokes paradox in fluid dynamics: no solution exists in an unbounded, two-dimensional domain that satisfies the condition at infinity. To circumvent this difficulty, we adopt a singular perturbation approach and construct local (inner) and global (outer) ex-pansions with appropriate length scales in each subdomain to resolve the non-uniformity at infinity. 64 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies We declare y to be the local variable and w(y;e) the local solution. In the global region, far from the body, we rescale the variables as x = ey and u(x;e) = w ( £ _ 1 x ; e ) and expand the local and global solutions as w(y; e) = a + fi{e)W^\y) + / 2 ( e ) ^ 2 ) ( y ) + • • • , (5.6) u(x; e) = 1 + P 1(e)^( 1)(x) + F2(s)U^2\x) + •••. (5.7) The logarithmic growth in W^(y) as |y| —>• oo requires that we choose the gauge functions, fj and Fj, as fJ{s) = F3{s)= [~^-£) , i = l , 2 , . . . . . (5.8) This reciprocal logarithmic form of the gauge functions, as in (2.14) with d = 1, enables us to match the local and global expansions without difficulty. In particular, since ( —1/loge}7' > £ for any j = 1,2,. . . , we can obtain from (5.3a) and (5.3b) that the local solution for the temperature is of the form °° / 1 V lu(y;e) = a + w^(y)J2^[-l—-) +••• ' = 1 V 8 J (5.9) Here, w^c\y) is a canonical local solution with a logarithmic far-field structure, w ( c ) ( y ) ~ l o g | y | + --- , | y | - » oo. (5.10) Also in (5.9), the cij for j > 1 are constant and A(e) is asymptotic to an infinite logarithmic series of the form A<e> ~ f> h ^ y (5-n) This asymptotic expression is the same as (2.19) with d = 1. Thus, we have shown that the problem in (5.3) is singular in nature and that its solution contains an infinite logarithmic expansion in terms of the small parameter e. Our goal is to avoid calculating the individual coefficients aj in (5.11) by formulating a problem that determines A(e) directly. We outline here the hybrid asymptotic-numerical method for this single body case, which we will extend in detail to the case of many bodies in the next section. Rescaling to the global variables in (5.3a) and (5.3c), and using (5.8), we can write the global solution as °° / 1 V ,(x;e) = ! + $ > ( - — j G(x) G(x)A(£) + i = i 1 log£ (5.12) 65 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Here, G(x) is the Green's function that satisfies AG = v-VG, x e f t 2\{0}, (5.13a) G -+ 0, |x | -> oo, (5.13b) G = log |x | + i2 + o(l), x ^ O , (5.13c) where R is the regular part of the Green's function. We define i/(e) = — l / l o g £ and, from the second term of the global expansion (5.12), we set u-l = u* where u* ~ u(e)G(x)A(e), £ 0. (5.14) We will formulate a related problem for u*, one that is straightforward to solve, and show that its solution contains the infinite logarithmic expansion A(e). Using (5.13c) in (5.14), we obtain, as x —s- 0, that the leading-order form for u* is u* ~ v(e)A(e) log |x| + v(e)A{e)R. (5.15) Using (5.10), the far-field form of the local solution, in terms of the global variable, is w(y;e) ~ a + i/(e)A(e)[log|x|+ i / - 1 (£ ) ] +••• , |y |-»• oo. (5.16) We require that the global solution as x —> 0 match to this far-field form of the local solution. Hence, we get that u* satisfies A u * = v W , x £ U 2\{0}, (5.17a) u* -»• 0, |x| -> oo, (5.17b) u* - y(e)A(£) log |x | + a + A{e) - 1 + o(l), x -+ 0. (5.17c) The first term in the right-hand sides of (5.15) and (5.17c) automatically agree, and the re-maining terms gives us the expression for A(s), w s 1 — ct . The expression for A(e) in this unbounded-domain problem is similar in form to the expression in (2.31) for a bounded domain. In formulating this related problem, we have avoided calculating the individual coefficients cij in (5.9) and (5.12) and, instead, have determined A(e) directly. 5.2 Matched Asymptotics for an Arbitrary Flow Field We now formulate the general convective heat transfer problem for an array of N small, cylin-drical bodies Df of arbitrary cross-sectional shape. The surrounding fluid has a known velocity field v = v(x) where x is the Oseen variable. The governing equation for the temperature u(x;e), the boundary condition at each body surface and the condition at infinity are A« - v • Vu = 0, x e TZ2\ U Df, (5.19a) £—- + Ki(u- cti) = 0, x e dDf, i = 1,... ,N, (5.19b) on u ~ 1, |x | = (x\ + 3 , - 2 )2 0 0 . (5.19c) 66 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Here, Df is the cross-section of the ith cylindrical body of "radius" 0(e) for i = 1,.. . , TV and Df is the scaled ith body, such that Df = eDf. Also, d/dn is the inward normal derivative to the body and dDf is the boundary of Df. In the boundary condition (5.19b), which represents imperfect Newtonian cooling on the body surfaces, a; > 0 and K ; > 0 are constants, where a; and K ; represent the temperature on the boundary and the surface heat transfer coefficient respectively. The condition at infinity in (5.19c) represents the ambient state of the temperature far from the bodies. In contrast to the approach of Leal [32], we study the convective heat transfer problem in the Oseen region where x is 0(1) far from the bodies of size 0(e). Although the small parameter e now appears as the measure of the body cross-section size and not in the governing equation as in the one body case of Section 5.1, the singular nature of the problem remains the same. We use the method of matched asymptotic expansions to construct a hybrid asymptotic-numerical method to solve (5.19) in the limit e 0. In the local region, close to the ith body located at £ i 5 we rescale the variables as y; = £ _ 1 ( x - £t-) and u>i(yi;e) = u(eyi + £,•;£)• These local variables_have the same form as (4.27). From (5.19a) and (5.19b), the temperature, Wi(yf,e), about the i th body satisfies Aywt - ev(eyi + £ ) -VyWi = 0, yt- £ D°, (5.20a) + Kt(Wi - ai) = 0, Yi € dD°i, (5.20b) where the subscript y indicates differentiation with respect to the local variable. We expand the local solution as w,-(yt-;e) = a< + W^0)(yl; vu ... ,vN) + eW^)(yl;v1,... ,vN) + ••• , (5.21) and take w/°^(y;; z^,... , vN) = AiViw\°\yi). Here, A,- = Ai(v\,... , vN) is arbitrary and analogous to the infinite logarithmic series that we defined in (5.11), V{ = Vi(edi) = — l/log(ed,-) and w\c\yi) is the canonical local solution that satisfies y, i A° , (5.22a) y,- G dDQi, (5.22b) |y,-| -> <x>. (5.22c) This, is a multi-body version of the canonical local problem in (2.10). In (5.22c), d; = di(ni) is a positive constant that depends on the shape of the ith body. When K,- = oo, d8- represents the logarithmic capacity (see Garabedian [14]) of Df. Section 2.5 in Chapter 2 contains details on the computation of d;. Also in (5.22c), p is a vector that depends on K{ and on the shape of the ith body. The far-field form of the local solution, to leading order in e and in terms of the global variable, is Wi(yt;e) ~ a; + A^-flog |x - &| + v'1} + • • • , |y,-| -+ oo. (5.23) Ayw\c) = 0, w i c ) ~ log |y«l dn - log di + p KiW] ' = 0, y«-67 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies In the global region, far from the bodies, we expand the solution as «(x; e) = 1 + c7(°)(x; vx, ...,vN) + £t / ( 1 ) (x; vN) + e2U^{x; vu ... , vN) + • • •. (5.24) We substitute (5.24) into (5.19a) and (5.19c) and require that match to the far-field, form of the local solution in (5.23) as x —> £{ for each i = 1,... , N. Hence, we get that f/t0) satisfies AUW - v • v*7(0) = 0, x e n2\ u {&} , (5.25a) i=l £/(°) -* 0, |x| -»• 0 0 , (5.25b) U {0) = Aji/.-log |x - &| + ai + A l - l + o(l) , x -» i = l,...,N. (5.25c) To solve for in (5.25), it is convenient to introduce the Green's function G'(x; £ ) that satisfies A G ' - v - V G = 0, x G K 2 \ { ( } , (5.26a) G - > 0 , |x | —* oo, (5.26b) G = l o g | x - C | + P + o(l), x ^ C - (5.26c) Here, R = R(C) is a constant that we determine from the solution to (5.26). If the velocity field v is independent of x, then by translational invariance R is independent of the location, £, of the singularity. Now, using the principle of superposition, the solution for is N (°>(x; « * , . . . , „ „ ) = £ AkvkG{x; (5.27) fc=i As x —» the leading-order form of C/(0' is AT [/(°) - A ^ t log |x - & | + A ^ + AkukG(£i, £fe). (5.28) Comparing the first term in the right-hand sides of (5.28) and (5.25c), we observe that they automatically agree. Comparison of the remaining terms in these expressions gives us the linear system of N equations for A ; , N [ViRi - l]Ai + AkukG(ti;Zk) = a, - 1. (5.29) k=l k^i In (5.28) and (5.29), Ri = R(£i) is the regular part of the Green's function that we can obtain from the solution to (5.26). We can also obtain the constants d,-, for i — 1,. . . , N, from the solution to the canonical local problem in (5.22). In some cases we can determine dt- analytically (see Section 2.5). Since z/,- = — l/Tog(£d;) —>• 0 as e —> 0, the system (5.29) is asymptotically 68 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies diagonally dominant and so is non-singular in this limit. To leading order, as e —*• 0, we have from (5.29) that Ai~l- a,. (5.30) This leading-order approximation, which is just the difference between the temperature at infinity and at the bodies, ignores the off-diagonal terms in (5.29) that represent the interaction between the bodies. In Section 5.5, we examine the asymptotic results for the special cases of N = 1 and N = 2 cylindrical bodies. The hybrid method that we have just described combines asymptotic analysis and simple nu-merics to solve (5.19). Using the hybrid method, a change in the cross-sectional shape of the i th cylindrical body requires only the computation of the shape-dependent parameter cZ; from (5.22). This process is independent of e, the order of magnitude of the size of the bodies. However, in a direct numerical solution to (5.19), restructuring the solution grid is necessary for a change in cross-sectional shape as well as for a change in e. Thus, the hybrid method has the advantage of avoiding such a computationally intensive procedure. In summary, after specifying the parameters: £ order of magnitude of the size of the bodies N number of bodies ai temperature on the ith body £j location of the ith body di body shape-dependent parameters Ri regular part of the Green's function, one can solve (5.29) for Ai and hence determine the global solution C/(°). Recall that Ai is analogous to the infinite logarithmic series in the small parameter e that we stated in (5.11). Hence, the hybrid method essentially sums the entire logarithmic expansion of the solution with an error \u — (1 + U^)\ = 0(e) away from the bodies. We calculate a further term in the asymptotic expansion in the next section. 5.3 A Higher-Order Term in the Expansion In Section 5.1, we showed that the solution for convective heat transfer from a small, cylindrical body involves an infinite sum of powers of l / l o g £ , where e represents the order of magnitude of the size of the bodies. Then in Section 5.2, we showed how the solution to a hybrid problem that is asymptotic to the desired solution essentially sums the entire infinite logarithmic series. In this section, for the special case of N cylinders of circular cross-section with radii p,-e for i = 1,.. . , /V, we show how to continue the asymptotic expansion of the solution to include the first term that is smaller than all positive powers of l / l o g £ . We will demonstrate how this term reveals the asymmetry of the solution near each cylindrical body of circular cross-section. In the local region, we substitute the expansion (5.21) into (5.20a) and (5.20b). As in Section 5.2, we take wf^ = AiUiw\°\ where w\ c^ is the canonical local solution to (5.22). For bodies of 69 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies circular cross-section, the vector p in (5.22c) vanishes and we can express analytically as W$0) = A ^ U o g \yt\ - log di). (5.31) Here, the coefficient Ai is the solution to (5.29) and the body shape-dependent parameter di (see Section 2.5) is di=-pie-1/Ki(,i. (5.32) Continuing to the next order in the expansion and using v ( £ + eyi) = v ( £ J + 0(e), we find that the solution W-- 1^ to the 0(e) local problem satisfies A „ w f > = v (&) • Vywl°\ yt £ A° , (5.33a) + nlWy ) = 0, yt- G dDl (5.33b) Substituting (5.31) into the right-hand side of (5.33a), we obtain that the solution w\ 1^ is =B,Vi[log |y,-| - logdt] + Q • (yt- - A y , / | y i | 2 ) ( \ \ (5.34) + AiViviti) • ( -y.-log |y,-| + E i y M * • To satisfy the boundary condition (5.33b), the constants Di and Ei are Di = pi ( , 2* = f ( ( l - ^ ) l ° 6 f t + l ) . (5.35) V KiPi + 1J 2 V 1 + KiPi ) In (5.34), we will determine the unknown constant vector C,- through a leading-order matching condition and the unknown Bi through a higher-order matching condition. Substituting (5.31) and (5.34) into (5.21), we obtain that the far-field form of the local solution, to leading order in e and in terms of the global variable, is u>i =cti + AiVi log |x - £,| + Ai + Ci • (x - f,•) + ^ [ l o g |x - £.| - loge] v (&) • (x - + • • • , ( 5 " 3 6 ) as |y;| —r oo. In the global region, far from the bodies, we expand the solution u as in (5.24), where the leading-order global solution is in (5.27) and the Green's function G satisfies (5.26). To match the global solution with the (x — ^ ) l o g | x — £j| term in (5.36), G(x;£j ) must have the form G - log |x - & | - R{ii) = ^ • (x - ii) log |x - t{\ + a(&) • (x - + • • • , (5.37) as x —> £ i , for some constant vector a ( £ t ) . As x —> £ i 5 we can decompose the solution to the leading-order global problem in (5.27) as the sum of a singular part representing the influence of only the i t h body and an analytic part representing the interaction with the other bodies, N [/(°) = AiViGfrti) + Y,AkVkG(x;tk). (5.38) fc=i 70 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Substituting (5.38) into (5.24) and using (5.37) and expanding the analytic part in a Taylor series, we obtain that the form of the global solution as x —> £ i 5 for i = 1,... , N, is u=l + AiVi log |x - £-| + A^R , Axv, v(&) • (x - &) log |x - &| + AiVi&tfi) • (x -iV JV + £ A ^ G ( f c ; + £ A ^ V G f o ; £fc) • (x - £,.) + • (5.39) k=l fe=l The leading-order matching condition requires that (5.36) and (5.39) agree. First, we observe that the (x — £ 8 ) log | x - £ J and log |x —£ J terms in these expressions automatically agree. Next, matching the 0(1) terms, we obtain the linear system for A; in (5.29) as we expect. Finally, in matching the (x — terms, we determine the constant vector C,- as (1 s N ^ p v ^ ) + a ( & ) J + £ AkvkVG(ti; £k) (5.40) k = l kjti where G is the solution to (5.26a) and (5.26b) with the singular form in (5.37). In order to determine fully the constant vector C,-, we must obtain the constant vector a(£,-). To do so, we introduce polar coordinates x = + (i\ cos 0, i\ sin 0) and we compute the Fourier coefficients of the left-hand side of (5.37) at r,- = 8 <C 1, G\c)(6) = - I V[G(rl,0)-logrl-Rl] K Jo G\S)(S) = - r[G(rl,9)-logrl-Rl n Jo r{=S cos 0 d0, Ti=S sin 0 d0. (5.41) (5.42) In general, the Fourier coefficients, G\ c\6) and G\s\s), depend on Using the right-hand side of (5.37) in (5.41) and (5.42), we obtain that G\C\S),G\S\S) Thus, the constant vector a (£ t ) is •Slog 8+ ^ 6 , 0. a ( £ i ) = J i m (5.43) (5.44) Next, we describe how to determine the constant B{ in (5.34). Substituting (5.31) and (5.34) into (5.21) and taking the far-field expansion up to 0(e), we find that the solution to the 0(e) global problem satisfies AU(1) - v VC/W = 0, x £ ft2\ U { £ } , OO, 1 7( 1 ) ^ 0 , | x _ UM ~ J B ^ [ l o g | x - e 4 | + ^ - 1 ] , x - & , i = l,...,N. (5.45a) (5.45b) (5.45c) 71 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Using the principle of superposition, the solution for is JV UM(x; vu...,vN) = Y, BkVkG(x; 4 ) . (5.46) k=i In analogy to the linear system for A,- in (5.29), we obtain the N equations for U ; , [uiRi - l]Bi + BkVkGfa tk) = 0, i = 1,.. . , N. (5.47) Since the system in (5.47) is asymptotically diagonally dominant as e —• 0, it is non-singular in this limit. This determines that B{ = 0 for i = 1,... , JV, and thus = 0. In summary, the global solution for an array of A circular cylinders is JV u= l + ^ A f c ^ G ( x ; ^ ) + 0 ( £ 2 ) , (5.48) k=i and the local solution near the i th cyUndrical body is u>i =oti + A,;/i[log |y;| - log dt] + £ | c , - • ( y i - A - y , 7 | y , f ) + AiUiw(d) • Q y » i o g M + ^ . -y .V ly i l 2 ) } (5-49) + 0(s2). Here, A{ is the solution to (5.29), the constant vector C,- is in (5.40) and the constants D{ and E{ are in (5.35). We can see that the 0(e) term in (5.49) involves y,- and not just the magnitude |y; | . This indicates that the temperature field is asymmetric near each cylindrical body, even for a body of circular cross-section. We will illustrate this effect in examples 5.5.2 and 5.5.3 in Section 5.5. 5.4 Two Specific Flow Fields 5.4.1 A r r a y of C y l i n d r i c a l B o d i e s i n U n i f o r m F l o w We examine the special case where the velocity field is a uniform flow in the positive a^-direction and where we fix the temperature on the boundary of each cylindrical body. Thus, in the general problem (5.19), we set the velocity field to v = (1, 0) and m = oo for i = 1,.. . , A . In this case, we determine analytically that the Green's function G(x; £), satisfying (5.26) with v = (1,0), is Xl - Cl A o CI <2(x; C) = - e x P where £ = (C15C2) and R'o(z) is a modified Bessel function. As x G(x; C) ~ log |x - C| - log 4 + 7 + • • (5.50) C, G(x; C) has the form (5.51) 72 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies where 7 = 0.5772... is Euler's constant. Comparing (5.51) and (5.26c), we obtain that R is independent of £ and is # = 7 - l o g 4 . (5.52) Substituting (5.50) into (5.27) and using (5.24), we have that the temperature in the global region satisfies N u(x; e) = 1 - £ A{i/i exp t=i Here, 671 is the a; 1-coordinate of and V{ = — 1/log(edj-). We obtain the coefficients A,-, for i = 1,.. . , N, by substituting (5.50) and (5.52) into (5.29) which yields the linear system N [vi(j - log 4) - l]Ai - 53 Ak"k exp x • 0(e). (5.53) 6,1 - 6c,1 Kn \ti-tk a,- — 1. (5.54) The off-diagonal terms in (5.54), containing exp[(£ t )i - &c,i)/2], reflect the asymmetry of the temperature field in the global region. To leading order, as e —• 0, we have Ai ~ 1 - a;. Thus, to leading order, the A» are independent of the locations of the small, cylindrical bodies and hence the asymmetry of the temperature field. In Section 5.5, we test the validity of this no-interaction limit for the specific example of two identical cylinders of elliptic cross-section. 5.4.2 Array of Cylindrical Bodies in Shear Flow Next, we consider the special case of N small, cylindrical bodies of arbitrary cross-section in a simple shear flow velocity field. Frankel and Acrivos [13] examined the case of one small cylinder of circular cross-section in their study of heat and mass transfer from spheres and cylinders in shear flow. As in Section 5.4.1, we fix the temperatures on the boundaries. In this case, u satisfies (5.19) with v = (X2, 0) and = 0 0 for i — 1,... , N. To calculate R, the regular part of the Green's function, we introduce the new variables X = x — £ and ip(X.) — G(x; £) into (5.26) with v = (^2,0), giving AiP - (X2 + (2, = 0, x e ft2\{o}, dXi ^ - > 0 , | X | -> 0 0 , ^ = l o g | X | + J? + o(l), X ^ O , (5.55a) (5.55b) (5.55c) where X = (Xi,X2). We can see from (5.55) that R depends only on ( 2 , the .T2-coordinate of £. The solution to (5.55) is (see Bretherton [3]), <Kx)=- r Jo (1 + ^ / 1 2 ) " It exp (Xi - \X2t - Q2tf (X2f At(l + t2/12) At dt. (5.56) 73 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Thus, in terms of x , the solution to (5.26) with v = (x2,0) is ( l + i 2 / 1 2 ) - 2 G(x;C) = -o 2t exp [(xi-Cl)-l(x2-C2)t-(2tf • ( 0 . - 2 - C 2 ) 2 4i(l + t2/12) At dt. (5.57) Figure 5.1: The regular part of the Green's function R versus £^2 for the simple shear flow case. To determine the form of G(x; C) as x —> C> from (5.26c) we write R= l i m { G ( x ; C ) - l o g | x - C | } . As in [3], we add the modified Bessel function A'o(|x — CI) a s R(C) = lim (G(x; C) + A 0 ( | x - C|) - log 2 + 7 ) . This makes the integral in (5.57) converge when x = C- Thus, (C 2 ) 2 i E ( ° = - r M ( i + i 2 / i 2 ) " e x p ( - 4 ( T + i 2/12) — e dt - log 2 + 7. (5.58) (5.59) (5.60) Since the velocity field v = (x2,0) for the simple shear flow depends on spatial position, the constant R{ = R(Xi) i s dependent on the location Here, R{ depends only on £ t i 2 , the x2-coordinate of Figure 5.1 contains the plot of R, that we computed numerically from (5.60), versus the x2-coordinate of the location of the ith cylindrical body. We obtain the linear system for Ai in the simple shear flow case by substituting the expressions for Rt = Rfc) from (5.60) and from (5.57) into (5.29). 74 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies A -0.3 -0.32 -0.34 b -0.36 -0.4 -0.42 1 1 1 1 1 Al £ = : A2 • • • -i i i i = 0.02, a = 1.5, d = 0.75 • r T ° — • i i i i i i i i i 0 50 100 150 200 250 300 350 400 450 500 r Figure 5.2: Test of no-interaction limit: A \ and A 2 versus the separation distance 7" for two identical, elliptic cylinders aligned with the uniform flow. 5.5 Results and Discussion We present the results of our study through various examples. In one particular example, we compare the asymptotic solution with the exact analytical solution. E x a m p l e 5.5.1 We consider two identical, elliptic cylinders with constant temperatures on the boundaries in a uniform flow with v = (1,0). The inset in Figure 5.2 shows the alignment of the cylinders in the uniform flow. The cylinders have semi-axes ea and eb, where a — 1.0 and b = 0.5, so the body shape-dependent parameter for both cylinders is d = 0.75 (see Section 2.5). Since the constant temperature a on each cylinder boundary is greater than that of the surrounding fluid, the two cylinders behave as heat sources. In Figure 5.2, for e = 0.02, we have plotted the values of Ai and A 2 from (5.54) versus r, the separation distance between the two cylinders. The figure indicates that there is still a significant interaction between the two bodies even at r — 500. This suggests that the no-interaction limit is not valid at this value of £. In Figure 5.3, for r — 5, we have plotted A \ , A 2 and the first approximation CLQ versus e. The first approximation cto is the solution to (5.54) ignoring the off-diagonal terms, and is a — 1 ^(7 - log 4) — 1 where v = —l/\og(ed) and 7 = 0.5772 . . . is Euler's constant. The results reveal that ao is a reasonable approximation for A \ , the coefficient associated with the upstream cylindrical body. 75 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies However, ao greatly underestimates A 2 , the coefficient for the downstream body that lies in the "wake" of the other. -0.2 -0.25 -0.3 -0.35 -0.4 -0.45 1 1 1 r a = 1.5, d = 0.75 > £b3—-"l i 1 1 r Ax ••• A2 ••• a 0 — 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Figure 5.3: Test of first approximation: Ai, A 2 and the first approximation ao versus the small parameter e for two identical, elliptic cylinders aligned with the uniform flow. In Figure 5.4, for e = 0.02 and r = 5, we show plots of A\ and A 2 versus 9, the alignment angle between the cylinders (see the inset in Figure 5.4). We observe that the coefficients for the two cylinders are equal at 9 = TT/2. This indicates that the first approximation a 0 that we plotted in Figure 5.3 would be a reasonable approximation for both A\ and A 2 in this case. E x a m p l e 5.5.2 We consider N circular cylinders of radius e with constant temperatures on the boundaries in a uniform flow. In this case, p; = 1, Ki — oo and v = (1,0) for i — 1,. . . , JV. We use (5.50) and (5.52) in (5.41) and (5.42) to obtain the Fourier sine and cosine coefficients, G {s) = 0, G<c) ~ - ^ ( l o g 4 - 7 ) + i t f l o g t f , S^Q. (5.62) Thus, from (5.44), the components of a = (ai , a2) are ^ = - ^ 4 - 7 ) , «2 = 0. (5.63) Note that for a uniform flow, a is independent off,-. Using (5.63) in (5.40), we determine that the vectors C ; = (Ct,i, Ci,2), for i = 1,... , N, are A N d = -f- (log ! + 7 , 0 ) + £ AkvkVG{^ik). (5.64) ft=i 76 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies -0.35 -0.39 e = 0.02, a = 1.5, d = 0.75 -0.36 *2 -0.37 -0.38 -0.4 \--0.41 -0.42 A2 TT /4 TT / 2 Figure 5.4: A \ and A 2 versus the alignment angle 0 for two identical elliptic cylinders. The summation term in (5.64) represents the interaction between the cylindrical bodies. We compare the local solution W{ in (5.49) for two circular cylinders (JV — 2) to that for one circular cylinder (JV = 1) in a uniform flow. We substitute (5.64) into (5.49) in which Di — 1, Ei = 0 and Ai is the solution to (5.54). In Figure 5.5, for cylinders of radius e = 0.02, we have plotted contour lines of constant temperature 1.3 for each of two circular cylinders a distance r — 1.5 apart and the corresponding contour line for one circular cylinder alone. We plotted each contour in reference to the origin (^,1,^,2) = (0,0) of the i th cylindrical body. In comparison to the contour for one cylinder alone, the figure indicates a substantial expansion of the w2 contour, corresponding to the downstream body, whereas the w\ contour shows that the upstream body is relatively insensitive to the presence of the other body. Hence, we can see that the main interaction between the two cylinders in the uniform flow is the influence of the upstream body on the "wake" of the one downstream. E x a m p l e 5.5.3 We show how to calculate even further terms in the asymptotic expansion of the solution for one circular cylinder of radius e with a constant temperature on its boundary in a uniform flow. In this case, JV = 1, p = 1, v = (1,0) and K — 0 0 , where we have dropped the subscript 1 for brevity. Thus, the body shape-dependent parameter in (5.32) simplifies to d — 1 and from (5.31), we can write Ai/ log |y | . Here, A = (1 - a.)/[v (log 4 — 7 ) + 1] is the solution to (5.54) and v = we obtain that D = 1 and E = 0. Hence, (5.34) simplifies to 1\ Av , QJ 2 cos 0, (5.65) -1/loge. From (5.35), (5.66) 77 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Vi,2 0 a = 1.5, d = 1.0, £ = 0.02, r = 1.5 w -i i i i i -15 -10 5 10 15 20 25 Figure 5.5: Local solution temperature contours at w = 1.3 for two circular cylinders (wi and iu2) and for one circular cylinder (w) of radius £ in a uniform flow. where y = (£> cos 0, £> sin 0) and g = |y | . Also, C\ = (a - l ) / 2 from (5.64). We solve for V F ^ 2 \ that satisfies (5.33a) and (5.33b) with v = (1,0) and K = 0 0 , and obtain T..(o\ (C\ Av\ 2 Av 2 l 7 1 Az/ C i ^(2) = ( T - T6 j 5+ T e loge" Mog^ + 16 - T + (5.67) cos 20. To determine the constants (5 and 6, we substitute (5.65) and (5.67) into (5.21) and express the local solution in terms of the global variable. A leading-order matching condition requires that P = (a — 1)/16. Then, retaining only terms of 0(e2), we obtain the singular form for the global solution f/(2) of the expansion in (5.24). Thus, the solution to the 0(e2) global problem satisfies (5.45a) and (5.45b) and, as r = |x — f | —> 0, has the singular form TT(7\ . , „ cos0 „ cos 20 b Av IJW ~ -felogr - C i C i — + — . r A v 16 The form of as r —>• 0 requires that we express the solution as U {2\r,9) = 6 exp r cos 9 K C i exp r cos 0 A i ( - ) cos 1 (5.68) (5.69) Here, Ko{z) a n ( 1 Ki(z) a r e modified Bessel functions. As x — f , or as r —> 0, the expression in (5.69) has the form TTO\ c o s # ,-1 cos2 9 ; / 1 . /y(2) ~ -fclogr - C i Cl~^2~ + H l o g 4 ~ 7)- (5.70) 78 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Matching requires that (5.68) and (5.70) agree. We observe that the first three terms in these expressions automatically agree. Comparing the remaining terms, we find that the unknown constant b is Av2 v2{l-o) 16[1 + / / ( l o g 4 - 7 ) ] 16[1 + i/(log4 - 7 ) ] 2 ' Thus, the global solution is of the form u = 1 + + e2U^2h (5.71) We can find the exact analytical solution to (5.19) for this special case. Following Philip et al. [39], we let i(x) = 1 - (1 - a)eXl/2H(x), (5.72) and introduce polar coordinates (x\ - £\,x2 - £2) = (r cos 0, r sin 0). Here, £ = (£1,62) is the location of the circular cylinder. Then H satisfies AH - i f / 4 = 0, r > £ H = e-rcose'2, r = e He rcos0/ 2^ O, r 0 0 . The solution for II in terms of the polar coordinates is / / M ) = ^ ^^/ 2) + 2E(-1)nM^ /^2)cos,f9, n=l Kn(e/2) (5.73a) (5.73b) (5.73c) (5.74) where In(z) and Kn(z) are modified Bessel functions of order n. Changing back to the original variables, we obtain that the exact analytical solution uE to (5.19) with K = 0 0 about a circular cylinder of radius e in a uniform flow is uE(x; e) =1 — (1 — a) exp x\ - £1 J 0 (£ /2 ) / - | x - £ | ^ 0 2 V r - i r J " ( £ / 2 ) K l ) Kn(s/2)An n=l K0(e/2)~ v V 2 ^ ^ cos n# (5.75) where tanfl = (a,-2 - £2)/(a-'i - £i)-We compare the asymptotic solution with the exact analytical solution in terms of the heat flux Q across the boundary. The heat flux is JdD on J0 2ir du 0|x| edO |x|=£ 2*" 0w 0 5|y| d0. |y|=i (5.76) Using the asymptotic solution to compute the heat flux, we substitute the local expansion (5.21), using (5.65)-(5.67), into (5.76). To leading order in e, the asymptotic heat flux Q\ is Q i = - 2 - K A V . (5.77) 79 Chapter 5. Convective Heat Transfer Past Small Cylindrical Bodies Figure 5.6: Asymptotic and exact heat flux Q versus cylinder radius e for one circular cylinder in a uniform flow. To 0(e 2 ) , the asymptotic heat flux Q2 is Q2 = Qi-2we2(C1/2-b), (5.78) where the constant b is in (5.71). The 0(e) term in the local expansion vanishes in the integration and hence, does not contribute to the asymptotic heat flux. For computing the exact heat flux, we substitute the first three terms of (5.75) into (5.76) and use a Taylor expansion to 0(e 4 ) for the exp[(xi — fi)/2] factor. The exact heat flux QE is then of the form e - ^ - M f ^ M ^ ) ( 5 7 9 ) where z = e/2. Figure 5.6 displays the leading-order and 0(e 2 ) asymptotic heat flux, Q\ and Q2, versus the exact heat flux QE. The leading-order asymptotic heat flux agrees reasonably well with QE for values of e, the radius of the cylindrical body, up to approximately 0.14. The agreement improves significantly after including the first of the higher-order terms in computing the asymp-totic heat flux. In the next chapter, we examine a second application of the hybrid method on an unbounded domain, but this time on a non-linear problem. The application of Chapter 6 models low Reynolds number fluid flow past a cylindrical body that is asymmetric to the uniform free-stream flow. In Chapter 7, we touch upon how to apply the hybrid method to a multi-body low Reynolds number fluid flow problem. 80 Chapter 6 Low Reynolds Number F lu id Flow Past an Asymmetr ic Cylinder In this chapter, we apply the hybrid method to a non-linear problem on an unbounded domain. We consider two-dimensional, steady, incompressible, viscous fluid flow at low Reynolds number about an arbitrarily shaped cylindrical body with a uniform free-stream velocity of magnitude f/co in the positive x\-direction. Low Reynolds number fluid flow can model the locomotion of micro-organisms (see Light hill [34]) with a Reynolds number in the range of 10~ 3 to 1. In measuring the force that the fluid exerts on the body, the dimensionless lift and drag coefficients are of particular interest. The equations for the velocity u = (ui, u2) and pressure p of the fluid flow are the Navier-Stokes equations, which are ( u - V ) u = - - V p + i/Au, (6.1a) V - u = 0. (6.1b) Here, u and p are functions of the spatial variable, x = (xi,X2). Also, p is the density, and v is the kinematic viscosity of the fluid. A l l of the quantities above are dimensional. The boundary conditions of the problem are u = 0, x<E<9D; u - ^ ( C / o o , 0 ) , |x | -* oo. (6.2) Here, dD is the boundary of the body. Typically, we express the problem in terms of dimen-sionless variables using a characteristic length scale L of the body and the magnitude of the uniform free-stream velocity, Ucx>- With this form of non-dimensionalization, we are able to identify the dimensionless quantity known as the Reynolds number, Re = XJ^L/u. The problem of slow, steady uniform flow past a stationary body is rich in history. In the mid nineteenth century, Stokes considered an approximation of the Navier-Stokes equations in which he neglected the effects of inertial forces. The "Stokes paradox" refers to his inability to find a solution to the resulting Stokes equations in two dimensions. About sixty years later, Oseen determined that the Stokes equations were not valid at infinity and constructed a linear, first approximation to the Navier-Stokes equations which he could solve in a region where the flow is nearly uniform. We can see that this low Reynolds number problem is a singular perturbation 81 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder problem, in which the Stokes approximation is valid in a local (inner) region, close to the body, and the Oseen approximation is valid in a global (outer) region, far from the body. In 1957, Kaplun [24] and Proudman & Pearson [44] used the method of matched asymptotic expansions to resolve the Stokes paradox in two-dimensional, steady, viscous flow past a circular cylinder. Proudman & Pearson formulated the problem in terms of the dimensionless stream function, whereas Kaplun used velocity and pressure in his formulation. In their separate studies, each were able to determine analytically the hrst two terms in an asymptotic expansion for the drag coefficient and, with some difficulty, Kaplun was able to determine the third term. Also, Kaplun remarked on how to obtain the form of an expansion for a cylinder of arbitrary cross-sectional geometry. In Chapter 2, we discussed the equivalence principle of Kaplun that asymptotically links the drag coefficients of cylinders of any cross-sectional shape. Twenty-five years later, Shintani, Umemura & Takano [49] applied the method of matched asymptotic expansions to determine the lift and drag coefficients of an elliptic cylinder in low Reynolds number fluid flow. They were able to obtain terms up to order (log R e ) - 2 in the inner expansion for the lift and drag forces acting on the cylinder. However, the truncated series for the drag and lift coefficients are only accurate for moderately small Reynolds number. Thus, further terms are necessary to provide reasonable accuracy for a wider range of low Reynolds numbers. For the special case of an elliptic .cross-section, we will compare the leading-order form of our lift coefficient result to theirs. Shortly thereafter, in 1986, Lee & Leal [33] numerically implemented the method of matched asymptotic expansions using velocity and pressure as variables in their study of low Reynolds number flow past cylinders of arbitrary cross-sectional shape. Like Shintani et al, they were able to determine expressions for the lift and drag force on the cylinders that were correct up to order ( l o g R e ) - 2 . We extend the analysis of Kropinski, Ward & Keller [29] who applied the hybrid method in calculating the drag coefficient, correct to all logarithmic terms, of a cylindrical body that is symmetric about the direction of the free-stream. We will refer to this as the symmetric case. In extending the work of Kropinski et ai, we allow the cylinder cross-section to be asymmetric with respect to the free-stream and hence, the body could have a non-zero lift force. For one cylindrical body, of arbitrary cross-sectional shape Do and asymmetric with respect to the free-stream, we construct an asymptotic solution for the lift and drag coefficients in the limit of Re —r 0. Applying the hybrid method enables us to sum all the logarithmic terms appearing in the expansions for the lift and drag forces, resulting in an error that is 0 (Re p ) instead of 0((log R e ) - 9 ) for some p and q. We introduce the stream function ip, in terms of the dimensionless fluid velocity components v = (wi,-y2), as Hence, the continuity equation (in dimensional form in (6.1b)) is automatically satisfied. As is dip dy v2 = -dip dx' 82 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder well-known, in terms of polar coordinates centred in the body, the stream function satisfies r £ D0, (6.3a) r e 8D0, (6.3b) r = (x2 + x2)1!2 -> oo. (6.3c) Here, e = Re = U^L/v <C 1 is the Reynolds number based on the length scale L of the cylinder cross-section Do, and Jr is the Jacobian, Jr(f,g) = r~1(frgg — grfg). In Sections 6.1 and 6.2, we outline the standard singular perturbation, analysis of (6.3) in the two regions of the solution domain; the "Stokes" (inner) and "Oseen" (outer) regions. Using the asymptotic structure that unfolds in the standard analysis, we apply the hybrid method to (6.3) in Section 6.3 and formulate a related problem for the stream function that we will solve numerically by extending the finite-difference code of Kropinski et al. [29]. The solution of the related problem contains the entire infinite logarithmic expansion of the flow field and the force coefficients. In Section 6.4, we describe certain details of numerically solving the hybrid related problem, including the necessary modifications to the symmetric case finite-difference code. We derive an asymptotic expression for the lift coefficient, CL, in Section 6.5, that is correct to all logarithmic terms. For the special case of an inclined elliptic cylinder, in Section 6.6, we determine an analytic expression for a body shape-dependent matrix M that is analogous to the parameter d of the symmetric case. In Section 6.7, we provide details of certain analytical formulae that we require in the numerical solution of the hybrid related problem. Finally, in Section 6.8, we illustrate the hybrid method results in terms of the lift coefficient, CL. A2ip + eJr(ip,Aip) = 0, on ip ~ r sin 0, 6.1 The Stokes Region In the Stokes (inner) region where r — 0(1), the stream function satisfies (6.3a) with e = 0 and (6.3b). We declare r and ip to be the Stokes variables and expand the stream function in the form CO iP(r,0;e) = Y,^(z)^(r,0) + ... , (6.4) where v(e) = - 1 / l o g e . Substituting (6.4) into (6.3a)-(6.3b), we find that ipj for j — 1,2,.. . satisfies A 2i>3 = 0, r$ D0, (6.5a) ^ = <tk = 0, re 8D0. (6.5b) on The asymptotic form of ipj as r —> oo involves linear combinations of { r 3 , r log r, r, r - 1 } cos# and {r3, rlogr, r, r - 1 } sin 0. To match with the Oseen expansion, the r3 terms must vanish. Thus, we can write the far-field form of ipj as ipj ~ a.j • \{r log r cos 0, r log r sin 9) + M ( r cos 0, r sin #)] + • • - , (6-6) 83 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder as r —r oo. We were inspired in our expression for the far-field form of the stream function ipj by the section on viscous flow problems of Hsiao & MacCamy [19]. In (6.6), &j = (acj,asj), for j > 1, are constant vectors that are independent of the Reynolds number, e. Also, M is a 2x2 matrix that depends on the cross-sectional shape of the body and on its orientation with respect to the free-stream. The matrix M is analogous to the body shape-dependent parameter d of the symmetric case. Thus, the Stokes expansion has the far-field behaviour CO tp ~ f'(e) SLJ • [(?- log ?' cos 9, r log ?' sin 8) 4- M ( r cos 8, r sin 9)], (6-7) as r —r oo. We can also write the far-field behaviour of the Stokes expansion as CO ^ ~ E ^ ( £ ) a r [ y l o g | y | + My], (6.8) as 7- = |y] —> oo. Here, y = (r cos 0, r sin 9). We compare this to the symmetric case that Kropinski et al. [29] studied in which, due to the symmetry of the flow field, it was only necessary to include sin 9 terms in the stream function expansion. 6.2 The Oseen Region In the Oseen (outer) region of the solution domain, where r = 0 ( £ _ 1 ) , we introduce new variables p = er, x = £y = (p cos9, p sin 9) and $ (p ,# ;£ ) = etp^pe -1, 9; s), and expand $ as CO Again, v = — l / l o g £ . Substituting (6.9) into (6.3a) and (6.3c) and matching $ as p —> 0 to the far-field form of the Stokes expansion in (6.7), we find that a x = (0,1) so that $ i satisfies (sin0 d d \ - C O S 0 — J A $ ! = 0, p>0, (6.10a) *x -> 0, p ^ o o , (6.10b) $ i ~ (logp + 7n 2 2 + a ^ p s i n f l + ( m 2 i + ac2)p cos8, p -> 0. (6.10c) We will elaborate later on the significance of the form of the first constant vector, a i . In (6.10a), Los is the linearized Oseen operator and $ i is the linearized Oseen solution. In (6.10c), m,-j is the entry in row i and column j of the matrix M . For j — 2, 3 , . . . , the functions $ j satisfy i - i Loa^j = -J2jp[^k,^3-k], P > 0 , (6.11a) fc=i * j -»• 0, p->oo, (6.11b) tfj ~ aj • [xlog|x | + M x ] + a j + 1 • x, p = |x | 0. (6.11c) 84 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder Again, the constant vectors a3 for j = 2 , 3 , . . . are independent of the Reynolds number, £. The solution to (6.11) recursively determines these constant vectors. For the symmetric body case, all the acj components of a.j vanish. As well for this case, we will see in Section 6.6 that the off-diagonal entries of the matrix M are zero. For £ —*• 0, the drag CD and lift coefficient CL for a cylinder of arbitrary cross-section are of the form [see Section 6.5] We begin the infinite sum in the lift coefficient expression at j = 2 since a c l = 0 (recall that a-i = (0,1)). At this stage, we see that with u = —1/loge, where £ is the Reynolds number, the coefficient of drag is 0 ( (£ log £ ) - 1 ) and the coefficient of lift is O ( £ _ 1 ( l o g £ ) - 2 ) . Kropinski et al. [29] noted that Kaplun's three-term expression for the drag coefficient provides a poor approximation of the experimental values unless £ is very small (see Van Dyke [56]). It is possible to compute numerically further coefficients in the series from the infinite sequence of partial differential equations. However, one would still have to truncate the series at some finite j. In Chapter 1, we demonstrated the poor accuracy of a five-term reciprocal logarithmic expansion at moderate values of the perturbation parameter. Instead of truncating the series for the lift and drag coefficients, we show that the hybrid asymptotic-numerical method allows us to sum all the terms in the infinite logarithmic series while avoiding the direct and tedious calculations of the individual coefficients in the asymptotic expansions. 6.3 The Hybrid Formulation We define the vector function A(e) as asymptotic to the infinite logarithmic series CO A(£) = ( A c ( £ ) , A s ( £ ) ) ~ ^ a ^ ' 7 " 1 ( £ ) ' £ ^ 0 - ^ 6 - 1 3 ) As in the previous section, the &j = (acj,asj) for j > 1 are constant vectors, and v = —1/loge where e is the Reynolds number. To obtain these constant vectors, it would be necessary to solve a recursive set of linearized, forced Oseen problems. In the symmetric case, Kaplun [24] was able to determine asj for j — 1,2,3. However, it is analytically intractable to calculate any more of the asj. For this reason, we will apply the hybrid method in order to find A(e) directly. We define the vector function ipc(p, 6) to be the canonical inner solution satisfying r i D0, (6.14a) r e 8D0, (6.14b) r = |y | oo. (6.14c) Again, M is a 2x2 matrix that depends on the shape of the body and is analogous to the shape-dependent parameter d in the symmetric case. The solution of the canonical inner problem - c - 0, dn ipc ~ y l o g | y | + M y 85 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder for a specific cylinder cross-section provides the matrix M for that shape. For certain cross-sectional shapes, such as an ellipse, it is possible to determine M analytically. In Section 6.6, we show how to determine M for an ellipse inclined at an angle a to the free-stream. The canonical problem in (6.14) is a vector analogue in terms of polar coordinates of the canonical local problem in (2.10). Using (6.13) and (6.14) together with (6.6), the Stokes expansion (6.4) is asymptotic to 1>(r, 8;e) = v(e)A(e) • ^ c ( r , $) + • • • . (6.15) Substituting (6.14c) into (6.15) and writing the result in terms of the Oseen variable x = ey, we get the far-field form ip ~ e _ 1 A ( e ) • [x + i / (£)xlog|x | + ; / (£)Mx] , |y | -» oo. (6.16) We now formulate the related problem for the stream function. The related problem for A(e) and the auxiliary stream function $ H = H(p,9; s) is A 2 $ f f + J P ( $ H , A $ H ) = 0, />= |x |>0 , (6.17a) # H ~psin6>, p —• oo, (6.17b) $ H ~ A(£) • [x + iv(e)xlog|x| + i / (e)Mx], p 0. (6.17c) The solution to the related problem will allow us to compute A(e). The related problem is a hybrid asymptotic-numerical formulation of the original problem (6.3) but in terms of the Oseen (outer) variables. In the related problem, we have replaced the boundary conditions on the cylindrical body in (6.3b) by the singularity structure (6.17c). We derived the form of the singularity through the far-field behaviour of the logarithmic expansion in the Stokes region. Applying the hybrid method reduces the problem to computing the solution to the parameter-dependent problem (6.17) instead of computing the solutions to the infinite sequence of outer problems in (6.11). In terms of A , the asymptotic formula for the lift and drag coefficients, valid to within all logarithmic terms, is (CL,-CD) = -^[v(e)A{e) + •• • ] , !/(£) = - 1 / l o g e . (6.18) We found that a i = (0,1), which we substitute into (6.13) to see that A c = 0(v) and that A s = 0(1). Thus, we see that the lift coefficient is smaller than the drag coefficient by a factor O f O ( ( l o g £ ) - 1 ) . 6.4 Numerical Solution of the Hybrid Related Problem in the Oseen Region In numerically solving the parameter-dependent hybrid related problem in (6.17), we first de-compose the solution as 0;e) = psm.e + v{e)A{e) • ( * o c , * o s ) + **(/»,0)- (6-19) 86 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder Using this expression, we will construct a problem to solve for that is regular as p —>• 0. In (6.19), $ o c and $ o s correspond to the cos n9 and sin?i# parts of the linearized Oseen solution respectively. From the form of (6.10), we have that $ o c and $ o s satisfy X o s ( * o c , $ o s ) = 0, p > 0 , (*oc,* O S)^0, p-±oo, $ o c ~./>log/9cos0, \P 0 S ~ plogpsinfl p -s- 0. Here, Los is the linearized Oseen operator as defined in (6.10a). The solution for tyos, from Proudman & Pearson [44], is (6.20a) (6.20b) (6.20c) 71 = 1 sin n9. (6.21) In this expression, Kn(z) and In(z) are modified Bessel functions of order n of the first and second kind. In a similar manner (see Section 6.7 for details), we determine $o c(/>, 0) to be Voc(p, 0) = -4j2 Ko($) U^) cos n9 (6.22) n=l For small z = p/2, the asymptotic behaviours as z —> 0 of the modified Bessel functions are Ii'o(z) ~ - log 7 + Jti(z) ~ - + I0(z) -! + •••, h(z) ~ - + (6.23a) (6.23b) In (6.23a), 7 = 0.5772 . . . is Euler's constant. Using (6.23) in (6.21) and (6.22), we obtain that the behaviours of \P 0 C and $ o s , as p —*• 0, are $oc ~ p log /O cos 6 + (7 — log 4)/9 cos $os ~ plogp sin0 + (7 — log 4 — l)/9 sin #, (6.24a) (6.24b) Substituting (6.19) into (6.17), and using (6.24), we have that is regular as p satisfies Los%* = -Jp[vA • (^oc^os) + uA • ( A * o c , A * o s ) + A $ * ] , $* 0, p 0 0 , $* - [A c + j / ( m n A c + m i 2 A s - A c ( 7 - log4))]p cosf? + [A s + u(m2\ Ac + m22As - A 5 (7 - log 4 - 1)) - l]p sin 0 + • • P>0, 0. *• 0 and (6.25a) (6.25b) (6.25c) Again, A = ( A c , A s ) and rn-ij is the i j th entry of the matrix M . In Section 6.7, we derive the linearized Oseen solution, $ o c , and determine analytical formulae for its various derivatives that we require to evaluate numerically the Jacobian Jp in (6.25a). The corresponding formulae for $ o s are in the paper of Kropinski et al. [29] for the symmetric case. 87 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder We can also write (6.25c) in the form $* = #( c)(p)cos0-r $ ( s )(p)sin6> + ••• , p ->• 0. (6.26) In (6.26), \p(c) and $( s) are the Fourier cosine and sine coefficients respectively, given by i r2-w $(<=) = _/ y*(p,0) cos 0d8, (6.27a) 1" Jo $(*) = ! / $*(p,fl)sin#d0. (6.27b) I" Jo Here, \1>'C' and depend on p <C 1 and on the vector parameter K = ( K C , K s ) = v A . We now outline how to determine the vector A(e), where £ is the Reynolds number. Comparison of the two expressions for as p —*• 0 in (6.25c) and (6.26) results in a 2 X 2 non-linear system for K of the form m 2 i v~x + m 2 2 - (7 - l°g 4 - 1) y-1 + ?Tlii - (7 - log 4) 77112 (6.28) For various K, we compute the solution from the parameter-dependent problem (6.25a) and (6.25b), noting that Ac — V~XKC and As = Z / - 1 K s . We fix p = S <C 1 and using (6.27), we compute a "table of values" for the right-hand side of the non-linear system involving the Fourier sine and cosine coefficients. In the specific case of low Reynolds number fluid flow past an elliptic cylinder, we can find the 7?7;J entries of the matrix M analytically. For a general cross-sectional shape, we can employ an integral equation method for biharmonic boundary value problems such as the one of Greengard et al. [15] to determine numerically the mt-7-. For a given Reynolds number, e, we have the value of v = — l / l o g £ . We can then solve the 2 x 2 non-linear system in (6.28) for n(e). We use a bi-cubic spline to interpolate values of the right-hand side of (6.28) at arbitrary n from our "table of values" of the Fourier sine and cosine coefficients. To solve (6.28), we employ Newton's method in which we compute the Jacobian with centred finite differences. Finally, A(e) = v~1n(e). With A(e) in hand, we calculate the coefficients of lift and drag from the asymptotic expression in (6.18). To compute the solution $* from (6.25), we extend the finite difference code of Kropinski et al. [29]. They based their code for the symmetric case on a stream function/vorticity formula-tion of the problem and solved the resulting non-linear system of equations for the unknowns $*(/?, 6) and w*(p, 0) = Aty*(p, 0) using Newton's iterations. They stretched the radial variable according to r = log(l-)-p) and applied a second-order centred discretization on a uniform polar grid to the equations in terms of the variables (r, 9). Exploiting the symmetry of the flow field, they solved for the unknowns $* and w* on the domain 0 < r < r c o , 0 < # < 7 r , where is an artificial far-field boundary. The main modifications to the symmetric version of the code were: to expand the solution domain to 0 < r < r ^ O < 0 < 2vr since the flow field is no longer symmetric in general; to determine the solution for various input vector parameters K = (KC,KS) instead of the scalar KC 1 + l i m ^ o lim p-»0 88 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder parameter « in the symmetric case; and using a Fourier cosine and sine expansion to produce a 2 x 2 non-linear system to solve for the vector A(s). In the symmetric case, Kropinski et al. [29] could exploit the symmetry of the flow field to restrict the solution domain to the upper half plane. The corresponding computational domain was then 0 < T < r o o , 0 < # < 7 r . For flow past an arbitrarily shaped cylinder, the flow field is, in general, asymmetric. Thus, we extend the code to solve for $* and w* on 0 < r < TCO?0 < 9 < 2ir and impose the periodicity conditions, $*(r, 0) = $*(r, 27r) and w*(r, 0) = OJ*(T, 2TT). We enforce periodicity in the computations through the addition of an artificial grid line in the direction. We add cos 9 terms to the linearized Oseen solution in the Jacobian term in (6.25a), which were not present in the symmetric case version of that equation, and compute the solution for various input vector parameters K — (KC, KS) using n values of KC and m values of KS. This is in contrast to the symmetric case where the parameter dependence was on a scalar K, corresponding to our KS. For a given KC<{, for i = 1,... , n, we compute $* for a range of KSJ, for j = 1,. . . , m, and at each stage, use the solution of the previous stage where K = ( K C ) ; , K S J _ I ) as the initial guess. When we step to the next value of KC, we use the solution from the stage where K = (KCJ, K S , I ) as the initial guess for the stage where K = ( K C I ; + I , K S J I ) . Kropinski et al. [29] matched the computed solution to its asymptotic behaviour near the origin in terms of its derivative, d^*/dr, which gave them a narrow range for the constants /1(e). In contrast, we follow the technique in Keller & Ward [25] of comparing the behaviour of the solution at the origin to an expansion in terms of its Fourier cosine and sine coefficients. This technique provides us with a 2 x 2 non-linear system to solve for K — vA, and in turn, to find A . In the next three sections, we present the details of determining the asymptotic expression for the lift and drag coefficients in terms of the vector A(e), the matrix M for an inclined ellipse, and analytic formulae for various derivatives of the linearized Oseen solution for use in our computations. In the last section of this chapter, we present graphical results of our study. 6.5 Calculating the Lift Coefficient CL In this section, we derive an asymptotic expression for the lift coefficient CL in terms of the vector A . At the same time, we link this expression to that of the drag coefficient from Kropinski et al. [29]. Following Imai [22], we have an expression for the drag force X and lift force Y, which is Here, p is the density of the surrounding fluid, is the magnitude of the uniform free-stream velocity, L is a characteristic length scale of the cylinder cross-section, p is the dynamic viscosity, C is a closed contour surrounding the cylinder, ip is the non-dimensional stream function, u is the dimensionless vorticity and z = x + iy, where (x,y) are the non-dimensional spatial coordinates. X-iY = - dz - ipU^L j (6.29) 89 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder To obtain the drag coefficient CD and the lift coefficient CL, we divide (6.29) by p U ^ L . Then, defining the Reynolds number to be £ = Re (6.30) we get that CD — iC'L = —2i / (-7r-\ dz - i i> LOZdtp — - <Jb ^zdz. Jc\dzJ Jc £ Jc dz (6.31) We convert this expression to polar coordinates using z = re and obtain, after a bit of algebra, that C, 2TT sin 9 and CL 2TT cos 9 The vorticity u is 8ip\2 sin9 (djj\2 + 2cos9dipdi} d ?- d9 r dr d9 2 r2ir d9 2TT f2v dtb rl f'iv dio r + r LO cos 0 - £ d0 / cos9—d9+-l u> cos 9 d8. (6.32) Jo 99 £ J0 dr £ J 0 (hp Ik cos9_fdipY 2sm9dipdip ~r^~ \~d9~ r2* _ r r I UJ sin 9-— d9 -\ /o 99 £ J o r dr d9 ,2 /•2TT d9 du> 2TT smv—d0 / ush\9d9. (6.33) dr e 1 to = -Aip. (6.34) Substituting the asymptotic behaviour of the Stokes (inner) expansion from (6.15) into (6.34), we see that the vorticity in the inner region is of the form u ~ -v(e)A(e) • Aipc(r, 9) + • (6.35) Here, tpc is the vector canonical inner solution satisfying (6.14). Using (6.14c) in (6.35), we get the far-field behaviour of the vorticity, which is 1 u ~ -2v(e)—(A(e) -y ) , as | y | o o . (6.36) Substituting (6.14a) and (6.15) into (6.32) and (6.33), we find that there are no contributions from the terms involving derivatives of the stream function tp. Considering only the expression for the lift coefficient, we can show that it reduces to 2 r2ir du r [2* cos9—-d9 + - / ucos9d9. 9r £ Jo (6.37) 90 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder Substitution of (6.36) into (6.37), gives us that the lift coefficient is CL = -^-v(e)Ac(e), (6.38) where Ac is the first component of the vector A(e) = (Ac(e), As(e)). Similarly, we obtain that the expression for the drag coefficient is 47T CD = —v(e)As(e). (6.39) Finally, we can write (CL,-CD) = -^-v(e)A(e). (6.40) Thus, we have linked the drag and lift coefficient in this asymptotic expression involving the vector A . In the next section, we show how to determine analytically the body shape-dependent matrix M for a cylinder with an elliptic cross-section, which is analogous to the parameter d of the symmetric case. 6.6 Determining M for an Ellipse For certain cross-sectional shapes, such as an ellipse, we can determine analytically the matrix M that we first introduced in (6.6). For non-elliptic shapes, we could employ an integral method technique such as the one of Greengard et al. [15] to determine numerically M from (6.14). We consider a uniform free-stream in the positive x-direction about a cylinder with an elliptic cross-section, such that the major axis makes an angle a with the free-stream (see Figure 6.1(a)). The ellipse has major semi-axis a and minor semi-axis b, where a > b. We rotate the ellipse in u (a) (b) Figure 6.1: (a) Inclined ellipse in the (x,y) reference frame, with polar coordinates x = rcosO, y = r sin 9. (b) Ellipse in the (w, v) reference frame, with polar coordinates u = r cos <fi, v = r sin 6. the (a;, y) reference frame by a positive angle a, using u\ = / cosa - s m a \ x\ vj ys ina cos a J \y 1 91 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder In the new (u, v) reference frame, which we show in Figure 6.1(b), the equation of the ellipse is (u/a)2+ (v/b)2 = 1. We will solve (6.5) for the special case where Do is an ellipse. Then, we will compare the far-field behaviour of the solution as r —s- oo with (6.6) to determine the matrix M . To solve (6.5) about an ellipse, it is convenient to convert (u, v) to the elliptic coordinates (£, ?/) using the definitions u = c cosh £ cos ?7, v = c sinh £ sin n, where c = (a2 — ft2)1/2. The boundary of the ellipse is the level line £ = £ 0 , where The solution ipj to (6.5a) in terms of elliptic coordinates is where f(0 = (f ~ fo) c o s h £ + sinh £o cosh £o cosh £ - cosh2 £o sinh £, = (£ ~~- £o) s m n £ - sinh £o cosh £o sinh £ + sinh 2 £ 0 cosh £. (6.42) (6.43) (6.44) (6.45a) (6.45b) To determine the behaviour of ipj as r —> oo, we first convert ( £ , 7 / ) to the polar coordinates (r, (p) in the far field. We note that u — r cos <p and v = r sin As £ —> oo (and hence, as ?' —* oo), we have that cosh£ ~ e^/2, s i n h £ ~ e ^ / 2 , £ ~ log(2p/c), 77 ~ (p. (6.46) Using (6.46) and (6.45) in (6.44), we obtain an expression for the far-field form of ipj, as r —> 0 0 , which is ^ ~ c { A j (log r ~ l o g ( ^ T ^ a + b r cos < — Bj I log 7' — log a + 6 Now, we revert to the polar coordinates (r, 0) of the original reference frame, where x = r cos 6 and y = rs in0 , using (6.41). As well, we define -(Aj cos a — B3 sin a), After a bit of algebra, we obtain that, as r —> 0 0 V'j ~ a c j r log r cos 0 + asj r log r sin 6 + \acj isj — —(Aj sin a + B3 cos a) (6.48) (b - a) cos2 a — b , fa + b log + < a, a + b (a — b) sin a cos a a -f b + asj (a - b) sin a cos a a + b (a — b) cos2 a — a , fa + b " T T T l o § r cos 0 rsinfl . (6.49) 92 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder Comparing (6.49) with (6.6), where a.j = (acj,asj) and is the i j th entry of the matrix M , we determine, for an ellipse with major semi-axis a and minor semi-axis b at an angle of attack a, that (b - a) cos2 a - b , fa + b\ m^ = -—a-Tb l o gbrJ' (6-50a) (a — b) sin a cos a m12 = TO2i = '—— , (6.50b) a 4 - b (a — b) cos2 a- a (a + b\ , m22=^ >— l o g ^ _ J . (6.50c) From (6.50), we see that for any angle of incline a, the matrix M for an ellipse is symmetric. For a = 0 (no angle of incline), M is diagonal and the 77122 entry must be equal to — log(de 1 / / 2), where d is the shape-dependent parameter of the symmetric case, given by a + b f b — a \ f/ = ^ e x p ( ^ T i o J - (6-5l) For no angle of incline (a = 0), 1 fa + b\ b m M = - l o g ( ^ _ J - — . (6.52) Indeed, this expression corresponds to — log(de 1 / ' 2) from the symmetric case. Thus, we have an analytical expression for the entries of the matrix M for an inclined ellipse. In the next section, we derive the linearized Oseen solution $ o c and certain analytical formulae of its derivatives that we require for solving the hybrid related problem numerically for the stream function. 6.7 The Linearized Oseen Solution, \& o c , and its Derivatives In order to evaluate numerically the Jacobian, Jp, in the hybrid formulation in (6.25a), we require analytical formulae for various derivatives of *&oc(p,0), which satisfies (6.20). In partic-ular, we require (i) dp^oc, (ii) <9#$0c, (iii) <9pA\P0C and (iv) deA^oc. At the end of this section, we will reproduce the formulae for the corresponding derivatives of $ o s ( / 3 , (?) that appear in the paper by Kropinski et al. [29]. First, we will outline the derivation of the solution $ o c and then use this solution to find the necessary analytical formulae of certain of its derivatives. To obtain the solution $ o c , we first define the negative vorticity u>oc — A $ o c , and let cooc be of the form uooc = e^ pcoseCJOC. Thus, (2>oc satisfies AOJOC — \ U J O C = 0. Using the technique of separation of variables, we let cooc — AnRc(p) cos nd. Substituting this form into the equation for u>oc gives that Rc(p) satisfies dz1 dz where z = p/2. This is the modified Bessel equation of order rt, which has solutions Rc(z) = I±n(z),Kn(z). (6.54) 93 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder To allow us the proper behaviour in the far-field, we choose Kn(z) since Kn{z) —• 0 as z —> oo, whereas In(z) —> o as z —• oo. Therefore, in general, LOoc = 53 BnKn C O S (6.55) n=0 Using (6.20c) in the definition tooc = A $ o c , we get that w o c ~ (2/p) cos# as p —>• 0, and so p -+0 . (6.56) -cos fle-*'008' 2 / 1 - cos 9 1 p cos ( We look at the asymptotic behaviour for Kn(z) for small z, which is T(n) (z\~n Kn{z) n > 0, (6.57) as z —• 0. The behaviour for iifo(-z) for small 2 is in (6.23a). Using this behaviour and the fact that we will need to integrate u>oc twice in obtaining $ o c , we set B N — 0 for n > 1 in (6.55). Wi th this, we have that uoc — BQKQ(P/2) + B\K\(p/2) cos0, and so B 0 K 0 ( £ ) cos i (6.58) To find $ o c , we need to integrate the previous expression. But first, we will rewrite the expres-sion a bit. Using a Fourier cosine series, we write e 2 ' BoKo(^) +B1K1(P-) cos9\ = J2K(p)cosn9, where bn(p) is bn(p)= - [\kp™t Jo i = 0 B 0 K 0 ( £ ) + 5 1iif 1(Q cos( cos n9 d,9. Using the identities [1], and we can write bn(p) as 2coszi cos22 = cos(zi — z 2 ) + cos(zi + £2), /„(*) = I y e z c o s 6 cos n9d9, bn(p) = 2B0K0 (0 4 (0 + 5 i K i /n - l ( | )+ /n+ l (^ (6.59) (6.60) (6.61) (6.62) (6.63) Next, we write $ o c(p,#) in the form CO (6.64) n=l 94 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder Here, we have started the summation at n = 1 since a term that is 0-independent cannot satisfy Los^oc — 0. From this form, we have that n=l 1 n2 y'n(p) + -y'n(p) ~ ~^Vn{p) C O S 710. (6.65) Thus, using (6.63), yn(p) satisfies 1 n2 Vn(p) + -V'n(p) ~ -jpVn = /„_i(|) + / „ + i ( | n > 0. (6.66) We let yn(p) = (p/n)fn(p/2) and z — p/2, and so fn(z) satisfies 3 1 — T ? 2 2r? /"(*) + -f'n(z) + = —{2?ltfl(*)[/n-l(*) + /n+l(*)] + 2B0K0(z)} Z z z The solution to (6.67), or equivalently (6.66), determines the solution for $ o c . We follow Shintani et al. [49] and try fn(z) of the form (3n fn(z) = —K0(z)In(z). Substituting this form into the left-hand side of (6.67), we obtain 2-^[K0{z)In{z) - Ki(z)I'n(z)}. 3.67) (6.68) (6.69) To agree with the right-hand side of (6.67), we require that (3 = 2Bo and (3 = —2B\. We accomplish this with B \ = 1, (3 = - 2 and BQ = - 1 . Thus, fn(z) = ( - 2 n / z )Ko ( z ) I n ( z ) and so ^ ( 2 ) ~n fn(2 Finally, we obtain that CO OC(/>, 0) = - 4 K0 (0 /„ (0 cos n0. (6.70) (6.71) n=l We need to verify that the solution for $ o c satisfies the far-field condition in (6.20b) and the behaviour for small p in (6.20c). For large p, the asymptotic behaviour of Ko(p/2) and In(p/2) is M , / 2 ) . ^ ( 1 _ _ L + ...} (6.72a) (6.72b) 95 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder as p —oo. Thus, we have that 1 f .3 ~ - < 1 h 2) p\ 4p + p oo, (6.73) which verifies that the solution satisfies the far-field condition in (6.20b). Using the behaviour of the modified Bessel functions for small p, in particular, that Ii'o(p/2) ~ - log p+log 4--y + 0(p2) and Ii(p/2) ~ p/4 + 0(p3), we verify that the solution satisfies the behaviour in (6.20c). Using (6.58) and (6.71), we obtain analytical expressions for certain derivatives of $ o c , which are d^oc = ex>2 * i ( % cos 6K + ^ (|)^(l)-^(f)^(f dg-$,-ir. = ex/2K0(^) psmt 2dpAVc -2K0 ( 0 cos 0 + Kx ( | ) cos2 6 + Kt (| 2T, (P^ - A i - ] cos I p V2> 2 0 « A * O C = - W 2 sm i -pA'o(0+/>A-1(0 cos6> + 2 A ' 1 ( | (6.74a) (6.74b) (6.74c) (6.74d) Here, x = pcosO. The corresponding analytical formulae for various derivatives of $os, from Kropinski et al. [29], are Op * OS 2dpA9os 2d0AVos -ex/2K0{'-)sml de^os = -pe x/2 + 2, exl2 sin i Kx [ - ) cos i ex/2K\^) [2 cos 0 - / 3 sin 2 9]. (6.75a) (6.75b) (6.75c) (6.75d) Again, x — pcos9. We require these expressions in numerically evaluating the Jacobian, Jp, in (6.25a). In the next section, we present the results of this application of the hybrid method in terms of the lift coefficient, CL. 6.8 Results and Discussion We present the results of our study through various examples using the specific cross-sectional shape of an inclined ellipse. The shape of the cylinder cross-section enters into the hybrid method solution through the matrix M , which we can determine analytically for an ellipse from (6.50). We use the numerical procedure of Section 6.4 to solve the hybrid related problem in (6.25) for the stream function \P* on a 100 X 100 grid with an artificial boundary condition corresponding to a radial variable of value p^ = 60. We use the solution for \P* to compute the Fourier sine and cosine coefficients in (6.27), which we require for solving the 2 x 2 non-linear system in (6.28) for A(e) = K(s)v~ l(e). We then use A(e) in (6.18) to obtain the coefficients of lift and drag, correct to all logarithmic terms. 96 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder C (hybrid) C D (Kaplun)| C L (hybrid) _J u .2 0.4 0.6 0.8 1 1.2 e (Reynolds number) 1.4 1.6 Figure 6.2: Drag coefficient, CD from the hybrid results and from Kaplun's three-term ex-pansion, and the lift coefficient, C'L, from the hybrid results, versus Reynolds number, e, of a circular cylinder (an ellipse with equal major and minor semi-axes a — b — 1). In examining the leading-order form of the asymptotic expansions for the lift and drag coeffi-cients in (6.12), we find that the drag coefficient is 0(vje) and that the lift is 0(u2/e), where v = — l / l o g £ . For an elliptic cylinder with major semi-axis a and minor semi-axis b inclined at an angle a to the free-stream at Reynolds number e, the leading-order form of our expression for the lift coefficient, CL, is s~ 4TT a — b . C i = - 7 -T; r s m a cos a. (6.76) e( log£)^a + 6 To obtain this leading-order form from (6.12), we need the first non-zero ac-component of the constant vectors, ELJ. In Section 6,2, we found that the first constant vector is a i = (0,1). We continue the matching procedure between the Stokes and Oseen regions to determine that the component ac2 of the second constant vector, a 2 , is ac2 = —m2\, where m2i is an entry of the body shape-dependent matrix M . Substituting the matrix M for an inclined ellipse in (6.50) and the form of ac2 into (6.12), we obtain our leading-order form for the lift coefficient of an inclined ellipse. We confirm the leading-order form of our result by comparing it with the leading-order expres-sion for the lift coefficient of Shintani et al. [49] in their study of low Reynolds number flow past an elliptic cylinder. Their expression for the lift coefficient, correct to 0 ( £ _ 1 ( l o g £ ) ~ 2 ) , can be written as 47T R(logR- t+)(logR- t-) sin 2a, (6.77) 97 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder where t± = - 7 + 41og(2) - log[l + b/a] ± 1 a — b cos 2a a + b Here, R = 2e and 7 = 0.5772 . . . is Euler's constant. The leading-order expression of Shintani et al. [49] in (6.77) is consistent with our leading-order form in (6.76). In their paper on flow past cylinders of arbitrary cross-sectional shape, Lee & Leal [33] numerically calculated the lift force, correct to 0 ( ( l o g £ ) - 2 ) , which corresponds to a lift coefficient of 0 ( £ - 1 ( l o g £ ) - 2 ) . They showed that their numerically calculated values for the lift agree with the analytical results of Shintani et al. [49] for an elliptic cross-section. Chester [4] examined the motion of an inclined elliptic cylinder at low Reynolds number and obtained an analytic expression for the lift force, correct to 0 ( ( l o g £ ) ~ 2 ) , that (once adjusted for the fluid in motion past a stationary body) is also consistent with our leading-order form. These results of previous researchers substantiate the form of the first term in our asymptotic expression for the lift coefficient as an infinite expansion of reciprocal logarithms. With the numerical solution of the hybrid related problem that we show in the next two examples, we obtain the lift coefficient correct to all logarithmic terms. 2 0 0 10 C L (hybrid) C L (Shintani et al.) C L (leading-order) 0.2 0.3 0.4 0.5 0.6 e (Reynolds number) 0.7 0.8 Figure 6.3: Lift coefficient, CL, versus Reynolds number, e, of an elliptic cylinder with major semi-axis a = 1 and minor semi-axis b = 0.5 at an angle of inclination, a — TT /4 , comparing the hybrid results with the leading-order form in (6.76) and the results from Shintani et al. in (6.77). Example: Varying Reynolds Number. First, we use the curve for the coefficient of drag as a check of the extended asymmetric code, to ensure that we compute the same solution as 98 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder the original symmetric code of Kropinski et al. [29]. For a circular cross-section of radius 1 (an ellipse with equal major and minor semi-axes, a = b — 1), Figure 6.2 displays the force coefficients versus the Reynolds number, £, showing a constant zero lift coefficient and the same coefficient of drag curve as produced in Kropinski et al. [29]. For the drag coefficient curve in Figure 6.2, we compare with the three-term expansion of Kaplun [24], which we can write as (CD)K ~ 4VTJ>[1 - 0.8669P2], 1 v = log 3.7027- log £ (6.78) In Figure 6.3, we plot the lift coefficient versus Reynolds number, e, for an elliptic cross-section with major semi-axis a = 1 and minor semi-axis b = 0.5 and at an angle of inclination of a = 7r /4. The figure contains the CL(e) curves from the hybrid method result correct to all logarithmic terms in (6.18), from the leading-order expression in (6.76), and from the result of Shintani et al. in (6.77). The plot shows reasonable agreement with the result of Shintani et al. [49], which is valid up to 0(v2/e). Figures 6.2 and 6.3 indicate that the three-term, expansion for CD in (6.78) and the leading-order expression for CL in (6.76) are valid only for a narrow range of e. This is analogous to the poor accuracy of the five-term reciprocal logarithmic expansion at moderate values of £ that we demonstrated in Chapter 1. Drag Coefficient Figure 6.4: Lift coefficient, CL, and drag coefficient, CD, versus angle of incline, a, of an elliptic cylinder with major semi-axis a — 1 and minor semi-axis b = {1.0, 0.5,0.2} at Reynolds number £ = 0.1. E x a m p l e : V a r y i n g A n g l e o f Incl ine. In Figure 6.4, we fix the Reynolds number, £ = 0.1, and vary the angle of incline, a, of the elliptic cylinder with major semi-axis a = 1 and minor semi-axis b = 1,0.5,0.2. In all cases, we see that at a = 0 and a = 7r /2 , when the ellipse is symmetric to the free-stream, the coefficient of lift is 0. For the case with minor semi-axis 99 Chapter 6. Low Reynolds Number Fluid Flow Past an Asymmetric Cylinder 6 = 1 , the cross-sectional shape is a circle, and as we expect, we see that the drag coefficient is constant and non-zero, and that the lift coefficient is 0. The figure also displays that the drag coefficient increases as the angle of incline increases, and decreases as the minor semi-axis 6 decreases. The graphs in Figure 6.4 are qualitatively similar in nature to those of Lee & Leal [33], in which they plotted the contributions to the force on an elliptic cylinder at 0((1/loge)2), where e is the Reynolds number, as a function of the angle of inclination, a. This low Reynolds number fluid flow problem was the second of two applications of the hybrid method on problems with unbounded domains. In the next chapter, we discuss some possible directions for future applications of the hybrid method. 100 Chapter 7 Other Applications of the Hybr id Method So far, we have shown the details of four major applications of the hybrid method: fully developed laminar flow in a straight pipe with a core (Chapter 3), oxygen transport from multiple capillaries to skeletal muscle tissue (Chapter 4), convective heat transfer past an array of cylindrical bodies (Chapter 5) and low Reynolds number fluid flow past a cylinder asymmetric to the free-stream (Chapter 6). A l l of these problems are independent of time, and, with the exception of the application in Chapter 6, linear. It is possible to extend the exploration of this thesis to a more general framework that would incorporate singular perturbation problems that arise in other disciplines, and that include unsteady problems and more general non-linearity. In this chapter, we discuss other possible applications of the hybrid method. One is anon-linear problem that could be viewed as an extension to the linear convective heat transfer problem of Chapter 5 as well as a rudimentary model of the non-linear low Reynolds number fluid flow problem of Chapter 6. We also comment upon extending the low Reynolds number application to fluid flow past an array of bodies symmetric to the free-stream. Our formulation of the application in Chapter 6 involved the biharmonic operator. We elaborate on another problem involving this operator, but this time on an eigenvalue problem that models the vibration of thin plates with small cutouts or concentrated masses. To end this chapter, we touch upon some possible extensions to the general framework of applicable problems that we presented in Chapter 2. 7.1 A Non-linear Model Problem To begin our discussion of future directions for the hybrid method, we describe a problem that is an extension to the convective heat transfer problem of Chapter 5. A n analogous problem, that we studied in Chapter 6, is low Reynolds number fluid flow past cylindrical bodies, although in that case, the governing equation is non-linear. Hsiao [20] considered a non-linear exterior Dirichlet problem that was a crude model for a class of problems in two dimensions, including stationary viscous incompressible flow past a cylinder. He considered 101 Chapter 7. Other Applications of the Hybrid Method problems for u = u(xi,X2',e) of the form On Au-eu—— = 0, x e f t 2 \ A ) , (7.1a) OXi u = a, x G dD0, (7.1b) u —> -a, |x| -» oo. (7-lc) Here, Do is a domain of radius 0(1) that is independent of e, 8DQ is its boundary, and a and a are constants. As a first step in constructing an asymptotic solution for the non-linear problem, he examined the singular nature of the linear problem 0, x e n2\D0, a, x G 8D0, a, |x | —> oo, where b is a constant. In essence, Hsiao's linear problem is the same as the convective heat transfer problem for one cylinder in a uniform flow. Hsiao used the linear analysis as a guide to construct an asymptotic expansion for the exact solution to (7.1), carrying the solution up to order ( l ogg ) - 2 . To this order, he rigorously justified the use of the method of matched asymptotic expansions on this problem. We could also apply the hybrid method to a class of non-linear model problems that includes (7.1) as a special case. This would be an intermediate step between the linear convective heat transfer problem of Chapter 5 and the non-linear low Reynolds number problem of Chapter 6. In contrast to Hsiao, applying the hybrid method to (7.1) would extend the solution further than the ( l oge ) - 2 term, since the method avoids the difficulty of having to calculate analytically each individual term in the infinite logarithmic expansion. To obtain an asymptotic solution to (7.1), we would employ the hybrid method as we outlined in Chapter 2, with the necessary modifications for a non-linear problem. Another similar problem is of the form x e TZ2\D0, (7.2a) x e dD0, (7.2b) |x| —oo. (7.2c) This is an extension of the analysis of heat convection past a small body to treat non-linear source terms to solve for the temperature «(x; e). For example, if F(u) = u4 - U^, this problem models radiative heating such as that for black-body radiation. 7.2 Low Reynolds Number Flow Past an Array of Cylindrical Bodies We could also consider low Reynolds number fluid flow past an array of N cylindrical bodies. In this case, we assume that the cross-sectional shapes are arbitrary although symmetric with Au — eb du dxx u u —> Au - eVv • V« = F(u), u = a, 102 Chapter 7. Other Applications of the Hybrid Method respect to the free-stream. We construct asymptotic expressions for the drag coefficient of each body. The governing equations are the steady-state Navier-Stokes equations, which, in terms of the dimensionless Stokes variables, are A v - V p = £(v • V ) v V • v = 0 ^ N x £ U Di with the boundary conditions v = 0, x G dDi, v —> i, p —¥ 1, |x | —> oo. Here, v(x ; e) is the velocity vector, p(x; e) is the pressure, e = Re is the Reynolds number based on a characteristic length L of the cross-section, and dDi is the boundary of the cylinder i, for i = 1,.. . A . Here, we have non-dimensionalized with respect to the length L, the magnitude of the free-stream velocity UQO and the free-stream pressure p^. Applying the hybrid method to these problems would enable us to determine asymptotic ex-pressions for the drag coefficient, correct to all logarithmic terms, and to measure the error in the approximation. In the method, we formulate a related problem whose solution contains the infinite logarithmic expansion. For this case involving multiple bodies, we would be able to examine the effect of the interaction between the bodies. Two different limits in terms of the Stokes variables that we may consider are: a "lumped-body" limit where the bodies of 0(1) radii are separated by 0(1) distances, and an "unlumped-body" limit in which the radii of the bodies have magnitude of 0(1) and are separated by distances of 0(1/e) where e <C 1 is the Reynolds number. In the first limit, we can lump the bodies together via one parameter d for the entire array of cylinders and in such a way, reduce the problem to the previous case that Kropinski et al. [29] examined for one cylinder. Greengard, Kropinski & Mayo [15] developed a boundary integral solver to perform the sophisticated numerical computation to determine this d. The lumped-body problem is analogous to the one-body problem in that one can consider the Stokes region containing all of the bodies. This was the approach of Lee & Leal [33] who used boundary integral methods to solve for low Reynolds number flow past multiple cylinders, assuming that they were close enough together so that the Stokes flow region encompassed both cylinders. Even though they were considering cylinders of simple cross-sectional shape, they only computed their solution to order ( log£)~ 2 . In the second limit, we cannot lump the array of cylinders together. This makes the "unlumped-body" problem analogous to the multi-body convective heat transfer problem of Chapter 5 in which we considered the long range interaction of the bodies. Each of the N cylinders in the array would have its own body-shape dependent parameter di associated with it, for i = 1,. . . , A . To treat the multi-body problem, it would be necessary to formulate the related problem of the hybrid method in terms of the velocity and pressure, instead of the stream function formulation that was possible in the one-cylinder case. There are further possible extensions for the low Reynolds number fluid flow application. One extension is to compute the higher-order terms, beyond all logarithmic terms, in the asymptotic 103 Chapter 7. Other Applications of the Hybrid Method expansion for a body of arbitrary cross-sectional shape. These terms, which are transcenden-tally small compared to the logarithmic terms, are necessary to reveal any asymmetry of the near-body flow pattern, such as the emergence of standing eddies, due to inertial terms. Skin-ner [50] showed how to determine a few of the transcendentally small terms for the case of a circular cylinder. Also, for one circular cylinder, Keller & Ward [25] extended the hybrid method to determine the solution beyond logarithmic orders. Another extension is to apply the hybrid method to slender-body theory, with possible applications in biofluiddynamics (see Lighthill [34]). Batchelor [2] considered slender-body theory for rigid bodies of arbitrary cross-section in Stokes flow. Khayat & Cox [27] examined the effect of small inertial terms on the motion of slender bodies using asymptotic expansions based on the ratio of cross-sectional ra-dius to body length, that they assumed to be small. They obtained up to three terms in the logarithmic series of the asymptotic solutions for the drag, lift and torque on the bodies. The slender-body problem can also be related to an extension of the capillary oxygen transport problem of Chapter 4 in which we would take into account the slow variation in the axial direction of the capillary. 7.3 A Biharmonic Eigenvalue Problem Another application of the hybrid method involves strongly localized perturbations of bihar-monic eigenvalue problems in a bounded, two-dimensional domain D. The strongly localized perturbations could be of two types: domain perturbations from the removal of N small sub-domains Df from D and imposing certain boundary conditions on the resulting holes, and domain perturbations from concentrations of mass in N small regions Df of the domain. For both problems, we are interested in the effects of the perturbations on the eigenvalue A of the biharmonic operator. The biharmonic eigenvalue problem models the free vibration of a uniform, thin plate where the eigenvalue is proportional to the square of the frequency of vibration. The theory of plate vibration has applications in structural engineering, such as the design of flat panels in machines and buildings. As well, plate vibration can represent certain percussion instruments that are essentially flat, metal plates (see Fletcher & Rossing [12]). We can contrast the biharmonic eigenvalue problem to the Laplacian eigenvalue problem, which models the vibration of an ideal membrane. The basic difference between a membrane and a plate is that the restoring force for a membrane is due primarily to the external applied tension, whereas for a plate the main restoring force is due to the stiffness of the material. The equation of motion for the small displacement U from equilibrium of a thin material with stiffness is of the form Utt = aAU - f3A2U. (7.3) Here, a = T/o~ and j3 = h2E/[12p(l — u2)] are positive constants depending on the tension T, the mass per unit area cr, the thickness h of the material, the Young's modulus E, the density p, and the Poisson's ratio v. For a typical thin plate, a is small, so that to a first approximation we can neglect the AU term. Similarly for a typical membrane, (3 is small and so we would neglect the A2U term in this case. We sketch out the application of the hybrid method on the first of the two biharmonic eigenvalue problems that we mentioned: that is, the perturbation of the domain D by removing N small 104 Chapter 7. Other Applications of the Hybrid Method subdomains Df and setting certain boundary conditions on the resulting holes. In (7.3), setting a = 0 and assuming harmonic vibrations with frequency ui so that U(x., t) = u{x)elLjt, we obtain that the governing equation for the perturbed linear biharmonic eigenvalue problem is o TV c A u - Au = 0, x G D\ U Df, (7.4) i=l where A = co 2 j'(3. We impose boundary conditions that represent clamped edges of the plate: u = = 0, x e 3D, x G 0£>f, i = 1,... , JV. (7.5) on Swanson [51] established the first few terms in asymptotic formulae for the eigenvalues of the biharmonic operator in bounded, two-dimensional domains. He showed a key feature in the solution to the clamped vibrating plate problem, (7.4) together with (7.5): in the limit of small subdomains, the perturbed eigenvalues do not converge in general to the classical eigenvalues for the unperturbed domain. Eastep & Hemmig [8] considered the vibration of a thin, uniform plate with a cutout whose boundaries deviated slightly from a circular shape. By approximating the boundary as a Fourier cosine series, they solved (7.4) with JV = 1 for the displacement u as a regular perturbation expansion and constructed an approximation for the fundamental frequency of the vibrating plate. The analogous membrane problem is one where the displacement u of the membrane satisfies Au + Xu = 0. This membrane problem (or Laplacian eigenvalue problem), with appropriate boundary condi-tions, has been the subject of much research. Ozawa [37] rigorously derived the leading-order behaviour as e —^ 0 of the eigenvalue X(e) for the Laplacian with a small hole of order e re-moved from a bounded, two-dimensional domain. Lange & Weinitschke [30] used the method of matched asymptotic expansions to determine a few terms in the logarithmic series of A(e), applying their results to the vibration of a rectangular membrane with one or two circular holes. Ward, Henshaw & Keller [57] were able to sum the entire logarithmic series of the eigenvalue by applying the hybrid method to the Laplacian eigenvalue problem. The second type of strongly localized biharmonic eigenvalue problem, which we will not discuss in detail here, involves concentrations of mass in JV small regions Df, each centred at of the domain. The problem is of the form A2u-Xpu = 0, xeD\uDf, (7.6) where p = p(x) satisfies K x ) = < ^ y ^ ) xeDf,i=l,...,N, 1 x £ D\ U Df. i=l (7.7) The solution of the problem would determine A(e) as e 0 for some range of m > 0 and for arbitrary pi. 105 Chapter 7. Other Applications of the Hybrid Method Ingber, Pate & Salazar [23] examined both experimentally and numerically the vibration of clamped plates with concentrated masses and spring attachments. It is possible to cast their formulation, in which they treated the masses and attachments as point sources, into the concentrated mass problem in (7.6) and (7.7). Doherty & Dowell [6] took an experimental approach to determine the response of a rectangular plate under an applied force, with and without point masses, taking into consideration the interaction between the masses. There is also an analogous concentrated mass problem for membranes. Leal k, Sanchez-Hubert [31] used the method of matched asymptotic expansions to examine the effect of a concentrated mass on the eigenmodes of a vibrating membrane in the limit as e —> 0, where e measured the size of the region of concentration. We now return to the first biharmonic eigenvalue problem. We outline the hybrid method applied to the biharmonic eigenvalue problem on a bounded, two-dimensional domain D with the removal of one small subdomains D£ from D. Our goal is to formulate an asymptotic expression for the eigenvalue A(e) for strongly localized perturbations of the bounded domain of the biharmonic eigenvalue problem. In some cases, the biharmonic eigenvalue problem is similar in nature to the stream function formulation of low Reynolds number fluid flow past one cylindrical body of Chapter 6. In these cases, we are able to apply the hybrid method since an infinite logarithmic expansion arises in the expansion for the eigenvalue A as it does in the expansion for the drag coefficient, and the canonical inner problem for both involve the biharmonic operator. In other cases, the asymptotic expansion does not involve reciprocal logarithmic terms and hence, we cannot apply the hybrid method and it would be necessary to establish a different asymptotic approach. The biharmonic eigenvalue problem to solve for the eigenfunction u(x;e) and eigenvalue A(e) is A2u-Xu = 0, xe D\D£ CH2, (7.8a) u = — = 0, x e dD, x G dD£, (7.8b) On J u2dx = 1. (7.8c) D\D£ Here, dD and dD£ are the boundaries of D and D£ respectively, and D£ is a small removed subdomain located at a point £ in D. The aim of applying the hybrid method to this problem is to determine the change in the simple eigenvalue An as a result of the perturbation. We assume that the position £ of the small subdomain lies on a nodal line, that is, the unperturbed solution UQ vanishes at x = ( . Without the nodal line assumption, the perturbed problem would have to have a point constraint at the centre of the removed subdomain. Although we assume that fio(C) = 0, we also require that VUQ does not vanish at £. As the first step of the hybrid method, we apply the method of matched asymptotic expansions to (7.8). We expand the eigenfunction u and eigenvalue A in the global region, where |x — £| = 0(1), away from the removed subdomain, in the form A(e) ~ A 0 + i/(e)Ai + • • • , (7.9a) u(x; e) - UQ(X) + v(s)ui(x) -) . (7.9b) 106 Chapter 7. Other Applications of the Hybrid Method Here, v = —1/ log e. In the local region, near the removed subdomain, we use the local variables y = e _ 1 (x - Q and v(y; e) = e~1u(ey + £; e), and we expand the solution in the local region as v(y; e) = K e > i ( y ) + W)fHy) + •••• (7-10) The scaling between the local and global variables in this biharmonic problem has the same form as the that between the Stokes and Oseen variables in the low Reynolds number application of Chapter 6. The local problems for Vj(y) in the biharmonic eigenvalue problem, where j > 1, are analogous to those of the Stokes region problems in (6.5)-(6.6) of the low Reynolds number flow application. The problems to solve for Vj(y), where j > 1, are A2v3 = 0, y £ D 0 , (7.11a) VJ = = 0, y e dD0, . (7.11b) o n VJ ~ a i • [ylog |y | + M y ] , | y | - » °o. (7.11c) Here, Do = z~xDe and M is the 2 x 2 matrix that we introduced in (6.6). We note again that M depends on the shape of Do-Using (7.11c) in (7.10), we find that the local expansion v(y\ s) has the same far-field, behaviour as ip in (6.8), as |y | —> oo. In the global region, we substitute (7.9) into (7.8), excluding the boundary condition on dD£. Thus, the governing equations and boundary conditions for Uj(x) , j > 0, are j — i A2UJ - X0Uj = (1 - e3o) Y ^j-kUk, x e -D\{C}) (7.12a) k=o Oil' u3 = —± = 0, x e 3D. (7.12b) on • ' Here, eik is the Kronecker delta function as defined in (4.33). In addition to the governing equations and boundary conditions above, the global functions tt j(x) must satisfy appropriate matching conditions to the local expansion and the normalization conditions / Y]ukUj-k «x = ej0. (7-13) M f O \f0 J For the global solution behaviour as x —> we expand the unperturbed solution uo(x) in a Taylor series about x = £ to get tto(x) = V 7 i 0 ( x ) | x = c • ( x - C) + • • • - (7.14) Here, we used our assumption that wo(C) = 0- To match between the global and local regions, we use (7.14) in (7.9b) as x —> £, and compare with the local solution in (7.10) as |y | —> oo using (7.11c). Noting that v scales like e^u, matching requires that a i = Vwo( x ) | x =c a n c l that the behaviour of ttj(x) as x —> for j > 1, is UJ ~ a j • [ ( x - C ) l o g | x - C | + M ( x - C ) ] + a J + 1 - ( x - C ) , x - ^ C - (7-15) 107 Chapter 7. Other Applications of the Hybrid Method This behaviour has the same form as (6.11c) in the low Reynolds number flow application that also involved the biharmonic operator. Considering only the problem for ^ i (x ) , we can write (7.12a) and (7.15) with j = 1 as A 2 « i - A 0 «i = AiUrj + 4-Trai • V<S(x - ( ) + • • • , x € D. (7.16) Here, V£(x - Q is the gradient of the Dirac delta function. Using (7.12b) and Green's identity, we require that < Lui, u0 > = < Lu0, ui >, (7-17) where L is the operator L = A 2 - A 0 . Since Lu0 = 0, and with Lu\ = XIUQ + 47rai • V<5(x - £) , we obtain the solvability condition J «o(x) [A x « 0 (x) + 47r a i • V5(x - Q] dx = 0. (7.18) D We solve this condition for A i , using that a x = Vuo|x=C: to obtain Ai = 47r|V7i0|x=c|2- (7.19) Hence, we have determined A i , the leading-order correction for the simple eigenvalue in (7.9a). To obtain the further corrections, \j for j > 2, the method of matched asymptotic expansions would require us solve an infinite set of such solvability conditions. We now present how the hybrid method will avoid determining each Aj individually. We introduce A*(e) as a function that is asymptotic to the sum A*(e) ~ A 0 + Ke)Ai + [^( £)] 2A 2 + • • • , s - 0. (7.20) In applying the hybrid method, we formulate a related problem to the original whose solution essentially sums the entire logarithmic series for A*. Using (7.11c) in (7.10), we can write the far-field behaviour of the local solution for the eigen-function as v(y;e) ~ u(e)A(s) • [ylog|y | + M y ] + ••• , | y | o o . (7.21) Here, A(e) is a constant vector that is asymptotic to an infinite sum as in (6.13) for the low Reynolds number fluid flow application. As well, we can express the local solution for f(y;e) in terms of a canonical local vector function as in (6.15). We formulate a related problem to the original perturbation problem. We solve for the auxiliary function uH(x.;e) which satisfies x G D\{Q, (7.22a) x e dD, (7.22b) (7.22c) D uH~ A ( £ ) - K e ) ( x - C ) l o g | x - C | + K £ ) M ( x - C ) + ( x - 0 ] , * - C- (7.22d) A2uH A uH dn (uH)2dx 0, 0, 1, 108 Chapter 7. Other Applications of the Hybrid Method Here, the boundary conditions at the removed subdomain are replaced, using the far-field form of the local solution, by the singularity structure in (7.22d). This singularity structure imposes two conditions; one on the cos# part and the other on the sin 6 part of the solution. We determine the matrix M by solving the local canonical problem in (6.14). To solve the related problem, it is convenient to remove the logarithmic singularity through a change of variables. We let % ( x ; e) = 4jri/(e)A(e) • V G ( x ; Q + e). (7.23) In (7.23), V G ( x ; C) is the gradient of the Green's function that satisfies A2G - XG = 0, x € TZ2\{Q, (7.24a) G ~ - ^ - | x - C | 2 l o g | x - C | , x ^ C (7.24b) The function t/*(x; v) in (7.23) is regular as x —*• C and satisfies A V - X*u* = 0, xeD, (7.25a) u* = -\-KVK • VG, x e dD, (7.25b) Bv* (9 — = -47TZ/A- — VG, x £ dD, (7.25c) on on u* ~ A(e)v(e) • M ( x - C), x -»• C- (7.25d) As well as (7.25), w*(x;i/) must satisfy the normalization condition J[u*]2dx=l. (7.26) The shape-dependent matrix M in (7.25d) is found from the solution to (6.14) for a given shape of Do. In Section 6.6, we showed how to determine M analytically when D0 is an ellipse. The Green's function G7(x;C) is G(x; 0 = -73^175 [ ^ o ( A 1 / 4 | x - C|) + 2 A 0 ( A 1 / 4 | x - fl)] . (7.27) Here, Y0(z) is a Bessel function and KQ(Z) is a modified Bessel function, both of order zero. We now outline how to determine u*(x;z/) and A* using the hybrid related problem. To solve (7.25), we decompose the solution as u* = A i u i + A2u2, A = ( A i , A2). (7.28) For /c = 1, 2; Ufc is a regular function as x —• C that satisfies A 2 u f e - X*uk = 0, x G D , (7.29a) «t = - 4 T T I / | ^ - , x € 51?, (7.29b) 109 Chapter 7. Other Applications of the Hybrid Method We solve for uk numerically, for k = 1,2, as a function of A*. In terms of this solution, we compute Ck and Sk, for k = 1,2, which are Ck = U m — / uk cos Od0, Sk = U r n — / uksmOdd. (7.30) irp J0 P-^O irp J0 Here, p = |x — CI and we note that Ck and for k = 1,2, depend on A*. Then, as x —• £, we have = A i (C ipcos0 + S1pslne) + A 2 (C 2 /9cos0 + S2p sin 6). (7.31) Comparing (7.31) with (7.25d), we have that Ai(C\p cos 0 + sin 0) + A2(C2p cos 0 + S2p sin 0) ~ Aii/(mii />cos0 + ?72i2/0sin0) -f A2v(m2\p cosO + m22psin 0). Then, equating the coefficients of cos# and sin (9, we obtain C\ — 7/77211 c 2 - 2/7722i "Ai" "0" S\ - 7/77212 5*2 - 7/77222 A 2 0 (7.32) (7.33) We require that the determinant of the matrix on the left-hand side is equal to 0, which provides the equation for A*. Noting again that Ck and Sk, for k = 1,2, depend on A*, the equation is ( C i - vmu)(S2 - 7/77222) - ( £ 1 - vf}\X2)(C2 - 7/7722i) = 0. (7.34) This equation gives A* as a function of v = — l / l o g £ and in terms of the body shape through the m,-j coefficients. We calculate the ratio of A i to A 2 from this system. Without loss of generality, we set A i = 1, and then, the normalization condition in (7.26) provides the final relation between A i and A 2 . In the next, and last section, of this chapter, we continue our description of extended applica-tions for the hybrid method. 7.4 Extensions to the General Framework We present briefly some further possible applications of the hybrid method that extend the general framework that we discussed in Chapter 2. One could apply the hybrid method to study linear diffusion problems with localized non-linear reactions. Peirce & Rabitz [38] considered a problem of this type in their study of the effect of defect structures on two-dimensional chemically active surfaces. The governing equation for the concentration 7t(x, t; e) of a chemical species on a bounded catalytic surface D is of the form fill — - Au + e'pV{e-1\x-<:\)F(u) = 0, x e f l . (7.35) Here, V is a potential that measures the finite range of the reaction that is located at x = £ (ie. V ( | y | ) —> 0 as |y | —> 0 0 ) and p measures the strength of the reaction. We note that if p — 2 in two-dimensions, then the reaction has the behaviour of a Dirac delta function. 110 Chapter 7. Other Applications of the Hybrid Method Another application of the hybrid method is to examine low frequency scattering of light by small cylindrical bodies Df of arbitrary cross-sectional shape. The scattering of a plane wave is modelled by the Helmholtz equation A $ + A;2$ = 0, x£D\(J?=1Df, (7.36) with boundary conditions and far-field condition of the form $ = 0, x E dDf ,i = 1,... , A , (7.37a) $ ~ e t k o c , |x| oo. (7.37b) Here, $(x;e) is the scattering field and k is the wave number. In the low frequency limit, ks <C 1 where £ is a measure of the cross-sectional size of the bodies. The form of the governing equation conforms to the second-order framework for applicable problems that we described in Chapter 2. In a series of articles (one of which is [55]), Twersky studied this problem, which is the optical version of the low Peclet number heat transfer application in a uniform flow of Chapter 5. Thus, mathematically, the solution involves reciprocal logarithms and so we can apply the hybrid method. A n unsteady application of the hybrid method is to consider a perturbed-boundary diffusion equation such as the diffusion of a chemical species out of a mostly impervious container. The problem to solve for the concentration c(x,t;e) of the chemical species would be of the form x E D, (7.38a) N C x E U dDf, (7.38b) x E dD\ U dDf, (7.38c) j=i with an initial concentration c(x, 0;e) = Co(x;e). The boundary dD is mostly insulating, but contains many small regions dDf,... ,dD£N where the material could leak out. The goal of applying the hybrid method would be to calculate the amount of material in the container after a large amount of time and to determine how long it takes it to leak out completely. This problem has applications in the study of controlled-release drug therapy. To solve this strongly localized domain perturbation problem, we would consider an eigenvalue expansion for c(x, t; s) of the form CO c(x,t;e) = Yaie~Xjt(f )j(^£)- - (7-39) i=o To obtain asymptotic estimates for large time t, we need to study the eigenvalues Aj. In particular, the resulting eigenvalue problem for the first eigenpair (</>o, Ao) is a slight variation on the second-order eigenvalue problem in Ward et al [57]. Thus, we can formulate this leaky container problem so that we may apply the hybrid method. The last application of the hybrid method that we mention here is in the field of electrochem-istry. Lucas et al. [35] studied a reaction-diffusion problem with a periodic array of circular microelectrodes on the surface of a fluid. Their interest was in calculating the steady-state c = 0, dc dn 111 Chapter 7. Other Applications of the Hybrid Method current flow due to chemical processes at the electrodes, which is proportional to the con-centration flux of the chemical species. The governing equation for the concentration of the chemical species in their model is the modified Helmholtz equation, having the form of (2.1a) with c = 1 and JV(x, $ , V $ ) = a 2 $ , for some positive constant a. The boundary conditions are of mixed type, which fit into the framework of (2.1b) and (2.1c). Thus, the hybrid method could be applied to a variation of this electrochemistry problem that would allow for arbitrary shape and location of the many microelectrodes on the fluid surface. We have called attention to some possible applications of the hybrid method to mathematical models that occur in many different fields. The problems that we mentioned include models of fluid flow at low Reynolds number, vibration of thin plates with small holes, small defect structures on chemically active surfaces and low frequency scattering of light. 112 Chapter 8 Conclusions The goal of this thesis was to demonstrate a hybrid asymptotic-numerical method for treating two-dimensional singular perturbation problems whose asymptotic solution involves reciprocal logarithms of the small perturbation parameter, e. In Chapter 1, we illustrated the difficulty of slow convergence of an infinite expansion S(e) of the form The purpose of this hybrid method is to treat the slow convergence problems of asymptotic expansions of this form (since, for the problems that we considered in this thesis, we believe there is sufficient evidence of convergence of these expansions for small enough e) and to im-prove the accuracy of approximate solutions. The hybrid method uses the method of matched asymptotic expansions to exploit the asymptotic structure to reduce the problem to one that is asymptotically related but easier to solve than the original. In general, one must solve this related problem numerically. The hybrid related problem contains the entire infinite logarithmic expansion in its solution, which circumvents the task of obtaining each coefficient in successive terms individually, as one would have to do using only the method of matched asymptotic expansions. Since the hybrid solution essentially sums the infinite expansion of reciprocal logarithms, the error of the approximation is smaller than any power of ( — 1/loge). A n important feature of the hybrid related problem is that it is non-stiff, which means that it does not suffer from the difficulty of applying full numerics to the original problem of resolving the rapidly varying scale structure. Another advantage of the hybrid method solution is that the parameter dependence of the problem is reduced from that of the original. Using the hybrid method solution, one can compute an e-curve of the solution for a given number and location of removed subdomains, where £ is a measure of the subdomain cross-sectional size. The shape of the subdomain enters into the hybrid method solution through a single parameter d — d(b), where b is the coefficient of the Dirichlet component of the subdomain boundary condition. In some problems, the parameter dependence is further reduced with e and d occurring only in terms of their product, (ed). This exploits Kaplun's equivalence principle, which states that there is an asymptotic equivalence between subdomains of different cross-sectional shape, based on an effective radius, (ed), of the cross-section. The reduction in parameter dependence means that the hybrid method solution is less computationally intensive than a full numerical solution, £ 0. (8.1) 113 Chapter 8. Conclusions in which one would have to restructure the solution grid for each change in size or shape of the subdomains. We have shown that singular perturbation problems containing infinite logarithmic expansions arise in a wide variety of contexts. We dedicated four chapters of this thesis to the detailed application of the hybrid method to such singular perturbation problems occurring in fluid flow in a straight pipe with a core, skeletal tissue oxygenation from capillary systems, heat transfer convected from small cylindrical objects, and low Reynolds number fluid flow past a cylinder that is asymmetric to the uniform free-stream. In the previous chapter, we remarked on possible extensions to the general framework of ap-plicable problems. We briefly discussed applications in black body radiation, multi-body low Reynolds number fluid flow, vibration of thin plates with small holes or concentrated masses, localized non-linear reactions on catalytic surfaces, low frequency scattering of light, diffusion of a chemical species out of an almost impervious container, and steady-state current flow from microelectrodes. A l l told, we have drawn attention to hybrid method applications on problems that are steady or unsteady, linear or non-linear, second-order or fourth-order, and eigenvalue problems. Although we are at the end of this thesis, there are many more applications of the hybrid asymptotic-numerical method for future studies. 114 Bibliography [1] M . Abramowitz and I. A . Stegun. Handbook of Mathematical Functions. Dover Publica-tions, Inc., New York, 1965. [2] G. Batchelor. Slender-body theory for particles of arbitrary cross-section in Stokes flow. Journal of Fluid Mechanics, 44:419-440, 1970. [3] F . P. Bretherton. Slow viscous motion round a cylinder in a simple shear. Journal of Fluid Mechanics, 12:591-613, 1962. [4] W . Chester. The matching condition for two-dimensional flow at low Reynolds number: with application to the motion of an elliptic cylinder. Proceedings of the Royal Society of London, Series A, 437:195-198, 1992. [5] P. A . Clark, S. P. Kennedy, and A . Clark Jr. Buffering of muscle tissue P O 2 levels by the superposition of the oxygen field from many capillaries. In Oxygen Transport to Tissue XI, volume 248 of Advances in Experimental Medicine and Biology, pages 165-174. Plenum Press, New York and London, 1989. [6] S. M . Doherty and E . H . Dowell. Experimental study of asymptotic modal analysis ap-plied to a rectangular plate with concentrated masses. Journal of Sound and Vibration, 170(5):671-681, 1994. [7] A . Dutta and A . S. Popel. A theoretical analysis of intracellular oxygen diffusion. Journal of Theoretical Biology, 176:433-445, 1995. [8] F . E . Eastep and F . G. Hemmig. Estimation of fundamental frequency of non-circular plates with free, circular cutouts. Journal of Sound and Vibration, 56(2):155-165, 1978. [9] M . L . Ellsworth, A . S. Popel, and R. N . Pittman. Assessment and impact of heterogeneities of convective oxygen transport parameters in capillaries of striated muscle: experimental and theoretical. Microvascular Research, 35:341-362, 1988. [10] J . E . Fletcher. Mathematical modeling of the microcirculation. Mathematical Biosciences, 38:159-202,1978. [11] J . E . Fletcher. On facilitated oxygen diffusion in muscle tissues. Biophysical Journal, 29:437-458,1980. [12] N . H . Fletcher and T. D. Rossing. The Physics of Musical Instruments. Springer-Verlag, New York, 1991. Chapter 3: two-dimensional systems: membranes and plates. [13] N . A . Frankel and A . Acrivos. Heat and mass transfer from small spheres and cylinders freely suspended in shear flow. The Physics of Fluids, 11(9):1913—1918, 1968. [14] P. R. Garabedian. Partial Differential Equations. John Wiley & Sons, Inc., New York, 1964. 115 Bibliography L. Greengard, M . C. Kropinski, and A . Mayo. Integral equation methods for Stokes flow and isotropic elasticity in the plane. Journal of Computational Physics, 125(2):403-414, 1996. P. Hemker and J . Miller, editors. Numerical Analysis of Singular Perturbation Problems, London, 1979. Academic Press Inc. L. Hoofd. Updating the Krogh model - assumptions and extensions. In Oxygen Transport in Biological Systems, Society for Experimental Biology Seminar Series 51, pages 197-229. Cambridge University Press, Cambridge, 1992. L . Hoofd. Calculation of oxygen pressures in tissue with anisotropic capillary orientation. I. Two-dimensional analytical solution for arbitrary capillary characteristics. Mathematical Biosciences, 129:1-23, 1995. G. Hsiao and R. C. MacCamy. Solution of boundary value problems by integral equations of the first kind. SIAM Review, 15(4):687-705, 1973. G . C. Hsiao. Singular perturbation of an exterior Dirichlet problem. SIAM Journal of Mathematical Analysis, 9(1):160-184, 1978. R. Hsu and T. W . Secomb. A Green's function method for analysis of oxygen delivery to tissue by microvascular networks. Mathematical Biosciences, 96:61-78, 1989. I. Imai. On the asymptotic behaviour of viscous fluid flow at a great distance from a cylindrical body, with special reference to Filon's paradox. Proceedings of the Royal Society, 208:487-516,1951. Al . S. Ingber, A . L . Pate, and J . M . Salazar. Vibration of a clamped plate with concentrated mass and spring attachments. Journal of Sound and Vibration, 153(1):143-166, 1992. S. Kaplun. Low Reynolds number flow past a circular cylinder. Journal of Mathematics and Mechanics, 6(5):595-603, 1957. J . B . Keller and M . J . Ward. Asymptotics beyond all orders for a low Reynolds number flow. Journal of Engineering Mathematics, 30:253-265, 1996. J . Kevorkian and J . Cole. Multiple Scale and Singular Perturbation Methods, volume 114 of Applied Mathematical Sciences. Springer, 1996. R. E . Khayat and R. G. Cox. Inertia effects on the motion of long slender bodies. Journal of Fluid Mechanics, 209:435-462, 1989. A . Krogh. The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue. Journal of Physiology (London), 52:409-415,1919. M . C. Kropinski, M . J . Ward, and J . B . Keller. A hybrid asymptotic-numerical method for calculating low Reynolds number flows past symmetric cylindrical bodies. SIAM Journal of Applied Mathematics, 55:1484-1510, 1995. 116 Bibliography C. G. Lange and H . J . Weinitschke. Singular perturbations of elliptic problems on domains with small holes. Studies in Applied Mathematics, 92:55-93, 1994. C. Leal and J . Sanchez-Hubert. Perturbation of the eigenvalues of a membrane with a concentrated mass. Quarterly of Applied Mathematics, 47(1):93-103, 1989. L . G . Leal. Laminar Flow and Convective Transport Processes: Scaling Principles and Asymptotic Analysis. Series in Chemical Engineering. Butterworth-Heinemann, Stoneham, M A , 1992. S. H . Lee and L. G. Leal. Low-Reynolds-number flow past cylindrical bodies of arbitrary cross-sectional shape. Journal of Fluid Mechanics, 164:401-427, 1986. S. J . Lighthill. Mathematical Biofluiddynamics. Regional conference series in applied mathematics. S I A M , Philadelphia, 1975. S. K . Lucas, R. Sipcic, and H . A . Stone. A n integral equation solution for the steady-state current at a periodic array of surface microelectrodes. SIAM Journal of Applied Mathematics, 1997. MathWorks, Inc. Partial Differential Equations Toolbox, User's Guide. M A T L A B . The MathWorks, Inc., Natick, Mass., 1996. S. Ozawa. Singular variation of domains and eigenvalues of the Laplacian. Duke Mathe-matical Journal, 48(4):767-778, 1981. A . B . Peirce and H . Rabitz. Modelling the effect of changes in defect geometry on chemically active surfaces by the boundary element technique. Surface Science, 202:32-57, 1988. J . R. Philip, J . H . Knight, and R. T. Waechter. Unsaturated seepage and subterranean holes: Conspectus, and exclusion problem for circular cylindrical cavities. Water Resources Research, 25(l):16-28, 1989. N . A . V . Piercy, M . S. Hooper, and H . F . Winny. Viscous flow through pipes with cores. Philosophical Magazine, 15(7):647-676, 1933. A . S. Popel. Analysis of capillary-tissue diffusion in multicapillary systems. Mathematical Biosciences, 39:187-211, 1978. A . S. Popel. Theory of oxygen transport to tissue. Critical Reviews in Biomedical Engi-neering, 17:257-321,1989. L . Prandtl. Uber Flussigkeiten bei sehr kleiner Reibung. International Math. Kongress, Heidelberg, 111:484-491, 1905. Motions of fluids with very little viscosity (English transla-tion), Tech. Memo., N . A . C . A . , No. 452, 1928. I. Proudman and J . R. A . Pearson. Expansions at small Reynolds numbers for the flow past a sphere and a circular cylinder. Journal of Fluid Mechanics, 2:237-262, 1957. T. Ransford. Potential Theory in the Complex Plane, volume 28 of London Mathematical Society Student Texts. Cambridge University Press, 1995. 117 Bibliography L. A . Romero. Low or high Peclet number flow past a sphere in a saturated porous medium. SIAM Journal of Applied Mathematics, 54(1):42-71, 1994. H . - G . Roos, M . Stynes, and L. Tobiska. Numerical Methods for Singularly Perturbed Dif-ferential Methods: Convection-Diffusion and Flow Problems, volume 24 of Computational Mathematics. Springer, 1996. M . H . Ross, E . J . Reith, and L. J . Romrell. Histology: A Text and Atlas, chapter 10, page 204. Williams & Wilkins, Baltimore, Maryland, 2 edition, 1989. K . Shintani, A . Umemura, and A . Takano. Low-Reynolds-number flow past an elliptic cylinder. Journal of Fluid Mechanics, 136:277-289, 1983. L . A . Skinner. Generalized expansions for slow flow past a cylinder. Quarterly Journal of Mechanical and Applied Mathematics, 28:333-340, 1975. C. A . Swanson. Domain perturbations of the biharmonic operator. Canadian Journal of Mathematics, 17:1053-1063, 1965. R. C. Tai and H. -K . Chang. Oxygen transport in heterogeneous tissue. Journal of Theo-retical Biology, 43:265-276, 1974. M . S. Titcombe and M . J . Ward. Convective heat transfer past small cylindrical bodies. Studies in Applied Mathematics, 91:81-105, 1997. M . S. Titcombe and M . J . Ward. A n asymptotic study of oxygen transport from multiple capillaries to skeletal muscle tissue. Submitted to SIAM Journal on Applied Mathematics, July 1998. V . Twersky. Acoustic bulk parameters in distributions of pair-correlated scatterers. Journal of the Acoustical Society of America, 64(6):1710-1719, December 1978. M . Van Dyke. Perturbation Methods in Fluid Mechanics. Parabolic Press, Stanford, 1975. M . J . Ward, W . D. Henshaw, and J . B . Keller. Summing logarithmic expansions for sin-gularly perturbed eigenvalue problems. SIAM Journal of Applied Mathematics, 53(3):799-828,1993. [58] A . J . Ward-Smith. Internal Fluid Flow. Clarendon Press, Oxford, 1980. 118
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Applications of a hybrid asymptotic-numerical method...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Applications of a hybrid asymptotic-numerical method for certain two-dimensional singular perturbation… Titcombe, Michele Susanne 1999
pdf
Page Metadata
Item Metadata
Title | Applications of a hybrid asymptotic-numerical method for certain two-dimensional singular perturbation problems |
Creator |
Titcombe, Michele Susanne |
Date Issued | 1999 |
Description | This thesis examines a hybrid asymptotic-numerical method for treating two-dimensional singular perturbation problems whose asymptotic solution involves reciprocal logarithms of the small perturbation parameter, ε, in the form [equation]. The purpose of this hybrid asymptotic-numerical method is to treat the slow convergence problems of asymptotic expansions of this form. For the applications that we consider in this thesis, we believe there is sufficient evidence of convergence of these expansions for small enough ε. The hybrid method uses the method of matched asymptotic expansions to exploit the asymptotic structure to reduce the problem to one that is asymptotically related to the original. In general, one must solve this related problem numerically. The hybrid related problem contains the entire infinite logarithmic expansion in its solution, thus removing the necessity of obtaining each coefficient in successive terms individually, as one would have to do using only the method of matched asymptotic expansions. The hybrid solution essentially sums the infinite expansion of reciprocal logarithms and in so doing, improves the accuracy of the solution since the error of the approximation is smaller than any power of ( -1/log ε). An important feature of the hybrid related problem is that it is non-stiff. Thus, it does not suffer from the difficulty in solving the original problem numerically of resolving the rapidly varying scale structure. Another advantage of the hybrid method solution is that the parameter dependence of the problem is reduced from that of the original. The reduction in parameter dependence means that the hybrid method solution is less computationally intensive than a full numerical solution. We show that singular perturbation problems containing infinite logarithmic expansions arise in a wide variety of contexts. Four chapters of this thesis are dedicated to the detailed application of the hybrid method to such singular perturbation problems occurring in fluid flow in a straight pipe with a core, skeletal tissue oxygenation from capillary systems, heat transfer convected from small cylindrical objects, and low Reynolds number fluid flow past a cylinder that is asymmetric to the uniform free-stream. Following the detailed analysis of these four problems, we remark on possible extensions to the general framework of applicable problems. For example, we discuss applications in black body radiation, multi-body low Reynolds number fluid flow, vibration of thin plates with small holes or concentrated masses, localized non-linear reactions on catalytic surfaces, low frequency scattering of light, diffusion of a chemical species out of an almost impervious container, and steady-state current flow from microelectrodes. |
Extent | 8275352 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080051 |
URI | http://hdl.handle.net/2429/10036 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1999-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1999-389901.pdf [ 7.89MB ]
- Metadata
- JSON: 831-1.0080051.json
- JSON-LD: 831-1.0080051-ld.json
- RDF/XML (Pretty): 831-1.0080051-rdf.xml
- RDF/JSON: 831-1.0080051-rdf.json
- Turtle: 831-1.0080051-turtle.txt
- N-Triples: 831-1.0080051-rdf-ntriples.txt
- Original Record: 831-1.0080051-source.json
- Full Text
- 831-1.0080051-fulltext.txt
- Citation
- 831-1.0080051.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080051/manifest