Trefftz-type FEM for solving orthotropic potential problems

A simple Trefftz-type finite element method (TFEM) is proposed for solving certain potential problems in orthotropic media. The “body force”, which is induced by internal sources or sinks, may produce domain integrals in the standard Trefftz finite element formulation. This will make the advantage “only-boundary integration” of TFEM lose entirely. To overcome this difficulty, the dual reciprocity method (DRM) is employed to transfer the original problem into a homogeneous one. Then, a particular solution (PS) Trefftz-type finite element model is established based on the modified functional. Three benchmark examples are investigated by the proposed approach and compared with the analytical solutions.

As highlighted by Jirousek and Venkatesh (1992); Qin (2000;2008), the TFEM couples the advantages of FEM and BEM.Due to the fact that no domain integrals are involved in the formulation of TFEM, the Trefftz-type elements are less sensitive to mesh distortion in practical applications.This feature has been investigated by Jirousek and Wroblewski (1995); Jirousek and Qin (1996); Choi et al. (2006); Cen et al. (2011); Wang et al. (2012) using different 4-node quadrilateral elements.On the other hand, the employment of two independent fields also makes the TFEM easier to generate arbitrary polygonal or polyhedral elements with or without inclusions, which are accurate, efficient and natural for micromechanical modeling of heterogeneous materials (Dong and Atluri, 2012a,b,c).And the special-purpose elements of TFEM, which may achieve a level of popularity unequalled nowadays, can efficiently capture the stress concentration (or high gradients) without mesh refinement.These special elements with embedded complexities, which use complete general solutions in the domain, can greatly reduce the computational burden and preprocessing effort.Among the researchers who made contributions to special elements, the work of Piltner (1985); Jirousek and Venkatesh (1992); Leconte et al. (2010); Wang and Qin (2011;2012a;2012b) should be mentioned.Besides, there is also another kind of special element which can capture the singularity at crack-tip (Freitas and Ji, 1996;Kaczmarczyk and Pearce, 2009).Zhao and Zhao (2011) recently proposed a hybrid finite element model for anisotropic potential problems.In their model, the fundamental solutions are used as trial functions for intraelement field.Wang et al. (2012) developed a novel Trefftz finite element model, whose intraelement interpolation functions can reflect varying properties, for simulating heat conduction in nonlinear functionally graded materials (FGMs).Wang et al. (2012) recently focused their attention on the orthotropic potential problems.The original problem is mapped into an equivalent isotropic one by coordinate transformation so that Trefftz functions are readily obtained from Laplace equation.In the meantime, the original boundary conditions must also be transformed before their imposition on the new domain.For a potential problem with "body force", a domain integral will also be required.This may cause one of the key features of TFEM, namely its "onlyboundary integration" formulation, to lose entirely.To smooth out this difficulty, a dual reciprocity method (DRM) developed by Wrobel et al. (1986); Nardini and Brebbia (1987) has usually been used for handling the potential problems with "body force" in isotropic bodies by Cao et al. (2012); Kita (2005); Wang et al. (2012); Balakrishnan and Ramachandran (1999).However, relatively few contributions applying TFEM to orthotropic potential problems with "body force" can be found in the literature.
Motivated by the strength of TFEM and the popularity of DRM, this paper makes an effort in dealing with potential problems with "body force" in orthotropic media.In this methodology, Trefftz functions proposed by Wang et al. (2012) are employed to construct the intra-element field.The "body force" is represented by a series expansion in terms of radial basis functions Latin American Journal of Solids and Structures 11 (2014) 2537-2554 (RBFs) for which particular solutions can be easily determined.To alleviate the inconvenience as addressed in the work (Wang et al., 2013), an alternative modified functional has been established for homogeneous orthotropic problems so that the solution process can be conducted in the original domain.As a benchmark, three examples are numerically investigated and a fair agreement is found in comparison with analytical solutions.

Basic equations
Let us consider a 2D well-posed, orthotropic potential problem ( ) subjected to the Dirichlet boundary condition and the Neumann boundary condition where u and q are the potential and its derivative in normal direction, f is the term of "body force" induced by internal sources or sinks, 1 k and 2 k are the horizontal and vertical material coefficients, respectively, and they are assumed to be parallel to the major axes of anisotropy, 1 n and 2 n are direction cosines of the outward normal to the boundary, G denotes a bounded domain in 2 Â space with boundary G , For the sake of clarity, Eq. ( 3) is rewritten in the matrix form as

Methodology of DRM
In order to solve the potential problems with "body force" by boundary-type methods, it is necessary to eliminate the right-hand side in Eq. ( 1).This can be done by decomposing the solution to Eq. (1) into two parts, namely a particular solution p u and a homogeneous solution h u , such that in which ( ) and together with the modified boundary conditions To obtain the particular solution p u and homogeneous solution h u , we can rewrite the differential operator in Eq. (1) as follows (Wang et al., 2012) where By virtue of Eqs. ( 6) and (10) the particular solution p u can be straightforwardly expressed as where 1 r and 2 r are arbitrary reference values.It is noted that p u dose not necessarily satisfy the boundary conditions (2) and (3) and is not unique.Besides, the exact expression for p u can be explicitly obtained using Eq. ( 12) if ( ) f X,Y is constant or in a simple form.However, it is often intractable or even impossible to get the analytical derivation of ( ) is a general function.Thus, the approximate particular solution becomes necessary.To determine the particular solution ( ) p u X,Y , the right-hand side term ( ) f X,Y in Eq. ( 1) is usually approximated in the manner ( ) ( ) where L stands for the number of reference points in the solution domain, k a are the unknown Latin American Journal of Solids and Structures 11 (2014) 2537-2554 coefficients, and ( ) k X, Y j denote the basis functions in which the RBFs are selected in this paper.Now, the problem of finding a particular solution is reduced to where . Subsequently, the approaximate particular solution ( ) p u X,Y can be expressed as follows (Qin and Wang, 2008) by solving Eq. ( 15) directly.Hence, we can rewrite Eq. ( 14) in the following form where . Performing integration twice analytically for Eq. ( 16), the corresponding

( )
k X, Y F can be readily determined.Here, the power spline-type RBFs in 2 Â space such that are chosen.Therefore, we arrive at ( ) Consequently, are obtained according to Eq. ( 16), we can solve Eq. ( 13) for determination of the unknown coefficients k a by means of the singular value decomposition (SVD).Then, the particular solution ( ) p u X,Y can be evaluated at any field point from Eq. ( 15).The corresponding particular heat flux ( ) p q X,Y may be readily expressed as Latin American Journal of Solids and Structures 11 (2014) 2537-2554 3 TREFFTZ-TYPE FINITE ELEMENT FORMULATION

Assumed fields
To solve the orthotropic potential problem governed by Eqs. ( 1)-(3) using TFEM approach, the solution domain W has to be divided into a number of elements as done in the conventional FEM.For each element e occupied by a sub-domain e W , two independent fields (Jirousek and Qin, 1996;Qin and Wang, 2008;Wang et al., 2012), i.e. intra-element potential field and frame potential field, are assumed as shown in Figure 1.
(i) Intra-element potential field where e N is a vector of interior trial functions in which a finite terms only of homogeneous solutions to Eq. ( 1) are retained and are called Trefftz functions, e c is of the unknown parameters ej c , e  N is of the conventional shape functions and e d is of nodal degrees of freedom (DOF) of the element.The symbol "~" allows the two fields to be distinguished and ( ) x,y is the local Cartesian coordinate system.The proper number m of truncated set of homogeneous solutions ej N is chosen in such a way that (Qin and Wang, 2008) to avoid spurious zero-energy modes.Here, d m denotes the number of nodal DOF of the element.It should be noted that Eq. ( 23) is only a necessary but not a sufficient condition.In practice, more Trefftz functions are usually required to guarantee the resulting element stiffness matrix with full rank and to stabilize the performance of the element (Jirousek and Venkatesh, 1992).In the framework of TFEM, the homogeneous solution h u is approximated as a superposition of Trefftz functions k L that exactly satisfies the governing differential equation ( 7).For the Laplace equation over a 2D bounded domain, a set of homogeneous solutions to Eq. ( 1) is given by (Wang et al., 2012) where Truncating m terms of homogeneous solutions for Eq. ( 24), the vector of Trefftz functions e N can be explicitly written as It can be observed from Eq. ( 26) that 0 1 e N = , which represents the rigid-body motion mode, is excluded from e N in generating the sequence ej N .As a result, m should always be an even number (Qin and Wang, 2008).
From Eq. ( 21), the corresponding outward normal derivative of e u may be readily deduced as where x y On the other hand, the vector of frame functions e  N can be generated in customary way as done in the conventional FEM.For instance, at the point locating the side consisting of nodes 3, 4 and 5, the vector e  N and the nodal DOF vector e d can be defined as Latin American Journal of Solids and Structures 11 (2014) 2537-2554 ( ) )

The finite element stiffness equation
To establish the linkage between the two independent fields ( 21) and ( 22), a modified variational functional, which includes boundary integrations only, is constructed based on the work of Wang and Qin (2008): ( ) e q e I me eh eh eh eh e ep eh eh eh eh q u q u q q q u q u where e e u e q e I , and eI G is the inter-element boundary between elements.
Further, Eq. ( 34) may be readily rewritten in a compact form ( ) me eh eh eh eh e ep eh q u q u q q u Then, substitution of Eqs. ( 21), ( 22) and ( 27 To ensure good numerical conditioning of e H and to prevent overflow or underflow in evaluating the inverse of e H , the introduction of a local non-dimensional coordinates system ( ) , x h is usually suggested such that (Jirousek and Venkatesh, 1992;Qin and Wang, 2008) ( )  may be obtained by assembling Eq. ( 42) for all individual elements.After the modified Dirichlet boundary condition ( 8) is introduced, the homogeneous potential values h u of all nodes will be evaluated simultaneously by solving Eq. ( 43).Then, the coefficient vector e c can be computed from Eq. (41).

Recovery of the lacking rigid-body motion
It is necessary to recover the lacking rigid-body motion term when calculating the intra-element field eh u of any element.The discarded term 0 u can be readily reintroduced by setting for the augmented intra-element potential field (Jirousek and Venkatesh, 1992;Qin and Wang, 2008;Wang et al., 2012)  which finally leads to ( ) Once the rigid-body motion term 0 u is determined by Eq. ( 46), the full potential field e u at any internal point can be evaluated in combination with Eqs. ( 5), ( 15), ( 21) and ( 46).

NUMERICAL EXAMPLES
The 2D steady-state particular solutions have been incorporated into an in-house standard TFEM code.To validate the numerical implementation, solutions to three test problems are presented below: In the first two, the domain is a simple rectangle; the second involves a curved geometry which may be more representative of an actual system.The analytical solutions are also provided for the purpose of a fair comparison.In each example, eight-node quadratic elements are invoked for discretization and four Gaussian points are utilized along each element side.The particular solutions related to "body force" are approximated using the method of RBFs.Following the work of Qin and Wang (2008), all nodes and elemental centroids are chosen as the reference points in the following numerical examples.

Example 1: Rectangular temperature field with linear "body force" term
In the first example, we consider a 2D steady-state temperature field over a rectangular domain with length 1 L = and width 0.8 W = as illustrated in Figure 2. The Dirichlet boundary conditions are prescribed on the left and right surfaces, while the Neumann boundary conditions are specified at the rest of the surfaces.The heat conductivity coefficients are given by 1 1 k = and 2 4 k = .In this example, the "body force" term is assumed to be of linear variation along the X direction such that ( ) f X,Y X = -.This problem admits the analytical solution in the form As mentioned in Subsection 3.1, the choice of number of Trefftz functions, m , which may affect the accuracy and convergence of the procedure, is an important factor in the practical computation.For this purpose 16 elements with 65 nodes are first used to model the entire domain.The results at selected points are listed against different values of m in Table 1, from which it can be seen that large error will appear once m arrives at 18.The reason for this phenomenon is that too many Trefftz functions can result in the numerical overflow when calculating the square matrix e H .For the best compromise between accuracy and computational effort, the number of Trefftz functions, 10 m = , is selected in this example.It should be noted that the number of 10 m = is also the optimal one for remaining two examples based on lots of numerical experiments.Besides, the convergent performance is also investigated using three meshing densities.As expected, improved numerical accuracy is observed in Table 2 along with an increase in the number of elements.Finally, the sensitivity to mesh distortion of PS-TFEM is illustrated using the distorted scheme defined in Figure 3.The scheme is implemented on a distorted 4×4 mesh, controlled by .The case of 0 l < indicates that elements 6, 7, 10 and 11 distort towards the corner while 0 l > indicates that elements 1, 4, 13 and 16 distort towards the center.As a measure of sensitivity to mesh distortion, the relative error defined as below is adopted Table 3 displays the relative errors for temperature u and its derivative u X ¶ ¶ with the distortion parameter l .Obviously, elements 1, 4, 13 and 16 will become eight-node triangles when 0.125 l = -and the similar trend is for elements 6, 7, 10 and 11 when 0.125 l = .Once 0.125 l > , the associated elements will collapse to concave quadrilaterals.As is known, it is very intractable or even impossible for conventional isoparametric elements to treat this case because Jacobian matrix will be less than zero when the internal angle of an element is equal to and greater than 180°.This limitation may be easily eliminated by taking the advantage of the "onlyboundary integration" Trefftz finite element formulation.It is obviously observed from Table 3 that the relative errors for the temperature u are very close to zero, which means the results are not too sensitive to mesh distortion.Although the maximum relative errors for the temperature derivative u X ¶ ¶ at points (0.50, 0.40), (0.25, 0.20) and (0.25, 0.60) reaches -3.425%, 4.647% and 4.647% respectively, the results still meet the request of engineering precision (within 5%).

Example 2: Rectangular temperature field with quadratic "body force" term
As a second example, once again, a 2D steady-state temperature field over a rectangular domain as demonstrated in Figure 4 is considered.The Dirichlet boundary conditions are prescribed on the left and right surfaces while the Neumann boundary conditions are specified at the rest of the surfaces.However, the geometry is 3 L = in length and 2 W = in width and the heat conductivity coefficients are assumed to be 1 4 k = and 2 9 k = .Moreover, the "body force" term is assumed Latin American Journal of Solids and Structures 11 (2014) 2537-2554 to be of quadratic variation in the X-direction such that ( ) . For reference, the exact solutions of temperature is given by   For the study of this example, this solution domain has been idealized by 24 elements with 93 nodes according to Figure 4.The distribution for temperature u on the top surface together with heat flux q (= u n ¶ ¶ ) on the right surface is depicted in Figure 5 and Figure 6, respectively.One can scarcely see differences in the results between the PS-TFEM and analytical solutions.The maximum error is only 0.0085% for the temperature u on the top surface while 1.25% for the heat flux q on the right surface.The model has demonstrated good comparison between analytical and numerical results.

Example 3: Torsion of an elliptic shaft
To demonstrate the efficacy of the proposed method for curved geometries, the pure torsion of an elliptic shaft as illustrated in Figure 7 is investigated.The values of semi-major and semi-minor axes are 10 a = and 5 b = .The Dirichlet boundary conditions are prescribed on the outer surface of the shaft.In the finite element solution presented, the "body force" term of ( ) is explored for simplicity.The cross section is represented by a material for which the reciprocal values of shear modulus along the X-and Y axes are taken to be 1 4 k = and 2 1 k = .This problem admits the analytical solution in the form where the stress function formulation is deduced using the representation for stresses One particular solution may be exactly expressed by means of Eq. ( 12) such that Figure 8 shows the distribution of stress function u along the X-and Y axes. Figure 9 shows the distribution of shear stresses ZX s and ZY s along the outer surface.For clear comparison, the results for shear stresses at selected points are also presented in Table 4.A slight difference was observed between the results based on RBFs and on Eq. ( 52).

CONCLUSIONS
In this paper, we presented a particular solution Trefftz-type finite element approach for solving certain potential problems with "body force" in plane orthotropic materials.The original problem under consideration is first transferred into a homogeneous one using DRM.And then, a modified functional is constructed for the homogeneous orthotropic problem so that the Trefftz-type finite element formulation is derived.In doing so, the advantage 'only-boundary formulation' of TFEM is well preserved.Three examples are presented to verify the methodology.It is seen that the approach described in this paper to solve potential problems is, indeed, extremely accurate and robust in comparison with analytical solutions.Although the implementation is carried out for orthotropic potential problems, the idea and development are applicable to anisotropic cases.This work is underway.

Figure 1 :
Figure 1: Intra-element field and frame field for a particular element e.
change at the kth node, and is the element stiffness matrix with symmetric and positive definite characteristics.The calculations for e H , e G and e p can resort to the popular Gaussian quadrature performed along the entire element boundary.Finally, the whole stiffness equation of the system rigid-body potential 0 u can be calculated using the least square matching of eh u and eh u  at nodes on the entire element boundary e G

Figure 2 :
Figure 2: A rectangular temperature field with linear "body force" term.

Figure 3 :
Figure 3: Scheme for mesh distortion analysis.

Figure 4 :
Figure 4: A rectangular temperature field with quadratic "body force" term.

Figure 5 :
Figure 5: Distribution of the temperature u along the top surface.

Figure 6 :
Figure 6: Distribution of the outward normal heat flux q along the right surface.

Figure 8 :Figure 9 :
Figure 8: Distribution of stress function u along the X-and Y axes.
X and o Y are global Cartesian coordinates of element centroid, l n is the number of element nodes and c a the average distance between centroid and nodes of the element.To enforce inter-element continuity on the common element boundary, the unknown vector e c should be expressed in terms of nodal DOF e d .The stationary condition of the functional me P with respect to e c and e d , respectively, yields the following formulae Latin American Journal ofSolids and Structures 11 (2014) 2537-2554where o

Table 1 :
Results at selected points for different number, m, of Trefftz functions.
1 data in the 1st row denote the temperature u and the 2nd row its derivative u X ¶ ¶ (similarly herein- after)

Table 3 :
Relative errors for different mesh distortion parameters.

Table 4 :
A comparison of shear stresses ZX s between PS-TFEM and exact solutions.