A New Strong Form Technique for Thermo-Electro-Mechanical Behaviors of Piezoelectric Solids

: Piezoelectric materials are widely fabricated and investigated for potential applications in microelectromechanical systems as direct converters between mechanical and electrical signals, where some show pyroelectric features involving thermo-electro-mechanical interactions. This study aimed to introduce a novel numerical technique to predict the thermo-electro-mechanical behaviors of piezoelectric structures, based on a strong-form numerical framework called the element differential method. In this method, the shape functions of the isoparametric element and their ﬁrst two derivatives were derived analytically by interpolating the temperature, displacement, and electric potentials. Then, a point collocation method based on node positions in the elements was proposed to generate the ﬁnal system of equations without any domain integrations. Thus, the coupled behaviors of thermal piezoelectric structures, including the pyroelectric features, can be simulated by the strong-form formulation of the governing equations. Several numerical examples, including the piezoelectric composites structures, are presented, and the coupled thermo-electro-mechanical responses have been analyzed to validate the proposed method.


Introduction
With the fast development of smart materials and manufacturing technologies, the thermo-electro-mechanical coupled responses of the smart structures have attracted instant attention from scientists. P piezoelectric materials have been broadly employed as actuators or sensors depending on the coupled mechanical and electrical characteristics. These devices can provide high-level features such as high force generation, response speed, and displacement accuracy [1,2]. Moreover, because thermo-piezoelectric materials possess the potential capability in aviation, transportation, and automating industries [3], they have been also widely selected in designing the devices for damage detecting, health inspecting, actuating, and load-carrying applications [4,5]. For instance, research indicates that the sensing devices based on pyroelectric and piezoelectric effects have potential applications in the design of active human-computer interaction and wearable robots [6,7].
The smart piezoelectric devices composed of thermo-electro-elastic properties often run at the different loading environments, such as cyclic electrical or mechanical, thermal, and or their combined conditions. In literature, a significant amount of research focuses on simulating and investigating the coupling structural responses of piezoelectric materials. Completely understanding the coupling behaviors or the physical laws will be beneficial in utilizing the thermo-piezoelectric materials in real engineering applications. The research [8,9] demonstrated that the polarization of the crystalline materials and structures could be affected by the temperature. Due to the temperature changes, the piezoelectric material may be polarized by establishing the electrical potentials [10].
On the other hand, a complete understanding of thermal stress is necessary to guarantee the reliabilities of piezoelectric structures in real applications, such as the aerospace industry, especially when piezoelectric structures are exposed to complex and changeable temperature environments [11,12]. Namli and Taya introduced a thermal energy piezo-SMA harvester subjected to temperature fluctuation in different frequencies [13]. Srivastava et al. designed a device to directly convert heat to electricity, in which energy conversion devices are designed according to the distributions of the electro-thermo-elastic coupling fields [14]. In addition, thermal effect analyses of the thermal piezoelectric structures could provide a theoretical basis of piezoelectric materials to predict the crack failure behaviors [15]. Consequently, the accurate analysis of the thermo-electro-mechanical behaviors of the piezoelectric materials under thermal environments is of great importance to understand the mechanism of interaction among different physical fields.
In this regard, recent literature has concentrated on developing novel numerical methods to analyze the multi-field coupling problems of the piezoelectric structures subjected to electrical, mechanical, and thermal loadings. Several illustrative numerical methods have been developed to analyze the thermo-electro-mechanical responses of piezoelectric structures, such as the finite element methods (FEM), the finite difference method (FDM), and the boundary element methods (BEM). Among them, the FEM is one of the most convenient computational tools. It has been widely utilized to calculate the boundary value problems of the piezoelectric medium subjected to complex loading conditions. For instance, the finite element formulation has been utilized to solve the crack problems under thermal loadings for piezoelectric structures with different size scales [16]. Rahman used FEM to analyze the free vibration analysis and steady harmonic analysis of a multilayer porous cantilever beam with a piezoelectric layer [17]. Jiang and Li introduced a finite element formulation to compute the energy harvesting problems of the piezoelectric composite beams under thermo-electro-mechanical loadings [18]. Amini et al. selected the FEM to calculate the coupled response of the piezoelectric harvesters made of functional graded materials in unimorph or bimorph arrangement [19].
Meanwhile, Yan and Jiang studied the dynamic fracture behaviors of piezoelectric materials based on the same numerical method [20]. The standard FEM to simulate the coupled problems of piezoelectric structures has been widely integrated into some of the commercial finite element codes, such as ABAQUS and ANSYS [21]. Although the overall success has been achieved through the finite element method to analyze piezoelectric problems, to the best of the authors' knowledge, even some of the commercial finite element codes still cannot simulate the pyroelectric behaviors directly due to the complexity of the problems. Therefore, it will be interesting to study the pyroelectric properties of piezoelectric cellular structures with more convenient methods. For instance, Lv et al. introduced an extended multiscale finite element method for thermoelectric coupling responses of piezoelectric structures with pyroelectric effects [22].
The numerical methods derived from strong formulations to calculate piezoelectric behaviors have received growing interest in recent years. Generally, the solutions based on the strong-form technique can directly generate the system of equations without any other supplementary methods [23,24]. Thus, the strong-form method shows unrivaled favorable circumstances to simulate complex engineering problems, especially for multiphysics problems with different governing equations [25]. For instance, Moradi-Dastjerdi et al. used a meshless method to investigate the damped dynamic motions of piezoelectric sandwich plates under thermo-electromechanical loads [26]. Liew proposed a Galerkin differential quadrature method (DQM) to simulate the nonlinear behaviors of the thermopiezoelectric structures [27]. Alashti and Khorsand also proposed a DQM to calculate the three-dimensional thermo-elastic problems of the cylindrical shells composed of piezoelectric layers under the thermo-electro-mechanical loads [28]. Recently, Lv et al. presented a type of strong form method to predict response of piezoelectric structures [29]. Their research indicated that the element differential method derived from the strong-form technique combines the benefits of standard FEMs, i.e., mesh grids can be generated adaptively for arbitrary configurations in real applications. At the same time, the strong form of governing equations can be directly utilized.
In this work, the research focused on the strong-form techniques to simulate the coupled responses of thermo-piezoelectric structures. A novel strong-form collocation method is presented to solve the coupled thermo-piezoelectric-mechanical behaviors. The isoparametric elements were introduced to assemble the governing partial differential equations of coupling problems. Then, according to the distributions of the nodes, a special collocation technique is used to construct the system equations for the coupled problems directly. The uncoupled procedure based on the standard Fourier heat conduction theory is utilized for thermo-electro-elastic behaviors. Furthermore, an additional innovation of the present strong-form method is that the pyroelectric effects of the piezoelectric structure are directly substituted into the system of equations through the collocation technique without any domain integrations. Finally, serval illustrative examples, including the complex piezoelectric composite structures, are proposed to validate the accuracy of the current method.
This paper's remaining is organized in the following manner: The main equations, including the governing equations and constitutive equations, are briefly introduced for the linear thermo-electro-elastic problems in Section 2. Section 3 presents the element differential formulations and the point collocation technique to construct the system equations to analyze the piezoelectric structures subjected to the thermo-electro-elastic loadings. The point collocation technique is separately introduced for the two types of problems. Section 4 gives some numerical examples for the code validation. The last section gives a summary and conclusions.

Unified Governing Equations for Thermo-Electro-Mechanical Problems
As indicated in Figure 1, a linear isotropic homogeneous piezoelectric domain subjected to the thermo-electro-mechanical loadings is considered. In this domain, the mechanical boundary conditions and prescribed electrical potentials are applied to the boundaries Γ u , Γ φ , respectively. At the same time, the traction t and the heat flux q act on boundaries Γ t , Γ q , respectively. piezoelectric structures [29]. Their research indicated that the element differential method derived from the strong-form technique combines the benefits of standard FEMs, i.e., mesh grids can be generated adaptively for arbitrary configurations in real applications. At the same time, the strong form of governing equations can be directly utilized. In this work, the research focused on the strong-form techniques to simulate the coupled responses of thermo-piezoelectric structures. A novel strong-form collocation method is presented to solve the coupled thermo-piezoelectric-mechanical behaviors. The isoparametric elements were introduced to assemble the governing partial differential equations of coupling problems. Then, according to the distributions of the nodes, a special collocation technique is used to construct the system equations for the coupled problems directly. The uncoupled procedure based on the standard Fourier heat conduction theory is utilized for thermo-electro-elastic behaviors. Furthermore, an additional innovation of the present strong-form method is that the pyroelectric effects of the piezoelectric structure are directly substituted into the system of equations through the collocation technique without any domain integrations. Finally, serval illustrative examples, including the complex piezoelectric composite structures, are proposed to validate the accuracy of the current method. This paper's remaining is organized in the following manner: The main equations, including the governing equations and constitutive equations, are briefly introduced for the linear thermo-electro-elastic problems in Section 2. Section 3 presents the element differential formulations and the point collocation technique to construct the system equations to analyze the piezoelectric structures subjected to the thermo-electro-elastic loadings. The point collocation technique is separately introduced for the two types of problems. Section 4 gives some numerical examples for the code validation. The last section gives a summary and conclusions.

Unified Governing Equations for Thermo-Electro-Mechanical Problems
As indicated in Figure 1, a linear isotropic homogeneous piezoelectric domain subjected to the thermo-electro-mechanical loadings is considered. In this domain, the mechanical boundary conditions and prescribed electrical potentials are applied to the boundaries Γ , Γ , respectively. At the same time, the traction t and the heat flux q act on boundaries Γ , Γ , respectively. The constitutive equations for thermo-electro-mechanical problems can be written as [22]: The constitutive equations for thermo-electro-mechanical problems can be written as [22]: where σ ij are the stress tensors, ε kl strain tensors, D i electric displacement vectors, E m electric field vectors, C ijkl elastic coefficient vectors, e mij piezoelectric stress coefficient vectors, λ ik dielectric coefficient vectors, β ij thermal-stress coefficient vectors and χ i pyroelectric coefficient vectors, T the temperature variations. In two-dimensional problems, the ranges for the repeated subscripts i and j are 2, while in three-dimensional problems, the ranges are 3. The parameters used in Equation (1) for three dimensional problems can be written in a unified way to simplify the expression of the constitutive equations for the thermopiezoelectric materials In the thermo-piezoelectric materials, the constitutive equations can be rewritten in an extended form σ I j = L I jMn ε Mn + π I j T In Equation (3), L I jMn and π I j denote the constitutive coefficient matrixes, which can be written as where, the ranges of the subscripts I and M is 3 for 2D and 4 for 3D problems. With the symbolic forms presented above, the prescribed boundary conditions in the mechanical and electrical fields of the piezoelectric structures can be rewritten as: where U i (x) denotes the fixed displacement and electrical potential values applied at the boundary Γ L , the character L M denotes the unified form about the traction vectors and surface charges. The outward normal vector n n refers to the boundaries Γ L . All the boundaries of the problem satisfy Therefore, based on the extended form of the parameters, the obtained second-order partial differential equations of the thermal piezoelectric structure can be written as: where F I denotes the body force. At the same time, the boundary conditions related to the tractions and surface charges can be expressed as On the other hand, for the thermal conduction problems of the piezoelectric structure, the governing equation about the steady-state heat conduction response can be written as where, k ij (x) is the heat conductivity varying with coordinates x, Q(x) is the heat generated rate. Three types of boundary conditions can be generally considered to solve heat transfer problems of the piezoelectric solids where Γ 1 ∪ Γ 2 ∪ Γ 3 = Γ, T(x) and T are the temperature and prescribed values, respectively. q is the specified heat flux. T ∞ is the temperature on the boundary Γ 3 . h denotes the heat transfer coefficients.

Elemental Differential Method for Thermo-Electro-Mechanical Problems
In this section, a novel strong-form method called the thermal piezoelectric elemental differential method (TP-EDM) is introduced to analyze the coupled responses of the structures under thermo-electro-mechanical loadings. The proposed method can be divided into three steps. In the first step, target domains are discretized into isoparametric elements, consistent with those standard finite element methods. In the second step, the final system of equations is generated by the point collocation method. These equations are generated node by node in this step, while domain integrations element by element construct the total stiffness matrix. In the final step, the final system of the equation is solved to calculate the unknowns, i.e., nodal displacements, electrical potentials, and temperatures.

Lagrange Isoparametric Elements
The standard C0-continuous 2D and 3D Lagrange isoparametric elements are used to discretize the computation domain [25]. As illustrated in Figure 2, a 2D quadrilateral element with nine nodes is considered here. The main characteristic of the elements is that an independent node is distributed inside the elements. Based on the isoparametric element, the geometrical and physical variables over the element can be approximated according to the nodal values of the element. Therefore, spatial coordinates, displacement, electrical potential, and temperature within an element can be calculated with the interpolation functions by (11) where, N α denotes the shape functions. The subscript α denotes the summation ranged from 1 to the number of the elemental nodes. It can be found that the mesh grids used in the standard finite element method can be directly employed for the computations of the current method.

Lagrange Isoparametric Elements
The standard C0-continuous 2D and 3D Lagrange isoparametric elements are used to discretize the computation domain [25]. As illustrated in Figure 2, a 2D quadrilateral element with nine nodes is considered here. The main characteristic of the elements is that an independent node is distributed inside the elements. Based on the isoparametric element, the geometrical and physical variables over the element can be approximated according to the nodal values of the element. The shape functions N α for the isoparametric elements can be explicitly derived with the intrinsic coordinates [25]. In this regard, the first two orders of the partial differential equations for the shape functions referring to the spatial coordinates can be derived explicitly. The transformations of the derivatives can be written as follows: where, the Jacobian matrix J denotes the mapping relationship between the Cartesian coordinate and the intrinsic coordinate systems. Based on the analytical expansion of the shape functions, the Jacobian matrix, its inverse matrix [J] −1 ik and derivative matrixes can be derived analytically, respectively [25].

Point Collocation Technique
This work assumes that the temperature distributions over the piezoelectric structures are independent of the electrical and elastic fields. The heat transfer problem is firstly calculated without considering the effects of the mechanical fields. Then the thermo-stress and the pyroelectric features of the piezoelectric structures are solved based on the obtained thermal fields. Although two different types of governing equations govern the two types of problems, the system equations of these boundary value problems can be constructed by the collocation technique in a unified way.

Point Collocation Technique for Thermal Conductions
As shown in Figure 3, the element nodes in the structures meshed with isoparametric elements can be clarified into three types according to their distributions, i.e., internal nodes located within the elements, interface nodes between elements, and outer surface nodes along boundaries.  For internal nodes located within an element, the set of equations can be directly created by governing equations of the thermal conduction problems. By combining Equations (7) and (12), the governing equations are expressed by the nodal element values of temperature [30] [ In the above equation, represents the intrinsic coordinates of the nodes, the value of at node . On the other hand, the heat flux should satisfy the equilibrium conditions at the other nodes distributed at element boundaries or interfaces shared by several elements. The equilibrium equations of the heat flux can be written as follows: where M denotes the quantity of element surfaces related to the interface nodes.
( ) is the heat flux of surface f; is the prescribed value of the heat flux. For the nodes shared by several elements, the relationship between heat flux and temperature gradients can be obtained by substituting Equation (9) into Equation (14) where is the intrinsic coordinates of the interface nodes; is the outward normal vector to surface f.
On the other hand, for the nodes on the outer boundaries of mesh grids, the heat flux equilibrium conditions should be satisfied [30] where, is the intrinsic coordinates of the outer boundary nodes, ( ) = ̅ at the boundary Γ 2 , and ( ) = ℎ( ( ) − ∞ ) at the boundary Γ 3 .

Point Collocation Techniques for the Thermal Elastic Problems
A similar way can generate the system of equations for the thermal elastic response of the thermal piezoelectric structure as those for heat transfer problems. The governing equations for the thermal stress problems should be satisfied at all nodes inside elements For internal nodes located within an element, the set of equations can be directly created by governing equations of the thermal conduction problems. By combining Equations (7) and (12), the governing equations are expressed by the nodal element values of temperature [30] In the above equation, ξ represents the intrinsic coordinates of the nodes, k β ij the value of k ij at node β.
On the other hand, the heat flux should satisfy the equilibrium conditions at the other nodes distributed at element boundaries or interfaces shared by several elements. The equilibrium equations of the heat flux can be written as follows: where M denotes the quantity of element surfaces related to the interface nodes. q f i (x) is the heat flux of surface f ; q i is the prescribed value of the heat flux.
For the nodes shared by several elements, the relationship between heat flux and temperature gradients can be obtained by substituting Equation (9) into Equation (14) where ξ I is the intrinsic coordinates of the interface nodes; n f i is the outward normal vector to surface f.
On the other hand, for the nodes on the outer boundaries of mesh grids, the heat flux equilibrium conditions should be satisfied [30] where, ξ b is the intrinsic coordinates of the outer boundary nodes, q ξ b = q at the boundary Γ 2 , and q ξ b = h T ξ b − T ∞ at the boundary Γ 3 .

Point Collocation Techniques for the Thermal Elastic Problems
A similar way can generate the system of equations for the thermal elastic response of the thermal piezoelectric structure as those for heat transfer problems. The governing equations for the thermal stress problems should be satisfied at all nodes inside elements where, ξ denotes intrinsic coordinates of the node. The thermal-stress coefficients π I j can be directly considered in the governing equation. In this method the temperature difference is directly utilized in the governing equations to create the system of equations. Thus, the complexity in dealing with the thermal coupling problems of the piezoelectric structures can be greatly reduced.
For the other nodes, traction equilibrium conditions should be satisfied.
In above equation, N represents the quantity of elements that includes the node considered, S e the quantity of surfaces of the element related to the node, L i (x) the surface tractions and charges applied.
Therefore, when the nodes are located at the element interfaces, the tractions for all the faces related to the considered nodes should be at an equilibrium state. The equation can be written as follows: where the symbol n j represents outward normal to the surface.

The Final System of Equations for the Coupled Problem
The system of equations for the heat conduction and thermal-mechanical problems of the piezoelectric structures can be obtained based on the collocation technique mentioned above. According to their distributions, each of these nodes can be used to assemble parts of the system of equations with governing equations or boundary conditions. The unknowns in the final system of the equation need to be numbered in a global sequence. For instance, the unknowns of the thermal conduction problem are related to nodal temperature. Each node can generate an equation containing all nodes, including the interface nodes. At the same time, the nodes at the outer boundaries can generate equations containing node temperatures and heat fluxes by substituting intrinsic coordinates ξ into Equations (15) and (16). Finally, once global nodal vectors of temperatures are numbered in a global sequence, then the final system of the equation can be written as [25]: where [A] denotes the coefficient matrix of the system of equation; {x} denotes solution vector and {b} right side vector. It should be noted that for the thermo-electro-mechanical problems, the global nodal vector is composed of the displacement and the electrical potentials. Therefore, the heat conduction and thermal, mechanical problems can be performed by solving the algebraic equations. Each node's stresses and surface charges for the thermalmechanical problems can be directly calculated by substituting the displacement and electrical potentials into the constitutive equations.

Numerical Examples
In this section, the computational performance of the proposed element differential method for thermo-electro-mechanical structures (TP-EDM) is validated through several numerical examples.

Piezoelectric Composite Structure Composed of Periodical Unit Cells
A piezoelectric composite structure composed periodically distributed unit cells (4 × 10) is considered in the first example. As indicated in Figure 4, the whole size of the square unit cell is 60 mm × 60 mm. The unit cell is made of two types of materials (Epoxy and PZT-7A). The diameter of the circular piezoelectric fiber is 30 mm. The left end of this structure is firmly fixed. The bottom temperature value is set to 100 • C, while the value for the top face is set to 500 • C. At the same time, the electrical potential at this part of the structure is set to zero. The material constants for the piezoceramic and non-piezoelectric matrix are listed as follow: where E is the elastic modulus of the materials, µ Poisson's ratio.

Piezoelectric Composite Structure Composed of Periodical Unit Cells
A piezoelectric composite structure composed periodically distributed unit cells (4 × 10) is considered in the first example. As indicated in Figure 4, the whole size of the square unit cell is 60 mm × 60 mm. The unit cell is made of two types of materials (Epoxy and PZT-7A). The diameter of the circular piezoelectric fiber is 30 mm. The left end of this structure is firmly fixed. The bottom temperature value is set to 100 °C , while the value for the top face is set to 500 °C . At the same time, the electrical potential at this part of the structure is set to zero. The material constants for the piezoceramic and non-piezoelectric matrix are listed as follow: PZT-7A: where E is the elastic modulus of the materials, Poisson's ratio. As shown in Figure 4, the piezoelectric composite structure is divided into 44,000 four-node quadrilateral elements. Both TP-EDM and a household standard finite element method developed are utilized to solve the coupled thermo-electro-mechanical problems. The displacement and electric potential results computed by these methods are plotted in Figures 5 and 6. It could be observed that a bending deformation will be generated under the uniform thermal temperature. The maximum displacement in the Y direction occurs at the top-right corner (Point A). At the same time, the electrical potential will be induced by the uniform expansion of the structure. As indicated in Figure 6, the maximum value of the electric potentials occurs at the right end of the structure, which is consistent with the deformation of the structure.
Furthermore, the maximum electrical potential and displacement at point A are listed in Table 1  As shown in Figure 4, the piezoelectric composite structure is divided into 44,000 four-node quadrilateral elements. Both TP-EDM and a household standard finite element method developed are utilized to solve the coupled thermo-electro-mechanical problems.
The displacement and electric potential results computed by these methods are plotted in Figures 5 and 6. It could be observed that a bending deformation will be generated under the uniform thermal temperature. The maximum displacement in the Y direction occurs at the top-right corner (Point A). At the same time, the electrical potential will be induced by the uniform expansion of the structure. As indicated in Figure 6, the maximum value of the electric potentials occurs at the right end of the structure, which is consistent with the deformation of the structure.
Furthermore, the maximum electrical potential and displacement at point A are listed in Table 1

Clamped Thermal Piezoelectric 3D Beam
In this example, the thermo-electro-mechanical response of a clamped thermal piezoelectric 3D beam is considered. The two ends of the piezoelectric beam are fixed, as shown in Figure 8. The size of the beam is 0.005 m× 0.1 m × 0.01 m. Since the commercial finite element method code cannot be directly used to model the piezoelectric structures with the pyroelectric behaviors, the pyroelectric coefficients for the 3D beam are set to zero here. The other material constants can be given as follows:

Clamped Thermal Piezoelectric 3D Beam
In this example, the thermo-electro-mechanical response of a clamped thermal piezoelectric 3D beam is considered. The two ends of the piezoelectric beam are fixed, as shown in Figure 8. The size of the beam is 0.005 m× 0.1 m × 0.01 m . Since the commercial finite element method code cannot be directly used to model the piezoelectric structures with the pyroelectric behaviors, the pyroelectric coefficients for the 3D beam are set to zero here. The other material constants can be given as follows: Two loading cases separately composed of thermal loading and the combined thermal-electric and mechanical loading are considered in this example, respectively, to validate the effectiveness of the presented method.
Case 1: The clamped thermal piezoelectric beam structure considered is subjected to the thermal loadings. A temperature load (T = 200 K) is subjected to the beam structure, while the electrical potential is set to zero at the bottom of the structure. The piezoelectric beam will expand under the temperatures. The beam structure meshes with 900 grids. The commercial finite element codes Abaqus with the C3D8P elements are applied to validate the proposed method TP-EDM. Since the same mesh grids are utilized, the results calculated by Abaqus can be seen as the reference results.
The deformation and the electrical potential responses predicted by both TP-EDM and Abaqus can be found in Figures 9 and 10, respectively. It can be found that a uniform expansion happens at most of the beam structures except the two ends, which are constrained by the fixed boundary conditions. Thus, thermal stress will be induced, which will lead to electrical potentials. At the same time, the maximum values about displacement and electrical potential induced by temperature variations in the structure can be found in Table 2. Furthermore, Figure 11 gives the displacement and electrical potential for the nodes at the bottom line. The results computed by the presented method are in good agreement with those of the standard finite element method.  Two loading cases separately composed of thermal loading and the combined thermalelectric and mechanical loading are considered in this example, respectively, to validate the effectiveness of the presented method.
Case 1: The clamped thermal piezoelectric beam structure considered is subjected to the thermal loadings. A temperature load (T = 200 K) is subjected to the beam structure, while the electrical potential is set to zero at the bottom of the structure. The piezoelectric beam will expand under the temperatures. The beam structure meshes with 900 grids. The commercial finite element codes Abaqus with the C3D8P elements are applied to validate the proposed method TP-EDM. Since the same mesh grids are utilized, the results calculated by Abaqus can be seen as the reference results.
The deformation and the electrical potential responses predicted by both TP-EDM and Abaqus can be found in Figures 9 and 10, respectively. It can be found that a uniform expansion happens at most of the beam structures except the two ends, which are constrained by the fixed boundary conditions. Thus, thermal stress will be induced, which will lead to electrical potentials. At the same time, the maximum values about displacement and electrical potential induced by temperature variations in the structure can be found in Table 2. Furthermore, Figure 11 gives the displacement and electrical potential for the nodes at the bottom line. The results computed by the presented method are in good agreement with those of the standard finite element method.
Case 2: In this case, the combined thermal and mechanical loadings are applied on the clamped piezoelectric beam. The upper surface of the piezoelectric beam is subjected to a uniform pressure (10MPa), while the whole model is still subjected to a temperature load of 200 K, as illustrated in Figure 8. Both of the two ends of the piezoelectric beam are fixed. At the same time, the electrical potential is set to zero at the bottom face.          In order to further verify the influence of the mesh generation type on the accuracy of the presented method, a detailed mesh refinement convergence study is performed. Three types of mesh grids are considered, consisting of 900, 5000, and 12,103 brick elements, respectively. In the last one, the irregular mesh grids with different types of hexahedral elements are considered. The deformation and the electrical potential generated by the thermal and mechanical environment can be found in Figures 12 and 13, respectively. It can be found that under the combined coupled loading, the bending deformation will be generated at the clamped piezoelectric beam structure. The electrical potential generated by the mechanical loading will gather at the top surface of the structure. to a uniform pressure (10MPa), while the whole model is still subjected to a temperature load of 200 K, as illustrated in Figure 8. Both of the two ends of the piezoelectric beam are fixed. At the same time, the electrical potential is set to zero at the bottom face. In order to further verify the influence of the mesh generation type on the accuracy of the presented method, a detailed mesh refinement convergence study is performed. Three types of mesh grids are considered, consisting of 900, 5000, and 12,103 brick elements, respectively. In the last one, the irregular mesh grids with different types of hexahedral elements are considered. The deformation and the electrical potential generated by the thermal and mechanical environment can be found in Figures 12 and 13, respectively. It can be found that under the combined coupled loading, the bending deformation will be generated at the clamped piezoelectric beam structure. The electrical potential generated by the mechanical loading will gather at the top surface of the structure.  Besides, it can be observed that the deformation and the electrical response predicated by TP-EDM are finely consistent with those calculated by Abaqus. The displacement and electrical potential values along the Y-axis of the upper surface of the 3D piezoelectric beam are plotted in Figure 14. At the same time, the results with different mesh grids can be found in Table 3. It can be seen from Figure 13 that the maximum electrical potential appeared at some part of the top surfaces, which ranges from 20mm to 80mm in the Y direction along the structure. Furthermore, the curve is approximately saturated in this range, shown in Figure 14b, while the maximum value happens at the center point of the structure. It can be observed that the maximum values of the displacement and electrical potential computed by the proposed method vary little with the number of the mesh grids (from 900 to 12103), which demonstrates a relatively high convergence rate of the proposed method. Thus, it can be concluded that the TP-EDM has to a uniform pressure (10MPa), while the whole model is still subjected to a temperature load of 200 K, as illustrated in Figure 8. Both of the two ends of the piezoelectric beam are fixed. At the same time, the electrical potential is set to zero at the bottom face. In order to further verify the influence of the mesh generation type on the accuracy of the presented method, a detailed mesh refinement convergence study is performed. Three types of mesh grids are considered, consisting of 900, 5000, and 12,103 brick elements, respectively. In the last one, the irregular mesh grids with different types of hexahedral elements are considered. The deformation and the electrical potential generated by the thermal and mechanical environment can be found in Figures 12 and 13, respectively. It can be found that under the combined coupled loading, the bending deformation will be generated at the clamped piezoelectric beam structure. The electrical potential generated by the mechanical loading will gather at the top surface of the structure.  Besides, it can be observed that the deformation and the electrical response predicated by TP-EDM are finely consistent with those calculated by Abaqus. The displacement and electrical potential values along the Y-axis of the upper surface of the 3D piezoelectric beam are plotted in Figure 14. At the same time, the results with different mesh grids can be found in Table 3. It can be seen from Figure 13 that the maximum electrical potential appeared at some part of the top surfaces, which ranges from 20mm to 80mm in the Y direction along the structure. Furthermore, the curve is approximately saturated in this range, shown in Figure 14b, while the maximum value happens at the center point of the structure. It can be observed that the maximum values of the displacement and electrical potential computed by the proposed method vary little with the number of the mesh grids (from 900 to 12103), which demonstrates a relatively high convergence rate of the proposed method. Thus, it can be concluded that the TP-EDM has Besides, it can be observed that the deformation and the electrical response predicated by TP-EDM are finely consistent with those calculated by Abaqus. The displacement and electrical potential values along the Y-axis of the upper surface of the 3D piezoelectric beam are plotted in Figure 14. At the same time, the results with different mesh grids can be found in Table 3. It can be seen from Figure 13 that the maximum electrical potential appeared at some part of the top surfaces, which ranges from 20 mm to 80 mm in the Y direction along the structure. Furthermore, the curve is approximately saturated in this range, shown in Figure 14b, while the maximum value happens at the center point of the structure. It can be observed that the maximum values of the displacement and electrical potential computed by the proposed method vary little with the number of the mesh grids (from 900 to 12103), which demonstrates a relatively high convergence rate of the proposed method. Thus, it can be concluded that the TP-EDM has satisfied numerical accuracy and stability to solve the coupled thermal-electrical-mechanical problems of the piezoelectric structure.   Table 3. Displacement in negative direction (m) and electrical potential(V) for the piezoelectric beam.

Cylindrical Structure Reinforced by Piezoelectric Fiber
In the final example, a composite structure reinforced by piezoelectric fiber under combined electro-thermo-mechanical loading is considered. This fiber-reinforced piezoelectric composite has received growing attention recently since it can provide novel properties and individual synergism between different types of materials [31]. As illustrated in Figure 15, the hollow cylindrical structure comprises a thin carbon shell, elastic core, and piezoelectric fibers. It can be found that the elastic core contains 16 sticks made of PZT. Since the model and the loading conditions are symmetric, only 1/4 model is selected to simulate the coupling behaviors. The detailed geometrical information of the structure can be found in Figure 15, in which the thickness of the carbon fiber film shell is 0.15 mm, the length of the cylinder L = 100 mm, the diameter of the piezoelectric fiber d = 10 mm, the radius of the hollow inclusion R c = 30 mm, while the thickness of the hollow core is Δ = R s − R c = 46.2 mm.

Cylindrical Structure Reinforced by Piezoelectric Fiber
In the final example, a composite structure reinforced by piezoelectric fiber under combined electro-thermo-mechanical loading is considered. This fiber-reinforced piezoelectric composite has received growing attention recently since it can provide novel properties and individual synergism between different types of materials [31]. As illustrated in Figure 15, the hollow cylindrical structure comprises a thin carbon shell, elastic core, and piezoelectric fibers. It can be found that the elastic core contains 16 sticks made of PZT. Since the model and the loading conditions are symmetric, only 1/4 model is selected to simulate the coupling behaviors. The detailed geometrical information of the structure can be found in Figure 15, in which the thickness of the carbon fiber film shell is 0.15 mm, the length of the cylinder L = 100 mm, the diameter of the piezoelectric fiber d = 10 mm, the radius of the hollow inclusion R c = 30 mm, while the thickness of the hollow core is ∆ = R s − R c = 46.2 mm. satisfied numerical accuracy and stability to solve the coupled thermal-electricalmechanical problems of the piezoelectric structure.

Cylindrical Structure Reinforced by Piezoelectric Fiber
In the final example, a composite structure reinforced by piezoelectric fiber under combined electro-thermo-mechanical loading is considered. This fiber-reinforced piezoelectric composite has received growing attention recently since it can provide novel properties and individual synergism between different types of materials [31]. As illustrated in Figure 15, the hollow cylindrical structure comprises a thin carbon shell, elastic core, and piezoelectric fibers. It can be found that the elastic core contains 16 sticks made of PZT. Since the model and the loading conditions are symmetric, only 1/4 model is selected to simulate the coupling behaviors. The detailed geometrical information of the structure can be found in Figure 15, in which the thickness of the carbon fiber film shell is 0.15 mm, the length of the cylinder L = 100 mm, the diameter of the piezoelectric fiber d = 10 mm, the radius of the hollow inclusion R c = 30 mm, while the thickness of the hollow core is Δ = R s − R c = 46.2 mm.
In order to simulate the coupled response of the piezoelectric fiber-reinforced composites, two loading cases are considered. In the first case, the loading and boundary conditions could be described as below: at the bottom of the structure (Z = 0), the electrical potential is set to zero, the normal displacement is constrained, the temperature is set to 20 • C. At the same time, the applied pressure is 100 MPa at the top face of the structure (Z = 100), while the temperature of the face is 100 • C. The boundary condition can be listed as follow: The piezoelectric fiber-reinforced composite cylindrical structure meshes with 23,902 elements. The displacement and electrical potential result of the 1/4 structure calculated by TP-EDM and Abaqus are given in Figures 16 and 17. The displacement and electrical potential data along the line at the position (X = 0, Y = 76.35 mm) are given in Figure 18. It can be seen from Figure 18 that the displacement and potential results curves calculated by TP-EDM are basically consistent as those calculated by Abaqus in the coupling problem of heat conduction and mechanical force. Moreover, the displacement and electrical potential created by Abaqus and TP-EDM of the one point (point B) on the structure are listed in Table 4. It was evident that the results obtained by these two methods are well-matched, although different types of materials are considered.      For the second case, to verify the accuracy and diversity of the method, the combined thermo-electro-mechanical boundary condition is considered in this example. The heat exchange boundary conditions are applied on the two end faces of the piezoelectric fibers.     For the second case, to verify the accuracy and diversity of the method, the combined thermo-electro-mechanical boundary condition is considered in this example. The heat exchange boundary conditions are applied on the two end faces of the piezoelectric fibers.  For the second case, to verify the accuracy and diversity of the method, the combined thermo-electro-mechanical boundary condition is considered in this example. The heat exchange boundary conditions are applied on the two end faces of the piezoelectric fibers. As indicated in Figure 15a, at one end of the cylinder, the heat flux of the PZT face is set to q = 0.2 Wm −2°C , while the reference temperature is 60 • C. At the other end, the heat flux is set to q = 0.1 Wm −2°C , while the environmental temperature is 200 • C. These boundary conditions applied on the target structure are listed as follow: PZT fibers: The mesh grid used in this case is the same as that for the first case. The temperature distribution of the cylindrical structure calculated by Abaqus and TP-EDM is shown in Figure 19. It can be observed that as the heat flux is only applied on the end of the piezoelectric fibers, the resulting temperature at corresponding areas is higher than that of the matrix. The temperature calculated by TP-EDM agrees well with that simulated by ABAQUS, which shows the accuracy of TP-EDM in solving the heat transfer problem.
The mesh grid used in this case is the same as that for the first case. The temperature distribution of the cylindrical structure calculated by Abaqus and TP-EDM is shown in Figure 19. It can be observed that as the heat flux is only applied on the end of the piezoelectric fibers, the resulting temperature at corresponding areas is higher than that of the matrix. The temperature calculated by TP-EDM agrees well with that simulated by ABAQUS, which shows the accuracy of TP-EDM in solving the heat transfer problem. Furthermore, once the calculated temperature is calculated for the piezoelectric composite structure, the other high order responses such as the electrical potential and stress can be further calculated by the proposed method. As illustrated in Figures 20 and  21, the left end with higher temperatures generates high electric potentials. At the same time, the von Mises stress results mainly happen at the shell surrounding the cylinder caused by the thermal expansion of the core. Moreover, the computed results by TP-EDM are in good agreement with those calculated by Abaqus. Therefore, we can conclude that TP-EDM can obtain the results with satisfying accuracy in dealing with the piezoelectric composites subjected to combined thermo-electrical-mechanical loadings. Further, since thermal loads are directly utilized in the governing equations to establish the system of equations without any domain integrations, the proposed method can reduce the complexity in dealing with the thermal coupling problems of the piezoelectric structures. Furthermore, once the calculated temperature is calculated for the piezoelectric composite structure, the other high order responses such as the electrical potential and stress can be further calculated by the proposed method. As illustrated in Figures 20 and 21, the left end with higher temperatures generates high electric potentials. At the same time, the von Mises stress results mainly happen at the shell surrounding the cylinder caused by the thermal expansion of the core. Moreover, the computed results by TP-EDM are in good agreement with those calculated by Abaqus. Therefore, we can conclude that TP-EDM can obtain the results with satisfying accuracy in dealing with the piezoelectric composites subjected to combined thermo-electrical-mechanical loadings. Further, since thermal loads are directly utilized in the governing equations to establish the system of equations without any domain integrations, the proposed method can reduce the complexity in dealing with the thermal coupling problems of the piezoelectric structures.

Conclusions
In this study, a new strong form method, the elemental differential method, has been proposed to simulate the coupled response of the piezoelectric structures subjected to

Conclusions
In this study, a new strong form method, the elemental differential method, has been proposed to simulate the coupled response of the piezoelectric structures subjected to

Conclusions
In this study, a new strong form method, the elemental differential method, has been proposed to simulate the coupled response of the piezoelectric structures subjected to combined thermal-electric-mechanical loadings. The main contribution of this method is to derive the analytical formulas for combining the governing equations of the piezoelectric problem and the heat conduction problem directly. Based on the proposed method, the thermal conduction and thermal, mechanical problems for the thermal-electric-mechanical problem can be resolved with the same strategy. Meanwhile, the pyroelectric effect considering the relationship between thermal and electrical fields is directly addressed.
Several numerical examples, including the piezoelectric composites, are utilized to evaluate the effectiveness of the proposed method. From the comparison of the results calculated by different methods based on the standard finite element framework, it can be found that the proposed method can extract the benefits from the point collocation methods (PCM) and standard finite element methods. Furthermore, the proposed method's final equations about elastic and electric field variables could be constructed from the isoparametric mesh grid without any domain integration. It also can be found that the system of equations for the heat transfer and thermal-elastic response of the piezoelectric problems can be generated in a unified way. At the same time, the effect of thermal expansion can be directly handled. Furthermore, in the proposed method, the pyroelectric effects of the piezoelectric structure can be directly substituted into the system of equations through the collocation technique without any domain integrations. This provides a convenient technique for future research on the pyroelectric effects in actual engineering applications.  Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.