**3.1 The Modified Local Green's Function Method - MLGFM**

The Modified Local Green's Function Method (MLGFM) is an integral technique that associates the Finite Element Method and the Boundary Element Method, solving the problem through an integral equations system. Unlike to the BEM and the Trefitz Methods, the MLGFM does not use a fundamental solution and/or a Green's function. The term "Local" indicates that the calculation of the GF projections can be done locally, that is, for each element.

Essentially, the MLGFM uses a transverse integration technique and reciprocity relations to determine, at a local level, the Green's Function, transforming the partial differential operator in an ordinary partial operator (Barcelos & Silva, 1987). The MLGFM uses finite elements at the domain with the purpose to create discrete projections of the Green's Function, corresponding to fundamental solutions, that are used later in the integral equations system associated to the boundary approximation. To analyze a continuum mechanics problem through the MLGFM, such as the plate indicated in figure 4, two meshes are necessary, one for the domain and other for the boundary. With domain elements, the method generates a set of domain equations, which are used to generate automatically the domain Green's projections. Later, a set of boundary equations are also generated and the boundary Green's projections can be determined with the domain projections developed before. At the end, the system is solved only for boundary equations, where the main variables are calculated. Domain values may be obtained once the boundary values are known after the solution of the boundary equation system.

Considering a conventional structural approximation by conventional Finite Element Method, the problem can be expressed by a system of algebraic equations, representing a

> 1 1 11 12 16 21 22 26 2 2 61 62 66 6 6

where *Kij* are the element stiffness matrix, (*dx*, *dy*, *dz*) are the components of the element displacement vector and {*Fa*} and {*Fd*} are the applied force vector and the element damage

The procedure used in this paper to obtain the expected results is a little different because it uses a different computational method known as Modified Local Green's Function Method (MLGFM), in witch the system defined in expression (22) is not directly applied in a conventional FEM. A detailed explanation of this procedure can be found in Barbieri *et al* (1998a,b) and Machado *et al* (2008). The MLGFM is an integral method that determines the unknowns on the boundary, similarly to the Boundary Element Method, but the fundamental solutions are generated automatically by projections of the Green's Functions developed from de field, as in the Finite Element Method. The matrix and the vectors indicated in (22) will be

The Modified Local Green's Function Method (MLGFM) is an integral technique that associates the Finite Element Method and the Boundary Element Method, solving the problem through an integral equations system. Unlike to the BEM and the Trefitz Methods, the MLGFM does not use a fundamental solution and/or a Green's function. The term "Local" indicates that the calculation of the GF projections can be done locally, that is, for

Essentially, the MLGFM uses a transverse integration technique and reciprocity relations to determine, at a local level, the Green's Function, transforming the partial differential operator in an ordinary partial operator (Barcelos & Silva, 1987). The MLGFM uses finite elements at the domain with the purpose to create discrete projections of the Green's Function, corresponding to fundamental solutions, that are used later in the integral equations system associated to the boundary approximation. To analyze a continuum mechanics problem through the MLGFM, such as the plate indicated in figure 4, two meshes are necessary, one for the domain and other for the boundary. With domain elements, the method generates a set of domain equations, which are used to generate automatically the domain Green's projections. Later, a set of boundary equations are also generated and the boundary Green's projections can be determined with the domain projections developed before. At the end, the system is solved only for boundary equations, where the main variables are calculated. Domain values may be obtained once the boundary values are

used to produce values at the boundary, as explained in the next topic.

**3.1 The Modified Local Green's Function Method - MLGFM**

known after the solution of the boundary equation system.

*KKK d F F KKK d F F KKK d F F* 

*x*

*y*

*a d*

*a d*

(22)

*a d <sup>z</sup>*

typical element:

force vector.

each element.

Fig. 4. Symmetric boundary conditions and plate for a 2x2 mesh (4 finite elements and 8 contour elements) by the MLGFM.

The most important steps of the MLGFM are detailed in the work sof Barbieri *et al* (1998a) and Machado *et al* (2008). It is possible to show that through the MLGFM two sets of equations are formed, the first one in the domain (equation (23)) and the other one on the boundary (equation (24)):

$$\mathbf{u}(\mathbf{Q}) = \int\_{\Omega} \left[ \mathbf{G}^{\mathrm{T}}(\mathbf{P}, \mathbf{Q}) \, \mathbf{a}(\mathbf{P}) \right] d\Omega + \int\_{\Gamma} \left[ \mathbf{G}^{\mathrm{T}}(\mathbf{p}, \mathbf{Q}) \, \mathbf{f}(\mathbf{p}) \right] d\Gamma \, ; \; \mathbf{P}, \mathbf{Q} \in \Omega \; ; \mathbf{p} \in \Gamma \tag{23}$$

$$\mathbf{u}(\mathbf{q}) = \int\_{\Omega} \left[ \mathbf{G}^{\mathrm{T}} \{ \mathbf{P}, \mathbf{q} \} \, \mathbf{a}(\mathbf{P}) \right] d\Omega + \int\_{\Gamma} \left[ \mathbf{G}^{\mathrm{T}} \{ \mathbf{p}, \mathbf{q} \} \, \mathbf{f}(\mathbf{p}) \right] d\Gamma \, ; \; \mathbf{P} \in \Omega \ ; \; \mathbf{p}, \mathbf{q} \in \Gamma \tag{24}$$

where Q, P are two points in the domain; q, p are other two points on the boundary; **a**(P) is the vector of independent terms for the original problem; f(p) is the vector associated to the fluxes on the boundary; **G**(i,j) are the Green's functions which may be understood as the generalized displacement in the point i in the direction of an unitary vector **ni**, when a generalized force is applied over the point j, in the direction of a unitary vector **n**j.

Equations (23) and (24) describe completely the problem. Since these equations involve domain and boundary integrals, two types of meshes are necessary, one in the domain and the other on the boundary, using FE and BE methods, respectively. The FE domain approximation is also used to develop the Green's functions which are associated to the matrices **G(**P,Q**), G(**p,Q**), G(**P,q**)** and **G(**p,q**).** 

Progressive Stiffness Loss Analysis

where

following functional **F**, which one depends on **G** or **G** :

**B** – is a bilinear form, developed to **G** or **G;**  and are constants whose values are: = 1 and = 0 to determine **G** = 0 and = 1 to determine **G**

 **B1(G,[])** =

 **B2(G,[])** =

(26), which are the main system of the MLGFM.

 **B3(G,G)** = <sup>1</sup>

Green's projections

**4. Damage evolution** 

where *Fa* is the failure criterion,

**B1**, **B2**, **B3** – are bilinear forms which can be written as:

**G** – corresponds to **G** or **G** depending on the case of interest;

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

In order to determine the Green's functions automatically, it must be considered the

 **F(G,G) = B(G,G) - B1(G,) - B2(G, ) + B3(G,G)** (37)

As in the variational approach of the FEM, the minimization of functional **F(G,G)** in Equation (37) results in a linear equation system which can be solved to determine the

where [**K**] is the global stiffness matrix, evaluated in the same way of the conventional finite element stiffness matrix; **A** and **D** are the matrices of Equations (27) and (30), respectively. In this way, the Green's projections are determined directly from Equation (41), and can be applied in Equations (28), (29), (31) and (32) to complete the matrices of equations (25) and

With purpose to quantify the damage accumulation due to a monotonic load increment,

*ij a ij F Z X* 

and *Z* is the failure value characteristic to each criterion. The rigidity degradation of a component occurs due to a progressive process during its serviceable life. It is important to note that the evolution of rigidity loss in a structure can be characterized by a single crack or by the occurrence of generalized cracks. Combinations of failure modes can act together

some failure criterion will be used. Generically, failure criteria can be considered as:

2  **<sup>G</sup>(**Q**) (**Q**)** d (38)

**<sup>G</sup>(**q**) (**q**)** d (39)

(42)

**N#(**q**) G (**q**). G (**q**)** d (40)

**KG G A D Ω Γ** (Q) (Q) (41)

*ij* are the local stress, *Xij* are the principal material strength

The approximation shape functions are the same as in the conventional FE and BE methods. For the present work, quadratic shape functions were employed to construct nine nodes lagrangean finite elements and three nodes boundary elements.

Developing discrete equations from the nodal values, two sets of linear equations are determined:

$$\mathbf{A \ u\_{i\bullet} = B \ f + C \ a \quad \text{(in the domain)}\tag{25}$$

$$\mathbf{D \ u\_{l} = E \ f + F \ a \quad \text{(on the boundary)}\tag{26}$$

where **u** and **u** are the domain and the boundary displacements, respectively, **a** and **f** are the independent and the fluxes variables vectors. The matrices **A**, **B**, **C**, **D**, **E** and **F** can be written as:

$$\mathbf{A} = \int\limits\_{\Omega} \Psi(\mathbf{Q})^{\mathrm{T}} \Psi(\mathbf{Q}) \mathrm{d}\Omega \tag{27}$$

$$\mathbf{B} = \int\limits\_{\Omega} \Psi(\mathbf{Q})^{\mathrm{T}} \mathbf{G}\_{\mathbf{T}}(\mathbf{Q}) d\Omega \tag{28}$$

$$\mathbf{C} = \int\limits\_{\Omega} \Psi(\mathbf{Q})^{\mathrm{T}} \, \mathbf{G}\_{\mathbf{D}}(\mathbf{Q}) d\Omega \tag{29}$$

$$\mathbf{D} = \int\_{\Gamma} \Phi(\mathbf{q})^{\mathbf{r}} \Phi(\mathbf{q}) d\Gamma \tag{30}$$

$$\mathbf{E} = \int\_{\Gamma} \Phi(\mathbf{q}) \mathbf{r} \, \mathbf{G}\_{\mathbf{T}}(\mathbf{q}) d\Gamma \tag{31}$$

$$\mathbf{F} = \int\limits\_{\Gamma} \Phi(\mathbf{q})^{\mathbf{T}} \mathbf{G}\_{\mathbf{d}\mathbf{d}}(\mathbf{q}) d\Gamma \tag{32}$$

where **(**Q**)** and **(**q**)** are matrices with the shape functions in the domain and on the boundary, respectively, and **G(**Q**), G(**Q**), G(**q**), G(**q**)** are the Green's function projections over the boundary and the domain , evaluated on the points Q and q. The Green's projections can be written as:

$$\mathbf{G\_{T}(Q) = \int\limits\_{\Gamma} G^{\tau}(p, Q) \,\Phi(p) \,\mathrm{d}\Gamma} \tag{33}$$

$$\mathbf{G\_{D}(Q) = \int\limits\_{\Omega} \mathbf{G^{T}(P,Q)} \,\Psi(P) \,\mathrm{d}\Omega \,\tag{34}$$

$$\mathbf{G\_{T(q)}} = \int\_{\Gamma\_{\parallel}} \mathbf{G\_{T(p,q)}} \, \Phi(\mathbf{p}) \, \mathrm{d}\Gamma \,\tag{35}$$

$$\mathbf{G\_{f2}(q) = \int\_{\Omega} \colon \mathbf{G^{T}(P,q)} \,\Psi(P) \,\mathrm{d}\Omega \,\tag{36}$$

In order to determine the Green's functions automatically, it must be considered the following functional **F**, which one depends on **G** or **G** :

$$\mathbf{F(G\_{D},G\_{T}) = B(G,G) - \alpha \ B\_{1}(G\_{D}, \mathbf{y}) - \beta \ B\_{2}(G\_{T}, \mathbf{d}) + B\_{3}(G, G)}\tag{37}$$

where

132 Continuum Mechanics – Progress in Fundamentals and Engineering Applications

The approximation shape functions are the same as in the conventional FE and BE methods. For the present work, quadratic shape functions were employed to construct nine nodes

Developing discrete equations from the nodal values, two sets of linear equations are

 **A u = B f + C a** (in the domain) (25)

 **D u = E f + F a** (on the boundary) (26) where **u** and **u** are the domain and the boundary displacements, respectively, **a** and **f** are the independent and the fluxes variables vectors. The matrices **A**, **B**, **C**, **D**, **E** and **F** can be

(Q)T (Q)d (27)

**(**Q**)T G(**Q**)**d (28)

**(**Q**)T G(**Q**)**d (29)

**(**q**)T (**q**)**d (30)

**(**q**)T G(**q**)**d (31)

**(**q**)T G(**q**)**d (32)

**GT(**p,Q**) (**p**)** d (33)

**GT(**P,Q**) (**P**)** d (34)

**GT(**p,q**) (**p**)** d (35)

**GT(**P,q**) (**P**)** d (36)

lagrangean finite elements and three nodes boundary elements.

**A** = 

**G(**Q**)** =

where **(**Q**)** and **(**q**)** are matrices with the shape functions in the domain and on the boundary, respectively, and **G(**Q**), G(**Q**), G(**q**), G(**q**)** are the Green's function projections over the boundary and the domain , evaluated on the points Q and q. The Green's

determined:

written as:

 **B =** 

 **C =** 

 **D =** 

 **E =** 

 **F =** 

 **G(**Q**) =** 

 **G(**q**) =** 

 **G(**q**) =** 

projections can be written as:

**G** – corresponds to **G** or **G** depending on the case of interest;

**B** – is a bilinear form, developed to **G** or **G;** 

and are constants whose values are:

= 1 and = 0 to determine **G**

= 0 and = 1 to determine **G**

**B1**, **B2**, **B3** – are bilinear forms which can be written as:

$$\mathbf{B}\_1(\mathbf{G}\_\mathbf{D}[\Psi]) = \int\_{\Omega} \mathbf{G}\_\mathbf{B}(\mathbf{Q}) \, \Psi(\mathbf{Q}) \, \mathrm{d}\Omega \tag{38}$$

$$\mathbf{B\_2(G\_{\sf T}[\Phi]) = \int\_{\sf T} G\_{\sf T}(\mathbf{q}) \, \boldsymbol{\Phi}(\mathbf{q}) \, d\sf{T} \tag{39}$$

$$\mathbf{B}\_3(\mathbf{G}, \mathbf{G}) = \frac{1}{2} \int\limits\_{\Gamma} \mathbf{N}^\bullet(\mathbf{q}) \, \mathbf{G} \text{ (q)}. \, \mathbf{G} \text{ (q)} \, \text{d}\Gamma \tag{40}$$

As in the variational approach of the FEM, the minimization of functional **F(G,G)** in Equation (37) results in a linear equation system which can be solved to determine the Green's projections

$$\begin{bmatrix} \mathbf{K} \end{bmatrix} \begin{bmatrix} \mathbf{G}\_{\Omega}(\mathbf{Q}) \ \vdots \ \mathbf{G}\_{\Gamma}(\mathbf{Q}) \end{bmatrix} = \begin{bmatrix} \mathbf{A} & \vdots & \mathbf{D} \end{bmatrix} \tag{41}$$

where [**K**] is the global stiffness matrix, evaluated in the same way of the conventional finite element stiffness matrix; **A** and **D** are the matrices of Equations (27) and (30), respectively. In this way, the Green's projections are determined directly from Equation (41), and can be applied in Equations (28), (29), (31) and (32) to complete the matrices of equations (25) and (26), which are the main system of the MLGFM.

#### **4. Damage evolution**

With purpose to quantify the damage accumulation due to a monotonic load increment, some failure criterion will be used. Generically, failure criteria can be considered as:

$$F\_a \left(\frac{\sigma\_{ij}}{X\_{ij}}\right) = Z$$

where *Fa* is the failure criterion, *ij* are the local stress, *Xij* are the principal material strength and *Z* is the failure value characteristic to each criterion. The rigidity degradation of a component occurs due to a progressive process during its serviceable life. It is important to note that the evolution of rigidity loss in a structure can be characterized by a single crack or by the occurrence of generalized cracks. Combinations of failure modes can act together

Progressive Stiffness Loss Analysis

plane stress condition, resulted in the expression (43):

**4.5 Strain Energy Release Rate Criterion** 

O

 

 

Fig. 5. Stiffness loss in composite laminates [0°l/90°m]s – Lim & Tay (1995).

The result is a reduction in the effective stiffness of the laminate in the loading direction, and this is represented by the OB segment in figure 5. Upon further loading, this reduction is verified by the segment BC. This process is repeated until line OF is reached. Note that

**4.3 Tsai-Hill Criterion** 

**4.4 Tsai-Wu Criterion** 

Tsai-Hill Criterion.

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

An adaptation made by Tsai in the Hill Criterion for transverse orthotropic laminate at

 

 

(43)

(44)

22 2 1 2 1 2 12 22 2 2 1 *XY X C*

A simple procedure was proposed by Tsai-Wu, changing the Tsai-Hill Criterion in equation (43). When the tensile and compression strength are similar, the expression (44) becomes the

> 22 2 1 2 1 2 12 22 2 1 *XY C XY*

The Strain Energy Release Rate Criterion to describe the damage evolution was applied by Lim and Tay (Lim & Tay, 1996). Consider a symmetric laminated composite of width *b* and length *L*, with the configuration [0°l/90°m]s, where l and m are integers. When the laminated is loaded uniaxially in tension, the stress-strain curve is linear until the failure criterion is reached for the first time, at point A (figure 5). A transverse crack is introduced in 90° layer.

C

E

F

D

H

A

B

G

causing changes in the material properties and in its local stress distribution. In this way, the main difficulty in this kind of analysis is the adoption of a failure criterion, *Fa* , that conveniently describes the damage evolution due to a failure mode.

The theories introduced to prevent the failure of an orthotropic laminate are adaptations of failure criterion for isotropic materials, modified for biaxial stress cases, such as, Maximum Stress Criterion, Maximum Strain Criterion, Tsai-Hill Criterion and Tsai-Wu Criterion (Reddy, 1997; Vasiliev & Morozov, 2001; Mendonça, 2005).

A new criterion is presented, based on the strain energy release rate, to evaluate the formation of a new micro crack (Anderssen et al., 1998; Ji et al., 1998; Kobayashi et al.,2000). The released energy is used because it is practically independent from the crack length (Anderssen et al., 1998). Some of these criterions are presented here.
