a formulation of optimal mesh for problems with

COMPUTATIONAL MECHANICS New Trends and Applications S. Idelsohn, E. Oñate and E. Dvorkin (Eds.) CIMNE, Barcelona, Spain 1998

A FORMULATION OF OPTIMAL MESH FOR PROBLEMS WITH MULTIPLE ANALYSIS CONDITIONS F. Javier Fuenmayor†, Jorge L. Restrepo‡, Juanjo J. Ródenas†, and José E. Tarancón† †

Departamento de Ingeniería Mecánica y de Materiales Universidad Politécnica de Valencia Camino de Vera, 14, 46022 Valencia, Spain e-mail: [email protected] ‡

Departamento de Ingeniería Mecánica Universidad EAFIT Carrera 49 No 7 Sur -50 Apartado Aéreo 3300 Medellín, Colombia e-mail: [email protected]

Key words: finite element method, optimal mesh, h-refinement. Abstract. Mesh adaptivity in Finite Element has been extensively studied during the last years. There are different approaches for the definition of the optimal mesh for problems with a single analysis condition or load case (LC). Traditionally the optimal mesh has been defined as the one in which the global error is distributed equally in each element. It has been recently proved that the minimization of the number of elements criterion is equivalent to the former one. h-adaptivity criteria are well defined when a single LC is considered, however, little research has been carried out relating to the problem of optimum mesh generation with multiple LC. Due to the fact that the error distribution depends on each LC, a mesh which is optimum for one of the LC will not generally be optimum for the rest of them. Therefore, the objective of this paper is to show a new general h-adaptivity method for multiple LC problems. The method is based on minimising the number of elements so that the global errors obtained are smaller or equal than the specified values. It should be noted that by means of the proposed method it is possible to obtain the simultaneous refinement for all the different LC as a combination of the necessary refinements for each of the LC individually considered. This allow us to evaluate the participation of each individual LC over the global refinement. The proposed method has been validated by applying the method to different test problems.

1

F. Javier Fuenmayor, Jorge L. Restrepo, Juan J. Ródenas, and José E. Tarancón

1 INTRODUCTION The aim of an h-adaptive process is to improve the finite element solution with the sole modification of the element size. In order to carry out this process, an estimation of the discretization error for every element is needed. In addition, an optimum mesh criterion must predict the new element size in the refined mesh to reduce the error to a previously specified value. Nowadays, a common practice is to consider that the optimum mesh is the one which gives the specified error with either the minimum number of degrees of freedom or minimum number of elements. The definition of optimum meshes for one LC has been fairly well studied during the last years. Fuenmayor and Oliver1 present five strategies for adaptive procedures, together with the equations employed to evaluate the element size of the new mesh and their features derived from numerical experience. However, not all of these strategies can be applied to problems with several LC. Therefore the definition of optimum the mesh in problems with more than one LC has not been sufficiently treated. The purpose of this work is to extend one of the strategies for adaptive processes defined by Fuenmayor and Oliver1, to problems with several LC. The procedure consists of a redefinition of the element size to achieve the prescribed error (or a lower one) for every LC, minimizing the number of degrees of freedom in the new mesh. Firstly, the method used to estimate the discretization error is reviewed and the basic equations related to global and local (i.e. at the element level) convergence are presented. Secondly, the procedure to get the optimum mesh for only one LC and its extension for several LC are exposed. Finally, the reliability of the proposed method was evaluated by its application to a two-dimensional problem. 2 ESTIMATION OF THE DISCRETIZATION ERROR The discretization error in the finite element solution is often quantified on the basis of the energy norm of the error in displacements. For the i-th LC, this norm can be expressed as: e ex( i )

2

= u ex( i ) − u fe( i )

2

=

∫ (σ

ex( i )

− σ fe ( i )

Ω

)

T

(

)

D−1 σ ex( i ) − σ fe( i ) dΩ (1)

where uex(i), ufe(i), σ ex(i) and σ fe(i) are the exact and finite element displacement fields and the exact and finite element stress fields corresponding to the i-th LC, D the stress-strain relation matrix, and Ω the domain on which the problem is defined. The exact relative error is defined as the ratio of the energy norm of the error to that of the exact solution: ηex( i ) =

2

e ex( i ) u ex( i )

(2)


The underlying principle behind the Zienkiewicz-Zhu error estimator is to use a smoothed stress field σ *(i) derived by post-processing of the finite element stress field, rather than the exact stress field σ ex(i), which is unknown. Thus, the discretization error estimator, for the i-th LC, can be defined as: 2

e es( i )

=

∫ (σ

*( i )

)

T

(

)

− σ fe( i ) D−1 σ *( i ) − σ fe ( i ) dΩ

Ω

(3)

The error estimator can also be defined as a function of the error estimation at element level: e es( i )

2

=

Ne

∑e

2 (e) es ( i )

=

e =1

∑ ∫(σ Ne

e =1

*( i )

− σ fe ( i )

Ωe

)

T

(

)

D −1 σ *( i ) − σ fe( i ) dΩ

(4)

where Ne represents the number of elements in the mesh and Ω e is the domain corresponding to element e. The estimation of the relative error can be calculated as: ηes( i ) =

e es( i ) u es( i )

e es( i )

= u fe( i )

2

+ e es( i )

2

(5)

The reliability of the error estimator can be assessed via the effectivity index, defined as the ratio of the estimated error to the exact error: e es( i )

θ ( i) =

(6)

e ex( i )

The error estimator is said to be asymptotically correct if the effectivity index tends to unity as the discretization error tends to zero. Obviously, the error estimator and its effectivity index will be conditioned by the way in which the smoothed stress field (σ σ *(i)) is obtained. The SPR2 technique was used, which is accurate for any type of elements3. 3 CONVERGENCE RELATED EQUATIONS For h-adaptive refinement, when the sequence of meshes is obtained through uniform refinement, the energy norm of the exact error is bounded (for bi-dimensional problems) by: e ex( i ) ≤ C1( i ) h c ≈ C2 ( i ) N

1 − c 2

(7)

where N is the number of degrees of freedom, h is the element size and C1(i) and C2(i) are positive constants independent of h or N. The exponent c is expressed as c = min(p,λ), where λ represents the ‘strength’of the singularities and p is the polynomial degree of displacements interpolation functions

3


being used. The exponent of N is called the asymptotic rate of convergence and depends on the finite element mesh and the smoothness of the exact solution. When the sequence of meshes is designed in such a way that the absolute error is approximately equally distributed over all the elements, the maximum asymptotic rate of convergence for h-adaptive refinement is obtained4, and the exact error is bounded by: e ex( i ) ≤ C 2 ( i ) N

1 − p 2

(8)

These equations can be evaluated locally, i.e. considering the domain defined by one element in the previous mesh and the elements contained in this domain in the new mesh. In general terms, the value of c depends on the element under consideration. If the exact solution inside the element is smooth then c = p, but if it is not smooth (element close to a singularity) then the value of c will tend to λ. In order to simplify the equations, it will be assumed that c = p. Assuming it is possible to write the previous expressions as a function of the estimated errors instead of the exact errors, the ratio of the error in the new mesh to the error in the previous mesh can be evaluated as: e (ese()i ) e (ese()pi ) where e (ese()i )

p

 h (pe()i )  ≈  ( e )p  = r((ie))  h n ( i )  c

p

n

( )

c

(9)

is the estimated error for element e in the previous mesh, having a size h(pe()i ) , and (e)

e (ese()ip) is the estimated error in all elements in the new mesh (with size h n ( i p) ), contained in the n

element e of the previous mesh. This equation defines the local refinement ratio r((ie) ) of each LC. The total error in the new mesh can be evaluated as a function of the local error estimation in the previous mesh and the local refinement ratio r((ie) ) by taking into account the former equation: e es( i )

2 n

=

Nn

∑ e=1

e (ese()i )

2 n

Np

= ∑ e (ese )( pi ) e =1

2 n

=

Np

∑  e e=1

2 (e) es( i ) p

(r ) (e) (i )

− 2c

 

(10)

where Nn and Np are the number of elements in the new and the previous meshes respectively. On the other hand, the number of elements in the new mesh contained in the element e, corresponding to the previous mesh, can be estimated as:

( )

N (ne()i p) ≈ r((ie))

2

(11)

Therefore, the total number of elements in the new mesh can be estimated as: Np

( )

N n ( i ) ≈ ∑ r((ie)) e =1

4

2

(12)


4 OPTIMAL MESH FOR A SINGLE LOAD CASE In order to define the optimal mesh the scheme described by Ladeveze et al.5,6,7 will be followed. In his works, the optimal mesh is defined as that which obtains the required error with a minimum number of elements. It has been recently proved that this criterion is equivalent to the traditional one8 in which the discretization error is equally distributed in each element of the new mesh1,9,10. Defining e ( i )

d

as the desired absolute error for the new mesh and for each LC, the problem of

finding the local size ratio would be:

∑ (r Np

Minimize: Nn ( i ) =

e =1

Subject to: ∑  e (ese()i ) e=1 Np

(e) ( i)

2 p

)

2

(r ) (e) (i )

− 2c

 − e (i ) 

(13) 2 d

=0

The solution to this problem can be evaluated by means of a multiplier α (i) for each of the LC so that the restriction is incorporated to the function to be minimized. By solving this problem, the size ratio which leads to the optimal mesh is defined as:

(r ) (e) (i )

2 ( c +1)

2

= c α( i ) e (ese )( i )

p

(14)

where:

c α( i )

 Np ( e ) c2+1   ∑ e es( i ) p  e=1  =  e 2  ( i) d    

c+ 1 c

(15)

The estimated errors at element level are already known. Therefore, and for each LC, the last expressions allow us to calculate how each of the elements in the previous mesh has to be refined in order to obtain the desired error in the new mesh. 5 OPTIMAL MESH FOR MULTIPLE LOAD CASES When the optimal mesh for multiple LC is to be found , the problem to be solved is that of minimizing the number of elements in the new mesh under the restrictions of obtaining the specified error for each of the LC. Hence, the optimal mesh is found by solving the following equations:

5


∑ (r Np

Minimize: N n =

(e)

e =1


)

2 p

2

(r ) (e)

− 2c

 − e (i ) 

(16) 2 d

=0

i = 1, .., N LC

where NLC is the number of LC under consideration. As in the previous case, the problem can be solved by considering a multiplier β (i) for each of the different LC. Once the problem has been solved, the size ratio is given by:

(r( e ) )

2 ( c+1)

=

∑ c β N LC

i =1

(i )

2 e (ese()i )  p

(17)

where the values of β (i) are obtained by solving the following system of non-linear equations: Np

∑ e =1

c  −  N LC c +1   2  (e) 2   (e) c β( j ) ees ( j )   e es( i ) p  ∑  = e ( i) p j =1  

2

i = 1, .., NLC

d

(18)

Let us define the global refinement contribution factor (ξ (i)) for the i-th LC as: ξ( i ) =

β( i ) α( i )

(19)

Taking into account equation (14) the local refinement ratio can be evaluated as:

(r ) (e)

2 ( c+ 1 )

=

N LC

∑ ξ i =1

( i)

(r ) (e) (i )

2 ( c+ 1 )

 

(20)

where the values for the global refinement contribution factors for each LC are obtained by solving the following system of non-linear equations: Np

∑ e =1

 N  ( e ) 2  LC e ξ( j ) r((je))  es( i ) p  ∑ j =1 

( )

2 ( c + 1)

  

−

c c+1

   = e( i ) 

2 d

i = 1, .., N LC

(21)

Therefore, according to equation (20), when multiple LC are considered, the ‘global’size ratio can be evaluated as a function of the size ratios evaluated for each of the different LC weighted by the global refinement contribution factors ξ (i). The values of the global refinement contribution factors ξ (i) can be positive, negative or zero. ξ (i) = 0 means that the i-th LC does not participate at all in the refinement, i.e. the refinement imposed by some other LC adequately refines the LC under consideration. A positive value of ξ (i) means that, in order to obtain the specified error for the i-th LC, the refinement defined by the LC under consideration has to be taken into account when calculating the global refinement. A negative value

6


for ξ (i) means that other LC are already refining the mesh in such a way that the error for the i-th LC would be less than its specified value. Therefore, the i-th LC would tend to un-refine the mesh so that the error takes the specified value. The last situation is shown in Fig. 1. Two different LC are taken into account. Both cases are analysed by a previous mesh with Np degrees of freedom, obtaining relative errors of ηp(1) and ηp(2) respectively. It is desired to obtain a new mesh with relative errors of ηd(1) and ηd(2). By refining considering only the first LC, the error reduction would follow the lines labelled as ‘LC1 ref.’. NLC1 degrees of freedom would be necessary to obtain an error equal or smaller than that specified for the first LC. With NLC1 degrees of freedom and for the second LC, the error in the new mesh would be smaller than that initially specified (ηn(2) < ηd(2)). Notwithstanding, by considering both LC, the errors in the new mesh would be equal to those specified at the expense of using a larger number of degrees of freedom (Nglobal > NLC1), i.e. un-refining according to the second LC. In this situation the value of the global refinement contribution factor for the second LC would be negative (ξ (2) < 0). As a consequence, in this example, the optimal mesh should be evaluated only through the first LC refinement.

Figure 1. Solution convergence as a function of the LC which controls the h-adaptive refinement.

Therefore the initial problem definition should be modified in order to consider inequality constraints in the optimization problem:

7


∑ (r Np

Minimize: N n =

(e)

e =1


2 p

)

2

(r ) (e)

− 2c

 − e (i ) 

(22) 2 d

≤0

i = 1, .., N LC

The optimum mesh could be found either by solving this optimization problem or by consecutively excluding from equation (22) those LC with negative values of the global refinement contribution factor. 6 NUMERICAL VERIFICATION The proposed scheme for the definition of the h-adaptive process has been applied for the resolution of a number of test problems. In this work we will show the results for only one example, as other examples would not show significative additional result. Figure 2 shows the geometry of the proposed problem, material data and LC under consideration. Each of the different LC is defined as a constant pressure applied over different contours. The problem has been solved using quadratic triangular elements. The original mesh is show in Fig 2.b. Four different h-adaptive processes have been applied. In the first three processes only one of the LC is used to guide the h-adaptive refinement. In the 4th process all of the LC are simultaneously considered in the h-adaptive process. The objective is that of obtaining a 2% discretization error by means of 5 consecutive refinements. Fig. 3 shows the different mesh sequence obtained for each of the defined processes. The element size distribution required by each of the LC can be observed in this figure. In this work the new meshes have been generated by using the Delaunay’s triangularization method, as described in11.

Figure 2.a. Problem definition.

Fig. 2.b.- Initial mesh.

8


Load case 1 refinement



Global refinement

Figure 3. Mesh sequence for each of the h-adaptive processes.

9



100

1.5


LC 2 LC 3

Efectivity Index

Estimated Error (%)

LC 1

10

1000

LC 1 0.5

0 100

1 100

1

10000


1.5

100

Efectivity Index

Estimated Error (%)

LC 1 LC 2 LC 3

10

1000

LC 1 LC 2 LC 3

1000

10000

d.o.f.



1.5

LC 1 LC 2 LC 3

Efectivity Index

Estimated Error (%)


0.5

0 100

10000

10

1

LC 1 LC 2 LC 3

0.5

1

0 100

1000

10000

100

d.o.f.

1000

Global refinement

Global refinement

1.5

Efectivity Index

LC 1 LC 2 LC 3 10

1

LC 1 LC 2 LC 3

0.5

1 100

10000

d.o.f.

100

Estimated Error (%)

10000

1

d.o.f.

100

1000

d.o.f.

d.o.f.

1 100

LC 2 LC 3

0 1000

10000

d.o.f.

100

1000

10000

d.o.f.

Figure 4. Convergence of error estimation for the different h-adaptive processes.

10

Figure 5. Effectivity index of the estimated discretization error for the different hadaptive processes


Fig. 4 shows the estimated discretization error vs. the total number of degrees of freedom for each of the four h-adaptive mesh control strategies. If mesh size is only controlled by LC 1, this LC will reach a 2% error while the other two LC wont. If mesh size is only controlled by LC 2 then both LC 2 and 3 will reach error values smaller than 2%; the error for LC 1 will be slightly higher than the required value. When the mesh refinement is only controlled by LC 3, the errors obtained for LC 1 and 2 are higher than 2%. As a consequence of these results it is clear that LC 3 should not be considered in the global h-adaptive refinement as the rest of them will already refine the mesh sufficiently for this LC. When a global h-adaptive refinement is carried out by using the proposed calculation scheme, LC 1 and 2 reach the desired error while for LC 3 the error obtained is smaller than 2%. This can be explained taking into account the participation of each of the LC over the global refinement, as shown in Table I. The contribution factor for LC 3 would always be negative, except for the case of mesh 1, therefore this LC has been excluded from the control of the h-adaptive global process.

LC1 LC2 LC3

mesh2 0.447 0.677 0.124

mesh3 0.357 0.807 0.000

mesh4 0.087 0.948 0.000

mesh5 0.255 0.837 0.000

mesh6 0.361 0.766 0.000

Table I: Global refinement contribution factor for each LC.

Fig. 5 shows the estimator error effectivity index for the mesh sequence of each h-adaptive process. The solution obtained with a very refined mesh was taken as the exact solution. For this mesh the estimated errors were, at least, one order of magnitude lower than the errors for the meshes under study. Fig 5 shows that the estimator effectivity is better for those LC with the highest participation in the refinement. Therefore, error estimation for those LC that do not participate in the refinement could be unreliable. Notwithstanding, in general terms, the effectivity values obtained are acceptable, as shown in the present example. 7 CONCLUSIONS In this work, a method to define an h-adaptive procedure for problems with several LC is proposed. The purpose of the procedure is to minimize the number of elements in the refined mesh, with a discretization error less than or equal to a previously specified value for each LC. This scheme allows to define the global refinement as a function of the refinements associated to each LC. Accordingly, the global refinement contribution factor for every LC is defined and a simple method to exclude those LC which are enough refined by other LC is also considered. The procedure is computationally affordable, since the system of equations required to solve the problem has as many unknowns as the number of LC considered. The proposed method has been numerically validated through different examples.

11


REFERENCES [1]

F. J. Fuenmayor and J. L. Oliver, ‘Criteria to achieve nearly optimal meshes in the h-adaptive Int. J. Numer. Meth. Eng, 39, 4039-4061 (1996). [2] O.C. Zienkiewicz and J.Z. Zhu, 'The Superconvergent Patch Recovery (SPR) and Adaptive Finite Element Refinement' Int. J. Numer. Meth. Eng, 101, 207-224, (1993). [3] I. Babuska, T. Strouboulis, C.S. Upadhyay, S.K. Gangaraj and K Copps, ´Validation of a posteriori error estimators by numerical approach´Int. J. Numer. Meth. Eng, 37, 1073-1123 (1994). [4] I. Babuska and B. Szabó, ‘On the rates of convergence of the finite element method’, Int. J. Num. Meth. Engng, 18, 323-341 (1982). [5] P. Ladeveze and D. Leguillon, ‘Error estimate procedure in the finite element method and SIAM J. Numer. Anal., 20(3), 485-509 (1983). [6] P. Ladeveze, P. Marin and J. P. Pelle, ‘Accuracy and optimal meshes in finite element computation for nearly incompressible materials’, Comput. Method Appl. M. Eng, 94, 303315 (1992). [7] P. Coorevits, P. Ladeveze and J. P. Pelle, ‘An automatic procedure with control of accuracy for finite element analysis in 2D elasticity’, Comput. Method Appl. M. Eng., 121, 91-120 (1995). [8] O. C. Zienkiewicz and J. Z. Zhu, ‘A simple error estimator and adaptive procedure for practical engineering analysis’, Int. J. Numer. Meth. Eng, 24, 337-357 (1987). [9] L. Y. Li, P. Bettes, J. Bull, T. Bomd and Y. Applegarth, ‘Theoretical formulations for adaptative finite element computations’, Commun. Numer. Meth. En., 11, 911-915 (1995). [10] L. Y. Li and P. Bettes, ‘Notes on mesh optimal criteria in adaptive finite element computations’, Commun. Numer. Meth. En., 39, 4039-4061 (1996). [11] F. J. Fuenmayor and J. L. Oliver, A ‘ plicación de la Triangulación de Delaunay en procesos adaptativos de elementos finitos’, Memoria I Congreso de Métodos Numéricos en Ingeniería, SEMNI, 698-705 (1990).

12