Advances in Mathematical Physics
Volume 2015 (2015), Article ID 916029, 38 pages
http://dx.doi.org/10.1155/2015/916029

Review Article

Hybrid Fundamental Solution Based Finite Element Method: Theory and Applications

Changyong Cao and Qing-Hua Qin

Research School of Engineering, The Australian National University, Acton, ACT 2601, Australia

Received 13 October 2014; Revised 23 December 2014; Accepted 24 December 2014

Academic Editor: Luigi C. Berselli

Copyright © 2015 Changyong Cao and Qing-Hua Qin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 -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 -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 Section 2, 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]

with the boundary conditions

where is the unknown field variable and represents the boundary flux, is the component of outward normal vector to the boundary , and and are specified functions on the related boundaries, respectively. The space derivatives are indicated by a subscript comma, that is, , and the subscript index takes values for two-dimensional and for three-dimensional problems. Additionally, the repeated subscript indices imply summation convention.

For convenience, (3) can be rewritten in matrix form as

with .

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 element occupying subdomain , 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):

where is undetermined coefficients, is the number of virtual sources, and is the fundamental solution to the partial differential equation:

as

Evidently, (5) analytically satisfies (1) due to the solution property of .

(a) 4-node 2D element

(b) 8-node 2D element

Figure 1: Intraelement field and frame field of a HFS-FEM element for 2D potential problems.

In 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 can be generated by means of the method employed in the MFS [32–35]

where is a dimensionless coefficient, is the point on the element boundary (the nodal point in this work), and is the geometrical centroid of the element (see Figure 1). Determination of was discussed in [31, 36] and is usually used in practice.

The corresponding outward normal derivative of on is

where

with

2.2.2. Auxiliary Frame Field

In order to enforce the conformity on the field variable , for instance, on of any two neighboring elements and , an auxiliary interelement frame field is used and expressed in terms of the same degrees of freedom (DOF) as those used in the conventional finite elements. In this case, is confined to the whole element boundary as

which is independently assumed along the element boundary in terms of nodal DOF , where 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

where 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 between intraelement fields of element “” and element “” (see Figure 3) [36, 37]:

Figure 3: Illustration of continuity between two adjacent elements “” and “.”

A modified variational functional is developed as follows:

where

in which the governing equation (1) is assumed to be satisfied, a priori, for deriving the HFS FE model. The boundary of a particular element consists of the following parts:

where represents the interelement boundary of the element “” 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:

from which the first-order variational yields

Using divergence theorem

we can obtain

For the displacement-based method, the potential conformity should be satisfied in advance

then (21) can be rewritten as

The Euler equation and boundary conditions can be obtained as

using the stationary condition .

As for the continuous requirement between two adjacent elements “” and “” given in (14), we can obtain it in the following way. When assembling elements “” and “,” we have

from which the vanishing variation of leads to the reciprocity condition on the interelement boundary .

If the following expression:

is uniformly positive (or negative) in the neighborhood of , where the displacement has such a value that and where stands for the stationary value of , we have

in which the relation that is identical on 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 and using the constrained conditions (26), we find

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 Figure 1), we can generate element stiffness equation by the variational functional derived in Section 2.3. Following the approach described in [21], the variational functional corresponding to a particular element of the present problem can be written as

Appling the divergence theorem (20) to the functional (29), we have the final functional for the HFS-FE model

Then, substituting (5), (9), and (12) into the functional (30) produces

in which

The symmetry of is obvious from the scalar definition (31) of variational functional .

To enforce interelement continuity on the common element boundary, the unknown vector should be expressed in terms of nodal DOF . The minimization of the functional with respect to and , respectively, yields

from which the optional relationship between and and the stiffness equation can be produced

where stands for the element stiffness matrix.

2.5. Numerical Integral for and 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 the matrix, one can express it as

by introducing the matrix function

Equation (36) can be further rewritten as

where

and is the Jacobean expressed as

where

Thus, the Gaussian numerical integration for matrix can be calculated by

where is the number of edges of the element and is the Gaussian sampling points employed in the Gaussian numerical integration. Similarly, we can calculate the matrix using

The calculation of vector 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 as

where the undetermined rigid-body motion parameter can be calculated using the least square matching of and at element nodes

which finally gives

in which and is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector can be evaluated from (34), and then is evaluated from (45). Finally, the potential field 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 (), the governing equations of a plane elastic body can be expressed as

If written as matrix form, it can be presented as

where is a stress vector, is a body force vector, and the differential operator matrix is given as

where is a strain vector and is a displacement vector.

The constitutive equations for the linear elasticity are given in matrix form as

where is the material coefficient matrix with constant components for isotropic materials, which can be expressed as follows:

where

The two different kinds of boundary conditions can be expressed as

where denotes the traction vector and is a transformation matrix related to the direction cosine of the outward normal

Substituting (49) and (50) into (47) yields the well-known Navier partial differential equations in terms of displacements:

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:

where 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 and the fundamental solution matrix are now in the form

in which and are, respectively, the field point and source point in the local coordinate system (). The components are the fundamental solution, that is, induced displacement component in -direction at the field point due to a unit point load applied in -direction at the source point , which are given by [40, 41]

where . 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):

where

As a consequence, the traction is written as

in which

The components for plane strain problems are given as

The unknown 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 with nodal displacement . Thus, the frame is defined as

where the symbol “” is used to specify that the field is defined on the element boundary only, is the matrix of shape functions, and 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, and can be expressed as

where , and can be expressed by natural coordinate as

3.3. Modified Functional for the Hybrid FEM

As in Section 2, 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 used for deriving the present HFS-FEM can be constructed as [22]

where and are the intraelement displacement field defined within the element and the frame displacement field defined on the element boundary, respectively. and are the element domain and element boundary, respectively.  , and stand, respectively, for the specified traction boundary, specified displacement boundary, and interelement boundary (). 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

Due to the satisfaction of the equilibrium equation with the constructed intraelement field, we have the following expression for HFS finite element model:

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 Section 2, the element stiffness equation can be generated by setting . Substituting (56), (64), and (61) into the functional of (69), we have

where

To enforce interelement continuity on the common element boundary, the unknown vector should be expressed in terms of nodal DOF . The stationary condition of the functional with respect to and yields, respectively, (33) and (34).

3.5. Recovery of Rigid-Body Motion

For the same reason stated in Section 2.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]

where the undetermined rigid-body motion parameter can be calculated using the least square matching of and at element nodes

which finally gives

where

in which and is the number of element nodes. Once the nodal field is determined by solving the final stiffness equation, the coefficient vector can be evaluated from (56), and then is evaluated from (74). Finally, the displacement field 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 () 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

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

where is the stress tensor, is the strain tensor, and is the Kronecker delta. Substituting (77) into (76), the equilibrium equation is rewritten as

For a well-posed boundary value problem, the boundary conditions are prescribed as follows:

where is the boundary of the solution domain and 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 term 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 is decomposed into two parts, a homogeneous solution and a particular solution

where the particular solution should satisfy the governing equation:

without any restriction of boundary condition. However, the homogeneous solution should satisfy

with the modified boundary conditions

From the above equations it can be seen that once the particular solution is known, the homogeneous solution 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 force , 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

where is the number of interpolation points, are the RBFs, and are the coefficients to be determined. Subsequently, the particular solution can be approximated by

where 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:

To solve this equation, the displacement is expressed in terms of the Galerkin-Papkovich vectors

Substituting (87) into (86), we can obtain the following biharmonic equation:

Taking the Spline Type RBF as an example, we get the following solutions:

where

and represents the Euclidean distance between a field point and a given point in the domain of interest. The corresponding particular solution of stresses can be obtained as

where . Substituting (90) into (92), we have

where

4.4. HFS-FEM for Homogeneous Solution

After obtaining the particular solution in Sections 4.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, , 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

where the matrix and unknown vector can be further written as

in which and are, respectively, the field point and source point in the local coordinate system (). The fundamental solution is given by [40]

where , is the number of source points. The source point can also be generated by means of the following method [36] as in two-dimensional cases

where is a dimensionless coefficient, is the point on the element boundary (the nodal point in this work), and is the geometrical centroid of the element (see Figure 5).

(a)

(b)

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).

According to (50) and (49), the corresponding stress fields can be expressed as

where

The components are given by

As a consequence, the traction can be written in the form

in which

To link the unknown and the nodal displacement , the frame is defined as

where the symbol “” is used to specify that the field is defined on the element boundary only, is the matrix of shape functions, and 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 and vector can be expressed as

where the shape functions are expressed as

where () can be expressed by natural coordinate

and is the natural coordinate of the -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 used for deriving the present HFS-FEM is written as [22]

By applying the Gaussian theorem, (108) can be simplified as

Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for HFS finite element model:

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

where

To enforce interelement continuity on the common element boundary, the unknown vector should be expressed in terms of nodal DOF . The stationary condition of the functional with respect to and yields again, respectively, (33) and (34).

4.4.4. Numerical Integral over Element

Considering a surface of the 3D hexahedron element, as shown in Figure 6, the vector normal to the surface can be obtained by

where and are the tangential vectors in the -direction and -direction, respectively, calculated by

where is the number of nodes of the surface and are the nodal coordinates. Thus the unit normal vector is given by

where

is the Jacobian of the transformation from Cartesian coordinates () to natural coordinates ().

For the matrix, we introduce the matrix function

Then we can get

and we rewrite it to the component form as

Using the relationship

and the Gaussian numerical integration, we can obtain

where and 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 matrix by

It should be mentioned that the calculation of vector 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 Section 2, 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

and using a least-square procedure to match and at the nodes of the element boundary

The above equation finally yields

where

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 () denote the coordinates in Cartesian coordinate system. The equilibrium governing equations of the thermoelasticity with the body force are expressed as

where is the stress tensor, is the body force vector, and . The generalized thermoelastic stress-strain relations and kinematical relation are given as

where is the strain tensor, is the displacement vector, is the temperature change, is the shear modulus, is the Poisson’s ratio, is the Kronecker delta, and

is the thermal constant with being the coefficient of linear thermal expansion. Substituting (128) into (127), the equilibrium equations may be rewritten as

For a well-posed boundary value problem, the following boundary conditions, either displacement or traction boundary condition, should be prescribed as

where is the boundary of the solution domain , and are the prescribed boundary values, and

is the boundary traction, in which 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 can be eliminated by employing the method of particular solution [16, 36, 44]. Using superposition principle, we decompose the displacement into two parts, the homogeneous solution and the particular solution as follows:

in which the particular solution should satisfy the governing equation

but 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

with modified boundary conditions

From above equations, it can be seen that once the particular solution is known, the homogeneous solution 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 force and the temperature field in order to obtain the particular solution. To implement this approximation, we may consider two different ways: one is to treat body force and the temperature field separately as in Tsai [54]. The other is to treat 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 force and temperature are assumed to be by the following two equations:

where is the number of interpolation points, are the basis functions, and and are the coefficients to be determined by collocation. Subsequently, the approximate particular solution can be written as follows:

where and 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:

To solve (140), the displacement is expressed in terms of the Galerkin-Papkovich vectors [43, 55–57]

Substituting (142) into (140), one can obtain the following biharmonic equation:

If taking the Spline Type RBF , one can get the following solutions:

where

for two-dimensional problem and

where

for three-dimensional problem, where represents the Euclidean distance of the given point from a fixed point in the domain of interest. The corresponding stress particular solution can be obtained by

where . Substituting (147) into (149), one can obtain

where

for two-dimensional problem, and substituting (147) into (149), one can obtain

where

for three-dimensional problem.

To solve (141), one can treat as the gradient of a scalar function

Substituting (154) into (141) obtains the Poisson’s equation

Thus, taking , its particular solution can be obtained [55] as follows:

Then from (154) we can get as follows:

The corresponding stress particular solution can be obtained by substituting (147) into

Then we have

5.3.2. Interpolating Temperature and Body Force Together

Considering the temperature gradient plays the role of body force, we can approximate together by the following equation:

Thus the approximate particular solution can be written as

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 and , 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 (), if we neglect the body force , the equilibrium equations, stress-strain laws, and strain-displacement equations for anisotropic elasticity are [61]

where , 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

The boundary conditions of the boundary value problem (163)–(165) are

where and are the prescribed boundary displacement and traction vector, respectively. In addition, is the unit outward normal to the boundary and is the boundary of the solution domain .

For the generalized two-dimensional deformation of anisotropic elasticity is assumed to depend on and only. Based on this assumption, the general solution to (165) can be written as [61, 62]

where is the displacement vector, is the stress function vector, and is a function vector composed of three holomorphic complex functions , which is an arbitrary function with argument 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, are the material eigenvalues with positive imaginary part, and and are 3 × 3 complex matrices formed by the material eigenvector associated with , which can be obtained by the following Eigen relations [61]:

where is a 6 × 6 foundational elasticity matrix and is a 6 × 1 column vector defined by

where , and the matrices and are 3 × 3 matrices extracted from as follows:

The stresses can be obtained from the derivation of stress functions as follows:

where

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) applied at an internal point far from the boundary. The boundary conditions of this problem can be written as

Thus, the Green’s function satisfying the above boundary conditions is found to be [72]

Therefore, fundamental solutions of the problem can be expressed as

The corresponding stress components can be obtained from stress function as

where are chosen to be , and , respectively, stands for the diagonal matrix corresponding to subscript , Im denotes the imagery part of a complex number, and superscript 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 frame , and of the structural element rather than to coordinates , and 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 Figure 7, the angle between the axis- and axis- 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

where the transformation matrix and its inverse matrix are defined as

with . Therefore, the constitutive relationship in the global coordinate system is given by

Figure 7: Schematic of the relationship between global coordinate system () and local material coordinate system ().

6.2. Formulations of HFS-FEM

6.2.1. Assumed Independent Fields

The intraelement displacement fields for a particular element is approximated in terms of a linear combination of fundamental solutions of the problem as

where the matrix and unknown vector can be further written as

in which is the number of source points, and are, respectively, the field point and source point in the coordinate system () local to the element under consideration. The fundamental solution 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

where

The components are given by (176) when is selected to be , , and , respectively. As a consequence, the traction can be written as

in which

The unknown 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

where the symbol “” is used to specify that the field is defined on the element boundary only, is the matrix of shape functions, and 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, and can be expressed as

where the shape functions are expressed as

and , and are expressed by natural coordinate

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 for a particular element is constructed as

where the boundary of the element is

and is the interelement boundary of element . Performing a variation of , one obtains

Applying Gaussian theorem

and the definitions of traction force

we obtain

Then, substituting (195) into (192) gives

Considering the fact that

we can finally obtain the following form:

Therefore, the Euler equations for (190) result in (165) and (166) because the quantities , and 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 and

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:

Substituting (180), (184), and (186) into the functional (200) yields the formulation as

where

The stationary condition of the functional with respect to and , 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 system () are given by

where is the stress tensor, 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

where is the elasticity tensor measured under zero electric field, and are, respectively, the piezoelectric tensor and dielectric tensor measured under zero strain. and are the elastic strain tensor and the electric field vector, respectively. The relation between the strain tensor and the displacement is given by

and the electric field component is related to the electric potential by

The boundary conditions of the boundary value problem (203)–(206) can be defined by

where , and are, respectively, the prescribed boundary displacement, the prescribed traction vector, the prescribed surface charge, and the prescribed electric potential. In addition, is the unit outward normal to the boundary and

is the boundary of the solution domain .

For the transversely isotropic material, if is taken as the isotropic plane, one can employ either or plane to study the plane electromechanical phenomenon. Thus, choosing the former and considering the plane strain conditions ( and ) (204) can be reduced to

For the plane stress piezoelectric problem ( and ), the constitutive equations can be obtained by replacing the coefficients , and in (212) as

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 field and an independent DEP frame field along element boundaries [21, 29].

7.2.1. Intraelement Field

The intraelement DEP field 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

where the fundamental solution matrix is now given by

in which and are, respectively, the field point (i.e., the nodal points of the element in this work) and source point in the local coordinate system (). The component is the induced displacement component () or electric potential () in -direction at the field point due to a unit point load () or point charge () applied in -direction at the source point . The fundamental solution is given as [76, 82, 84]

where and is the three different roots of the characteristic equation . The source point can be generated by the following method [36]:

Making use of (205) and the expression of intraelement DEP field in (214), the corresponding stress and electric displacement in (212) can be written as

where and

in which denotes the corresponding stress components () or electric displacement () along -direction at the field point due to a unit point load () or a unit point charge () applied in -direction at the source point and can be derived from (216) and are listed below, which are derived by substituting the fundamental solutions into constitutive equations [85]:

in which the coefficients , and are defined as in literature [84].

From (204), (208), and (209), the generalized traction forces and electric displacement are given as

where

7.2.2. Auxiliary Frame Field

For the two-dimensional piezoelectric problem under consideration, the frame field is assumed as

where 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 and nodal vector can be given in the form

where the shape functions are expressed by natural coordinate

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]:

where “” and “” 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 for a particular element is constructed as

where

and the boundary of the element is

and is the interelement boundary of element . 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 , one obtains [86]

in which the first term is

Applying Gaussian theorem

and the definitions of traction force and electrical displacement

we obtain

Then, substituting (235) into (231) gives

Considering the fact that

we finally get the following form:

Therefore, the Euler equations for (238) result in (203)–(210) and (226) because the quantities , 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 and :

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 . 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

Due to the satisfaction of the equilibrium equation with the constructed intraelement fields, we have the following expression for the HT finite element model:

Substituting (214), (223), and (221) into the above functional (241) yields the formulation as

where

The stationary condition of the functional with respect to and , 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 as

where the undetermined rigid-body motion parameter can be calculated using the least square matching of and at element nodes

which finally gives

in which and is the number of element nodes. As a consequence, can be calculated by (246) once the nodal DEP fields and the interpolation coefficients are, respectively, determined by (214) and (223). Then the complete DEP fields 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 10¹⁹ 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 , , , and , 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 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: , , , , and .

Displacement

Electric Potential

Stress

Electric induction

Compliance

Impermeability

Electric field

Piezoelectric strain constant

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 Figure 8, 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 and the stress component at point 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 and 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.

(a)

(b)

Figure 9: Cubic block with body force under uniform distributed load: convergent study of (a) displacement and (b) stress .

(a)

(b)

(c)

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)

Figure 11: Contour plots of (a) displacement and (b) stress of the cube.

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 temperature 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 , , , , , and . The two approaches listed in Section 4 to approximate the body force and temperature are discussed and analyzed in this example.

(a)

(b)

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).

Figure 13 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.

(a)

(b)

Figure 13: (a) Radial thermal stresses and (b) circumferential thermal stresses with the cylinder radius.

Figure 14 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.

(a)

(b)

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).

8.3. 3D Cube under Arbitrary Temperature and Body Force

As shown in Figure 15, a 3D cube of with center located at () is considered in this example. The material properties of the cube are Young’s modulus , Poisson’s ratio , and linear thermal expansion coefficient . The bottom surface is fixed on the ground and the temperature distribution and body force are assumed to be

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.

(a)

(b)

(c)

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).

Figure 16 presents the displacement and stress along one edge of the cube which is coinciding with 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.

(a)

(b)

Figure 16: (a) Displacement and (b) stress along one cube edge when subjected to arbitrary temperature and body force.

8.4. Orthotropic Composite Plate with an Elliptic Hole under Tension

A finite composite plate containing an elliptical hole (Figure 17) is investigated in this example. A uniform tension of is applied in direction. The material parameters of the orthotropic plate are taken as , , , and . The length and width of the plate are and ; the major and minor lengths of the ellipses and 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.

(a)

(b)

(c)

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.

Figure 18 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 , , and 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 .

(a)

(b)

(c)

Figure 19: Contour plots of stress components around the elliptic hole in the composite plate: (a) , (b) , and (c) .

Figure 20 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 , whereas the smallest appears at . 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 Figure 21, an isotropic plate containing multianisotropic inclusions of square geometry (edge length ) is considered. The distance between any two inclusions is assumed to be . The material parameters for the inclusions are chosen as , , , , and , . The material parameters for isotropic matrix are elastic modules and poison’s ratio . The mesh configuration of the plate for HFS-FEM is shown in Figure 21, which uses 272 quadratic general elements.

(a)

(b)

Figure 21: (a) Schematic and (b) mesh configuration of an isotropic plate with multianisotropic inclusions.

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 eigenvalues [59]. However, a small perturbation of the material constants, such as , , and , can be applied to make the eigenvalues distinct and the results can be applied conveniently.

Table 2 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., ) 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 and along the right edge () 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.
A
0.04322
0.04318
0.04335

B
0.03719
0.03721
0.03744

C
0.03062
0.03076
0.03091

Stress
A
10.0446
9.9219
9.9992

B
9.8585
9.8304
9.9976

C
20.6453
13.7302
23.6625

(a)

(b)

Figure 22: The variation of displacement component (a) and (b) along the right edge of the plate () by HFS-FEM and ABAQUS.

8.6. Infinite Piezoelectric Medium with Hole

Consider an infinite piezoelectric plane with a circular hole as shown in Figure 23. The material parameters are given in Table 3. Suppose that mechanical load parallel to the 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 and is employed in the analysis.

Table 3: Properties of the material PZT-4 used for the model.

Parameters
Values
Parameters
Values

(a)

(b)

(c)

Figure 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.

Figure 24 shows the distribution of hoop stress and radial stress along the line for remote loading and along the line for remote loading , respectively. Figure 25 presents the variations of the normalized stress and the normalized electric displacement 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].

(a)

(b)

Figure 24: Distribution of (a) hoop stress and (b) radial stress along the line when subjected to remote mechanical load and along the line when subjected to remote mechanical load .

(a)

(b)

Figure 25: Variation of (a) normalized stress and (b) electrical displacement along the hole boundary under remote mechanical loading.

It can be seen from Figure 24 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 when increases toward infinity. Compared with the hoop stress , it is obvious that the radial stress 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 for case of and at and for case of loading , both of which agree well with the analytical solution from Sosa [74]. The electric displacement produced by and is nearly the same and is symmetrical with respect to the axis. It is found that the maximum values of appear at and , 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

H. Wang and Q. H. Qin, “FE approach with Green's function as internal trial function for simulating bioheat transfer in the human eye,” Archives of Mechanics, vol. 62, no. 6, pp. 493–510, 2010.

C. Cao, A. Yu, and Q. H. Qin, “Evaluation of effective thermal conductivity of fiber-reinforced composites,” International Journal of Architecture, Engineering and Construction, vol. 1, no. 1, pp. 14–29, 2012.

W. J. Santos, J. A. F. Santiago, and J. C. F. Telles, “Optimal positioning of anodes and virtual sources in the design of cathodic protection systems using the method of fundamental solutions,” Engineering Analysis with Boundary Elements, vol. 46, pp. 67–74, 2014.

A. Karageorghis, D. Lesnic, and L. Marin, “The method of fundamental solutions for an inverse boundary value problem in static thermo-elasticity,” Computers and Structures, vol. 135, pp. 32–39, 2014.

C. Y. Cao, Q. H. Qin, and A. B. Yu, “Micromechanical analysis of heterogeneous composites using hybrid trefftz FEM and hybrid fundamental solution based FEM,” Journal of Mechanics, vol. 29, no. 4, pp. 661–674, 2013.

K. Y. Wang, Q. H. Qin, Y. L. Kang, J. S. Wang, and C. Y. Qu, “A direct constraint-Trefftz FEM for analysing elastic contact problems,” International Journal for Numerical Methods in Engineering, vol. 63, no. 12, pp. 1694–1718, 2005.

F. L. S. Bussamra, E. Lucena Neto, and W. M. Ponciano, “Simulation of stress concentration problems by Hexahedral Hybrid-Trefftz finite element models,” Computer Modeling in Engineering & Sciences, vol. 99, no. 3, pp. 255–272, 2014.

C. J. Pearce, G. Edwards, and L. Kaczmarczyk, “3D cohesive crack propagation using hybrid-Trefftz finite elements,” in Computational Modelling of Concrete Structures, N. Bicanic, H. Mang, G. Meschke, and R. de Borst, Eds., vol. 1, CRC Press, New York, NY, USA, 2014.

H. Wang and Q. H. Qin, “Special fiber elements for thermal analysis of fiber-reinforced composites,” Engineering Computations, vol. 28, no. 8, pp. 1079–1097, 2011.

Q. H. Qin, “Nonlinear analysis of reissner plates on an elastic foundation by the BEM,” International Journal of Solids and Structures, vol. 30, no. 22, pp. 3101–3111, 1993.

O. Atak, S. Jonckheere, E. Deckers, D. Huybrechs, B. Pluymers, and W. Desmet, “A hybrid Boundary Element-Wave BASed Method for an efficient solution of bounded acoustic problems with inclusions,” Computer Methods in Applied Mechanics and Engineering, vol. 283, pp. 1260–1277, 2015.

S. Natarajan, J. Wang, C. Song, and C. Birk, “Isogeometric analysis enhanced by the scaled boundary finite element method,” Computer Methods in Applied Mechanics and Engineering, vol. 283, pp. 733–762, 2015.

H. T. Wang, Z. H. Yao, and P. B. Wang, “On the preconditioners for fast multipole boundary element methods for 2D multi-domain elastostatics,” Engineering Analysis with Boundary Elements, vol. 29, no. 7, pp. 673–688, 2005.

Q. H. Qin and S. Diao, “Nonlinear analysis of thick plates on an elastic foundation by HT FE with p-extension capabilities,” International Journal of Solids and Structures, vol. 33, no. 30, pp. 4583–4604, 1996.

J. Jirousek, A. Wróblewski, Q. H. Qin, and X. Q. He, “A family of quadrilateral hybrid-Trefftz p-elements for thick plate analysis,” Computer Methods in Applied Mechanics and Engineering, vol. 127, no. 1–4, pp. 315–344, 1995.

H. Wang, Q. H. Qin, and D. Arounsavat, “Application of hybrid Trefftz finite element method to non-linear problems of minimal surface,” International Journal for Numerical Methods in Engineering, vol. 69, no. 6, pp. 1262–1277, 2007.

S.-W. Gao, Y.-S. Wang, Z.-M. Zhang, and X.-R. Ma, “Dual reciprocity boundary element method for flexural waves in thin plate with cutout,” Applied Mathematics and Mechanics, vol. 26, no. 12, pp. 1564–1573, 2005.

M. Dhanasekar, J. J. Han, and Q. H. Qin, “A hybrid-Trefftz element containing an elliptic hole,” Finite Elements in Analysis and Design, vol. 42, no. 14-15, pp. 1314–1323, 2006.

J. A. T. de Freitas and Z.-Y. Ji, “Hybrid-trefftz equilibrium model for crack problems,” International Journal for Numerical Methods in Engineering, vol. 39, no. 4, pp. 569–584, 1996.

C. Hennuyer, N. Leconte, B. Langrand, and E. Markiewicz, “Interpolation functions of a hybrid-Trefftz perforated super-element featuring nodes on the hole boundary,” Finite Elements in Analysis and Design, vol. 91, pp. 40–47, 2014.

Q. H. Qin, The Trefftz Finite and Boundary Element Method, WIT Press, Southampton, UK, 2000.

Q. H. Qin and H. Wang, MATLAB and C Programming for Trefftz Finite Element Methods, CRC Press, New York, NY, USA, 2008.

J. A. de Freitas and M. Toma, “Hybrid-Trefftz stress elements for incompressible biphasic media,” International Journal for Numerical Methods in Engineering, vol. 79, no. 2, pp. 205–238, 2009.

J. Jirousek and Q. H. Qin, “Application of hybrid-Trefftz element approach to transient heat conduction analysis,” Computers & Structures, vol. 58, no. 1, pp. 195–201, 1996.

N. Leconte, B. Langrand, and E. Markiewicz, “On some features of a plate hybrid-Trefftz displacement element containing a hole,” Finite Elements in Analysis and Design, vol. 46, no. 10, pp. 819–828, 2010.

I. D. Moldovan, T. D. Cao, and J. A. T. de Freitas, “Hybrid-Trefftz finite elements for biphasic elastostatics,” Finite Elements in Analysis and Design, vol. 66, pp. 68–82, 2013.

I. D. Moldovan, T. D. Cao, and J. A. de Freitas, “Hybrid-Trefftz displacement finite elements for elastic unsaturated soils,” International Journal of Computational Methods, vol. 11, no. 2, Article ID 1342005, 31 pages, 2014.

Q. H. Qin, “Hybrid Trefftz finite-element approach for plate bending on an elastic foundation,” Applied Mathematical Modelling, vol. 18, no. 6, pp. 334–339, 1994.

Q. H. Qin, “Variational formulations for TFEM of piezoelectricity,” International Journal of Solids and Structures, vol. 40, no. 23, pp. 6335–6346, 2003.

Q. H. Qin, “Trefftz finite element method and its applications,” Applied Mechanics Reviews, vol. 58, no. 5, pp. 316–337, 2005.

H. Wang and Q. H. Qin, “Hybrid FEM with fundamental solutions as trial functions for heat conduction simulation,” Acta Mechanica Solida Sinica, vol. 22, no. 5, pp. 487–498, 2009.

H. Wang, Q. H. Qin, and Y.-L. Kang, “A meshless model for transient heat conduction in functionally graded materials,” Computational Mechanics, vol. 38, no. 1, pp. 51–60, 2006.

W. Hui and Q. H. Qin, “Some problems with the method of fundamental solution using radial basis functions,” Acta Mechanica Solida Sinica, vol. 20, no. 1, pp. 21–29, 2007.

H. Wang, Q. H. Qin, and Y. L. Kang, “A new meshless method for steady-state heat conduction problems in anisotropic and inhomogeneous media,” Archive of Applied Mechanics, vol. 74, no. 8, pp. 563–579, 2005.

M. A. Golberg, C. S. Chen, and H. Bowman, “Some recent results and proposals for the use of radial basis functions in the BEM,” Engineering Analysis with Boundary Elements, vol. 23, no. 4, pp. 285–296, 1999.

H. Wang and Q. H. Qin, “Fundamental-solution-based finite element model for plane orthotropic elastic bodies,” European Journal of Mechanics, A/Solids, vol. 29, no. 5, pp. 801–809, 2010.

H. Wang, Q. H. Qin, and X.-P. Liang, “Solving the nonlinear Poisson-type problems with F-Trefftz hybrid finite element model,” Engineering Analysis with Boundary Elements, vol. 36, no. 1, pp. 39–46, 2012.

H. C. Simpson and S. J. Spector, “On the positivity of the second variation in finite elasticity,” Archive for Rational Mechanics and Analysis, vol. 98, no. 1, pp. 1–30, 1987.

S. P. Timoshenko and J. N. Goodier, Theory of Elasticity, McGraw-Hill, New York, NY, USA, 2nd edition, 1951.

S. A. Sauter and C. Schwab, Boundary Element Methods, Springer, Berlin, Germany, 2010.

Q. H. Qin and Y. Y. Huang, “BEM of postbuckling analysis of thin plates,” Applied Mathematical Modelling, vol. 14, no. 10, pp. 544–548, 1990.

C. Cao, Q. H. Qin, and A. Yu, “A new hybrid finite element approach for three-dimensional elastic problems,” Archive of Mechanics, vol. 64, no. 3, pp. 261–292, 2012.

A. H.-D. Cheng, C. S. Chen, M. A. Golberg, and Y. F. Rashed, “BEM for theomoelasticity and elasticity with body force—a revisit,” Engineering Analysis with Boundary Elements, vol. 25, no. 4-5, pp. 377–387, 2001.

D. P. Henry and P. K. Banerjee, “A new boundary element formulation for two- and three-dimensional thermoelasticity using particular integrals,” International Journal for Numerical Methods in Engineering, vol. 26, no. 9, pp. 2061–2077, 1988.

D. Aubry, D. Lucas, and B. Tie, “Adaptive strategy for transient/coupled problems. Applications to thermoelasticity and elastodynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 176, no. 1–4, pp. 41–50, 1999.

R. Carrazedo and H. B. Coda, “Alternative positional FEM applied to thermomechanical impact of truss structures,” Finite Elements in Analysis and Design, vol. 46, no. 11, pp. 1008–1016, 2010.

R. J. Yang, “Shape design sensitivity analysis of thermoelasticity problems,” Computer Methods in Applied Mechanics and Engineering, vol. 102, no. 1, pp. 41–60, 1993.

A. Chaudouet, “Three-dimensional transient thermo-elastic analyses by the BIE method,” International Journal for Numerical Methods in Engineering, vol. 24, no. 1, pp. 25–45, 1987.

L. L. Cao, Q. H. Qin, and N. Zhao, “An RBF-MFS model for analysing thermal behaviour of skin tissues,” International Journal of Heat and Mass Transfer, vol. 53, no. 7-8, pp. 1298–1307, 2010.

Z.-W. Zhang, H. Wang, and Q. H. Qin, “Transient bioheat simulation of the laser-tissue interaction in human skin using hybrid finite element formulation,” Molecular & Cellular Biomechanics, vol. 9, no. 1, pp. 31–54, 2012.

Q. H. Qin and J. Q. Ye, “Thermoelectroelastic solutions for internal bone remodeling under axial and transverse loads,” International Journal of Solids and Structures, vol. 41, no. 9-10, pp. 2447–2460, 2004.

Q. H. Qin, “Thermoelectroelastic analysis of cracks in piezoelectric half-plane by BEM,” Computational Mechanics, vol. 23, no. 4, pp. 353–360, 1999.

C. Cao, Q. H. Qin, and A. Yu, “A novel boundary-integral based finite element method for 2D and 3D thermo-elasticity problems,” Journal of Thermal Stresses, vol. 35, no. 10, pp. 849–876, 2012.

C. C. Tsai, “The method of fundamental solutions with dual reciprocity for three-dimensional thermoelasticity under arbitrary body forces,” Engineering Computations (Swansea, Wales), vol. 26, no. 3, pp. 229–244, 2009.

A. H.-D. Cheng, “Particular solutions of Laplacian, Helmholtz-type, and polyharmonic operators involving higher order radial basis functions,” Engineering Analysis with Boundary Elements, vol. 24, no. 7-8, pp. 531–538, 2000.

X.-R. Fu, M.-W. Yuan, S. Cen, and G. Tian, “Characteristic equation solution strategy for deriving fundamental analytical solutions of 3D isotropic elasticity,” Applied Mathematics and Mechanics, vol. 33, no. 10, pp. 1253–1264, 2012.

D. Palaniappan, “A general solution of equations of equilibrium in linear elasticity,” Applied Mathematical Modelling, vol. 35, no. 12, pp. 5494–5499, 2011.

V. V. Vasiliev and E. V. Morozov, Advanced Mechanics of Composite Materials, Elsevier Science, 2007.

S. G. Lekhnitskii, Theory of Elasticity of An Anisotropic Body, Mir Publishers, Moscow, Russia, 1981.

Q. H. Qin, “New solution for thermopiezoelectric solid with an insulated elliptic hole,” Acta Mechanica Sinica, vol. 14, no. 2, pp. 157–170, 1998.

T. C. T. Ting, Anisotropic Elasticity: Theory and Applications, vol. 45 of Oxford Science Publications, Oxford University Press, New York, NY, USA, 1996.

A. N. Stroh, “Dislocations and cracks in anisotropic elasticity,” Philosophical Magazine, vol. 3, no. 30, pp. 625–646, 1958.

Q. H. Qin, Y.-W. Mai, and S.-W. Yu, “Some problems in plane thermopiezoelectric materials with holes,” International Journal of Solids and Structures, vol. 36, no. 3, pp. 427–439, 1999.

Q. H. Qin and Y. W. Mai, “Crack growth prediction of an inclined crack in a half-plane thermopiezoelectric solid,” Theoretical and Applied Fracture Mechanics, vol. 26, no. 3, pp. 185–191, 1997.

Q. H. Qin and M. Lu, “BEM for crack-inclusion problems of plane thermopiezoelectric solids,” International Journal for Numerical Methods in Engineering, vol. 48, no. 7, pp. 1071–1088, 2000.

Q. H. Qin, “Thermoelectroelastic Green's function for a piezoelectric plate containing an elliptic hole,” Mechanics of Materials, vol. 30, no. 1, pp. 21–29, 1998.

Q. H. Qin, “2D Green's functions of defective magnetoelectroelastic solids under thermal loading,” Engineering Analysis with Boundary Elements, vol. 29, no. 6, pp. 577–585, 2005.

J. Wang, S. G. Mogilevskaya, and S. L. Crouch, “A numerical procedure for multiple circular holes and elastic inclusions in a finite domain with a circular boundary,” Computational Mechanics, vol. 32, no. 4–6, pp. 250–258, 2003.

K. N. Rajesh and B. N. Rao, “Two-dimensional analysis of anisotropic crack problems using coupled meshless and fractal finite element method,” International Journal of Fracture, vol. 164, no. 2, pp. 285–318, 2010.

Q. H. Qin and Y. W. Mai, “BEM for crack-hole problems in thermopiezoelectric materials,” Engineering Fracture Mechanics, vol. 69, no. 5, pp. 577–588, 2002.

C. Cao, A. Yu, and Q. H. Qin, “A novel hybrid finite element model for modeling anisotropic composites,” Finite Elements in Analysis and Design, vol. 64, pp. 36–47, 2013.

Q. H. Qin, Fracture Mechanics of Piezoelectric Materials, WIT Press, Southampton, UK, 2001.

Q. H. Qin, “Solving anti-plane problems of piezoelectric materials by the Trefftz finite element approach,” Computational Mechanics, vol. 31, no. 6, pp. 461–468, 2003.

H. Sosa, “Plane problems in piezoelectric media with defects,” International Journal of Solids and Structures, vol. 28, no. 4, pp. 491–505, 1991.

Q. H. Qin and S.-W. Yu, “An arbitrarily-oriented plane crack terminating at the interface between dissimilar piezoelectric materials,” International Journal of Solids and Structures, vol. 34, no. 5, pp. 581–590, 1997.

Q. H. Qin, Green's Function and Boundary Elements of Multifield Materials, Elsevier Science, Oxford, UK, 2007.

S. W. Yu and Q. H. Qin, “Damage analysis of thermopiezoelectric properties: part II. Effective crack model,” Theoretical and Applied Fracture Mechanics, vol. 25, no. 3, pp. 279–288, 1996.

Q. H. Qin and Y.-W. Mai, “Thermoelectroelastic Green's function and its application for bimaterial of piezoelectric materials,” Archive of Applied Mechanics, vol. 68, no. 6, pp. 433–444, 1998.

Q. H. Qin, “General solutions for thermopiezoelectrics with various holes under thermal loading,” International Journal of Solids and Structures, vol. 37, no. 39, pp. 5561–5578, 2000.

Q. H. Qin and X. Zhang, “Crack deflection at an interface between dissimilar piezoelectric materials,” International Journal of Fracture, vol. 102, no. 4, pp. 355–370, 2000.

S. W. Yu and Q. H. Qin, “Damage analysis of thermopiezoelectric properties: part I—crack tip singularities,” Theoretical and Applied Fracture Mechanics, vol. 25, no. 3, pp. 263–277, 1996.

C. Cao, Q. H. Qin, and A. Yu, “Hybrid fundamental-solution-based FEM for piezoelectric materials,” Computational Mechanics, vol. 50, no. 4, pp. 397–412, 2012.

C. Cao, A. Yu, and Q. H. Qin, “A new hybrid finite element approach for plane piezoelectricity with defects,” Acta Mechanica, vol. 224, no. 1, pp. 41–61, 2013.

H. Ding, G. Wang, and J. Liang, “General and fundamental solutions of plane piezoelectroelastic problem,” Acta Mechanica Sinica, vol. 28, pp. 441–448, 1996.

W. Yao and H. Wang, “Virtual boundary element integral method for 2-D piezoelectric media,” Finite Elements in Analysis and Design, vol. 41, no. 9-10, pp. 875–891, 2005.

C. Cao, Micro-Macro Modeling of Advanced Materials by Hybrid Finite Element Method, College of Engineering and Computer Science, Australian National University, Canberra, Australia, 2013.

W. G. Jin, N. Sheng, K. Y. Sze, and J. Li, “Trefftz indirect methods for plane piezoelectricity,” International Journal for Numerical Methods in Engineering, vol. 63, no. 1, pp. 139–158, 2005.

G. S. Fam and Y. F. Rashed, “The method of fundamental solutions applied to 3D structures with body forces using particular solutions,” Computational Mechanics, vol. 36, no. 4, pp. 245–254, 2005.

H. J. Al-Gahtani and F. M. Mukhtar, “RBF-based meshless method for the free vibration of beams on elastic foundations,” Applied Mathematics and Computation, vol. 249, pp. 198–208, 2014.

F. Toja-Silva, J. Favier, and A. Pinelli, “Radial basis function (RBF)-based interpolation and spreading for the immersed boundary method,” Computers & Fluids, vol. 105, pp. 66–75, 2014.