- Email: [email protected]
A flexible strategy for numerical path integration
PII:
ELSEVIER
Advances in Engineering Software 28 (1997) 203-210 0 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 0965-9978/97/$17.00
SO965-9978(96)00047-6
Technical Note A flexible strategy for numerical path integration Gang Xue, Xucheng Wang & Qinghua Du Department
of Engineering Mechanics, Tsinghua University, Beijing 100084, P.R. China
(IReceived 11 April 1995; received in revised form 24 July 1996; accepted 10 October 1996) In this paper, a flexible strategy for numerical path integration is proposed using the J integral as an example to describe the algorithm that can be extended to the evaluation of other J-like integrals without difficulty. A program called PathInt is developed to run in expert mode or automatic mode. When in expert mode, the program may be controlled in a more detailed way to meet the needsof proficient researchers. When in automatic mode, the program can determine several parameters automatically according to the default settings. A numerical example of a typical elastic fracture problem demonstrates the effectiveness of the algorithm. The applicability of several finite element meshes with or without special crack-tip elements is also investigated. It is found that eight-node collapsed elements are of great potential usefulness since the results from the meshesemploying them at the crack tip are satisfactory in both elastic and plastic problems. 0 1997 Elsevier Science Limited. All rights reserved. Key
words:
numerical path integration, J-like integral, expert mode, automatic
1 INTRODUCTION
when it is applied to the practical problem. The limitation is mainly on the nodal numbering thus also
Since Rice’s’ introduction of the path-independent J integral, many other similar path integrals have been proposed as fracture characteristic parameters, such as Jv integral2 for viscoelastic media, i integral3 for incremental plasticity media, C(tj4 and Ct5 integrals for transient creep media both of which can degenerate into the C* integral6 for steady-state creep media. In most practical fracture problems it is too complex to find the analytical solutions, so numerical techniques, particularly focused on the finite-element method (FEM), play an important role in the solution procedure. Nowadays, there are many general purpose programs for finite-element analysis from which only the solutions of displacement and stress fields can be obtained in most cases. Thus, a posit-processor is needed to evaluate the path integrals deiined above. Concerning numerical methods for fracture mechanics, Owen & Fawkes7 have developed a FORTRAN subroutine named as JevaM to evaluate J integral, but there are some severe restrictions
on the local coordinate system of the element. Nowadays there are so many structured or unstructured, semiautomatic or automatic and interactive mesh generation programs by which it is inconvenient for the user to know exactly the local numbering sequence of the nodes in each element. In this paper, a more flexible strategy is proposed for numerical path integration on the basis of JevaI6. More
general considerations are discussedfor path integration in an individual element. According to our new program named PathInt,
developed with extended FORTRAN
77 language, only the numbering of the elements along the integral path is needed as input data while the program is running in automatic mode for the numerical evaluation of J-like path integrals without thinking about the axis directions of the local coordinate system. In addition, the program has the capability for proficient researchers to control the integration path in detail while the program is running in expert mode. To describe the algorithm, we use, as an example, the 203
G. Xue, X. Wang, Q. Du
(4
one of the local coordinate directions of the element. With respect to this, Bakker* gives the following two important reasons:
rl 5
2
1
4 6
y 8
0
4 ~ 3
t 4
7 I
4
Gaussian
0
Element
Quadrstare
Points
Nodes
04
-5
+
Gaussian
0
Element
Quadrature
Points
Nodes
Fig. 1. (a) Non-corner element on the integration contour.
(b) Corner elementon the integration contour. J integral on contour l? defined as’
where U is the strain energy density, ti is the traction vector and ui is the displacement vector. The numerical example of a typical elastic fracture problem, the single-notched plate under a pin load, is employed to demonstrate the effectiveness and ease of use of the algorithm. The applicability of several finiteelement meshes with or without special crack-tip elements is also investigated by the example; it is found that eight-node collapsed elements are of great potential usefulness since the results from the meshes applying them at the crack tip are satisfactory in both elastic and plastic problems.
(1) A finite-element analysis program normally only produces the stress at the integration points that are considered optimal stress points.’ The path chosen above can avoid the need for inte~olation procedures which usually induce some errors. (2) In elastic-plastic or other nonlinear analyses, the integration points are the only locations where material behavior coincides with the given constitutive law exactly. Below, we discuss in detail the general considerations for an element on the integration path. These considerations relate to the shape, the location and the direction of the path, and are general but not too complicated to handle. The elements intersecting the integration contour can be divided into two categories as non-corner elements and corner elements. In the non-corner element shown in Fig. l(a), the integration path is along one of the local coordinate axes without turning. In the corner element shown in Fig. l(b), the path changes its direction, namely there exists at least one corner in it. There are two sets of non-corner elements in which the integration path is along the n direction with < constant or along the 5 direction with v constant, respectively. In fact, there may exist several corners in a corner element especially for high Gaussian quadrature order. For simplicity, it is assumed that the path is chosen so that at most only one corner exists in the element. Thus there are only two types of corner element in which the integration path is changed from the q direction to the t direction or vice versa. So altogether there are four kinds of element intersecting the integration contour. In the program, an integer variable ikind can take different values according to the respective situations. Thus ’ 1 non-corner element in which path is along the n direction 2 non-corner element in which path is along the 6 direction ikind = 4
2 GENERAL CONSIDERATION FOR PATH INTEGRATION IN AN ELEMENT In computational fracture problems, quadrilateral elements are widely used. Since triangular elements can be considered as degenerated quadrilateral elements, we only deal with quadrilateral elements in this paper. For a path-independ~t integral such as the J integral, the integration contour would be chosen freely in principle. But for the finite-element method, the path is restricted just through element integration points (assuming with NG x NG Gaussian quadrature order), and parallel with
\
3 corner element in which path changes from the Q to the < direction 4 corner element in which path changes from the 5 to the q direction
Furthermore, the location and the direction of the path in an element must be thought about. For each pattern of non-corner element, say, for the case of ikind = 1, there are NG types of integration path with t held are the index and constant at
Numerical path integration
205
and ipeta are the indices of the Gaussian quadrature points so that < held constant as <$xsi on the half-path along the 77direction, and q held constant as qipetaon the other half-path along the E direction, satisfying that 1 < ipxsi, ipeta 5 NG. It is similar to this when ikind = 4. Namely, ipxsi
-5
iab = ideta x 10 + idxsi
and icd = ipeta x 10 + ipxsi @
Gaussian
@
Element
Quadrature
Points
Nodes
Fig. 2. Gaussian quadrature points in an eight-mode quad-
rilateral elementwhen NC = 2. may be 1,2 . . IVG, satisfying c1 < & < . . . < &, namely, the points are ordered according to their
and icd = ipxsi x lOi+ @eta
where 1 along the negative direction of 77-axis (with E held constant)
idxsi =
I ideta =
2 along the positive direction of q-axis (with 5 held constant) 1 along the negative direction of E-axis (with 7 held constant) 2 along the positive direction of {-axis (with 17held constant)
To summarize, two integer variables, ikind and icode, are employed to represent all the information about the path in an individual element. The representation is easy to handle and capable of describing a rather complex situation. For example, for the integration path shown in Fig. l(a) it is that ikind = 1 and icode = 2, and for the case shown in Fig. l(b) it is that ikind = 3 and icode = 2122. The path integration is assembled from all the contributions of each individual element on the contour. The integration in the non-corner element is simple. For example, in the case of integrating along the 77 direction,7 the contribution from an individual element is Jqelem =
where U, and UP are the elastic and plastic work, IAand w are displacements, (T,, aY, and 7XVare stresses.The expression of Jpelemthat corresponds to the case along the < direction is similar. It is a little complicated for the corner element. As shown above, the path in the corner element may be divided into two half-paths. One is along the Qdirection and the other along the < direction. It can be treated in an approximate way.7 As shown in Fig. 3, it is scaled so that J elem &dem+~Jy-m L t
where L is the total range of the integration along the local axis, L = 2 if the local coordinate arranges from - 1 to + 1, and Lf or L, is the actual range corresponding to the direction along the c or Qaxis such that LE and L, are obtained from
and L, = Irlr - ‘%ideI
in which (Et,Q) is the local coordinate vector of the turning point, and [sideor vsidecorresponds to the local coordinate which is held constant on the element side on which the path gets into or gets out of the element. The
206
G. he,
X. Wang, Q. Du
ilocal =2
ie&=+1
I
l l
Gaussian
Quadrature
Element
Nodes
Point
i&l =2 iedge=-l
Fig. 3. The scalingprocessin a corner element.
[sideand qside,- 1 or 1, depend on the pattern of the corner element and the directions of the two halfpaths. In the case shown in Fig. 3,
VdLleS Of
‘%ide
=
-1.
Obviously, some numerical inaccuracy is introduced by this scaling process, but since the contribution of such corner elements to the complete J integral will be small, the error involved should be acceptable, which can be verified by the numerical example shown later. From the above, most integration path situations in an element can be defined by two integer variables. This can meet the demands of the scientist who intends to control the program in a wide range. When the program is running under this kind of thorough control by the user, it is said to be in expert mode. But not every researcher needs to control the path integration completely in this way. In this case, it is convenient to set some default values such that only the most necessary data need to be prepared by the user when the program is running in automatic mode.
3 AUTOMATIC
MODE
In the expert mode described above, the user should know exactly the local information of the element and have complete control to define the path in each individual element intersecting the integration contour in sequence. Hand calculation may be inconvenient, sometimes even nearly impossible, for most practical problems whose finite-element meshesare composed of a large number of nodes and elements, particularly for those whose meshes are generated automatically by some mesh generator. Next, the algorithm determines ikind and icode for each element automatically based on the sequenceof elements on the path when the program is running in automatic mode. Then the user only need to prepare a few necessarydata. A subroutine named GetCode is developed to determine the path-control parameters for each element. It is assumed that the numbering of the previous, current and next elements on the contour are expressed
Fig. 4. The values of ilocal and iedge on each side of the
element. by three integer variables ieprev, iecur, ienext, respectively. By comparing the nodal numberings of the two neighboring elements ieprev and iecur, it can be found which side of them is common. This common side in element iecur may be expressedby two integer variables ilocal and iedge whose values are defined in Fig. 4. Similarly, the common side of elements iecur and ienext can be expressed by another two integer variables jlocal and jedge. If element iecur is the first or last element on the contour, the common side is on the boundary and the values of these integer variables can also be obtained, which will not influence the algorithm afterwards. If ilocal is equal to jZoca1, it is a non-corner element, otherwise a corner element. According to the values of the four variables defined above, the pattern and the direction of the integration path in the element can be determined, but the location, namely the value(s) of ipxsi or/and ipeta are still unknown. Since the consecutive intersections of the contour crossing with the elements are continuous, the starting point of the path in element iecur, at which point the contour intersects the common side of elements iecur and ieprev, is the sameas the endpoint of the path in element ieprev. In the noncorner element, the endpoint can be obtained easily if the starting point of this element is known. In the corner element, with the assumption that the path gets into and out of the element in a symmetric manner, namely LE = L, just as the situation shown in Fig. 3, the endpoint can be obtained too if the starting point in this element is known. Thus the location of the path in each element can be determined once a starting point of the contour is defined. So from a starting point and the numbering sequence of the element series on the contour, the ikind and icode of each element can be obtained. In the program, it is important to treat some particular elements carefully, such as degenerate and collapsed elements. 4 STRUCTURE
OF THE PROGRAM
Since the efficiency of the path integration program is
207
Numerical path integration
Lr’ main
dimea
calculatesvariables
associated
with dynamic dimensioning.
input reads input data that includes definition and definition of integration paths.
of mesh, results of finite element analysis,
SetCaussConst
t
‘E
evaluates
ikind and
icode for an element
BoundEdge -
judges whether the element is on structure sides on boundary if there are.
boundary,
and determines
-
judges whether two elements are adjacent, and determines common side if there is.
the
AdEdge which is the
GetGpInt gets the index of the Gaussian point on the path
GetIntersectCoord -
gets the coordinates of the intersection neighboring elements.
Outputs the coordinates
on the common
side of two
of the Gaussian points on the path.
Fig. 5. The program structure.
not so crucial as that for finite-element analysis, the program PathInt j.s developed inclined to its understandability and expandability. The program is designed modularly, consisting of about 20 modules. The structure of the program is shown schematically in Fig. 5. The two integer variables ikind and icode are the most important co:ntrol data of the program on which the function of the major subroutines are concentrated. According to whether they are input by the user or not, the program, is classified to be running in expert mode or automatic mode. The automatic mode is recommended because of its easy use. For expert mode, the subroutine of CheckCode is necessary to verify the legality of the two control data prepared by the user. In any mode, the subroutine WrPoint is called to output the coordinates of the Gaussian quadrature points on the contour so as to display the contour graphically further.
The program is written in extended FORTRAN 77 language of about 2000 source lines of which about 3040% are comments. The following features are extended based on standard FORTRAN 77 to ensure correctness, understandability and expandability: (1) lowercase letters for keywords and variable names; (2) names longer than six characters; (3) use of ‘!’ to delimit comments; (4) the constructs of ‘do. end do’ and ‘do while. . . end do’. Fortunately, many modern FORTRAN compilers will accept all of these deviations from the standard. 5 NUMERICAL EXAMPLE To demonstrate the algorithm developed in this paper, a single-edge cracked plate is considered under tension of pin loading P as a linear elastic fracture problem. The
208
G. he, X. Wang, Q. Du
(4
(4
t
-crack
crack tip 2
---
tip
integration contours element edges
Fig. 6. Finite element mesh (if: (a) overview, (b) eight integration contours surrounding the notch tip orderedfrom inside to outside. thickness of the plate is assumed to be unity. The width of the plate is symbolized as t, the length of the plate, L, taken as 4t, and the crack length, a, taken as 0*2t. The plate is considered for plane stress and plane strain cases,respectively. The Young modulus and the Poisson ratio of material properties assumed for analysis are taken as E = 10’ MPa and v = O-3.The width is taken as t = 1m, and the load is taken as P = lo7 N. Five finite-element meshes,termed as (i), (ii), (iii), (iv) and (v), respectively, shown in Figs 6-9, are applied in the analysis. Mesh (i) is a single-edge notched plate and the others are single-edge cracked plates. Most of the elements making up these meshes are standard eightnode isoparametric quadrilateral elements. But meshes (iii), (iv) and (v) all contain some particular elements at the crack tip, where mesh (iii) includes some six-node isoparametric triangular elements, degenerated from eight-node quadrilateral elements, which are non-
-crack tip
crack tip-
~
integration contours element edges
Fig. 7. Finite element mesh (ii): (a) overview, {b) eight integration contours surrounding the crack tip orderedfrom inside to outside. singular; mesh (iv) includes some six-node quarterpoint elements which are well known to model Y-‘/~ singularity; and mesh (v) includes some eight-node collapsed elements which can model r-’ singularity. For each eight-node collapsed element there are three nodes at the crack tip location that are free to displace independently, and the mid-side nodes are left at their usual midpoints. Meshes (iii), (iv) and (v) look the same, but mesh (v) has more nodes, and mesh (iv) has some quarter-point nodes. A detailed comparison of these meshesis listed in Table 1. For each mesh, eight integral paths are employed to evaluate the J integral which surround the crack tip from inside to outside; furthermore, to determine the stress intensity factor lu, according to
KI= m
209
Numerical path integration
crack tip2 Fig. 9.
ccack
tip
t-
integration contours elementedges
8. Finite elems:ntmesh (iii): (a) overview, (b) eight integration contours surrounding the crack tip ordered from inside to outside.
Fig.
where El =
E E/( 1 - v*)
for plane stress for plane strain
For all five mesh[es,automatic mode is adopted, since it is much easier to handle than expert mode. In fact, for such dense meshes,preparing the data for expert mode should be a rather dull task. The results of J and KI from these five meshes are listed in Tables 2 and 3 for plane stress and plane strain cases, respectively, and compared with the analytic solution of KI where J,, Js, Js and 5, are the J integral on the contours 1, 3, 5 and 7 surrounding the crack tip ordered from inside to outside respectively, Jlim is the limit value of J to the outside one, Ki = dm, and the error item is the deviation of K, from the analytic solutioni’ Kpnalytlc= 1.087 x lo7 N/m312. According to this comparison, it is clear that the algorithm given in the present
Near crack tip details of finite elementmesh(iv).
paper is effective and the scaling process does not cause serious error since it is found that each mesh can give a fairly good solution. It is noted that the results of KI from plane strain cases are all larger than those from plane stress cases by about 0.02 x 10’ N/m3’*. The reason may be associated with discretization error. The result from mesh (i) deviates a little more from the analytical one than those from others, which may be easily understood since there exists a little geometrical difference between the notched plate and the cracked one. Meshes (iv) with six-node quarter-point elements and mesh (v) with eight-node collapsed elements are of perfect accuracy for the plane stress case. It is well known that six-node quarter-pointer elements are important for the elastic fracture problem. Theoretically, a crack-tip field obtained from eight-node collapsed elements is appropriate only for a perfectly plastic body whose crack tip is of r-l singularity, but in the practical computation with the finite-element method, the element arrangement with this type of singular element may obtain satisfactory crack-tip fields for strain-hardening plastic analysis.” According to the comparison shown above, it is clear that eight-node collapsed elements are as applicable as six-node quarterpoint elements for the elastic fracture problem, at least for calculating the stress intensity factor as above, although they just contain r-’ singularity in principle. Thus this type of singular element is of great potential usage for the elasto-plastic problem. 6 CONCLUSION
In this paper, a flexible strategy for numerical path integration is proposed such that the J integral is taken as the example to describe the algorithm that can be extended to the evaluation of other J-like integrals without difficulty. The program PathInt is developed to run in expert mode or automatic mode. The design technique can be applied in other programs since this is a rather good compromise of detailed control and easy use. The numerical example has demonstrated the
210
G. Xue, X. Wang, Q. Du Table 1. The five finite-element meshes Mesh
Geometry
Element type
Element number
Node number
6) (ii) (iii)
notched cracked cracked
123 32 34
418 437 435
(iv)
cracked
34
435
(-4
cracked
standard eight-node standard eight-node standard eight-node non-singular six-node standard eight-node quarter-point six-node standard eight-node collapsed eight-node
134
457
Table 2. Comp~n
of J integral paths with exact solution; redts for plane stress case
Mesh (103&m) 0) (ii) (iii) (iv> 69
1.203 1.089 1,186 1.176 1.211
(103%/m)
(103%/m)
(103%/m)
1.210 1.136 1.164 1.177 1.187
1.211 1.150 1,162 1,175 1.186
1.211 1.160 1,162 1.175 1,185
(lO’$m) 1.211 1.160 1,162 l-175 1.185
KI = qQ&z ( lo7 N/mV-2) 1.100 1.077 1.080 1.084 1.089
Error (“/I 1.20 -0.92 -0.64 -0.28 0.18
Table 3. Comparison of J integral paths with exact solution; results for plane strain case ~~03~,rn~
~~03~,rn~
~103~,rn~
~~03~,rn~
(lO~~,rn~
KI = gc&JF (10’ N/m%‘r)
8;)
1.045 1.157
1.081 1.157
1.091 1.153
1.090 1.153
1.090 1.153
1.126 1.094
Error (%I 0.64 3.59
(iii) (iv) iv)
1.124 1.118 1.153
1.103 1.116 I.127
1.101 1.113 1.124
1.100 I.113 I.123
1.100 I.113 1,123
1.099 1.106 1.111
1.10 l-7.5 2.21
Mesh
effectivenessof the algorithm and has shown that several finite-element meshes with or without special crack-tip elements are applicable for the elastic fracture problem, and particularly, the eight-node collapsed element is of great potential usage since it is appropriate for both elastic and plastic problems.
The work of this paper is sponsored by Natural Science Foundation of China and Opening Foundation of the State Key Laboratory of Mechanical Structural Strength and Vibration of China. REFERENCES 1. Rice, J. R., A path independent integral and the approximate analysis of strain concentrations by notches and cracks. J. appl. Mech., 1968, 35, 379-86. 2. Schapery, R. A., On some path inde~ndent integrals and their use in fracture of nonlinear viscoelastic media. Int. J. Fracture, 1990, 42, 189-207.
3. Carpenter, W. C., Read, D. T. & Dodds, R. H. Jr., Comparison of several path independent integrals including plasticity effects. Znt. J. Fracture, 1986, 31, 303-23. 4. Riedel, H. & Rice, J. R., Tensile cracks in creeping solids. ASTM STP 700, American Society for Testing and Materials, Philadelphia, pp. 112-130, 1980. 5. Saxena, A., Creep crack growth under non-steady-state conditions. ASTM STP 905, American Society for Testing and Materials, Philadelphia, pp. 185-20 1, 1986. 6. Landes, J. D. & Begley, J. A., A fracture mechanics approach to creep crack growth. ASTM STP 595, American Society for Testing and Materials, Phjladelphia, pp. 128-148, 1976. 7. Owen, D. R. J. & Fawkes, A. J., Engineering Fracture Mechanics: Numerical Methods and Applications. Pineridge Press, Swansea, UK, 1983. 8. Bakker, A., An analysis of the numerical path dependence of the J-integral. Znt. J. Press. Vess. Piping, 1983, 14, 15379. 9. Barlow, J., Optimal stress locations in finite element models. Znt. J. Numer. Meth. Engng., 1976, 10, 243-51. 10. Ewalds, H. L. & Wanhill, R. J. H., Fracture Mechanics. Edward Arnold, London, Delftse Uitgevers Maatschappij, Delft, 1984. 11. Shih, C. F. & German, M. D., Requirements for a oneparameter characterization of crack tip fields by the HRR singularity. Znt. J. Fracture, 1981, 17, 27-43.
ELSEVIER
Advances in Engineering Software 28 (1997) 203-210 0 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 0965-9978/97/$17.00
SO965-9978(96)00047-6
Technical Note A flexible strategy for numerical path integration Gang Xue, Xucheng Wang & Qinghua Du Department
of Engineering Mechanics, Tsinghua University, Beijing 100084, P.R. China
(IReceived 11 April 1995; received in revised form 24 July 1996; accepted 10 October 1996) In this paper, a flexible strategy for numerical path integration is proposed using the J integral as an example to describe the algorithm that can be extended to the evaluation of other J-like integrals without difficulty. A program called PathInt is developed to run in expert mode or automatic mode. When in expert mode, the program may be controlled in a more detailed way to meet the needsof proficient researchers. When in automatic mode, the program can determine several parameters automatically according to the default settings. A numerical example of a typical elastic fracture problem demonstrates the effectiveness of the algorithm. The applicability of several finite element meshes with or without special crack-tip elements is also investigated. It is found that eight-node collapsed elements are of great potential usefulness since the results from the meshesemploying them at the crack tip are satisfactory in both elastic and plastic problems. 0 1997 Elsevier Science Limited. All rights reserved. Key
words:
numerical path integration, J-like integral, expert mode, automatic
1 INTRODUCTION
when it is applied to the practical problem. The limitation is mainly on the nodal numbering thus also
Since Rice’s’ introduction of the path-independent J integral, many other similar path integrals have been proposed as fracture characteristic parameters, such as Jv integral2 for viscoelastic media, i integral3 for incremental plasticity media, C(tj4 and Ct5 integrals for transient creep media both of which can degenerate into the C* integral6 for steady-state creep media. In most practical fracture problems it is too complex to find the analytical solutions, so numerical techniques, particularly focused on the finite-element method (FEM), play an important role in the solution procedure. Nowadays, there are many general purpose programs for finite-element analysis from which only the solutions of displacement and stress fields can be obtained in most cases. Thus, a posit-processor is needed to evaluate the path integrals deiined above. Concerning numerical methods for fracture mechanics, Owen & Fawkes7 have developed a FORTRAN subroutine named as JevaM to evaluate J integral, but there are some severe restrictions
on the local coordinate system of the element. Nowadays there are so many structured or unstructured, semiautomatic or automatic and interactive mesh generation programs by which it is inconvenient for the user to know exactly the local numbering sequence of the nodes in each element. In this paper, a more flexible strategy is proposed for numerical path integration on the basis of JevaI6. More
general considerations are discussedfor path integration in an individual element. According to our new program named PathInt,
developed with extended FORTRAN
77 language, only the numbering of the elements along the integral path is needed as input data while the program is running in automatic mode for the numerical evaluation of J-like path integrals without thinking about the axis directions of the local coordinate system. In addition, the program has the capability for proficient researchers to control the integration path in detail while the program is running in expert mode. To describe the algorithm, we use, as an example, the 203
G. Xue, X. Wang, Q. Du
(4
one of the local coordinate directions of the element. With respect to this, Bakker* gives the following two important reasons:
rl 5
2
1
4 6
y 8
0
4 ~ 3
t 4
7 I
4
Gaussian
0
Element
Quadrstare
Points
Nodes
04
-5
+
Gaussian
0
Element
Quadrature
Points
Nodes
Fig. 1. (a) Non-corner element on the integration contour.
(b) Corner elementon the integration contour. J integral on contour l? defined as’
where U is the strain energy density, ti is the traction vector and ui is the displacement vector. The numerical example of a typical elastic fracture problem, the single-notched plate under a pin load, is employed to demonstrate the effectiveness and ease of use of the algorithm. The applicability of several finiteelement meshes with or without special crack-tip elements is also investigated by the example; it is found that eight-node collapsed elements are of great potential usefulness since the results from the meshes applying them at the crack tip are satisfactory in both elastic and plastic problems.
(1) A finite-element analysis program normally only produces the stress at the integration points that are considered optimal stress points.’ The path chosen above can avoid the need for inte~olation procedures which usually induce some errors. (2) In elastic-plastic or other nonlinear analyses, the integration points are the only locations where material behavior coincides with the given constitutive law exactly. Below, we discuss in detail the general considerations for an element on the integration path. These considerations relate to the shape, the location and the direction of the path, and are general but not too complicated to handle. The elements intersecting the integration contour can be divided into two categories as non-corner elements and corner elements. In the non-corner element shown in Fig. l(a), the integration path is along one of the local coordinate axes without turning. In the corner element shown in Fig. l(b), the path changes its direction, namely there exists at least one corner in it. There are two sets of non-corner elements in which the integration path is along the n direction with < constant or along the 5 direction with v constant, respectively. In fact, there may exist several corners in a corner element especially for high Gaussian quadrature order. For simplicity, it is assumed that the path is chosen so that at most only one corner exists in the element. Thus there are only two types of corner element in which the integration path is changed from the q direction to the t direction or vice versa. So altogether there are four kinds of element intersecting the integration contour. In the program, an integer variable ikind can take different values according to the respective situations. Thus ’ 1 non-corner element in which path is along the n direction 2 non-corner element in which path is along the 6 direction ikind = 4
2 GENERAL CONSIDERATION FOR PATH INTEGRATION IN AN ELEMENT In computational fracture problems, quadrilateral elements are widely used. Since triangular elements can be considered as degenerated quadrilateral elements, we only deal with quadrilateral elements in this paper. For a path-independ~t integral such as the J integral, the integration contour would be chosen freely in principle. But for the finite-element method, the path is restricted just through element integration points (assuming with NG x NG Gaussian quadrature order), and parallel with
\
3 corner element in which path changes from the Q to the < direction 4 corner element in which path changes from the 5 to the q direction
Furthermore, the location and the direction of the path in an element must be thought about. For each pattern of non-corner element, say, for the case of ikind = 1, there are NG types of integration path with t held are the index and constant at
Numerical path integration
205
and ipeta are the indices of the Gaussian quadrature points so that < held constant as <$xsi on the half-path along the 77direction, and q held constant as qipetaon the other half-path along the E direction, satisfying that 1 < ipxsi, ipeta 5 NG. It is similar to this when ikind = 4. Namely, ipxsi
-5
iab = ideta x 10 + idxsi
and icd = ipeta x 10 + ipxsi @
Gaussian
@
Element
Quadrature
Points
Nodes
Fig. 2. Gaussian quadrature points in an eight-mode quad-
rilateral elementwhen NC = 2. may be 1,2 . . IVG, satisfying c1 < & < . . . < &, namely, the points are ordered according to their
and icd = ipxsi x lOi+ @eta
where 1 along the negative direction of 77-axis (with E held constant)
idxsi =
I ideta =
2 along the positive direction of q-axis (with 5 held constant) 1 along the negative direction of E-axis (with 7 held constant) 2 along the positive direction of {-axis (with 17held constant)
To summarize, two integer variables, ikind and icode, are employed to represent all the information about the path in an individual element. The representation is easy to handle and capable of describing a rather complex situation. For example, for the integration path shown in Fig. l(a) it is that ikind = 1 and icode = 2, and for the case shown in Fig. l(b) it is that ikind = 3 and icode = 2122. The path integration is assembled from all the contributions of each individual element on the contour. The integration in the non-corner element is simple. For example, in the case of integrating along the 77 direction,7 the contribution from an individual element is Jqelem =
where U, and UP are the elastic and plastic work, IAand w are displacements, (T,, aY, and 7XVare stresses.The expression of Jpelemthat corresponds to the case along the < direction is similar. It is a little complicated for the corner element. As shown above, the path in the corner element may be divided into two half-paths. One is along the Qdirection and the other along the < direction. It can be treated in an approximate way.7 As shown in Fig. 3, it is scaled so that J elem &dem+~Jy-m L t
where L is the total range of the integration along the local axis, L = 2 if the local coordinate arranges from - 1 to + 1, and Lf or L, is the actual range corresponding to the direction along the c or Qaxis such that LE and L, are obtained from
and L, = Irlr - ‘%ideI
in which (Et,Q) is the local coordinate vector of the turning point, and [sideor vsidecorresponds to the local coordinate which is held constant on the element side on which the path gets into or gets out of the element. The
206
G. he,
X. Wang, Q. Du
ilocal =2
ie&=+1
I
l l
Gaussian
Quadrature
Element
Nodes
Point
i&l =2 iedge=-l
Fig. 3. The scalingprocessin a corner element.
[sideand qside,- 1 or 1, depend on the pattern of the corner element and the directions of the two halfpaths. In the case shown in Fig. 3,
VdLleS Of
‘%ide
=
-1.
Obviously, some numerical inaccuracy is introduced by this scaling process, but since the contribution of such corner elements to the complete J integral will be small, the error involved should be acceptable, which can be verified by the numerical example shown later. From the above, most integration path situations in an element can be defined by two integer variables. This can meet the demands of the scientist who intends to control the program in a wide range. When the program is running under this kind of thorough control by the user, it is said to be in expert mode. But not every researcher needs to control the path integration completely in this way. In this case, it is convenient to set some default values such that only the most necessary data need to be prepared by the user when the program is running in automatic mode.
3 AUTOMATIC
MODE
In the expert mode described above, the user should know exactly the local information of the element and have complete control to define the path in each individual element intersecting the integration contour in sequence. Hand calculation may be inconvenient, sometimes even nearly impossible, for most practical problems whose finite-element meshesare composed of a large number of nodes and elements, particularly for those whose meshes are generated automatically by some mesh generator. Next, the algorithm determines ikind and icode for each element automatically based on the sequenceof elements on the path when the program is running in automatic mode. Then the user only need to prepare a few necessarydata. A subroutine named GetCode is developed to determine the path-control parameters for each element. It is assumed that the numbering of the previous, current and next elements on the contour are expressed
Fig. 4. The values of ilocal and iedge on each side of the
element. by three integer variables ieprev, iecur, ienext, respectively. By comparing the nodal numberings of the two neighboring elements ieprev and iecur, it can be found which side of them is common. This common side in element iecur may be expressedby two integer variables ilocal and iedge whose values are defined in Fig. 4. Similarly, the common side of elements iecur and ienext can be expressed by another two integer variables jlocal and jedge. If element iecur is the first or last element on the contour, the common side is on the boundary and the values of these integer variables can also be obtained, which will not influence the algorithm afterwards. If ilocal is equal to jZoca1, it is a non-corner element, otherwise a corner element. According to the values of the four variables defined above, the pattern and the direction of the integration path in the element can be determined, but the location, namely the value(s) of ipxsi or/and ipeta are still unknown. Since the consecutive intersections of the contour crossing with the elements are continuous, the starting point of the path in element iecur, at which point the contour intersects the common side of elements iecur and ieprev, is the sameas the endpoint of the path in element ieprev. In the noncorner element, the endpoint can be obtained easily if the starting point of this element is known. In the corner element, with the assumption that the path gets into and out of the element in a symmetric manner, namely LE = L, just as the situation shown in Fig. 3, the endpoint can be obtained too if the starting point in this element is known. Thus the location of the path in each element can be determined once a starting point of the contour is defined. So from a starting point and the numbering sequence of the element series on the contour, the ikind and icode of each element can be obtained. In the program, it is important to treat some particular elements carefully, such as degenerate and collapsed elements. 4 STRUCTURE
OF THE PROGRAM
Since the efficiency of the path integration program is
207
Numerical path integration
Lr’ main
dimea
calculatesvariables
associated
with dynamic dimensioning.
input reads input data that includes definition and definition of integration paths.
of mesh, results of finite element analysis,
SetCaussConst
t
‘E
evaluates
ikind and
icode for an element
BoundEdge -
judges whether the element is on structure sides on boundary if there are.
boundary,
and determines
-
judges whether two elements are adjacent, and determines common side if there is.
the
AdEdge which is the
GetGpInt gets the index of the Gaussian point on the path
GetIntersectCoord -
gets the coordinates of the intersection neighboring elements.
Outputs the coordinates
on the common
side of two
of the Gaussian points on the path.
Fig. 5. The program structure.
not so crucial as that for finite-element analysis, the program PathInt j.s developed inclined to its understandability and expandability. The program is designed modularly, consisting of about 20 modules. The structure of the program is shown schematically in Fig. 5. The two integer variables ikind and icode are the most important co:ntrol data of the program on which the function of the major subroutines are concentrated. According to whether they are input by the user or not, the program, is classified to be running in expert mode or automatic mode. The automatic mode is recommended because of its easy use. For expert mode, the subroutine of CheckCode is necessary to verify the legality of the two control data prepared by the user. In any mode, the subroutine WrPoint is called to output the coordinates of the Gaussian quadrature points on the contour so as to display the contour graphically further.
The program is written in extended FORTRAN 77 language of about 2000 source lines of which about 3040% are comments. The following features are extended based on standard FORTRAN 77 to ensure correctness, understandability and expandability: (1) lowercase letters for keywords and variable names; (2) names longer than six characters; (3) use of ‘!’ to delimit comments; (4) the constructs of ‘do. end do’ and ‘do while. . . end do’. Fortunately, many modern FORTRAN compilers will accept all of these deviations from the standard. 5 NUMERICAL EXAMPLE To demonstrate the algorithm developed in this paper, a single-edge cracked plate is considered under tension of pin loading P as a linear elastic fracture problem. The
208
G. he, X. Wang, Q. Du
(4
(4
t
-crack
crack tip 2
---
tip
integration contours element edges
Fig. 6. Finite element mesh (if: (a) overview, (b) eight integration contours surrounding the notch tip orderedfrom inside to outside. thickness of the plate is assumed to be unity. The width of the plate is symbolized as t, the length of the plate, L, taken as 4t, and the crack length, a, taken as 0*2t. The plate is considered for plane stress and plane strain cases,respectively. The Young modulus and the Poisson ratio of material properties assumed for analysis are taken as E = 10’ MPa and v = O-3.The width is taken as t = 1m, and the load is taken as P = lo7 N. Five finite-element meshes,termed as (i), (ii), (iii), (iv) and (v), respectively, shown in Figs 6-9, are applied in the analysis. Mesh (i) is a single-edge notched plate and the others are single-edge cracked plates. Most of the elements making up these meshes are standard eightnode isoparametric quadrilateral elements. But meshes (iii), (iv) and (v) all contain some particular elements at the crack tip, where mesh (iii) includes some six-node isoparametric triangular elements, degenerated from eight-node quadrilateral elements, which are non-
-crack tip
crack tip-
~
integration contours element edges
Fig. 7. Finite element mesh (ii): (a) overview, {b) eight integration contours surrounding the crack tip orderedfrom inside to outside. singular; mesh (iv) includes some six-node quarterpoint elements which are well known to model Y-‘/~ singularity; and mesh (v) includes some eight-node collapsed elements which can model r-’ singularity. For each eight-node collapsed element there are three nodes at the crack tip location that are free to displace independently, and the mid-side nodes are left at their usual midpoints. Meshes (iii), (iv) and (v) look the same, but mesh (v) has more nodes, and mesh (iv) has some quarter-point nodes. A detailed comparison of these meshesis listed in Table 1. For each mesh, eight integral paths are employed to evaluate the J integral which surround the crack tip from inside to outside; furthermore, to determine the stress intensity factor lu, according to
KI= m
209
Numerical path integration
crack tip2 Fig. 9.
ccack
tip
t-
integration contours elementedges
8. Finite elems:ntmesh (iii): (a) overview, (b) eight integration contours surrounding the crack tip ordered from inside to outside.
Fig.
where El =
E E/( 1 - v*)
for plane stress for plane strain
For all five mesh[es,automatic mode is adopted, since it is much easier to handle than expert mode. In fact, for such dense meshes,preparing the data for expert mode should be a rather dull task. The results of J and KI from these five meshes are listed in Tables 2 and 3 for plane stress and plane strain cases, respectively, and compared with the analytic solution of KI where J,, Js, Js and 5, are the J integral on the contours 1, 3, 5 and 7 surrounding the crack tip ordered from inside to outside respectively, Jlim is the limit value of J to the outside one, Ki = dm, and the error item is the deviation of K, from the analytic solutioni’ Kpnalytlc= 1.087 x lo7 N/m312. According to this comparison, it is clear that the algorithm given in the present
Near crack tip details of finite elementmesh(iv).
paper is effective and the scaling process does not cause serious error since it is found that each mesh can give a fairly good solution. It is noted that the results of KI from plane strain cases are all larger than those from plane stress cases by about 0.02 x 10’ N/m3’*. The reason may be associated with discretization error. The result from mesh (i) deviates a little more from the analytical one than those from others, which may be easily understood since there exists a little geometrical difference between the notched plate and the cracked one. Meshes (iv) with six-node quarter-point elements and mesh (v) with eight-node collapsed elements are of perfect accuracy for the plane stress case. It is well known that six-node quarter-pointer elements are important for the elastic fracture problem. Theoretically, a crack-tip field obtained from eight-node collapsed elements is appropriate only for a perfectly plastic body whose crack tip is of r-l singularity, but in the practical computation with the finite-element method, the element arrangement with this type of singular element may obtain satisfactory crack-tip fields for strain-hardening plastic analysis.” According to the comparison shown above, it is clear that eight-node collapsed elements are as applicable as six-node quarterpoint elements for the elastic fracture problem, at least for calculating the stress intensity factor as above, although they just contain r-’ singularity in principle. Thus this type of singular element is of great potential usage for the elasto-plastic problem. 6 CONCLUSION
In this paper, a flexible strategy for numerical path integration is proposed such that the J integral is taken as the example to describe the algorithm that can be extended to the evaluation of other J-like integrals without difficulty. The program PathInt is developed to run in expert mode or automatic mode. The design technique can be applied in other programs since this is a rather good compromise of detailed control and easy use. The numerical example has demonstrated the
210
G. Xue, X. Wang, Q. Du Table 1. The five finite-element meshes Mesh
Geometry
Element type
Element number
Node number
6) (ii) (iii)
notched cracked cracked
123 32 34
418 437 435
(iv)
cracked
34
435
(-4
cracked
standard eight-node standard eight-node standard eight-node non-singular six-node standard eight-node quarter-point six-node standard eight-node collapsed eight-node
134
457
Table 2. Comp~n
of J integral paths with exact solution; redts for plane stress case
Mesh (103&m) 0) (ii) (iii) (iv> 69
1.203 1.089 1,186 1.176 1.211
(103%/m)
(103%/m)
(103%/m)
1.210 1.136 1.164 1.177 1.187
1.211 1.150 1,162 1,175 1.186
1.211 1.160 1,162 1.175 1,185
(lO’$m) 1.211 1.160 1,162 l-175 1.185
KI = qQ&z ( lo7 N/mV-2) 1.100 1.077 1.080 1.084 1.089
Error (“/I 1.20 -0.92 -0.64 -0.28 0.18
Table 3. Comparison of J integral paths with exact solution; results for plane strain case ~~03~,rn~
~~03~,rn~
~103~,rn~
~~03~,rn~
(lO~~,rn~
KI = gc&JF (10’ N/m%‘r)
8;)
1.045 1.157
1.081 1.157
1.091 1.153
1.090 1.153
1.090 1.153
1.126 1.094
Error (%I 0.64 3.59
(iii) (iv) iv)
1.124 1.118 1.153
1.103 1.116 I.127
1.101 1.113 1.124
1.100 I.113 I.123
1.100 I.113 1,123
1.099 1.106 1.111
1.10 l-7.5 2.21
Mesh
effectivenessof the algorithm and has shown that several finite-element meshes with or without special crack-tip elements are applicable for the elastic fracture problem, and particularly, the eight-node collapsed element is of great potential usage since it is appropriate for both elastic and plastic problems.
The work of this paper is sponsored by Natural Science Foundation of China and Opening Foundation of the State Key Laboratory of Mechanical Structural Strength and Vibration of China. REFERENCES 1. Rice, J. R., A path independent integral and the approximate analysis of strain concentrations by notches and cracks. J. appl. Mech., 1968, 35, 379-86. 2. Schapery, R. A., On some path inde~ndent integrals and their use in fracture of nonlinear viscoelastic media. Int. J. Fracture, 1990, 42, 189-207.
3. Carpenter, W. C., Read, D. T. & Dodds, R. H. Jr., Comparison of several path independent integrals including plasticity effects. Znt. J. Fracture, 1986, 31, 303-23. 4. Riedel, H. & Rice, J. R., Tensile cracks in creeping solids. ASTM STP 700, American Society for Testing and Materials, Philadelphia, pp. 112-130, 1980. 5. Saxena, A., Creep crack growth under non-steady-state conditions. ASTM STP 905, American Society for Testing and Materials, Philadelphia, pp. 185-20 1, 1986. 6. Landes, J. D. & Begley, J. A., A fracture mechanics approach to creep crack growth. ASTM STP 595, American Society for Testing and Materials, Phjladelphia, pp. 128-148, 1976. 7. Owen, D. R. J. & Fawkes, A. J., Engineering Fracture Mechanics: Numerical Methods and Applications. Pineridge Press, Swansea, UK, 1983. 8. Bakker, A., An analysis of the numerical path dependence of the J-integral. Znt. J. Press. Vess. Piping, 1983, 14, 15379. 9. Barlow, J., Optimal stress locations in finite element models. Znt. J. Numer. Meth. Engng., 1976, 10, 243-51. 10. Ewalds, H. L. & Wanhill, R. J. H., Fracture Mechanics. Edward Arnold, London, Delftse Uitgevers Maatschappij, Delft, 1984. 11. Shih, C. F. & German, M. D., Requirements for a oneparameter characterization of crack tip fields by the HRR singularity. Znt. J. Fracture, 1981, 17, 27-43.