Nonisoparametric formulations for the three-dimensional boundary element method

Nonisoparametric formulations for the threedimensional boundary element method M. H. ALIABADI

Computational Mechanics Institute, Ashurst, Southampton S04 2AA, UK W. S. HALL

Department of Mathematics and Statistics, Teesside Polytechnic, Middlesbrough, Cleveland TS1 3BA, UK The paper investigates the inaccuracies arising from discretization procedures in the threedimensional Boundary Element Method for potential problems, particularly from the representation of the surface and unknown function by interpolating shape functions. Nonisoparametric formulations are considered; that is, sub-parametric in which the surface is not as well represented as the unknown function, and super-parametric in which the surface is better represented than the unknown function. These nonisoparametric formulations are based on 8node and 9-node quadrilaterals which are also used to provide two different levels of isoparametric formulations. Numerical results based on the various formulations are presented for some problems for which exact solutions are known. Firstly, Dirichlet problems of finding the surface charge for a unit sphere with given uniform and nonuniform potentials corresponding to first and second order harmonics; secondly, a Neuman problem of finding the velocity potential on a sphere placed in a uniform stream. The results show that there is an advantage in accuracy and reduced calculation time for the super-parametric formulation compared to iso-parametric formulations. The sub-parametric formulation in which the unknown function is better represented than the surface is less accurate and offers no time reduction. It also presents special difficulties in defining the position of any extra unknown function node not corresponding to a surface node. Key Words: boundary element method, shape functions, potential problems, sub-parametric, super-parametric, iso-parametric, quadrilateral elements.

INTRODUCTION Inaccuracies in the Boundary Element Method arise from numerical integration and discretization procedures. Numerical integration procedures have had much attention, for example by Lachat and Watson x, Jeng and Wexler2, Pina, Fernandes and Brebbia a and Aliabadi, Hall and Phemister4, however in this paper it is inaccuracies associated with discretization procedures which are investigated. Such discretization inaccuracies may be attributed to the representation of both the unknown function and the curved surface using interpolating functions. By refining these representations, a higher degree of accuracy is achievable, although in practice computational time and effort are as important as the accuracy. For the Boundary Element Method, since the problem is defined on the boundary, surface representation would appear to be more important than say in the Finite Element Method. Iso-parametric, sub-parametric and super-parametric formulations based on 8- and 9-node quadrilateral elements are considered in an attempt to see whether Accepted May 1988. Discussion closes May 1989.

198

Enoineerin9 Analysis, 1988, Vol. 5, No. 4

representation of the unknown function or of the surface geometry is the most important factor affecting the accuracy of the solution. ACCURACY OF ELEMENTS IN REPRESENTING CURVED SURFACES Before going on to the larger question of the effect of different representations on the full Boundary Element formulations it is possible simply to examine how well 8node and 9-node quadrilaterals, both of which have quadratic behaviour on their boundaries, represent surface geometry. Two types of geometrical representation error can be tested for a known surface. The first is the difference between known surface points and the points derived from the shape function interpolation. The second is the amount by which the elements fail to join smoothly, measured by the angle between the tangents at their boundaries. The shape functions for 8-node and 9-node elements are tested by fitting the surface of a unit sphere. This is a demanding surface for a polynomial fit and there is also a simple measure of geometric error which is the difference

© 1988ComputationalMechanics Publications

Nonisoparametric formulations: M. H. Aliabadi and I4/. S. Hall from unity of the radius to any interpolated point. The surface of the sphere is divided into six identical elements, each of the element sides being part of a great circle on the sphere. A further subdivision of each element into four is made, giving 24 elements. This process is described fully in Aliabadi 5. In Table 1 both types of error are shown for the elements investigated. The superiority of the nine-node quadrilateral is obvious, both for radius and for angle between two surfaces. It is almost four times better than the eight-node quadrilateral with only 30 % more nodes in total for the finer subdivision. It is worth noting that AliabadP also found that the nine-node quadrilateral is three times better than a sixnode triangular element with the same total number of nodes. NONISOPARAMETRIC BOUNDARY ELEMENT FORMULATIONS

The bounary element formulation in discretized form for potential problems may be written as 5'6 ,=1 ~=1 [t~nJ

1

1 J(¢,tl)Ml(~,q)G[x(~,q),x '] dCdr/

= ~(x')~(x')+ d?l n=1

J(~, q)M'(~, q)

1=1

[x(~, t/),x'] d~ dq

1 , -1

(1) in which the unknown function q~ and its normal derivative Od?/Onhave been represented over each element by L

4,(¢, q) = ~ M'(¢, n)4,'

(2)

1=1

c3n-(~' q) = ~ Mt(~' q)

(3)

1=1

and N is the total number of elements. G is the potential kernel, J is the Jacobian of the transformation K

x(¢, q)= ~ Nk(~, q)x k

(4)

k=l

which is also the representation of the curved surface over each element. M(~,q) and N(¢,q) (dropping the superscripts) are polynominal shape functions associated with the curvilinear coordinates (~, q). If L = K and M and N are the same, the element is called iso-parametric. If M is a higher order polynominal than N, the element is called sub-parametric and if N is of higher order than M then the element is called super-parametric. It can be seen from equation (1) that a higher order of approximation to the surface may be achieved with a small increase in overall computation cost since the associated increase in K only has an effect of adding one term in equation (4) to the computation of points x(~, q) on the approximated surface and on the computation of the Jacobian J(¢, q). On the other hand a higher order of approximation of ~bor Odp/Onis associated with additional

unknown function values. The corresponding larger set of linear algebraic equations must eventually be solved at a cost which may vary with the cube of the number of unknowns.

NUMERICAL TESTS

A series of simple but systematic tests, or numerical experiments, has been carried out involving shape function representations for 8-node and 9-node quadrilateral elements, in the combinations given in Table 2. The various formulations tabulated are compared with the exactly known solutions for test problems of both Dirichlet and Neumann type. These problems are: (a) Surface charge for a unit sphere with uniform potential of 4n, which has a uniform solution of unity5. (b) Surface charge for a unit sphere with nonuniform potential, corresponding to first and second order spherical harmonics 5. For the first order spherical harmonics the given surface potential of cos ~b has a surface charge distribution 3cos~b/4rc. For the second order harmonics surface potentials of 3 cos q~ sin q~ cos 0, (3 cos 2 ~b- 1)/2 and 3 sin 2 q~ sin 2tk have surface charge distributions obtained by multiplying the appropriate potential by 5n/4. In these expressions ~b is the elevation angle and 0 the circumferential angle on the sphere. These problems with varying surface distributions provide a more demanding test than the uniform distributions of problem (a). (c) Flow past a sphere 7,s is a Neumann problem, in which the given potential flux - c o s 0 has a known corresponding potential - i / 2 c o s 0 as given by Lamb 9. In the three problems described above, the sphere is discretized in the same way as before. That is, into 6 equal elements each divided into four to give 24 elements each of which are further divided into four to give 96 elements. Once again further details of this discretization may be found in Aliabadi 5. Table 3 gives the number of nodes associated with the different levels of discretization for 8and 9-node elements. Results are now given for the three problems and the formulations A, C and D. The subparametric case, B, is considered separately later since it presents special difficulties. For problem (a), Fig. 1 shows the maximum error between the calculated surface charge and the exact surface charge at a node for 6, 24 and 96 elements on the sphere. There are marked improvements as more elements are used to define the surface of the sphere. Generally the results show the super-parametric formulation C being nearly as good as the 9-node isoparametric formulation D. For problem (b), Figs 2-5 show the error using formulations A, C and D for the first and second order spherical harmonics. The maximum error of the surface distribution, normalized with respect to the largest value of the surface distribution, is presented for 6, 24 and 96 elements on the sphere. There are marked improvements as more elements are used to define the surface of the sphere. Generally the results show a similar pattern to the

Enoineering Analysis, 1988, Vol. 5, No. 4

199

Nonisoparametric .formulations: M. H. Aliabadi and W. S. Hall

Isoparametric

-I 0

(8-node)

Superparametric Isoparametric

(9-node)

© D @

-2

0

Log10E

e -3

I

u

24"

G

I

!

i

96

number" o£ elements Fig. 1.

Maximum nodal error of the surface charge on a sphere. Uniform given potential

-I

Isoparametric

0 D

(8-node

Superparametric Isoparametric

-2

0 F1

@

(9-node)

O

e

Log10E

-3~

l

l

G

|

J

2"tnumber oF elements (a)

Fig. 2.

Maximum nodal error of the surface distribution on a sphere. Given potential cos c~

-I

Isoparametric

0 [3

(8-node)

Superparametric 0

Isoparametric

0 [2

(9-node)

0

-2 0

log10E

-3

!

'6 Fig. 3.

200

24

' number oF elements

Maximum nodal error of the surface distribution on a sphere. Given potential 3 cos (a sin c~ cos 0

Engineering Analysis, 1988, Vol. 5, No. 4

'$6

Nonisoparametric formulations: M. H. Aliabadi and W. S. Hall 0 D

-I

0

Isopar~metric (8-node)

0

Superparametric

[]

Isoparametric (9-node)

<>

-2 0

log10E

-3

I

!

G

II

!

l

24

9G number oF elements

Fig. 4.

Maximum nodal error of the surface distribution on a sphere, Given potential (3 cos2 gS- I)/2

-I

0

Isoparametric (8-node)

[]

Superparametric

0 []

isoparametric (9-node) -2

log 1@E -3

I

|

G

L

1

2÷

number oF elements Fig. 5.

Maximum nodal error of the surface distribution on a sphere. Given potential 3 sin 2 dpsin 20

-I.

0

Isoparametric (8-node)

0

Superparametric

r-I

Isoparametric (9-node)

<>

--2.

0

log10E

-3

!

I

S

'

~

'

I

2~

!

9G

number oF elements (a) Fig. 6.

Maximum nodal error for the potential on a sphere in a uniform flow

Engineering Analysis, 1988, Vol. 5, No. 4

201

Nonisoparametric formulations: M. H. Aliabadi and W. S. Hall Table I. Errors of difjerent elements in fittim3 the surface of a unit sphere 8-node

9-node

No. of elements

No. of nodes

Greatest error in radius

Greatest error in angle in degrees

No. of nodes

Greatest error in radius

Greatest error in angle in degrees

6 24

20 74

0.163 0.010

49.69 7.34

26 98

0.015 0.003

I0.73 3.64

Maximum error at the centre

0

Maximum error at the centre

0

Maximum error at node

[]

Maximum error at node 0

[]

[] -1

[q

-1

-2

-2

fl

log10E

l ° g 10 E -3

-3

number

oF

number

elements

Fig. 7. Comparisons of maximum error between the centre values and nodal values for an 8-node iso-parametric formulation. Uniform (liven potential

-1

'G

24

'6

Maximum error at the centre

0

Maximum error at node

[]

oF

elements

Fig. 9. Comparisons of maximum error between the centre values and nodal values for an 8-node iso-parametric formulation. Given potential 3 cos 49 sin 49 cos 0

Maximum error at the centre Maximum error at node

0 []

0 D

o

n

-I

o E1 o 13

-2

l°gloE

-2 l ° g IO E

-3

-3.

'6 number

~4 oF

elements

'6 number

'2÷ oF'

elements

Fig. 8. Comparisons of maximum error between the centre values and nodal values for an 8-node iso-parametric formulation. Given potential cos dp

Fig. 10. Comparisons of maximum error between the centre values and nodal values for an 8-node iso-parametric formulation. Given potential (3 cos 2 49- 1)/2

uniform p o t e n t i a l test results with the s u p e r - p a r a m e t r i c f o r m u l a t i o n C being nearly as g o o d as the 9 - n o d e i s o p a r a m e t r i c f o r m u l a t i o n D. O n e difference however is for the sphere discretized into 6 elements for which the

s u p e r p a r a m e t r i c f o r m u l a t i o n C gives results close to the 8-node i s o - p a r a m e t r i c f o r m u l a t i o n A. P r e s u m a b l y in this case there are so few d e m e n t s t h a t the function a p p r o x i m a t i o n is now insufficient to represent its

202

Engineering Analysis, 1988, Vol. 5, No. 4

Nonisoparametric formulations: M. H. Aliabadi and W. S. Hall Table 2.

Element types for cases studied

Case

Surface

Unknown function

Comment

A B C D

8-node 8-node 9-node 9-node

8-node 9-node 8-node 9-node

iso sub super iso

Maximum error at the centre

DIFFICULTIES OF THE SUB-PARAMETRIC FORMULATION

0

In the sub-parametric formulation B for curved boundaries, difficulties arise in placing the extra nodes

Maximum error at node -1

solutions. Fig. 6 shows the maximum error at nodal points. In the same way as for the Dirichlet problems there is a clear improvement coming from a better surface representation, and little difference in accuracy between the 9-node isoparametric formulation D and the superparametric formulation C.

9th node on the sphere

o

0

9th node at the centre of 8-node element

0 ta

-2

0

log 10 E -I

-3

0

-2

'6 number oF elements

24

log10E -3

Fig. 11. Comparisons of maximum error between the centre values and nodal values for an 8-node iso-parametric formulation. Given potential 3 sin 2 c~sin 2 0

9th node on the sphere

'6 " 24 number oF e L e m e n t ~ Fig. 13. Comparisons of errors for different position of the 9th node in the sub-parametric formulation. Given potential cos 49

0

9th node at the centre of 8-node element

-1

0

9th node on the sphere

0

m

0 9th node at the centre of 8-node element

0 0

-2

[]

log10E -I

-3

0 ta

-2

'6 number

Log10E oF

elements

~Fi9. 12. Comparisons of errors for different position of the 9th node in the sub-parametric formulation

variations. Nevertheless it can be seen that good surface representation is generally more important than function representation. For problem (c), the Boundary Element formulations A, C and D are again used to provide alternative

-3

24 number

oF

elements

Fig. 14. Comparisons of error for different position of the 9th node in the sub-parametric formulation. Given potential 3 cos dpsin cpcos 0

Enyineerin9 Analysis, 1988, Vol. 5, No. 4

203

Nonisoparametric )brmulations: M. H. Aliabadi and W. S. Hall Numbers of nodes and elements on the sphere

Table 3.

No. of elements

No. of nodes 8-node elements 20 74 290

6 24 96

9-node elements 26 98 386

9th node on the sphere

0

9th node at the centre of 8-node element

-I.

Q --2.

log 1 oE -3.

'6 number

d24 oF e l e m e n t s

Fig. 15. Comparison of error for different position of the 9th node in the sub-parametric formulation. Given potential (3 cos 2 ~p- 1)/2

9th node on the sphere

F o r the Dirichlet and N e u m a n n potential problems considered, it is seen that representation of the surface geometry is more important than the u n k n o w n function. Hence there is an advantage in a super-parametric formulation which here gives accuracies approaching that of a higher order iso-parametric formulation with a 30% reduction of calculation nodes on a sphere.

REFERENCES

-2 log10E

-3

's

~4 oF e l e m e n t s

Fig. 16. Comparisons of error for different position of the 9th node in the sub-parametric formulation. Given potential 3 sin 2 ~ sin 20

204

CONCLUSIONS

0

9th node at the centre of 8-node element

number

required by the higher order approximation of the function. In the examples being considered the geometry of the problem is k n o w n exactly (i.e., a sphere), therefore there are two possible positions for the extra node required by the 9-node representation of the function. The first possible position being the centre of 8-node element, obtained by taking ((, r/) = (0, 0), so that the extra node is on the modelled surface. The second possibility is to move this centre node along its radius so that it is placed on the actual surface defined by the spherical co-ordinates. As the result the extra node will not lie on the modelled surface. It is worth noting that the m a x i m u m error on the surface of the sphere represented by the 8-node shape function is expected to be higher than that obtained at the nodal points. Indeed the m a x i m u m error on the surface represented by the 8-node shape function is expected to lie at the centre of the element as shown by Aliabadi 6. Figs 711 c o m p a r e the errors at the centre to those found at the nodal points using the 8-node iso-parametric formulation and confirm that the nodal error is slightly lower. Figs 12-16 c o m p a r e the m a x i m u m nodal error obtained by the two possible placings of the extra node. It seems that the most suitable placing is that of the centre of 8-node element. However both results are similar to those obtained at the centre of 8-node iso-parametric formulation A. Moreover the increase in the computational effort is significant, since there is an immediate increase in the size of the coefficient matrix because of the extra node required for the function evaluation.

Engineering Analysis, 1988, Vol. 5, No. 4

1 Lachat, J. C. and Watson, J. O. Effective numerical treatment of boundary integral equations: a formulation for elastostatics, Int. J. Numer. Methods Eng., 1976, 10, 991-1005 2 Jeng, G. and Wexler, A. Iso-parametric, finite element, variational solution of integral equations for three-dimensional fields, Int. J. Numer. Methods Eng., 1977, ll, 1455-1471 3 Pina, H. L. G., Fernandes, J. L. M. and Brebbia, C. A. Some numerical integration formulae over triangles and squares with a 1/R singularity, App. Mathl Modell., 1981, 5, 209-211 4 Aliabadi, M. H., Hall, W. S. and Phemister, T. G. Taylor expansions for singular kernels in the Boundary Element Method, Int. J. Numer. Methods Eng., 1985, 21, 2221-2236 5 Aliabadi, M. H. The treatment of problems in the three-dimensional Boundary Element Method, MPhil Thesis, CNAA, London, 1981 6 Aliabadi, M. H. Singular Kernel integration and surface representation in the three-dimensional Boundary Element Method, PhD Thesis, CNAA, London, 1985 7 Jaswon, M. A. and Symm, G. T. Integral Equation Methods in Potential Theory and Elastostatics, Academic Press, 1977 8 Aliabadi, M. H., Hall, W. S. and Phemister, T. G. Taylor expansions in the Boundary Element Method for Neumann problems in Boundary Elements VII, (Eds C.A. Brebbia and G. Maier), Proceedings of the 7th Int. Conf. Como, Italy, Computational Mechanics Publications, 1985, 12-31 to 1~39 9 Lamb, H. Hydrodynamics, Cambridge University Press, Cambridge, 1945