Studies of fiber-matrix debonding

Abstract

In this paper, the debonding of a single fiber-matrix system of carbon fiber reinforced composite (CFRP) AS4/Epson 828 material is studied using Cohesive Zone Model (CZM). The effect of parameters namely, maximum tangential contact stress, tangential slip distance and artificial damping coefficient on the debonding length at the interface of the fiber-matrix is analyzed. Contact elements used in the CZM are coupled based on a bilinear stress-strain curve. Load is applied on the matrix, tangential to the interface. Hence, debonding is observed primarily in Mode II. Wide range of values are considered to study the inter-dependency of the parameters and its effect on debonding length. Out of the three parameters mentioned, artificial damping coefficient and tangential slip distance significantly affect debonding length. A thorough investigation is recommended for case wise interface debonding analysis, to estimate the optimal parametric values while using CZM.

Keywords

single fibre / cohesive zone model / interface debonding / carbon fiber reinforced composite (CFRP)

Introduction

Carbon fibers are extensively used in advanced applications like aerospace [1−4], formula one [5−8] oil and gas industry [9], to name a few. With the advancement in processing composite fibers the cost of manufacturing is reduced, which leads to diverse applications of composites [10−13]. Because of its low weight, high strength, electro-magnetic resistance, thermal resistance and high wear and tear, carbon fiber reinforced composite (CFRP) is desirable in the majority of advanced applications [1,14].

CFRP composites fail mainly in three ways; 1) fiber damage, 2) matrix damage and 3) fiber-matrix interface debonding [4]. In the present work, focus is mainly on failure by interface debonding. In fiber reinforced composite materials, the interface is defined as the entire ‘shared’ surface between the fiber and the matrix [15]. The interface between carbon fibers and matrix plays a critical role in controlling the overall properties of composites, such as mechanical performance, fracture toughness and environment stabilities. Interfacial characteristics determine the way loads can be transferred from the polymer to the fiber [16]. To tailor composites having desirable properties, it is necessary to know mechanisms of the polymer−fiber adhesive contact formation and its behavior under mechanical loading [17].

The interface contains air bubbles or voids, and impurities absorbed during the fabrication of the composite. The cooling process involved during fabrication induces micro-cracks that develop due to the thermal residual stresses. The failure mechanism at the interface, induced through applied loading, primarily determines the overall mechanical response and strength of the composite. Interfacial failure is a multi-stage progression in which each stage is characterized by certain parameters. The process is fundamentally nonlinear due to the effects of matrix plasticity and interface crack initiation and growth, together with the associated frictional stress transfer. These characteristic features are essential in determining the reinforcement capability of the fiber. Thus traditional assumptions of linear-elastic behavior and perfect adhesion between the fiber and matrix are no longer applicable and can lead to erroneous results [18].

Experimental methods of fracture analysis of CFRP are summarized in [19]. It is widely recognized that an interfacial region exists between the fiber and the matrix in fiber reinforced composites. The importance of the interfacial region has been under extensive study in recent years both experimentally and analytically. It is believed that the properties of the interface play a critical role in the performance of the composite materials, especially fiber-matrix adhesion [20]. Since, the fiber is embedded in the matrix, the load is primarily applied on the matrix. One of the major roles of the interface is to transfer the stress from the lower modulus matrix to the reinforcing higher modulus fibers. Hence, interface ultimately determines the maximum load that can be applied to the reinforcing fiber.

Interfacial region in [20] is a bonded layer with thickness 0.16µm. In the present study, model of a single carbon fiber surrounded by a matrix with interface layer reduced to zero thickness is considered. Figure 1 shows a schematic of the fiber-matrix composite material at different scales. The use of single fiber-matrix composite can be viewed as a starting point in understanding carbon composite properties at a micro level. With the use of multiscale methodologies [21−37], such an analysis can be expanded to fit the full scale CRFP structures.

Several models are proposed in Ref. [38] to study the debonding/delamination of a composite fiber. These estimate the delamination growth of composites with an initial crack based on the Linear Elastic Fracture Mechanics (LEFM) approach. In the virtual crack closure technique (VCCT) [39−43], J-integral method [44], virtual crack extension technique [40,45] and the stiffness derivative method [46], the delamination growth is predicted based on the energy release rate [47]. However, some issues are noticed when above techniques are implemented based on the finite element (FE) codes. For example, calculation of fracture parameters like, stress intensity factor/energy release rate requires the topology information of the nodes ahead and behind the crack front. Such computations are easy for a stationary crack, but are very difficult for a progressive crack [48−50]. On the other hand, several techniques exists in the literature to model the fracture at continuum level [51−56].

The arrangement of the article is as follows. Cohesive Zone Model (CZM) is presented in Section 2. Numerical simulations of the single fiber-matrix debonding analyses are presented in Section 3. The effect of maximum tangential contact stress, tangential slip distance and artificial damping coefficient on debonding length is studied in Section 4. Section 5 concludes this study.

Cohesive zone model (CZM)

The origin of the cohesive crack model goes back to Dugdale [57] who introduced the concept that stresses in the material are limited by the yield stress and that a thin plastic zone is generated in front of the notch. Barenblatt [58] introduced cohesive forces on a molecular scale in order to solve the problem of equilibrium in elastic bodies with cracks. Cohesive zone model is based on the principle that the energy required to separate the interface surfaces (critical fracture energy) can be estimated knowing the traction separation laws. The interfacial separation is defined in terms of the contact gap and the tangential slip distance. The interface surfaces are modeled by the interface/contact elements. The computation of contact gap and tangential slip is based on the type of contact element and the location of contact detection point.

CZM relates tractions to displacement jumps at an interface where a crack may occur. Damage initiation is related to the interfacial strength, which is estimated as the maximum traction on the traction−displacement relation. When the area under the traction−displacement curve is equal to the fracture toughness, the traction is reduced to zero and new crack surfaces are formed. The advantages of cohesive zone models are their simplicity and the unification of crack initiation and growth within a single model. Although the cohesive damage models cannot be considered non-local damage models [59], they allow a mesh-independent representation of material softening, provided that the mesh is sufficiently refined [38]. CZM is used for bonded contact with Augmented Lagrangian method.

In the present study, Mode II debonding is studied based on the bilinear behavior. The bilinear cohesive zone model is based on the model proposed by Alfano and Crisfield[60]. Mode II debonding defines a mode of separation of the interface surfaces where tangential slip dominates the separation normal to the interface. The tangential contact stress ({\tau }_{\mathrm{t}}) and tangential slip distance (u_t) are related as

(1){\tau }_{\mathrm{t}}=K_{\mathrm{t}}{u}_{\mathrm{t}}(\mathrm{1}-d_{\mathrm{t}}),

where K_{\mathrm{t}} = tangential contact stiffness and d_t = tangential slip distance at the completion of debonding, which is given by

(2)d_t=\left({\displaystyle \frac{u_t-{\bar{u}}_t}{u_t}}\right)\left({\displaystyle \frac{u_t^c}{u_t^c-{\bar{u}}_t}}\right),

where d_t=\mathrm{0} for {\mathrm{\Delta }}_t\le \mathrm{1} and for {\mathrm{\Delta }}_{\mathrm{t}}>\mathrm{1}. {\bar{u}}_t is the tangential slip distance at the maximum tangential contact stress and u_t^c is the tangential slip distance at the completion of debonding. u_t^c is one of the three parameters studied in relation to its effect on the debonding length. The time step is defined as

(3){\mathrm{\Delta }}_t={\displaystyle \frac{u_t}{{\bar{u}}_t}}.

For the 3-D stress state an “isotropic” behavior is assumed and the debonding parameter is computed using an equivalent tangential slip distance

(4)u_t=\sqrt{u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}}},

where u_{\mathrm{1}} and u_{\mathrm{2}} are the slip distances in the two principal directions in the tangent plane. The components of the tangential contact stress are defined as

(5){\tau }_{\mathrm{1}}=K_t{u}_{\mathrm{1}}(\mathrm{1}-d_t),

and

(6){\tau }_{\mathrm{2}}=K_t{u}_{\mathrm{2}}(\mathrm{1}-d_t).

The tangential critical fracture energy is computed as

(7)G_{\mathrm{ct}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\tau }_{\max }{u}_t^c,

where {\tau }_{\max } is the maximum tangential contact stress before debonding initiates. In the present work, a wide range of values of {\tau }_{\max } have been considered to analyze its influence on interface debonding length of fiber-matrix. {\tau }_{\max } is used as the reference parameter for comparison of inter-dependency of the other two parameters.

Numerical simulations

Consider a single fiber-matrix composite body as shown in Fig. 2(a). Because of the symmetry, only a quarter fiber as represented in Fig. 2(b) is chosen for the simulation. The symmetry plane is selected as the plane with positive y axis and negative x axis. Figure 2(c) shows the transverse section of the quarter body along with the boundary conditions. Faces at one end of fiber and matrix are constrained in all degrees of freedom. A displacement load of 8 µm is prescribed along the positive x axis on the other end of the matrix face.

The aspect ratio of the single fiber-matrix body is 126.5. Total length of the fiber is 506 µm and that of the matrix is 507 µm. The outer matrix diameter of the model is 10.8 times the fiber diameter, large enough to ensure that the displacement and stress fields approach an undisturbed uniaxial stress state in the matrix. Fiber diameter is 8 µm. The material properties for fiber and matrix are given in Table 1.

Linear solid elements are used, both for fiber and matrix during all the numerical simulations. The orientation of the solid elements is parallel to the global Cartesian coordinate system.

Analysis settings

Contact zone between fiber and matrix is bonded contact. In this study, CZM code is applied in the contact zone between fiber and matrix. The initial settings of this code are set at reasonably high input values for the variables, so that no debonding is observed at fiber-matrix interface. The initial set up values of CZM used in present work are listed in Table 2.

Meshing of the model is performed with the mesh element size at 9µm. Total number of elements in the model are 2166 and number of nodes are 11776. Convergence study performed is presented in Table 3. Let the debonding length estimated considering 2166 elements be dbl_ref. The % deviation of debonding length is defined as the “change in ratio of debonding length with respect to dbl_ref and dbl_ref, multiplied by 100.” Convergence error occurred for majority of the u_t^c and μ analyses when more than 2166 elements is used. Furthermore, debonding length deviation is only 2.4% at a high increase in number of elements of 3400. With these two considerations, 2166 elements is used for all analyses.

A displacement load of 8 µm is applied to one face of the matrix. The other end of the model is constrained on both fiber and matrix face. The load is applied parallel to the fiber-matrix interface, such that Mode II debonding at the interface can be observed after the values in the above CZM are gradually reduced as described in detail in the following section.

Results and discussion

The objective of this study is to understand the debonding behavior at fiber-matrix interface with reference to three parameters. Out of the six parameters of the cohesive zone model as listed in Table 2, maximum tangential contact stress, tangential slip distance and artificial damping coefficient are selected in the present study to understand their influence on the debonding length.

A force reaction probe is applied on the loading face of the matrix for correlation of debonding length. Debonding length and maximum tangential contact stress are inversely proportional with each other. Since, the load is in the form of displacement, a force reaction plot gives an estimate on the amount of force required for the same effect of debonding, which is directly proportional to the maximum tangential contact stress.

Effect of the maximum tangential contact stress on debonding length

The initial estimate of the maximum tangential contact stress ({\tau }_{\max }) is 500MPa at which zero debonding length is observed. To understand the effect of {\tau }_{\max } on the debonding length, {\tau }_{\max } is reduced in steps of 25 MPa, the effect of which is observed in Fig. 3. As {\tau }_{\max } decreases, the maximum tangential contact stress required for initiation of debonding also decreases. The maximum tangential contact stress value required for the initiation of debonding at fiber-matrix interface is observed at 150 MPa. As {\tau }_{\max } is further reduced tending to zero, in steps of 25 MPa, a sudden increase in the debonding length is observed. Constant load is maintained throughout the process. As the tangential stress limit of the interface is progressively reduced, the load experienced by the interface progressively increases resulting in higher debonding lengths. The reduction in {\tau }_{\max } indicates a weaker structure resulting in a weaker interface. A sudden increase in the debonding length is observed when {\tau }_{\max } drops below 75 MPa. Therefore, after the initiation of debonding at the interface, the fiber-matrix system remains considerably strong in a certain region. After that it collapses suddenly, leading to the complete failure.

Figure 4 shows the decrease in force reaction as {\tau }_{\max } decreases. The force reaction calculations work from summing the internal forces on the underlying elements under a contact region. Figure 4 is in inverse to Fig. 3 where an increase in debonding length as {\tau }_{\max } decreases is observed. Such a contrary relation concludes that, as debonding length increases for decreasing {\tau }_{\max }, the force required for the respective decreasing {\tau }_{\max } also decreases. This trend can be particularly observed between {\tau }_{\max } values of 50MPa and 25MPa where a sudden increase in debonding length in Fig. 3 corresponds to a decrease in force reaction between the same region in Fig. 4.

Effect of tangential slip distance on debonding length

When applied tangential force is increased between two surfaces pressed together, relative micro-displacement can be seen before macro-slip occurs. The displacement generally causes a nonlinear force−displacement relation and this further leads to a hysteresis loss at the surface contacts of elements under a cyclic tangential force. This micro-displacement has also strong influences on the stiffness of contact interface, the stick-slip behavior of sliding table, fretting and other many engineering problems [61]. In Ref. [62], tangential slip distance is studied at a micro level to understand the initial tangential slip distance required to shift from sticking to sliding condition. It is observed that at some critical tangential force, the displacement reaches its maximum value and then complete sliding starts. This displacement can be interpreted as the characteristic length of transition from sticking to sliding.

u_t^c is the characteristic length of transition from sticking to sliding or tangential slip distance at the completion of debonding. For analyses in section 4.1, u_t^c value is 0.001. However, to study the effect of u_t^c on the debonding length, analyses are performed with values of u_t^c ranging from u_t^c = 0.1 to 0.000005 with constant {\tau }_{\max } = 100 MPa. The results of four analyses are presented in Fig. 5. It is observed that as u_t^c value decreases, the debonding length increases.

Comparison of the change in debonding length with respect to the decreasing maximum tangential contact stress at two particular values of u_t^c (0.001 and 0.0005) are plotted in Fig. 6, which also reflects the interdependency of u_t^c and {\tau }_{\max }. When both values of u_t^c and {\tau }_{\max } are reduced, the resultant debonding length is much higher than that observed in Fig. 3. The higher values of debonding lengths observed in Fig. 6 compared to Fig. 3 is because of the reduction of u_t^c from 0.001 to 0.0005. This is an important relation supporting that a slight decrease in tangential slip distance can cause high increase in debonding length. This is in relation that, decreasing tangential slip distance decreases the area under the bilinear CZM-curve, which is proportional to the energy release rate and to the failure energy. Therefore, less energy is required to enforce failure for one CZM-element.

Effect of artificial damping coefficient (μ) on the debonding length

The reduction of the computational effort by reducing the simulation time is of considerable interest. Reduction of the computational effort is especially desirable because many calculations often are necessary during the analysis of nonlinear systems, due to the fact that a number of values of the main design parameters must be accommodated (e.g., rpm, load conditions, valve characteristics, pipe lengths etc.). The reduction of the simulation time for the determination of the steady-state of a nonlinear system can be achieved by adding a controlled “artificial” damping. Artificial damping can be applied to a nonlinear system with multiple degrees of freedom and nonlinearities [63]. It has been numerically observed that the artificial damping, even if slight, makes the coarse grid problem much easier to solve, without deteriorating the overall convergence rate. For most problems, 2-6 times speed up have been observed [64].

This study presents artificial damping not in terms of its effect on convergence rate, rather its effect on the debonding length at fiber-matrix interface. As will be seen, the results are particularly important for further study of interface using CZM model. For the analyses in Section 4.1, the artificial damping coefficient value is 0.001. To study the effect of artificial damping coefficient, analyses are performed ranging from 0.1 to 0.000005 with a constant value of {\tau }_{\max } equal to 100 MPa. The results of analyses are presented in Fig. 7. It is observed that as μ decreases, debonding length increases. However, it is not as direct an increase in debonding length as seen in u_t^c analyses in section 4.2. The “bottom” graph in Fig. 7 shows the steep increase in debonding length. From Fig. 7 a small interval of 0.001 to 0.00005 results in a significantly large change in debonding length. Figure 7 clearly reflects the importance of selection of an appropriate μ value for a particular debonding analysis.

Figure 8 shows the effect of artificial damping coefficient on the reaction force of contact elements. Figure 8 is in inverse correlation to the “bottom” graph of Fig. 7. Similar trend is observed between Figs. 3 and 4. The sharp increase in debonding length in the “bottom” graph of Fig. 7 represents the sudden weakening of the interface between the fiber and matrix when μ is decreased, especially after 0.001. This weakening interface relates to the weak contact elements that fail with decreasing μ. As μ decreases, the contact elements become weak and hence the force required for failure also decreases. This phenomenon can be observed in Fig. 8, where the force reaction takes a similar trend as plotted in the “bottom” picture of Fig. 7. Hence, this value of μ = 0.001 is considered as an optimal value, for the considered problem after which the structure fails.

Figure 9 shows a comparison of the debonding length with decreasing {\tau }_{\max } and a lower value of μ = 0.0005 than in Fig. 3. The value of μ is 0.0005 to check the effect on debonding length for decreasing values of {\tau }_{\max } as plotted. The effect of {\tau }_{\max } and u_t^c on debonding length, is found to be in similar both Figs. 6 and 9. For a lower initial value of artificial damping coefficient μ = 0.0005, at the first point of the graph ({\tau }_{\max } = 150 MPa), there is no significant increase in debonding length. But as {\tau }_{\max } decreases further, rate of increase of debonding length increases. From the trend observed in Fig. 9, it can be stated that drastically higher debonding length would be observed if value of μ is chosen as 0.0005.

Conclusions

Maximum tangential contact stress, tangential slip distance and artificial damping coefficient and their interdependency have a significant effect on the debonding length at fiber-matrix interface. Each of these three parameters initiate debonding at three different values at varying lengths of debonding. These values differ for different problems and hence, a case by case study for finding these minimum values is recommended. Artificial damping coefficient is of particular importance as its effect on debonding length is more drastic at significantly lower damping values. In contrast, a small decrease in tangential slip distance results in a higher debonding compared to artificial damping at similar values. Hence, wide range of artificial damping and tangential slip values should be studied along with their interdependency due to their inverse effect on debonding. Case wise detailed study of these parameters is recommended as they play a critical role in determining the stability of the composite.

References

[1] Pattabhi Ramaiah B, Rammohan B, Vijay Kumar S, Satish Babu D, Raghuatnhan R. Aero-elastic analysis of stiffened composite wing structure. Advances in Vibration Engineering, 2009, 8(3): 255–264

[2] Sudhir Sastry Y B, Budarapu P R, Madhavi N, Krishna Y. Buckling analysis of thin wall stiffened composite panels. Computational Materials Science, 2015, 96(B): 459–471

[3] Frolov V A. Strength of a composite material for structural applications. Mechanics of Composite Materials, 1987, 23(2): 148–154

[4] Sudhir Sastry Y B, Budarapu P R, Krishna Y, Devaraj S. Studies on ballistic impact of the composite panels. Theoretical and Applied Fracture Mechanics, 2014, 72: 2–12

[5] BPO'Rourke. The uses of composite materials in the design and manufacture of formula 1 racing cars. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering, 1990, 204: 41–48

[6] Parga-Landa B, Hernández-Olivares F. Analytical model to predict behaviour of soft armours. International Journal of Impact Engineering, 1995, 16(3): 455–466

[7] Anderson C E Jr, Bodner S R. Ballistic impact: The status of analytical and material modeling. International Journal of Impact Engineering, 1988, 7(1): 9–35

[8] Naik N K, Shrirao P. Composite structures under ballistic impact. Composite Structures, 2004, 66(1-4): 579–590

[9] James C. Leslie, The use of composite material increase the availability of oil and gas and reduce the cost of drilling and production operations. Proc. SPIE 6531, Nondestructive Characterization for Composite Materials, Aerospace Engineering, Civil Infrastructure, and Homeland Security, 65310E, 2007,

[10] Anderson C E Jr, Bodner S R. Ballistic impact: The status of analytical and material modeling. International Journal of Impact Engineering, 1988, 7(1): 9–35

[11] Naik N K, Shrirao P. Composite structures under ballistic impact. Composite Structures, 2004, 66(1−4): 579–590

[12] Katz S, Grossman E, Gouzman I, Murat M, Wiesel E, Wagner H D. Response of composite materials to hypervelocity impact. International Journal of Impact Engineering, 2008, 35(12): 1606–1611

[13] Budarapu P R, Narayana T S S, Rammohan B, Rabczuk T. Directionality of sound radiation from rectangular panels. Applied Acoustics, 2015, 89: 128–140

[14] Scoutis C. Carbon fiber reinforced plastics in aircraft construction. Materials Science and Engineering: A, 2005, 412(1−2): 171–176

[15] Paipetis A S. Room vs. temperature studies of model composites: modes of failure of carbon fibre/epoxy interfaces. Composite Interfaces, 2012, 19(2): 135–158

[16] Xu Z, Li J, Wu X, Huang Y, Chen L, Zhang G. Effect of kidney-type and circular cross sections on carbon fiber surface and composite interface. Composites. Part A, Applied Science and Manufacturing, 2008, 39(2): 301–307

[17] Zhao J, Ho K K C, Shamsuddin S R, Bismarck A, Dutschk V. A comparative study of fibre/matrix interface in glass fibre reinforced polyvinylidene fluoride composites. Colloids Surfaces A Physicochem, Eng Asp, 2012, 413: 58–64

[18] Nath R B, Fenner D N, Galiotis C. Progressional approach to interfacial failure in carbon reinforced composites: Elasto-plastic finite element modelling of interface cracks. Composites. Part A, Applied Science and Manufacturing, 2000, 31(9): 929–943

[19] Kim B W, Nairn J A. Observations of fiber fracture and interfacial debonding phenomena using the fragmentation test in single fiber composites. Journal of Composite Materials, 2002, 36(15): 1825–1858

[20] Ho H, Drzal L T. Non-linear numerical study of the single-fiber fragmentation test Part I: Test mechanics. Composites Engineering, 1995, 5(10−11): 1231–1244

[21] Yang S W, Budarapu P R, Roy Mahapatra D, Bordas S, Rabczuk T. A meshless adaptive multiscale method for fracture. Computational Materials Science, 2015, 96(B): 382–395

[22] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modelling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[23] Budarapu P R, Gracie R, Bordas S P A, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[24] Sudhir Sastry Y B, Krishna Y, Budarapu P R. Parametric studies on buckling of thin walled channel beams. Computational Materials Science, 2015, 96(B): 416–424

[25] Budarapu P R, Yb S S, Javvaji B, Mahapatra D R. Vibration analysis of multi-walled carbon nanotubes embedded in elastic medium. Frontiers of Structural and Civil Engineering, 2014, 8(2): 151–159

[26] Xiao S P, Belytschko T. A bridging domain method for coupling continua with molecular dynamics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(17−20): 1645–1669

[27] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2008, 197(6−8): 641–658

[28] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid−structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[29] Robert G, Belytschko T. Concurrently coupled atomistic and XFEM models for dislocations and cracks. International Journal for Numerical Methods in Engineering, 2008, 78: 354–378

[30] Rabczuk T, Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[31] Rabczuk T, Bordas S P, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[32] Bordas S P, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[33] Gracie R, Belytschko T. Adaptive continuum-atomistic simulations of dislocation dynamics. International Journal for Numerical Methods in Engineering, 2011, 86(4−5): 575–597

[34] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29−30): 2777–2799

[35] Rabczuk T, Areias P M A, Belytschko T. A simplified meshfree method for shear bands with cohesive surfaces. International Journal for Numerical Methods in Engineering, 2007, 69(5): 993–1021

[36] Rabczuk T, Song J H, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics, 2009, 76(6): 730–741

[37] Rabczuk T, Zi G, Bordas S P, Nguyen-Xuan H. A simple and robust threedimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37−40): 2437–2455

[38] Turon A, Camanho P P, Costa J, Dávila C G. A damage model for the simulation of delamination in advanced composites under variable-mode loading, Mechanics of Materials, 2006, 38(11): 1072–1089 

[39] Irwin G R. Fracture. Handbuch der Physik. Flugge S, ed. Berlin: Springer, 1958, VI: 551–590

[40] Rybicki E F, Kanninen M F. A finite element calculation of stress intensity factors by a modified crack closure integral. Engineering Fracture Mechanics, 1977, 9(4): 931–938

[41] Shivakumar K N. Tan,PW. Newman,  JC. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. International Journal of Fracture, 1988, 36: R43–R50

[42] Krueger R. Virtual crack closure technique: History, approach, and applications. Applied Mechanics Reviews, 2004, 57(2): 109–143

[43] Raju I S. Calculation of strain-energy release rates with higher order and singular finite elements. Engineering Fracture Mechanics, 1987, 28(3): 251–274

[44] Rice J R. A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics, 1968, 35(2): 379–386

[45] Hellen T K. On the method of virtual crack extensions. International Journal for Numerical Methods in Engineering, 1975, 9(1): 187–207

[46] Parks D M. A stiffness derivative finite element technique of crack tip stress intensity factors. International Journal of Fracture, 1974, 10(4): 487–502

[47] Griffith A A. The phenomenon of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1921, 221(582−593): 163–198

[48] Turon A, Costa J, Camanho P P, Davila, C G. Simulation of delamination propagation in composites under high-cycle fatigue by means of cohesive-zone models, NASA/TM 214532 1− 28, 2006

[49] Turon A, Camanho P P, Costa J, Dávila C G. Mechanics of Materials, 2006, 38: 1072–1089

[50] Turon A, Costa J, Camanho P P, Dávila C G. Composites. Part A, Applied science and manufacturing, 2008, 38: 2270–2282

[51] Kumar S, Singh I V, Mishra B K. A homogenized XFEM approach to simulate fatigue crack growth problems. Computers & Structures, 2015, 150: 1–22

[52] Kumar S, Singh I V, Mishra B K, Rabczuk T. Modeling and simulation of kinked cracks by virtual node XFEM. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 1425–1466

[53] Kumar S, Singh I V, Mishra B K. A multigrid coupled (FE-EFG) approach to simulate fatigue crack growth in heterogeneous materials. Theoretical and Applied Fracture Mechanics, 2014, 72: 121–135

[54] Agarwal A, Singh I V, Mishra B K. Numerical prediction of elasto-plastic behaviour of interpenetrating phase composites by EFGM. Composites. Part B, Engineering, 2013, 51: 327–336

[55] Bhardwaj G, Singh I V, Mishra B K, Bui T Q. Numerical simulation of functionally graded cracked plates using NURBS based XIGA under different loads and boundary conditions. Composite Structures, 2015, 126: 347–359

[56] Bhardwaj G, Singh I V, Mishra B K. Stochastic fatigue crack growth simulation of interfacial crack in bi-layered FGMs using XIGA. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 186–229

[57] Dugdale D S. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, 1960, 8(2): 100–104

[58] Barenblatt G I. The mathematical theory of equilibrium cracks in brittle f<?Pub Caret?>racture. Advances in Applied Mechanics, 1962, 7(C): 55–129

[59] Bažant Z P, Jirásek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. Journal of Engineering Mechanics, 2002, 128(11): 1119–1149

[60] Alfano G, Crisfield M A. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. International Journal for Numerical Methods in Engineering, 2001, 50(7): 1701–1736

[61] Fujimoto T, Kagami J, Kawaguchi T, Hatazawa T. Micro-displacement characteristics under tangential force. Wear, 2000, 241(2): 136–142

[62] Grzemba B, Pohrt R, Teidelt E, Popov V L. Maximum micro-slip in tangential contact of randomly rough self-affine surfaces. Wear, 2014, 309(1−2): 256–258

[63] Design M, Section C S, Technical N, Introduction I. An artificial damping method for the harmonically excited non-linear systems. 1988, 120: 597–608

[64] Kim S. Artificial damping in multigrid methods. Applied Mathematics Letters, 2001, 14(3): 359–364

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg