An Hp-Adaptive Hierarchical Formulation for the Boundary Element Method Applied to Elasticity in Two Dimensions

An Hp-Adaptive Hierarchical Formulation for the Boundary Element Method Applied to Elasticity in Two Dimensions

R. B. V. Pessolani

Universidade Federal Fluminense

Departamento de Engenharia Mecânica

Av. Ari Parreiras, 660

24230-323 Niterói, RJ. Brazil [e-mail]

This paper presents an HP-Adaptive Procedure with Hierarchical formulation for the Boundary Element Method in 2-D Elasticity problems. Firstly, H, P and HP formulations are defined. Then, the hierarchical concept, which allows a substantial reduction in the dimension of equation system, is introduced. The error estimator used is based on the residual computation over each node inside an element. Finally, the HP strategy is defined and applied to two examples.

Keywords: Numerical methods, boundary element methods, adaptive procedures elasticity

Introduction

The self-adaptive programs represent an area of great interest for researchers and analysts in the Boundary Element Method (BEM) field. Those techniques can be seen as ways of automation of calculation procedures, seeking for and minimizing errors. For such, the user provides the adaptive program the minimum necessary information to define the geometry, the loads and the boundary conditions, and the program tries to generate the best element mesh that gives the most accurate calculation within a certain tolerance. This procedure has the advantage of reducing the processing time, reducing also considerably the dimension of the equation system and resulting in a great computational economy. At the same time, the adaptive strategy helps the engineering in the problem analysis, avoiding occasional evaluation errors in the mesh distribution.

There are three types of Adaptive Formulations: The P Adaptive Technique increases the degree of the interpolation function in the critical areas, whereas H refines the mesh of elements, and HP is a combination of H and P techniques.

In BEM, P technique was introduced by Alarcón and Reverter (1986) and developed by Cerrolaza et al. (1988), Parreira (1987) and Umetami (1986) and others who showed that its convergence is asymptotic in areas where the solution does not present singular points.

The information obtained with P refinement is hierarchical, that is, stored from one iteration to the following, incorporating only the new refinement computation. For those, the families of Lagrangean's functions cannot be used. Although, Peano (1976) hierarchical functions were initially the most used, nowadays other families of functions such as polynomials of Chebischev, Legendre and Integrated Legendre are also employed.

The H Technique was studied previously by Rank (1989), and Rencis and Mullen (1986). These authors checked the good operation even in the presence of singularities. At first, H technique was not hierarchical. That means that in a new iteration it was necessary to recompute all the coefficients for the old and the new collocation points. Guigianni (1992) and Parreira and Dong (1992) presented also error estimators adopting a technique that allows the hierarchization of the elements, bringing large computational savings.

Babuska and Guo (1986) studied the convergence of the P technique for the Finite Element Method (FEM) with different meshes. They drew curves representing the error as a function of the number of degrees of freedom (dof) for different initial meshes, and noticed that the format is that of an inverted S, as represented in figure 1. In the first stage, there is an exponential convergence where a small increase of the number of dof produces a great reduction error. Starting from a critical zone, the curve tends to stabilize, like an asymptotic curve, or, in other words, the increment of dof does not produce substantial increase of precision. This critical point depends on the number of singularities in the solution. The fact that different meshes produce different curves leads to the verification that precision of the P version depends on the initial distribution of the mesh. The ideal curve would be the one that coincides in first phase with mesh (a) until the intersection with the curve of mesh (b) and then it would slide until the intersection with mesh (c). For that, it is necessary to include changes during the processing, which results in the HP refinement. In a different study, Zienkiewicz and Zhu (1989) proved that H Technique leads to a good precision for quadratic interpolation functions. However, for a good accuracy, it is also necessary to use the P technique with the inclusion of functions of higher degrees. Even so, in the presence of singularities, H technique is the most suitable.

Rank (1989) verified the validity of those concepts for BEM, and concluded that a combination of mesh refinement and increase of the interpolation order is the best procedure, combining the advantages of P and H Techniques. The solution converges asymptotically, even in presence of singularities. He also defined an implementation strategy for an HP Adaptive program, consisting of the following steps:

i) Choice of a basic mesh, enough to describe the geometry and boundary conditions.

ii) Separation of the elements in two categories: the one that is expected to have a smooth solution, and the other that is adjacent to points of singularity solution, which can be corners or points of boundary condition changes.

iii) Choice of initial interpolation degrees for all the elements.

iv) Solution with the BEM

v) Computing the errors for each element.

vi) If the error is smaller than the tolerance, the execution is finalized. Otherwise, if it is larger than the maximum specified, it is necessary to increase the degree of the interpolation function in the non-critical areas (P) or subdivide the elements in the areas of singularity (H). Then return to step iv.

The present work was first published as a PHD Thesis (Pessolani (1995)) and develops an HP Adaptive Program applied to the BEM, combining the hierarchization concepts of H and P, and using the error indicator of Parreira (1992).

Classical BEM Formulation

The direct BEM formulation (Brebbia, 1978) for linear two-dimensional elasticity is given by:

with boundary conditions .

Usually, inside the element, variables are interpolated according to Lagrange functions. Applying boundary conditions, u_j= in G_u and t_j= in G_t, and rearranging the equation, the system 2N_ex2N_e is set,

The system is solved by Gauss elimination or by iterative solvers, especially when there is a larger number of dof (Pessolani, 1999).

The resolution of the integral is usually made by the Gauss Quadrature procedure. When the collocation point is integrated inside the element, it is necessary to use expressions developed to calculate the singular integral such as those given by Guiggiani and Lombardi (1992) or Pessolani (1997b).

H-Hierarchical Formulation

Using Lagrange Interpolation Functions, any element subdivision or any increase of the number of collocations points requires that the system be reassembled with the calculation of all the coefficients in matrices A and b. When the process is adaptive, several discretizations are necessary until meeting the correct solution. For this reason, the classical method is not convenient to be applied to the adaptive procedure, because each new iteration requires a new computation of the whole system.

An alternative approach introduced in FEM by Peano (1976), and reported by Zienkiewicz, Gago and Kelly (1983) uses hierarchical families of interpolation functions. Those functions are accumulative, that is, the higher order functions are generated without modifying lower order contributions. In the rediscretization, the new system is set up, increasing the previous system with submatrices A₁₂, A₁₂, A₂₂ and the vector b₂, corresponding to the new collocation points.

In BEM, the H-hierarchical formulation was applied by Parreira (1992). The idea is to include intermediary collocation points in the element with interpolations functions of the same order of the previous, that can be the conventional Langrangeans.

For a quadratic interpolation is defined a set of interpolations functions N₂, that is:

where Dh is half the element to be refined in natural coordinates, k the number of requested divisions, h₁₁ is the natural coordinate of the nodes used in the linear interpolation and h ₂₁ is the natural node coordinate used in the quadratic interpolation.

For each degree l , 2^l-1 Bubbles functions are generated with the corresponding collocation points arranged symmetrically, equally spaced along the element. For example, the second hierarchization level would include collocation points according to the figure 2.

The approximate solution j on the element is expressed in terms of the previous and the new values on the nodes:

Generically, let be the iteration l , L the number of hierarchical levels and P_l the list of integers indicating which pair of hierarchical functions P are being applied in the mesh (for each level l there are 2^{l -1} pairs). The approximate solution is given by

Where in iteration l , it is:

Where:

t_l ,n is the natural coordinate mapped in the interval [a_l ,_n,b_l ,_n], defined by:

To accelerate the convergence rate it is quite convenient that the starting point of the adaptive scheme be a refinement that takes into account the points of singularity of the problem.

P-Hierarchical Formulation

According to Zienkiewicz and Zhu (1989), in FEM the H formulation gives for linear and quadratic functions an accuracy from 5 to 10% measured in the global energy norm, being especially efficient for p=2 (quadratic functions). However, when a higher accuracy is required, the P or H-P formulation is necessary.

So, in the present work the Hierarchical P Formulation was developed for, associated with the H procedure, optimize the resolution of the problem.

It was adopted a base of functions formed by the family of polynomials Integrated Legendre polynomials that are used in FEM, and have series of properties that gives a good numeric stability for elliptical problems.

Let the Legendre polynomials be:

where f i is defined in term of the Legendre polynomials P_i-1:

Which integrated yields:

The Integrated Legendre polynomials are:

Therefore, the interpolation functions of the Integrated Polynomials will be:

Polynomials of Legendre of higher order can be used, but for the purposes of the present work, these are enough.

To increase the order of interpolation function, a new collocation point for each added function is and placed inside the element.

HP-Hierarchical Formulation

The idea is to combine the H and P procedures in a single program. The strategy of HP formulation can be better visualized in figure 3. The initial system is composed by Lagrangean quadratic functions. After solving the system, the procedures are applied and an error indicator is applied to evaluate the solution accuracy. If a new refinement is requested, collocation points are generated with new interpolation functions, and added to the original system matrix. This process is repeated until one obtain the requested precision or a bad condition appears.

Error Indicators

The manager of the adaptive process is the error indicator that provides information about size and error distribution. If it is local, this information is called indicator. The estimator shows the measure of the global error. In the present study, Parreira's (1992) indicator and estimator was adopted, based on the calculation of the residue inside the element.

Definition of Residue

Considering the body forces bi the equation (1) can be written as:

Removing the functions arguments and isolating b_i, it comes:

To obtain solution to the equation (20), a discretization approaching the values of t and u in terms of nodal values is made:

To simplify the explanation of adaptive concepts, the volume forces will not be considered. It will be only necessary to calculate the boundary integrals.

Introducing the approximate solution in equation (20), the residue can be expressed as:

The Residue Interpolation Technique

Recently, with the H-hierarchical formulation, Parreira and Dong (1992) developed a calculation procedure involving nondimensional residues. First, they defined an expression called the residue characteristic:

This parameter is calculated in the first iteration, starting from the components of the vector b_i and is be used as an invariant during the process. Defining also the characteristic estimator to obtain a nondimensional estimator:

Being E computed according to the L2 norm:

Where Ne is the number of Boundary Elements

Let's the system of the first refinement written in the form:

In the hierarchical formulation, for the second iteration, where new degrees of freedom have been introduced, the system becomes:

where A11 and b1 are kept unchanged from eq. (26) and the solution is improved.

The sum of the residues corresponding to the increased points in the element Gk, in the direction i is:

and the total over the boundary is:

Considering that l _i^k depends on the length of the element, Parreira and Dong (1992) define the error density in the element k in the form:

where Lk is the length of the element k.

For the whole boundary the global density error is defined as:

This expression has to be always smaller than an arbitrated maximum tolerance. The estimator can, therefore, also be written in the form:

where .

To obtain the error indicator, Parreira and Dong (1992) make use of an interpolation of the residue:

where Nk are the Lagrangean interpolation functions associated to the m+1 nodes of the element Gk.

According to the Collocation Method, for a degree n the residue is null in the n+1 collocation points. So the residue is calculated for the new collocation points:

The degree functions between Nn+1 and Nm are increased in the hierarchical process. To calculate the residue in each hierarchical collocation point, Peano method is used:

That expression calculates the residue in the collocation points associated with the hierarchical functions. Notice that this coefficients are always computed by the hierarchical process.

Finally, the error indicator for each element is calculated according to the expression (25):

where ,

or ,

in the norms L2 and H1 respectively. The whole calculation process is systematized in the figure 4.

The technique utilizes the own hierarchical formulation to compute the Residue. The residue calculation is done taking for base the coefficients calculated in the hierarchical formulation. This is a great advantage and its performance will be verified.

Examples

Square Plate with a Hole in the Center

As a first example, a plate with a hole in the center, as shown in figure 5, was tested.

This plate has a solution in terms of displacements and tractions, shown in figure 6.a and 6.b. As it can be noticed, the solution is smooth in sides BC, CD and EA. There is a zone of concentration of tensions near vertex A (0,0) and E (0,-1).

To capture those particularities the refinement used in the adaptive process is illustrated in figure 7a.. To help the discretization of the hole area, a geometric processor of Pessolani(1997) was used, to evaluate the difference between the original geometry and the discretization. The resulting mesh with three elements is indicated in figure 7b. Notice that the elements near the two vertex are smaller than the central one.

To visualize the effectiveness of the formulation HP, the problem was solved with the formulations H, P and HP, and compared the final meshes obtained and the global accuracy. (Fig 8).

For the HP Technique, the strategy H was used in the element 1 – since the presence of the singularity - and in the other areas P technique was used. It can be observed how well the error indicator works, accumulating collocation points in areas where the solution requires a larger refinement, especially in the singular area.

The error percentage of the element 1 is shown in figure 9. It can be observed the better convergence of H technique compared with the strategy P. It confirms the concept that the H strategy works better in elements of singular solution. The HP-Technique has the same rate of the H technique, with a smaller number of degrees of freedom.

Otherwise, for the element 7, where the solution is smooth, the convergence is shown in figure 10. The solution of P and HP technique are the same. It can be noticed that the P strategy has a greater precision, confirming the argument that the P formulation is indicated for elements with smooth solution.

The HP Formulation gives more precision with less collocation points both for smooth or singular solutions.

Example 2 – 2D Stress Analysis with Several Singularities

For example 2, the structure described in figure 11 will be submitted by the adaptive analysis.

The solution in terms of displacement and surface forces is indicated in figure 12.

As it can be observed, the structure has a critical point at vertex D. Also on the side FA and near points A and F there is an increase of the gradient of the support reactions, characterizing the need of a larger number of collocation points.

Considering that the initial mesh is fundamental to the speed of adaptive process convergence, after some attempts was adopted the refinement shown in figure 13 to start the process. Notice the smaller elements near the corners A, D and F.

The results for the H and P formulation are shown in figures 14 (a) and 14 (b). The convergence can be observed in figures ^{15 (a)} and ^{15 (b)} compared with the real solution. Notice the sensibility of the indicator for the singular areas.

For the HP formulation, the elements 3, 7, 8, 11 and 15 utilized the H strategy and the others the P strategy. The final discretization is shown in figure 14 (c) and the convergence is in figure 15 (c). In the last iteration a non-convergence process appear. That fact is attributed to the appearance of special singular integrals, due to the presence of collocation points very close to neighboring elements. This integrals have to be calculated with a more precision technique like the one developed by Guiggiani and Lombardi (1992) or Pessolani (1997b).

Conclusions

The Adaptive Hierarchical HP Formulation evidenced the applicability of adaptive programs that reach a high degree of accuracy with a smaller number of collocation points. Using it, the engineer has the analysis simplified. Those programs can be executed in personal computers of small capacity, since it requests less memory and computational effort than the conventional analysis.

The hierarchical formulation was shown as a rational way to do the adaptive analysis, optimizing the use of the machine, reaching the same degree of accuracy that can be reached using the classical formulation, in less time.

The HP formulation proved to be efficient and complete, combining the advantage of H and P Techniques, reaching the best solution with less collocation points, even in the presence of singularities.

The integration process, however, should be reviewed, due to the appearance of Quasi-Singular Integrals during the refinement. It is also necessary to extend this formulation for Fleury's studies (1993), Guiggiani (1995) and Guimarães(1994).

The Parreira error indicator and estimator worked well, capturing the regions where it was over the tolerance, increasing the number of collocation points or the order of the interpolation function. The calculation of the error has the advantage of using terms that are calculated in the hierarchical formulation.

That Formulation has to be extended to other fields, like Potential, Sub-Regions and other types of problems.

Manuscript received: July 2000, Technical Editor: José Roberto de França Arruda.

• ALARCÓN, E.; REVERTER, A. (1986), P-Adaptive Boundary Elements, Int. J. Num. Meth. Engng, 23, 801-829.
• BABUSKA, I.; GUO, B., (1986), The H-P Version of the Finite Element Method. Part I: The Basic Aproximation Results. Part 2: General Results and Applications, Comp. Mech., 1, 21-41.
• BREBBIA, C.A. (1978), The Boundary Element Method for Engineers, Pentech Press, London.
• CERROLAZA, m.; GONG-LERA, M.S.; ALARCÓN, E., (1988), Elastostatics P-Adaptive B.E. for Micros, Software for Engng Workstation, January, vol. 4.
• FLEURY JR. P., (1993), A Study on the Hyper-singular Formulation for The Boundary Element Method for Potential problems. Tangential Derivative Calculation, Ph. D. Qualifying Dissertation COPPE/UFRJ, Rio de Janeiro, Brasil
• GUIGGIANI M.; LOMBARDI F (1992), Self-Adaptive Boundary Elements with h-Hierarchical Shape Functions, Advances in Engng. Software, pags 269 a 278.
• GUIGGIANI M.(1992); A General Algorithm for the Numerical Solution of Hyper-singular Boundary Integral Equations, Journal of Applied Mechanics, Transactions of the ASME, vol. 59, 604-614.
• GUIMARÃES, S. and TELLES, J.C.F. (1994); On hyper-singular Boundary-element Formulation for Fracture-mechanics Applications, Engineering Analysis with \boundary Elements, 13, 353-363.
• PARREIRA, P.G., (1987), Error Analysis in the Boundary Element Method in Elasticity, (in portuguese) Master Thesis, Univ. Technique of Lisbon, Portugal.
• PARREIRA P.; DONG Y.F (1992), Adaptive Hierarchical Boundary Elements, Advances in Engng. Software, pags 249 to 259.
• PEANO, A.G. (1976), Hierarchics of Conforming Finite Elements goes Elasticity and Plate Bending, Comput it Planes. and Maths. with Appl., 2, No. 3-4.
• PESSOLANI, R.B.V. (1995), An HP Adaptive Hierachical Formulation for Elasticity in the Boundary Element Methods, (in portuguese) ,Ph.D. Thesis, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil, 195pp.
• PESSOLANI, R.B.V. (1997), An Adaptive Procedure for Geometry Defined by B-Splines applied in Boundary Element Method, (in Portuguese) RBCM, Revista Brasileira de Ciências Mecânicas, Vol XIX, No. 4.
• PESSOLANI, R.B.V. (1997), Semi-Analytical Singular Integration for the Boundary Element Method In Elasticity, Proc. of the Brazilian Congress of Mechanical Engineering, COBEM-97, Baurú, São Paulo.
• PESSOLANI, R.B.V. (1999), Utilization of Iterative Solvers in Adaptive Hierarchical programs for the Boundary Element Method, (in Portuguese) Proc. of the Brazilian Congress of Mechanical Engineering, COBEM-99, Aguas de Lindóia, São Paulo.
• RANK, E. (1989), Adaptive h -, p - and h-p Versions goes Integral Boundary Elements Methods, Int. J. In a. Meth. in Engng., vol 28, 1335-1389.
• RENCIS, J.J.; MULLEN R.L. (1986), Solution of Elasticity Problems by Self Adaptive Mesh Refinement Technique in the Boundary Element Computation, Int. J. In a. Meth in Engng, vol 23, 1509-1527.
• ZIENKIEWICZ, O.C.; GAGO, J.P. OF S.R.; AND KELLY, D.W. (1983), The Hierarchical Concept in Finite Element Analisys, Comp. and Struct., 16, In the 1-4.
• ZIENKIEWICZ, O.C.; ZHU, J.Z.; GONG, N.G. (1989), Effective and Practical H-P Version Adaptive Analysis Procedures goes the FEM, Int. J. In a Meth in Engng, vol 28, 879-891.

Publication Dates

• Publication in this collection
  18 Sept 2002
• Date of issue
  Mar 2002

History

• Received
  July 2000