Chihyu Liu, Chengyu Ku, , Jingen Xiao and Weichung Yeih,
Abstract: In this article, a meshless method using the spacetime collocation for solving the two-dimensional backward heat conduction problem (BHCP) is proposed. The spacetime collocation meshless method (SCMM) is to derive the general solutions as the basis functions for the two-dimensional transient heat equation using the separation of variables. Numerical solutions of the heat conduction problem are expressed as a series using the addition theorem. Because the basis functions are the general solutions of the governing equation, the boundary points may be collocated on the spacetime boundary of the domain. The proposed method is verified by conducting several heat conduction problems. We also carry out numerical applications to compare the SCMM with other meshless methods. The results show that the SCMM is accurate and efficient.Furthermore, it is found that the recovered boundary data on inaccessible boundary can be obtained with high accuracy even though the over specified data are provided only at a 1/6 portion of the spacetime boundary.
Keywords: Spacetime, collocation, meshless method, backward heat conduction problem,basis functions.
There are many backward heat conduction problems (BHCP) encountered in various engineering and industrial practices. Since the initial value of the temperature field of heat distribution is unknown in the BHCP, it is of importance that the BHCP is to recover the elapsed temperature as well as to obtain the initial temperature [Mera, Elliott, Ingham et al. (2002)]. Though boundary conditions for such a system are given, the absent initial temperature renders BHCP into the inverse problems [Movahedian, Boroomand and Soghrati (2013); Liu (2016); Liu, Qu and Zhang (2018)]. Since the unknown temperatures have to be determined from final time and boundary data which may be contaminated by measurement error, the BHCP is categorized into one of the ill-posed problems where small perturbation on the measured data may amplify error in the numerical solutions [Marin (2008)]. As a result, the approach for obtaining the numerical solutions is challenging because of the ill-conditioned system of the BHCP. For these reasons, the BHCPs remain difficult to solve.
Several mesh-based approaches have already been studied for solving the BHCPs, such as the boundary element method [Han, Ingham and Yuan (1995); Mera, Elliott, Ingham et al. (2001)], the iterative boundary element method [Lesnic, Elliott and Ingham (1998)],the finite element method [Wang and Mai (2005)], the finite difference method [Le Dret and Lucquin (2016)] , the local Petrov-Galerkin collocation method [Sladek, Sladek and Atluri (2004); Wu, Shen and Tao (2007)] or the lattice-free high-order finite difference method [Iijima (2004); Iijima and Onishi (2007)]. Besides, the group preserving scheme[Liu (2004); Liu, Chang and Chang (2006)], the Lie-group shooting method in the time direction [Chang, Liu and Chang (2007); Chang, Liu and Chang (2009); Liu (2010)], the fictitious time integration method [Chang and Liu (2010)], and the modified polynomial expansion algorithm [Chang, Liu and Wang (2018)] have also been used in solving the BHCP. In the past, the finite and boundary element methods play an important role as numerical tools for the numerical solutions of boundary value problems [Atalay, Aydin and Aydin (2004); Ferronato, Mazzia and Pini (2010)]. However, the mesh-based methods require complicated mesh generation in the domain and usually more computational time.
Recently published work demonstrates that the meshless methods such as the Trefftz method, the reproducing kernel particle method (RKPM) [Cheng and Liew (2009); Mu,He and Dong (2015)], the method of fundamental solutions (MFS) [Mera (2005); Hon and Wei (2005); Johansson and Lesnic (2008); Shi, Wang and Wei (2012)], the collocation method [Dai, Yue and Liu (2014)], the time domain collocation method [Dai,Schnoor and Atluri (2012); Dai, Yue, Yuan et al. (2014)], the smooth particle hydrodynamic method (SPH) [Jeong, Jhon, Halow et al. (2003); Shadloo, Oger and Le Touzé (2016); Alshaer, Rogers and Li (2017)], and the radial basis function (RBF)[Dehghan and Mohammadi (2014); Hanoglu and Sarler (2015)] have been developed and attracted much attention to solve heat conduction problems. Johansson and Lesnic introduced an application of the MFS defining the solution of the heat conduction problems as linear combinations of fundamental solutions to solve time-dependent heat conduction problems [Johansson and Lesnic (2008)]. Movahedian et al. [Movahedian,Boroomand and Soghrati (2013)] applied the Trefftz method based on adopting the exponential basis functions to solve the inverse heat conduction problem. Arghand et al.[Arghand and Amirfakhrian (2015)] proposed a meshless method based on the MFS and the RBF to solve the inverse heat conduction problem. Moreover, Liu solved the highdimensional inverse heat conduction problem using a multiple/scale/direction polynomial Trefftz method [Liu (2016)]. The indirect Trefftz method, also called the collocation Trefftz method (CTM) which can be traced back to 1926, is a boundary-type solution procedure. The CTM belongs to the meshless family in which numerical solutions for dealing with boundary value problems are approximated as basis functions satisfying exactly the governing equation [Kita and Kamiya (1995); Chen, Wu, Lee et al. (2007);Fan and Chan (2011)]. The CTM is originally proposed to solve boundary value problems in Euclidean space. Consequently, the adoption of the CTM for solving transient problems has never been reported [Karageorghis, Lesnic and Marin (2014); Ku,Xiao, Liu et al. (2016); Liu, Kuo and Jhao (2016)]. Recently, a novel advanced meshless method based on the spacetime coordinate system has been developed which can provide promising numerical solution for the transient flow problems [Ku, Liu, Xiao et al. (2017);Ku, Liu, Su et al. (2018)].
In this study, we present a novel spacetime collocation meshless method (SCMM) for solving the two-dimensional BHCP. The SCMM is to derive the general solutions as the basis functions for the two-dimensional transient heat equation using the separation of variables. Numerical solutions of the heat conduction problem are expressed as a series using the addition theorem. Since the basis functions are the general solutions of the governing equation, the boundary points may be placed on the spacetime boundary of the domain. The proposed method is verified by conducting several heat conduction problems. We also carry out numerical applications to compare the SCMM with other meshless methods. The remainder of this study is arranged as follows. Section 2 describes the formulation of the heat equation. In Section 3, we introduce the SCMM to explain the two-dimensional numerical solutions to the heat equation for heat conduction problems. In Section 4, we conducted several numerical examples of heat conduction problems to show the stability, convergence and high accuracy of the SCMM. Finally,findings and conclusions are summarized in Section 5.
The two-dimensional heat conduction problem in the spatial and time domain is described as follows.
where α2is the heat conductivity, u is the heat distribution, t is the time, L is the length of the space time domain, and T is the final elapsed time. For the two-dimensional direct heat conduction problem (DHCP) as Eq. (1), the initial and boundary conditions are required as follows.


where ?0(x,y,t=0)is the initial temperature and h(x,y,t)is the boundary data. Eq. (1)is the governing equation of the DHCP with given initial and boundary conditions for modeling the heat transfer process.
For the inverse heat conduction problem (IHCP), the initial or boundary conditions may not be provided. Typically, there are two types of IHCPs. The first type of IHCPs is the case where the initial value of the temperature field of heat distribution is unknown. The distribution of the temperature in the domain from the specified data on the final time boundary is determined by the BHCP [Movahedian, Boroomand and Soghrati (2013)].The elapsed temperature can be recovered by the BHCP. The governing equation of the BHCP is expressed as Eq. (1). The boundary condition is expressed as Eq. (3). The final time condition is required as follows.

where ?f(x ,y,t =T)is the temperature for the final time condition. Solving the governing equation as depicted in Eq. (1) with the given boundary and final time conditions is called the BHCP. The BHCP, which is an inverse problem, is to determine the value of h(x,y,t )for 0≤t≤Tfrom the data h (x,y,t)and ?f(x ,y,t =T). The final time heat distributions?f(x ,y,t =T)are unknowns based on the measurement.
The other type of the IHCP is to recover the temperature field only from the measured data on partial boundary. In addition, white noise for the measured data may also be considered [Movahedian, Boroomand and Soghrati (2013)]. A small random noise in the measured data?f(x ,y,t =T)may cause arbitrarily large error in the solution of h(x,y,t)because of the ill-posedness. In this study, the noised data on accessible boundary and final time data are represented as follows.

where h (x,t)is the exact boundary data in Eq. (5),?f(x,t =T)is the exact final time data in Eq. (6),is the noise data on accessible boundary,(x,t =T)is the noise data on final time data, rand is the random number generated by the uniform distribution in the range of[-1,1], and sis the noise level. The number of random number generated in Eqs. (5) and (6) are the same with the number of boundary collocation points.To deal with the BHCP, it is of importance to examine the accuracy and stability of the SCMM with the consideration of different levels of random noise in the measured data.
The spacetime coordinate system is adopted to deal with the transient modeling. Recently published work indicates that time can be defined as an independent variable based on the spacetime coordinate system [Ku, Liu, Xiao et al. (2017); Ku, Liu, Su et al. (2018)]. For example, we conduct a three-dimensional coordinate which includes one-dimensional in time and two-dimensional in space for solving a two-dimensional BHCP. Consequently,the original two-dimensional BHCP can be transformed into three-dimensional inverse boundary value problem (IBVP). It is significant that the temperature values at the final time and at space boundary may be given on the spacetime boundary of the domain.
Considering a two-dimensional heat conduction problem, there is two-dimensional simply-connected domain in space and one-dimensional in time. The two-dimensional heat conduction problem is expressed as

whereρdenotes the radius, andθdenotes the polar angle in the polar coordinate system.To formulate transient solutions of two-dimensional heat equation, the method of separation of variables is applied by assuming the following solution.

For simplicity, we consider the following equation.

Inserting the above equations into Eq. (7), we may obtain

Dividing by X(ρ ) Y (θ)T(t)on both sides in the above equation, Eq. (10) is rewritten as following equations.

whereλand?are unknown constant. To obtain the eigenvalue of the above equations,we adopt the constant v andkto ensure positive or negative constant. The detailed formulations including six cases are described as follows.
(1) The first case:?=0and λ=0
Considering the first case,?=0and λ=0, we obtain the solutions as follows.

where C1,C2,C3,C4and C5are constants. Using the boundary conditions of Y(ρ,θ=0,t)=Y(ρ,θ=2π,t), we may find that C3=0. Inserting Eq. (14) into Eq. (7),we may obtain

(2) The second case:?=0and λ=k2
Considering the second case,?=0and λ=k2, we may obtain the solutions as follows.

where I0is the modified Bessel function of the first kind,K0is the modified Bessel function of the second kind [Abramowitz and Stegun (1965); Watson (1995)],C6,C7,C8,C9and C10are constants. Substituting Eq. (16) into Eq. (7), we may have

(3) The third case:?=0and λ=-k2
Considering the third case,?=0and λ=-k2, we may obtain the solutions as follows.

where J0is the Bessel function of the first kind,Y0is the Bessel function of the second kind [Abramowitz and Stegun (1965); Watson (1995)],C11,C12,C13,C14and C15are constants. Substituting Eq. (18) into Eq. (7) may yield the following equation.

(4) The fourth case:?=-v2and λ=0
Considering the fourth case,?=-v2and λ=0, we may obtain the solutions as follows.

where C16,C17,C18,C19and C20are constants. Substituting Eq. (20) into Eq. (7) may obtain the following equation.

(5) The fifth case:?=-v2and λ=k2
Considering the fifth case,?=-v2and λ=k2, we may obtain the solutions as follows.

where C21,C22,C23,C24and C25are constants. Substituting Eq. (22) into Eq. (7) may have the following equation.

where Ivis the modified Bessel function of the first kind of the v order,Kvis the modified Bessel function of the second kind of the v order [Abramowitz and Stegun(1965); Watson (1995)],,,andare constants.
(6) The sixth case:?=-v2and λ=-k2
Finally, we consider the last case,?=-v2and λ=-k2. We may obtain the solutions as follows.

where Jvis the Bessel function of the first kind of the v order,Yvis the Bessel function of the second kind of the v order [Abramowitz and Stegun (1965); Watson (1995)],C26,C27,C28,C29and C30are constants. Substituting Eq. (24) into Eq. (7), we have


where wis the order of the basis functions for approximating the solution. For a simply connected domain, only positive basis functions need to be considered. Therefore, the above equation can be simplified as following equation.

Using the final time and boundary conditions given by Eqs. (4) and (3), Eq. (27) may be discretized at numerous boundary collocation points on the spacetime boundary.Consequently, a system of linear equations may be obtained as following equation.

C=,A is a matrix (size of p×q) of the basis functions,Cis a vector (size of q×1) of unknown coefficients to be determined,Bis a vector (size of p×1) of given Dirichlet boundary data at boundary collocation points,pis the number of boundary collocation points,qis the number of the terms related to the order of the basis function as depicted in Eq. (28), which can be defined asq =4w +4w2+1,ρ1, ρ2,…,ρpare radiuses,θ1, θ2,…,θpare polar angles,,,, …,are unknown coefficients to be evaluated, and bb1, b b2, bb3, …,bbpare Dirichlet boundary data. The unknown coefficients for the three-dimensional spacetime domain may be determined for solving Eq. (28).To obtain the temperature field of the heat distribution for the spacetime domain, the inner collocation points in threedimensional spacetime domain have to be placed. The heat distribution at inner collocation points can be computed by Eq. (28) for the three-dimensional spacetime domain.
The scenario investigated is the modeling of two-dimensional DHCP. The timedependent heat transfer in an isotropic and homogeneous space is expressed by the heat equation given by Eq. (1). The initial condition is given as follows.

The Dirichlet boundary conditions are given using the following exact solution.

In this example, the length and the width of the space domain is πm, the heat conductivity is 1 W m-1K-1and the final elapsed time is 0.25 h. The initial condition is applied on the bottom side of the spacetime domain and the boundary conditions are applied on both left and right sides of the three-dimensional spacetime domain, as shown in Fig. 1(a).

Figure 1: Illustration of the two-dimensional DHCP and BHCP
For a simply connected domain, it is necessary to locate the source point inside the domain and the number of source point is only one for in the SCMM. The source point is the reference point for computing the radius and polar angle. There are 1280 boundary collocation points. The numbers of boundary collocation points in space domain and in time domain are considered to be 256 and 1024, respectively. The Dirichlet boundary values obtained from the exact solution are specified on boundary collocation points which placed on four vertical sides of the spacetime domain. The order of the basis function is 20.
To show the accuracy of the computed result from the proposed method and that from other meshless methods, we conduct a comparison example with the consideration of DHCP. We place 4056 inner collocation points which placed uniformly inside the three-dimensional spacetime domain to obtain the computed results of the heat distribution. Then, we select the profiles of the computed results at different time to compare with the exact solution.
Fig. 2 demonstrates the computed results from the proposed method and the exact solution.The computed heat distribution from the proposed method are compared with those of the RKPM and Trefftz method based on EBFs, as depicted in Tab. 1. For the RKPM, 121 scattered nodes inside the domain were considered [Cheng and Liew (2009)]. As for the Trefftz method based on EBFs, 80 points on the boundary and 100 points as the initial nodes inside the domain were considered [Movahedian, Boroomand and Soghrati (2013)]. It is found that the computed results using the SCMM agree very well with the exact solution.We then compare the maximum absolute error of the proposed method with those of the RKPM and Trefftz method based on EBFs. The maximum absolute error of the RKPM and Trefftz method based on EBFs reach only to the order of 10-3and 10-5[Movahedian,Boroomand and Soghrati (2013)]. On the other hand, the maximum absolute error of the SCMM can reach to the order of 10-10, as depicted in Fig. 3. It is significant that the SCMM is able to yield highly accurate results.

Table 1: The computed heat distribution from the proposed method

Figure 2: Comparison of results for solving two-dimensional DHCP

Figure 3: Comparison of the absolute error of the computed results and the Trefftz method based on EBFs at final time
The second example is the modeling of two-dimensional ill-posed BHCP. Because the BHCP is ill-posed in nature, it is quite challenging for solving the BHCP with the absent initial conditions. The governing equation is given as Eq. (1). We assume the boundary conditions to be the Dirichlet boundary conditions. The Dirichlet boundary data are given as follows.

The final time condition is given by

In this example, the length and the width of the space domain is 1 m, the heat conductivity is 1 W m-1K-1and the final elapsed time is assumed to be 0.25 h. The numerical procedure for the two-dimensional BHCP is similar with the previous two-dimensional DHCP. The final time and boundary conditions are specified on the top and four vertical sides of the three-dimensional spacetime domain, as depicted in Fig. 1(b). We collocate 1680 boundary collocation points and a source point. The numbers of boundary collocation points in space domain and in time domain are considered to be 336 and 1344, respectively. The Dirichlet boundary data using the exact solution are applied on boundary collocation points and are placed on four vertical sides of the three-dimensional spacetime domain. The order of the basis function for the numerical analysis is 20.
We select the profiles of the numerical solution on different time to view the results clearly.To obtain the temperature values of the heat distribution, 4056 inner points are collocated uniformly inside the three-dimensional spacetime domain. The absolute error of the recovered numerical solution is demonstrated in Fig. 4(a). The maximum absolute error is found in the order of 10-6. The global space-time multiquadric (MQ) method as well as the Trefftz method based on EBFs have also been applied to solve this BHCP. For the MQ method, 121 regular grid points in space domain were considered [Li and Mao (2011)]. As for the Trefftz method based on EBFs, 16 uniform discrete data at the final time were considered [Movahedian, Boroomand and Soghrati (2013)]. We compare the absolute error of the proposed method with those of the MQ method and Trefftz method based on EBFs.Based on their results, the root mean square error (RMSE) using the MQ method and Trefftz method based on EBFs can only reach to the order of 10-3and 10-6. However, the RMSE using the SCMM can reach to the order of 10-7. It is significant that the initial heat distribution can be recovered with high accuracy using the proposed method.
Since the proposed method is advantageous in solving problems of irregular geometry,the following example is the two-dimensional BHCP enclosed by irregular boundary.With a two-dimensional simply connected domain, the governing equation can be expressed as Eq. (1). The two-dimensional boundary, as shown in Fig. 5(a), is defined as

where ρ(θ)=sec(3θ)sin6θ, 0≤ θ≤2π.
We assume that the boundary conditions are from the analytical solutions. The boundary data are given by following exact solution.

The final time condition is given by

In this example, the heat conductivity is 1 W m-1K-1, and the final elapsed time is 0.5 h.To examine the effectiveness and stability of the SCMM, we consider that the boundary and final time data are contaminated by random noise. The noisy data as depicted in Eqs.(5) and (6) are used in the calculation of two-dimensional BHCP. The two-dimensional BHCP is one-dimensional in time and two-dimensional in space, as depicted in Fig. 5(a).We therefore transform the domain of the two-dimensional BHCP into a threedimensional spacetime domain, as depicted in Fig. 5(b). The final time conditions are given on the top side of the spacetime domain. The boundary conditions are given on the circumferential boundary of the spacetime domain. The two-dimensional BHCP is transformed into the three-dimensional IBVP in which the bottom side boundary values,also known as the initial conditions, are not assigned.
We adopt 486 boundary collocation points and one source point. The numbers of boundary collocation points in space domain and in time domain are considered to be 236 and 250, respectively. The order of the basis function for the analysis is 8.
To obtain the computed results at different times, 1090 inner points are placed. The computed results indicate that the numerical solutions agree well with the exact solution.Fig. 6(a) shows the absolute error of the recovered solutions. The accuracy of the error is found in the order of 10-11and highly accurate recovered solutions can be obtained using the SCMM for this two-dimensional BHCP enclosed by irregular boundary. In addition,three different noise level with s =0.01,s =0.05, and s=0.1are considered to investigate the application of the SCMM for dealing with the two-dimensional BHCP.The results illustrate that accurate recovered solutions with the accuracy in the order of 10-7, 10-6and 10-5can be obtained for the noise level s =0.01,s =0.05, and s=0.1,respectively, as depicted in Fig. 6(b), Fig. 6(c) and Fig. 6(d).

Figure 4: Comparison of the absolute error of the computed results and the Trefftz method based on EBFs at initial time

Figure 5: Illustration of the collocation scheme for the two-dimensional BHCP

Figure 6: Absolute error of the computed results with exact solutions with the consideration of noise level
Previous examples, including two-dimensional DHCP and BHCP, illustrate that the SCMM is appropriately verified using the harmonic boundary conditions. In this particular example, we conduct an example with non-harmonic initial and boundary conditions. The governing equation can be expressed as Eq. (1). The two-dimensional boundary is defined as following equation.

The non-harmonic boundary and initial data are given using the following equation.

We consider the final elapsed time to be 0.5 h. It is significant that the above equation is not the exact solution of the two-dimensional BHCP. All the input data are the same with those for the previous example of two-dimensional BHCP except the boundary shape,non-harmonic initial and boundary conditions. Since this example does not have an exact solution to verify the correctness, we first carry out a forward analysis of DHCP to obtain the heat distribution at the final time. After that, a BHCP is conducted where the computed heat distribution at the final time are assumed to be provided. Through the BHCP, the initial heat distribution at bottom side of the three-dimensional spacetime boundary can then be obtained. The algorithm for the numerical implementation is represented as a flow chart depicted in Fig. 7.
In this example, we adopt 3905 boundary collocation points and a source point. The order of the basis function is 10. To validate the example, we compare the given non-harmonic initial condition with the computed initial temperature from the BHCP. The computed results with the given non-harmonic initial condition are illustrated in Fig. 8. Results obtained show that the computed initial temperature from the BHCP agree well with the given non-harmonic initial condition, as depicted in Fig. 8. The difference between given non-harmonic initial condition and computed initial temperature from the BHCP can reach to the order of 10-5, as depicted in Fig. 9.

Figure 7: The flow chart of the validation procedure for the two-dimensional BHCP

Figure 8: Comparison of the results computed with the given non-harmonic initial data

Figure 9: The difference between given non-harmonic initial condition and computed initial temperature
The last example is to investigate the effectiveness and stability of the SCMM to solve the two-dimensional inverse BHCP in a rectangular domain. In this example, the boundary and final time data contaminated by random noise are adopted. The timedependent heat transfer in an isotropic and homogeneous space is expressed by the heat equation given by Eq. (1). The spacetime boundary data using the Dirichlet boundary condition are applied from the following exact solution.

In this example, we consider the heat conductivity to be 1 m-1K-1, length and width of the space domain to be 1 m, and final elapsed time to be 0.25 h. The two-dimensional inverse BHCP in this study is defined that not only the initial conditions are absent, but the partial data of the final time and the spacetime boundary conditions are also not specified. In this example, four cases with the consideration of the different combinations of missing initial and boundary conditions are considered, as depicted in Fig. 10. In case A, we consider absent initial condition which is the traditional BHCP, as depicted in Fig. 10(a). In case B, we consider partial data of the boundary conditions are absent in which the accessible boundary values are provided at 1/2 portion of the overall spacetime boundary, as shown in Fig. 10(b).In case C, we consider that the accessible boundary values are specified at 1/3 portion of the overall spacetime boundary, as shown in Fig. 10(c). In case D, we further consider that the accessible boundary values are provided at only 1/6 portion of the overall spacetime boundary,as shown in Fig. 10(d). It is significant that case D is the most challenging and difficult example. We collocate 800 boundary collocation points and a source point. The order of the basis function is 10. The maximum absolute error for case A to case D are illustrated in Fig.11. Results obtained show that the maximum absolute error for case A to case D are 10-10,10-7, 10-6and 10-5, respectively. It is found that the absent boundary data can be recovered with very high accuracy using the proposed method, as depicted in Tab. 2. To show the accuracy of the SCMM, we consider that the measurements of the specified data are polluted by random noise where the noise level is s =10-3. The maximum absolute error for case A to case D with the consideration of noise level s =10-3is illustrated in Fig. 12. It is found that recovered boundary data with the accuracy at the order of 10-2on inaccessible boundary can be obtained even 5/6 portion of the overall spacetime boundary are absent.

Figure 10: Four cases with the consideration of missing initial and boundary conditions

Table 2: The maximum absolute error for two-dimensional inverse BHCP

Figure 11: Absolute error of the computed results for case A to D

Figure 12: Absolute error of the computed results for case A to D with noise level s=10-3
A novel spacetime collocation meshless method for solving the two-dimensional backward heat conduction problem is proposed. We may conclude the following findings.For the modeling of the two-dimensional BHCP, we propose an innovative concept that one may collocate boundary points in the spacetime coordinate system. Both the boundary and initial conditions can then be regarded as boundary conditions on the spacetime domain boundary. We therefore transform the domain of the two-dimensional BHCP into a three-dimensional spacetime domain. As a result, we solve a threedimensional IBVP instead of the original two-dimensional BHCP by using the SCMM.Moreover, we demonstrate that recovered boundary data with the accuracy in the order of 10-2on inaccessible boundaries can be obtained even the accessible data are only specified at 1/6 portion of the whole spacetime boundary.
Acknowledgement:This study was partially supported by the National Science Council under project grant MOST 106-2621-M-019-004-MY2. The authors thank the National Science Council for the generous financial support. The authors declare that there is no conflict of interest regarding the publication of this paper.
Computer Modeling In Engineering&Sciences2019年1期