UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
EXTENSIONS AND REFINEMENTS TO THE COMPLEX VARIABLE BOUNDARY ELEMENT METHOD INCLUDING ITS APPLICATION TO NUMERICAL GRID GENERATION By ROBERT TIMOTHY BAILEY A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1991 ACKNOWLEDGMENTS I would first like to thank my committee chair, Dr. C. K. Hsieh, for his continual guidance and encouragement during the course of my doctoral work. His dedication to excellence in both the classroom and the laboratory has set an example that I can only hope to one day follow. Much appreciation goes out to the four remaining members of my supervisory committee, not only for their part in this, my last academic effort as a student, but also for their continual influence during my time here at the University. Dr. Gater's particular brand of Thermodynamics and his "master scratch sheet" have followed me throughout my undergraduate and graduate studies. Dr. Selfridge's unique approach to computer programming has permanently altered my way of thinking about problems (although I'm still not an APL convertmuch to his chagrin). I have enjoyed many conversations with Dr. Lear on a wide range of topics, and I hope that I can someday have the same kind of rapport with students. Dr. Roan has earned my admiration for his superior teaching and my gratitude for the opportunity of working with him on a very interesting ejector analysis project. Additional thanks goes to my fellow students, especially Don Kalinich and Roy Johannesen, for their biting wit and sparkling lunchtime conversations. My parents, family, and friends deserve thanks for their continued interest in my progress (even though I will not miss answering "What exactly is it that you're doing?" for the 1200th time). Finally, my biggest thanks and my love goes to my wife, Patti, without whom none of this would have been possible. She has supported me (both emotionally and financially) through much of my Ph.D. work, and I hope that I can someday return the favor. TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................... ii ABSTRACT ........................................... vi CHAPTERS I INTRODUCTION AND LITERATURE REVIEW .......... 1 1.1 The Corer Problem ................... .......... 3 1.2 HigherOrder Elements ........................ 11 1.3 Numerical Grid Generation ........................ 11 1.4 Objectives ................................. 16 1.5 Overview .................................. 16 II DESCRIPTION OF THE LINEARELEMENT CVBEM ..... 18 2.1 Boundary Conditions ............................ 23 2.2 Solution Methods ............................ 25 2.3 Calculation of the Normal Gradient of the Potential ...... 33 2.4 Additional Comments and the Psi Reference Node ....... 34 III FORMULATION OF THE CVBEM USING QUADRATIC ELEMENTS ....................... 36 3.1 Boundary Conditions and Solution Methods ............ 42 3.2 Additional Comments ......................... 48 IV TREATMENT OF CORNERS IN THE CVBEM ........... 49 4.1 BCI Nodes ........................... 49 4.2 Corner Nodes and DoublySpecified BCI Nodes ......... 56 4.3 Rules for the Nine Corner Cases .. ................ 59 V NUMERICAL GRID GENERATION USING THE CVBEM ................................. 65 5.1 Derivation of the Inverted CVBEM Equations ......... 65 5.2 The CVBEM Numerical Grid Generation Algorithm ...... 69 5.3 Additional Comments on the CVBEGGM ............ 82 VI EXAMPLES AND RESULTS ......................... 85 6.1 Quadratic Element Examples and Results ............. 85 6.2 Corner Treatment Examples and Results ............. 91 6.3 CVBEGGM Examples and Results .................. 105 6.4 Hardware and Software .......................... 112 VII CONCLUSIONS AND RECOMMENDATIONS ........... 117 APPENDICES A FORMATION OF THE CVBEM MATRIX EQUATIONS ........................ 121 A.1 Explanation of the Matrix Equation Variables .......... 121 A.2 The Explicit Method Matrix Equation ................ 123 A.3 The Implicit Method Matrix Equation ............... 125 A.4 The Hybrid Method Matrix Equation ................ 126 B DERIVATIONS FOR THE QUADRATIC ELEMENT CVBEM .................................... 129 B.1 Derivation of the QuadraticElement CVBEM InteriorPoint Equation ........................ 129 B.2 Derivation of the QuadraticElement CVBEM NodalPoint Equations ......................... 132 REFERENCES ........................................... 146 BIOGRAPHICAL SKETCH ................................. 151 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy EXTENSIONS AND REFINEMENTS TO THE COMPLEX VARIABLE BOUNDARY ELEMENT METHOD INCLUDING ITS APPLICATION TO NUMERICAL GRID GENERATION By Robert Timothy Bailey May, 1991 Chairman: Dr. C. K. Hsieh Major Department: Mechanical Engineering The complex variable boundary element method (CVBEM) for simply connected domains is extended in three new areas. A complete formulation of the CVBEM using quadratic elements is presented. The derivation follows the format for linear elements given in the literature, and both the nodal and interiorpoint equations are given in detail. Second, a methodology for the treatment of covers in the CVBEM is developed for all possible combinations of Dirichlet, Neumann, and Robin boundary conditions at the corner. It is shown that no special augmentation of the CVBEM matrix equations is necessary when one follows the guidelines set forth in this work. Finally, an implicit, iterative variation of the linearelement CVBEM is applied to numerical grid generation in twodimensional, simplyconnected spatial domains. The quadraticelement formulation is successfully tested by solving example problems with available analytical solutions. Analytical solutions are also used to verify the corner handling methodology, and the CVBEM results are compared with those from a real variable boundary element method for accuracy. The application of the CVBEM to numerical grid generation is demonstrated by generating grids for several irregular geometries. The quality of the resulting grids is assessed by visual inspection and by examining the distribution of the metrics. CHAPTER I INTRODUCTION AND LITERATURE REVIEW Over the past ten years, a new class of numerical methods for solving partial differential equations (PDEs) has risen to challenge the wellestablished finitedifference methods (FDMs) and finiteelement methods (FEMs) [1,2]. Unlike the FDMs and FEMs which require that the entire domain be broken down into grid points or elements (discretized), these relatively new methods, known as boundary element methods (BEMs), only involve discretization of the boundaries of the domain. In this way, the BEMs effectively reduce the dimension of the discretization by one, and in the process, significantly decrease the number of simultaneous equations which need to be solved. These benefits do not come free, however, as the matrices associated with these systems of equations tend to be nearly fully populated in the BEMs, while they are comparably sparse in the FEMs and the implicit FDMs. Even with this tradeoff, it has been shown that the BEMs are often much more computationally efficient than the FDMs or FEMs [16]. One way of classifying the different kinds of BEMs is according to the variables used in their formulation. Traditionally, the BEMs have been formulated using real variables (RVBEMs), but a recent and powerful advance involves the use of complex variables, the result being a method known as the complex variable boundary element method (CVBEM) [7]. This new method is limited in application to Laplace's and Poisson's equations, and is further limited to twodimensional (2d) domains, but it does possess the following two significant advantages over the RVBEMs: 1) approximations are made only at the boundary of the domain, and 2) all integration are carried out analytically without the need for numerical integration. These two factors combine to give the CVBEM excellent potential for both accuracy and efficiency. The origins of the CVBEM can be traced to a paper by Hunt and Isaacs [8] where it was applied to the solution of groundwater flow problems. Hromadka and Guymon [9] subsequently used the method to predict freezing fronts in soils and then formalized the rather loose development into the CVBEM [10]. The method was also shown by Hromadka [11] to be a generalization of the analytic function method of Van Der Veer [12]. In the process of refining the CVBEM, Hromadka [13] developed a technique for error visualization at the domain boundary and determined relative error bounds for the method [14]. He also considered the use of variable trial functions for improving accuracy [15] and investigated the proper placement of collocation points on the domain boundary [16]. Various physical phenomena have been modeled by the CVBEM including advective and groundwater contaminant transfer [8,17,18], conduction heat transfer [16], prediction of freezing fronts in soil [9,19], and stratified flows [20]. All of these applications and developments have been detailed in a book on the CVBEM by Hromadka and Lai [7]; however, more recent developments in the CVBEM (post 1987) are not covered in this book. The CVBEM has since been extended from simplyconnected to doubly and multiplyconnected domains by Kassab [21], Kassab and Hsieh [22], and Hsieh and Kassab [23]. It was found that the complex potential along cuts in the domain does not cancel out but results in a complex stream function that plays the role of a perturbation in the nodal equations. About the same time, Harryman et al. [24] and Hromadka [25] applied the CVBEM to specific multiplyconnected domains where such perturbations did not appear due to a zeroflux condition on the interior boundaries. Very recently, Mokry [26] has applied the CVBEM to external potential flows. Examples were given for flows over airfoils whose exact solutions are known. In all cases, the agreement of the CVBEM results with the theory was excellent. This dissertation is concerned with extending the CVBEM in three areas, namely, boundary corner treatment, higherorder elements, and numerical grid generation. Each of these topics will be discussed in the next three sections. 1.1 The Corner Problem An interesting problem common to many BEMs is the proper treatment of the corners or edges in the boundaries of the domain. A corner (in 2d) or edge (in 3d) is defined as a point or locus of points where the tangent to the boundary possesses a sharp discontinuity. This discontinuity must be due solely to the problem geometry and not due to any discretization scheme associated with the BEM. At such a corner or edge, quantities associated with the direction normal to the boundary are doublevalued, and the conventional RVBEMs can produce systems of equations which are underconstrained. In 2d potential problems, the three traditional boundary conditions are 1) the Dirichlet condition where 0 is known and is unknown, 2) the Neumann On condition where is known and 0 is unknown, and 3) the Robin (or convective) condition where both 6 and are unknown. The symbol 4 represents the on potential function, while n represents the direction normal to the boundary. The application of these three conditions logically leads to nine possible boundary condition combinations at any corner; see Table 11. Here, the subscripts U and D refer to the "upstream" and "downstream" sides of the corner, respectively, with the direction of travel along the boundary always placing the domain to the TABLE 11. The nine traditional boundary condition combinations at a corner node in a potential problem. Upstream Nodes Condition Ouantity Specified Dirichlet Neumann (')u Neumann (a)u Dirichlet 0 Robin none Robin none Neumann (a)u Robin none Dirichlet 0 Downstream Nodes Condition Quantity Specified Neumann ()D Dirichlet Neumann ( D)O Robin none Dirichlet 4 Robin none Robin none Neumann ()D Dirichlet Case 1 2 3 4 5 6 7 8 9 I left. Also, the symbol s will be used to represent the direction tangent to the boundary. For the optimum situation where both the upstream and downstream conditions are available at the corner node, the unknowns and corresponding number of equations which can be applied in ordinary RVBEM formulations are given in Table 12. In all nine cases, one boundary element nodal equation is always available for application. Thus, with only one unknown, Cases 1 through 3 pose no problem. Cases 4 and 5 have two unknowns, but in both cases an extra equation is available from the Robin condition. Case 6 has three unknowns, but two additional equations are available from the Robin conditions. Cases 7 and 8 have two unknowns each, but once again, the Robin condition can be used as an extra equation in either case. Finally, Case 9 has two unknowns, but no additional equation is available for application. This renders Case 9 underconstrained, and it is here that investigators have been forced to seek methods for special treatment of corners in the RVBEMs. Probably the simplest way to treat Case 9 is to assume that Uj and (jD are equal at the corner. This approach was used by Lachat and Watson [27] in the early treatment of problems in 3d elastostatics. They found that the resulting errors were confined mostly to the region of the domain near the corner. Brebbia and Dominguez [28] also investigated this method with regard to 2d potential problems but found that, for their cases, the errors in the gradients near the corners were unacceptable. They then attempted to use two nodal points in close proximity, one on either side of the corner (binodes), concluding that the results using this approach were generally better but were still subject to large errors. Ricardella [29] used essentially the same "binode" approach in his treatment of problems in elasticity, but he considered only step changes in tractions (Case 3), or traction/displacement combinations (Cases 1 and 2), avoiding Case 9 altogether. TABLE 12. The unknowns and corresponding number of equations available in the standard RVBEMs for the nine possible corer node boundary condition combinations listed in Table 11. Number of Case Unknowns 4an. Available 1 (1 2 ()D 1 3 1 4 ( )v and ()D 2 5 ( )U and ()D 2 6 (n)U, ( )D and 2 3 7 ( )D, and 2 8 ()U, and 2 9 ( ) and ( )D 1 Alarcon et al. [30] proposed a different approach based on the fact that, in principle, knowledge of the potential distribution along the boundary allows calculation of the gradient there. By assuming a linear variation of L and aL over the segments connecting the corner node with its adjacent nodes, each of the normal gradients is expressed in terms of a single "new" unknown, V . After solution, the values of ( and (D are calculated from V The drawback with this method is that as the angle turned by the tangent to the boundary decreases, the approximation to V 4 becomes less and less accurate. Paris et al. [31] later proposed calculating the values of F and 'JD directly from the values of the potential at the corner and coreradjacent nodes. This is done before the BEM solution process is even begun. The equations associated with the corner nodes are thus completely eliminated from the BEM system of equations; however, such direct approximation of the gradients is subject to truncation errors and, as before, the approximations to the normal gradients become unreliable as the corer becomes less sharp. Walker and Fenner [32] carried the "gradient approximation" line of thought in a slightly different direction by combining approximate expressions for and into a single equation relating the two. This equation is then used as the "extra" equation needed to properly constrain Case 9. The drawback associated with this method is that the equation relating the normal gradients is sometimes inaccurate, with special treatment being necessary when the tangent to the boundary turns through an angle near ninety degrees or when the tangential gradients are nearly zero and the angle turned by the tangent is less than thirty degrees. Grilli et al. [33] have reported good results using a similar "gradient approximation" approach in their investigation of nonlinear waves, as have Bruch and Lejeune [34] in their study of 2 and 3d potential problems. The latter's treatment is somewhat problem dependent, however, since they state that "the additional equation (for Case 9) will vary from problem to problem." Yet another approach was investigated by Gray [35], who used the so called hypersingular equation (obtained by differentiating the traditional boundary integral equation) as the "extra" equation needed in Case 9. Reasonable results were obtained for practical problems involving the simulation of a 3d electrochemical plating process. The drawback associated with this method is the careful treatment required for proper evaluation of the integrals in the hypersingular equation. A completely different line of reasoning involves the use of discontinuous (or nonconforming) boundary elements at a corner. A discontinuous element is one in which the collocation points are removed from the element perimeter and do not coincide with the geometric nodes used to prescribe the domain geometry. Using this approach, a collocation point is placed on either side of a corner with the customary single geometric node still located at the corner (see Fig. 11). The corner geometry is thus retained, while the use of two collocation points permits the writing of two boundary element nodal equations. This allows the normal gradient of the potential to be doublevalued at the corner. One drawback associated with this approach is that since no collocation node actually exists at the corner, the 0 and values must be interpolated from values at the interior collocation points. This introduces some inaccuracies at the corners. Also, the number of simultaneous equations is increased by the use of discontinuous elements since two nodes are used near each corner as opposed to one node for continuous elements. Patterson and Sheikh [36] successfully applied this method to 2d harmonic and 3d potential problems, and Danson et al. [37] have made use of discontinuous elements in their general purpose BEASY boundary element program with good results. However, Manolis and Banerjee 9 separate collocation points ordinary node corner geometric node FIGURE 11. Using "discontinuous" elements at a boundary corner. [38] compared the use of continuous (conforming) and discontinuous (non conforming) boundary elements and found that, in many cases, the use of discontinuous elements lead to large errors in the computed solution. They attributed these errors to the fact that discontinuous element collocation nodes often do not lie on the actual domain boundary, concluding that, "In general, conforming (continuous) elements yield more accurate results than non conforming (discontinuous) elements." This conclusion was subsequently challenged by Brebbia and Niku [39] who suggested that the errors attributed to the use of discontinuous elements were actually the result of poor numerical methodsa charge subsequently denied by Manolis and Banerjee. It appears that the use of discontinuous elements is still under debate in BEM circles; nevertheless, the BEASY code continues to use them with success. Recently, Detournay [40] has investigated the handling of corer singularities using a complex variable integral equation method similar to the CVBEM. In that paper, special corer elements characterized by two intersecting straight line segments were used to handle problems where the gradient of the potential had a constant direction along each of the segments. By using these elements, it is intended that the special behavior of the potential function in the vicinity of a corner can be more accurately modeled. A drawback with this approach is that, in many cases, such an element yields a nonlinear nodal equation. By formulating 2d potential problems with complex variables instead of real variables, the CVBEM avoids the problem of gradient discontinuity at boundary corners. Instead of solving for the doublevalued gradients at a corner node, one now solves for a singlevalued stream function. Still, Neumann and Robin conditions can lead to special circumstances at covers, and a methodology for the treatment of covers in the CVBEM is needed. 1.2 HigherOrder Elements The earliest BEM formulations approximated the domain boundary as a series of linear segments over which the value of the variables of interest (temperature, displacement, etc.) were constant [41,42]. These segments are known as constant elements (see Fig 12a). The collocation points where the unknown values are considered in the BEM solution are known as nodes or nodal points, and, for constant elements, they are located right at the center of each element. Constant elements provided reasonable solutions in many cases, but when greater accuracy was required, investigators used more refined boundary approximations such as linear, quadratic, or cubic elements [1,2,28,29]. Linear elements have nodes at both ends (see Fig 12b), and the variables of interest are assumed to vary linearly over each element. Quadratic and cubic elements are referred to as higherorder elements, and they allow the curved structure of boundaries to be retained. As their names imply, quadratic elements pass a seconddegree polynomial through three successive nodal points, while cubic elements involve thirddegree polynomials and four nodal points (see Figs. 12c and 12d). Hromadka [7] has laid the groundwork for the use of higherorder elements in the CVBEM via generalized proofs, but linear elements remain the highest order elements that have been reported to date in the CVBEM literature. 1.3 Numerical Grid Generation It was mentioned earlier that finitedifference methods (FDMs) are widely used in the solution of PDEs, especially in the field of computational fluid dynamics (CFD) [43,44]. All FDMs require that the continuous domain of interest be replaced by a discrete domain composed of a collection of points distributed throughout the original continuous domain. These points are known .nodes constant element linear element nodes a) constant elements quadratic element b) linear elements middle nodes end nodes c) quadratic elements d) cubic elements FIGURE 12. Different types of boundary elements. cubic element nodes as grid points, and the collection of grid points is known as a grid system. Proper placement of these grid points is essential if accurate solutions are to be obtained. Often, nonuniform distributions are necessary in order to accommodate irregular geometries or complex flow conditions. Unsteady problems may even require grid points which move over time. Unfortunately, it is extremely difficult and impractical to derive finitedifference equations (FDEs) with nonuniformly distributed and moving grid points. For this reason, the grid points in the "physical" spatial domain are mapped onto a "transformed" domain where they are stationary and uniformly distributed (see Fig. 13). In practice, it is easier to perform the mapping in the opposite direction, and this process is known as grid generation, i.e., the generation of grid points in the "physical" spatial domain [45]. Methods for generating grid systems are traditionally divided into two major classesdifferential equation methods (including conformal mappings) and algebraic methods. Differential equation methods generate grids by solving one or more PDEs that describe the transformation between the "physical" spatial domain and the "transformed" domain. This usually requires a significant computational effort, but it produces grid systems whose grid lines are smooth and nonoverlapping. On the other hand, algebraic methods generate grid systems by interpolating between the boundaries of the "physical" spatial domain. Since the solution of PDEs is not involved, the algebraic methods tend to be much more computationally efficient than the differential equation methods. However, grid systems generated by algebraic methods may have grid lines which are nonsmooth and overlappingcharacteristics which must be corrected if meaningful finitedifference solutions are to be obtained. Thompson et al. [46] have written a book which thoroughly describes the techniques used in numerical grid generation. Also, two comprehensive survey !grid points grid lines a) the physical domain b) the transformed domain FIGURE 13. The physical and transformed domains used in numerical grid generation (2d case). papers on grid generation are listed as References 47 and 48 for the interested reader. Here, only those grid generation methods which are in some way connected to the CVBEM are reviewed. Since the CVBEM utilizes complex variables, it shares common ground with the conformal mapping techniques. Generally, these techniques have been restricted to particular geometries, but some of the more flexible approaches include those of Anderson [49], Davis [50], Halsey [51], and Harrington [52]. The CVBEM is set apart from even these more advanced conformal mapping methods in terms of its ability to handle irregular geometries due to its boundary element nature. For this reason, references which involve the application of BEMs to grid generation will now be reviewed. The integral equation method of Symm [42], mentioned earlier, has been used for mapping a simplyconnected 2d spatial domain onto a unit disk and for mapping doublyconnected region onto annular regions. Hayes et al. [53] subsequently improved this approach and were able to obtain more accurate results. In each of these applications, the primary concern was to obtain the mapping rather than to use the mapping for generating grid systems. Indeed, the only application of a BEM directly to grid generation appears to be that described by Adamczyk [54], who used an electrostatic analog to generate grids between bodies in a cascade. Force and potential lines calculated using a panel method solution of an integral equation were used to define the resulting grid system. Also, a doubly infinite line of alternating charges separated by a finite distance was used as the fundamental solution. It appears that the CVBEM has never been applied to numerical grid generation; however, its accuracy and efficiency make it an excellent candidate for bridging the gap between existing algebraic and differential equation methods. Since it involves the solution of the Laplace equation, it would be properly classified as a differential equation method (although the term integral equation method would be more appropriate), but it is expected to be much faster than the existing elliptic solvers. 1.4 Objectives Based on the background information presented in the previous three sections, the objectives of this dissertation are as follows: 1) Develop a formulation of the CVBEM using quadratic elements. Previously, linear elements were the highestorder elements used. 2) Devise a methodology for the proper treatment of boundary corners in the CVBEM. This topic has been treated extensively in the RVBEMs but has not been addressed in the CVBEM. 3) Derive a variation of the CVBEM which can be used for numerical grid generation. It is noted that the CVBEM has never been used for this purpose. 1.5 Overview In the next chapter, a description of the linear element CVBEM as formulated by Hromadka is presented together with some observations and suggestions regarding some of the method's more subtle points. Chapter III then extends the CVBEM via the use of quadratic elements. Chapter IV goes on to discuss the proper treatment of corners in both the linear and quadraticelement CVBEMs, while Chapter V presents a novel application of the CVBEM in the area of numerical grid generation. Examples and results which demonstrate the new techniques of the previous three chapters are contained in Chapter VI. Finally, in Chapter VII, conclusions are stated and recommendations for further work are given. 17 To present the material in a concise manner, all derivations are relegated to the Appendix. CHAPTER II DESCRIPTION OF THE LINEARELEMENT CVBEM When solving 2d potential problems, it is often convenient to make use of the theory of complex variables. Using this approach, one can construct a complex potential, w(z) = (zx,y) + io(z,y), where is the potential function and i) is the stream function. This complex potential is analytic so that 0 and f satisfy the CauchyRiemann conditions, 80 al an s8, (2.1) where s represents the streamline coordinate and n represents the coordinate normal to it; furthermore, one can apply the wellknown Cauchy integral formula which states that (zo) = ( zo e 1 zo P r (2.2) This expression relates the value of complex potential w at point zo located inside the kconnected Jordan domain, Q, to a contour integral (containing w) along the boundary, F. The direction of travel for the contour integral is such that the interior of the domain is always to the left. The linearelement formulation of the CVBEM transforms the Cauchy integral formula into a BEM by using two major approximations. First, boundary F is discretized into N finitelength segments (elements), rj, the end points of which are referred to as nodal points or nodes. Domain boundary r is then taken as the union of these elements as shown in Fig. 21, i.e., N r=u r, (2.3) j=1 Second, the function w(z) is approximated by a linear global trial function, GI(z), given by N G,(z) = E N,(z) wj (2.4) j=1 where w is the value of w at nodal point zj, and Nj(z) is a continuous basis function weighting the effect of wj over elements F I and F,. It is necessary to explain the notation that will be used throughout the remainder of this dissertation. The subscript e.ct will identify a quantity whose value is known exactly. The bar, , will refer to a quantity whose value is specified, such as in the boundary conditions. Finally, the hat, ^, will represent a quantity whose value is treated as an unknown in the solution by the CVBEM. Returning now to the description of the linear element CVBEM, basis function Nj(z) is taken as a first degree polynomial of the following form: N,(z) = ( +1Z) / (z+1z) z E ri 0 z r,_1 U r, By substituting G1(C) for w(C) in the right hand side of the Cauchy integral formula (Eqn. (2.2)), the firstorder CVBEM approximation of w can be FIGURE 21. 8(j+1,j;O) 8(i+1,i j z0 z z z 3 z z 2 3 .X Boundary discretization and angle definitions in the linear element CVBEM. expressed as (zo) = I 0C zo E zo r (2.5) r After much simplification (see [7]), the contour integral can be eliminated, and the above equation reduces to 1N hJ o(Zo) = 2 [.j+1(zo zj) Wj(Zo zj+,)] (2.6a) 2i = (zi +1 zj) where hj = In L(zj+ zo) = In (zi1z) + iO(j +1, j; 0) (2.6b) L (z o) J (z zo) This formula forms the basis for the linearelement CVBEM. (See Fig. 21 for the definition of O(j + 1, j; 0).) Eqn. (2.6), being expressed in complex variables, actually embodies two equationsone for the real part and one for the imaginary part. If the values of 0 and k (and thus w) are known at each boundary node, Eqn. (2.6) can be used to calculate and at any interior point, z0. In most potential problems, however, boundary conditions specify either 0 0, or neither of them explicitly. To solve for the unknown values of 0 and i, it is necessary to derive an extended version of Eqn. (2.6) by moving zo to the position of zj on the boundary. In this effort, one cannot simply plug in z3 for zo in Eqn. (2.6) because (zj zo) appears in the denominator of the natural log term. Instead, one takes the limit of Eqn. (2.5) as zo approaches zj. Hromadka [7] has performed this somewhat involved task, with the result zWi'l + 0(j+1 ^1 wj= wj wIn + (j + 1, j 1; j) + L [wi+1(z z ) j(z,zi+ )] hi/(zi+, z.) (2.7a) i,i+1 i j where h = In + j ( + i(i+1, i; j) (2.7b) F i (z Iz) The angles O(j + 1, j 1; j) and B(i + 1, i; j) are shown in Fig. 21. Eqn. (2.7) can be applied at any boundary node, but like Eqn (2.6), Eqn. (2.7) has real and imaginary parts. Two equations can thus be derived for arbitrary boundary node j as 1 Zf i+lzj, 27= 1 znx J + OjO(j+1, j1;j) + E {,+CA + i+ Ci ,,C4 iC3}} (2.8) i,i+lij and = In ; z(j+1,j1; j) + E { ,+1C i,+iC2 C, + AC, (2.9) where ii+ 1 C, = [(x x,)C (y, yi)D] (2.10a) C2 = [(xj x)D + (y y)C] (2.10b) C3 = [(xixi+i)C (yjy;+1)D] (2.10c) C4 = [(xjxi+)D + (y,yi+,)C] (2.10d) C = [A(xi+, xi) + B(y+,1 yi)] / F (2.10e) D = [B(xzi+,x) A(y,+1y)] / F (2.10f) F = (xi +1 )2 + (Yi +iYi)2 (2.10g) A = In ) (2.10h) I(z z') B = O(i + 1,i; j) (2.10i) Henceforth, Eqns. (2.8) and (2.9) will be referred to as the phi nodal and psi nodal equations, respectively. 2.1 Boundary Conditions Having generated the desired equations for calculating q and i at the boundary nodes, a question arises as to how one uses Eqns. (2.8) and (2.9) to solve for the unknown boundary values of and 1. Before this question can be answered, consideration must be given to the conditions that can exist at the boundary nodes. As it turns out, each of the three "classic" boundary conditions of the potential theory yields a different equation as detailed below [21]. Dirichlet Condition: The potential at node j is known and specified, i.e., j = e,,xtc, j (2.11) Notice that for such a condition tj is unknown. In keeping with the notation defined previously, the specified j becomes Neumann Condition: The normal gradient of the potential at node j is known. The stream function there is then related to the normal gradient of the potential as += 21f(1 + }1Izz .1 (2.12) Here, it is assumed that ( varies linearly from node j 1 to node j [21]. For this condition, both ,j and ,j are unknown and are referred to as fj and ,j. Robin Condition: The normal gradient of the potential and the potential itself are related by = H(O 0,). The stream function there is then related to the potential as j = i + {(H(O)), + (H(,)),_} \ zz.il\ (2.13) Here, it is assumed that HO and H,. vary linearly from node j 1 to node j [21]. As in the Neumann condition, both qj and k, are unknown and again become 4j and 7j. The use of complex variables in the CVBEM also allows for the occurrence of the stream function condition as follows: Stream Function Condition: The stream function at node j is known and specified, i.e., i = .ctj (2.14) For this condition, Oj is unknown. The specified kj is renamed 4', and the unknown 4, becomes j,. 2.2 Solution Methods Eqns. (2.8) and (2.9) will be used to estimate the unknown values of 0 and ? at the boundary nodes, but prior to the presentation of the methodology for doing so, a close examination of these equations is in order. It is clear from their format that qj and j on the left can be calculated by using the known values of 0 and 0 at all boundary nodes whose indices appear on the right. As discussed previously, some of the values of 4 and i at the boundary nodes are not specified by the boundary conditions. The methods for estimating the unspecified values of 0 and I (designated q and f) thus hinge on how these quantities are related to the specified quantities (4 and b) and on which of the equations ((2.8) or (2.9)) is used in the construction of the matrix equation for solution of the problem. At nodal point j where a Dirichlet or stream function condition is imposed, there are three methods of solving for the estimated nodal values of 0j or Sy. Explicit Method: For a Dirichlet condition imposed at node j, Oj is known (as 4j) but Cj is not. One thus sets aj = j = ', and 0i = si. The first setting is governed by the fact that a Dirichlet condition is specified; no estimation is thus needed for Oj. The second setting is made in order to estimate the unknown value of Oj. Without such a setting, there will be two unknowns (0j and jl) at the same nodal point, a situation which does not allow for a unique solution. Since the field theory predicts that the estimated ij, if error free, will be exactly equal to that imposed, equating ij and ij is certainly reasonable. In the explicit method, an effort is made to keep all of the unknowns on the right side of the equation. It is therefore impossible to use Eqn. (2.9) since ij on the left side is unknown for the Dirichlet condition specified. Thus, Eqn. (2.8) is used as S= y{ln + j (j +1, j1; j) + E {,i+c2 + +iC1 0,4 ,iC} (2.15) i=1 i,i+1 Ai This equation is referred to as the explicit phi nodal equation. For the explicit method, it is sufficient to use Eqn (2.15) to solve for l,. Eqn. (2.9) is dropped. For a stream function condition at node j, Cj is known (as 4j) but Oj is not. One thus sets j = i = 4i, and Oj = Cj. This time, for the explicit method of solution, Eqn (2.9) is used, and the nodal equation is = In/ z 0(j+l, j1;J) + {i,, C3 + Pc4} (2.16) i,i+ij This equation is referred to as the explicit psi nodal equation. Eqn. (2.8) is dropped for this node in the solution. For the explicit method, solution is effected by applying Eqn. (2.15) at each Dirichlet node and Eqn. (2.16) at each stream function node to form a system of simultaneous equations which is solved for the unknown values of 4 and 4. Implicit Method: For a Dirichlet condition imposed at node j, one sets ij = 4i, and ij = ij. However, unlike the explicit method, the implicit method involves a nodal equation where unknowns appears on both sides. Eqn. (2.9) is thus used as I='j '{in z + O1(jl, j1;j) + O {i+i i+xCi ,c3 + AC4 (2.17) i= 1 i, i +1 This equation is referred to as the implicit psi nodal equation. Eqn. (2.8) is dropped. Similarly, for a stream function condition specified at node j, one sets j = i, and ,i = j,. Eqn. (2.8) is now used as In = t lzIn z, + Ou(j+1, J1; j) + I i+1C2 + O+1C i.C4 iC3} (2.18) i,i + 1 This equation is referred to as the implicit phi nodal equation. Eqn. (2.9) is dropped. In a manner similar to the explicit method, solution by the implicit method is effected by applying Eqn. (2.17) at each Dirichlet node and Eqn. (2.18) at each stream function node to form a system of simultaneous equations which is solved for the unknown values of 4 and 6. Hybrid Method: For a Dirichlet condition imposed at node j, one sets 4j = 4j on the left side of Eqn. (2.8) and 4j' = 4j and #1 = Oj on the right sides of Eqns. (2.8) and (2.9) to obtain i= 0{ in zI Z + 1o(j+1, j1;j) N + i {+A + ,+ACc iC4 AiC} (2.19) i,i+1j j1 Iz z1 2r il zj U +(j +1, 1; j) + O {,+1C i+c2 ,iC3 + i4 (2.20) i,i+1j These equations are referred to as the Dirichlet hybrid phi nodal and psi nodal equations, respectively, since Eqn. (2.19) is explicit while Eqn. (2.20) is implicit. Notice that both Oj and Oi are unknown. For a stream function condition imposed at node j, one sets Oj = 0, on the left side of Eqn. (2.9) and 1 = ^j and Sj = 1 j ln/zj S= {,ln Zj + 1 0(j+ 1, j 1; j) + E {i,+iC2 + V,+ICI ,C, AC, (2.21) i,i+1 j 2 2 i1 z z O(j+1, j1; j) + f {.+1c1 i+1C2 OiCa + .C4} (2.22) i,i +l J These equations are referred to as the stream function hybrid phi nodal and psi nodal equations, respectively, since Eqn. (2.21) is implicit while Eqn. (2.22) is explicit. Again, both S, and tk are unknown. In the hybrid method, solution is effected by applying Eqns. (2.19) and (2.20) at each Dirichlet node and Eqns. (2.21) and (2.22) at each stream function node. The number of equations solved simultaneously for unknowns ( and t is thus doubled in the hybrid method as compared to the explicit and implicit methods. Table 21 gives a summary of the solution methods described above based on a node imposed with a Dirichlet condition. A close examination of this table reveals the strengths and weaknesses associated with each method. Attention is first directed to the column headed "Quantities Equated" in the table. It appears that the only seemingly minor difference between the explicit method and the implicit method is the relation between fj and ,j, i.e., the two are set equal in the explicit method but not in the implicit method. In fact, Eqn. (2.8) for calculating (j is dropped in the implicit method when solving for Oj (see the column headed "Equation Dropped"). However, this equation can still be used to advantage for estimating the overall accuracy of the solution. In this effort, qj is calculated using Eqn. (2.8), and the estimated Oj is compared with the specified Oj for error. A "small" error gives an indication of high accuracy in the overall solution. The above advantage does not apply to the explicit method. For a Dirichlet condition, the explicit method requires that Eqn. (2.9) be dropped. If after solution, Eqn. (2.9) is subsequently used to estimate the unknown Oj, no specified iji exists to compare it to. Worse still, errors in ik at the boundary nodes lead to errors in the calculated values of since is found by On' in TABLE 21. Summary of solution methods based on a node imposed with a Dirichlet condition. Method Known Unknown Quantities Equated Equation Used Equation Dropped Explicit j, j = j = ect (2.8) (2.9) Implicit j ij = =exc,,i (2.9) (2.8) Hybrid qj = tj = ,,ex<,j (2.8) and (2.9) none on the left side of Eqn. (2.8); j = j and = j on the right side of Eqns. (2.8) and (2.9) numerically differentiating 0 as will be discussed later. Experience has shown that the implicit method generally yields more accurate values of 1 than the explicit method. The hybrid method is essentially the combination of the explicit and implicit methods. As such, one might expect that it is the most accurate of the three methods developed, and this often (but not always) turns out to be the case. Table 21 shows that, in the hybrid method j = j on the left side of Eqn. (2.8); j = s on the right sides of Eqns. (2.8) and (2.9); j = i on the right sides of Eqns. (2.8) and (2.9). The fact that 4i is equated to 'j on the left side of the Dirichlet hybrid phi nodal equation (Eqn. (2.19)) renders this equation very similar to Eqn. (2.15) in the explicit method. Likewise, the Dirichlet hybrid psi nodal equation (Eqn. (2.20)) is much like Eqn. (2.17) in the implicit method. The differences lie in the fact that unknowns 4' and j are both present in the hybrid method, with Eqns. (2.19) and (2.20) both being applied at Dirichlet node j. Since is also specified, it may be compared after solution with the estimated si, thus offering the potential for error detection afforded by the implicit method. Unfortunately, as just noted, the hybrid method requires the application of two equations at each Dirichlet or stream function node, while the explicit and implicit methods require only one. This results in a doubling of the number of simultaneous equations associated with Dirichlet and stream function nodes. Be mindful that all of the above comments concerning the explicit, implicit, and hybrid methods have been made for the Dirichlet case; however, they are just as valid for the stream function case with the exchange of 4 with i. For nodes where Neumann or Robin conditions are specified, Eqn. (2.12) or Eqn. (2.13) will be used in conjunction with the implicit phi nodal equation, Eqn. (2.18), in the solution of unknowns. By applying Eqns. (2.12) and (2.18) at all Neumann nodes, Eqns. (2.13) and (2.18) at all Robin nodes, and either the explicit, implicit, or hybrid method equations at all Dirichlet and stream function nodes, a system of linear algebraic equations is obtained. This system can be represented in a matrix form as cf{I = R (2.23) Here, C is the square matrix of coefficients on the unknown values of 4 and ? and R is a vector of known constants. A detailed description of how C and R are formed is given in Appendix A. The system can be solved by any direct method such as Gaussian elimination. Once the unknown values of 4 and are found at each boundary node, one can use these values, together with the specified values, 4 and as the boundary nodal values needed by Eqn. (2.6) for calculating C at any interior point. The matrix formulation of Eqn. (2.23) provides another perspective as to the preference of the implicit method over the explicit and hybrid methods. In the explicit method, the nodal equations produce a coefficient matrix with zeros on the main diagonala very poor situation, numerically. By contrast, the implicit method produces a matrix with the largest element in each row located on the main diagonal, a situation conducive to obtaining good numerical results. The hybrid method suffers from the same numerical pitfall as the explicit method (zeros on the main diagonal) but to a lesser degree since implicit equations are also used to form the hybrid coefficient matrix. The primary reason for choosing the implicit method over the hybrid method is the undesirable increase in matrix size in the hybrid method due to the need for two equations at each Dirichlet or stream function node instead of just one. 2.3 Calculation of the Normal Gradient of the Potential After 0 and 0 have been evaluated at the boundary nodes and interior points, the values of b may be used to determine at the boundary nodes. Use is made of the CauchyRiemann conditions, an Fs' where n is now taken as the outward normal to r with s as the tangential coordinate along F, assumed positive in the direction of the contour integral in Eqn. (2.2). The partial differential (and thus ) can be approximated by using finitedifference formulae. For example, a secondorder accurate, three (go point backward finitedifference approximation to at node j is given by (a tyj[(sjsj 2) (SjSji)] ) j_ti(,s2j)2 s (s sj_)(s ii2)( i_ 2) Sj2(S  j_ 1)(2.24) (Si_Sj_1)(SSj_2)(Sj1_Si_2) The values of the tangential arc length, s, can be estimated as the cumulative sum of the lengths of the linear segments connecting each successive nodal point, or a more sophisticated polynomial or spline procedure may be used. Onesided formulae should be used at corner nodes and at any others where the normal gradients are different on either side of the node. Centraldifferencing can be used at all remaining nodes. Equations such as (2.24) can be greatly simplified if the nodal points are equally spaced along the boundary, and References 43 and 44 provide lists of some common finitedifference formulae which can be used when this is the case. 2.4 Additional Comments and the Psi Reference Node Two more points are worthy of note. First, nodal points which are located at corners can be considered as a special case which will be studied in great detail in Chapter IV. Second, no matter what boundary conditions exist, one must specify a reference value of 0 at some nodal point along the boundary in order to serve as the constant of integration in Eqn. (2.5). Numerically speaking, the associated matrices become singular if no 0 value is provided. This anchor node will be referred to as the "psi reference node". Numbering it as node j, one has 4'j = '"exact,j The implicit phi nodal equation (Eqn. (2.18)) can then be used to determine the unknown ~j as was discussed earlier in the description of the implicit method. If, however, the value of .exact is also known at node j (i.e., a Dirichlet condition), one can specify both ) and k at this node as j = kezact,j and 4' = ezxact,j The complex potential w is thus fully specified, and no nodal equation needs to be applied. This node will be called a "completely specified node". In some circumstances, the position of the psi reference node can be important to the application of the CVBEM. To properly discuss this phenomenon, a detailed examination of the way that Neumann and Robin boundary conditions are implemented in the CVBEM is necessary. Since such an examination is also necessary when considering treatment of corners, further discussion concerning placement of the psi reference node will be delayed until Chapter IV where the treatment of corners is covered. This chapter has laid the foundation for the development of the CVBEM. A linearelement model has been presented here; a more refined model follows in the next chapter. CHAPTER III FORMULATION OF THE CVBEM USING QUADRATIC ELEMENTS As with the linearelement CVBEM, the formulation of the quadratic element CVBEM begins with Cauchy's integral formula, W(zo) = Z (C o E f zo r (3.1) 27ri r C( zo This time, however, boundary F is discretized into M finitelength curved elements, ri. These elements are formed by passing a quadratic polynomial through three successive nodal points. The domain boundary is taken as the union of these elements as shown in Fig. 31, i.e., M r=U ri (3.2) j=1 Notice that the element j contains nodes k, k +1, and k + 2, so that the indices j and k are related by k = 2j 1. Further, M elements are formed from N = 2M nodal points. The function w(z) is now approximated by a quadratic global trial function, G2(z), given by M G,(z) = .Mk(z) Wk + Mk+l(Z) Wk+1 + Mk+2(z) Wk+2 (3.3) j=1 k = 2j1 37 r Zk2 8(k+2,k2;k) Zk1 \ 1+2 zk 8(1+2,1;k) ri 8(k+2,k;0) Zk+2 F2 z 5 iy zN z3 ZN z X FIGURE 31. Boundary discretization and angle definitions in the quadraticelement CVBEM. Here, Mk(z), Mk +(z), and Mk + 2(z) are continuous basis functions representing the effects of wk, k + 1, and k + 2 over elements rj and Fj. These basis functions are taken as seconddegree Lagrange polynomials of the form ( k 2)( k 1) (Zk Zk 2)(Zk Zk  (zk+1Z) (Zk+2Z) (Zk 2)(Zk k 1) 0 z e r 1 z E r, z g i_, U rI By substituting G2(() for w(() in the right hand side of the Cauchy integral formula (Eqn. (3.1)), the secondorder CVBEM approximation of w can be expressed as (Zo) = I j zo r zo E Q zo o r As in the linearelement CVBEM, the contour integral can be eliminated (see Appendix B), and the above equation reduces to M O(z0) = 2 Ck wLk + Ck+ k+1 + Ck+2 Wk+2 j=1 k = 2j 1 (3.6a) where Ck{ (Zk +2 2k + 1)(Zk + 2 Zk + )hk D (zk + 2 + k+ )(k + 2 zk + )[(zk +2 zk) + zhk] D (zk+2 k+ ){(k+2z)+o(+2 )+ Dk ZO[(Zk + 2 Zk) + z0ohk + D I Mk(z) = (3.4) (3.5) (3.6b) Ck+ (k + 2k)(k + 2 k)hk Ck+i = D S(Zk + 2 + Zk)(Zk + 2 Zk)[(Zk + 2 Zk) + zohk] D (zo z}k) (Zk + 2 k) k + 22 + Z[(zk + 2 ) + hk D (3.6c) Ck+2 (k + lZk)(Zk + 1 zk)hk (zk+1 + Zk)(Zk +1 Zk)[(Zk+ 2 Zk) + 0hk] D (Zk +1 Zk){( k+ + zo[(zk + 2 Zk) + } (3.6d) + D2 (3.6d) D = (k +2zk)(zk + Zk)(Zk+2 Zk +) (3.6e) hk = I(k + 2 In (k + 2 + i0(k + 2, k; 0) (3.6f) (Zk ZO) (Zk 0Z) See Fig. 31 for the definition of 8(k + 2, k; 0). Eqn. (3.6) is analogous to Eqn. (2.6) in the linearelement CVBEM, and, like Eqn. (2.6), it embodies two equationsone for O(zo) and one for O(zo). Close inspection of Eqn. (3.6) reveals that it may be used to calculate the value of 0 at any location z0o zk, provided the values of w are known at all nodes. Thus, one can apply Eqn. (3.6) not only at interior points, but also at the "middle" node of any boundary element, r (zo = zk +1). The quadraticelement analog to Eqn. (2.7) which can be applied at the "end" nodes (zo = zk) is derived in Appendix B, and is given by Fk2 k2 + FkI k1 + Fkwk + Fk+l Wk+1 + Fk+2 k+2 M + C CW + Ct+Ii w+ + i=1 i,i+1/j k=2j 1 I = 2i1 where F z fk2zk_ +Zk2 F2 I 2(zk I1 Zk 2) Fk ( Zk2 ) 2(zk: )(Z1 + 2) Fk = { (Zk+2zk) + io(k +2z, k2; k) (zk 2 Zk) S3z 2zk k_ 2 2(zk zk 1) 3Zk 2k+ zk+2k 2(k k + 1) k +Zk Zk +2 2 S= 2(Zk k + )(Zk +1 Zk + 2) Fk+2 2(zk2k+z k+2 Cl = { (zl+2Z +1)(Zl +2 z2 + )h C= D D + D ( + 22+ !ZI +zk~l Cl+2 w+2} (3.7a) (3.7b) (3.7c) (3.7d) (3.7e) (3.7f) (3.7g) k = 21{ C (z (z+2I)(ZI+2 z)hI (zZ + 2 + z1)(Z1 + 2 1)[(z + 2 Z) + khl] D ((2 Z) + Zk[( 2 ZI) + Zh] 1 2 D (3.7h) C { (z1+lzt)(z+ I zl)h, C+2 D (Zl +1 + ZI)(ZI +1 I)[(z + 2 ZI) + Zkhl] ( Z ) +{.2 z [( z)) + ]} + D (3.7i) D = (zI+2zl)(ZI+Iz1)(z1+2zg+l) (3.7j) _(__ _I ) (zI+2zk) h, = In (T +2 zk) in ( Zk) +iO( + 2, ; k) (3.7k) (ZI Zk) I (Z Zk) See Fig. 31 for the definition of O(k + 2, k 2; k) and 0(1 + 2, 1; k). As with Eqn. (3.6), this complicated expression can be broken down into two equationsone for qk and one for fk. This breakdown must be performed so that the explicit, implicit, and hybrid methods of solving for the unknown values of 0 and 0 at the boundary nodes may be implemented just as in the linear element CVBEM. The desired nodal equations as derived in Appendix B are S2= F,_2 k2 2 + F2k2 + F I1 Ok1 k+ F1 k1 + Fk Ok + Fk Ok + FIk+1 k+1 + F+1 4k+1 + F +2 k+2 + F+2 Ok+2 M + E [Cf I + CfR + c+, +1 + Cf+1 +1 i= I 1,2+1 #kR (+.+) k=2j1 + C2 +2 1+ '+2 0+2 (3.8) 1= 2i1 and k Fk_ k2 Fi2 Ok2 Ok = F kr2  + F_, 1ki F_1 _k1 +F F'k f k + Fk+1 k+1 F+1 7k+1 + Fk+2 k+2 Fk+2 Ok+2 M + E [cl cf C + C+1 ,+1 cJ+iV +1 i=1 i +i+ 2 + 2 C1+ 2 1 + (3.9) k = 2j i + Of+2 91+2 1+21+2 1=2i1 Here, the superscripts R and I refer to the real and imaginary parts of the superscripted quantities, respectively. Complete expressions for these real and imaginary parts are given in Appendix B. 3.1 Boundary Conditions and Solution Methods The boundary conditions and solution methods discussed in Sections 2.1 and 2.2 in conjunction with the linearelement CVBEM are equally applicable when quadratic elements are used, with some minor modifications. The four types of boundary conditionsDirichlet, Neumann, Robin, and stream functionand their associated equations (Eqns. (2.11) through (2.14)) remain unchanged except for the replacement of nodal index j with index k. The derivation of the explicit, implicit, and hybrid method equations require only that Eqns. (3.8) and (3.9) be used in place of Eqns. (2.8) and (2.9) as detailed below. Explicit Method: For a Dirichlet condition imposed at node k, one sets Ok = Ok = qk and Ok = Ck in Eqn. (3.8) to obtain F2 k2 + F2 k2 Fi, + F R F+ ik + F R Fk+21 k+l + Fk+1 k+1 F/+2 Ok+2 + Fk+2 k+2 M iI =1 S= 2i 1 (3.10) This equation is referred to as the explicit phi nodal equation. Eqn. (3.9) is dropped. For a stream function condition imposed at node k, one sets 1,k = Ok = O and Ck = k in Eqn. (3.9) to obtain Fk2 k2 F2 k2 Fk1 kI FI _k F, k F k FR+1 +k+ F+2 k+2 F2+2j =k { k = { M + E [Cf Q Cf 4, + CR+, 1 +1 c+, +1 k=2j1 + \ +2 1+2 O(+2 01+2J (3.11) 1= 2i1 This equation is referred to as the explicit psi nodal equation. Eqn. (3.8) is dropped. For the explicit method, solution is effected by applying Eqn. (3.10) at each Dirichlet node and Eqn. (3.11) at each stream function node to form a system of simultaneous equations which is solved for the unknown values of and 4. Implicit Method: For a Dirichlet condition imposed at node k, one sets k = k and Ok = k in Eqn. (3.9) to yield k 2= F~k2k2 k2 k2 + F Ck1 F[,1  R + Fk~ Fk O, + F Oqk+1 Fk+1 k+1 + FR+i +i F+ + F k+2 k+2 Fk+2 Ok+2 M + E [c~ c + cR+, ,+1 c+, + =21 + C+2 + C+2 1+2] (3.12) S= 2i 1 This equation is referred to as the implicit psi nodal equation. Eqn. (3.8) is dropped. Similarly, for a stream function condition imposed at node k, one sets Ck = 1k and Ok = k in Eqn. (3.8) to obtain F2 k2 k 2 k2 Fk1 k1 k1 OkI Fk g + Fk 4k F+1 k+1 + F+1 k+i FI+2 k+ + 2 +2k+2 M + [Cf 0, + CR , i=1 k=2j1 I = 2i1 + cfI+++ + CRf+l t+] + Cf2 1+2 + C+2 01+21 (3.13) This equation is referred to as the implicit phi nodal equation. Eqn. (3.9) is dropped. In a manner similar to the explicit method, solution by the implicit method is effected by applying Eqn. (3.12) at each Dirichlet node and Eqn. (3.13) at each stream function node to form a system of simultaneous equations which is solved for the unknown values of q and b. Hybrid Method: For a Dirichlet condition imposed at node k, one sets Ok = k on the left side of Eqn. (3.8) and Ok = Ok and Ok = k on the right sides of Eqns. (3.8) and (3.9) to obtain Fk2I + FR k2 Ok2 k 2 Fi + F R Fk1 k1 + F1 Fkq Sk + F'k Fi+1 k+1 + Fk+l Fk+2 Ok+2 + Fk+2 Ok+1 OkC+2 1 k2n 1 S= 2 { M + M [cIf + cf + cf+,1 0+1, + cf+ 0i+1 S+= 1 + C2 k=2j1 1+2 0'1+2 1+2 1+2 (3.14) 1=2i1 I k=  Fk2 k2 F2k2 + FkR_1 k1 F1 _k1 + Fk qk Fk + F+ +lk+I Fi+I k+ + Fk+2a k+2 Fk+2 k+2 M + E [CR I CI 0 + CC + +I +1 :,+11 + (3+1) k 2j + C ,I+2 1+2 C1+2 I+21 (3.15) S=2i 1 These equations are referred to as the Dirichlet hybrid phi nodal and psi nodal equations, respectively, since Eqn. (3.14) is explicit while Eqn. (3.15) is implicit. For a stream function condition imposed at node k, one sets k = k on the left side of Eqn. (3.9) and Ok = k and Ok = "k on the right sides of Eqns. (3.8) and (3.9) to yield F2 2k2 + Fk2 k_2 F1 Ok1 + Fk1 OkI Fi k + FkR F[+k+i + F+Il+l Ff+2 1k+2 I+ F +1 2 k+2 FT k+2 + 2 Ok+2 M + o [cf, + cfR + CI/+f 1+I + cr+01 + ,+1 1 1 +1 1+1 C+11+ 1 k=211 + 2+ 21+2 + 0, +2 (3.16) S= 2 1 k = 2  k 2= F k2k2 F2 k2 + Fi1 k1 FL1 k + F k F k + Fk+1 k+1 FI+1 Ck+1 + Fk+2 k+2 Fk+2 ?k+2 M + E [Cf 1 CI, + CR+I 1+1 C+ +1 i=1 k = 1 + +2 01+2 C2 +2 (3.17 1= 2i1 These equations are referred to as the stream function hybrid phi nodal and psi nodal equations, respectively, since Eqn. (3.16) is implicit while Eqn. (3.17) is explicit. In the hybrid method, solution is effected by applying Eqns. (3.14) and (3.15) at each Dirichlet node and Eqns. (3.16) and (3.17) at each stream function node to form a system of simultaneous equations which is solved for the unknown values of j and '. For nodes where Neumann or Robin conditions are specified, Eqn. (2.12) or Eqn. (2.13) (with j replaced by k) will be used in conjunction with the implicit phi nodal equation, Eqn. (3.13), in the solution of unknowns. By applying Eqns. (2.12) and (3.13) at all Neumann nodes, Eqns. (2.13) and (3.13) at all Robin nodes, and either the explicit, implicit, or hybrid method equations at all Dirichlet and stream function nodes, a system of linear algebraic equations of the same form as that for the linearelement CVBEM is obtained. As mentioned previously, a detailed description of how the system is formed is given in Appendix A. The system can be solved by any direct method such as Gaussian elimination. Once the unknown values of q and i are found, one can use these values, together with the specified values (< and b), as the boundary nodal values needed by Eqn. (3.6) for calculating W at any interior point. 3.2 Additional Comments Here, a few remaining comments regarding the quadratic element CVBEM are made. As with the linearelement CVBEM, the implicit method is recommended over the explicit and hybrid methods for the same reasons mentioned in Section 2.2. Also, the calculation of at the boundary nodes proceeds just as described in Section 2.3. Finally, corner nodes are still a special case, as will be explained in Chapter IV. The description of the quadratic element CVBEM is now complete. In the next chapter, a methodology for the proper treatment of boundary covers in either the linearelement or quadraticelement CVBEM is presented. CHAPTER IV TREATMENT OF CORNERS IN THE CVBEM It was mentioned in Chapter I that the covers in a domain boundary can cause problems for the RVBEMs since is doublevalued at a corner node. n The CVBEM avoids this problem by solving for 0 instead of Since 0 is both continuous and singlevalued at a corner (just like 0), corner nodes can generally be handled more easily than in the RVBEMs. A different "corner" problem arises in the CVBEM, however, due to the way that Neumann and Robin boundary conditions are implemented. In fact, this problem is not limited to corer nodes but exists at nodes which lie at any interface between two different boundary conditions as well. This chapter is concerned with the treatment of corner nodes in the CVBEM, but due to the nature of the method, it is instructive to first consider how nodes which lie at boundary condition interfaces (BCI nodes) should be handled. The principles developed for dealing with BCI nodes are then used to establish a methodology for the treatment of corner nodes. 4.1 BCI Nodes Before considering how BCI nodes should be treated in the CVBEM, it is necessary to clearly define what a BCI node is. A BCI node is defined as the first node imposed with a different type of boundary condition than the node or stretch of nodes immediately preceding it. Here, the direction of travel is assumed to be the same as that of the contour integral in Eqn. (2.2). Each time the boundary condition changes, a BCI node is encountered. Thus, in the sample domain of Fig. 41, there are four BCI nodes along the boundary. (The letters D, N, and R are used to represent Dirichlet, Neumann, and Robin conditions, respectively, imposed at the nodes.) In this way, the "new" boundary condition is assumed to begin along the boundary just after the node preceding the BCI node. Thus, only one condition is actually imposed at the BCI node. A "doublyspecified" BCI node where two boundary conditions are available for specification is a different case and will be discussed in the next section along with corner nodes. Proper treatment of BCI nodes in the CVBEM requires examination of the equations for implementing the boundary conditions. Dirichlet BCI nodes pose no problem whatsoever as the potential is directly specified via Eqn. (2.11). In contrast, the equations for implementing Neumann or Robin boundary conditions (Eqns. (2.12) and (2.13)) must be applied with caution since they are "upstream" equations, meaning that they require information from the previous "upstream" node. For example, suppose that node m is the BCI node which starts a Neumann or Robin stretch of boundary. To properly apply Eqn. (2.12) or (2.13), the values of n or [H( O 0)] must be known at node m 1. If node m 1 is a Neumann or Robin node, then these values will be known, and BCI node m can be treated like any other Neumann or Robin node; however, if node m 1 is a Dirichlet node, then these values will not be known, and neither Eqn. (2.12) nor Eqn. (2.13) can be applied directly at node m. One way to remedy this deficiency is to treat node m as the psi reference node mentioned in Section 2.4. Equation (2.12) or (2.13) would then be applied at all of the "downstream" nodes along the remainder of the Neumann or Robin stretch of boundary (nodes m + 1, m + 2, etc.). Unfortunately, only one psi reference node is allowed along the boundary since the constant of integration can only be specified once. Consider the case shown in Fig. 42. Node 3 is the FIGURE 41. Example showing the location of the BCI nodes. n is not known here D N 6 N N 3 10 11 second Neumann BCI node S2 \ psi reference node D 2 (first Neumann BCI node) FIGURE 42. Example showing a section of boundary containing two Neumann BCI nodes. BCI node which starts the first Neumann stretch of boundary, so it becomes the psi reference node. Node 8, which is the BCI node starting the second Neumann stretch of boundary, cannot then also be made the psi reference node because the reference value for 0 has already been specified at node 3. Neither can the Neumann equation for psi (Eqn. (2.12)) be applied at node 8 since at node 7 is not known. There are two possible remedies for the situation: (1) both and 4 must ,On be specified at node 7 or (2) a must be estimated at node 7 using the Cauchy Riemann conditions and a finitedifference expression, e.g., 5: = 07 0z6 This leads to the following altered version of Eqn. (2.12) which can be applied in lieu of that equation at node 8: = =7 + + } 0 '6 1Z8 Z71 Fn Z7 Z6l In a more general form, for node j, the above equation becomes = ik 1 + 2{()i 2 Other more accurate finitedifference expressions can certainly be used to approximate ( 1 when possible. If node 8 had started a Robin stretch of boundary, a similar situation would occur, and the two possible remedies would be that (1) both Ho, and ( must be specified at node 7 or (2) H(4 0) must be estimated at node 7 using the Robin equation, the CauchyRiemann conditions, and a finitedifference expression, e.g., = I7 6 This leads to the following altered version of Eqn. (2.13) which can be applied in lieu of that equation at node 8: 08 = t + { H( 0)s + it7 6 }Z8 71 2 z IZ7 Z61 In a more general form, for node j, the above equation becomes 01 = 0C, + {H(1)j+ i  2 j 2 1 Iz1Zj21 Once again, other more accurate finitedifference expressions can be used to approximate H(0 ~oo) 1 where appropriate. There is yet another reason why a BCI node which begins a Neumann or Robin stretch of boundary should be treated as the psi reference node when possible. Consider a case where a Neumann or Robin condition exists along some part of the boundary, say from node j to node j + k; see Fig. 43. Equations (2.12) and (2.13) show that, in either case, the value of 0 at nodes j + 1 through j + k depend intimately upon the value of 0,. In fact, any error in 0, will be carried over as errors in 0i + through 0j+k. The errors in 4 along the Neumann or Robin stretch of boundary will, in turn, cause errors in the estimated values of 4 and 0 at all of other nodes along the boundary due to the coupled nature of the equations in matrix equation (2.23). If the Neumann or Robin part of the boundary comprises a significant portion of the total boundary length, this can result in a large amount of error. Thus, it is imperative that Neumann or Robin condition imposed here / j+k BCI node FIGURE 43. A Neumann or Robin stretch of boundary. the value of 0 at the first node of a Neumann or Robin stretch of boundary be as accurate as possible. With this in mind, the following guidelines are provided for the choice of the psi reference nodal location: 1) If portions of the boundary (or the entire boundary) are imposed with Neumann or Robin conditions, let the psi reference node be the BCI node of the longest Neumann or Robin stretch following the direction of contour integration of Eqn. (2.2). 2) If the entire boundary is imposed with a Dirichlet condition, no BCI nodes exist, and the optimum location of the psi reference node can be determined by trial and error. Experience has shown that its placement in this "all Dirichlet" case makes little difference in the final overall solution. The position is not crucial to accuracy. To summarize, Dirichlet BCI nodes require no special treatment and can be handled just like any other Dirichlet node. On the other hand, Neumann and Robin BCI nodes should be treated using the remedies and guidelines detailed above, subject to the limitation that only one psi reference node is allowed. 4.2 Corner Nodes And DoublySpecified BCI Nodes When using the CVBEM to solve potential problems with corners in the boundary, the same nine traditional boundarycondition combinations in Table 11 are possible. As in Chapter I, the optimal case where both the upstream and downstream conditions are known at the corner is considered. Nonoptimalcase corner nodes (i.e., corner nodes which are imposed with only one boundary condition) can be treated in the same manner as the BCI nodes described in the preceding section; however, corner nodes where two boundary conditions are imposed simultaneously can be treated specially in the CVBEM. For such "doublyspecified" nodes, the unknowns in Cases 1 through 9 are listed in Table 41. Since the CVBEM corerr" problem is due to the way that Neumann and Robin boundary conditions are implemented and is not due to the doublevalued nature of the normal gradient of the potential (as in the RVBEM), "doubly specified" BCI nodes (i.e., BCI nodes where two boundary conditions are available for specification) may also be treated as "corner nodes". From this point on, no distinction will be made between "doublyspecified" corer nodes and "doublyspecified" BCI nodes, and the term corer node will be understood to encompass both. Earlier in Section 2.1, four boundary conditions were discussed, namely, Dirichlet, Neumann, Robin, and stream function. A nomenclature is now introduced whereby a particular node is identified by the boundary condition imposed there. Thus, four types of nodes are identifiedDirichlet nodes, Neumann nodes, Robin nodes, and stream function nodes. A fifth type of node the completely specified nodewas discussed previously in Section 2.3. Now, a question arises as to whether in each of the nine different corer cases, the corner node can be treated as one of the five nodal types above. In order to address this question, the equations associated with the different boundary conditions must again be examined as they were for BCI nodes. The application of these equations at arbitrary corner node m will now be investigated. As mentioned in the Section 4.1, the Neumann and Robin equations (Eqns. (2.12) and (2.13)) are "upstream" equations for calculating 0,, meaning that they require information from the previous upstream node, m 1. Due to this "upstream" nature, applying either Eqn. (2.12) or (2.13) at a "doublyspecified" corer node which is the first node of a Neumann or Robin stretch of boundary is TABLE 41. The unknowns and corresponding available in the implicit CVBEM "doublyspecified" corner node combinations listed in Table 11. number of equations for the nine possible boundary condition Unknowns Number of ns. Available 1 4 2 2 4 2 3 q and 3 4 4 2 5 4 2 6 and 3 7 and 3 8 and 4 3 9 4 1 Case undesirable since it forces the downstream Neumann or Robin condition back to the previous node m 1 where the condition is not imposed (see Fig. 44). Notice that this situation is different from that described in the previous section for BCI nodes. The BCI node boundary condition was actually assumed to begin just after the node immediately preceding the BCI node. Here, for the "doubly specified" case, the upstream and downstream boundary conditions are assumed to terminate and begin, respectively, at the corer node itself (see Fig. 45). Thus, in Cases 1, 3, 4, 6, 7, and 8, although two conditions are available, the corner node should be treated using the upstream condition rather than the downstream Neumann or Robin condition. In contrast to the Neumann and Robin equations, the equation for imposing a Dirichlet boundary condition (Eqn. (2.11)) involves direct specification of 0,m. Therefore, Case 9 is easily handled as a Dirichlet node. This leaves Cases 2 and 5 which are grouped together as an upstream equation for b,, followed by a direct specification of 0k,. For these cases, ,,, can be specified directly using Eqn (2.11), and the upstream Neumann or Robin equation (Eqn (2.12) or (2.13)) can then be applied at node m. Clearly, this is not one of the five nodal types described above. This new type of node will be referred to as a NeumannRobin/Dirichlet node. To summarize, it has been shown that, of the nine corner cases, seven can be handled using previously defined nodal types. Only twoCases 2 and 5 require the addition of a new type of node. 4.3 Rules for the Nine Corer Cases A list is now presented which gives the proper treatment for each of the nine cases in the implicit CVBEM at arbitrary corner node m. The equation numbers that are given in parentheses refer to the appropriate nodal equations downstream Neumann or Robin condition . The downstream Neumann or Robin condition will be forced back to this node if Eqn.(2.12) or (2.13) is applied at node m. upstream boundary condition Diagram showing the effect of applying the downstream Neumann or Robin condition at corer node m. FIGURE 44. Neumann condition begins along boundary just after N this node N 0 BCI node D Dirichlet condition terminates at this node. Neumann condition begins at this node. J N D FIGURE 45. Diagrams showing the location where a boundary condition change is assumed to take place for a BCI node and a "doublyspecified" corner node. for both the linear and the quadraticelement formulations of the CVBEM. Case 1: Treat the corner node as a Dirichlet node; 1m is then the only unknown. Apply the implicit nodal equation for 0., (Eqn. (2.17) or Eqn. (3.12)). Case 2: Treat the corner node as a NeumannRobin/Dirichlet node; 1m is then the only unknown. Apply the Neumann equation for 1m (Eqn. (2.12)). Case 3: Treat the corner node as an upstream Neumann node. Both , and 1,k are unknown. Apply the implicit nodal equation for m, (Eqn. (2.18) or Eqn. (3.13)) and the upstream Neumann equation for km (Eqn. (2.12)). Case 4: Treat the corner node as a Dirichlet node; 1.m is then the only unknown. Apply the implicit nodal equation for 1,m (Eqn. (2.17) or Eqn. (3.12)). Case 5: Treat the corner node as a NeumannRobin/Dirichlet node; 1, is then the only unknown. Apply the Robin equation for 1, (Eqn. (2.13)). Case 6: Treat the corner node as an upstream Robin node. Both .m and 1m, are unknown. Apply the implicit nodal equation for 0m (Eqn. (2.18) or Eqn. (3.13)) and the upstream Robin equation for 1m (Eqn. (2.13)). Cae : Treat the corner node as an upstream Neumann node. Both , and 1,m are unknown. Apply the implicit nodal equation for qm (Eqn. (2.18) or Eqn. (3.13)) and the upstream Neumann equation for 1km (Eqn. (2.12)). Case a: Treat the corner node as an upstream Robin node. Both 0, and ,m are unknown. Apply the implicit nodal equation for 0m (Eqn. (2.18) or Eqn. (3.13)) and the upstream Robin equation for Cm (Eqn. (2.13)). Case 9: Treat the corner node as a Dirichlet node; 1,m is the only unknown. Apply the implicit nodal equation for 1,, (Eqn. (2.17) or Eqn. (3.12)). An exception to the above nine rules involves the psi reference node described in Section 2.3. This exception ties in with the comments made in Section 4.1 for BCI nodes and is listed below. Exception: The first nodal point of the longest stretch of boundary imposed with a Neumann or Robin boundary condition (traveling from the initial point along the boundary, keeping the interior of the domain on the left) should be treated as the psi reference node described in Chapter II. This important exception can supersede the rules listed for Cases 1, 3, 4, 6, 7, and 8. This completes the discussion of the methodology for handling corners and boundary condition interfaces in the CVBEM. Examples of solutions obtained using the corner methodology will be given in Section 6.1 of Chapter VI. In the 64 chapter that follows, a novel application of the CVBEM in the area of numerical grid generation is presented. CHAPTER V NUMERICAL GRID GENERATION USING THE CVBEM As mentioned in Chapter I, grid generation methods can be divided into two main classesalgebraic methods and differential equation methods. The algebraic methods, being much faster, tend to be the methods of choice for unsteady, deformingboundary problems where a new grid is required at each time step. Unfortunately, unlike the differential equation methods, the algebraic methods often do not produce the smooth, nonoverlapping grid lines essential to obtaining accurate FDM solutions. In this chapter, a variation of the linearelement CVBEM which can be used to generate grids for 2d, simplyconnected spatial domains is presented. Since the solution of Laplace's equation is involved, this method can be classified as a differential equation method; however, it is anticipated that the method could prove to be computationally efficient enough to bridge the gap between the existing algebraic and differential equation methods. 5.1 Derivation of the Inverted CVBEM Equations The linear element CVBEM interior point equations (represented collectively by Eqn. (2.6)) will first be rearranged so that if one knows the values of and ? at an interior point, one can solve for the point's location, Zo = x0 + iyo, provided that the values of ( and 4 at the boundary nodes are known. These "inverted" equations will form an integral part of the CVBEM grid generation method. To begin, Eqn. (2.6) is split into its real and imaginary parts as N1 ,(zo)= _ + 2 E 3=1 j+I[(xoxj)D + j[(xo zxj + )D + +,[(xox1j+l)C  CJ+I[(xoo XJ)C  Oj + I[(xo zj)D  Oj[(z o Xj+ )C + Oj[(xo xj + )D (yo Y)C] (yo y,)DI (Yo yj + )C] (yoy +i+)D]  (yo yj)D] + (yo y)C]  (yoyj+,)D] + (yo yji+)C] C = [A(x+1x,) + B(yj+,yUj)] / F D = [B(xz+lI.) A(yj+ly)] / F F = (xj+lxj)2 + (yj+iy))2 A = in (zi zo) B = (j + 1,j; O) Eqns. (5.1) and (5.2) can be rewritten as N 27r(zo) = j=1 j[zxoC xj+,C yoD + yji+D] Oj[xoD x.+,D + yoC yi+C] +,j+l[xoC xC yoD + yD] + j+1[xoD xjD + yoC yC]) and (5.1) where (5.2) (5.3a) (5.3b) (5.3c) (5.3d) (5.3e) and 27r(zo)= E{ j[xoD xj+D + yoC yj+iC] j=1 j[xoC xj+iC yoD + yj+iD] j+ xzoD zxD + yoC yC] + j+l[zoC xjC yoD + yD] } which can be further rearranged as 27r(zo) = j=1 and  2r(zo) = N j=1 xo[C( J+1 Oi) 4 yo[D(O i+ 1) 4 [ j(X +,C + Oj(xj+,D  +oj+1(xjC + + j+1(zjD  xo[D(j j + ) + yo[C(4j ,j +1) + [ Oj(x5+,D j(xi+iC j+i(xjD +,( zXD + j+ 1( xC  D(j+ 1 j)] C(j +I 0] yi+ D) y,D) yj + C) YjD) y,C)] } + C(Oj +1 j)] + D(j +,1)]  y,+iC) + Yj+ID)  yC) + yjD)] Terms in the braces can be expressed compactly as Cl, = C(Oj,+1O) + D(j+I,) C2,j = D(Oj j+1) + C(Oj+,1j) C3,, = j(x+iC + yj+D) Oj(xj+lD yj+C) + i+1(xC + yjD) + 4j+l(xiD yC) (5.4) (5.5) (5.6a) (5.6b) (5.6c) C4,j = Oj(xj+lD yj+lC) Oj(xj+lC + yj+,D) Oj+l(xjD yC) + j+l1(zjC + yD) (5.6d) Substituting Eqns. (5.6a) through (5.6d) into Eqns. (5.4) and (5.5) yields N N N xoE C,, + yoEC2,, = 27r(zo) EC3,, (5.7) j=1 j=1 j=1 and N N N XEC2,j yo3CI,j = 2r(zo) EC4, (5.8) j=1 j=1 j=1 With the domain geometry specified and the values of q and b known at interior point zo (the hats can then be dropped) and at the boundary nodes, Eqns. (5.7) and (5.8) become two simultaneous equations in two unknownsao and yo. Unfortunately, these equations are nonlinear due to the presence of terms A and B in the summations. Recall that A and B contain zo = xo + iyo (see Eqns. (5.3d) and (5.3e)). Since Eqns. (5.7) and (5.8) are nonlinear, a singlepoint iterative solution process will be used to solve the equations simultaneously for xo and yo. To facilitate the implementation of this process, Eqns (5.7) and (5.8) are combined to yield N N 27ro(zo) E~C3, 2 (zo) E C4,j j=1 j=1 N + N clC, I1C J zi Yo = __ __ (5.9) SN Cl, E C2,j j=1 j and N N 2r0(zo) EC,,j EC2,, 0 = (5.10) EC xj EC1,,j j=1 j=1 The process of solving for ox and Yo proceeds as follows: 1) Guess initial values for xo and Yo. 2) Calculate Cj,, C2,j, C3,j, and C4,j using the values of xo and yo. 3) Calculate a new value of yo using Eqn. (5.9). 4) Calculate a new value of zx using Eqn. (5.10); the values of C1,j, C2,j, C3,j, and C4, are taken from Step 2, while the new value of yo from Step 3 is used. 5) Repeat Steps 2 through 5 until the percent difference between the new and previous values of xz and Yo is less than a predetermined constant. These five steps comprise the desired iterative process for calculating the location of an interior point when the values of 4 and 0 are known at the boundary nodes and at the point itself. 5.2 The CVBEM Numerical Grid Generation Algorithm The iterative process described above is now combined with the linear element CVBEM to create a method for generating grid systems in 2d, simply connected spatial domains. The method will henceforth be referred to as the complex variable boundary element grid generation method or CVBEGGM for short. The fundamental concept behind the CVBEGGM is the use of potential lines and streamlines as grid lines, and it is intended for domains whose boundaries can be divided into four separate, continuous, smooth or nonsmooth curves. Examples of such domains are given in Fig. 51. The CVBEGGM consists of the following eight steps: Step 1 Specify the Domain Geometry. Prescribe the locations of the nodal points along the domain boundary. The number of nodal points will have a direct effect on the accuracy of the solution (i.e., the quality of the grid). Generally, as more nodal points are used, the grid quality improves, but the amount of time required to generate the grid increases as well. There is therefore a compromise to be made, and a rule of thumb is that one should use as many nodal points as possible for the computational time constraints imposed. Certainly, the number of nodal points should be great enough to adequately describe the boundary contour itself but need not be equal to the number of grid points ultimately desired along the boundaries. Step 2 Define the Mapping. The CVBEGGM is intended to map a 2d, simplyconnected, physical spatial domain onto a 2d, rectangular transformed domain. The boundary of the physical domain will be mapped onetoone to the four sides of the transformed domain, forming a socalled bodyfitted coordinate system [46]. Since the grid lines of opposite families connect opposing sides of the rectangular transformed domain (see Fig. 52b), the boundary of the physical domain needs to be broken down into four separate boundary curves (see Fig. 5 2a). Generally, there will be a "natural" way to perform this breakup, but consideration should be given to the type of grid desired. Reference 46 gives a good discussion of various breakup considerations for simplyconnected domains. Step 3 Apply the Boundary Conditions. The boundary conditions which should be imposed to insure that the grid lines intersect the boundaries orthogonally are shown in Fig. 53. These boundary conditions force boundary FIGURE 51. Simplyconnected domains whose boundaries can be broken up into four separate boundary curves. boundary curve 3 curve 2 boundary curve 1 a) The physical spatial domain. lines of constant 3 grid lines of constant 4 5 b) The transformed domain. FIGURE 52. The physical and transformed domains. boundary curve 4 0.01 0.0 ao = 0 a 0 On~ (boundary curve 3) (b=1 (boundary curve 2) (boundary curve 4 ) y S=o0 (boundary curve 1) FIGURE 53. Boundary conditions for the CVBEGGM. curves 1 and 3 to be streamlines and boundary curves 2 and 4 to be potential lines. Streamlines in the domain will thus connect and intersect orthogonally boundary curves 2 and 4, while potential lines in the domain will connect and intersect orthogonally boundary curves 1 and 3. The streamlines and potential lines will become grid lines and will be mapped to lines of constant ( and ry, respectively, in the transformed domain. These grid lines will intersect both each other and the domain boundaries orthogonallyexcellent characteristics for a grid system to possess [46]. Corner nodes should be treated following the rules set forth in Chapter IV, which are given as follows: a) the corner node between boundary curves 4 and 1 completely specified node (0 = 0, 0 = 0), b) the corer node between boundary curves 1 and 2 completely specified node (4 = 1, ip = 0), c) the corner node between boundary curves 2 and 3 Dirichlet node (4 = 1), and d) the corner node between curves 3 and 4 is treated Robin/Dirichlet node. is treated as a is treated as a is treated as a as a Neumann Step 4 Solve for the Unknown Values of 4 and 0 at the Boundary Nodes. This is accomplished using the linear element CVBEM as described in Chapter II. The implicit formulation is recommended. Step 5 Calculate the Boundary Grid Point Locations. There are a number of ways of specifying the boundary grid point locations, but generally, one first prescribes equally spaced grid points along boundary curves 1 and 4. The locations of these grid points are calculated by approximating the arc length along the boundary as the cumulative sum of the lengths of the linear segments connecting the nodal points. Linear spline interpolation is then used to determine the equallyspaced positions. In order to facilitate the presentation of the grid generation process, some new notation is in order; see Figs. 54 and 55. Previously in the linearelement CVBEM, the nodal points were numbered from 1 to N, proceeding around the boundary of the domain from an arbitrary starting position in the direction of the contour integral of Eqn. (2.2). Now, the nodal points associated with each of the four boundary curves will be numbered from 1 to N,, with N, representing the number of nodal points along boundary curve m (m being either 1, 2, 3 or 4). The location of nodal point k along boundary curve m will be identified as (Xm,k,Ym,k); see the example in Fig. 54. Values of the potential and stream function at this nodal point will be designated as ,m,k and 0m,k, respectively. Approximate arc lengths along the boundary at nodal point k along boundary curve m will be represented by Sm,k; its reference point for zero length will be (X,,1,Ym,,). Note that k takes on values between 1 and Nm on each boundary curve and that each corner point is associated with two boundary curves. In addition to the nodal points, there will be grid points located along each of the boundary curves. In general, these grid points will not coincide with the nodal points except at the corners and will always be referenced separately. The number of grid points along streamlines will be designated as IL, while the number of grid points along potential lines will be designated as JL. Grid point locations will be given by ((i, j), y(i,j)) where i runs from 1 to IL and j runs from 1 to JL; see the example in the physical domain shown in Fig. 55. Values of the potential and the stream function at such grid points will be designated as 0(i, j) and I(i, j), respectively. Finally, approximate arc lengths along the boundary at the grid points will be identified as s(i,j) where i and j will correspond to the i and j indices of the point's ordered pair designation, N2 2 N4 (X1,3' 1,3) 2 1 2 N S FIGURE 54. The notation used for the nodal points in the CVBEGGM. (x(8,JL),y(8,JL) I. IL 2~~Y ~ FIGURE 55. The notation used for the grid points in the CVBEGGM. ((i, j), y(i, j)); the reference point for zero arc length along each boundary curve will be given later. Using the variables defined above, the locations of the grid points along boundary curve 1 are given by x(i, 1) = s(i, ) (Xi,+X) + XI, (5.11) 1(1x, + (5.11) Sl,t+l~3lIl and y(i, 1) = S') S (Y1,+I Y1,) + Y1,1 (5.12) SI1+1 S ,I +) where i = 1, 2, ..., IL. Similarly, the locations of the grid points along boundary curve 4 are given by x(1, j) = S(1 S4" (X4,n+i X4,) + X4,,. (5.13) and y(l, j) = S, S4, (Y4, + Y4,,) + Y,, (5.14) S4,nn+1 S4, ,n where j = 1, 2, ..., JL. Here, indices 1 and n identify the nodes along boundary curves 1 and 4 that immediately precede boundary grid points (x(i, 1), y(i, 1)) and ((1, j), y(l,j)), respectively. The values of I and n are found via a simple ordered table search of the approximate arc lengths along the boundary curves at the grid points and nodal points. The approximate arc lengths for the nodal points are calculated by setting Sm,1 = 0 (5.15a) Then k Sm,k= E (XmiXmi) + (m,.i_Y,,.i_)2 (5.15b) i=2 where m = 1, 2, 3, 4, and k = 2, 3, ..., Nm. Similarly, by setting s(1,1) =0 (5.16a) the approximate arc lengths for the remaining grid points along boundary curve 1 can be computed using s(i,l) = (x(k, 1) x(k 1,1))2 + (y(k, 1) y(k 1,1))2 (5.16b) k=2 and s(i, JL) = { (x(k,JL) x(k 1,JL))2 k=2 + (y(k,JL) y(k 1,JL))2 }05 (5.16c) where i = 2, 3, ..., IL. Further, the approximate arc lengths at the remaining grid points along boundary curve 4 can be computed using s(1,j) = E (x(1,k)x(1,k1))2 + (y(1,k) y(1,k1))2 (5.17b) k=2 and s(IL,j) = ((IL,k) x(IL,k 1))2 + (y(IL,k) y(IL,k 1))2 }) (5.17c) where j = 2, 3, ..., JL. Once the grid point locations along boundary curves 1 and 4 are determined, the values of 4 and i at the grid points on these boundary curves can be calculated. It was found that linear interpolation between the values of and 4 at the boundary nodes produced satisfactory results for 4 and 0, provided that Eqns. (2.8) and (2.9) were used to calculate the boundary node values. The equations used for interpolation are given below. (i, 1) = 1) S (,+ 1,) + 1,, (5.18a) si, 1+ i i, l and s(i, 1) S1t  (i, 1) = SI (1, + ,) + 1, (5.18b) where i = 1, 2, ..., IL along boundary curve 1, and s(1, j) S4, 0(1, ) = (4",n+x 4,n) + 04,n (5.19a) S4,n+1 4,n and s(1, j) S4, n 0(1, j) = , S4 (4, +104,.) + 04,n (5.19b) where j = 1, 2, ..., JL along boundary curve 4. The values of I and n have been defined earlier, and the approximate arc lengths are the same as those calculated previously. Finally, the values of 0 and 0 at the grid points along boundary curves 1 and 4 are used to calculate the grid point locations along boundary curves 2 and 3. Because of the boundary conditions imposed, each grid point on boundary curve 1 has a corresponding grid point on boundary curve 3 which has the same value of 4. Likewise, each grid point on boundary curve 4 has a corresponding grid point on boundary curve 2 which has the same value of 0. For each grid point on boundary curves 1 and 4, the corresponding grid point on boundary curves 2 and 3 is located by using linear interpolation between the values at the nodal points. The equations used to perform this operation are given below. 0(i, 1) 3,, x(i, JL) = (X3,t+I X3,1) + X3,1 (5.20a) and y(i, JL) = ( 1 (Y3, I+Y3,,) + Y3,I (5.20b) 03, +1 3,1 where i = 1, 2, ..., IL along boundary curve 3, and x(IL, j) = (4," (X2,n1 X2,) + X2,, (5.21a) 02,n+1 2,n and y(IL, j) = 4 (Y2, n+ Y2, ) + Y2,n (5.21b) 02,n+1 2,n where j = 1, 2, ..., JL along boundary curve 2. This time, the values of I and n are found via a simple ordered table search of the values of 4 and 0 at the grid points and nodal points along boundary curves 3 and 4, respectively. Once again, Eqns. (2.8) and (2.9) should be used to calculate the boundary node values of < and 5. Step 6 Calculate the Initial Interior Grid Point Locations. There are several ways which an initial distribution of interior grid points can be specified; here it is recommended that transfinite bilinear interpolation be used [45]. This method is quite computationally efficient and has been found to provide acceptable initial guesses for the interior grid point locations. The equations for calculating the initial interior grid point locations using this method are x(i,j) = (1 i ) x(i, 1) + 77 x(i, JL) + (1 C) x(1, j) + C x(IL, j) [ (l )(l 7) x(1, 1) + ( (1 ) x(IL, 1) (1 ) (1, JL) + 77 x(IL, JL) and y(i,j) = (1 9) y(i, 1) + 7 y(i, JL) + (1 ) y(1, j) + C y(IL, j) [ (1 )(1 7) y(l, 1) + C (1 q) y(IL, 1) (1 ) y(1, JL) + I y(IL, JL) where i1 j1 = = 1 JL1 and i = 2, 3, ..., IL 1 j = 2, 3, ... JL 1 Step 7 Calculate the Final Interior Grid Point Locations. This is accomplished by using the iterative procedure described in Section 5.1. Each interior grid point location is calculated using this procedure. Step 8 Calculate the Grid Point Locations in the Transformed Domain. Having calculated the locations of the grid points in the physical domain, it is now necessary to calculate their mapped locations in the transformed domain. This operation could have been performed earlier, but it is sufficient to do so now. The locations of the grid points in the transformed domain are given by the ordered pairs (&,, j) where & = (i1) A i = 1,2,...,IL j. = (jl) A j = 1,2,..., JL and 6 = 1 J1 = IL1 J 5.3 Additional Comments On The CVBEGGM This section is included in order to mention some of the more subtle points associated with the CVBEGGM. It was stated in Step 1 of the method that the number of nodal points used can affect the quality of the resulting grid. If the boundaries are composed mostly of flat segments, then very few nodal points may be required to describe the domain geometry. This does not mean that few nodal points are required to produce an acceptable grid, however, since the values of 4 and 0 may deviate substantially from linearity over a geometrically linear boundary. Since the CVBEGGM is a numerical process, some irregularities in smoothness and orthogonality are to be expected no matter how many nodal points are used. (One might view the CVBEGGM as an approximate numerical conformal mapping.) If no errors were present along the boundary, then only roundoff error would occur, but generally there will be errors at the boundary due to the linear approximation. Increasing the number of nodal points usually increases the accuracy of the solution, thus improving smoothness and nearorthogonality. One word of caution, however. Adding nodal points also increases the number of simultaneous equations which must be solved. The recommended direct method of solution (Gaussian elimination with scaled partial pivoting) will be subject to roundoff errors when the number of equations (nodes) becomes "large." ("Large" will be machine and language dependent.) Thus, some care should be exercised when increasing the number of nodal points to insure that the solution is actually improved. Another point worthy of note involves the distribution of grid points in the physical domain. Traditional elliptic PDE grid generation methods employ Poisson's equation, utilizing the source term to control the distribution of grid lines and grid points [46]. Since the CVBEM is capable of solving Poisson's equation in certain cases, the use of this approach in the CVBEGGM is currently being investigated; however, some measure of grid line control can be achieved by simply varying the boundary grid point distribution. Instead of prescribing equallyspaced grid points along boundaries 1 and 4 (as described in Step 5), one can vary the spacing using onedimensional stretching functions. A list of such functions is given in Reference 43. Unfortunately, because of the smoothing characteristics of the Laplacian, grid lines tend to be equally spaced away from the boundaries regardless of the stretching function used [46]. Still, variation of 84 the boundary point distribution does provide a simple albeit limited way of controlling the location of the grid points in the physical domain. This concludes the description of the CVBEGGM. Examples of grids generated using this method will be included in the next chapter, along with examples demonstrating the treatment of covers and the use of quadratic elements in the CVBEM. CHAPTER VI EXAMPLES AND RESULTS In this chapter, the techniques and methods developed in the preceding three chapters are verified and tested. First, the quadratic element formulation of the CVBEM derived in Chapter III is used to solve two example problems, and the results are compared with those obtained using linear elements. Then, the corer handling methodology developed in Chapter IV is applied to several problems with varied geometries. The results obtained with the CVBEM are compared to available analytical solutions and to those of other investigators using RVBEMs where appropriate. Finally, the CVBEGGM of Chapter V is used to generate grids for four simplyconnected domains with irregular boundaries. 6.1 Quadratic Element Examples and Results The quadratic element formulation of the CVBEM was tested by application to the circular domain shown in Figure 61. Two separate complex potential distributions were assumed, namely w = z2 (6.1a) which results in X = y2 (6.1b) = 2xy (6.1c) y 9 3.0 c b\ 2.0 17 1 32 d e 1.0 25 0.0 1 I x 0.0 1.0 2.0 3.0 FIGURE 61. The circular domain used to test the quadraticelement CVBEM. a0 = 2[a cos(g) y sin(9)] (6.1d) and w = e= (6.2a) which results in S= exp(x) cos(y) (6.2b) 4 = exp(x) sin(y) (6.2c) S= exp(x)[cos(y) cos(0) sin(y) sin(g)] (6.2d) For both distributions, Dirichlet boundary conditions were specified at all 32 boundary nodes (including node 1, which was treated as a completely specified node with 4 taken to be zero). The stream function and normal gradient of the potential were then calculated at each node using the quadraticelement CVBEM. Also, 4 and 4 were calculated at five internal points labeled "a" through "e" in Fig. 61. The distribution given by Eqn. (6.1) was chosen specifically for its quadratic 4 and 0. Therefore, the solutions obtained using quadratic elements should be exact for them (except for round off errors). This was indeed the case, with the values of 4 and 4 at both the boundary nodes and interior points being 04 accurate to at least fourteen decimal places. The results for at the boundary, while not exact, were also good, with a maximum percent error of 0.59 at nodal point 17 and a mean error of 0.18 percent over all of the boundary nodes. The exponential distribution given by Eqn. (6.2) was chosen to provide for a comparison between the solutions obtained using the linearelement and quadraticelement formulations of the CVBEM. The results for the percent error in at the boundary nodes are given in Fig. 62. Here, it is noted that the percent error in at the boundary is a good indication of the overall error in On the computation for two reasons. First, in the CVBEM, must be calculated i04 from the nodal values of 4 using finitedifference formulae. This makes  generally less accurate than either 4 or 4'. Second, quantities at the boundary which are not specified by the boundary conditions tend to be less accurate than those calculated at interior points. In fact, in the solution of potential problems by the CVBEM, the maximum error in the estimated values of 0 and 4I occur on the boundary [7]. As shown in Fig. 62, the errors at each of the nodal points were reduced significantly by the use of quadratic elements, with the maximum percent error in being 13.74 for linear elements and 3.79 for quadratic elements, an almost threefold improvement. The results at the interior points are given in Table 61, and the same improvement in the accuracy due to the quadratic elements is apparent. The computational times for the linear and quadraticelement solutions referred to in Fig. 62 were 19 and 35 seconds, respectively, on an 80286/80287 based personal computer. However, it was found that a quadraticelement solution using only 16 equallyspaced nodes still possessed a smaller maximum error in than the 32node linearelement solution. The runtime for this 16 On node quadraticelement solution was 11 seconds. Additional examples of solutions generated using the quadraticelement CVBEM are presented in the next section for domains with covers in their boundaries. 9 11 13 15 17 19 21 23 25 27 29 Nodal Point Number + Linear Elements 6 Quadratic Elements FIGURE 62. Percent error in the normal gradient of the potential at the boundary nodes for the CVBEM and RVBEM solutions on the domain shown in Fig. 61 imposed with the complex potential given by Eqn. (6.2). TABLE 61. The percent error in 4 and 0b at the interior points for the CVBEM solution on the domain of Fig. 61 imposed with the complex potential of Eqn. (6.2). Linear Elements % Error in % % Error in i 0.00 0.98 0.18 1.23 0.18 2.11 3.46 1.45 3.48 0.56 Quadratic Elements % Error in % Error in 0.00 0.142 0.002 0.123 0.005 0.354 0.054 0.216 0.044 0.075 Interior Point a b c d e 6.2 Corner Treatment Examples and Results In order to assess the accuracy of the CVBEM in solving problems with covers in the domain boundary, several examples were investigated. The ones shown in Figs. 63 and 64 were taken from the paper by Walker and Fenner [32] referred to in Chapter II. For Corner Example 1 (Fig. 63), the potential and stream function distributions over the domain were taken to be S= sin(x)cosh(y) and S= cos(x)sinh(y) while for Comer Example 2 (Fig. 64), the somewhat more complicated expressions S= ln{[(x2)2 + (y2)2]0.5} + x2y2, and tan ( I 1 + 2xy (x 2) were assumed. In both examples, Dirichlet boundary conditions were imposed along all parts of the boundary (except at one node, where both 0 and 0 were specified). Sixteen nodes were used in Corner Example 1, while 32 nodes were used in Comer Example 2. Corner Examples 3 and 4 involve square domains, and they differ only in the boundary conditions imposed as shown in Figs. 65 and 66. Here, the potential and stream function distributions over the domains were assumed to be 0 = exp(z) cos(y) 13 14 15 16 Y 0.5 0.25. 0.0 0.0 135 0.25 FIGURE 63. The geometry and nodal discretization for Corner Example 1. 0 1I 0.0 1.0 FIGURE 64. The geometry and nodal discretization for Corer Example 2. 0.51 0.0. 0.5 1 .0 1.0 ' 