XFEM for Fracture Analysis of Centrally Cracked Laminated Plates Subjected to Biaxial Loads

P. Shailesh ^a*, A. Lal ^b

^a Department of Mechanical Engineering, SRE’s Sanjivani College of Engineering, Affiliated to Savitribai Phule Pune University, Pune, Kopargaon, 432603, India

^b Department of Mechanical Engineering, S.V. National Institute of Technology, Surat, 395007, India

1.     Introduction

In the recent year’s fiber reinforced composites (FRC) are being increasingly used in the aerospace, space craft, marine, automotive, and bio-medical engineering. This is because FRC has several advantages such as mass production, low fabrication cost, good resistance to impact and high ability to meet multi-functional needs (wear, corrosion, thermal resistance and toughness).

During manufacturing and fabrications, due to improper design of the constituent components, imperfection in production and manufacturing procedures and inaccurate measurement of curing parameters these structures are often subjected to different types of internal and external discontinuities in the form of cracks, inclusions, holes, and voids. Since most of these materials exhibit brittle fracture failure with little or no ductility as offered by metals, the behavior of these structures must be understood, and analysis to predict the fracture failure needs to be performed so that the structural safety and reliability of FRC structures can be increased up to a significant level.

In most of the applications, FRC components are often subjected to complex service loading conditions which may lead to development of cracks in these structures. The complex service loading conditions involved in various sensitive applications of FRC materials are rarely uniaxial, typically involving biaxial or even triaxial states of stress. The presence and growth of cracks in composite structures may lead to the initiation and catastrophic failure of the structures. Hence, accurate predictions of fracture parameters like stress intensity factors (SIFs) and their effect on unstable crack propagation under the action of biaxially acting tensile, shear or even the combination of both these stresses becomes necessary.

Several numerical methods have been suggested in the last decades in order to predict and analyze fracture characteristics in composites. Due to simplicity and stability in handling material heterogeneity, non-linearity and complex boundary condition, the finite element method (FEM) is used extensively as a numerical tool in this area. However, the generation of conforming meshes with the evolving discontinuities is an expensive task using conventional FEM. Hence, after several efforts of many researchers, significant improvements in crack modelling using conventional FEM are realized with the development of a partition of unity-based enrichment method.

In earlier days Melenk and Babuska [1] and Duarte and Oden [2] proposed and combined the partition of unity method with the conventional FEM, which was then implemented in actual sense by Mo¨es et al. [3] for modelling and evolving cracks. Further, this method is known as an extended finite element method (XFEM). In XFEM, for modelling the crack, the generalized Heaviside step functions are utilized to account for crack face discontinuity and the two-dimensional linear elastic asymptotic crack-tip displacement functions are used for the crack tip discontinuity.

Several efforts have been made by many researchers in the early days to increase and check the accuracies of XFEM [4-6].  Belytschko et al. [7] and Sukumar et al. [8] described the state of the art for the problems in material science and reviewed the extended, generalized, moving least squares and meshless methods with an emphasis on their applications to various problems in damage mechanics.

To improve the capabilities of XFEM many researchers used and modified this approach to solve more complicated problems [9-13]. For the implementation of the XFEM within the commercial FE code ABAQUS, Giner et al. [14] presented a procedure based on UEL and explained the procedures that interact with ABAQUS. Chatzi et al. [15] proposed improvements in XFEM–GA algorithm by using a weighted average approach.  Abdelaziz and Hamouine [16] presented an overview and recent progress of the XFEM in the analysis of crack growth modelling. 

For use XFEM to carry out fracture analysis of orthotropic materials, particular treatments and more efficient enrichment functions are required. Hattori et al. [17] proposed the new set of enrichment functions for the crack-tip by means of the Stroh’s formulation for fully anisotropic materials. In this direction, the progressive work has been carried out by Asadpoure et al. [18-19] and Asadpoure and Mohammadi [20], they proposed a small change in crack face and crack tip enrichment functions. Ebrahimi et al. [21] by using meshless based ideas, modified isotropic enrichment functions to enable modelling orthotropic problems. Further to increase the solution accuracy of the XFEM in orthotropic materials, Ghorashi et al. [22] adopted orthotropic enrichment functions along with sub-triangle technique and described a new approach for modelling discrete cracks in two-dimensional orthotropic media by the element free Galerkin method. Bhardwaj et al. [23] by modelling crack face and crack tip enrichment functions, performed XIGA and simulated the through thickness cracked FGM and composite plates using FSDT under different types of loading and boundary conditions.  Chen et al. [24] developed adaptive extended (XIGA) based on LR B-splines and by enhancing signed-distance and orthotropic crack-tip enrichment functions for the simulation of arbitrary holes in orthotropic materials. Further, this approach and enrichment functions are modified by Gu et al. [25-26] and Yu et al. [27] for carrying out simulations of through-cracked Mindlin-Reissner plates and also to perform crack growth in isotropic and orthotropic media. Fang et al. [28] by combining XFEM, developed a new enriched technique based on sign function which can eliminate the blending element problems.

In the last decade, XFEM has become one of the important tools to carry out analysis of structures with discontinuities, various researchers have developed different approaches by utilizing XFEM. Liu et al. [29] developed an accurate extended 3-node triangular plate element in the context of XFEM by integrating the DSG to eliminate shear-locking and presented numerical results of buckling failure analysis of cracked composite functionally graded plates subjected to uniaxial and biaxial compression loads. Wang et al. [30-31] developed a novel local mesh refinement approach and described an adaptive numerical framework for modeling arbitrary inclusions and holes in terms of the XFEM. Yin et al. [32] combined XIGA based on the NURBS and Reissner- Mindlin plate theory, they have evaluated the critical buckling parameters and natural frequencies of defective FGM plates with internal cracks or voids. Ma et al. [33] introduced a novel detection scheme suitable for multiple complex flaw clusters in large structures by adopting XFEM. Yu and Bui [34] described a new adaptive local mesh refinement XFEM approach which combines a posteriori error estimation algorithm, a local non-conformal mesh connection strategy. They have presented some introductory studies of the developed method in two dimensional problems of strong and weak discontinuities. Bui et al. [35] for carrying out the analysis of transient DSIFs of two-dimensional fracture problems of FGMs, introduced and improved the extended moving Kriging based meshfree method by introducing three new different types of correlation function which fully eliminate the user parameter for MK shape functions. Ding et al. [36-37] proposed and enhanced a Matlab object-oriented implementation of the variable-node XFEM, an easy-to-use local mesh refinement scheme, they have applied this scheme for modeling strong and weak discontinuities such as cracks, inclusions, voids and for simulations of multiple crack growth in brittle materials.

For carrying out fatigue crack growth analysis, [38-40] proposed and carried out the simulations of three-dimensional fatigue crack growth of 3-D planar, non-planar and arbitrary shape cracks under mechanical and cyclic thermal loads by using XFEM. By applying XFEM, crack interactions have been studied by Pathak [41] in the heterogeneous FGM under mixed mode mechanical and thermal loading environment and Mishra et al. [42] studied the behavior of piezoelectric components in the presence of multiple cracks under thermo-electro-mechanical loading environment. Patil et al. [43] integrated phase field method with multiscale XFEM and simulated the crack growth in highly heterogeneous materials. For 3-D analysis of heterogeneous unit cell, Bansal et al. [44] proposed a multi-split XFEM approach, to reduce the computational cost, a parallel algorithm is developed for large degrees of freedom systems. Kumar et al. [45] proposed a homogenized multigrid XFEM approach for carrying out the crack growth simulations in ductile materials by modelling the geometric and material nonlinearity. They have investigated the load carrying capacity of the structures with microstructural defects.

During fracture toughness tests significant amounts of crack-tip plasticity is exhibited by most of the structural materials, in such situations only biaxially applied stresses are expected to exert a considerable influence on their fracture characteristics along with various biaxial load factors [46-48]. Therefore, there is considerable incentive to perform fracture analysis in which the biaxial stresses along with various biaxial load factors varied over a considerable range of values.

In this direction, earlier work by Liebowitz et al. [49] explained the importance of the second term in the series representation of an exact analytical solution for the crack tip stresses for an infinite center-cracked specimen subjected to symmetric and uniform biaxial loadings. Oladimeji [50] employed the principle of linear superposition by combining the series type analytical solution around the crack tip for studying cracks emanating from a circular hole in a finite sheet under biaxial loading. To simulate crack initiation and propagation, from an initially inclined crack in isotropic and orthotropic material, Theocaris and Papadopoulos [51] proposed a procedure based on either of the two-term S2-criterion (S2) and the T-criterion (T).  Lim et al. [52] suggested to consider the evaluation of singular expression for the stresses and the corresponding displacements near a crack tip and their effect on the predicted crack growth direction in an anisotropic solid subjected to a biaxial loading. [53-54] developed complex variables formulation for the derivation of the complex variable expressions of the elastic field for the orthotropic materials, and studied the elasto-static fracture response of an orthotropic cracked plate. The influence of load biaxiality on the stress intensity factors, crack growth, and on the local stress components in an orthotropic and isotropic medium have been studied by [55-58]. Meek and Ainsworth [59] performed upper and lower bound limit load analyses and finite element analyses of a center-cracked plate under a wide range of biaxial loadings for a wide range of crack sizes and plate lengths.

It is evident from the literature, the various numerical and experimental results indicate that biaxial loadings have shown considerable influence on fracture resistance, crack driving force and load carrying capacity. Considering the heterogeneous nature of composite materials, strengths obtained using uniaxial testing procedures for these materials often does not give the true insights about multi-axially loaded conditions. This puts significant limitations for wide scale, efficient usage of composites in various sensitive engineering structures.

To the author’s knowledge, in contrast to the isotropic materials, lesser attention has been paid to the study of effects of the multi-axial loads in cracked composite structures. Additionally, the ability to successfully model and simulate the multi-axial behavior of composite laminates largely depends on proper choice of numerical algorithm/approach that eventually ensures wide applicability of the respective fracture failure model in describing the corresponding material behavior under variety of complex loading conditions and for different crack configurations and their propagation.

This paper therefore presents a systematic study of the effects of biaxially applied loadings (tensile, shear and tensile and shear combined) on the MSIFs of an inclined through thickness cracked finite angle ply laminated composite plates along with crack propagation path analysis by implementing the well-established, simple and accurate XFEM algorithm. The MSIFs for different biaxial load factors, crack angles, lamination angles, location of the crack along the X and or Y axis by varying crack length to plate width ratio are evaluated. In addition, the detailed investigation of the crack propagation path analysis by applying global tracking crack growth algorithm is carried out in the framework of XFEM. The effect of different types of biaxial stresses on the MSIFs with each step of the crack growth and its effect on the crack propagation direction are investigated.

An excellent agreement of the present results with the literature will be helpful for actual manufacturing and determining the fracture behavior of these materials.

1. FEM Formulation for Fracture Analysis of Composite Structures

To evaluate the MSIFs, essential XFEM formulation for fracture analysis of composite materials developed by Mohammadi [60, 61], applied and explained by Asadpoure et al. [18], Asadpoure and Mohammadi [20] is used in the present analysis.

In XFEM, to model the crack problems a numerical model is generated by dividing the model into two parts; first generate a mesh for the total domain geometry (neglecting the existence of any discontinuities/cracks) by the classical finite element method. Then in the second part by considering the location of discontinuities/ cracks and by using appropriate functions, enrichment of the selected nodes near to the discontinuities/ cracks is carried out by adding a few degrees of freedom to the classical finite element. Then, by defining a transition zone between enriched and non-enriched domains the compatibility is achieved. Thus, the total geometry is divided into three sub domains: the standard FEM domain, elements cut by the crack face are enriched by the crack face enrichment functions and the elements cut by crack tip are enriched by crack tip enrichment functions.

Assume a cracked orthotropic body with global Cartesian coordinates (X, Y) and the axes of elastic symmetry coincide with the local Cartesian coordinate axes (x₁, y₁), the local polar coordinates (r, θ) are defined on the crack tip subjected to arbitrary forces with general displacement and traction boundary conditions as illustrated in Fig.1.

The equilibrium equations and boundary conditions for a cracked body can be written as

    with the following boundary conditions,

where Γ_t , Γ_u and Γ_c are external traction applied on the body, prescribed displacement and traction free crack boundaries, respectively,

σ is the stress tensor and,are the body force and external traction vectors, respectively.

In the XFEM, the displacement field of any point x located within the cracked domain, following approximation along with the classical finite element method (CFEM) is used        

Therefore, the Eq. (6) will take the form,

where u_j is the nodal degrees of freedom in CFEM, a_k is the added set of degrees of freedom to CFEM model, L_j and L_k are the shape functions associated with nodes j and k; and ϕ(x) is the discontinuous enrichment function defined for the set of nodes those have the crack in their influence/support domain. The enrichment function ϕ(x) can be chosen by applying appropriate analytical solutions according to the type of crack/ discontinuity.

The displacement field can be approximated by rearranging Eq. (7) as below

where three parentheses in Eq. (8) represents linear, discontinuous and tip enrichment parts of an approximation respectively. Nodes that belong to two crack-tips are enriched with the two enrichment functions  respectively, for each crack tip, and the nodes which contain the crack surface within their support domain are enriched with the Heaviside function H (ξ).

In Eq. (8) the cf is the set of nodes that have the crack face in their support domain. t₁ and t₂ are the sets of nodes associated with crack tips 1 and 2 in their influential domain, respectively. u_j are the nodal displacements (standard degrees of freedom),  a_k,  are the vectors of additional degrees of freedom for the nodes located on crack face and the two crack tips respectively. The functions  represent the crack tip enrichment functions for both the crack tips. The necessary crack tip enrichment functions for an orthotropic body used in this study are explained in the succeeding section 2.1.1.  The Heaviside function H (ξ(x)) is defined as a step function or the signed distance function used for modelling strong discontinuity/crack face.

It should be noted that, in this study; the value of +1 if the point is on the positive side of the crack face and –1, otherwise is used.

Fig.1. An arbitrary orthotropic body with crack, having global Cartesian coordinates (X, Y), local polar co-ordinate (r, θ) defined at the crack-tip surrounded by contour ᴦ and its interior area A with arbitrary boundary conditions subjected to traction f ^t

2.1. General Displacement, Stress Field and Enrichment Functions for an Orthotropic Cracked Body 

The general form of Hooke’s law for the stress- strain relationship of an orthotropic body can be defined as

where C_ij (i, j =1, 2, 6) are the relevant compliance coefficients of the orthotropic material along the local Cartesian coordinate axes (x₁, y₁) as shown in Fig. 1. Using the basic theories of elasticity and applying the stress function, the governing fourth-order partial differential characteristic equation of an anisotropic body by applying equilibrium and compatibility conditions can be given as

It can be shown that the roots of Eq. (11) are in the complex form, these roots always occur in conjugate pairs as  and having the general form and are derived as [62]. For numerical purposes, only those roots that have either positive imaginary coefficient or negative one is selected.

The roots μ_k of Eq. (12) are purely imaginary, which can be given as , [63]

For an orthotropic case in the plane problems, the Eq. (11) is reduced to the following simplified equation [63]:

Sih et al. [64] derived the two-dimensional asymptotic displacement and stress fields in the vicinity of the crack tip by means of analytical functions and complex variables.

The stress components for pure mode I are defined as,

The displacement and stress components for pure mode II can be defined in the similar way and can be given as.

and

where K_I and K_II  are stress intensity factors (SIFs) for mode I and mode II, respectively, R_e denotes the real part of the complex function, and g_k,  p_k , q_k  can be defined as

2.1.1. Orthotropic Crack Tip Enrichment Functions

By using the general displacement and stress conditions given in the Eqs. (15-26), the necessary crack tip enrichment functions can be extracted by the evaluation of near tip displacement fields for a general traction free crack within an infinite orthotropic body. Asadpoure et al. [18, 19] and Asadpoure and Mohammadi [20] investigated various types of crack tip enrichment functions for an orthotropic media.

In the most general case, the orthotropic crack tip enrichment functions which can cover the whole range of orthotropic media are presented in [20].

The crack tip enrichment functions in the crack tip local polar coordinate system (r, θ) are defined as

where μ_{x or y} have been defined in Eq. (13).

2.2. Evaluation of Stress Intensity Factors (SIFs)

In the XFEM, the global form of linear equilibrium equations for a discrete system by discretizing Eq. (8) can be written as,

where K is the stiffness matrix, U is the vector of degrees of freedom for nodes (for both CFEM and enriched ones by XFEM) and F is the vector of external forces. The global matrix and vectors are calculated by assembling matrices and vectors of each element.

K and F for each element are defined as

and U is the vector of nodal parameters

In Eq. (36), D is the material stiffness matrix and B is the matrix of shape function derivatives.  The B matrix can be as given 

The parameter D for orthotropic and composite material can be defined as;  

(a) For orthotropic material

(b) For laminated composite material

where NL is the number of layers of the laminate, with

and

where β is the fiber orientation angle of the lamina and 

The standard path independent J-integral method for homogeneous orthotropic materials as explained by Mohammadi [61] has been used in this study to calculate the MSIFs. This method can be represented by adding assumptions and exploiting divergence theorem in the form of domain integral approach, the general form of this method can be explained as

where A is the area surrounding the crack tip and would be the interior region of Γ as shown in Fig. 1, n_j is j^th component of the outward unit normal to A, W is the strain energy density,

Using the equivalent domain integral, Eq. (49) can be transformed to

where, Δ_1j is the Kronecker delta and q is a simple function varying linearly from q=1 at the crack tip to q=0 at the exterior boundary ᴦ as shown in Fig. 1.

The interaction integral is employed to compute mode I and II SIFs. It is based on the superimposition of auxiliary and actual fields. The J integral can be obtained by superposition of the actual and auxiliary states of J integrals as:

where M is the interaction integral and can be given as

For linear elastic conditions

The auxiliary stress and displacement fields, which are defined by asymptotic fields near the crack-tip, can be given as

and

where, µ₁ and µ₂ are the roots obtained from Eq. (13),  g_k, p_k and q_k are defined in Eqs. (27- 29), and superscript aux stands for the auxiliary state. The strain of the auxiliary field can be chosen by either the strain–stress relationship, or selected so as to satisfy the equilibrium as well as traction free condition.

It has been proved that for the two superimposed fields [65]

By setting states as  and leads to a system of linear algebraic equations which are used to calculate the MSIFs.

The MSIFs associated with auxiliary and actual states can be evaluated by calculating M from both the Eqs. (53) and (60) and solving a system of linear algebraic equations as.

3. Crack Growth and Propagation Paths with MSIFs Variation

The accuracy and reliability of the analysis of a structure with evolving cracks are primarily depends upon the accurate determination and the continuity of the crack path. It is worthwhile to mention here that in the case of CFEM this requirement becomes computationally expensive and burdensome, as the mesh is required to be updated at each step of crack increment. XFEM gives an elegant way of modelling discontinuities, where the discontinuities arbitrarily aligned with the mesh could be modelled. This alleviates the need of a conforming mesh and hence no mesh update is required, which in the case of CFEM might result in the loss of accuracy as the data is transferred from one mesh to the other.

It is therefore very much essential to select the crack growth criteria very carefully.  Some of the commonly used crack growth criteria/ algorithms are:

I. Minimum strain energy density criteria proposed by Sih [66]

II. Maximum energy release rate criteria proposed by Nuismer [67]

III. Maximum hoop stress or maximum principal stress criteria proposed by Erdogan and Sih [68]

IV. The global tracking algorithm proposed by Oliver et al. [69].

The first three criteria/ algorithms are widely used in CFEM. In these algorithms, the basic requirement is that the algorithm is needed to be evaluated for each individual crack growth segment, which is an expensive and may result in inaccuracy as discussed earlier.

In contrast to these tracking algorithms, Global tracking crack growth algorithm proposed by Oliver et al. [69- 70] traces all discontinuity paths at once and does not need to be evaluated at each individual crack propagation step. This algorithm shows good results [71- 72] in predicting crack paths and can be easily and elegantly incorporated in the XFEM algorithm. The basic idea is to construct a scalar function ψ whose isolines run perpendicular to the direction of principal stresses in all integration points of the investigated structure. However, this comes at the cost of solving additional global system of equations with one degree of freedom per node. An Iso-line is then defined as

If crack growth is indicated by the crack propagation criterion, the value of ψ is calculated at the crack tip. The crack is then extended following the isoline corresponding to this value of ψ. The crack growth criterion which is used for present analyses is an energy-based crack propagation criterion proposed by Xie and Gerstle [73]. According to this criterion the crack propagation condition can be given as follows,

where A_cr is the area of a new crack segment.

4. Numerical Examples

In this section, several numerical examples are presented by developing a computer code in MATLAB. From the literature available it is being observed that, in the investigation of the effects of biaxially applied stresses, MSIFs are considered among the important parameters for fracture analysis and also for the determination of the path of crack propagation. Therefore, firstly to show the accuracy of MATLAB code developed for implementation of XFEM algorithm to handle various fracture problems, the MSIFs are calculated and compared with the exact solution and those available in the literature, and then in the further study the effect of different types of biaxially applied stresses like tensile, shear and tensile and shear combined on the MSIFs by varying the load and crack parameters along with crack propagation paths is studied.

In all the examples considered in this study, the following average material parameters which are being the functions of independent engineering constants (E_ij, υ_ij, G_ij, i, j = 1, 2) are used for all orthotropic and laminated composite plates as suggested by Mohammadi [60]

where E'  is the efficient Young’s modulus, ν' is the effective Poisson’s ratio, δ⁴ is the stiffness ratio, and κ' is the shear parameter.

The following normalized MSIFs, K_I and K_II are used in the present analysis unless otherwise mentioned, and

where   and   are the first and second mode SIFs, the parameter σ, τ and ϛ are the tensile, shear and combined (tensile + shear) stresses respectively, which are assumed as unity, as the thickness of the plate considered is unity for the present analysis (unless otherwise mentioned) and the parameter a is a semi crack length.

4.1. Convergence and Validation of Present Study

For validating the present computer program developed for XFEM algorithm, the program is firstly applied to standard test problem. In this example, a finite plate of isotropic material having material parameters as E = 200 GPa and υ = 0.3 with an inclined central crack under biaxial stress as explained by Tabarraei and Sukumar [10] is considered. A square plate with height H (=5 mm), width W (=5 mm) and unit thickness with centrally located inclined crack with the crack angle α measured with the Y-axis and the semi crack length a=1mm is considered. The plate is subjected to uniformly applied tensile stress of Kσ_X =1 MPa with K (=0.5) and σ_Y = 2 MPa in the X- and Y-directions respectively is shown in Fig. 2 (a). The factor K is the biaxial load factor, which is defined as the ratio of the load applied parallel to the crack to that applied perpendicular to it, where the crack is aligned with the X or Y coordinate axis. The meshing, modelling of the plate and elements selected for the enrichment of the crack face and crack tips are shown in Fig. 2 (b-c).

For the accuracy point of view, convergence studies of various structured meshes are carried out. It is evident from the Table 1 that the present computations of MSIFs for the plate with an inclined center crack with crack angle (α =40°) as explained in the literature [10] are converged well for the 4 nodes quadrilateral elements with mesh size of 40 x 40 elements. Therefore, for further computations, 40 x 40 number of element mesh is used. The interaction integral with semi crack length a_e is calculated within the domain of size as the values of the SIFs become nearly independent of the domain size [61].

The exact stress intensity factors which are the functions of the crack angle α and as given by Tabarraei and Sukumar [10] are

It should be noted here that the exact values of stress intensity factors mentioned by Tabarraei and Sukumar [10] are for an infinite plate, and as the problem at hand has the plate dimensions quite large as compared to the crack length used for analysis, therefore the numerical solution can be compared with the exact solution and results given in [10].

The MSIFs are determined and compared for different values of crack angle α, are presented in Fig. 3 (a-b). The present XFEM results, the results obtained for Quadrilateral meshes in [10] and the results of exact solution are in excellent agreement.

In this study, the effect of the biaxially applied tensile stress on an orthotropic plate with an inclined central crack is studied by varying the biaxial load factor (K) and crack angle (α). Analysis of variation of MSIFs along with crack propagation path is also carried out.

Table 1. Convergence study of normalized MSIFs of the isotropic plate with an inclined center crack with crack angle (α=40°)

(a)

(b)

(c)

Fig.2. Plate subjected to biaxial tensile stress (a) geometry and loading condition (b) X-FEM meshing of the plate with center inclined crack (c) crack tip and crack face enrichment by X-FEM

                                                    (a)

(b)

Fig.3. Comparison of present MSIFs results

4.2. Effect of Variation of Biaxial Load Factor (K) on an Orthotropic Plate with an Inclined Centre Crack Subjected to Biaxial Tensile Stress and Crack Propagation Path Analysis.

In this example, a square plate with a centrally located inclined crack is considered. The dimensions of the plate and crack are as explained in the previous case. The material properties used are E₁₁ =144.8 GPa, E₂₂=11.7 GPa, G₁₂=9.66 GPa and υ₁₂ =0.21. The plate is subjected to uniformly applied biaxial tensile stresses of Kσ_X and σ_Y in the X- and Y-directions respectively.

The obtained results are presented in Fig.4 (a-b), which shows the variation of MSIFs for various biaxial load factors (K= 0, 1, 2, 4) and crack angle (α =0, 10, 20, 30 to 180°).

Fig.4 (a) shows the variation of normalized mode I SIF K_I for various values of K, for K=1 the values of SIF K_I are almost equal and independent of the crack angle while for K=0, 2 and 4 the values of SIF K_I are symmetric about a vertical line passing through crack angle 90° and shows the sinusoidal variation with maxima at α= 90°. Fig.4 (b) shows the variation of normalized SIF K_II for various values of K, it is seen that for K=1, the value of K_II is independent of the crack angle and is equal to zero. Except for K=1, irrespective of the value of K, SIF K_II is zero for α= 0°, 90° and 180° and shows the sinusoidal variation with maxima/minima occurring at α= 45°/135°.

             (a)

       (b)

Fig.4. Effect of variation of K on MSIFs of an orthotropic square plate with a center inclined crack subjected to biaxial tensile stress

In order to examine further the effects of biaxially applied stress, we have predicted the crack extension/propagation path in this case. From the results obtained, it is observed that the values of K_I are higher when K =4 and K_II shows the sinusoidal variation with maxima/minima occurring at α= 45°/135° therefore, the case of an inclined crack with α (= 45°) and K (=4) is considered for crack propagation.

As explained previously, in the present study the global tracking algorithm is used to find the path of crack propagation by defining an ISO-line, crack increment length (Δ_incr) and number of crack increment step. In the present case Δ_incr = 1mm and the number of steps equal to 5 is used. The crack is propagated through the domain till it cut the whole body into two halves. The results are presented in Fig.5 (a-c).

Fig.5 (a-b) shows the path followed by the upper and lower crack tips. The obtained propagating paths are representing the crack growth pattern for the cracks subjected to tensile stress in plates. Both crack tips follow the similar propagation pattern, after first two steps crack follows a straight line perpendicular to X-direction. This indicates the mode I failure of the plate which can be confirmed from Fig.5 (c) which shows the values of K_Iincrease while the values of K_II decrease linearly in each step after the first two steps. This may be due to the fact that, in that zone the plate is in compression because of biaxial stress in Y-direction and after that, as K is large, the crack tips move almost perpendicular to X-direction which confirms the mode I failure of the plate.

4.3. Effect of Fiber Angle (± β°) and Crack Angle (α) on the MSIFs of Centrally Inclined Crack in an Orthotropic Square Plate Subjected to Tensile, Shear and Combined Stress

In this numerical example, the effects of fiber angle and crack angle on the MSIFs are evaluated by applying the different types of stresses like tensile, shear and combined i.e., tensile and shear. The plate geometry and loading conditions for the plate subjected to biaxial shear stress and combined stresses are shown in Figs. 6 and 7, respectively.

4.3.1. Effect of Fiber Angle

In this case, the dimensions of the plate and material properties used are as explained in the previous example, the crack angle is fixed at α=45° and K =2 is considered.

For orthotropic material under consideration the fiber orientation is defined as the angle of the counterclockwise rotation from the crack axis to the fiber axis as shown in Fig. 8.

(a)

(b)

(c)

Fig.5. Crack propagation path for center inclined crack (α = 45°) in an orthotropic plate with K =4, subjected to biaxial tensile stress (a) upper crack tip (b) lower crack tip, (c) K_I and K_II variation

Fig.6. Plate geometry and loading condition subjected to biaxial shear stress

Fig.7. Plate geometry and loading condition subjected to biaxial combined stress

Several fiber orientations (β = 0°, ±10°, ±20°, ±30°, ±40°, ±50°, ±60°, ±70°, ±80°, 90°) are used for the evaluation of normalized MSIFs by applying uniform biaxial tensile, shear and combined stress. The results obtained for MSIFs are presented in Fig. 9 (a-b).

From Fig.9 (a-b) it is observed that biaxially applied tensile, shear and combined stresses follow the same trend for K_I and K_II. The values of SIFs are decreasing up to β = ±30°, and then increases and are maximum at β = ±40° and again, these values decrease up to β = ±80°, for β = 0° and 90° the values are same. Values of K_II are very small in comparison with K_I which shows that there is no effect on the mode- II SIF which is also observed in uniaxial loading. The MSIFs K_I and K_II   follow almost sinusoidal variation about 40 ° fiber angle for all types of biaxial loads considered in this study.

4.3.2. Effect of Crack Angle

As the results of the previous example reveals that MSIFs shows almost a sinusoidal variation about 40 ° fiber angle. Therefore, in the present case the effect of variation of crack angle α on the MSIFs of an orthotropic plate with fiber angle β = ±40° is studied. The plate is subjected to different types of biaxially applied stresses as explained in the preceding section. The MSIFs are evaluated by keeping biaxial load factor K=2.

Fig.8. Positive rotation of principal material axes from material axes

      (a)

     (b)

Fig. 9. Effect of orthotropic axis ± β ° on the MSIFs of centrally inclined crack (α =45°) in orthotropic plate subjected to biaxial tensile, shear and combined (tensile+ shear) stress with (K=2)

(a)

(b)

Fig. 10.  Effect of crack angle α on the MSIFs of centrally inclined crack in an orthotropic plate subjected to biaxial tensile, shear and combined (tensile+ shear) stress with (K=2)

The material properties, plate and crack dimensions are kept same as in the preceding example. The numerical results obtained for K_I and K_II are presented in Fig.10 (a) and (b), respectively.

From the obtained results it is observed that when orthotropic plate is subjected to biaxial tensile stress, values of K_I start increasing up to the crack angle α =90° and then decreasing up to 180°, the values are same for α =0° and 180°, while the values of K_II show the symmetric sinusoidal variation. The values of K_II vary symmetrically about the mean line with K_II = 0, and are zero for the crack angle α =0°, 90° and 180° with maxima at α =50° and 130°.

In the present case, the crack propagation paths along with MSIFs variation in each crack increment step for the tensile, shear and combined stress are determined; the results are presented in Fig. 11-13 (a-c). For this analysis, we considered the case with crack angle α (= 45°), K (=2), Δ_incr (= 1 mm) and the number of steps equal to 5.

From the Fig. 11 (a-b) it is observed that when tensile stress is applied, both the crack tips propagate in the similar manner, and Fig. 11 (c) show the linear variation of MSIFs with K_I greater than K_II, this confirms the mode I failure of the plate. This may be due to the crack tips start to open almost perpendicular to the direction of application of stress as the value of K is kept equal to two.

From the Fig. 12 (a), it can be observed that when shear stress is applied, the crack tips start sliding and show the mode- II failure. The Fig. 12 (b) shows the variation of MSIFs, it is observed that KII increases linearly while KI decrease after every two steps, which indicate that the crack becomes unstable and propagate in failure mode- II.

Fig. 13 (a) shows the mixed mode crack propagation paths when combined stress I applied. In the first step crack propagates in mode-I and further it shows oscillations which indicate that the crack direction becomes unstable.

This is due to the fact that, after the first step, because of combined stresses the crack tips, slides and opens simultaneously and starts oscillating and may get fracture failure in the mixed mode manner. From the Fig. 13 (b) it is observed that the values of MSIFs obtained are greater than those obtained when tensile and shear stresses are applied. Again, it is observed that the K_I is greater than K_II, which indicates that though the crack tips show mixed mode failure, still K_I shows dominant effects than K_II in the fracture failure of the plate.

(a)

(b)

(c)

Fig.11. Crack propagation path for center inclined crack (α = 45°) in orthotropic plate with K =2subjected to biaxial tensile stress (a) upper crack tip (b) lower crack tip, (c) K_I and K_II variation

4.4. Effect of variation of K on MSIFs

The results of the previous example, show that, when an orthotropic plate is subjected to biaxial stresses the value of K, crack angle along with the fiber angle has considerable effects on the fracture failure of orthotropic materials. Therefore, to assess the effect of K and fiber angle, in this example, MSIFs of symmetric angle ply [0°/β°/ β°/0°] laminated composite plate by keeping β (=15°, 30°, 45º, 60°, 75º, and 90°) and keeping K (= 2, 1, 0.5, and 0.25) are assessed. The geometry of the crack, plate dimensions and material properties are kept same as in the previous example, the crack angle is fixed at α=45°.

(a)

(b)

Fig. 12. Crack propagation path for center inclined crack (α = 45°) in orthotropic plate with K =2subjected to biaxial shear stress (a) crack tip (b) K_I and K_II variation

(a)

(b)

Fig. 13. Crack propagation path for center inclined crack (α = 45°) in orthotropic plate with K =2subjected to biaxial combined stress (a) crack tip, (b) K_I and K_II variation

The MSIFs are evaluated by applying uniform biaxial tensile, shear and combined stress, and results are presented in Table 2.

The obtained results reveal that, for tensile stress the value of MSIF is increased up to a fiber angle β = 45° and then decreasing up to 90°.

The MSIFs obtained for all values of K for [0°/45°/ 45°/0°] laminated composite plate are higher among the lamination schemes considered for tensile stress.

For shear and combined stresses, the value of MSIFs is continuously increasing from a fiber angle β = 15° up to 90° and MSIFs are maximum for [0°/90°/ 90°/0°] laminated composite plate for all the values of K.

From Table 2 it is observed that for K =2, and 0.5, the MSIFs obtained are exactly same for all types of stresses considered in this study. This reveals that for the same value of total applied stress there is no effect of changing the direction of stress from the X-axis to Y-axis or from Y-axis to X-axis on the MSIFs when α=45° of an inclined crack. For K=0.25 the MSIFs are higher than the MSIFs of other values of K, this indicates that for an inclined crack with α=45°, the stresses applied along the Y-axis have a dominant effect in comparison with the stresses applied along X-axis.

From Table 2, again, it is observed that for        K =2, 0.5, and 0.25, the MSIFs obtained for shear stress are higher while for the tensile stress are lower and vice versa for K=1. This is in contradiction with the results obtained in section 4.3.2; this indicates that there is considerable effect of fiber angle, lamina configuration and variation of K on the fracture behavior of the laminated composite plate when biaxial stresses are applied.

In order to examine the effect of lamination scheme and the value of K on crack extension/propagation path, the [0°/45°/ 45°/0°] laminated composite plate with an inclined crack with α (= 45°) and K (=0.25) is considered in the analysis. For this analysis the Δ_incr and the number of steps is kept same as in the previous example.

Table 2. Effect of biaxial load factor K on SIFs of centrally located inclined crack (α=45°) in symmetric angle ply laminated composite square plate subjected to biaxial tensile, shear, combined (Tensile + shear) stress

(a)

                                                          (b)

Fig. 14. Crack propagation path for center inclined crack (α=45°) in [0°/45°/45°/0°] laminated composite plate with K =0.25subjected to biaxial tensile stress (a) crack tip, (b) K_I and K_II variation

From Fig. 14 (a) it is observed that when the tensile stress with K (= 0.25) is applied, the crack tips move almost perpendicular to the Y-axis and follow a straight path, which show mode-I failure along the Y-axis. From Fig.14 (b) it is observed that, initially when the plate is under the influence of biaxial stress with less stress along X-axis, the values of K_II are larger than K_I for the first two steps and crack shows some kinks. After this, in the further steps, the crack tips get stabilized and start moving in a straight path and the values of K_I start increasing and are larger than K_II.

From Fig. 15 (a) it is observed that for shear stress, the crack follows the same mode-II propagation pattern as observed in the previous section 4.3, this can be confirmed from Fig. 15 (b).

It is observed from Fig. 16 (a), when the crack is subjected to combined stress; the crack tips follow a zigzag trajectory, forming successive kinks. These phenomena occur as K_II begins to increase and the ratio of K_Ito K_II begins to decrease, this can be observed from the results of K_I and K_II shown in Fig.16 (b). In such cases after certain steps, this ratio becomes too large and the crack path is susceptible to oscillation and becomes unstable as identified by Belytchko and Flemming [53].

4.5. Effect of Crack Eccentricity in Y-Direction, Crack Length and Biaxial Load Factor on MSIFs of Laminated Composite [0°/45°/45°/0°] Plate

In this example an eccentric crack symmetric about Y-axis and subjected to uniformly applied biaxial tensile, shear and combined stresses as shown in Fig.17 (a-c) is considered. The [0°/45°/45°/0°] laminated composite plate is used to evaluate the effect of crack length and variation of K on the MSIFs. The material properties and plate dimensions used are as in the previous example. The eccentricity of the crack is defined by the ratio e_y/H, where e_y is the distance of the crack along the Y-axis and H is the semi height of the plate as shown in Fig. 17 (a-c).

To verify the effect of variation of crack length, the ratio of crack length to semi plate width (a/W) is used. The MSIFs are evaluated by keeping e_y/H (= 0. 2, 0.4 and 0.8), a/W (=0. 1, 0.2, 0.3, 0.4 and 0.5) and K (= 0.5 and 2). The results are presented in Table 3.

From the results it is observed that for an eccentric crack with eccentricity along Y-axis, for tensile stress the values of K_Iare greater than K_II, for shear stress the values of K_II are greater than K_I, but for combined stress the values of K_II are slightly greater than K_I for K=0. 5 and vice versa for K=2.

(a)

(b)

Fig. 15. Crack propagation path for center inclined crack (α = 45°) in [0°/45°/45°/0°] laminated composite plate with K =0.25subjected to biaxial shear stress (a) crack tip, (b) K_I  and K_II variation

Table 3. Effect of eccentricity e_y/H and a/W ratio on the laminated composite [0°/45°/45°/0°] square plate with eccentric crack subjected to bi-axial tensile, shear, combined (Tensile + shear) stress with different K