MESHLESS SIMULATION OF HEAT CONDUCTION IN FUNCTIONALLY GRADED MATERIALS SUBJECTED TO TEMPERATURE-DEPENDENT HEAT SOURCES

– In this paper, an efficient meshless local B-spline based finite difference (FD) method for analysis of functionally graded materials (FGMs) subjected to temperature-dependent heat sources is presented. The favourable properties of B-spline basis functions in having arbitrary degree for resolution of solution, partition of unity and the Kronecker delta properties are combined with high accuracy and low computational effort of differential quadrature method in approximating shape functions and their derivatives. In this study, the FGMs are assumed to have temperature-dependent materials properties that vary as a function of radial distance. The homogenized properties are evaluated with power-law mixture rule. The nonlinearities from material properties and heat source terms are handled by the predictor-corrector method along with the Crank-Nicolson scheme for time integration. Case of nonlinear 2D heat conduction in FGM due to temperature dependent-heat sources is examined. The method is shown to be accurate and efficient for complex thermal analysis of FGMs taking into account temperature dependency of material properties and heat sources.


Introduction
Devising more accurate and efficient numerical schemes for thermal analysis of functionally graded materials (FGM) is of great interest. In the last decades, mesh reduction techniques so called meshless methods have emerged as attractive numerical techniques in engineering and science. The distinct feature of meshless methods is that the problems of discrete material property and discontinuity of solution are not encountered in the meshless methods due to the absence of mesh. Such an interesting feature is beneficial in handling material variation/gradation or in capturing thermal shock in the materials. Meshless methods hence have attracted great attention in FGM modelling (Sladek et al., 2003;Qian and Batra, 2005;Wang et al., 2006;Ching and Yen, 2006;Sladek et al., 2008;Skouras et al., 2011, Zhang et al., 2013Hidayat et al., 2015;Hidayat et al., 2017;Hidayat, 2019). However, only few reports on meshless analysis and investigation of FGMs accounting for temperature dependency of material properties can be found (Ching and Chen, 2007;Khosravifard et al., 2011;Mehrabadi and Aragh, 2013). The motivation and objective of the present work is to develop and present an efficient and effective meshless collocation approach for thermal analysis of FGMs taking into account temperature dependency of material properties heat sources.
In this paper, an efficient meshless local B-spline based finite difference (FD) method for analysis of FGMs subjected to temperature-dependent heat sources is presented. The favourable properties of B-spline basis functions in having arbitrary order/degree for resolution of solution, partition of unity and the Kronecker delta properties are combined with high accuracy, simplicity and low computational effort of differential quadrature in approximating shape functions and their derivatives. The method is thus truly meshless. In this study, the FGMs are assumed to have temperature-dependent materials properties that vary as a function of radial distance. The FGM homogenized properties are evaluated with simple power-law mixture rule. The nonlinearities from the temperature-dependent material properties and heat source terms are handled by using the predictor-corrector method along with the Crank-Nicolson implicit scheme for time integration.

Meshless Local B-spline based FD method
In meshless methods, a geometrical domain is discretized by a set of scattered nodes rather than by mesh connectivity. A local approximation is achieved by overlapped local domains covering the problem domain.  be an increasing sequence of real numbers. We will refer to the i τ as knots and to T as the knot vector. An interval [ ] 1 0 ,τ τ , which is allowed to be empty, will be called knot span. The i th B-spline basis function of order k (or degree , is defined recursively as (Cox, 1972;de Boor, 1972;de Boor, 2001;Piegl and Tiller, 1995;Farin, 2002): Here, t and s are parametric coordinates.
The interpolation condition is enforced at the ns collocation points in the subset i χ as: The interpolation matrix B has entries of the B-spline basis functions and u is the function values at the collocation points.
The function ( ) x u derivative can be approximated as: Evaluating the derivative at j x of the subset i χ and substituting where: L is the differential operator, ( ) ( ) then the meshless local B-FD collocation approach for the approximation of derivative can be written as: with weights vector w : Note that the local matrix B is non-square matrix. QR factorization or SVD techniques is employed to solve Eq. (10) in the least squares sense, from which the weight values for the node stencils can be obtained. In the next step, the vectors of weights w for all stencils are assembled together to form the differentiation matrix m D to give the integrated solutions ( ) x u over the problem domain Ω as: where: u is the vector of global unknowns and f is the vector of external forces. Eq. (11) is nothing but Eq. (8) in the form of global matrix assembled from the local systems of stencil.

Governing equation for general transient heat conduction in an isotropic solid body occupying a closed domain
Ω bounded by Γ is given by: The boundary conditions are: The initial condition is: where ( ) t T , x represents the temperature field, 0 T denotes the initial temperatures, t denotes time, ( ) x is the heat generation/source per unit volume at the internal point x and 1 Γ , 2 Γ and 3 Γ denote the three kinds of boundary conditions of first kind (specified temperature), second kind (specified heat flux) and third kind (convection heat transfer), respectively. Moreover, _ T and _ q are the prescribed temperature and the prescribed heat flux on the corresponding boundaries, respectively, h is the convection heat transfer coefficient, f T is the environmental temperature, and n is the unit normal vector outward to the boundary.
Temperature function ( ) t x T h , approximation and its time derivative by the B-spline shape functions (Hidayat, et al., 2015) are given as: For time integration, the following two-point difference time stepping technique or θ-method is used: Equation (17) then can be written as: where: Δt is the time step and q denotes the time level (i.e. Due to the temperature dependent material properties, a linearization scheme is required to solve Eq. (18). In this work, a predictor-corrector scheme based on direct substitution has been employed as the linearization scheme. The scheme has predictor and corrector steps as follows (Lewis and Roberts, 1987;Thakur et al. 2009):

Results and Discussion
Case of FGM made of zirconium oxide (ZrO 2 ) and titanium alloy (Ti-6Al-4V) due to with temperature dependent heat source is considered in this study. The nonlinear heat conduction problem is taken from Khosravifard et al. (2011). The problem geometry, boundary and initial conditions are depicted in Fig. 2. Further, Fig. 3 depicts the FE reference mesh (3356 nodes and 6464 elements) employed in this study. The numerical results obtained are compared with those of FEM with fine mesh as benchmark results. FEM is shown to be an effective method for thermal problems (Syahroni and Hidayat, 2011).
Effective property of the FGM at a point is computed by a simple mixture rule as follows: where: p 1 and p 2 are the material properties of Material 1 and Material 2, respectively, and v 1 and v 2 are the corresponding volume fractions for each material in the FGM. Volume fractions for the FGM constituents are: where: . Using the volume fractions relationship, the FGM property in the geometry is assumed to vary linearly with respect to the radial distance, from Ti-6Al-4V to ZrO 2 . Material properties of FGM used in this study are shown in Table 1.   (Touloukian, 1976 The heat source profile is depicted in Fig. 3. A quite sharp transition is observed in the heat source profile, thus inducing stronger nonlinearity in the analysis.
Temperature profile at several time instants with respect to different number of nodes are depicted in Fig. 5a, in which k = 7 is employed in the simulation. In Fig. 5b, the temperature variation along line P-Q in the FGM with respect to time step sizes are further shown. Temperature histories in the FGM is shown in Fig. 5c. Finally, the accuracies of the present method are given in Table 2.

Conclusions
In this study, meshless local B-spline based finite difference method has been presented for analysis of nonlinear heat conduction problem in FGM. Good accuracy and convergence are demonstrated in predicting temperature variation in the FGM. The method is shown to be accurate and efficient for complex thermal analysis of FGM taking into account temperature dependency of material properties and heat sources.