Fundamental solutions and boundary element methods

Fundamental solutions and boundary element methods J. C. W U

Professor of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA 30342, USA

INTRODUCTION

The majority of problems of modern engineering are concerned with physical phenomenon describable in terms of fields. Field problems are familiarly described by partial differential equations, either homogeneous or inhomogeneous, subject to certain prescribed boundary conditions and, for some problems, also initial conditions. Historically, classical mathematical analyses were the only effective avenue of obtaining solutions to field problems. In recent decades, however, the advent of high-speed digital computers has opened up a new passage to these solutions. Indeed, at the present, a preponderance of engineering analyses is performed via computations rather than classical mathematics. It is well known, however, that computational methods can, and in many cases must, benefit from classical mathematical analyses. In the present paper, the role of 'fundamental solutions in boundary element methods' is discussed. It needs to be emphasized, even at the risk of appearing trivial and repetitive, that fields are not merely solutions of differential equations. No meaningful solution of a physical problem is possible without an adequate knowledge of the boundary conditions. Improper or imprecise numerical treatment of boundary conditions invariably leads to unreliable or unacceptable solutions of field problems. With the boundary element approach, the differential equations describing the fields are reformulated as integral relations containing boundary integrals. Numerical procedures are then developed to process boundary information utilizing the boundary integrals. This distinguishing feature of the boundary element approach, that it processes boundary information utilizing the boundary integral, offers an opportunity for developing numerical procedures that possess superior solution efficiency and accuracy. 1'2 It has been found that the fundamental solutions, in addition to their mathematical attributes, offer important physical insights which are invaluable to the development of boundary element procedures. In the present paper, the flow of a viscous fluid is chosen as a focal problem to describe physical perceptions associated with the boundary element approach. In particular, elliptic differential equations are examined in the present paper. Issues discussed here are, of course, common to other applications such as elastostatics and electromagnetism.

that, under quite general circumstances, the flow problem can be partitioned into a kinematic aspect and a kinetic aspect. The kinematic aspect of the problem consists of the definition of derived field variables such as the rotation field, or vorticity, te and the dilitation field, or source, g, viz. V.v=g (1) V'v = ~

(2)

where v is the velocity field. Equations (1) and (2) are valid for compressible and incompressible flows, steady and unsteady flows, turbulent and laminar flows, internal and external flows. It is well known that the velocity field can be decomposed into a solenoidal part and an irrotational part. The solenoidal part possesses a vector potential and the irrotational part possesses a scalar potential. Both the scalar potential and the vector potential satisfy Poisson's equations. Equations (1) and (2), therefore, is a set of elliptic vector equations. The kinetic aspect of the flow problem is concerned with the various kinetic processes -convection, diffusion, production, dissipation, etc. - that govern the distribution of field variables such as the momentum, the energy, the dilatation and the vorticity. This kinetic aspect is described by differential equations known as transport equation. For brevity, only the vorticity transport equation for a steady incompressible laminar is studied in the present paper, although references are made from time to time to unsteady flows. Since the transport of vorticity is similar to the transport of other field variables, the vorticity transport equation is similar to the transport equations for other field variables. The discussions of the present paper are, therefore, useful under quite general circumstances. Through the use of fundamental solutions, Wu3 reexpressed the kinematic and kinetic differential equations in the form of integral representations for the field variables, velocity, vorticity, etc. These integral representations as well as their applications to various types of flows, are summarized in a series of earlier articles. o '-S A general expression for the velocity field 3 is

~v(r)=-fgoFo ~o- I ~0x FodR0 R

R

+ ~(vo'no) F dBo--~(Vo x no) x Fo dBo FLOW EQUATIONS IN INTEGRAL FORM

The motion of a viscous fluid is described by the continuity, the Navier-Stokes and the energy equations which are well known in their differential form. Wu3 has shown Accepted August 1986. Discussion closesMay 1987.

2 Engineering Analysis, 1987, Vol. 4, No. i

B

(3)

B

where R is a region of interest, B is the boundary of R, n is a unit normal vector directed outward from the region R, r is a position vector, the subscript '0' indicates a variable, or a differentiation or an integration in the ro space, e.g. go = g(ro), fl is the angle of the infinitesimal region sur0264-682X/87/010002-05 $2.00 © 1987 Computational MechanicsPublications

Fundamental solutions and boundary element methods: J. C ICu rounding the point r. Thus, for a three-dimensional problem /3 = 4rt if r is an interior point in R and/3 = 2rr if r is on a smooth boundary. For a two-dimensional problem, ~ = 27r if r is an interior point in R and ~ = ir if r is on a smooth boundary. The function F 0 is the gradient, in the re space, of the fundamental solution of scalar elliptic equations multiplied by 41r or 2rr, depending on whether the problem is three-dimensional or two-dimensional. Specifically, r0 - r Fo - - Iro--rl a

(4)

where d is the dimensionality of the problem, i.e. d = 3 or 2 depending on whether the problem is three-dimensional or two-dimensional. The function Fo in fact may be viewed as the fundamental solution of the elliptic set of vector equations (1) and (2). It is easy to show that the differentiations of equation (3) in the r space give equations (1) and (2). Equation (3), therefore, expresses the kinematics of the general flow problem in the form of an integral representation for the velocity vector. Equation (3) expresses the velocity field as the sum of four contributions. Each of these contributions is represented by an integral and is a velocity field by itself. The first integral in equation (3) is the contribution of the source field to the velocity field. This contribution is an irrotational velocity field. The second integral in equation (3) represents the contribution of the vorticity field to the velocity field. The integral is a generalized statement of the well-known law of Biot-Savart. The Biot-Savart integral gives a solenoidal velocity field. The source and Biot-Savart integrals are integrals over the region R. The form of these two integrals shows that the contributions of the source and vorticity fields to the velocity field are not dependent upon the geometry of the boundary B. These contributions are the same whether or not the boundary is present. In consequence, these contributions may be computed as ifR is an infinite unlimited region. Relatively simple numerical procedure can, therefore, be utilized to evaluate these contributions. The third and fourth integrals in equation (3) are boundary integrals representing respectively the contributions of the normal component and the tangential component of the velocity boundary condition to the velocity field in R. Each of these contributors gives a solenoidal and irrotational field in R. The integrand of the third integral in equation (3) is similar in form to that of the first integral. In consequence, the contribution of the normal velocity boundary condition to the velocity field in the region R is equivalent to that of a sheet of concentrated source, of strength --v.n lying on the boundary B. The negative sign for this source strength is due to our sign convention, the normal direction is considered positive when directed away from R. The integrand of the fourth integral in equation (3) is similar in form to that of the second integral. The contribution of the tangential velocity boundary condition to the velocity field in R is, therefore, equivalent to that of a vortex sheet of strength v x n lying on B. The boundary integrals are thus equivalent to the source and Biot-Savart integrals whose physical interpretations are more readily recognized. Wu and Wahbah 6 presented an integral representation for the vorticity vector. This integral representation, given below, is equivalent to the well-known vorticity transport equation describing the kinetics of a steady, incompressible and laminar flow.

t3ca(r)

=

- -

1H f ( nv ° x°t o °x) x FF° d R° ° + d~ B R

°

v

B

+ ~ (cao "no)F0 dBo--~ (too x n0) x F0 dBo (5) B

B

where H is the total head of the flow divided by the kinematic viscosity v. The first integral in equation (5) is a domain integral representing the contribution of the convection process to the vorticity field. The form of this integral shows that this contribution is a Biot-Savart contribution. In other words, the effects of the quantity v x to on the vorticity field in a steady flow is similar to the effect of w on the velocity field. Since the vorticity field is solenoidal, the counterpart to the source integral of equation (3) is absent in equation (5). The third and fourth integrals are boundary integrals representing the contributions of the vorticity boundary conditions to the vorticity field. The second integral in equation (5) is a boundary integral representing the contribution of the total head to the vorticity field. This integral has a special role in the boundary integral procedure. Equations (3) and (5), with g = 0, together with prescribed velocity boundary condition % on the closed boundary B of a region R, completely describe the steady, laminar, incompressible flow problem. These two equations describe respectively the kinematics and the kinetic aspect of the problem.

BOUNDARY ELEMENTCONCEPTION Consider a well-posed problem described by equations (1) and (2). If g and ca are completely known in a singularlyconnected region R, then the correct boundary conditions for the problem is either the normal velocity component v.n, or the tangential velocity v x n, or a linear combination of the two over the entire boundary B. 7 Equation (3), therefore, contains an unknown value of either v ' n or v x n at every point on B. In order to use equation (3) to compute the field variable v in R, these unknown boundary values must first be determined. With a boundary element procedure, one determines the unknown values of v on B by letting r = rb, a boundary point, and thus specializing equation (3) to the boundary B. This specialization yields a boundary integral equation containing unknown values of v on B. The solution of this boundary integral equation yield values of these unknowns which can then be used along with the known, prescribed, values of v on B to compute the field variable v in R, using equation (3). Consider the following equations which are components of equation (3) as applied to the boundary point rb: /~vb- ~(vob" no)Fob dBo + ~ (Vob x no) x Fob dBo = Gb B

B

(6) where the subscript 'b' is used to designate quantities on the boundary point rb, i.e. vb = v(r~), etc., the function Gb represents the contributions of the region integrals in equation (3) at the point rb. Since g and ca are known in R, Gb is considered a known function of rb. For the case where the normal velocity n~.v b is prescribed at every point on

Engineering Analysis, 1987, Vol. 4, No. 1

3

Fundamental solutions and boundary element methods: J. C lfu B, the first integral in equation (6) is known. One therefore has, from equation (6) n. (~ (Vob X no) X Fob dBo = Hb

(7)

B

/3n x vb + n x ~ (Vob x no) x Fob dBo = Jb

(8)

where Hb and Jb are known functions of rb. Equation (7) is the normal component of equation (6), and equation (8) is the tangential component. For three-dimensional problems, equation (8) has two scalar components. The three scalar equations of equations (7) and (8) contain the two components of v b x n as unknown functions. It can be shown that only two of the three scalar equations are independent. The solution of any two of the three equations yields vb x n at every point rb on B. For two-dimensional problems, equation (7) has only one component, in the direction perpendicular to the plane of the field v. Thus equations (7) and (8) form a set of two scalar equations. These two scalar equations are equivalent to one another and either equation can be solved to yield the tangential velocity component vb x n at every point rb on B. Equation (7), being a Fredholm equation of the first kind, leads to matrix equations with coefficient matrices containing zero principal diagonal elements. For this reason, it is generally preferable to employ equation (8), and not equation (7), in the computation of vb x n. For the case where the tangential velocity v b x n is prescribed on B, equation (6) yields two scalar equations for three-dimensional problems and two scalar equations for two-dimensional problems containing v b -n as the unknown function. Only one of the two equations is independent. In this case it is preferable to solve for vb" n using the normal component of equation (6), since it is a Fredholm equation of the second kind in this case. With a singly-connected region R, equations (7) and (8), or similar equations obtained from equation (6) permit the missing boundary information to be uniquely determined. For problems involving multiply-connected regions, with g and co known in R, in addition to the correct boundary conditions for the problem stated earlier, the field v is uniquely determined only if the circulations or the fluxes around the independent circuits in R are also known. 7 Information about the circulations and fluxes form auxiliary conditions which, together with boundary integral equations obtained from equation (6), permit the missing boundary information to be uniquely determined. The boundary element procedure outlined above utilizes an integral representation obtained on the basis of the fundamental solution concept and contains two major steps: a step in which missing information about the field variable on the boundary of the region of interest is determined and a step in which the field variable is computed in the interior of the region. These two major steps are the distinguishing feature of boundary element methods. They are, with appropriate modifications for specific problems involved, common to all boundary element procedures.

both the normal and the tangential components of the velocity vector on tile boundary are established from the physics of the problem. If the vorticity field is completely known in R and on B, then the problem is overspecified. If one uses known distributions of ~a and g and known boundary values of v~ to compute a velocity field based on equation (3), then this computed velocity field in general will have values on B that differ from the input v b values. The computed v b values agree with the input v b values only if the dilatation aad vorticity fields in R are suitably restricted. For incompressible flows, the dilatation is everywhere zero and the vorticity field is restricted. In computing unsteady incompressible flows, it is convenient to follow the development of the vorticity field with time. A computation loop to advance the solution from an old time level to a new time level begins with a kinetic part, where new vorticity values are computed, and follows with a kinematic part, where corresponding new velocity values are computed. In the kinetic part, new vorticity values are obtained at all interior points of the region R. The boundary information about the vorticity field, however, is not known. If one wishes to use equation (3), to compute v values, then the boundary information about co needs to be established first. Since vb is completely prescribed, the two boundary integrals in equation (3) are known. Therefore, applying equation (3) to the boundary points r b, one obtains a vector integral containing the boundary vorticity cob as the unknown function. The integral equation can be solved to determine the needed boundary vorticity information. At this point, all the integrands in integrals of equation (3) are known. Equation (3) can be used to compute v values inR. The procedure just outlined is a boundary element procedure. This procedure is similar to that described earlier. In the present procedure, since no velocity boundary information is missing, one does not compute any boundary velocity values. Instead, the missing vorticity information is computed. As discussed earlier, the boundary integral in equation (3) containing the tangential velocity is a Biot-Savart integral. Since the present computation of the missing boundary vorticity information involves a Biot-Savart integral, the comments presented earlier in connection with the computation of the missing tangential boundary velocity are directly applicable to the computation of the boundary vorticity values. To demonstrate this conclusion, consider a thin layer adjacent to the boundary. The vorticity in this thin layer is approximated by a vortex sheet adjacent to the boundary. The strength 7 of this vortex sheet is the integrated value of e~ across the layer. This strength contains the unknown boundary vorticity information and needs to be established. One 7 is computed, the total vorticity represented by 7 may be suitably distributed across the layer to yield the boundary vorticity information. One obtains the following two equations by specializing equation (3) to the boundary point rb and writing the resulting equation in its component form :

n- + 70 x Fob dBo = Kb KINEMATICS In the foregoing discussions, it is assumed that the vorticity field is known in R and on B and that the velocity boundary information is not overspecified. In most flow problems,

4

EngineeringAnalysis, 1987, Vol. 4, No. 1

(9)

B ÷

n x (~ 7o x Fob dBo = Lb t/

B ÷

(10)

Fundamental solutions and boundary element methods: J. C. Wu where Ko and L 0 are known functions of r 0 and the integration is over B ÷, which represents the location of the vortex sheet 3' and is adjacent to the boundary B but not coinciding with B. Equations (9) and (10) are identical to equations (7) and (8) except that the unknown function in the present case is 7 in place of v x n and that the known function Lb includes the quantity fin x v 0 which is known. Earlier discussions in connection with equations (7) and (8) concerning the inter-dependency of the component equations, the preferability of the particular equation(s) to use in the numerical procedure and the issues relating the dimensionality of the problem as well as the connectedness of the region involved are, therefore, all applicable to equations (9) and (10). It is worthy of note that equation (10) is a Fredholm equation of the second kind and is the preferred equation even though the unknown function 7 is not explicitly present outside of the integral. The fact that equation (10) is a Fredholm equation of the second kind can be demonstrated by letting B + and B coincide and consider the integral in equation (10) in the Cauchy principal value sense. Instead of employing the concept of a vortex sheet adjacent to the boundary, equation (3) can be discretized, yielding a set of algebraic equations containing velocity values and boundary vorticity values as unknowns. A subset of this set of equations involving boundary points contains only the boundary vorticity values as unknown. This subset can be solved to determine boundary vorticity values. The discussions above concerns the unsteady flow and view the boundary element procedure as a means to simulate the process of generation of vorticity on the boundary of a flowfield. Since the kinematics of both the unsteady flow and the steady flow are described by the same equations, namely equations (1) and (2), these discussions are equally applicable to steady flow problems. It is more appropriate, however, to view the boundary element procedure, when used in connection with the kinematics of the steady flow, as a part of an iteration loop. KINETICS OF STEADY FLOW

Consider the kinetic part of the iteration loop. At the beginning of each computation loop, values of v and t~ at the old iteration level are known. Boundary values of H, which are also needed in order to compute new vorticity values in R using equation (5), however, are not known. By specializing equation (5) to boundary points rb and using old iteration level to o values in the left side of the specialized equation, one obtains a boundary integral equation containing H as the unknown function. This boundary integral equation is solved to determine the quantity H on the boundary B. The quantity H thus determined is then placed into the right side of equation (5) along with old values of v and ~ to compute new values of ¢~ inR. The kinetic part of the computational loop just outlined constitutes another boundary element procedure in which some missing boundary information is first computed and the computed boundary information is then utilized in the computation of the field variable of interest. The overall problem of steady flow contains, therefore, two boundary element procedures, a kinetic one and a kinematic one. The two procedures are similar to one another and they are linked through the boundary conditions for the vorticity field. This loop is summarized in the next section of this paper.

NUMERICAL PROCEDURE

A boundary element procedure for solving two-dimensional steady incompressible flow is as follows. The boundary B is divided into boundary elements and the region R is divided into field elements, each with a number of associated nodes. Integrals in equations (5) and (3) are replaced by sums of element integrals. Eelement interpolar functions are utilized to evaluate the element integrals analytically and to express them in terms of nodal values of field variables. One obtains as a result the following sets of algebraic equations approximating equations (5) and (3). 1

1

F

11

{oa} = - [ A I {uoa} + ' [ B ] {voa} + [C] {H0) + [D] {~%} (11)

{u} = --[B] {oa} + [D] {Uo} -- [C] {Vo}

(12)

{v} = [A] {oa} + [Cl {uo} + [DI {vo}

(13)

where u, v, w, uw and vw are M x 1 matrices, M being the total number of nodes in R and on B; u and v are two Cartesian components of v, co is the magnitude of t~, which has only one component for the present two-dimensional problem. Ho, W0, ub and Vo are N x 1 matrices, N being the number of boundary nodes on B; [A] and [B] are M x M matrices and [C] and [D] are M x N matrices. The matrices [A], [B], [C] and [D] are geometric coefficient matrices. They are determinate once the element configurations in R and on B are geometrically specified. These coefficient matrices are shared by the kinetic equation, equation (11), and the kinematic equations, equations (12) and (13). Equation (11) possesses a sub-set of algebraic equation: 1 1 {oa0} = - [ A o l {wu} +--[Bo] {oar} + [C01 {Ho} 1)

P

+ [Do] (oat,)

(14)

where [Ab] and [B0] are N × M matrices and [Co] and [Db] are N × N matrices, with the subscript b denoting boundary nodes. Equations (12) and (13) possess sub-sets of algebraic equations which can be combined to give

{vt} = [A o sin ~--Bo cos c~] {co)+ [Do] (vt}-- [Col {Vn)

(15) where a is the slope of the boundary B relative to u, v t and vn are respectively the tangential and the normal velocity components of v0. An equation similar to equation (I 5) is also obtainable expressing {Vn} in terms of {co}, etc. This equation for {Vn) is equivalent to equation (15). Values of {u}, {v}, {oa), {uoa} and {voa) are assigned to begin an iteration process. At the beginning of each succeeding iteration, these values are known from the preceding iteration. These values contain the boundary values {oao}. The values of {ub} and (vo} or equivalently {vt} and {Vn} are prescribed. Each iteration loop contains the following steps: (a) Using known values of {oao}, {uoa) and {z~.o}, the matrix {Ho} is computed by solving equation (14). (b) Computed (Hb) are used together with known {uoa), {voa) and {oao}in equation (11) to compute new oa values at nodes in R away from B. There are (M--N) interior values of co computed in this step.

Engineering Analysis, 1987, Vol. 4, No. 1

5

Fundamental solutions and boundary element methods: J. C. lCu (c) Placing the computed interior values of ~ and prescribed values of (v t) and {v n) in equation (15), one obtains a matrix equation containing new (~o0) as unknown. This equation is solved to determine new ~wo) values. (d) Computed {o~b) and interior t~ values, together with prescribed {u0) and ~v0) values, are placed in equations (12) and (13) to compute new { u ) a n d {v) in values. Steps (a) and (b) constitute a boundary element procedure treating the kinetic aspect of the steady flow problem. Steps (c) and (d) constitute a boundary element procedure treating the kinematic aspect of the problem. Steps (a) and (c) require the inversion of the matrices [C0] and [Do] respectively. This inversion can be carried out using either an iterative procedure s or a direct procedure through a conformal transformation technique. 9 It should be noted that [C0] and [Do] are fully occupied matrices and they are N x N matrices, N being the number of boundary nodes. Steps (b) and (d) are explicit steps and involve no matrix inversion. The procedure just outlined has been implemented by Wahbah s who presented analytic expressions for all the geometric matrices in equations (11)-(15)using polynomial element interpolation functions. Wang 9 used a Fourier series method in combination with a conformal transformation technique to obtain a remarkably efficient and accurate procedure. Numerical results for various steady flow problems are presented by Wahbah and by Wang.

6

Engineering Analysis, 1987, Vol. 4, No. 1

ACKNOWLEDGEMENTS This research is supported by the Air Force Office of Scientific Research under Grant No. 82-0108.

REFERENCES 1 Wu, J. C. Boundary element methods and inhomogeneouselliptic differential equations, Proe. 7th lnt. Conf. on Boundary Element Methods in Engineering, 1985a 2 Wu, J. C. Boundary element methods and inhomogeneousparabolic differential equations, Proe. 1st Boundary Element Technology Conference, Springer-Verlag, 1985b 3 Wu, J. C. Problems of general viscous flow, Chapter 4 of Developments in Boundary Element Methods, Vol. 2, Elsevier Applied Science Publishers, London, 1982 4 Wu, J. C. Fundamental solutions and numerical methods for flow problems, lnt. J. for Numerical Methods in Fluids, 1984, 4, 185 5 Wu, J. C. Problems of time-dependent Navier-Stokes flow, Chapter 6 of Developments in Boundary Element Methods, Vol. 3, Elsevier Applied Science Publishers, London, 1985c 6 Wu, J. C. and Wahbah,M. M. Numerical solution of viscousflow equations using integral representations, Lecture Series in Physics, Springer-Verlag, 1976 7 Wu, J. C. Numerical boundary conditions for viscousflow problems,AIAAJ. 1976, 14 (8), 423 8 Wahbah, M. M. Computation of internal flows with arbitrary boundaries using the integral representation method, Georgia Institute of Technology Report, 1978 9 Wang,C. M. Computation of steady internal flows using integral representation methods with finite series expansions, Georgia Institute of Technology PhD Thesis, 1983