# Reduction of Free Edge Peeling Stress of Laminated Composites Using Active Piezoelectric Layers

Authors: Bin Huang; Heung Soo Kim
Journal: The Scientific World Journal
 (2014)

---

## Abstract

An analytical approach is proposed in the reduction of free edge peeling stresses of laminated composites using active piezoelectric layers. The approach is the extended Kantorovich method which is an iterative method. Multiterms of trial function are employed and governing equations are derived by taking the principle of complementary virtual work. The solutions are obtained by solving a generalized eigenvalue problem. By this approach, the stresses automatically satisfy not only the traction-free boundary conditions, but also the free edge boundary conditions. Through the iteration processes, the free edge stresses converge very quickly. It is found that the peeling stresses generated by mechanical loadings are significantly reduced by applying a proper electric field to the piezoelectric actuators.


---

## Body

## 1. Introduction

Laminated composite materials [1] have been increasingly attracted in the engineering industries, due to their high stiffness-to-weight ratio, high strength-to-weight ratio, high corrosion resistance, high chemical stability, and so on. They have been popularized in the aircraft manufacturing, vehicle industries, medical fields, and even sport equipment industries. With the improvement of manufacturing process and productivity, they have potential in replacing the conventional metal alloys and nonmetal materials massively in the future. However, it is well known that due to the mismatch of elastic properties between the adjacent layers, there exist large undesirable interlaminar peeling stresses. Because of these highly concentrated peeling stresses, failures such as delamination or crack may occur near the free edges and further lead to the damage of structures.The concentrations of peeling stress often occur within a short distance from the free edges, known as the free edge boundary effect [2, 3]. The stresses decrease and vanish in the interior region of the laminates. Also, even an in-plane load applied on the laminate will cause 3D state of stress distributions near the free edge. Thus, some theories are not capable of predicting the free edge peeling stress distributions which assume the plane stress condition, such as classical laminate theory (CLT) [4, 5].The concentrated free edge peeling stresses are critical to the fatigues of laminated structures. Some passive control methods have been developed either by adjusting the stacking sequence of laminates [6] or by changing the edge shape of fibers [7] to prevent the onset of delamination. Nevertheless, due to their passive nature, they cannot work effectively and other desired goals cannot be achieved. Thus, in this study, we will investigate the peeling stress reduction by using piezoelectric actuators. Owing to their electromechanical coupling nature, piezoelectric actuators have potential to generate stresses which have opposite signs of those peeling stresses generated by mechanical loadings and, in return, to reduce the magnitudes of those peeling stresses.Since the last decades, a large volume of literature sources has introduced the modeling of smart composite laminated plates and beams [8–12]. Chen and Qiao investigated the electromechanical behavior and fracture of piezoelectric bilayer beams [13–15]. They also investigated the interlaminar stresses of plated beam with adhesive layers [16]. They found that the adhesive layers can affect the local deformations and in return alleviate the stress concentrations. While predicting the interlaminar stresses of laminated structures, most approaches neglect the adhesive layer effect. Numerous approximate approaches, without considering the adhesive layer effect, have been developed since there is no method of obtaining exact 3D elastic solutions. The developed theories vary from the displacement field based approaches [17] to stress field based approaches [18, 19]. All of these theories have their own characteristics and specific applications.For the displacement field based equivalent single layer theories, such as classical laminate theory, first order shear deformation theory [20, 21], and higher order shear deformation theory [22–26], these theories are efficient in predicting the global response of laminates in bending, buckling, and vibration analysis, while they are not capable of predicting the localized free edge stress concentrations. The stress components are usually calculated from the constitutive relation or the equilibrium equations. Thus, the free edge stress concentration at the layer interface cannot be revealed accurately. For the stress function based equivalent single layer theories [27, 28], the localized free edge stress concentrations can be accurately predicted, while, on the other hand, displacement information cannot be achieved. Moreover, prescribed boundary conditions are mandatory in those equivalent single layer theories either for the preassumed displacement trial functions or for the stress trial functions. The layerwise theories [29] are much more accurate in analyzing the stress concentration at the ply level. For the displacement based layerwise theories, they allow the in-plane displacement to vary continuously through the laminate thickness, while their derivatives with respect to the thickness coordinate may be discontinuous through the thickness, thereby allowing for the possibility of continuous transverse stresses at the ply interfaces. For the stress based layerwise theories [30, 31], the stress fields are assumed layerwisely and transverse stress continuity conditions have to be enforced at the ply interfaces. Both the displacement and stress based theories can accurately predict the stress concentration at the free edges. However, the layerwise theories are much more computationally expensive than the equivalent single layer theories.Among the analytical methods, another approximate method named as the extended Kantorovich method [32, 33] can be referred to. The extended Kantorovich method is an iterative method for solving partial differential equations. Although the initial assumption of admissible functions influences the convergence time a lot, the solutions can be improved gradually by repeating solving the governing equations and applying boundary and continuity conditions during the iteration processes. The assumed admissible functions with prescribed boundary conditions are of priority and any types of boundary condition problems can be solved by using this approach. The interlaminar stresses at the straight free edges of laminates under extension loads were investigated by Cho and Yoon [34]. And later Cho and Kim [35] analyzed the free edge interlaminar stresses of composite laminates under extension, bending, twisting, and thermal loads by using the extended Kantorovich method. In their study, they preferred using the stress trial functions, and traction-free boundary conditions and free edge boundary conditions were satisfied automatically during the iterative processes. Ungbhakorn and Singhatanadgid [36] proposed an extended Kantorovich method to investigate the buckling problem of rectangular laminates with various boundary conditions. The initial trial functions were used arbitrarily and the principle of minimum total potential energy was used to derive governing equations. Andakhshideh and Tahani [37] investigated the interlaminar stresses near free edges of finite length general composite laminates subjected to axial and shearing loads using three-dimensional extended Kantorovich method. They provided convergent 3D elasticity solutions with four terms of trial functions. Due to their arbitrary choice of initial trial functions and various boundary conditions, the extended Kantorovich method could be a good methodology in the interlaminar stress analysis of laminated composite structures.In the present study, a stress function based extended Kantorovich method is proposed to investigate the reduction of free edge peeling stresses using piezoelectric actuators. Symmetrically layered laminates will be taken into consideration as symmetric laminates are most commonly used to avoid coupling of bending and twisting modes. Most studies start from the initial assumption that the strain or stress is independent of the longitudinal coordinate, and this study also adopts this point. The three-dimensional stress solutions are obtained by iterative processes involving multiterms of initially assumed trial functions. By the extended Kantorovich method, the influence of initial assumptions can be well eliminated and the accuracy of stress distributions can be improved a lot. Finally, convergent results are given as numerical examples and compared with those obtained by finite element method. It is expected to obtain significant peeling stress reduction by applying proper electric fields to piezoelectric actuators.


## 2. Mathematical Modeling

### 2.1. Derivation of Governing Equation

For the given geometry of laminated composite plate, as shown in Figure1, two piezoelectric actuators are symmetrically surface bonded on the top and bottom surfaces of the laminates without considering the adhesive layer affect. This kind of structure is usually considered to be smart composite laminates [38]. Two piezoelectric layers are used to control the free edge peeling stresses that are generated by mechanical loadings. The longitudinal dimension of the laminates is considered long enough so that the stress fields are independent of the x-axis. Therefore, generalized plane strain states are assumed to be in the x direction.Figure 1
Geometry of laminated composite plate with surface bonded piezoelectric actuators under combined electromechanical loadings.The extended Kantorovich method requires iterative calculation methodology. By preassumption of trial functions and solving a set of ordinary differential equations, the solution of the first process can be obtained. Due to the rough initial assumption, it cannot guarantee the accuracy of solution. Thus, to overcome this weak point, iterative processes can be conducted to improve the accuracy of solution. The governing equations can be derived by using either Galerkin’s method or the principle of virtual work. In this work, the principle of complementary virtual work is employed. The complementary strain energy of the laminated composite plate is calculated by the following equation:(1)
δ
U
=
∬
s
ε
i
δ
σ
i
d
A
i
=
∬
s
ε
2
δ
σ
2
+
ε
3
δ
σ
3
+
ε
4
δ
σ
4
+
ε
5
δ
σ
5
+
ε
6
δ
σ
6
d
A
=
0
.Before substituting strain and stress components into the complementary strain energy equation, an initially assumed stress field should be introduced. Instead of using displacement assumption, stress field assumption is preferred. A good choice for stress field is Lekhnitskii stress functions [39] which can express the three-dimensional stress state of piezo-bonded composite laminates. Lekhnitskii stress functions can automatically satisfy the pointwise equilibrium equations as follows:
(2)
σ
2
=
∂
2
F
∂
η
2
,
σ
3
=
∂
2
F
∂
ξ
2
,
σ
4
=
-
∂
2
F
∂
ξ
∂
η
σ
5
=
-
∂
ψ
∂
ξ
,
σ
6
=
∂
ψ
∂
η
,
where η
=
z
/
H and ξ
=
y
/
H are nondimensionalized coordinates and H is the total thickness of the laminates. The solutions F and ψ are series combinations and assumed to be a combination of in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)) and out-of-plane stress functions (g
i
(
η
)) using separation of variables as follows:
(3)
F
=
∑
i
=
1
n
f
i
(
ξ
)
g
i
(
η
)
,
ψ
=
∑
i
=
1
n
p
i
(
ξ
)
g
i
I
(
η
)
.With the initial assumption of stress field, the strain tensor{
ε
i
} can be found from the linear elastic constitutive relation. This relation is assumed for each individual layer with piezo induced strain term:
(4)
{
ε
i
}
(
k
)
=
[
S
-
]
(
k
)
{
σ
i
}
(
k
)
+
[
d
]
(
k
)
{
E
}
(
k
)
,
(
i
=
1,2
,
…
,
6
)
,
where [
S
-
] is the generalized compliance matrix for orthotropic materials. [
d
]  is the piezoelectric constant matrix and [
E
] is the applied electric field.As the mechanical extension load is applied in thex direction, all other strain components can be expressed using ε
1 and applied piezoelectric fields. Consider
(5)
ε
i
=
S
^
i
j
σ
j
+
S
i
1
S
11
ε
1
+
d
^
i
3
E
3
,
(
i
=
2,3
,
…
,
6
,
j
=
2,3
,
…
,
6
)
,
where
(6)
S
^
i
j
=
S
i
j
-
S
i
1
S
1
j
S
11
,
d
^
i
3
=
d
i
3
-
S
i
1
S
11
d
13
.
Substituting (2) and (5) into (1), the complementary strain energy becomes
(7)
δ
U
=
∬
s
ε
i
δ
σ
i
d
A
i
=
∬
s
[
(
S
^
i
j
σ
j
+
S
i
1
S
11
ε
1
+
d
^
3
i
E
3
)
δ
σ
i
]
d
ξ
d
η
=
0
,
mmmmmmmmmmmmmmmm
(
i
,
j
=
2,3
,
…
,
6
)
.To simplify the above equation, integration by parts and free edge boundary conditions (σ
2
=
σ
4
=
σ
6
=
0, at ξ
=
0) are applied. Then, the double integration is simplified to the following equation and the governing equations can be obtained by using the variational principle. Consider
(8)
∫
[
a
i
j
(
4
)
f
j
IV
+
a
i
j
(
2
)
f
j
II
+
a
i
j
(
0
)
f
j
+
b
i
j
(
2
)
p
j
II
+
b
i
j
(
0
)
p
j
+
r
i
]
δ
f
i
d
ξ
+
∫
[
d
i
j
(
2
)
f
j
II
+
d
i
j
(
0
)
f
j
+
c
i
j
(
2
)
p
j
II
+
c
i
j
(
0
)
p
j
+
s
i
]
δ
p
i
d
ξ
=
0
,
mmmmmmmmiimmmmmmm
(
i
,
j
=
1,2
,
…
,
n
)
,
where
(9)
a
i
j
(
4
)
=
∫
S
^
33
g
i
g
j
d
η
,
a
i
j
(
2
)
=
∫
S
^
23
(
g
i
II
g
j
+
g
i
g
j
II
)
d
η
-
∫
S
^
44
g
i
I
g
j
I
d
η
,
a
i
j
(
0
)
=
∫
S
^
22
g
i
II
g
j
II
d
η
,
b
i
j
(
2
)
=
∫
S
^
36
g
i
g
j
II
d
η
-
∫
S
^
45
g
i
I
g
j
I
d
η
,
b
i
j
(
0
)
=
∫
S
^
26
g
i
II
g
j
II
d
η
,
c
i
j
(
2
)
=
-
∫
S
^
55
g
i
I
g
j
I
d
η
,
c
i
j
(
0
)
=
∫
S
^
66
g
i
II
g
j
II
d
η
,
d
i
j
(
2
)
=
∫
S
^
36
g
i
II
g
j
d
η
-
∫
S
^
45
g
i
I
g
j
I
d
η
,
d
i
j
(
0
)
=
∫
S
^
26
g
i
II
g
j
II
d
η
,
r
i
=
∫
(
S
21
S
11
ε
1
+
d
^
23
E
3
)
g
i
II
d
η
,
s
i
=
∫
(
S
61
S
11
ε
1
+
d
^
63
E
3
)
g
i
II
d
η
.


### 2.2. Solutions of Governing Equations and Generalized Eigenvalue Problem

To solve the fourth-order and second-order coupled ordinary differential equations, the trial functionsg
i
(
η
)  are specified first in (10). Although prescribed boundary conditions are not necessarily satisfied for trial functions when using the extended Kantorovich method, the mode shapes of the clamped-clamped beam are used as the initial out-of-plane trail functions to satisfy the traction-free boundary conditions automatically (σ
3
=
σ
4
=
σ
5
=
0
at
η
=
±
1
/
2). Consider
(10)
g
i
(
η
)
=
cos
⁡
(
β
i
η
)
+
cos
⁡
(
β
i
/
2
)
cosh
⁡
(
β
i
/
2
)
cosh
⁡
(
β
i
η
)
,
where β
i are the roots of the following characteristic equation. Consider
(11)
tan
(
β
2
)
+
tanh
(
β
2
)
=
0
.Thus, all coefficients in the governing equations can be obtained and the general solutions of governing equations can be assumed as the sum of homogeneous solutions(
f
i
(
H
)
and
p
i
(
H
)
) and particular solutions (
f
i
(
P
)
and
p
i
(
P
)
). The homogeneous solutions of f
i
(
H
) and p
i
(
H
) are assumed to be the following form:
(12)
f
i
(
H
)
=
v
i
f
e
λ
ξ
,
p
i
(
H
)
=
v
i
p
e
λ
ξ
.Substituting (12) into the governing equations and defining v
j
f
II
=
λ
2
v
j
f, the governing equations can be rewritten in the matrix form and a standard eigenvalue problem is formed. Consider
(13)
[
0
-
a
i
j
(
4
)
-
b
i
j
(
2
)
0
0
c
i
j
(
0
)
I
0
0
]
-
1
[
a
i
j
(
0
)
a
i
j
(
2
)
b
i
j
(
0
)
d
i
j
(
0
)
d
i
j
(
2
)
b
i
j
(
0
)
0
I
0
]
[
v
j
f
v
j
f
II
v
j
p
]
=
λ
2
[
v
j
f
v
j
f
II
v
j
p
]
,
(
i
,
j
=
1,2
,
…
,
n
)
,
where λ
2 are eigenvalues of the above equation. Since the interlaminar stresses decay in the interior region of piezo-bonded composite laminates, only negative roots of λ
2 are selected. Equation (13) results in 3n eigenvalues and eigenvectors, so that the homogeneous solutions consist of 3n terms linear combination as follows:
(14)
f
i
(
H
)
=
t
j
v
i
j
f
e
-
λ
j
ξ
,
p
i
(
H
)
=
t
j
v
i
j
p
e
-
λ
j
ξ
,
(
i
=
1,2
,
…
,
n
)
,
(
j
=
1,2
,
…
,
3
n
)
,
where t
j  are constants that will be determined by the prescribed boundary conditions.The particular solutions are obtained by consideringf
i
(
P
) and p
i
(
P
) as constants; thus all their derivatives are zeros in the governing equations and the remaining parts can be rewritten in the following equation:
(15)
[
f
j
(
P
)
p
j
(
P
)
]
=
-
[
a
i
j
(
0
)
b
i
j
(
0
)
b
i
j
(
0
)
c
i
j
(
0
)
]
-
1
[
r
i
s
i
]
(
i
,
j
=
1,2
,
…
,
n
)
.The final step of solving the governing equations is to obtain the constantst
j in the homogeneous solutions. They can be calculated by enforcing the free edge boundary conditions (σ
2
=
σ
4
=
σ
6
=
0
at
ξ
=
0), and this will lead f
i,  f
i
I, and p
i to be zeros at the free edges. Consider
(16)
t
j
v
i
j
f
e
-
λ
j
ξ
|
ξ
=
0
+
f
i
(
P
)
=
0
,
-
t
j
λ
j
v
i
j
f
e
-
λ
j
ξ
|
ξ
=
0
=
0
,
t
j
v
i
j
p
e
-
λ
j
ξ
|
ξ
=
0
+
p
i
(
P
)
=
0
.After solving the governing equations, all stress components can be obtained by substituting the solutions into the Lekhnitskii stress functions and they have the following expressions:(17)
σ
2
=
∑
i
=
1
n
∑
j
=
1
3
n
(
t
j
v
i
j
f
e
-
λ
j
ξ
+
f
i
(
P
)
)
g
i
II
,
σ
3
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
i
2
t
j
v
i
j
f
e
-
λ
j
ξ
g
i
,
σ
4
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
j
t
j
v
i
j
f
e
-
λ
j
ξ
g
i
I
,
σ
5
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
j
t
j
v
i
j
p
e
-
λ
j
ξ
g
i
I
,
σ
6
=
∑
i
=
1
n
∑
j
=
1
3
n
(
t
j
v
i
j
p
e
-
λ
j
ξ
+
p
i
(
P
)
)
g
i
II
.


### 2.3. Iterative Procedures

The accuracy of the first step largely depends on the initial assumption of trail functions and the number of series terms involved in the trail functions. The results may also be influenced by the choice of trail functions which may cause singularity problem, harmonic distributions of stresses, and so forth. Thus, it is important to overcome the weaknesses appeared in the first process. We prefer using the extended Kantorovich method which is an iterative process. Through the iterative processes, the influence of the initial assumption can be eliminated and the convergent results can be obtained.In the second step, the solutions of the first step (f
i
(
ξ
) and p
i
(
ξ
)) are considered to be the known functions and the principle of complementary virtual work is used again to derive governing equations. Here, the unknown out-of-plane stress functions (g
i
(
η
)) are assumed to be layer dependent and the continuity conditions are enforced at the layer interface. Thus, from the principle of complementary virtual work in (7) and integration by parts, the governing equation for the out-of-plane stress functions can be obtained as follows:
(18)
∫
[
m
i
j
(
4
)
(
k
)
g
j
(
k
)
(
IV
)
+
m
i
j
(
2
)
(
k
)
g
j
(
k
)
(
II
)
+
m
i
j
(
0
)
(
k
)
g
j
(
k
)
+
x
i
(
k
)
]
×
δ
g
j
(
k
)
d
η
+
Λ
=
0
(
i
,
j
=
1,2
,
…
,
n
)
,
mmmmmmmmmmm
(
k
=
1,2
,
…
,
N
)
,
where Λ represents the boundary terms induced by integration by parts and the superscript (k) refers to the kth ply. The layerwise coefficients are listed as follows:
(19)
m
i
j
(
4
)
(
k
)
=
∫
[
S
^
22
(
k
)
f
i
f
j
+
S
^
26
(
k
)
(
f
i
p
j
+
p
i
f
j
)
+
S
^
66
(
k
)
p
i
p
j
]
d
ξ
,
m
i
j
(
2
)
(
k
)
=
∫
[
S
^
23
(
k
)
(
f
i
f
j
I
I
+
f
i
I
I
f
j
)
+
S
^
36
(
k
)
(
f
i
I
I
p
j
+
p
i
f
j
I
I
)
m
m
m
m
m
i
-
S
^
44
(
k
)
f
i
I
f
j
I
-
S
^
45
(
k
)
(
f
i
I
p
j
I
+
p
i
I
f
j
I
)
-
S
^
55
(
k
)
p
i
I
p
j
I
]
d
ξ
,
m
i
j
(
0
)
(
k
)
=
∫
S
^
33
(
k
)
f
i
I
I
f
j
I
I
d
ξ
,
x
i
(
k
)
=
∫
(
S
13
(
k
)
S
11
(
k
)
ε
1
+
d
^
33
E
3
)
f
i
I
I
d
ξ
.The solution procedure is the same as in the first step. A generalized eigenvalue problem is solved and the general solution is calculated by adding the homogeneous solutionsg
i
(
k
)
(
H
) and the particular solutions g
i
(
k
)
(
P
):
(20)
g
i
(
k
)
=
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
+
g
i
(
k
)
(
P
)
,
(
k
=
1,2
,
…
,
N
)
,
(
i
=
1,2
,
…
,
n
)
,
where μ
(
k
) are the roots of eigenvalues and v
i
g
(
k
) are the eigenvectors.The constantsb
j
(
k
) in the homogeneous solutions can be determined from the traction-free boundary conditions on the top and bottom surfaces (σ
i
=
0
at
η
=
±
1
/
2
,
i
=
3,4
,
5), interlaminar stress continuity conditions at the ply interfaces (σ
i
+
(
k
)
=
σ
i
-
(
k
+
1
)
,
i
=
3,4
,
5), and boundary terms Λ set to be zero. These conditions lead to the following equations:
(21)
g
i
(
k
)
=
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
±
1
/
2
+
g
i
(
k
)
(
P
)
=
0
,
(
k
=
1
,
N
)
,
g
i
(
k
)
I
=
b
j
(
k
)
μ
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
±
1
/
2
=
0
,
(
k
=
1
,
N
)
,
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
η
k
+
+
g
i
(
k
)
(
P
)
=
b
j
(
k
+
1
)
v
i
g
(
k
+
1
)
e
μ
(
k
+
1
)
η
|
η
=
η
k
+
1
-
hhhhhhhhhhhhhhhhhhhhhhh
+
g
i
(
k
+
1
)
(
P
)
,
(
k
=
1,2
,
…
,
N
-
1
)
,
b
j
(
k
)
μ
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
η
k
+
=
b
j
(
k
+
1
)
μ
(
k
+
1
)
v
i
g
(
k
+
1
)
e
μ
(
k
+
1
)
η
|
η
=
η
k
+
1
-
,
hhhhhhhhhhhhhhhhhhhhhh
(
k
=
1,2
,
…
,
N
-
1
)
.
Thus, after obtaining the solutions of governing equation, all stress components can be obtained and are listed as follows:
(22)
σ
2
=
∑
i
=
1
n
∑
j
=
1
4
n
f
i
(
ξ
)
[
μ
j
(
k
)
2
v
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
3
=
∑
i
=
1
n
∑
j
=
1
4
n
f
i
II
(
ξ
)
[
v
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
+
g
i
(
k
)
(
P
)
]
,
σ
4
=
∑
i
=
1
n
∑
j
=
1
4
n
-
f
i
I
(
ξ
)
[
μ
j
(
k
)
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
5
=
∑
i
=
1
n
∑
j
=
1
4
n
-
p
i
I
(
ξ
)
[
μ
j
(
k
)
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
6
=
∑
i
=
1
n
∑
j
=
1
4
n
p
i
(
ξ
)
[
μ
j
(
k
)
2
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
.After conducting the second step, if the results are still not convergent, more processes can be conducted. The iterative algorithm can be concluded as follows.(1)
Assume an initial set of out-of-plane stress functions (g
i
(
η
)) and then take the principle of complementary virtual work to solve the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)).
(2)
Use the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)) obtained from the first step as the known functions and repeat taking the principle of complementary virtual work to solve the out-of-plane stress functions (g
i
(
η
)).
(3)
Use the out-of-plane stress functions (g
i
(
η
)) obtained from the second step as the known functions and repeat taking the principle of complementary virtual work to solve the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)).Three steps are narrated in the above algorithm and more processes can be conducted, if required. However, one essential difference should be emphasized. In the first step, the initially assumed out-of-plane stress functions represent the global state of stress distributions which implies that the first step is a stress based equivalent single layer theory, while from the second step, the out-of-plane stress functions are assumed layerwisely and continuity conditions should be enforced at the layer interfaces. All stress components satisfy the boundary conditions automatically when using the extended Kantorovich method.


## 2.1. Derivation of Governing Equation

For the given geometry of laminated composite plate, as shown in Figure1, two piezoelectric actuators are symmetrically surface bonded on the top and bottom surfaces of the laminates without considering the adhesive layer affect. This kind of structure is usually considered to be smart composite laminates [38]. Two piezoelectric layers are used to control the free edge peeling stresses that are generated by mechanical loadings. The longitudinal dimension of the laminates is considered long enough so that the stress fields are independent of the x-axis. Therefore, generalized plane strain states are assumed to be in the x direction.Figure 1
Geometry of laminated composite plate with surface bonded piezoelectric actuators under combined electromechanical loadings.The extended Kantorovich method requires iterative calculation methodology. By preassumption of trial functions and solving a set of ordinary differential equations, the solution of the first process can be obtained. Due to the rough initial assumption, it cannot guarantee the accuracy of solution. Thus, to overcome this weak point, iterative processes can be conducted to improve the accuracy of solution. The governing equations can be derived by using either Galerkin’s method or the principle of virtual work. In this work, the principle of complementary virtual work is employed. The complementary strain energy of the laminated composite plate is calculated by the following equation:(1)
δ
U
=
∬
s
ε
i
δ
σ
i
d
A
i
=
∬
s
ε
2
δ
σ
2
+
ε
3
δ
σ
3
+
ε
4
δ
σ
4
+
ε
5
δ
σ
5
+
ε
6
δ
σ
6
d
A
=
0
.Before substituting strain and stress components into the complementary strain energy equation, an initially assumed stress field should be introduced. Instead of using displacement assumption, stress field assumption is preferred. A good choice for stress field is Lekhnitskii stress functions [39] which can express the three-dimensional stress state of piezo-bonded composite laminates. Lekhnitskii stress functions can automatically satisfy the pointwise equilibrium equations as follows:
(2)
σ
2
=
∂
2
F
∂
η
2
,
σ
3
=
∂
2
F
∂
ξ
2
,
σ
4
=
-
∂
2
F
∂
ξ
∂
η
σ
5
=
-
∂
ψ
∂
ξ
,
σ
6
=
∂
ψ
∂
η
,
where η
=
z
/
H and ξ
=
y
/
H are nondimensionalized coordinates and H is the total thickness of the laminates. The solutions F and ψ are series combinations and assumed to be a combination of in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)) and out-of-plane stress functions (g
i
(
η
)) using separation of variables as follows:
(3)
F
=
∑
i
=
1
n
f
i
(
ξ
)
g
i
(
η
)
,
ψ
=
∑
i
=
1
n
p
i
(
ξ
)
g
i
I
(
η
)
.With the initial assumption of stress field, the strain tensor{
ε
i
} can be found from the linear elastic constitutive relation. This relation is assumed for each individual layer with piezo induced strain term:
(4)
{
ε
i
}
(
k
)
=
[
S
-
]
(
k
)
{
σ
i
}
(
k
)
+
[
d
]
(
k
)
{
E
}
(
k
)
,
(
i
=
1,2
,
…
,
6
)
,
where [
S
-
] is the generalized compliance matrix for orthotropic materials. [
d
]  is the piezoelectric constant matrix and [
E
] is the applied electric field.As the mechanical extension load is applied in thex direction, all other strain components can be expressed using ε
1 and applied piezoelectric fields. Consider
(5)
ε
i
=
S
^
i
j
σ
j
+
S
i
1
S
11
ε
1
+
d
^
i
3
E
3
,
(
i
=
2,3
,
…
,
6
,
j
=
2,3
,
…
,
6
)
,
where
(6)
S
^
i
j
=
S
i
j
-
S
i
1
S
1
j
S
11
,
d
^
i
3
=
d
i
3
-
S
i
1
S
11
d
13
.
Substituting (2) and (5) into (1), the complementary strain energy becomes
(7)
δ
U
=
∬
s
ε
i
δ
σ
i
d
A
i
=
∬
s
[
(
S
^
i
j
σ
j
+
S
i
1
S
11
ε
1
+
d
^
3
i
E
3
)
δ
σ
i
]
d
ξ
d
η
=
0
,
mmmmmmmmmmmmmmmm
(
i
,
j
=
2,3
,
…
,
6
)
.To simplify the above equation, integration by parts and free edge boundary conditions (σ
2
=
σ
4
=
σ
6
=
0, at ξ
=
0) are applied. Then, the double integration is simplified to the following equation and the governing equations can be obtained by using the variational principle. Consider
(8)
∫
[
a
i
j
(
4
)
f
j
IV
+
a
i
j
(
2
)
f
j
II
+
a
i
j
(
0
)
f
j
+
b
i
j
(
2
)
p
j
II
+
b
i
j
(
0
)
p
j
+
r
i
]
δ
f
i
d
ξ
+
∫
[
d
i
j
(
2
)
f
j
II
+
d
i
j
(
0
)
f
j
+
c
i
j
(
2
)
p
j
II
+
c
i
j
(
0
)
p
j
+
s
i
]
δ
p
i
d
ξ
=
0
,
mmmmmmmmiimmmmmmm
(
i
,
j
=
1,2
,
…
,
n
)
,
where
(9)
a
i
j
(
4
)
=
∫
S
^
33
g
i
g
j
d
η
,
a
i
j
(
2
)
=
∫
S
^
23
(
g
i
II
g
j
+
g
i
g
j
II
)
d
η
-
∫
S
^
44
g
i
I
g
j
I
d
η
,
a
i
j
(
0
)
=
∫
S
^
22
g
i
II
g
j
II
d
η
,
b
i
j
(
2
)
=
∫
S
^
36
g
i
g
j
II
d
η
-
∫
S
^
45
g
i
I
g
j
I
d
η
,
b
i
j
(
0
)
=
∫
S
^
26
g
i
II
g
j
II
d
η
,
c
i
j
(
2
)
=
-
∫
S
^
55
g
i
I
g
j
I
d
η
,
c
i
j
(
0
)
=
∫
S
^
66
g
i
II
g
j
II
d
η
,
d
i
j
(
2
)
=
∫
S
^
36
g
i
II
g
j
d
η
-
∫
S
^
45
g
i
I
g
j
I
d
η
,
d
i
j
(
0
)
=
∫
S
^
26
g
i
II
g
j
II
d
η
,
r
i
=
∫
(
S
21
S
11
ε
1
+
d
^
23
E
3
)
g
i
II
d
η
,
s
i
=
∫
(
S
61
S
11
ε
1
+
d
^
63
E
3
)
g
i
II
d
η
.


## 2.2. Solutions of Governing Equations and Generalized Eigenvalue Problem

To solve the fourth-order and second-order coupled ordinary differential equations, the trial functionsg
i
(
η
)  are specified first in (10). Although prescribed boundary conditions are not necessarily satisfied for trial functions when using the extended Kantorovich method, the mode shapes of the clamped-clamped beam are used as the initial out-of-plane trail functions to satisfy the traction-free boundary conditions automatically (σ
3
=
σ
4
=
σ
5
=
0
at
η
=
±
1
/
2). Consider
(10)
g
i
(
η
)
=
cos
⁡
(
β
i
η
)
+
cos
⁡
(
β
i
/
2
)
cosh
⁡
(
β
i
/
2
)
cosh
⁡
(
β
i
η
)
,
where β
i are the roots of the following characteristic equation. Consider
(11)
tan
(
β
2
)
+
tanh
(
β
2
)
=
0
.Thus, all coefficients in the governing equations can be obtained and the general solutions of governing equations can be assumed as the sum of homogeneous solutions(
f
i
(
H
)
and
p
i
(
H
)
) and particular solutions (
f
i
(
P
)
and
p
i
(
P
)
). The homogeneous solutions of f
i
(
H
) and p
i
(
H
) are assumed to be the following form:
(12)
f
i
(
H
)
=
v
i
f
e
λ
ξ
,
p
i
(
H
)
=
v
i
p
e
λ
ξ
.Substituting (12) into the governing equations and defining v
j
f
II
=
λ
2
v
j
f, the governing equations can be rewritten in the matrix form and a standard eigenvalue problem is formed. Consider
(13)
[
0
-
a
i
j
(
4
)
-
b
i
j
(
2
)
0
0
c
i
j
(
0
)
I
0
0
]
-
1
[
a
i
j
(
0
)
a
i
j
(
2
)
b
i
j
(
0
)
d
i
j
(
0
)
d
i
j
(
2
)
b
i
j
(
0
)
0
I
0
]
[
v
j
f
v
j
f
II
v
j
p
]
=
λ
2
[
v
j
f
v
j
f
II
v
j
p
]
,
(
i
,
j
=
1,2
,
…
,
n
)
,
where λ
2 are eigenvalues of the above equation. Since the interlaminar stresses decay in the interior region of piezo-bonded composite laminates, only negative roots of λ
2 are selected. Equation (13) results in 3n eigenvalues and eigenvectors, so that the homogeneous solutions consist of 3n terms linear combination as follows:
(14)
f
i
(
H
)
=
t
j
v
i
j
f
e
-
λ
j
ξ
,
p
i
(
H
)
=
t
j
v
i
j
p
e
-
λ
j
ξ
,
(
i
=
1,2
,
…
,
n
)
,
(
j
=
1,2
,
…
,
3
n
)
,
where t
j  are constants that will be determined by the prescribed boundary conditions.The particular solutions are obtained by consideringf
i
(
P
) and p
i
(
P
) as constants; thus all their derivatives are zeros in the governing equations and the remaining parts can be rewritten in the following equation:
(15)
[
f
j
(
P
)
p
j
(
P
)
]
=
-
[
a
i
j
(
0
)
b
i
j
(
0
)
b
i
j
(
0
)
c
i
j
(
0
)
]
-
1
[
r
i
s
i
]
(
i
,
j
=
1,2
,
…
,
n
)
.The final step of solving the governing equations is to obtain the constantst
j in the homogeneous solutions. They can be calculated by enforcing the free edge boundary conditions (σ
2
=
σ
4
=
σ
6
=
0
at
ξ
=
0), and this will lead f
i,  f
i
I, and p
i to be zeros at the free edges. Consider
(16)
t
j
v
i
j
f
e
-
λ
j
ξ
|
ξ
=
0
+
f
i
(
P
)
=
0
,
-
t
j
λ
j
v
i
j
f
e
-
λ
j
ξ
|
ξ
=
0
=
0
,
t
j
v
i
j
p
e
-
λ
j
ξ
|
ξ
=
0
+
p
i
(
P
)
=
0
.After solving the governing equations, all stress components can be obtained by substituting the solutions into the Lekhnitskii stress functions and they have the following expressions:(17)
σ
2
=
∑
i
=
1
n
∑
j
=
1
3
n
(
t
j
v
i
j
f
e
-
λ
j
ξ
+
f
i
(
P
)
)
g
i
II
,
σ
3
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
i
2
t
j
v
i
j
f
e
-
λ
j
ξ
g
i
,
σ
4
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
j
t
j
v
i
j
f
e
-
λ
j
ξ
g
i
I
,
σ
5
=
∑
i
=
1
n
∑
j
=
1
3
n
λ
j
t
j
v
i
j
p
e
-
λ
j
ξ
g
i
I
,
σ
6
=
∑
i
=
1
n
∑
j
=
1
3
n
(
t
j
v
i
j
p
e
-
λ
j
ξ
+
p
i
(
P
)
)
g
i
II
.


## 2.3. Iterative Procedures

The accuracy of the first step largely depends on the initial assumption of trail functions and the number of series terms involved in the trail functions. The results may also be influenced by the choice of trail functions which may cause singularity problem, harmonic distributions of stresses, and so forth. Thus, it is important to overcome the weaknesses appeared in the first process. We prefer using the extended Kantorovich method which is an iterative process. Through the iterative processes, the influence of the initial assumption can be eliminated and the convergent results can be obtained.In the second step, the solutions of the first step (f
i
(
ξ
) and p
i
(
ξ
)) are considered to be the known functions and the principle of complementary virtual work is used again to derive governing equations. Here, the unknown out-of-plane stress functions (g
i
(
η
)) are assumed to be layer dependent and the continuity conditions are enforced at the layer interface. Thus, from the principle of complementary virtual work in (7) and integration by parts, the governing equation for the out-of-plane stress functions can be obtained as follows:
(18)
∫
[
m
i
j
(
4
)
(
k
)
g
j
(
k
)
(
IV
)
+
m
i
j
(
2
)
(
k
)
g
j
(
k
)
(
II
)
+
m
i
j
(
0
)
(
k
)
g
j
(
k
)
+
x
i
(
k
)
]
×
δ
g
j
(
k
)
d
η
+
Λ
=
0
(
i
,
j
=
1,2
,
…
,
n
)
,
mmmmmmmmmmm
(
k
=
1,2
,
…
,
N
)
,
where Λ represents the boundary terms induced by integration by parts and the superscript (k) refers to the kth ply. The layerwise coefficients are listed as follows:
(19)
m
i
j
(
4
)
(
k
)
=
∫
[
S
^
22
(
k
)
f
i
f
j
+
S
^
26
(
k
)
(
f
i
p
j
+
p
i
f
j
)
+
S
^
66
(
k
)
p
i
p
j
]
d
ξ
,
m
i
j
(
2
)
(
k
)
=
∫
[
S
^
23
(
k
)
(
f
i
f
j
I
I
+
f
i
I
I
f
j
)
+
S
^
36
(
k
)
(
f
i
I
I
p
j
+
p
i
f
j
I
I
)
m
m
m
m
m
i
-
S
^
44
(
k
)
f
i
I
f
j
I
-
S
^
45
(
k
)
(
f
i
I
p
j
I
+
p
i
I
f
j
I
)
-
S
^
55
(
k
)
p
i
I
p
j
I
]
d
ξ
,
m
i
j
(
0
)
(
k
)
=
∫
S
^
33
(
k
)
f
i
I
I
f
j
I
I
d
ξ
,
x
i
(
k
)
=
∫
(
S
13
(
k
)
S
11
(
k
)
ε
1
+
d
^
33
E
3
)
f
i
I
I
d
ξ
.The solution procedure is the same as in the first step. A generalized eigenvalue problem is solved and the general solution is calculated by adding the homogeneous solutionsg
i
(
k
)
(
H
) and the particular solutions g
i
(
k
)
(
P
):
(20)
g
i
(
k
)
=
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
+
g
i
(
k
)
(
P
)
,
(
k
=
1,2
,
…
,
N
)
,
(
i
=
1,2
,
…
,
n
)
,
where μ
(
k
) are the roots of eigenvalues and v
i
g
(
k
) are the eigenvectors.The constantsb
j
(
k
) in the homogeneous solutions can be determined from the traction-free boundary conditions on the top and bottom surfaces (σ
i
=
0
at
η
=
±
1
/
2
,
i
=
3,4
,
5), interlaminar stress continuity conditions at the ply interfaces (σ
i
+
(
k
)
=
σ
i
-
(
k
+
1
)
,
i
=
3,4
,
5), and boundary terms Λ set to be zero. These conditions lead to the following equations:
(21)
g
i
(
k
)
=
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
±
1
/
2
+
g
i
(
k
)
(
P
)
=
0
,
(
k
=
1
,
N
)
,
g
i
(
k
)
I
=
b
j
(
k
)
μ
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
±
1
/
2
=
0
,
(
k
=
1
,
N
)
,
b
j
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
η
k
+
+
g
i
(
k
)
(
P
)
=
b
j
(
k
+
1
)
v
i
g
(
k
+
1
)
e
μ
(
k
+
1
)
η
|
η
=
η
k
+
1
-
hhhhhhhhhhhhhhhhhhhhhhh
+
g
i
(
k
+
1
)
(
P
)
,
(
k
=
1,2
,
…
,
N
-
1
)
,
b
j
(
k
)
μ
(
k
)
v
i
g
(
k
)
e
μ
(
k
)
η
|
η
=
η
k
+
=
b
j
(
k
+
1
)
μ
(
k
+
1
)
v
i
g
(
k
+
1
)
e
μ
(
k
+
1
)
η
|
η
=
η
k
+
1
-
,
hhhhhhhhhhhhhhhhhhhhhh
(
k
=
1,2
,
…
,
N
-
1
)
.
Thus, after obtaining the solutions of governing equation, all stress components can be obtained and are listed as follows:
(22)
σ
2
=
∑
i
=
1
n
∑
j
=
1
4
n
f
i
(
ξ
)
[
μ
j
(
k
)
2
v
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
3
=
∑
i
=
1
n
∑
j
=
1
4
n
f
i
II
(
ξ
)
[
v
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
+
g
i
(
k
)
(
P
)
]
,
σ
4
=
∑
i
=
1
n
∑
j
=
1
4
n
-
f
i
I
(
ξ
)
[
μ
j
(
k
)
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
5
=
∑
i
=
1
n
∑
j
=
1
4
n
-
p
i
I
(
ξ
)
[
μ
j
(
k
)
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
,
σ
6
=
∑
i
=
1
n
∑
j
=
1
4
n
p
i
(
ξ
)
[
μ
j
(
k
)
2
ν
i
j
g
(
k
)
b
j
(
k
)
e
μ
j
(
k
)
η
]
.After conducting the second step, if the results are still not convergent, more processes can be conducted. The iterative algorithm can be concluded as follows.(1)
Assume an initial set of out-of-plane stress functions (g
i
(
η
)) and then take the principle of complementary virtual work to solve the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)).
(2)
Use the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)) obtained from the first step as the known functions and repeat taking the principle of complementary virtual work to solve the out-of-plane stress functions (g
i
(
η
)).
(3)
Use the out-of-plane stress functions (g
i
(
η
)) obtained from the second step as the known functions and repeat taking the principle of complementary virtual work to solve the in-plane stress functions (f
i
(
ξ
) and p
i
(
ξ
)).Three steps are narrated in the above algorithm and more processes can be conducted, if required. However, one essential difference should be emphasized. In the first step, the initially assumed out-of-plane stress functions represent the global state of stress distributions which implies that the first step is a stress based equivalent single layer theory, while from the second step, the out-of-plane stress functions are assumed layerwisely and continuity conditions should be enforced at the layer interfaces. All stress components satisfy the boundary conditions automatically when using the extended Kantorovich method.


## 3. Numerical Results

Based on the theoretical formulations discussed in the previous section, numerical examples are investigated. The laminated composite structures are subjected to uniaxial extension and electric load. The material properties of graphite/epoxy composite laminates [40] and piezoelectric actuator (PZT-5H) are given in Table 1. The lamina thickness is 0.125 mm, which is almost one-fourth of the PZT thickness (0.5 mm). For convenience, 4 laminae, which have the same fiber orientations, are considered one lamina, so that each lamina thickness is then the same as the PZT thickness. The width of the laminate b is equal to four times the total thickness H.Table 1
Material properties of composite lamina and PZT-5H actuator.
Composite
Actuator
E
1
=
138 GPa
E
1
=
E
2
=
62 GPa
E
2
=
E
3
=
14.5 GPa
E
3
=
48 GPa
G
12
=
G
13
=
G
23
=
5.9 GPa
G
12
=
23.5 GPa, G
13
=
G
23
=
23.5 GPa
ν
12
=
ν
23
=
ν
21
=
0.21
ν
12
=
ν
23
=
ν
21
=
0.29
d
31
=
d
32
=
274
×
10
-
12 m/V
d
33
=
-
593
×
10
-
12 m/V, d
15
=
-
741
×
10
-
12 m/V


### 3.1. Convergence Study

The laminates are subjected to uniaxial tension only in this part. To investigate the convergence of the proposed method, different terms of initially assumed stress functions are used and several iterative processes are conducted. To validate the proposed modeling, 3D finite element analysis is also conducted using commercial software ANSYS. The Solid 64 element is used for composite laminates and Solid 5 is used for PZTs. To obtain accurate interlaminar stress distributions, fine mesh is required and it should be refined at the layer interface. Thus, we used 20 elements along the thickness direction for a single layer. The total degree of freedom is 640,281 that is much larger than our proposed method. The three-dimensional solid laminate has a relatively long length, compared with the width and height, to express the generalized plane strain state effectively.Figure2 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge of [PZT/04/904]s laminate under 0.1% uniaxial strain obtained from the first process. The result is compared with the FEM result, and the influence of term number used in the first process is investigated. It is found that the accuracy of the stress distribution largely depends on the number of terms in the first process. With 8 terms of initially assumed stress functions, the result cannot converge well and undesired oscillation appears due to the initial assumption of harmonic functions. Moreover, the stress concentration at the layer interfaces also cannot be predicted accurately. As we expected the peak values should be exactly located at the PZT/0 layer interface (z
/
H
=
1
/
3) and 0/90 layer interface (z
/
H
=
1
/
6), while it cannot be recovered by using 8 terms. Only by increasing the number of terms, the locations can be predicted more precisely. Until using 40 terms of initially assumed stress functions, the stress concentrations can be found exactly located at the layer interfaces, but undesired oscillations still exist. The proposed result can be verified by the FEM result. Both of them show the same tendency and magnitude, whereas the result obtained by finite element method shows difficulties in predicting the stress concentrations at the layer interfaces even using fine mesh. A proper choice of the term number not only guarantees the accuracy of results, but also improves the computational ability. It is to be noted that the computational time using 40 terms for solving the first process only costs seconds, while FEM takes 10–20 minutes to run a single case.Figure 2
σ
3 distribution of [PZT/04/904]s laminate at the free edge, using different terms, under 0.1% uniaxial strain.The results of iterative process for [PZT/04/904]s laminates under 0.1% uniaxial strain are given in Figures 3–5. These results are obtained by only 8 terms and different iterative processes. Figure 3 shows the interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. From the graph, it is found that the iterative processes can significantly improve the accuracy of the stress prediction. The locations of stress concentration can be exactly predicted at the layer interfaces and the oscillations appeared in the first iteration disappeared since the third iteration. Figures 4 and 5 show interlaminar normal stress distribution (σ
3) and interlaminar shear stress distribution (σ
4) through the in-plane direction, respectively. First iteration with 8 terms of initially assumed stress functions cannot guarantee the convergence of the stresses, while the third and fifth iterations provide almost converged stresses. The stresses obtained by FEM show great coincidence with those obtained by the proposed method but finite element analysis underestimates stress concentrations at the free edge. Interlaminar shear stress (σ
4) also cannot satisfy the free edge boundary condition exactly in finite element analysis. It is concluded that the iterative processes can provide the converged stresses with small number of terms. Not only the influence of initial assumption can be eliminated, but also the stress concentration locations can be predicted accurately by the iterative processes. In the following work, the results of the fifth iteration with 8 terms of initially assumed that the stress functions are considered to investigate the reduction of peeling stresses by piezoelectric excitation.Figure 3
σ
3 distribution of [PZT/04/904]s laminate at the free edge from different iterations, under 0.1% uniaxial strain.Figure 4
σ
3 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6 from different iterations, under 0.1% uniaxial strain.Figure 5
σ
4 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6 from different iterations, under 0.1% uniaxial strain.


### 3.2. Reduction of Peeling Stresses Using Piezoelectric Excitations

The existence of peeling stresses generated by the mechanical loadings can be reduced by applying proper electric fields to the actuators. Due to the electromechanical coupling phenomenon, the induced strain or stress generated by the electric loadings has the possibility of compensating the strain or stress generated by the mechanical loadings, in turn, to reduce the peeling stresses and to expand the service life of laminated structures. To investigate the feasibility of peeling stress reduction using piezoelectric actuators, two general lay-ups are considered and numerical examples are discussed in this part.The first example is given by cross ply composite laminate ([PZT/04/904]s) which is subjected to 0.1% uniaxial strain. Moreover, different electric voltages are applied to the actuators that are expected to reduce the free edge peeling stresses. Figure 6 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. It can be clearly found that the peeling magnitude of σ
3 is significantly reduced at the 0/90 layer interface by applying a positive electric voltage to the actuators. Without electric potential, the concentrated peeling stress is 8.3 MPa. It is reduced to 4.9 MPa, when 300 V electric potential is applied. This is of great meaningfulness to the reduction of positive peeling stress at the interface. Figures 7 and 8 show interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) distribution through the in-plane direction at the 0/90 interface, where z
/
H
=
1
/
6. Our ambition is to reduce the peeling stress at the composite layer interface by using piezoelectric actuator. Thus, this interface is what we concerned. σ
3 and σ
4 show large positive values at the free edge and near the free edge under the mechanical loading. With the increasing of the electric potential, the magnitudes of peeling stresses are reduced dramatically. The peeling stress σ
3 can be reduced to almost 50.2 percent when 300 V electric potential is applied which is summarized in Table 2. For interlaminar shear stress (σ
4), 17.3 percent decrement is observed when 300 V electric potential is applied.Table 2
Comparison of peak stresses for cross ply composite laminate ([PZT/04/904]s) at the 0/90 layer interface, where z
/
H
=
1
/
6, under 0.1% uniaxial strain, fifth iteration, n
=
8.
Applied electricpotential (V)
σ
3
σ
4
(MPa)
Reduction (%)
(MPa)
Reduction (%)
0
6.62
0
2.37
0
100
5.51
16.8%
2.19
7.6%
200
4.41
33.4%
2.06
13.1%
300
3.30
50.2%
1.96
17.3%Figure 6
σ
3 distribution of [PZT/04/904]s laminate at the free edge, under 0.1% uniaxial strain.Figure 7
σ
3 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6, under 0.1% uniaxial strain.Figure 8
σ
4 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6, under 0.1% uniaxial strain.The second example is given by angle ply composite laminate ([PZT/454/−454]s) which is subjected to −0.1% uniaxial strain load. The laminate is subjected to compressive load while it sustains large peeling stresses at the 45/−45 layer interfaces. Figure 9 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. It is found that the peeling stress is severely concentrated at the 45/−45 layer interface which could lead to the onset of delamination of composite layers. However, this peeling stress can also be significantly reduced in angle ply laminate, when a positive electric potential is applied to the actuators. The peak value is reduced from 11.83 MPa to 8.52 MPa. The interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) distribution through the in-plane direction at the 45/−45 layer interface are given in Figures 10 and 11, respectively. The interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) well match with those obtained by the finite element analysis for the mechanical loading case. However, finite element analysis still underestimates the stress concentration at the free edge. The detailed reductions of peeling stresses can be found in Table 3. Since angle ply composite laminates sustain larger shear stresses than cross ply composite laminates, the reduction of σ
3 is smaller than that of cross ply, while the reduction of σ
4 is larger than that of the cross ply. The proposed reduction algorithm also works for the in-plane stress (σ
2) which is obtained at the interior of the laminate, shown in Figure 12. It is interesting to find that the magnitude of σ
2 also decreases when a positive electric potential is applied. This result can be verified by using classical laminate theory (CLT). The fifth result can accurately recover the CLT solution. Thus, the proposed reduction approach is well suited to the angle ply composite laminate.Table 3
Comparison of peak stresses for angle ply laminate ([PZT/454/−454]s) at the 45/−45 layer interface, where z
/
H
=
1
/
6, under −0.1% uniaxial strain, fifth iteration, n
=
8.
Applied electricpotential (V)
σ
3
σ
4
(MPa)
Reduction (%)
(MPa)
Reduction (%)
0
11.28
0
1.81
0
100
10.73
9.3%
1.63
10.0%
200
9.63
18.6%
1.46
19.3%
300
8.52
24.5%
1.30
28.2%Figure 9
σ
3 distribution of [PZT/454/−454]s laminate at the free edge, under −0.1% uniaxial strain.Figure 10
σ
3 distribution of [PZT/454/−454]s laminate at the interface of z
/
H
=
1
/
6, under −0.1% uniaxial strain.Figure 11
σ
4 distribution of [PZT/454/−454]s laminate at the interface of z
/
H
=
1
/
6, under −0.1% uniaxial strain.Figure 12
σ
2 distribution of [PZT/454/−454]s laminate at the interior, under −0.1% uniaxial strain.


## 3.1. Convergence Study

The laminates are subjected to uniaxial tension only in this part. To investigate the convergence of the proposed method, different terms of initially assumed stress functions are used and several iterative processes are conducted. To validate the proposed modeling, 3D finite element analysis is also conducted using commercial software ANSYS. The Solid 64 element is used for composite laminates and Solid 5 is used for PZTs. To obtain accurate interlaminar stress distributions, fine mesh is required and it should be refined at the layer interface. Thus, we used 20 elements along the thickness direction for a single layer. The total degree of freedom is 640,281 that is much larger than our proposed method. The three-dimensional solid laminate has a relatively long length, compared with the width and height, to express the generalized plane strain state effectively.Figure2 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge of [PZT/04/904]s laminate under 0.1% uniaxial strain obtained from the first process. The result is compared with the FEM result, and the influence of term number used in the first process is investigated. It is found that the accuracy of the stress distribution largely depends on the number of terms in the first process. With 8 terms of initially assumed stress functions, the result cannot converge well and undesired oscillation appears due to the initial assumption of harmonic functions. Moreover, the stress concentration at the layer interfaces also cannot be predicted accurately. As we expected the peak values should be exactly located at the PZT/0 layer interface (z
/
H
=
1
/
3) and 0/90 layer interface (z
/
H
=
1
/
6), while it cannot be recovered by using 8 terms. Only by increasing the number of terms, the locations can be predicted more precisely. Until using 40 terms of initially assumed stress functions, the stress concentrations can be found exactly located at the layer interfaces, but undesired oscillations still exist. The proposed result can be verified by the FEM result. Both of them show the same tendency and magnitude, whereas the result obtained by finite element method shows difficulties in predicting the stress concentrations at the layer interfaces even using fine mesh. A proper choice of the term number not only guarantees the accuracy of results, but also improves the computational ability. It is to be noted that the computational time using 40 terms for solving the first process only costs seconds, while FEM takes 10–20 minutes to run a single case.Figure 2
σ
3 distribution of [PZT/04/904]s laminate at the free edge, using different terms, under 0.1% uniaxial strain.The results of iterative process for [PZT/04/904]s laminates under 0.1% uniaxial strain are given in Figures 3–5. These results are obtained by only 8 terms and different iterative processes. Figure 3 shows the interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. From the graph, it is found that the iterative processes can significantly improve the accuracy of the stress prediction. The locations of stress concentration can be exactly predicted at the layer interfaces and the oscillations appeared in the first iteration disappeared since the third iteration. Figures 4 and 5 show interlaminar normal stress distribution (σ
3) and interlaminar shear stress distribution (σ
4) through the in-plane direction, respectively. First iteration with 8 terms of initially assumed stress functions cannot guarantee the convergence of the stresses, while the third and fifth iterations provide almost converged stresses. The stresses obtained by FEM show great coincidence with those obtained by the proposed method but finite element analysis underestimates stress concentrations at the free edge. Interlaminar shear stress (σ
4) also cannot satisfy the free edge boundary condition exactly in finite element analysis. It is concluded that the iterative processes can provide the converged stresses with small number of terms. Not only the influence of initial assumption can be eliminated, but also the stress concentration locations can be predicted accurately by the iterative processes. In the following work, the results of the fifth iteration with 8 terms of initially assumed that the stress functions are considered to investigate the reduction of peeling stresses by piezoelectric excitation.Figure 3
σ
3 distribution of [PZT/04/904]s laminate at the free edge from different iterations, under 0.1% uniaxial strain.Figure 4
σ
3 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6 from different iterations, under 0.1% uniaxial strain.Figure 5
σ
4 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6 from different iterations, under 0.1% uniaxial strain.


## 3.2. Reduction of Peeling Stresses Using Piezoelectric Excitations

The existence of peeling stresses generated by the mechanical loadings can be reduced by applying proper electric fields to the actuators. Due to the electromechanical coupling phenomenon, the induced strain or stress generated by the electric loadings has the possibility of compensating the strain or stress generated by the mechanical loadings, in turn, to reduce the peeling stresses and to expand the service life of laminated structures. To investigate the feasibility of peeling stress reduction using piezoelectric actuators, two general lay-ups are considered and numerical examples are discussed in this part.The first example is given by cross ply composite laminate ([PZT/04/904]s) which is subjected to 0.1% uniaxial strain. Moreover, different electric voltages are applied to the actuators that are expected to reduce the free edge peeling stresses. Figure 6 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. It can be clearly found that the peeling magnitude of σ
3 is significantly reduced at the 0/90 layer interface by applying a positive electric voltage to the actuators. Without electric potential, the concentrated peeling stress is 8.3 MPa. It is reduced to 4.9 MPa, when 300 V electric potential is applied. This is of great meaningfulness to the reduction of positive peeling stress at the interface. Figures 7 and 8 show interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) distribution through the in-plane direction at the 0/90 interface, where z
/
H
=
1
/
6. Our ambition is to reduce the peeling stress at the composite layer interface by using piezoelectric actuator. Thus, this interface is what we concerned. σ
3 and σ
4 show large positive values at the free edge and near the free edge under the mechanical loading. With the increasing of the electric potential, the magnitudes of peeling stresses are reduced dramatically. The peeling stress σ
3 can be reduced to almost 50.2 percent when 300 V electric potential is applied which is summarized in Table 2. For interlaminar shear stress (σ
4), 17.3 percent decrement is observed when 300 V electric potential is applied.Table 2
Comparison of peak stresses for cross ply composite laminate ([PZT/04/904]s) at the 0/90 layer interface, where z
/
H
=
1
/
6, under 0.1% uniaxial strain, fifth iteration, n
=
8.
Applied electricpotential (V)
σ
3
σ
4
(MPa)
Reduction (%)
(MPa)
Reduction (%)
0
6.62
0
2.37
0
100
5.51
16.8%
2.19
7.6%
200
4.41
33.4%
2.06
13.1%
300
3.30
50.2%
1.96
17.3%Figure 6
σ
3 distribution of [PZT/04/904]s laminate at the free edge, under 0.1% uniaxial strain.Figure 7
σ
3 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6, under 0.1% uniaxial strain.Figure 8
σ
4 distribution of [PZT/04/904]s laminate at the interface of z
/
H
=
1
/
6, under 0.1% uniaxial strain.The second example is given by angle ply composite laminate ([PZT/454/−454]s) which is subjected to −0.1% uniaxial strain load. The laminate is subjected to compressive load while it sustains large peeling stresses at the 45/−45 layer interfaces. Figure 9 shows interlaminar normal stress (σ
3) distribution through the thickness direction at the free edge. It is found that the peeling stress is severely concentrated at the 45/−45 layer interface which could lead to the onset of delamination of composite layers. However, this peeling stress can also be significantly reduced in angle ply laminate, when a positive electric potential is applied to the actuators. The peak value is reduced from 11.83 MPa to 8.52 MPa. The interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) distribution through the in-plane direction at the 45/−45 layer interface are given in Figures 10 and 11, respectively. The interlaminar normal stress (σ
3) and interlaminar shear stress (σ
4) well match with those obtained by the finite element analysis for the mechanical loading case. However, finite element analysis still underestimates the stress concentration at the free edge. The detailed reductions of peeling stresses can be found in Table 3. Since angle ply composite laminates sustain larger shear stresses than cross ply composite laminates, the reduction of σ
3 is smaller than that of cross ply, while the reduction of σ
4 is larger than that of the cross ply. The proposed reduction algorithm also works for the in-plane stress (σ
2) which is obtained at the interior of the laminate, shown in Figure 12. It is interesting to find that the magnitude of σ
2 also decreases when a positive electric potential is applied. This result can be verified by using classical laminate theory (CLT). The fifth result can accurately recover the CLT solution. Thus, the proposed reduction approach is well suited to the angle ply composite laminate.Table 3
Comparison of peak stresses for angle ply laminate ([PZT/454/−454]s) at the 45/−45 layer interface, where z
/
H
=
1
/
6, under −0.1% uniaxial strain, fifth iteration, n
=
8.
Applied electricpotential (V)
σ
3
σ
4
(MPa)
Reduction (%)
(MPa)
Reduction (%)
0
11.28
0
1.81
0
100
10.73
9.3%
1.63
10.0%
200
9.63
18.6%
1.46
19.3%
300
8.52
24.5%
1.30
28.2%Figure 9
σ
3 distribution of [PZT/454/−454]s laminate at the free edge, under −0.1% uniaxial strain.Figure 10
σ
3 distribution of [PZT/454/−454]s laminate at the interface of z
/
H
=
1
/
6, under −0.1% uniaxial strain.Figure 11
σ
4 distribution of [PZT/454/−454]s laminate at the interface of z
/
H
=
1
/
6, under −0.1% uniaxial strain.Figure 12
σ
2 distribution of [PZT/454/−454]s laminate at the interior, under −0.1% uniaxial strain.


## 4. Conclusions

In this study, the reduction of free edge peeling stresses of laminated composite structures using piezoelectric actuators was analyzed by using the extended Kantorovich method. Multiterms of the initially assumed stress functions were involved, and three-dimensional stress distributions were obtained. The accuracy and effectiveness of the present approach were demonstrated by comparing the stresses obtained by ANSYS. The provided results showed excellent agreements with those obtained by FEM, while it showed accuracy and great efficiency in computation. It was found that small terms of trail functions could lead to accurate solutions. The present method is general and applicable to various lay-up configurations. The peeling stresses generated by mechanical load at the free edges were significantly reduced by applying a proper electric potential. The present method could serve as an efficient tool for reducing the free edge peeling stresses and extending the service life of laminated composite structures.



---

*Source: 439492-2014-06-16.xml*