Interface fracture in orthotropic composite plates using second-order shear deformation theory

Abstract

The second-order shear deformation theory is used in this study to calculate the stresses and the energy release rates in orthotropic composite plates. A novel double-plate system is utilized with the imposition of the proper kinematic constraints in the interface plane of a double-plate system. The governing equations of the system were derived and as a demonstrative example a simply-supported plate subjected to a point force was analyzed. Using Lévy plate formulation, the plate problem was solved by a state-space model, incorporating four different regions. The distribution of the stress resultants and the interlaminar stresses in the uncracked part were also determined. Moreover, the distributions of the mode-II and mode-III energy release rates along the crack front were calculated by the J-integral. The 3D finite element model of the plate was created providing reference data for the analytical model. The results show reasonably good agreement between the analytical and numerical results. Also, the present model eliminates the physical inconsistency of previous models and reveals that under mixed-mode II/III condition, the energy release rate is not contributed by the bending, twisting moments and shear forces at all.

Keywords

Laminated plate delamination J-integral mixed-mode II/III fracture second-order plate theory interface constraint

Introduction

Nowadays laminated composite materials play a very significant role in the industrial practice as well as in the scientific life. The numerous application fields (e.g.: bodywork construction, cars, airplanes and aircrafts) show that the development of these heterogeneous materials continues probably very intensively in the future. One of the many damage modes of composites (Oskay and Pal, 2010) is interlaminar fracture or delamination (Anderson, 2005; Lundmark and Varna, 2011; Nagarajan et al., 2012). The resistance against the delamination under static and cyclic (e.g., Gornet and Ijaz, 2011) loading is characterized by the energy release rate (ERR). The basic concept of the linear elastic fracture mechanics (LEFM) is that the crack initiates/propagates if the critical energy release rate (CERR) is reached in the material (Anderson, 2005). The three basic fracture modes are the mode-I (opening), mode-II (in-plane shear) and mode-III (anti-plane shear). Under these fracture modes, the CERR is determined through standard or nonstandard test methods. Apparently, the major part of the literature deals with mode-I (e.g., Jumel et al., 2011; Peng et al., 2011), mode-II (e.g., Argüelles et al., 2011; Plain and Tong, 2011; Wang and Qiao, 2004a) and mixed-mode I/II (e.g., Davidson and Sundararaman, 1996; Plain and Tong, 2011; Reeder and Crews, 1990) fractures. However, during the last decades more and more attention was favored to the mode-III fracture (Browning et al., 2010, 2011; Lee, 1993; Li and O’Brien, 1996; Mehrabadi and Khosravan, 2013; Szekrényes, 2009a). In contrast, the investigation of the mixed-mode I/III (Pereira and de Morais, 2009; Szekrényes, 2009b), II/III (de Morais and Pereira, 2008; Kondo et al., 2010, 2011; Mehrabadi, 2012; Miura et al., 2012; Suemasu et al., 2009, 2010; Szekrényes, 2007, 2012a) and I/II/III (Davidson and Sediles, 2011; Davidson et al., 2010; Szekrényes, 2011) delamination fracture of composite materials is related only to the 21st century. The first attempts include beam and plate specimens. These tests work more or less fine, however, much more effort is necessary to develop so simple and effective test methods like those for mode-I and mode-II. Compared to mode-I and mode-II, the mode-III fracture involves significant difficulties: pure mode-III fracture does not exist; the geometry of the samples is also a critical point. Beam specimens are in general very stiff, plate-like specimens are much more difficult to manufacture.

Independently of the fact whether we use beam or plate-like specimens, an analytical solution – in general – makes the data reduction relatively simple. For beam-like specimens many improved models have been developed (e.g., Jumel et al., 2011; Wang and Qiao, 2004b; Yazdi and Rezaeepazhand, 2012) based on elastic foundation beams, crack tip shear deformation and similar considerations. For plates the analytical solution is much more difficult to obtain (Sriram and Armanios, 1993; Tian and Fu, 2010): such solutions exist only for some relatively simple systems, like the edge crack torsion (ECT) specimen that involves simple loading conditions and analytical solution (de Moura et al., 2009; Lee, 1993). In the last few years, several plate-like specimens were developed for the mode-III (de Morais and Pereira, 2009), mixed-mode I/III (Pereira and de Morais, 2009) and II/III (de Morais and Pereira, 2008) fracture testing of laminated composite materials. Without any exception the data reduction is always made by the finite element (FE) method incorporating the virtual-crack closure technique (VCCT) and cohesive zone model (CZM) applications (e.g., Omiya and Kishimoto, 2010). For delaminated plates, Davidson et al. (2000) applied shell elements to calculate the ERRs in plate-like structures, Sankar and Sonik (1995) performed similar computations. The crack tip force method (CTFM; Park and Sankar, 2002) is a similar solution to the VCCT, utilizing the crack tip forces to calculate and separate the ERRs. However, its result does not differ from that of a VCCT analysis. The main problem of the FE models is that a 3D model is necessary to construct and the VCCT is not available as a built-in command in most of the FE packages.

The main aim of this paper is to present the application of second-order plate theory (SSDT) to analyze delaminated plates and to eliminate the physical inconsistency of previous beam and plate solutions. Wang and Qiao (2004a, 2004b), Qiao and Wang (2004), Qiao and Chen (2011) and Chen (2011) applied several flexible joint models mainly for beams. The continuity of the displacement along the interface of a double beam system was ensured by considering the interface peel and shear stresses. This model was successfully applied to fracture and vibration problems too (Qiao and Chen, 2012), although its physical inconsistency is evident. First, the interface shear compliance is defined arbitrarily to obtain acceptable results. Moreover, the basic equations of Euler-Bernoulli or Timoshenko beams are utilized, but in the displacement field a second-order term is assumed apart from the constant and linear one. This model does not conform to the basic equations of the Euler-Bernoulli and Timoshenko beams; furthermore, a contradiction takes place when we apply the equations of linear elasticity. The flexible joint model was later extended to analyze delaminated plates (Szekrényes, 2012b, 2013), although its extension to asymmetrically delaminated plates is limited.

In this work, the SSDT (Khdeir and Reddy, 1999) is utilized to develop a mechanical model for the fracture of delaminated orthotropic composite plates with symmetrical lay-up and straight delamination front. First, the displacement field is formulated by imposing the interface constraints. Second, the basic equations of linear elasticity are applied to derive the strain and stress fields in elastic orthotropic composite plates. The present formulation does not incorporate arbitrarily defined parameters and it is shown that the developed model is physically consistent with the equations of linear elasticity. As an example a simply-supported plate subjected to a point force is analyzed applying the state-space model (Reddy, 2004). The distribution of the interlaminar stresses is calculated, moreover, the J-integral (Cherepanov, 1997; Rice, 1968) is utilized to determine the distribution of the mode-II and mode-III ERRs along the crack front. A FE model is also created and the numerical results are compared to those obtained from the analytical model. The good agreement obtained shows the usefulness of plate theories with convenient interface constraints in the delamination analysis.

Second-order plate theory – general formulation

The plate theory presented in this section is utilized to capture the displacement and stress fields in the delaminated portion of an elastic laminated orthotropic plate with symmetric lay-up presented in Figure 1. The thickness of the top and bottom plates is t. The assumed displacement field based on SSDT for elastic plates can be written as (Khdeir and Reddy, 1999; Shahrjerdi and Mustapha, 2011):

u_{top} (x, y, z) = u_{0} (x, y) + θ_{x} (x, y) \cdot z + φ_{x} (x, y) \cdot z^{2} v_{top} (x, y, z) = v_{0} (x, y) + θ_{y} (x, y) \cdot z + φ_{y} (x, y) \cdot z^{2}

(1)

for the top plate element of the delaminated portion in Figure 1(b), and:

u_{bot} (x, y, z) = - u_{0} (x, y) + θ_{x} (x, y) \cdot z - φ_{x} (x, y) \cdot z^{2} v_{bot} (x, y, z) = - v_{0} (x, y) + θ_{y} (x, y) \cdot z - φ_{y} (x, y) \cdot z^{2}

(2)

for the bottom plate element of the delaminated portion in Figure 1(b). It should be mentioned that the SSDT applies unsymmetrical distribution for the in-plane displacements u and v over the plate thickness of the delaminated plate portions. In this paper, only mixed-mode II/III problems are considered, i.e. the crack opening is zero. Therefore, the transverse deflection, w = w(x,y), is the same for both the top and bottom plates. Applying the basic equations of linear elasticity (Chou and Pagano, 1967) for an orthotropic plate we obtain the strain and stress fields, then by integrating the stresses over the thickness the stress resultants in terms of strains and shear strains become:

{{N} {M} {L}} = [[A] [0] [D] [0] [D] [0] [D] [0] [F]] {{ɛ^{(0)}} {ɛ^{(1)}} {ɛ^{(2)}}}, {{Q} {R}} = [[A] [0] [0] [D]] {{γ^{(0)}} {γ^{(1)}}}

(3)

where {N}^T= {N_x, N_y, N_xy} is the vector of in-plane normal and shear forces, {M}^T= {M_x, M_y, M_xy} is the vector of bending and twisting moments, {Q}^T= {Q_x, Q_y} is the vector of shear forces, furthermore {L}^T= {L_x, L_y, L_xy} and {R}^T= {R_x, R_y} are the vectors of higher order stress resultants (Khdeir and Reddy, 1999; Shahrjerdi and Mustapha, 2011). The stress resultants can be calculated from the integration of the stresses in the through-thickness direction (Reddy, 2004):

{N_{α β} M_{α β} L_{α β}} = \int_{- t / 2}^{t / 2} σ_{α β} {1 z z^{2}} dz, {Q_{α} R_{α}} = \int_{- t / 2}^{t / 2} σ_{α z} {1 z} dz

(4)

where α and β take the symbols x and y. The stiffness parameters are defined as (Kollár and Springer, 2003; Reddy, 2004):

A_{ij}, D_{ij}, F_{ij} = \sum_{k = 1}^{N_{l}} \int_{z_{k}}^{z_{k + 1}} {\bar{C}}_{ij}^{(k)} (1, z^{2}, z^{4}) dz

(5)

where the coordinates z_k and z_k+1 give the position of the k^th layer. In Figure 1, the relevant coordinates are indicated. The vectors of strains including constant, linear and quadratic terms are written as (Reddy, 2004):

{ɛ_{x} ɛ_{y} γ_{xy}} = {ɛ_{x}^{(0)} ɛ_{y}^{(0)} γ_{xy}^{(0)}} + z {ɛ_{x}^{(1)} ɛ_{y}^{(1)} γ_{xy}^{(1)}} + z^{2} {ɛ_{x}^{(2)} ɛ_{y}^{(2)} γ_{xy}^{(2)}}, {γ_{xz} γ_{yz}} = {γ_{xz}^{(0)} γ_{yz}^{(0)}} + z {γ_{xz}^{(1)} γ_{yz}^{(1)}}

(6)

where:

{ɛ_{x}^{(0)} ɛ_{y}^{(0)} γ_{xy}^{(0)}} = {\frac{\partial u_{0}}{\partial x} \frac{\partial v_{0}}{\partial x} \frac{\partial u_{0}}{\partial y} + \frac{\partial v_{0}}{\partial x}}, {ɛ_{x}^{(1)} ɛ_{y}^{(1)} γ_{xy}^{(1)}} = {\frac{\partial θ_{x}}{\partial x} \frac{\partial θ_{y}}{\partial x} \frac{\partial θ_{x}}{\partial y} + \frac{\partial θ_{y}}{\partial x}}, {ɛ_{x}^{(2)} ɛ_{y}^{(2)} γ_{xy}^{(2)}} = {\frac{\partial φ_{x}}{\partial x} \frac{\partial φ_{y}}{\partial x} \frac{\partial φ_{x}}{\partial y} + \frac{\partial φ_{y}}{\partial x}}

(7)

and:

{γ_{xz}^{(0)} γ_{yz}^{(0)}} = {θ_{x} + \frac{\partial w}{\partial x} θ_{y} + \frac{\partial w}{\partial y}}, {γ_{xz}^{(1)} γ_{yz}^{(1)}} = 2 {φ_{x} φ_{y}}

(8)

Formulating the total potential energy for the elastic plate subjected to transverse load and deriving the Euler-Lagrange equations, we obtain the equilibrium equations of the second-order plate in the following form (Shahrjerdi and Mustapha, 2011):

A_{11} \frac{\partial^{2} u_{0}}{\partial x^{2}} + A_{66} \frac{\partial^{2} u_{0}}{\partial y^{2}} + (A_{12} + A_{66}) \frac{\partial^{2} v_{0}}{\partial x \partial y} + D_{11} \frac{\partial^{2} φ_{x}}{\partial x^{2}} + D_{66} \frac{\partial^{2} φ_{x}}{\partial y^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} φ_{y}}{\partial x \partial y} = 0

(9)

A_{22} \frac{\partial^{2} v_{0}}{\partial y^{2}} + A_{66} \frac{\partial^{2} v_{0}}{\partial x^{2}} + (A_{12} + A_{66}) \frac{\partial^{2} u_{0}}{\partial x \partial y} + D_{22} \frac{\partial^{2} φ_{y}}{\partial y^{2}} + D_{66} \frac{\partial^{2} φ_{y}}{\partial x^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} φ_{x}}{\partial x \partial y} = 0

(10)

D_{11} \frac{\partial^{2} θ_{x}}{\partial x^{2}} + D_{66} \frac{\partial^{2} θ_{x}}{\partial y^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} θ_{y}}{\partial x \partial y} - A_{55} (θ_{x} + \frac{\partial w}{\partial x}) = 0

(11)

D_{22} \frac{\partial^{2} θ_{y}}{\partial y^{2}} + D_{66} \frac{\partial^{2} θ_{y}}{\partial x^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} θ_{y}}{\partial x \partial y} - A_{44} (θ_{y} + \frac{\partial w}{\partial y}) = 0

(12)

D_{11} \frac{\partial^{2} u_{0}}{\partial x^{2}} + D_{66} \frac{\partial^{2} u_{0}}{\partial y^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} v_{0}}{\partial x \partial y} + F_{11} \frac{\partial^{2} φ_{x}}{\partial x^{2}} + F_{66} \frac{\partial^{2} φ_{x}}{\partial y^{2}} + (F_{12} + F_{66}) \frac{\partial^{2} φ_{y}}{\partial x \partial y} - 4 D_{55} φ_{x} = 0

(13)

D_{22} \frac{\partial^{2} v_{0}}{\partial y^{2}} + D_{66} \frac{\partial^{2} v_{0}}{\partial x^{2}} + (D_{12} + D_{66}) \frac{\partial^{2} u_{0}}{\partial x \partial y} + F_{22} \frac{\partial^{2} φ_{y}}{\partial y^{2}} + F_{66} \frac{\partial^{2} φ_{y}}{\partial x^{2}} + (F_{12} + F_{66}) \frac{\partial^{2} φ_{x}}{\partial x \partial y} - 4 D_{44} φ_{y} = 0

(14)

A_{55} (\frac{\partial θ_{x}}{\partial x} + \frac{\partial^{2} w}{\partial x^{2}}) + A_{44} (\frac{\partial θ_{y}}{\partial y} + \frac{\partial^{2} w}{\partial y^{2}}) = - q

(15)

where A_ij, D_ij and F_ij are stiffness parameters and q = q(x,y) is the external load of the plate.

Second-order plate theory with interface constraint

The second-order plate theory utilized in this section is based on an assumed displacement field including an interface constraint to formulate the model of the uncracked region of a delaminated orthotropic composite plate. The mathematical form of in-plane displacement components are the same as those given by equations (1) and (2), however we have to ensure the displacement continuity between the top and bottom plate elements, as it is shown by Figure 1. The interface constraint equations of the uncracked portion are:

u_{top} |_{z = - \frac{t}{2}} = u_{bot} |_{z = \frac{t}{2}}, v_{top} |_{z = - \frac{t}{2}} = v_{bot} |_{z = \frac{t}{2}}

(16)

These conditions make it possible to express φ_x and φ_y in terms of the remaining parameters:

φ_{x} = - \frac{4}{t^{2}} (u_{0} - \frac{t}{2} θ_{x}), φ_{y} = - \frac{4}{t^{2}} (v_{0} - \frac{t}{2} θ_{y})

(17)

Taking these back into equations (1) and (2), we obtain the displacement field satisfying the interface constraint conditions:

u_{top} (x, y, z) = u_{0} + θ_{x} \cdot z - \frac{4}{t^{2}} (u_{0} - \frac{t}{2} θ_{x}) \cdot z^{2}, v_{top} (x, y, z) = v_{0} + θ_{y} \cdot z - \frac{4}{t^{2}} (v_{0} - \frac{t}{2} θ_{y}) \cdot z^{2}

(18)

For the bottom plate, similar expressions can be derived. Due to the symmetric lay-up with respect to the x − y plane, we analyze only the top plate in the sequel. Based on equation (18), the strains and shear strains become:

{ɛ_{x}^{(0)} ɛ_{y}^{(0)} γ_{xy}^{(0)}} = {\frac{\partial u_{0}}{\partial x} \frac{\partial v_{0}}{\partial x} \frac{\partial u_{0}}{\partial y} + \frac{\partial v_{0}}{\partial x}}, {ɛ_{x}^{(1)} ɛ_{y}^{(1)} γ_{xy}^{(1)}} = {\frac{\partial θ_{x}}{\partial x} \frac{\partial θ_{y}}{\partial x} \frac{\partial θ_{x}}{\partial y} + \frac{\partial θ_{y}}{\partial x}}

(19)

and:

{ɛ_{x}^{(2)} ɛ_{y}^{(2)} γ_{xy}^{(2)}} = - c_{1} {\frac{\partial u_{0}}{\partial x} - \frac{t}{2} \frac{\partial θ_{x}}{\partial x} \frac{\partial v_{0}}{\partial y} - \frac{t}{2} \frac{\partial θ_{y}}{\partial y} \frac{\partial u_{0}}{\partial y} + \frac{\partial v_{0}}{\partial x} - \frac{t}{2} (\frac{\partial θ_{x}}{\partial y} + \frac{\partial θ_{y}}{\partial x})}

(20)

moreover:

{γ_{xz}^{(0)} γ_{yz}^{(0)}} = {θ_{x} + \frac{\partial w}{\partial x} θ_{y} + \frac{\partial w}{\partial y}}, {γ_{xz}^{(1)} γ_{yz}^{(1)}} = - 2 c_{1} {u_{0} - \frac{t}{2} θ_{x} v_{0} - \frac{t}{2} θ_{y}}

(21)

where c₁ = 4/t². The system of equilibrium equations consists of only five equations in this case:

\frac{\partial N_{x}}{\partial x} + \frac{\partial N_{xy}}{\partial y} - c_{1} (\frac{\partial L_{x}}{\partial x} + \frac{\partial L_{xy}}{\partial y}) + 2 c_{1} R_{x} = 0,

(22)

\frac{\partial N_{xy}}{\partial x} + \frac{\partial N_{y}}{\partial y} - c_{1} (\frac{\partial L_{xy}}{\partial x} + \frac{\partial L_{y}}{\partial y}) + 2 c_{1} R_{y} = 0,

(23)

\frac{\partial M_{x}}{\partial x} + \frac{\partial M_{xy}}{\partial y} + \frac{c_{1} t}{3} (\frac{\partial L_{x}}{\partial x} + \frac{\partial L_{xy}}{\partial y}) - Q_{x} - {tc}_{1} R_{x} = 0,

(24)

\frac{\partial M_{xy}}{\partial x} + \frac{\partial M_{y}}{\partial y} + \frac{c_{1} t}{3} (\frac{\partial L_{xy}}{\partial x} + \frac{\partial L_{y}}{\partial y}) - Q_{y} - {tc}_{1} R_{y} = 0,

(25)

\frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y} + q = 0

(26)

Calculating the stress resultants (given by equation (3)) in terms of the basic parameters of the displacement field and taking them back into the equilibrium equations (equations (22)–(26)), the following system of equations is obtained:

a_{1} \frac{\partial^{2} u_{0}}{\partial x^{2}} + a_{2} \frac{\partial^{2} u_{0}}{\partial y^{2}} + a_{3} u_{0} + a_{4} \frac{\partial^{2} v_{0}}{\partial x \partial y} + a_{5} \frac{\partial^{2} θ_{x}}{\partial x^{2}} + a_{6} \frac{\partial^{2} θ_{x}}{\partial y^{2}} + a_{7} θ_{x} + a_{8} \frac{\partial^{2} θ_{y}}{\partial x \partial y} = 0

(27)

b_{1} \frac{\partial^{2} u_{0}}{\partial x \partial y} + b_{2} \frac{\partial^{2} v_{0}}{\partial x^{2}} + b_{3} \frac{\partial^{2} v_{0}}{\partial y^{2}} + b_{4} v_{0} + b_{5} \frac{\partial^{2} θ_{x}}{\partial x \partial y} + b_{6} \frac{\partial^{2} θ_{y}}{\partial x^{2}} + b_{7} \frac{\partial^{2} θ_{y}}{\partial y^{2}} + b_{8} θ_{y} = 0

(28)

c_{1} \frac{\partial^{2} u_{0}}{\partial x^{2}} + c_{2} \frac{\partial^{2} u_{0}}{\partial y^{2}} + c_{3} u_{0} + c_{4} \frac{\partial^{2} v_{0}}{\partial x \partial y} + c_{5} \frac{\partial^{2} θ_{x}}{\partial x^{2}} + c_{6} \frac{\partial^{2} θ_{x}}{\partial y^{2}} + c_{7} θ_{x} + c_{8} \frac{\partial^{2} θ_{y}}{\partial x \partial y} + c_{9} \frac{\partial w}{\partial x} = 0

(29)

d_{1} \frac{\partial^{2} u_{0}}{\partial x \partial y} + d_{2} \frac{\partial^{2} v_{0}}{\partial x^{2}} + d_{3} \frac{\partial^{2} v_{0}}{\partial y^{2}} + d_{4} v_{0} + d_{5} \frac{\partial^{2} θ_{x}}{\partial x \partial y} + d_{6} \frac{\partial^{2} θ_{y}}{\partial x^{2}} + d_{7} \frac{\partial^{2} θ_{y}}{\partial y^{2}} + d_{8} θ_{y} + d_{9} w = 0

(30)

e_{1} \frac{\partial θ_{x}}{\partial x} + e_{2} \frac{\partial θ_{y}}{\partial y} + e_{3} \frac{\partial^{2} w}{\partial x^{2}} + e_{4} \frac{\partial^{2} w}{\partial y^{2}} + q = 0

(31)

where the coefficients depend on the stiffness parameters and the thickness of the plate, and are collected in Appendix 1. It is important to highlight that the reference plane of the top and bottom plates of the uncracked portion are the individual midplanes (see Figure 1), and consequently the stiffness parameters in equation (5) for the delaminated and uncracked plate portions are the same. The different deformation of the two portions is included through the interface constraint conditions in equation (16).

In the next section, we solve a simply-supported plate subjected to a point force by the state-space model. The displacement and stress fields are calculated and the J-integral is utilized to calculate the energy release rate distributions along the crack front.

Example – simply-supported plate, Lévy plate formulation

In this section, we apply the state-space model (Reddy, 2004) to solve the system of equations for a delaminated plate subjected to a point force, shown by Figure 2. The governing partial differential equation (PDE) system is different for the delaminated and uncracked parts; therefore, the state-space models are developed separately.

Delaminated portion

In accordance with Lévy plate formulation, the displacement components and the external load for simply-supported second-order plates (Khdeir and Reddy, 1999) can be written as:

{u_{0} (x, y) v_{0} (x, y) θ_{x} (x, y) θ_{y} (x, y) φ_{x} (x, y) φ_{y} (x, y) w (x, y)} = \sum_{n = 1}^{\infty} {U_{0 n} (x) \sin β y V_{0 n} (x) \cos β y X_{n} (x) \sin β y Y_{n} (x) \cos β y T_{xn} (x) \sin β y T_{yn} (x) \cos β y W_{n} (x) \sin β y}, q = \sum_{n = 1}^{\infty} Q_{n} \sin β y

(32)

Taking these solutions back into the system of governing equations of the delaminated portion (equations (9)–(15)), we can derive the state-space model, in general form it becomes:

\underline{Z}'^{^{(d)} = {\underline{\underline{T}}}^{(d)} {\underline{Z}}^{(d)} + {\underline{F}}^{(d)}}

(33)

where the superscript “(d)” refers to the delaminated plate portion and Z is the state vector:

{\underline{Z}}^{(d)} = [U_{0 n} U'_{0 n} V_{0 n} V'_{0 n} X_{n} X'_{n} Y_{n} Y'_{n} T_{xn} T'_{xn} T_{yn} T'_{yn} W_{n} W'_{n}]^{T}

(34)

moreover,

\underline{\underline{T}}

is the system matrix, F contains the inhomogeneity of the governing PDE system, respectively. In an expanded form equation (33) becomes:

[U'_{0 n} U ″_{0 n} V'_{0 n} V ″_{0 n} X'_{n} X ″_{n} Y'_{n} Y ″_{n} T'_{xn} T ″_{xn} T'_{yn} T ″_{yn} W'_{n} W ″_{n}] = [01000000000000 k_{1} 00 k_{2} 0000 k_{3} 00 k_{4} 00000100000000000 m_{1} m_{2} 000000 m_{3} m_{4} 000000001000000000000 n_{1} 00 n_{2} 00000 n_{3} 0000000100000000000 p_{1} p_{2} 00000 p_{3} 000000000010000 r_{1} 00 r_{1} 0000 r_{3} 00 r_{4} 00000000000001000 s_{1} s_{2} 000000 s_{3} s_{4} 0000000000000000100000 t_{1} t_{2} 00000 t_{3} 0] [U_{0 n} U'_{0 n} V_{0 n} V'_{0 n} X_{n} X'_{n} Y_{n} Y'_{n} T_{xn} T'_{xn} T_{yn} T'_{yn} W_{n} W'_{n}] + [0000000000000 t_{4} Q_{n}]

(35)

where the coefficients in the 14 × 14 system matrix

\underline{\underline{T}}

are collected in Appendix 2. The general solution of equation (33) is (Reddy, 2004):

{\underline{Z}}^{(d)} (x) = e^{{\underline{\underline{T}}}^{(d)} x} [{\underline{K}}^{(d)} + \int_{x_{0}}^{x} e^{- {\underline{\underline{T}}}^{(d)} ξ} {\underline{F}}^{(d)} (ξ) d ξ] = {\underline{\underline{G}}}^{(d)} (x) {\underline{K}}^{(d)} + {\underline{H}}^{(d)} (x)

(36)

where K is the vector of constants (14), H is the vector of particular solutions.

Uncracked portion

It has been shown that because of the kinematic constraints, two of the displacement parameters can be eliminated, therefore for the uncracked part we have:

{u_{0} (x, y) v_{0} (x, y) θ_{x} (x, y) θ_{y} (x, y) w (x, y)} = \sum_{n = 1}^{\infty} {U_{0 n} (x) \sin β y V_{0 n} (x) \cos β y X_{n} (x) \sin β y Y_{n} (x) \cos β y W_{n} (x) \sin β y}, q = \sum_{n = 1}^{\infty} Q_{n} \sin β y

(37)

Utilizing equations (27)–(31), the state space model can be derived as:

\underline{Z}' (ud) = {\underline{\underline{T}}}^{(ud)} {\underline{Z}}^{(ud)} + {\underline{F}}^{(ud)}

(38)

where the superscript “(ud)” refers to the undelaminated region. The state vector in this case is:

{\underline{Z}}^{(ud)} = [U_{0 n} U'_{0 n} V_{0 n} V'_{0 n} X_{n} X'_{n} Y_{n} Y'_{n} W W'_{n}]^{T}

(39)

The expanded state space model becomes:

[U'_{0 n} U ″_{0 n} V'_{0 n} V ″_{0 n} X'_{n} X ″_{n} Y'_{n} Y ″_{n} W'_{n} W ″_{n}] = [0100000000 f_{1} 00 f_{2} f_{3} 00 f_{4} 0 f_{5} 00010000000 g_{1} g_{2} 00 g_{3} g_{4} 0 g_{5} 00000010000 h_{1} 00 h_{2} h_{3} 00 h_{4} 0 h_{5} 00000001000 j_{1} j_{2} 00 j_{3} j_{4} 0 j_{5} 0000000000100000 - \frac{e_{1}}{e_{3}} \frac{β e_{2}}{e_{3}} 0 \frac{β^{2} e_{4}}{e_{3}} 0] [U_{0 n} U'_{0 n} V_{0 n} V'_{0 n} X_{n} X'_{n} Y_{n} Y'_{n} W_{n} W'_{n}] + [000000000 - Q_{n} / e_{3}]

(40)

which involves a 10 × 10 system matrix with coefficients defined in Appendix 3. Considering the fact that there is no external load at the uncracked plate portion, the solution of equation (38) is:

{\underline{Z}}^{(ud)} (x) = e^{{\underline{\underline{T}}}^{(ud)} x} {\underline{K}}^{(ud)} = {\underline{\underline{G}}}^{(ud)} (x) {\underline{K}}^{(ud)}

(41)

Boundary and continuity conditions

The elements of the state vectors in equations (36) and (41) can be referred to as:

Z_{i}^{(d)} = \sum_{j = 1}^{14} G_{ij}^{(d)} K_{j}^{(d)} + H_{j}^{(d)}, Z_{i}^{(ud)} = \sum_{j = 1}^{10} G_{ij}^{(ud)} K_{j}^{(ud)} + H_{j}^{(ud)}

(42)

In accordance with Figure 2, we have four different plate portions. The point force causes singularity in the PDEs, therefore a plate portion loaded by a constant line force was applied, the length d₀ was a very small value compared to the plate dimensions. In this case, Q_n = 2q₀/b·sin(βy₀) (Reddy, 2004). Thus, the four parts are denoted by “1a”, “1q”, “1” for the delaminated portion and “2” for the undelaminated region. Consequently, the state-space model in section ‘Delaminated portion' is utilized for the “1a”, “1q” and “1” portions, while the one in section ‘Uncracked portion' was used to model the undelaminated “2” region. The boundary conditions (B.C.s) are formulated through the displacement parameters and the stress resultants. The latter ones can be expressed in the following forms:

{N_{x} N_{y} M_{x} M_{y}} = \sum_{n = 1}^{\infty} {n_{xn} n_{yn} m_{xn} m_{yn}} \sin β y, {L_{x} L_{y} Q_{x} R_{x}} = \sum_{n = 1}^{\infty} {l_{xn} l_{yn} q_{xn} r_{xn}} \sin β y, {N_{xy} M_{xy} Q_{y} R_{y}} = \sum_{n = 1}^{\infty} {n_{xyn} m_{xyn} q_{yn} r_{yn}} \cos β y

(43)

i.e. n_xn is the function coefficient in the Fourier series of N_x, etc. For the present problem, 54 conditions need be formulated. The B.C.s are:

W_{n}^{(1 a)} (a) = 0, V_{0 n}^{(1 a)} (a) = 0, Y_{n}^{(1 a)} (a) = 0, T_{yn}^{(1 a)} (a) = 0, n_{xn}^{(1 a)} (a) = 0, m_{xn}^{(1 a)} (a) = 0, l_{xn}^{(1 a)} (a) = 0, W_{n}^{(2)} (- c) = 0, V_{0 n}^{(2)} (a) = 0, Y_{n}^{(2)} (- c) = 0, n_{xn}^{(2)} (a) = 0, m_{xn}^{(2)} (- c) = 0

(44)

The continuity conditions between regions “1” and “2” are:

U_{0 n}^{(1)} (0) = U_{0 n}^{(2)} (0), X_{n}^{(1)} (0) = X_{n}^{(2)} (0), T_{xn}^{(1)} (0) = - \frac{4}{t^{2}} (U_{0 n}^{(2)} - \frac{t}{2} X_{n}^{(2)}) |_{x = 0} V_{0 n}^{(1)} (0) = V_{0 n}^{(2)} (0), Y_{n}^{(1)} (0) = Y_{n}^{(2)} (0), T_{yn}^{(1)} (0) = - \frac{4}{t^{2}} (V_{0 n}^{(2)} - \frac{t}{2} Y_{n}^{(2)}) |_{x = 0} W_{n}^{(1)} (0) = W_{n}^{(2)} (0), W {' (1)}_{n}^{^{}} (0) = W {' (2)}_{n} (0) n_{xn}^{(1)} (0) = n_{xn}^{(2)} (0), n_{xyn}^{(1)} (0) = n_{xyn}^{(2)} (0), m_{xn}^{(1)} (0) = m_{xn}^{(2)} (0), m_{xyn}^{(1)} (0) = m_{xyn}^{(2)} (0)

(45)

The continuity conditions between regions “1q” and “1” are:

U_{0 n}^{(1 q)} (x_{01}) = U_{0 n}^{(1)} (x_{01}), X_{n}^{(1 q)} (x_{01}) = X_{n}^{(1)} (x_{01}), T_{xn}^{(1 q)} (x_{01}) = T_{xn}^{(1)} (x_{01}) V_{0 n}^{(1 q)} (0) = V_{0 n}^{(1)} (0), Y_{n}^{(1 q)} (0) = Y_{n}^{(1)} (0), T_{yn}^{(1 q)} (x_{01}) = T_{yn}^{(1)} (x_{01}) W_{n}^{(1 q)} (x_{01}) = W_{n}^{(1)} (x_{01}), W {' (1 q)}_{n} (x_{01}) = W {' (1)}_{n} (x_{01}) n_{xn}^{(1 q)} (x_{01}) = n_{xn}^{(1)} (x_{01}), n_{xyn}^{(1 q)} (x_{01}) = n_{xyn}^{(1)} (x_{01}), m_{xn}^{(1 q)} (x_{01}) = m_{xn}^{(1)} (x_{01}), m_{xyn}^{(1 q)} (x_{01}) = m_{xyn}^{(1)} (x_{01}), l_{xn}^{(1 q)} (x_{01}) = l_{xn}^{(1)} (x_{01}), l_{xyn}^{(1 q)} (x_{01}) = l_{xyn}^{(1)} (x_{01})

(46)

where x₀₁ = x_Q − d₀. Moreover, 14 continuity conditions can be formulated between regions “1a” and “1q” at x₀₂ = x_Q + d₀, but these are similar to those in equation (46), therefore they are not necessary to present. Totally, the system involves 54 equations, which is equal to the number of unknowns. The problem was solved in the code MAPLE 12. Based on the continuity conditions, the following stress resultants are continuous across the delamination front: M_x, M_y, M_xy, Q_x, Q_y, R_x and R_y.

Calculation of the J-integral

In the general 3D case, the J-integral is defined as (Rigby and Aliabadi, 1998; Shivakumar and Raju, 1992):

J_{k} = \int_{C} ({Wn}_{k} - σ_{ij} u_{i, k} n_{j}) ds - \int_{A} (σ_{i 3} u_{i, k})_{, 3} d A, k = 1, 2, 3

(47)

where W is the strain energy density:

W = \int_{0}^{ɛ_{ij}} σ_{ij} d ɛ_{ij} = \frac{1}{2} (σ_{x} ɛ_{x} + σ_{y} ɛ_{y} + τ_{xy} γ_{xy} + τ_{xz} γ_{xz} + τ_{yz} γ_{yz})

(48)

Moreover, based on Figure 3, n_k is the outward normal vector of the contour C, σ_ij is the stress tensor, u_i is the displacement vector, A is the area enclosed by the contour C. The separation of the modes is possible by using a direct method (Rigby and Aliabadi, 1998; Shivakumar and Raju, 1992):

J_{1} = J_{I} + J_{II} + J_{III}, J_{2} = - 2 \sqrt{J_{I} J_{II}}, J_{3} = J_{III}

(49)

In our problem, x₁ = x, x₂ = z and x₃ = y. For the calculation, we apply a zero-area path around the crack tip (Szekrényes, 2012a). This way the surface integral in equation (47) becomes zero. The layerwise stress–strain relations in laminated composite plates are (Kollár and Springer, 2003):

[σ_{x} σ_{x} τ_{xy}]^{(k)} = [{\bar{C}}_{11} {\bar{C}}_{12} 0 {\bar{C}}_{21} {\bar{C}}_{22} 000 {\bar{C}}_{66}]^{(k)} [ɛ_{x} ɛ_{y} γ_{xy}], [τ_{xz} τ_{yz}]^{(k)} = [{\bar{C}}_{55} 00 {\bar{C}}_{44}]^{(k)} [γ_{xz} γ_{yz}]

(50)

where {

{\bar{C}}_{ij}

} is the stiffness matrix of the k^th layer. Taking the stresses and strains back into J₁ in equation (47), we obtain an expression including the product of stress resultants and strain components in the form of:

J_{1} = \frac{1}{2} (N_{x 1} ɛ_{x 1}^{(0)} |_{\binom{x = -}{+ 0}} - N_{x 2} ɛ_{x 2}^{(0)} |_{x = - 0}) - \frac{1}{2} (N_{y 1} ɛ_{y 1}^{(0)} |_{x = + 0} - N_{y 2} ɛ_{y 2}^{(0)} |_{x = - 0}) + \frac{1}{2} (M_{x 1} ɛ_{x 1}^{(2)} |_{x = + 0} - M_{x 2} ɛ_{x 2}^{(1)} |_{x = - 0}) - \frac{1}{2} (M_{y 1} ɛ_{y 1}^{(1)} |_{x = + 0} - M_{y 2} ɛ_{y 2}^{(1)} |_{x = - 0}) + \frac{1}{2} (L_{x 1} ɛ_{x 1}^{(2)} |_{x = + 0} - L_{x 2} ɛ_{x 2}^{(2)} |_{x = - 0}) - \frac{1}{2} (L_{y 1} ɛ_{y 1}^{(2)} |_{x = + 0} - L_{y 2} ɛ_{y 2}^{(2)} |_{x = - 0}) - \frac{1}{2} (N_{xy 1} {\overset{\land}{γ}}_{xy 1}^{(0)} |_{x = + 0} - N_{xy 2} {\overset{\land}{γ}}_{xy 2}^{(0)} |_{x = - 0}) - \frac{1}{2} (M_{xy 1} {\overset{\land}{γ}}_{xy 1}^{(1)} |_{x = + 0} - M_{xy 2} {\overset{\land}{γ}}_{xy 2}^{(1)} |_{x = - 0}) - \frac{1}{2} (L_{xy 1} {\overset{\land}{γ}}_{xy 1}^{(2)} |_{x = + 0} - L_{xy 2} {\overset{\land}{γ}}_{xy 2}^{(2)} |_{x = - 0}) + (Q_{x 1} (\frac{1}{2} γ_{xz 1}^{(0)} - \frac{\partial w_{1}}{\partial x}) |_{x = + 0} - Q_{x 2} (\frac{1}{2} γ_{xz 2}^{(0)} - \frac{\partial w_{2}}{\partial x}) |_{x = - 0}) + \frac{1}{2} (Q_{y 1} γ_{yz 1}^{(0)} |_{x = + 0} - Q_{y 2} γ_{yz 2}^{(0)} |_{x = - 0}) + \frac{1}{2} (R_{x 1} γ_{xz 1}^{(1)} |_{x = + 0} - R_{x 2} γ_{xz 2}^{(1)} |_{x = - 0}) + \frac{1}{2} (R_{y 1} γ_{yz 1}^{(1)} |_{x = + 0} - R_{y 2} γ_{yz 2}^{(1)} |_{x = - 0})

(51)

where e.g. N_x1 is the in-plane normal force for the delaminated portion,

ɛ_{x 1}^{(0)}

is the constant part of the strain in the x direction as it is defined in equation (6), etc. Moreover, the shear strains including the hat are

{\overset{\land}{γ}}_{xy}^{(0)} = \frac{\partial u_{0}}{\partial y} - \frac{\partial v_{0}}{\partial x}, {\overset{\land}{γ}}_{xy}^{(1)} = \frac{\partial θ_{x}}{\partial y} - \frac{\partial θ_{y}}{\partial x}, {\overset{\land}{γ}}_{xy}^{(2)} = \frac{\partial φ_{x}}{\partial y} - \frac{\partial φ_{y}}{\partial x}

(52)

Considering the fact that M_x, M_y, M_xy, Q_x, Q_y, R_x, R_y and the corresponding strains are continuous across portions “1” and “2,” a significant part of the J-integral vanishes. The remaining part can be separated based on the direct method (refer to equation (49)) or by simply separating the terms including the sin (mode-II) and cos (mode-III) functions, leading to:

J_{II} = \frac{1}{2} (N_{x 1} ɛ_{x 1}^{(0)} |_{x = + 0} - N_{x 2} ɛ_{x 2}^{(0)} |_{x = - 0}) - \frac{1}{2} (N_{y 1} ɛ_{y 1}^{(0)} |_{x = + 0} - N_{y 2} ɛ_{y 2}^{(0)} |_{x = - 0}) + \frac{1}{2} (L_{x 1} ɛ_{x 1}^{(2)} |_{x = + 0} - L_{x 2} ɛ_{x 2}^{(2)} |_{x = - 0}) - \frac{1}{2} (L_{y 1} ɛ_{y 1}^{(2)} |_{x = + 0} - L_{y 2} ɛ_{y 2}^{(2)} |_{x = - 0})

(53)

J_{III} = \frac{1}{2} (N_{xy 1} {\overset{\land}{γ}}_{xy 1}^{(0)} |_{x = + 0} - N_{xy 2} {\overset{\land}{γ}}_{xy 2}^{(0)} |_{x = - 0}) - \frac{1}{2} (L_{xy 1} {\overset{\land}{γ}}_{xy 1}^{(2)} |_{x = + 0} - L_{xy 2} {\overset{\land}{γ}}_{xy 2}^{(2)} |_{x = - 0})

(54)

Consequently, under mixed-mode II/III fracture condition, the ERR is not contributed by the bending and twisting moments, shear forces, as well as L_xy, R_x, R_y (higher order stress resultants) at all. Also, the present formulation does not incorporate any physically inconsistent parameters (e.g. shear compliances like in Qiao and Wang, 2004; Szekrényes, 2012a; Wang and Qiao, 2004b), it is based on an entirely exact formulation including the material law of orthotropic solids.

Results and discussion

The properties of the analyzed simply-supported plate were (refer to Figure 2): a = 105 mm (crack length), c = 45 mm (uncracked length), b = 100 mm (plate width), t = 2 mm (half plate thickness), Q₀ = 1000 N (point force magnitude), x_Q = 31 mm, y_Q = 50 mm (point of action coordinates of Q₀, d₀ = 0.1 mm and q₀ = Q₀/2d₀. The plate is made of a carbon/epoxy material, the lay-up of the plate was [

\pm 45_{2}^{f}

; 0₁₂;

\pm 45_{2}^{f}

] for the delaminated and [

\pm 45_{2}^{f}

; 0₁₂;

\pm 45_{2}^{f}

]_s for the uncracked part. The superscript “f” refers to the fact that the cross-ply laminate is a woven fabric panel. The properties of the individual laminae are given by Table 1 (Kollár and Springer, 2003). The computation was performed in the code MAPLE 12 (Garvan, 2002) in accordance with the following points. The stiffness matrices of each single layer, the stiffness matrices ({A}, {D} and {F}) of the plate, were determined based on the elastic properties of the laminates given in Table 1. The problem in Figure 3 was solved varying the number of Fourier series terms (N) by creating a for-do cycle. Based on the displacement parameters, the stress resultants and the stresses were calculated, while the ERRs were calculated using the J-integral. The convergence of the results was analyzed and it was found that after the 13^th Fourier term there was no change in the displacement field, stresses, stress resultants and ERRs.

Table 1.

Elastic properties of single carbon/epoxy laminates.

E_x (Gpa)	E_y (Gpa)	E_z (Gpa)	G_yz (Gpa)	G_xz (Gpa)	G_xy (Gpa)	ν _yz	ν _xz	ν _xy
±45°F	16.39	16.39	16.4	5.46	16.4	5.46	0.5	0.3	0.5
0°	148	9.65	9.65	4.91	3.71	4.66	0.27	0.3	0.25

Finite element model

In order to verify the analytical results, a FE analysis was carried out. The 3D FE model of the plate was created in the code ANSYS 12 using 8-node solid elements. Similar 3D models are documented by de Morais and Pereira (2008) and Pereira and de Morais (2009), therefore the model is not shown here: 50, 78 and 10 elements were applied along the plate width (y), length (x) and thickness (z), respectively. The global element size was 2 mm × 2 mm × 0.4 mm. In the vicinity of the crack tip, a refined mesh was constructed including trapezoid shape elements (Davidson and Sundararaman, 1996). The displacements in the z direction of the contact nodes over the delaminated surface were imposed to be the same. The mode-II and mode-III ERRs were calculated by the VCCT (e.g., de Morais and Pereira, 2008), the size of the crack tip elements were Δx = 0.2 mm, Δy = 0.2 mm and Δz = 2 mm. For the determination of G_II and G_III along the delamination front, a so-called 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 node along the crack front.

Analytical and numerical results

Figure 4 shows the in-plane displacements at two points of the plate: u(0,b/2,z) and v(0,b,z), i.e. each point lies in the crack front. It is seen that the analytical solution agrees excellently with the FE results. Although the nonlinearity of in-plane displacements is not so significant, the SSDT captures this change in the through-thickness direction very well.

The distribution of the normal stresses, σ_x and σ_y at x = 0, y = b/2, are demonstrated in Figure 5. The distributions were determined layerwise using the stress–strain relations given by equation (50). Although there are some differences compared to the FE results, the overall agreement is reasonable, especially in the case of σ_y. It has to be mentioned that in the FE model the displacement and stress continuity is ensured, but in the analytical model only the continuity of displacement parameters and some of the stress resultants (N_x, M_x, M_y, M_xy, Q_x, Q_y) can be realized. Consequently the stresses are not continuous in accordance with the SSDT, however, the discrepancies in stresses in the transition between the delaminated and undelaminated plate portions are not so significant.

The distributions of transverse shear stresses are plotted in Figure 6. The FE solution, the solution by SSDT (piecewise linear, red line) as well as the solution calculated by the 3D equilibrium equations are equally presented. The latter ones were calculated by the $\underline{\underline{σ}} \nabla = \underline{0}$ equation (Reddy, 2004). The mismatch of the FE and analytical solutions are somewhat higher than in the case of the normal stresses, however the overall agreement is acceptable. From another point of view, the SSDT and 3D equilibrium solutions agree well. The distribution of the in-plane normal force, N_x(x,y) and the shear force N_xy(x,y) are plotted in Figure 7.

Although these stress resultants could be calculated by the FE model too (by integrating the normal and in-plane shear stresses in the through thickness direction), this would be a very long process, therefore, in this case only analytical results are presented. It is seen that near the delamination front (x = 0) these stress resultants change suddenly and significantly. Moreover, N_x involves the sin, while N_xy involves the cos function (refer to Equation (43)) leading to the completely different nature of these stress resultants. One of the advantages of plate theory over the FE model is that these plots can be obtained very simply in the analytical way.

The interlaminar shear stress (=transverse shear in accordance with the duality of shear stresses) distributions along the global midplane of the undelaminated plate portion (top plate, refer to Figure 1) can be calculated as:

τ_{xz} |_{z = - \frac{t}{2}} = G_{13}^{\pm 45} γ_{xz} |_{z = - \frac{t}{2}} = G_{13}^{\pm 45} (\frac{4}{t} u_{02} - θ_{x 2} + \frac{\partial w_{2}}{\partial x}), τ_{yz} |_{z = - \frac{t}{2}} = G_{23}^{\pm 45} γ_{yz} |_{z = - \frac{t}{2}} = G_{23}^{\pm 45} (\frac{4}{t} v_{02} - θ_{y 2} + \frac{\partial w_{2}}{\partial y})

(55)

where the shear moduli are defined in Table 1. The function plots are given by Figure 8. Both stresses change significantly near the delamination front and decays suddenly by getting far from it. It is important to note that τ_xz does not vanish at the x = −c boundary, on the other hand τ_yz vanishes entirely at the same location.

The mode-II, mode-III ERRs and the mode ratios along the delamination front are plotted in Figure 9. The symbols show the results of the VCCT, while the curves represent the results by SSDT. The main conclusion is that except for the relatively small regions at the edges of the plate, the SSDT agrees excellently with the results by VCCT. Also, it is clear that near the edges the agreement becomes not so good. In Figure 9, the dashed blue line shows that after this line towards the plate edge, the difference between the SDDT and the VCCT in the case of G_III becomes higher than 10%. The corresponding distance from the edge of the plate is again approximately 10% of the plate width. In Figure 9(b) G_T = G_II + G_III is the total ERR. In the case of the ratios of G_II/G_T and G_III/G_T, the agreement is excellent, however at the edges some discrepancies appear again. These differences are attributed to the distinctions in the boundary conditions. The plate theory assumes that the midplane of both the top and bottom plates is simply-supported. On the other hand, in the FE model – relating to practical conditions – only the contour lines of the bottom plates are supported in the z direction. In spite of the distinctions in the stress distributions by FEM and analyses in Figures 5 and 6, it is seen that the energy release rates agree well. From equations (53)–(54) it is clear that the ERRs depend on the product of stress resultants and strain components. The stress resultants are calculated as the integration of the stress distributions over the thickness, and in Figures 5 and 6 the area under the curves are approximately the same for the FE and plate theory solutions. That is why the two solutions match well in the case of the ERRs.

The overall agreement between the SSDT and VCCT is fairly good. It is important to mention that the VCCT method is mesh-sensitive to a certain degree and the investigation of the effect of mesh refinement was outside the scope of this paper. It has to be also mentioned that the present model does not take the effect of possible nonlinearites into account, like fiber-bridgings (e.g., Tamuzs et al., 2001) in the delaminated area and the effect of the so-called fracture process zone (FPZ; e.g., Amrutharaja et al., 1995; Tsouvalis and Anyfantis, 2012).

Conclusions

The second order shear deformation plate theory is utilized in this work to develop a double-plate system for delaminated orthotropic composite plates. The model is based on the continuity of the displacement field across the delamination front by imposing the interface constraint along the interface. A simply-supported delaminated plate subjected to a point force was analyzed using Lévy plate formulation, the stresses and the energy release rates were calculated. The results were compared to those of a 3D FE model and very good agreement was obtained.

The present model eliminates the physically inconsistent shear compliance of the flexible joint models and includes the effect of interface deformation based on the equations of linear elasticity and the material law of orthotropic solids. It was shown that although the displacement components are continuous across the delamination, there are stress resultants, which remain discontinuous. Moreover, a significant amount of terms vanish in the J-integral, and the mode-II and mode-III energy release rates are defined using the stress resultants and strains around the delamination front. However, it must be mentioned that only mixed-mode II/III fracture problems were considered in this study and the transverse deflection of the top and bottom plates was considered to be the same. It has been shown that the difference between the FE and plate theory solutions differs moderately at the edges of the plate, which can be explained by the differences in the boundary conditions.

Considering the available methods for the calculation of the ERR in plates, the first alternative is in general the VCCT. However, for the 3D FE model the computation could be lengthy, especially if the model has relatively large dimensions. Furthermore, in the crack tip a refined mesh should be constructed to obtain accurate ERR values. Finally, in most of the commercial FE packages the VCCT has not yet been implemented. The present work provides another possibility for the calculation of the ERR in plates subjected to bending. The possible application field of the presented method is the fracture mechanics of composite materials. In the last few years, fracture test methods including plate-like specimens have been developed to characterize the mode-III, mixed-mode II/III and mixed-mode I/III fracture behavior of laminated materials. By preparing a detailed user-friendly worksheet in MAPLE, it is possible to provide a data reduction scheme for delaminated plates for the experimentalists. Also, the application to asymmetrically delaminated orthotropic and angle-ply laminated plates as well as sandwich panels need to be investigated. These tasks will be carried out in the near future.

Footnotes

Acknowledgements

This work was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences. This work is connected to the scientific program of the “Development of quality-oriented and harmonized R + D + I strategy and functional model at BME” project. This project is supported by the New Hungary Development Plan (Project ID: TÁMOP-4.2.1/B-09/1/KMR-2010-0002).

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

References

Amrutharaja

Lama

Cotterell

(1995) Fracture process zone concept and delamination of composite laminates. Theoretical and Applied Fracture Mechanics 24(1): 57–64.

Anderson

(2005) Fracture Mechanics – Fundamentals and Applications, 3rd edn. Boca Raton, London, New York, and Singapore: CRC Press, Taylor & Francis Group.

Argüelles

Viña

Canteli

(2011) Influence of resin type on the delamination behavior of carbon fiber reinforced composites under mode-II loading. International Journal of Damage Mechanics 20(7): 963–978.

Browning

Carlsson

Ratcliffe

(2010) Redesign of the ECT test for mode III delamination testing. Part I: finite element analysis. Journal of Composite Materials 44: 1867–1881.

Browning

Carlsson

Ratcliffe

(2011) Modification of the edge crack torsion specimen for mode III delamination testing. Part II-experimental study. Journal of Composite Materials 45: 2633–2640.

Chen F (2011) Interface mechanics of layered composite beam-type structures. PhD Thesis, Department of Civil and Environmental Engineering, Washington State University.

Cherepanov

(1997) Methods of Fracture Mechanics: Solid Matter Physics, Dordrecht, Boston, MA, and London: Kluwer Academic Publishers.

Chou

Pagano

(1967) Elasticity – Tensor, Dyadic, and Engineering Approaches, Princeton, NJ, Toronto, and London: D. Van Nostrand Company Inc.

Davidson

Sundararaman

(1996) A single leg bending test for interfacial fracture toughness determination. International Journal of Fracture 78: 193–210.

10.

Davidson

Sediles

(2011) Mixed-mode I-II-III delamination toughness determination via a shear-torsion-bending test. Composites Part A: Applied Science and Manufacturing 42: 589–603.

11.

Davidson BD, Sediles FO and Humphreys KD (2010) A shear-torsion-bending test for mixed-mode I-II-III delamination toughness determination. In: 25th technical conference of the American society for composites and 14th US-Japan conference on composite materials, Dayton, Ohio, USA, 2010, pp.1001–1020.

12.

Davidson

(2000) Determination of energy release rate and mode mix in three-dimensional layered structures using plate theory. International Journal of Fracture 105: 81–104.

13.

de Morais

Pereira

(2008) Mixed mode II + III interlaminar fracture of carbon/epoxy laminates. Composites Science and Technology 68: 2022–2027.

14.

de Morais

Pereira

(2009) Mode III interlaminar fracture of carbon/epoxy laminates using a four-point bending plate test. Composites Part A - Applied Science and Manufacturing 40: 1741–1746.

15.

de Moura

MFSF

Fernandez

MVC

de Morais

(2009) Numerical analysis of the edge crack torsion test for mode III interlaminar fracture of composite laminates. Engineering Fracture Mechanics 76: 469–478.

16.

Garvan

(2002) The Maple Book, Boca Raton, London, New York, and Washington, DC: Chapman & Hall/CRC.

17.

Gornet

Ijaz

(2011) High cycle fatigue damage model for delamination crack growth in CF/epoxy composite laminates. International Journal of Damage Mechanics 20(5): 783–807.

18.

Jumel

Budzik

Shanahan

MER

(2011) Beam on elastic foundation with anticlastic curvature: application to analysis of mode I fracture tests. Engineering Fracture Mechanics 78: 3253–3269.

19.

Khdeir

Reddy

(1999) Free vibrations of laminated composite plates using second-order shear deformation theory. Computers & Structures 71: 617–626.

20.

Kollár

Springer

(2003) Mechanics of Composite Structures, Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, and São Paolo: Cambridge University Press.

21.

Kondo

Sato

Suemasu

(2010) Characterization of fracture resistance of carbon/epoxy composite laminates during mixed-mode II and III stable damage propagation. Journal of the Japan Society for Composite Materials 36(5): 179–188.

22.

Kondo

Sato

Suemasu

(2011) Fracture resistance of carbon/epoxy composite laminates under mixed-mode II and III failure and its dependence on fracture morphology. Advanced Composite Materials 20(5): 405–418.

23.

Lee

(1993) Fall. An edge crack torsion method for mode III delamination fracture testing. Journal of Composite Technology and Research 15: 193–201.

24.

O’Brien

(1996) Simplified data reduction methods for the ECT test for mode III interlaminar fracture toughness. Journal of Composite Technology and Research 18: 96–101.

25.

Lundmark

Varna

(2011) Stiffness reduction in laminates at high intralaminar crack density: effect of crack interaction. International Journal of Damage Mechanics 20(2): 279–297.

26.

Mehrabadi

(2012) Analysis of pure mode III and mixed mode (III + II) interlaminar crack growth in polymeric woven fabrics. Materials & Design 44: 429–437.

27.

Mehrabadi FA and Khosravan M (2013) Mode III interlaminar fracture and damage characterization in woven fabric-reinforced glass/epoxy composite laminates. Journal of Composite Materials 47(13): 1583--1592.

28.

Miura

Shindo

Takeda

(2012) Interlaminar fracture characterization of woven glass/epoxy composites under mixed-mode II/III loading conditions at cryogenic temperatures. Engineering Fracture Mechanics 96: 615–625.

29.

Nagarajan

Sundaram

Thyagarajan

(2012) Refined delamination factor failure characterization of composite wind turbine blade. International Journal of Damage Mechanics 21: 1227–1244.

30.

Omiya

Kishimoto

(2010) Damage-based cohesive zone model for rate-depend interfacial fracture. International Journal of Damage Mechanics 19(4): 397–420.

31.

Oskay

Pal

(2010) A multiscale failure model for analysis of thin heterogeneous plates. International Journal of Damage Mechanics 19(5): 575–610.

32.

Park

Sankar

(2002) Crack-tip force method for computing energy release rate in delaminated plates. Composite Structures 55: 429–434.

33.

Peng

Zhang

Zhao

(2011) Mode I delamination growth of multidirectional composite laminates under fatigue loading. Journal of Composite Materials 45: 1077–1090.

34.

Pereira

de Morais

(2009) Mixed-mode I + III interlaminar fracture of carbon/epoxy laminates. Composites Part A: Applied Science and Manufacturing 40: 518–523.

35.

Plain

Tong

(2011) An experimental study on mode I and II fracture toughness of laminates stitched with a one-sided stitching technique. Composites Part A: Applied Science and Manufacturing 42: 203–210.

36.

Qiao

Chen

(2011) On the compliance and energy release rate of generically-unified beam-type fracture specimens. Journal of Composite Materials 45(1): 65–101.

37.

Qiao

Chen

(2012) On the improved dynamic analysis of delaminated beams. Journal of Sound and Vibration 331: 1143–1163.

38.

Qiao

Wang

(2004) Mechanics and fracture of crack-tip deformable bimaterial interface. International Journal of Solids and Structures 41: 7423–7444.

39.

Reddy

(2004) Mechanics of laminated composite plates and shells – theory and analysis, Boca Raton, London, New York, and Washington, DC: CRC Press.

40.

Reeder

Crews

(1990) Mixed-mode bending method for delamination testing. AIAA Journal 28: 1270–1276.

41.

Rice

(1968) A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics 35: 379–386.

42.

Rigby

Aliabadi

(1998) Decomposition of the mixed-mode J-integral – revisited. International Journal of Solids and Structures 35: 2073–2099.

43.

Sankar

Sonik

(1995) Pointwise energy release rate in delaminated plates. AIAA Journal 33: 1312–1318.

44.

Shahrjerdi A and Mustapha F (2011) Second order shear deformation theory (SSDT) for free vibration analysis on a functionally graded quadrangle plate. In: Baddour N (ed) Recent advances in vibration analysis. Croatia, Europe; Shanghai, China: Intech, pp. 59–78.

45.

Shivakumar

Raju

(1992) An equivalent domain integral method for three-dimensional mixed-mode fracture problems. Engineering Fracture Mechanics 42: 935–959.

46.

Sriram

Armanios

(1993) A shear deformation model for transverse cracking in composite laminates. International Journal of Damage Mechanics 2: 73–91.

47.

Suemasu H, Kondo A, Gozu K, et al. (2009) Double notched split cantilever test method to measure the mixed mode II and III interlaminar toughness. In: ICCM-17 17th international conference on composite materials. Paper ID. F9.6. 27–31 July 2009. Edinburgh, UK: Edinburgh International Convention Centre (EICC).

48.

Suemasu

Kondo

Gozu

(2010) Novel test method for mixed mode II and III interlaminar fracture toughness. Advanced Composite Materials 19(4): 349–361.

49.

Szekrényes

(2007) Delamination fracture analysis in the G_II-G_III plane using prestressed composite beams. International Journal of Solids and Structures 44: 3359–3378.

50.

Szekrényes

(2009a) Improved analysis of the modified split-cantilever beam for mode III fracture. International Journal of Mechanical Sciences 51: 682–693.

51.

Szekrényes

(2009b) Interlaminar fracture analysis in the G_I-G_III plane using prestressed composite beams. Composites Part A – Applied Science and Manufacturing 40: 1621–1631.

52.

Szekrényes

(2011) Interlaminar fracture analysis in the G_I-G_II-G_III space using prestressed composite beams. Journal of Reinforced Plastics and Composites 30: 1655–1669.

53.

Szekrényes

(2012a) Interlaminar fracture analysis in the G_II-G_III plane using prestressed composite beams. Composites Part A – Applied Science and Manufacturing 43: 95–103.

54.

Szekrényes

(2012b) Interlaminar stresses and energy release rates in delaminated orthotropic composite plates. International Journal of Solids and Structures 49: 2460–2470.

55.

Szekrényes A (2013) Interface crack between isotropic Kirchhoff plates. Meccanica. DOI: 10.1007/s11012-012-9611-9.

56.

Tamuzs

Tarasovs

Vilks

(2001) Progressive delamination and fiber bridging modeling in double cantilever beam composite specimens. Engineering Fracture Mechanics 68(5): 513–525.

57.

Tian

(2010) Modeling of composite laminated plates with interfacial damage evolution. International Journal of Damage Mechanics 20(3): 370–399.

58.

Tsouvalis

Anyfantis

(2012) Determination of the fracture process zone under mode I fracture in glass fiber composites. Journal of Composite Materials 46(1): 27–41.

59.

Wang

Qiao

(2004a) Novel beam analysis of the end notched flexure specimen for mode-II fracture. Engineering Fracture Mechanics 71: 219–231.

60.

Wang

Qiao

(2004b) Interface crack between two shear deformable elastic layers. Journal of the Mechanics and Physics of Solids 52: 891–905.

61.

Yazdi AA and Rezaeepazhand J (2012) Applicability of small-scale models in prediction flutter pressure of delaminated composite beam-plates. International Journal of Damage Mechanics (to appear).