# Hybrid Fundamental Solution Based Finite Element Method: Theory and Applications

Authors: Changyong Cao; Qing-Hua Qin
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

An overview on the development of hybrid fundamental solution based finite element method (HFS-FEM) and its application in engineering problems is presented in this paper. The framework and formulations of HFS-FEM for potential problem, plane elasticity, three-dimensional elasticity, thermoelasticity, anisotropic elasticity, and plane piezoelectricity are presented. In this method, two independent assumed fields (intraelement filed and auxiliary frame field) are employed. The formulations for all cases are derived from the modified variational functionals and the fundamental solutions to a given problem. Generation of elemental stiffness equations from the modified variational principle is also described. Typical numerical examples are given to demonstrate the validity and performance of the HFS-FEM. Finally, a brief summary of the approach is provided and future trends in this field are identified.


---

## Body

## 1. Introduction

A novel hybrid finite element formulation, called the hybrid fundamental solution based FEM (HFS-FEM), was recently developed based on the framework of hybrid Trefftz finite element method (HT-FEM) and the idea of the method of fundamental solution (MFS) [1–5]. In this method, two independent assumed fields (intraelement filed and auxiliary frame field) are employed and the domain integrals in the variational functional can be directly converted to boundary integrals without any appreciable increase in computational effort as in HT-FEM [6–8]. It should be mentioned that the intraelement field of HFS-FEM is approximated by the linear combination of fundamental solutions analytically satisfying the related governing equation, instead of T-complete functions as in HT-FEM. The resulting system of equations from the modified variational functional is expressed in terms of symmetric stiffness matrix and nodal displacements only, which is easy to implement into the standard FEM. It is noted that no singular integrals are involved in the HFS-FEM by locating the source point outside the element of interest and do not overlap with field point during the computation [9].The HFS-FEM mentioned above inherits all the advantages of HT-FEM over the traditional FEM and the boundary element method (BEM), namely, domain decomposition and boundary integral expressions, while avoiding the major weaknesses of BEM [10–12], that is, singular element boundary integral and loss of symmetry and sparsity [13]. The employment of two independent fields also makes the HFS-FEM easier to generate arbitrary polygonal or even curve-sided elements. It also obviates the difficulties that occur in HT-FEM [14, 15] in deriving T-complete functions for certain complex or new physical problems [16]. The HFS-FEM has simpler expressions of interpolation functions for intraelement fields (fundamental solutions) and avoids the coordinate transformation procedure required in the HT-FEM to keep the matrix inversion stable. Moreover, this approach also has the potential to achieve high accuracy using coarse meshes of high-degree elements, to enhance insensitivity to mesh distortion, to give great liberty in element shape, and to accurately represent various local effects (such as hole, crack, and inclusions) without troublesome mesh adjustment [17–20]. Additionally, HFS-FEM makes it possible for a more flexible element material definition which is important in dealing with multimaterial problems, rather than the material definition being the same in the entire domain in BEM. However, we noticed that there are also some limitations of HFS-FEM, for example, determining the positions of source points used for approximation interpolations. It is also known that fundamental solution based approximations can perform remarkably well in smooth problems but tend to deteriorate when high-gradient stress fields are presented.This paper is organized as follows: in Section2, the basic idea and formulations of the HFS-FEM are presented through a simple potential problem. Then, plane elasticity problems are described in Section 3. Section 4 extends the 2D formulations of the HFS-FEM to general three-dimensional (3D) elasticity problems. The method of particular solution and radial basis function approximation are shown to deal with body force in this Section. In Section 5 we extend the HFS-FEM to thermoelastic problems with arbitrary body force and temperature change. In Section 6, the HFS-FEM for 2D anisotropic elastic materials is described based on the powerful Stroh formalism. Plane piezoelectric problem is discussed in Section 7. Finally, typical numerical examples are presented in Section 8 to illustrate applications and performance of the HFS-FEM. Concluding remarks and future development are discussed at the end of this paper.


## 2. Potential Problems

### 2.1. Basic Equations of Potential Problems

The Laplace equation of a well-posed potential problem (e.g., heat conduction) in a general plane domainΩ can be expressed as [21, 22](1)∇2ux=0∀x∈Ωwith the boundary conditions(2)u=u-on⁡Γu,(3)q=u,ini=q-on⁡Γq,where u is the unknown field variable and q represents the boundary flux, ni is the ith component of outward normal vector to the boundary Γ=Γu∪Γq, and u- and q- are specified functions on the related boundaries, respectively. The space derivatives are indicated by a subscript comma, that is, u,i=∂u/∂xi, and the subscript index i takes values (1,2) for two-dimensional and (1,2,3) for three-dimensional problems. Additionally, the repeated subscript indices imply summation convention.For convenience, (3) can be rewritten in matrix form as(4)q=Au,1u,2=q-with A=[n1n2].


### 2.2. Assumed Independent Fields

In this section, the procedure for developing a hybrid finite element model with fundamental solution as interior trial function is described based on the boundary value problem defined by (1)–(3). Similar to the conventional FEM and HT-FEM, the domain under consideration is divided into a series of elements [15, 16, 21, 23–30]. In each element, two independent fields are assumed in the way as described in [31] and are given in Section 2.2.


#### 2.2.1. Intraelement Field

Similar to the method of fundamental solution (MFS) in removing singularities of fundamental solution, for a particular elemente occupying subdomain Ωe, we assume that the field variable defined in the element domain is extracted from a linear combination of fundamental solutions centered at different source points located outside the element (see Figure 1):(5)uex=∑j=1nsNex,yjcej=Nexcekkkkkkkikkkk∀x∈Ωe,yj∉Ωe,where cej is undetermined coefficients, ns is the number of virtual sources, and Ne(x,yj) is the fundamental solution to the partial differential equation:(6)∇2Nex,y+δx,y=0∀x,y∈R2as(7)Nex,y=-12πln⁡⁡rx,y.Evidently, (5) analytically satisfies (1) due to the solution property of Ne(x,yj).Figure 1
Intraelement field and frame field of a HFS-FEM element for 2D potential problems.
(a)
4-node 2D element
(b)
8-node 2D elementIn implementation, the number of source points is taken to be the same as the number of element nodes, which is free of spurious energy modes and can keep the stiffness equations in full rank, as indicated in [21]. The source point ysj(j=1,2,…,ns) can be generated by means of the method employed in the MFS [32–35](8)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 1). Determination of γ was discussed in [31, 36] and γ=5–10 is usually used in practice.The corresponding outward normal derivative ofue on Γe is(9)qe=∂ue∂n=Qece,where(10)Qe=∂Ne∂n=ATewith(11)Te=∂Ne∂x1∂Ne∂x2T.


#### 2.2.2. Auxiliary Frame Field

In order to enforce the conformity on the field variableu, for instance, ue=uf on Γe∩Γf of any two neighboring elements e and f, an auxiliary interelement frame field u~ is used and expressed in terms of the same degrees of freedom (DOF) d as those used in the conventional finite elements. In this case, u~ is confined to the whole element boundary as(12)u~ex=N~exdewhich is independently assumed along the element boundary in terms of nodal DOF de, where N~ex represents the conventional finite element interpolating functions. For example, a simple interpolation of the frame field on a side with three nodes of a particular element can be given in the form(13)u~=N~1u1+N~2u2+N~3u3,where N~i(i=1,2,3) stands for shape functions in terms of natural coordinate ξ defined in Figure 2.Figure 2
Typical quadratic interpolation for frame field.


### 2.3. Modified Variational Principle

For the boundary value problem defined in (1)–(3) and (5), since the stationary conditions of the traditional potential or complementary variational functional cannot guarantee the interelement continuity condition required in the proposed HFS FE model, as in the HT FEM [21, 26], a variational functional corresponding to the new trial functions should be constructed to assure the additional continuity across the common boundaries ΓIef between intraelement fields of element “e” and element “f” (see Figure 3) [36, 37]:(14)kkiue=ufconformityqe+qf=0(reciprocity)kkikkkkon⁡ΓIef=Γe∩Γf.Figure 3
Illustration of continuity between two adjacent elements “e” and “f.”A modified variational functional is developed as follows:(15)Πm=∑eΠme=∑eΠe+∫Γeu~-uqdΓ,where(16)Πe=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓin which the governing equation (1) is assumed to be satisfied, apriori, for deriving the HFS FE model. The boundary Γe of a particular element consists of the following parts:(17)Γe=Γue∪Γqe∪ΓIe,where ΓIe represents the interelement boundary of the element “e” shown in Figure 3.To show that the stationary condition of the functional (15) leads to the governing equation (Euler equation), boundary conditions, and continuity conditions, invoking (16) and (15) gives the following functional for the problem domain:(18)Πme=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓfrom which the first-order variational yields(19)δΠme=∫Ωeu,iδu,idΩ-∫Γqeq-δu~dΓ+∫Γeδu~-δuqdΓ+∫Γeu~-uδqdΓ.Using divergence theorem(20)∫Ωf,ih,idΩ=∫Γhf,inidΓ-∫Ωh∇2fdΩwe can obtain (21)δΠme=∫Ωeu,iiδudΩ-∫Γqeq--qδu~dΓ+∫Γueqδu~dΓ+∫ΓIeqδu~dΓ+∫Γeu~-uδqdΓ.For the displacement-based method, the potential conformity should be satisfied in advance(22)δu~=0on⁡Γue∵u~=u-δu~e=δu~fon⁡ΓIef∵u~e=u~fthen (21) can be rewritten as(23)δΠme=∫Ωeu,iiδudΩ-∫Γqeq--qδu~dΓ+∫ΓIeqδu~dΓ+∫Γeu~-uδqdΓ.The Euler equation and boundary conditions can be obtained as(24)u,ii=0inΩeq=q-on⁡Γqeu~=uon⁡Γeusing the stationary condition δΠme=0.As for the continuous requirement between two adjacent elements “e” and “f” given in (14), we can obtain it in the following way. When assembling elements “e” and “f,” we have(25)δΠm(e+f)=∫Ωe+fu,iiδudΩ-∫Γqe+Γqfq--qδu~dΓ+∫Γeu~-uδqdΓ+∫Γfu~-uδqdΓ+∫ΓIefqe+qfδu~efdΓ+⋯from which the vanishing variation of Πm(e+f) leads to the reciprocity condition qe+qf=0 on the interelement boundary ΓIef.If the following expression:(26)∫Γqδqδu~ds-∑e∫ΓIeδqeδu~eds+∫Γeδqeδu~e-uedsis uniformly positive (or negative) in the neighborhood of {u}0, where the displacement {u}0 has such a value that Πm({u}0)=(Πm)0 and where (Πm)0 stands for the stationary value of Πm, we have (27)Πm≥Πm0orΠm≤Πm0in which the relation that {u~}e={u~}f is identical on Γe∩Γf has been used. This is due to the definition in (14) in Section 2.3.Proof.
For the proof of the theorem on the existence of extremum, we may complete it by the so-called “second variational approach” [38]. In doing this, performing variation of δΠm and using the constrained conditions (26), we find(28)δ2Πm=∫Γqδqδu~ds-∑e∫ΓIeδqeδu~eds+∫Γeδqeδu~e-ueds.Therefore the theorem has been proved from the sufficient condition of the existence of a local extreme of a functional [38].


### 2.4. Element Stiffness Equation

With the intraelement field and frame field independently defined in a particular element (see Figure1), we can generate element stiffness equation by the variational functional derived in Section 2.3. Following the approach described in [21], the variational functional Πe corresponding to a particular element e of the present problem can be written as(29)Πe=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓ.Appling the divergence theorem (20) to the functional (29), we have the final functional for the HFS-FE model(30)Πe=12∫ΓequdΓ+∫Ωeuk∇2udΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓ=-12∫ΓequdΓ-∫Γqeq-u~dΓ+∫Γequ~dΓ.Then, substituting (5), (9), and (12) into the functional (30) produces(31)Πe=-12ceTHece-deTge+ceTGedein which(32)He=∫ΓeQeTNedΓ=∫ΓeNeTQedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓqeN~eTq-dΓ.The symmetry of He is obvious from the scalar definition (31) of variational functional Πe.To enforce interelement continuity on the common element boundary, the unknown vectorce should be expressed in terms of nodal DOF de. The minimization of the functional Πe with respect to ce and de, respectively, yields(33)∂Πe∂ceT=-Hece+Gede=0,∂Πe∂deT=GeTce-ge=0,from which the optional relationship between ce and de and the stiffness equation can be produced(34)ce=He-1Gede,Kede=ge,where Ke=GeTHe-1Ge stands for the element stiffness matrix.


### 2.5. Numerical Integral forH and G Matrix

Generally, it is difficult to obtain the analytical expression of the integral in (32) and numerical integration along the element boundary is required. Herein the widely used Gaussian integration is employed [22].For theHe matrix, one can express it as(35)He=∫ΓeQeTNedΓ=∫ΓeFxdΓby introducing the matrix function (36)Fx=Fijxm×m=QeTNe.Equation (36) can be further rewritten as(37)Hij=∫ΓeFijxdΓ=∑l=1ne∫ΓelFijxdΓ,where(38)dΓ=dx12+dx22=dx1dξ2+dx2dξ2dξ=Jdξand J is the Jacobean expressed as (39)J=dx1dξ2+dx2dξ2,where(40)dx1dξ,dx2dξT=∑i=1nodNiξdξx1ix2i.Thus, the Gaussian numerical integration for H matrix can be calculated by(41)Hij=∑l=1ne∫-1+1FijxξJξdξ≈∑l=1ne∑k=1npwkFijxξkJξk,where ne is the number of edges of the element and np is the Gaussian sampling points employed in the Gaussian numerical integration. Similarly, we can calculate the Ge matrix using (42)Gij=∑l=1ne∫-11F~ijxξJξdξ≈∑l=1ne∑k=1npwkF~ijxξkJξk.The calculation of vectorge in (32) is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program. Besides, the flux is directly computed from (9). The boundary DOF can be directly computed from (12) while the unknown variable at interior points of the element can be determined from (5) plus the recovered rigid-body modes in each element, which is discussed in the following section.


### 2.6. Recovery of Rigid-Body Motion

Considering the physical definition of the fundamental solution, it is necessary to recover the missing rigid-body motion modes from the above results. Following the method presented in [21], the missing rigid-body motion can be recovered by writing the internal potential field of a particular element e as(43)ue=Nece+c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(44)∑i=1nNece+c0-u~e2Nece+c0-u~e2nodei=min⁡which finally gives(45)c0=1n∑i=1nΔueiin which Δuei=u~e-Necenodei and n is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector ce can be evaluated from (34), and then c0 is evaluated from (45). Finally, the potential field u at any internal point in an element can be obtained by means of (43).


## 2.1. Basic Equations of Potential Problems

The Laplace equation of a well-posed potential problem (e.g., heat conduction) in a general plane domainΩ can be expressed as [21, 22](1)∇2ux=0∀x∈Ωwith the boundary conditions(2)u=u-on⁡Γu,(3)q=u,ini=q-on⁡Γq,where u is the unknown field variable and q represents the boundary flux, ni is the ith component of outward normal vector to the boundary Γ=Γu∪Γq, and u- and q- are specified functions on the related boundaries, respectively. The space derivatives are indicated by a subscript comma, that is, u,i=∂u/∂xi, and the subscript index i takes values (1,2) for two-dimensional and (1,2,3) for three-dimensional problems. Additionally, the repeated subscript indices imply summation convention.For convenience, (3) can be rewritten in matrix form as(4)q=Au,1u,2=q-with A=[n1n2].


## 2.2. Assumed Independent Fields

In this section, the procedure for developing a hybrid finite element model with fundamental solution as interior trial function is described based on the boundary value problem defined by (1)–(3). Similar to the conventional FEM and HT-FEM, the domain under consideration is divided into a series of elements [15, 16, 21, 23–30]. In each element, two independent fields are assumed in the way as described in [31] and are given in Section 2.2.


### 2.2.1. Intraelement Field

Similar to the method of fundamental solution (MFS) in removing singularities of fundamental solution, for a particular elemente occupying subdomain Ωe, we assume that the field variable defined in the element domain is extracted from a linear combination of fundamental solutions centered at different source points located outside the element (see Figure 1):(5)uex=∑j=1nsNex,yjcej=Nexcekkkkkkkikkkk∀x∈Ωe,yj∉Ωe,where cej is undetermined coefficients, ns is the number of virtual sources, and Ne(x,yj) is the fundamental solution to the partial differential equation:(6)∇2Nex,y+δx,y=0∀x,y∈R2as(7)Nex,y=-12πln⁡⁡rx,y.Evidently, (5) analytically satisfies (1) due to the solution property of Ne(x,yj).Figure 1
Intraelement field and frame field of a HFS-FEM element for 2D potential problems.
(a)
4-node 2D element
(b)
8-node 2D elementIn implementation, the number of source points is taken to be the same as the number of element nodes, which is free of spurious energy modes and can keep the stiffness equations in full rank, as indicated in [21]. The source point ysj(j=1,2,…,ns) can be generated by means of the method employed in the MFS [32–35](8)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 1). Determination of γ was discussed in [31, 36] and γ=5–10 is usually used in practice.The corresponding outward normal derivative ofue on Γe is(9)qe=∂ue∂n=Qece,where(10)Qe=∂Ne∂n=ATewith(11)Te=∂Ne∂x1∂Ne∂x2T.


### 2.2.2. Auxiliary Frame Field

In order to enforce the conformity on the field variableu, for instance, ue=uf on Γe∩Γf of any two neighboring elements e and f, an auxiliary interelement frame field u~ is used and expressed in terms of the same degrees of freedom (DOF) d as those used in the conventional finite elements. In this case, u~ is confined to the whole element boundary as(12)u~ex=N~exdewhich is independently assumed along the element boundary in terms of nodal DOF de, where N~ex represents the conventional finite element interpolating functions. For example, a simple interpolation of the frame field on a side with three nodes of a particular element can be given in the form(13)u~=N~1u1+N~2u2+N~3u3,where N~i(i=1,2,3) stands for shape functions in terms of natural coordinate ξ defined in Figure 2.Figure 2
Typical quadratic interpolation for frame field.


## 2.2.1. Intraelement Field

Similar to the method of fundamental solution (MFS) in removing singularities of fundamental solution, for a particular elemente occupying subdomain Ωe, we assume that the field variable defined in the element domain is extracted from a linear combination of fundamental solutions centered at different source points located outside the element (see Figure 1):(5)uex=∑j=1nsNex,yjcej=Nexcekkkkkkkikkkk∀x∈Ωe,yj∉Ωe,where cej is undetermined coefficients, ns is the number of virtual sources, and Ne(x,yj) is the fundamental solution to the partial differential equation:(6)∇2Nex,y+δx,y=0∀x,y∈R2as(7)Nex,y=-12πln⁡⁡rx,y.Evidently, (5) analytically satisfies (1) due to the solution property of Ne(x,yj).Figure 1
Intraelement field and frame field of a HFS-FEM element for 2D potential problems.
(a)
4-node 2D element
(b)
8-node 2D elementIn implementation, the number of source points is taken to be the same as the number of element nodes, which is free of spurious energy modes and can keep the stiffness equations in full rank, as indicated in [21]. The source point ysj(j=1,2,…,ns) can be generated by means of the method employed in the MFS [32–35](8)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 1). Determination of γ was discussed in [31, 36] and γ=5–10 is usually used in practice.The corresponding outward normal derivative ofue on Γe is(9)qe=∂ue∂n=Qece,where(10)Qe=∂Ne∂n=ATewith(11)Te=∂Ne∂x1∂Ne∂x2T.


## 2.2.2. Auxiliary Frame Field

In order to enforce the conformity on the field variableu, for instance, ue=uf on Γe∩Γf of any two neighboring elements e and f, an auxiliary interelement frame field u~ is used and expressed in terms of the same degrees of freedom (DOF) d as those used in the conventional finite elements. In this case, u~ is confined to the whole element boundary as(12)u~ex=N~exdewhich is independently assumed along the element boundary in terms of nodal DOF de, where N~ex represents the conventional finite element interpolating functions. For example, a simple interpolation of the frame field on a side with three nodes of a particular element can be given in the form(13)u~=N~1u1+N~2u2+N~3u3,where N~i(i=1,2,3) stands for shape functions in terms of natural coordinate ξ defined in Figure 2.Figure 2
Typical quadratic interpolation for frame field.


## 2.3. Modified Variational Principle

For the boundary value problem defined in (1)–(3) and (5), since the stationary conditions of the traditional potential or complementary variational functional cannot guarantee the interelement continuity condition required in the proposed HFS FE model, as in the HT FEM [21, 26], a variational functional corresponding to the new trial functions should be constructed to assure the additional continuity across the common boundaries ΓIef between intraelement fields of element “e” and element “f” (see Figure 3) [36, 37]:(14)kkiue=ufconformityqe+qf=0(reciprocity)kkikkkkon⁡ΓIef=Γe∩Γf.Figure 3
Illustration of continuity between two adjacent elements “e” and “f.”A modified variational functional is developed as follows:(15)Πm=∑eΠme=∑eΠe+∫Γeu~-uqdΓ,where(16)Πe=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓin which the governing equation (1) is assumed to be satisfied, apriori, for deriving the HFS FE model. The boundary Γe of a particular element consists of the following parts:(17)Γe=Γue∪Γqe∪ΓIe,where ΓIe represents the interelement boundary of the element “e” shown in Figure 3.To show that the stationary condition of the functional (15) leads to the governing equation (Euler equation), boundary conditions, and continuity conditions, invoking (16) and (15) gives the following functional for the problem domain:(18)Πme=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓfrom which the first-order variational yields(19)δΠme=∫Ωeu,iδu,idΩ-∫Γqeq-δu~dΓ+∫Γeδu~-δuqdΓ+∫Γeu~-uδqdΓ.Using divergence theorem(20)∫Ωf,ih,idΩ=∫Γhf,inidΓ-∫Ωh∇2fdΩwe can obtain (21)δΠme=∫Ωeu,iiδudΩ-∫Γqeq--qδu~dΓ+∫Γueqδu~dΓ+∫ΓIeqδu~dΓ+∫Γeu~-uδqdΓ.For the displacement-based method, the potential conformity should be satisfied in advance(22)δu~=0on⁡Γue∵u~=u-δu~e=δu~fon⁡ΓIef∵u~e=u~fthen (21) can be rewritten as(23)δΠme=∫Ωeu,iiδudΩ-∫Γqeq--qδu~dΓ+∫ΓIeqδu~dΓ+∫Γeu~-uδqdΓ.The Euler equation and boundary conditions can be obtained as(24)u,ii=0inΩeq=q-on⁡Γqeu~=uon⁡Γeusing the stationary condition δΠme=0.As for the continuous requirement between two adjacent elements “e” and “f” given in (14), we can obtain it in the following way. When assembling elements “e” and “f,” we have(25)δΠm(e+f)=∫Ωe+fu,iiδudΩ-∫Γqe+Γqfq--qδu~dΓ+∫Γeu~-uδqdΓ+∫Γfu~-uδqdΓ+∫ΓIefqe+qfδu~efdΓ+⋯from which the vanishing variation of Πm(e+f) leads to the reciprocity condition qe+qf=0 on the interelement boundary ΓIef.If the following expression:(26)∫Γqδqδu~ds-∑e∫ΓIeδqeδu~eds+∫Γeδqeδu~e-uedsis uniformly positive (or negative) in the neighborhood of {u}0, where the displacement {u}0 has such a value that Πm({u}0)=(Πm)0 and where (Πm)0 stands for the stationary value of Πm, we have (27)Πm≥Πm0orΠm≤Πm0in which the relation that {u~}e={u~}f is identical on Γe∩Γf has been used. This is due to the definition in (14) in Section 2.3.Proof.
For the proof of the theorem on the existence of extremum, we may complete it by the so-called “second variational approach” [38]. In doing this, performing variation of δΠm and using the constrained conditions (26), we find(28)δ2Πm=∫Γqδqδu~ds-∑e∫ΓIeδqeδu~eds+∫Γeδqeδu~e-ueds.Therefore the theorem has been proved from the sufficient condition of the existence of a local extreme of a functional [38].


## 2.4. Element Stiffness Equation

With the intraelement field and frame field independently defined in a particular element (see Figure1), we can generate element stiffness equation by the variational functional derived in Section 2.3. Following the approach described in [21], the variational functional Πe corresponding to a particular element e of the present problem can be written as(29)Πe=12∫Ωeu,iu,idΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓ.Appling the divergence theorem (20) to the functional (29), we have the final functional for the HFS-FE model(30)Πe=12∫ΓequdΓ+∫Ωeuk∇2udΩ-∫Γqeq-u~dΓ+∫Γequ~-udΓ=-12∫ΓequdΓ-∫Γqeq-u~dΓ+∫Γequ~dΓ.Then, substituting (5), (9), and (12) into the functional (30) produces(31)Πe=-12ceTHece-deTge+ceTGedein which(32)He=∫ΓeQeTNedΓ=∫ΓeNeTQedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓqeN~eTq-dΓ.The symmetry of He is obvious from the scalar definition (31) of variational functional Πe.To enforce interelement continuity on the common element boundary, the unknown vectorce should be expressed in terms of nodal DOF de. The minimization of the functional Πe with respect to ce and de, respectively, yields(33)∂Πe∂ceT=-Hece+Gede=0,∂Πe∂deT=GeTce-ge=0,from which the optional relationship between ce and de and the stiffness equation can be produced(34)ce=He-1Gede,Kede=ge,where Ke=GeTHe-1Ge stands for the element stiffness matrix.


## 2.5. Numerical Integral forH and G Matrix

Generally, it is difficult to obtain the analytical expression of the integral in (32) and numerical integration along the element boundary is required. Herein the widely used Gaussian integration is employed [22].For theHe matrix, one can express it as(35)He=∫ΓeQeTNedΓ=∫ΓeFxdΓby introducing the matrix function (36)Fx=Fijxm×m=QeTNe.Equation (36) can be further rewritten as(37)Hij=∫ΓeFijxdΓ=∑l=1ne∫ΓelFijxdΓ,where(38)dΓ=dx12+dx22=dx1dξ2+dx2dξ2dξ=Jdξand J is the Jacobean expressed as (39)J=dx1dξ2+dx2dξ2,where(40)dx1dξ,dx2dξT=∑i=1nodNiξdξx1ix2i.Thus, the Gaussian numerical integration for H matrix can be calculated by(41)Hij=∑l=1ne∫-1+1FijxξJξdξ≈∑l=1ne∑k=1npwkFijxξkJξk,where ne is the number of edges of the element and np is the Gaussian sampling points employed in the Gaussian numerical integration. Similarly, we can calculate the Ge matrix using (42)Gij=∑l=1ne∫-11F~ijxξJξdξ≈∑l=1ne∑k=1npwkF~ijxξkJξk.The calculation of vectorge in (32) is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program. Besides, the flux is directly computed from (9). The boundary DOF can be directly computed from (12) while the unknown variable at interior points of the element can be determined from (5) plus the recovered rigid-body modes in each element, which is discussed in the following section.


## 2.6. Recovery of Rigid-Body Motion

Considering the physical definition of the fundamental solution, it is necessary to recover the missing rigid-body motion modes from the above results. Following the method presented in [21], the missing rigid-body motion can be recovered by writing the internal potential field of a particular element e as(43)ue=Nece+c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(44)∑i=1nNece+c0-u~e2Nece+c0-u~e2nodei=min⁡which finally gives(45)c0=1n∑i=1nΔueiin which Δuei=u~e-Necenodei and n is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector ce can be evaluated from (34), and then c0 is evaluated from (45). Finally, the potential field u at any internal point in an element can be obtained by means of (43).


## 3. Plane Elasticity Problems

### 3.1. Linear Theory of Plane Elasticity

In linear elastic theory, the strain displacement relations can be used and equilibrium equations refer to the undeformed geometry [39]. In the rectangular Cartesian coordinates (X1,X2), the governing equations of a plane elastic body can be expressed as(46)σij,j=bi,i,j=1,2.If written as matrix form, it can be presented as(47)Lσ=b,where σ=[σ11σ22σ12]T is a stress vector, b=[b1,b2]T is a body force vector, and the differential operator matrix L is given as(48)L=∂∂x10∂∂x20∂∂x2∂∂x1,(49)ε=LTu,where ε=[ε11ε22ε12]T is a strain vector and u=[u1,u2]T is a displacement vector.The constitutive equations for the linear elasticity are given in matrix form as(50)σ=Dε,where D is the material coefficient matrix with constant components for isotropic materials, which can be expressed as follows: (51)D=λ~+2Gλ~0λ~λ~+2G000G,where(52)λ~=2ν~1-2ν~G,G=E21+ν,ν~=νforplanestrainν1+νforplanestress.The two different kinds of boundary conditions can be expressed as(53)u=u-on⁡Γu,t=Aσ=t-on⁡Γt,where t=[t1t2]T denotes the traction vector and A is a transformation matrix related to the direction cosine of the outward normal (54)A=n10n20n2n1.Substituting (49) and (50) into (47) yields the well-known Navier partial differential equations in terms of displacements:(55)LDLTu=b.


### 3.2. Assumed Independent Fields

For elasticity problem, two different assumed fields are employed as in potential problems: intraelement and frame field [1, 31, 36]. The intraelement continuity is enforced on nonconforming internal displacement field chosen as the fundamental solution of the problem [36]. The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem of interest:(56)ux=u1xu2x=∑j=1nsu11*x,ysju12*x,ysju21*x,ysju22*x,ysjc1jc2j=Necex∈Ωe,ysj∉Ωe,where ns is again the number of source points outside the element domain, which is equal to the number of nodes of an element in the present research based on the generation approach of the source points [31]. The vector ce and the fundamental solution matrix Ne are now in the form (57)Ne=u11*x,ys1u12*x,ys1⋯u11*x,ysnsu12*x,ysnsu21*x,ys1u22*x,ys1⋯u21*x,ysnsu22*x,ysns,ce=c11c21⋯c1nc2nTin which x and ysj are, respectively, the field point and source point in the local coordinate system (X1,X2). The components uij*(x,ysj) are the fundamental solution, that is, induced displacement component in i-direction at the field point x due to a unit point load applied in j-direction at the source point ysj, which are given by [40, 41](58)uij*x,ysj=-18π1-νG3-4νδijln⁡⁡r-r,ir,j,where ri=xi-xis,r=r12+r22. The virtual source points for elasticity problems are generated in the same manner as that in potential problems described in Section 2.With the assumption of intraelement field in (56), the corresponding stress fields can be obtained by the constitutive equation (50):(59)σx=σ11σ22σ12T=Tece,where (60)Te=σ111*x,y1σ211*x,y1⋯σ111*x,ynsσ211*x,ynsσ122*x,y1σ222*x,y1⋯σ122*x,ynsσ222*x,ynsσ112*x,y1σ212*x,y1⋯σ112*x,ynsσ212*x,yns.As a consequence, the traction is written as (61)t1t2=nσ=Qecein which (62)Qe=nTe,n=n10n20n2n1.The components σijk*(x,y) for plane strain problems are given as (63)σijk*x,y=-14π1-νr·1-2νr,kδij+r,jδki-r,iδjk+2r,ir,jr,k.The unknownce in (56) is calculated by a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as that in conventional FEM (see Figure 1). This means that, in the HFS-FEM, a conforming displacement field should be independently defined on the element boundary to enforce the field continuity between elements and also to link the unknown c with nodal displacement de. Thus, the frame is defined as (64)u~x=u~1u~2=N~1N~2de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N~e and de can be expressed as (65)N~e=0⋯0N~10N~20N~300⋯00⋯0︸40N~10N~20N~30⋯0︸6,de=u11u21u12u22⋯u18u28T,where N~1,N~2, and N~3 can be expressed by natural coordinate as(66)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkξ∈-1,1.


### 3.3. Modified Functional for the Hybrid FEM

As in Section2, HFS-FE formulation for a plane elastic problem can also be established by the variational approach [36]. In the absence of body forces, the hybrid functional Πme used for deriving the present HFS-FEM can be constructed as [22](67)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ,where u~i and ui are the intraelement displacement field defined within the element and the frame displacement field defined on the element boundary, respectively. Ωe and Γe are the element domain and element boundary, respectively.  Γt,Γu, and ΓI stand, respectively, for the specified traction boundary, specified displacement boundary, and interelement boundary (Γe=Γt+Γu+ΓI). Compared to the functional employed in the conventional FEM, the present variational functional is constructed by adding a hybrid integral term related to the intraelement and element frame displacement fields to guarantee the satisfaction of displacement and traction continuity conditions on the common boundary of two adjacent elements.By applying the Gaussian theorem, (67) can be simplified to(68)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement field, we have the following expression for HFS finite element model:(69)Πme=12∫ΓetiuidΓ-∫Γtt-iu~idΓ+∫Γeti(u~i-ui)dΓ=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.The variational functional in (69) contains boundary integrals only and will be used to derive HFS-FEM formulation for the plane isotropic elastic problem.


### 3.4. Element Stiffness Matrix

As in Section2, the element stiffness equation can be generated by setting δΠme=0. Substituting (56), (64), and (61) into the functional of (69), we have(70)Πme=-12ceTHece+ceTGede-deTge,where (71)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.To enforce interelement continuity on the common element boundary, the unknown vector ce should be expressed in terms of nodal DOF de. The stationary condition of the functional Πme with respect to ce and de yields, respectively, (33) and (34).


### 3.5. Recovery of Rigid-Body Motion

For the same reason stated in Section2.6, it is necessary to reintroduce the discarded rigid-body motion terms after we have obtained the internal field of an element. The least squares method can be employed for this purpose and the missing terms can be recovered easily by setting for the augmented internal field [22](72)ue=Nece+10x201-x1c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(73)∑i=1nu1i-u~1i2+u2i-u~2i2=min⁡which finally gives(74)c0=Re-1re,where(75)Re=∑i=1n10x2i01-x1ix2i-x1ix1i2+x2i2,re=∑i=1nΔue1iΔue2iΔue1ix2i-Δue2ix1iin which Δueji=u~eji-u^eji(j=1,2) and n is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector ce can be evaluated from (56), and then c0 is evaluated from (74). Finally, the displacement field ue at any internal point in an element can be obtained by (72).


## 3.1. Linear Theory of Plane Elasticity

In linear elastic theory, the strain displacement relations can be used and equilibrium equations refer to the undeformed geometry [39]. In the rectangular Cartesian coordinates (X1,X2), the governing equations of a plane elastic body can be expressed as(46)σij,j=bi,i,j=1,2.If written as matrix form, it can be presented as(47)Lσ=b,where σ=[σ11σ22σ12]T is a stress vector, b=[b1,b2]T is a body force vector, and the differential operator matrix L is given as(48)L=∂∂x10∂∂x20∂∂x2∂∂x1,(49)ε=LTu,where ε=[ε11ε22ε12]T is a strain vector and u=[u1,u2]T is a displacement vector.The constitutive equations for the linear elasticity are given in matrix form as(50)σ=Dε,where D is the material coefficient matrix with constant components for isotropic materials, which can be expressed as follows: (51)D=λ~+2Gλ~0λ~λ~+2G000G,where(52)λ~=2ν~1-2ν~G,G=E21+ν,ν~=νforplanestrainν1+νforplanestress.The two different kinds of boundary conditions can be expressed as(53)u=u-on⁡Γu,t=Aσ=t-on⁡Γt,where t=[t1t2]T denotes the traction vector and A is a transformation matrix related to the direction cosine of the outward normal (54)A=n10n20n2n1.Substituting (49) and (50) into (47) yields the well-known Navier partial differential equations in terms of displacements:(55)LDLTu=b.


## 3.2. Assumed Independent Fields

For elasticity problem, two different assumed fields are employed as in potential problems: intraelement and frame field [1, 31, 36]. The intraelement continuity is enforced on nonconforming internal displacement field chosen as the fundamental solution of the problem [36]. The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem of interest:(56)ux=u1xu2x=∑j=1nsu11*x,ysju12*x,ysju21*x,ysju22*x,ysjc1jc2j=Necex∈Ωe,ysj∉Ωe,where ns is again the number of source points outside the element domain, which is equal to the number of nodes of an element in the present research based on the generation approach of the source points [31]. The vector ce and the fundamental solution matrix Ne are now in the form (57)Ne=u11*x,ys1u12*x,ys1⋯u11*x,ysnsu12*x,ysnsu21*x,ys1u22*x,ys1⋯u21*x,ysnsu22*x,ysns,ce=c11c21⋯c1nc2nTin which x and ysj are, respectively, the field point and source point in the local coordinate system (X1,X2). The components uij*(x,ysj) are the fundamental solution, that is, induced displacement component in i-direction at the field point x due to a unit point load applied in j-direction at the source point ysj, which are given by [40, 41](58)uij*x,ysj=-18π1-νG3-4νδijln⁡⁡r-r,ir,j,where ri=xi-xis,r=r12+r22. The virtual source points for elasticity problems are generated in the same manner as that in potential problems described in Section 2.With the assumption of intraelement field in (56), the corresponding stress fields can be obtained by the constitutive equation (50):(59)σx=σ11σ22σ12T=Tece,where (60)Te=σ111*x,y1σ211*x,y1⋯σ111*x,ynsσ211*x,ynsσ122*x,y1σ222*x,y1⋯σ122*x,ynsσ222*x,ynsσ112*x,y1σ212*x,y1⋯σ112*x,ynsσ212*x,yns.As a consequence, the traction is written as (61)t1t2=nσ=Qecein which (62)Qe=nTe,n=n10n20n2n1.The components σijk*(x,y) for plane strain problems are given as (63)σijk*x,y=-14π1-νr·1-2νr,kδij+r,jδki-r,iδjk+2r,ir,jr,k.The unknownce in (56) is calculated by a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as that in conventional FEM (see Figure 1). This means that, in the HFS-FEM, a conforming displacement field should be independently defined on the element boundary to enforce the field continuity between elements and also to link the unknown c with nodal displacement de. Thus, the frame is defined as (64)u~x=u~1u~2=N~1N~2de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N~e and de can be expressed as (65)N~e=0⋯0N~10N~20N~300⋯00⋯0︸40N~10N~20N~30⋯0︸6,de=u11u21u12u22⋯u18u28T,where N~1,N~2, and N~3 can be expressed by natural coordinate as(66)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkξ∈-1,1.


## 3.3. Modified Functional for the Hybrid FEM

As in Section2, HFS-FE formulation for a plane elastic problem can also be established by the variational approach [36]. In the absence of body forces, the hybrid functional Πme used for deriving the present HFS-FEM can be constructed as [22](67)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ,where u~i and ui are the intraelement displacement field defined within the element and the frame displacement field defined on the element boundary, respectively. Ωe and Γe are the element domain and element boundary, respectively.  Γt,Γu, and ΓI stand, respectively, for the specified traction boundary, specified displacement boundary, and interelement boundary (Γe=Γt+Γu+ΓI). Compared to the functional employed in the conventional FEM, the present variational functional is constructed by adding a hybrid integral term related to the intraelement and element frame displacement fields to guarantee the satisfaction of displacement and traction continuity conditions on the common boundary of two adjacent elements.By applying the Gaussian theorem, (67) can be simplified to(68)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement field, we have the following expression for HFS finite element model:(69)Πme=12∫ΓetiuidΓ-∫Γtt-iu~idΓ+∫Γeti(u~i-ui)dΓ=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.The variational functional in (69) contains boundary integrals only and will be used to derive HFS-FEM formulation for the plane isotropic elastic problem.


## 3.4. Element Stiffness Matrix

As in Section2, the element stiffness equation can be generated by setting δΠme=0. Substituting (56), (64), and (61) into the functional of (69), we have(70)Πme=-12ceTHece+ceTGede-deTge,where (71)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.To enforce interelement continuity on the common element boundary, the unknown vector ce should be expressed in terms of nodal DOF de. The stationary condition of the functional Πme with respect to ce and de yields, respectively, (33) and (34).


## 3.5. Recovery of Rigid-Body Motion

For the same reason stated in Section2.6, it is necessary to reintroduce the discarded rigid-body motion terms after we have obtained the internal field of an element. The least squares method can be employed for this purpose and the missing terms can be recovered easily by setting for the augmented internal field [22](72)ue=Nece+10x201-x1c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(73)∑i=1nu1i-u~1i2+u2i-u~2i2=min⁡which finally gives(74)c0=Re-1re,where(75)Re=∑i=1n10x2i01-x1ix2i-x1ix1i2+x2i2,re=∑i=1nΔue1iΔue2iΔue1ix2i-Δue2ix1iin which Δueji=u~eji-u^eji(j=1,2) and n is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector ce can be evaluated from (56), and then c0 is evaluated from (74). Finally, the displacement field ue at any internal point in an element can be obtained by (72).


## 4. Three-Dimensional Elastic Problems

In this section, the HFS-FEM approach is extended to three-dimensional (3D) elastic problem with/without body force. The detailed 3D formulations of HFS-FEM are firstly derived for elastic problems by ignoring body forces, and then a procedure based on the method of particular solution and radial basis function approximation are introduced to deal with the body force [42]. As a consequence, the homogeneous solution is obtained by using the HFS-FEM and the particular solution associated with body force is approximated by using the strong form of basis function interpolation.


### 4.1. Governing Equations and Boundary Conditions

Let (x1,x2,x3) denote the coordinates in a Cartesian coordinate system and consider a finite isotropic body occupying the domain Ω, as shown in Figure 4. The equilibrium equation for this finite body with body force can be expressed as (76)σij,j=-bii,j=1,2,3.Figure 4
Geometrical definitions and boundary conditions for a general 3D solid.The constitutive equations for linear elasticity and the kinematical relation are given as(77)σij=2Gv1-2vδijekk+2Geij,eij=12ui,j+uj,i,where σij is the stress tensor, eij is the strain tensor, and δij is the Kronecker delta. Substituting (77) into (76), the equilibrium equation is rewritten as(78)Gui,jj+G1-2vuj,ji=-bi.For a well-posed boundary value problem, the boundary conditions are prescribed as follows:(79)ui=u-ion⁡Γu,ti=t-ion⁡Γt,where Γu∪Γt=Γ is the boundary of the solution domain Ω,u-i and t-i are the prescribed boundary values. In the following parts, we will present the procedure for handling the body force appearing in (78).


### 4.2. The Method of Particular Solution

The inhomogeneous termbi associated with the body force in (78) can be effectively handled by means of the method of particular solution presented in [22]. In this approach, the displacement ui is decomposed into two parts, a homogeneous solution uih and a particular solution uip(80)ui=uih+uip,where the particular solution uip should satisfy the governing equation: (81)Gui,jjp+G1-2vuj,jip=-biwithout any restriction of boundary condition. However, the homogeneous solution should satisfy(82)Gui,jjh+G1-2vuj,jih=0with the modified boundary conditions(83)uih=u-i-uipon⁡Γu,tih=t-i-tipon⁡Γt.From the above equations it can be seen that once the particular solutionuip is known, the homogeneous solution uih in (82) and (83) can be obtained using HFS-FEM. The final solution can then be given by (80). In the next section, radial basis function approximation is introduced to obtain the particular solution, and the HFS-FEM is presented for solving (82) and (83).


### 4.3. Radial Basis Function Approximation

For body forcebi, it is generally impossible to find an analytical solution which enables us to convert the domain integral into a boundary one. So we must approximate it by a combination of basis (trial) functions or other methods with the HFS-FEM. Here, we use radial basis function (RBF) [33, 43] to interpolate the body force. We assume(84)bi≈∑j=1Nαijφj,where N is the number of interpolation points, φj are the RBFs, and αij are the coefficients to be determined. Subsequently, the particular solution can be approximated by(85)uip=∑j=1NαijΦikj,where Φikj is the approximated particular solution kernel of displacement satisfying (86) below. Once the RBFs are selected, the problem of finding a particular solution is reduced to solve the following equation:(86)GΦil,kk+G1-2νΦkl,ki=-δilφ.To solve this equation, the displacement is expressed in terms of the Galerkin-Papkovich vectors(87)Φik=1-νGFik,mm-12GFmk,mi.Substituting (87) into (86), we can obtain the following biharmonic equation:(88)∇4Fil=-11-νδilφ.Taking the Spline Type RBF φ=r2n-1 as an example, we get the following solutions:(89)Fli=-δli1-vr2n+32n+12n+22n+32n+4,(90)Φli=A0A1δli+A2r,ir,l,where(91)A0=-18G1-vr2n+1n+1n+22n+1,A1=7+4n-4vn+2,A2=-2n+1,and rj represents the Euclidean distance between a field point (x1,x2,x3) and a given point (x1j,x2j,x3j) in the domain of interest. The corresponding particular solution of stresses can be obtained as(92)Slij=GΦli,j+Φlj,i+λδijΦlk,k,where λ=(2v/(1-2v))G. Substituting (90) into (92), we have(93)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,l,where (94)B0=-141-vr2nn+1n+2,B1=3+2n-2vn+2,B2=2vn+2-1,B3=1-2n.


### 4.4. HFS-FEM for Homogeneous Solution

After obtaining the particular solution in Sections4.2 and 4.3, we can determine the modified boundary conditions (83). Finally, we can treat the 3D problem as a homogeneous problem governed by (82) and (83) by using the HFS-FEM presented below. It is clear that once the particular and homogeneous solutions for displacement and stress components at nodal points are determined, the distribution of displacement and stress fields at any point in the domain can be calculated using the element interpolation function. However, for 3D elasticity problems in the absent of body force, that is, bi=0, the procedures in Sections 4.2 and 4.3 will become unnecessary and we can employ the following procedures to find the solution directly.


#### 4.4.1. Assumed Intraelement and Auxiliary Frame Fields

The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem as(95)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as(96)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysnsce=c11c21c31⋯c1nc2nc3nTin which x and ysj are, respectively, the field point and source point in the local coordinate system (x1,x2). The fundamental solution uij*(x,ysj) is given by [40](97)uij*x,ysj=116π1-νGr3-4νδij+r,ir,j,where ri=xi-xis,r=r12+r22+r32, ns is the number of source points. The source point ysj(j=1,2,…,ns) can also be generated by means of the following method [36] as in two-dimensional cases (98)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 5).Figure 5
Intraelement field and frame field of a hexahedron HFS-FEM element for 3D elastic problem (the source points and centroid of 20-node element are omitted in the figure for clarity and clear view, which is similar to that of the 8-node element).
(a)
(b)According to (50) and (49), the corresponding stress fields can be expressed as (99)σx=σ11σ22σ33σ23σ31σ12T=Tece,where (100)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ133*x,y1σ233*x,y1σ333*x,y1⋯σ133*x,ynsσ233*x,ynsσ333*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ212*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by(101)σijk*x,y=-18π1-νr2·1-2νr,kδij+r,jδki-r,iδjk+3r,ir,jr,k.As a consequence, the traction can be written in the form (102)t1t2t3=nσ=Qecein which (103)Qe=nTe,n=n1000n3n20n20n30n100n3n2n10.To link the unknownce and the nodal displacement de, the frame is defined as (104)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the surface 2-3-7-6 of a particular 8-node brick element (see Figure 5) as an example, matrix N~e and vector de can be expressed as (105)N~e=0N-1N-200N-4N-30,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (106)N-i=N~i000N~i000N~i,0=000000000,where N~i (i=1–4) can be expressed by natural coordinate ξ,η∈[-1,1](107)N~1=1+ξ1+η4,N~2=1-ξ1+η4,N~3=1-ξ1-η4,N~4=1+ξ1-η4and ξi,ηi is the natural coordinate of the i-node of the element (Figure 6).Figure 6
Typical linear interpolation for the frame fields of 3D brick elements.


#### 4.4.2. Modified Functional for Hybrid Finite Element Method

In the absence of body forces, the hybrid functionalΠme used for deriving the present HFS-FEM is written as [22](108)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.By applying the Gaussian theorem, (108) can be simplified as(109)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for HFS finite element model:(110)Πme=12∫ΓetiuidΓ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.The functional (110) contains only boundary integrals of the element and will be used to derive HFS-FEM formulation for the three-dimensional elastic problem in the following section.


#### 4.4.3. Element Stiffness Matrix

Substituting (95), (102), and (104) into the functional (110), we have (111)Πme=-12ceTHece+ceTGede-deTge,where (112)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.To enforce interelement continuity on the common element boundary, the unknown vector ce should be expressed in terms of nodal DOF de. The stationary condition of the functional Πme with respect to ce and de yields again, respectively, (33) and (34).


#### 4.4.4. Numerical Integral over Element

Considering a surface of the 3D hexahedron element, as shown in Figure6, the vector normal to the surface can be obtained by (113)vn=vξ×vη=vnxvnyvnz=dxdξdydξdzdξ×dxdηdydηdzdη=dydξdzdη-dydηdzdξdzdξdxdη-dzdηdxdξdxdξdydη-dxdηdydξ,where vξ and vη are the tangential vectors in the ξ-direction and η-direction, respectively, calculated by(114)vξ=dxdξdydξdzdξ=∑i=1nd∂Niξ,η∂ξxiyizi,vη=dxdηdydηdzdη=∑i=1nd∂Niξ,η∂ηxiyizi,where nd is the number of nodes of the surface and xi,yi,zi are the nodal coordinates. Thus the unit normal vector is given by(115)n=vnvn,where(116)Jξ,η=vn=vnx2+vny2+vnz2is the Jacobian of the transformation from Cartesian coordinates (x,y) to natural coordinates (ξ,η).For theH matrix, we introduce the matrix function (117)Fx,y=Fijx,ym×m=QeTNe.Then we can get(118)He=∫ΓeQeTNedΓ=∫ΓeFx,ydΓand we rewrite it to the component form as (119)Hij=∫ΓeFijx,ydS=∑l=1nf∫ΓelFijx,ydS.Using the relationship (120)dS=Jξ,ηdξdηand the Gaussian numerical integration, we can obtain(121)Hij=∑l=1nf∫-11Fijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtFijxξs,ηt,yξs,ηtJξs,ηt,where nf and np are, respectively, the number of surface of the 3D element and the number of Gaussian integral points in each direction of the element surface. Similarly, we can calculate the Ge matrix by(122)Gij=∑l=1nf∫-11F~ijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtF~ijxξs,ηt,yξs,ηtJξs,ηt.It should be mentioned that the calculation of vector ge in equation is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program. Besides, the stress and traction estimations are directly computed from (99) and (100), respectively. The boundary displacements can be directly computed from (104) while the displacements at interior points of element can be determined from (95) plus the recovered rigid-body modes in each element, which is introduced in the following section.


#### 4.4.5. Recovery of Rigid-Body Motion Terms

As in Section2, the least square method is employed to recover the missing terms of rigid-body motions. The missing terms can be recovered by setting for the augmented internal field(123)ue=Nece+1000x3-x2010-x30x1001x2-x10c0and using a least-square procedure to match ueh and u~eh at the nodes of the element boundary(124)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+u3i-u~3i2.The above equation finally yields(125)c0=Re-1re,where(126)Re=∑i=1n1000x3i-x2i010-x3i0x1i001x2i-x1i00-x3ix2ix2i2+x3i2-x1ix2i-x1ix3ix3i0-x1i-x1ix2ix1i2+x3i2-x2ix3i-x2ix1i0-x1ix3i-x2ix3ix1i2+x2i2,re=∑i=1nΔue1iΔue2iΔue3iΔue3ix2i-Δue2ix3iΔue1ix3i-Δue3ix1iΔue2ix1i-Δue1ix2i.


## 4.1. Governing Equations and Boundary Conditions

Let (x1,x2,x3) denote the coordinates in a Cartesian coordinate system and consider a finite isotropic body occupying the domain Ω, as shown in Figure 4. The equilibrium equation for this finite body with body force can be expressed as (76)σij,j=-bii,j=1,2,3.Figure 4
Geometrical definitions and boundary conditions for a general 3D solid.The constitutive equations for linear elasticity and the kinematical relation are given as(77)σij=2Gv1-2vδijekk+2Geij,eij=12ui,j+uj,i,where σij is the stress tensor, eij is the strain tensor, and δij is the Kronecker delta. Substituting (77) into (76), the equilibrium equation is rewritten as(78)Gui,jj+G1-2vuj,ji=-bi.For a well-posed boundary value problem, the boundary conditions are prescribed as follows:(79)ui=u-ion⁡Γu,ti=t-ion⁡Γt,where Γu∪Γt=Γ is the boundary of the solution domain Ω,u-i and t-i are the prescribed boundary values. In the following parts, we will present the procedure for handling the body force appearing in (78).


## 4.2. The Method of Particular Solution

The inhomogeneous termbi associated with the body force in (78) can be effectively handled by means of the method of particular solution presented in [22]. In this approach, the displacement ui is decomposed into two parts, a homogeneous solution uih and a particular solution uip(80)ui=uih+uip,where the particular solution uip should satisfy the governing equation: (81)Gui,jjp+G1-2vuj,jip=-biwithout any restriction of boundary condition. However, the homogeneous solution should satisfy(82)Gui,jjh+G1-2vuj,jih=0with the modified boundary conditions(83)uih=u-i-uipon⁡Γu,tih=t-i-tipon⁡Γt.From the above equations it can be seen that once the particular solutionuip is known, the homogeneous solution uih in (82) and (83) can be obtained using HFS-FEM. The final solution can then be given by (80). In the next section, radial basis function approximation is introduced to obtain the particular solution, and the HFS-FEM is presented for solving (82) and (83).


## 4.3. Radial Basis Function Approximation

For body forcebi, it is generally impossible to find an analytical solution which enables us to convert the domain integral into a boundary one. So we must approximate it by a combination of basis (trial) functions or other methods with the HFS-FEM. Here, we use radial basis function (RBF) [33, 43] to interpolate the body force. We assume(84)bi≈∑j=1Nαijφj,where N is the number of interpolation points, φj are the RBFs, and αij are the coefficients to be determined. Subsequently, the particular solution can be approximated by(85)uip=∑j=1NαijΦikj,where Φikj is the approximated particular solution kernel of displacement satisfying (86) below. Once the RBFs are selected, the problem of finding a particular solution is reduced to solve the following equation:(86)GΦil,kk+G1-2νΦkl,ki=-δilφ.To solve this equation, the displacement is expressed in terms of the Galerkin-Papkovich vectors(87)Φik=1-νGFik,mm-12GFmk,mi.Substituting (87) into (86), we can obtain the following biharmonic equation:(88)∇4Fil=-11-νδilφ.Taking the Spline Type RBF φ=r2n-1 as an example, we get the following solutions:(89)Fli=-δli1-vr2n+32n+12n+22n+32n+4,(90)Φli=A0A1δli+A2r,ir,l,where(91)A0=-18G1-vr2n+1n+1n+22n+1,A1=7+4n-4vn+2,A2=-2n+1,and rj represents the Euclidean distance between a field point (x1,x2,x3) and a given point (x1j,x2j,x3j) in the domain of interest. The corresponding particular solution of stresses can be obtained as(92)Slij=GΦli,j+Φlj,i+λδijΦlk,k,where λ=(2v/(1-2v))G. Substituting (90) into (92), we have(93)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,l,where (94)B0=-141-vr2nn+1n+2,B1=3+2n-2vn+2,B2=2vn+2-1,B3=1-2n.


## 4.4. HFS-FEM for Homogeneous Solution

After obtaining the particular solution in Sections4.2 and 4.3, we can determine the modified boundary conditions (83). Finally, we can treat the 3D problem as a homogeneous problem governed by (82) and (83) by using the HFS-FEM presented below. It is clear that once the particular and homogeneous solutions for displacement and stress components at nodal points are determined, the distribution of displacement and stress fields at any point in the domain can be calculated using the element interpolation function. However, for 3D elasticity problems in the absent of body force, that is, bi=0, the procedures in Sections 4.2 and 4.3 will become unnecessary and we can employ the following procedures to find the solution directly.


### 4.4.1. Assumed Intraelement and Auxiliary Frame Fields

The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem as(95)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as(96)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysnsce=c11c21c31⋯c1nc2nc3nTin which x and ysj are, respectively, the field point and source point in the local coordinate system (x1,x2). The fundamental solution uij*(x,ysj) is given by [40](97)uij*x,ysj=116π1-νGr3-4νδij+r,ir,j,where ri=xi-xis,r=r12+r22+r32, ns is the number of source points. The source point ysj(j=1,2,…,ns) can also be generated by means of the following method [36] as in two-dimensional cases (98)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 5).Figure 5
Intraelement field and frame field of a hexahedron HFS-FEM element for 3D elastic problem (the source points and centroid of 20-node element are omitted in the figure for clarity and clear view, which is similar to that of the 8-node element).
(a)
(b)According to (50) and (49), the corresponding stress fields can be expressed as (99)σx=σ11σ22σ33σ23σ31σ12T=Tece,where (100)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ133*x,y1σ233*x,y1σ333*x,y1⋯σ133*x,ynsσ233*x,ynsσ333*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ212*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by(101)σijk*x,y=-18π1-νr2·1-2νr,kδij+r,jδki-r,iδjk+3r,ir,jr,k.As a consequence, the traction can be written in the form (102)t1t2t3=nσ=Qecein which (103)Qe=nTe,n=n1000n3n20n20n30n100n3n2n10.To link the unknownce and the nodal displacement de, the frame is defined as (104)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the surface 2-3-7-6 of a particular 8-node brick element (see Figure 5) as an example, matrix N~e and vector de can be expressed as (105)N~e=0N-1N-200N-4N-30,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (106)N-i=N~i000N~i000N~i,0=000000000,where N~i (i=1–4) can be expressed by natural coordinate ξ,η∈[-1,1](107)N~1=1+ξ1+η4,N~2=1-ξ1+η4,N~3=1-ξ1-η4,N~4=1+ξ1-η4and ξi,ηi is the natural coordinate of the i-node of the element (Figure 6).Figure 6
Typical linear interpolation for the frame fields of 3D brick elements.


### 4.4.2. Modified Functional for Hybrid Finite Element Method

In the absence of body forces, the hybrid functionalΠme used for deriving the present HFS-FEM is written as [22](108)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.By applying the Gaussian theorem, (108) can be simplified as(109)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for HFS finite element model:(110)Πme=12∫ΓetiuidΓ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.The functional (110) contains only boundary integrals of the element and will be used to derive HFS-FEM formulation for the three-dimensional elastic problem in the following section.


### 4.4.3. Element Stiffness Matrix

Substituting (95), (102), and (104) into the functional (110), we have (111)Πme=-12ceTHece+ceTGede-deTge,where (112)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.To enforce interelement continuity on the common element boundary, the unknown vector ce should be expressed in terms of nodal DOF de. The stationary condition of the functional Πme with respect to ce and de yields again, respectively, (33) and (34).


### 4.4.4. Numerical Integral over Element

Considering a surface of the 3D hexahedron element, as shown in Figure6, the vector normal to the surface can be obtained by (113)vn=vξ×vη=vnxvnyvnz=dxdξdydξdzdξ×dxdηdydηdzdη=dydξdzdη-dydηdzdξdzdξdxdη-dzdηdxdξdxdξdydη-dxdηdydξ,where vξ and vη are the tangential vectors in the ξ-direction and η-direction, respectively, calculated by(114)vξ=dxdξdydξdzdξ=∑i=1nd∂Niξ,η∂ξxiyizi,vη=dxdηdydηdzdη=∑i=1nd∂Niξ,η∂ηxiyizi,where nd is the number of nodes of the surface and xi,yi,zi are the nodal coordinates. Thus the unit normal vector is given by(115)n=vnvn,where(116)Jξ,η=vn=vnx2+vny2+vnz2is the Jacobian of the transformation from Cartesian coordinates (x,y) to natural coordinates (ξ,η).For theH matrix, we introduce the matrix function (117)Fx,y=Fijx,ym×m=QeTNe.Then we can get(118)He=∫ΓeQeTNedΓ=∫ΓeFx,ydΓand we rewrite it to the component form as (119)Hij=∫ΓeFijx,ydS=∑l=1nf∫ΓelFijx,ydS.Using the relationship (120)dS=Jξ,ηdξdηand the Gaussian numerical integration, we can obtain(121)Hij=∑l=1nf∫-11Fijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtFijxξs,ηt,yξs,ηtJξs,ηt,where nf and np are, respectively, the number of surface of the 3D element and the number of Gaussian integral points in each direction of the element surface. Similarly, we can calculate the Ge matrix by(122)Gij=∑l=1nf∫-11F~ijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtF~ijxξs,ηt,yξs,ηtJξs,ηt.It should be mentioned that the calculation of vector ge in equation is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program. Besides, the stress and traction estimations are directly computed from (99) and (100), respectively. The boundary displacements can be directly computed from (104) while the displacements at interior points of element can be determined from (95) plus the recovered rigid-body modes in each element, which is introduced in the following section.


### 4.4.5. Recovery of Rigid-Body Motion Terms

As in Section2, the least square method is employed to recover the missing terms of rigid-body motions. The missing terms can be recovered by setting for the augmented internal field(123)ue=Nece+1000x3-x2010-x30x1001x2-x10c0and using a least-square procedure to match ueh and u~eh at the nodes of the element boundary(124)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+u3i-u~3i2.The above equation finally yields(125)c0=Re-1re,where(126)Re=∑i=1n1000x3i-x2i010-x3i0x1i001x2i-x1i00-x3ix2ix2i2+x3i2-x1ix2i-x1ix3ix3i0-x1i-x1ix2ix1i2+x3i2-x2ix3i-x2ix1i0-x1ix3i-x2ix3ix1i2+x2i2,re=∑i=1nΔue1iΔue2iΔue3iΔue3ix2i-Δue2ix3iΔue1ix3i-Δue3ix1iΔue2ix1i-Δue1ix2i.


## 4.4.1. Assumed Intraelement and Auxiliary Frame Fields

The intraelement displacement fields are approximated in terms of a linear combination of fundamental solutions of the problem as(95)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as(96)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysnsce=c11c21c31⋯c1nc2nc3nTin which x and ysj are, respectively, the field point and source point in the local coordinate system (x1,x2). The fundamental solution uij*(x,ysj) is given by [40](97)uij*x,ysj=116π1-νGr3-4νδij+r,ir,j,where ri=xi-xis,r=r12+r22+r32, ns is the number of source points. The source point ysj(j=1,2,…,ns) can also be generated by means of the following method [36] as in two-dimensional cases (98)ys=x0+γx0-xc,where γ is a dimensionless coefficient, x0 is the point on the element boundary (the nodal point in this work), and xc is the geometrical centroid of the element (see Figure 5).Figure 5
Intraelement field and frame field of a hexahedron HFS-FEM element for 3D elastic problem (the source points and centroid of 20-node element are omitted in the figure for clarity and clear view, which is similar to that of the 8-node element).
(a)
(b)According to (50) and (49), the corresponding stress fields can be expressed as (99)σx=σ11σ22σ33σ23σ31σ12T=Tece,where (100)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ133*x,y1σ233*x,y1σ333*x,y1⋯σ133*x,ynsσ233*x,ynsσ333*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ212*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by(101)σijk*x,y=-18π1-νr2·1-2νr,kδij+r,jδki-r,iδjk+3r,ir,jr,k.As a consequence, the traction can be written in the form (102)t1t2t3=nσ=Qecein which (103)Qe=nTe,n=n1000n3n20n20n30n100n3n2n10.To link the unknownce and the nodal displacement de, the frame is defined as (104)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the surface 2-3-7-6 of a particular 8-node brick element (see Figure 5) as an example, matrix N~e and vector de can be expressed as (105)N~e=0N-1N-200N-4N-30,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (106)N-i=N~i000N~i000N~i,0=000000000,where N~i (i=1–4) can be expressed by natural coordinate ξ,η∈[-1,1](107)N~1=1+ξ1+η4,N~2=1-ξ1+η4,N~3=1-ξ1-η4,N~4=1+ξ1-η4and ξi,ηi is the natural coordinate of the i-node of the element (Figure 6).Figure 6
Typical linear interpolation for the frame fields of 3D brick elements.


## 4.4.2. Modified Functional for Hybrid Finite Element Method

In the absence of body forces, the hybrid functionalΠme used for deriving the present HFS-FEM is written as [22](108)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.By applying the Gaussian theorem, (108) can be simplified as(109)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for HFS finite element model:(110)Πme=12∫ΓetiuidΓ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.The functional (110) contains only boundary integrals of the element and will be used to derive HFS-FEM formulation for the three-dimensional elastic problem in the following section.


## 4.4.3. Element Stiffness Matrix

Substituting (95), (102), and (104) into the functional (110), we have (111)Πme=-12ceTHece+ceTGede-deTge,where (112)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.To enforce interelement continuity on the common element boundary, the unknown vector ce should be expressed in terms of nodal DOF de. The stationary condition of the functional Πme with respect to ce and de yields again, respectively, (33) and (34).


## 4.4.4. Numerical Integral over Element

Considering a surface of the 3D hexahedron element, as shown in Figure6, the vector normal to the surface can be obtained by (113)vn=vξ×vη=vnxvnyvnz=dxdξdydξdzdξ×dxdηdydηdzdη=dydξdzdη-dydηdzdξdzdξdxdη-dzdηdxdξdxdξdydη-dxdηdydξ,where vξ and vη are the tangential vectors in the ξ-direction and η-direction, respectively, calculated by(114)vξ=dxdξdydξdzdξ=∑i=1nd∂Niξ,η∂ξxiyizi,vη=dxdηdydηdzdη=∑i=1nd∂Niξ,η∂ηxiyizi,where nd is the number of nodes of the surface and xi,yi,zi are the nodal coordinates. Thus the unit normal vector is given by(115)n=vnvn,where(116)Jξ,η=vn=vnx2+vny2+vnz2is the Jacobian of the transformation from Cartesian coordinates (x,y) to natural coordinates (ξ,η).For theH matrix, we introduce the matrix function (117)Fx,y=Fijx,ym×m=QeTNe.Then we can get(118)He=∫ΓeQeTNedΓ=∫ΓeFx,ydΓand we rewrite it to the component form as (119)Hij=∫ΓeFijx,ydS=∑l=1nf∫ΓelFijx,ydS.Using the relationship (120)dS=Jξ,ηdξdηand the Gaussian numerical integration, we can obtain(121)Hij=∑l=1nf∫-11Fijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtFijxξs,ηt,yξs,ηtJξs,ηt,where nf and np are, respectively, the number of surface of the 3D element and the number of Gaussian integral points in each direction of the element surface. Similarly, we can calculate the Ge matrix by(122)Gij=∑l=1nf∫-11F~ijxξ,η,yξ,ηJξ,ηdξdη≈∑l=1nf∑s=1np∑t=1npwswtF~ijxξs,ηt,yξs,ηtJξs,ηt.It should be mentioned that the calculation of vector ge in equation is the same as that in the conventional FEM, so it is convenient to incorporate the proposed HFS-FEM into the standard FEM program. Besides, the stress and traction estimations are directly computed from (99) and (100), respectively. The boundary displacements can be directly computed from (104) while the displacements at interior points of element can be determined from (95) plus the recovered rigid-body modes in each element, which is introduced in the following section.


## 4.4.5. Recovery of Rigid-Body Motion Terms

As in Section2, the least square method is employed to recover the missing terms of rigid-body motions. The missing terms can be recovered by setting for the augmented internal field(123)ue=Nece+1000x3-x2010-x30x1001x2-x10c0and using a least-square procedure to match ueh and u~eh at the nodes of the element boundary(124)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+u3i-u~3i2.The above equation finally yields(125)c0=Re-1re,where(126)Re=∑i=1n1000x3i-x2i010-x3i0x1i001x2i-x1i00-x3ix2ix2i2+x3i2-x1ix2i-x1ix3ix3i0-x1i-x1ix2ix1i2+x3i2-x2ix3i-x2ix1i0-x1ix3i-x2ix3ix1i2+x2i2,re=∑i=1nΔue1iΔue2iΔue3iΔue3ix2i-Δue2ix3iΔue1ix3i-Δue3ix1iΔue2ix1i-Δue1ix2i.


## 5. Thermoelasticity Problems

Thermoelasticity problems arise in many practical designs such as steam and gas turbines, jet engines, rocket motors, and nuclear reactors. Thermal stress induced in these structures is one of the major concerns in product design and analysis. The general thermoelasticity is governed by two time-dependent coupled differential equations: the heat conduction equation and the Navier equation with thermal body force [44]. In many engineering applications, the coupling term of the heat equation and the inertia term in Navier equation are generally negligible [44]. As a consequence, most of the analyses are employing the uncoupled thermoelasticity theory which is adopted in this topic [45–52]. In this section, the HFS-FEM is presented to solve 2D and 3D thermoelastic problems with considering arbitrary body forces and temperature changes [53]. The method used herein is similar to that in Section 4.


### 5.1. Basic Equations for Thermoelasticity

Consider an isotropic material in a finite domainΩ and let (x1,x2,x3) denote the coordinates in Cartesian coordinate system. The equilibrium governing equations of the thermoelasticity with the body force are expressed as (127)σij,j=-bi,where σij is the stress tensor, bi is the body force vector, and i,j=1,2,3. The generalized thermoelastic stress-strain relations and kinematical relation are given as (128)σij,j=2Gν1-2νδije+2Geij-mδijT,eij=12ui,j+uj,i,where eij is the strain tensor, ui is the displacement vector, T is the temperature change, G is the shear modulus, ν is the Poisson’s ratio, δij is the Kronecker delta, and (129)m=2Gα1+v1-2vis the thermal constant with α being the coefficient of linear thermal expansion. Substituting (128) into (127), the equilibrium equations may be rewritten as (130)Gui,jj+G1-2νuj,ji=mT,i-bi.For a well-posed boundary value problem, the following boundary conditions, either displacement or traction boundary condition, should be prescribed as(131)ui=u-ion⁡Γu,ti=t-ion⁡Γt,where Γu∪Γt=Γ is the boundary of the solution domain Ω, u-i and t-i are the prescribed boundary values, and (132)ti=σijnjis the boundary traction, in which nj denotes the boundary outward normal.


### 5.2. The Method of Particular Solution

For the governing equation (130) given in the previous section, the inhomogeneous term mT,i-bi can be eliminated by employing the method of particular solution [16, 36, 44]. Using superposition principle, we decompose the displacement ui into two parts, the homogeneous solution uih and the particular solution uip as follows:(133)ui=uih+uipin which the particular solution uip should satisfy the governing equation (134)Gui,jjp+G1-2νuj,jip=mT,i-bibut does not necessarily satisfy any boundary condition. It should be pointed out that its solution is not unique and can be obtained by various numerical techniques. However, the homogeneous solution should satisfy(135)Gui,jjh+G1-2νuj,jih=0with modified boundary conditions(136)uih=u-i-uip,on⁡Γu,(137)tih=t-i+mTni-tip,on⁡Γt.From above equations, it can be seen that once the particular solution is known, the homogeneous solution uih in (135)–(137) can be solved by (135). In the following section, RBF approximation is described to illustrate the particular solution procedure, and the HFS-FEM is presented which can be used for solving (135)–(137).


### 5.3. Radial Basis Function Approximation

RBF is to be used to approximate the body forcebi and the temperature field T in order to obtain the particular solution. To implement this approximation, we may consider two different ways: one is to treat body force bi and the temperature field T separately as in Tsai [54]. The other is to treat mT,i-bi as a whole [53]. Here we demonstrated that the performance of the latter one is usually better than the former one.


#### 5.3.1. Interpolating Temperature and Body Force Separately

The body forcebi and temperature T are assumed to be by the following two equations: (138)bi≈∑j=1Nαijφji=1,2inR2,i=1,2,3inR3,T≈∑j=1Nβjφj,where N is the number of interpolation points, φj are the basis functions, and αij and βj are the coefficients to be determined by collocation. Subsequently, the approximate particular solution can be written as follows: (139)uip=∑k=13∑j=1NαijΦikj+∑j=1NβjΨij,where Φikj and Ψij are the approximated particular solution kernels. Once the RBF is selected, the problem of finding a particular solution will be reduced to solve the following equations:(140)GΦil,kk+G1-2νΦkl,ki=-δilφ,(141)GΨi,kk+G1-2νΨk,ki=mφ,i.To solve (140), the displacement is expressed in terms of the Galerkin-Papkovich vectors [43, 55–57](142)Φik=1-νGFik,mm-12GFmk,mi.Substituting (142) into (140), one can obtain the following biharmonic equation:(143)∇4Fil=-11-νδilφ.If taking the Spline Type RBF φ=r2n-1, one can get the following solutions: (144)Fli=-δli1-vr2n+32n+122n+32R2forn=1,2,3,…,Φli=A0A1δli+A2r,ir,lR2forn=1,2,3,…,where(145)A0=-12G1-vr2n+12n+122n+3,A1=5+4n-2v2n+3,A2=-2n+1for two-dimensional problem and (146)Fli=-δli1-vr2n+32n+12n+22n+32n+4R3forn=1,2,3,…,(147)Φli=A0A1δli+A2r,ir,lR3forn=1,2,3,…,where(148)A0=-18G1-vr2n+1n+1n+22n+1,A1=7+4n-4vn+2,A2=-2n+1for three-dimensional problem, where rj represents the Euclidean distance of the given point (x1,x2) from a fixed point (x1j,x2j) in the domain of interest. The corresponding stress particular solution can be obtained by(149)Slij=GΦli,j+Φlj,i+λδijΦlk,k,where λ=(2v/1-2v)G. Substituting (147) into (149), one can obtain(150)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR2forn=1,2,3,…,where(151)B0=-11-vr2n2n+12n+3,B1=2n+2-v2n+3,B2=v2n+3-1,B3=1-2nfor two-dimensional problem, and substituting (147) into (149), one can obtain(152)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR3forn=1,2,3,…,where(153)B0=-141-vr2nn+1n+2,B1=3+2n-2vn+2,B2=2v(n+2)-1,B3=1-2nfor three-dimensional problem.To solve (141), one can treat Ψi as the gradient of a scalar function (154)Ψi=U,i.Substituting (154) into (141) obtains the Poisson’s equation(155)∇2U=m1-2ν2G1-νφ.Thus, taking φ=r2n-1, its particular solution can be obtained [55] as follows:(156)U=m1-2ν2G1-νr2n+12n+12R2forn=1,2,3,…,U=m1-2ν2G1-νr2n+12n+12n+2R3forn=1,2,3,….Then from (154) we can get Ψi as follows:(157)Ψi=m1-2ν2G1-νr,ir2n2n+1R2forn=1,2,3,…,Ψi=m1-2ν2G1-νr,ir2n2n+2R3forn=1,2,3,….The corresponding stress particular solution can be obtained by substituting (147) into(158)Sij=GΨi,j+Ψj,i+λδijΨk,k.Then we have(159)Sij=mr2n-11-v2n+1(1+2nv)δij+1-2v2n-1r,ir,jR2forn=1,2,3,…,Sij=mr2n-11-v2n+21-2ν2n-1r,ir,j+δij1+2nvR3forn=1,2,3,….


#### 5.3.2. Interpolating Temperature and Body Force Together

Considering the temperature gradient plays the role of body force, we can approximatemT,i-bi together by the following equation: (160)mT,i-bi≈∑j=1Nαijφj.Thus the approximate particular solution can be written as (161)uip=∑k=13∑j=1NαijΦikj.Consequently, we can follow the same way as that for body force in Method 1 and employing (140) and (142)–(150) to obtain the desired Φikj and Slij, which are the same as those for body force case only. It is noted that Method 2 has a relatively smaller number of equations to solve for the coefficients and the condition number of the coefficient matrix is smaller as well, which will be beneficial to the solution.Once we have obtained the particular solutions of (133), we can use them to get the modified boundary conditions in (136) to obtain the homogeneous solution by considering (135). Then, we can employ the HFS-FEM described in Section 3 for 2D problem and Section 4.4 for 3D problem to obtain the homogeneous solutions.


## 5.1. Basic Equations for Thermoelasticity

Consider an isotropic material in a finite domainΩ and let (x1,x2,x3) denote the coordinates in Cartesian coordinate system. The equilibrium governing equations of the thermoelasticity with the body force are expressed as (127)σij,j=-bi,where σij is the stress tensor, bi is the body force vector, and i,j=1,2,3. The generalized thermoelastic stress-strain relations and kinematical relation are given as (128)σij,j=2Gν1-2νδije+2Geij-mδijT,eij=12ui,j+uj,i,where eij is the strain tensor, ui is the displacement vector, T is the temperature change, G is the shear modulus, ν is the Poisson’s ratio, δij is the Kronecker delta, and (129)m=2Gα1+v1-2vis the thermal constant with α being the coefficient of linear thermal expansion. Substituting (128) into (127), the equilibrium equations may be rewritten as (130)Gui,jj+G1-2νuj,ji=mT,i-bi.For a well-posed boundary value problem, the following boundary conditions, either displacement or traction boundary condition, should be prescribed as(131)ui=u-ion⁡Γu,ti=t-ion⁡Γt,where Γu∪Γt=Γ is the boundary of the solution domain Ω, u-i and t-i are the prescribed boundary values, and (132)ti=σijnjis the boundary traction, in which nj denotes the boundary outward normal.


## 5.2. The Method of Particular Solution

For the governing equation (130) given in the previous section, the inhomogeneous term mT,i-bi can be eliminated by employing the method of particular solution [16, 36, 44]. Using superposition principle, we decompose the displacement ui into two parts, the homogeneous solution uih and the particular solution uip as follows:(133)ui=uih+uipin which the particular solution uip should satisfy the governing equation (134)Gui,jjp+G1-2νuj,jip=mT,i-bibut does not necessarily satisfy any boundary condition. It should be pointed out that its solution is not unique and can be obtained by various numerical techniques. However, the homogeneous solution should satisfy(135)Gui,jjh+G1-2νuj,jih=0with modified boundary conditions(136)uih=u-i-uip,on⁡Γu,(137)tih=t-i+mTni-tip,on⁡Γt.From above equations, it can be seen that once the particular solution is known, the homogeneous solution uih in (135)–(137) can be solved by (135). In the following section, RBF approximation is described to illustrate the particular solution procedure, and the HFS-FEM is presented which can be used for solving (135)–(137).


## 5.3. Radial Basis Function Approximation

RBF is to be used to approximate the body forcebi and the temperature field T in order to obtain the particular solution. To implement this approximation, we may consider two different ways: one is to treat body force bi and the temperature field T separately as in Tsai [54]. The other is to treat mT,i-bi as a whole [53]. Here we demonstrated that the performance of the latter one is usually better than the former one.


### 5.3.1. Interpolating Temperature and Body Force Separately

The body forcebi and temperature T are assumed to be by the following two equations: (138)bi≈∑j=1Nαijφji=1,2inR2,i=1,2,3inR3,T≈∑j=1Nβjφj,where N is the number of interpolation points, φj are the basis functions, and αij and βj are the coefficients to be determined by collocation. Subsequently, the approximate particular solution can be written as follows: (139)uip=∑k=13∑j=1NαijΦikj+∑j=1NβjΨij,where Φikj and Ψij are the approximated particular solution kernels. Once the RBF is selected, the problem of finding a particular solution will be reduced to solve the following equations:(140)GΦil,kk+G1-2νΦkl,ki=-δilφ,(141)GΨi,kk+G1-2νΨk,ki=mφ,i.To solve (140), the displacement is expressed in terms of the Galerkin-Papkovich vectors [43, 55–57](142)Φik=1-νGFik,mm-12GFmk,mi.Substituting (142) into (140), one can obtain the following biharmonic equation:(143)∇4Fil=-11-νδilφ.If taking the Spline Type RBF φ=r2n-1, one can get the following solutions: (144)Fli=-δli1-vr2n+32n+122n+32R2forn=1,2,3,…,Φli=A0A1δli+A2r,ir,lR2forn=1,2,3,…,where(145)A0=-12G1-vr2n+12n+122n+3,A1=5+4n-2v2n+3,A2=-2n+1for two-dimensional problem and (146)Fli=-δli1-vr2n+32n+12n+22n+32n+4R3forn=1,2,3,…,(147)Φli=A0A1δli+A2r,ir,lR3forn=1,2,3,…,where(148)A0=-18G1-vr2n+1n+1n+22n+1,A1=7+4n-4vn+2,A2=-2n+1for three-dimensional problem, where rj represents the Euclidean distance of the given point (x1,x2) from a fixed point (x1j,x2j) in the domain of interest. The corresponding stress particular solution can be obtained by(149)Slij=GΦli,j+Φlj,i+λδijΦlk,k,where λ=(2v/1-2v)G. Substituting (147) into (149), one can obtain(150)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR2forn=1,2,3,…,where(151)B0=-11-vr2n2n+12n+3,B1=2n+2-v2n+3,B2=v2n+3-1,B3=1-2nfor two-dimensional problem, and substituting (147) into (149), one can obtain(152)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR3forn=1,2,3,…,where(153)B0=-141-vr2nn+1n+2,B1=3+2n-2vn+2,B2=2v(n+2)-1,B3=1-2nfor three-dimensional problem.To solve (141), one can treat Ψi as the gradient of a scalar function (154)Ψi=U,i.Substituting (154) into (141) obtains the Poisson’s equation(155)∇2U=m1-2ν2G1-νφ.Thus, taking φ=r2n-1, its particular solution can be obtained [55] as follows:(156)U=m1-2ν2G1-νr2n+12n+12R2forn=1,2,3,…,U=m1-2ν2G1-νr2n+12n+12n+2R3forn=1,2,3,….Then from (154) we can get Ψi as follows:(157)Ψi=m1-2ν2G1-νr,ir2n2n+1R2forn=1,2,3,…,Ψi=m1-2ν2G1-νr,ir2n2n+2R3forn=1,2,3,….The corresponding stress particular solution can be obtained by substituting (147) into(158)Sij=GΨi,j+Ψj,i+λδijΨk,k.Then we have(159)Sij=mr2n-11-v2n+1(1+2nv)δij+1-2v2n-1r,ir,jR2forn=1,2,3,…,Sij=mr2n-11-v2n+21-2ν2n-1r,ir,j+δij1+2nvR3forn=1,2,3,….


### 5.3.2. Interpolating Temperature and Body Force Together

Considering the temperature gradient plays the role of body force, we can approximatemT,i-bi together by the following equation: (160)mT,i-bi≈∑j=1Nαijφj.Thus the approximate particular solution can be written as (161)uip=∑k=13∑j=1NαijΦikj.Consequently, we can follow the same way as that for body force in Method 1 and employing (140) and (142)–(150) to obtain the desired Φikj and Slij, which are the same as those for body force case only. It is noted that Method 2 has a relatively smaller number of equations to solve for the coefficients and the condition number of the coefficient matrix is smaller as well, which will be beneficial to the solution.Once we have obtained the particular solutions of (133), we can use them to get the modified boundary conditions in (136) to obtain the homogeneous solution by considering (135). Then, we can employ the HFS-FEM described in Section 3 for 2D problem and Section 4.4 for 3D problem to obtain the homogeneous solutions.


## 5.3.1. Interpolating Temperature and Body Force Separately

The body forcebi and temperature T are assumed to be by the following two equations: (138)bi≈∑j=1Nαijφji=1,2inR2,i=1,2,3inR3,T≈∑j=1Nβjφj,where N is the number of interpolation points, φj are the basis functions, and αij and βj are the coefficients to be determined by collocation. Subsequently, the approximate particular solution can be written as follows: (139)uip=∑k=13∑j=1NαijΦikj+∑j=1NβjΨij,where Φikj and Ψij are the approximated particular solution kernels. Once the RBF is selected, the problem of finding a particular solution will be reduced to solve the following equations:(140)GΦil,kk+G1-2νΦkl,ki=-δilφ,(141)GΨi,kk+G1-2νΨk,ki=mφ,i.To solve (140), the displacement is expressed in terms of the Galerkin-Papkovich vectors [43, 55–57](142)Φik=1-νGFik,mm-12GFmk,mi.Substituting (142) into (140), one can obtain the following biharmonic equation:(143)∇4Fil=-11-νδilφ.If taking the Spline Type RBF φ=r2n-1, one can get the following solutions: (144)Fli=-δli1-vr2n+32n+122n+32R2forn=1,2,3,…,Φli=A0A1δli+A2r,ir,lR2forn=1,2,3,…,where(145)A0=-12G1-vr2n+12n+122n+3,A1=5+4n-2v2n+3,A2=-2n+1for two-dimensional problem and (146)Fli=-δli1-vr2n+32n+12n+22n+32n+4R3forn=1,2,3,…,(147)Φli=A0A1δli+A2r,ir,lR3forn=1,2,3,…,where(148)A0=-18G1-vr2n+1n+1n+22n+1,A1=7+4n-4vn+2,A2=-2n+1for three-dimensional problem, where rj represents the Euclidean distance of the given point (x1,x2) from a fixed point (x1j,x2j) in the domain of interest. The corresponding stress particular solution can be obtained by(149)Slij=GΦli,j+Φlj,i+λδijΦlk,k,where λ=(2v/1-2v)G. Substituting (147) into (149), one can obtain(150)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR2forn=1,2,3,…,where(151)B0=-11-vr2n2n+12n+3,B1=2n+2-v2n+3,B2=v2n+3-1,B3=1-2nfor two-dimensional problem, and substituting (147) into (149), one can obtain(152)Slij=B0B1r,jδli+r,iδjl+B2δijr,l+B3r,ir,jr,lR3forn=1,2,3,…,where(153)B0=-141-vr2nn+1n+2,B1=3+2n-2vn+2,B2=2v(n+2)-1,B3=1-2nfor three-dimensional problem.To solve (141), one can treat Ψi as the gradient of a scalar function (154)Ψi=U,i.Substituting (154) into (141) obtains the Poisson’s equation(155)∇2U=m1-2ν2G1-νφ.Thus, taking φ=r2n-1, its particular solution can be obtained [55] as follows:(156)U=m1-2ν2G1-νr2n+12n+12R2forn=1,2,3,…,U=m1-2ν2G1-νr2n+12n+12n+2R3forn=1,2,3,….Then from (154) we can get Ψi as follows:(157)Ψi=m1-2ν2G1-νr,ir2n2n+1R2forn=1,2,3,…,Ψi=m1-2ν2G1-νr,ir2n2n+2R3forn=1,2,3,….The corresponding stress particular solution can be obtained by substituting (147) into(158)Sij=GΨi,j+Ψj,i+λδijΨk,k.Then we have(159)Sij=mr2n-11-v2n+1(1+2nv)δij+1-2v2n-1r,ir,jR2forn=1,2,3,…,Sij=mr2n-11-v2n+21-2ν2n-1r,ir,j+δij1+2nvR3forn=1,2,3,….


## 5.3.2. Interpolating Temperature and Body Force Together

Considering the temperature gradient plays the role of body force, we can approximatemT,i-bi together by the following equation: (160)mT,i-bi≈∑j=1Nαijφj.Thus the approximate particular solution can be written as (161)uip=∑k=13∑j=1NαijΦikj.Consequently, we can follow the same way as that for body force in Method 1 and employing (140) and (142)–(150) to obtain the desired Φikj and Slij, which are the same as those for body force case only. It is noted that Method 2 has a relatively smaller number of equations to solve for the coefficients and the condition number of the coefficient matrix is smaller as well, which will be beneficial to the solution.Once we have obtained the particular solutions of (133), we can use them to get the modified boundary conditions in (136) to obtain the homogeneous solution by considering (135). Then, we can employ the HFS-FEM described in Section 3 for 2D problem and Section 4.4 for 3D problem to obtain the homogeneous solutions.


## 6. Anisotropic Composite Materials

In materials science, composite laminates are usually assemblies of layers of fibrous composite materials which can be joined together to provide required engineering properties, such as specified in-plane stiffness, bending stiffness, strength, and coefficient of thermal expansion [58]. Individual layers (or laminas) of the laminates consist of high-modulus, high-strength fibers in a polymeric, metallic, or ceramic matrix material. From the viewpoint of micromechanics, the fiber and matrix in each lamina can be treated as the inclusion and matrix, respectively. On the other hand, from the viewpoint of macromechanics, both a lamina and the whole laminate can be viewed as a general anisotropic body by classical lamination theory. Hence, the analysis of anisotropic bodies is important for understanding of the micro- or macromechanical behavior of composites [58].In the literature, there are two main approaches dealing with generalized two-dimensional anisotropic elastic problems. One is Lekhnitskii formalism [59, 60] which begins with stresses as basic variables and the other is Stroh formalism [61, 62] which starts with displacements as basic variables. Both of them are formulated in terms of complex variable functions. The Stroh formalism, which has been shown to be elegant and powerful, is used to find the analytical solutions for the corresponding infinite bodies [61, 63]. The formalism is also widely employed in the derivation of the inclusion or crack problems of anisotropic materials [64, 65]. Because of the limitations of the analytical solutions which are only available for some problems with simple geometry and boundary conditions [66, 67], numerical methods such as finite element method (FEM), boundary element method (BEM), mesh free method (MFM), and Hybrid-Trefftz (HT) FEM are usually resorted to solve more complex problems with complicated boundary constraints and loading conditions [68–70]. In this section, we presented the HFS-FEM for analyzing anisotropic composite materials based on the associated fundamental solutions in terms of Stroh formalism [71].


### 6.1. Linear Anisotropic Elasticity

#### 6.1.1. Basic Equations and Stroh Formalism

In the Cartesian coordinate system (x1,x2,x3), if we neglect the body force bi, the equilibrium equations, stress-strain laws, and strain-displacement equations for anisotropic elasticity are [61] (162)σij,j=0,(163)σij=Cijklekl,(164)eij=12ui,j+uj,i,where i,j=1,2,3, Cijkl is the fourth-rank anisotropic elasticity tensor. The equilibrium equations can be rewritten in terms of displacements by substituting (163) and (164) into (162) as (165)Cijkluk,jl=0.The boundary conditions of the boundary value problem (163)–(165) are (166)ui=u-ion⁡Γu,ti=σijnj=t-ion⁡Γt,where u-i and t-i are the prescribed boundary displacement and traction vector, respectively. In addition, ni is the unit outward normal to the boundary and Γ=Γu+Γt is the boundary of the solution domain Ω.For the generalized two-dimensional deformation of anisotropic elasticityui is assumed to depend on x1 and x2 only. Based on this assumption, the general solution to (165) can be written as [61, 62](167)u=2ReAfz,φ=2ReBfz,where u=(u1,u2,u3)T is the displacement vector, φ=(φ1,φ2,φ3)T is the stress function vector, and f(z)=[f1(z1),f2(z2),f3(z3)]T is a function vector composed of three holomorphic complex functions fα(zα),α=1,2,3, which is an arbitrary function with argument zα=x1+pαx2 and will be determined by satisfying the boundary and loading conditions of a given problem. In (167), Re stands for the real part of a complex number, pα are the material eigenvalues with positive imaginary part, and A=[a1,a2,a3] and B=[b1,b2,b3] are 3 × 3 complex matrices formed by the material eigenvector associated with pα, which can be obtained by the following Eigen relations [61]:(168)Nξ=pξ,where N is a 6 × 6 foundational elasticity matrix and ξ is a 6 × 1 column vector defined by(169)N=N1N2N3N1T,ξ=ab,where N1=-T-1RT,N2=T-1,N3=RT-1RT-Q, and the matrices Q,R and T are 3 × 3 matrices extracted from Cijkl as follows:(170)Qik=Ci1k1,Rik=Ci1k2,Tik=Ci2k2.The stresses can be obtained from the derivation of stress functions φ as follows:(171)σi1=2ReLf′z,σi2=2ReBf′z,where(172)L=-p1b1,-p2b2,-p3b3,-p4b4.


#### 6.1.2. Foundational Solutions

To find the fundamental solution needed in our analysis, we have to first derive the Green’s function of the problem: an infinite homogeneous anisotropic elastic medium loaded by a concentrated point force (or line force for two-dimensional problems)p^=(p^1,p^2,p^3) applied at an internal point x^=(x^1,x^2) far from the boundary. The boundary conditions of this problem can be written as (173)∫Cdϕ=p^foranyclosedcurveCenclosingx^,∫Cdu=p^foranyclosedcurveC,limx→∞σij=0.Thus, the Green’s function satisfying the above boundary conditions is found to be [72] (174)fz=12πiln⁡zα-z^αATp^.Therefore, fundamental solutions of the problem can be expressed as(175)u=1πImAln⁡⁡zα-z^αATp^,ϕ=1πImBln⁡⁡zα-z^αATp^.The corresponding stress components can be obtained from stress function ϕ as(176)σi1*=-ϕ,2=-1πImBpαzα-z^αATp^,σi2*=ϕ,1=1πImB1zα-z^αATp^,where p^ are chosen to be (1,0,0)T,(0,1,0)T, and (0,0,1)T, respectively, · stands for the diagonal matrix corresponding to subscript α, Im denotes the imagery part of a complex number, and superscript T denotes the matrix transpose.


#### 6.1.3. Coordinate Transformation

A typical composite laminate consists of individual layers, which are usually made of unidirectional plies with the same or regularly alternating orientation. A layer is generally referred to the global coordinate framex1,x2, and x3 of the structural element rather than to coordinates x,y, and z associated with the ply orientation. So it is necessary to transform the constitutive relationship of each layer from the material coordinate frame to the uniform global coordinate frame.For the two coordinate systems mentioned in Figure7, the angle between the axis-x1 and axis-x is denoted by φ, which is positive along the anticlockwise direction, and the relationship for transformation of stress and strain components between the local material coordinates and global coordinates is given by(177)σ11,σ22,σ23,σ31,σ12T=T-1σxx,σyy,σyz,σzx,σxyT,εxx,εyy,εyz,εzx,εxyT=T-1Tε11,ε22,ε23,ε31,ε12T,where the transformation matrix T and its inverse matrix are defined as(178)T=c2s2002css2c200-2cs00c-s000sc0-cscs00c2-s2,T-1=c2s200-2css2c2002cs00cs000-sc0cs-cs00c2-s2with c=cos⁡(φ),s=sin(φ). Therefore, the constitutive relationship in the global coordinate system is given by (179)σ11,σ22,σ23,σ31,σ12T=T-1CT-1Tε11,ε22,ε23,ε31,ε12T.Figure 7
Schematic of the relationship between global coordinate system (x1,x2) and local material coordinate system (x,y).


### 6.2. Formulations of HFS-FEM

#### 6.2.1. Assumed Independent Fields

The intraelement displacement fields for a particular elemente is approximated in terms of a linear combination of fundamental solutions of the problem as (180)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as (181)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nc2nc3nTin which ns is the number of source points, x and ysj are, respectively, the field point and source point in the coordinate system (x1,x2) local to the element under consideration. The fundamental solution uij*(x,ysj) is given by (175) for general elements. γ is a dimensionless coefficient for determining source points as described in previous sections.The corresponding stress fields can be expressed as(182)σx=σ11σ22σ23σ31σ12T=Tece,where (183)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ231*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by (176) when p^i is selected to be (1,0,0)T, (0,1,0)T, and (0,0,1)T, respectively. As a consequence, the traction can be written as (184)t1t2t3=nσ=Qecein which (185)Qe=nTe,n=n100n3n20n2n30n100n2n10.The unknown ce in (180) and (182) may be calculated using a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as in conventional FEM (see Figure 1). Thus, the frame field is defined as (186)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N~e and de can be expressed as (187)N~e=00N-1N-2N-3000,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (188)N-i=N~i000N~i000N~i,0=000000000and N~1,N~2, and N~3 are expressed by natural coordinate ξ∈[-1,1](189)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


#### 6.2.2. Modified Functional for HFS-FEM

With the assumption of two distinct intraelement field and frame field for elements, we can establish the modified variational principle based on (165) and (166) for the hybrid finite element method of anisotropic materials [21, 73]. In the absence of the body forces, the hybrid variational functional Πme for a particular element e is constructed as(190)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ,where the boundary Γe of the element e is (191)Γe=Γeu∪Γet∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γtand ΓeI is the interelement boundary of element e. Performing a variation of Πm, one obtains(192)δΠme=∬Ωeσijδui,jdΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~i-δuidΓ.Applying Gaussian theorem (193)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force (194)ti=σijnjwe obtain (195)δΠe=-∬Ωeσij,jδuidΩ+∫ΓetiδuidΓ-∫Γett-iδu~idΓ.Then, substituting (195) into (192) gives(196)δΠme=-∬Ωeσij,jδuidΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~idΓ.Considering the fact that(197)∫Γeutiδu~idΓ=0we can finally obtain the following form:(198)δΠme=-∬Ωeσij,jδuidΩ+∫Γetti-t-iδu~idΓ+∫Γeu~i-uiδtidΓ+∫ΓItiδu~idΓ.Therefore, the Euler equations for (190) result in (165) and (166) because the quantities δui,δti, and δu~i may be arbitrary. As for the continuity condition between elements, it can be easily seen from the following variational of two adjacent elements such as e and f(199)δΠm(e∪f)=-∬Ωe∪Ωfσij,jδuidΩ+∫Γet+Γetti-t-iδu~idΓ+∫Γe+Γfu~i-uiδtidΓ+∫ΓefItie+tifδu~i.This indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (190) can be easily proved by the so-called “second variational approach” as well, which indicates functional (190) has a local extreme. Therefore, we can conclude that the variational functional (190) can be used for deriving hybrid finite element formulations.


#### 6.2.3. Element Stiffness Equation

Using Gaussian theorem and equilibrium equations, all domain integrals in (190) can be converted into boundary integrals as follows:(200)Πme=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.Substituting (180), (184), and (186) into the functional (200) yields the formulation as (201)Πme=-12ceTHece+ceTGede-deTge,where(202)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and  (34).


## 6.1. Linear Anisotropic Elasticity

### 6.1.1. Basic Equations and Stroh Formalism

In the Cartesian coordinate system (x1,x2,x3), if we neglect the body force bi, the equilibrium equations, stress-strain laws, and strain-displacement equations for anisotropic elasticity are [61] (162)σij,j=0,(163)σij=Cijklekl,(164)eij=12ui,j+uj,i,where i,j=1,2,3, Cijkl is the fourth-rank anisotropic elasticity tensor. The equilibrium equations can be rewritten in terms of displacements by substituting (163) and (164) into (162) as (165)Cijkluk,jl=0.The boundary conditions of the boundary value problem (163)–(165) are (166)ui=u-ion⁡Γu,ti=σijnj=t-ion⁡Γt,where u-i and t-i are the prescribed boundary displacement and traction vector, respectively. In addition, ni is the unit outward normal to the boundary and Γ=Γu+Γt is the boundary of the solution domain Ω.For the generalized two-dimensional deformation of anisotropic elasticityui is assumed to depend on x1 and x2 only. Based on this assumption, the general solution to (165) can be written as [61, 62](167)u=2ReAfz,φ=2ReBfz,where u=(u1,u2,u3)T is the displacement vector, φ=(φ1,φ2,φ3)T is the stress function vector, and f(z)=[f1(z1),f2(z2),f3(z3)]T is a function vector composed of three holomorphic complex functions fα(zα),α=1,2,3, which is an arbitrary function with argument zα=x1+pαx2 and will be determined by satisfying the boundary and loading conditions of a given problem. In (167), Re stands for the real part of a complex number, pα are the material eigenvalues with positive imaginary part, and A=[a1,a2,a3] and B=[b1,b2,b3] are 3 × 3 complex matrices formed by the material eigenvector associated with pα, which can be obtained by the following Eigen relations [61]:(168)Nξ=pξ,where N is a 6 × 6 foundational elasticity matrix and ξ is a 6 × 1 column vector defined by(169)N=N1N2N3N1T,ξ=ab,where N1=-T-1RT,N2=T-1,N3=RT-1RT-Q, and the matrices Q,R and T are 3 × 3 matrices extracted from Cijkl as follows:(170)Qik=Ci1k1,Rik=Ci1k2,Tik=Ci2k2.The stresses can be obtained from the derivation of stress functions φ as follows:(171)σi1=2ReLf′z,σi2=2ReBf′z,where(172)L=-p1b1,-p2b2,-p3b3,-p4b4.


### 6.1.2. Foundational Solutions

To find the fundamental solution needed in our analysis, we have to first derive the Green’s function of the problem: an infinite homogeneous anisotropic elastic medium loaded by a concentrated point force (or line force for two-dimensional problems)p^=(p^1,p^2,p^3) applied at an internal point x^=(x^1,x^2) far from the boundary. The boundary conditions of this problem can be written as (173)∫Cdϕ=p^foranyclosedcurveCenclosingx^,∫Cdu=p^foranyclosedcurveC,limx→∞σij=0.Thus, the Green’s function satisfying the above boundary conditions is found to be [72] (174)fz=12πiln⁡zα-z^αATp^.Therefore, fundamental solutions of the problem can be expressed as(175)u=1πImAln⁡⁡zα-z^αATp^,ϕ=1πImBln⁡⁡zα-z^αATp^.The corresponding stress components can be obtained from stress function ϕ as(176)σi1*=-ϕ,2=-1πImBpαzα-z^αATp^,σi2*=ϕ,1=1πImB1zα-z^αATp^,where p^ are chosen to be (1,0,0)T,(0,1,0)T, and (0,0,1)T, respectively, · stands for the diagonal matrix corresponding to subscript α, Im denotes the imagery part of a complex number, and superscript T denotes the matrix transpose.


### 6.1.3. Coordinate Transformation

A typical composite laminate consists of individual layers, which are usually made of unidirectional plies with the same or regularly alternating orientation. A layer is generally referred to the global coordinate framex1,x2, and x3 of the structural element rather than to coordinates x,y, and z associated with the ply orientation. So it is necessary to transform the constitutive relationship of each layer from the material coordinate frame to the uniform global coordinate frame.For the two coordinate systems mentioned in Figure7, the angle between the axis-x1 and axis-x is denoted by φ, which is positive along the anticlockwise direction, and the relationship for transformation of stress and strain components between the local material coordinates and global coordinates is given by(177)σ11,σ22,σ23,σ31,σ12T=T-1σxx,σyy,σyz,σzx,σxyT,εxx,εyy,εyz,εzx,εxyT=T-1Tε11,ε22,ε23,ε31,ε12T,where the transformation matrix T and its inverse matrix are defined as(178)T=c2s2002css2c200-2cs00c-s000sc0-cscs00c2-s2,T-1=c2s200-2css2c2002cs00cs000-sc0cs-cs00c2-s2with c=cos⁡(φ),s=sin(φ). Therefore, the constitutive relationship in the global coordinate system is given by (179)σ11,σ22,σ23,σ31,σ12T=T-1CT-1Tε11,ε22,ε23,ε31,ε12T.Figure 7
Schematic of the relationship between global coordinate system (x1,x2) and local material coordinate system (x,y).


## 6.1.1. Basic Equations and Stroh Formalism

In the Cartesian coordinate system (x1,x2,x3), if we neglect the body force bi, the equilibrium equations, stress-strain laws, and strain-displacement equations for anisotropic elasticity are [61] (162)σij,j=0,(163)σij=Cijklekl,(164)eij=12ui,j+uj,i,where i,j=1,2,3, Cijkl is the fourth-rank anisotropic elasticity tensor. The equilibrium equations can be rewritten in terms of displacements by substituting (163) and (164) into (162) as (165)Cijkluk,jl=0.The boundary conditions of the boundary value problem (163)–(165) are (166)ui=u-ion⁡Γu,ti=σijnj=t-ion⁡Γt,where u-i and t-i are the prescribed boundary displacement and traction vector, respectively. In addition, ni is the unit outward normal to the boundary and Γ=Γu+Γt is the boundary of the solution domain Ω.For the generalized two-dimensional deformation of anisotropic elasticityui is assumed to depend on x1 and x2 only. Based on this assumption, the general solution to (165) can be written as [61, 62](167)u=2ReAfz,φ=2ReBfz,where u=(u1,u2,u3)T is the displacement vector, φ=(φ1,φ2,φ3)T is the stress function vector, and f(z)=[f1(z1),f2(z2),f3(z3)]T is a function vector composed of three holomorphic complex functions fα(zα),α=1,2,3, which is an arbitrary function with argument zα=x1+pαx2 and will be determined by satisfying the boundary and loading conditions of a given problem. In (167), Re stands for the real part of a complex number, pα are the material eigenvalues with positive imaginary part, and A=[a1,a2,a3] and B=[b1,b2,b3] are 3 × 3 complex matrices formed by the material eigenvector associated with pα, which can be obtained by the following Eigen relations [61]:(168)Nξ=pξ,where N is a 6 × 6 foundational elasticity matrix and ξ is a 6 × 1 column vector defined by(169)N=N1N2N3N1T,ξ=ab,where N1=-T-1RT,N2=T-1,N3=RT-1RT-Q, and the matrices Q,R and T are 3 × 3 matrices extracted from Cijkl as follows:(170)Qik=Ci1k1,Rik=Ci1k2,Tik=Ci2k2.The stresses can be obtained from the derivation of stress functions φ as follows:(171)σi1=2ReLf′z,σi2=2ReBf′z,where(172)L=-p1b1,-p2b2,-p3b3,-p4b4.


## 6.1.2. Foundational Solutions

To find the fundamental solution needed in our analysis, we have to first derive the Green’s function of the problem: an infinite homogeneous anisotropic elastic medium loaded by a concentrated point force (or line force for two-dimensional problems)p^=(p^1,p^2,p^3) applied at an internal point x^=(x^1,x^2) far from the boundary. The boundary conditions of this problem can be written as (173)∫Cdϕ=p^foranyclosedcurveCenclosingx^,∫Cdu=p^foranyclosedcurveC,limx→∞σij=0.Thus, the Green’s function satisfying the above boundary conditions is found to be [72] (174)fz=12πiln⁡zα-z^αATp^.Therefore, fundamental solutions of the problem can be expressed as(175)u=1πImAln⁡⁡zα-z^αATp^,ϕ=1πImBln⁡⁡zα-z^αATp^.The corresponding stress components can be obtained from stress function ϕ as(176)σi1*=-ϕ,2=-1πImBpαzα-z^αATp^,σi2*=ϕ,1=1πImB1zα-z^αATp^,where p^ are chosen to be (1,0,0)T,(0,1,0)T, and (0,0,1)T, respectively, · stands for the diagonal matrix corresponding to subscript α, Im denotes the imagery part of a complex number, and superscript T denotes the matrix transpose.


## 6.1.3. Coordinate Transformation

A typical composite laminate consists of individual layers, which are usually made of unidirectional plies with the same or regularly alternating orientation. A layer is generally referred to the global coordinate framex1,x2, and x3 of the structural element rather than to coordinates x,y, and z associated with the ply orientation. So it is necessary to transform the constitutive relationship of each layer from the material coordinate frame to the uniform global coordinate frame.For the two coordinate systems mentioned in Figure7, the angle between the axis-x1 and axis-x is denoted by φ, which is positive along the anticlockwise direction, and the relationship for transformation of stress and strain components between the local material coordinates and global coordinates is given by(177)σ11,σ22,σ23,σ31,σ12T=T-1σxx,σyy,σyz,σzx,σxyT,εxx,εyy,εyz,εzx,εxyT=T-1Tε11,ε22,ε23,ε31,ε12T,where the transformation matrix T and its inverse matrix are defined as(178)T=c2s2002css2c200-2cs00c-s000sc0-cscs00c2-s2,T-1=c2s200-2css2c2002cs00cs000-sc0cs-cs00c2-s2with c=cos⁡(φ),s=sin(φ). Therefore, the constitutive relationship in the global coordinate system is given by (179)σ11,σ22,σ23,σ31,σ12T=T-1CT-1Tε11,ε22,ε23,ε31,ε12T.Figure 7
Schematic of the relationship between global coordinate system (x1,x2) and local material coordinate system (x,y).


## 6.2. Formulations of HFS-FEM

### 6.2.1. Assumed Independent Fields

The intraelement displacement fields for a particular elemente is approximated in terms of a linear combination of fundamental solutions of the problem as (180)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as (181)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nc2nc3nTin which ns is the number of source points, x and ysj are, respectively, the field point and source point in the coordinate system (x1,x2) local to the element under consideration. The fundamental solution uij*(x,ysj) is given by (175) for general elements. γ is a dimensionless coefficient for determining source points as described in previous sections.The corresponding stress fields can be expressed as(182)σx=σ11σ22σ23σ31σ12T=Tece,where (183)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ231*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by (176) when p^i is selected to be (1,0,0)T, (0,1,0)T, and (0,0,1)T, respectively. As a consequence, the traction can be written as (184)t1t2t3=nσ=Qecein which (185)Qe=nTe,n=n100n3n20n2n30n100n2n10.The unknown ce in (180) and (182) may be calculated using a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as in conventional FEM (see Figure 1). Thus, the frame field is defined as (186)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N~e and de can be expressed as (187)N~e=00N-1N-2N-3000,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (188)N-i=N~i000N~i000N~i,0=000000000and N~1,N~2, and N~3 are expressed by natural coordinate ξ∈[-1,1](189)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


### 6.2.2. Modified Functional for HFS-FEM

With the assumption of two distinct intraelement field and frame field for elements, we can establish the modified variational principle based on (165) and (166) for the hybrid finite element method of anisotropic materials [21, 73]. In the absence of the body forces, the hybrid variational functional Πme for a particular element e is constructed as(190)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ,where the boundary Γe of the element e is (191)Γe=Γeu∪Γet∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γtand ΓeI is the interelement boundary of element e. Performing a variation of Πm, one obtains(192)δΠme=∬Ωeσijδui,jdΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~i-δuidΓ.Applying Gaussian theorem (193)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force (194)ti=σijnjwe obtain (195)δΠe=-∬Ωeσij,jδuidΩ+∫ΓetiδuidΓ-∫Γett-iδu~idΓ.Then, substituting (195) into (192) gives(196)δΠme=-∬Ωeσij,jδuidΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~idΓ.Considering the fact that(197)∫Γeutiδu~idΓ=0we can finally obtain the following form:(198)δΠme=-∬Ωeσij,jδuidΩ+∫Γetti-t-iδu~idΓ+∫Γeu~i-uiδtidΓ+∫ΓItiδu~idΓ.Therefore, the Euler equations for (190) result in (165) and (166) because the quantities δui,δti, and δu~i may be arbitrary. As for the continuity condition between elements, it can be easily seen from the following variational of two adjacent elements such as e and f(199)δΠm(e∪f)=-∬Ωe∪Ωfσij,jδuidΩ+∫Γet+Γetti-t-iδu~idΓ+∫Γe+Γfu~i-uiδtidΓ+∫ΓefItie+tifδu~i.This indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (190) can be easily proved by the so-called “second variational approach” as well, which indicates functional (190) has a local extreme. Therefore, we can conclude that the variational functional (190) can be used for deriving hybrid finite element formulations.


### 6.2.3. Element Stiffness Equation

Using Gaussian theorem and equilibrium equations, all domain integrals in (190) can be converted into boundary integrals as follows:(200)Πme=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.Substituting (180), (184), and (186) into the functional (200) yields the formulation as (201)Πme=-12ceTHece+ceTGede-deTge,where(202)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and  (34).


## 6.2.1. Assumed Independent Fields

The intraelement displacement fields for a particular elemente is approximated in terms of a linear combination of fundamental solutions of the problem as (180)ux=u1xu2xu3x=Necex∈Ωe,ysj∉Ωe,where the matrix Ne and unknown vector ce can be further written as (181)Ne=u11*x,ys1u12*x,ys1u13*x,ys1⋯u11*x,ysnsu12*x,ysnsu13*x,ysnsu12*x,ys1u22*x,ys1u23*x,ys1⋯u12*x,ysnsu22*x,ysnsu23*x,ysnsu13*x,ys1u32*x,ys1u33*x,ys1⋯u13*x,ysnsu32*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nc2nc3nTin which ns is the number of source points, x and ysj are, respectively, the field point and source point in the coordinate system (x1,x2) local to the element under consideration. The fundamental solution uij*(x,ysj) is given by (175) for general elements. γ is a dimensionless coefficient for determining source points as described in previous sections.The corresponding stress fields can be expressed as(182)σx=σ11σ22σ23σ31σ12T=Tece,where (183)Te=σ111*x,y1σ211*x,y1σ311*x,y1⋯σ111*x,ynsσ211*x,ynsσ311*x,ynsσ122*x,y1σ222*x,y1σ322*x,y1⋯σ122*x,ynsσ222*x,ynsσ322*x,ynsσ123*x,y1σ223*x,y1σ323*x,y1⋯σ123*x,ynsσ223*x,ynsσ323*x,ynsσ131*x,y1σ231*x,y1σ331*x,y1⋯σ131*x,ynsσ231*x,ynsσ331*x,ynsσ112*x,y1σ231*x,y1σ312*x,y1⋯σ112*x,ynsσ212*x,ynsσ312*x,yns.The components σijk*(x,y) are given by (176) when p^i is selected to be (1,0,0)T, (0,1,0)T, and (0,0,1)T, respectively. As a consequence, the traction can be written as (184)t1t2t3=nσ=Qecein which (185)Qe=nTe,n=n100n3n20n2n30n100n2n10.The unknown ce in (180) and (182) may be calculated using a hybrid technique [31], in which the elements are linked through an auxiliary conforming displacement frame which has the same form as in conventional FEM (see Figure 1). Thus, the frame field is defined as (186)u~x=u~1u~2u~3=N~1N~2N~3de=N~ede,x∈Γe,where the symbol “~” is used to specify that the field is defined on the element boundary only, N~e is the matrix of shape functions, and de is the nodal displacements of elements. Taking the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N~e and de can be expressed as (187)N~e=00N-1N-2N-3000,de=u11u21u31u12u22u32⋯u18u28u38T,where the shape functions are expressed as (188)N-i=N~i000N~i000N~i,0=000000000and N~1,N~2, and N~3 are expressed by natural coordinate ξ∈[-1,1](189)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


## 6.2.2. Modified Functional for HFS-FEM

With the assumption of two distinct intraelement field and frame field for elements, we can establish the modified variational principle based on (165) and (166) for the hybrid finite element method of anisotropic materials [21, 73]. In the absence of the body forces, the hybrid variational functional Πme for a particular element e is constructed as(190)Πme=12∬ΩeσijεijdΩ-∫Γtt-iu~idΓ+∫Γetiu~i-uidΓ,where the boundary Γe of the element e is (191)Γe=Γeu∪Γet∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γtand ΓeI is the interelement boundary of element e. Performing a variation of Πm, one obtains(192)δΠme=∬Ωeσijδui,jdΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~i-δuidΓ.Applying Gaussian theorem (193)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force (194)ti=σijnjwe obtain (195)δΠe=-∬Ωeσij,jδuidΩ+∫ΓetiδuidΓ-∫Γett-iδu~idΓ.Then, substituting (195) into (192) gives(196)δΠme=-∬Ωeσij,jδuidΩ-∫Γett-iδu~idΓ+∫Γeu~i-uiδti+tiδu~idΓ.Considering the fact that(197)∫Γeutiδu~idΓ=0we can finally obtain the following form:(198)δΠme=-∬Ωeσij,jδuidΩ+∫Γetti-t-iδu~idΓ+∫Γeu~i-uiδtidΓ+∫ΓItiδu~idΓ.Therefore, the Euler equations for (190) result in (165) and (166) because the quantities δui,δti, and δu~i may be arbitrary. As for the continuity condition between elements, it can be easily seen from the following variational of two adjacent elements such as e and f(199)δΠm(e∪f)=-∬Ωe∪Ωfσij,jδuidΩ+∫Γet+Γetti-t-iδu~idΓ+∫Γe+Γfu~i-uiδtidΓ+∫ΓefItie+tifδu~i.This indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (190) can be easily proved by the so-called “second variational approach” as well, which indicates functional (190) has a local extreme. Therefore, we can conclude that the variational functional (190) can be used for deriving hybrid finite element formulations.


## 6.2.3. Element Stiffness Equation

Using Gaussian theorem and equilibrium equations, all domain integrals in (190) can be converted into boundary integrals as follows:(200)Πme=-12∫ΓetiuidΓ+∫Γetiu~idΓ-∫Γtt-iu~idΓ.Substituting (180), (184), and (186) into the functional (200) yields the formulation as (201)Πme=-12ceTHece+ceTGede-deTge,where(202)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and  (34).


## 7. Piezoelectric Materials

Piezoelectric materials have the property of converting electrical energy into mechanical energy and vice versa. This reciprocity in the energy conversion makes them very attractive for using in electromechanical devices, such as sensors, actuators, transducers, and frequency generators. To enhance understanding of the electromechanical coupling mechanism in piezoelectric materials and to explore their potential applications in practical engineering, numerous investigations, either analytically or numerically, have been reported over the past decades [73–81]. In this section, the HFS-FEM is developed for modeling two-dimensional piezoelectric material [82, 83]. The detailed formulations and procedures are given in the following sections.


### 7.1. Basic Equations for Piezoelectric Materials

For a linear piezoelectric material in absence of body forces and electric charge density, the differential governing equations in the Cartesian coordinate systemxi (i=1,2,3) are given by (203)σij,j=0,Di,i=0inΩ,where σij is the stress tensor, Di is the electric displacement vector, and Ω is the solution domain. With strain and electric field as the independent variables, the constitutive equations are written as (204)σij=cijklεkl-ekijEk,Di=eiklεkl+κikEk,where cijkl is the elasticity tensor measured under zero electric field, eikl and κij are, respectively, the piezoelectric tensor and dielectric tensor measured under zero strain. εij and Ei are the elastic strain tensor and the electric field vector, respectively. The relation between the strain tensor εij and the displacement ui is given by(205)εij=12ui,j+uj,iand the electric field component Ei is related to the electric potential ϕ by(206)Ei=-ϕ,i.The boundary conditions of the boundary value problem (203)–(206) can be defined by(207)ui=u-ion⁡Γu,(208)ti=σijnj=t-ion⁡Γt,(209)Dn=Dini=-q-n=D-non⁡ΓD,(210)ϕ=ϕ-on⁡Γϕ,where u-i,t-i,q-n, and ϕ- are, respectively, the prescribed boundary displacement, the prescribed traction vector, the prescribed surface charge, and the prescribed electric potential. In addition, ni is the unit outward normal to the boundary and (211)Γ=Γu+Γt=ΓD+Γϕis the boundary of the solution domain Ω.For the transversely isotropic material, ifx1-x2 is taken as the isotropic plane, one can employ either x1-x3 or x2-x3 plane to study the plane electromechanical phenomenon. Thus, choosing the former and considering the plane strain conditions (ε22=ε32=ε12=0 and E2=0) (204) can be reduced to(212)σ11σ33σ13=c11c130c13c33000c44ε11ε332ε13-0e310e33e150E1E3,D1D3=00e15e31e330ε11ε332ε13+κ1100κ33E1E3.For the plane stress piezoelectric problem (σ22=σ32=σ12=0 and D2=0), the constitutive equations can be obtained by replacing the coefficients c11,c13,c33,c44,e15,e31,e33,κ11, and κ33 in (212) as(213)c11*=c11-c122c11,c13*=c13-c12c13c11,c33*=c33-c132c11,c44*=c44,e15*=e15,e31*=e31-c12e31c11,e33*=e33-c13e31c11,κ11*=κ11,κ33*=κ33+e312c11.


### 7.2. Assumed Independent Fields

For the piezoelectric problems, HFS-FEM is based on assuming two distinct displacement and electric potential (DEP) fields: intraelement DEP fieldu and an independent DEP frame field u~ along element boundaries [21, 29].


#### 7.2.1. Intraelement Field

The intraelement DEP fieldu identically fulfills the governing differential equations (203) and is approximated by a linear combination of foundational solutions at different source points located outside of the element domain (214)ue=u1u2ϕ=∑j=1nsu11*x,ysju21*x,ysju31*x,ysju12*x,ysju22*x,ysju32*x,ysju13*x,ysju23*x,ysju33*x,ysj·c1jc2jc3j=Necex∈Ωe,ysj∉Ωe,where the fundamental solution matrix Ne is now given by (215)Ne=u11*x,ys1u21*x,ys1u31*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu12*x,ys1u22*x,ys1u32*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu13*x,ys1u23*x,ys1u33*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nsc2nsc3nsTin which x and ysj are, respectively, the field point (i.e., the nodal points of the element in this work) and source point in the local coordinate system (x1,x2). The component uij*(x,ysj) is the induced displacement component (i=1,2) or electric potential (i=3) in i-direction at the field point x due to a unit point load (j=1,2) or point charge (j=3) applied in j-direction at the source point ysj. The fundamental solution uij*(x,ysj) is given as [76, 82, 84](216)u11*=1πM11∑j=13sj1t2j11ln⁡⁡rj,u12*=1πM12∑j=13sj2t2j21arctgx1-xs1sjx3-xs3,u13*=1πM13∑j=13sj3t2j31arctgx1-xs1sjx3-xs3,u21*=1πM11∑j=13dj1t2j11arctgx1-xs1sjx3-xs3,u22*=1πM12∑j=13dj2t2j21ln⁡⁡rj,u23*=1πM13∑j=13dj3t2j31ln⁡⁡rj,u31*=1πM11∑j=13gj1t2j11arctgx1-xs1sjx3-xs3,u32*=1πM12∑j=13gj2t2j21ln⁡⁡rj,u33*=1πM13∑j=13gj3t2j31ln⁡⁡rj,where rj=(x1-xs1)2+sj2(x3-xs3)2 and sj is the three different roots of the characteristic equation asi6-bsi4+csi3-d=0. The source point ysj(j=1,2,…,ns) can be generated by the following method [36]: (217)ys=x0+γx0-xc.Making use of (205) and the expression of intraelement DEP field u in (214), the corresponding stress and electric displacement in (212) can be written as(218)σ=Tece,where σ=[σ11σ22σ12D1D2]T and (219)Te=σ11*x,ys1σ12*x,ys1σ13*x,ys1⋯σ11*x,ysnsσ12*x,ysnsσ13*x,ysnsσ21*x,ys1σ22*x,ys1σ23*x,ys1⋯σ21*x,ysnsσ22*x,ysnsσ23*x,ysnsσ31*x,ys1σ32*x,ys1σ33*x,ys1⋯σ31*x,ysnsσ32*x,ysnsσ33*x,ysnsσ41*x,ys1σ42*x,ys1σ43*x,ys1⋯σ41*x,ysnsσ42*x,ysnsσ43*x,ysnsσ51*x,ys1σ52*x,ys1σ53*x,ys1⋯σ51*x,ysnsσ52*x,ysnsσ53*x,ysnsin which σij*(x,ysj) denotes the corresponding stress components (i=1,2,3) or electric displacement (i=4,5) along i-direction at the field point x due to a unit point load (j=1,2) or a unit point charge (j=3) applied in j-direction at the source point ys and can be derived from (216) and are listed below, which are derived by substituting the fundamental solutions into constitutive equations [85]: (220)σ11*=1πM11∑j=13c11sj1-c13dj1sj-e31gj1sjx1-xs1rj2kkkkkkkkkkki·t2j11x1-xs1rj2,σ12*=1πM12∑j=13c11sj2+c13dj2sj+e31gj2sjsjx3-xs3rj2kkkkkkkkkki·t2j21sjx3-xs3rj2,σ13*=1πM13∑j=13c11sj3+c13dj3sj+e31gj3sjsjx3-xs3rj2kkkkkkkkkkk·t2j31sjx3-xs3rj2,σ21*=1πM11∑j=13c13sj1-c33dj1sj-e33gj1sjx1-xs1rj2kkkkkkkkkkkkkk·t2j11x1-xs1rj2,σ22*=1πM12∑j=13c13sj2+c33dj2sj+e33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ23*=-1πM13∑j=13c13sj3+c33dj3sj+e33gj3sjsjx3-xs3rj2kkkkkkkkkkkki·t2j31sjx3-xs3rj2,σ31*=1πM11∑j=13c44sj1sj+c44dj1+e15gj1sjx3-xs3rj2kkkkkkkkkkkkkk·t2j11sjx3-xs3rj2,σ32*=1πM12∑j=13-c44sj2sj+c44dj2+e15gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ33*=-1πM13∑j=13-c44sj3sj+c44dj3+e15gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ41*=1πM11∑j=13e15sj1sj+e15dj1-λ11gj1sjx3-xs3rj2kkkkkkkkkkk·t2j11sjx3-xs3rj2,σ42*=1πM12∑j=13-e15sj2sj+e15dj2-λ11gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ43*=-1πM13∑j=13-e15sj3sj+e15dj3-λ11gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ51*=1πM11∑j=13e31sj1-e33dj1sj+λ33gj1sjx1-xs1rj2kkkkkkkkkkkk·t2j11x1-xs1rj2,σ52*=1πM12∑j=13e31sj2+e33dj2sj-λ33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ53*=-1πM13∑j=13e31sj3+e33dj3sj-λ33gj3sjsjx3-xs3rj2kkkkkkkkkkkkk·t2j31sjx3-xs3rj2in which the coefficients sij,dij,gij,tij,M11,M12, and M13 are defined as in literature [84].From (204), (208), and (209), the generalized traction forces and electric displacement are given as (221)t1t2Dn=Q1Q2Q3ce=Qece,where(222)Qe=nTe,n=n10n2000n2n100000n1n2.


#### 7.2.2. Auxiliary Frame Field

For the two-dimensional piezoelectric problem under consideration, the frame field is assumed as(223)u~x=u~1u~2ϕ~=N~1N~2N~3de=N~ede,x∈Γe,where N~e is a matrix of the corresponding shape functions. For the side 3-4-5 of a particular quadratic element as shown in Figure 1, the shape function matrix N~e and nodal vector de can be given in the form (224)N~e=0⋯0N~100N~200N~3000⋯00⋯00N~100N~200N~300⋯00⋯0︸600N~100N~200N~30⋯0︸9,de=u11u21ϕ1⋯u14u24ϕ4⋯u18u28ϕ8T,where the shape functions N~e are expressed by natural coordinate ξ(225)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


### 7.3. HFS-FEM Formulations

#### 7.3.1. Variational Principles

Based on the assumption of two distinct DEP fields, the Euler equations of the proposed variational functional should also satisfy the following interelement continuity requirements in addition to (203)–(210) [82, 83]:(226)uie=uifϕe=ϕfon⁡Γe∩Γf,conformity,(227)tie+tif=0Dne+Dnf=0on⁡Γe∩Γf,reciprocity,where “e” and “f” stand for any two neighboring elements. Equations (203)–(210) together with (226) and (227) can now be taken as the basis to establish the modified variational principle for the hybrid finite element method of piezoelectric materials [21, 29].Since the stationary conditions of the traditional potential or complementary variational functional cannot satisfy the interelement continuity condition required in the proposed HFS-FEM, new modified variational functional should be developed. In the absence of the body forces and electric charge density, the hybrid functionalΠme for a particular element e is constructed as(228)Πme=Πe+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ,where(229)Πe=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓand the boundary Γe of the element e is (230)Γe=Γeu∪Γet∪ΓeI=Γeϕ∪ΓeD∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γt,Γeϕ=Γe∩Γϕ,ΓeD=Γe∩ΓDand ΓeI is the interelement boundary of element e. Compared to the functional employed in the conventional FEM, the present hybrid functional is constructed by adding two integral terms related to the intraelement and element frame DEP fields to guarantee the satisfaction of displacement and electrical potential continuity condition on the common boundary of two adjacent elements.It can be proved that the stationary conditions of the above functional (228) lead to (203)–(210). To this end, performing a variation of Πm, one obtains [86](231)δΠme=δΠe+∫Γeu~i-uiδti+tiδu~i-δuidΓ+∫Γeϕ~-ϕδDn+Dnδϕ~-δϕdΓin which the first term is(232)δΠe=∬ΩeσijδεijdΩ+∬ΩeDiδEidΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ=∬Ωeσijδui,jdΩ+∬ΩeDiδϕ,idΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Applying Gaussian theorem (233)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force and electrical displacement(234)ti=σijnj,Dn=Diniwe obtain (235)δΠe=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫ΓetiδuidΓ+∫ΓeDnδϕdΓ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Then, substituting (235) into (231) gives(236)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ+∫Γeu~i-uiδti+tiδu~idΓ+∫Γeϕ~-ϕδDn+Dnδϕ~dΓ.Considering the fact that(237)∫Γeutiδu~idΓ=0,∫ΓeϕDnδϕ~dΓ=0we finally get the following form:(238)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫Γet(ti-t-i)δu~idΓ+∫ΓeDDn-D-nδϕ~dΓ+∫Γeu~i-uiδtidΓ+∫Γeϕ~-ϕδDndΓ+∫ΓItiδu~idΓ+∫ΓIDnδϕ~dΓ.Therefore, the Euler equations for (238) result in (203)–(210) and (226) because the quantities δui,δti,δϕ,δDn,δu~i, and δϕ~ may be arbitrary. As for the continuity condition of (227), it can easily be seen from the following variation of two adjacent elements such as e and f:(239)δΠme∪f=-∬Ωe∪Ωfσij,jδuidΩ-∬Ωe∪ΩfDi,iδϕdΩ+∫Γet+Γetti-t-iδu~idΓ+∫ΓeD+ΓfDDn-D-nδϕ~dΓ+∫Γe+Γfu~i-uiδtidΓ+∫Γe+Γfϕ~-ϕδDndΓ+∫ΓefItie+tifδu~i+∫ΓefIDne+Dnfδϕ~dΓwhich indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (228) can be easily proved by the “second variational approach” as well, which indicates functional (228) has a local extreme.


#### 7.3.2. Element Stiffness Equation

The element stiffness equation can be generated by settingδΠme=0. To simplify the derivation, we first transform all domain integrals in (228) into boundary ones. With Gaussian theorem, the functional in (228) may be simplified as(240)Πme=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ+12∫ΓeDnϕdΓ-∬ΩeDi,iϕdΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for the HT finite element model:(241)Πme=12∫ΓetiuidΓ+12∫ΓeDnϕdΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=-12∫Γetiui+DnϕdΓ+∫Γe(tiu~i+Dnϕ~)dΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ.Substituting (214), (223), and (221) into the above functional (241) yields the formulation as(242)Πme=-12ceTHece+ceTGede-deTge,where (243)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ+∫ΓDN~eTD-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and (34). Following the procedure described in [21, 22], the missing rigid-body motion can be recovered by setting the augmented internal field of a particular element e as(244)ue=Nece+10x2001-x100001c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(245)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+ϕi-ϕ~i2which finally gives(246)c0=Re-1re,(247)Re=∑i=1n10x2i001-x1i0x2i-x1ix1i2+x2i200001,(248)re=∑i=1nΔue1iΔue2iΔue1ix2i-Δue2ix1iΔϕeiin which Δuei=(u~e-u^e)nodei and n is the number of element nodes. As a consequence, c0 can be calculated by (246) once the nodal DEP fields de and the interpolation coefficients ce are, respectively, determined by (214) and (223). Then the complete DEP fields ue can be obtained from (244).


### 7.4. Normalization

The order of magnitudes of the material constants and the corresponding field variables in piezoelectricity have a wide spectrum as large as 1019 in SI unit. This will lead to ill-conditioned matrix of the system [72, 87]. To resolve this problem, normalization of each quantity by its reference value is employed. The reference values for the stiffness, piezoelectric stress constant, dielectric constants, and strain are selected to be c0=1011N/m2, e0=101(N/mV), k0=10-9C/mV, and ε0=10-3V/m, respectively. The reference values of other quantities, as shown in Table 1, are determined in terms of these four fundamental reference variables and the characteristic length x0=100m of the problem so that the normalized governing equations remain in exactly the same form as the original equations.Table 1
Reference values for material constants and field variables in piezoelectricity derived from basic reference variables:c0, e0, k0, ε0, and x0.
Displacement
u
0
=
x
0
ε
0
=
1
0
-
3
(
m
)
Electric Potential
ϕ
0
=
x
0
E
0
=
1
0
7
(
V
)
Stress
σ
0
=
c
0
ε
0
=
1
0
8
(
N
/
m
2
)
Electric induction
D
0
=
k
0
E
0
=
1
0
-
2
(
C
/
m
2
)
Compliance
s
0
=
ε
0
/
σ
0
=
1
0
-
11
(
m
2
/
N
)
Impermeability
β
0
=
E
0
/
D
0
=
1
0
9
(
mV
/
C
)
Electric field
E
0
=
σ
0
/
e
0
=
1
0
7
(
V
/
m
)
Piezoelectric strain constant
g
0
=
E
0
/
σ
0
=
1
0
-
1
(
mV
/
N
)


## 7.1. Basic Equations for Piezoelectric Materials

For a linear piezoelectric material in absence of body forces and electric charge density, the differential governing equations in the Cartesian coordinate systemxi (i=1,2,3) are given by (203)σij,j=0,Di,i=0inΩ,where σij is the stress tensor, Di is the electric displacement vector, and Ω is the solution domain. With strain and electric field as the independent variables, the constitutive equations are written as (204)σij=cijklεkl-ekijEk,Di=eiklεkl+κikEk,where cijkl is the elasticity tensor measured under zero electric field, eikl and κij are, respectively, the piezoelectric tensor and dielectric tensor measured under zero strain. εij and Ei are the elastic strain tensor and the electric field vector, respectively. The relation between the strain tensor εij and the displacement ui is given by(205)εij=12ui,j+uj,iand the electric field component Ei is related to the electric potential ϕ by(206)Ei=-ϕ,i.The boundary conditions of the boundary value problem (203)–(206) can be defined by(207)ui=u-ion⁡Γu,(208)ti=σijnj=t-ion⁡Γt,(209)Dn=Dini=-q-n=D-non⁡ΓD,(210)ϕ=ϕ-on⁡Γϕ,where u-i,t-i,q-n, and ϕ- are, respectively, the prescribed boundary displacement, the prescribed traction vector, the prescribed surface charge, and the prescribed electric potential. In addition, ni is the unit outward normal to the boundary and (211)Γ=Γu+Γt=ΓD+Γϕis the boundary of the solution domain Ω.For the transversely isotropic material, ifx1-x2 is taken as the isotropic plane, one can employ either x1-x3 or x2-x3 plane to study the plane electromechanical phenomenon. Thus, choosing the former and considering the plane strain conditions (ε22=ε32=ε12=0 and E2=0) (204) can be reduced to(212)σ11σ33σ13=c11c130c13c33000c44ε11ε332ε13-0e310e33e150E1E3,D1D3=00e15e31e330ε11ε332ε13+κ1100κ33E1E3.For the plane stress piezoelectric problem (σ22=σ32=σ12=0 and D2=0), the constitutive equations can be obtained by replacing the coefficients c11,c13,c33,c44,e15,e31,e33,κ11, and κ33 in (212) as(213)c11*=c11-c122c11,c13*=c13-c12c13c11,c33*=c33-c132c11,c44*=c44,e15*=e15,e31*=e31-c12e31c11,e33*=e33-c13e31c11,κ11*=κ11,κ33*=κ33+e312c11.


## 7.2. Assumed Independent Fields

For the piezoelectric problems, HFS-FEM is based on assuming two distinct displacement and electric potential (DEP) fields: intraelement DEP fieldu and an independent DEP frame field u~ along element boundaries [21, 29].


### 7.2.1. Intraelement Field

The intraelement DEP fieldu identically fulfills the governing differential equations (203) and is approximated by a linear combination of foundational solutions at different source points located outside of the element domain (214)ue=u1u2ϕ=∑j=1nsu11*x,ysju21*x,ysju31*x,ysju12*x,ysju22*x,ysju32*x,ysju13*x,ysju23*x,ysju33*x,ysj·c1jc2jc3j=Necex∈Ωe,ysj∉Ωe,where the fundamental solution matrix Ne is now given by (215)Ne=u11*x,ys1u21*x,ys1u31*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu12*x,ys1u22*x,ys1u32*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu13*x,ys1u23*x,ys1u33*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nsc2nsc3nsTin which x and ysj are, respectively, the field point (i.e., the nodal points of the element in this work) and source point in the local coordinate system (x1,x2). The component uij*(x,ysj) is the induced displacement component (i=1,2) or electric potential (i=3) in i-direction at the field point x due to a unit point load (j=1,2) or point charge (j=3) applied in j-direction at the source point ysj. The fundamental solution uij*(x,ysj) is given as [76, 82, 84](216)u11*=1πM11∑j=13sj1t2j11ln⁡⁡rj,u12*=1πM12∑j=13sj2t2j21arctgx1-xs1sjx3-xs3,u13*=1πM13∑j=13sj3t2j31arctgx1-xs1sjx3-xs3,u21*=1πM11∑j=13dj1t2j11arctgx1-xs1sjx3-xs3,u22*=1πM12∑j=13dj2t2j21ln⁡⁡rj,u23*=1πM13∑j=13dj3t2j31ln⁡⁡rj,u31*=1πM11∑j=13gj1t2j11arctgx1-xs1sjx3-xs3,u32*=1πM12∑j=13gj2t2j21ln⁡⁡rj,u33*=1πM13∑j=13gj3t2j31ln⁡⁡rj,where rj=(x1-xs1)2+sj2(x3-xs3)2 and sj is the three different roots of the characteristic equation asi6-bsi4+csi3-d=0. The source point ysj(j=1,2,…,ns) can be generated by the following method [36]: (217)ys=x0+γx0-xc.Making use of (205) and the expression of intraelement DEP field u in (214), the corresponding stress and electric displacement in (212) can be written as(218)σ=Tece,where σ=[σ11σ22σ12D1D2]T and (219)Te=σ11*x,ys1σ12*x,ys1σ13*x,ys1⋯σ11*x,ysnsσ12*x,ysnsσ13*x,ysnsσ21*x,ys1σ22*x,ys1σ23*x,ys1⋯σ21*x,ysnsσ22*x,ysnsσ23*x,ysnsσ31*x,ys1σ32*x,ys1σ33*x,ys1⋯σ31*x,ysnsσ32*x,ysnsσ33*x,ysnsσ41*x,ys1σ42*x,ys1σ43*x,ys1⋯σ41*x,ysnsσ42*x,ysnsσ43*x,ysnsσ51*x,ys1σ52*x,ys1σ53*x,ys1⋯σ51*x,ysnsσ52*x,ysnsσ53*x,ysnsin which σij*(x,ysj) denotes the corresponding stress components (i=1,2,3) or electric displacement (i=4,5) along i-direction at the field point x due to a unit point load (j=1,2) or a unit point charge (j=3) applied in j-direction at the source point ys and can be derived from (216) and are listed below, which are derived by substituting the fundamental solutions into constitutive equations [85]: (220)σ11*=1πM11∑j=13c11sj1-c13dj1sj-e31gj1sjx1-xs1rj2kkkkkkkkkkki·t2j11x1-xs1rj2,σ12*=1πM12∑j=13c11sj2+c13dj2sj+e31gj2sjsjx3-xs3rj2kkkkkkkkkki·t2j21sjx3-xs3rj2,σ13*=1πM13∑j=13c11sj3+c13dj3sj+e31gj3sjsjx3-xs3rj2kkkkkkkkkkk·t2j31sjx3-xs3rj2,σ21*=1πM11∑j=13c13sj1-c33dj1sj-e33gj1sjx1-xs1rj2kkkkkkkkkkkkkk·t2j11x1-xs1rj2,σ22*=1πM12∑j=13c13sj2+c33dj2sj+e33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ23*=-1πM13∑j=13c13sj3+c33dj3sj+e33gj3sjsjx3-xs3rj2kkkkkkkkkkkki·t2j31sjx3-xs3rj2,σ31*=1πM11∑j=13c44sj1sj+c44dj1+e15gj1sjx3-xs3rj2kkkkkkkkkkkkkk·t2j11sjx3-xs3rj2,σ32*=1πM12∑j=13-c44sj2sj+c44dj2+e15gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ33*=-1πM13∑j=13-c44sj3sj+c44dj3+e15gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ41*=1πM11∑j=13e15sj1sj+e15dj1-λ11gj1sjx3-xs3rj2kkkkkkkkkkk·t2j11sjx3-xs3rj2,σ42*=1πM12∑j=13-e15sj2sj+e15dj2-λ11gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ43*=-1πM13∑j=13-e15sj3sj+e15dj3-λ11gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ51*=1πM11∑j=13e31sj1-e33dj1sj+λ33gj1sjx1-xs1rj2kkkkkkkkkkkk·t2j11x1-xs1rj2,σ52*=1πM12∑j=13e31sj2+e33dj2sj-λ33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ53*=-1πM13∑j=13e31sj3+e33dj3sj-λ33gj3sjsjx3-xs3rj2kkkkkkkkkkkkk·t2j31sjx3-xs3rj2in which the coefficients sij,dij,gij,tij,M11,M12, and M13 are defined as in literature [84].From (204), (208), and (209), the generalized traction forces and electric displacement are given as (221)t1t2Dn=Q1Q2Q3ce=Qece,where(222)Qe=nTe,n=n10n2000n2n100000n1n2.


### 7.2.2. Auxiliary Frame Field

For the two-dimensional piezoelectric problem under consideration, the frame field is assumed as(223)u~x=u~1u~2ϕ~=N~1N~2N~3de=N~ede,x∈Γe,where N~e is a matrix of the corresponding shape functions. For the side 3-4-5 of a particular quadratic element as shown in Figure 1, the shape function matrix N~e and nodal vector de can be given in the form (224)N~e=0⋯0N~100N~200N~3000⋯00⋯00N~100N~200N~300⋯00⋯0︸600N~100N~200N~30⋯0︸9,de=u11u21ϕ1⋯u14u24ϕ4⋯u18u28ϕ8T,where the shape functions N~e are expressed by natural coordinate ξ(225)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


## 7.2.1. Intraelement Field

The intraelement DEP fieldu identically fulfills the governing differential equations (203) and is approximated by a linear combination of foundational solutions at different source points located outside of the element domain (214)ue=u1u2ϕ=∑j=1nsu11*x,ysju21*x,ysju31*x,ysju12*x,ysju22*x,ysju32*x,ysju13*x,ysju23*x,ysju33*x,ysj·c1jc2jc3j=Necex∈Ωe,ysj∉Ωe,where the fundamental solution matrix Ne is now given by (215)Ne=u11*x,ys1u21*x,ys1u31*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu12*x,ys1u22*x,ys1u32*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysnsu13*x,ys1u23*x,ys1u33*x,ys1⋯u33*x,ysnsu33*x,ysnsu33*x,ysns,ce=c11c21c31⋯c1nsc2nsc3nsTin which x and ysj are, respectively, the field point (i.e., the nodal points of the element in this work) and source point in the local coordinate system (x1,x2). The component uij*(x,ysj) is the induced displacement component (i=1,2) or electric potential (i=3) in i-direction at the field point x due to a unit point load (j=1,2) or point charge (j=3) applied in j-direction at the source point ysj. The fundamental solution uij*(x,ysj) is given as [76, 82, 84](216)u11*=1πM11∑j=13sj1t2j11ln⁡⁡rj,u12*=1πM12∑j=13sj2t2j21arctgx1-xs1sjx3-xs3,u13*=1πM13∑j=13sj3t2j31arctgx1-xs1sjx3-xs3,u21*=1πM11∑j=13dj1t2j11arctgx1-xs1sjx3-xs3,u22*=1πM12∑j=13dj2t2j21ln⁡⁡rj,u23*=1πM13∑j=13dj3t2j31ln⁡⁡rj,u31*=1πM11∑j=13gj1t2j11arctgx1-xs1sjx3-xs3,u32*=1πM12∑j=13gj2t2j21ln⁡⁡rj,u33*=1πM13∑j=13gj3t2j31ln⁡⁡rj,where rj=(x1-xs1)2+sj2(x3-xs3)2 and sj is the three different roots of the characteristic equation asi6-bsi4+csi3-d=0. The source point ysj(j=1,2,…,ns) can be generated by the following method [36]: (217)ys=x0+γx0-xc.Making use of (205) and the expression of intraelement DEP field u in (214), the corresponding stress and electric displacement in (212) can be written as(218)σ=Tece,where σ=[σ11σ22σ12D1D2]T and (219)Te=σ11*x,ys1σ12*x,ys1σ13*x,ys1⋯σ11*x,ysnsσ12*x,ysnsσ13*x,ysnsσ21*x,ys1σ22*x,ys1σ23*x,ys1⋯σ21*x,ysnsσ22*x,ysnsσ23*x,ysnsσ31*x,ys1σ32*x,ys1σ33*x,ys1⋯σ31*x,ysnsσ32*x,ysnsσ33*x,ysnsσ41*x,ys1σ42*x,ys1σ43*x,ys1⋯σ41*x,ysnsσ42*x,ysnsσ43*x,ysnsσ51*x,ys1σ52*x,ys1σ53*x,ys1⋯σ51*x,ysnsσ52*x,ysnsσ53*x,ysnsin which σij*(x,ysj) denotes the corresponding stress components (i=1,2,3) or electric displacement (i=4,5) along i-direction at the field point x due to a unit point load (j=1,2) or a unit point charge (j=3) applied in j-direction at the source point ys and can be derived from (216) and are listed below, which are derived by substituting the fundamental solutions into constitutive equations [85]: (220)σ11*=1πM11∑j=13c11sj1-c13dj1sj-e31gj1sjx1-xs1rj2kkkkkkkkkkki·t2j11x1-xs1rj2,σ12*=1πM12∑j=13c11sj2+c13dj2sj+e31gj2sjsjx3-xs3rj2kkkkkkkkkki·t2j21sjx3-xs3rj2,σ13*=1πM13∑j=13c11sj3+c13dj3sj+e31gj3sjsjx3-xs3rj2kkkkkkkkkkk·t2j31sjx3-xs3rj2,σ21*=1πM11∑j=13c13sj1-c33dj1sj-e33gj1sjx1-xs1rj2kkkkkkkkkkkkkk·t2j11x1-xs1rj2,σ22*=1πM12∑j=13c13sj2+c33dj2sj+e33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ23*=-1πM13∑j=13c13sj3+c33dj3sj+e33gj3sjsjx3-xs3rj2kkkkkkkkkkkki·t2j31sjx3-xs3rj2,σ31*=1πM11∑j=13c44sj1sj+c44dj1+e15gj1sjx3-xs3rj2kkkkkkkkkkkkkk·t2j11sjx3-xs3rj2,σ32*=1πM12∑j=13-c44sj2sj+c44dj2+e15gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ33*=-1πM13∑j=13-c44sj3sj+c44dj3+e15gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ41*=1πM11∑j=13e15sj1sj+e15dj1-λ11gj1sjx3-xs3rj2kkkkkkkkkkk·t2j11sjx3-xs3rj2,σ42*=1πM12∑j=13-e15sj2sj+e15dj2-λ11gj2x1-xs1rj2kkkkkkkkkkk·t2j21x1-xs1rj2,σ43*=-1πM13∑j=13-e15sj3sj+e15dj3-λ11gj3x1-xs1rj2kkkkkkkkkkkkk·t2j31x1-xs1rj2,σ51*=1πM11∑j=13e31sj1-e33dj1sj+λ33gj1sjx1-xs1rj2kkkkkkkkkkkk·t2j11x1-xs1rj2,σ52*=1πM12∑j=13e31sj2+e33dj2sj-λ33gj2sjsjx3-xs3rj2kkkkkkkkkkk·t2j21sjx3-xs3rj2,σ53*=-1πM13∑j=13e31sj3+e33dj3sj-λ33gj3sjsjx3-xs3rj2kkkkkkkkkkkkk·t2j31sjx3-xs3rj2in which the coefficients sij,dij,gij,tij,M11,M12, and M13 are defined as in literature [84].From (204), (208), and (209), the generalized traction forces and electric displacement are given as (221)t1t2Dn=Q1Q2Q3ce=Qece,where(222)Qe=nTe,n=n10n2000n2n100000n1n2.


## 7.2.2. Auxiliary Frame Field

For the two-dimensional piezoelectric problem under consideration, the frame field is assumed as(223)u~x=u~1u~2ϕ~=N~1N~2N~3de=N~ede,x∈Γe,where N~e is a matrix of the corresponding shape functions. For the side 3-4-5 of a particular quadratic element as shown in Figure 1, the shape function matrix N~e and nodal vector de can be given in the form (224)N~e=0⋯0N~100N~200N~3000⋯00⋯00N~100N~200N~300⋯00⋯0︸600N~100N~200N~30⋯0︸9,de=u11u21ϕ1⋯u14u24ϕ4⋯u18u28ϕ8T,where the shape functions N~e are expressed by natural coordinate ξ(225)N~1=-ξ1-ξ2,N~2=1-ξ2,N~3=ξ1+ξ2ξ∈-1,1.


## 7.3. HFS-FEM Formulations

### 7.3.1. Variational Principles

Based on the assumption of two distinct DEP fields, the Euler equations of the proposed variational functional should also satisfy the following interelement continuity requirements in addition to (203)–(210) [82, 83]:(226)uie=uifϕe=ϕfon⁡Γe∩Γf,conformity,(227)tie+tif=0Dne+Dnf=0on⁡Γe∩Γf,reciprocity,where “e” and “f” stand for any two neighboring elements. Equations (203)–(210) together with (226) and (227) can now be taken as the basis to establish the modified variational principle for the hybrid finite element method of piezoelectric materials [21, 29].Since the stationary conditions of the traditional potential or complementary variational functional cannot satisfy the interelement continuity condition required in the proposed HFS-FEM, new modified variational functional should be developed. In the absence of the body forces and electric charge density, the hybrid functionalΠme for a particular element e is constructed as(228)Πme=Πe+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ,where(229)Πe=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓand the boundary Γe of the element e is (230)Γe=Γeu∪Γet∪ΓeI=Γeϕ∪ΓeD∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γt,Γeϕ=Γe∩Γϕ,ΓeD=Γe∩ΓDand ΓeI is the interelement boundary of element e. Compared to the functional employed in the conventional FEM, the present hybrid functional is constructed by adding two integral terms related to the intraelement and element frame DEP fields to guarantee the satisfaction of displacement and electrical potential continuity condition on the common boundary of two adjacent elements.It can be proved that the stationary conditions of the above functional (228) lead to (203)–(210). To this end, performing a variation of Πm, one obtains [86](231)δΠme=δΠe+∫Γeu~i-uiδti+tiδu~i-δuidΓ+∫Γeϕ~-ϕδDn+Dnδϕ~-δϕdΓin which the first term is(232)δΠe=∬ΩeσijδεijdΩ+∬ΩeDiδEidΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ=∬Ωeσijδui,jdΩ+∬ΩeDiδϕ,idΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Applying Gaussian theorem (233)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force and electrical displacement(234)ti=σijnj,Dn=Diniwe obtain (235)δΠe=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫ΓetiδuidΓ+∫ΓeDnδϕdΓ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Then, substituting (235) into (231) gives(236)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ+∫Γeu~i-uiδti+tiδu~idΓ+∫Γeϕ~-ϕδDn+Dnδϕ~dΓ.Considering the fact that(237)∫Γeutiδu~idΓ=0,∫ΓeϕDnδϕ~dΓ=0we finally get the following form:(238)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫Γet(ti-t-i)δu~idΓ+∫ΓeDDn-D-nδϕ~dΓ+∫Γeu~i-uiδtidΓ+∫Γeϕ~-ϕδDndΓ+∫ΓItiδu~idΓ+∫ΓIDnδϕ~dΓ.Therefore, the Euler equations for (238) result in (203)–(210) and (226) because the quantities δui,δti,δϕ,δDn,δu~i, and δϕ~ may be arbitrary. As for the continuity condition of (227), it can easily be seen from the following variation of two adjacent elements such as e and f:(239)δΠme∪f=-∬Ωe∪Ωfσij,jδuidΩ-∬Ωe∪ΩfDi,iδϕdΩ+∫Γet+Γetti-t-iδu~idΓ+∫ΓeD+ΓfDDn-D-nδϕ~dΓ+∫Γe+Γfu~i-uiδtidΓ+∫Γe+Γfϕ~-ϕδDndΓ+∫ΓefItie+tifδu~i+∫ΓefIDne+Dnfδϕ~dΓwhich indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (228) can be easily proved by the “second variational approach” as well, which indicates functional (228) has a local extreme.


### 7.3.2. Element Stiffness Equation

The element stiffness equation can be generated by settingδΠme=0. To simplify the derivation, we first transform all domain integrals in (228) into boundary ones. With Gaussian theorem, the functional in (228) may be simplified as(240)Πme=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ+12∫ΓeDnϕdΓ-∬ΩeDi,iϕdΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for the HT finite element model:(241)Πme=12∫ΓetiuidΓ+12∫ΓeDnϕdΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=-12∫Γetiui+DnϕdΓ+∫Γe(tiu~i+Dnϕ~)dΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ.Substituting (214), (223), and (221) into the above functional (241) yields the formulation as(242)Πme=-12ceTHece+ceTGede-deTge,where (243)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ+∫ΓDN~eTD-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and (34). Following the procedure described in [21, 22], the missing rigid-body motion can be recovered by setting the augmented internal field of a particular element e as(244)ue=Nece+10x2001-x100001c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(245)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+ϕi-ϕ~i2which finally gives(246)c0=Re-1re,(247)Re=∑i=1n10x2i001-x1i0x2i-x1ix1i2+x2i200001,(248)re=∑i=1nΔue1iΔue2iΔue1ix2i-Δue2ix1iΔϕeiin which Δuei=(u~e-u^e)nodei and n is the number of element nodes. As a consequence, c0 can be calculated by (246) once the nodal DEP fields de and the interpolation coefficients ce are, respectively, determined by (214) and (223). Then the complete DEP fields ue can be obtained from (244).


## 7.3.1. Variational Principles

Based on the assumption of two distinct DEP fields, the Euler equations of the proposed variational functional should also satisfy the following interelement continuity requirements in addition to (203)–(210) [82, 83]:(226)uie=uifϕe=ϕfon⁡Γe∩Γf,conformity,(227)tie+tif=0Dne+Dnf=0on⁡Γe∩Γf,reciprocity,where “e” and “f” stand for any two neighboring elements. Equations (203)–(210) together with (226) and (227) can now be taken as the basis to establish the modified variational principle for the hybrid finite element method of piezoelectric materials [21, 29].Since the stationary conditions of the traditional potential or complementary variational functional cannot satisfy the interelement continuity condition required in the proposed HFS-FEM, new modified variational functional should be developed. In the absence of the body forces and electric charge density, the hybrid functionalΠme for a particular element e is constructed as(228)Πme=Πe+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ,where(229)Πe=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓand the boundary Γe of the element e is (230)Γe=Γeu∪Γet∪ΓeI=Γeϕ∪ΓeD∪ΓeI,Γeu=Γe∩Γu,Γet=Γe∩Γt,Γeϕ=Γe∩Γϕ,ΓeD=Γe∩ΓDand ΓeI is the interelement boundary of element e. Compared to the functional employed in the conventional FEM, the present hybrid functional is constructed by adding two integral terms related to the intraelement and element frame DEP fields to guarantee the satisfaction of displacement and electrical potential continuity condition on the common boundary of two adjacent elements.It can be proved that the stationary conditions of the above functional (228) lead to (203)–(210). To this end, performing a variation of Πm, one obtains [86](231)δΠme=δΠe+∫Γeu~i-uiδti+tiδu~i-δuidΓ+∫Γeϕ~-ϕδDn+Dnδϕ~-δϕdΓin which the first term is(232)δΠe=∬ΩeσijδεijdΩ+∬ΩeDiδEidΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ=∬Ωeσijδui,jdΩ+∬ΩeDiδϕ,idΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Applying Gaussian theorem (233)∬Ωef,idΩ=∫Γef·nidΓand the definitions of traction force and electrical displacement(234)ti=σijnj,Dn=Diniwe obtain (235)δΠe=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫ΓetiδuidΓ+∫ΓeDnδϕdΓ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ.Then, substituting (235) into (231) gives(236)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ-∫Γett-iδu~idΓ-∫ΓeDD-nδϕ~dΓ+∫Γeu~i-uiδti+tiδu~idΓ+∫Γeϕ~-ϕδDn+Dnδϕ~dΓ.Considering the fact that(237)∫Γeutiδu~idΓ=0,∫ΓeϕDnδϕ~dΓ=0we finally get the following form:(238)δΠme=-∬Ωeσij,jδuidΩ-∬ΩeDi,iδϕdΩ+∫Γet(ti-t-i)δu~idΓ+∫ΓeDDn-D-nδϕ~dΓ+∫Γeu~i-uiδtidΓ+∫Γeϕ~-ϕδDndΓ+∫ΓItiδu~idΓ+∫ΓIDnδϕ~dΓ.Therefore, the Euler equations for (238) result in (203)–(210) and (226) because the quantities δui,δti,δϕ,δDn,δu~i, and δϕ~ may be arbitrary. As for the continuity condition of (227), it can easily be seen from the following variation of two adjacent elements such as e and f:(239)δΠme∪f=-∬Ωe∪Ωfσij,jδuidΩ-∬Ωe∪ΩfDi,iδϕdΩ+∫Γet+Γetti-t-iδu~idΓ+∫ΓeD+ΓfDDn-D-nδϕ~dΓ+∫Γe+Γfu~i-uiδtidΓ+∫Γe+Γfϕ~-ϕδDndΓ+∫ΓefItie+tifδu~i+∫ΓefIDne+Dnfδϕ~dΓwhich indicates that the stationary condition of the functional satisfies both the required boundary and interelement continuity equations. In addition, the existence of extremum of functional (228) can be easily proved by the “second variational approach” as well, which indicates functional (228) has a local extreme.


## 7.3.2. Element Stiffness Equation

The element stiffness equation can be generated by settingδΠme=0. To simplify the derivation, we first transform all domain integrals in (228) into boundary ones. With Gaussian theorem, the functional in (228) may be simplified as(240)Πme=12∬Ωeσijεij+DiEidΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=12∫ΓetiuidΓ-∬Ωeσij,juidΩ+12∫ΓeDnϕdΓ-∬ΩeDi,iϕdΩ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ.Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for the HT finite element model:(241)Πme=12∫ΓetiuidΓ+12∫ΓeDnϕdΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ+∫Γetiu~i-uidΓ+∫ΓeDnϕ~-ϕdΓ=-12∫Γetiui+DnϕdΓ+∫Γe(tiu~i+Dnϕ~)dΓ-∫Γtt-iu~idΓ-∫ΓDD-nϕ~dΓ.Substituting (214), (223), and (221) into the above functional (241) yields the formulation as(242)Πme=-12ceTHece+ceTGede-deTge,where (243)He=∫ΓeQeTNedΓ,Ge=∫ΓeQeTN~edΓ,ge=∫ΓtN~eTt-dΓ+∫ΓDN~eTD-dΓ.The stationary condition of the functional Πme with respect to ce and de, respectively, yields (33) and (34). Following the procedure described in [21, 22], the missing rigid-body motion can be recovered by setting the augmented internal field of a particular element e as(244)ue=Nece+10x2001-x100001c0,where the undetermined rigid-body motion parameter c0 can be calculated using the least square matching of ue and u~e at element nodes(245)min⁡=∑i=1nu1i-u~1i2+u2i-u~2i2+ϕi-ϕ~i2which finally gives(246)c0=Re-1re,(247)Re=∑i=1n10x2i001-x1i0x2i-x1ix1i2+x2i200001,(248)re=∑i=1nΔue1iΔue2iΔue1ix2i-Δue2ix1iΔϕeiin which Δuei=(u~e-u^e)nodei and n is the number of element nodes. As a consequence, c0 can be calculated by (246) once the nodal DEP fields de and the interpolation coefficients ce are, respectively, determined by (214) and (223). Then the complete DEP fields ue can be obtained from (244).


## 7.4. Normalization

The order of magnitudes of the material constants and the corresponding field variables in piezoelectricity have a wide spectrum as large as 1019 in SI unit. This will lead to ill-conditioned matrix of the system [72, 87]. To resolve this problem, normalization of each quantity by its reference value is employed. The reference values for the stiffness, piezoelectric stress constant, dielectric constants, and strain are selected to be c0=1011N/m2, e0=101(N/mV), k0=10-9C/mV, and ε0=10-3V/m, respectively. The reference values of other quantities, as shown in Table 1, are determined in terms of these four fundamental reference variables and the characteristic length x0=100m of the problem so that the normalized governing equations remain in exactly the same form as the original equations.Table 1
Reference values for material constants and field variables in piezoelectricity derived from basic reference variables:c0, e0, k0, ε0, and x0.
Displacement
u
0
=
x
0
ε
0
=
1
0
-
3
(
m
)
Electric Potential
ϕ
0
=
x
0
E
0
=
1
0
7
(
V
)
Stress
σ
0
=
c
0
ε
0
=
1
0
8
(
N
/
m
2
)
Electric induction
D
0
=
k
0
E
0
=
1
0
-
2
(
C
/
m
2
)
Compliance
s
0
=
ε
0
/
σ
0
=
1
0
-
11
(
m
2
/
N
)
Impermeability
β
0
=
E
0
/
D
0
=
1
0
9
(
mV
/
C
)
Electric field
E
0
=
σ
0
/
e
0
=
1
0
7
(
V
/
m
)
Piezoelectric strain constant
g
0
=
E
0
/
σ
0
=
1
0
-
1
(
mV
/
N
)


## 8. Numerical Examples

Several numerical examples are presented in this section to illustrate the application of the HFS-FEM and to demonstrate its effectiveness and accuracy. Unless otherwise indicated, mesh convergence tests were conducted for the reference solutions obtained from ABAQUS in the following examples.


### 8.1. Cubic Block under Uniform Tension and Body Force

An isotropic cubic block, with dimension 10 × 10 × 10 and subject to a uniform tension as shown in Figure8, is considered in this example. A constant body force of 10 Mpa and uniform distributed tension of 100 MPa are applied to the cube. Figure 9 presents the displacement component u1 and the stress component σ11 at point A of the block, which are calculated by the HFS-FEM on the three meshes with distorted 8-node brick elements (Figure 10). The results calculated by ABAQUS with a very fine mesh (with 40 × 40 × 40 C3D8 and EAS element with 68921 nodes) are taken as a reference value for comparison. The results from C3D8 and EAS elements in ABAQUS are also presented in Figure 9 for comparison. It can be seen from these figures that the results obtained from both the HFS-FEM and ABAQUS converge to the benchmark value with the number of degree of freedom (DOF) increasing. For Mesh 1, the hybrid EAS element has the best performance while for Mesh 2 and Mesh 3 it can be seen that HFS-FEM with 8-node brick elements exhibits better accuracy in both displacement and stress compared with EAS element in traditional FEM. Contour plots of u1 and σ11 obtained by HFS-FEM on Mesh 3 are also presented in Figure 11. It should be noted that for problems involving body forces the accuracy of the RBF interpolation has to be considered for a satisfactory solution. The details on the RBF interpolation can be found in previous literatures [43, 88–90].Figure 8
Cubic block under uniform tension and body force: geometry, boundary condition, and loading.Figure 9
Cubic block with body force under uniform distributed load: convergent study of (a) displacementu1 and (b) stress σ11.
(a)
(b)Figure 10
Cubic block under uniform tension and body force: (a) mesh 1 (4 × 4 × 4 elements), (b) mesh 2 (6 × 6 × 6 elements), and (c) mesh 3 (10 × 10 × 10 elements).
(a)
(b)
(c)Figure 11
Contour plots of (a) displacementu1 and (b) stress σ11 of the cube.
(a)
(b)


### 8.2. Circular Cylinder with Axisymmetric Temperature Change

In this example, a long circular cylinder with axisymmetric temperature change in domain is considered. Both inside and outside surfaces of the cylinder are assumed to be free from traction. The temperatureT changes logarithmically along the radial direction. With the symmetry condition of the problem, only one quarter of the cylinder is modeled. The configurations of geometry and the boundary conditions are shown in Figure 12. In our computation, the parameters a=5, b=20, E=1000, ν=0.3, α=0.001, and T0=10. The two approaches listed in Section 4 to approximate the body force and temperature are discussed and analyzed in this example.Figure 12
(a) Geometry and boundary conditions of the long cylinder with axisymmetric temperature change and (b) mesh configurations of a quarter of circular cylinder (128 eight-node elements).
(a)
(b)Figure13 presents the variation of the radial and the circumferential thermal stresses with the cylinder radius, respectively, in which the theoretical values are given for comparison [39]. It is seen from Figure 13 that the results from Method 2 are much better than those obtained from Method 1 for both radial and circumferential stress. It can be inferred that the error may be to a large extent due to the RBF interpolation, for which the number of interpolation points has a significant influence on its accuracy.Figure 13
(a) Radial thermal stresses and (b) circumferential thermal stresses with the cylinder radius.
(a)
(b)Figure14 displays the contour plots of (a) radial and (b) circumferential thermal stresses (the meshes used for contour plot are different from that for calculation due to using quadratic elements). It demonstrates that treating temperature gradient and body force together is more superior to dealing with them separately. However, we have to rely on Method 1 when the temperature change is discrete distribution or the gradient of the temperature field is not available.Figure 14
Contour plots of (a) radial and (b) circumferential thermal stresses (the background mesh in (a) and (b) is used for plots of the calculated results in postprocessing).
(a)
(b)


### 8.3. 3D Cube under Arbitrary Temperature and Body Force

As shown in Figure15, a 3D cube of 1×1×1 with center located at (0.5,0.5,0.5) is considered in this example. The material properties of the cube are Young’s modulus E=5000, Poisson’s ratio ν=0.3, and linear thermal expansion coefficient α=0.001. The bottom surface is fixed on the ground and the temperature distribution and body force are assumed to be(249)T=30x3,b1=0,b2=0,b3=-2000x-0.52+y-0.52.Because there is no analytical solution available, the results from ABAQUS herein are employed for comparison. The meshes used by HFS-FEM and ABAQUS are given in Figure 15.Figure 15
(a) Schematic of the 3D cube under arbitrary temperature and body force, (b) mesh used by HFS-FEM (125 20-node brick elements), and (c) mesh for ABAQUS (8000 C3D20R elements).
(a)
(b)
(c)Figure16 presents the displacement u3 and stress σ33 along one edge of the cube which is coinciding with x3 axis. It can be seen that the results from HFS-FEM again agree very well with those by ABAQUS. It is demonstrated that the proposed HFS-FEM is able to predict the response of 3D thermoelastic problems under arbitrary temperature and body force. It is also shown that HFS-FEM with RBF interpolation can give satisfactory results using coarse meshes.Figure 16
(a) Displacementu3 and (b) stress σ33 along one cube edge when subjected to arbitrary temperature and body force.
(a)
(b)


### 8.4. Orthotropic Composite Plate with an Elliptic Hole under Tension

A finite composite plate containing an elliptical hole (Figure17) is investigated in this example. A uniform tension of σ0=1GPa is applied in x2 direction. The material parameters of the orthotropic plate are taken as E1=113GPa, E2=52.7GPa, G12=28.5GPa, and v12=0.45. The length and width of the plate are L=20mm and W=20mm; the major and minor lengths of the ellipses a and b are, respectively, 2 mm and 1 mm. In the computation, plane stress condition is used. The mesh configurations used by HFS-FEM and ABAQUS are given in Figure 17.Figure 17
(a) Schematic of an orthotropic composite plate with an elliptic hole under uniform tension, and its mesh configurations for (b) HFS-FEM, 1515 quadratic elements, and (c) ABAQUS, 9498 quadratic elements.
(a)
(b)
(c)Figure18 shows the corresponding variation of the hoop stress along the rim of the elliptical hole when the orientation angle φ of reinforced fibers is equal to 0°, 45°, and 90°, respectively. It is found from Figure 18 that the results from HFS-FEM have a good agreement with the reference solutions from ABAQUS. This indicates that the proposed method is able to capture the variations of the hoop stress induced by the elliptical hole in the plate. The contour plots of stress components σ11, σ22, and σ12 around the elliptic hole in the composite plate for several different fiber angles are shown in Figure 19.Figure 18
Variation of hoop stresses along the rim of the elliptical hole for different fiber orientationφ.Figure 19
Contour plots of stress components around the elliptic hole in the composite plate: (a)φ=0°, (b) φ=45°, and (c) φ=90°.
(a)
(b)
(c)Figure20 shows the stress concentration factor (SCF) along with the inclined angle φ of the reinforced fibers, which exhibits a good agreement with the solutions from ABAQUS. It is obvious that the SCF of the punched plate rises with the increasing fiber angle φ. It is found from Figure 20 that the largest SCF occurs at φ=90°, whereas the smallest appears at φ=0°. It indicates the effectiveness of the proposed method in predicting the SCF for anisotropic composites as well.Figure 20
Variation of SCF with the lamina angleφ.


### 8.5. Isotropic Plate with Multianisotropic Inclusions

In this example, a multi-inclusion problem is investigated to show the capability of the HFS-FEM to deal with both isotropic and anisotropic materials in a unified way. As shown in Figure21, an isotropic plate containing multianisotropic inclusions of square geometry (edge length a=2) is considered. The distance between any two inclusions is assumed to be b=3. The material parameters for the inclusions are chosen as E1=134.45GPa, E2=E3=11.03GPa, G23=2.98GPa, G31=G12=28.5GPa, and v23=0.49, v31=v12=0.301. The material parameters for isotropic matrix are elastic modules E=2.8Gpa and poison’s ratio ν=0.3. The mesh configuration of the plate for HFS-FEM is shown in Figure 21, which uses 272 quadratic general elements.Figure 21
(a) Schematic and (b) mesh configuration of an isotropic plate with multianisotropic inclusions.
(a)
(b)In general, the Stroh formalism is suitable for the anisotropic material with distinct material eigenvalues, and it fails for the degenerated materials like isotropic material with repeated eigenvaluespα=i,(α=1,2,3) [59]. However, a small perturbation of the material constants, such as p1=(1-0.004)i, p2=i, and p3=(1+0.004)i, can be applied to make the eigenvalues distinct and the results can be applied conveniently.Table2 shows the displacement and stresses at points A, B, and C as indicated in Figure 21. It is observed that there is a good agreement between the results by the HFS-FEM and those from ABAQUS using very fine mesh, in which the maximum relative error for displacement and stress by HFS-FEM occur at Point B (i.e., x2=0) and are 0.7% and 1.3%, respectively. Additionally, it can be found that the results from HFS-FEM are better than those from ABAQUS using the same mesh. Although the stress at the vicinity of the inclusions (Point C) has a little degradation due to the high stress concentration and stress contrast at the adjacent elements, the displacement still agrees well with the reference solution. The variations of displacement components u1 and u2 along the right edge (x=8) by HFS-FEM are shown in Figure 22.Table 2
Comparison of displacement and stress at points A and B.
Items
Points
HFS-FEM (272 elements)
ABAQUS (272 elements)
ABAQUS (30471 elements )
Disp.u1
A
0.04322
0.04318
0.04335
B
0.03719
0.03721
0.03744
C
0.03062
0.03076
0.03091
Stressσ11
A
10.0446
9.9219
9.9992
B
9.8585
9.8304
9.9976
C
20.6453
13.7302
23.6625Figure 22
The variation of displacement component (a)u1 and (b) u2 along the right edge of the plate (x=8) by HFS-FEM and ABAQUS.
(a)
(b)


### 8.6. Infinite Piezoelectric Medium with Hole

Consider an infinite piezoelectric plane with a circular hole as shown in Figure23. The material parameters are given in Table 3. Suppose that mechanical load σxx∞=σ0=10 parallel to the x axis is imposed at infinity with traction and electric charge-free at the boundary of the hole. In our calculation, the radius of the hole is set to be r=1 and L/r=20 is employed in the analysis.Table 3
Properties of the material PZT-4 used for the model.
Parameters
Values
Parameters
Values
c
11
13.9
×
1
0
10
N
m
-
2
e
15
13.44
C
m
-
2
c
12
7.78
×
1
0
10
N
m
-
2
e
31
-
6.98
C
m
-
2
c
13
7.43
×
1
0
10
N
m
-
2
e
33
13.84
C
m
-
2
c
33
11.3
×
1
0
10
N
m
-
2
κ
11
6.0
×
1
0
-
9
C
/
Nm
c
44
2.56
×
1
0
10
N
m
-
2
κ
33
5.47
×
1
0
-
9
C
/
NmFigure 23
(a) An infinite piezoelectric plate with a circular hole subjected to remote stress. (b) Mesh used by HFS-FEM. (c) Mesh used by ABAQUS.
(a)
(b)
(c)Figure24 shows the distribution of hoop stress σθ and radial stress σr along the line z=0 for remote loading σxx∞ and along the line x=0 for remote loading σzz∞, respectively. Figure 25 presents the variations of the normalized stress σθ/σ0 and the normalized electric displacement Dθ/σ0×1010 along the hole edge under remote mechanical loading. It is obvious that the results obtained from HFS-FEM agree well with the results from ABAQUS and Sosa [74].Figure 24
Distribution of (a) hoop stressσθ and (b) radial stress σr along the line z=0 when subjected to remote mechanical load σxx∞ and along the line x=0 when subjected to remote mechanical load σzz∞.
(a)
(b)Figure 25
Variation of (a) normalized stressσθ/σ0 and (b) electrical displacement Dθ/σ0×1010 along the hole boundary under remote mechanical loading.
(a)
(b)It can be seen from Figure24 that hoop stress σθ has maximum value on the rim of the hole and it decreases dramatically with the increase of the distance from the hole edge. It is also shown that σθ tends to be equal to the remote applied load σ0 when r increases toward infinity. Compared with the hoop stress σθ, it is obvious that the radial stress σr is much smaller and usually does not need to be considered. It is obvious that loading along the poling direction will produce smaller stress concentration due to coupling effect. The maximum values of σθ appear at θ=90∘ for case of σxx∞ and at θ=0∘ and θ=180∘ for case of loading σzz∞, both of which agree well with the analytical solution from Sosa [74]. The electric displacement Dθ/σ0×1010 produced by σxx∞ and σzz∞ is nearly the same and is symmetrical with respect to the x axis. It is found that the maximum values of Dθ appear at θ=65∘ and θ=115∘, which also agrees well with the analytical solution.


## 8.1. Cubic Block under Uniform Tension and Body Force

An isotropic cubic block, with dimension 10 × 10 × 10 and subject to a uniform tension as shown in Figure8, is considered in this example. A constant body force of 10 Mpa and uniform distributed tension of 100 MPa are applied to the cube. Figure 9 presents the displacement component u1 and the stress component σ11 at point A of the block, which are calculated by the HFS-FEM on the three meshes with distorted 8-node brick elements (Figure 10). The results calculated by ABAQUS with a very fine mesh (with 40 × 40 × 40 C3D8 and EAS element with 68921 nodes) are taken as a reference value for comparison. The results from C3D8 and EAS elements in ABAQUS are also presented in Figure 9 for comparison. It can be seen from these figures that the results obtained from both the HFS-FEM and ABAQUS converge to the benchmark value with the number of degree of freedom (DOF) increasing. For Mesh 1, the hybrid EAS element has the best performance while for Mesh 2 and Mesh 3 it can be seen that HFS-FEM with 8-node brick elements exhibits better accuracy in both displacement and stress compared with EAS element in traditional FEM. Contour plots of u1 and σ11 obtained by HFS-FEM on Mesh 3 are also presented in Figure 11. It should be noted that for problems involving body forces the accuracy of the RBF interpolation has to be considered for a satisfactory solution. The details on the RBF interpolation can be found in previous literatures [43, 88–90].Figure 8
Cubic block under uniform tension and body force: geometry, boundary condition, and loading.Figure 9
Cubic block with body force under uniform distributed load: convergent study of (a) displacementu1 and (b) stress σ11.
(a)
(b)Figure 10
Cubic block under uniform tension and body force: (a) mesh 1 (4 × 4 × 4 elements), (b) mesh 2 (6 × 6 × 6 elements), and (c) mesh 3 (10 × 10 × 10 elements).
(a)
(b)
(c)Figure 11
Contour plots of (a) displacementu1 and (b) stress σ11 of the cube.
(a)
(b)


## 8.2. Circular Cylinder with Axisymmetric Temperature Change

In this example, a long circular cylinder with axisymmetric temperature change in domain is considered. Both inside and outside surfaces of the cylinder are assumed to be free from traction. The temperatureT changes logarithmically along the radial direction. With the symmetry condition of the problem, only one quarter of the cylinder is modeled. The configurations of geometry and the boundary conditions are shown in Figure 12. In our computation, the parameters a=5, b=20, E=1000, ν=0.3, α=0.001, and T0=10. The two approaches listed in Section 4 to approximate the body force and temperature are discussed and analyzed in this example.Figure 12
(a) Geometry and boundary conditions of the long cylinder with axisymmetric temperature change and (b) mesh configurations of a quarter of circular cylinder (128 eight-node elements).
(a)
(b)Figure13 presents the variation of the radial and the circumferential thermal stresses with the cylinder radius, respectively, in which the theoretical values are given for comparison [39]. It is seen from Figure 13 that the results from Method 2 are much better than those obtained from Method 1 for both radial and circumferential stress. It can be inferred that the error may be to a large extent due to the RBF interpolation, for which the number of interpolation points has a significant influence on its accuracy.Figure 13
(a) Radial thermal stresses and (b) circumferential thermal stresses with the cylinder radius.
(a)
(b)Figure14 displays the contour plots of (a) radial and (b) circumferential thermal stresses (the meshes used for contour plot are different from that for calculation due to using quadratic elements). It demonstrates that treating temperature gradient and body force together is more superior to dealing with them separately. However, we have to rely on Method 1 when the temperature change is discrete distribution or the gradient of the temperature field is not available.Figure 14
Contour plots of (a) radial and (b) circumferential thermal stresses (the background mesh in (a) and (b) is used for plots of the calculated results in postprocessing).
(a)
(b)


## 8.3. 3D Cube under Arbitrary Temperature and Body Force

As shown in Figure15, a 3D cube of 1×1×1 with center located at (0.5,0.5,0.5) is considered in this example. The material properties of the cube are Young’s modulus E=5000, Poisson’s ratio ν=0.3, and linear thermal expansion coefficient α=0.001. The bottom surface is fixed on the ground and the temperature distribution and body force are assumed to be(249)T=30x3,b1=0,b2=0,b3=-2000x-0.52+y-0.52.Because there is no analytical solution available, the results from ABAQUS herein are employed for comparison. The meshes used by HFS-FEM and ABAQUS are given in Figure 15.Figure 15
(a) Schematic of the 3D cube under arbitrary temperature and body force, (b) mesh used by HFS-FEM (125 20-node brick elements), and (c) mesh for ABAQUS (8000 C3D20R elements).
(a)
(b)
(c)Figure16 presents the displacement u3 and stress σ33 along one edge of the cube which is coinciding with x3 axis. It can be seen that the results from HFS-FEM again agree very well with those by ABAQUS. It is demonstrated that the proposed HFS-FEM is able to predict the response of 3D thermoelastic problems under arbitrary temperature and body force. It is also shown that HFS-FEM with RBF interpolation can give satisfactory results using coarse meshes.Figure 16
(a) Displacementu3 and (b) stress σ33 along one cube edge when subjected to arbitrary temperature and body force.
(a)
(b)


## 8.4. Orthotropic Composite Plate with an Elliptic Hole under Tension

A finite composite plate containing an elliptical hole (Figure17) is investigated in this example. A uniform tension of σ0=1GPa is applied in x2 direction. The material parameters of the orthotropic plate are taken as E1=113GPa, E2=52.7GPa, G12=28.5GPa, and v12=0.45. The length and width of the plate are L=20mm and W=20mm; the major and minor lengths of the ellipses a and b are, respectively, 2 mm and 1 mm. In the computation, plane stress condition is used. The mesh configurations used by HFS-FEM and ABAQUS are given in Figure 17.Figure 17
(a) Schematic of an orthotropic composite plate with an elliptic hole under uniform tension, and its mesh configurations for (b) HFS-FEM, 1515 quadratic elements, and (c) ABAQUS, 9498 quadratic elements.
(a)
(b)
(c)Figure18 shows the corresponding variation of the hoop stress along the rim of the elliptical hole when the orientation angle φ of reinforced fibers is equal to 0°, 45°, and 90°, respectively. It is found from Figure 18 that the results from HFS-FEM have a good agreement with the reference solutions from ABAQUS. This indicates that the proposed method is able to capture the variations of the hoop stress induced by the elliptical hole in the plate. The contour plots of stress components σ11, σ22, and σ12 around the elliptic hole in the composite plate for several different fiber angles are shown in Figure 19.Figure 18
Variation of hoop stresses along the rim of the elliptical hole for different fiber orientationφ.Figure 19
Contour plots of stress components around the elliptic hole in the composite plate: (a)φ=0°, (b) φ=45°, and (c) φ=90°.
(a)
(b)
(c)Figure20 shows the stress concentration factor (SCF) along with the inclined angle φ of the reinforced fibers, which exhibits a good agreement with the solutions from ABAQUS. It is obvious that the SCF of the punched plate rises with the increasing fiber angle φ. It is found from Figure 20 that the largest SCF occurs at φ=90°, whereas the smallest appears at φ=0°. It indicates the effectiveness of the proposed method in predicting the SCF for anisotropic composites as well.Figure 20
Variation of SCF with the lamina angleφ.


## 8.5. Isotropic Plate with Multianisotropic Inclusions

In this example, a multi-inclusion problem is investigated to show the capability of the HFS-FEM to deal with both isotropic and anisotropic materials in a unified way. As shown in Figure21, an isotropic plate containing multianisotropic inclusions of square geometry (edge length a=2) is considered. The distance between any two inclusions is assumed to be b=3. The material parameters for the inclusions are chosen as E1=134.45GPa, E2=E3=11.03GPa, G23=2.98GPa, G31=G12=28.5GPa, and v23=0.49, v31=v12=0.301. The material parameters for isotropic matrix are elastic modules E=2.8Gpa and poison’s ratio ν=0.3. The mesh configuration of the plate for HFS-FEM is shown in Figure 21, which uses 272 quadratic general elements.Figure 21
(a) Schematic and (b) mesh configuration of an isotropic plate with multianisotropic inclusions.
(a)
(b)In general, the Stroh formalism is suitable for the anisotropic material with distinct material eigenvalues, and it fails for the degenerated materials like isotropic material with repeated eigenvaluespα=i,(α=1,2,3) [59]. However, a small perturbation of the material constants, such as p1=(1-0.004)i, p2=i, and p3=(1+0.004)i, can be applied to make the eigenvalues distinct and the results can be applied conveniently.Table2 shows the displacement and stresses at points A, B, and C as indicated in Figure 21. It is observed that there is a good agreement between the results by the HFS-FEM and those from ABAQUS using very fine mesh, in which the maximum relative error for displacement and stress by HFS-FEM occur at Point B (i.e., x2=0) and are 0.7% and 1.3%, respectively. Additionally, it can be found that the results from HFS-FEM are better than those from ABAQUS using the same mesh. Although the stress at the vicinity of the inclusions (Point C) has a little degradation due to the high stress concentration and stress contrast at the adjacent elements, the displacement still agrees well with the reference solution. The variations of displacement components u1 and u2 along the right edge (x=8) by HFS-FEM are shown in Figure 22.Table 2
Comparison of displacement and stress at points A and B.
Items
Points
HFS-FEM (272 elements)
ABAQUS (272 elements)
ABAQUS (30471 elements )
Disp.u1
A
0.04322
0.04318
0.04335
B
0.03719
0.03721
0.03744
C
0.03062
0.03076
0.03091
Stressσ11
A
10.0446
9.9219
9.9992
B
9.8585
9.8304
9.9976
C
20.6453
13.7302
23.6625Figure 22
The variation of displacement component (a)u1 and (b) u2 along the right edge of the plate (x=8) by HFS-FEM and ABAQUS.
(a)
(b)


## 8.6. Infinite Piezoelectric Medium with Hole

Consider an infinite piezoelectric plane with a circular hole as shown in Figure23. The material parameters are given in Table 3. Suppose that mechanical load σxx∞=σ0=10 parallel to the x axis is imposed at infinity with traction and electric charge-free at the boundary of the hole. In our calculation, the radius of the hole is set to be r=1 and L/r=20 is employed in the analysis.Table 3
Properties of the material PZT-4 used for the model.
Parameters
Values
Parameters
Values
c
11
13.9
×
1
0
10
N
m
-
2
e
15
13.44
C
m
-
2
c
12
7.78
×
1
0
10
N
m
-
2
e
31
-
6.98
C
m
-
2
c
13
7.43
×
1
0
10
N
m
-
2
e
33
13.84
C
m
-
2
c
33
11.3
×
1
0
10
N
m
-
2
κ
11
6.0
×
1
0
-
9
C
/
Nm
c
44
2.56
×
1
0
10
N
m
-
2
κ
33
5.47
×
1
0
-
9
C
/
NmFigure 23
(a) An infinite piezoelectric plate with a circular hole subjected to remote stress. (b) Mesh used by HFS-FEM. (c) Mesh used by ABAQUS.
(a)
(b)
(c)Figure24 shows the distribution of hoop stress σθ and radial stress σr along the line z=0 for remote loading σxx∞ and along the line x=0 for remote loading σzz∞, respectively. Figure 25 presents the variations of the normalized stress σθ/σ0 and the normalized electric displacement Dθ/σ0×1010 along the hole edge under remote mechanical loading. It is obvious that the results obtained from HFS-FEM agree well with the results from ABAQUS and Sosa [74].Figure 24
Distribution of (a) hoop stressσθ and (b) radial stress σr along the line z=0 when subjected to remote mechanical load σxx∞ and along the line x=0 when subjected to remote mechanical load σzz∞.
(a)
(b)Figure 25
Variation of (a) normalized stressσθ/σ0 and (b) electrical displacement Dθ/σ0×1010 along the hole boundary under remote mechanical loading.
(a)
(b)It can be seen from Figure24 that hoop stress σθ has maximum value on the rim of the hole and it decreases dramatically with the increase of the distance from the hole edge. It is also shown that σθ tends to be equal to the remote applied load σ0 when r increases toward infinity. Compared with the hoop stress σθ, it is obvious that the radial stress σr is much smaller and usually does not need to be considered. It is obvious that loading along the poling direction will produce smaller stress concentration due to coupling effect. The maximum values of σθ appear at θ=90∘ for case of σxx∞ and at θ=0∘ and θ=180∘ for case of loading σzz∞, both of which agree well with the analytical solution from Sosa [74]. The electric displacement Dθ/σ0×1010 produced by σxx∞ and σzz∞ is nearly the same and is symmetrical with respect to the x axis. It is found that the maximum values of Dθ appear at θ=65∘ and θ=115∘, which also agrees well with the analytical solution.


## 9. Conclusions

In this paper, we have reviewed the HFS-FEM and its application in engineering applications. The HFS-FEM is a promising numerical method for solving complex engineering problems. The main advantages of this method include integration along the element boundaries only, easily adopting arbitrary polygonal or even curve-sided elements and symmetric and sparse stiffness matrix, and avoiding the singularity integral problem as encountered in BEM. Moreover, as in HT-FEM, this method offers the attractive possibility to develop accurate crack singular, corner, or perforated elements, simply by using appropriate special fundamental solutions as the trial functions of the intraelement displacements. It is noted that the HFS-FEM has attracted more attention of researchers in computational mechanics in the past few years, and good progress has been made in the field of potential problems, plane elasticity, piezoelectric problems, and so on. However, there are still many possible extensions and areas in need of further development in the future:(1)
to develop various special-purpose elements to effectively handle singularities attributable to local geometrical or load effects (holes, cracks, inclusions, interface, corner, and load singularities), with the special-purpose functions warranting that excellent results are obtained at minimal computational cost and without local mesh refinement,(2)
to extend the HFS-FEM to elastodynamics, fluid flow, thin and thick plate bending, and fracture mechanics,(3)
to develop efficient schemes for complex engineering structures and improve the related general purpose computer codes with good preprocessing and postprocessing capabilities,(4)
to extend this method to the case of multifield problems such as thermoelastic-piezoelectric materials and thermomagnetic-electric-mechanical materials, and to develop multiscale framework across from continuum to micro- and nanoscales for modeling heterogeneous materials.



---

*Source: 916029-2015-03-17.xml*