A New Meshless "Fragile Points Method" and A Local Variational Iteration Method for General Transient Heat Conduction in Anisotropic Nonhomogeneous Media
Yue Guan, Rade Grujicic, Xuechuan Wang, Leiting Dong, Satya N. Atluri
aa r X i v : . [ c s . C E ] J a n A New Meshless “Fragile Points Method” and A Local VariationalIteration Method for General Transient Heat Conduction inAnisotropic Nonhomogeneous Media
Yue Guan a, ∗ , Rade Grujicic b , Xuechuan Wang c , Leiting Dong d , Satya N. Atluri a a Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 79415, United States b Faculty of Mechanical Engineering, University of Montenegro, 81000 Podgorica, Montenegro c School of Astronautics, Northwestern Polytchnical University, Xi’an 710072, China d School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China
Abstract
A new and effective computational approach is presented for analyzing transient heat conduc-tion problems. The approach consists of a meshless Fragile Points Method (FPM) being utilizedfor spatial discretization, and a Local Variational Iteration (LVI) scheme for time discretization.Anisotropy and nonhomogeneity do not give rise to any difficulties in the present implementation.The meshless FPM is based on a Galerkin weak-form formulation and thus leads to symmetric ma-trices. Local, very simple, polynomial and discontinuous trial and test functions are employed. Inthe meshless FPM, Interior Penalty Numerical Fluxes are introduced to ensure the consistency of themethod. The LVIM in the time domain is generated as a combination of the Variational IterationMethod (VIM) applied over a large time interval and numerical algorithms. A set of collocationnodes are employed in each finitely large time interval. The FPM + LVIM approach is capable ofsolving transient heat transfer problems in complex geometries with mixed boundary conditions,including pre-existing cracks. Numerical examples are presented in 2D and 3D domains. Bothfunctionally graded materials and composite materials are considered. It is shown that, with suit-able computational parameters, the FPM + LVIM approach is not only accurate, but also efficient,and has reliable stability under relatively large time intervals. The present methodology represents aconsiderable improvement to the current state of science in computational transient heat conductionin anisotropic nonhomogeneous media.
Keywords:
Fragile Points Method, Numerical Flux Corrections, Local Variational Iteration ∗ Corresponding author.
Email address: [email protected] (Yue Guan)
Preprint submitted to January 27, 2020 ethod, Collocation Method, transient heat conduction, anisotropy, nonhomogeneity.
1. Introduction
Transient heat conduction analysis in nonhomogeneous media is of great interest in researchand engineering applications [1–3]. Typical nonhomogeneous media include conventional com-posite materials and functionally graded materials (FGMs), etc. In FGMs, the material propertiesvary gradually in space [4, 5]. These materials have been found widespread prospective applicationsin aerospace industry, computer circuit industry, etc[6]. Therefore, there is an increasing demandfor a reliable, accurate and efficient numerical approach for heat conduction problems in nonho-mogeneous materials. Due to various processing techniques, the FGM or composite material mayexhibit isotropic, orthotropic or anisotropic properties [4, 7]. Hence the anisotropy should also betaken into account.Numerical methods for solving transient heat conduction problems usually have two discretiza-tion stages [8, 9]. First, spatial discretization is employed. This step reduces the original partialdifferential equation (PDE) to a set of ordinary differential equations (ODEs) in time. These ODEs,known as “semi-discrete equations”, are then discretized in time using some standard ODE solvers.The most commonly used discretization methods are the Finite Element Method (FEM) [10] in thespatial domain and the finite difference schemes [11] in the time domain. However, in commercialFEM codes, the material properties are usually considered to be uniform in each element. It also hasdifficulties in analyzing systems with fragmentation, e.g., crack propagation in thermally shockedbrittle materials. The finite difference methods in the time domain, on the other hand, may havestability and accuracy problems while using large time steps. In general, the current commercialcodes are far from perfect.Apart from the FEM, other mesh-based numerical methods such as the Finite Volume Method(FVM) [12] and Boundary Element Method (BEM) [13, 14] can also be utilized in spatial discretiza-tion in heat transfer problems. In BEM, the non-availability of fundamental solutions in nonhomo-geneous anisotropic media is a serious limitation which is often insurmountable. The same as theFEM, the FVM and BEM approaches also have drawbacks in solving fragmentation problems suchas thermal-shock induced crack propagation in brittle solids. Their accuracy is also threatened whenmesh distortion occurs. Another category of methods, known as “meshless methods”, are partly or2ompletely free of mesh discretization. As a result, the human and computer cost in generating ahigh-quality contiguous mesh can be eliminated or reduced. This is a great improvement especiallyin 3D problems involving complex geometries. The Smoothed Particle Hydrodynamics (SPH) pro-posed by Randles and Libersky [15] is one of the earliest meshless methods. Though its originalformulation has a problem of stability and particle deficiency on and near the boundaries, a numberof improved methods based on the SPH have been carried out, including the modified SPH methodby Randles and Libersky [15], the Reproducing Kernel Particle Method (RKPM) by Liu et al. [16]and the Corrective Smoothed Particle Method (CSPM) by Chen et al. [17], etc. These methodsare extensively used in thermal analysis and fluid and solid mechanics. The SPH and its improvedmethods are based on a strong form, making it difficult to study their stability. Yet numerical testshave implied that the SPH method may turn unstable under random point distributions [18].Another category of meshless methods is weak-form-based, namely, employing the variationalprinciple to minimize the weighted residual of the governing differential equations. The DiffuseElement Method (DEM) [19] was initially introduced as a generalization of the FEM by remov-ing some limitations related to the trial functions and mesh generations. The original formulationfails in passing the patch test. Nevertheless, Krongauz and Belytschko [20] established an improvedDEM based on Petrov-Galerkin formulation (PG DEM) which satisfies the patch test but resultingin asymmetric matrices. Based on the DEM, the Element-Free Galerkin (EFG) method was carriedout by Belytschko et al. [21] in which the shape functions are developed by Moving Least Squares(MLS) or Radial Basis Function (RBF) approximations. The MLS approximation does not havedelta-function properties and hence imposition of essential boundary conditions is a serious lim-itation [22]. In addition, the integration of the Galerkin functional in EFG requires back-groundmeshes and is tedious. Furthermore, the EFG is not necessarily objective when the back-groundmeshes are rotated. After that, a Local Boundary Integral Equation (LBIE) Method is introducedby Zhu et al. [23]. While analogous to the BEM, the LBIE method circumvents the problem ofglobal fundamental solutions in nonhomogeneous anisotropic materials in BEM. It uses local fun-damental solutions (assuming locally homogeneous material properties) or Heaviside functions astest functions. However, the method still has drawbacks in solving fragmentation problems and hasasymmetric matrices. Finally, Atluri and Zhu [24] proposed the Meshless Local Petrov-Galerkin(MLPG) approach in 1998. Compared with the EFG method, the MLPG approach employs the lo-3al Petrov-Galerkin weak-formulation instead of the global Galerkin weak-formulation. The MLPGmethod is a truly meshless method and has shown its capability and accuracy in 2D and 3D transientheat conduction analysis involving anisotropy, nonhomogeneity and temperature-dependent mate-rial properties [4, 25]. Yet the MLPG method still has its limitations: the matrices are asymmetric;the computation of integration in the Petrov-Galerkin weak-form over local subdomains is compli-cated, as a result of the complex shape functions; and (the same as the EFG method) the essentialboundary conditions cannot be imposed directly, since the MLS approximation usually does notpass its corresponding data points. A modified collocation method has to be applied to enforce theessential boundary conditions [22].In contrast to the complex, global continuous MLS approximated shape function used in the EFGand MLPG, local, simple, polynomial, piecewise continuous trial functions are applied in generatingthe Fragile Points Method (FPM) in [26]. Nevertheless, the method would become inconsistent ifthe Galerkin weak-form is employed directly with these discontinuous trial and test functions. TheNumerical Flux Correction, which is widely used in Discontinuous Galerkin (DG) methods [27, 28],is introduced to remedy the problem. Whereas on the other hand, the inherent discontinuity can alsobe a benefit. Since it is convenient to relax the continuity requirement between neighboring points,the FPM has a great potential in analyzing systems involving cracks, ruptures and fragmentations,such as in problems of thermal shock in brittle solids. The method has already shown its stability,accuracy and efficiency in solving 1D and 2D Poisson equations [26] and elasticity problems [29].In the current work, it is further extended to 2D and 3D heat transfer analysis for discretizationin space. The Galerkin functional in the FPM can be integrated quite simply, and finally leads tosymmetric matrices. Thus, the FPM based on a Galerkin weak-form is far more efficient than eitherthe EFG or the MLPG method.After the spatial discretization is carried out, a semi-discrete system is achieved. The final com-puting accuracy and efficiency are also significantly affected by the ODE solver employed in thetime domain. Though nonlinearity is not emphasized in this paper, here we focus on generalizedtime discretization methods that can deal with linear as well as nonlinear ODEs. These methodscan roughly be divided into two categories: 1). The finite difference method, which is simple andthe most widely used, and 2). weak-form method based on weighted residual approximations [30].The most well-known finite difference schemes include: the central, forward and backward differ-4nce schemes [31], Runge-Kutta method [32], Newmark- β method [33] and Hilber-Hughes-Taylor(HHT)- α method [34]. The Houbolt’s method [35] based on a third-order interpolation is also fa-mous in dynamic analysis. Ode45 , the most popular, highly optimized built-in ODE solver in MAT-LAB, is based on an explicit Runge-Kutta formula [36, 37] and has excellent performance in solvinglow-dimensional nonlinear ODEs. However, we aim here at a high-dimensional semi-discrete sys-tem generated with the 2D or 3D FPM. The efficiency of the classic finite difference methods maynot be sufficient then. Moreover, many finite difference methods may encounter stability problemswhen considering nonlinearity or under large time steps [33, 34, 38].On the other hand, the weak-form methods, though somewhat hard to implement, have a poten-tial in analyzing high dimensional and nonlinear systems more efficiently. For periodic systems, theHarmonic Balance (HB) method [39] and the Spectral Time Domain Collocation (TDC) method[40] were developed. However, periodic responses are rarely seen in heat conduction problems.Considering more generalized transient solutions, a series of analytical or semi-analytical asymp-totic methods are introduced. He [41] proposed the Variational Iteration Method (VIM) which canbe seen as an extension of the Newton-Raphson method to nonlinear algebraic equations (NAEs).The Adomian Decomposition Method (ADM) is then developed by Adomian [42], in which theinitial guess is corrected step by step by adding components of an Adomian polynomial. Followingthat, the Picard Iteration Method (PIM) is carried out by Fukushima [43] and modified by Wool-lands et al. [44]. The method is also based on an initial guess and correctional iterative formula. Theformula is relatively concise, yet computing the integral of nonlinear terms in each time step is achallenge. Though developed independently, the previous VIM, ADM and PIM approaches can beunified using a generalized Lagrange multiplier [45]. Based on that, Wang et al. [46] employed theVIM over a finitely large time interval, along with a collocation method and numerical discretiza-tion, leading to the Local Variational Iteration Method (LVIM). Unlike the VIM, the LVIM is anumerical method, applicable to digital computation and has the potential in being implementedwith parallel processing. The initial guess can be constructed in a relatively simple form. Themethod possesses excellent efficiency in solving nonlinear ODEs in fluid mechanics, structural me-chanics, and astrophysics, etc [46]. Several approximated algorithms can be generated to furtherimprove the computing efficiency. In the current work, in order to maintain the best stability, onlythe classic LVIM based on the first kind of Chebyshev polynomials is considered. This method is5lso named as Chebyshev Local Iterative Collocation - 1 (CLIC-1) algorithm in [45].In this paper, we focus on transient heat conduction problems in anisotropic nonhomogeneousmedia. The system is discretized by using the FPM in space, and the LVIM in the time domain.Section 2 presents the governing equations for heat conduction in anisotropic nonhomogeneousmedia, and the formulation of the FPM. The LVIM and its numerical implementations are introducedin section 3. Numerous 2D and 3D examples, solved with the proposed FPM + LVIM approach,are carried out in section 4, followed by a discussion on the computational parameters, and a briefconcluding section.
2. Fragile Points Method (FPM) based on a Galerkin Weak-Form and Point “Stiffness” Ma-trices
Consider a transient heat conduction problem in a continuously anisotropic nonhomogeneousmedium, which is governed by the following partial differential equation [4, 47]: ρ ( x ) c ( x ) ∂u∂t ( x , t ) = ∇ · [ k ∇ u ( x , t )] + Q ( x , t ) , x ∈ Ω , t ∈ [0 , T ] , (1)where Ω is the entire 2D or 3D domain under study, the coordinate vector x = [ x y ] in 2D or x = [ x y z ] in 3D, u ( x , t ) is the temperature field, and Q ( x , t ) is the density of heat sources. ∇ isthe gradient operator. The thermal conductivity tensor k ( x ) , mass density ρ ( x ) and specific heatcapacity c ( x ) are dependent on the spatial coordinates in nonhomogeneous media. The thermalconductivity tensor components k ij can also be directionally dependent for anisotropic materials.Three main kinds of boundary conditions are considered:1). Dirichlet bc : u ( x , t ) = e u D ( x , t ) , on Γ D , : q ( x , t ) = k ∇ u ( x , t ) · n = e q N ( x , t ) , on Γ N , : q ( x , t ) = h ( x ) [ e u R ( x ) − u ( x , t )] , on Γ R , (2)where the global boundary ∂ Ω = Γ D ∪ Γ N ∪ Γ R , n is the unit outward normal of ∂ Ω , h ( x ) is the heattransfer coefficient, and e u R ( x ) is the temperature of the medium outside the convective boundary.The initial condition is assumed as: u ( x , t ) | t =0 = u ( x , , in Ω ∪ ∂ Ω . (3)6 .2. Local, polynomial, point-based discontinues trial and test functions In the domain Ω , a set of random points are introduced. The global domain can then be par-titioned into several confirming and nonoverlapping subdomains, with only one point in each sub-domain (as shown in Fig. 1). The subdomains could be of arbitrary geometric shapes. And thepartition is not unique. For simplicity, in this paper, the Voronoi Diagram partition [48] is applied.Unlike the FEM and other element-based methods, the shape and trial functions in the present Frag-ile Points Method (FPM) are totally established based on the random points and are independentof the domain partition . As a result, the temperature field could be discontinuous between subdo-mains, as well as the shape and trial functions. In FEM on the other hand, the trial and test functionsare element-based, and are continuous at the interelement boundaries.In each subdomain, we define the simple, local, polynomial trial function in terms of temperature u and its gradient at the internal point. For instance, the trial function u h in subdomain E whichcontains an internal point P can be written as: u h ( x ) = u + ( x − x ) · ∇ u (cid:12)(cid:12)(cid:12) P , x ∈ E (4)where u is the value of u h at P , and x is the coordinate vector of P .The gradient of temperature ∇ u at point P remains unknown. In this paper, we employ theGeneralized Finite Difference (GFD) method [49] to estimate ∇ u (cid:12)(cid:12) P in terms of the value u h atseveral neighboring points of P . Unlike the common definition of the support of P which includesall the points P ∈ { P ( x ) | ( x − x ) ≤ r } (as shown in Fig. 2(a), where r is a constant radius), inthe present work, the support of P is defined to involve all the nearest neighboring points of P insubdomains sharing boundaries with E in the Voronoi partition (shown in Fig. 2(b)). The pointsare named as P , P , · · · , P m .In order to estimate the gradient ∇ u (cid:12)(cid:12) P , we minimize the following weighted discrete L norm J : J = m X i =1 (cid:20) ( x i − x ) · ∇ u (cid:12)(cid:12)(cid:12) P − ( u i − u ) (cid:21) w i , (5)where x i donates the coordinate vector of P i , u i is the value of u h at P i , and w i is the value ofweight function at P i ( i = 1 , , · · · , m ). For convenience, we assume constant weight functions in7 (cid:2868) (cid:1842) (cid:2869) (cid:1842) (cid:2870) (cid:1842) (cid:2871) (cid:1842) (cid:2872) (cid:1842) (cid:2873) (cid:1842) (cid:2874) (a) (cid:1842) (cid:2868) (cid:1842) (cid:2869) (cid:1842) (cid:2870) (cid:1842) (cid:2871) (cid:1842) (cid:2872) (cid:1842) (cid:2873) (cid:1842) (cid:2874) (cid:1842) (cid:2875) (b)(c)Figure 1: The domain Ω and its partitions. (a) 2D domain with points distributed inside it ( P ∈ Ω ). (b) 2D domainwith points distributed inside it and on its boundary ( P ∈ Ω ∪ ∂ Ω ). (c) 3D domain with points distributed inside it( P ∈ Ω ). a) (b)Figure 2: Two kinds of support of a point in 3D. this paper. Hence, the temperature gradient at P is solved as: ∇ u (cid:12)(cid:12)(cid:12) P = Bu E , (6)where u E = h u u u · · · u m i T , B = ( A T A ) − A T h I I i , I = (cid:18)h − − · · · − i × m (cid:19) T , I = · · ·
00 1 . . . ...... . . . . . . · · · m × m , A = x − x x − x · · · x m − x . Therefore, the relation between u h and u E can be obtained: u h ( x ) = Nu E , x ∈ E (7)where N is called the shape function of u h in E : N = [ x − x ] B + h · · · i × ( m +1) (8)9hus the shape function is defined independently in each subdomain, and no continuity require-ment exists at the internal boundaries. Fig. 3 shows the graphs of shape functions in 2D and 3Ddomains respectively. The trial function u h can be derived in each subdomain E i ∈ Ω by the sameprocess. Thus u h can also be discontinuous at the internal boundaries. For instance, the trial func-tions simulating an exponential function u a = e − | x − x a | is shown in Fig. 4. As can be seen, the trialfunction is a simple local polynomial and just piecewise-continuous in the entire domain. The testfunction v h in the Galerkin weak-form in FPM is prescribed to have the same piecewise-continuousshape as u h .Unfortunately, the discontinuous trial and test functions will lead to an inconsistent and inaccu-rate result under the traditional Galerkin weak-form. To resolve that problem, we introduce Numer-ical Flux Corrections to the present FPM. (a) (b)Figure 3: The shape functions. (a) 2D domain with 49 points. (b) 3D domain with 343 points ( x, y, z ≤ . ). The Numerical Fluxes are widely used in Discontinuous Galerkin Methods [28] to help improv-ing their accuracy and stability. In our work, the Interior Penalty (IP) Numerical Flux Correctionsare employed.The governing equation Eqn. 1 can be written in the local weak-form with test function v in each10 a) (b)Figure 4: The trial functions for u a = e − | x − x a | . (a) 2D domain with 49 points, x a = [0 . . . (b) 3D domain with343 points ( x, y, z ≤ . ), x a = [0 . . . . subdomain E ∈ Ω . After applying the Gauss divergence theorem, we can get: Z E ρcv ∂u∂t dΩ + Z E ∇ v T k ∇ u dΩ = Z E Qv dΩ + Z ∂E v n T k ∇ u dΓ , (9)where ∂E is the boundary of the subdomain, and n is the unit vector outward to ∂E .Let Γ donate the set of all internal and external boundaries, i.e., Γ = Γ h + ∂ Ω = Γ h + Γ D +Γ N + Γ R , where Γ h is the set of all internal boundaries. We sum Eqn. 9 over all subdomains andrewrite it with the jump operator [] and average operator {} : X E ∈ Ω Z E ρcv ∂u∂t dΩ + X E ∈ Ω Z E ∇ v T k ∇ u dΩ = Z Ω Qv dΩ + X e ∈ Γ h (cid:18)Z e { v } (cid:2) n T k ∇ u (cid:3) + [ v ] (cid:8) n T k ∇ u (cid:9)(cid:19) dΓ + X e ∈ Γ D ∪ Γ N ∪ Γ R Z e [ v ] (cid:8) n T k ∇ u (cid:9) dΓ , (10)where the jump operator [] and average operator {} are defined as (for ∀ w ∈ R ): [ w ] = w (cid:12)(cid:12)(cid:12) E e − w (cid:12)(cid:12)(cid:12) E e e ∈ Γ h w (cid:12)(cid:12)(cid:12) e e ∈ ∂ Ω , { w } = (cid:18) w (cid:12)(cid:12)(cid:12) E e + w (cid:12)(cid:12)(cid:12) E e (cid:19) e ∈ Γ h w (cid:12)(cid:12)(cid:12) e e ∈ ∂ Ω . When e ∈ Γ h ( e ∈ ∂E ∩ ∂E ), n ej is a unit vector normal to e and pointing outward from E j (seeFig. 5). 11 (cid:2869) (cid:1842) (cid:2870) (cid:2196) (cid:2869) (cid:2196) (cid:2870) (a) (b)Figure 5: The inner boundary and normal vectors. (a) 2D case. (b) 3D case. Substituting the boundary conditions (Eqn. 2) into the last term in Eqn. 10, we obtain: X e ∈ Γ N Z e [ v ] (cid:8) n T k ∇ u (cid:9) dΓ = X e ∈ Γ N Z e v e q N dΓ , X e ∈ Γ R Z e [ v ] (cid:8) n T k ∇ u (cid:9) dΓ = X e ∈ Γ R Z e hv e u R dΓ − X e ∈ Γ R Z e hvu dΓ . (11)When u is the exact solution, since there is no âĂŸjumpâĂŹ on the internal boundaries, (cid:2) n T k ∇ u (cid:3) =0 . Similarly, [ u ] = 0 . This leads to (cid:8) n T k ∇ v (cid:9) [ u ] = 0 . Hence, we can replace the term { v } (cid:2) n T k ∇ u (cid:3) in Eqn. 10 by (cid:8) n T k ∇ v (cid:9) [ u ] without influencing the accuracy of the formula.Two Internal Penalty Numerical Fluxes are applied on Γ h and Γ D with different penalty param-eters η and η respectively. The formula of the FPM with IP Numerical Flux Corrections can thenbe achieved: X E ∈ Ω Z E ρcv ∂u∂t dΩ + X E ∈ Ω Z E ∇ v T k ∇ u dΩ − X e ∈ Γ h ∪ Γ D Z e (cid:0)(cid:8) n T k ∇ u (cid:9) [ v ] + (cid:8) n T k ∇ v (cid:9) [ u ] (cid:1) dΓ+ X e ∈ Γ R Z e hvu dΓ + X e ∈ Γ h η h e Z e [ u ] [ v ] dΓ + X e ∈ Γ D η h e Z e uv dΓ= X E ∈ Ω Z E Qv dΩ + X e ∈ Γ D Z e (cid:18) η h e v − n T k ∇ v (cid:19) e u D dΓ + X e ∈ Γ N Z e v e q N dΓ + X e ∈ Γ R Z e hv e u R dΓ , (12)where h e is a boundary-dependent parameter with the unit of length. For instance, h e can be definedas the length of the boundary (in 2D), the square root of the boundary area (in 3D), or the distance12etween the points in subdomains sharing the boundary. The penalty parameters η , η are positivenumbers having the same unit of k and independent of the boundary size. The method is onlystable when the penalty parameters are large enough. However, on the other hand, an excessivelylarge η is harmful for the accuracy and may cause a condition number problem. A discussion onrecommended values of these penalty parameters is presented in section 4.3. The IP NumericalFlux Correction terms vanish when u h equals to the exact solution, that is, when there is no jumpon internal boundaries and the Dirichlet boundary conditions are well satisfied.There are two ways to impose the Dirichlet boundary conditions in practice. If there are nopoints distributed on the boundaries (as shown in Fig. 1(a)), the Interior Penalty terms on Γ D areresponsible for the boundary conditions. The approach is analogous to the collocation method in-troduced in [22]. Alternatively, if boundary points are employed (as shown in Fig. 1(b)), we canalso impose u = e u D strongly at the boundary points and thus the corresponding Internal Penaltyterms in Eqn. 12 vanish. The formula of the FPM can be written in the matrix form finally: C ˙ u + Ku = q , (13)where C and K are the global heat capacity and thermal conductivity matrices respectively, u is theunknown vector with nodal temperatures, q is the heat flux vector.Substituting the shape function N for u h and v , B for ∇ u and ∇ v , the point heat capacitymatrix C E , point thermal conductivity matrix K E and the boundary thermal conductivity matrices13 h , K D , K R can be written as: C E = Z E ρc N T N dΩ , E ∈ Ω , K E = Z E B T kB dΩ , E ∈ Ω , K h = − Z e (cid:0) N T1 n T1 kB + B T1 k T n N (cid:1) dΓ + η h e Z e N T1 N dΓ − Z e (cid:0) N T2 n T2 kB + B T2 k T n N (cid:1) dΓ + η h e Z e N T2 N dΓ − Z e (cid:0) N T1 n T1 kB + B T1 k T n N (cid:1) dΓ − η h e Z e N T1 N dΓ − Z e (cid:0) N T2 n T2 kB + B T2 k T n N (cid:1) dΓ − η h e Z e N T2 N dΓ , e ∈ ∂E ∩ ∂E , K D = − Z e (cid:0) N T n T kB + B T k T nN (cid:1) dΓ + η h e Z e N T N dΓ , e ∈ Γ D , K R = Z e h N T N dΓ , e ∈ Γ D , (14)The global heat capacity and thermal conductivity matrices can be established by assemblingall the submatrices. The process is the same as the FEM. Similarly, the point and boundary heatflux vectors are developed: q E = Z E N T Q dΩ , E ∈ Ω , q D = − Z e B T k T n e u D dΓ + η h e Z e N T e u D dΓ , e ∈ Γ D , q N = Z e N T e q N dΓ , e ∈ Γ N , q R = Z e h N T e u R dΓ , e ∈ Γ R . (15)The global heat flux vector is assembled in the same way. Eventually, a set of discretized ODEs(Eqn. 13) with sparse and symmetric matrices are achieved.
3. Local Variational Iteration Method (LVIM)
Eqn. 13 can be rewritten as a system of standard first-order ODEs: ˙ u = g ( u , t ) = − C − Ku + C − q ( t ) , t ∈ [0 , T ] . (16)14he unknown temperature vector u = [ u , u , · · · , u L ] T , where L is the number of points used inthe FPM.In a finitely large time interval [ t i , t i +1 ] ⊂ [0 , T ] , with a given initial approximation u ( τ ) , theLocal Variational Iteration Method (LVIM) approximates the exact solution u at any time t with thefollowing correctional iterative formula [46]: u n +1 ( t ) = u n ( t ) + Z tt i λ ( τ ) R ( u n , τ ) d τ, (17)where the error residual R ( u n , τ ) is defined as: R ( u n , τ ) = ˙ u n ( τ ) − g ( u n , τ ) , (18) λ ( τ ) is a matrix of Lagrange multipliers which are yet to be determined.Eqn. 17 can also be regarded as a correctional iteration based on an optimally weighted errorresidual in time interval [ t i , t ] , where λ ( τ ) is the set of optimal weighting functions.By making the right side of Eqn. 17 stationary, we obtain the following constraints for λ ( τ ) : I + λ ( τ ) (cid:12)(cid:12)(cid:12) τ = t = ˙ λ ( τ ) = λ ( τ ) J ( u n , τ ) , τ ∈ [ t i , t ] , (19)where J ( u n , τ ) = ∂ g ( u n , τ ) ∂ u n = − C − K (20)is the Jacobian matrix. I is the unit matrix.Using the theory of Magnus series [50], it can be proved that [51]: λ ( t ) = − I ∂ λ ( t ) ∂t = − J ( u n , t ) λ ( t ) , t ∈ [ τ, t i +1 ] , (21)Differentiating Eqn. 17 and substituting Eqn. 21 into it: ˙ u n +1 ( t ) = ˙ u n ( t ) + λ ( t ) R ( u n , t ) + Z tt i ∂ λ ( t ) ∂t R ( u n , τ ) d τ = g ( u n , t ) − J ( u n , t ) Z tt i λ ( τ ) R ( u n , τ ) d τ = g ( u n , t ) − J ( u n , t ) [ u n +1 ( t ) − u n ( t )] . (22)15herefore, we obtain the recursive formula: ˙ u n +1 ( t ) + J ( u n , t ) u n +1 ( t ) = g ( u n , t ) + J ( u n , t ) u n ( t ) , t ∈ [ t i , t i +1 ] . (23) Eqn. 23 can be written in the weak-form in the time interval [ t i , t i +1 ] with a matrix of test func-tions v ( t ) : Z t i +1 t i v ( t ) [ ˙ u n +1 ( t ) + J ( u n , t ) u n +1 ( t )] d t = Z t i +1 t i v ( t ) [ g ( u n , t ) + J ( u n , t ) u n ( t )] d t. (24)Let v ( t ) = diag ([ v, v, · · · , v ]) , where v is the Dirac Delta function for a set of collocation nodes t , t , · · · , t M ∈ [ t i , t i +1 ] , that is: v = δ ( t − t m ) , t m ∈ [ t i , t i +1 ] , m = 1 , , · · · , M. (25)The weak-form formula leads to: ˙ u n +1 ( t m ) + J ( u n , t m ) u n +1 ( t m ) = g ( u n , t m ) + J ( u n , t m ) u n ( t m ) ,t m ∈ [ t i , t i +1 ] , m = 1 , , · · · , M. (26)A set of orthogonal basis functions Φ = { φ , φ , · · · , φ N } are used to construct the trial function u e : u e ( t ) = N X n =0 a e,n φ n ( t ) , (27)where u e ( t )( e = 1 , , · · · , L ) are elements of the solution vector u ( t ) . A number of types of basisfunctions can be used in the collocation method, including harmonics, polynomials, Radial BasisFunctions (RBFs), etc. In this paper, we employed the first kind of Chebyshev polynomials [52] asan example. The collocation nodes are selected as Chebyshev-Gauss-Lobatto points. From Eqn. 27,we can get: U e = QA e , ˙ U e = ( LQ ) A e (28)where U e = [ u e ( t ) , u e ( t ) , · · · , u e ( t M )] T , A e = [ a e, , a e, , · · · , a e,N ] T , Q = φ ( t ) φ ( t ) · · · φ N ( t ) φ ( t ) φ ( t ) · · · φ N ( t ) ... ... . . . ... φ ( t M ) φ ( t M ) · · · φ N ( t M ) ( N +1) × M , LQ = ˙ φ ( t ) ˙ φ ( t ) · · · ˙ φ N ( t )˙ φ ( t ) ˙ φ ( t ) · · · ˙ φ N ( t ) ... ... . . . ... ˙ φ ( t M ) ˙ φ ( t M ) · · · ˙ φ N ( t M ) ( N +1) × M . M = N + 1 . Thus, we achieve the relation between U e and its derivative: ˙ U e = ( LQ ) Q − U e . (29)Finally, substituting the relation into Eqn. 26 and rearranging the sequence of the collocationequations: (cid:16) e E + e J (cid:17) e U n +1 = (cid:16) e E + e J (cid:17) e U n − e R , (30)where e U = (cid:2) U T1 , U T2 , · · · , U T L (cid:3) T , e E = I L × L ⊗ (cid:2) ( LQ ) Q − (cid:3) , e J = J (cid:2) diag (ˆ t ) (cid:3) = (cid:0) − C T K (cid:1) ⊗ I M × M , ˆ t = [ t , t , · · · , t M ] T , e C = C ⊗ I M × M , e K = K ⊗ I M × M , e q = q (ˆ t ) , e R = e E e U n + e C − e K e U n − e C − e q , here ⊗ denotes the Kronecker product.The LVIM can usually achieve good estimates with very simple initial guess functions, e.g.,linear functions. In the current approach, we simply assume e U = u ( t i ) ⊗ [1 , , · · · , T , that is, aconstant function as the initial condition at all time steps. To apply the initial conditions, we usuallyselect the first collocation point at the initial boundary, i.e., t = t i . However, this would makeEqn. 30 overdetermined. To solve that problem, the collocation equations at the initial boundaryneed to be eliminated. Thus, the final iteration formula in LVIM can be written as: e U r n +1 = e U r n − (cid:16) e E r + e J r (cid:17) − e R r , (31)where [] r stands for the remained vector (or matrix) after eliminating the ( pM + 1 )th rows (andcolumns), p = 0 , , · · · , L − .There are some other modifications of LVIM in which the matrix of Lagrange multipliers λ ( τ ) is approximated in different ways [45]. Some of these modifications have potentials in further im-proving the computing efficiency by avoiding the inversion of the Jacobian matrix, especially forsystems dominated by a few eigenvalues. Yet in this paper, we just concentrate on the basic LVIMscheme shown in Eqn. 31. 17 . Numerical results and discussion In this section, a number of 2D and 3D numerical examples are carried out to illustrate theimplementation and effectiveness of our approach. Both steady-state and transient heat conductionproblems are presented. The FPM is employed for spatial discretization in all the examples witheither uniform or random points. The LVIM is applied in transient examples and the results arecompared with explicit and implicit Euler schemes. Anisotropic nonhomogeneous materials areconsidered. Some complex and practical examples are solved after that, followed by a discussionon the penalty parameters in the FPM and the number of collocation nodes in LVIM. The relativeerrors r and r used in this section are defined as: r = (cid:13)(cid:13) u h − u (cid:13)(cid:13) L k u k L , r = (cid:13)(cid:13) ∇ u h − ∇ u (cid:13)(cid:13) L k∇ u k L (32)where k u k L = (cid:18)Z Ω u dΩ (cid:19) / , k∇ u k L = (cid:18)Z Ω |∇ u | dΩ (cid:19) / . (33) In the first example (Ex. (1.1)), a circular isotropic and homogenous domain is considered. With-out loss of generality, we assume the material properties ρ = 1 , c = 1 , thermal conductivity tensorcomponents k = k = 1 , k = k = 0 . The body source density Q is absent. The simplifiedgoverning equation can be written as: ˙ u ( x, y, t ) = ∇ u ( x, y, t ) . (34)We consider a postulated analytical solution: u ( x, y, t ) = e x + y cos( x + y + 4 t ) , ( x, y ) ∈ (cid:8) ( x, y ) | x + y ≤ (cid:9) , (35)Dirichlet boundary conditions are prescribed on the circumference, corresponding to the given pos-tulated solution. A total of 601 points are distributed uniformly or randomly in the domain, 30 ofwhich are on the boundary ( x + y = 1 ). The Dirichlet boundary condition is applied directly.Hence the penalty parameter η is eliminated. The solutions based on the FPM + LVIM / BackwardEuler scheme and their relative errors at t = 0 . are presented in Fig. 6 and Table 1, in which η M is the number of collocation points in each timeinterval, and tol is the error tolerance in stopping criteria in the LVIM. As can be seen, the FPMcan be incorporated with different ODE solvers and achieve highly accurate solutions. Whereas theLVIM in the time domain reduces the computational cost significantly. Table 1: Relative errors and computational time of FPM + LVIM / backward Euler approach in solving Ex. (1.1).
Method Computationalparameters Time step Relative errors Computationaltime (s)FPM + LVIM(601 uniform points) η = 2 , M = 5 , tol = 10 − ∆ t = 0 . r = 6 . × − r = 1 . × − η = 2 , M = 5 , tol = 10 − ∆ t = 0 . r = 5 . × − r = 1 . × − η = 2 ∆ t = 0 . r = 7 . × − r = 1 . × − η = 2 ∆ t = 0 . r = 5 . × − r = 1 . × − u ( x, y, t ) = √ e − π t/ h cos( πx − π πy − π i , ( x, y ) ∈ { ( x, y ) | x ∈ [0 , , y ∈ [0 , } . (36)Neumann boundary condition consistent with the postulated solution is applied on x = 1 , while theother sides are under Dirichlet boundary conditions. 144 uniform or random points are utilized, ofwhich 44 points are on the boundaries. The computed solutions are shown in Table 2 and Fig. 7.Our current FPM + LVIM approach presents significantly high accuracy for the mixed boundaryvalue problem. The computational speed is ten times higher than the forward and backward Eulerschemes. While the forward Euler scheme may become unstable and result in divergent results witha large time step, the LVIM shows its reliability under relatively large time intervals.19 a) (b)Figure 6: Ex. (1.1) - The computed solution when t = 0 . . (a) 601 uniform points. (b) 601 random points.Table 2: Relative errors and computational time of FPM + LVIM / forward Euler / backward Euler approach in solvingEx. (1.2). Method Computationalparameters Time step Relative errors Computationaltime (s)FPM + LVIM(144 uniform points) η = 2 , M = 5 , tol = 10 − ∆ t = 0 . r = 5 . × − r = 1 . × − η = 2 , M = 5 , tol = 10 − ∆ t = 0 . r = 1 . × − r = 2 . × − η = 2 ∆ t = 1 × − r = 6 . × − r = 1 . × − η = 2 ∆ t = 1 × − r = 1 . × − r = 2 . × − η = 2 ∆ t = 5 × − r = 6 . × − r = 1 . × − η = 2 ∆ t = 1 × − r = 1 . × − r = 2 . × − a) (b)Figure 7: Ex. (1.2) - The computed solution when t = 1 . (a) 144 uniform points. (b) 144 random points. In the following four examples, a benchmark mixed boundary value problem in anisotropic non-homogeneous materials is considered. The tested domain is a L × L square with Dirichlet boundaryconditions on y = 0 and y = L . Symmetric boundary conditions are applied on the lateral sides. Forisotropic problems, symmetry is equivalent to Neumann boundary condition with e q N = 0 . Whereasfor anisotropic problems, an additional boundary thermal conductivity matrix has to be employed: K S = − Z e (cid:2) N T n T ( k − k I ) B (cid:3) dΓ , e ∈ Γ S , (37)where k is the first diagonal element of the thermal conductivity tensor k . Γ S stands for thesymmetric boundaries. Clearly, the matrix vanishes in isotropic domain. The initial, boundaryconditions and material properties are given as: u ( x, , t ) = u , u ( x, L, t ) = u L , u ( x, y,
0) = u ,ρ ( x, y ) = 1 , c ( x, y ) = f ( y ) , k ( x, y ) = f ( y ) ˆ k ˆ k ˆ k ˆ k , (38)where u , u L , ˆ k ij ( i, j = 1 , are constant. In isotropic case, ˆ k ij = δ ij . Whereas in anisotropiccase, ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 . The body source density Q = const = 0 . It turns out that theresulting temperature distribution is not dependent on x , i.e., the example can be equivalent to a 1Dheat conduction problem. 21n Ex. (1.3), u = 1 , u L = 20 , the gradation function f ( y ) = exp( δy/L ) . The exact solutionis obtained and given in [1]. When δ = 0 , the material is homogenous. The computed solution forisotropic homogenous, isotropic nonhomogeneous, and anisotropic nonhomogeneous materials arepresented and compared with exact solutions in Fig. 8. With only 121 ( × ) uniform points inthe domain, the result shows great agreement with the exact solution. It is also consistent with theresults shown in [1] based on meshless point interpolation method (PIM) and Laplace-transform(LT) approach. The time cost and average errors of the present FPM + LVIM approach is listed inTable 3, as well as the backward Euler scheme. The average error r is defined as the average valueof r in time interval [0 , . . It should be pointed out that in order to get a continuous solutionin the entire domain, the FPM with random points usually requires a larger penalty parameter η .Unfortunately, the accuracy drops down as η increases.As can be seen from Fig. 8 and Table 3, while the nonhomogeneity and anisotropy of the materialhave a significant influence on the temperature distribution, they do not give rise to any difficulties inthe present computing method. As the solution achieves steady state before t = 0 . , the advantageof LVIM approach in computational time is not distinct, especially when comparing with Ex. (1.1)and (1.2) in which the temperature solution varies violently. Yet the LVIM approach still savesapproximately one half of the computing time. t u ( L / , L / , t ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (a) y / L u ( x , y , . ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (b)Figure 8: Ex. (1.3) - The computed solution with different material properties. (a) transient temperature solution atthe midpoint of the domain in time scope [0 , . . (b) vertical temperature distribution when t = 0 . . In Ex. (1.4) âĂŞ (1.6), we consider the same initial boundary value problem as shown in Ex. (1.3).22 able 3: Relative errors and computational time of FPM + LVIM / backward Euler approach in solving Ex. (1.3).
Method Computationalparameters Time step Average errors Computationaltime (s)Homogenous isotropic ( δ = 0 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM(144 uniform points) η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 5 . × − η = 20 , M = 5 , tol = 10 − ∆ t = 0 . r = 4 . × − η = 10 ∆ t = 0 . r = 2 . × − η = 20 ∆ t = 0 . r = 2 . × − δ = 3 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM(144 uniform points) η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 7 . × − η = 20 , M = 5 , tol = 10 − ∆ t = 0 . r = 4 . × − η = 10 ∆ t = 0 . r = 7 . × − η = 20 ∆ t = 0 . r = 2 . × − δ = 3 ; ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 )FPM + LVIM(144 uniform points) η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 7 . × − η = 20 , M = 5 , tol = 10 − ∆ t = 0 . r = 4 . × − η = 10 ∆ t = 0 . r = 7 . × − η = 20 ∆ t = 0 . r = 4 . × − f ( y ) and boundary values are given as:Ex. (1.4): exponential : f ( y ) = [exp( δy/L ) + exp( − δy/L )] , δ = 2 , u = 1 , u L = 20; Ex. (1.5): trigonometric : f ( y ) = [cos( δy/L ) + 5sin( δy/L )] , δ = 2 , u = 0 , u L = 100; Ex. (1.6): power-law : f ( y ) = (1 + δy/L ) , δ = 3 , u = 1 , u L = 20; The computed solutions of these three examples are shown in Fig. 9, 10 and 11 respectively.121 uniform points are utilized. The results achieve great agreement with the analytical solutions,confirming that the nonhomogeneity and anisotropy do not cause any difficulties in the FPM + LVIMapproach. The corresponding relative errors and computational times are shown in Table 4, 5 and6. The LVIM approach cuts the computing time approximately by a half and does not cause anystability problems. All these results are consistent with the numerical example solutions in [1]. t u ( L / , L / , t ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (a) y / L u ( x , y , . ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (b)Figure 9: Ex. (1.4) - The computed solution with different material properties. (a) transient temperature solution intime scope [0 , . . (b) vertical temperature distribution when t = 0 . . In Ex. (1.6), we can also replace the symmetric boundary conditions by Neumann boundaryconditions with heat flux vanishing on the sides. In the anisotropic case, as a result, the temperaturevariation in x - direction is no longer constant. The computed 2D temperature distribution basedon 144 random points is shown in Fig. 12. In Fig. 12(a), 44 of the points are distributed on theboundaries, hence the Dirichlet boundary condition is imposed directly. Whereas in Fig. 12(b), nopoints are on the boundary. A collocation method based on integral terms on the boundaries isemployed to enforce the essential boundary conditions. That is, in Eqn. 12, penalty parameter η is24 t u ( L / , L / , t ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (a) y / L u ( x , y , . ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (b)Figure 10: Ex. (1.5) - The computed solution with different material properties. (a) transient temperature solution intime scope [0 , . . (b) vertical temperature distribution when t = 0 . . t u ( L / , L / , t ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (a) y / L u ( x , y , . ) homog.isotr. :Exacthomog.isotr. :FPMnonhom.isotr. :Exactnonhom.isotr. :FPMnonhom.anis. :Exactnonhom.anis. :FPM (b)Figure 11: Ex. (1.6) - The computed solution with different material properties. (a) transient temperature solution intime scope [0 , . . (b) vertical temperature distribution when t = 0 . . able 4: Relative errors and computational time of FPM + LVIM / backward Euler approach (121 uniform points) insolving Ex. (1.4). Method Computationalparameters Time step Average errors Computationaltime (s)Homogenous isotropic ( δ = 0 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 6 . × − η = 10 ∆ t = 0 . r = 5 . × − δ = 2 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 7 . × − η = 10 ∆ t = 0 . r = 5 . × − δ = 2 ; ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 8 . × − η = 10 ∆ t = 0 . r = 8 . × − Table 5: Relative errors and computational time of FPM + LVIM / backward Euler approach (121 uniform points) insolving Ex. (1.5).
Method Computationalparameters Time step Average errors Computationaltime (s)Homogenous isotropic ( δ = 0 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 6 . × − η = 10 ∆ t = 0 . r = 6 . × − δ = 2 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 2 . × − η = 10 ∆ t = 0 . r = 2 . × − δ = 2 ; ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 9 . × − η = 10 ∆ t = 0 . r = 1 . × − able 6: Relative errors and computational time of FPM + LVIM / backward Euler approach (121 uniform points) insolving Ex. (1.6). Method Computationalparameters Time step Average errors Computationaltime (s)Homogenous isotropic ( δ = 0 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 6 . × − η = 10 ∆ t = 0 . r = 5 . × − δ = 2 ; ˆ k = ˆ k = 1 , ˆ k = ˆ k = 0 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 1 . × − η = 10 ∆ t = 0 . r = 1 . × − δ = 2 ; ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 )FPM + LVIM η = 10 , M = 5 , tol = 10 − ∆ t = 0 . r = 7 . × − η = 10 ∆ t = 0 . r = 7 . × − Ex. (1.7) is still in a square domain. However, the material property is no longer continuous.As shown in Fig. 13(a), in the top half of the domain ( y >
50 m ), the medium is isotropic andhas a thermal conductivity k = 2W / (m ◦ C) , while in the bottom half ( y <
50 m ), the isotropicthermal conductivity is k = 1W / (m ◦ C) . An adiabatic crack emanates on the midline of the domain(
25 m < x <
75 m , y = 50 m ). Dirichlet boundary condition is applied on all the externalboundaries. On the bottom and lateral sides, e u D = 0 ◦ C , while on the top side, e u D = 100 ◦ C . Forsimplicity, only the steady-state solution is considered in this example.In FPM, the subdomain boundaries shared by two points on either side of the crack are regardedas external boundaries. In this example, Neumann (adiabatic) boundary condition is applied. Thecomputed steady-state temperature distribution is presented in Fig. 13. The result based on 100uniform points (Fig. 13(b)) shows a very good accuracy and is consistent with the numerical result27 a) (b)Figure 12: Ex. (1.6) - The computed solution with vanishing heat fluxes on the lateral sides. (a) 44 points on theboundaries, η = 10 . (b) no points on the boundaries, η = 10 , η = 20 . given in [53], since the crack is just on top of some subdomain boundaries. However, in a partitionwith random distributed points, the crack may not coincide with the internal boundaries (as can beseen in Fig. 13(a)). Yet the FPM can still get a good approximation of the temperature distribution inthe entire domain, especially outside the vicinity of the crack. When the number of points increases(see Fig. 13(d)), the computed result approaches the exact solution gradually.Such a result shows the potential of the FPM in solving thermal-shock problems with crackpropagation in brittle materials. Without knowing the exact geometry of the cracks, an approximatesolution can be obtained by simply shifting some internal subdomain boundaries from Γ h to ∂ Ω inwhere the thermal stress is above the yield stress. Other than the adiabatic crack, multiple thermalcrack models can be incorporated with the FPM.In Ex. (1.8), a L-shaped orthotropic domain is considered. As shown in Fig. 14(a), the tem-perature is fixed to ◦ C on the left and bottom sides. The other sides are Neumann boundarieswith e q N = 0 on the black sides and e q N = 12W / m on the red sides. The orthotropic material hasthermal conductivity coefficients k = 4 W / m ◦ C , k = 7 W / m ◦ C , and k = k = 0 . Dirichletboundary conditions are applied by the collocation method. The penalty parameters η = 5 √ k k , η = 20 √ k k . The steady-state results are shown in Fig. 14. The FPM solution agrees well withthe FEM solution achieved by ABAQUS with 310 linear quadrilateral elements (341 nodes).In the last 2D example (Ex. (1.9)), we consider the transient heat conduction in a semi-infinite28 .5 m0.5 m ° C100 ° Ck = 2 W/(m ° C)k = 2 W/(m ° C) (a) x (m) y ( m ) (b)(c) (d)Figure 13: Ex. (1.7) âĂŞ The adiabatic crack and computed solutions ( η = 5 √ k k , η = 20 √ k k ). (a) the points,partition and adiabatic crack. (b) 100 uniform points. (c) 100 random points. (d) 2000 random points. m2 m2 m1 m ° C (a) NT11+6.623e+00+6.905e+00+7.186e+00+7.468e+00+7.749e+00+8.030e+00+8.312e+00+8.593e+00+8.874e+00+9.156e+00+9.437e+00+9.719e+00+1.000e+01Step: Step−1Increment 1: Step Time = 1.000Primary Var: NT11ODB: EX8−02.odb Abaqus/Standard Student Edition 2018 Tue Dec 31 08:19:08 Central Standard Time 2019
XYZ (b)(c) (d)Figure 14: Ex. (1.8) - The boundary conditions and computed solutions. (a) the problem domain and boundaryconditions. (b) ABAQUS solution with 310 DC2D4 elements (341 nodes) (c) FPM solution with 48 uniform points.(d) FPM solution with 100 random points.
12 m × domain is considered. According tothe symmetry, we only compute one half of the domain. As shown in Fig. 15(a), the pipe wall witha radius of .
45 m is modeled as a Dirichlet boundary with e u D = 20 ◦ C . The infinite boundary isapplied as e u D = 10 ◦ C on the right and bottom sides ( x = 6 m and y = − ). The left side ( x = 0 )is symmetric, and the top side ( y = 0 ) is adiabatic. The two boundary conditions are equivalent heresince the material is isotropic. A number of points are distributed in the domain. As the variation oftemperature is more violent, more points are distributed in the vicinity of the pipe wall. Generally,the number of points in a unit area decreases exponentially with the distance from the pipe center.The material properties are given as: ρ = 2620 kg / m , c = 900 J / kg ◦ C , k = 2 .
92 W / m ◦ C . Theinitial condition is u ( x, y,
0) = const = 10 ◦ C .In FPM, the essential boundary conditions are imposed by the collocation method. The penaltyparameters are set as: η = 5 k , η = 20 k . The computed time-variation of temperature at fourrepresentative points on the adiabatic side is shown in Fig. 15(b). The results present great consis-tency with the FEM result achieved by ABAQUS with 661 DC2D4 elements (715 nodes) and anexplicit solver. The number of time steps is 100 in ABAQUS and 10 in the LVIM approach. Thetemperature distribution results at t = In this section, we consider a number of 3D heat conduction examples in a cubic domain
Ω = { ( x, y, z ) | x, y, z ∈ [0 , L ] } . Various boundary conditions and material properties are tested. Theheat source density Q vanishes in all the following examples.First, a steady-state problem with homogenous anisotropic material is considered. The thermalconductivity tensor components k = k = k = 1 × − , k = 0 . × − , k = k = 0 . Apostulated analytical solution is considered: u ( x, t, z ) = y + y − yz + xz. (39)31 m 10 ° C(infinite soil)8 msym1.27 mr=0.45 m20 ° C(pipe) (a) t (hours) u ( x , , t ) ( ° C ) FPM, LVIM: x = 0.75 mFPM, LVIM: x = 1.5 mFPM, LVIM: x = 3.0 mFPM, LVIM: x = 4.5 m ABAQUS: x = 0.75 mABAQUS: x = 1.5 mABAQUS: x = 3.0 mABAQUS: x = 4.5 m (b)(c)(d)Figure 15: Ex. (1.9) - The problem and computed solutions. (a) the problem domain and boundary conditions. (b)transient temperature solution. (c) temperature distribution when t = t = z = 0 . L is shown in Fig. 16. In Fig. 16(a), × × points are distributed uniformlyin the cube, while in Fig. 16(b) the points are distributed randomly. The penalty parameters are: η = 5 k for the uniform points, and η = 10 k , η = 20 k for the random points. Both resultsmatch well with the exact solution. The relative errors r are . × − and . × − respectively. y / L -120-100-80-60-40-200204060 u ( x , y , . ) x/L=0.5: Exactx/L=0.9: Exactx/L=0.5: FPMx/L=0.9: FPM (a) y / L -120-100-80-60-40-200204060 u ( x , y , . ) x/L=0.5: Exactx/L=0.9: Exactx/L=0.5: FPMx/L=0.9: FPM (b)Figure 16: Ex. (2.1) - The computed solutions at z = 0 . L . (a) 1000 uniform points. (b) 1000 random points. Next, a transient heat conduction example is considered. The material properties are given as: ρ = 1 , c = 1 , and k = k = k = 1 , k = k = k = 0 . The boundary condition on the topsurface ( z = L ) is prescribed as a thermal shock e u D = H ( t − , where H is the Heaviside timestep function. The bottom boundary condition on z = 0 is given as e u D = 0 . And all the lateralsurfaces ( x, y = 0 , L ) have vanishing heat fluxes. The initial condition is u ( x, y, z,
0) = 0 . Theside length L = 10 . It turns out that the temperature distribution in this example is not dependenton x and y coordinates. As a result, the problem can be analyzed equivalently in 2D. The transienttemperatures at z = 0 . L , z = 0 . L and z = 0 . L are computed by the 2D and 3D FPM andpresented in Fig. 17 respectively. The computational times cost by the LVIM approach and thebackward Euler scheme are shown in Table 7. As the time-variation of temperature is smooth inthis case, the LVIM approach can only improve the computing efficiency slightly.In Ex. (2.3), we consider a similar initial boundary condition problem as Ex. (2.2) in an isotropicmedium. The thermal conductivity tensor components are: k = k = 1 , k = 1 . , k = 0 . ,33 = k = 0 . Symmetric boundary conditions are given on the left and right surfaces ( x = 0 , L )instead of the Neumann boundary conditions. The temperature distribution is then independent of x coordinate, i.e., the example can also be equivalent to a 2D problem. Fig. 18 compares the computedsteady-state temperature distribution on y = 0 and y = L analyzed by 2D and 3D FPM ( η = 10 k , η = 20 k in both cases). Very good agreement can be observed. The transient result is shown asa comparison of Ex. (2.4) in the following Fig. 19(a). t u ( x , y , z , t ) z/L= 0.1(3D: 1000 points)z/L= 0.5(3D: 1000 points)z/L= 0.8(3D: 1000 points)z/L= 0.1(2D: 100 points)z/L= 0.8(2D: 100 points)z/L= 0.5(2D: 100 points) Figure 17: Ex. (2.2) - The computed transienttemperature solution. z / L u ( x , y , z ) y/L=0: FPM (2D:100 points)y/L=1: FPM (2D:100 points)y/L=0: FPM (3D:512 points)y/L=1: FPM (3D:512 points) Figure 18: Ex. (2.3) - The computed steady-state result.Table 7: Computational time of 3D FPM + LVIM / backward Euler approach (1000 points) in solving Ex. (2.2).
Method Computationalparameters Time step Computationaltime (s)FPM + LVIM η = 10 , η = 20 , M = 3 , tol = 10 − ∆ t = 5 . η = 10 , η = 20 ∆ t = 0 . ρ and heat capacity c remain constant in the whole domain.Whereas the thermal conductivity tensor is prescribed as: k ( z ) = 1 + z/L , k = 1 , k = 1 . , k = 0 . , k = k = 0 . The side length L = 10 . The example can also be analyzed in 2D.A comparison of the transient 2D and 3D computed temperatures on z = 0 . L is presented inFig. 19(a). The homogenous result (Ex. (2.3)) is also shown as a comparison. Table 8 shows the34omputational times for the LVIM approach and backward Euler scheme when achieving the sameaccuracy. As can be seen, the nonhomogeneity has a considerable influence on the temperature dis-tribution, while it has no influence on the accuracy or efficiency of the FPM + LVIM approach. Thetransient temperature solution approaches the steady-state result, as shown in Fig. 19(b), gradually. t u ( x , y , z , t ) y/L= 1 ; z/L= 0.2(3D: 1331 points)y/L= 0 ; z/L= 0.2(3D: 1331 points)y/L= 1 ; z/L= 0.2(2D: 121 points)y/L= 0 ; z/L= 0.2(2D: 121 points)y/L= 1 ; z/L= 0.2(homogeneous)y/L= 0 ; z/L= 0.2(homogeneous) (a) z / L u ( x , y , z ) y/L=0: FPM (2D:100 points)y/L=1: FPM (2D:100 points)y/L=1: FPM (3D:512 points)y/L=0: FPM (3D:512 points)y/L=0: homogeneousy/L=1: homogeneous (b)Figure 19: Ex. (2.4) - The computed solutions. (a) transient temperature solution. (b) steady-state result.Table 8: Computational time of 3D FPM + LVIM / backward Euler approach (1331 points) in solving Ex. (2.4). Method Computationalparameters Time step Computationaltime (s)FPM + LVIM η = 10 , η = 20 , M = 4 , tol = 10 − ∆ t = 25 η = 10 , η = 20 ∆ t = 0 . L × L × L cube with vanishing flux on all the lateral surfaces. The boundaryconditions on the top and bottom surfaces are given as: e u D = H ( t − , for z = L ; and e u D =0 , for z = 0 . The homogenous anisotropic thermal conductivity coefficients: k = k = 1 , k = 1 . , k = k = k = 0 . . The other conditions are the same as the previous examples.The computed solution is compared with FEM result achieved by ABAQUS with 1000 linear heattransfer elements (DC3D8) and shown in Fig. 20(a). The homogenous solution (Ex. (2.2)) is alsoshown for comparison. As can be seen, a good agreement is observed between the FPM + LVIM35nd ABAQUS results. As time goes on, the transient solution keeps approaching the steady state.The computed temperature distributions on the four lateral sides of the domain ( x, y = 0 , L ) areshown in Fig. 20(b), as well as the ABAQUS results. It is clear that the solution is dependent onall x , y and z coordinates, and cannot be simplified as a 2D problem. The penalty parameters andcomputational times are listed in Table 9, confirming that the FPM + LVIM approach can work withconsiderable large time intervals and achieving accurate solutions. t u ( x , . , . , t ) FEM: anisotropic: x/L= 0FEM: anisotropic: x/L= 1FEM: isotropicABAQUS: anisotropic: x/L= 0ABAQUS: anisotropic: x/L= 1ABAQUS: isotropic (a) z / L u ( x , y , z ) ABAQUS: x/L=1; y/L=1ABAQUS: x/L=1; y/L=0ABAQUS: x/L=0; y/L=1ABAQUS: x/L=0; y/L=0FPM: x/L=1; y/L=1FPM: x/L=1; y/L=0FPM: x/L=0; y/L=1FPM: x/L=0; y/L=0 (b)Figure 20: Ex. (2.5) - The computed solutions. (a) transient temperature solution. (b) steady-state result.Table 9: Computational time of 3D FPM + LVIM / backward Euler approach (1000 points) in solving Ex. (2.5).
Method Computationalparameters Time step Computationaltime (s)FPM + LVIM η = 5 , η = 10 , M = 5 , tol = 10 − ∆ t = 25 η = 5 , η = 10 ∆ t = 0 . k ( z ) = 1 + z/L , k = 1 , k = 1 . , k = k = k = 0 . . All the boundary conditions are the same as Ex. (2.5). Fig. 21 presents thecomputed steady-state solution obtained by the FPM. The solution, as well as all the previous solu-tions in Ex. (2.1) âĂŞ Ex. (2.5), are consistent with the computed solutions achieved by Sladek et al.[4] with the Meshless Local Petrov-Galerkin (MLPG) method and Laplace-transform technique.36n Ex. (2.7), a transient heat conduction example with Robin boundary condition is tested. Thematerial is homogenous and isotropic: ρ = 1 , c = 1 , k = k = k = 1 , k = k = k = 0 .The top surface has a heat transfer coefficient h = 1 . . And the temperature outside the top surfaceis prescribed as e u R = H ( t − . All the lateral surfaces and bottom surface have heat fluxes e q N = 0 .Started from an initial condition u ( x, y, z,
0) = 0 , the temperature distribution depends only on z coordinate and the time. The analytical solution can be written as [4]: u ( x, y, z, t ) = u ( z, t ) = 1 − m ∞ X i =1 sin β i cos (cid:0) β i zL (cid:1) exp (cid:16) − β i k tρcL (cid:17) β i (cid:0) m + sin β i (cid:1) , where β i are roots of the transcendental equation: β sin β cos β − m = 0 , where m = hLk . (40)Let L = 10 . The computed time-variations of temperature on the bottom and midsurface of thecube ( z = 0 , . L ) are shown in Fig. 22, in which an excellent agreement is observed between theFPM + LVIM solution and the analytical result. z / L u ( , y , z ) y/L=0: FPM (3D:1000 points)y/L=1: FPM (3D:1000 points)y/L=0: homogeneousy/L=1: homogeneous Figure 21: Ex. (2.6) - The computed steady-state result. t u ( x , y , z , t ) z/L= 0: FPM (3D:1000 nodes)z/L= 0.5: FPM (3D:1000 nodes)z/L= 0: Exactz/L= 0.5: Exact Figure 22: Ex. (2.7) - The computed transienttemperature solution.
Finally, two practical examples with multiple materials and complicated geometries are consid-ered. Ex. (2.8) shows the heat conduction in a wall with crossed U-girders. The example is presentedin [54]. As shown in Fig. 23(a), the wall is consisted of two gypsum wallboards, two steel crossed37-girders and insulation materials (the insulation material is not presented in the sketch). The U-girders are separated by
300 mm . Thus, we can only focus on a
300 mm ×
300 mm ×
262 mm cellof the wall. The material properties are listed in Table 10. The indoor ( z = 262 mm ) and outdoor( z = 0 mm ) surfaces are under convection boundary conditions. The corresponding heat transfercoefficients h and the temperatures outside the surfaces are shown in Table 11. All the other lateralsurfaces are symmetric, i.e., e q N = 0 in this case. The initial condition is ◦ C in the whole domain.A total of 2880 points are used in the FPM analysis. Notice that though the insulation materialis not shown in the sketch, there are still points distributed in it. As a result of the uneven varia-tion of material properties, the density of points used in the gypsum and steel are higher than theinsulation. It should be pointed out that when the point distribution is extremely uneven, as in thisexample, it is highly recommended to define the boundary-dependent parameter h e in Eqn. 12 as thedistance of the two points sharing the subdomain boundary. Fig. 23(b) presents the time-variationof temperatures on three representative points A, B, and C (shown in Fig. 23(a)) in 10 hours. TheFPM + LVIM solution shows an excellent consistency with the ABAQUS result obtained with 9702DC3D8 elements (11132 nodes). The computed temperature distribution in the gypsum wallboardsand U-girders when t = 1 hours and t ≥
10 hours (steady-state) are exhibited in Fig. 23(c) and23(d). The results also agree well with ABAQUS. The computational parameters and times areshown in Table 12. As can be seen, the LVIM approach helps to save approximately one half ofthe total computing time. Ex. (2.8) demonstrates the accuracy and efficiency of the FPM + LVIMapproach in solving complicated 3D transient heat conduction problems with multiple materials andhighly uneven point distributions.
Table 10: Material properties in Ex. (2.8).
Material ρ (( kg/m )) c ( × J / (kg ◦ C)) k (W / (m ◦ C)) gypsum 2300 1.09 0.22steel 7800 0.50 60insulation 1.29 1.01 0.036
Table 11: Boundary conditionsin Ex. (2.8). bc e u R ( ◦ C) h (W / (m ◦ C)) outdoor 20 25indoor 30 7.7Ex. (2.9) is also given in [54]. In this example, the heat transfer through a wall corner is studied.38 a) t (hours) u ( ° C ) FPM: AFPM: BFPM: CABAQUS: AABAQUS: BABAQUS: C (b)(c) (d)Figure 23: Ex. (2.8) âĂŞ The problem and computed solutions. (a) the problem domain and boundary conditions. (b)transient temperature solution. (c) temperature distribution when t = 1 hours. (d) steady-state result.Table 12: Computational time of FPM + LVIM / backward Euler approach (2880 points) in solving Ex. (2.8). Method Computationalparameters Time step Computationaltime (s)FPM + LVIM η = 11 , M = 3 , tol = 10 − ∆ t = 2 h η = 11 ∆ t = 0 . h δ stands for adiabatic boundaries, while α , β and γ are allconvection boundaries. Table 14 presents their heat transfer coefficients and surface temperatures.The initial condition is ◦ C in the whole domain.First, we concentrate on the temperatures of four representative points (A, B, C, D) as shownin Fig. 24(a). The time-variation of temperatures on these points is presented in Fig. 24(b). Theresult approaches steady state as time increases. Table 15 illustrates how the number of pointsused in the FPM influences the steady-state solution. The results are compared with data in theEuropean standards (CEN, 1995) [54]. As can be seen, when the number of points rises, the solutionapproaches the reference solution gradually. With more than 6288 points, the result keeps stable andhas no more than . ◦ C error compared with the CEN solution. Fig. 24(c) and Fig. 24(d) presentthe temperature distribution in the corner when t = 6 hours and after steady-state. Table 15 showsthe computational parameters and times comparing with the backward Euler scheme. Similar withthe previous examples, the LVIM approach works well with large time intervals and has no stabilityproblem. Table 13: Material properties in Ex. (2.9).
Material ρ (( kg/m )) c ( × J / (kg ◦ C)) k (W / (m ◦ C))
M1 849 0.9 0.7M2 80 0.84 0.04M3 2000 0.8 1.0M4 2711 0.88 2.5M5 2400 0.96 1.0
Table 14: Boundary conditionsin Ex. (2.9). bc e u R ( ◦ C) h (W / (m ◦ C)) α
20 5 β
15 5 γ δ – 0 (adiabatic) As have been stated in the previous sections, the penalty parameters η and η have a significantinfluence on the accuracy and stability of the FPM. For example, if η is too small, the method40 a) t (hours) u ( ° C ) ABCD (b)(c) (d)Figure 24: Ex. (2.9) âĂŞ The problem and computed solutions. (a) the problem domain and boundary conditions. (b)transient temperature solution. (c) temperature distribution when t = 6 hours hours. (d) steady-state result.Table 15: The computed steady-state temperatures ( ◦ C ) obtained by the FPM with different numbers of points ( L ) -Ex. (2.9). Point L = 1120 L = 3006 L = 6288 L = 11350 L = 28350 CEN [54]A 12.7 12.8 12.7 12.6 12.6 12.6B 10.9 11.1 11.1 11.0 11.0 11.1C 12.7 14.6 15.1 15.2 15.2 15.3D 15.9 16.4 16.5 16.4 16.4 16.441 able 16: Computational time of FPM + LVIM / backward Euler approach (3006 points) in solving Ex. (2.9).
Method Computationalparameters Time step Computationaltime (s)FPM + LVIM η = 10 , M = 3 , tol = 10 − ∆ t = 30 h η = 10 ∆ t = 1 . h η is too small, the Dirichlet boundaryconditions may not be satisfied. On the contrary, if η is very large, small jumps of temperatureon the internal boundaries can be expected, but the accuracy of the solution is doubtable. In thissection, parametric studies on η and η are carried out on both 2D and 3D examples.First, the steady-state solution of 2D example Ex. (2.3) is considered. We concentrate on thenonhomogeneous anisotropic case, i.e., δ = 3 , ˆ k = ˆ k = 2 , ˆ k = ˆ k = 1 . A total of 225 pointsare used in the FPM. Fig. 25(a) shows the influence of the penalty parameters on the relative errors r and r . The penalty parameters are nondimensionalized by k = (cid:16) ˆ k ˆ k (cid:17) / = 2 . As can beseen, in order to get a stable and accurate solution, it is recommended to define η in the range of . k to k , and η larger than k . The best choice in this example is η = 5 k and η > k .The accuracy decreases dramatically when η is too large or η is too small. Yet there is no upperlimit of the recommended range of η . Notice that in homogenous or isotropic case, the effectiverange of η and η can be much larger.In 3D case, the anisotropic steady-state example Ex. (2.1) is considered. With 1000 points dis-tributed uniformly in the domain, the relative errors r and r of the computed FPM solution undervarying η and η are shown in Fig. 25(b), in which the penalty parameters are nondimensionalizedby k = ( k k k ) / = 1 × − . To get a continuous and accurate computed solution, the penaltyparameters should be defined in the range . k < η < k , and k < η < k . In this ex-ample, the best choice is η = 15 k and η = 5000 k . However, the best choice varies under differentpoint distributions. As can be seen from the parametric studies, the relative errors shoot up when η or η is too small, as the continuity or essential boundary conditions may not be satisfied then.An excessively large η should also be avoided. Whereas a large η is still acceptable since it doesnot do much harm to the accuracy. 42n general, the recommended values of η and η are proportional to the thermal conductivity k .The approximate effective ranges of η and η are k < η < k , and k < η < × k , where k = ( k k ) / in 2D case and k = ( k k k ) / in 3D case. Notice that the range may vary underdifferent definitions of h e and different point distributions. Generally, the best choice of η shouldbe slightly larger than η since a small discontinuity of temperature on the internal boundaries isacceptable, while the essential boundary conditions should be satisfied strictly. Homogenous andisotropic problem usually has less requirement on the effective penalty parameters. Next, the recommended value of the number of collocation nodes in each time interval ( M ) in theLVIM is discussed. Take Ex. (1.1) as an example, Fig. 26 presents the relationship of computationaltime and average relative error r achieved by the LVIM approach with different M and backwardEuler scheme. The example is discretized with 601 uniform points in the FPM with η = 5 . Andthe error r is defined as the average value of relative error between the computed solution and theconverged solution (obtained with an extremely small time step) in time scope [0 , . As can beseen, though the backward Euler scheme and LVIM approach with M = 3 have an advantage incomputational time under low accuracy requirement, e.g., r ≤ . , their computational timesincrease rapidly when the required relative error decreases. As a result, large M has a benefitin achieving relatively accurate solution, while small M is more suitable for exploring a roughapproximation. Notice that the backward Euler scheme is equivalent to the LVIM approach with M = 2 , and follows the same tendency of accuracy and computational time for the LVIM.Table 17 shows the computational times required for the LVIM approach and backward Eulerscheme ( M = 2 ) when obtaining the same relative errors. To get a solution with r ≈ × − , theLVIM approach with M = 3 costs the least computational time, which is approximately one third ofthe computational time of the backward Euler scheme. On the other hand, in order to achieve higheraccuracy, e.g. r ≈ × − , the best choice would become M = 5 . Comparing with the backwardEuler scheme, the LVIM approach shows extraordinary efficiency under high accuracy requirement.Since the computational time rises rapidly with M , too many collocation nodes (e.g., M > ) arenot recommended. We usually apply M in the range of 3 to 5. For problems with lower accuracyrequirement and higher numbers of nodes, a small M is recommended. Whereas for problems with43 -5 / k -3 -2 -1 R e l a t i v e E rr o r s / k = 1000 r r -2 / k -2 -1 R e l a t i v e E rr o r s / k = 200 r r / k -3 -2 -1 R e l a t i v e E rr o r s / k = 5 r r (a) / k -2 -1 R e l a t i v e E rr o r s / k = 20 r r (b)Figure 25: Parametric studies on η and η . (a) 2D case: Ex. (1.3). (b) 3D case: Ex. (2.1). M could be more beneficial. Computational time (s) -4 -3 -2 -1 A v e r age E rr o r Backward EulerLVIM: M=3LVIM: M=4LVIM: M=5LVIM: M=6LVIM: M=7
Figure 26: Parametric study on the number of collocation nodes in each time interval ( M ) - Ex. (1.1).
5. Conclusion
A new computational approach is developed for analyzing 2D and 3D transient heat conduc-tion problems in complex anisotropic nonhomogeneous media. The truly meshless Fragile PointsMethod (FPM) based on Galerkin weak-form formulation is employed for spatial discretization,while the Local Variational Iteration (LVI) scheme is used to achieve the solution in the time do-main. The meshless FPM is a significant advancement over either the Element-Free Galerkin (EFG)Method or the Meshless Local Petrov-Galerkin (MLPG) Method. The EFG is based on GlobalGalerkin weak-form and requires back-ground cells to integrate the weak-form terms. The integra-tion becomes tedious while using the meshless Moving Least Squares (MLS) approximations. Alsowhen the mesh of back-ground cells is rotated, the EFG may not be an objective method. The MLPGis a truly meshless method, based on a local Petrov-Galerkin weak-form, and the integration of theweak-form is complicated when MLS approximations are used as trial functions and test functionsare different from the trial functions. The FPM is also a truly meshless method based on a Galerkinweak-form, uses very simple polynomial discontinuous trial and test functions and the integrationof the weak-form is simple. The imposition of essential boundary conditions in the FPM is simi-45 able 17: Computational time of FPM + LVIM / backward Euler approach under varying number of collocation nodesin each interval ( M ) in solving Ex. (1.1). Method M Time step ∆ t Computational time (s)Average relative error r ≈ × − FPM + backward Euler 2 0.013 16FPM + LVIM 3 0.27 54 0.53 75 0.80 86 1.33 127 1.60 16Average relative error r ≈ × − FPM + backward Euler 2 0.0016 144FPM + LVIM 3 0.08 174 0.27 135 0.53 126 0.8 207 1.1 2446ar to that in EFG and MLPG. The FPM leads to sparse symmetric matrices. The time integrationscheme LVIM is considerablely superior to the finite difference methods. Thus, the FPM + LVIMmethod for transient heat conduction in anisotropic nonhomogeneous solids presented in this paperis a superior meshless method as compared to those in earlier literatures.The FPM is generated by local, simple, polynomial, point-based (as opposed to element-basedin the FEM) and piecewise-continuous trial and test functions. Numerical Flux Corrections areemployed in terms of internal penalty functions. With large enough penalty parameters, the methodpresents its consistency and accuracy with both regularly and randomly distributed points. A simpledomain partition is still required, but just for integral computation. Symmetric and sparse matricescan be achieved in most cases. In the time domain, the highly efficient LVIM is introduced. As acombination of the VIM and a collocation method in each time interval, the LVIM shows excellentaccuracy and efficiency in solving nonlinear ODEs.Plenty of numerical examples are presented both in 2D and 3D. Mixed boundary conditions areinvolved, including Dirichlet, Neumann, Robin, and purely symmetric boundary conditions. Bothfunctionally graded materials and composite materials are considered. The computed solutions arecompared with analytical results, equivalent 1D or 2D results, and FEM solutions obtained by acommercial software. The forward and backward Euler schemes are used together with the FPM asa comparison to the LVIM. The FPM + LVIM approach exhibits great accuracy and efficiency andhas no stability problem under relatively large time intervals. The anisotropy and nonhomogeneitygive rise to no difficulties in the current approach. The computing efficiency is extraordinary whenthe response varies dramatically, or a high accuracy is required. The approach is also capable of an-alyzing systems with preexisting cracks, even if the domain partition does not coincide on the crackgeometry. This implies the further potential of the FPM + LVIM approach in solving crack propa-gation problems. At last, a recommended range of the computational parameters is given. We canconclude that, with suitable computational parameters, the FPM + LVIM approach shows excellentperformance in analyzing transient heat conduction systems with anisotropy and nonhomogeneity.
References [1] V. Sladek, J. Sladek, M. Tanaka, and C. Zhang. Transient heat conduction in anisotropic andfunctionally graded media by local integral equations.
Engineering Analysis with Boundary lements , 29(11):1047–1065, 2005.[2] K. J. Quint, S. Hartmann, S. Rothe, N. Saba, and K. Steinhoff. Experimental validation ofhigh-order time integration for non-linear heat transfer problems. Computational Mechanics ,48(1):81–96, 2011.[3] J. Zhang and S. Chauhan. Fast explicit dynamics finite element algorithm for transient heattransfer.
International Journal of Thermal Sciences , 139:160–175, 2019.[4] J. Sladek, V. Sladek, C. L. Tan, and S. N. Atluri. Analysis of transient heat conduction in 3Danisotropic functionally graded solids, by the MLPG method.
CMES - Computer Modeling inEngineering and Sciences , 32(3):161–174, 2008.[5] Y. Miyamoto, W. A. Kaysser, B. H. Rabin, A. Kawasaki, and R. G. Ford.
Functionally gradedmaterials: design, processing and applications , volume 5. Springer Science & Business Me-dia, 2013.[6] M. Simsek. Static analysis of a functionally graded beam under a uniformly distributed loadby Ritz method.
International Journal of Engineering & Applied Sciences , 1(3):1–11, 2009.[7] L. Chen, W. Lengauer, P. Ettmayer, K. Dreyer, H. W. Daub, and D. Kassel. Fundamentalsof liquid phase sintering for modern cermets and functionally graded cemented carbonitrides(FGCC).
International Journal of Refractory Metals and Hard Materials , 18(6):307–322,2000.[8] R. J. LeVeque and R. J. Leveque.
Numerical methods for conservation laws , volume 132.Springer, 1992.[9] M. Afrasiabi, M. Roethlin, and K. Wegener. Contemporary Meshfree Methods for ThreeDimensional Heat Conduction Problems.
Archives of Computational Methods in Engineering ,pages 1–35, 2019.[10] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu.
The finite element method: its basis andfundamentals . Elsevier, 2005. 4811] B. Shen, A. J. Shih, and G. Xiao. A heat transfer model based on finite difference method forgrinding.
Journal of Manufacturing Science and Engineering , 133(3):31001, 2011.[12] J. C. Chai, H. S. Lee, and S. V. Patankar. Finite volume method for radiation heat transfer.
Journal of thermophysics and heat transfer , 8(3):419–425, 1994.[13] L. C. Wrobel and A. J. Kassab. Boundary element method, volume 1: Applications in thermo-fluids and acoustics.
Appl. Mech. Rev. , 56(2):B17–B17, 2003.[14] J. Wen and M. M. Khonsari. Transient heat conduction in rolling/sliding components by adual reciprocity boundary element method.
International Journal of Heat and Mass Transfer ,52(5-6):1600–1607, 2009.[15] P. W. Randles and L. D. Libersky. Smoothed particle hydrodynamics: Some recent improve-ments and applications.
Computer Methods in Applied Mechanics and Engineering , 139(1-4):375–408, 1996.[16] W. K. Liu, S. Jun, and Y. F. Zhang. Reproducing kernel particle methods.
International journalfor numerical methods in fluids , 20(8-9):1081–1106, 1995.[17] J. K. Chen, J. E. Beraun, and T. C. Carney. A corrective smoothed particle method for boundaryvalue problems in heat conduction.
International Journal for Numerical Methods in Engineer-ing , 46(2):231–252, 1999.[18] D. I. Graham and J. P. Hughes. Accuracy of SPH viscous flow models.
International journalfor numerical methods in fluids , 56(8):1261–1269, 2008.[19] B. Nayroles, G. Touzot, and P. Villon. Generalizing the finite element method: diffuse approx-imation and diffuse elements.
Computational mechanics , 10(5):307–318, 1992.[20] Y. Krongauz and T. Belytschko. A Petrov-Galerkin diffuse element method (PG DEM) and itscomparison to EFG.
Computational Mechanics , 19(4):327–333, 1997.[21] T. Belytschko, Y. Y. Lu, and L. Gu. Element-free Galerkin methods.
International journal fornumerical methods in engineering , 37(2):229–256, 1994.4922] T. Zhu and S. N. Atluri. A modified collocation method and a penalty formulation for enforc-ing the essential boundary conditions in the element free Galerkin method.
ComputationalMechanics , 21(3):211–222, 1998.[23] T. Zhu, J. Zhang, and S. N. Atluri. Meshless numerical method based on the local boundaryintegral equation (LBIE) to solve linear and non-linear boundary value problems.
EngineeringAnalysis with Boundary Elements , 23(5):375–389, 1999.[24] S. N. Atluri and T. Zhu. A new Meshless Local Petrov-Galerkin (MLPG) approach in compu-tational mechanics.
Computational Mechanics , 22(2):117–127, 1998.[25] M. Shibahara and S. N. Atluri. The meshless local Petrov-Galerkin method for the analysis ofheat conduction due to a moving heat source, in welding.
International Journal of ThermalSciences , 50(6):984–992, 2011.[26] L. Dong, T. Yang, K. Wang, and S. N. Atluri. A new Fragile Points Method (FPM) in computa-tional mechanics, based on the concepts of Point Stiffnesses and Numerical Flux Corrections,Engineering Analysis with Boundary Elements.
Engineering Analysis with Boundary Ele-ments , 107:124–133, 2019.[27] D. N. Arnold, F. Brezzi, B. Cockburn, and L. Donatella Marini. Unified analysis of discon-tinuous Galerkin methods for elliptic problems.
SIAM Journal on Numerical Analysis , 39(5):1749–1779, 2001.[28] I. Mozolevski, E. S¨uli, and P. R. B¨osing. hp-version a priori error analysis of interior penaltydiscontinuous Galerkin finite element approximations to the biharmonic equation.
Journal ofScientific Computing , 30(3):465–491, 2007.[29] T. Yang, L. Dong, and S. N. Atluri. An Elementarily Simple Galerkin Meshless Method: theFragile Points Method (FPM) Using Point Stiffness Matrices, for 2D Elasticity Problems inComplex Domains. arXiv preprint arXiv:1909.04149 , 2019.[30] X. Wang, W. Pei, and S. N. Atluri. Bifurcation & chaos in nonlinear structural dynamics: Novel& highly efficient optimal-feedback accelerated Picard iteration algorithms.
Communicationsin Nonlinear Science and Numerical Simulation , 65:54–69, 2018.5031] G. D. Smith, G. D. Smith, and G. D. S. Smith.
Numerical solution of partial differentialequations: finite difference methods . Oxford university press, 1985.[32] E. Fehlberg. Low-order classical Runge-Kutta formulas with stepsize control and their appli-cation to some heat transfer problems. 1969.[33] N. M. Newmark. A method of computation for structural dynamics. American Society ofCivil Engineers, 1959.[34] H. M. Hilber, T. J. R. Hughes, and R. L. Taylor. Improved numerical dissipation for time inte-gration algorithms in structural dynamics.
Earthquake Engineering & Structural Dynamics ,5(3):283–292, 1977.[35] J. C. Houbolt. A recurrence matrix solution for the dynamic response of elastic aircraft.
Jour-nal of the Aeronautical Sciences , 17(9):540–550, 1950.[36] H. Vie and R. A. Miller. Estimation by limiting dilution analysis of human IL 2-secretingT cells: Detection of IL 2 produced by single lymphokine-secreting T cells.
Journal of Im-munology , 136(9):3292–3297, 1986.[37] J. R. Dormand and P. J. Prince. A reconsideration of some embedded Runge-Kutta formulae.
Journal of Computational and Applied Mathematics , 15(2):203–211, 1986.[38] X. Wang, Q. Xu, and S. N. Atluri. A Simple Local Variational Iteration Method and RelatedAlgorithm for Nonlinear Science and Engineering. arXiv preprint arXiv:1904.11021 , 2019.[39] J. P. Thomas, C. H. Custer, E. H. Dowell, K. C. Hall, and C. Corre. Compact implementationstrategy for a harmonic balance method within implicit flow solvers.
AIAA journal , 51(6):1374–1381, 2013.[40] T. A. Elgohary, L. Dong, J. L. Junkins, and S. N. Atluri. Time domain inverse problems innonlinear systems using collocation & radial basis functions.
CMES - Computer Modeling inEngineering and Sciences , 100(1):59–84, 2014.[41] J. H. He. Variational iteration method - A kind of non-linear analytical technique: Someexamples.
International Journal of Non-Linear Mechanics , 34(4):699–708, 1999.5142] G. Adomian. A review of the decomposition method in applied mathematics.
Journal ofMathematical Analysis and Applications , 135(2):501–544, 1988.[43] T. Fukushima. Picard iteration method, Chebyshev polynomial approximation, and global nu-merical integration of dynamical motions.
The Astronomical Journal , 113:1909–1914, 1997.[44] R. M. Woollands, A. Bani Younes, and J. L. Junkins. New solutions for the perturbed lam-bert problem using regularization and picard iteration.
Journal of Guidance, Control, andDynamics , 38(9):1548–1562, 2015.[45] X. Wang.
Optimized Picard Iteration Methods for Nonlinear Dynamical Systems with Non-Smooth Nonlinearities, and Orbital Mechanics . PhD thesis, Texas Tech University, 2019.[46] X. Wang, Q. Xu, and S. N. Atluri. Combination of the variational iteration method and nu-merical algorithms for nonlinear problems.
Applied Mathematical Modelling , 79:243–259,2020.[47] D. W. Mackowski. Conduction heat transfer: Notes for MECH 7210.
Mechanical EngineeringDepartment, Auburn University , 2011.[48] G. Voronoi. Nouvelles applications des param`etres continus `a la th´eorie des formes quadra-tiques. Deuxi`eme m´emoire. Recherches sur les parall´ello`edres primitifs.
Journal f¨ur die reineund angewandte Mathematik , 134:198–287, 1908.[49] T. Liszka and J. Orkisz. The finite difference method at arbitrary irregular grids and its appli-cation in applied mechanics.
Computers & Structures , 11(1-2):83–95, 1980.[50] S. Blanes, F. Casas, J. A. Oteo, and J. Ros. The Magnus expansion and some of its applications.
Physics Reports , 470(5-6):151–238, 2009.[51] X. Wang and S. N. Atluri. A novel class of highly efficient and accurate time-integrators innonlinear computational mechanics.
Computational Mechanics , 59(5):861–876, 2017.[52] J. C. Mason and D. C. Handscomb.
Chebyshev polynomials . Chapman and Hall/CRC, 2002.5253] S. Liu, G. Fang, B. Wang, M. Fu, and J. Liang. Study of Thermal Conduction Problem Us-ing Coupled Peridynamics and Finite Element Method.
Chinese Journal of Theoretical andApplied Mechanics , 50:339–348, 2018.[54] T. Blomberg. Heat conduction in two and three dimensions.