Numerical Simulation of Convection-Diffusion Phenomena by Four Inverse-Quadratic-RBF Domain-Meshfree Schemes

The convection-diffusion type of PDEs is numerically solved by four numerical methods in this work. These four comparatively young numerical approaches are categorized as ‘domain-meshfree’ as they require no internal meshing but rely only on the collocation process amongst nodes via. the inverse-quadratic radial basis function (IQ-RBF). They are; the wellknown Kansa Collocation Method (KCM), the Hermite Collocation Method (HCM), the Radial Point Interpolation Method (RPIM), and the Dual Reciprocity Boundary Element Method (DRBEM). The work aims to demonstrate the use of IQ-RBF as well as to compare the practical use of the methods. Moreover, engineering senses of criteria judging the quality of the methods are considered. It is found in this work that while KCM is the simplest to construct and deploys, it is highly sensitive to the number of nodes and the IQ-RBF shape parameter. The asymmetric and populated matrix problem are alleviated when HCM or DRBEM are in use yet more CPU-time and storage seem to be the price to pay, particularly HCM. However, it is actually RPIM that has appeared to be an optimal choice under all the criteria imposed.


INTRODUCTION
Convection diffusion problems are known to be governed PDE mathematical models and they are found to appear in many branches of sciences and engineering such as biological, physical chemical, physical in fluid mechanics, astrophysics, meteorology, and multiphase flow in oil reservoirs, polymer flow and many other areas [1] .The unsteady-state form of the equation can be represented by; In the above definition, ̇ represents the transient term, is the gradient operator,  2 is the Laplacian operator, and  is the sum of additional source and/or sink terms.Very often the dimensionless parameters,  and  , that measure the relative strength of the diffusion to the convection is quite small leading to situations where thin boundary and interior layers are presented and singular perturbation problems arises [2].
In the study of this type of problem, the numerical solution in many cases has been tackled using traditional method such as the finite difference method (FDM), the finite element method (FEM), or the finite volume method (FVM), see the book by Zienkiewicz and Taylor [3] and references herein.In these conventional schemes, before the computing process can be performed, a mesh generation or meshing process over both the domain  and boundary  is required to take place.This inevitably makes the methods difficult and time-consuming particularly when solving complicated geometries.To improve the situation, affords may be put to automatically allocating the mesh-grid in hope to obtain optimal solutions and this is known as 'grid-adaptation method' [4][5].
Due to this undesirable aspect of mesh-generated numerical methods, together with a few more (please see [6]), some alternatives have been proposed over the past decade and one of which is those based on the use of a multivariate function called 'Radial Basis Function (RBF)'.The methods under this category is known as 'Meshfree/Meshless Method' [7].
A radial basis function (RBF) can be defined as a function of the distance of the point to the origin.That is,  is a RBF if () = (‖‖) , for  is a vector in ℝ  and ‖⋅‖ is a chosen norm.In this work, the main attention is paid to a specific type of inverse multiquadric and it is called 'Inverse-Quadratic (IQ)', defined as follows; With ‖⋅‖ 2 being the Euclidean norm.The appearance of the parameter  > 0 is known to have a great effect on the solution of the problem at hand.FIG. 1 demonstrates the IQ-RBF profiles influenced by different nonnegative shape  .Under the theory of interpolation, this RBF has been proved to be 'Globally Supported, Strictly Positive Definite Functions' and this function is  ∞ at the origin [6].While other kinds of RBF mentioned above have been receiving a huge interest from researchers [8][9][10], it is interesting to see this is not the case for IQ and, therefore, it deserves more investigation.
Numerical methods which require no mesh generation, at least in the domain, have recently been more popular.With the use of RBF, one of the pioneers of this idea is the work nicely done my Kansa in 1990 [11].In his work, a multiquadric type of RBF was introduced and applied to the interpolation problem and was extended to solving PDEs and the method is named "Kansa Method' ever since; some applications include groundwater contaminant transport [12], convection-diffusion problems [13], plate and shell analysis [14], microelectromechanical system analysis [15] and many more.Li [19] concluded that the accuracy of the RBF method was more superior than FEM.
The traditional version of Kansa Collocation Method (KCM), nevertheless, is known to severely suffer from having a fully asymmetric and populated collocation matrix, increasing the risk of being ill-conditioned.Attempts to avoid these drawbacks have been proposed and developed, see the book by Hua and Shantanu [17].
One of the studies aiming to alleviate the problems encountered in the use of KCM is that proposed by Fasshauer [18] in 1997.Fasshauer [18] described the scattered Hermite interpolation for the solution of elliptic partial differential equations by using collocation.An example of comparison studies is that done by Kazerm et al. [19].Another attempt is that proposed in 2002 by Wang and Liu [20] and it has been known as 'The Radial Point Interpolation Method (RPIM)'.One of the nice numerical experiments on applying this method is that carried out by Liu [21] in 2011, where it was concluded that the singularity problem can be improved (see also Bozkurt et.al. [22]).For some more recently and nicely documented, the interested reader is referred to [23] and [24].
Another comparatively new method that does not require any domain meshing is called 'Boundary Element Method (BEM)' [25].Nevertheless, an undesirable aspect is the need to discretize the domain into a series of internal cells to deal with the terms taken to the boundary by application of the fundamental solution, inevitably destroying some attraction of the method.To remedy this drawback, one of the attempts is the so-called 'Dual Reciprocity Boundary Element Method (DRBEM)', by Nardini and Brebbia [26].The idea has extensively been developed by many researchers such as the perturbation DRBEM [27], the separation of variables DRBEM [28], and the Laplace transform DRBEM [29].
Over the decade, even though some numerical comparison studies are available, see [30] for instance, the only criterion adopted appears to be the accuracy.This is rather impractical when there are other factors deserve to be taken into consideration as criteria as well.
In this work, four comparatively young domain-meshless numerical methods are studied and they are; the traditional or conventional Kansa Collocation Method (KCM), the Hermite Collocation Method (HCM), the Radial Point Interpolation Method (RPIM), and the Dual Reciprocity Boundary Element Method (DRBEM).The main objective is to shed more light into their overall practical aspects that should cover at least four criteria, listed below; 1.The overall accuracy and the error growth in time.2. The CUP-time requirement in computing process.3. The sensitivity to factors involved such as the density of nodes and shape values.4. The simplicity to setup and deploy.
The problem chosen for this task is one of the classical convection-diffusion PDE in two dimensions.

MATHEMATICAL BACKGROUND
Where   ,   are convection coefficients, and   ,   are diffusion coefficients.The last two terms  and the source term (, , ) are additional and needed only in specific cases.
For the test cases studied in this work, it is set that   =   =  the it leads to the following expression; Leading to; Subject to the initial condition (, , 0) =  1 (, ) and the boundary condition (, , ) =  2 (, , ) with  > 0 and  is a domain of the problem,  is its boundary,  1 ,  2 are known functions.

The Kansa Collocation Method (KCM)
The collocation scheme starts with considering the following elliptical partial differential equation defined on a bounded and connected domain ; for  ∈ Ω ⊂ ℝ  (6) Where ∂Ω is the domain boundary containing two non-overlap sections; Γ 1 and Γ 2 , with Γ 1 ∩ Γ 2 =  .These differential operators; Φ, and  1 ,  2 are applied on the domain, and the two boundary sections respectively.Three known functions (), (), ℎ() can be dependent The collocation scheme writes the approximate solution,  �(), as a linear combination of the basis function {(: −)}   , shown in the following form; Where   are coefficients and ‖. ‖ 2 being the Euclidean norm.The basis function used now is the inverse-quadratic radial type as defined previously.
Applying the operators Φ, and  1 ,  2 to on both domain and boundary sections, satisfying the governing system of equations, allows the system to arrive at; Where  = [  1  2 . . .  ] , the known  vector is as follows; ( ) ( ) And by setting  to be a matrix with entries   =  ��  −   � 2 �.Hence, for ,  = 1,2, . . .,  we have; Once   are obtained, the approximate solution are straightforward yielded.This method is known as 'Kansa', in honor of a great mathematician Prof. Edward Kansa who discovered the idea in 1990.The method has been applied to a wide range of problem ever since [31].It, nevertheless, is not of no shortcoming where it is known to suffer the problem of asymmetric interpolation matrix, , and the rigorous mathematical proof of its solvability is still not available [32].It is also very often produces low quality results particularly in boundaryadjacent region [33].
The Laplacian form, LHS of equation ( 5), can now be expressed as follows; Where the other terms in the equation can also be replaced by derivatives expressed above accordingly.

The Hermite Collocation Method (HCM)
In 1997, Fasshauer [21] proposed a way of interpolation by applying the self-adjoint operators  * , and  1 * ,  2 * to the governing system of equations and rewrite the approximate solution as; An application of this Hermite type of collocation method to elastostatic problem was done by Leitao [34].Some interesting implementations of the scheme to the transient and nonlinear plate problems can be found in [35,36] and [37].
In order to implement the Hermite concept, it is necessary that the Laplacian be applied twice resulting in the following so-called fourth-order biharmonic form; Hence, the inverse-quadratic type of RBF used in this work and its first four orders of derivatives can be respectively expressed as follows; Where   = � −   � 2 , for each pair of center nodes;   ,   .

The Radial Point Interpolation Method (RPIM)
The method was proposed by Wang and Liu [23], and it writes the approximate solution for a given PDE,  �() , as the linear combination of the basis function and monomials   (), shown in the following form; With  representing the number of polynomial basis (usually,  < ).The polynomial function can be chosen from Pascal's triangle which, for 2D problems, as; To ensure the unique solution of the system, additional  equations can be added as the constraint conditions, as follows; for  = 1,2,3, . . ., .This leads to the following form of matrix equations describing the interpolation process for all centers all over the domain.
Leading to Where Substituting these back into the collocation equation, yielding; By setting the shape functions matrix,  ̑(), as; Then the previous equation can be re-written as; Where the approximate solutions at each center can now be obtained.

The Dual Reciprocity Boundary Element Method (DRBEM)
The mathematical construction of the dual reciprocity boundary element method can start with the Poisson equation as follows; Which as its equivalent integral form, given by [29], as; Where  * is the fundamental solution and the term  �  is defined as , where  is the unit outward normal to Γ, and can be written as; With  and  being the number of boundary and internal nodes respectively,  can be now approximated by; Here, the function  is the radial basis function, then; For some particular solution  �  .By applying Green's theorem, the boundary element approximation to then it becomes, at a node i th ; Where the definition of the terms   and   can be found in [28].The index  is used for the boundary nodes which are the field points.Using a collocation technique, the above equation can be compactly expressed in matrix form as follows; By substituting  =  −1  into equation ( 35) making the right-hand side) a known vector.Therefore, it can be rewritten as; Where  = � � −  � � −1 .
Applying boundary condition(s) to equation (36) then it can be seen as the simple form as follows; Where  contains  unknown boundary values of ' and ' and the resulting linear system can be solved by Gaussian elimination scheme.Therefore, the solution is obtained via; When dealing with convection-diffusion forms of PDEs, it can be done by setting; Where ̇ is the transient term which will be discretized using RK4 as detailed in the following section.

COMPUTATION SETUPS
All numerical solutions obtained from the whole experiment are validated mainly by comparing to the exact or analytical solutions.For this, the following error norms are used; 2. Root-Mean-Square Error (  ); The transient term of the governing equation is tackled numerically by using the wellknown Runge-Kutta (RK4) method for all schemes under investigation.The process begins by setting; For (,   ) =   ,  ∈ ℝ  .Then let the following; Then, the time-increased  +1 is obtained by; (50) With a set of  computational nodes uniformly-distributed over the unit square (i.e. the summation of internal and boundary nodes).The interpolation function,  ̃(  ) at the i-th center node   , is defined as a linear combination of the inverse quadratic radial basis function as; With  ̃(  ) = (  ) for all  = 1,2, . . ., .By imposing this function on all nodes, it leads to the linear system expressed as; Functions' and the coefficient matrix;  = [  1  2 . . .  ]  , with the known vector function ...
. By the well-known Gauss-Seidel method, the coefficient matrix  is easily obtained and it will then be used to interpolate the value of (, ) at a new set of interpolation nodes via. the same linear summation as defined above.
With starting the interpolation process by utilizing  ctr = 10 × 10 centers, FIG. 2. shows the RMS error produced at a wide range of shape parameter when using   = 7 × 7 interpolation nodes.It is observed that the error is clearly varies with the shape and the best interpolation solution is seen to take place when  ∈ (2.0,5.0).At  ≈ 3.5 in particular, the results obtained are in good agreement with the exacts as illustrated in FIG. 3. The comparison of absolute error (..)produced when using two values of shape,  = 3.5 and  = 10.0, is provided in Tab. 1.It can be seen from this Table.that the shape has great effect on the final numerical solution, both locally and globally, of the domain.Nodes (0.333,0.000), (1.000,0.333),and (1.000,0.666)are some of the locations where the highest growths in ..are found; from ≈ 9.00 − 06 up to ≈ 1.00 − 01 .
This is defined on the domain with an elliptical boundary, expressed as; Where the boundary condition is taken directly from the exact solution which is expressed as follows; In this example, the shape values expected to lead to reasonably good results are adopted from the findings in example 1, i.e.  ∈ (2.0,5.0).The main focus is paid to the effect of the density of nodes contained in the domain and for this, at least 4 levels of nodes density are investigated.FIG. 4 displays two levels of density where the number of boundary nodes is kept constant at 40 for all cases.
Both RMS and  ∞ -error norms are carefully monitored, and the results are in Tab. 2.
As the number of nodes increases, it is found from both error norms that all the four numerical schemes under investigation produce results with more accuracy.Nevertheless, RPIM is seen to provide results with the comparatively lowest error at every level of nodes density.At the highest density of nodes, while KCM, HCM, and DRBEM have  ∞ > 1.0 − 03, it is RPIM that produces solutions with only  ∞ ≈ 3.25 − 04 .However, all errors revealed from all schemes in this example and at all levels of density of nodes, are found to remain under 5.00 − 03 which is reasonably well.An example of solutions obtained from DRBEM is plotted against the exact solution and is depicted in FIG. 5. Tab. 2. RMS-and  ∞ -error, produced at 4 levels of densities of internal nodes.Tab. 3 compares RMS errors measured at chosen 7 internal nodes obtained at two time levels,  = 0.5 and  = 1.0 , using  = 0.001 and 12 × 12 computational nodes.At the smaller  = 0.5, KCM and DRBEM are found to have approximately the same error magnitude where the other two, HCM and RPIM are close to each other with comparatively lower error range,≈ 5.00 − 04 to 9.00 − 04 .When the time level increases to  = 1.0, it is interesting to have found that while KCM, HCM, and DRBEM methods are seen to produce solutions with noticeably lower accuracy, it is RPIM that has been only slightly affected, with the RMS remaining in the range, i.e.< 8.00 − 04 .
When time increases even further, see FIG. 7, it can be clearly seen that all four domainmeshfree methods have growing their own  ∞ − error norm.The highest magnitude in error is found to have been produced by KCM with above  ∞ ≈ 1.00 − 01 at  = 2.00.DRBEM is also found to have approximately the same growth rate in error as KCM, but is only slightly lower in overall.On the other hand, HCM and RPIM are found to remain at lower level of this type of error norm with the lowest value of  ∞ − error is found to be obtained from RPIM, i.e. ∞ < 1.00 − 02 .The governing equation is set to contain = 0 and the source term is given by; The computational domain is (, ) ∈  = {(, ) ∈ ℝ 2 |0 < ,  < 2},  ∈ 0,  and all boundary conditions are taken from the analytical solution which is provided as follows; On this regular domain, it is set that the total number of centers or computational nodes is defined as the sum of the number of internal nodes and the boundary nodes or  =   +   .In order to cover as wide aspect; accuracy, computational time, node density and parameter effect, as possible, a large amount of numerical experiments were carried out and the main results are listed in Tab. 4. In this Table .,only comparatively best results for each case produced by a certain value of shape parameter,   , are presented for each node density level.
In terms of effect due to the shape value, it is clearly seen that it is very much varying for each numerical method and no noticeable correlation can be established.When compare between DRBEM and KCM method, for instance, while the optimal values of shape used in DRBEM remain under 1.00 − 01 , those found for KCM seem to increase beyond 1.00 + 01 as   = �  −   � 2 is reduced, for all   =   =  = 10, 10 2 and 10 3 .Moreover, the optimal shape does not always increase when the number of nodes creases and this aspect can be seen in the case of HCM and RPIM.When setting   =   = 5,000,   =   = 1 and with  = 19 × 19, FIG. 8 clearly shows that KCM is highly sensitive to the change of shape parameter value.In this FIGURE., the good quality of solutions produced by HCM, RPIM, and DRBEM can still be expected while utilizing the same shape values as listed in Tab. 4, yet this is not the case for KCM.Another evidence supporting this assumption is what is displayed in FIG. 9.It shows that the solution quality obtained from KCM is significantly improved when an optimal,   , is found and used for that particular case.
The exact solution is as follows; The initial and boundary conditions are then imposed using the above exact forms; The Dirichlet boundary conditions; At   =   = 1.0,   =   = 0.8 and time  = 0.05, Tab. 5 contains the numerical solutions at 8 points over the domain, obtained from all 4 schemes.
When time increases, nevertheless, from  = 0.05 to  = 0.5 , all four methods have been found to have lost their capability as clearly displayed in FIG. 10 .KCM is found to strongly be affected by the time-increment while HCM and RPIM are much less affected.When time has been even further increased, it is very interesting to have seen that it is actually DRBEM that can provide good results.FIG.11. confirms this argument where DRBEM leads to approximately the same quality of accuracy as those produced by RPIM.As the computation process continues from  = 1.0 to  = 2.0, however, the only scheme that can remain desirable solution accuracy while keeping the same   is actually RPIM and the evidence of this is shown in FIG.12.
So far, it has been clearly seen that RPIM is capable of producing the best quality of numerical results and is least sensitive to the change of the shape parameter.The last experiment is on the effect on nodes density in the computational domain.FIG. 13 depicts the RMS error generated by all four method when the nodes increase, for the case of   =   = 10,   =   = 0.8 at  = 1.5 when using  = 0.005.Once again, RPIM is seen to be slightly affected by the distance between nodes while strong fluctuations in error clearly appear in the case of DRBEM and KCM.When comparing RPIM with HCM, the results show that HCM is even less sensitive to nodes density particularly during 50 <  < 250 .Nevertheless, in terms of the accuracy, RPIM is noticeably more accurate with RMS as low as ≈ 1.00 − 05.

CONCLUSION
In this study, one of the classical PDEs namely convection-diffusion is numerical solved by four numerical domain-meshfree approaches.The common and crucial component of all the methods is radial basis function and the inverse-quadratic (IQ) type is chosen.To reach the main purposes; to shed more light into the practical use of IQ-RBF and to compare the effectiveness of the four methods, a large series of numerical experiments were carries out and 4 aspects of each method are measured.The computed results have revealed the following conclusions; 1.In terms of the overall accuracy and error growth in time; all four schemes are found to have approximately the same error growth rate but different in magnitude, i.e.KCM and DRBEM are roughly one-order of magnitude higher than the other two.2. In the time-consuming aspect, it appears that the Hermite type of collocation requires comparatively noticeably higher CPU-time.KCM is, on the other hand, found to be the fastest approach.3. The impact caused by the number of nodes, , and the change of shape parameter,  , are seen to be very high under the context of KCM.DRBEM and RPIM are found the be much less sensitive to those factors.4. In terms of the simplicity to construct and deploy, this can easily be judged by the way of approximating the solution  �(  ) .Based on this, it is concluded in this work that KCM is the simplest while HCM is the most complicated one followed by DRBEM and RPIM respectively.
All in all, if a good shape parameter can be achieved, this work has proved that KCM is the most suitable.for most practical use.Under the situation where internal nodes are less needed, DRBEM is an obvious choice where only boundary nodes play the major rules.For large and multidimensional domains, while KCM might encounter the problem of asymmetric and populated collocation matrix and while DRBEM can no longer handle properly, HCM can well be another alternative (provided that CPU-time has no limitation).In summary, it has been revealed that RPIM is the optimal approach in terms of all aspects and criteria mentioned above.

FIG. 4 .
FIG. 4. Two levels of nodes density with a fixed number of boundary nodes of 40 and with; (a) 51 internal nodes, and (b) 217 internal nodes.

Example 4 . 5 :
Transient with Dirichlet boundary conditionsThe special form of the problem is considered by setting;
2.1.The Governing Equation of Convection-Diffusion PDEThis work aims to numerically solve two-dimensional convection-diffusion problems that are modelled and governed by the following partial differential equation; Example 4.3: With Zero Source TermIn this example, the governing equation explained in Section 2.1 is solved with the main objective of investigating the error accumulated in time where the best shape parameter is still believed to follow the findings from example 4.1.In this case, we set   =   = 0.8 and   =   = 0.01 with zero sink and zero source terms, i.e.  = () = 0 .The governing equation is of the form as shown below; Tab. 3. RMS error comparison at two time levels with Δ = 0.001 and using 12 × 12 nodes.