**2.1 Constitutive relations**

The models developed by Talreja & Boehler (1990), Allen et. al. (1987 a,b) and Lim and Tay (1996) to describe the damaged composite laminates were based on the Continuum Damage Mechanics using internal state variables. In the presente paper, the model proposed by Allen et al (1987 a,b) will be used, which describes the damage through a set of internal state variables. The final result of the distributed damage is built in the constitutive equations through these variables. Thus, the stress-strain relationship of the representative volume of a damaged material at the level of a lamina is assumed as:

$$
\sigma\_{ij} = \mathbb{C}\_{ijkl}\varepsilon\_{kl} + I\_{ijkl}^{\eta} \alpha\_{kl}^{\eta} \tag{1}
$$

Progressive Stiffness Loss Analysis

internal state variable of the problem.

**2.2 Determination of internal variables** 

of Symmetric Laminated Plates due to Transverse Cracks Using the MLGFM 127

*k k k k*

*k k k k*

*x y xy <sup>t</sup>*

*x y xy <sup>t</sup>*

1

1

22 12

 

*<sup>k</sup> <sup>k</sup>* 

where *P* is the total number of damage models being considered and 22 and 12 are the

In spite of the random characteristic, as can be found in the work of Silberschmidt (2005), the transverse cracks are assumed to be uniformly distributed. In this way, the laminate behavior can be adequately represented by a representative unit volume of material containing a transverse crack, as shown in figure 2. In the particular case of symmetric laminates, the damage models are simplified and incorporate only two types of fracture, namely, Mode I (crack opening) and Mode I coupled with Mode III (shear out of plane). They are represented respectively by the internal variables 22 e 12. As only symmetric

The internal variable, equation (15), proposed by Allen (Allen et al., 1987 a,b) can be determined by a computational analysis based on Finite Element Method. The representative volume is modeled, as shown in figure 3, for the symmetric laminate [0o/90o/0o]. A uniform displacement is imposed in one side of the element to determine the opening of the crack. The size and shape of the representative volume depend on the thickness of the different plies and the crack density (the number of cracks per unit of

> 22 2 2 1

*V*

*Sc*

*u n dS*

(15)

 

1

1

(8)

(9)

*dz* (10)

*dz* (11)

(12)

(13)

(14)

\_ 2 2

\_ 3 3

 /2 <sup>0</sup> /2 *t*

 /2 <sup>0</sup> /2 *t*

 

\_ \_

\_ \_ 2 2

*k k k k k*

*D zz Q*

*k k k k k*

1

1

*D zz Q*

1 0 *P*

1 2 *N*

The internal state variables vector has two components and is expressed by:

laminates are analyzed in this paper, just the *α*22 variable will be developed.

volume). Then, the internal variable can be determined by:

*N*

*D zz Q* 

*B zz Q* 

1

1

1 2 *N*

1 3 *N*

*N*

*M*

where ij is the applied stress tensor, Cijkl is the constitutive relation tensor of the undamaged material, *kl* are the strain tensor, *ijkl I* are the elements of the damage matrix, *kl* are the internal state variables, and = 1, 2, 3, …, refers to the damage modes. As suggested by Allen *et al* (1987 a,b), a first simplification can be made considering that the tensor Iijkl is the actual tensor of constitutive relationships, as shown in Equation (2).

$$I\_{ijkl}^{\eta} = \mathbb{C}\_{ijkl} \tag{2}$$

However, it is important to emphasize that equation (1) does not provide any information on how the damage state has been attained, that is, the history of damage accumulation. Thus, it is necessary to turn to Fracture Mechanics in search of a suitable criterion to evaluate the damage growth. Thereby, equation (1) is sufficiently general to permit the use of Classical Laminate Theory to determine the composite laminate constitutive relations with transverse cracks in the matrix.

Supposing the representation of the laminated plate by plane elements located in its middle surface, the loads in a certain point inside this surface can be evaluated by the following expressions:

$$\{N\} = \int\_{-t/2}^{t/2} \begin{pmatrix} \sigma\_x & \sigma\_y & \tau\_{xy} \end{pmatrix} dz \tag{3}$$

$$\mathbb{E}\left\{M\right\} = \int\_{-t/2}^{t/2} \begin{pmatrix} \sigma\_x & \sigma\_y & \tau\_{xy} \end{pmatrix} zdz\tag{4}$$

where {*N*}e {*M*} are, respectively, the force and moment resultants vectors, *x, y, xy* are the stresses in the plane of the lamina and *t* is the thickness of the laminate. Taking to account that {*ε0*} e {*κ0*} are the strain and bending vectors in the middle surface of the plate, [*A*], [*B*] and [*D*] are the laminate extensional stiffness matrix, coupling stiffness matrix and bending stiffness matrix, respectively, {*DN*}e {*DM*} are the damage vectors related to the force and moment resultants, the expressions (3) and (4) can be transformed to:

$$\mathbb{E}\{N\} = \left[A\right]\left\{\boldsymbol{\varepsilon}\_{0}\right\} + \left[B\right]\left\{\boldsymbol{\kappa}\_{0}\right\} + \left\{\boldsymbol{D}^{N}\right\} \tag{5}$$

$$\{M\} = \left[B\right]\{\varepsilon\_0\} + \left[D\right]\{\kappa\_0\} + \left\{D^M\right\} \tag{6}$$

Assuming that zk-1 e zk are the corresponding distances from the middle surface to the inner and outer surfaces of the *k*th lamina, respectively, [ ] *Q <sup>k</sup>* and { } *<sup>k</sup>* are the transformed reduced material stiffness matrix and the transformed vector of the internal state variables (expressed in global coordinates), respectively. The matrix and vectors presents in equations (5) and (6) can be expressed by:

$$\mathbb{E}\left[A\right] = \sum\_{k=1}^{N} \left(z\_k - z\_{k-1}\right) \left[\bar{Q}\right]\_k \tag{7}$$

*ij ijkl kl C Iijkl kl* 

where ij is the applied stress tensor, Cijkl is the constitutive relation tensor of the

suggested by Allen *et al* (1987 a,b), a first simplification can be made considering that the

*ijkl ijkl I C* 

However, it is important to emphasize that equation (1) does not provide any information on how the damage state has been attained, that is, the history of damage accumulation. Thus, it is necessary to turn to Fracture Mechanics in search of a suitable criterion to evaluate the damage growth. Thereby, equation (1) is sufficiently general to permit the use of Classical Laminate Theory to determine the composite laminate constitutive relations

Supposing the representation of the laminated plate by plane elements located in its middle surface, the loads in a certain point inside this surface can be evaluated by the following

> /2 /2 *t*

 /2 /2 *t*

*x y xy <sup>t</sup> <sup>M</sup>* 

where {*N*}e {*M*} are, respectively, the force and moment resultants vectors,

moment resultants, the expressions (3) and (4) can be transformed to:

and outer surfaces of the *k*th lamina, respectively, [ ] *Q <sup>k</sup>* and { }

*x y xy <sup>t</sup> N d* 

stresses in the plane of the lamina and *t* is the thickness of the laminate. Taking to account that {*ε0*} e {*κ0*} are the strain and bending vectors in the middle surface of the plate, [*A*], [*B*] and [*D*] are the laminate extensional stiffness matrix, coupling stiffness matrix and bending stiffness matrix, respectively, {*DN*}e {*DM*} are the damage vectors related to the force and

> 0 0 *<sup>N</sup> NA B D*

 0 0 *MB D D<sup>M</sup>* 

Assuming that zk-1 e zk are the corresponding distances from the middle surface to the inner

reduced material stiffness matrix and the transformed vector of the internal state variables (expressed in global coordinates), respectively. The matrix and vectors presents in equations

\_

*k k k k A zz Q* 

1

*N*

 

 

1

  (1)

are the elements of the damage matrix,

= 1, 2, 3, …, refers to the damage modes. As

(2)

*<sup>z</sup>* (3)

*zdz* (4)

(7)

*x, y,* 

(5)

(6)

*<sup>k</sup>* are the transformed

*xy* are the

 

tensor Iijkl is the actual tensor of constitutive relationships, as shown in Equation (2).

are the strain tensor, *ijkl I*

undamaged material, *kl*

*kl* 

expressions:

are the internal state variables, and

with transverse cracks in the matrix.

(5) and (6) can be expressed by:

$$\mathbb{E}\left[B\right] = \frac{1}{2} \sum\_{k=1}^{N} \left(z\_k^2 - z\_{k-1}^2\right) \left[\bar{Q}\right]\_k \tag{8}$$

$$\mathbb{E}\left[D\right] = \frac{1}{\mathfrak{Z}} \sum\_{k=1}^{N} \left(z\_k^3 - z\_{k-1}^3\right) \left[\bar{\mathcal{Q}}\right]\_k \tag{9}$$

$$\mathbb{E}\left\{\boldsymbol{\varepsilon}^{0}\right\} = \int\_{-t/2}^{t/2} \begin{Bmatrix} \boldsymbol{\varepsilon}\_{x} & \boldsymbol{\varepsilon}\_{y} & \boldsymbol{\gamma}\_{xy} \end{Bmatrix} d\boldsymbol{z} \tag{10}$$

$$\mathbb{E}\left\{\boldsymbol{\kappa}^{0}\right\} = \int\_{-t/2}^{t/2} \begin{pmatrix} \kappa\_x & \kappa\_y & \kappa\_{xy} \end{pmatrix} d\boldsymbol{z} \tag{11}$$

$$\mathbb{E}\left\{\boldsymbol{D}^{N}\right\}=\sum\_{k=1}^{N}\left(\boldsymbol{z}\_{k}-\boldsymbol{z}\_{k-1}\right)\left[\bar{\boldsymbol{Q}}\right]\_{k}\left\{\bar{\boldsymbol{a}}\right\}\_{k}\tag{12}$$

$$\mathbb{E}\left\{\boldsymbol{D}^{\mathcal{M}}\right\} = \frac{1}{2} \sum\_{k=1}^{N} \left(\boldsymbol{z}\_{k}^{2} - \boldsymbol{z}\_{k-1}^{2}\right) \left[\bar{\boldsymbol{Q}}\right]\_{k} \left\{\bar{\boldsymbol{a}}\right\}\_{k} \tag{13}$$

The internal state variables vector has two components and is expressed by:

$$\left\{ \alpha \right\}\_{k} = \sum\_{\eta=1}^{p} \begin{Bmatrix} 0 & \alpha\_{22}^{\eta} & \alpha\_{12}^{\eta} \end{Bmatrix}\_{k} \tag{14}$$

where *P* is the total number of damage models being considered and 22 and 12 are the internal state variable of the problem.

#### **2.2 Determination of internal variables**

In spite of the random characteristic, as can be found in the work of Silberschmidt (2005), the transverse cracks are assumed to be uniformly distributed. In this way, the laminate behavior can be adequately represented by a representative unit volume of material containing a transverse crack, as shown in figure 2. In the particular case of symmetric laminates, the damage models are simplified and incorporate only two types of fracture, namely, Mode I (crack opening) and Mode I coupled with Mode III (shear out of plane). They are represented respectively by the internal variables 22 e 12. As only symmetric laminates are analyzed in this paper, just the *α*22 variable will be developed.

The internal variable, equation (15), proposed by Allen (Allen et al., 1987 a,b) can be determined by a computational analysis based on Finite Element Method. The representative volume is modeled, as shown in figure 3, for the symmetric laminate [0o/90o/0o]. A uniform displacement is imposed in one side of the element to determine the opening of the crack. The size and shape of the representative volume depend on the thickness of the different plies and the crack density (the number of cracks per unit of volume). Then, the internal variable can be determined by:

$$
\alpha\_{22} = \frac{1}{V} \int\_{S\_c} u\_2 n\_2 dS \tag{15}
$$

Progressive Stiffness Loss Analysis

of Symmetric Laminated Plates due to Transverse Cracks Using the MLGFM 129

non-dimensionalized crack density (the quantity of cracks per unit of length), *δ* is the nondimensionalized maximum crack opening displacement, *u2* is the maximum crack opening displacement, *ψ* is a normalized function of the crack opening profile and is the normalized distance between the cracks center, the following expressions can be defined:

> 2*t L*

2 2 *u t* 

> 2 *u*( ) *u*

The maximum crack opening displacement *u2* can be determined by a simple finite element analysis, considering the boundary conditions specified in figure 3. An arbitrary displacement is imposed. The value of *δ* is determined by the equation (18), and the displacement *u2* is determined by MEF. The non-dimensionalized maximum crack opening

> 1 1 123 *a b ce ce c*

The constants *a1, b1, c1, c2* e *c3* in the expression (20) depends on the type of the used material. The table 1 exposes the value of these constants for the laminate glass/epoxi (Gl/Ep).

Glass/Epóxi 1.03 -0.81 2.28E-2 0.94 1.00

As the internal variable used in this problem depends on the maximum crack opening

Mode I: 22 2

In this paper, the maximum crack opening displacement *u2* is determined by the expressions (18) and (20). The crack density *ρ* depends on the distance between two adjacent transverse cracks *l*, and its values are successively modified by the verification of the composite

*u*

displacement according to equation (15), and the crack density is calculated by

can be shown that the state variable associated to the Mode I becomes:

 

**Formulation in terms of** *ρ* c1 c2 c3 a1 b1

> 8 5

(20)

(21)

1 /*L* , it

displacement *δ* can be obtained using *ρ*, as shown in the expression (20):

Table 1. Coefficients for the expression (20) (Lim & Tay, 1996)

**3. Approximation by computational methods**

material rigidity loss.

**Material** 

1 2 *tt t* (16)

(17)

(18)

(19)

where *u2* is the crack opening displacement, *n2* is the unitary vector normal to the crack surface, *V* is the representative element volume and *Sc* is the crack surface.

Fig. 2. A [0°/90°/0o] laminated plate with generalized cracks: definition of parameters and the representative volume (Machado et al, 2008).

Fig. 3. Boundary conditions and finite element mesh to evaluate the crack opening in a representative volume [0°/90°/0o] – (as suggested by Lim & Tay, 1996)

Considering that *t1* and *t2* are the thickness of the 0° and 90° plies, respectively, *t* is the total thickness of the laminate, *l* is the distance between two adjacent transverse cracks, *ρ* is the

where *u2* is the crack opening displacement, *n2* is the unitary vector normal to the crack

*t1*

*t2*

*L/2 <sup>X</sup>*

*u2*

*0o Ply*

*Y* 

*Representative Volume* 

Fig. 2. A [0°/90°/0o] laminated plate with generalized cracks: definition of parameters and

**t1**

Opening of the crack

Fig. 3. Boundary conditions and finite element mesh to evaluate the crack opening in a

Considering that *t1* and *t2* are the thickness of the 0° and 90° plies, respectively, *t* is the total thickness of the laminate, *l* is the distance between two adjacent transverse cracks, *ρ* is the

representative volume [0°/90°/0o] – (as suggested by Lim & Tay, 1996)

**L/2 u2**

the representative volume (Machado et al, 2008).

*L* 

**90o t2**

 **0o** 

*0o Ply*

 *90o Ply*

Free edge of the crack

surface, *V* is the representative element volume and *Sc* is the crack surface.

 *0o/90o/0o Laminated Plate* 

non-dimensionalized crack density (the quantity of cracks per unit of length), *δ* is the nondimensionalized maximum crack opening displacement, *u2* is the maximum crack opening displacement, *ψ* is a normalized function of the crack opening profile and is the normalized distance between the cracks center, the following expressions can be defined:

$$t = t\_1 + t\_2 \tag{16}$$

$$
\rho = \frac{2t}{L} \tag{17}
$$

$$
\delta = \frac{\mu\_2}{t\_2} \tag{18}
$$

$$\psi = \frac{\mu(\xi)}{\mu\_2} \tag{19}$$

The maximum crack opening displacement *u2* can be determined by a simple finite element analysis, considering the boundary conditions specified in figure 3. An arbitrary displacement is imposed. The value of *δ* is determined by the equation (18), and the displacement *u2* is determined by MEF. The non-dimensionalized maximum crack opening displacement *δ* can be obtained using *ρ*, as shown in the expression (20):

$$\mathcal{S} = \mathcal{c}\_1 \left( e^{-a\_1 \rho} \right) + \mathcal{c}\_2 \left( e^{-b\_1 \rho} \right) + \mathcal{c}\_3 \tag{20}$$

The constants *a1, b1, c1, c2* e *c3* in the expression (20) depends on the type of the used material. The table 1 exposes the value of these constants for the laminate glass/epoxi (Gl/Ep).


Table 1. Coefficients for the expression (20) (Lim & Tay, 1996)

As the internal variable used in this problem depends on the maximum crack opening displacement according to equation (15), and the crack density is calculated by 1 /*L* , it can be shown that the state variable associated to the Mode I becomes:

$$\text{Mode I:} \qquad a\_{22} = \frac{8}{5} \mu\_2 \zeta \tag{21}$$
