- Email: [email protected]
Singular shape functions for an automatic mesh generation around singularities
Engineering
Fracture
Mechanics
Printed in Great Britain.
Vol. 36, No. 3, pp. 451-458,
$3.00+ 0.00 0013-7944/90 0 1990 Pergamon Press plc.
1990
SINGULAR SHAPE FUNCTIONS FOR AN AUTOMATIC MESH GENERATION AROUND SINGULARITIES G. TSAMASPHYROS
and M. PAPAIOANNOU
Section of Mechanics, Department of Engineering Sciences, The National Technical University of Athens, 5, Heroes of Polytechnion, GR 157 74, Greece Abstract-A
proper gradual condensation of the finite element mesh around singularities is obtained by using a mapping technique. This technique is based on the construction of special singular shape functions. An advantage of the method is that it is entirely automatic and has not the drawbacks of a similar previous method.
INTRODUCTION LAST YEAR’S finite element method became one of the most powerful procedures for a numerical approximate solution of boundary value problems. One of the inherent difficulties of the method is the need of a great number of input data related to the topology of the mesh. In the case of problems without singularities the problem has already been solved by the use of various mesh generation schemes[l+. But the standard automatic mesh generation schemes cannot confront problems where a stress singularity appears. In this case a concentration of elements is needed around the singular point even in the case where special elements are used[5-81 (see also [lO-141). Thus the finite element mesh must present a gradual transition from small elements, around the singular points, to large elements, in the areas where regular stresses are expected. This is a very difficult task since the elements must be as regular as possible[9] and the optimum size or the rate of the gradual change of the size of the elements are not a priori known. In the present paper a method for an automatic generation of optimum mesh is presented. The method allows to confront singularities of any order and it is based on the use of special shape functions which maps a regular mesh to a mesh with an optimum concentration of elements. In fact the basic idea of the present method is the same with the previous ones[13-141, i.e. we use an approximate transformation which removes the singularity from the solution of the problem in the transformed domain.
MAPPING PROCEDURE Let K?: (Fig. la) equal a finite element. This element using isoparametric shape functions Mj(%, H) is mapped to a regular shape SQ(Fig. lb) in the parametric plane. Any field of parameters L can be described in the Q: domain by: L = f M,(E, H)L, j= I
where the (Lj}/“_I are the values of the parameter L at the m nodes of the element. In the case of linear elements (triangular or quadrilateral) the quantity: pi=l--zMi
(2)
represents the non dimensional radius from node i. This assertion remains true in the more general case of the element of the Fig. l(b) only for the corner node 1 since along the edges l-m and 1-2, M, behaves linearly. In the general case the quantity pi can be considered in the vicinity of the node i, as a measure of the non dimensional radius. Consider now a plane body (Fig. 2) having a singularity at A. We subdivide this body to a number of “super elements” 0; and let node 1 of each super element coincide with the singular point A. If we denote by r the radius from A, then the stresses and strains behave like O(r”-‘). Obviously the same (2 - 1) singularity will exist in each super element and the same singularity 451
452
G. TSAMASPHYROS
and M. PAPAIOANNOU
3 H
n-l-1
2
LE
5 mcl
1
m
(b) Parametric plane or
(a) Initial plane
mapped plane Fig. 1. (a) Initial super element Q:, (b) element q into the parametric plane obtained from Q using singular mapping.
[O(p”-I)] will exist after the mapping into the parametric plane. (This will be explained by the linearity of the mapping.) The above singular behavior is due to terms like p which exist in the displacement field. If a change of variable, like: Pi= Pt
(3)
is performed, then the stress field in the transformed p-plane would be regular around the node i. That means that the strain and the stress field and consequently the displacements field is reproducible by polynomials in the new plane. In this way we can use conventional polynomial shape functions {(Mi):_ 1 in order to describe the displacements u = (uX, uY} in the transformed plane: u=Mq’
(4)
q’ being the vector of nodal displacements of the super element Q. Another immediate consequence is that no special mesh is needed in the transformed CL q )-plane. Thus a special mesh can be constructed from a regular one in the mapped plane by using the inverse transformation. This fundamental observation has already been used[13-141 by taking the conformal mapping as the transformation which could perform a change of variables like the one described by the relation (3). The essential disadvantage of this transformation is that it distorts the boundaries of each super-element &2:. That means the boundaries of Qt (which is the image of a:) cannot be
Fig. 2. Plane body with singularity at A. Subdivision into “super elements” q.
Singular shape functions
453
represented by polynomials. This creates some difficulties in the exact representation of the boundaries of the initial domain. Methods to overcome these difficulties are described in the above papers[l3-141, but a method without this drawback is desirable. In the method which will be described here the initial super element Pz will be mapped onto the unit square (Fig. lb). This mapping is based on the previous described possibility to represent the normalized radius from node 1 through the shape function, see relation (2). To this end we will represent by: JJ,=l-M,
(5)
the non dimensional radius from node 1 (in the transformed plane). The shape functions N which could describe the singular behavior of strains in the initial plane must be singular. Nevertheless we consider again the quantity: p, = 1 -N,
(6)
as a measure of the normalized radius from node 1. Using (3) we have: (1 -M,)
= (1 - iv,)“.
(7)
That means: Ni = 1 - (1 -M,)“? Now taking into consideration
(8)
that Ni must satisfy the relations: Ni(
qj)
=
$,
i, j = l(l)m
6,;
N(419V)=
(9)
l
(10)
we can easily construct the other shape functions Ni (i = 2(l)m). To this end we take into account
that M satisfies relations similar to the relations (9) and (10). Nj=(l
5 MA. j=2(l)m
-N,)Mj/
(11)
k==2
These relations (8) and (11) have been used locally by Okabe et al.[ 15-l 61 in order to construct singular elements. The possibility to map a part of the whole body has not been explored. The previous idea can be generalized by using appropriate weight coefficients aj. Thus Nj=(l
-N,)ajMj/
2
(akitfk).
j =
2(l)m
(12)
k=2
The mesh obtained by the singular mapping as well as the influence of parameters aj in the distortion of the mesh is depicted in Fig. 3 (cases a-c). The above analysis is based on the representation of the dimensional radius in terms of 1 - M, . But the dimensional radius pl can be represented by the relation: pI = maxE
+ 1)/2, (tt +
WI.
(13)
If the above relation is substituted into relation (3) expressions similar to relations (8) and (12) could be obtained. Another possibility is to use a combination of (2) and (13) in order to represent the dimensional radius; for example: PI = PZP,p,8*
a+/3=1
(14)
where p,. is given by (2) and pls is given by (13). In the case where the body has many singular points it is possible to use successively the previous transformations for each singular point in order to obtain the final mapping. Consecutive transformations can be used in the case where Mi is not a linear function and the edges l-m and l-2 are curvilinear. In that case using the parametric transformation (l), the initial super element Pz is mapped to a square in the parametric plane. We map afterwards the square using the singular transformation (8) and (12) but supposing that the shape functions Mi
G. TSAMASPHYROS
454
and M. PAPAIOANNOU
(b)
(4 Fig.
3. Singular
mesh obtained from relations (8) and (12) with: and (c) a,=a,=O.S,a, = 4 = a2 = a, = f.O,a,=OS;
(a) a,=~,
=a4= I; (b)
1.0.
(1 = l(l)m) involved in these formulas are the shape functions M; (i = 1(1)4) of the four node quad~latera1. The above examples concerned plane problems but the extension to three-dimensional problems is obvious. On the other hand we are only interested for problems where the singularity is at described points. We consider now the case where the singularity is along a particular line S of the body. In that case one or more elements into which the body is subdivided, must have one side along S. The normalized radial distance from S is p = 1 -(Nj+Nj+l)
(15)
where nodes i and (i + 1) are along S. If the transformation (3) is performed we would have: 1 - (N, + Mz) = (1 - [N, + NJ)“.
(16)
In the case of a four-node element (16) gives: 043 + A44)= (A$+ N4)1:
(17)
and a possible choice of A$ and A$ will be: N, = Mi(M3 + M#‘“-1;
i = 3,4
(18)
Singular shape functions
455
and consequently: N, = M,(l - (M3 + M,))““-‘.
i = 1,2
(19)
An illustration of the distortion of the regular mesh produced by the above formulas (18) and (19) (with 1 = l/2) is given in Fig. 4(a). We use now the previous transformation (15) not for a distance from an arbitrary line S but from the Ox and 0~ axis. In this way we obtain the normalized coordinates x and y which are given by: x=1-(N,+N,)=N*+N,,
JJ=l-(N,+N,)=N,+N,
(20)
in the case of a four node quadrilateral. Using the transformation (17) for x and y we have: &+M3=(Nz+N,)i,
M,+M,=(N,+N,)”
(21)
01
N2 + N3 = (M2 + h4p,
N, + N2 = (M, +
M*)‘?
(22)
Consequently, N, - NS = (M, - M,)[(M, + M,)“‘l-
(M3 + M,)“V(M,
- M,)
(23)
and thus we can take:
w’il/(~3 + 4) N, = A43 [(M,+ A&)“”- w3 + W’“lI(~3 + W) N, = M, [(M, + &)“A - (MS +
(24)
N2 = (M2 + &)I’*: - N,
I
N,=l-N,-Nz-N,.
In the special case of a crack (1 = l/2) the previous formulas can be simplified. An illustration of this transformation is given in Fig. 4(b). The above examples concerned plane problems but the extension to three-dimensional problems is obvious. NUMERICAL
EXAMPLES
In order to check the efficiency of the proposed method a rectangular sheet of width 26, height 2h, containing a straight crack of length 2u was considered. The plate was subjected to uniform
(a)
(b)
Fig. 4. Singular mesh obtained (a) from formulas (18) and (18) or (b) from formulas (24) with I = l/2.
G. TSAMASPHYROS
456
and M. PAPAIOANNOU
Table 1. Percentage error of the stress intensity factor for a central crack with a/b = 0.5 using the singular mappin,g of Fig. 3(a) and varying number of nodes hr, of coefficients C, and collocation points n, of the formula (25) Total number of nodes (N)
Number of nodes along the crack (n)
21 65 121 133 173 181 199 245 309
3 5 9 7 9 11 9 11 11
Percentage difference between exact value(21] and
formula (25) 13.29 -3.98
- 12.08 -5.46 -8.39 - 12.43 -7.03 -9.33 -7.06
special collocation formula[ 141 -2.82 -0.85 - 0.05 0.10 0.30 0.25 0.38 0.47 0.57
tension a, along the parallel to the crack edges. Three values (1.0,0.7,0.4) of the ratio h/b and various values of the ratio a/b were considered in order to study the behavior of the procedure when the loaded boundary comes near the crack. In order to reduce the computational effort the symmetry was exploited, thus a quarter of the sheet is descritized by the above described method and a singular finite element mesh is obtained. After that a classical finite element program can be used. Since, any of our singular finite element meshes, has curvilinear elements, isoparametric elements must be utilized. We have selected quad~lateral elements and especially the 8-node ones. Since, in the case of cracked body the order of sin~a~ty 1 is equal to (l/2) we can easily see from formula (7) that the crack tip element is automatically reduced to a quarter point singular element[ 17j. The stress intensity is calculated by using the special collocation formula described in [14]. The stress intensity factor can also be obtained by the formula: K, - ZK,,=
E,/%(4AU,-AU,)
(25) 4&K=% where AU, (or AU, ) is the difference of complex displacements of the opposite points 2 and 2’ (or 1 and l’), with 2’ (or 1’) being nodes of the lower lip of the crack and 2 (or 1) nodes of the upper lip of the crack. The nodes 1 and 2 or the nodes 1’ and 2’ belongs to the same crack tip element. Node 0 coincides with the crack tip. This formula is given in [18] and inde~nden~y in 119,201. First the influence of the number N of nodes is studied. To this end the case a/b = 0.5 and h/b = 1.0 is considered with a varying number of total nodes iV and a varying number of nodes m along the crack lips (m > n), The results figure in Table 1 where the first column corresponds to the number of nodes N, the second one to the number n of nodes along the crack lips, The third Table 2. Percentage error of the stress intensity factor for a central crack with a/b = 0.5 using an arbitrary mesh, and a varying number of nodes N, of coefficients C,, and collocation points n Total number of nodes (N) 53 ;; 121 133 199 245 309
Number of nodes along the crack (u)
Percentage difference between exact vaIue(21] and formula (25)
special collocation fO~~~I4]
7 7 7 9 7 9 11 11
81.81 106.80 124.99 142.21 148.40 181.41 199.88 218.56
29.25 34.93 37.88 41.13 41.17 45.20 47.12 48.35
451
Singular shape functions Table 3. Percentage error of the stress intensity factor for a central crack with h/b = 1.Oand various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 200)
alb 0.3 0.4 0.5 0.6 0.7 0.8
Number of nodes along the crack (n) 1 9 9 11 11 11
Percentage difference between exact value[ll] and
formula (25) - 8.94 - 12.44 -7.03 -5.14 -5.33 -2.60
special collocation formula[ 141 -0.86 -0.41 0.38 -0.85 - 10.12 0.91
column gives the difference between exact value[21] of the stress intensity factor and the one obtained from the formula (13). As we can see the results of the collocation procedure are very satisfactory even in the case of a very coarse mesh. On the contrary the accuracy of the formula (13) is not satisfactory and this contradicts with the claims of Ingraffea et a1.[20]. In order to have an idea of the improvement obtained by the use of the above described singular finite element mesh (see formula 7 and 10) a regular finite element mesh is constructed with the same number of nodes. The results of the above mesh are shown in Table 2. As we can see the improvement is spectacular. For h/b = 1.O,0.7,0.4 and a varying crack length the difference of the calculated stress intensity factor from the exact one is given in Tables 3,4,5, respectively. The results of Table 3 are obtained from finite element patterns with approximately N = 200 nodes. In the case of Tables 4 and 5 the number of nodes was N = 150 and N = 130, respectively.
Table 4. Percentage error of the stress intensity factor for a central crack with = 0.7 and various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 150)
h/b
alb
Number of nodes along the crack (n)
Percentage difference between exact value[21] and
formula (25)
0.3 0.4 0.5 0.6 0.7
5 7 9 9 11
-1.58 -5.15 - 5.47 -4.38 -3.41
special collocation formula[ 141 -0.71 0.46 -1.38 - 1.43 -1.18
Table 5. Percentage error of the stress intensity factor for a central crack with h lb = 0.4 and various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 130)
a/b 0.3 0.4 0.5 0.6 0.1
Number of nodes along the crack (n) I 9 11 11 11
Percentage difference between exact value[21] and formula (25) -4.82 -5.65 -6.00 -5.44 -5.34
special collocation formula[l4] 1.64 -2.12 -1.79 -2.06 -2.18
458
G. TSAMASPHYROS
and M. PAPAIOANNOU
REFERENCES [I] W. P. Buell and A. A. Bush, Mesh generation-a survey. J. uppf. Mech. 332-338 (1973). [2] R. D. Shaw and R. G. Pithen, Modifications to Suhara-Fukuda method of network generation. Inr. J. namer. Merh. Engng 12, 93-99 (1973). (31 S. E. Umanskij, Triangulation algorithm and program two-dimensional domains of any shape. Srrengrh Marer. 9, 700-706 (1978). [4] I. Imafucu, Y. Kodera and M. Sayawaki, A generalized automatic mesh generation scheme for finite element method. Int. J. numer. Meth. Engng 15, 713-731 (1980). [S] 0. C. Zienkiewicz and D. W. Philips, An automatic mesh generation scheme for plane and curved surfaces by “isoparametric” co-ordinates. hr. J. numer. Me&. Engng 3, 519-528 (1971). [6] R. H. Gallagher, A Review of Finite EIemenr Techniques in Fracture Mechanics (Edited by A. R. Luxmoore and D. R. J. Owen), pp l-25. University of Wales, Swansea (1978). [7] K. D. Henshell, Crack tip fmite elements are unnecessary. Inr. J. numer. Merh. Engng 9, 495-507 (1975). (81 0. C. Zienkiewicz, The Finire Element Method, 2nd edn. McGraw-Hill, New York (1977). [9] M. Zlamal, A finite element procedure of the second order of accuracy. Numer. Murh. 14, 394Xt2 (1970). [IO] J. J. Guydish and J. F. Fleming, Optimization of the finite element mesh for the solution of fractural problems. Engng Fracture Mech. 10, 31-42 (1978). [I I] J. J. Guydish, J. F. Fleming, J. R. Pentx and C. E. Runnion, The finite element method vs the edge function method for linear fractural analysis. Engng Fracrure Mech. 13, 43-55 (1980). [12] P. R. Brown, A non-interactive method for the automatic generation of finite element meshes using the Schwartz-Christoffel transformation. Compur. iUerh. appl. Mech. Engng 25, 101-126 (1978). [ 131 G. Tsamasphyros and A. E. Giannakopoulos, The mapped elements for the solution of cracked bodies. Compur. Merh. appl. Mech. Engng 49, 331-342 (1985). [ 141 G. Tsamasphyros and M. Papaioannou, Mapping technique for the automatic mesh generation around singularities. Compur. Srrucr. 29, 815-823 (1988). [IS] Y. Yamada, Y. Ezawa, I. Nishiquchi and M. Okabe, Reconsiderations on singularity or crack tip elements. Inr. J. numer. Merh. Engng 14, 1525-1544 (1979). [16] M. Okabe, Fundamental theory of the semi-radial singularity mapping with applications to fracture mechanics. Compur. Merh. appl. Mech. Engng 26, 53-73 (1981). [17] R. S. Barsoum, Application of quatratic isoparametric finite elements in linear mechanics. Inr. J. Fracrure 10,603-605 (1974). [18] G. Tsamasphyros and A. E. Giannakopoulos, Automatic optimum mesh around singularities using conformal mapping. Engng Fracture Mech. 23, 507-520 (1986). [19] C. F. Shih, H. G. de Lorenzi and M. D. German, Crack extension modeling with singular quadratic isoparametric elements. Inr. J. Fracrure 12, 647451 (1976). [20] A. R. Iggrafea, Discrete fracture propagation in rock laboratory tests and finite elements analysis. PhD Thesis, University of Colorado (1977). [21] M. Isida. Stress intensity factors of internally cracked plates. Inr. J. Frucrure Mech. 7, 301-316 (1971). (Received 12 April 1989)
Fracture
Mechanics
Printed in Great Britain.
Vol. 36, No. 3, pp. 451-458,
$3.00+ 0.00 0013-7944/90 0 1990 Pergamon Press plc.
1990
SINGULAR SHAPE FUNCTIONS FOR AN AUTOMATIC MESH GENERATION AROUND SINGULARITIES G. TSAMASPHYROS
and M. PAPAIOANNOU
Section of Mechanics, Department of Engineering Sciences, The National Technical University of Athens, 5, Heroes of Polytechnion, GR 157 74, Greece Abstract-A
proper gradual condensation of the finite element mesh around singularities is obtained by using a mapping technique. This technique is based on the construction of special singular shape functions. An advantage of the method is that it is entirely automatic and has not the drawbacks of a similar previous method.
INTRODUCTION LAST YEAR’S finite element method became one of the most powerful procedures for a numerical approximate solution of boundary value problems. One of the inherent difficulties of the method is the need of a great number of input data related to the topology of the mesh. In the case of problems without singularities the problem has already been solved by the use of various mesh generation schemes[l+. But the standard automatic mesh generation schemes cannot confront problems where a stress singularity appears. In this case a concentration of elements is needed around the singular point even in the case where special elements are used[5-81 (see also [lO-141). Thus the finite element mesh must present a gradual transition from small elements, around the singular points, to large elements, in the areas where regular stresses are expected. This is a very difficult task since the elements must be as regular as possible[9] and the optimum size or the rate of the gradual change of the size of the elements are not a priori known. In the present paper a method for an automatic generation of optimum mesh is presented. The method allows to confront singularities of any order and it is based on the use of special shape functions which maps a regular mesh to a mesh with an optimum concentration of elements. In fact the basic idea of the present method is the same with the previous ones[13-141, i.e. we use an approximate transformation which removes the singularity from the solution of the problem in the transformed domain.
MAPPING PROCEDURE Let K?: (Fig. la) equal a finite element. This element using isoparametric shape functions Mj(%, H) is mapped to a regular shape SQ(Fig. lb) in the parametric plane. Any field of parameters L can be described in the Q: domain by: L = f M,(E, H)L, j= I
where the (Lj}/“_I are the values of the parameter L at the m nodes of the element. In the case of linear elements (triangular or quadrilateral) the quantity: pi=l--zMi
(2)
represents the non dimensional radius from node i. This assertion remains true in the more general case of the element of the Fig. l(b) only for the corner node 1 since along the edges l-m and 1-2, M, behaves linearly. In the general case the quantity pi can be considered in the vicinity of the node i, as a measure of the non dimensional radius. Consider now a plane body (Fig. 2) having a singularity at A. We subdivide this body to a number of “super elements” 0; and let node 1 of each super element coincide with the singular point A. If we denote by r the radius from A, then the stresses and strains behave like O(r”-‘). Obviously the same (2 - 1) singularity will exist in each super element and the same singularity 451
452
G. TSAMASPHYROS
and M. PAPAIOANNOU
3 H
n-l-1
2
LE
5 mcl
1
m
(b) Parametric plane or
(a) Initial plane
mapped plane Fig. 1. (a) Initial super element Q:, (b) element q into the parametric plane obtained from Q using singular mapping.
[O(p”-I)] will exist after the mapping into the parametric plane. (This will be explained by the linearity of the mapping.) The above singular behavior is due to terms like p which exist in the displacement field. If a change of variable, like: Pi= Pt
(3)
is performed, then the stress field in the transformed p-plane would be regular around the node i. That means that the strain and the stress field and consequently the displacements field is reproducible by polynomials in the new plane. In this way we can use conventional polynomial shape functions {(Mi):_ 1 in order to describe the displacements u = (uX, uY} in the transformed plane: u=Mq’
(4)
q’ being the vector of nodal displacements of the super element Q. Another immediate consequence is that no special mesh is needed in the transformed CL q )-plane. Thus a special mesh can be constructed from a regular one in the mapped plane by using the inverse transformation. This fundamental observation has already been used[13-141 by taking the conformal mapping as the transformation which could perform a change of variables like the one described by the relation (3). The essential disadvantage of this transformation is that it distorts the boundaries of each super-element &2:. That means the boundaries of Qt (which is the image of a:) cannot be
Fig. 2. Plane body with singularity at A. Subdivision into “super elements” q.
Singular shape functions
453
represented by polynomials. This creates some difficulties in the exact representation of the boundaries of the initial domain. Methods to overcome these difficulties are described in the above papers[l3-141, but a method without this drawback is desirable. In the method which will be described here the initial super element Pz will be mapped onto the unit square (Fig. lb). This mapping is based on the previous described possibility to represent the normalized radius from node 1 through the shape function, see relation (2). To this end we will represent by: JJ,=l-M,
(5)
the non dimensional radius from node 1 (in the transformed plane). The shape functions N which could describe the singular behavior of strains in the initial plane must be singular. Nevertheless we consider again the quantity: p, = 1 -N,
(6)
as a measure of the normalized radius from node 1. Using (3) we have: (1 -M,)
= (1 - iv,)“.
(7)
That means: Ni = 1 - (1 -M,)“? Now taking into consideration
(8)
that Ni must satisfy the relations: Ni(
qj)
=
$,
i, j = l(l)m
6,;
N(419V)=
(9)
l
(10)
we can easily construct the other shape functions Ni (i = 2(l)m). To this end we take into account
that M satisfies relations similar to the relations (9) and (10). Nj=(l
5 MA. j=2(l)m
-N,)Mj/
(11)
k==2
These relations (8) and (11) have been used locally by Okabe et al.[ 15-l 61 in order to construct singular elements. The possibility to map a part of the whole body has not been explored. The previous idea can be generalized by using appropriate weight coefficients aj. Thus Nj=(l
-N,)ajMj/
2
(akitfk).
j =
2(l)m
(12)
k=2
The mesh obtained by the singular mapping as well as the influence of parameters aj in the distortion of the mesh is depicted in Fig. 3 (cases a-c). The above analysis is based on the representation of the dimensional radius in terms of 1 - M, . But the dimensional radius pl can be represented by the relation: pI = maxE
+ 1)/2, (tt +
WI.
(13)
If the above relation is substituted into relation (3) expressions similar to relations (8) and (12) could be obtained. Another possibility is to use a combination of (2) and (13) in order to represent the dimensional radius; for example: PI = PZP,p,8*
a+/3=1
(14)
where p,. is given by (2) and pls is given by (13). In the case where the body has many singular points it is possible to use successively the previous transformations for each singular point in order to obtain the final mapping. Consecutive transformations can be used in the case where Mi is not a linear function and the edges l-m and l-2 are curvilinear. In that case using the parametric transformation (l), the initial super element Pz is mapped to a square in the parametric plane. We map afterwards the square using the singular transformation (8) and (12) but supposing that the shape functions Mi
G. TSAMASPHYROS
454
and M. PAPAIOANNOU
(b)
(4 Fig.
3. Singular
mesh obtained from relations (8) and (12) with: and (c) a,=a,=O.S,a, = 4 = a2 = a, = f.O,a,=OS;
(a) a,=~,
=a4= I; (b)
1.0.
(1 = l(l)m) involved in these formulas are the shape functions M; (i = 1(1)4) of the four node quad~latera1. The above examples concerned plane problems but the extension to three-dimensional problems is obvious. On the other hand we are only interested for problems where the singularity is at described points. We consider now the case where the singularity is along a particular line S of the body. In that case one or more elements into which the body is subdivided, must have one side along S. The normalized radial distance from S is p = 1 -(Nj+Nj+l)
(15)
where nodes i and (i + 1) are along S. If the transformation (3) is performed we would have: 1 - (N, + Mz) = (1 - [N, + NJ)“.
(16)
In the case of a four-node element (16) gives: 043 + A44)= (A$+ N4)1:
(17)
and a possible choice of A$ and A$ will be: N, = Mi(M3 + M#‘“-1;
i = 3,4
(18)
Singular shape functions
455
and consequently: N, = M,(l - (M3 + M,))““-‘.
i = 1,2
(19)
An illustration of the distortion of the regular mesh produced by the above formulas (18) and (19) (with 1 = l/2) is given in Fig. 4(a). We use now the previous transformation (15) not for a distance from an arbitrary line S but from the Ox and 0~ axis. In this way we obtain the normalized coordinates x and y which are given by: x=1-(N,+N,)=N*+N,,
JJ=l-(N,+N,)=N,+N,
(20)
in the case of a four node quadrilateral. Using the transformation (17) for x and y we have: &+M3=(Nz+N,)i,
M,+M,=(N,+N,)”
(21)
01
N2 + N3 = (M2 + h4p,
N, + N2 = (M, +
M*)‘?
(22)
Consequently, N, - NS = (M, - M,)[(M, + M,)“‘l-
(M3 + M,)“V(M,
- M,)
(23)
and thus we can take:
w’il/(~3 + 4) N, = A43 [(M,+ A&)“”- w3 + W’“lI(~3 + W) N, = M, [(M, + &)“A - (MS +
(24)
N2 = (M2 + &)I’*: - N,
I
N,=l-N,-Nz-N,.
In the special case of a crack (1 = l/2) the previous formulas can be simplified. An illustration of this transformation is given in Fig. 4(b). The above examples concerned plane problems but the extension to three-dimensional problems is obvious. NUMERICAL
EXAMPLES
In order to check the efficiency of the proposed method a rectangular sheet of width 26, height 2h, containing a straight crack of length 2u was considered. The plate was subjected to uniform
(a)
(b)
Fig. 4. Singular mesh obtained (a) from formulas (18) and (18) or (b) from formulas (24) with I = l/2.
G. TSAMASPHYROS
456
and M. PAPAIOANNOU
Table 1. Percentage error of the stress intensity factor for a central crack with a/b = 0.5 using the singular mappin,g of Fig. 3(a) and varying number of nodes hr, of coefficients C, and collocation points n, of the formula (25) Total number of nodes (N)
Number of nodes along the crack (n)
21 65 121 133 173 181 199 245 309
3 5 9 7 9 11 9 11 11
Percentage difference between exact value(21] and
formula (25) 13.29 -3.98
- 12.08 -5.46 -8.39 - 12.43 -7.03 -9.33 -7.06
special collocation formula[ 141 -2.82 -0.85 - 0.05 0.10 0.30 0.25 0.38 0.47 0.57
tension a, along the parallel to the crack edges. Three values (1.0,0.7,0.4) of the ratio h/b and various values of the ratio a/b were considered in order to study the behavior of the procedure when the loaded boundary comes near the crack. In order to reduce the computational effort the symmetry was exploited, thus a quarter of the sheet is descritized by the above described method and a singular finite element mesh is obtained. After that a classical finite element program can be used. Since, any of our singular finite element meshes, has curvilinear elements, isoparametric elements must be utilized. We have selected quad~lateral elements and especially the 8-node ones. Since, in the case of cracked body the order of sin~a~ty 1 is equal to (l/2) we can easily see from formula (7) that the crack tip element is automatically reduced to a quarter point singular element[ 17j. The stress intensity is calculated by using the special collocation formula described in [14]. The stress intensity factor can also be obtained by the formula: K, - ZK,,=
E,/%(4AU,-AU,)
(25) 4&K=% where AU, (or AU, ) is the difference of complex displacements of the opposite points 2 and 2’ (or 1 and l’), with 2’ (or 1’) being nodes of the lower lip of the crack and 2 (or 1) nodes of the upper lip of the crack. The nodes 1 and 2 or the nodes 1’ and 2’ belongs to the same crack tip element. Node 0 coincides with the crack tip. This formula is given in [18] and inde~nden~y in 119,201. First the influence of the number N of nodes is studied. To this end the case a/b = 0.5 and h/b = 1.0 is considered with a varying number of total nodes iV and a varying number of nodes m along the crack lips (m > n), The results figure in Table 1 where the first column corresponds to the number of nodes N, the second one to the number n of nodes along the crack lips, The third Table 2. Percentage error of the stress intensity factor for a central crack with a/b = 0.5 using an arbitrary mesh, and a varying number of nodes N, of coefficients C,, and collocation points n Total number of nodes (N) 53 ;; 121 133 199 245 309
Number of nodes along the crack (u)
Percentage difference between exact vaIue(21] and formula (25)
special collocation fO~~~I4]
7 7 7 9 7 9 11 11
81.81 106.80 124.99 142.21 148.40 181.41 199.88 218.56
29.25 34.93 37.88 41.13 41.17 45.20 47.12 48.35
451
Singular shape functions Table 3. Percentage error of the stress intensity factor for a central crack with h/b = 1.Oand various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 200)
alb 0.3 0.4 0.5 0.6 0.7 0.8
Number of nodes along the crack (n) 1 9 9 11 11 11
Percentage difference between exact value[ll] and
formula (25) - 8.94 - 12.44 -7.03 -5.14 -5.33 -2.60
special collocation formula[ 141 -0.86 -0.41 0.38 -0.85 - 10.12 0.91
column gives the difference between exact value[21] of the stress intensity factor and the one obtained from the formula (13). As we can see the results of the collocation procedure are very satisfactory even in the case of a very coarse mesh. On the contrary the accuracy of the formula (13) is not satisfactory and this contradicts with the claims of Ingraffea et a1.[20]. In order to have an idea of the improvement obtained by the use of the above described singular finite element mesh (see formula 7 and 10) a regular finite element mesh is constructed with the same number of nodes. The results of the above mesh are shown in Table 2. As we can see the improvement is spectacular. For h/b = 1.O,0.7,0.4 and a varying crack length the difference of the calculated stress intensity factor from the exact one is given in Tables 3,4,5, respectively. The results of Table 3 are obtained from finite element patterns with approximately N = 200 nodes. In the case of Tables 4 and 5 the number of nodes was N = 150 and N = 130, respectively.
Table 4. Percentage error of the stress intensity factor for a central crack with = 0.7 and various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 150)
h/b
alb
Number of nodes along the crack (n)
Percentage difference between exact value[21] and
formula (25)
0.3 0.4 0.5 0.6 0.7
5 7 9 9 11
-1.58 -5.15 - 5.47 -4.38 -3.41
special collocation formula[ 141 -0.71 0.46 -1.38 - 1.43 -1.18
Table 5. Percentage error of the stress intensity factor for a central crack with h lb = 0.4 and various values of a/b using the singular mapping and a varying number of coefficients C,. (Number of nodes N = 130)
a/b 0.3 0.4 0.5 0.6 0.1
Number of nodes along the crack (n) I 9 11 11 11
Percentage difference between exact value[21] and formula (25) -4.82 -5.65 -6.00 -5.44 -5.34
special collocation formula[l4] 1.64 -2.12 -1.79 -2.06 -2.18
458
G. TSAMASPHYROS
and M. PAPAIOANNOU
REFERENCES [I] W. P. Buell and A. A. Bush, Mesh generation-a survey. J. uppf. Mech. 332-338 (1973). [2] R. D. Shaw and R. G. Pithen, Modifications to Suhara-Fukuda method of network generation. Inr. J. namer. Merh. Engng 12, 93-99 (1973). (31 S. E. Umanskij, Triangulation algorithm and program two-dimensional domains of any shape. Srrengrh Marer. 9, 700-706 (1978). [4] I. Imafucu, Y. Kodera and M. Sayawaki, A generalized automatic mesh generation scheme for finite element method. Int. J. numer. Meth. Engng 15, 713-731 (1980). [S] 0. C. Zienkiewicz and D. W. Philips, An automatic mesh generation scheme for plane and curved surfaces by “isoparametric” co-ordinates. hr. J. numer. Me&. Engng 3, 519-528 (1971). [6] R. H. Gallagher, A Review of Finite EIemenr Techniques in Fracture Mechanics (Edited by A. R. Luxmoore and D. R. J. Owen), pp l-25. University of Wales, Swansea (1978). [7] K. D. Henshell, Crack tip fmite elements are unnecessary. Inr. J. numer. Merh. Engng 9, 495-507 (1975). (81 0. C. Zienkiewicz, The Finire Element Method, 2nd edn. McGraw-Hill, New York (1977). [9] M. Zlamal, A finite element procedure of the second order of accuracy. Numer. Murh. 14, 394Xt2 (1970). [IO] J. J. Guydish and J. F. Fleming, Optimization of the finite element mesh for the solution of fractural problems. Engng Fracture Mech. 10, 31-42 (1978). [I I] J. J. Guydish, J. F. Fleming, J. R. Pentx and C. E. Runnion, The finite element method vs the edge function method for linear fractural analysis. Engng Fracrure Mech. 13, 43-55 (1980). [12] P. R. Brown, A non-interactive method for the automatic generation of finite element meshes using the Schwartz-Christoffel transformation. Compur. iUerh. appl. Mech. Engng 25, 101-126 (1978). [ 131 G. Tsamasphyros and A. E. Giannakopoulos, The mapped elements for the solution of cracked bodies. Compur. Merh. appl. Mech. Engng 49, 331-342 (1985). [ 141 G. Tsamasphyros and M. Papaioannou, Mapping technique for the automatic mesh generation around singularities. Compur. Srrucr. 29, 815-823 (1988). [IS] Y. Yamada, Y. Ezawa, I. Nishiquchi and M. Okabe, Reconsiderations on singularity or crack tip elements. Inr. J. numer. Merh. Engng 14, 1525-1544 (1979). [16] M. Okabe, Fundamental theory of the semi-radial singularity mapping with applications to fracture mechanics. Compur. Merh. appl. Mech. Engng 26, 53-73 (1981). [17] R. S. Barsoum, Application of quatratic isoparametric finite elements in linear mechanics. Inr. J. Fracrure 10,603-605 (1974). [18] G. Tsamasphyros and A. E. Giannakopoulos, Automatic optimum mesh around singularities using conformal mapping. Engng Fracture Mech. 23, 507-520 (1986). [19] C. F. Shih, H. G. de Lorenzi and M. D. German, Crack extension modeling with singular quadratic isoparametric elements. Inr. J. Fracrure 12, 647451 (1976). [20] A. R. Iggrafea, Discrete fracture propagation in rock laboratory tests and finite elements analysis. PhD Thesis, University of Colorado (1977). [21] M. Isida. Stress intensity factors of internally cracked plates. Inr. J. Frucrure Mech. 7, 301-316 (1971). (Received 12 April 1989)