Extensions and refinements to the complex variable boundary element method including its application to numerical grid g...

MISSING IMAGE

Material Information

Title:
Extensions and refinements to the complex variable boundary element method including its application to numerical grid generation
Physical Description:
vii, 151 leaves : ill. ; 28 cm.
Language:
English
Creator:
Bailey, Robert Timothy, 1963-
Publication Date:

Subjects

Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1991.
Bibliography:
Includes bibliographical references (leaves 146-150).
Statement of Responsibility:
by Robert Timothy Bailey.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001683418
notis - AHZ5389
oclc - 25034722
System ID:
AA00003290:00001

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 convert-much 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 Higher-Order Elements ........................ 11
1.3 Numerical Grid Generation ........................ 11
1.4 Objectives ................................. 16
1.5 Overview .................................. 16

II DESCRIPTION OF THE LINEAR-ELEMENT 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 Doubly-Specified 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 Quadratic-Element CVBEM
Interior-Point Equation ........................ 129
B.2 Derivation of the Quadratic-Element CVBEM
Nodal-Point 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

interior-point 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 linear-element CVBEM is applied to numerical grid generation

in two-dimensional, simply-connected spatial domains.

The quadratic-element 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 well-established

finite-difference methods (FDMs) and finite-element 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 [1-6].

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 two-dimensional (2-d) 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 simply-connected to doubly- and multiply-connected 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 multiply-connected domains where such









perturbations did not appear due to a zero-flux 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, higher-order 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 2-d) or edge

(in 3-d) 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 double-valued, and the conventional

RVBEMs can produce systems of equations which are underconstrained.

In 2-d 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 1-1. 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 1-1. 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 1-2. 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 3-d 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 2-d

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 1-2.


The unknowns and corresponding number of
equations available in the standard RVBEMs for
the nine possible corer node boundary condition
combinations listed in Table 1-1.


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 corer-adjacent 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 non-linear waves, as have
Bruch and Lejeune [34] in their study of 2- and 3-d 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 3-d 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 non-conforming) 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. 1-1).

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 double-valued 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 2-d harmonic and 3-d 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 1-1. 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

methods-a 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 non-linear

nodal equation.

By formulating 2-d 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 double-valued gradients at a corner

node, one now solves for a single-valued 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 Higher-Order 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 1-2a). 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 1-2b), and the variables of interest are

assumed to vary linearly over each element. Quadratic and cubic elements are

referred to as higher-order elements, and they allow the curved structure of

boundaries to be retained. As their names imply, quadratic elements pass a

second-degree polynomial through three successive nodal points, while cubic

elements involve third-degree polynomials and four nodal points (see Figs. 1-2c

and 1-2d).

Hromadka [7] has laid the groundwork for the use of higher-order 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 finite-difference 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 1-2. 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 finite-difference equations (FDEs)

with non-uniformly 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. 1-3). 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 classes-differential 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 non-overlapping. 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 non-smooth and overlapping-characteristics which must be

corrected if meaningful finite-difference 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 1-3.


The physical and transformed domains used in numerical
grid generation (2-d 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 simply-connected 2-d spatial domain onto a unit disk and for

mapping doubly-connected 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 highest-order 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 quadratic-element

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 LINEAR-ELEMENT CVBEM


When solving 2-d 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 Cauchy-Riemann 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 well-known 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 k-connected 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 linear-element formulation of the CVBEM transforms the Cauchy
integral formula into a BEM by using two major approximations. First,

boundary F is discretized into N finite-length 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. 2-1, 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) = ( +1-Z) / (z+1-z) 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 first-order CVBEM approximation of w can be

























































FIGURE 2-1.


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 (zi1-z) + iO(j +1, j; 0) (2.6b)
L (z o) J (z -zo)

This formula forms the basis for the linear-element CVBEM. (See Fig. 2-1 for
the definition of O(j + 1, j; 0).)
Eqn. (2.6), being expressed in complex variables, actually embodies two
equations-one 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. 2-1.
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+l-zj,
27= 1- zn-x J + OjO(j+1, j-1;j)


+ E {,+CA + i+ Ci ,,C4 iC3}} (2.8)
i,i+lij

and

= In ; -z(j+1,j-1; 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 = [(xi-xi+i)C (yj-y;+1)D] (2.10c)

C4 = [(xj--xi+)D + (y,-yi+,)C] (2.10d)

C = [A(xi+, xi) + B(y+,1 yi)] / F (2.10e)









D = [B(xzi+,-x) A(y,+1-y)] / 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 j, becomes "j.


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









+= 2-1f(1 + }1Iz-z .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(,-)),_} \ z-z.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, j-1; 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, j-1;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, j-1;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, J-1; 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, j-1;j)

N
+ i {+A + ,+ACc iC4 AiC} (2.19)
i,i+1-j


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 = Eqns. (2.8) and (2.9) to yield


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 i-1 z -z O(j+1, j-1; 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 2-1 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 2-1. 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 2-1 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 diagonal-a 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 Cauchy-Riemann 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 finite-difference formulae. For example, a second-order accurate, three-
(go
point backward finite-difference approximation to at node j is given by



(a tyj[(sj-sj- 2) (Sj--Sji)] -) j_t-i(,-s2j)2
s (s -sj_-)(s i-i-2)( i-_ 2)


Sj2(S -- j_ 1)(2.24)
(Si-_Sj_-1)(S-Sj-_2)(Sj-1_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. One-sided
formulae should be used at corner nodes and at any others where the normal
gradients are different on either side of the node. Central-differencing 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 finite-difference 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 linear-element 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 linear-element 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 finite-length 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. 3-1, 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 = 2j-1






37








r




Zk-2
8(k+2,k-2;k) Zk-1 \ 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 3-1. Boundary discretization and angle definitions in the
quadratic-element 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 second-degree Lagrange polynomials of the form


( k 2)( k -1)
(Zk- Zk 2)(Zk Zk -

(zk+1--Z) (Zk+2-Z)
(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 second-order CVBEM approximation of w can be

expressed as


(Zo) = I j zo
r


zo E Q zo o r


As in the linear-element 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 +2-zk)(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. 3-1 for the definition of 8(k + 2, k; 0).
Eqn. (3.6) is analogous to Eqn. (2.6) in the linear-element CVBEM, and,
like Eqn. (2.6), it embodies two equations-one 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 quadratic-element analog to Eqn. (2.7) which can be applied at the
"end" nodes (zo = zk) is derived in Appendix B, and is given by












Fk-2 k-2 + Fk-I k-1 + 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 = 2i-1


where

F z fk-2zk_- +Zk-2
F-2 I 2(zk -I1 Zk 2)


Fk- ( --Zk-2 )
2(zk: --)(Z-1 -+ -2)




Fk = { (Zk+2-zk) + io(k +2z, k-2; k)
(zk 2 Zk)


S3z 2zk -k_ 2
2(zk zk 1)


3Zk 2k+ z-k+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
( + 2--2--------+ -!-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+2-zl)(ZI+I-z1)(z1+2-zg+l) (3.7j)
_(__ _I- ) (zI+2-zk)
h, = In (T +2 zk) in ( Zk) +iO( + 2, ; k) (3.7k)
(ZI Zk) I (Z Zk)

See Fig. 3-1 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 equations-one 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 k-2 2 + F-2k2
+ F I1 Ok-1 k+ F-1 k-1









+ 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=2j-1 + C2 +2 1+ '+2 0+2 (3.8)
1= 2i-1


and

k -Fk_ k-2 Fi-2 Ok-2
Ok = F kr2 -
+ F-_, 1k-i F_-1 _k-1

+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=2i-1


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 linear-element CVBEM are equally applicable

when quadratic elements are used, with some minor modifications. The four

types of boundary conditions-Dirichlet, Neumann, Robin, and stream

function-and 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


F-2 k-2 + F-2 k-2


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


Fk-2 k-2 F-2 k-2

Fk-1 k-I 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=2j-1 + \ +2 1+2 O(+2 01+2J (3.11)
1= 2i-1


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~k-2k-2 k-2 k-2

+ F Ck-1 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+, +

=2-1 + 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











F-2 k-2 k -2 k-2
Fk-1 k-1 k1 Ok-I
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=2j-1
I = 2i-1


+ 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


Fk2-I + FR-
k-2 Ok-2 k -2


Fi + F R
Fk-1 k-1 -+ 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=2j-1 1+2 0'1+2 1+2 1+2 (3.14)
1=2i-1 I


k= -


Fk-2 k-2 F-2k-2
+ FkR_1 k-1 F-1 _k-1

+ 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


F-2 2k-2 + Fk-2 k_-2

F-1 Ok-1 + Fk-1 Ok-I
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=21-1 + 2+ 21+2 + 0, +2 (3.16)
S= 2 -1


k = 2 -











k 2= F k-2k-2 F-2 k-2
+ Fi1 k-1 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= 2i-1


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 linear-element 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 linear-element 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 linear-element or quadratic-element 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 double-valued at a corner node.
n
The CVBEM avoids this problem by solving for 0 instead of Since 0 is both

continuous and single-valued 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. 4-1, 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 "doubly-specified" 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. 4-2. Node 3 is the





















































FIGURE 4-1. 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 4-2.


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 finite-difference expression, e.g.,


5: -= 07- 0z-6


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 finite-difference 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 Cauchy-Riemann conditions, and a finite-difference









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 Iz-1-Zj-21


Once again, other more accurate finite-difference 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. 4-3. 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 4-3. 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 Doubly-Specified BCI Nodes

When using the CVBEM to solve potential problems with corners in the

boundary, the same nine traditional boundary-condition combinations in Table

1-1 are possible. As in Chapter I, the optimal case where both the upstream and

downstream conditions are known at the corner is considered. Non-optimal-case

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

"doubly-specified" nodes, the unknowns in Cases 1 through 9 are listed in Table

4-1.

Since the CVBEM corerr" problem is due to the way that Neumann and

Robin boundary conditions are implemented and is not due to the double-valued

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 "doubly-specified" corer nodes

and "doubly-specified" 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 identified-Dirichlet nodes,

Neumann nodes, Robin nodes, and stream function nodes. A fifth type of node-

the completely specified node-was 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 "doubly-specified"

corer node which is the first node of a Neumann or Robin stretch of boundary is















TABLE 4-1.


The unknowns and corresponding
available in the implicit CVBEM
"doubly-specified" corner node
combinations listed in Table 1-1.


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. 4-4).

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. 4-5).

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 Neumann-Robin/Dirichlet node.

To summarize, it has been shown that, of the nine corner cases, seven can

be handled using previously defined nodal types. Only two-Cases 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 4-4.
















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 4-5. Diagrams showing the location where a boundary condition
change is assumed to take place for a BCI node and a
"doubly-specified" corner node.









for both the linear- and the quadratic-element 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 Neumann-Robin/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 Neumann-Robin/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 classes-algebraic methods and differential equation methods. The

algebraic methods, being much faster, tend to be the methods of choice for

unsteady, deforming-boundary 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, non-overlapping grid lines essential to

obtaining accurate FDM solutions.

In this chapter, a variation of the linear-element CVBEM which can be

used to generate grids for 2-d, simply-connected 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[(xo-xj)D +



j[(xo zxj + )D +
+,[(xo-x1j+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]

(yo-y +i+)D]



- (yo- yj)D]

+ (yo- y)C]

- (yo-yj+,)D]

+ (yo yji+)C]


C = [A(x+1-x,) + B(yj+,-yUj)] / F

D = [B(xz+lI-.) A(yj+l-y)] / F

F = (xj+l-xj)2 + (yj+i-y))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,+1-O) + D(j+I-,)

C2,j = D(Oj- j+1) + C(Oj+,1-j)
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 unknowns-a-o

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 single-point 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 2-d, 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 non-smooth

curves. Examples of such domains are given in Fig. 5-1.

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 2-d,

simply-connected, physical spatial domain onto a 2-d, rectangular transformed

domain. The boundary of the physical domain will be mapped one-to-one to the

four sides of the transformed domain, forming a so-called body-fitted coordinate

system [46]. Since the grid lines of opposite families connect opposing sides of

the rectangular transformed domain (see Fig. 5-2b), 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 simply-connected 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. 5-3. These boundary conditions force boundary






















































FIGURE 5-1.


Simply-connected 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 5-2. 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 5-3. 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 orthogonally-excellent 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 equally-spaced positions.

In order to facilitate the presentation of the grid generation process, some
new notation is in order; see Figs. 5-4 and 5-5. Previously in the linear-element
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. 5-4. 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. 5-5. 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 5-4. The notation used for the nodal points in the CVBEGGM.



















(x(8,JL),y(8,JL)


I. IL
2~~Y ~


FIGURE 5-5. 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 -(Xmi-Xmi-) + (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,k-1))2 + (y(1,k) -y(1,k-1))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, +1-04,.) + 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
i-1 j-1
= = -1 JL-1









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


& = (i-1) A i = 1,2,...,IL
j. = (j-l) A j = 1,2,..., JL
and
6 = 1 J1
= IL-1 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 round-off

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 near-orthogonality.

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

equally-spaced grid points along boundaries 1 and 4 (as described in Step 5), one

can vary the spacing using one-dimensional 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 simply-connected 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 6-1. 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 6-1. The circular domain used to test the quadratic-element
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 quadratic-element
CVBEM. Also, 4 and 4 were calculated at five internal points labeled "a"
through "e" in Fig. 6-1.

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 linear-element and










quadratic-element formulations of the CVBEM. The results for the percent error

in at the boundary nodes are given in Fig. 6-2. 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 finite-difference 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. 6-2, 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 6-1,

and the same improvement in the accuracy due to the quadratic elements is

apparent.

The computational times for the linear- and quadratic-element solutions

referred to in Fig. 6-2 were 19 and 35 seconds, respectively, on an 80286/80287

based personal computer. However, it was found that a quadratic-element

solution using only 16 equally-spaced nodes still possessed a smaller maximum

error in than the 32-node linear-element solution. The run-time for this 16-
On
node quadratic-element solution was 11 seconds.

Additional examples of solutions generated using the quadratic-element

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 6-2. 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. 6-1 imposed with the complex
potential given by Eqn. (6.2).















TABLE 6-1. The percent error in 4 and 0b at the interior points for the
CVBEM solution on the domain of Fig. 6-1 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. 6-3 and 6-4 were taken from the paper by Walker and Fenner [32]

referred to in Chapter II. For Corner Example 1 (Fig. 6-3), 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. 6-4), the somewhat more complicated

expressions


S= ln{[(x-2)2 + (y-2)2]-0.5} + x2-y2,

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. 6-5 and 6-6. 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 6-3. The geometry and nodal discretization for Corner Example 1.
























































0 1I
0.0 1.0


FIGURE 6-4. The geometry and nodal discretization for Corer Example 2.


0.51


0.0.




-0.5-


-1
-----.0
-1.0


'