Global-element-based free element method for solving non-linear and inhomogeneous heat conduction problems

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Global-element-based free element method for solving non-linear and inhomogeneous heat conduction problems Xiao-Wei Gao∗, Jin-Xing Ding, Miao Cui, Kai Yang Dalian University of Technology, State Key Laboratory of Structural Analysis for Industrial Equipment, China

a r t i c l e

i n f o

Keywords: Free element method FREM Global element Global free element method Finite element method Mesh free method

a b s t r a c t In this paper, a family of global elements (GEs) are constructed for modeling geometries and representing physical variables, based on a set of complete basis functions formulated in terms of normalized global coordinates. The main benefits of using GEs are that the elemental nodes can be distributed and numbered in an arbitrary manner and the global spatial partial derivatives of geometries and physical variables appearing in the governing equations of engineering problems can be directly derived. Based on the constructed GEs and their spatial derivatives of global shape functions, a simple and robust new numerical method, called as the Global-Element-based FRee Element Method (GEFREM), is proposed for solving general two-dimensional heat conduction problems. GEFREM inherits the advantages of the finite element method, mesh free method and free element method. A detailed description of GEFREM for solving general non-linear and inhomogeneous heat conduction problems will be presented in the paper and a number examples are given to verify the correctness and demonstrate the potential of the proposed method.

1. Introduction In the numerical analysis of engineering problems, the frequently used numerical method is the finite element method (FEM) [1–3], attributed to the features that different material properties can be easily imposed through elements and very stable results can be obtained due to the use of elemental nodes distributed in certain rule. In FEM, both the number and distribution of the elemental nodes are subjected to some constraints, the basic two of which are (1) the number of required nodes for an element must make the shape functions of the element complete, so that any problems with different orientation geometries can be modeled in a unified basis function; and (2) the nods defined over an element should distribute in certain manner to guarantee the stiffness matrix health. In one hand, these constraints can guarantee the results of FEM very stable. On the other hand, these constraints make FEM awkward to discrete irregular geometries [4] and model shock waves [5] and thermal concentrations [6,7]. In the standard FEM, the shape functions of elements are formulated in terms of selected parameters with fixed varying ranges, usually from −1 to +1, (they are also called as intrinsic coordinates), and they are used to represent both the problem geometry and related physical variables and, therefore, they are called as the isoparametric elements [1–3]. Since the parameters can also be directly used as the integration variables when utilizing the Gaussian quadrature to evaluate ∗

element stiffness matrixes, the isoparametric elements are extensively used in FEM. In the conventional FEM, the nodes in the intrinsic coordinate system are equally-spaced. However, it has been found that non-equally spaced nodes may give more accurate results [8]. Although the non-equally spaced elements can be constructed by taking the Lobatto points as nodes as done in [8], the intrinsic coordinates at the Lobatto points are still limited to the values solved from the equation of Legendre polynomials [9]. It is not easy to derive the shape functions for elements with arbitrarily spaced nodes which are very useful when reducing the number of elements at some locations where local refinement is required. In this aspect, the use of global elements (GEs) [10,11] can naturally overcome this difficulty, since the shape functions of the GEs are directly set up from the global coordinate system. The initial use of GEs was performed by Delves et al. [12,13], where the method was called the global element method (GEM) and can be viewed as a variable order, non-conforming finite element method. Later, GEM was applied by Krason and Szczecin to solve static field problems [14]. This method splits the given region into a number of elements which then are mapped on to a standard shape of region [−1,1] × [−1,1] and the solution is approximated by a product expansion in the mapped variables of each element. The basis functions of this type GEM are the Chebyshev polynomials and the solution scheme is based on the variational principle. Between the neighboring elements, the usual interface conditions are imposed.

Corresponding author. E-mail address: [email protected] (X.-W. Gao).

https://doi.org/10.1016/j.enganabound.2019.09.018 Received 4 June 2019; Received in revised form 17 July 2019; Accepted 24 September 2019 Available online 9 October 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

The above mentioned FEM and GEM are weak-form algorithms, in which the integration over each element is required. The most important feature of the weak-form algorithms is that very stable solutions can be obtained. The drawbacks of it are that a variational or energy principle is needed for a specific problem and some violent change phenomena such as the stress concentration in solid mechanics and shock waves in fluid mechanics are not easy to be captured. This may be the main reason that the dominant numerical methods in computational fluid mechanics are the finite different method [15] and the finite volume method [16], rather than FEM. Recently, the strong-form FEMs were proposed for solving engineering problems [17–19]. In these methods, the variational or energy principles are not required, nor the integration. Solution schemes are established by the governing equations and elemental interface conditions. The advantages of these methods are that the stress-concentration phenomenon can be more obviously embodied and a unified algorithm for solving multi-field problems can be more easily established than the conventional FEM [20]. To make up for the drawback of FEM in discretization of complicated and irregularly-shaped geometries, the mesh free method (MFM) was proposed and has been extensively used in solving various engineering problems [21–26]. The most attractive feature of MFM is that no elements with certain regularly distributed nodes are required to solve an engineering problem governed by partial differential equations (PDEs) and the mesh for a computational point can be independently formed by selecting a cluster of points around the computational point [27,28]. This is quite convenient to solve problems with complicated geometries or variable boundaries. An attractive and easy to implement MFM is the point interpolation method (PIM) applied to the strong form differential equations [24,25]. In this approach, polynomial [24] or radial basis functions [25] can be adopted to create the shape functions. The former has the Kronecker delta property and the latter is more stable. PIM has the advantages of very easy to use and requiring least number of nodes, hence it is most efficient. The drawback of the PIM is that its accuracy deteriorates when dealing with natural boundary conditions [29]. To improve the accuracy of the method, a weak strong form approach was proposed by Liu et al. [30]. The other drawbacks in PIM are that volume integrals may be required, the moment matrix could be singular, the construction process of shape functions may fail [24], and therefore special techniques are needed to overcome the singularity issue [29]. In addition, the approximation function may not be continuous over the global problem domain. To absorb the advantages of FEM and MFM, Gao et al. proposed a new type of collocation mesh free method, called as the FRee Element Method (FREM) [31], which employs a set of element differential formulations derived for general isoparametric elements [19,20] to compute the spatial partial derivatives. The distinct feature of FREM is that only one independent local element is needed for each collocation point, and the element is completely free, that is, the element not only can be formed by freely selecting certain surround nodes, but also is not connected to other points’ elements. This feature not only make the mesh generation much easy, but also give us the convenience to form the global system without need to consider the C1 -continuity of physical variables across the elements’ interfaces as required in the standard FEM and PIM. The free element method developed in [31] works based on the isoparametric elements as used in the standard FEM. A drawback exists in this type of elements, that is, the nodes of an element is equally spaced in the parameter space, which usually can give rise to a big computational error for elements with non-equally spaced nodes in the physical space. It is quite difficult to create the shape functions for high order isoparametric elements with nodes non-equally spaced in the parameter space. To absorb the advantages and overcome the drawbacks appearing in the above mentioned methods, a new strong-form mesh free method, called as the Global-Element-based FRee Element Method (GEFREM),

is proposed in the paper. The novelty and advantages of GEFREM over other existing methods can be identified through the following comparison. • Comparing to FEM, the proposed GEFREM uses only one free element for each point, which can overcome the drawback of FEM in need of considering the connection between the current element and other points’ elements when employing a variational principle or an energy principle to set up the solution scheme. And, due to the use of elements, GEFREM can inherit the advantage of FEM that stable results can be achieved from the inherent characteristic of elements that the variation of physical variables is consistent through the element’s nodes. • Comparing to PIM, the GEFREM also uses the complete basis functions to create the shape functions for each collocation point. However, GEFREM uses elements with fixed forms as used in the standard FEM, but PIM uses a cluster of points without fixed form. The benefit from the use of the fix-formed elements is that the different material properties can be easily applied at the different locations of the problem by specifying the properties on different elements, and the Neumann boundary conditions can be easily imposed to the boundary of the problem by specifying the boundary condition on the surface of the element having Neumann boundary conditions. Moreover, PIM directly uses the global coordinates to create the shape functions, which may cause the computation failed in some problems with huge geometry aspect ratio. In the proposed GEFREM, the normalized coordinates are used in the creation of nodal shape functions for an element, which can guarantee the computation successful. • Comparing to GEM, the GEFREM is a strong-form formulation, without integrations being involved and only one free element is used for each point. The existing GEM is a weak-form scheme, whole solution procedure of which is similar to the standard FEM, except for the approximation of physical variables in terms of the global coordinates. To fulfill the integration over an element, GEM employs the Chebyshev polynominals to approximate the physical variables, which is convenient to map the integration region into [−1,1] × [−1,1]. However, the proposed GEFREM uses the normalized global coordinates to represent the physical variables, and no any variational principles are required to set up the solution scheme. • Comparing to FREM, the GEFREM is an algorithm working in the global coordinate system space, while FREM works in the intrinsic (parameter) coordinate system space. The current one can automatically model the elements with non-equally-spaced nodes, and therefore is more stable and accurate than the FREM. In the following sections, the detailed description about the GEFREM will be given. Although the focus of the paper is paid on the treatment of the two-dimensional (2D) inhomogeneous and non-linear heat conduction problems, the theory can be easily applied to solve other kind of engineering problems. 2. Global elements formulated in normalized coordinates As in FEM, PIM and GEM, the variation of the temperature T with respect to the global coordinates (x, y) over an element can be expressed in terms of a complete set of basis functions 𝑓𝐼 (𝑥, 𝑦), i.e., ∑ 𝑇 (𝑥, 𝑦) = 𝑛𝐼=1 𝑎𝐼 𝑓𝐼 (𝑥, 𝑦) where n is the number of element nodes and 𝑎𝐼 are the coefficients to be determined through a simple collocation scheme as described late. In FEM [1–3] and PIM [24,25], the polynomials in (x, y) are used and in GEM the Chebyshev polynominals are used. After numerical investigations, it was found that when the aspect ratio of the problem is big or the element is distorted, the direct use of the global coordinates (x, y) may result in an inaccurate approximation and sometimes the computation may fail when determining the coefficients 𝑎𝐼 . For this reason, the normalized coordinates (𝑥̃ , 𝑦̃) are used in the paper, which have been proved very stable. Thus, the following function 118

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

13

15

16

9

10

11

12

5

6

7

8

1

(a) 9-node quadratic element

14

2

3

4

(b) 16-node cubic element

(c) 25-node quartic element

Fig. 1. Structured elements with different numbers of nodes.

approximation can be written. 𝑇 (𝑥, 𝑦) =

𝑛 ∑ 𝐼=1

𝑎𝐼 𝑓𝐼 (𝑥̃ (𝑥, 𝑦), 𝑦̃(𝑥, 𝑦))

2.1.2. Unstructured elements When the number of element nodes is not 𝑛 = 𝑛𝑥̃ × 𝑛𝑦̃ , the element is called as unstructured elements. In the sense of the isoparametric elements used in FEM, some of this type of elements are called as the Serendipity elements [1–3]. For example, the commonly used 8-node rectangular element is a kind unstructured element. According to the number of element nodes, the following types of unstructured elements are proposed, which have been successfully used during this study. Below each type of figures is followed by the approximation function which can be uniformly represent this type of elements and has been validated very stable for different distributions of element nodes. 4-node elements: ( ) 𝑇 (𝑥, 𝑦) = 𝑎1 + 𝑎2 𝑥̃ + 𝑎3 𝑦̃ + 𝑎4 𝑥̃ 2 + 𝑦̃2 + 𝑥̃ 𝑦̃ (5)

(1)

where, 𝑥̃ and 𝑦̃ are the normalized coordinates and the following expressions are used in the paper. 𝑥̃ (𝑥, 𝑦) =

𝑦 − 𝑦𝑐 𝑥 − 𝑥𝑐 , 𝑦̃(𝑥, 𝑦) = 𝑙𝑥 𝑙𝑦

(2)

in which, 𝑥𝑐 and 𝑦𝑐 are the average coordinates of all element nodes, and lx and ly are the reference lengths of the element under considered in x and y directions, respectively, which in this paper are taken as the differences between the 𝑥𝑐 and 𝑦𝑐 and the minimum coordinates of all the element nodes.

6-node elements:

( ) ( ) 𝑇 (𝑥, 𝑦) = 𝑎1 + 𝑎2 𝑥̃ + 𝑎3 𝑦̃ + 𝑎4 𝑥̃ 𝑦̃ + 𝑎5 𝑥̃ 2 + 𝑦̃2 + 𝑎6 𝑥̃ 2 𝑦̃ + 𝑦̃2 𝑥̃

2.1. Basis functions for different elements The basis functions fI used in Eq. (1) are determine by the number of nodes of the used element. In this study, two types of elements are introduced, which are called the structured and unstructured elements.

(6)

7-node elements:

( ) 𝑇 (𝑥, 𝑦) = 𝑎1 + 𝑎2 𝑥̃ + 𝑎3 𝑦̃ + 𝑎4 𝑥̃ 𝑦̃ + 𝑎5 𝑥̃ 2 + 𝑎6 𝑦̃2 + 𝑎7 𝑥̃ 2 𝑦̃ + 𝑥̃ 𝑦̃2 .

(7)

8-node elements: 2.1.1. Structured elements The structured elements refer here to the elements with structured nodal distributions, which are also called the Lagrange elements in FEM [1], EDM [19], FBM [32], and SFEM [33]. Fig. 1 shows three types of the structured elements with different numbers of nodes. For the structured elements, the function approximation can be written in a unified formulation as follows. 𝑇 (𝑥, 𝑦) =

𝑛𝑥̃ −1 𝑛𝑦̃ −1

∑ ∑

𝑏𝑖𝑗 𝑥̃ 𝑖 𝑦̃𝑗 =

𝑛 ∑

𝑎𝐼 𝑓𝐼 (𝑥̃ , 𝑦̃)

(3)

𝑓𝐼 (𝑥̃ , 𝑦̃) = 𝑥̃ 𝑖 𝑦̃𝑗 , 𝐼 = 1 + 𝑖 + 𝑗 ⋅ 𝑛𝑥̃ ( ) 𝑖 = 0, 1, ⋯ , 𝑛𝑥̃ − 1, 𝑗 = 0, 1, ⋯ , 𝑛𝑦̃ − 1

(4)

𝑖=0 𝑗=0

𝐼=1

𝑇 (𝑥, 𝑦) = 𝑎1 + 𝑎2 𝑥̃ + 𝑎3 𝑦̃ + 𝑎4 𝑥̃ 𝑦̃ + 𝑎5 𝑥̃ 2 + 𝑎6 𝑦̃2 + 𝑎7 𝑥̃ 2 𝑦̃ + 𝑎8 𝑥̃ 𝑦̃2

(8)

From Eqs. (5) to (8), it can be seen that the basis functions appearing in some terms of the function approximation expressions are not the monomials of the normalized coordinates. The benefit from the use of this type of basis functions is that the triangular, quadrilateral, and polygonal elements with or without internal nodes can be approximated by these functions in a unified way, regardless of the node distribution manner. As described later, in GEFREM, each element is served for collocation purpose at one point. The nodes numbers given in the above elements are only used for imposing the Neumann boundary conditions at boundary collocation points. In this case, the side numbers also should be assigned for the element over which the Neumann boundary conditions are imposed. For elements used for internal collocation points, actually no nodes numbers need to be defined over them, and only the nodal coordinates are required for internal elements. Since the free elements used in this paper are independent each other, they can be easily generated by writing up a simple pre-process code. However, for problems with complicated geometries, the preprocesses in existing commercial softwares such as ANSYS and ABAQUS can be used to generate the initial computational mesh, and then a simple transfer code can be written to transform the initial mesh into the model used in the current GEFREM. Both these two methods are used in this study. In principle, all the elements with 4, 5 or 6 edges as shown in Figs. 1–5 can be used in the presented method. However, the adequate

where 𝑎𝐼 = 𝑏𝑖𝑗 , and

where 𝑛𝑥̃ and 𝑛𝑦̃ are the numbers of nodes along the normalized directions 𝑥̃ and 𝑦̃, respectively. For example, they are 3, 4 and 5 for the three types of the elements shown in Fig. 1. Obviously, the total number of nodes shown in Eq. (3) is 𝑛 = 𝑛𝑥̃ × 𝑛𝑦̃ . It is noted that although the element shapes shown in Fig. 1 are regular rectangles, in the paper they can be arbitrarily oriented and their nodes can be arbitrarily distributed even distorted, provided that no any two nodes are too close each other. Also, since the proposed method belongs to the strong-form mesh free collocation method, it has been proved that even order basis functions have no the order discrepancy issue [34]. This is a reminder of using higher order elements. 119

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 2. The 4-node elements.

Fig. 5. The 8-node elements.

description of generating the shape functions for the above presented elements is given in the following. To determine the coefficients 𝑎𝐼 included in Eq. (1), collocating Eq. (1) at every node of the element under consideration, for instance, at node J results in 𝑛 ( ) ∑ ( ( ) ( )) (9) 𝑇 𝑥𝐽 , 𝑦𝐽 = 𝑎𝐼 𝑓𝐼 𝑥̃ 𝑥𝐽 , 𝑦𝐽 , 𝑦̃ 𝑥𝐽 , 𝑦𝐽 𝐼=1

Considering all n element nodes, Eq. (9) can be written as the following matrix form:

Fig. 3. The 6-node elements.

𝑭𝒂 = 𝑻 one should be the one which can well fit the geometry of the analyzed problem, especially at a corner node. Usually, the higher the order of the used element is, the higher the accuracy of the computational results is, of cause, the higher the difficulty is to generate the computational model.

(10)

where a = [a1 , a2 , a3 , … , an ⎡𝑓 (𝑥̃ 1 , 𝑦̃1 ) ⎢𝑓1 (𝑥̃ , 𝑦̃ ) 𝑭 =⎢ 1 2 2 ⎢ ⋮ ⎢𝑓 (𝑥̃ , 𝑦̃ ) ⎣ 1 𝑛 𝑛

𝑓2 (𝑥̃ 1 , 𝑦̃1 ) 𝑓2 (𝑥̃ 2 , 𝑦̃2 ) ⋮ 𝑓2 (𝑥̃ 𝑛 , 𝑦̃𝑛 )

]T ,

T = [T1 , T2 , T3 , … , Tn

⋯ ⋯ ⋮ ⋯

𝑓𝑛 (𝑥̃ 1 , 𝑦̃1 )⎤ 𝑓𝑛 (𝑥̃ 2 , 𝑦̃2 )⎥⎥ ⎥ ⋮ 𝑓𝑛 (𝑥̃ 𝑛 , 𝑦̃𝑛 )⎥⎦

]T ,

and

(11)

in which, TJ = T(xJ ,yJ ) and 𝑓𝐼 (𝑥̃ 𝐽 , 𝑦̃𝐽 ) represents 𝑓𝐼 (𝑥̃ (𝑥𝐽 , 𝑦𝐽 ), 𝑦̃(𝑥𝐽 , 𝑦𝐽 )). From Eq. (10), the coefficients 𝑎𝐼 can be solved as follows

2.2. Shape functions and derivatives of different elements

𝒂 = 𝑭 −1 𝑻

The most important task to develop a new numerical method based on the concept of elements is to generate shape functions for the proposed elements. The generated shape functions should be very stable (without singularity) for arbitrary distributions of the element nodes. The shape functions generated based on the above proposed function approximation expressions can achieve this purpose. The detailed

(12)

Thus, in terms of Eq. (1), the physical variable, temperature T, varying over the element can be expressed as [ ( )] 𝑛 𝑛 ∑ ∑ −1 𝑇 (𝑥, 𝑦) = 𝑓𝐼 𝑥̃ (𝑥, 𝑦), 𝑦̃(𝑥, 𝑦) 𝐹𝐼𝛼 𝑇𝛼 (13) 𝐼=1

Fig. 4. The 7-node elements. 120

𝛼=1

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 6. Free elements for different collocation nodes.

Fig. 9. Computed temperature using different orders of elements.

Fig. 7. Dimension and boundary conditions of the plate.

Fig. 10. Convergence of the residual Log10 (R) in iteration.

The shape functions formed in the above manner have the following Kronecker delta property: 𝑁𝛼 (𝑥𝐽 , 𝑦𝐽 )=1 𝑤ℎ𝑒𝑛 𝛼 = 𝐽 𝑁𝛼 (𝑥𝐽 , 𝑦𝐽 )=0 𝑤ℎ𝑒𝑛 𝛼 ≠ 𝐽 𝑛 ∑ 𝑁𝛼 (𝑥, 𝑦) = 1 𝛼=1

Attributed to this property, the temperature boundary conditions can be directly imposed [24]. In this study, the first and second order spatial derivatives of the shape functions are also needed. Thus, from Eq. (15), it can be easily derived that

Fig. 8. Model of 425 nodes used in GEFREM analysis of the plate.

Rearranging Eq. (13), we can express it as follows 𝑇 (𝑥, 𝑦) = 𝑁𝛼 (𝑥, 𝑦)𝑇𝛼

𝑛 𝜕 𝑁𝛼 (𝑥, 𝑦) 1 ∑ 𝜕 𝑓𝐼 (𝑥̃ , 𝑦̃) −1 = 𝐹𝐼𝛼 𝜕𝑥 𝑙𝑥 𝐼=1 𝜕 𝑥̃

(14)

𝑛 𝜕 𝑁𝛼 (𝑥, 𝑦) 1 ∑ 𝜕 𝑓𝐼 (𝑥̃ , 𝑦̃) −1 = 𝐹𝐼𝛼 𝜕𝑦 𝑙𝑦 𝐼=1 𝜕 𝑦̃

where 𝑁𝛼 (𝑥, 𝑦) =

𝑛 ∑ 𝐼=1

( ) −1 𝑓𝐼 𝑥̃ (𝑥, 𝑦), 𝑦̃(𝑥, 𝑦)𝐹𝐼𝛼

(15)

(16)

for the first order derivatives, and 𝜕 2 𝑁𝛼 (𝑥, 𝑦)

For simplicity, in Eq. (14) the repeated subscript 𝛼 represents the summation of all terms formed by 𝛼=1, 2, … , n and 𝑇𝛼 is the value of T at −1 is the Ith row and 𝛼th column value of node 𝛼, and in Eq. (15) the 𝐹𝐼𝛼 the inverse matrix of F shown in Eq. (11). The N𝛼 (x,y) shown in Eq. (15) is the shape functions for the element with n nodes. It not only can represent physical variables as shown in Eq. (14), but also can represent coordinates in a similar expression. It is noted that there are two cases which can make the inverse process of the matrix F shown in Eq. (12) failure, i.e., F is singular. One is the case that there are two element nodes coincided and the other is that most element nodes lie on a line.

𝜕𝑥2 𝜕 2 𝑁𝛼 (𝑥, 𝑦) 𝜕𝑦2 𝜕 2 𝑁𝛼 (𝑥, 𝑦) 𝜕 𝑥𝜕 𝑦

= = =

𝑛 2 1 ∑ 𝜕 𝑓𝐼 (𝑥̃ , 𝑦̃) −1 𝐹𝐼𝛼 2 𝑙𝑥 𝐼=1 𝜕 𝑥̃ 2 𝑛 2 1 ∑ 𝜕 𝑓𝐼 (𝑥̃ , 𝑦̃) −1 𝐹𝐼𝛼 2 𝑙𝑦 𝐼=1 𝜕 𝑦̃2 𝑛 2 1 ∑ 𝜕 𝑓𝐼 (𝑥̃ , 𝑦̃) −1 𝐹𝐼𝛼 𝑙𝑥 𝑙𝑦 𝐼=1 𝜕 𝑥̃ 𝜕 𝑦̃

(17)

for the second order derivatives. In terms of Eq. (1), the derivatives of fI with respect to the normalized coordinates 𝑥̃ and𝑦̃ appearing in Eqs. (16) and (17) can be easily obtained for a specific expression given by Eqs. (4)–(8). 121

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 12. Variation of L2 error with respect to nodal space in x-direction.

(a) Model of 1785 nodes with nodal space of 0.125cm

Fig. 13. Model with irregularly distributed nodes.

The following four types of boundary conditions related to the governing Eq. (18) are considered in this study.

(b) Model of 6969 nodes with nodal space of 0.0625cm

𝑇 (𝒙) = 𝑇̄ (𝒙), 𝒙 ∈ Γ1

Fig. 11. GEFREM models with 1785 and 6969 nodes.

3. Global free element method for non-linear heat conduction problems Using the function approximation Eq. (14) and the related first and second derivatives Eqs. (16) and (17) for the global elements shown in Figs. 1–5, the governing equations and related Neumann boundary conditions can be directly used to generate the system of equations for non-linear and nonhomogeneous heat conduction problems by employing the free element method.

The governing partial differential equation (PDE) for heat conduction problems can be expressed as [7,19] (

𝜆𝑖𝑗 (𝑇 , 𝒙) 𝜕𝑇𝜕 𝑥(𝒙) 𝑗

) + 𝑄 ( 𝒙) = 0 (18)

𝑜𝑟 2 𝜆𝑖𝑗 (𝑇 , 𝒙) 𝜕𝜕 𝑥𝑇𝜕(𝑥𝒙) 𝑖 𝑗

+

𝜕 𝜆𝑖𝑗 (𝑇 ,𝒙) 𝜕𝑇 (𝒙) 𝜕 𝑥𝑖 𝜕 𝑥𝑗

−𝜆𝑖𝑗 (𝑇 , 𝒙)

𝜕𝑇 (𝒙) 𝑛 (𝒙) = 𝑞̄(𝒙), 𝒙 ∈ Γ2 𝜕 𝑥𝑗 𝑖

(20a)

−𝜆𝑖𝑗 (𝑇 , 𝒙)

( ) 𝜕𝑇 (𝒙) 𝑛 (𝒙) = ℎ(𝒙) 𝑇 (𝒙) − 𝑇∞ , 𝒙 ∈ Γ3 𝜕 𝑥𝑗 𝑖

(20b)

−𝜆𝑖𝑗 (𝑇 , 𝒙)

( ) 𝜕𝑇 (𝒙) 𝑛 (𝒙) = 𝜎𝜀 𝑇 4 (𝒙) − 𝑇𝑠4 , 𝒙 ∈ Γ4 𝜕 𝑥𝑗 𝑖

(20c)

in which, Γ1 , Γ2 , Γ3 , and Γ4 are the boundaries specified respectively by temperature, heat flux, heat convection, and heat radiation boundary conditions, with Γ1 ∪Γ2 ∪Γ3 ∪Γ4 = Γ, Γ = 𝜕Ω, n is the outward normal to the boundary Γ, h is the heat transfer coefficient; 𝑇̄ ,𝑞̄ and T∞ are the prescribed temperature, heat flux and environmental temperature on the boundary, respectively; 𝜎 = 5.67051 × 10−8 W/(m2 K4 ) is the Stefan– Boltzman constant, and 𝜀 is the radiation interchange factor between the irradiated boundary Γ4 and the environment having a temperature 𝑇𝑠 . The repeated subscripts i and j in Eqs. (18) and (19) represent the summation through 1 to 2. The temperature boundary condition shown in Eq. (19) is called the Dirichlet boundary condition, and the flux boundary conditions shown in Eqs. (20a)–(20c) are called the Neumann, Robin, and radiation boundary conditions, respectively.

3.1. Governing equations for non-linear heat conduction problems

𝜕 𝜕 𝑥𝑖

(19)

+ 𝑄(𝒙) = 0

3.2. Collocation nodes and free elements

where Q is the heat source, and 𝜆ij is the heat conductivity which is the function of the temperature T for non-linear problems and/or coordinates x for nonhomogeneous problems.

From the governing Eq. (18), it can be seen that the PDE for heat conduction problems is the second order PDE of the physical variable, 122

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 14. Contour plot of computed temperature over the irregular model.

where x represents x and y, subscripts i and j take values from 1 to 2 with notation of x1 and x2 being x and y, respectively. The first and second order spatial derivatives of the shape functions in above equations are computed using Eqs. (16) and (17). As in the mesh free method, to solve the heat conduction problem, the computational domain Ω is discretized into a serious of arbitrarily distributed nodes in the GEFREM analysis and the system of equations is generated node by node. The node at which the equation will be generated for the system is called as the collocation node. For each collocation node, a free element is formed by the node and some selected surround nodes following the nodes arrangement shown in Figs. 1–5. It is noted that although the nodes to form the free element can be arbitrarily distributed, to achieve a high computational accuracy, the selected nodes should be as regular as possible. According to the different geometric positions of the collocation nodes located, three kinds of collocation nodes and corresponding free elements can be identified. Fig. 6 is a schematic example, in which the free elements are formed for a corner node (marked with red color), two smooth boundary nodes (marked with blue color), and two internal nodes (marked with black color). It is emphasized that for an internal node, the free element should have at least one internal node which can be used by the collocation node to establish the system of equations from the governing equations. However, for the corner or smooth boundary nodes, any types of isoparametric elements shown in this paper and those used in the standard FEM can be served as the free elements.

Fig. 15. Relative error of all models to the analytical solution.

temperature, with the global coordinates, and the related boundary conditions shown in Eqs. (19) and (20) involve the first order derivatives of the temperature. Therefore, the key point of solving this type of boundary value problems is to set up a scheme to compute the first and second order partial derivatives of the physical variables with respective to the global coordinates. The shape functions generated in this paper are formulated in terms of the normalized coordinates which have the linear relationship with the global coordinates (see Eq. (2)). This provide us with a great convenience to obtain the derivatives of the temperature with respect to the global coordinates. From Eq. (14), the following expressions can be easily derived: 𝜕 𝑁 𝛼 ( 𝒙) 𝜕𝑇 (𝒙) = 𝑇𝛼 𝜕 𝑥𝑖 𝜕 𝑥𝑖 𝜕 2 𝑁 𝛼 ( 𝒙) 𝜕 2 𝑇 ( 𝒙) = 𝑇 𝜕 𝑥𝑖 𝜕 𝑥𝑗 𝜕 𝑥𝑖 𝜕 𝑥𝑗 𝛼

3.3. Generation of the system of equations from the global free element method

(21)

Using the above described free elements, the system of equations can be easily generated for internal nodes and boundary nodes.

(22)

3.3.1. Generation of the system of equations for internal nodes The most important feature of the free element method is that the derived partial derivative expressions of the shape functions with 123

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 16. A half tube with a sector under heat convection and radiation boundary conditions.

respect to global coordinates can be directly substituted into the governing differential equations to generate the system of equations. To do this, substituting Eqs. (22) and (21) into the governing Eq. (18) yields [ ] 𝜕 2 𝑁𝛼 (𝒙) 𝜕 𝜆𝑖𝑗 (𝑇 , 𝒙) 𝜕 𝑁𝛼 (𝒙) 𝜆𝑖𝑗 (𝑇 , 𝒙) + 𝑇𝛼 + 𝑄(𝒙) = 0 (23) 𝜕 𝑥𝑖 𝜕 𝑥𝑗 𝜕 𝑥𝑖 𝜕 𝑥𝑗

where B is a known vector formed by the given heat fluxes and temperatures at all nodes, A is a square matrix formed by the coefficients of Eqs. (23) and (25), and E is a sparse matrix formed by the emission term (heat radiation) of Eq. (25c) . It is noted that for the nodes at which the temperature is specified, the column of the matrix A corresponding to these nodes are taken as 0 except for the diagonal ones which are taken as 1. In this case, the given temperatures are multiplied with their corresponding coefficients in Eqs. (23) and (25) and have been included in the right-hand side term B of Eq. (26). It is noted that for each node, only coefficients corresponding to the nodes of the element including the node are not zero, therefore, the matrix A in Eq. (26) is highly sparse. On the other hand, since the coefficient matrix A is the function of the unknown variable, temperature T, for non-linear problems, the system Eq. (26) needs to be solved through an iterative process. In this study, the Newton–Raphson iterative technique is employed to solve this non-linear equation set. In this technique, the temperature at each node is updated in the iterative process from a given initial values. We assume that the temperature vector is Tk after kth iteration. Based on the Tk , the residual of Eq. (26) can be written as

where x represents the coordinate at the collocation node which should be one of the element nodes, and the term 𝜕 𝜆ij /𝜕 xi can be calculated using the following formulation: 𝜕 𝜆𝑖𝑗 (𝑇 , 𝒙) 𝜕 𝑥𝑖

=

𝜕 𝑁 𝛽 ( 𝒙) 𝜕 𝑥𝑖

𝜆𝛽𝑖𝑗 (𝑇 )

(24)

in which, 𝜆𝛽𝑖𝑗 (𝑇 ) is the value of the conductivity 𝜆𝑖𝑗 at node 𝛽 when temperature is T, and the repeated subscript 𝛽 represents summation through all element nodes. 3.3.2. Generation of the system of equations for boundary nodes For nodes located on the outer boundary of the problem, the Dirichlet boundary condition shown in Eq. (19) can be directly used as an equation for a node where the temperature is specified. The flux boundary conditions shown in Eq. (20) are related to the first order spatial derivative of the physical variable. The system of equations for this type of boundary conditions can be generated by substituting Eq. (21) into Eq. (20). This can generate the following equations: −𝜆𝑖𝑗 (𝑇 , 𝒙)𝑛𝑖 (𝒙)

𝜕 𝑁𝛼 (𝒙) 𝑇𝛼 = 𝑞̄(𝒙), 𝒙 ∈ Γ2 𝜕 𝑥𝑗

( ) 𝜕 𝑁𝛼 (𝒙) −𝜆𝑖𝑗 (𝑇 , 𝒙)𝑛𝑖 (𝒙) 𝑇𝛼 = ℎ(𝒙) 𝑇 (𝒙) − 𝑇∞ , 𝒙 ∈ Γ3 𝜕 𝑥𝑗 −𝜆𝑖𝑗 (𝑇 , 𝒙)𝑛𝑖 (𝒙)

( ) 𝜕 𝑁𝛼 (𝒙) 𝑇𝛼 = 𝜎𝜀 𝑇 4 (𝒙) − 𝑇𝑠4 , 𝒙 ∈ Γ4 𝜕 𝑥𝑗

𝑹𝑘 = 𝑩 𝑘 − 𝑨𝑘 (𝑇 𝑘 )𝑻 𝑘 − 𝑬 (𝑻 𝑘 )4 To reduce the residual to a specified tolerance, forcing next iteration to be zero produces: 0 = 𝐑𝑘+1 = 𝐑𝑘 +

(25a)

in the

(28)

(25b)

where ΔT is the change in temperature. From above equation, we can obtain [ ] 3 (29) 𝑨𝑘 (𝑇 𝑘 ) + 4𝑬 (𝑻 𝑘 ) Δ𝑻 = 𝐑𝑘

(25c)

Solving above equation for ΔT, we can obtain the changes in temperature. Then, the temperature vector can be updated by 𝑻 𝑘+1 =𝑻 𝑘 + 𝑤Δ𝑻

In Eqs. (25b) and (25c), T(x) on the right-hand sides is the temperature at the collocation node x, which must be one of the nodal values of 𝑇𝛼 .

(30)

where w is a relaxation factor, which can be determined by an optimization process for a best convergence value [35]. In this paper, it takes the value of 1. The computational process from Eq. (27) to Eq. (30) is repeated for the next iteration, until convergence is achieved, i.e., until the following inequation holds true: √ 𝑇 𝑅 = (𝑹 𝑘 ) 𝑹 𝑘 ∕𝑁 < 𝑒 (31)

3.4. Solving the non-linear system using Newton–Raphson iteration method After applying the specified boundary conditions to Eqs. (23) and (25), the following system equations in the matrix can be obtained: 𝑨(𝑇 )𝑻 + 𝑬 𝑻 4 = 𝑩

𝜕 𝐑𝑘 Δ𝑻 = 𝐑𝑘 − 𝑨𝑘 (𝑇 𝑘 )Δ𝑻 − 4𝑬 (𝑻 𝑘 )3 Δ𝑻 𝜕𝑻 𝑘

(27) Rk + 1

(26) 124

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 17. Nodes and elements in GEFREM computation.

The top and bottom sides of the plate are insulated, while the left and right sides are imposed with the temperature condition of 𝑇̄ and the flux condition of 𝑞̄, respectively. Fig. 7 shows the dimension and boundary conditions of the plate. For this problem, the analytical solution can be derived as follows: [ ] ( ) 1 𝐻𝑐 𝑥3 𝐻𝑐 𝑥4 𝑐𝑥 𝐻𝐿 ̄ 𝑇 (𝑥) = − 𝐿𝑜𝑔 − − − 𝑞̄ + 𝑒−𝑐 𝑇 𝑐 6𝐿𝑘 𝑘 6 12𝐿2 𝑘

where N is the number of the total collocation nodes of the problem and e is a specified small tolerance. Once the values of temperature at all nodes are solved using the above iterative process, unknown heat fluxes at boundary nodes can be easily evaluated using Eqs. (25) and (21). 4. Numerical examples

where L = 6, c = 0.001, k = 20, 𝑇̄ = 800, 𝑞̄ = 400, H = 1.

A code has been developed for using the proposed method GEFREM to solver nonlinear heat conduction problems and two numerical examples are given to verify the correctness of the proposed method and examine the behavior of the proposed elements.

In the simulation using GEFREM, the plate is discretized into 425 nodes as shown in Fig. 8 with the space of adjacent nodes in x-direction being 0.25 cm. To investigate the behavior of GEFREM on the use of different orders of elements, the 25-node and 9-node structured elements as shown in Fig. 1, and 8, 7, 6- and 4-node unstructured elements as shown in Figs. 5(d), 4(b), 3(c) and 2(a), respectively, are respectively used at all internal node. Fig. 9 shows the computed results along the line of y = 2 cm using different orders of elements. By comparing to the Analytical result, it can be seen that the most elements can give good results, but the result

Example 1. A rectangular plate with non-linear conductivity and varying heat source The first example is concerned with a rectangular plate of 6 cm × 4 cm as shown in Fig. 7. The heat conductivity and source use here are as follows: ( ) 𝑥 𝑥 𝜆𝑖𝑗 (𝒙) = 𝜆(𝒙)𝛿𝑖𝑗 with 𝜆(𝒙) = 𝑘𝑒−𝑐𝑇 , 𝑄(𝒙) = 𝐻 1− 𝐿 𝐿 125

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

Fig. 18. Contours of the computed temperatures.

marked by GEFREM4 using the 4-node element has a bigger discrepancy with the analytical result. To show the convergence process of the residual R defined by Eq. (31) in the iterative computation, Fig. 10 gives the variation curve of the Log10 (R) with the iteration time. It can be seen that after 10 iterations, the residual R already decreased below 10−8 , which is the specified tolerance e in Eq. (31). To examine the mesh convergence performance of GEFREM to the number and distribution of discretized nodes, two finer computational models are also simulated, which are shown in Fig. 11. The numbers of nodes of the two models are 1785 and 6969, respectively, and the corresponding adjacent nodal spaces in x-direction are 0.125 cm and 0.0625 cm, respectively. Fig. 12 shows the variation of the L2 error with √ respect to the adjacent nodal space, which is defined [31] as 𝐿2 = (𝑻 − 𝑻̄ ) ⋅ (𝑻 − 𝑻̄ )∕𝑛 where T and 𝑻̄ are the computed temperature vector and the analytical solution vector at all nodes, with n being the number of total nodes. From Fig. 12, it can be seen that the proposed method has a good mesh convergence.

To investigate the stability of the proposed method to the irregular nodes distribution, a computational model with irregularly distributed nodes is generated, in which the nodal space of adjacent nodes in ydirection is not equally spaced and around four locations some nodes are deliberately distorted as shown in Fig. 13 or more obviously in the final overlapped elements show in Fig. 14. In this model, the number of nodes is 1785, which is the same as the model shown in Fig. 11(a). Fig. 14 also shows the contour of the computed temperature using the irregular model. From Fig. 14, it can be seen that all the contour lines along y-direction are straight, even if around the distorted locations. This indicates that the results of the present method are less dependent of nodes distribution, which proves that the stability of the method is good. To further examine the accuracy and stability of GEFREM, Fig. 15 plots the relative errors of all the four computational models shown in Figs. 8, 11 and 13, in which the Model2 (Non-uniform) refers to the irregular model shown in Fig. 13. Fig. 15 shows that all relative

126

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

overlapped elements are formed for all nodes, which are shown in Fig. 17. For comparison, the problem is also computed using the finite element method commercial software ANSYS, the model of which is constituted by 14,079 quadratic elements and 43,322 nodes. Fig. 18 is the contour plots of the computed temperatures using the commercial software ANSYS and current method GEFREM, and Figs. 19 and 20 show the variation curves of the temperature along line AB and BC, respectively, as marked in Fig. 17(a). From Fig. 18 it can be seen that the contour lines in both figures are very smooth and the pattern of GEFREM is the same as that of ANSYS. This indicates that the overall behavior of the proposed method is good. On the other hand, Figs. 19 and 20 show that the temperatures computed using GEFREM are very close to those using ANSYS. The maximum relative error is within 0.2%. These validates the correctness of the proposed method. Since the radiation boundary condition is applied on the outer sides of the sector, the problem becomes a nonlnear one (see Eq. (26)). To examine the iterative convergence of the proposed method in the radiation-induced non-linearity, Fig. 21 gives the variation of the residual with the number of iterations. It can be seen that only 3 iterations are needed to achieve the convergence under the specified tolerance e = 10−10 in Eq. (31).

Fig. 19. Variation curves of the computed temperature along line AB.

5. Concluding discussions A simple and robust new mesh free method has been presented in the paper for solving general nonlinear and nonhomogeneous heat conduction problems. The following conclusions and discussions can be drawn. (1) The shape functions of various orders of elements can be formulated in terms of only normalized global coordinates, not needing to use the radial basis functions (RBFs) to enhance the stability of the function approximation. (2) The two different types of nonlinear numerical examples, one being variable conductivity induced nonlinearity and other being the nonlinear boundary condition induced nonlinearity, indicate that the proposed method can effectively solve the general nonlinear heat conduction problems with high computational accuracy and excellent convergence. (3) The numerical Example 2 used in the paper has a very complicated geometry and boundary conditions. The correct computational results of this example indicate that the constructed global elements and the proposed free element method have the ability of treating problems with convex corners, concave corners, circular and triangular holes. This demonstrates the potential of the proposed method to solve complicated engineering problems. (4) A new representation form of the nonlinear system has been proposed, which explicitly include the radiation effect (see Eq. (26)). As a result, very fast iteration convergence can be achieved for the radiation induced nonlinear problems (see Fig. 16). (5) Although the 2D heat conduction problem is treated in the paper, the proposed elements and solution scheme can be easily used to solve other types of engineering problems.

Fig. 20. Variation curves of the computed temperature along line BC.

Fig. 21. Convergence of the residual Log10 (R) in iteration.

errors are within 0.001% and the finer the model is, the smaller the relative error becomes. Example 2. A half tube with a radiation sector The second example to be treated is a half tube with a radiation sector (see Fig. 16). The inner and outer radii of the tube are 240 cm and 300 cm, respectively. The boundary of the circular hole with the radius of 40 cm embedded in the sector is imposed with a heat convection condition of h = 20 W/(m2 K) and T∞ = 600 K and two straight outer sides of the sector are applied with the heat radiation boundary condition of 𝜎 = 5.67 × 10−12 W∕(cm2 K 4 ), 𝜀 = 0.8 𝑎𝑛𝑑 𝑇𝑠 = 300 K. The conductivity 𝜆 = 100 is used for both the tube and the sector.

Acknowledgments The supports of this investigation by the National Natural Science Foundation of China (11672061, 51576026, 11572075) and the Fundamental Research Funds for the Central Universities (DUT18GF105) are gratefully acknowledged. References [1] Zienkiewicz OC, Taylor RL. The finite element method. London, UK: McGraw-Hill; 1989. [2] Hughes TJR. The finite element method: linear static and dynamic finite element analysis. New Jersey: Prentice-Hall, Englewood Cliffs; 1987. [3] Belytschko T, Liu WK, Moran B, Elkhodary K. Nonlinear finite elements for continua and structures. 2nd updated and extended ed. Chichester: John Wiley & Sons Ltd; 2000.

The computational region of this problem is discretized into 7919 nodes, and the 25-node structured element as shown in Fig. 1(c) is used for most nodes; the 6-node element as shown in Fig. 3(a) is used at the tip node A marked in Fig. 17(a), and the 8-node element as shown in Fig. 5(b) is used at the concave nodes B and D. Eventually, 127

X.-W. Gao, J.-X. Ding and M. Cui et al.

Engineering Analysis with Boundary Elements 109 (2019) 117–128

[4] Yang K, Feng WZ, Wang J, Gao XW. RIBEM for 2D and 3D nonlinear heat conduction with temperature dependent conductivity. Eng Anal Bound Elem 2018;87:1–8. [5] Gao XW, Liu HY, Cui M, Yang K, Peng HF. Free element method and its application in CFD. Eng Comput 2019. https://doi.org/10.1108/EC-10-2018-0471. [6] Sladeka V, Sladeka J, Tanakab M, Zhang C. Transient heat conduction in anisotropic and functionally graded media by local integral equations. Eng Anal Bound Elem 2005;29:1047–65. [7] Divo E, Kassab AJ. Boundary element method for heat conduction: with applications in non-homogenous media. Southampton: WIT Press; 2003. [8] Wang DD, Li XW, Pan FX. A unified quadrature-based superconvergent finite element formulation for eigenvalue computation of wave equations. Comput Mech 2017;59:37–72. [9] Hu YC, Sze KY, Zhou YX. Stabilized plane and axisymmetric Lobatto finite element models. Comput Mech 2015;56:879–903. [10] Gao XW. Strong-form global element method and its application in mechanical analysis of functionally graded materials. ICM13, June 11–14, Melbourne, Australia; 2019. [11] Gao XW. Global-element-based free element method. ICCM2019, July 9–13, Singapore; 2019. [12] Delves LM, Hall CA. An implicit matching procedure for global element calculations. J Inst Maths Applies 1979;23:223–34. [13] Delves LM, Phillips C. A fast implementation of the global element method. JIMA 1980;25:177–97. [14] Krason P. Szczecin. Static field analysis by means of global method. Archiv fur Elektrotechnik 1986;69:93–6. [15] Liszka T, Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11(1–2):83–95. [16] Tao WQ, He YL, Wang QW, Qu ZG, Song FQ. A unified analysis on enhancing single phase convective heat transfer with field synergy principle. Int J Heat Mass Transf 2002;45:4871–9. [17] Wen PH, Cao P, Korakianitis T. Finite block method in elasticity. Eng Anal Bound Elem 2014;46:116–25. [18] Fantuzzi N, Tornabene F, Viola E, Ferreira AJM. A strong formulation finite element method (SFEM) based on RBF and GDQ techniques for the static and dynamic analyses of laminated plates of arbitrary shape. Meccanica 2014;49:2503–42. [19] Gao XW, Huang SZ, Cui M, Ruan B, Zhu QH, Yang K, Lv J, Peng HF. Element differential method for solving general heat conduction problems. Int J Heat Mass Transfer 2017;115:882–94.

[20] Gao XW, Li ZY, Yang K, Lv J, Peng HF, Cui M, Ruan B, Zhu QH. Element differential method and its application in thermal-mechanical problems. Int J Numer Methods Eng 2018;113(1):82–108. [21] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [22] Liu WK, Li SF, Belytschko T. Moving least-square reproducing kernel methods: I methodology and convergence. Comput Methods Appl Mech Eng 1997;143:113–54. [23] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [24] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. Int J Numer Methods Eng 2001;50:937–51. [25] Liu GR, Wang JG. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54:1623–48. [26] Zhang X, Liu XH, Song KZ, Lu MW. Least-square collocation meshless method. Int J Numer Methods Eng 2001;51(9):1089–100. [27] Sladek J, Sladek V, Atluri SN. Meshless local Petrov-Galerkin method for heat conduction problem in an anisotropic medium. Comput Model Eng Sci 2004;6(3):309–18. [28] Wang D, Wu J. An efficient nesting sub-domain gradient smoothing integration algorithm with quadratic exactness for Galerkin meshfree methods. Comput Methods Appl M 2016;298:485–519. [29] Liu GR. An overview on meshfree methods: for computational solid mechanics. Int J Comput Methods 2016;13(5):1630001. [30] Liu GR, Gu YT. A meshfree method: meshfree weak-strong (MWS) form method for 2-D solids. Comput Mech 2003;33(1):2–14. [31] Gao XW, Gao LF, Zhang Y, Cui M, Lv J. Free element collocation method: a new method combining advantages of finite element and mesh free methods. Comput Struct 2019;215:10–26. [32] Wen PH, Cao P, Korakianitis T. Finite block method in elasticity. Eng Anal Bound Elem 2014;46:116–25. [33] Fantuzzi N, Tornabene F, Viola E, Ferreira AJM. A strong formulation finite element method (SFEM) based on RBF and GDQ techniques for the static and dynamic analyses of laminated plates of arbitrary shape. Meccanica 2014;49:2503–42. [34] Wang DD, Wang JR, Wu JC. Superconvergent gradient smoothing meshfree collocation method. Comput Methods Appl Mech Eng 2018;340:728–66. [35] Cui M, Gao XW, Zhang JB. A new approach for the estimation of temperature-dependent thermal properties by solving transient inverse heat conduction problems. Int J Thermal Sci 2014;58:113–19.

128