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 conditionswhere 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 aswith .
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:asEvidently, (5) analytically satisfies (1) due to the solution property of .
(a) 4-node 2D element
(b) 8-node 2D element
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 iswherewith
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 aswhich 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 formwhere stands for shape functions in terms of natural coordinate defined in Figure 2.
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]:
A modified variational functional is developed as follows:wherein 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 yieldsUsing divergence theoremwe can obtain For the displacement-based method, the potential conformity should be satisfied in advancethen (21) can be rewritten asThe Euler equation and boundary conditions can be obtained asusing 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 havefrom 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 findTherefore 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 asAppling the divergence theorem (20) to the functional (29), we have the final functional for the HFS-FE modelThen, substituting (5), (9), and (12) into the functional (30) producesin whichThe 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, yieldsfrom which the optional relationship between and and the stiffness equation can be producedwhere 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 asby introducing the matrix function Equation (36) can be further rewritten aswhereand is the Jacobean expressed as whereThus, the Gaussian numerical integration for matrix can be calculated bywhere 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 aswhere the undetermined rigid-body motion parameter can be calculated using the least square matching of and at element nodeswhich finally givesin 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 asIf written as matrix form, it can be presented aswhere is a stress vector, is a body force vector, and the differential operator matrix is given aswhere is a strain vector and is a displacement vector.
The constitutive equations for the linear elasticity are given in matrix form aswhere is the material coefficient matrix with constant components for isotropic materials, which can be expressed as follows: whereThe two different kinds of boundary conditions can be expressed aswhere 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 toDue 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 havewhere 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 nodeswhich finally giveswherein 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
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 satisfywith 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 assumewhere is the number of interpolation points, are the RBFs, and are the coefficients to be determined. Subsequently, the particular solution can be approximated bywhere 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 vectorsSubstituting (87) into (86), we can obtain the following biharmonic equation:Taking the Spline Type RBF as an example, we get the following solutions:whereand 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 aswhere . Substituting (90) into (92), we havewhere
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 asin 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)
According to (50) and (49), the corresponding stress fields can be expressed as where The components are given byAs 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).
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 asDue 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 bywhere is the number of nodes of the surface and are the nodal coordinates. Thus the unit normal vector is given bywhereis the Jacobian of the transformation from Cartesian coordinates () to natural coordinates ().
For the matrix, we introduce the matrix function Then we can getand we rewrite it to the component form as Using the relationship and the Gaussian numerical integration, we can obtainwhere 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 byIt 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 fieldand using a least-square procedure to match and at the nodes of the element boundaryThe above equation finally yieldswhere
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 satisfywith modified boundary conditionsFrom 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: wherefor two-dimensional problem and wherefor 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 bywhere . Substituting (147) into (149), one can obtainwherefor two-dimensional problem, and substituting (147) into (149), one can obtainwherefor three-dimensional problem.
To solve (141), one can treat as the gradient of a scalar function Substituting (154) into (141) obtains the Poisson’s equationThus, 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) intoThen 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 bywhere , 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 aswhere 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 bywhere the transformation matrix and its inverse matrix are defined aswith . Therefore, the constitutive relationship in the global coordinate system is given by
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 aswhere the boundary of the element is and is the interelement boundary of element . Performing a variation of , one obtainsApplying Gaussian theorem and the definitions of traction force we obtain Then, substituting (195) into (192) givesConsidering the fact thatwe 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 whereThe 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 byand the electric field component is related to the electric potential byThe boundary conditions of the boundary value problem (203)–(206) can be defined bywhere , 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 toFor 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 aswhere 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 aswhere 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 aswhereand 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 isApplying Gaussian theorem and the definitions of traction force and electrical displacementwe obtain Then, substituting (235) into (231) givesConsidering the fact thatwe 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 asDue 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 aswhere 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 aswhere the undetermined rigid-body motion parameter can be calculated using the least square matching of and at element nodeswhich finally givesin 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 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 , , , 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.
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].
(a)
(b)
(a)
(b)
(c)
(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 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 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 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)
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 beBecause 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 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)
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 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.
(a)
(b)
(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.
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)
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.
(a)
(b)
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.
(a)
(b)
(c)
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)
(a)
(b)
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.