A novel numerical approach for fracture analysis in orthotropic media

This paper presents a novel approach for fracture analysis in two-dimensional orthotropic domain. The proposed method is based on consecutive-interpolation procedure (CIP) and enrichment functions. The CIP were recently introduced as an improvement for standard Finite Element method, such that higher-accurate and higher-continuous solution can be obtained without smoothing operation and without increasing the number of degrees of freedom. To avoid re-meshing, the enrichment functions are employed to mathematically describe the jump in displacement fields and the singularity of stress near crack tip. The accuracy of the method for analysis of cracked body made of orthotropic materials is studied. For that purpose, various examples with different geometries and boundary conditions are considered. The results of stress intensity factors, a key quantity in fracture analysis, are validated by comparing with analytical solutions and numerical solutions available in literatures. Index Terms — consecutive-interpolation procedure, crack analysis, enrichment functions, orthotropic materials, stress intensity factor


INTRODUCTION
hanks to its specific high strength and stiffness per unit weight, orthotropic composite materials have been widely used in many modern engineering applications such as automobile industries, shipbuilding and aerospace components.Due to the demand to improve the durability of those structures, studies on fracture behavior of orthotropic media has arisen as an important and indispensible task.Analytical investigation on fracture mechanics of orthotropic materials have been presented for some particular problems with relatively simple geometries and boundary conditions [1,2,3].For more complicated problems, as usually encountered in engineering structures, numerical approach is more suitable.
Currently, the finite element method is the most popular which is widely used in both academic and industrial communities due to its simplicity and efficiency.To avoid the cumbersome task of remeshing in modelling cracks, the extended finite element method (XFEM) was proposed [4].In XFEM, cracks are not directly modelled as geometric discontinuities.Instead, additional terms, namely enriched terms, are introduced into the approximated displacement formulation to mathematically describe the discontinuities.Usually, the enriched functions are proposed based on knowledge of closed form asymptotic fields at crack tip, see [4] for isotropic material and [5] for orthotropic materials.Although the enriched functions proposed by [5] for orthotropic materials have been later employed by many authors [6, 7, 8, 9], they are not sufficient in the special case of isotropic materials.In the last few years, a new type of enriched functions, namely ramp functions, was proposed [10], which is not based on analytical solution.However, the ramp-type function is currently not suitable for orthotropic media as there is no information on material orientation incorporated.
Despite of popularity, XFEM still contains the inherent drawbacks of FEM.For example, the gradient fields calculated by FEM (as well as A novel numerical approach for fracture analysis in orthotropic media Nguyen Ngoc Minh, Nguyen Thanh Nha, Bui Quoc Tinh, Truong Tich Thien T XFEM), e.g.strain and stress, are non-physically discontinuous across nodes.Recently, the consecutive-interpolation procedure (CIP) has been introduced as an improvement for FEM [11,12,13], such that the conventional FEM formulatin is enhanced by averaged nodal gradient.Desirable properties of CIP include the smooth stress fields and the higher accuracy of field variables due to refined interpolation.In this paper, the CIP is combined with the enriched functions to model behavior of twodimensional cracked solids.A slight modification is proposed for the enriched functions originally developed by [5], as an attempt to clear the gap between orthotropic materials and isotropic materials.
The outline of the paper is as follows.A brief on CIP formulation on a particular case of 4-node quadrilateral element is reported in Section 2. Section 3 presents the application of CIP in linear elastic fracture mechanics with the aid of enriched functions.Several numerical examples are investigated in Section 4, in order to demonstrate the accuracy of the proposed approach.Conclusions and remarks are given in Section 5.

BRIEF ON FORMULATION OF THE CONSECUTIVE -INTERPOLATION 4-NODE QUADRILATERAL ELEMENT (CQ4)
Details on the formulation of the CQ4 element was previously described by the authors [12].In this paper, the consecutive-interpolation procedure (CIP) is briefly presented for the sake of completeness.Consider a 2D body in the domain Ω bounded by Г, which is discretized into nonoverlapping sub-domains Ωe, namely finite elements.Any function u(x) defined in Ω can then be approximated using the CIP as where n is the number of nodes I u ˆis the nodal value of function u(x) and RI is the consecutiveinterpolation shape function associated with node I (global index).The vector of shape functions R is determined by Here, Si is the the set of elements interconnected at node i, and e being the area of element e.
It is important to highlight that the auxiliary functions ϕI, ϕIx, ϕIy have to be developed for each type of elements [11,12], which is actually not a trivial task.Fortunately, a general formulation to determine auxiliary functions for a wide range of finite elements from one to three dimensions has been recently proposed by [9], resolving the bottleneck.For the sake of completeness, the general formulation of auxiliary functions is shown in the followings.
Given an element e with ne number of nodes, the auxiliary functions associated with the local i th node where N is the Lagrange shape functions and Σ1 and Σ2 are determined by Replacing x-coordinate by y-coordinate in (8), the function ϕiy is obtained.
Fig. 1 illustrates the application of CIP approach into Q4 element described particularly in an irregular finite element mesh.The point of interest x is located inside a 4-node quadrilateral element, where the four local nodes are subsequently denoted as i, j, k, m.The four sets Si, Sj, Sk, Sm are established by collecting the elements share the node i, j, k, m, respectively.Once the sets Si, Sj, Sk, Sm are determined, the consecutive-interpolation shape functions can be calculated through (2).As shown in Fig. 1, the set of nodes that support a point of interest x used in CIP is in any cases larger than that of the conventional FEM, because it includes not only the nodes of the element containing the point x but also the nodes of the adjacent elements.

Governing equations
The governing equation for static equilibrium in a domain Ω bounded by Г, assuming small displacements and small strains, is given by 0 = b σ , (9) where b is the body force and σ denotes the Cauchy stress tensor.The stress-strain relation is given by Hook's law: ) in which C is the fourth-order tensor of material property and the strain tensor ε is determined by The associated boundary conditions are as follows u u = on Γu: prescribed displacement, (12) t n σ = on Γt: prescribed traction, (13) with Г = Гu + Гt and Гu ∩ Гt = {ø}.A crack existing in Ω is denoted by Гc, which is assumed to be traction free.

Enriched formulations
In order to mathematically describe the discontinuity, enriched approach [4] is usually employed.Recently, [14] introduced the extended consecutive-interpolation 4-node quadrilateral element (XCQ4), which incorporates the enriched formulation into CQ4 element, such that the approximated displacement field in Eq. ( 1) is rewritten by (14) In ( 14), J split is the set of nodes belong to elements completely cut by crack and K tip is the set of nodes belong to the elements containing the crack tips.The employment of enriched terms lead to additional DOFs aj and bk.Function H(x) is the Heaviside step function, describing the jump in displacement fied across the crack, while the four branch functions F α (α=1,...,4) are crack-tip enrichments, capturing the singularities of asymptotic stress fields.For 2D linear isotropic elasticity problems, the four branch functions are given in the local polar coordinate (r, θ) defined at the crack tip by Functions gq(θ) and θq (q =1,2) are defined as where Cij is the components of the tensor of material property as defined in (10).
When the material is isotropic, i.e. s1x = s2x = 0 and s1y = s2y = 1, the branch functions calculated using (16,17,18) degenerate into , which are clearly not sufficient as a set of basis functions.Thus, a modification for functions g1 and g2 in ( 17) is proposed as follows With equation (20), the enriched functions degenerates exactly into (15) in the special case of isotropic material.

Computation of Stress Intensity Factors
Stress Intensity Factors (SIFs) are important parameters reflecting the singular fields near the crack tip in linear elastic fracture mechanics.
Numerically, SIFs can be determined by the following relation: in which KI, KII are the mode I and mode II SIFs, respectively.Subscript (1) denotes the present state of the cracked body, while subscript (2) denotes an auxiliary state, which can be chosen as the asymptotic fields of pure mode I (i.e., KI (2) = 1 and K II (2) = 0) or pure mode II (i.e., K I (2) = 0 and K II (2) = 1), see [4,   (1,2) is a path-independent integral, namely interactive integral calculated as follows [6, 7] Here, Γ is an arbitrary contour surround the crack tip, which encloses no other types of discontinuities, and nj is the j th component of the outward unit vector normal to Γ.

NUMERICAL EXAMPLES
Four isotropic and orthotropic problems are examined in this section to assess the accuracy and performance of the proposed method.All the numerical examples are listed in the following for clarity: Finite rectangular isotropic plate with an edge crack Finite rectangular orthotropic plate with an edge crack Finite rectangular orthotropic plate with a central slanted crack An inclined center crack in an orthotropic disk subjected to point loads The standard extended 4-node quadrilateral element is denoted by XQ4 and XCQ4 denotes the extended consecutive-interpolation 4-node quadrilateral element.Note that the modified enriched functions in (20) are used by default, otherwise, it is stated clearly whether equation ( 20) or (17) are used.

Finite isotropic rectangular plate with an edge crack
An isotropic plate with an edge crack under uniform tensile loading σ0 = 1 is considered in this problem, see Fig. 2. The purpose of this example is to demonstrate that the modified branch functions using (20) performs better than the ones originally proposed by Asadpoure and Mohammadi [5] in case of isotropic material.The plate is determined by L = 2W = 16 and a crack length a. Material parameters are given by: Young's modulus E = 1000 and Poisson's ratio ν = 0.3.This problem is pure mode I, in which the closed form stress intensity factor KI is given by [4] 1 is achieved by XCQ4-(20), demonstrating that XCQ4 element, with the enhanced consecutive-interpolation, outperforms the XQ4 element.Thanks to the consecutive-interpolation procedure, the stress fields evaluated by XCQ4 elements are smooth across element nodes (except for the regions containing crack), which is physically more appropriate than the non-smooth stress provided by XQ4 elements, as depicted in Fig. 3 for the normal stress component σyy.This is the reason that higher accuracy for SIF values, which are based on stress components, is obtained when XCQ4 elements are used.

Finite orthotropic rectangular plate with an edge crack
A rectangular orthotropic plate with an edge crack subjected to distributed load, as shown in Fig. 4, is investigated in this problem.The material is made from graphic-epoxy with the following properties: E1 = 114.8GPa, E2 = 11.7 GPa, G12 = 9.66 GPa, ν12 = 0.21.The crack length is determined by a/W = 0.5.The same mesh of 25 x 49 quadrilateral elements as in Example 4.1 is used for discretization of the problem domain.
Effects of the material orthotropic angle β on mixed-mode stress intensity factors are shown in Fig. 5.The present approach is in good agreement with reference results [5,6,15].KI tends to increase from β = 0 o to β = 45 o and decreases from β = 90 o .
For KII, the peak value is reached at about β = 30 o .In Fig. 5, the SIFs are normalized by The problem domain is discretized by a mesh of 45 x 91 quadrilateral elements (i.e.4232 nodes).Fig. 7 depicts the variation of normalized mode I and mode II SIFs with respect to the inclined angle of crack.As the inclined angle increased from 0 o to 90 o , mode I SIF, I K ~, gradually decreases from 1 to 0, while the mode II SIF, II K ~, increases to the peak value at 45 o and then decreases to 0. Good agreement with results using meshfree method (4560 nodes) reported in [6] is observed.Largest discrepancy in Fig. 7 is recorded between the curves of II K ~at slanted angle α = 45 o .Thus, further comparison is conducted and reported in Table 2, showing the consistency between present approach (XCQ4 -(20)) and literatures.

CONCLUSIONS AND OUTLOOKS
In this paper, the XCQ4 element has been successfully extended for modelling cracks in twodimensional orthotropic problems.The accuracy and performance of the present formulation has been verified through a series of numerical examples.Preliminary results indicate that the present approach is in good agreement with other authors.Furthermore, XCQ4 element is observed to perform better than its XFEM counterpart, the XQ4 element, such that higher accuracy of SIFs is achieved.As SIFs are key quantities to numerically determine the propagating direction during crack advancement, the approach is promising to be extended to problems involving crack growth.
The higher accuracy of XCQ4 over XQ4 is possibly due to the enhanced interpolation by CIP, by which the erroneous non-smooth stress fields in XQ4 can be overcome by XCQ4.It is important to emphasize that no extra degrees of freedom is required for CIP.Although in this work, only the quadrilateral element is investigated, the approach is possible for other types of element.With the aid of general formulation for auxiliary functions, see [13], CIP can be integrated into a wide range of existing finite elements without difficulties The set of crack-tip enriched functions proposed by [5] is shown to be not well-chosen.Thus, a modified version of the enriched functions is presented, which properly degenerates into those proposed by [4] for the special case of isotropic material.The new set of enriched functions outperfoms the set by [5] when material is isotropic.For orthotropic material, the new set of enriched functions is consistent with references available in literatures.
of Lagrange shape functions with respect to x-and y-directions,

Figure 1 .
Figure 1.Schematic sketch of CQ4 element 3 APPLICATION OF CIP IN LINEAR ELASTIC FRACTURE PROBLEMS For 2D linear orthotropic materials, the crack-tip enrichments are introduced as in[5, 11]

Figure 2 .
Figure 2. Example 4.1 Isotropic rectangular plate with an edge crack under uniform tensile loading For numerical calculation, a mesh of 25 x 49 quadrilateral elements (1300 nodes) is used to discretize the problem domain.The values of KI evaluated for different ratios a/W are presented in Table 1, in which a comparison between XQ4-(17), XQ4-(20), XCQ4-(20).Analytical solution is used as reference to assess the accuracy of three approaches.Results indicate that the evaluation of KI by XQ4-(20) is closer than XQ4-(17), evidently showing the appropriateness of the modified enriched function in (20).The highest accuracy inTable 1 is achieved by XCQ4-(20), demonstrating

Figure 4 .Figure 5 .
Figure 4. Example 4.2: Orthotropic rectangular plate with an edge crack under uniform tensile loading

Figure 7 .
Figure 7. Example 4.3: Normalized mode I and mode II SIFs computed according to the crack inclined angle

TABLE 1 .
EXAMPLE 4.1: VALUES OF KI CALCULATED WITH

TABLE 2 .
EXAMPLE 4.3: COMPARISON OF NORMALIZED MODE I AND MODE II SIFS AT INCLINED CRACK ANGLE 45 O