A Hybrid Interpolating Meshless Method for 3D Advection–Diffusion Problems

: A hybrid interpolating meshless (HIM) method is established for dealing with three-dimensional (3D) advection–diffusion equations. To improve computational efﬁciency, a 3D equation is changed into correlative two-dimensional (2D) equations. The improved interpolating moving least-squares (IIMLS) method is applied in 2D subdomains to obtain the required approximation function with interpolation property. The ﬁnite difference method (FDM) is utilized in time domain and the splitting direction. Setting diagonal elements to one in the coefﬁcient matrix is chosen to directly impose Dirichlet boundary conditions. Using the HIM method, difﬁculties created by the singularity of the weight functions, such as truncation error and calculation inconvenience, are overcome. To prove the advantages of the new method, some advection–diffusion equations are selected and solved by HIM, dimension splitting element-free Galerkin (DSEFG), and improved element-free Galerkin (IEFG) methods. Comparing and analyzing the calculation results of the three methods, it can be shown that the HIM method effectively improves computation speed and precision. In addition, the effectiveness of the HIM method in the nonlinear problem is veriﬁed by solving a 3D Richards’ equation.


Introduction
In practical engineering problems, many partial differential equations are hard to solve using analytic methods, so numerical methods have become an important tool to obtain approximate solutions. The finite element method (FEM) [1][2][3] divides the continuum into finite elements, and the junction nodes of the elements are taken as the discrete points. The boundary element method (BEM) [4] reduces the dimension of the problem. Under the condition of the same discretization accuracy, the accuracy of the BEM is higher than that of the FEM. However, when dealing with elastic-plastic problems or large finite deformation problems, the advantage of BEM disappears. The meshless method [5][6][7] is a typical method. Since it is implemented on the basis of multiple groups of discrete nodes, the construction of a trial function is not constrained by elements or grids, and there is no need for mesh reconstruction. When simulating complex problems such as large deformation, the meshless method has great adaptivity [8,9].
Advection-diffusion equations are mainly used to describe the variation law of some physical quantities carried by fluid particles, such as the distribution of pollutants in nuclear waste pollution, heat conduction in fluid, and other physical phenomena [10,11]. The numerical method for advection-diffusion problems can be widely applied in the fields of hydromechanics, aerodynamics, environmental science and energy development. Although the FEM and BEM can be used in solving advection-diffusion problems, the numerical solutions obtained by these methods usually cause numerical oscillation for advectiondominated problems. The advantages of the meshless method include its great precision and its ability to overcome the mesh-dependence phenomenon. Therefore, utilizing the meshless method to calculate solutions of advection-diffusion equations has important research significance.
Moving least square (MLS) approximation [12][13][14] is a common way to fit discrete data using the meshless method. When exploring the meshless method, the element-free Galerkin (EFG) method [15], which is extensively researched and applied in various fields, provides a new direction for development. Adopting the improvement of MLS approximation [16], an improved EFG (IEFG) method that reduces the amount of computation is presented to solve diffusion equations, etc. [17][18][19][20]. For the convenience of applying boundary conditions, an interpolating EFG method [21,22] is proposed according to the interpolating MLS method [23]. By modifying the basic functions, an improved interpolating MLS (IIMLS) method [24] with a simpler format is proposed. From this, the improved interpolating EFG (IIEFG) method [25][26][27][28] is established for some classical problems.
To overcome the disadvantages of the IEFG method, such as its slow speed when calculating 3D equations, the dimension splitting method (DSM) [29,30] is introduced to further develop the meshless method. A dimension splitting EFG (DSEFG) method [31][32][33], which greatly improves the computation speed, is proposed and applied to calculate the numerical solutions of 3D equations. The hybrid complex variable meshless method [34][35][36] is presented according to the improved complex variable EFG method [37] and DSM. Moreover, a reducing dimension interpolating EFG method [38,39], which is highly accurate, is proposed. However, this method has difficulties, such as truncation error and calculation inconvenience, created by singular weight functions. To solve the above difficulties, an efficient interpolating EFG method [40] is further proposed and successfully used to solve 3D wave equations.
In this paper, a hybrid interpolating meshless (HIM) method is proposed for dealing with 3D advection-diffusion equations. To improve computational efficiency, a threedimensional equation is changed into correlative two-dimensional equations on the basis of DSM. The IIEFG method is applied on two-dimensional planes. Utilizing the IIMLS method, the approximation function with interpolation property is constructed. The FDM (finite difference method) is utilized in the time domain and axis after splitting. Dirichlet boundary conditions are introduced directly by setting diagonal elements in the coefficient matrix to one. Then, dealing with the solvable algebraic equations, the numerical solutions can be obtained. To prove the advantages of the HIM method, some examples are selected and solved by the IEFG, HIM and DSEFG methods. The influence of time step ∆t, node distributions, scaling parameters of influence domain d max , temporal interval t 0 and splitting directions on the numerical results is discussed. Comparing results obtained by the above three methods, it is demonstrated that the HIM method significantly increases calculation efficiency and accuracy. Moreover, by solving a 3D Richards' equation [41,42], the effectiveness of the HIM method in the nonlinear problem is verified.

Splitting Process of 3D Advection-Diffusion Problems
Consider a 3D advection-diffusion equation as follows.
When α ∈ Γ u , it has the essential boundary condition When α ∈ Γ q , it has the natural boundary condition The initial condition is where α = (x, y, z), u(α, t) is the field function used to describe physical quantities carried per unit volume of a fluid, t is time, f (α, t) is a source term representing the unsteady term that cannot be included in the governing equation, β x , β y , β z are diffusion coefficients which can be constants, γ x , γ y , γ z are advection coefficients in direction x, y, z, respectively. The 3D solution domain Ω gives the boundary Γ = Γ u ∪ Γ q with Γ u ∩ Γ q = Ø. u and q(α, t) are known functions on Γ u and Γ q , respectively. η x , η y , η z are the projections on the direction x, y, z of n = η x , η y , η z ) , which is a unit outer normal vector of Γ q . u 0 is a prescribed function at t = 0. The DSM was adopted to reduce the dimension of the 3D problem. Generally, a splitting direction that is easy to calculate and program is selected according to the form of control equation or boundary condition. For complex problems, we tend to choose the direction with a uniform cross section as the splitting direction. Assume Ω is split along direction x, and L − 1 2D planes are inserted in Ω with equal spacing, then a group of subdomains Ω (k) (k = 0, 1, · · · , L) can be obtained. The distance between every two planes is denoted by ∆x, and the 3D solution domain Ω can be given by If the values of x (k) and t are fixed, u, ∂u ∂x and ∂ 2 u ∂x 2 can be thought of as functions of y and z. Similarly, u = u(x, y, z) can be considered as a function in subdomains with the new boundary with the new form of boundary conditions To obtain its Galerkin weak form after splitting, applying the variational principle to Equation (8), we can obtain Using the integration aspect of the last integral term in Equation (11), we can obtain Form (9) and (10), on the essential boundary Γ (k) u , we have δu = 0; on the natural boundary Γ Therefore, Equation (12) becomes Equation (15) is the Galerkin weak form for the original 3D problem described by Equations (1)-(3) under the DSM.

The IIMLS Method on the Basis of Nonsingular Weight Functions
The IIMLS approximation [24] is utilized in subdomains Ω (k) to acquire the shape function, and the approximating function with interpolation property can be further formed.
To obtain the desired shape function, we choose the following basis function.
The following approximation of u(α (k) ) can be obtained according to the combination of Equations (21), (22) and (27).
The following formula can also express function u h (α (k) ) in Equation (34).
Compared with singular weight functions given by the fixed format in the IMLS method, a variety of weight functions can be tried for different problems in the IIMLS method. Moreover, difficulties created by singular weight functions, such as truncation error and calculation inconvenience, can be overcome.
The weight function in the calculation process of this paper was chosen as follows. where I and the point closest to it. d max is a custom scaling parameter that determines the size of the influence domain.

The HIM Method for 3D Advection-Diffusion Equations
According to Equation (36), ∇u, ∂u ∂t , ∂u ∂x and ∂ 2 u ∂x 2 are written as ∂u ∂t ∂u ∂x Substitute (36) and (40)- (43) into (15), then we obtain To obtain the linear discrete system, integral terms in (47) are discussed as follows. where Substitute (48)-(53) into (47), then we obtain Because δu T is arbitrary, we obtain where To obtain the discrete form of (59), the approximation of u on planes x = x (k) , (k = 0, 1, · · · , L) when the temporal interval takes t is denoted as Utilizing the FDM, the approximate values of . u, u and u on intermediate planes From (62), Equation (59) is expressed as . . .

Nu
then Equation (67) becomes The method of setting diagonal elements in the coefficient matrix to one was chosen to impose Dirichlet boundary conditions. For a discrete system of equations as if the unsolved displacement u r in row r lies on the boundary Γ u , it means u r can be known by the essential boundary condition u = u. Replace the principal diagonal element a rr of corresponding coefficient matrix with 1. Change the other elements and corresponding right-hand vector elements in row r to 0 and u r , respectively. Then, (71) becomes a 11 a 12 · · · a 1r · · · a 1n a 21 a 22 · · · a 2r · · · a 2n . . . . . . . . . . . .
Considering the row r of (72), we can obtain u r = u r . In particular, if the essential boundary condition is u = u = 0, the elements in column r other than the main diagonal element should also be changed to 0, i.e., Considering the row i, (i = 1, 2, · · · , n, i = r) of (73), we can obtain a i1 u 1 + a i2 u 2 + · · · + a i(r−1) u r−1 + a i(r+1) u r+1 + · · · + a in u n = v i .
The calculation of Equation (74) avoids the participation of u r .
Solving Equation (70), we can acquire solutions on planes x = x (k) , (k = 0, 1, · · · , L) at time ∆t. Using the linear interpolation method, we can obtain the approximate value of u(α, t) at any time t. Hence, u(x, y, z, t) in the 3D solving domain Ω is expressed as

Numerical Examples
To demonstrate the advantages of the HIM method, some specific examples were selected. Each problem was solved using three methods: the DSEFG, HIM, and IEFG methods. The calculated results were analyzed with different time steps, scaling parameters of influence domain, splitting directions, node distributions and time interval. Moreover, a 3D Richards' equation was chosen to verify the application of the HIM method in the nonlinear problem.
In this section, 4 × 4 points Gaussian integration is adopted to obtain numerical integration on the background grids. The calculation program is implemented using MATLAB. For comparison, the relative error δ is an important indicator to describe the accuracy. δ is defined by
The initial condition is The expression of the corresponding exact solution is given by Solving this problem using the HIM method, the splitting direction is z, d max = 1.35, t 0 = 0.1 s, ∆t = 0.001 s. The node distribution takes 11 × 11 × 11. Finally, the relative error δ = 1.9780 × 10 −4 . Computing time takes 2.90 s.
Comparing data given by these three methods, the HIM method is about three-times faster than the DSEFG method. Additionally, it is about 100-times faster than the IEFG method. Moreover, the precision of the HIM method is one order of magnitude higher than that of the DSEFG method. Additionally, it is two orders of magnitude higher than that of the IEFG method. Table 1 presents a comparison of relative error and computing time when the three methods take different scaling parameters of influence domain d max . All methods have the same ∆t = 0.001 s, t 0 = 0.1 s, and node distribution 11 × 11 × 11. We can see that it is necessary to find appropriate d max for different methods. To obtain the high precision results under current parameters, we let d max = 1.35 in the HIM method, while the IEFG method takes d max = 1.5 and the DSEFG method takes d max = 1.3.  Table 2 shows the comparison of relative error and computing time when all three methods choose different ∆t. The distribution of discrete nodes takes 11 × 11 × 11, t 0 = 0.1 s, d max = 1.3. Compared with the other two methods, the HIM method extremely improves computing speed and accuracy. We can also draw a conclusion that the HIM method is convergent with ∆t.  Figures 1-3 are the comparison of calculation results along axis x, y and z when the IEFG method, HIM method and DSEFG method take different t 0 . The node distribution takes 11 × 11 × 11, ∆t = 0.001 s. The figures indicate that the calculation results of the above methods agree with analytical solutions well. Figure 4 provides a 3D diagram of the calculation results given by the HIM method on layers z = k/10, (k = 0, 1, · · · , 10). ∆t = 0.001 s, the node distribution takes 11 × 11 × 11, t 0 = 0.1 s. the IEFG method, HIM method and DSEFG method take different 0 t . The node distribution takes 11 × 11 × 11, 0.001 s t Δ = . The figures indicate that the calculation results of the above methods agree with analytical solutions well. Figure 4 provides a 3D diagram of the calculation results given by the HIM method on layers / 10 z k = , ( 0,1, ,10) k =    the IEFG method, HIM method and DSEFG method take different 0 t . The node distribution takes 11 × 11 × 11, 0.001 s t Δ = . The figures indicate that the calculation results of the above methods agree with analytical solutions well. Figure 4 provides a 3D diagram of the calculation results given by the HIM method on layers / 10 z k = , ( 0,1, ,10) k = 

An Advection-Diffusion Equation with the Source Term
The governing equation is The boundary conditions are (0, , , ) y z t u y z t e + − =
The initial condition is u(x, y, z, 0) = e x +y +z .
The expression of the corresponding exact solution is given by u(α, t) = e x +y +z−t .
Compared to the results of all three methods, the computing speed and precision of the HIM method are obviously improved. Table 3 is the comparison of relative error and computing time when the above three methods take different t 0 with the distribution of nodes 11 × 11 × 11. Time step ∆t = 0.001 s. The data in Table 3 show that the HIM method takes less running time and has higher precision than the other two methods.  Table 4 shows the comparison of relative error and computing time when the above three methods take different node distributions. The HIM and DSEFG methods choose different splitting layers, and the IEFG method distributes different nodes along direction z. All methods have the same ∆t = 0.001 s, t 0 = 0.05 s. From the data, when there are a few nodes distributed in the discrete area, the HIM method still has high accuracy while maintaining fast calculation speed. With the increase in discrete nodes, the calculation accuracy of all three methods is improved gradually. Compared with the other two methods, the HIM method improves calculation precision and reduces running time. Moreover, the HIM method is convergent with the number of splitting layers.  Figures 5-7 show the calculation results along axis x, y, z when the node distribution of these three methods takes 11 × 11 × 11, ∆t = 0.001 s. It demonstrates that the calculation results given by all methods fit analytical solutions quite well according to the figures. Figure 8 provides a 3D diagram of the calculation results on planes z = k/10, (k = 0, 1, · · · , 10) given by the HIM method with d max = 1.4, 11 × 11 × 11 nodes, ∆t = 0.001 s, t 0 = 0.1 s.

An Advection-Diffusion Equation with Dirichlet Boundary Conditions
The equation is
The expression of the corresponding exact solution is given by where For the sake of simplicity, parameters in the above equations are β x = 1.4, β y = 1.7, Solving this problem by the IEFG method, d max = 1.8, ∆t = 0.01 s, t 0 = 0.1 s. Penalty factor α = 1.1 × 10 4 . The distribution of discrete nodes takes 11 × 11 × 11. Finally, the relative error δ = 0.0012. The computing time is 64.46 s.
Solving this problem using the DSEFG method or the HIM method, we discuss three cases under different splitting directions when the scaling parameter of influence domain d max = 1.8, the node distribution takes 11 × 11 × 11, ∆t = 0.01 s, t 0 = 0.1 s. Table 5 shows the comparison of relative error and computing time under splitting directions x, y and z given by the DSEFG method and HIM method. It demonstrates that these two methods are more accurate when the splitting direction is z. In addition, the HIM method has less computing time and higher precision than the DSEFG method according to the data.  Table 6 shows the comparison of relative error and computing time at time interval t 0 = 1.0 s when the above three numerical methods take different ∆t. With the decrease in time step ∆t, the calculation precision of the HIM method and DSEFG method is improved, which indicates that the two methods are convergent with the time step. Furthermore, although the HIM method has slight improvement in accuracy compared with the DSEFG method, it takes much less running time than the other two methods. Figures 9-11 are calculation results along axis x, y and z when the IEFG method, HIM method and DSEFG method take different t 0 with the node distribution 11 × 11 × 11, ∆t = 0.01 s. From the figures, the calculation results of all three methods agree with well with the exact ones. Figure 12 provides a 3D diagram of the calculation results at (x, y, k/10, 0.1 s), (k = 0, 1, · · · , 10) given by the HIM method. The distribution of discrete nodes is 11 × 11 × 11, ∆t = 0.01 s.

A 3D Richards' Equation for Unsaturated Flow in an Initially Dry Rectangular Block of Soil
Consider the 3D Richards' equation where the direction of gravity is (0,0,1) (the vertical coordinate is positive down), the pressure head < 0 h , the relative permeability λ = h r k e , the water content θ is given Figure 12. Calculation results on layers z = k/10, (k = 0, 1, · · · , 10) given by the HIM method.

A 3D Richards' Equation for Unsaturated Flow in an Initially Dry Rectangular Block of Soil
Consider the 3D Richards' equation where the direction of gravity is (0, 0, 1) (the vertical coordinate is positive down), the pressure head h < 0, the relative permeability k r = e λh , the water content θ is given by θ = θ r + (θ s − θ r )e λh with saturated water content parameters θ r , θ s and λ, K s is the saturated hydraulic conductivity related to the physical property of the medium. The boundary conditions are The initial condition is h(x, y, z, 0) = h r .
The expression of corresponding exact solution is given by where According to [43,44], (96) can be transformed into a linear equation by using the Kirchhoff integral transformation and the Gardner relation Substituting (104) and (105) Let a = b = L = 15.24 m, h r = −15.24 m, K s = 1 h/m, θ r = 0.15, θ s = 0.45, λ = 0.5. Solving this problem by the HIM method, the splitting direction is z, d max = 1.1, ∆t = 0.05 h, t 0 = 5 h. The distribution of discrete nodes takes 11 × 11 × 11. Finally, the relative error δ = 0.0239. Computing time takes 3.19 s. Table 7 is the comparison of relative error and computing time of the HIM method when ∆t is different. It shows that the method is convergent with time step ∆t when dealing with the 3D Richards' equation.

Conclusions
A hybrid interpolating meshless method is proposed for solving 3D advectiondiffusion equations. To improve computational efficiency, a three-dimensional equation is changed into correlative two-dimensional equations. After the dimension reduction by the DSM, the IIEFG method can be utilized in 2D subdomains. On the basis of the IIMLS method, the approximation function with interpolation property can be formed. For obtaining solvable algebraic equations, FDM is utilized in the time domain and axis after splitting. Setting diagonal elements in the coefficient matrix to one is chosen to impose Dirichlet boundary conditions directly. Then, the discrete system of equations can be solved, and the numerical solutions can be acquired.
To prove the advantages of the HIM method, some advection-diffusion equations are selected and solved by the IEFG method, HIM method, and DSEFG method. The influence of time steps t Δ , node distributions, scaling parameters of influence domain max d , temporal interval 0 t and splitting directions on the calculation results is discussed and analyzed. Collecting the data and comparing the numerical results obtained by the above three methods, it can be shown that the HIM method significantly increases the computing speed and precision. Moreover, by solving a 3D Richards' equation, the effectiveness of the HIM method in the nonlinear problem is verified.

Conclusions
A hybrid interpolating meshless method is proposed for solving 3D advection-diffusion equations. To improve computational efficiency, a three-dimensional equation is changed into correlative two-dimensional equations. After the dimension reduction by the DSM, the IIEFG method can be utilized in 2D subdomains. On the basis of the IIMLS method, the approximation function with interpolation property can be formed. For obtaining solvable algebraic equations, FDM is utilized in the time domain and axis after splitting. Setting diagonal elements in the coefficient matrix to one is chosen to impose Dirichlet boundary conditions directly. Then, the discrete system of equations can be solved, and the numerical solutions can be acquired.
To prove the advantages of the HIM method, some advection-diffusion equations are selected and solved by the IEFG method, HIM method, and DSEFG method. The influence of time steps ∆t, node distributions, scaling parameters of influence domain d max , temporal interval t 0 and splitting directions on the calculation results is discussed and analyzed. Collecting the data and comparing the numerical results obtained by the above three methods, it can be shown that the HIM method significantly increases the computing speed and precision. Moreover, by solving a 3D Richards' equation, the effectiveness of the HIM method in the nonlinear problem is verified.
In addition, it can be found that parameters need to be selected in an appropriate range for obtaining numerical solutions with high precision. In order to explain the variation in error and convergence speed when the HIM method takes different parameter settings, further theoretical research on error estimation and stability analysis will be carried out in the future.