Numerical methods for estimating J integral in models with regular rectangular meshes

Cracks and delaminations are the common structural degradation mechanisms studied recently using numerous methods and techniques. Among them, numerical methods based on FEM analyses are in widespread commercial use. The scope of these methods has focused i.e. on energetic approach to linear elastic fracture mechanics (LEFM) theory, encompassing such quantities as the J-integral and the energy release rate G. This approach enables to introduce damage criteria of analyzed structures without dealing with the details of the physical singularities occurring at the crack tip. In this paper, two numerical methods based on LEFM are used to analyze both isotropic and orthotropic specimens and the results are compared with well-known analytical solutions as well as (in some cases) VCCT results. These methods are optimized for industrial use with simple, rectangular meshes. The verification is made based on two dimensional mode partitioning.


Introduction
It has been commonly known that in the scope of linear elastic fracture mechanics the J integral remains a measure of resistance to crack growth equivalent to the energy release rate (ERR) G [1,2]. Moreover, in ideally elastic brittle, as well as ductile-brittle materials with plastic zone at the crack tip, the J integral is highly dependent on the loading mode. Calculating the different mode components of J integral with good precision is therefore vital for correct design against fracture of structures, including laminates. Despite the existence of analytic solutions in terms of both the J integral J [2,3] and the ERR G [4,5], the numerical calculations based on FEM are still considered as the primary efficient method for estimating values of J and G [6].
An important problem in the field of determining resistance to crack growth is estimation of mode mixity for a given body-crack configuration. Complex two-dimensional cases of loading modes are usually described using a relative percentage of the mode II (shearing mode) to the total ERR (G II / G), the so-called mode mix, where: The estimation of mode mix (the so-called decomposition or mode partitioning) has been the subject of numerous scientific works of analytic character [2][3][4][5]. Subsequently, many works of experimental nature have been conducted with the aim of verifying and comparing both methods, as well as discovering their limits of applicability [3,5]. This paper is aimed at presenting and comparing numerical techniques which allow to calculate the J integral and its mode components. The results obtained using these techniques will be compared against those acquired using analytic methods. For detailed description of these methods, the readers are referred to [2,3,4,6].

Numerical mode partitioningthe J integral method
Numerical methods have found application in fracture mechanics analyses since FEM software came into common usage. The J integral method uses the energetic contour path integral described as [1]: The J integral is a path integral over a scalar field along an arbitrary curve in n-dimensional space.
In fracture mechanics problems the integral defined in two and three-dimensional Euclidean space is used. Another primary characteristic of the J integral is that it disappears when a path contour is closed. Consequently, J may be considered a measure of a crack's existence inside the analyzed contour. Additionally, J is an oriented integralthis last property together with the disappearance of J along a closed contour imply that the J integral is path independent.
Under the assumptions of LEFM, equality of the two quantities: J and G is proven. This fact together with the path independence property of the J integral cause J to be useful in analyzing nonhomogeneous structures with interfacial cracks. The energetic approach with the use of the path independent J integral enables to bypass the singularity dominated zone in the vicinity of the crack tip and the issues related to its existence [5,7], particularly in terms of mode partitioning accuracy.
J, similarly to G, may be partitioned into mode I and II (in a two-dimensional loading case), corresponding to opening and shearing fracture. The quantities J I and J II are obtained by partitioning the terms under the integral sign into symmetric and antisymmetric components respectively [8]. The equation (1.1) holds true also in case of the J integral.

Implementation of the J integral calculations
In the following section, two independent numerical techniques of calculating the J integral will be introduced. The details of computational implementation will be given in terms of the complete flowcharts of the methods' algorithms. Although the algorithms may already be used in commercial software environments, they have not yet been published to the best of the author's knowledge. The flowcharts are designed to be used as the basis for development of macros in typical parametric languages, such as ANSYS APDL, however general purpose programming languages may also be used if a FEM solution is known.

Direct integration method
The name direct integration was given to the algorithm in which an integration path is selected as a sorted set of nodes of finite element mesh from a rectangular contour of predefined thickness t (see Fig. 1), where t ≥ 0. In practice, the cases with t equalling zero occur only when a mesh is regular.
Formulas describing the J integral and its components on separate respective segments of Γ can be obtained using the definition (2.1) and the conventions depicted in Fig. 1. The trapezoidal rule is used for calculating Riemann sums of the integrated functions between the nodes (endpoints of respective partitions) of the integration contour (the (*) box, Fig. 3): IOP Publishing IOP Conf. Series: Materials Science and Engineering 175 (2017) 012062 doi:10.1088/1757-899X/175/1/012062 While describing the direct integration algorithm it is worth noticing that, due to specifics of many FEM software packages, it is possible to switch to the nearest node neighbouring the selected node during the iteration process within the inner loop (which iterates over the nodes subscripted with the letter k). This action is described by the condition (**): In general-purpose programming languages, this condition could be realized as a switch between two points (array elements) sorted in ascending order of distance from the starting point P 0 .

EDI integration method
Equivalent Domain Integral (EDI) is the name of a method, which invokes the divergence theorem to convert a line integral of a vector field around a simple closed curve to a domain integral of this field's divergence over the plane region bounded by this curve. By that means the J integral is calculated over a certain compact domain, which consists of finite elements which form a ring around the crack tip. The EDI method also introduces operations and symbols characteristic of a discrete space spread by the finite element mesh [8,9]. Fig. 4 describes the EDI algorithm as implemented in MATLAB environment. The algorithm is executed as a post-processing phase after FEM calculations carried out in ANSYS Mechanical.
In Fig. 2 the use of a parent element coordinate system Oξη is depicted. This system is defined in general as curvilinear, with its origin at an element's geometric centre. The coordinates (ξ, η) are selected (box (*), Fig. 4) so that the corner nodes are positioned in the following way (see Fig. 2): The coordinates of the even-numbered nodes are arithmetic means of the corresponding coordinates of two neighbouring nodes. The local system utilizes Gaussian points, which are the roots of a Gauss-Legendre polynomial, which in turn approximates the integrated function. The Gaussian points are numbered using Roman numerals according to the convention in Fig. 2.
The sum of the integrated function's values at the Gaussian points forms a Gaussian quadrature of the J integral's increment in an element [8,9]. Using the Einstein notation this can be written as:   The S i values are set so that in the inner nodes of the contour ring (nodes numbered 3, 4, 5 according to Fig. 2) the value of S is 1, while in the outer nodes (numbered 1, 7, 8) S is 0. In nodes 2 and 6 the S function takes arbitrary values from the interval [0, 1]. These values are selected in the inner loop of the algorithm in Fig. 4 so that the condition (1.1) expressed in terms of J integral and its components is met with predefined accuracy.
For the purposes of the EDI algorithm and its application, a type of mesh with two geometrical parameters has been developed. Thanks to its simplicity, a derivative of any quantity in the direction x at every Gaussian point may be calculated using the central finite difference formula. Therefore, the values of derivatives are approximated with good accuracy and, at the same time, the distance between the contour ring and the crack tip, as well as the number of elements forming the ring, may be controlled. Obtaining the values of various quantities at Gaussian points may be achieved in most FEM environments by turning off the extrapolation of the Gaussian point results to the nodes.

Verification and application of the J integral analyses
Verification of results of numerical analyses was conducted basing on loading tests (specimens) used for evaluation of the critical ERRs in various modes of loading. All the specimens have geometries of slender beams and remain under plane stress condition (width equalling 25.4 mm). The specimens' models were built with the use of ANSYS Mechanical software. Material properties of the models correspond to isotropic (steel) and unidirectional orthotropic (graphite/epoxy) material as described in Table 1. Values of loads have been selected so that every test case remains in elastic strain regime.   .933 J/m 2 (EDI) respectively with the element's size l e ranging between 0.08 mm and 0.24 mm for isotropic, as well as 0.08 mm and 0.32 mm for orthotropic case. A similar test has been carried out on the End Notched Flexure (ENF) specimen with a non-linear simulation of contact between elements on the opposite cracked legs. ENF solutions also appeared to be of little dependence on l e .
The results of the two tests displayed pure mode I (DCB) and pure mode II (ENF) loading as expected, with roughly 1% and 99% of mode mix respectively. The result J II /J = 94% for the ENF specimen calculated using the direct method can be seen as underestimated. The Single Leg Bending (SLB) test was performed to verify the mode prediction ability of the tested methods in complex loading configuration. A mode mix of approx. 40% was calculated using analytic methods. Fig. 7 illustrates changes of mode mix with increasing distance between the contour path and the crack tip (increment of the parameter k) for constant l e . For both materials examined, the J II /J(k) curves appear to converge except for small oscillations. An exception is the EDI method, which fails when 4-node elements are used for meshing. A big improvement comes after using elements with mid-side nodes. The EDI solution for orthotropic specimen converges on the value of 41%, which lies between the two analytically predicted mode mixes. The solution obtained using direct integration is lower, amounting to approximately 38%. It may be concluded that both numerical methods estimate mode mixes accurately. Also, both numerical methods estimate the total value of J to be 7.5% ÷ 12% higher than analytic solutions. Similar results obtained using the well-known VCCT method suggest a possible influence of FEM modelling on the analytic-numerical discrepancies.

Concluding remarks
In the present paper, the comparison between different methods which allow for the calculation of energy release rate (ERR) and its mode partitioning was presented. The methods used were of both analytic, semi-analytic and numerical origin. In the first part of the paper, the description of the two numerical techniques was provided, which use the J integral approach to calculate the ERR. The matter of numerical implementation in a typical FEM software of these techniques was also presented in detail. In the second part of the paper, the results obtained using the two numerical techniques, known as the direct integration method and the Equivalent Domain Integral method, were compared using three typical types of specimens: the DCB, ENF and SLB. The analytic, semi-analytic and (in some cases) VCCT solutions were considered as comparison criteria. The scope of the analysis (the loads applied, the maximum deflection of the specimens) never exceeded the linear elastic fracture mechanics (LEFM) criteria. The analysis allow for positive validation of the two numerical techniques compared. The direct integration method seems to diminish the J II /J values in the whole scope of integration contours considered. The direct integration results, however, are of less deviation and converge faster than EDI counterparts. The matter of optimizing the integration algorithms for future use in simple rectangular meshes remains forthcoming.