# A Coupling Electromechanical Cell-Based Smoothed Finite Element Method Based on Micromechanics for Dynamic Characteristics of Piezoelectric Composite Materials

Authors: Jianxiao Zheng; Zhishan Duan; Liming Zhou
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

Coupling electromechanical cell-based smoothed finite element method (CSFEM) with the asymptotic homogenization method (AHM) is presented to overcome the overstiffness of FEM. This method could accurately simulate the dynamic responses and electromechanical coupling effects of piezoelectric composite material (PCM) structures. Firstly, the efficient performances for active compounds of round cross-section fibers are calculated based on AHM. Secondly, in the CSFEM, electromechanical multi-physic-field FEM is coupled with gradient smoothing technique. CSFEM returns the nearly exact stiffness of continuum structures, which auto discretes the elements in complex areas more readily and thus remarkably reduces the numerical errors. Static and dynamic characteristics of four PCM structures are investigated using CSFEM with AHM. Results are compared with analytical solution and those of FEM, which proves that CSFEM with AHM is more accurate and reliable than the standard FEM when solving problems of complex structures. Additionally, CSFEM could provide results of higher accuracy even using distorted meshes. Therefore, such method is a robust tool for analyzing mechanical properties of PCM structures.


---

## Body

## 1. Introduction

Piezoelectric composite materials (PCMs) could confirm between mechanical energy and electrical energy [1]. PCMs are made of piezoceramics and passive non-piezoelectric polymers [2, 3]. These composites possess superior properties owing to the most promising characteristics of components as well as various structures.Because of the electromechanical effects, PCMs are more often used as sensors or actuators into noise, vibration, precision position control, energy harvesting, and structural health monitoring [4–6]. Thus, the electric-mechanical properties of PCMs were investigated widely. The heterogeneous media can be generally characterized by micromechanical models. Under such scenario, numerical or analytical ways were used to electromechanically characterize PCMs.Some analytical approaches were presented to investigate PCMs, but the methods were limited by boundary conditions and loading cases [7, 8]. Moreover, single inclusion was applied into a piezoelectric material based on the micromechanical theory with the electromechanical solution [9–13]. Also, the asymptotic homogenization method (AHM) was developed to solve the effective coefficients of PCMs with the square fibers distribution [14, 15]. AHM can be used to calculate effective performances of structures made of hexagonal symmetrical and randomly distributed fibers [16, 17]. The effective elastic coefficients of periodic composites could be solved by an asymptotic method [18]. Otero et al. expressed the effective properties of reinforced PCMs in the closed form [19]. de Medeiros et al. studied the effective coefficients of PCMs made of circular or squared cross-sectional fibers based on AHM [20]. Viaño et al. derived a high-order asymptotic expanding model of piezoelectric rods [21]. Le built an exact 2D theory based on the variation asymptotic method for functionally graded PCMs [22]. Fantoni et al. proposed a multifield AHM to analyze the periodic microstructured PCMs exposed to body force, charge density, and heat sources [23]. However, closed-form solutions which are very a valuable benchmark to issues concerning PCMs are often inaccessible except for relatively simple boundary and geometry conditions.Among numerical approaches, researchers usually use the finite element model (FEM) to develop a specific representative volume element (RVE) for various physical problems, including piezoelectricity. Kar-Gupta and Venkatesh developed an FEM to evaluate how the fiber form affected the general piezoelectric characteristic of PCMs [24]. Jin et al. developed a micromechanical model based on linear stress strain relations and the stress amplification factor to determine the microstress under various mechanical loadings [25]. de Medeiros et al. compared the analytical model and FEM with RVEs in calculating the effective coefficients of PCMs [20]. Würkner et al. established a numerical model to estimate the effective indexes of unidirectional fiber-reinforced PCMs with rhombic fiber arrangement [26]. Mishra et al. extended FEM to evaluate the effective properties of PCMs with SU8 photoresist as the matrix reinforced by the vertically arrayed ZnO nanowires [27]. Bowen et al. discussed the effective properties of new PCMs based on ferroelectric PCR-7M ceramic [28]. Though most of the FEMs can interpolate displacement and electric potential as kinematic field indicators with satisfactory compatibility equations, these methods are often limited by overly stiffness, inaccuracy, and sensitivity to mesh distortion. These limitations can be overcome by hybrid and mixed finite elements [29–32]. Other contributions include drilling degree-of-freedom piezoelectric elements [33–35]. For details and review on FEM establishment for PCM analysis and modeling, refer to Reference [36]. So far, much research has been conducted to develop new special elements [12, 37]. However, the application of FEM is yet limited by the occurrence of mesh distortion.To solve the problem of mesh distortion in FEM, Liu et al. presented a new smoothed FEM by incorporating the gradient smoothing technique (GST) into FEM for solid mechanics problems, based on the nodal integrated meshless methods [38–42]. The FEM employs standard Galerkin formulation in which the stiffness matrix is stiffened by the presumed displacement field [43]. In the smoothed FEM, the smoothed Galerkin weak form is used where GST softens the stiffness matrix through the softening effect. Theoretically, the smoothed FEM in the energy norm often creates a softer stiffness matrix than the FEM with the same background meshes [44, 45]. This unique ability endows smoothed FEM with many critical characteristics [46], such as the easier modeling, upper bound solution property [47], and even nearly perfect solutions in a norm [43, 48–52]. The softening effects enable smoothed FEM to process highly distorted meshes and n-sided polygonal elements [45]. The smoothed FEM interpolates displacement according to the same lower standard mesh and evaluates the weak form based on the smoothing zones.An ultraprecise hybrid smoothed FEM for the piezoelectric problem was designed using with the simplest three- and four-sided elements in 2D and 3D, respectively [53]. An edge-based smoothed FEM was proposed for static eigenvalue assessment of 2D piezoelectric structures [54]. An effective numerical method was presented to optimize and maximize the basic frequency of functionally graded carbon-nanotube-strengthened four-sided sheets [55]. On the commercial software ABAQUS, Bhowmick and Liu established a phase-field model based on cell-based smoothed FEM (CSFEM) to address the brittle fracturing in solids [56]. Pramod et al. developed the CSFEM with the Reissner’s mixed variational theorem to analyze the static and free vibration of cross-ply sheet plates [57].In this work, the dynamic characteristics on PCM structures were studied by using the technique based on the effective CSFEM with AHM. Longitudinal/transversal elastic, piezoelectric, and dielectric effective parameters of a piezoceramic fiber with an O-shaped geometrical section buried in a non-piezoelectric material were computed by AHM based on micromechanics. The CSFEM of PCM structures was presented by applying GST into the exit FEM for an electromechanical coupling field. Then, the equations of the dynamic responses under the multifield for PCM instruments were deduced. Finally, a bilayered PCM actuator and a PCM energy harvester based on micromechanics were calculated by CSFEM. Results of CSFEM were compared with those of FEM, which validated that CSFEM possessed the advantages of accuracy, convergence, and efficiency. Besides, CSFEM was insensitive to mesh distortion, which was very useful to analyze complex structures or large deformation problems. Therefore, CSFEM with AHM performed better than FEM.This paper is organized as follows: Section2 introduces the basic formulations for PCMs. Sections 3 and 4 describe the AHM based on micromechanics and CSFEM, respectively. Section 5 puts forward the modified Wilson-θ method. In Section 6, the bilayered PCM actuator and the PCMs energy harvester were elaborated. Section 7 gives the conclusions.


## 2. Basic Equations for PCMs

### 2.1. Geometry and Coordinate System

The cross section of a PCM beam was a rectangle with lengthL, breadth b, and height h. Here, the continuum theory was used. The PCMs, cylindrical piezoceramic fibers, were buried in the epoxy substrate. The Cartesian coordinate system (x1, x2, x3) and geometric parameters are illustrated in Figure 1.Figure 1
A PCMs beam: geometric parameters and Cartesian coordinate system.


### 2.2. Constitutive Equations

The 3D linear constitutive equations, which can very accurately model the electromechanical coupling behavior of PCMs, are polarized along the global coordinates as follows:(1)σ11σ22σ33σ23σ13σ12=c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66ε11ε22ε332ε232ε132ε12−00e3100e3200e330e240e1500000E1E2E3,(2)D1D2D3=0000e150000e2400e31e32e33000ε11ε22ε332ε232ε132ε12+χ11000χ22000χ33E1E2E3,where σij and εkl are the stress and infinitesimal strain tensors, respectively; Ei and Di are electric field and electric displacement vector components, respectively; eik, ckl, and χij are the piezoelectric, elastic, and dielectric material constants, respectively.In the PCM beam, equations (1) and (2) could be written in the matrix form:(3)σ=Cε−eE,D=eTε+χE,where(4)σ=σ11σ33σ13,C=c11c130c31c33000c55,e=0e310e33e150,E=E1E3,ε=ε11ε33ε13,D=D1D3,χ=χ1100χ33.The balance law and Gauss’ law underlie the description of the coupled field responses, and in the case of no body force and for quasistatic electromechanics, equations can be expressed as(5)σji,j=0,σji=σij,Di,i=0.The boundary conditions are as follows:For the mechanical fields,(6)uit=u¯itonΓu,σijnj=β¯itonΓβ,Γ=Γu∪Γβ.For the electric fields,(7)Φt=Φ¯tonΓp,Dini=Q¯tonΓb,Γ=Γa∪Γb,where Γ, u¯it, and Φ¯t are the global boundary, displacement, and electrical potential, respectively; Γu, Γβ, Γp, and Γq are the prescribed displacement, traction vector, normal component of the electric displacement vector, and electrical potential of Γ, respectively; β¯it and Q¯t are the surface load and surface density of free charge, respectively.


## 2.1. Geometry and Coordinate System

The cross section of a PCM beam was a rectangle with lengthL, breadth b, and height h. Here, the continuum theory was used. The PCMs, cylindrical piezoceramic fibers, were buried in the epoxy substrate. The Cartesian coordinate system (x1, x2, x3) and geometric parameters are illustrated in Figure 1.Figure 1
A PCMs beam: geometric parameters and Cartesian coordinate system.


## 2.2. Constitutive Equations

The 3D linear constitutive equations, which can very accurately model the electromechanical coupling behavior of PCMs, are polarized along the global coordinates as follows:(1)σ11σ22σ33σ23σ13σ12=c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66ε11ε22ε332ε232ε132ε12−00e3100e3200e330e240e1500000E1E2E3,(2)D1D2D3=0000e150000e2400e31e32e33000ε11ε22ε332ε232ε132ε12+χ11000χ22000χ33E1E2E3,where σij and εkl are the stress and infinitesimal strain tensors, respectively; Ei and Di are electric field and electric displacement vector components, respectively; eik, ckl, and χij are the piezoelectric, elastic, and dielectric material constants, respectively.In the PCM beam, equations (1) and (2) could be written in the matrix form:(3)σ=Cε−eE,D=eTε+χE,where(4)σ=σ11σ33σ13,C=c11c130c31c33000c55,e=0e310e33e150,E=E1E3,ε=ε11ε33ε13,D=D1D3,χ=χ1100χ33.The balance law and Gauss’ law underlie the description of the coupled field responses, and in the case of no body force and for quasistatic electromechanics, equations can be expressed as(5)σji,j=0,σji=σij,Di,i=0.The boundary conditions are as follows:For the mechanical fields,(6)uit=u¯itonΓu,σijnj=β¯itonΓβ,Γ=Γu∪Γβ.For the electric fields,(7)Φt=Φ¯tonΓp,Dini=Q¯tonΓb,Γ=Γa∪Γb,where Γ, u¯it, and Φ¯t are the global boundary, displacement, and electrical potential, respectively; Γu, Γβ, Γp, and Γq are the prescribed displacement, traction vector, normal component of the electric displacement vector, and electrical potential of Γ, respectively; β¯it and Q¯t are the surface load and surface density of free charge, respectively.


## 3. AHM

A fiber and matrix two-phase structure (Figure2(a)) is defined by a composite body Ω and the axis x3 parallel to the fibers. Straight reinforcement is periodically arranged. Hence, RVE could be illustrated in Figure 2(b). The periodic cell S has a squared transversal , which includes a circle of radius R (Figure 2(c)).Figure 2
The geometric mode of PCMs: (a) structure, (b) representative volume element (RVE), and (c) transversal cross-section of the RVE.The phases are supposed to be homogeneous and electroelastically linear. Between the fiber and the matrix is the interface Γ, which is regarded as perfect. In the composite, the dimensionsl and L are the fiber diameter and the center-to-center distance between two adjacent cylinders (fibers), respectively. Since ε = l/L is far smaller than the dimension of RVE, two space scales the slow variable x and the fast variable y = x/ε, can be separated. The stress and strain tensors, electrical displacement vector, and electric potential are continuous across Γ or namely between phases. Thus, Ckl, eik, and χij have piecewise constant roles in the periodic cells.The first set of math issues has new physical relationships over Ω with mechanical (C¯), piezoelectric (e¯), and dielectric (χ¯) constant coefficients, which reflect the characteristics of a uniform medium Ω and are called the effective characteristics of PCMs (Figure 2). The coefficients c¯ijpq, e¯ipq, and χ¯ip could be computed using the equations in Reference [14].


## 4. Electromechanical CSFEM

### 4.1. Weak Formulation

The virtual work for a piezoelectric material with volume Ω and regular boundary surface Γ can be expressed as(8)δΠ=δU−δW=−∫Ωδε¯TσdΩ+∫ΓδuTβ¯dΓ+∫ΩρδuTu¨dΩ+∫ΩδE¯TDdΩ+∫ΓQ¯δΦdΓ=0,where ρ and δ are mass density and virtual quantity, respectively. Ω is divided into np elements, which contains Nn nodes, and the approximation displacement u¯ and electrical potential Φ¯ for an FGPM problem can be expressed as(9)u¯=∑i=1npNiuui=Nuu,Φ¯=∑i=1npNiΦΦi=NΦΦ,where Nu and NΦ are the shape functions of electromechanical CSFEM displacement and electrical potential, respectively.A 4-node element is separated into four smoothing subdomainsΩik. The nodes of field and edge/center smoothing, the edge Gaussian points, outer normal vector distribution, and the shape functions are shown in Figure 3.Figure 3
Smoothing subcells and the values of shape functions.At any pointxk in Ωik, the smoothed form of strain ε¯xk and electric field E¯xk are(10)ε¯xk=∫Ωikεxrx−xkdΩ,(11)E¯xk=∫ΩikExrx−xkdΩ,where rx−xk is the constant function:(12)rx−xk=1/Aik,x∈Ωik,0,x∉Ωik,where(13)Aik=∫ΩikdΩ.Substituting equation (12) in equations (10) and (11), we have(14)ε¯xk=1Aik∫ΓiknukudΓ,(15)E¯xk=1Aik∫ΓiknΦkΦdΓ,where Γik is the boundary of Ωik and nuk and nΦk are the outer normal vector matrices of the smoothing domain boundary:(16)nuk=nx1k00nx3knx3knx1k,nΦk=nx1knx3k.Equations (14) and (15) can be newly expressed as(17)ε¯xk=∑i=1neB¯uixkui,E¯xk=−∑i=1neB¯ΦixkΦi,where ne is the number of smoothing elements.(18)B¯uixk=1Aik∫ΓkNiunx1k00Niunx3kNiunx3kNiunx1kdΓ,(19)B¯Φixk=1Aik∫ΓkNiΦnx1kNiΦnx3kdΓ,At the Gaussian pointxbG, equations (18) and (19) become(20)B¯uixk=1Aik∑b=1nbNiuxbGnx1k00NiuxbGnx3kNiuxbGnx3kNiuxbGnx1klbk,B¯Φixk=1Aik∑b=1nbNiΦxbGnx1kNiΦxbGnx3klbk,where lbk is the length of the smoothing boundary and nb is the total number of boundaries in a subdomain. As the shape function varies linearly along each side of the subdomain, one Gauss point is enough for precise boundary integration [38].The essential distinction between CSFEM and FEM is that FEM needs to construct the shape function matrix of the element, while CSFEM only needs to use the function at the Gaussian point of the smoothing element boundary and avoids the function derivatives, which reduces the continuity requirement of the shape function and improves the accuracy and convergence.The dynamic model of the PCM electromechanical system can be deduced from the Hamilton rule as follows:(21)M¯q¨+K¯q=F¯,where(22)M¯=Muu000,q¨=u¨Φ¨,K¯=KuuKuΦKuΦTKΦΦ,q=uΦ,F¯=FQ,Muu=∑eMuue,Muue=diagm1,m1,m2,m2,m3,m3,m4,m4,where mi=piT, Aik (i = 1, 2, 3, 4) is the mass of the i-th smoothing element related to node i, T is the smoothing element thickness, and pi is the density of the Gaussian integration point of the i-th smoothing subdomain.(23)Kuu=∑i=1ncB¯uiTC¯B¯uiAik,KuΦ=∑i=1ncB¯uiTe¯B¯ΦiAik,KΦΦ=−∑i=1ncB¯ΦiTχ¯B¯ΦiAik,F=∫ΓqNuTβdΓ,Q=∫ΩNΦTQ¯dΩ,where nc=np×ne.


## 4.1. Weak Formulation

The virtual work for a piezoelectric material with volume Ω and regular boundary surface Γ can be expressed as(8)δΠ=δU−δW=−∫Ωδε¯TσdΩ+∫ΓδuTβ¯dΓ+∫ΩρδuTu¨dΩ+∫ΩδE¯TDdΩ+∫ΓQ¯δΦdΓ=0,where ρ and δ are mass density and virtual quantity, respectively. Ω is divided into np elements, which contains Nn nodes, and the approximation displacement u¯ and electrical potential Φ¯ for an FGPM problem can be expressed as(9)u¯=∑i=1npNiuui=Nuu,Φ¯=∑i=1npNiΦΦi=NΦΦ,where Nu and NΦ are the shape functions of electromechanical CSFEM displacement and electrical potential, respectively.A 4-node element is separated into four smoothing subdomainsΩik. The nodes of field and edge/center smoothing, the edge Gaussian points, outer normal vector distribution, and the shape functions are shown in Figure 3.Figure 3
Smoothing subcells and the values of shape functions.At any pointxk in Ωik, the smoothed form of strain ε¯xk and electric field E¯xk are(10)ε¯xk=∫Ωikεxrx−xkdΩ,(11)E¯xk=∫ΩikExrx−xkdΩ,where rx−xk is the constant function:(12)rx−xk=1/Aik,x∈Ωik,0,x∉Ωik,where(13)Aik=∫ΩikdΩ.Substituting equation (12) in equations (10) and (11), we have(14)ε¯xk=1Aik∫ΓiknukudΓ,(15)E¯xk=1Aik∫ΓiknΦkΦdΓ,where Γik is the boundary of Ωik and nuk and nΦk are the outer normal vector matrices of the smoothing domain boundary:(16)nuk=nx1k00nx3knx3knx1k,nΦk=nx1knx3k.Equations (14) and (15) can be newly expressed as(17)ε¯xk=∑i=1neB¯uixkui,E¯xk=−∑i=1neB¯ΦixkΦi,where ne is the number of smoothing elements.(18)B¯uixk=1Aik∫ΓkNiunx1k00Niunx3kNiunx3kNiunx1kdΓ,(19)B¯Φixk=1Aik∫ΓkNiΦnx1kNiΦnx3kdΓ,At the Gaussian pointxbG, equations (18) and (19) become(20)B¯uixk=1Aik∑b=1nbNiuxbGnx1k00NiuxbGnx3kNiuxbGnx3kNiuxbGnx1klbk,B¯Φixk=1Aik∑b=1nbNiΦxbGnx1kNiΦxbGnx3klbk,where lbk is the length of the smoothing boundary and nb is the total number of boundaries in a subdomain. As the shape function varies linearly along each side of the subdomain, one Gauss point is enough for precise boundary integration [38].The essential distinction between CSFEM and FEM is that FEM needs to construct the shape function matrix of the element, while CSFEM only needs to use the function at the Gaussian point of the smoothing element boundary and avoids the function derivatives, which reduces the continuity requirement of the shape function and improves the accuracy and convergence.The dynamic model of the PCM electromechanical system can be deduced from the Hamilton rule as follows:(21)M¯q¨+K¯q=F¯,where(22)M¯=Muu000,q¨=u¨Φ¨,K¯=KuuKuΦKuΦTKΦΦ,q=uΦ,F¯=FQ,Muu=∑eMuue,Muue=diagm1,m1,m2,m2,m3,m3,m4,m4,where mi=piT, Aik (i = 1, 2, 3, 4) is the mass of the i-th smoothing element related to node i, T is the smoothing element thickness, and pi is the density of the Gaussian integration point of the i-th smoothing subdomain.(23)Kuu=∑i=1ncB¯uiTC¯B¯uiAik,KuΦ=∑i=1ncB¯uiTe¯B¯ΦiAik,KΦΦ=−∑i=1ncB¯ΦiTχ¯B¯ΦiAik,F=∫ΓqNuTβdΓ,Q=∫ΩNΦTQ¯dΩ,where nc=np×ne.


## 5. Modified Wilson-θ Method

The modified Wilson-θ method is an important scheme and an implicit integral way to solve the dynamic system equations. If θ > 1.37, the solution is unconditionally stable. The detailed procedures are shown as follows.


### 5.1. Initial Calculation

(1)
Formulate generalized stiffness matrixK¯, mass matrix M¯
(2)
Calculate initial values ofu, u˙, u¨
(3)
Select the time step Δt and the integral constant θ (θ = 1.4)(24)a0=6θΔt2,a1=3θΔt,a2=2a1,a3=θΔt2,a4=a0θ,a5=−a2θ,a6=1−3θ,a7=Δt2,a7=Δt2,a8=−Δt26.
(4)
Formulate an effective generalized stiffness matrixK˜: K˜=K¯+a0M¯.


### 5.2. For Each Time Step

(1)
Calculate the payload at timet + θΔt:(25)F¯t+θΔt=F¯t+θF¯t+Δt−F¯t+M¯a0qt+a2q˙t+2q¨t.
(2)
Calculate the generalized displacement at timet + θΔt:(26)K˜ut+θΔt=F¯t+θΔt.
(3)
Calculate the generalized acceleration, generalized speed, and generalized displacement at timet + Δt:(27)u¨t+Δt=a4ut+θΔt−ut+a5u˙t+a6u¨t,u˙t+Δt=u˙t+a7u¨t+Δt+u¨t,ut+Δt=ut+Δtu˙t+a8u¨t+Δt+2u¨t.


## 5.1. Initial Calculation

(1)
Formulate generalized stiffness matrixK¯, mass matrix M¯
(2)
Calculate initial values ofu, u˙, u¨
(3)
Select the time step Δt and the integral constant θ (θ = 1.4)(24)a0=6θΔt2,a1=3θΔt,a2=2a1,a3=θΔt2,a4=a0θ,a5=−a2θ,a6=1−3θ,a7=Δt2,a7=Δt2,a8=−Δt26.
(4)
Formulate an effective generalized stiffness matrixK˜: K˜=K¯+a0M¯.


## 5.2. For Each Time Step

(1)
Calculate the payload at timet + θΔt:(25)F¯t+θΔt=F¯t+θF¯t+Δt−F¯t+M¯a0qt+a2q˙t+2q¨t.
(2)
Calculate the generalized displacement at timet + θΔt:(26)K˜ut+θΔt=F¯t+θΔt.
(3)
Calculate the generalized acceleration, generalized speed, and generalized displacement at timet + Δt:(27)u¨t+Δt=a4ut+θΔt−ut+a5u˙t+a6u¨t,u˙t+Δt=u˙t+a7u¨t+Δt+u¨t,ut+Δt=ut+Δtu˙t+a8u¨t+Δt+2u¨t.


## 6. Numerical Results

### 6.1. A Singer-Player PCMs Strip

In this example to test the precision of CSFEM under mechanical and electrical boundary conditions, we used the shear deformation of a piezoelectric strip (1 × 1 mm2, thickness t = 10 μm) under compressive stress σ0 = 5 N·mm−2 and applied voltage V0 = 1000 V (Figure 4). The PCMs, piezoceramic fibers (PZT-5A), were buried in the non-piezoelectric epoxy matrix. The material properties are shown in Table 1 [20]. A block of 10000 square cells (each cell 10 × 10 × 10 μm3) was used in the composite z with the specific fiber volume fraction of 55.55%. One part of the PCM structure (60 × 10 × 40 μm3) is shown in Figure 5. The materials were polarized under the electric field set to the left and right edges, which resulted in shear strain. The boundary conditions were applied to the strip edges:(28)Φ,zx,z=±h=0,σzzx,z=±h=σ0,σxzx=L,z=0,σxzx,z=±h=0,Φx=L,z=−V0,σxxx=L,z=0,Φx=0,z=V0,ux=0,z=0,vx=0,z=0=0.Figure 4
A singer-player PCMs strip under shear deformation.Table 1
Material constants for constituents of the PZT-5A and epoxy matrix.
Material constants
c
11 (GPa)
c
12 (GPa)
c
13 (GPa)
c
33 (GPa)
c
44 (GPa)
c
66 (GPa)
e
13 (C/m2)
e
15 (C/m2)
e
33 (C/m2)
χ
11 (F/m)
χ
33 (F/m)
ρ (kg/m3)
Epoxy matrix
3.86
2.57
2.57
3.86
0.64
0.64
—
—
—
0.0797
0.0797
7500
PZt-5A
121.0
75.4
75.2
111.0
21.1
2.28
−5.4
12.3
15.8
8.11
7.35
7500Figure 5
A part homogeneous model of the singer-player PCMs strip by AHM.The analytical solution to this problem is shown below [58]:(29)u=s13σ0x,v=d15V0xh+s33σ0z,Φ=V01−2xL.It should be noted that all 11 effective coefficients in Table2 were determined using AHM. CSFEM, FEM, and analytical solution were all based on AHM. Three meshes (each 10 × 10, including one with uniform elements and two with distorted elements) were analyzed by CSFEM (Figure 6). The αir is the irregularity factor that is assigned between 0.0 and 0.5 [38]. Meanwhile, 60 × 60 uniform elements were used by FEM. Figure 7 illustrates the CSFEM, FEM, and analytical horizontal displacement ux at the central line (z = 0), while the vertical displacement uz at z = 0 of the single-layer PCM strip is shown in Figure 8. The distributions of the electric potential Φ at z = 0 with FEM and the analytical solutions are demonstrated in Figure 9. Clearly, the displacements and electric potential computed by CSFEM match well with the analytical solutions and outperform those estimated by FEM for all three meshes, which suggested that CSFEM can recreate the linear behavior of the analytical solutions.Table 2
The material properties of the singer-player PCMs strip were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
0.5555
9.74292
5.54232
5.95557
35.11748
2.17285
2.1003
0.021904
10.879245
−0.251715
0.28018
4.27248Figure 6
Typical irregular elements with irregularity factorαir = (a) 0.3 (I1) and (b) 0.5 (I2).
(a)
(b)Figure 7
Variation of horizontal displacementux at the central line (z = 0) of the singer-player PCMs strip.Figure 8
Variation of vertical displacementuz at the central line (z = 0) of the singer-player PCMs strip.Figure 9
Variation of electric potential Φ at the central line (z = 0) of the singer-player PCMs strip.Moreover, the computation cost of CSFEM over PCMs in a homogenized structure surpasses that of FEM with a heterogeneous structure (Figure10), but FEM costs more to reach the same precision. Thus, CSFEM clearly enhances the calculation efficiency.Figure 10
Computational cost by CSFEM and FEM.A comparison of costs (computation time for the same accuracy) for the homogenized PCM structures indicates CSFEM takes much lower cost than traditional FEM.


### 6.2. A Cantilever PCMs Beam

The free vibrations on a cantilever PCM beam were calculated by CSFEM under the geometrical parameters of lengthL = 20 mm, width H = 2 mm, and plane stress (Figure 11). The PCMs, piezoceramic fibers (PZT-5A), were buried in the non-piezoelectric epoxy matrix. The material properties are shown in Table 3. A block of 16000 square cells was used in the composite z with the specific fiber volume fraction of 66.66%. The relevant boundary condition was ux = uz = Φ = 0 at the clamped end.Figure 11
A cantilever PCMs beam.Table 3
The material properties of the cantilever PCMs beam were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
ρ (kg/m3)
0.6666
13.27801
7.01985
7.87201
42.31662
3.18867
3.12908
0.052676
12.984637
-0.39420
0.42031
5.1058
7500The 11 effective coefficients estimated by AHM are shown in Table3. The first 10 natural frequencies of the cantilever PCM beam with 100 × 10 uniform elements calculated by CSFEM (Figure 12) agree well with the solutions by FEM using 200 × 20 uniform elements. The CSFEM has higher accuracy than FEM. Figure 13 plots the first sixth-order modal shapes, which verify the correctness and validity of CSFEM.Figure 12
First 10 natural frequencies of the cantilever PCMs beam using CSFEM and FEM.Figure 13
First 6 mode shapes of the cantilever PCMs beam by CSFEM and FEM. (a) Mode shape 1. (b) Mode shape 2. (c) Mode shape 3. (d) Mode shape 4. (e) Mode shape 5. (f) Mode shape 6.
(a)
(b)
(c)
(d)
(e)
(f)


### 6.3. A Bilayered PCMs Actuator

In this example, the transient responses of a clamp-free bilayered PCM actuator (L = 20 mm, b = 2 mm, and layer width = 1 mm) exposed to triangular-wave load F (time period T = 8 s, and F0 = 5 N) at point A were investigated (Figures 14 and 15). The actuator consisted of a lower PCM layer and an upper red copper layer, which were supposed to be well bonded. The PCMs (PZT-5A) were buried in the epoxy matrix. A block of 8000 square cells was used in the composite z with the specific fiber volume fraction of 44.44%. The red copper was featured by Young’s modulus (E) = 108 GPa, Poisson ratio (ν) = 0.32, and volume density (ρ) = 8900 kg/m3 (Table 2). The boundary condition of the actuator was ux = uz = Φ = 0 (at the clamped end). The dynamic system equations were addressed by using the implicit integral Newmark scheme.Figure 14
Geometry of the bilayered PCMs actuator.Figure 15
Sine-wave load at pointA.The 12 effective coefficients calculated using AHM are shown in Table4 (Δt = 0.002 s). The generalized displacements at points A and B calculated by CSFEM and FEM are shown in Figures 16–21. Clearly, the high agreement of the simulation results confirms the accuracy of CSFEM. CSFEM achieves higher precision by using fewer elements than FEM (40 × 4 VS 160 × 16 meshes). Thus, CSFEM is well feasible for dynamic analysis of PCM structures.Table 4
The material properties of the bilayered PCMs actuator were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
ρ (kg/m3)
0.4444
7.61822
4.56846
4.77091
28.4779
1.61506
1.52487
0.010859
8.73225
-0.16364
0.20584
3.43606
7500Figure 16
The variation of the displacementux at point A with respect to time.Figure 17
The variation of the displacementuz at point A with respect to time.Figure 18
The variation of the electrical potential Ф at pointA with respect to time.Figure 19
The variation of the displacementux at point B with respect to time.Figure 20
The variation of the displacementuz at point A with respect to time.Figure 21
The variation of the electrical potential Ф at pointA with respect to time.


### 6.4. A PCMs Energy Harvester

The transient responses of a typical PCM energy harvester were investigated at pointsA and B (Figure 22). The mild steel was featured by E = 210 GPa, ρ = 7800.0 kg·m−3, and ν = 0.3, while red copper by E = 108 GPa, ρ = 8900 kg·m−3, and ν = 0.32 (Table 2). The PCMs (PZT-5A) were embedded in the epoxy matrix. A block of 8000 square cells was used in the composite z with the specific fiber volume fraction of 55.55%. The boundary and initial conditions were the same as stated above. The sine-wave payload with 8 s and time step Δt = 0.02 s was applied (Figure 23). The background cells of quadrilateral elements for this harvester were first discretized (Figure 24). The implicit integral Newmark scheme was used to address dynamic system equations.Figure 22
Geometry of PCMs energy harvester.Figure 23
Sine-wave load at top surface of steel.Figure 24
Discrete meshes of PCMs energy harvester. (a) Elements of CSFEM. (b) Elements of FEM.
(a)
(b)The 12 effective coefficients calculated using AHM are shown in Table5. The dynamic responses of the PCM energy harvester were calculated by CSFEM in comparison with FEM (700 VS 5600 quadrilateral elements). The ux, uz, and Ф at points A and B are shown in Figures 25–30. Clearly, the CSFEM produced the generalized displacements closer to real solutions and thereby was proved to reduce the number of meshes and return more-accurate results.Table 5
The material properties of the PCMs energy harvester were calculated using AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
0.5555
9.74292
5.54232
5.95557
35.11748
2.17285
2.1003
0.021904
10.879245
−0.251715
0.28018
4.27248Figure 25
The variation of the displacementux of the PCMs energy harvester at point A with respect to time.Figure 26
The variation of the displacementuz of the PCMs energy harvester at point A with respect to time.Figure 27
The variation of the electrical potential Ф of the PCMs energy harvester at pointA with respect to time.Figure 28
The variation of the displacementux of the PCMs energy harvester at point B with respect to time.Figure 29
The variation of the displacementuz of the PCMs energy harvester at point B with respect to time.Figure 30
The variation of the electrical potential Ф of the PCMs energy harvester at pointB with respect to time.


## 6.1. A Singer-Player PCMs Strip

In this example to test the precision of CSFEM under mechanical and electrical boundary conditions, we used the shear deformation of a piezoelectric strip (1 × 1 mm2, thickness t = 10 μm) under compressive stress σ0 = 5 N·mm−2 and applied voltage V0 = 1000 V (Figure 4). The PCMs, piezoceramic fibers (PZT-5A), were buried in the non-piezoelectric epoxy matrix. The material properties are shown in Table 1 [20]. A block of 10000 square cells (each cell 10 × 10 × 10 μm3) was used in the composite z with the specific fiber volume fraction of 55.55%. One part of the PCM structure (60 × 10 × 40 μm3) is shown in Figure 5. The materials were polarized under the electric field set to the left and right edges, which resulted in shear strain. The boundary conditions were applied to the strip edges:(28)Φ,zx,z=±h=0,σzzx,z=±h=σ0,σxzx=L,z=0,σxzx,z=±h=0,Φx=L,z=−V0,σxxx=L,z=0,Φx=0,z=V0,ux=0,z=0,vx=0,z=0=0.Figure 4
A singer-player PCMs strip under shear deformation.Table 1
Material constants for constituents of the PZT-5A and epoxy matrix.
Material constants
c
11 (GPa)
c
12 (GPa)
c
13 (GPa)
c
33 (GPa)
c
44 (GPa)
c
66 (GPa)
e
13 (C/m2)
e
15 (C/m2)
e
33 (C/m2)
χ
11 (F/m)
χ
33 (F/m)
ρ (kg/m3)
Epoxy matrix
3.86
2.57
2.57
3.86
0.64
0.64
—
—
—
0.0797
0.0797
7500
PZt-5A
121.0
75.4
75.2
111.0
21.1
2.28
−5.4
12.3
15.8
8.11
7.35
7500Figure 5
A part homogeneous model of the singer-player PCMs strip by AHM.The analytical solution to this problem is shown below [58]:(29)u=s13σ0x,v=d15V0xh+s33σ0z,Φ=V01−2xL.It should be noted that all 11 effective coefficients in Table2 were determined using AHM. CSFEM, FEM, and analytical solution were all based on AHM. Three meshes (each 10 × 10, including one with uniform elements and two with distorted elements) were analyzed by CSFEM (Figure 6). The αir is the irregularity factor that is assigned between 0.0 and 0.5 [38]. Meanwhile, 60 × 60 uniform elements were used by FEM. Figure 7 illustrates the CSFEM, FEM, and analytical horizontal displacement ux at the central line (z = 0), while the vertical displacement uz at z = 0 of the single-layer PCM strip is shown in Figure 8. The distributions of the electric potential Φ at z = 0 with FEM and the analytical solutions are demonstrated in Figure 9. Clearly, the displacements and electric potential computed by CSFEM match well with the analytical solutions and outperform those estimated by FEM for all three meshes, which suggested that CSFEM can recreate the linear behavior of the analytical solutions.Table 2
The material properties of the singer-player PCMs strip were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
0.5555
9.74292
5.54232
5.95557
35.11748
2.17285
2.1003
0.021904
10.879245
−0.251715
0.28018
4.27248Figure 6
Typical irregular elements with irregularity factorαir = (a) 0.3 (I1) and (b) 0.5 (I2).
(a)
(b)Figure 7
Variation of horizontal displacementux at the central line (z = 0) of the singer-player PCMs strip.Figure 8
Variation of vertical displacementuz at the central line (z = 0) of the singer-player PCMs strip.Figure 9
Variation of electric potential Φ at the central line (z = 0) of the singer-player PCMs strip.Moreover, the computation cost of CSFEM over PCMs in a homogenized structure surpasses that of FEM with a heterogeneous structure (Figure10), but FEM costs more to reach the same precision. Thus, CSFEM clearly enhances the calculation efficiency.Figure 10
Computational cost by CSFEM and FEM.A comparison of costs (computation time for the same accuracy) for the homogenized PCM structures indicates CSFEM takes much lower cost than traditional FEM.


## 6.2. A Cantilever PCMs Beam

The free vibrations on a cantilever PCM beam were calculated by CSFEM under the geometrical parameters of lengthL = 20 mm, width H = 2 mm, and plane stress (Figure 11). The PCMs, piezoceramic fibers (PZT-5A), were buried in the non-piezoelectric epoxy matrix. The material properties are shown in Table 3. A block of 16000 square cells was used in the composite z with the specific fiber volume fraction of 66.66%. The relevant boundary condition was ux = uz = Φ = 0 at the clamped end.Figure 11
A cantilever PCMs beam.Table 3
The material properties of the cantilever PCMs beam were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
ρ (kg/m3)
0.6666
13.27801
7.01985
7.87201
42.31662
3.18867
3.12908
0.052676
12.984637
-0.39420
0.42031
5.1058
7500The 11 effective coefficients estimated by AHM are shown in Table3. The first 10 natural frequencies of the cantilever PCM beam with 100 × 10 uniform elements calculated by CSFEM (Figure 12) agree well with the solutions by FEM using 200 × 20 uniform elements. The CSFEM has higher accuracy than FEM. Figure 13 plots the first sixth-order modal shapes, which verify the correctness and validity of CSFEM.Figure 12
First 10 natural frequencies of the cantilever PCMs beam using CSFEM and FEM.Figure 13
First 6 mode shapes of the cantilever PCMs beam by CSFEM and FEM. (a) Mode shape 1. (b) Mode shape 2. (c) Mode shape 3. (d) Mode shape 4. (e) Mode shape 5. (f) Mode shape 6.
(a)
(b)
(c)
(d)
(e)
(f)


## 6.3. A Bilayered PCMs Actuator

In this example, the transient responses of a clamp-free bilayered PCM actuator (L = 20 mm, b = 2 mm, and layer width = 1 mm) exposed to triangular-wave load F (time period T = 8 s, and F0 = 5 N) at point A were investigated (Figures 14 and 15). The actuator consisted of a lower PCM layer and an upper red copper layer, which were supposed to be well bonded. The PCMs (PZT-5A) were buried in the epoxy matrix. A block of 8000 square cells was used in the composite z with the specific fiber volume fraction of 44.44%. The red copper was featured by Young’s modulus (E) = 108 GPa, Poisson ratio (ν) = 0.32, and volume density (ρ) = 8900 kg/m3 (Table 2). The boundary condition of the actuator was ux = uz = Φ = 0 (at the clamped end). The dynamic system equations were addressed by using the implicit integral Newmark scheme.Figure 14
Geometry of the bilayered PCMs actuator.Figure 15
Sine-wave load at pointA.The 12 effective coefficients calculated using AHM are shown in Table4 (Δt = 0.002 s). The generalized displacements at points A and B calculated by CSFEM and FEM are shown in Figures 16–21. Clearly, the high agreement of the simulation results confirms the accuracy of CSFEM. CSFEM achieves higher precision by using fewer elements than FEM (40 × 4 VS 160 × 16 meshes). Thus, CSFEM is well feasible for dynamic analysis of PCM structures.Table 4
The material properties of the bilayered PCMs actuator were calculated by AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
ρ (kg/m3)
0.4444
7.61822
4.56846
4.77091
28.4779
1.61506
1.52487
0.010859
8.73225
-0.16364
0.20584
3.43606
7500Figure 16
The variation of the displacementux at point A with respect to time.Figure 17
The variation of the displacementuz at point A with respect to time.Figure 18
The variation of the electrical potential Ф at pointA with respect to time.Figure 19
The variation of the displacementux at point B with respect to time.Figure 20
The variation of the displacementuz at point A with respect to time.Figure 21
The variation of the electrical potential Ф at pointA with respect to time.


## 6.4. A PCMs Energy Harvester

The transient responses of a typical PCM energy harvester were investigated at pointsA and B (Figure 22). The mild steel was featured by E = 210 GPa, ρ = 7800.0 kg·m−3, and ν = 0.3, while red copper by E = 108 GPa, ρ = 8900 kg·m−3, and ν = 0.32 (Table 2). The PCMs (PZT-5A) were embedded in the epoxy matrix. A block of 8000 square cells was used in the composite z with the specific fiber volume fraction of 55.55%. The boundary and initial conditions were the same as stated above. The sine-wave payload with 8 s and time step Δt = 0.02 s was applied (Figure 23). The background cells of quadrilateral elements for this harvester were first discretized (Figure 24). The implicit integral Newmark scheme was used to address dynamic system equations.Figure 22
Geometry of PCMs energy harvester.Figure 23
Sine-wave load at top surface of steel.Figure 24
Discrete meshes of PCMs energy harvester. (a) Elements of CSFEM. (b) Elements of FEM.
(a)
(b)The 12 effective coefficients calculated using AHM are shown in Table5. The dynamic responses of the PCM energy harvester were calculated by CSFEM in comparison with FEM (700 VS 5600 quadrilateral elements). The ux, uz, and Ф at points A and B are shown in Figures 25–30. Clearly, the CSFEM produced the generalized displacements closer to real solutions and thereby was proved to reduce the number of meshes and return more-accurate results.Table 5
The material properties of the PCMs energy harvester were calculated using AHM.
Volume fraction
c
¯
11 (GPa)
c
¯
12 (GPa)
c
¯
13 (GPa)
c
¯
33 (GPa)
c
¯
44 (GPa)
c
¯
66 (GPa)
e
¯
15 (C/m2)
e
¯
33 (C/m2)
e
¯
13 (C/m2)
χ
¯
11 (nF/m)
χ
¯
33 (nF/m)
0.5555
9.74292
5.54232
5.95557
35.11748
2.17285
2.1003
0.021904
10.879245
−0.251715
0.28018
4.27248Figure 25
The variation of the displacementux of the PCMs energy harvester at point A with respect to time.Figure 26
The variation of the displacementuz of the PCMs energy harvester at point A with respect to time.Figure 27
The variation of the electrical potential Ф of the PCMs energy harvester at pointA with respect to time.Figure 28
The variation of the displacementux of the PCMs energy harvester at point B with respect to time.Figure 29
The variation of the displacementuz of the PCMs energy harvester at point B with respect to time.Figure 30
The variation of the electrical potential Ф of the PCMs energy harvester at pointB with respect to time.


## 7. Conclusions

The dynamic responses of PCM structures were investigated. Firstly, the effective properties of PCMs were calculated by AHM. Secondly, the effective CSFEM was established by applying gradient smoothing to calculate the electromechanical coupling field of FEM. Then, the equations for dynamic response computation over the multiphysics coupling field of PCMs were derived. Finally, a bilayered PCM actuator and a PCM energy harvester were calculated by both CSFEM and FEM.(i)
CSFEM reduced the systematic stiffness of FEM and thereby enhanced the accuracy under the same element number. CSFEM took less computation time than FEM under the same accuracy.(ii)
CSFEM avoided the derivation of shape functions and reduced the requirement of form function continuity by simply converting area integral to boundary integral in the solution domain.(iii)
The practical bilayered PCM actuator and the PCM energy harvester were modeled by CSFEM, which outperformed FEM in calculating the transient responses by reducing the number of meshes.



---

*Source: 4913784-2019-05-28.xml*