J-integral for delaminated beam and plate models

In this work the J-integral is applied to calculate the energy release rate in simple beam and plate models. Four examples are considered: the mode-I double-cantilever beam, the modeII end-notched flexure, the mixed-mode I/II single-leg bending specimens and a delaminated plate with simply-supported edges, respectively. In each example the details of calculation are given, in the latter case the distribution of the energy release rate along the crack front is calculated. While for delaminated beams the literature presents the energy release rates in numerous studies, the J-integral calculation is not trivial in theses cases. Moreover for delaminated bent plates the application of the J-integral is not documented. It is shown that considering the classical plate theory there are serious limitations to calculate representative energy release rates.


Introduction
The J-integral was developed in 1968 by Rice [1] to characterize the strain concentration around cracks and notches.The original definition of the integral is: where the parameters will be defined later.Essentially, 2D problems were presented including elastic and elastic-plastic field theories.Later, the method was extended for orthotropic composite materials [2,3] and 3D problems too [4][5][6].In the latter case the so-called J k vector was defined: where according to Fig. 1 n k is the outward normal vector of  the contour C, δ i j is the Kronecker tensor, σ i j is the stress tensor (σ i j n j is the traction vector), u i is the displacement vector, A is the area enclosed by contour C. The contour C contains the crack tip and the integration is carried out in the counterclockwise direction (see Fig. 1).Under static conditions the J-integral is equivalent to the energy release rate, G, which is also known as the crack driving force [7].The original brittle fracture model of Griffith [7] applies the relation to determine G, which is useful if the dependence of the strain energy, U and work of external forces, W F are known in the function of crack length a (dA = bda, where b is the width).However, for problems, where the mentioned relationships with the crack length are difficult to determine the J-integral is very effective.The energy release rate can be obtained using field theories, i.e. the stress and strain states are to be calculated analytically or numerically.In this paper some basic problems are analyzed, which are in fact simple, but -to telling the truthnot available in the literature and could be useful for those who face to the J-integral application for the first time.Four examples are considered here: the mode-I double-cantilever beam (DCB) [8][9][10][11][12][13][14][15][16], the mode-II end-notched flexure (ENF) specimen [14,15,17,18] and the mixed-mode I/II single-leg bending (SLB) system [14,15,17,18].Finally, a simply-supported delaminated elastic plate subjected to a point force is analyzed.In the latter problem a three dimensional analysis is necessary and some nontrivial aspects of the calculation are presented.

Example 1. -mode-I double-cantilever beam
The DCB specimen is a very well known, common specimen [11][12][13][14][15][16], schematically shown in Fig. 2. We analyze two loading schemes: DCB loaded in pure bending and subjected to a point force.Although mainly the second case is involved in practical tests, even the first case has practical importance, see e.g.: [10].

DCB in pure bending
Fig. 2 shows the DCB subjected to pure bending.The contour consists of four parts.The corresponding stress and displacement components can be obtained by simple beam theory [7]: (1) : (2) : (3) : (4) : where b is the specimen width, h is the thickness, I is the area moment of inertia, E is the flexural modulus.The strain energy density, W is nonzero only for contour No. 3: The integral of the strain energy density becomes: where the arc length coordinate is s = −y.The second term in the integral is nonzero again only for contour No. 3: The sum of the latter two expressions gives: Considering the fact that we have analyzed only a half model, we have: which is well-known from the literature [11][12][13][14][15][16].

DCB loaded by concentrated force
The DCB loaded by a wedge force is considered in Fig. 2b.The stress and displacement components for the contour parts based on simple beam theory are: (1) : (2) : (3) : (4) : The integral of strain energy density along the through-thickness direction is zero even for contour No. 3: which can be explained by the fact that transverse shear deformation is not considered.The second term in Eq. (2) becomes: Again, taking two times the result of the the half model, we have: Per. Pol.Mech.Eng.

64
András Szekrényes    The ENF specimen is a standard beam for the measurement of the mode-II toughness of composites [14][15][16][17][18].The specimen geometry is depicted by Fig. 3.For the calculation we define a zero-area path [23] around the crack tip.The dashed lines show the stress and displacement fields on both sides of the tip.The problem is symmetric with respect to the x-axis, therefore we consider the top half of the beam only.The outward normal vectors of contours No. 1 and 2 are:  The stress tensor, traction vector, displacement vector field and its derivative for the right and left contour parts, respectively are: (1) : (2) : where v=v(x) the transverse deflection of the beam.The strain energy densities on both sides of the crack tip are given by: where are the moments reduced to the crack tip.
The integral of strain energy densities in the thickness direction results in: The second part of the J-integral is: The sum of Eqs. ( 26) and ( 27) leads to: Taking Eq. ( 28) two times the energy release rate (ERR) becomes: Apparently, this simple formula is well-known in the literature [14,15,17,18].

Example 3. -the single-leg bending specimen
As a mixed-mode I/II fracture system we consider the SLB specimen [15,[19][20][21][22], shown by Fig. 4. As it is seen, the bottom arm is unloaded.It is reasonable to calculate the mode-I and mode-II ERRs separately.A mode decomposition method has been proposed by Shivakumar and Raju [5] which is based  on the separation of the displacement and stress components into symmetric and antisymmetric parts.Later, it was shown by Rigby and Aliabadi [6] that the stress decomposition in [5] was partly incorrect and the method has been revisited, which was applied later by numerous authors, e.g.: [24,25].Fig. 4 shows that in the case of the SLB specimen the stress state can be decomposed by superimposing a DCB loaded by a bending moment Pa/2 and an ENF specimen, shown by Fig. 3.For this simple beam problem the method of Shivakumar and Raju [5] is equivalent to the global method by Williams [26].Thus we have, based on Eqs.( 13) and ( 29) the followings: which has already been published based on strength of materials analysis [15,[19][20][21][22].

Example 4. -a simply-supported delaminated elastic plate
In this section we present the solution for a simply-supported delaminated plate subjected to a point force.The plate is symmetrically delaminated, hence -since there is no delamination opening -the plate is under mixed-mode II/III fracture condition.Let us consider now Fig. 5, which is the original problem.In fact this problem can be treated as a stepped thickness plate [27,28] as it is shown by Fig. 5b.Due to the fact that the bending stiffness is piecewise continuous the only alternative is Lévy plate formulation [29][30][31][32] to solve the problem.

Diplacement field -Lévy solution
To solve the problem we adopt Lévy plate formulation [29,32].The deflection functions and area moments of inertia of the uncracked and delaminated parts are given by: In the y direction the deflections are approximated by sin functions, i.e. we have: The governing equation of isotropic, elastic plates is [29]: where p is the surface pressure and E 1 = E/(1 − ν 2 ).Taking Eqs.( 31) and (32) back into (39) yields: where [32]: The solution functions must satisfy the following boundary conditions: Per. Pol.Mech.Eng.

J-integral evaluation
To calculate the J-integral -similarly to the ENF system -we define a zero-area path around the crack tip in accordance with Fig. 6.The problem is symmetric with respect to the midplane, therefore we analyze only the top half of the plate.The outward normals are given by Eq.( 21).This problem involves the 3D Jintegral, therefore, Eqs.( 2) and (3) must be used.However, due to the zero-area path, the surface integral term in Eq.(2) becomes zero.The stress tensor, traction vector, displacement vector and its derivative for contours No.  (1) : (2) : J-integral for delaminated beam and plate models 67 2012 56 1 Thus, the J 1 integral becomes: Since there is no crack opening, we have: Finally, J 3 is: In Eqs.(50) and (52) the terms related to σ y and y in J 1 and J 3 should be ignored.This can be justified by the incompatible displacement field illustrated in Fig. 7. Theoretically, if there is no crack opening then at the point where the cracked and uncracked parts are connected to each other the strain in the y direction and the distribution of the stress σ y are the same, in other words there is no discountinuity between the cracked and uncracked parts.Therefore, in the ideal case the integration of the stress, σ y producted by the strain in the y direction over the thickness would lead to zero.Consequently, the Kirchhoff plate model predicts erroneously the stress state in the transition zone, which can be counteracted only by ignoring the terms in question.Moreover, in accordance with Fig. 6 it can be seen that M xy is a negative stress couple, while M yx is positive, otherwise they have the same magnitude i.e.: M yx = −M xy .It can also be noticed in Fig. 6, that at an actual point they amplify each other.Thus, the J-integrals for the delaminated plate are: Next, we need to separate J 1 and J 3 into mode-II and mode-III J-integrals.There is a direct decomposition method [5,6], which provides the mode-I, mode-II and mode-III J-integrals in the form of: i.e. the J III -integral is directly obtained from J 3 .In our case the mode-I component is zero, thus we have: The result of Kirchhoff plate theory is compared to finite element results carried out by ANSYS 12, detailed in next section.is shown by Fig. 8.The model was created in the commercial ANSYS 12 package using 8 node solid elements.In the vicinity of the crack tip a refined mesh was constructed including trapezoid shape elements [19].The z displacements of contact nodes over the delaminated surface were imposed to be the same.The mode-II and mode-III ERRs were calculated by the virtual crack-closure technique (VCCT, e.g.: [34]), the size of the crack tip elements were ∆x = 0.25 mm, ∆y = 0.25 mm and ∆z = 2 mm.For the determination of G II and G III a socalled MACRO was written in the ANSYS Design and Parametric Language (ADPL).The MACRO gets the nodal forces and displacements at the crack tip and at each pair of nodes, respectively, then by defining the size of crack tip elements it determines and plots the ERRs at each point along the crack front.Apart from the VCCT G II and G III were determined also by the the J-integral, which is available in ANSYS 12 as a built-in command.

Comparison of analytical and numerical results
The mode-II, mode-III ERRs and the mode ratio along the crack front are shown by Fig. 9.The symbols show the result of the VCCT and the J-integral by FEM, while the curves represent the purely analytical plate theory solutions.It is seen that the mode-II component is significantly underestimated by the plate theory solution.On the contrary the mode-III component is overestimated.While the VCCT predicts that the mode-III ERR decays suddenly near the edges, there is not any decay in accordance with the plate solution.From the practical point of view this difference is insignificant.The crack is expected to initiate/propagate at the points where the highest ERR takes place.This point is in fact almost identical predicted by the two methods.The ANSYS' J-integral underpredicts both ERR components in comparison with the VCCT.In this respect, the plate theory solution agrees more or less with the mode-II ERR, on the contrary the mode-III component by the numerical J-integral is significantly underpredicted and neither the plate theory nor the VCCT shows good agreement with it.In Fig. 9a the curve marked by G My represents that part of the J-integral, which was calculated by the bending moment M y and the corresponding strains in the y direction.This curve shows extremely large discrepancy compared to the FE results.This confirms that Kirchhoff plate theory is insufficient in this respect.The mode ratio (G II /G III ) is depicted in Figure 9b.It shows that the VCCT and the numerical J-integral is in a reasonably good agreement with each other.The plate theory solution slightly underestimates the mode ratio, however the nature of the curves matches well with the numerical results.The overall agreement between the VCCT and plate theory methods is acceptable.The difference between the VCCT and plate theory solution can be attributed to the transverse and interlaminar shear effects.Further work is necessary to improve the plate theory solution.
Considering the possible application of plate theory solution it must be noted, that nowadays more and more attention is favored to the mode-III and mixed-mode II/III or I/III fracture of composite materials.Recently, de Morais and Pereira published many studies for carbon/epoxy systems involving the bending of plate specimens.The mode-III 4-point bending plate [35], the mixed-mode I/III 8-point bending [36] and the mixed-mode II/III 6-point bending plate [37] specimens involve the bending of delaminated plates.Analytical reduction techniques are not yet available for these systems.As a future work, by using a proper plate theory this problem is planned to be solved.

Conclusions
The application of J-integral has been demonstrated including common fracture specimens and a delaminated simply-supported isotropic plate.While for cracked beams the Jintegral is relatively simple to be calculated, in the plate theory solution some inconsistencies have been elaborated.First, plate have curvatures in the y direction too, and it has been shown that the transition zone between the cracked and uncracked parts are erroneously predicted by the plate model.This can be counteracted by ignoring the stress and displacement components in the y direction when we calculate the J-integral.Second, in the mode-III component it must be considered, that the twisting moments in the different planes amplify each other, consequently in the J-integral their sum appears.Otherwise, based on comparison with finite element calculations plate theory seems to be suitable to calculate energy release rate in delaminated plates, however higher order plate theories may help to achieve better accuracy.

Figure 1 :
Figure 1: Reference system for the 3D J-integral.

Figure 2 :Figure 3 :
Figure 2: Double-cantilever beam loaded in pure bending (a) and by a wedge force (b).

Fig. 2 .
Fig. 2. Double-cantilever beam loaded in pure bending (a) and by a wedge force (b).

Figure 1 :
Figure 1: Reference system for the 3D J-integral.

Figure 4 :Figure 5 :
Figure 4: Integration path and stress decomposition for the single-leg bending specimen.

Fig. 4 .
Fig. 4. Integration path and stress decomposition for the single-leg bending specimen.

Figure 6 :Figure 7 : 1 crackFigure 8 :
Figure 6: Integration path and stress couples around the crack tip for an elastic p

Fig. 6 .
Fig. 6.Integration path and stress couples around the crack tip for an elastic plate.

Figure 6 :Figure 7 :Fig. 7 .
Figure 6: Integration path and stress couples around the crack tip for an elastic plate.

Figure 7 : 1 crackFigure 8 :
Figure 7: Incompatible displacement field due to M y in a delaminated elastic plate.

Fig. 9 .
Fig. 9. Energy release rate (a) and mode ratio (b) for an elastic delaminated plate.