- Email: [email protected]
A methodology for crack tip mesh design
Engineering Fracture Mechank's Vol. 50, No. 5/6, pp. 713-726, 1995
Pergamon
Copyright © 1995 ElsevierScienceLtd Printed in Great Britain. All rights reserved 0013-7944/95 $9.50+ 0.00
0013-7944(94)E0056-M
A METHODOLOGY
FOR CRACK
TIP MESH DESIGN
FERNANDO C. M. MENANDRO, E. T. MOYER, JR and H. LIEBOWITZ School of Engineering and Applied Science, The George Washington University, Washington, DC 20052, U.S.A. Abstraet--A methodology for crack tip mesh design is developed which consists of comparing the mesh geometric parameters against the accuracy of the finite element solution. By successive changes in the mesh parameters a near optimal mesh can be obtained. This was done here for two-dimensional linear elastic single mode problems. The direct displacement extrapolation method for stress intensity factor estimation is used.
1. INTRODUCTION ThE FINITEELEMENTMETHOD(FEM) has been in use for about 30 years and is now established as a powerful numerical method for solution of partial differential equations. It is widely used for fracture problems, although some numerical difficulties occur. By analyzing the FEM response of two-dimensional linear elastic fracture problems, a methodology for mesh generation was developed. Significant amounts of research have been performed in automatic grid generation, as well as adaptive meshing schemes. There are, however, very few mathematical developments of these subjects for fracture problems since the theory is based on smooth functions. This led to this work's approach: to determine which geometric characteristics dominate the FEM mesh behavior and generate a near optimal mesh based on those parameters [1]. FEM work was carried out by analyzing some problems with known analytical solutions with a starting mesh and by successive changes in the mesh searching for a near optimum one. These analyses were performed with the MARC [2] software on an ALLIANT computer. The method adopted for the stress intensity factor (SIF) estimation was direct displacement extrapolation (DDE) [3]. Background in FEM analysis for linear elastic fracture mechanics (LEFM) is given only in relation to this work. The problem itself was attacked in three basic steps: (a) one-parameter studies; (b) two-parameter studies; (c) angular studies. Example problems were then tested for the desired accuracy [3]. 2. FEM FOR LEFM
The method adopted for SIF estimation in this work is the DDE method of Chan et al. [4]. The method consists of calculating the effective stress intensities at a given direction from the crack tip (usually 90 ° and 180°) and asymptotically extrapolating their values to radius zero (crack tip). Two-dimensional parabolic isoparametric quadrilateral elements were adopted, with the exception at the crack tip, where degenerated quadrilateral quarter-point elements were used. 3. ONE-PARAMETER STUDIES
In this section the accuracy of the crack tip mesh for problems in which the only operative parameter is the crack length was studied. Only two-dimensional problems in LEFM were studied. The approach chosen was to define a model problem (with known analytical solution) and study the convergence of the FEM solution due to systematic changes in the mesh. Meshes with 713
714
F . C . M . MENANDRO et al.
errors equal to or smaller than 5%, 3%, 2% and 1%, respectively, were obtained with the model problem and tested for other problems (to see if the accuracy was maintained).
3.1. Model adoption The model problem defined was the semi-infinite body with an edge crack subjected to normal stress perpendicular to the crack. The analytical values for the SIFs are shown below: K~ = 1.1215ax/(zra ) -~ 1.987807
(1)
KH = 0
(2)
for o- = 1 (applied stress) and a = 1 (crack length). The error in the analytical solution due to the finite dimensions of the model adopted (Fig. !) is of small order when compared with FEM error (i.e. for a finite dimension specimen, a power series correction factor should be introduced [4]). Since, for the dimensions adopted (width = 40, a = 1), this factor was smaller than 0.05% and the error of the FEM solution was an order of magnitude higher (~>0.1%), the model was considered to represent the semi-infinite body accurately. The results for K~ are presented as the ratio between the FEM solution and the analytical value of 1.987807 given previously. Near the crack tip, a mesh of circular rings of elements is needed. The number of rings of elements was considered as one of the parameters of the crack mesh to be studied, as well as the outer radius (r) of the last ring, as shown in Fig. 1. The third parameter was the distance (l), measured along the crack, between the crack tip and the edge of the small rectangle of refined mesh (crack tip superelement region), also shown in the figure. A rectangular superelement region was adopted since it is the most natural for straight crack problems and gives results as accurate as a circular one [6]. To use dimensionless parameters the actual parameters adopted were the number of circular rings of elements, r/l and l/a (a is the crack length). The number of elements in each ring was eight (for 180°), which have been proved to be optimal [7].
Rings of elements
Fig. 1. Finite element mesh of model problem.
A methodology for crack tip mesh design Table
1.
Kt(FEM~/giAnalytical(l/a=0.40 r/I = 0.625)
and
Table 2.
715
Kt(FEM)/Klanalytical(lla=0.40 r/l = 0.500)
Rings
180 °
90 °
Rings
180°
90 °
1 2 3 4 5 6 7 8 9 10
0.96007 0.97564 0.98225 0.99280 0.99306 0.99347 0.99349 0.99358 0.99367 0.99366
1.04953 1.04657 1.04206 1.02853 1.01372 1.01027 1.00848 1.00633 1.00403 1.00401
1 2 3 4 5 6 7 8 9 10
0.96530 0.97944 0.98526 0.99437 0.99320 0.99362 0.99361 0.99370 0.99379 0.99378
1.05197 1.04617 1.04111 1.03005 1.00975 1.00707 1.00547 1.00357 1.00148 1.00147
and
3.2. Number of circular rings tests The starting mesh selected is shown in Fig. 1 and has an edge of the superelement region over crack length (l/a) parameter equal to 0.40 and a radius of the circular mesh over edge of superelement region (r/l) parameter equal to 0.625. The results are shown for the 90 ° and 180 ° lines, for a number of rings from one to 10 (Fig. 1), in Table 1. It should be noted that the results from the 90 ° line converge to the correct solution from a higher value, while the 180 ° line gives lower bound results; this fact was used later to determine the optimal mesh. The same problem was solved with the r/l parameter equal to 0.500 and 0.750; the results are shown in Tables 2 and 3. The results along 0 = 180 ° demonstrated an invalid characteristic, namely that four circular rings of elements yield the best results. The D D E method was implemented to take 10 nodes along a radius from the crack tip, calculate the effective stress intensity factors and by linear regression extrapolate the value for the crack tip. For five circular rings of elements all the nodal displacements used were in the refined circular mesh, while for less than five there were nodes outside of this refined mesh (for one and two rings there will be less than 10 values for the linear regression in the 180 ° line, respectively six and eight) where the equations utilized no longer held. Since for most crack problems the geometry is such that the total mesh has many more elements than the superelement region alone, the change in the problem size in relation to the number of circular rings of elements is small (e.g. for this model problem approximately 1% more nodes for each additional ring). Therefore, little additional cost is introduced when using a more refined mesh. As expected, the accuracy showed improvement with more rings (and more elements), which was due to nothing more than convergence of the FEM (after the singularity was well represented). Table 4. Seven rings
Table 3. Kl(FEM)/KiAnalytica I r/l = 0.750)
(l/a
= 0.40 and
Rings
180 °
90 °
1 2 3 4 5 6 7 8 9 10
0.95509 0.97215 0.97959 0.99165 0.99285 0.99318 0.99330 0.99345 0.99356 0.99356
1.04671 1.04637 1.04224 1.02661 1.01736 1.01314 1.01109 1.00885 1.00629 1.00629
EFM 50/5-6~1
(l/a
= 0.40)
r/l
180 °
90 °
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03639 1.01996 1.01045 1.00463 1.00093 0.99853 0.99693 0.99586 0.99513 0.99465 0.99430 0.99406 0.99389 0.99377 0.99371 0.99364 0.99361 0.99357 0.99355 0.99353 0.99351 0.99349
1.03514 1.01967 1.01110 1.00623 1.00341 1.00192 1.00120 1.00100 1.00111 1.00147 1.00190 1.00241 1.00296 1.00356 1.00423 1.00483 1.00547 1.00609 1.00669 1.00731 1.00790 1.00848
716
F.C.M. Table 5. Six rings (l/a = 0.40)
M E N A N D R O et al. Table 6. Five rings (l/a = 0.40)
r/l
180 °
90 °
r/I
180°
90 "~
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03692 1.02035 1.01075 1.00485 1.00113 0.99869 0.99704 0.99595 0.99520 0.99471 0.99434 0.99410 0.99392 0.99381 0.99373 0.99367 0.99362 0.99359 0.99354 0.99353 0.99350 0.99347
1.03610 1.02059 1.01207 1.00717 1.00452 1.00308 1.00234 1.00220 1.00236 1.00276 1.00319 1.00380 1.00441 1.00508 1.00478 1.00542 1.00707 1.00677 1.00737 1.00804 1.00867 1.01027
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03497 1.01867 1.00929 1.00353 0.99993 0.99764 0.99611 0.99509 0.99443 0.99400 0.99371 0.99351 0.99339 0.99332 0.99327 0.99322 0.99320 0.99318 0.99316 0.99313 0.99310 0.99306
1.03576 1.02059 1.01240 1.00767 1.00519 1.00403 1.00354 1.00352 1.00387 1.00442 1.00506 1.00575 1.00654 1.00734 1.00816 1.00893 1.00975 1.01056 1.01135 1.01214 1.01293 1.01372
The rate of convergence might depend on the length parameters, but it can be proved that once the singular behavior is well represented the convergence itself is not dependent on them [8]. Comparing the results for various r/l, it was noted that the results improved for smaller values of the parameter (for meshes with five to 10 rings). This was understood as a need for a refined mesh region near the crack tip, but, as was discovered later, there is an optimum value beyond which the refined mesh is no longer sufficient to represent the asymptotic behavior. Comparing the results from the three tables, it was observed that for the same number of rings the accuracy increased for smaller values of the r/l parameter. This is shown in the next section. 3.3.
Tests for r/1
From the previous discussion it was decided to test the problem for small changes in the parameter r/l (Tables 4-6). Refinement of more than seven rings of elements was considered unnecessary and was not pursued. It is important to note that the 90 ° line showed what seemed to be a stationary point (minimum) where the solution was optimal, while for the 180 ° line the convergence was not clear (there was no stationary point in the function of r/l). The finding of a minimum was in agreement with the observations in the previous section, i.e. that the 90 ° line gives upper bound solutions. In the three tables the best result for the 90 ° line was at r/l = 0.275 (which gives a minimum) and for the 180 ° line at r/l = 0.200. For the 180 ° line, a maximum should be expected and could not be found; thus, the best result considered was the one which best approximated the analytical value. 3.4.
Tests for l/a
The
Continuing the search for an optimal solution, the l/a parameter was now analyzed (Table 7). r/l parameter was kept constant along with the number of circular rings of elements. From Table 7. r/I = 0.500
l/a
180 °
90 °
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.40
1.00219 0.99993 0.99825 0.99696 0.99598 0.99522 0.99464 0.99422 0.99361
1.00788 1.00637 1.00545 1.00486 1.00457 1.00448 1.00456 1.00481 1.00547
Table 8
l/a
r/l
180 ~
90"
0.24 0.24 0.32 0.32
0.200 0.275 0.200 0.275
1.00721 1.00213 1.00253 0.99745
1.00632 1.00324 1.00339 1.00067
A methodology for crack tip mesh design
717
Table 10. r/l = 0.275
Table 9. l/a = 0.32
r/l
180 °
90 °
I/a
18if'
90 ~'
0.250 0.275 0.300
0.99852 0.99745 0.99673
1.00092 1.00067 1.00068
0.30 0.32 0.34 0.36 0.38 0.40
0.99821 0.99745 0.99691 0.99646 0.99613 0.99586
1.00093 1.00067 1.00067 1.00065 1.00080 1.00100
then on, the tests were performed only for seven rings, which had shown good accuracy and convergence. Again, a stationary point was found for the 90 ° line, while the 180 ° line did not show the same behavior. The most accurate solutions were for l/a equal to 0.24 and 0.32 for the 180 ° and the 90 ° lines, respectively. 3.5. Combined tests Combining the best results for both parameters, the desired mesh should be obtained (Table 8). Since the most accurate mesh parameters might not be independent, it was important to test the neighborhood of this solution (by separate small changes in each of the parameters). Tables 9-11 show these test results. Since the 90 ° line had always given upper bound solutions and the 180 ° line had not shown the same consistency (the results seemed to be more affected by changes in the mesh parameters, sometimes diverging), the mesh with l/a = 0.32 and r/l = 0.275 was considered to be the most accurate (Table 8). The mesh selected in the previous paragraph was obtained by combining independently studied results for r/l a n d / / a . Since these parameters might have some correlation, there was no guarantee that this mesh was really the optimal one. To address this problem, the SIF for the 90 ° line was analyzed as a function of r/l and l/a and the neighborhood of the selected mesh was tested (searching for a minimum). As expected, the parameters were interrelated. The best result was actually for l/a = 0.36. As mentioned in Section 3.2, the FEM solution converges after the singular behavior is represented: this led to the use of a fixed number of circular rings of elements (seven) for the rest of the tests. Once the optimal length parameters have been found, the convergence due to the number of elements (or the number of rings) should be studied in order to define meshes for the desired accuracies. This is shown in Table 12. Calculating the percentage difference between the analytical solution and the FEM (90 ° line) asymptotic displacement solution, the error was 3.188% for four rings and 0.300% for five rings. This drop in accuracy had to do with the D D E method as mentioned before, and gave a warning about using less than five circular rings of elements. The most accurate acceptable mesh for an edge crack was thus the one with the following parameters: I - = 0.36 a r = 0.275
and at least five rings of elements. Table 12. I/a = 0.36, r/I = 0.275
Table 11. l/a = 0.36
r/1
180 '~
90 ~
0.250 0.275 0.300
0.99752 0.99646 0.99573
1.00088 1.00065 1.00073
Rings
180'
90 ~
1 2 3 4 5 6 7
0.97714 0.98858 0.99308 1.00011 0.99568 0.99656 0.99646
1.05441 1.04419 1.03871 1.03188 1.00300 1.00181 1.00065
718
F.C.M. MENANDRO et al.
Since the five-ring mesh had already given results with better than 1% accuracy, this was selected as the desired mesh (for 10%, 5%, 3%, 2% and 1%), which should be tested for some other problems. 3.6. Final tests Once the desired mesh had been obtained, testing for other one-parameter problems could either support the chosen mesh or disregard it, suggesting changes. The problems used to test the chosen mesh were the well known centrally cracked strip subjected to tension and the three-point bend specimen. For the centrally cracked strip (Fig. 2), Isida [9] developed (based on Laurent series expansions of the complex stress potentials) a solution which for our purposes could be considered exact. The mode I and II SIFs are given by: Kt = a x / ( ~ a ) F ( 2 )
(3)
Kn = 0,
(4)
where tr is the applied stress, a is the half crack length, 2 = 2a/W, W is the plate width and F(2) is a 36-term power series of 2 given in ref. [9]. For the three-point bend specimen (Fig. 3), a solution considered to be accurate up to 0.5% was obtained by Gross and Srawley in 1965 [10]. This solution was used only for comparison of results, since the reported accuracy was below this study's requirements. The SIFs, as presented by Sih [11], are:
K~ = 6_.~h,~ __M./-aa F(2)
(5)
KH = 0,
(6)
where M is the applied moment, a is the crack length, t is the thickness of the beam, h is the height of the beam, 2 = a/h and F(2) is a five-term power series of 2 given in ref. [11]. The results of the tests, already divided by the adopted solution, are shown in Table 13. The error was approximately 1% for every problem, which indicated that the mesh of the previous section might not be as accurate as expected. Some values greater than 1% were found where the expected errors should be less than 1%. All this evidence led to the application of the same procedure used in Section 3.5 to search for a better mesh.
tttttttttttt
v
b
J ' 2 . J'
b --
W
~
P
s T
Fig. 2. Centrallycracked strip.
Fig. 3. Three-pointbend specimen.
719
A methodology for crack tip mesh design Table 14. W = 5a, gt(FEMl/gtAnalyticaI
Table 13. KI~FEM)/KIAaalyti~I Model Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Three-point bending
W ,> 2a W/2a = W/2a = W/2a = W/2a =
5 4 3 2
180 °
90 °
l/a
r/l
180 °
90 °
0.99031 0.99044 0.99040 0.99012 0.98817 0.98686
! .01159 1.01192 1.01201 1.01206 1.01141 0.98879
0.38 0.38 0.38 0.38
0.275 0.250 0.225 0.200
0.99002 0.99091 0.99233 0.99449
1.01254 1.01163 1.01122 1.01143
0.36
0.300
0.98986
1.01302
0.36
0.275
0.99044
1.01192
0.36 0.36 0.36
0.250 0.225 0.200
0.99136 0.99273 0.99493
1.01114 1.01066 1.01101
0.34 0.34 0.34 0.34 0.34
0.300 0.275 0.250 0.225 0.200
0.99038 0.99098 0.99189 0.99329 0.99545
1.01242 1.01143 1.01066 1.01031 1.01062
0.32 0.32 0.32 0.32 0.32
0.300 0.275 0.250 0.225 0.200
0.99107 0.99165 0.99255 0.99397 0.99613
1.01197 1.01098 1.01025 1.01004 1.01043
0.30 0.30 0.30 0.30 0.30
0.300 0.275 0.250 0.225 0.200
0.99197 0.99253 0.99343 0.99481 0.99703
1.01177 1.01079 1.01013 1.00983 1.01044
0.28 0.28 0.28 0.28 0.28
0.300 0.275 0.250 0.225 0.200
0.99305 0.99363 0.99453 0.99593 0.99810
1.01163 1.01076 1.01013 1.00997 1.01048
0.26 0.26 0.26 0.26 0.26
0.300 0.275 0.250 0.225 0.200
0.99452 0.99508 0.99598 0.99737 0.99956
1.01194 1.01108 1.01052 1.01040 1.01103
Again, the neighborhood of the solution was tested, this time for the centrally cracked strip problem. The width over length (W/2a) ratio adopted was 5. The results are presented in Table 14: the result in italics was the starting mesh and the final result is given in bold. These results were obtained in the same way as in section 3.5, the final result being selected as a minimum for the 90 ° line. The error for the final result was less than 1% for both directions (90 ° and 180°). Again, the results for the 180 ° line were not considered due to the already discussed argument. To finish this analysis it was necessary to utilize this final mesh on the previous problems in order to check their accuracy. The results, divided by the adopted theoretical solution, were as follows: Edge crack (model problem) 180°:
900:
KI(FEM)=0.99988 KI(Th.)
KItFEM)=1.00359 KI(Th.)
Cracked infinite plate (modeled as strip with W = 108a) 180°:
KI~FEM)= 0.99472 KI(Th.)
90o:
KI(FEM) -- 1 . 0 0 9 5 7 gl(Th.)
720
F.C.M.
M E N A N D R O et al. 12-5-1990 m MENTAT
ORTHO V.P.O.O00e+O00.O00e+00
0.100e+O1
X J
Fig. 4. Crack tip mesh.
Three-point bend specimen 180°:
KI(FEM)=
0.99078
gl(Th.)
900: KI~FEM)=0.98521. gl(Th.)
The three-point bend problem continued to show up the error detected previously, but only for one direction. If the SIF was calculated as the mean of the results in both directions, the accuracy of the solution would improve. The explanation of this discrepancy is probably the lack of precision of the adopted theoretical solution, as mentioned previously. F o r the other two problems, the error in the solution was, as expected, less than 1% and therefore this mesh was adopted as optimal. The result for the edge crack seemed to be less precise than the one obtained in the previous section, but as mentioned earlier the second decimal digit was beyond this study's precision. Thus, the optimal mesh for problems in which the only operative parameter is the crack length is: l - = 0.30 a r
-! = 0.225 and at least five rings of elements. This mesh is shown in Fig. 4. The next section studies this mesh for problems in which more than one length parameter must be considered to define the problem geometry. 4. T W O - P A R A M E T E R S T U D I E S
Two-parameter problems constitute ones in which the crack length is no longer the only length parameter that determines the required crack tip mesh. This occurs when the other dimensions of the specimen are of the same order of magnitude as the crack length. Rules that take into account geometric parameters other than the crack length were developed for mapping the crack tip superelement mesh.
A methodology for crack tip mesh design
721
4.1. Model definition The adopted model problem was again the centrally cracked strip (Fig. 2) of Section 3.6, now with different values for the width over crack length (W/2a) parameter. The theoretical solution for this problem is repeated below for modes I and II: KI = a ~/(na )F(2)
(7)
K,, = O,
(8)
where e is the applied stress, a is the half crack length, 2 = 2a/W, W is the strip width and F(2) is a 36-term power series of 2 given in ref. [9]. 4.2. Tests Making the strip width (W) in a centrally cracked strip (Fig. 2) smaller in relation to the crack length (2a) increases the contribution of ligament 6 and, at some point, it starts to dominate. It should be expected that the smallest parameter dominates, which was confirmed for the tests run for g varying from four times the half crack length to one-quarter of it. The same mesh defined in the previous section was adopted for these problems with the difference that when 6 was smaller than a the length parameter used for the mapping was A These results divided by the respective analytical solutions are shown in Table 15. The results confirmed the suspected behavior and it can be seen that all were better than 1%. In the following section, another problem that occurred in two-parameter problems is discussed. 4.3. Tests for a _~ d A real two-parameter problem appears when the half crack length and the ligament have similar values, more specifically when 0,5 < g- < 2,
(9)
a
where 6 is the ligament and a is the half crack length (as shown in Fig. 2). When mapping the crack tip mesh for problems with 6/a agreeing with the inequality (9), some difficulties appeared in the transition from the superelement region to the regular domain mesh. This problem was studied using four different meshes, as shown in Fig. 5. All four meshes had the domain discretized the same way, with element sides roughly equal to the crack size. The transition between the superelement region and the regular domain elements was accomplished in four different ways from the following rules. To check the assumption that the dominant parameter was the smallest, one of the meshes used the other parameter for mapping. The transition adopted in this case seemed to be the most reasonable one (case 1). Using the smallest parameter for mapping created some ambiguity on what transition mesh should be adopted, the first option (case 2) being logically to use only four elements around the superelement to cover the entire extension (a + 6). If two elements are used to make the transition to the bigger side, one might opt to smooth these elements, or have the outer one rectangular (cases 3 and 4, respectively). By smoothing it is meant that the nodes are moved so that the distance between them is roughly the same [12]. The results of these analyses are shown in Table 16. The results for case 4 were usually the best, but since no significant difference existed between cases 3 and 4 the most straightforward was accepted, which was case 3. Table 15. Two parameters: for {~ >~a < a
I/a=0.30 I/3 = 0.30
3/a
180°
90 °
4 3 2 1 1/2 1/3 1/4
0.99481 0.99475 0.99442 0.99225 0.99064 0.99028 0.99009
1.00983 1.00988 1.00982 1.00870 0.99966 0.99663 0.99507
F . C . M . MENANDRO et al.
722
\ Case 1
Case 2
:
i
i
:
Case 3
Case 4 Fig. 5. Transition mesh cases.
M o r e tests s h o w e d t h a t all m e s h e s p e r f o r m well f o r a 10% difference b e t w e e n a a n d ~. In o t h e r w o r d s , f o r ~/a b e t w e e n 0.91 a n d 1.10, all m e s h e s g a v e a c c u r a c y to w i t h i n 1%. C a s e s 3 a n d 4 a d o p t the c r a c k l e n g t h as t h e m a p p i n g p a r a m e t e r , cases 3' a n d 4' the l i g a m e n t . T h e s e results a r e s h o w n in T a b l e 17. T h e p e r f o r m a n c e o f cases 3 a n d 4 for ~ = a was n o t t h a t desired d u e to the i n f l u e n c e o f a s p e c t ratios. It was s u g g e s t e d , as a g e n e r a l rule, t h a t f o r 0.91 < ~/a < 1.10 case 2 s h o u l d be a d o p t e d , case 3 r e m a i n i n g f o r t h e rest. A s a n o t h e r g e n e r a l rule it was n e c e s s a r y t h a t n o n e o f the t r a n s i t i o n e l e m e n t s ' side l e n g t h s be g r e a t e r t h a n t h e c r a c k (or l i g a m e n t ) length.
Table 16. 0.5 < ~/a < 2
d/a
Case
2/3
1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
3/4
4/3
3/2
180° 0.98633 1.00490 0.99646 0.99648 0.98946 0.99914 0.99460 0.99444 0.98996 0.99469 0.99133 0.99143 0.98937 0.99588 0.99247 0.99261
90 ° 1.00378 1.01721 1.00803 1.00799 1.00558 1.01238 1.00752 1.00730 1.01098 1.01069 1.00691 1.00701 1.01263 1.01174 1.00793 1.00807
Table 17. a ~
8/a 1.00
0.91
1.10
Case 2 3 3' 4 4' 2 3 4 2 3 4
180° 0.99225 0.98853 0.98986 0.98821 0.98882 0.99374 0.99016 0.98959 0.99310 0.99043 0.99031
90 ° 1.00870 1.00438 1.00661 1.00400 1.00550 1.00898 1.00552 1.00490 1.00938 1.00616 1.00601
A methodology for crack tip mesh design 9-12-1990 ca MENTAT
9-12-1990 oa MENTAT
ORTHO V.P. 0.000e+00 0.000e+00 0.100e+01
X )RTttO V.P. 0.000e+00 0.000e+00 0.100e+01
9-12-1990 oa MENTAT
ORTHO V.P. 0.000e+00 0.000e+00 0.100e+01
723
9-12-1990 ¢o MENTAT
x OR"~HO V.P. 0.000e+00 0.000e+00 0.100e+01
Fig. 6. Crack tip superelement regions. 5. A N G U L A R S T U D I E S Angular studies are the studies o f the influence o f angular refinement a r o u n d the crack tip on the accuracy. The accuracy o f the adopted mesh with eight elements a r o u n d the half circle was tested against six or 10 elements a r o u n d the half circle. The accuracy and the displacement fields were examined and compared. The choice o f eight elements was justified in view of the other options.
5.1. Model adoption In this section the refinement a r o u n d the crack tip in the angular direction was studied. The adopted solution o f eight elements a r o u n d the semi-circle was c o m p a r e d with options o f six and 10. The problem analyzed was the centrally cracked strip, as this is a well studied reference solution. The accuracy o f the solution, as well as the displacement fields, was studied to establish the optimal solution. The crack tip superelement regions are shown in Fig. 6. Meshes 6', and 10 and 10' were proposed and mesh 6 was obtained from 6' by a smoothing algorithm used at the non-circular region. The solutions o f the centrally cracked strip divided by the analytical one are shown in Tables 18-21. The non-dimensional parameters were taken close to the ones o f the previous analyses. It can be seen that although mesh 6 had a better behavior than 6', the same cannot be guaranteed for 10 elements. To improve the obtained solutions the accuracy against r/l was studied for cases 6, 10 and 10' (6' was obviously worse). Table 19. Mesh 6 (lla = 0.300) Table 18. Circular refinement test (l/a = 0.300; r/l = 0.200) Mesh 180° 90° 6' 0 . 9 8 7 7 9 0.99252 6 0 . 9 9 1 6 0 0.99561 10' 0 . 9 9 4 6 9 1.01659 10 0 . 9 9 3 6 9 1.01512
r/l 0.250 0.225 0.200 0.175 0.150 0.125 0.100
180° 0.98839 0.98964 0.99160 0.99473 0.99978 1.00819 1.02290
90° 0.99662 0.99584 0.99561 0.99635 0.99876 1.00420 1.01551
724
F.C.M.
M E N A N D R O et al.
Table 20. Mesh I0 (l/a = 0.300)
Table 21. Mesh 10' (l/a = 0.300)
r/I
180 °
90 °
r/l
180 °
90 °
0.250 0.225 0.200 0.175 0.150
0.99260 0.99299 0.99369 0.99501 0.99750
1.01697 1.01595 1.01512 1.01487 1.01581
0.250 0.225 0.200 0.175 0.150
0.99370 0.99405 0.99469 0.99592 0.99823
1.01849 1.01747 1.01659 1.01626 1.01704
The result for r/l equal to 0.150 was the best and, although a complete analysis was not performed, it would probably end up as the optimal. The best solution for mesh 10 was the one with r/l = 0.175 (it was a minimum for 90°). Unfortunately, it still did not give less than 1% in both directions. Originally, mesh 10' had the same r/l as mesh 10 and a similar inaccuracy in the results. For clarity, the best solutions for all meshes analyzed are shown together in Table 22. Undoubtedly the two best solutions are the ones for six and eight elements around the semi-circle. These meshes were thus compared thoroughly. Eight or six elements around the semi-circle were then utilized to compare displacement fields for the centrally cracked strip and the edge crack. The results of ux and uy are shown against 0 in Figs 7-10. The comparisons were taken at the outer border of the first and second rings of elements, which do not lie in the same distance for the different meshes. For the superposed charts, the solutions were scaled so that they were equal. It is worth remembering that the displacements in the x-direction are off by a constant and the zero was arbitrarily taken at 0 equal to 180 °. From Figs 7 and 8, one might assume that mesh 6 is better; however, since for the edge crack the pattern does not hold and references give eight elements as an optimum, it will stay as the selected mesh. Another reason for this decision is that the SIF results for the edge crack were slightly better and mesh 6 did not seem to present the 90 ° and 180 ° lines as upper and lower bounds, respectively, as experienced previously. The SIFs divided by the reference solution are shown in Table 23. For the purpose of this research, eight elements around the semi-circle were adopted as being precise enough. More research in this direction (0, angular) is advised before anything is concluded in problems that involve more than just LEFM. 6. CONCLUSIONS A methodology for crack tip mesh design was described which consists in determining what parameters influence the mesh behavior. By extensive analysis of FEM solutions the accuracy is checked against the mesh parameters. A near optimal crack tip mesh was obtained for single mode
o Ux - 8 elements
~,,,~¢,~,m---~
m Uy - 8 elements
t~
Uy/¢"
o Ux - 6 elements LxUy - 6 elements
Table 22. Angular comparison Mesh
180 °
90 °
8 6 10 10'
0.99481 0.99978 0.99501 0.99592
1.00983 0.99876 1.01487 1.01626
~/
,I
~ r
.o
Ux
:~..~'9~
~
,
.~ "~ °
Z Table 23. Angular test
t , ~ A'"
Mesh
Problem
180 °
90 °
6 8 6 8
Strip Strip Edge Edge
0.99978 0.99481 1.00584 0.99988
0.99876 1.00983 0.99618 1.00359
0
, 30
~ 60
J 90 Angle
~ 120
~ 150
\ N~ 180
Fig. 7. Displacements for centrally cracked strip. Six elements around crack tip, radius 0.009. Eight elements around crack tip, radius 0.0135.
A methodology for crack tip mesh design
725 • ~/)'m~'i
o Ux - 8 elements e~
Uy
-
y ~
8 elements
Uy ~(,/,f/
o Ux - 6 elements
~:
t, Uy - 6 elements
j
J
eO
O
z
0
30
60
90 Angle
120
150
180
Fig. 8. Displacements for centrally cracked strip. Six elements around crack tip, radius 0.018. Eight elements around crack tip, radius 0.027.
t3 Ux - 8 elements =
• Uy - 8 elements
St//
.z./
o Ux - 6 elements zx Uy - 6 elements
--
e~ @
0
30
60
90 Angle
120
150
180
Fig. 9. Displacements for edge crack specimen. Six elements around crack tip, radius 0.009. Eight elements around crack tip, radius 0.0135.
uy:J'~"~ J /
n Ux - 8 elements m Uy - 8 elements o Ux - 6 elements
8
~,
• Uy - 6 elements
/
i¢, Z 13 \
O
30
60
90
120
150
180
Angle Fig. 10. Displacements for edge crack specimen. Six elements around crack tip, radius 0.019. Eight elements around crack tip, radius 0.027.
726
F . C . M . MENANDRO et al.
c r a c k p r o b l e m s . It was s h o w n t h a t a l t h o u g h the s o l u t i o n s for the D D E m e t h o d are b e t t e r for the crack face (180°), the 90 ° d i r e c t i o n gives a b e t t e r p a r a m e t e r for analysis. It was also s h o w n t h a t for single m o d e L E F M p r o b l e m s o n l y the crack a n d l i g a m e n t (smallest d i s t a n c e b e t w e e n crack tip a n d the free edge) l e n g t h s i n f l u e n c e the F E M m e s h b e h a v i o r .
REFERENCES [I] F. C. M. Menandro, A knowledge based approach for crack tip mesh generation. M.Sc. Thesis, The George Washington University, Washington, DC (1991). [2] MARC general purpose finite element program, Marc Software International, Revision K3 (1988). [3] F. C. M. Menandro, E. T. Moyer, Jr and H. Liebowitz, A near optimum crack tip mesh. Engng Fracture Mech. 50, 703-711 (1995). [4] S. K. Chan, I. S. Tuba and W. K. Wilson, On the finite element method in linear fracture mechanics. Engng Fracture Mech. 2, 1-17 (1970). [5] H. L. Ewald and R. J. H. Wanhill, Fracture Mechanics. Edward Arnold and Delftse Uitgevers Maatschappij (1984). [6] E. T. Moyer, Jr and H. Liebowitz, Creep and fracture characteristics of materials and structures at elevated temperatures. ONR Final Report, no. N00014-84-K-00027 (1986). [7] E. T. Moyer, Jr and H. Liebowitz, Finite element mesh design--a knowledge based approach. Proc. Sixth SAS Worm Conf--Structural Analysis and Optimization (FEMCAD-89) (Edited by H. Liebowitz and G. A. O. Davies), pp. 223-231. I.I.T.T. International (1989). [8] J. T. Oden and G. F. Carey, Finite Elements, Vol. IV. Mathematical Aspects. Prentice-Hall, Englewood Cliffs, NJ (1983). [9] M. Isida, Analysis of stress intensity factors for the tension of a centrally cracked strip with stiffened edges. Engng Fracture Mech. 5, 647-665 (1973). [10] B. Gross and J. E. Srawley, Stress intensity factors for three-point bend specimens by boundary collocation. NASA Technical Note TN D-3092 (1965). [11] G. C. Sih, Handbook of Stress Intensity Factors. Institute of Fracture and Solid Mechanics, Lehigh University, Bethlehem, PA (1973). [12] R. A. Ludwig et al., Adaptive solutions of the Euler equations using finite quadtree and octree grids. Comput. Structures 30, 327 336 (1988).
Pergamon
Copyright © 1995 ElsevierScienceLtd Printed in Great Britain. All rights reserved 0013-7944/95 $9.50+ 0.00
0013-7944(94)E0056-M
A METHODOLOGY
FOR CRACK
TIP MESH DESIGN
FERNANDO C. M. MENANDRO, E. T. MOYER, JR and H. LIEBOWITZ School of Engineering and Applied Science, The George Washington University, Washington, DC 20052, U.S.A. Abstraet--A methodology for crack tip mesh design is developed which consists of comparing the mesh geometric parameters against the accuracy of the finite element solution. By successive changes in the mesh parameters a near optimal mesh can be obtained. This was done here for two-dimensional linear elastic single mode problems. The direct displacement extrapolation method for stress intensity factor estimation is used.
1. INTRODUCTION ThE FINITEELEMENTMETHOD(FEM) has been in use for about 30 years and is now established as a powerful numerical method for solution of partial differential equations. It is widely used for fracture problems, although some numerical difficulties occur. By analyzing the FEM response of two-dimensional linear elastic fracture problems, a methodology for mesh generation was developed. Significant amounts of research have been performed in automatic grid generation, as well as adaptive meshing schemes. There are, however, very few mathematical developments of these subjects for fracture problems since the theory is based on smooth functions. This led to this work's approach: to determine which geometric characteristics dominate the FEM mesh behavior and generate a near optimal mesh based on those parameters [1]. FEM work was carried out by analyzing some problems with known analytical solutions with a starting mesh and by successive changes in the mesh searching for a near optimum one. These analyses were performed with the MARC [2] software on an ALLIANT computer. The method adopted for the stress intensity factor (SIF) estimation was direct displacement extrapolation (DDE) [3]. Background in FEM analysis for linear elastic fracture mechanics (LEFM) is given only in relation to this work. The problem itself was attacked in three basic steps: (a) one-parameter studies; (b) two-parameter studies; (c) angular studies. Example problems were then tested for the desired accuracy [3]. 2. FEM FOR LEFM
The method adopted for SIF estimation in this work is the DDE method of Chan et al. [4]. The method consists of calculating the effective stress intensities at a given direction from the crack tip (usually 90 ° and 180°) and asymptotically extrapolating their values to radius zero (crack tip). Two-dimensional parabolic isoparametric quadrilateral elements were adopted, with the exception at the crack tip, where degenerated quadrilateral quarter-point elements were used. 3. ONE-PARAMETER STUDIES
In this section the accuracy of the crack tip mesh for problems in which the only operative parameter is the crack length was studied. Only two-dimensional problems in LEFM were studied. The approach chosen was to define a model problem (with known analytical solution) and study the convergence of the FEM solution due to systematic changes in the mesh. Meshes with 713
714
F . C . M . MENANDRO et al.
errors equal to or smaller than 5%, 3%, 2% and 1%, respectively, were obtained with the model problem and tested for other problems (to see if the accuracy was maintained).
3.1. Model adoption The model problem defined was the semi-infinite body with an edge crack subjected to normal stress perpendicular to the crack. The analytical values for the SIFs are shown below: K~ = 1.1215ax/(zra ) -~ 1.987807
(1)
KH = 0
(2)
for o- = 1 (applied stress) and a = 1 (crack length). The error in the analytical solution due to the finite dimensions of the model adopted (Fig. !) is of small order when compared with FEM error (i.e. for a finite dimension specimen, a power series correction factor should be introduced [4]). Since, for the dimensions adopted (width = 40, a = 1), this factor was smaller than 0.05% and the error of the FEM solution was an order of magnitude higher (~>0.1%), the model was considered to represent the semi-infinite body accurately. The results for K~ are presented as the ratio between the FEM solution and the analytical value of 1.987807 given previously. Near the crack tip, a mesh of circular rings of elements is needed. The number of rings of elements was considered as one of the parameters of the crack mesh to be studied, as well as the outer radius (r) of the last ring, as shown in Fig. 1. The third parameter was the distance (l), measured along the crack, between the crack tip and the edge of the small rectangle of refined mesh (crack tip superelement region), also shown in the figure. A rectangular superelement region was adopted since it is the most natural for straight crack problems and gives results as accurate as a circular one [6]. To use dimensionless parameters the actual parameters adopted were the number of circular rings of elements, r/l and l/a (a is the crack length). The number of elements in each ring was eight (for 180°), which have been proved to be optimal [7].
Rings of elements
Fig. 1. Finite element mesh of model problem.
A methodology for crack tip mesh design Table
1.
Kt(FEM~/giAnalytical(l/a=0.40 r/I = 0.625)
and
Table 2.
715
Kt(FEM)/Klanalytical(lla=0.40 r/l = 0.500)
Rings
180 °
90 °
Rings
180°
90 °
1 2 3 4 5 6 7 8 9 10
0.96007 0.97564 0.98225 0.99280 0.99306 0.99347 0.99349 0.99358 0.99367 0.99366
1.04953 1.04657 1.04206 1.02853 1.01372 1.01027 1.00848 1.00633 1.00403 1.00401
1 2 3 4 5 6 7 8 9 10
0.96530 0.97944 0.98526 0.99437 0.99320 0.99362 0.99361 0.99370 0.99379 0.99378
1.05197 1.04617 1.04111 1.03005 1.00975 1.00707 1.00547 1.00357 1.00148 1.00147
and
3.2. Number of circular rings tests The starting mesh selected is shown in Fig. 1 and has an edge of the superelement region over crack length (l/a) parameter equal to 0.40 and a radius of the circular mesh over edge of superelement region (r/l) parameter equal to 0.625. The results are shown for the 90 ° and 180 ° lines, for a number of rings from one to 10 (Fig. 1), in Table 1. It should be noted that the results from the 90 ° line converge to the correct solution from a higher value, while the 180 ° line gives lower bound results; this fact was used later to determine the optimal mesh. The same problem was solved with the r/l parameter equal to 0.500 and 0.750; the results are shown in Tables 2 and 3. The results along 0 = 180 ° demonstrated an invalid characteristic, namely that four circular rings of elements yield the best results. The D D E method was implemented to take 10 nodes along a radius from the crack tip, calculate the effective stress intensity factors and by linear regression extrapolate the value for the crack tip. For five circular rings of elements all the nodal displacements used were in the refined circular mesh, while for less than five there were nodes outside of this refined mesh (for one and two rings there will be less than 10 values for the linear regression in the 180 ° line, respectively six and eight) where the equations utilized no longer held. Since for most crack problems the geometry is such that the total mesh has many more elements than the superelement region alone, the change in the problem size in relation to the number of circular rings of elements is small (e.g. for this model problem approximately 1% more nodes for each additional ring). Therefore, little additional cost is introduced when using a more refined mesh. As expected, the accuracy showed improvement with more rings (and more elements), which was due to nothing more than convergence of the FEM (after the singularity was well represented). Table 4. Seven rings
Table 3. Kl(FEM)/KiAnalytica I r/l = 0.750)
(l/a
= 0.40 and
Rings
180 °
90 °
1 2 3 4 5 6 7 8 9 10
0.95509 0.97215 0.97959 0.99165 0.99285 0.99318 0.99330 0.99345 0.99356 0.99356
1.04671 1.04637 1.04224 1.02661 1.01736 1.01314 1.01109 1.00885 1.00629 1.00629
EFM 50/5-6~1
(l/a
= 0.40)
r/l
180 °
90 °
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03639 1.01996 1.01045 1.00463 1.00093 0.99853 0.99693 0.99586 0.99513 0.99465 0.99430 0.99406 0.99389 0.99377 0.99371 0.99364 0.99361 0.99357 0.99355 0.99353 0.99351 0.99349
1.03514 1.01967 1.01110 1.00623 1.00341 1.00192 1.00120 1.00100 1.00111 1.00147 1.00190 1.00241 1.00296 1.00356 1.00423 1.00483 1.00547 1.00609 1.00669 1.00731 1.00790 1.00848
716
F.C.M. Table 5. Six rings (l/a = 0.40)
M E N A N D R O et al. Table 6. Five rings (l/a = 0.40)
r/l
180 °
90 °
r/I
180°
90 "~
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03692 1.02035 1.01075 1.00485 1.00113 0.99869 0.99704 0.99595 0.99520 0.99471 0.99434 0.99410 0.99392 0.99381 0.99373 0.99367 0.99362 0.99359 0.99354 0.99353 0.99350 0.99347
1.03610 1.02059 1.01207 1.00717 1.00452 1.00308 1.00234 1.00220 1.00236 1.00276 1.00319 1.00380 1.00441 1.00508 1.00478 1.00542 1.00707 1.00677 1.00737 1.00804 1.00867 1.01027
0.100 0.125 0.150 0.175 0.200 0.225 0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625
1.03497 1.01867 1.00929 1.00353 0.99993 0.99764 0.99611 0.99509 0.99443 0.99400 0.99371 0.99351 0.99339 0.99332 0.99327 0.99322 0.99320 0.99318 0.99316 0.99313 0.99310 0.99306
1.03576 1.02059 1.01240 1.00767 1.00519 1.00403 1.00354 1.00352 1.00387 1.00442 1.00506 1.00575 1.00654 1.00734 1.00816 1.00893 1.00975 1.01056 1.01135 1.01214 1.01293 1.01372
The rate of convergence might depend on the length parameters, but it can be proved that once the singular behavior is well represented the convergence itself is not dependent on them [8]. Comparing the results for various r/l, it was noted that the results improved for smaller values of the parameter (for meshes with five to 10 rings). This was understood as a need for a refined mesh region near the crack tip, but, as was discovered later, there is an optimum value beyond which the refined mesh is no longer sufficient to represent the asymptotic behavior. Comparing the results from the three tables, it was observed that for the same number of rings the accuracy increased for smaller values of the r/l parameter. This is shown in the next section. 3.3.
Tests for r/1
From the previous discussion it was decided to test the problem for small changes in the parameter r/l (Tables 4-6). Refinement of more than seven rings of elements was considered unnecessary and was not pursued. It is important to note that the 90 ° line showed what seemed to be a stationary point (minimum) where the solution was optimal, while for the 180 ° line the convergence was not clear (there was no stationary point in the function of r/l). The finding of a minimum was in agreement with the observations in the previous section, i.e. that the 90 ° line gives upper bound solutions. In the three tables the best result for the 90 ° line was at r/l = 0.275 (which gives a minimum) and for the 180 ° line at r/l = 0.200. For the 180 ° line, a maximum should be expected and could not be found; thus, the best result considered was the one which best approximated the analytical value. 3.4.
Tests for l/a
The
Continuing the search for an optimal solution, the l/a parameter was now analyzed (Table 7). r/l parameter was kept constant along with the number of circular rings of elements. From Table 7. r/I = 0.500
l/a
180 °
90 °
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.40
1.00219 0.99993 0.99825 0.99696 0.99598 0.99522 0.99464 0.99422 0.99361
1.00788 1.00637 1.00545 1.00486 1.00457 1.00448 1.00456 1.00481 1.00547
Table 8
l/a
r/l
180 ~
90"
0.24 0.24 0.32 0.32
0.200 0.275 0.200 0.275
1.00721 1.00213 1.00253 0.99745
1.00632 1.00324 1.00339 1.00067
A methodology for crack tip mesh design
717
Table 10. r/l = 0.275
Table 9. l/a = 0.32
r/l
180 °
90 °
I/a
18if'
90 ~'
0.250 0.275 0.300
0.99852 0.99745 0.99673
1.00092 1.00067 1.00068
0.30 0.32 0.34 0.36 0.38 0.40
0.99821 0.99745 0.99691 0.99646 0.99613 0.99586
1.00093 1.00067 1.00067 1.00065 1.00080 1.00100
then on, the tests were performed only for seven rings, which had shown good accuracy and convergence. Again, a stationary point was found for the 90 ° line, while the 180 ° line did not show the same behavior. The most accurate solutions were for l/a equal to 0.24 and 0.32 for the 180 ° and the 90 ° lines, respectively. 3.5. Combined tests Combining the best results for both parameters, the desired mesh should be obtained (Table 8). Since the most accurate mesh parameters might not be independent, it was important to test the neighborhood of this solution (by separate small changes in each of the parameters). Tables 9-11 show these test results. Since the 90 ° line had always given upper bound solutions and the 180 ° line had not shown the same consistency (the results seemed to be more affected by changes in the mesh parameters, sometimes diverging), the mesh with l/a = 0.32 and r/l = 0.275 was considered to be the most accurate (Table 8). The mesh selected in the previous paragraph was obtained by combining independently studied results for r/l a n d / / a . Since these parameters might have some correlation, there was no guarantee that this mesh was really the optimal one. To address this problem, the SIF for the 90 ° line was analyzed as a function of r/l and l/a and the neighborhood of the selected mesh was tested (searching for a minimum). As expected, the parameters were interrelated. The best result was actually for l/a = 0.36. As mentioned in Section 3.2, the FEM solution converges after the singular behavior is represented: this led to the use of a fixed number of circular rings of elements (seven) for the rest of the tests. Once the optimal length parameters have been found, the convergence due to the number of elements (or the number of rings) should be studied in order to define meshes for the desired accuracies. This is shown in Table 12. Calculating the percentage difference between the analytical solution and the FEM (90 ° line) asymptotic displacement solution, the error was 3.188% for four rings and 0.300% for five rings. This drop in accuracy had to do with the D D E method as mentioned before, and gave a warning about using less than five circular rings of elements. The most accurate acceptable mesh for an edge crack was thus the one with the following parameters: I - = 0.36 a r = 0.275
and at least five rings of elements. Table 12. I/a = 0.36, r/I = 0.275
Table 11. l/a = 0.36
r/1
180 '~
90 ~
0.250 0.275 0.300
0.99752 0.99646 0.99573
1.00088 1.00065 1.00073
Rings
180'
90 ~
1 2 3 4 5 6 7
0.97714 0.98858 0.99308 1.00011 0.99568 0.99656 0.99646
1.05441 1.04419 1.03871 1.03188 1.00300 1.00181 1.00065
718
F.C.M. MENANDRO et al.
Since the five-ring mesh had already given results with better than 1% accuracy, this was selected as the desired mesh (for 10%, 5%, 3%, 2% and 1%), which should be tested for some other problems. 3.6. Final tests Once the desired mesh had been obtained, testing for other one-parameter problems could either support the chosen mesh or disregard it, suggesting changes. The problems used to test the chosen mesh were the well known centrally cracked strip subjected to tension and the three-point bend specimen. For the centrally cracked strip (Fig. 2), Isida [9] developed (based on Laurent series expansions of the complex stress potentials) a solution which for our purposes could be considered exact. The mode I and II SIFs are given by: Kt = a x / ( ~ a ) F ( 2 )
(3)
Kn = 0,
(4)
where tr is the applied stress, a is the half crack length, 2 = 2a/W, W is the plate width and F(2) is a 36-term power series of 2 given in ref. [9]. For the three-point bend specimen (Fig. 3), a solution considered to be accurate up to 0.5% was obtained by Gross and Srawley in 1965 [10]. This solution was used only for comparison of results, since the reported accuracy was below this study's requirements. The SIFs, as presented by Sih [11], are:
K~ = 6_.~h,~ __M./-aa F(2)
(5)
KH = 0,
(6)
where M is the applied moment, a is the crack length, t is the thickness of the beam, h is the height of the beam, 2 = a/h and F(2) is a five-term power series of 2 given in ref. [11]. The results of the tests, already divided by the adopted solution, are shown in Table 13. The error was approximately 1% for every problem, which indicated that the mesh of the previous section might not be as accurate as expected. Some values greater than 1% were found where the expected errors should be less than 1%. All this evidence led to the application of the same procedure used in Section 3.5 to search for a better mesh.
tttttttttttt
v
b
J ' 2 . J'
b --
W
~
P
s T
Fig. 2. Centrallycracked strip.
Fig. 3. Three-pointbend specimen.
719
A methodology for crack tip mesh design Table 14. W = 5a, gt(FEMl/gtAnalyticaI
Table 13. KI~FEM)/KIAaalyti~I Model Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Centrally cracked strip, Three-point bending
W ,> 2a W/2a = W/2a = W/2a = W/2a =
5 4 3 2
180 °
90 °
l/a
r/l
180 °
90 °
0.99031 0.99044 0.99040 0.99012 0.98817 0.98686
! .01159 1.01192 1.01201 1.01206 1.01141 0.98879
0.38 0.38 0.38 0.38
0.275 0.250 0.225 0.200
0.99002 0.99091 0.99233 0.99449
1.01254 1.01163 1.01122 1.01143
0.36
0.300
0.98986
1.01302
0.36
0.275
0.99044
1.01192
0.36 0.36 0.36
0.250 0.225 0.200
0.99136 0.99273 0.99493
1.01114 1.01066 1.01101
0.34 0.34 0.34 0.34 0.34
0.300 0.275 0.250 0.225 0.200
0.99038 0.99098 0.99189 0.99329 0.99545
1.01242 1.01143 1.01066 1.01031 1.01062
0.32 0.32 0.32 0.32 0.32
0.300 0.275 0.250 0.225 0.200
0.99107 0.99165 0.99255 0.99397 0.99613
1.01197 1.01098 1.01025 1.01004 1.01043
0.30 0.30 0.30 0.30 0.30
0.300 0.275 0.250 0.225 0.200
0.99197 0.99253 0.99343 0.99481 0.99703
1.01177 1.01079 1.01013 1.00983 1.01044
0.28 0.28 0.28 0.28 0.28
0.300 0.275 0.250 0.225 0.200
0.99305 0.99363 0.99453 0.99593 0.99810
1.01163 1.01076 1.01013 1.00997 1.01048
0.26 0.26 0.26 0.26 0.26
0.300 0.275 0.250 0.225 0.200
0.99452 0.99508 0.99598 0.99737 0.99956
1.01194 1.01108 1.01052 1.01040 1.01103
Again, the neighborhood of the solution was tested, this time for the centrally cracked strip problem. The width over length (W/2a) ratio adopted was 5. The results are presented in Table 14: the result in italics was the starting mesh and the final result is given in bold. These results were obtained in the same way as in section 3.5, the final result being selected as a minimum for the 90 ° line. The error for the final result was less than 1% for both directions (90 ° and 180°). Again, the results for the 180 ° line were not considered due to the already discussed argument. To finish this analysis it was necessary to utilize this final mesh on the previous problems in order to check their accuracy. The results, divided by the adopted theoretical solution, were as follows: Edge crack (model problem) 180°:
900:
KI(FEM)=0.99988 KI(Th.)
KItFEM)=1.00359 KI(Th.)
Cracked infinite plate (modeled as strip with W = 108a) 180°:
KI~FEM)= 0.99472 KI(Th.)
90o:
KI(FEM) -- 1 . 0 0 9 5 7 gl(Th.)
720
F.C.M.
M E N A N D R O et al. 12-5-1990 m MENTAT
ORTHO V.P.O.O00e+O00.O00e+00
0.100e+O1
X J
Fig. 4. Crack tip mesh.
Three-point bend specimen 180°:
KI(FEM)=
0.99078
gl(Th.)
900: KI~FEM)=0.98521. gl(Th.)
The three-point bend problem continued to show up the error detected previously, but only for one direction. If the SIF was calculated as the mean of the results in both directions, the accuracy of the solution would improve. The explanation of this discrepancy is probably the lack of precision of the adopted theoretical solution, as mentioned previously. F o r the other two problems, the error in the solution was, as expected, less than 1% and therefore this mesh was adopted as optimal. The result for the edge crack seemed to be less precise than the one obtained in the previous section, but as mentioned earlier the second decimal digit was beyond this study's precision. Thus, the optimal mesh for problems in which the only operative parameter is the crack length is: l - = 0.30 a r
-! = 0.225 and at least five rings of elements. This mesh is shown in Fig. 4. The next section studies this mesh for problems in which more than one length parameter must be considered to define the problem geometry. 4. T W O - P A R A M E T E R S T U D I E S
Two-parameter problems constitute ones in which the crack length is no longer the only length parameter that determines the required crack tip mesh. This occurs when the other dimensions of the specimen are of the same order of magnitude as the crack length. Rules that take into account geometric parameters other than the crack length were developed for mapping the crack tip superelement mesh.
A methodology for crack tip mesh design
721
4.1. Model definition The adopted model problem was again the centrally cracked strip (Fig. 2) of Section 3.6, now with different values for the width over crack length (W/2a) parameter. The theoretical solution for this problem is repeated below for modes I and II: KI = a ~/(na )F(2)
(7)
K,, = O,
(8)
where e is the applied stress, a is the half crack length, 2 = 2a/W, W is the strip width and F(2) is a 36-term power series of 2 given in ref. [9]. 4.2. Tests Making the strip width (W) in a centrally cracked strip (Fig. 2) smaller in relation to the crack length (2a) increases the contribution of ligament 6 and, at some point, it starts to dominate. It should be expected that the smallest parameter dominates, which was confirmed for the tests run for g varying from four times the half crack length to one-quarter of it. The same mesh defined in the previous section was adopted for these problems with the difference that when 6 was smaller than a the length parameter used for the mapping was A These results divided by the respective analytical solutions are shown in Table 15. The results confirmed the suspected behavior and it can be seen that all were better than 1%. In the following section, another problem that occurred in two-parameter problems is discussed. 4.3. Tests for a _~ d A real two-parameter problem appears when the half crack length and the ligament have similar values, more specifically when 0,5 < g- < 2,
(9)
a
where 6 is the ligament and a is the half crack length (as shown in Fig. 2). When mapping the crack tip mesh for problems with 6/a agreeing with the inequality (9), some difficulties appeared in the transition from the superelement region to the regular domain mesh. This problem was studied using four different meshes, as shown in Fig. 5. All four meshes had the domain discretized the same way, with element sides roughly equal to the crack size. The transition between the superelement region and the regular domain elements was accomplished in four different ways from the following rules. To check the assumption that the dominant parameter was the smallest, one of the meshes used the other parameter for mapping. The transition adopted in this case seemed to be the most reasonable one (case 1). Using the smallest parameter for mapping created some ambiguity on what transition mesh should be adopted, the first option (case 2) being logically to use only four elements around the superelement to cover the entire extension (a + 6). If two elements are used to make the transition to the bigger side, one might opt to smooth these elements, or have the outer one rectangular (cases 3 and 4, respectively). By smoothing it is meant that the nodes are moved so that the distance between them is roughly the same [12]. The results of these analyses are shown in Table 16. The results for case 4 were usually the best, but since no significant difference existed between cases 3 and 4 the most straightforward was accepted, which was case 3. Table 15. Two parameters: for {~ >~a < a
I/a=0.30 I/3 = 0.30
3/a
180°
90 °
4 3 2 1 1/2 1/3 1/4
0.99481 0.99475 0.99442 0.99225 0.99064 0.99028 0.99009
1.00983 1.00988 1.00982 1.00870 0.99966 0.99663 0.99507
F . C . M . MENANDRO et al.
722
\ Case 1
Case 2
:
i
i
:
Case 3
Case 4 Fig. 5. Transition mesh cases.
M o r e tests s h o w e d t h a t all m e s h e s p e r f o r m well f o r a 10% difference b e t w e e n a a n d ~. In o t h e r w o r d s , f o r ~/a b e t w e e n 0.91 a n d 1.10, all m e s h e s g a v e a c c u r a c y to w i t h i n 1%. C a s e s 3 a n d 4 a d o p t the c r a c k l e n g t h as t h e m a p p i n g p a r a m e t e r , cases 3' a n d 4' the l i g a m e n t . T h e s e results a r e s h o w n in T a b l e 17. T h e p e r f o r m a n c e o f cases 3 a n d 4 for ~ = a was n o t t h a t desired d u e to the i n f l u e n c e o f a s p e c t ratios. It was s u g g e s t e d , as a g e n e r a l rule, t h a t f o r 0.91 < ~/a < 1.10 case 2 s h o u l d be a d o p t e d , case 3 r e m a i n i n g f o r t h e rest. A s a n o t h e r g e n e r a l rule it was n e c e s s a r y t h a t n o n e o f the t r a n s i t i o n e l e m e n t s ' side l e n g t h s be g r e a t e r t h a n t h e c r a c k (or l i g a m e n t ) length.
Table 16. 0.5 < ~/a < 2
d/a
Case
2/3
1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
3/4
4/3
3/2
180° 0.98633 1.00490 0.99646 0.99648 0.98946 0.99914 0.99460 0.99444 0.98996 0.99469 0.99133 0.99143 0.98937 0.99588 0.99247 0.99261
90 ° 1.00378 1.01721 1.00803 1.00799 1.00558 1.01238 1.00752 1.00730 1.01098 1.01069 1.00691 1.00701 1.01263 1.01174 1.00793 1.00807
Table 17. a ~
8/a 1.00
0.91
1.10
Case 2 3 3' 4 4' 2 3 4 2 3 4
180° 0.99225 0.98853 0.98986 0.98821 0.98882 0.99374 0.99016 0.98959 0.99310 0.99043 0.99031
90 ° 1.00870 1.00438 1.00661 1.00400 1.00550 1.00898 1.00552 1.00490 1.00938 1.00616 1.00601
A methodology for crack tip mesh design 9-12-1990 ca MENTAT
9-12-1990 oa MENTAT
ORTHO V.P. 0.000e+00 0.000e+00 0.100e+01
X )RTttO V.P. 0.000e+00 0.000e+00 0.100e+01
9-12-1990 oa MENTAT
ORTHO V.P. 0.000e+00 0.000e+00 0.100e+01
723
9-12-1990 ¢o MENTAT
x OR"~HO V.P. 0.000e+00 0.000e+00 0.100e+01
Fig. 6. Crack tip superelement regions. 5. A N G U L A R S T U D I E S Angular studies are the studies o f the influence o f angular refinement a r o u n d the crack tip on the accuracy. The accuracy o f the adopted mesh with eight elements a r o u n d the half circle was tested against six or 10 elements a r o u n d the half circle. The accuracy and the displacement fields were examined and compared. The choice o f eight elements was justified in view of the other options.
5.1. Model adoption In this section the refinement a r o u n d the crack tip in the angular direction was studied. The adopted solution o f eight elements a r o u n d the semi-circle was c o m p a r e d with options o f six and 10. The problem analyzed was the centrally cracked strip, as this is a well studied reference solution. The accuracy o f the solution, as well as the displacement fields, was studied to establish the optimal solution. The crack tip superelement regions are shown in Fig. 6. Meshes 6', and 10 and 10' were proposed and mesh 6 was obtained from 6' by a smoothing algorithm used at the non-circular region. The solutions o f the centrally cracked strip divided by the analytical one are shown in Tables 18-21. The non-dimensional parameters were taken close to the ones o f the previous analyses. It can be seen that although mesh 6 had a better behavior than 6', the same cannot be guaranteed for 10 elements. To improve the obtained solutions the accuracy against r/l was studied for cases 6, 10 and 10' (6' was obviously worse). Table 19. Mesh 6 (lla = 0.300) Table 18. Circular refinement test (l/a = 0.300; r/l = 0.200) Mesh 180° 90° 6' 0 . 9 8 7 7 9 0.99252 6 0 . 9 9 1 6 0 0.99561 10' 0 . 9 9 4 6 9 1.01659 10 0 . 9 9 3 6 9 1.01512
r/l 0.250 0.225 0.200 0.175 0.150 0.125 0.100
180° 0.98839 0.98964 0.99160 0.99473 0.99978 1.00819 1.02290
90° 0.99662 0.99584 0.99561 0.99635 0.99876 1.00420 1.01551
724
F.C.M.
M E N A N D R O et al.
Table 20. Mesh I0 (l/a = 0.300)
Table 21. Mesh 10' (l/a = 0.300)
r/I
180 °
90 °
r/l
180 °
90 °
0.250 0.225 0.200 0.175 0.150
0.99260 0.99299 0.99369 0.99501 0.99750
1.01697 1.01595 1.01512 1.01487 1.01581
0.250 0.225 0.200 0.175 0.150
0.99370 0.99405 0.99469 0.99592 0.99823
1.01849 1.01747 1.01659 1.01626 1.01704
The result for r/l equal to 0.150 was the best and, although a complete analysis was not performed, it would probably end up as the optimal. The best solution for mesh 10 was the one with r/l = 0.175 (it was a minimum for 90°). Unfortunately, it still did not give less than 1% in both directions. Originally, mesh 10' had the same r/l as mesh 10 and a similar inaccuracy in the results. For clarity, the best solutions for all meshes analyzed are shown together in Table 22. Undoubtedly the two best solutions are the ones for six and eight elements around the semi-circle. These meshes were thus compared thoroughly. Eight or six elements around the semi-circle were then utilized to compare displacement fields for the centrally cracked strip and the edge crack. The results of ux and uy are shown against 0 in Figs 7-10. The comparisons were taken at the outer border of the first and second rings of elements, which do not lie in the same distance for the different meshes. For the superposed charts, the solutions were scaled so that they were equal. It is worth remembering that the displacements in the x-direction are off by a constant and the zero was arbitrarily taken at 0 equal to 180 °. From Figs 7 and 8, one might assume that mesh 6 is better; however, since for the edge crack the pattern does not hold and references give eight elements as an optimum, it will stay as the selected mesh. Another reason for this decision is that the SIF results for the edge crack were slightly better and mesh 6 did not seem to present the 90 ° and 180 ° lines as upper and lower bounds, respectively, as experienced previously. The SIFs divided by the reference solution are shown in Table 23. For the purpose of this research, eight elements around the semi-circle were adopted as being precise enough. More research in this direction (0, angular) is advised before anything is concluded in problems that involve more than just LEFM. 6. CONCLUSIONS A methodology for crack tip mesh design was described which consists in determining what parameters influence the mesh behavior. By extensive analysis of FEM solutions the accuracy is checked against the mesh parameters. A near optimal crack tip mesh was obtained for single mode
o Ux - 8 elements
~,,,~¢,~,m---~
m Uy - 8 elements
t~
Uy/¢"
o Ux - 6 elements LxUy - 6 elements
Table 22. Angular comparison Mesh
180 °
90 °
8 6 10 10'
0.99481 0.99978 0.99501 0.99592
1.00983 0.99876 1.01487 1.01626
~/
,I
~ r
.o
Ux
:~..~'9~
~
,
.~ "~ °
Z Table 23. Angular test
t , ~ A'"
Mesh
Problem
180 °
90 °
6 8 6 8
Strip Strip Edge Edge
0.99978 0.99481 1.00584 0.99988
0.99876 1.00983 0.99618 1.00359
0
, 30
~ 60
J 90 Angle
~ 120
~ 150
\ N~ 180
Fig. 7. Displacements for centrally cracked strip. Six elements around crack tip, radius 0.009. Eight elements around crack tip, radius 0.0135.
A methodology for crack tip mesh design
725 • ~/)'m~'i
o Ux - 8 elements e~
Uy
-
y ~
8 elements
Uy ~(,/,f/
o Ux - 6 elements
~:
t, Uy - 6 elements
j
J
eO
O
z
0
30
60
90 Angle
120
150
180
Fig. 8. Displacements for centrally cracked strip. Six elements around crack tip, radius 0.018. Eight elements around crack tip, radius 0.027.
t3 Ux - 8 elements =
• Uy - 8 elements
St//
.z./
o Ux - 6 elements zx Uy - 6 elements
--
e~ @
0
30
60
90 Angle
120
150
180
Fig. 9. Displacements for edge crack specimen. Six elements around crack tip, radius 0.009. Eight elements around crack tip, radius 0.0135.
uy:J'~"~ J /
n Ux - 8 elements m Uy - 8 elements o Ux - 6 elements
8
~,
• Uy - 6 elements
/
i¢, Z 13 \
O
30
60
90
120
150
180
Angle Fig. 10. Displacements for edge crack specimen. Six elements around crack tip, radius 0.019. Eight elements around crack tip, radius 0.027.
726
F . C . M . MENANDRO et al.
c r a c k p r o b l e m s . It was s h o w n t h a t a l t h o u g h the s o l u t i o n s for the D D E m e t h o d are b e t t e r for the crack face (180°), the 90 ° d i r e c t i o n gives a b e t t e r p a r a m e t e r for analysis. It was also s h o w n t h a t for single m o d e L E F M p r o b l e m s o n l y the crack a n d l i g a m e n t (smallest d i s t a n c e b e t w e e n crack tip a n d the free edge) l e n g t h s i n f l u e n c e the F E M m e s h b e h a v i o r .
REFERENCES [I] F. C. M. Menandro, A knowledge based approach for crack tip mesh generation. M.Sc. Thesis, The George Washington University, Washington, DC (1991). [2] MARC general purpose finite element program, Marc Software International, Revision K3 (1988). [3] F. C. M. Menandro, E. T. Moyer, Jr and H. Liebowitz, A near optimum crack tip mesh. Engng Fracture Mech. 50, 703-711 (1995). [4] S. K. Chan, I. S. Tuba and W. K. Wilson, On the finite element method in linear fracture mechanics. Engng Fracture Mech. 2, 1-17 (1970). [5] H. L. Ewald and R. J. H. Wanhill, Fracture Mechanics. Edward Arnold and Delftse Uitgevers Maatschappij (1984). [6] E. T. Moyer, Jr and H. Liebowitz, Creep and fracture characteristics of materials and structures at elevated temperatures. ONR Final Report, no. N00014-84-K-00027 (1986). [7] E. T. Moyer, Jr and H. Liebowitz, Finite element mesh design--a knowledge based approach. Proc. Sixth SAS Worm Conf--Structural Analysis and Optimization (FEMCAD-89) (Edited by H. Liebowitz and G. A. O. Davies), pp. 223-231. I.I.T.T. International (1989). [8] J. T. Oden and G. F. Carey, Finite Elements, Vol. IV. Mathematical Aspects. Prentice-Hall, Englewood Cliffs, NJ (1983). [9] M. Isida, Analysis of stress intensity factors for the tension of a centrally cracked strip with stiffened edges. Engng Fracture Mech. 5, 647-665 (1973). [10] B. Gross and J. E. Srawley, Stress intensity factors for three-point bend specimens by boundary collocation. NASA Technical Note TN D-3092 (1965). [11] G. C. Sih, Handbook of Stress Intensity Factors. Institute of Fracture and Solid Mechanics, Lehigh University, Bethlehem, PA (1973). [12] R. A. Ludwig et al., Adaptive solutions of the Euler equations using finite quadtree and octree grids. Comput. Structures 30, 327 336 (1988).