1. science
2. mathematics
3. algebra

HYBRID FUNDAMENTAL SOLUTION BASED FINITE
ELEMENT METHOD: THEORY AND APPLICATIONS
Changyong Cao and Qing-Hua Qin
Research School of Engineering, The Australian National University, Canberra, ACT 2601,
Australia
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, threedimensional elasticity, thermo-elasticity, anisotropic elasticity, and plane piezoelectricity are
presented. In this method, two independent assumed fields (intra-element 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-12]. In this method, two independent assumed fields (intra-element 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 [13-16]. It should be mentioned that the intra-element field of HFSFEM is approximated by the linear combination of fundamental solutions analytically

Corresponding author: [email protected]
1
satisfying the related governing equation, instead of T-complete functions as in HT-FEM.
The resulting system of equations from the modified variational functional is expressed in
terms of symmetric stiffness matrix and nodal displacements only, which is easy to
implement into the standard FEM. It is noted that no singular integrals are involved in the
HFS-FEM by locating the source point outside the element of interest and not overlap with
field point during the computation [17].
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 [18, 19],
i.e. singular element boundary integral and loss of symmetry and sparsity [20, 21]. 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 HTFEM [22-25] in deriving T-complete functions for certain complex or new physical problems
[22, 26]. The HFS-FEM has simpler expressions of interpolation functions for intra-element
fields (fundamental solutions) and avoids the coordinate transformation procedure required in
the HT-FEM to keep the matrix inversion stable [22]. 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 accurately
representing various local effects (such as hole, crack and inclusions) without troublesome
mesh adjustment [7, 27-29]. Additionally, HFS-FEM makes it possible for a more flexible
element material definition which is important in dealing with multi-material 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 tends 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 is 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
2
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 thermo-elastic 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
 2u  x   0
x 
(1)
with the boundary conditions
u u
on u
q  u,i ni  q
on q
(2)
(3)
where u is the unknown field variable and q represents the boundary flux, ni is the ith
component of outward normal vector to the boundary   u  q , and u and q are
specified functions on the related boundaries, respectively. The space derivatives are
indicated by a subscript comma, i.e. u,i  u / xi , and the subscript index i takes values (1,2)
for two dimensional and (1,2,3) for three dimensional problems. Additionally, the repeated
subscript indices imply summation convention.
For convenience, Eq. (3) can be rewritten in matrix form as
 u,1 
q  A   q
u,2 
(4)
with A   n1 n2  .
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
3
problem defined by Eqs. (1)-(3). Similar to the conventional FEM and HT-FEM, the domain
under consideration is divided into a series of elements. In each element, two independent
fields are assumed in the way as described in [12] and are given in Section 2.2.
2.2.1. Intra-element field
Similar to the method of fundamental solution (MFS) in removing singularities of
fundamental solution, for a particular element e occupying sub-domain  e , we assume that
the field variable defined in the element domain is extracted from a linear combination of
fundamental solutions centered at different source points located outside the element (see
Figure 1)
ns
ue  x    Ne  x, y j  cej  Ne  x  ce
j 1
x e , y j e
(5)
where cej is undetermined coefficients, ns is the number of virtual sources and Ne (x, y j ) is
the fundamental solution to the partial differential equation
2 Ne  x, y     x, y  =0
x, y 2
(6)
as
Ne  x, y   
1
ln r  x, y 
2
(7)
Evidently, Eq. (5) analytically satisfies Eq. (1) due to the solution property of Ne (x, y j ) .
Frame field u(x)  N ed e
Frame field u(x)  N ed e
xs
xs
4
e
x0
3
X2
xc
7
e
x0
6
5
X2
8
X1
4
X1
e
e
1
2
1
xc
Source
2
3
Node
Centroid
Intra-element field
Figure 1
problems
u  N ece
Intra-element field u  Nece
(a) 4-node 2D element
(b) 8-node 2D element
Intra-element field and frame field of a HFS-FEM element for 2D potential
4
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 [22]. The source point y sj ( j  1, 2,
, ns ) can be
generated by means of the method employed in the MFS [30-32]
y s  x 0   (x 0  x c )
(8)
where  is a dimensionless coefficient, x 0 is the point on the element boundary (the nodal
point in this work)and x c the geometrical centroid of the element (see Figure 1).
Determination of  was discussed in [2, 33] and =5-10 is usually used in practice.
The corresponding outward normal derivative of ue on  e is
u
qe  e  Qece
n
where
Qe 
Ne
 ATe
n
(9)
(10)
with
 N
Te   e
 x1
N e 

x2 
T
(11)
2.2.2. Auxiliary frame field
In order to enforce the conformity on the field variable u , for instance, ue  u f on
e   f of any two neighboring elements e and f, an auxiliary inter-element frame field u is
used and expressed in terms of the same degrees of freedom (DOF) d as those used in the
conventional finite elements. In this case, u is confined to the whole element boundary as
(12)
ue  x   Ne  x  de
which is independently assumed along the element boundary in terms of nodal DOF d e ,
where Ne  x  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
u  N1u1  N2u2  N3u3
(13)
where N i ( i  1, 2,3 ) stands for shape functions in terms of natural coordinate  defined in
Figure 2.
5
Figure 2
Typical quadratic interpolation for frame field
2.3. MODIFIED VARIATIONAL PRINCIPLE
For the boundary value problem defined in Eqs. (1)-(3) and (5), since the stationary
conditions of the traditional potential or complementary variational functional cannot
guarantee the inter-element continuity condition required in the proposed HFS FE model, as
in the HT FEM, a variational functional corresponding to the new trial functions should be
constructed to assure the additional continuity across the common boundaries  Ief between
intra-element fields of element ‘e’ and element ‘f’ (see Figure 3):
ue  u f (conformity) 
 on  Ief  e   f
qe  q f  0 (reciprocity) 
(14)
 Ief
e
Figure 3
f
Illustration of continuity between two adjacent elements ‘e’ and ‘f’
A modified variational functional is developed as follows

 m    me    e  
e
e
where
6
e

 u  u  qd
(15)
e 
1
u,i u,i d   qud
qe
2 e
(16)
in which the governing equation (1) is assumed to be satisfied, a priori, for deriving the HFS
FE model. The boundary  e of a particular element consists of the following parts
e  ue  qe   Ie
(17)
where  Ie represents the inter-element boundary of the element ‘e’ shown in Figure 3.
To show that the stationary condition of the functional (15) leads to the governing
equation (Euler equation), boundary conditions and continuity conditions, invoking (16) and
(15) gives the following functional for the problem domain
1
 me   u,i u,i d   qud   q  u  u  d
qe
e
2 e
from which the first-order variational yields
 me   u,i u,i d   q ud  
e
qe
e
(18)
 u   u  qd   u  u   qd
(19)
e
Using divergence theorem


we can obtain
 me   u,ii ud  
e
qe
f,i h,i d   hf,i ni d   h2 fd

(20)

 q  q   ud   q ud   q ud   u  u   qd
ue
Ie
(21)
e
For the displacement-based method, the potential conformity should be satisfied in advance
u  0
on ue ( u  u )
(22)
 u e   u f on  Ief
( ue  u f )
then, Eq. (21) can be rewritten as
 me   u,ii ud  
e
qe
 q  q   ud   q ud   u  u   qd
Ie
(23)
e
The Euler equation and boundary conditions can be obtained as
u,ii  0
in e
qq
on  qe
u u
on  e
(24)
using the stationary condition  me  0 .
As for the continuous requirement between two adjacent elements ‘e’ and ‘f’ given in
Eq. (14), we can obtain it in the following way. When assembling elements ‘e’ and ‘f’, we
have
7
 m ( e f )  
e  f

u,ii ud  
f
 qe  qf
 q  q   ud    u  u   qd
 u  u   qd  
e
Ief
q
e
 q f   u ef d  ...
(25)
from which the vanishing variation of  m( e f ) leads to the reciprocity condition qe  q f  0
on the inter-element boundary  Ief .
If the following expression
  q uds  [
q
e
 Ie
 qe ue ds  
e
 qe (ue  ue )ds]
(26)
is uniformly positive (or negative) in the neighborhood of {u}0 , where the displacement {u}0
has such a value that  m ({u}0 )  ( m )0 , and where ( m )0 stands for the stationary value of
 m , we have
 m  ( m )0
or
 m  ( m )0
(27)
in which the relation that {u}e  {u} f is identical on e   f has been used. This is due to the
definition in Eq. (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” [34]. In doing this, performing variation of
 m and using the constrained conditions (26), we find
 2 m  
q
 q uds  [
e
 Ie
 qe ue ds  
e
 qe (ue  ue )ds]
(28)
Therefore the theorem has been proved from the sufficient condition of the existence of a
local extreme of a functional [34].
2.4. ELEMENT STIFFNESS EQUATION
With the intra-element 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 reference [22], the
variational functional  e corresponding to a particular element e of the present problem can
be written as
e 
1
u,i u,i d   qud   q  u  u  d
qe
e
2 e
(29)
Appling the divergence theorem (20) to the functional (29), we have the final functional for
the HFS-FE model
8
1
qud   uk 2ud    qud   q  u  u  d



e
e
e

 qe
2
1
   qud   qud   qud
 qe
e
2 e
e 
Then, substituting Eqs. (5), (9) and (12) into the functional (30) produces
1
 e   ceT Hece  deT g e  ceTG ede
2
(30)
(31)
in which
He   QeT Ne d   NeTQe d, G e   QeT Ne d, g e   NeT qd
e
e
e
qe
(32)
The symmetry of H e is obvious from the scalar definition (31) of variational functional  e .
To enforce inter-element continuity on the common element boundary, the unknown
vector c e should be expressed in terms of nodal DOF d e . The minimization of the functional
 e with respect to c e and d e , respectively, yields
 e
 H ece  G ed e  0
ce T
 e
 G Te ce  g e  0
d e T
(33)
from which the optional relationship between c e and d e , and the stiffness equation can be
produced
ce = He1G ede
and
K ed e = g e
(34)
where K e = G eT He1G e stands for the element stiffness matrix.
2.5. NUMERICAL INTEGRAL FOR H AND G MATRIX
Generally, it is difficult to obtain the analytical expression of the integral in Eq. (32)
and numerical integration along the element boundary is required. Herein the widely used
Gaussian integration is employed.
For the He matrix, one can express it as
He   QTe Ne d   F(x)d
e
e
by introducing the matrix function
F(x)   Fij (x) 
 QeT Ne
mm
Eq. (36) can be further rewritten as
9
(35)
(36)
ne
H ij   Fij  x  d    Fij  x  d
e
l 1
el
(37)
where
d 
 dx1    dx2 
2
2
2
2
 dx   dx 
  1    2  d  Jd
 d   d 
(38)
and J is the Jacobean expressed as
2
2
 dx   dx 
J   1   2  .
 d   d 
(39)
no
dNi    x1i 
 dx1 dx2 
,

 

 d d 
d  x2i 
i 1


(40)
where
T
Thus, the Gaussian numerical integration for H matrix can be calculated by
ne
ne
 np

1
H ij     Fij  x    J   d     wk Fij  x  k   J  k  
 1
 l 1  k 1
l 1 

(41)
where ne is the number of edges of the element and n p is the Gaussian sampling points
employed in the Gaussian numerical integration. Similarly, we can calculate the G e matrix
using
ne
ne
 np

1


Gij     Fij  x( ) J ( )d     wk Fij  x( k ) J ( k ) 
 1
 l 1 
l 1 
 k 1


(42)
The calculation of vector g e in Eq. (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 Eqs. (9). The boundary DOF can be directly
computed from Eq. (12) while the unknown variable at interior points of the element can be
determined from Eq. (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’s necessary to
recover the missing rigid-body motion modes from the above results. Following the method
presented in [22], the missing rigid-body motion can be recovered by writing the internal
potential field of a particular element e as
ue  Nece  c0
10
(43)
where the undetermined rigid-body motion parameter c0 can be calculated using the least
square matching of ue and ue at element nodes
n
N c
i 1
e e
 c0  ue 
 min
2
(44)
node i
which finally gives
c0 
1 n
 uei
n i 1
(45)
in which uei   ue  Nece  node i and n is the number of element nodes. Once the nodal field
is determined by solving the final stiffness equation, the coefficient vector c e can be
evaluated from Eq. (34), and then c0 is evaluated from Eq. (45). Finally, the potential field u
at any internal point in an element can be obtained by means of Eq. (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 [35]. In the rectangular Cartesian coordinates
(X1, X2), the governing equations of a plane elasticitc body can be expressed as
 ij , j  bi ,
i, j  1, 2
(46)
If written as matrix form, it can be presented as
Lσ  b
(47)
where σ  11  22 12  is a stress vector, b  b1 , b2 
T
T
differential operator matrix L is given as
 
 x
1
L

 0

ε  11  22 12 
T
a body force vector, and the
 
x2 

 
x1 
0

x2
ε = LTu
(48)
(49)
is a strain vector, and u  u1 , u2  a displacement vector.
T
The constitutive equations for the linear elasticity are given in matrix form as
σ  Dε
11
(50)
where D is the material coefficient matrix with constant components for isotropic materials,
which can be expressed as follow
  2G


D 
  2G
 0
0

where

2
G
1  2
G
,
0

0
G 
(51)
E
2 1  
(52)
and
for plane strain
v

  
1  for plane stress
The two different kinds of boundary conditions can be expressed as
uu
on u
t  Aσ  t
on t
(53)
(54)
where t  [t1 t2 ]T denotes the traction vector and A is a transformation matrix related to the
direction cosine of the outward normal
n
A 1
0
0
n2
n2 
n1 
(55)
Substituting Eqs. (49) and (50) into Eq. (47) yields the well-known Navier partial differential
equations in terms of displacements
LDLTu  b
(56)
3.2. ASSUMED INDEPENDENT FIELDS
For elasticity problem, two different assumed fields are employed as in potential
problems: intra-element and frame field. The intra-element continuity is enforced on
nonconforming internal displacement field chosen as the fundamental solution of the problem
[2]. The intra-element displacement fields is approximated in terms of a linear combination of
fundamental solutions of the problem of interest
*
*
u1 (x)  ns  u11 (x, y sj ) u12 (x, y sj )   c1 j 
u(x)  
    N e ce
   *
*
u2 (x)  j 1 u21 (x, y sj ) u22 (x, y sj )  c2 j 
(x e ,y sj e ) (57)
where ns is again the number of source points outside the element domain, which is equal to
the number of nodes of an element in the present research based on the generation approach
12
of the source points [33], The vector ce and the fundamental solution matrix N e , are now in
the form
*
 u11
(x, y s1 ) u12* (x, y s1 ) ... u11* (x, y sns ) u12* (x, y sns ) 
Ne   *

*
*
*
u21 (x, y s1 ) u22 (x, y s1 ) ... u21 (x, y sns ) u22 (x, y sns ) 
(58)
ce  [c11 c21 ...... c1n c2 n ]T
(59)
in which x and y sj are respectively the field point and source point in the local coordinate
system (X1, X2). The components uij* (x, y sj ) is the fundamental solution, i.e. induced
displacement component in i-direction at the field point x due to a unit point load applied in jdirection at the source point y sj , which are given by [19, 36]
1
uij* (x, y sj ) 
(3  4 )ij ln r  r,i r, j 
8 (1  )G
(60)
where ri  xi  xis , r  r12  r22 . The virtual source points for elasticity problems is generated
in same manner as that in potential problems described in Section 2.
With the assumption of intra-element field in Eq. (57), the corresponding stress fields
can be obtained by the constitutive Eq. (50)
σ(x)  11  22 12   Tece
(61)
*
*
*
*
 111
(x, y1 )  211
(x, y1 ) ......  111
(x, y ns )  211
(x, y ns ) 
 *

*
*
*
Te   122 (x, y1 )  222
(x, y1 ) ......  122
(x, y ns )  222
(x, y ns ) 
 *

*
*
*
 112 (x, y1 )  212 (x, y1 ) ......  112 (x, y ns )  212 (x, y ns ) 
(62)
T
where
As a consequence, the traction is written as
t1 
   nσ  Qece
t2 
(63)
in which
n
Qe  nTe , n   1
0
0
n2
n2 
n1 
*
The components  ijk
(x, y ) for plane strain problems is given as
1
*
(1  2 )(r,k  ij  r, j ki  r,i jk )  2r,i r, j r,k 
 ijk
(x, y ) 
4 (1  )r 
(64)
(65)
The unknown ce in Eq. (57) is calculated by a hybrid technique [33], 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
13
conforming displacement field should be independently defined on the element boundary to
enforce the field continuity between elements and also to link the unknown c and with nodal
displacement d e . Thus, the frame is defined as
N1 

u1  
u ( x)       d e  N e d e ,
N 2 

u2  
(x e )
(66)
where the symbol “~” is used to specify that the field is defined on the element boundary
only, N e is the matrix of shape functions, d e is the nodal displacements of elements. Taking
the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N e
and d e can be expressed as
0
N e  0


0
N1
0
N2
0
N3
0
0
0
0
N1
0
N2
0
N3
0
4
6
de  u11 u21 u12 u22
0

0

u18 u28 
T
where N1 , N 2 and N 3 can be expressed by natural coordinate as
 1   
 1   
N1  
, N 2  1   2 , N3 
  1,1
2
2
(67)
(68)
(69)
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 [2]. In the absence of body forces, the hybrid
functional  me used for deriving the present HFS-FEM can be constructed as [21]
1
 me    ij ij d   ti ui d    ti (ui  ui )d 
t
e
2 e
(70)
where ui and ui are the intra-element displacement field defined within the element and the
frame displacement field defined on the element boundary, respectively. e and  e are the
element domain and element boundary, respectively.  t ,  u and  I stands respectively for
the specified traction boundary, specified displacement boundary and inter-element boundary
( e  t  u   I ). Compared to the functional employed in the conventional FEM, the
present variational functional is constructed by adding a hybrid integral term related to the
intra-element and element frame displacement fields to guarantee the satisfaction of
displacement and traction continuity conditions on the common boundary of two adjacent
elements.
14
By applying the Gaussian theorem, Eq. (70) can be simplified to
1
 me    ij ij d    ti ui d    ti (ui  ui )d 
t
e
2 e
1
 (  ti ui d     ij , j ui d )   tiui d    ti (ui  ui )d 
e
t
e
2 e
(71)
Due to the satisfaction of the equilibrium equation with the constructed intra-element field,
we have the following expression for HFS finite element model
1
 me   ti ui d    ti ui d    ti (ui  ui )d 
t
e
2 e
1
   ti ui d    ti ui d    ti ui d 
e
t
2 e
(72)
The variational functional in Eq. (72) 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  me  0 .
Substituting Eqs. (57), (66) and (63) into the functional of Eq. (72), we have
1
 me   ce T Hece  ceTG ede  de T g e
2
(73)
where
He   QeT Ne d,
e
G e   QeT Ne d,
e
ge   NeT t d 
t
(74)
To enforce inter-element continuity on the common element boundary, the unknown vector c e
should be expressed in terms of nodal DOF d e . The stationary condition of the functional
 me with respect to c e and d e yields, respectively, Eqs. (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 [21]
1 0 x2 
u e  N ece  
 c0
0 1  x1 
(75)
where the undetermined rigid-body motion parameter c0 can be calculated using the least
square matching of ue and ue at element nodes
15
n
  u
i 1
1i
2
2
 u1i    u2i  u2i    min

(76)
which finally gives
c0 = R e1re
(77)
where
0
x2i 
1
R e =   0
1
 x1i 
i 1
 x2i  x1i x12i  x22i 
ue1i


n


re   
ue 2i

i 1 

ue1i x2i  ue 2i x1i 
n

in which ueji  ueji  uˆeji

(78)
(79)
(j  1, 2) and n is the number of element nodes. Once the
nodal field is determined by solving the final stiffness equation, the coefficient vector c e can
be evaluated from Eq. (57), and then c0 is evaluated from Eq. (77). Finally, the displacement
field u e at any internal point in an element can be obtained by Eq. (75).
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, then a procedure based on the method
of particular solution and radial basis function approximation are introduced to deal with the
body force. As a consequence, the homogeneous solution is obtained by using the HFS-FEM
and the particular solution associated with body force is approximated by using the strong
form of basis function interpolation.
4.1. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS
Let (x1, x2, x3) denote the coordinates in a Cartesian coordinate system and consider a
finite isotropic body occupying the domain  , as shown in Figure 4. The equilibrium
equation for this finite body with body force can be expressed as
 ij , j  bi
i, j  1, 2,3
16
(80)
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
2Gv
 ij ekk  2Geij
1  2v
1
eij  (ui , j  u j ,i )
2
 ij 
(81)
(82)
where  ij is the stress tensor, eij the strain tensor,  ij the Kronecker delta. Substituting Eqs.
(81) and (82) into Eq. (80), the equilibrium equation is rewritten as
G
Gui , jj 
u j , ji  bi
1  2v
(83)
For a well-posed boundary value problem, following boundary conditions, are
prescribed as
ui  ui on u ,
ti  ti on t ,
where u
(84)
(85)
t   is the boundary of the solution domain  , ui and ti are the prescribed
boundary values. In the following parts, we will present the procedure for handling the body
force appearing in Eq. (83).
4.2. THE METHOD OF PARTICULAR SOLUTION
The inhomogeneous term bi associated with the body force in Eq. (83) can be
effectively handled by means of the method of particular solution presented in [21]. In this
approach, the displacement ui is decomposed into two parts, a homogeneous solution uih and
a particular solution uip
17
ui  uih  uip
(86)
where the particular solution uip should satisfy the governing equation
G
Guip, jj 
u jp, ji  bi
1  2v
(87)
without any restriction of boundary condition. However, the homogeneous solution should
satisfy
Guih, jj 
G h
u j , ji  0
1  2v
(88)
with the modified boundary conditions
_
uih  ui  uip
_
tih  ti  tip
on u
(89)
on t
(90)
From above equations it can be seen that once the particular solution uip is known, the
homogeneous
solution
uih in
Eqs.
(88)Error!
Reference
source
not
found.-
Error! Reference source not found.(90) can be obtained using HFS-FEM. The final
solution can then be given by Eq. (86)Error! Reference source not found.. In the next
section, radial basis function approximation is introduced to obtain the particular solution,
and
the
HFS-FEM
is
presented
for
solving
Eqs.
(88)
-
(90)
Error! Reference source not found..Error! Reference source not found.
4.3. RADIAL BASIS FUNCTION APPROXIMATION
For body force bi , it is generally impossible to find an analytical solution which
enable 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) [31, 37] to interpolate the body force. We assume
N
bi    i j j
(91)
j 1
where N is the number of interpolation points,  j are the RBFs and  i j are the coefficients to
be determined. Subsequently, the particular solution can be approximated by
N
uip    i j ikj
j 1
18
(92)
where  ikj is the approximated particular solution kernel of displacement satisfying Eq. (93)
below. Once the RBFs are selected, the problem of finding a particular solution is reduced to
solve the following equation
Gil ,kk 
G
 kl ,ki   il
1  2
(93)
To solve this equation, the displacement is expressed in terms of the Galerkin Papkovich vectors
1 
1
Fik ,mm 
Fmk ,mi
(94)
G
2G
Substituting Eq. Error! Reference source not found.(94) into Eq. (93)
Error! Reference source not found., we can obtain the following bi-harmonic equation
1
4 Fil  
 il
(95)
1 
ik 
Taking the Spline Type RBF   r 2n1 as an example, we get the following solutions
Fli  
 li
r 2 n 3
1  v  2n  1 2n  2  2n  3 2n  4 
li  A0  A1 li  A2 r,i r,l 
(96)
(97)
where
A0  
1
r 2 n 1
8G 1  v   n  1 n  2  2n  1
A1  7  4n  4v(n  2)
(98)
A2  (2n  1)
and rj represents the Euclidean distance between a field point ( x1 , x2 , x3 ) and a given point
( x1 j , x2 j , x3 j ) in the domain of interest. The corresponding particular solution of stresses can
be obtained as
where  
Slij  G(li , j  lj ,i )  ij lk ,k
(99)
2v
G . Substituting Eq.Error! Reference source not found. (97) into
1  2v
Eq.Error! Reference source not found. (99), we have

Slij  B0 B1  r, j li  r,i jl   B2 ij r,l  B3r,i r, j r,l
where
19

(100)
B0  
1
r 2n
4 1  v   n  1 n  2 
B1  3  2n  2v(n  2)
(101)
B2  2v(n  2)  1
B3  1  2n
4.4. HFS-FEM FOR HOMOGENEOUS SOLUTION
After obtained the particular solution in sections 4.2 and 4.3, we can determine the
modified boundary conditions Error! Reference source not found.(89) and (90)
Error! Reference source not found.. Finally, we can treat the 3D problem as a
homogeneous problem governed by Eqs. (88)-(90) 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, i.e. bi  0 , the procedures in
Sections 4.2 and 4.3 will become unnecessary and we can employ the following procedures
to find the solution directly.
4.4.1. Assumed intra-element and auxiliary frame fields
The intra-element displacement fields is approximated in terms of a linear
combination of fundamental solutions of the problem as
 u1 (x) 


u(x)  u2 (x)   Nece (x  e ,y sj  e )
u ( x) 
 3 
(102)
where the matrix N e and unknown vector ce can be further written as
*
*
*
u11
(x, y s1 ) u12
(x, y s1 ) u13
(x, y s1 )
 *
*
*
Ne  u12 (x, y s1 ) u22
(x, y s1 ) u23
(x, y s1 )
 *
*
*
u13 (x, y s1 ) u32 (x, y s1 ) u33 ( x, y s1 )
*
u13
( x, y sns ) 

*
*
*
u12
( x, y sns ) u22
( x, y sns ) u23
( x, y sns ) 

*
*
*
u13
( x, y sns ) u32
( x, y sns ) u33
(x, y sns ) 
(103)
c1n c2 n c3n ]T
(104)
*
*
u11
( x, y sns ) u12
( x, y sns )
ce  [c11 c21 c31
in which x and y sj are respectively the field point and source point in the local coordinate
system (x1, x2). The fundamental solution uij* (x, y sj ) are given by [36]
20
ui*j (x, y sj ) 
1
(3  4 )ij  r,i r, j 
16 (1  )Gr
(105)
where ri  xi  xis , r  r12  r22  r32 , ns is the number of source points. The source point
y sj ( j  1, 2,
, ns ) can also be generated by means of the following method [2] as in two-
dimensional cases
y s  x 0   (x 0  x c )
(106)
where  is a dimensionless coefficient, x 0 the point on the element boundary (the nodal point
in this work) and x c the geometrical centroid of the element (see Figure 5).
Intra-element field
u  N ece
Source points
Node
Centroid
8
19
7
e
5
Frame field
u ( x)  N e d e
xs
3
4
X2
11
10
5
8
X3
X2
4
1
2
3
2
20-node 3D element
8-node 3D element
Figure 5
6
7
9
15
X1
X1
1 x0
17
16
12
xc
e
e
14
13
X3
6
18
20
Intra-element 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 Eqs. (50) and (49), the corresponding stress fields can be expressed as
σ(x)  11  22  33  23  31 12   Tece
T
(107)
where
*
 111
(x, y1 )
 *
 122 (x, y1 )
 *
 133 (x, y1 )
Te   *
 (x, y )
1
 123
*
 131 (x, y1 )
 *
 112 (x, y1 )
*
*
 211
(x, y1 )  311
(x, y1 )
*
*
*
 111
(x, y n )  211
(x, y n )  311
( x, y n ) 

 (x, y n ) 
*
222
(x, y1 ) 
*
322
s
*
122
(x, y1 )
s
s
*
222
*
*
 233
(x, y1 )  333
(x, y1 )
*
*
 133
( x, y n )  233
( x, y n )

(x, y1 )
 (x, y n )  ( x, y n )
*
331
 (x, y1 )  (x, y1 )
 (x, y n )  (x, y n )

*
312
 (x, y n ) 
*
223
(x, y1 ) 
*
231
*
212
(x, y1 ) 
*
323
s
*
123
*
131
*
112
(x, y1 )
s
*
223
s
s
s
*
231
s
*
212
s
(x, y ns )
(108)
*
The components  ijk
(x, y ) is given by
21

( x, y ns ) 

*
 333
(x, y ns ) 
*
 323
( x, y ns ) 

*
 331
( x, y ns ) 

*
 312
(x, y ns ) 
(x, y ns ) 
s
*
322
*
 ijk
(x, y ) 
1
(1  2 )(r,kij  r, j ki  r,i jk )  3r,i r, j r,k 
8 (1  )r 2
As a consequence, the traction can be written in the form
t1 
 
t2   nσ  Qece
t 
 3
(109)
(110)
in which
 n1
Qe  nTe , n   0

 0
0
0
0
n3
n2
0
n3
0
0
n3
n2
n1
n2 
n1 
0 
(111)
To link the unknown c e and the nodal displacement d e , the frame is defined as
 u1  N1 
   
u(x)  u2   N 2  d e  N ed e ,
u   
 3  N3 
(x   e )
(112)
where the symbol “~” is used to specify that the field is defined on the element boundary
only, N e is the matrix of shape functions, d e is the nodal displacements of elements. Taking
the surface 2-3-7-6 of a particular 8-node brick element (see Figure 5) as an example, matrix
N e and vector d e can be expressed as
Ne  0 N1 N2 0 0 N4
de = u11 u21 u31 u12 u22 u32
where the shape functions are expressed as
 Ni 0 0 


Ni   0 Ni 0  ,
 0 0 Ni 


N3 0
(113)
u18 u28 u38 
T
0 0 0 
0  0 0 0 
0 0 0 
(114)
(115)
where N i (i=1-4) can be expressed by natural coordinate  , [1,1]
N1 
N3
1   1    ,
N2 
4
1   1   ,

N4
4
1   1   
4
1   1  

4
and i ,i  is the natural coordinate of the i-node of the element (Figure 6).
22
(116)

(1,1) 4
3 (1,1)
o
1
2
(1, 1)
Figure 6

(1, 1)
Typical linear interpolation for the frame fields of 3D brick elements
4.4.2. Modified Functional for Hybrid Finite Element Method
In the absence of body forces, the hybrid functional  me used for deriving the present
HFS-FEM is written as [21]
 me 
1
 ij ij d   tiui d    ti (ui  ui )d 
t
e
2 e
By applying the Gaussian theorem, Eq. (117) can be simplified as
1
 me    ij ij d    ti ui d    ti (ui  ui )d 
t
e
2 e
1
 (  ti ui d     ij , j ui d )   tiui d    ti (ui  ui )d 
e
t
e
2 e
(117)
(118)
Due to the satisfaction of the equilibrium equation with the constructed intra-element fields,
we have the following expression for HFS finite element model
1
 me   ti ui d    ti ui d    ti (ui  ui )d 
t
e
2 e
1
   ti ui d    ti ui d    ti ui d 
e
t
2 e
(119)
The functional (119) 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.
23
4.4.3. Element stiffness matrix
Substituting Eqs. (102), (110) and (112) into the functional (119), we have
1
 me   ce T Hece  ceTG ede  de T g e
(120)
2
where
He   QeT Ne d,
e
G e   QeT Ne d,
e
ge   NeT t d 
t
(121)
To enforce inter-element continuity on the common element boundary, the unknown vector
c e should be expressed in terms of nodal DOF d e . The stationary condition of the functional
 me with respect to c e and d e yields again, respectively, Eqs. (33) and (34).
4.4.4. Numerical integral over element
Consider a surface of the 3D hexahedron element, as shown in Figure 6, the vector
normal to the surface can be obtained by
 dx   dx   dy
 d   d   d
vnx      
   dy   dy   dz
vn  v  v  vny         
   d   d   d
vnz   dz   dz   dx
    
 d   d   d
dz dy dz 

d d d 
dx dz dx 


d d d 
dy dx dy 


d d d 
(122)
where v and v are the tangential vectors in the  -direction and  -direction, respectively,
calculated by
 dx 
 dx 
 d 
 d 
 
 
x


 xi 
i
 dy  nd Ni ( , )  
 dy  nd Ni ( , )  
v     
 yi  , v     
 yi 


 d  i 1
z 
 d  i 1
z 
 i
 i
 dz 
 dz 
 
 
 d 
 d 
(123)
where nd is the number of nodes of the surface,  xi , yi , zi  are the nodal coordinates. Thus
the unit normal vector is given by
n
vn
vn
(124)
where
J ( , )  vn  vnx2  vny2  vnz2
24
(125)
is the Jacobian of the transformation from Cartesian coordinates (x, y) to natural coordinates (
 , ).
For the H matrix, we introduce the matrix function
F(x, y)   Fij ( x, y) 
 QeT Ne
mm
(126)
Then we can get
He   QTe Ne d   F(x, y)d
e
e
(127)
and we rewrite it to the component form as
nf
H ij   Fij ( x, y)dS    Fij ( x, y)dS
(128)
dS=J ( , )d d
(129)
e
Using the relationship
l 1
el
and the Gaussian numerical integration, we can obtain
nf
H ij    Fij  x( , ), y ( , )  J ( , )d d
l 1
1
1
n
n
 p p

   ws wt Fij  x( s ,t ), y ( s ,t )  J ( s ,t ) 
l 1 
 s 1 t 1

nf
(130)
where n f and n p 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 G e matrix by
nf
Gij    Fij  x( , ), y ( , )  J ( , )d d
l 1
1
1
n
n
 p p

   ws wt Fij  x( s ,t ), y ( s ,t )  J ( s ,t ) 
l 1 

 s 1 t 1
nf
(131)
It should be mentioned that the calculation of vector g e in Eq. 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 Eqs. (107) and (108), respectively. The boundary displacements can be directly
computed from Eq. (112) while the displacements at interior points of element can be
determined from Eq.(102) plus the recovered rigid-body modes in each element, which is
introduced in the following section.
25
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
1 0 0 0
ue  N ece  0 1 0  x3
0 0 1 x2
 x2 
x1  c0
0 
x3
0
 x1
(132)
and using a least-square procedure to match ueh and ueh at the nodes of the element boundary
n
2
2
2
min    u1i  u1i    u2i  u2i    u3i  u3i  


i 1
(133)
The above equation finally yields
c0 = R e1re
(134)
where
 1
 0

n  0
Re =  
i 1  0
 x3i

  x2i
0
0
0
x3i
1
0
 x3i
0
0
1
x2i
 x1i
 x3i
x2i
x22i  x32i
 x1i x2i
0
 x1i
 x1i x2i
x12i  x32i
x1i
0
 x1i x3i
 x2i x3i
ue1i




ue 2i


n 

ue3i
re   

i 1  ue 3i x2 i  ue 2 i x3i 
 ue1i x3i  ue3i x1i 


 ue 2i x1i  ue1i x2i 
 x2i 
x1i 
0 

 x1i x3i 
 x2i x3i 

x12i  x22i 
(135)
(136)
5. Thermo-elasticity 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 [38]. In many
engineering applications, the coupling term of the heat equation and the inertia term in Navier
equation are generally negligible [38]. As a consequence, most of the analyses are employing
26
the uncoupled thermoelasticity theory which is adopted in this work [39-46]. In this section,
the HFS-FEM is presneted to solve 2D and 3D thermoelastic problems with considering
arbitrary body forces and temperature changes. The method used herein is similar to that in
Section 4.
5.1. BASIC EQUATIONS FOR THERMOELASTICITY
Consider an isotropic material in a finite domain  , and let (x1, x2, x3) denote the
coordinates in Cartesian coordinate system. The equilibrium governing equations of the
thermoelasticity with the body force are expressed as
 ij , j  bi
(137)
where  ij is the stress tensor, bi the body force vector and i, j  1, 2,3 . The generalized
thermoelastic stress-strain relations and kinematical relation are given as
2G
 ij , j 
 ij e  2Geij  m ijT
1  2
1
eij  (ui , j  u j ,i )
2
(138)
(139)
where eij is the strain tensor, ui the displacement vector, T the temperature change, G the
shear modulus,  the Poisson’s ratio,  ij the Kronecker delta, and
m  2G (1  v) / (1  2v)
(140)
is the thermal constant with  being the coefficient of linear thermal expansion. Substituting
Eqs. (138) and (139) into Eq. (137), the equilibrium equations may be rewritten as
G
Gui , jj 
u j , ji  mT,i  bi
(141)
1  2
For a well-posed boundary value problem, the following boundary conditions, either
displacement or traction boundary condition, should be prescribed as
ui  ui
on u ,
ti  ti
where u
on t ,
(142)
(143)
t   is the boundary of the solution domain  , ui and ti are the prescribed
boundary values, and
ti   ij n j ,
is the boundary traction, in which n j denotes the boundary outward normal.
27
(144)
5.2. THE METHOD OF PARTICULAR SOLUTION
For the governing equation (141) given in the previous section, the inhomogeneous
term mT,i  bi can be eliminated by employing the method of particular solution [2, 26, 38].
Using superposition principle, we decompose the displacement ui into two parts, the
homogeneous solution uih and the particular solution uip as follows
ui  uih  uip
(145)
in which the particular solution uip should satisfy the governing equation
G
Guip, jj 
u jp, ji  mT,i  bi
1  2
(146)
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
Guih, jj 
G
u hj , ji  0
1  2
(147)
with modified boundary conditions
_
uih  ui  uip ,
on u
_
tih  ti  mTni  tip ,
(148)
on t
(149)
From above equations, it can be seen that once the particular solution is known, the
homogeneous solution uih
Error!
Reference
source
in Eqs. (147)Error! Reference source
not
found.(149)
can
be
solved
by
not found.Eq.
(147)
Error! Reference source not found.. 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 the Eqs. (147)-(149)Error! Reference source not found..
5.3. RADIAL BASIS FUNCTION APPROXIMATION
RBF is to be used to approximate the body force bi and the temperature field T in
order to obtain the particular solution. To implement this approximation, we may consider
two different ways: one is to treat body force bi and the temperature field T separately as in
by Tsai [47]. The other is to treat mT,i  bi as a whole [3]. Here we demonstrated that the
performance of the latter one is usually better than the former one.
28
5.3.1. Interpolating temperature and body force separately
The body force bi and temperature T are assumed to be by the following two
equations
N
bi  i j j
(i =1,2 in 2 and i=1,2,3 in 3 )
(150)
j 1
and
N
T    j j
(151)
j 1
where N is the number of interpolation points,  j are the basis functions, and  i j and  j are
the coefficients to be determined by collocation. Subsequently, the approximate particular
solution can be written as follows
3
N
N
uip    i j ikj    j  ij
k 1 j 1
(152)
j 1
where  ikj and  ij is 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
G
Gil ,kk 
 kl ,ki   il
(153)
1  2
G
G i ,kk 
 k ,ki  m,i
(154)
1  2
To solve Eq.(153), the displacement is expressed in terms of the Galerkin -Papkovich
vectors [37, 48]
ik 
1 
1
Fik ,mm 
Fmk ,mi
G
2G
(155)
Substituting Eq. (155) into Eq. (153), one can obtain the following bi-harmonic equation
1
4 Fil  
 il
(156)
1 
If taking the Spline Type RBF   r 2n1 , one can get the following solutions
Fli  
 li
r 2 n 3
1  v  2n  1  2n  3
2
2
(2 ) for n=1,2,3
li  A0  A1 li  A2 r,i r,l  (2 ) for n=1,2,3
(157)
(158)
where
A0  
1
r 2 n 1
2G 1  v   2n  12  2n  3
A1  5  4n  2v(2n  3)
A2  (2n  1)
29
(159)
for two-dimensional problem and
Fli  
 li
r 2 n 3
(3 ) for n=1,2,3
1  v  2n  1 2n  2  2n  3 2n  4 
li  A0  A1 li  A2 r,i r,l  (3 ) for n=1,2,3
(160)
(161)
where
1
r 2 n 1
A0  
8G 1  v   n  1 n  2  2n  1
A1  7  4n  4v(n  2)
(162)
A2  (2n  1)
for three-dimensional problem, where rj represents the Euclidean distance of the given point
( x1 , x2 ) from a fixed point ( x1 j , x2 j ) in the domain of interest. The corresponding stress
particular solution can be obtained by
Slij  G(li , j  lj ,i )  ij lk ,k
where  
(163)
2v
G . Substituting Eq. (161) into (163)Error! Reference source not found.,
1  2v
one can obtain

Slij  B0 B1  r, j li  r,i jl   B2 ij r,l  B3r,i r, j r,l

(2 ) for n=1,2,3
(164)
where
B0  
1
r 2n
1  v   2n  1 2n  3
B1  (2n  2)  v(2n  3)
(165)
B2  v(2n  3)  1
B3  1  2n
for
two-dimensional
problem,
and
substituting
Error! Reference source not found., one can obtain

Slij  B0 B1  r, j li  r,i jl   B2 ij r,l  B3r,i r, j r,l

Eq.
(161)
(3 ) for n=1,2,3
into
(163)
(166)
where
1
r 2n
B0  
4 1  v   n  1 n  2 
B1  3  2n  2v(n  2)
B2  2v(n  2)  1
B3  1  2n
for three-dimensional problem.
30
(167)
To solve Eq. (154), one can treat  i as the gradient of a scalar function
 i  U ,i
Substituting Eq. (168) into Eq. (154) obtains the Poisson’s equation
m(1  2 )
 2U 

2G(1  )
(168)
(169)
Thus, taking   r 2n1 , its particular solution can be obtained [48]
m(1  2 ) r 2 n1
(2 ) for n=1,2,3
2
2G(1  ) (2n  1)
m(1  2 )
r 2 n1
U
(3 ) for n=1,2,3
2G(1  ) (2n  1)(2n  2)
U
(170)
(171)
Then from Eq. (168) we can get  i as follows
2n
m(1  2 ) r,i r
(2 ) for n=1,2,3
2G(1  ) 2n  1
2n
m(1  2 ) r,i r
i 
(3 ) for n=1,2,3
2G(1  ) (2n  2)
i 
(172)
(173)
The corresponding stress particular solution can be obtained by substituting Eq. (161) into
Sij  G(i , j   j ,i )  ij k ,k
(174)
Then we have
mr 2 n1
(1  2nv) ij  (1  2v)(2n  1)r,i r, j  (2 ) for n=1,2,3

(1  v)  2n  1
(175)
mr 2 n1
(1  2 )(2n  1)r,i r, j   ij (1  2nv) (3 ) for n=1,2,3

1  v  2n  2 
(176)
Sij 
Sij 
5.3.2. Interpolating temperature and body force together
Considering the temperature gradient plays the role of body force, we can
approximate mT,i  bi together by the following equation
N
mT,i  bi    i j j
(177)
j 1
Thus the approximate particular solution can be written as
3
N
uip    i j ikj
(178)
k 1 j 1
Consequently, we can follow the same way as that for body force in Method 1 and employing
Eq. (153) and Eqs. (155)-(164) to obtain the desired  ikj and Slij , which are the same as those
for body force case only. It is noted that Method 2 has a relatively smaller number of
31
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 Eq. (145), we can use them to get
the modified boundary conditions in Eq.Error! Reference source not found. (148) to obtain
the homogeneous solution by considering Eq.Error! Reference source not found. (147).
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 [49]. Individual layers (or laminas) of the laminates consist of high-modulus, highstrength 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 macro-mechanical behavior of composites [49].
In the literature, there are two main approaches dealing with generalized twodimensional anisotropic elastic problems. One is Lekhnitskii formalism [50, 51] which begins
with stresses as basic variables and the other is Stroh formalism [52, 53] 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 [52, 54]. The formalism is
also widely employed in the derivation of the inclusion or crack problems of anisotropic
materials [55, 56]. Because of the limitations of the analytical solutions which are only
available for some problems with simple geometry and boundary conditions [57, 58],
32
numerical methods such as finite element method (FEM), boundary element method (BEM),
mesh free method (MFM), Hybrid-Trefftz (HT) FEM are usually resorted to solve more
complex problems with complicated boundary constraints and loading conditions [59-61]. In
this section, we presented the HFS-FEM for analyzing anisotropic composite materials based
on the associated fundamental solutions in terms of Stroh formalism.
6.1. LINEAR ANISOTROPIC ELASTICITY
6.1.1. Basic equations and Stroh Formalism
In the Cartesian coordinate system (x1, x2, x3), if we neglect the body force bi , the
equilibrium equations, stress-strain laws and strain-displacement equations for anisotropic
elasticity are [52]
 ij , j  0
(179)
 ij  Cijkl ekl
1
eij  (ui , j  u j ,i )
2
(180)
(181)
where i, j  1, 2,3 , Cijkl the fourth-rank anisotropic elasticity tensor. The equilibrium
equations can be rewritten in terms of displacements by substituting Eqs. (180) and (181) into
Eq. (179) as
Cijkl uk , jl  0
The boundary conditions of the boundary value problem (180)-(182) are
on u
ui  ui
ti   ij n j  ti
on t
(182)
(183)
(184)
where ui and ti are the prescribed boundary displacement and traction vector, respectively. In
addition, ni is the unit outward normal to the boundary and =u+ t is the boundary of the
solution domain .
For the generalized two-dimensional deformation of anisotropic elasticity ui is
assumed to depend on x1 and x2 only. Based on this assumption, the general solution to
(182) can be written as [52, 53]
u  2Re{Af ( z)}, φ  2Re{Bf ( z)}
33
(185)
where u =  u1 , u2 , u3  is the displacement vector, φ  1 , 2 , 3  the stress function vector,
T
T
f ( z)  [ f1 ( z1 ), f 2 ( z2 ), f3 ( z3 )]T a function vector composed of three holomorphic complex
function f ( z ) ,   1, 2,3 , which is an arbitrary function with argument z  x1  p x2 and
will be determined by satisfying the boundary and loading conditions of a given problem. In
Eq. (185), Re stands for the real part of a complex number, p are the material eigenvalues
with positive imaginary part, A = a1 ,a2 ,a3  and B = b1 ,b2 ,b3  are 3×3 complex matrices
formed by the material eigenvector associated with p , which can be obtained by the
following eigen-relations [52]
Nξ  pξ
(186)
where N is a 6×6 foundational elasticity matrix and ξ is a 6×1 column vector defined by
a 
 N N2 
, ξ 
(187)
N 1
T
b 
 N 3 N1 
where N1 = -T-1RT , N2 = T-1 , N3 = RT-1RT - Q and the matrices Q , R and T are 3×3 matrices
extracted from Cijkl as follow
Qik  Ci1k1 , Rik  Ci1k 2 , Tik  Ci 2k 2
The stresses can be obtained from the derivation of stress functions φ as follow
 i1  2 Re{Lf ( z)},  i 2   2 Re{Bf ( z)}
where
L    p1b1 ,  p2b2 ,  p3b3 ,  p4b4 
(188)
(189)
(190)
6.1.2. Foundational solutions
To find the fundamental solution needed in our analysis, we have to first derive the Green’s
function of the problem: an infinite homogeneous anisotropic elastic medium loaded by a
concentrated point force (or line force for two-dimensional problems) pˆ  ( pˆ1 , pˆ 2 , pˆ 3 ) applied
at an internal point xˆ  ( xˆ1 , xˆ2 ) far from the boundary. The boundary conditions of this
problem can be written as


C
C
d   pˆ
for any closed curve C enclosing xˆ
du  pˆ
for any closed curve C
lim  ij  0
x 
Thus, the Green’s function satisfying above boundary conditions is found to be [62]
34
(191)
f ( z) 
1
ln( z  zˆ ) AT pˆ
2 i
(192)
Therefore, fundamental solutions of the problem can be expressed as
1
u  Im A ln  z  zˆ  A T pˆ

1
  Im B ln  z  zˆ  A T pˆ





The corresponding stress components can be obtained from stress function  as
1
 i*1  ,2   Im B p /  z  zˆ  A T pˆ

1
 i*2  ,1  Im B 1/  z  zˆ  A T pˆ





(193)
(194)
where pˆ are chosen to be (1,0,0)T ,(0,1,0)T ,(0,0,1)T , respectively,    stands for the
diagonal matrix corresponding to subscript  , Im denotes the imagery part of a complex
number, and superscript T denotes the matrix transpose.
6.1.3. Coordinate transformation
A typical composite laminate consists of individual layers, which are usually made of
unidirectional plies with the same or regularly alternating orientation. A layer is generally
referred to the global coordinate frame x1 , x2 , and x3 of the structural element rather than to
coordinates x, y, and z associated with the ply orientation. So it is necessary to transform the
constitutive relationship of each layer from the material coordinate frame to the uniform
global coordinate frame.
35
Figure 7
Schematic of the relationship between global coordinate system ( x1 , x2 ) and
local material coordinate system ( x, y ).
For the two coordinate systems mentioned above (see Figure 7), the angle between the axis-
x1 and axis-x is denoted by  , which is positive along the anti-clockwise direction, the
relationship for transformation of stress and strain components between the local material
coordinates and global coordinates is given by
11, 22 , 23 , 31,12 
T
and
 T1  xx ,  yy ,  yz ,  zx ,  xy 
T
(195)
 xx ,  yy ,  yz ,  zx ,  xy    T1  11 ,  22 ,  23 ,  31 , 12 
T
T
T
(196)
where the transformation matrix T and its inverse matrix are defined as
 c2 s2 0 0
c 2 s 2 0 0 -2cs 
2cs 
 2

 2

c2 0 0
2cs 
c2 0 0
2cs 
s
s
1
T 0
0 c s
0  T 0 0 c s
0 




0 s c
0 
0 0 -s c
0 
 0

,
 cs cs 0 0 c 2  s 2 
 cs -cs 0 0 c 2  s 2 




(197)
with c  cos( ), s  sin( ) . Therefore, the constitutive relationship in the global coordinate
system is given by
11, 22 , 23 , 31,12 
T
 T1C  T1  11 ,  22 ,  23 ,  31 , 12 
T
36
T
(198)
6.2. FORMULATIONS OF HFS-FEM
6.2.1. Assumed independent fields
The intra-element displacement fields for a particular element e is approximated in terms of a
linear combination of fundamental solutions of the problem as
 u1 (x) 


u(x)  u2 (x)   Nece (x  e ,y sj  e )
u ( x) 
 3 
(199)
where the matrix N e and unknown vector ce can be further written as
*
*
*
u11
(x, y s1 ) u12
(x, y s1 ) u13
(x, y s1 )
 *
*
*
Ne  u12 (x, y s1 ) u22
(x, y s1 ) u23
(x, y s1 )
 *
*
*
u13 (x, y s1 ) u32 (x, y s1 ) u33 ( x, y s1 )
ce  [c11 c21 c31
*
u13
( x, y sns ) 

*
*
*
u12
( x, y sns ) u22
( x, y sns ) u23
( x, y sns )  (200)

*
*
*
u13
( x, y sns ) u32
( x, y sns ) u33
(x, y sns ) 
c1n c2 n c3n ]T
(201)
*
*
u11
( x, y sns ) u12
( x, y sns )
in which ns is the number of source points, x and y sj are respectively the field point and
source point in the coordinate system (x1, x2) local to the element under consideration. The
fundamental solution uij* (x, y sj ) is given by Eq.(193) for general elements. γ is a
dimensionless coefficient for determining source points as described in previous sections.
The corresponding stress fields can be expressed as
σ(x)  11  22  23  31 12   Tece
T
(202)
where
*
 111
(x, y1 )
 *
 122 (x, y1 )
 *
(x, y1 )
Te   123
 * (x, y )
1
 131
*
 112 (x, y1 )
*
*
 211
(x, y1 )  311
(x, y1 )
*
*
*
 111
(x, y n )  211
(x, y n )  311
( x, y n ) 

 (x, y n ) 

(
x
,
y
)
ns 
s

*
*
*
*
*
 223
(x, y1 )  323
(x, y1 )
 123
( x, y ns )  223
( x, y ns )  323
(x, y ns ) 
*
*
*
*
*
 231
(x, y1 )  331
(x, y1 )
 131
(x, y ns )  231
(x, y ns )  331
( x, y ns ) 

*
*
*
*
*
 231
(x, y1 )  312
(x, y1 )
 112
(x, y ns )  212
(x, y ns )  312
( x, y ns ) 
(203)
*
222
(x, y1 ) 
*
322
s
*
122
(x, y1 )
s
*
222
(x, y ns ) 
s
*
322
*
The components  ijk
(x, y ) is given by Eq. (194) when pˆ i is selected to be (1, 0, 0)T , (0,1, 0)T
and (0, 0,1)T , respectively. As a consequence, the traction can be written as
t1 
 
t2   nσ  Qece
t 
 3
(204)
in which
Qe  nTe
37
(205)
 n1
n   0
 0
0
0
n3
n2
n3
0
0
n2
n1
n2 
n1 
0 
(206)
The unknown c e in Eq. (199) and Eq. (202) may be calculated using a hybrid technique [33],
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
 u1  N1 
   
(207)
u(x)  u2   N 2  d e  N ed e ,
(x   e )
u   
 3  N3 
where the symbol “~” is used to specify that the field is defined on the element boundary
only, N e is the matrix of shape functions, d e is the nodal displacements of elements. Taking
the side 3-4-5 of a particular 8-node quadrilateral element (see Figure 1) as an example, N e
and d e can be expressed as
Ne  0 0 N1 N2
N3 0 0 0
de = u11 u21 u31 u12 u22 u32
where the shape functions are expressed as
 Ni 0

N i   0 Ni
0
0

(208)
u18 u28 u38 
T
0
0 0 0 

0  , 0  0 0 0 
0 0 0 
Ni 
and N1 , N 2 and N 3 are expressed by natural coordinate  [1,1]
 1   
 1   
N1  
, N 2  1   2 , N3 
  1,1
2
2
(209)
(210)
(211)
6.2.2. Modified functional for HFS-FEM
With the assumption of two distinct intra-element field and frame field for elements, we can
establish the modified variational principle based on Eqs.(182)-(184) for the hybrid finite
element method of anisotropic materials [22, 28]. In the absence of the body forces, the
hybrid variational functional  me for a particular element e is constructed as
1
 me    ij ij d   ti ui d    ti (ui  ui )d 
t
e
2 e
where the boundary  e of the element e is
e  eu  et  eI
and
38
(212)
(213)
eu  e u , et  e t
(214)
and  eI is the inter-element boundary of element e. Performing a variation of  m , one
obtains
 me    ij ui , j d    ti ui d    [(ui  ui ) ti  ti ( ui   ui )]d 
e
et
e
(215)
Applying Gaussian theorem

f,i d    f  ni d 
e
(216)
e
and the definitions of traction force
ti   ij n j
(217)
we obtain
e    ij , j ui d    ti ui d    ti ui d 
e
e
(218)
et
Then, substituting Eq.(218) into Eq.(215) gives
 me    ij , j ui d    ti ui d    [(ui  ui ) ti  ti ui ]d 
e
et
e
(219)
Considering the fact that

t  ui d   0
(220)
eu i
we can finally obtain the following form
 me    ij , j ui d    (ti  ti ) ui d    (ui  ui ) ti d    ti ui d  (221)
e
et
e
I
Therefore, the Euler equations for Eq. (212) result in Eqs. (182)-(184) because the quantities
 ui ,  ti and  ui may be arbitrary. As for the continuity condition between elements, it can
be easily seen from the following variational of two adjacent elements such as e and f
 m ( e f )     ij , j ui d    (ti  ti ) ui d 
e  f


e  f
et et
(ui  ui ) ti d   
efI
tie  tif   ui
(222)
This indicates that the stationary condition of the functional satisfies both the required
boundary and inter-element continuity equations. In addition, the existence of extremum of
functional (212) can be easily proved by the so-called “second variational approach” as well,
which indicates functional (212) have a local extreme. Therefore, we can conclude that the
variational functional (212) can be used for deriving hybrid finite element formulations.
39
6.2.3. Element Stiffness Equation
Using Gaussian theorem and equilibrium equations, all domain integrals in Eq. (212) can be
converted into boundary integrals as follows
1
 me    ti ui d    ti ui d    tiui d 
e
t
2 e
(223)
Substituting Eqs. (199), (204) and (207) into the functional (223) yields the formulation as
1
 me   ceT Hece  ceT G ede  deT g e
(224)
2
where
He   QeT Ne d , G e   QeT Ne d , g e   NeT t d 
e
e
t
(225)
The stationary condition of the functional  me with respect to ce and de, respectively, yields,
Eqs. (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 [5, 28, 63-70]. In this section, the HFS-FEM is developed for modeling
two-dimensional piezoelectric material. 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 xi (i=1, 2, 3)
are given by
 ij , j  0,
Di,i  0 in 
40
(226)
where  ij is the stress tensor, Di the electric displacement vector,  the solution domain.
With strain and electric field as the independent variables, the constitutive equations are
written as
 ij  cijkl  kl  ekij Ek , Di  eikl  kl  ik Ek
(227)
where cijkl is the elasticity tensor measured under zero electric field, eikl and  ij are
respectively the piezoelectric tensor and dielectric tensor measured under zero strain.
 ij and Ei are the elastic strain tensor and the electric field vector, respectively. The relation
between the strain tensor  ij and the displacement ui is given by
1
2
 ij  (ui , j  u j ,i )
(228)
and the electric field component Ei is related to the electric potential  by
Ei  ,i
(229)
The boundary conditions of the boundary value problem (226)-(229) can be defined by:
on  u
(230)
ui  ui
ti   ij n j  ti
on  t
(231)
Dn  Di ni  qn  Dn on  D
 
on 
(232)
(233)
where ui , ti , qn and  are respectively the prescribed boundary displacement, the prescribed
traction vector, the prescribed surface charge and the prescribed electric potential. In
addition, ni is the unit outward normal to the boundary and
  u  t   D  
(234)
is the boundary of the solution domain .
For the transversely isotropic material, if x1 - x2 is taken as the isotropic plane, one can
employ either x1 - x3 or x2 - x3 plane to study the plane electromechanical phenomenon. Thus,
choosing the former and considering the plane strain conditions (  22   32  12  0 and E2  0
) Eq.(227) can be reduced to
 11   c11
  
 33   c13
   0
 13  
 D1   0
 
 D3  e31
c13
c33
0
0
e33
0   11   0


0    33    0
2  e
c44  
13 
 15
e31 
E 
e33   1 
E
0   3 
 
e15   11  11 0   E1 
  
 
0   33   0  33   E3 
213 
41
(235)
For the plane stress piezoelectric problem (  22   32  12  0 and D2  0 ), the constitutive
equations can be obtained by replacing the coefficients c11 , c13 , c33 , c44 , e15 , e31 , e33 , 11 ,  33 in
Eq. (235) as
*
c11
 c11  c122 / c11 , c13*  c13  c12c13 / c11 ,
*
*
c33
 c33  c132 / c11 , c44
 c44 , e15*  e15 ,
*
*
e31
 e31  c12e31 / c11 , e33
 e33  c13e31 / c11 ,
(236)
11*  11 ,  33*   33  e312 / c11
7.2. ASSUMED INDEPENDENT FIELDS
For the piezoelectric problems, HFS-FEM is based on assuming two distinct
displacement and electric potential (DEP) fields: intra-element DEP field u and an
independent DEP frame field u along element boundaries [16, 22].
7.2.1. Intra-element field
The intra-element DEP field u identically fulfills the governing differential equations
(226) and is approximated by a linear combination of foundational solutions at different
source points located outside of the element domain
*
*
*
 u1  ns u11 (x, y sj ) u21 (x, y sj ) u31 (x, y sj )   c1 j 
 *
 
 
*
*
ue  u2    u12
(x, y sj ) u22
(x, y sj ) u32
(x, y sj )  c2 j   Nece (x  e , y sj  e ) (237)
   j 1 u* (x, y ) u* (x, y ) u* (x, y )  c 
sj
23
sj
33
sj   3 j 
 
 13
where the fundamental solution matrix N e is now given by
*
*
*
u11
(x, y s1 ) u21
(x, y s1 ) u31
(x, y s1 )
 *
*
*
Ne  u12 (x, y s1 ) u22 (x, y s1 ) u32 (x, y s1 )
 *
*
*
u13 (x, y s1 ) u23 (x, y s1 ) u33 (x, y s1 )
ce  [c11 c21 c31
*
*
*
... u33
( x, y sns ) u33
( x, y sns ) u33
( x, y sns ) 

*
*
*
... u33
(x, y sns ) u33
( x, y sns ) u33
( x, y sns )  (238)

*
*
*
... u33
( x, y sns ) u33
( x, y sns ) u33
(x, y sns ) 
(239)
... c1ns c2 ns c3ns ]T
in which x and y sj are respectively the field point (i.e. the nodal points of the element in this
work) and source point in the local coordinate system (x1, x2). The component uij* (x, y sj ) is
the induced displacement component (i=1, 2) or electric potential (i=3) in i-direction at the
field point x due to a unit point load (j=1, 2) or point charge (j=3) applied in j-direction at the
source point y sj . The fundamental solution uij* (x, y sj ) is given as [5, 65, 71]
42
 
1 3
s j1t12 j 1 ln rj
u11 

 M 11 j 1


x1  xs1
1 3

s j 2t12 j 2 arctg
u12 

 M 12 j 1
s j  x3  xs 3 


x1  xs1
1 3
u13 
s j 3t12 j 3arctg

 M 13 j 1
s j  x3  xs 3 

x1  xs1
 
1 3
1
u

 21  M  d j1t 2 j 1arctg s x  x
11 j 1
j  3
s3 

3

1
 
d j 2t12 j 2 ln rj
u22 

 M 12 j 1

 
1 3
d j 3t12 j 3 ln rj
u23 

 M 13 j 1


(240)
(241)
x1  xs1
 
1 3
1
u

 31  M  g j1t 2 j 1arctg s x  x
11 j 1
j  3
s3 


1 3
 
u

g j 2t12 j 2 ln rj
 32


M
12 j 1

 
1 3
g j 3t12 j 3 ln rj
u33 

 M 13 j 1


where rj 
 x1  xs1 
2
(242)
 s 2j  x3  xs 3  , and s j is the three different roots of the characteristic
2
equation asi6  bsi4  csi3  d  0 . The source point y sj ( j  1, 2,
following method [2]
, ns ) can be generated by the
y s  x 0   (x 0  x c )
(243)
Making use of Eq.(228) and the expression of intra-element DEP field u in Eq.(237),
the corresponding stress and electric displacement in Eq. (235) can be written as
σ  Tece
where σ  11  22 12
 11* (x, y s1 )
 *
 21 (x, y s1 )
 *
(x, y s1 )
Te   31
*
 (x, y )
s1
 41
*
 51 (x, y s1 )
D1
(244)
D2  and
T
 12* (x, y s1 )  13* (x, y s1 ) ...  11* (x, y sn )  12* (x, y sn )  13* (x, y sn ) 
s
s
s

*
*
*
*
*
 22
(x, y s1 )  23
(x, y s1 ) ...  21
(x, y sn )  22
(x, y sn )  23
( x, y sn ) 

 32* (x, y s1 )  33* (x, y s1 ) ...  31* (x, y sn )  32* ( x, y sn )  33* (x, y sn ) 
s
s
s
s
s
s
 (x, y s1 )  (x, y s1 ) ...  (x, y sn )  (x, y sn )  (x, y sn ) 
*
42
*
43
*
41
s
*
42
s
*
43
s

 (x, y s1 )  (x, y s1 ) ...  (x, y sn )  (x, y sn )  (x, y sn ) 
*
52
*
53
*
51
s
*
52
s
*
53
s
(245)
in which  (x, y sj ) denotes the corresponding stress components (i=1,2,3) or electric
*
ij
displacement (i=4,5) along i-direction at the field point x due to a unit point load (j=1,2) or a
43
unit point charge (j=3) applied in j-direction at the source point ys and can be derived from
Eqs.(240)-(242), and are listed below, which are derived by substituting the fundamental
solutions into constitutive equations [72].
 
x1  xs1 
1 3 
1
 11 
 c11s j1  c13d j1s j  e31 g j1s j  t 2 j 1


 M 11 j 1 
rj2 


s j  x3  xs 3  
1 3 
 
1


 c11s j 2  c13d j 2 s j  e31 g j 2 s j  t 2 j 2

 12

 M 12 j 1 
rj2



3 
s j  x3  xs 3  
   1
c11s j 3  c13d j 3 s j  e31 g j 3 s j  t12 j 3




13

 M 13 j 1 
rj2


 
x1  xs1 
1 3 
1
 21 
 c13 s j1  c33 d j1s j  e33 g j1s j  t 2 j 1


 M 11 j 1 
rj2 


s j  x3  xs 3  
1 3 
 
1
 c13 s j 2  c33 d j 2 s j  e33 g j 2 s j  t 2 j 2

 22 

 M 12 j 1 
rj2



3 
s j  x3  xs 3  
    1
c13 s j 3  c33d j 3 s j  e33 g j 3 s j  t12 j 3




23
2


M
r
j 1 

13
j


 
s j  x3  xs 3  
1 3 
1
 31 
 c44 s j1s j  c44 d j1  e15 g j1  t 2 j 1


 M 11 j 1 
rj2



x1  xs1 
1 3 
 
1


 c44 s j 2 s j  c44 d j 2  e15 g j 2  t 2 j 2

 32

 M 12 j 1 
rj2 


3 
x x 
    1
c44 s j 3 s j  c44 d j 3  e15 g j 3  t12 j 3 1 2 s1 



33

 M 13 j 1 
rj 

 
s j  x3  xs 3  
1 3 
1
 41 
 e15 s j1s j  e15 d j1  11 g j1  t 2 j 1


 M 11 j 1 
rj2



x1  xs1 
1 3 
 
1
 e15 s j 2 s j  e15 d j 2  11 g j 2  t 2 j 2

 42 

 M 12 j 1 
rj2 


3 
x x 
    1
e15 s j 3 s j  e15 d j 3  11 g j 3  t12 j 3 1 2 s1 



43

 M 13 j 1 
rj 

44
(246)
(247)
(248)
(249)
 
x1  xs1 
1 3 
1
 51 
 e31s j1  e33d j1s j  33 g j1s j  t 2 j 1


 M 11 j 1 
rj2 


s j  x3  xs 3  
1 3 
 
1


 e31s j 2  e33 d j 2 s j  33 g j 2 s j  t 2 j 2

 52

 M 12 j 1 
rj2



3 
s j  x3  xs 3  
    1
e31s j 3  e33d j 3 s j  33 g j 3 s j  t12 j 3




53

 M 13 j 1 
rj2


(250)
in which the coefficients sij , dij , gij , tij , M11 , M12 , M13 are defined as in literature [71].
From Eqs.(227), (231) and (232), the generalized traction forces and electric
displacement are given as
 t1   Q1 
   
 t2   Q 2  ce  Qece
 D  Q 
 n  3
(251)
Qe  nTe
(252)
where
 n1
n   0
 0
0
n2
0
n2
n1
0
0
0
n1
0
0 
n2 
(253)
7.2.2. Auxiliary frame field
For the two-dimensional piezoelectric problem under consideration, the frame field is
assumed as
 u1   N1 
   
u(x)  u2   N 2  d e  N ed e ,
   N 
   3
(x   e )
(254)
where N e is a matrix of the corresponding shape functions. For the side 3-4-5 of a particular
quadratic element as shown in Figure 1, the shape function matrix N e and nodal vector d e
can be given in the form

0 0 N1
N e  0 0 0

0 0 0
 6
0
0
N2
0
0
N3
0
0
0
N1
0
0
N2
0
0
N3
0
0
0
N1
0
0
N2
0
0
N3
0
de  u11 u21 1
9
u14 u24 4
u18 u28 8 
where the shape functions N e are expressed by natural coordinate 
45
T

0
0

0

(255)
(256)
N1  
 1   
2
, N 2  1   2 , N3 
 1   
2
  1,1
(257)
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 inter-element continuity
requirements in addition to Eqs.(226)-(233):
uie  uif e   f
(on e   f , conformity)
tie  tif  0
Dne  Dnf  0 (on e   f , reciprocity)
(258)
(259)
where ‘ e ’ and ‘ f ’ stand for any two neighboring elements. Eqs.(226)-(233) together with
Eq.(258) and (259) can now be taken as the basis to establish the modified variational
principle for the hybrid finite element method of piezoelectric materials [16, 22].
Since the stationary conditions of the traditional potential or complementary
variational functional cannot satisfy the inter-element continuity condition required in the
proposed HFS-FEM, new modified variational functional should be developed. In the
absence of the body forces and electric charge density, the hybrid functional  me for a
particular element e is constructed as
 me   e   ti (ui  ui )d    Dn (   )d 
e
e
(260)
where
e 
1
( ij ij  Di Ei )d    tiui d    Dn d 
t
D
2 e
(261)
and the boundary  e of the element e is
e  eu et eI  e eD eI
(262)
eu  e u , et  e t , e  e  , eD  e  D
(263)
and
and  eI is the inter-element boundary of element e . Compared to the functional employed in
the conventional FEM, the present hybrid functional is constructed by adding two integral
terms related to the intra-element and element frame DEP fields to guarantee the satisfaction
of displacement and electrical potential continuity condition on the common boundary of two
adjacent elements.
46
It can be proved that the stationary conditions of the above functional (260) leads to
Eqs.(226)-(233). To this end, performing a variation of  m , one obtains
 me  e   [(ui  ui ) ti  ti ( ui   ui )]d    [(   ) Dn  Dn (   )]d  (264)
e
e
in which the first term is
 e    ij ij d )   Di Ei d    ti ui d   
Dn d 
   ij ui , j d    Di,i d    ti ui d   
Dn d 
e
e
e
et
e
eD
et
eD
(265)
Applying Gaussian theorem

e
f,i d    f  ni d 
(266)
e
and the definitions of traction force and electrical displacement
ti   ij n j , Dn  Di ni
(267)
we obtain
 e     ij , j ui d    Di ,i d    ti ui d    Dn d 
e
e
  ti ui d   
et
eD
e
e
Then, substituting Eq. (268) into Eq. (264) gives
 me     ij , j ui d    Di ,i d    ti ui d   
e
(268)
Dn d 
e
et
eD
Dn d 
(269)
  [(ui  ui ) ti  ti ui ]d    [(   ) Dn  Dn ]d 
e
e
Considering the fact that

t  ui d   0,
eu i

e
Dn d   0
(270)
we finally get the following form
 me     ij , j ui d    Di ,i d    (ti  ti ) ui d   
e
e
et
eD
( Dn  Dn ) d 
  (ui  ui ) ti d    (   ) Dn d    ti ui d    Dn d 
e
e
I
(271)
I
Therefore, the Euler equations for Eq. (271) result in Eqs. (226)-(233) and Eq. (258) because
the quantities  ui ,  ti ,  ,  Dn ,  ui , and  may be arbitrary. As for the continuity
condition of Eq.(259), it can easily be seen from the following variational of two adjacent
elements such as e and f

 m ( e f )  
e  f


e  f
 ij , j ui d  

Di ,i d  
e  f
(ui  ui ) ti d  

e  f

(ti  ti ) ui d  
et et
(   ) Dn d   
efI
47

( Dn  Dn ) d 
eD  fD
t
ie
 tif   ui  
efI
D
ne
 Dnf   d 
(272)
which indicates that the stationary condition of the functional satisfies both the required
boundary and inter-element continuity equations. In addition, the existence of extremum of
functional (260) can be easily proved by the “second variational approach” as well, which
indicates functional (260) have a local extreme.
7.3.2. Element Stiffness Equation
The element stiffness equation can be generated by setting  me  0 . To simplify the
derivation, we first transform all domain integrals in Eq. (260) into boundary ones. With
Gaussian theorem, the functional in Eq. (260) may be simplified as
1
 me   ( ij  ij  Di Ei )d    ti ui d    Dn d    ti (ui  ui )d    Dn (   )d 
t
D
e
e
2 e
1
1
 (  ti ui d     ij , j ui d )  (  Dn d    Di ,i d )   tiui d    Dn d 


e
t
D
e
2 e
2 e
  ti (ui  ui )d    Dn (   )d 
e
e
(273)
Due to the satisfaction of the equilibrium equation with the constructed intra-element fields,
we have the following expression for the HT finite element model
1
1
 me   ti ui d    Dn d    ti ui d    Dn d    ti (ui  ui )d    Dn (   )d 

t
D
e
e
2 e
2 e
1
   (ti ui  Dn )d    (tiui  Dn )d    tiui d    Dn d 
e
t
D
2 e
(274)
Substituting Eqs. (237), (254) and (251) into the above functional (274) yields the
formulation as
1
 me   ceT Hece  ceT G ede  deT g e
2
(275)
where
He   QeT Ne d
e
,
G e   QeT Ne d
e
,
ge   NeT t d    NeT Dd 
t
D
.
The stationary condition of the functional  me with respect to ce and de, respectively, yields
Eqs. (33) and (34). Following the procedure described in [21, 22], the missing rigid-body
motion can be recovered by setting the augmented internal field of a particular element e as
1 0 x2 0 
(276)
u e  N ece  0 1  x1 0  c0
0 0 0 1 
48
where the undetermined rigid-body motion parameter c 0 can be calculated using the least
square matching of u e and u e at element nodes

n

2
2
2
min    u1i  u1i    u2i  u2i   i  i 


i 1 
(277)
which finally gives
c0 = R e1re
0
x2i
1

n
0
1
 x1i
Re =  
2
 x1i x1i  x22i
i 1  x2 i

0
0
0
ue1i




n
ue 2i


re  



u
x


u
x
i 1
e1i 2 i
e 2 i 1i


ei


(278)
0
0 
0

1
(279)
(280)
in which Δuei =  ue - uˆ e  node i and n is the number of element nodes. As a consequence, c 0
can be calculated by Eq. (278) once the nodal DEP fields d e and the interpolation
coefficients c e are respectively determined by Eqs. (237) and (254). Then the complete DEP
fields u e can be obtained from Eq. (276).
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 [62, 73]. To resolve this problem, normalization of each
quantity by its reference value is employed. The reference values for the stiffness,
piezoelectric stress constant, dielectric constants and strain are selected to be
c0  1011  N / m2  , e0  101 ( N / mV ) , k0  109  C / mV  ,  0  103 V / m  , respectively.
The reference values of other quantities, as shown in Table 1, are determined in terms of
these four fundamental reference variables and the characteristic length x0  100  m  of the
problem so that the normalized governing equations remain in exactly the same form as the
original equations.
49
Table 1
Reference values for material constants and field variables in piezoelectricity
derived from basic reference variables: c0 , e0 , k0 ,  0 and x0 .
Displacement
u0  x0 0  103 (m)
Electric Potential
0  x0 E0  107 (V )
Stress
 0  c0 0  108 ( N / m2 )

s0  0  1011 (m2 / N )
0
0
7
Electric induction
D0  k0 E0  102 (C / m2 )
E
0  0  109 (mV / C )
D0
E
g0  0  101 (mV / N )
0
Compliance
Electric field
E0 
e0
 10 (V / m)
Impermeability
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 10Mpa and
uniform distributed tension of 100MPa are applied to the cube. Figure 10 present the
displacement component u1 and the stress component  11 at Point A of the block, which are
calculated by the HFS-FEM on the three meshes with distorted 8-node brick elements (Figure
9). 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 10 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 compare with EAS element in traditional FEM. Contour plots
of u1 and  11 obtained by HFS-FEM on Mesh 3 are also presented in Figure 11. It should be
50
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 [37, 74].
W  10
A
P  100
H  10
x3
L  10
x1
x2
Figure 8
Cubic block under uniform tension and body force: geometry, boundary
condition and loading.
Figure 9
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).
51
Figure 10
Cubic block with body force under uniform distributed load: Convergent study
of (a) displacement u1 and (b) stress  11 .
52
Figure 11
Contour plots of (a) displacement u1 and (b) stress  11 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 T changes logarithmically along the radial direction. With the
symmetry condition of the problem, only one quarter of the cylinder is modeled. The
configurations of geometry and the boundary conditions are shown in Figure 12. In our
computation, the parameters a=5, b=20, E=1000, ν =0.3, α =0.001, and T0 =10. The two
approaches listed in Section 4 to approximate the body force and temperature are discussed
and analyzed in this example.
53
Figure 12
(a) Geometry and boundary conditions of the long cylinder with axisymmetric
temperature change, (b) Mesh configurations of a quarter of circular cylinder
(128 eight-node elements)
Figure 13 present the variation of the radial and the circumferential thermal stresses
with the cylinder radius, respectively, in which the theoretical values are given for
comparison [35]. 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.
54
Figure 13
(a) Radial thermal stresses and (b)circumferential thermal stresses with the
cylinder radius.
55
Figure 14
Contour plots of (a) radial and (b) circumferential thermal stresses (The
background mesh in the figures (a) and (b) is used for plots the calculated results
in post processing)
Figure 14 displays the contour plots of (a) radial and (b) circumferential thermal
stresses (The meshes used for contour plot is 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 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.
8.3. 3D CUBE UNDER ARBITRARY TEMPERATURE AND BODY FORCE
As shown in Figure 15, a 3D cube of 111 with center located at (0.5, 0.5, 0.5) is
considered in this example. The material properties of the cube are Young’s modulus
E  5000 , Poisson’s ratio v  0.3 , and linear thermal expansion coefficient   0.001 . The
56
bottom surface is fixed on the ground and the temperature distribution and body force are
assumed to be
T  30 x3 , b1  0, b2  0, b3  2000 ( x  0.5)2  ( y  0.5)2 
Because there is no analytical solution available, the results from ABAQUS herein are
employed for comparison. The meshes used by HFS-FEM and ABAQUS are given in Figure
15.
Figure 15
(a) Schematic of the 3D cube under arbitrary temperature and body force, (b)
mesh used by HFS-FEM (125 20-node brick elements), (c) mesh for ABAQUS
(8000 C3D20R elements).
Figure 16 presents the displacement u3 and stress  33 along one edge of the cube
which is coinciding with x3 axis. It can be seen that the results from HFS-FEM are 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
57
body force. It is also shown that HFS-FEM with RBF interpolation can give satisfactory
results using coarse meshes.
Figure 16
(a) Displacement u3 and (b) stress  33 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  0 =1 GPa is applied in x2 direction. The material parameters
of the orthotropic plate are taken as El=113 GPa, E2=52.7 GPa, G12=28.5 GPa, v12= 0.45.
The length and width of the plate are L=20 mm and W=20 mm; the major and minor lengths
of the ellipse a and b are respectively 2 mm and 1 mm. In the computation, plane stress
condition is used. The mesh configurations used by HFS-FEM and ABAQUS are given in
Figure 17.
58
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.
7
0
=0 (ABAQUS)
6
0
=30 (ABAQUS)
Hoop stress  (GPa)
5
0
=45 (ABAQUS)
0
4
=60 (ABAQUS)
3
=0 (HFS-FEM)
0
=90 (ABAQUS)
0
0
=30 (HFS-FEM)
2
0
=45 (HFS-FEM)
0
=60 (HFS-FEM)
1
0
=90 (HFS-FEM)
0
-1
-2
0
10
20
30
40
50
60
70
80
90
Angle  (Degree)
Figure 18
Variation of hoop stresses along the rim of the elliptical hole for different fiber
orientation  .
59
Figure 19
Contour plots of stress components around the elliptic hole in the composite
plate: (a)   00 , (b)   450 , (c)   900 .
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  11 ,  22 ,  12 around the elliptic hole in the
composite plate for several different fiber angles are shown in Figure 19.
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.
60
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   900 whereas the smallest appears at
  0o . It indicates the effectiveness of the proposed method in predicting the SCF for
anisotropic composites as well.
7.0
ABAQUS
HFS-FEM
6.5
SCF (/0)
6.0
5.5
5.0
4.5
4.0
0
10
20
30
40
50
60
70
80
90
Angle  (Degree)
Figure 20
Variation of SCF with the lamina angle  .
8.5. ISOTROPIC PLATE WITH MULTI-ANISOTROPIC 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 multi-anisotropic inclusions of square
geometry (Edge length a=2) are considered. The distance between any two inclusions is
assumed to be b=3. The material parameters for the inclusions are chosen as El=134.45 GPa,
E2= E3=11.03 GPa, G23=2.98 GPa, G31= G12=28.5 GPa, v23=0.49, v31= v12=0.301. The
material parameters for isotropic matrix are elastic modules E=2.8 Gpa and poison’s ratio
v=0.3. The mesh configuration of the plate for HFS-FEM is shown in Figure 21, which uses
272 quadratic general elements.
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 p  i,(  1, 2,3) (Lekhnitskii, 1968). However, a small perturbation
of the material constants, such as p1= (1-0.004)i, p2=i, p3= (1+0.004)i, can be applied to
make the eigenvalues be distinct and the results can be applied conveniently.
61
Figure 21
(a) Schematic and (b) mesh configuration of an isotropic plate with multianisotropic inclusions.
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. x2=0) and are 0.7% and 1.3%,
respectively. Additionally, it can be found that the results from HFS-FEM are better than
those from ABAQUS using the same mesh. Although the stress at the vicinity of the
inclusions (Point C) has a little degradation due to the high stress concentration and stress
contrast at the adjacent elements, the displacement still agrees well with the reference
solution. The variations of displacement components u1 and u2 along the right edge (x=8) by
HFS-FEM is shown in Figure 22.
62
Table 2
Comparison of displacement and stress at points A and B.
Items
Disp. u1
Stress  11
Figure 22
Points
HFS-FEM
ABAQUS
ABAQUS
(272 elements) (272 elements) (30471elements )
A
0.04322
0.04318
0.04335
B
0.03719
0.03721
0.03744
C
0.03062
0.03076
0.03091
A
10.0446
9.9219
9.9992
B
9.8585
9.8304
9.9976
C
20.6453
13.7302
23.6625
The variation of displacement component (a) u1 and (b) u2 along the right edge
of the plate (x=8) by HFS-FEM and ABAQUS.
63
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  xx   0  10
parallel to the x axis be imposed at infinity with traction- and electric charge-free at the
boundary of the hole. In our calculation, the radius of the hole is set to be r = 1 and L/r =20 is
employed in the analysis.
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.
Table 3
Properties of the material PZT-4 used for the model
Parameters
Values
Parameters
Values
c11
13.9 1010 Nm2
e15
13.44Cm2
64
c12
7.78 1010 Nm2
e31
6.98Cm2
c13
7.43 1010 Nm2
e33
13.84Cm2
c33
11.3 1010 Nm2
11
6.0 109 C / Nm
c44
2.56 1010 Nm2
 33
5.47 109 C / Nm
Figure 24 shows the distribution of hoop stress   and radial stress  r along the line
z=0 for remote loading  xx and along the line x=0 for remote loading  zz , respectively.
Figure 25 presents the variations of the normalized stress   /  0 and the normalized electric
displacement D /  0 1010 along the hole edge under remote mechanical loading. It is
obvious that the results obtained from HFS-FEM agree well with the results from ABAQUS
and Sosa [63].
Figure 24
Distribution of (a) hoop stress   and (b) radial stress  r along the line z=0
when subjected to remote mechanical load  xx and along the line x=0 when
subjected to remote mechanical load  zz .
65
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  0 when r increasing
toward infinity. Compared with the hoop stress   , it is obvious that the radial stress  r is
much smaller and usually does not need to be considered. It is obvious that loading along the
poling direction will produce smaller stress concentration due to coupling effect. The
maximum values of   appears at   90 for case of  xx and at   0 and   180 for case
of loading  zz , both of which agree well with the analytical solution from Sosa [63]. The
electric displacement D /  0 1010 produced by  xx and  zz are nearly the same, and are
symmetrical with respect to the x axis. It is found that the maximum values of D appears at
  65 and   115 , which also agrees well with the analytical solution.
Figure 25
Variation of (a) normalized stress   /  0 and (b) electrical displacement
D /  0 1010 along the hole boundary under remote mechanical loading
66
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 of arbitrary polygonal or even curve-sided
elements, 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 intra-element
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, and 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). The special-purpose functions warrant 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 multi-field problems such as thermo-elasticpiezoelectric materials, thermo-magnetic-electric-mechanical materials, and to develop
multiscale framework across from continuum to micro- and nano-scales for modeling
heterogeneous materials.
67
References
[1] 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(6): pp.
493-510, 2010.
[2] 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(5): pp.
801-809, 2010.
[3] C. Cao, Q.H. Qin, and A. Yu, "A novel Boundary-Integral based fnite element method
for 2D and 3D thermo-elasticity problems," Journal of Thermal Stresses., vol. 35(10): pp.
849-876, 2012.
[4] C. Cao, Q.H. Qin, and A. Yu, "New hybrid finite element approach for three-dimensional
elastic problems," Archive of Mechanics, vol. 64(3): pp. 261-292, 2012.
[5] C. Cao, Q.H. Qin, and A. Yu, "Hybrid fundamental-solution-based FEM for piezoelectric
materials," Computational Mechanics, vol. 50(4): pp. 397-412, 2012.
[6] C. Cao, A. Yu, and Q.H. Qin, "Evaluation of effective thermal conductivity of fiberreinforced composites," International Journal of Architecture, Engineering and
Construction, vol. 1(1): pp. 14-29, 2012.
[7] H. Wang and Q.H. Qin, "A new special element for stress concentration analysis of a
plate with elliptical holes," Acta Mechanica, vol. 223(6): pp. 1323-1340, 2012.
[8] 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(1): pp. 39-46, 2012.
[9] H. Wang and Q.H. Qin, "Boundary integral based graded element for elastic analysis of
2D functionally graded plates," European Journal of Mechanics-A/Solids, vol. 33: pp. 1223, 2012.
[10] 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.
[11] 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(04): pp. 661-674, 2013.
[12] C. Cao, A. Yu, and Q.-H. Qin, "A new hybrid finite element approach for plane
piezoelectricity with defects," Acta Mechanica, vol. 224(1): pp. 41-61, 2013.
[13] Q.H. Qin, "Hybrid Trefftz finite-element approach for plate bending on an elastic
foundation," Applied Mathematical Modelling, vol. 18(6): pp. 334-339, 1994.
[14] Q.-H. Qin, "Trefftz finite element method and its applications," Applied Mechanics
Reviews, vol. 58(5): pp. 316-337, 2005.
68
[15] K. Wang, Q. Qin, Y. Kang, J. Wang, and C. Qu, "A direct constraint‐Trefftz FEM
for analysing elastic contact problems," International journal for numerical methods in
engineering, vol. 63(12): pp. 1694-1718, 2005.
[16] Q.H. Qin, "Variational formulations for TFEM of piezoelectricity," International
Journal of Solids and Structures, vol. 40(23): pp. 6335-6346, 2003.
[17] H. Wang and Q.H. Qin, "Special fiber elements for thermal analysis of fiberreinforced composites," Engineering Computations, vol. 28(8): pp. 1079-1097, 2011.
[18] Q.H. Qin, "Nonlinear analysis of reissner plates on an elastic-foundation by the
BEM," International Journal of Solids and Structures, vol. 30(22): pp. 3101-3111, 1993.
[19] Q.H. Qin and Y.Y. Huang, "BEM of postbuckling analysis of thin plates," Applied
Mathematical Modelling, vol. 14(10): pp. 544-548, 1990.
[20] J.A.T. de Freitas and M. Toma, "Hybrid-Trefftz stress elements for incompressible
biphasic media," International journal for numerical methods in engineering, vol. 79(2):
pp. 205-238, 2009.
[21] Q.H. Qin and H. Wang, Matlab and C programming for Trefftz finite element
methods. New York: CRC Press, 2008.
[22] Q.H. Qin, The Trefftz finite and boundary element method. Southampton WIT Press,
2000.
[23] J. Jirousek and Q.H. Qin, "Application of hybrid Trefftz element approach to transient
heat-conduction analysis," Computers & Structures, vol. 58(1): pp. 195-201, 1996.
[24] 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(30): pp. 4583-4604, 1996.
[25] J. Jirousek, A. Wroblewski, 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(1): pp. 315-344, 1995.
[26] 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(6): pp. 1262-1277, 2007.
[27] 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-English Edition, vol. 26(12): pp. 1564-1573, 2005.
[28] Q.H. Qin, "Solving anti-plane problems of piezoelectric materials by the Trefftz finite
element approach," Computational Mechanics, vol. 31(6): pp. 461-468, 2003.
[29] 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(14-15): pp. 1314-1323,
2006.
69
[30] H. Wang, Q.H. Qin, and Y. Kang, "A meshless model for transient heat conduction in
functionally graded materials," Computational Mechanics, vol. 38(1): pp. 51-60, 2006.
[31] H. Wang and Q.H. Qin, "Some problems with the method of fundamental solution
using radial basis functions," Acta Mechanica Solida Sinica, vol. 20(1): pp. 21-29, 2007.
[32] H. Wang, Q.H. Qin, and Y. Kang, "A new meshless method for steady-state heat
conduction problems in anisotropic and inhomogeneous media," Archive of Applied
Mechanics, vol. 74(8): pp. 563-579, 2005.
[33] H. Wang and Q.H. Qin, "Hybrid FEM with fundamental solutions as trial functions
for heat conduction simulation," Acta Mechanica Solida Sinica, vol. 22(5): pp. 487-498,
2009.
[34] 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(1): pp. 1-30, 1987.
[35] S.P. Timoshenko and J.N. Goodier, Theory of Elasticity (second edition). New York:
McGraw-Hill, 1951.
[36]
S.A. Sauter and C. Schwab, Boundary Element Methods. Berlin: Springer, 2010.
[37] 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(4-5): pp. 377-387, 2001.
[38] 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(9): pp. 2061-2077, 1988.
[39] 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(1-4): pp. 41-50, 1999.
[40] R. Carrazedo and H.B. Coda, "Alternative positional FEM applied to
thermomechanical impact of truss structures," Finite Elements in Analysis and Design,
vol. 46(11): pp. 1008-1016, 2010.
[41] R.J. Yang, "Shape design sensitivity analysis of thermoelasticity problems,"
Computer Methods in Applied Mechanics and Engineering, vol. 102(1): pp. 41-60, 1993.
[42] A. Chaudouet, "Three-dimensional transient thermo-elastic analyses by the BIE
method," International journal for numerical methods in engineering, vol. 24(1): pp. 2545, 1987.
[43] 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(78): pp. 1298-1307, 2010.
[44] Z.W. Zhang, H. Wang, and Q.H. Qin, "Transient Bioheat Simulation of the LaserTissue Interaction in Human Skin Using Hybrid Finite Element Formulation," MCB:
Molecular & Cellular Biomechanics, vol. 9(1): pp. 31-54, 2012.
70
[45] 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(9): pp. 2447-2460, 2004.
[46] Q.H. Qin, "Thermoelectroelastic analysis of cracks in piezoelectric half-plane by
BEM," Computational mechanics, vol. 23(4): pp. 353-360, 1999.
[47] C. Tsai, "The method of fundamental solutions with dual reciprocity for threedimensional thermoelasticity under arbitrary body forces," Engineering Computations,
vol. 26(3): pp. 229-244, 2009.
[48] 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(7-8): pp. 531-538, 2000.
[49] V.V. Vasiliev and E.V. Morozov, Advanced mechanics of composite materials.
Elsevier Science, 2007.
[50] S.G. Lekhnitskii, Theory of Elasticity of an Anisotropic Body. Moscow: Mir
Publishers, 1981.
[51] Q.H. Qin, "A new solution for thermopiezoelectric solid with an insulated elliptic
hole," Acta Mechanica Sinica, vol. 14(2): pp. 157-170, 1998.
[52] T.C.T. Ting, Anisotropic Elasticity: Theory and Applications NewYork: Oxford
Science Publications, 1996.
[53] A.N. Stroh, "Dislocations and cracks in anisotropic elasticity," Philosophical
Magazine, vol. 3(30): pp. 625-646, 1958.
[54] 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(3): pp. 427439, 1999.
[55] Q.H. Qin and Y.W. Mai, "Crack growth prediction of an inclined crack in a half-plane
thermopiezoelectric solid," Theor Appl Fract Mech, vol. 26(3): pp. 185-191, 1997.
[56] Q.H. Qin and M. Lu, "BEM for crack-inclusion problems of plane
thermopiezoelectric solids," International Journal for Numerical Methods in
Engineering, vol. 48(7): pp. 1071-1088, 2000.
[57] Q.H. Qin, "Thermoelectroelastic Green's function for a piezoelectric plate containing
an elliptic hole," Mechanics of materials, vol. 30(1): pp. 21-29, 1998.
[58] Q.H. Qin, "2D Green's functions of defective magnetoelectroelastic solids under
thermal loading," Engineering analysis with boundary elements, vol. 29(6): pp. 577-585,
2005.
[59] 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(4): pp. 250-258, 2003.
71
[60] K. Rajesh and B. Rao, "Two-dimensional analysis of anisotropic crack problems
using coupled meshless and fractal finite element method," International Journal of
Fracture, vol. 164(2): pp. 285-318, 2010.
[61] Q.H. Qin and Y.W. Mai, "BEM for crack-hole problems in thermopiezoelectric
materials," Engineering fracture mechanics, vol. 69(5): pp. 577-588, 2002.
[62] Q.H. Qin, Fracture mechanics of piezoelectric materials. Southampton: WIT press
2001.
[63] H. Sosa, "Plane problems in piezoelectric media with defects," International Journal
of Solids and Structures, vol. 28(4): pp. 491-505, 1991.
[64] 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(5): pp. 581-590, 1997.
[65] Q.H. Qin, Green's function and boundary elements of multifield materials. Oxford:
Elsevier Science, 2007.
[66] S.W. Yu and Q.H. Qin, "Damage analysis of thermopiezoelectric properties: Part II.
Effective crack model," Theor Appl Fract Mech, vol. 25(3): pp. 279-288, 1996.
[67] 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(6): pp.
433-444, 1998.
[68] Q.H. Qin, "General solutions for thermopiezoelectrics with various holes under
thermal loading," International journal of solids and structures, vol. 37(39): pp. 55615578, 2000.
[69] Q.H. Qin and X. Zhang, "Crack deflection at an interface between dissimilar
piezoelectric materials," International Journal of Fracture, vol. 102(4): pp. 355-370,
2000.
[70] S.W. Yu and Q.H. Qin, "Damage analysis of thermopiezoelectric properties: Part I—
crack tip singularities," Theor Appl Fract Mech, vol. 25(3): pp. 263-277, 1996.
[71] H. Ding, G. Wang, and J. Liang, "General and fundamental solutions of plane
piezoelectroelastic problem," Acta Mechanica Sinica (Chinese Edition), vol. 28: pp. 441448, 1996.
[72] W. Yao and H. Wang, "Virtual boundary element integral method for 2-D
piezoelectric media," Finite elements in analysis and design, vol. 41(9-10): pp. 875-891,
2005.
[73] 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(1):
pp. 139-158, 2005.
72
[74] G.S.A. Fam and Y.F. Rashed, "The Method of Fundamental Solutions applied to 3D
structures with body forces using particular solutions," Computational Mechanics, vol.
36(4): pp. 245-254, 2005.
73