Series Solution-Based Approach for the Interlaminar Stress Analysis of Smart Composites under Thermo-Electro-Mechanical Loading

Abstract

This paper introduces a new loading condition considering the combined thermo-electro-mechanical coupling effect in a series solution-based approach to analyze the free-edge interlaminar stresses in smart composite laminates. The governing equations are developed using the principle of complementary virtual work. The assumed stress fields satisfy the traction-free and free-edge boundary conditions. The accurate stress states of the composite structures are acquired through the procedure of generalized eigenvalue problems. The uniform temperature is employed throughout the laminate, and the electric field loading is applied to the symmetric piezo-bonded actuators to examine the combined effect of thermal and electrical stresses on the overall deformation of smart composite laminates. It was observed that the magnitude of the peeling stresses generated by mechanical loading was reduced by the combined thermal and electric excitation loading (up to 25.3%), which in turn resulted in expanding the service life of the smart composite structures. The proposed approach is implemented on three different layup configurations. The efficiency of the current methodology is confirmed by comparing the results with the 3D finite element (FEM) solution computed by ABAQUS.

1. Introduction

Recently, smart composite laminates are being used in numerous fields due to their potential applications in energy harvesting [1], vibration suppression [2], and modern control engineering [3]. Compared to other metals and nonmetals, fiber-reinforced laminated composites offer many physical advantages, such as high strength, light weight, and high stiffness. These composite laminates can be suitable candidates for substrate materials in smart structures. Lately, different finite element modeling techniques have been developed to analyze composite structures [4,5]. Sayed et al. [6] presented a novel approach for modeling the dust-loaded electrostatic precipitators by developing a finite difference method integrated with a full multi-grid method to resolve Poisson and continuity equations on one fine computational domain. Still, the accurate prediction of delamination is a significant issue due to interlaminar stresses in the laminated composites [7]. The high free-edge stresses between the composite laminates and the smart elements (sensors/actuators) ultimately result in the element’s debonding failure from the composite structure.

The free-edge effect is a phenomenon that occurs due to the existence of materials and geometric discontinuities at the free edges. Due to the discontinuities, stresses are locally concentrated to achieve a stress equilibrium state at the interface of each lamina [8]. The stress concentration near the free edge is a type of boundary layer solution that is usually discriminated from an interior solution. The boundary layer solution (i.e., self-equilibrated stress field) decays rapidly at the inward part of the solids according to Saint Venant’s principle [9]. In composite laminates, the decaying distance is much longer than the conventional isotropic solid, since the material anisotropy and geometrical data determine the distance. The local substantial stress concentrations cause delamination, which is the first sign of the failure of composite laminates. After the delamination, continued stress concentration leads to the propagation of cracks, diminishing the overall strength and stiffness of the composite laminates. Thus, when analyzing or designing mechanical and aerospace structures or structural components using composite laminates, the free-edge effect should never be overlooked.

After Pipes and Pagano’s [10] pioneering work, numerous approximation methods have been developed. Comprehensive reviews and assessments related to the free-edge effect can be found in the literature [8,11,12]. Theoretical approaches of the interlaminar stress analysis are divided into displacement-based and stress-based formulations. Because of the convenience of deriving the finite element formulation, displacement-based techniques are widely studied. The theoretical models are primarily developed using equivalent single-layer [13,14,15] or layerwise models [16,17]. The stress-based approaches [18,19,20] have also been popularly developed. Compared to the displacement-based approaches, these approaches can precisely fulfill the traction-free boundary condition.

Yin et al. [19,20] adopted piecewise polynomial approximation using the layerwise theory (LWT) to predict interlaminar stress distribution. However, compared to the equivalent single-layer theory (ESLT), the LWT has many degrees of freedom (DOF). Therefore, from the computational aspect, it is more efficient to use ESLT over LWT. Flanagan et al. [21] proposed an efficient ESLT approach based on the series expansion of eigenmode shapes. Although Flanagan obtained accurate in-plane stresses, the out-of-plane stresses showed oscillation behavior. The extended Kantorovich method [22,23] was proposed to eliminate the oscillation behavior of interlaminar stresses, and the converged stresses were obtained. Kim and co-workers applied this methodology to the strength analyses of composite laminates [24], calculation of the stress field of composite laminates containing an internal ply-drop [25], and the stress analysis of the composite patch for repairing [26]. Layup optimization of composite laminates was carried out under different loading conditions by combining the methodology with a genetic algorithm [27,28,29]. Nosier et al. [30] used the first-order shear deformation theory of plates and Reddy’s layerwise theory to analyze the edge effect in the laminates subjected to axial strain and developed a reduced form of the displacement field for long antisymmetric angle-ply composite laminates subjected to extensional and/or torsional loads. Similarly, Lo et al. [31] developed a higher-order theory to analyze free-edge problems in symmetric and unsymmetric laminates under extension, bending, and thermal loading.

Meanwhile, various delamination prevention methods have been proposed. These methods include varying the assembling sequence of layers and changing the fiber edge shape. However, these techniques could not achieve the expected designed goals. Kim et al. [32,33] employed thermal loading through the thickness direction to effectively control the interlaminar stresses produced by the extension load. The proposed approach significantly reduced interlaminar stresses. Further, Huang et al. [34] successfully employed the appropriate electric fields to piezoelectric actuators to minimize the interlaminar stresses produced by mechanical loading. However, to the author’s knowledge, no research has yet been conducted on analyzing the free-edge stresses under the combined thermo-electro-mechanical loading by using the series solution-based approach.

The stress concentration occurs at the free edge of smart composites due to the actuating force of the piezoelectric actuators. The debonding and delamination failure of the piezoelectric actuator initiates due to the stress concentration. The previous studies were focused on the global response of smart laminates, such as health monitoring [35] and energy harvesting [36]. Huang et al. [37] proposed the stress-function–based approach to analyze the stress concentration in smart laminates by considering the electro-mechanical coupling effect. The accurate interlaminar stresses were obtained, and the stress-function-based approach was proven to be efficient and accurate in predicting the interlaminar stresses in smart composites. This study introduces a stress-function-based approach to analyze the free edge interlaminar stresses in smart composite laminates under the combined thermo-electro-mechanical loading. The governing equations are developed using the principle of complementary virtual work. The assumed stress fields satisfy the traction-free and free-edge boundary conditions. The accurate stress states of the composite structures are acquired through the procedure of generalized eigenvalue problems. The proposed approach is implemented on three different layup configurations. The efficiency of the current methodology is confirmed by comparing the results with the 3D finite element (FEM) solution computed by ABAQUS.

2. Theoretical Formulation

A schematic of the piezo-bonded smart composite laminate is shown in Figure 1. It is composed of piezoelectric (PZT) actuators attached on the upper and lower surfaces of the laminates. The laminate width is represented by ‘2b’, which is four times the thickness of ‘h’. First, the tensile load is applied on the laminate’s free edges (0, 2b) along the x-direction, generating significant interlaminar stresses. These stresses are reduced by using the combined electric and thermal fields. The constitutive equations [37,38] used for the modeling of the smart composite laminated structure are discussed in this section.

The constitutive equations of the kth layer with the applied thermal and electric loading are described as follows:

$$\left\{{\epsilon }_i\right\}={\left[d\right]}^{\left(k\right)}{\left\{E\right\}}^{\left(k\right)}+{\left[\overline{S}\right]}^{\left(k\right)}\left\{{\sigma }_i\right\}+{\left[\alpha \right]}^{\left(k\right)}{\left\{\Delta T\right\}}^{\left(k\right)}$$

(1)

where $\left[d\right]$ represents the piezoelectric constant; $\left\{E\right\}$ and $\left\{\Delta T\right\}$ are the applied electric field and thermal loading, respectively; $\left\{{\sigma }_i\right\}$ and $\left\{{\epsilon }_i\right\}$ denote the stress and strain vectors, respectively; and $\left[\overline{S}\right]$ represents the global compliance matrix for the orthotropic materials. “i” indicates the number of plies from 1 to 6, and “k” represents the kth ply in the laminate.

To achieve the static equilibrium conditions, the assumed stress field is made to be independent of the x direction. The normal axial stress from Equation (1) is computed as follows:

$${\sigma }_1^{\left(k\right)}=\frac{{\epsilon }_1^{\left(k\right)}-S_{1j}^{\left(k\right)}{\sigma }_j^{\left(k\right)}-d_{13}{}^{\left(k\right)}{E}_3{}^{\left(k\right)}-{\alpha }_1^{\left(k\right)}\Delta T^{\left(k\right)}}{S_{11}^{\left(k\right)}},\hspace{1em}\left(j=2,3,6\right)$$

(2)

The strain components are obtained by substituting Equation (3) into Equation (1), as follows:

$${\epsilon }_i^{\left(k\right)}={\widehat{S}}_{ij}^{\left(k\right)}{\sigma }_j^{\left(k\right)}+\frac{S_{1i}^{\left(k\right)}}{S_{11}^{\left(k\right)}}{\epsilon }_1^{\left(k\right)}+{\widehat{\alpha }}_i^{\left(k\right)}\Delta T^{\left(k\right)}+{\widehat{d}}_{13}^{\left(k\right)}{E}_3{}^{\left(k\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i,j=2,3,\cdots ,6\right)$$

(3)

where

$${\widehat{S}}_{ij}^{\left(k\right)}=S_{ij}^{\left(k\right)}-\frac{S_{1i}^{\left(k\right)}{S}_{1j}^{\left(k\right)}}{S_{11}^{\left(k\right)}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\widehat{\alpha }}_i^{\left(k\right)}={\alpha }_i^{\left(k\right)}-\frac{{\alpha }_1^{\left(k\right)}{S}_{1i}^{\left(k\right)}}{S_{11}^{\left(k\right)}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\widehat{d}}_{i3}^{\left(k\right)}=d_{i3}^{\left(k\right)}-\frac{d_{13}^{\left(k\right)}{S}_{1i}^{\left(k\right)}}{S_{11}^{\left(k\right)}}$$

(4)

Under the uniaxial extension load, the displacement fields presenting the generalized plane strain state can be expressed as:

$$\begin{array}{c}u\hspace{0.17em}\left(x,y,z\right)=Ax+U\left(y,z\right)+{\omega }_{1}z-{\omega }_{2}y+u_0\\ v\hspace{0.17em}\left(x,y,z\right)=V\left(y,z\right)+{\omega }_{2}x-{\omega }_{3}z+v_0\\ w\hspace{0.17em}\left(x,y,z\right)=W\left(y,z\right)+{\omega }_{3}y-{\omega }_{1}x+w_0\end{array}$$

(5)

where $A$ represents the extension of material along x-axis. $U_{y,z}$, $V_{y,z}$, and $W_{y,z}$ denote the deformations in YZ plane. The constants (${\omega }_1$, ${\omega }_2$, ${\omega }_3$, $u_0$, $v_0$, and $w_0$) represent the displacements of the rigid body.

$$\xi =\frac{y}{h} , \eta =\frac{z}{h}$$

(6)

where ‘’h’’ denotes the thickness of the laminate. Equation (6) nondimensionalizes the coordinates. For this study, pointwise equilibrium equations are satisfied using Lekhnitskii stress functions [39]. Stress functions are separated into “ξ”-dependent and “η”-dependent functions, respectively. The separation is carried out using the variables separation method:

$$\frac{{\partial }^{2}F}{\partial {\eta }^2}={\sigma }_2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{{\partial }^{2}F}{\partial {\xi }^2}={\sigma }_3,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{{\partial }^{2}F}{\partial \eta \partial \xi }=-{\sigma }_4,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial G}{\partial \xi }=-{\sigma }_5,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial G}{\partial \eta }={\sigma }_6$$

(7)

where

$$F={\displaystyle \sum }_{i=1}^n{f}_i\left(\xi \right)g_i\left(\eta \right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}G={\displaystyle \sum }_{i=1}^n{p}_i\left(\xi \right)g_{i,\eta }\left(\eta \right)$$

(8)

where $n$ is the number of stress functions which represents the eigenmodes of the present system. The subscript $i$ denotes the differentiation of the function with respect to $\eta$. In the previous studies, Kim et al. and Lee et al. [19,20] have already shown the convergence according to n. In this paper, the authors follow the previous works. Here, total complementary strain energy expression is used to get the governing equations:

$$\begin{array}{c}\delta U={{\displaystyle \iint }}^{ }\epsilon \delta \sigma dydz\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={{\displaystyle \iint }}^{ }\left[{\sigma }_i{\widehat{S}}_{ij}\delta {\sigma }_j+\left({\epsilon }_1\frac{S_{i1}}{S_{11}}+{\widehat{\alpha }}_i\Delta T\right)\delta {\sigma }_j\right]d\xi d\eta =0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i,j=2,3,\cdots ,6\right)\end{array}$$

(9)

Replacing Equations (7) and (8) into Equation (9), and repeating the integration by part, gives:

$$\begin{array}{c}{{\displaystyle \int }}^{ }\left[a_{ij}^{\left(4\right)}{f}_{j,\xi \xi \xi \xi }+a_{ij}^{\left(2\right)}{f}_{j,\xi \xi }+a_{ij}^{\left(0\right)}{f}_j+b_{ij}^{\left(2\right)}{p}_{j,\xi \xi }+b_{ij}^{\left(0\right)}{p}_{j,\xi \xi }+r_i\right]\delta f_{i}d\xi \\ +{{\displaystyle \int }}^{ }\left[c_{ij}^{\left(2\right)}{p}_{j,\xi \xi }+c_{ij}^{\left(0\right)}{p}_j+d_{ij}^{\left(2\right)}{f}_j^{II}+d_{ij}^{\left(0\right)}{f}_j+s_i\right]\delta p_{i}d\xi =0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i,j=1,2,\cdots ,n\right)\end{array}$$

(10)

where

$$\begin{array}{c}{a}_{ij}^{\left(4\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{33}{g}_i{g}_{j}d\eta \\ a_{ij}^{\left(2\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{23}\left(g_{i,\eta \eta }{g}_j+g_i{g}_{j,\eta \eta }\right)d\eta -{{\displaystyle \int }}^{ }{\widehat{S}}_{44}{g}_{i,\eta }{g}_{j,\eta }d\eta \\ a_{ij}^{\left(0\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{22}{g}_{i,\eta \eta }{g}_{j,\eta \eta }d\eta \\ b_{ij}^{\left(2\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{36}{g}_i{g}_{j,\eta \eta }d\eta -{{\displaystyle \int }}^{ }{\widehat{S}}_{45}{g}_{i,\eta }{g}_{j,\eta }d\eta \\ b_{ij}^{\left(0\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{26}{g}_{i,\eta \eta }{g}_{j,\eta \eta }d\eta \\ c_{ij}^{\left(2\right)}=-{{\displaystyle \int }}^{ }{\widehat{S}}_{55}{g}_{i,\eta }{g}_{j,\eta }d\eta \\ c_{ij}^{\left(0\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{66}{g}_{i,\eta \eta }{g}_{j,\eta \eta }d\eta \\ d_{ij}^{\left(2\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{36}{g}_{i,\eta \eta }{g}_{j}d\eta -{{\displaystyle \int }}^{ }{\widehat{S}}_{45}{g}_{i,\eta }{g}_{j,\eta }d\eta \\ d_{ij}^{\left(0\right)}={{\displaystyle \int }}^{ }{\widehat{S}}_{26}{g}_{i,\eta \eta }{g}_{j,\eta \eta }d\eta \\ r_i={{\displaystyle \int }}^{ }\left(\frac{S_{12}{}^{\left(k\right)}}{S_{11}{}^{\left(k\right)}}{\epsilon }_1+{\widehat{\alpha }}_2{}^{\left(k\right)}\Delta T+{\widehat{d}}_{23}{}^{\left(k\right)}{E}_3{}^{\left(k\right)}\right)g_{i,\eta \eta }d\eta \\ s_i={{\displaystyle \int }}^{ }\left(\frac{S_{16}{}^{\left(k\right)}}{S_{11}{}^{\left(k\right)}}{\epsilon }_1+{\widehat{\alpha }}_6{}^{\left(k\right)}\Delta T+{\widehat{d}}_{63}{}^{\left(k\right)}{E}_3{}^{\left(k\right)}\right)g_{i,\eta \eta }d\eta \end{array}$$

(11)

In Equation (10), f(ξ) and p(ξ), both the in-plane stress functions, are unknown; however, the solution begins with the calculation of “g(η)”as the out-of-plane stress function, so the generalized configuration of “g(η)” is expressed as follows:

$$g_i\left(\eta \right)=\cos \mathrm{h}{k}_i{\eta }^{*}-\cos k_i{\eta }^{*}\hspace{0.17em}-K_i\left(\sin \mathrm{h}{k}_i{\eta }^{*}-\sin k_i{\eta }^{*}\right)$$

(12)

where

$$K_i=\frac{\cos \mathrm{h}{k}_{i}h-\cos k_{i}h}{\sin \mathrm{h}{k}_{i}h-\sin k_{i}h},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\eta }^{*}=\eta +\frac{1}{2}$$

(13)

After determining the coefficients of Equation (10), the governing equations are solved by assuming exponential functions:

$$f_i(\xi )={\upsilon }_i^f{e}^{\lambda \xi },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_i(\xi )={\upsilon }_i^p{e}^{\lambda \xi }$$

(14)

Substituting Equation (14) into Equation (10) gives:

$$\begin{array}{l}{a}_{ij}^{(0)}{\upsilon }_j^f+\left(a_{ij}^{(2)}+{\lambda }^2{a}_{ij}^{(4)}\right){\upsilon }_{j,\xi \xi }^f+\left(b_{ij}^{(0)}+{\lambda }^2{b}_{ij}^{(2)}\right){\upsilon }_j^p=0\\ d_{ij}^{(0)}{\upsilon }_j^f+d_{ij}^{(2)}{\upsilon }_{j,\xi \xi }^f+\left(c_{ij}^{(0)}+{\lambda }^2{c}_{ij}^{(2)}\right){\upsilon }_j^p=0\\ {\lambda }^2{\upsilon }_j^f-{\upsilon }_{j,\xi \xi }^f=0\end{array}$$

(15)

The third part of Equation (15), an auxiliary equation, transforms the original problem to the generalized eigenvalue form:

$$\left[\begin{array}{ccc}{a}_{ij}^{\left(0\right)}& a_{ij}^{\left(2\right)}& b_{ij}^{\left(0\right)}\\ d_{ij}^{\left(0\right)}& d_{ij}^{\left(2\right)}& c_{ij}^{\left(0\right)}\\ 0& -I& 0\end{array}\right]\left[\begin{array}{c}{\upsilon }_j^f\\ {\upsilon }_{j,\xi \xi }^f\\ {\upsilon }_j^p\end{array}\right]={\lambda }^2\left[\begin{array}{ccc}0& -a_{ij}^{\left(4\right)}& -b_{ij}^{\left(2\right)}\\ 0& 0& -c_{ij}^{\left(0\right)}\\ -I& 0& 0\end{array}\right]\left[\begin{array}{c}{\upsilon }_j^f\\ {\upsilon }_{j,\xi \xi }^f\\ {\upsilon }_j^p\end{array}\right]$$

(16)

Considering the characteristics of interlaminar stress, the negative roots of ${\lambda }^2$ are chosen in Equation (16). The solution of this problem gives 3n eigen values, and a linear combination of the 3n-term also determines the homogeneous solutions of in-plane stress functions.

$$\begin{array}{l}{f}_i^{(H)}={\upsilon }_{ij}^f{t}_j{e}^{-{\lambda }_j\xi }\\ p_i^{(H)}={\upsilon }_{ij}^p{t}_j{e}^{-{\lambda }_j\xi },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i=1,\hspace{0.17em}2,\hspace{0.17em}\hspace{0.17em}\cdots ,\hspace{0.17em}n,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\hspace{0.17em}2,\hspace{0.17em}\cdots ,\hspace{0.17em}3n\right)\end{array}$$

(17)

The particular solution can be obtained when f(ξ) and p(ξ) are constants in Equation (11).

$$\begin{array}{l}{a}_{ij}^{(0)}{f}_j^{(P)}+b_{ij}^{(0)}{p}_j^{(P)}=-r_i\\ d_{ij}^{(0)}{f}_j^{(P)}+c_{ij}^{(0)}{p}_j^{(P)}=-s_i\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i,j=1,2,\cdots ,n\right)\end{array}$$

(18)

To compute the constants, the free-edge boundary conditions of σ₄ and σ₆ are employed, which leads to:

$$\begin{array}{l}{\upsilon }_{ij}^f{t}_j=-f_i^{(P)}\\ {\lambda }_j{\upsilon }_{ij}^f{t}_j=0\\ {\upsilon }_{ij}^p{t}_j=-p_i^{(P)},\hspace{1em}\hspace{1em}\left(i=1,2,\cdots ,n,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,\cdots ,3n\right)\end{array}$$

(19)

By replacing the obtained in-plane stress functions into Equations (8) and (9), the final interlaminar stresses can be determined as follows:

$$\begin{array}{l}{\sigma }_2=\left({\upsilon }_{ij}^f{t}_j{e}^{-{\lambda }_j\xi }+f_j^{(P)}\right)g_{i,\eta \eta }\\ {\sigma }_3={\lambda }_j^2{\upsilon }_{ij}^f{t}_j{e}^{-{\lambda }_j\xi }{g}_i\\ {\sigma }_4={\lambda }_j{\upsilon }_{ij}^f{t}_j{e}^{-{\lambda }_j\xi }{g}_{i,\eta }\\ {\sigma }_5={\lambda }_j^2{\upsilon }_{ij}^p{t}_j{e}^{-{\lambda }_j\xi }{g}_i\\ {\sigma }_6=\left({\upsilon }_{ij}^p{t}_j{e}^{-{\lambda }_j\xi }+p_j^{(P)}\right)g_{i,\eta }\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(i=1,2,\cdots ,n,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,\cdots ,3n\right)\end{array}$$

(20)

3. Computational Results

This section presents the computational results of the proposed methodology.

Table 1. Material properties of the PZT-5H layer and composite lamina.

PZT-5H Layer [40]	Composite Lamina [14]
E1 = E2 = 62 × 103 MPa	E1 = 138 × 103 MPa
E3 = 48 × 103 MPa	E2 = E3 = 14.5 × 103 MPa
G12 = 23.5 × 103 MPa	G12 = G23 = G13 = 5.9 × 103 MPa
G13 = G23 = 5.9 × 103 MPa	ν12 = ν13 = ν23 = 0.21
ν12 = ν13 = ν23 = 0.49	α1 = 1.2 × 10−7 (0C−1)
d31 = d32 = −2.74 × 10−10 m/N	α2 = α3 = 8.4 × 10−6 (0C−1)
d33 = 5.93 × 10−10 m/N
d42 = d51 = 7.41 × 10−10 m/N
α1 = α2 = α3 = 3.5 × 10−6 (0C−1)

shows the material properties of the PZT-5H actuator [40] and the graphite/epoxy composite [9] used:

Although the stress-function-based approach for interlaminar stress analysis is effective, an FE-based solution is computed in ABAQUS to validate the present analysis results. The thickness of each lamina is 0.0125 cm and PZT layer is four times the thickness of the composite lamina. Four laminae with the same fiber orientation are assumed as one lamina to match the thickness of the PZT layer. The obtained interlaminar stresses are compared with the 3D finite element solutions computed from the ABAQUS software [41]. Three-dimensional ‘HEX’ (C3D8) is chosen in ABAQUS, with a total of 449,280 DOF. The length of the solid laminate follows the plane strain condition, i.e., the length is ten times larger than the width of the laminate (L = 10 × 2b = 40H). Axis-symmetric boundary conditions are employed, and the quarter model is used to reduce the computational cost, as shown in Figure 2. The results are extracted at the middle part of the laminate, where the boundary layer effects along the x-axis decayed. To reduce the computational time, the element size is refined at the middle region where the results are extracted, while the relatively coarse mesh was applied in the other region. Global and local coordinates were defined to describe the laminate and material behavior correctly.

This study’s primary purpose is to show how to reduce the interlaminar stresses generated by extension load by applying thermal and electrical loading. Therefore, the phenomenon of interlaminar stress reduction in smart composite laminates is demonstrated by combined thermo-mechanical and thermo-electro-mechanical loadings.

3.1. Uniaxial Extension and Thermal Loading

This section consists of piezo-bonded composite laminates under combined 0.1% uniaxial tension and unit temperature rise thermal loading (ΔT = 1). Three layup configurations are analyzed. A clamped boundary condition is adopted at one end along the x-axis, while at the other, a 0.1% longitudinal tensile strain load is applied. Unit temperature rise is employed throughout the composite laminates. Figure 3a,b show the free-edge interlaminar stress (σ₃) distribution under tensile and thermal loadings. Stress concentrations occurred at the interfaces of the layers. The result of the present method correlates well with the 3D finite element result, but the FE-based solution underestimates the stress concentrations at the free edge. All the graphs show that the present method’s solutions are greater than the 3D finite element solution, because the 3D solid elements popularly used in FE analysis are based on the displacement field; the solution cannot guarantee the statically admissible stress field, even though the displacement field is exact.

Figure 4a,b signify the normal (σ₃) and shear stresses (σ₄) at the PZT/0 and PZT/45 interfaces of the laminates under 0.1% uniaxial extension and unit thermal loading (ΔT = 1). The maximum normal stress (σ₃) concentration in both cross-ply and angle-ply laminates occurs at the free edge that can cause the PZT layer’s delamination failure. The present analysis results of free-edge of cross-ply and angle-ply indicate that the large shear stress (σ₄) represents the satisfaction of the traction-free boundary condition and becomes null at the interior region. The results of the 3D-FEM solution show good agreement with the present analysis results.

3.2. Reduction of Interlaminar Stress by Combined Thermo-Mechanical Loading

Thermal loading has the capability to reduce the interlaminar stresses generated by the extension load. The stresses or strains caused by thermal loading can compensate for the stresses and strains induced by the mechanical loading, due to the thermal-mechanical coupling phenomenon. This, in turn, reduces the peeling stresses and enhances the lifespan of the composite structures. Figure 5 shows an example of free-edge interlaminar stress σ₃ of cross-ply laminate under the combined 0.1% uniaxial extension loading, and the variable rise in the thermal loading throughout the laminate. The interlaminar stress distribution’s magnitude is reduced at the PZT/0 layer and 0/90 layer interface by increasing the laminate’s temperature. Without the thermal loading, the magnitude of peeling stress is 8.3 MPa; however, as the temperature loading of ΔT = 30 is applied, peeling stress is reduced to 2.9 MPa.

Similarly, Figure 6a,b show the (σ₃) and (σ₄) distribution through the in-plane direction at the PZT/0 interface. The magnitude of the peeling stresses decreased significantly with the increase in thermal loading:

3.3. Piezoelectric Excitation

This section shows interlaminar stress analysis results under the electric field loading. 2 × 10⁵ V/m of electric excitation is applied on both surfaces of the piezoelectric actuators. Figure 7 shows the free edge σ₃ distribution. The maximum stress concentration occurs at the PZT/0 and PZT/45 layer interfaces. The FEM results show a good correlation with the present analysis results.

Figure 8 shows the normal stress (σ₃) and the shear stresses (σ₄) at the PZT/0 and PZT/45 interfaces of the studied laminates. For all laminates the maximum normal stress (σ₃) concentration appears at the free edge that can cause the PZT layer’s delamination failure. On the other hand, FEM results indicate that the traction-free boundary conditions are not satisfied as the FEM follows the displacement-based approach. Therefore, the present analysis results provide more accurate results as compared to the FEM-based results.

3.4. Reduction in Interlaminar Stress by Combined Thermo-Electro-Mechanical Loading

This section provides the stress distributions at the free edge and the interfaces of the laminas under the combined thermo-electro-mechanical loading. The induced stresses and strains generated by mechanical loading can be reduced by the applied thermal and electric excitation loading due to the combined coupling effect. Thus, the extent of free-edge interlaminar stresses can be lower down, which results in expanding the service life of the smart composite structures. Figure 9a shows the interlaminar stress σ₃ of the free edge of the cross-ply laminate under the combined thermo-electro-mechanical loading with 0.1% uniaxial tension, temperature loading of ΔT = 10, and the variable electric field of (E3 = 0 V/m, E3 = 1 × 10⁵ V/m, and E3 = 2 × 10⁵ V/m). It can be observed that, without applying the electric field, the concentration of peeling stress is 6.3 MPa, and upon applying the excitation of 2 × 10⁵ V/m, it is reduced to 4.7 MPa (up to 25.3% reduction). Similarly, at the PZT/0 layer in the in-plane direction, the peak value of the concentrated peeling stress drops with increased electric field loading.

In the case of angle-ply laminate, Figure 10a shows the interlaminar stress σ₃ at the free edge under the combined thermo-electro-mechanical loading with 0.1% uniaxial tension, temperature loading of ΔT = 10, and the variable electric field of (E3 = 0 V/m, E3 = 1 × 10⁵ V/m, and E3 = 2 × 10⁵ V/m). Large stress concentrations can be seen to occur on the free edge of the 45/−45 layer interface. When an electric excitation at the piezoelectric layer is applied, the peeling stress in angle-ply laminate could be significantly reduced.

In the literature, different approaches were introduced to analyze the free-edge interlaminar stresses.

Table 2. Comparative analysis of the proposed approach with the previous literature approaches.

Approach	Contribution	Limitation
Layerwise theory [19]	Prediction of interlaminar stress in a multilayered strip of a laminate subjected to extension, bending, and twisting loading	Computationally inefficient, many degrees of freedom
Equivalent single-layer theory [21]	Accurate prediction of inplane stresses	Out-of-plane stresses showed oscillation behavior
Series solution based approach [37]	Minimized interlaminar stresses by electro-mechanical coupling phenomenon	Missing thermo-electro-mechanical coupling effect
Proposed approach	Accurate prediction of interlaminar stresses under thermo-electro-mechanical loading	Only applicable for the simple geometry

summarizes the literature’s approaches and their limitations compared to the currently proposed methodology. The proposed work solved the limitations of the previous studies and provided an accurate interlaminar stress distribution analysis under thermo-electro-mechanical loading. The drawback of the proposed approach is that it is only applicable for simple geometries.

4. Conclusions

This study provides a simple and efficient approach to estimate the free-edge interlaminar stress distribution in smart composites under the combined thermo-electro-mechanical loading. Accurate interlaminar stresses are obtained by a series solution-based approach. The assumed stress fields satisfy the traction-free and free-edge boundary conditions. The proposed approach is implemented on three different composite layups. The obtained results are compared with the 3D-FEM solutions computed by ABAQUS. The present methodology results provide good agreement with the 3D-FEM results and deliver more accurate results with superior computational efficiency. The proposed approach is independent of the number-stacking layers, making it a more efficient approach than the displacement-based approach. The interlaminar stresses produced by the mechanical load are substantially decreased at the free edge and the layer interfaces by applying the appropriate thermal and electric loading due to an efficient thermo-electro-mechanical coupling effect. It was observed that the magnitude of the peeling stresses generated by mechanical loading was reduced by the combined thermal and electric excitation loading (up to 25.3%), which in turn resulted in expanding the service life of the smart composite structures. A limitation of this approach is that only simple geometry can be analyzed; thus, it can be used as a reference approach to finite element analysis and the initial design phase of smart composite structures.