Research Article  Open Access
Jiao Wang, Lijun Su, Xinqiang Qin, "Hermite Radial Basis Collocation Method for Unsaturated Soil Water Movement Equation", Mathematical Problems in Engineering, vol. 2018, Article ID 8298915, 19 pages, 2018. https://doi.org/10.1155/2018/8298915
Hermite Radial Basis Collocation Method for Unsaturated Soil Water Movement Equation
Abstract
Due to the nonlinear diffusion term, it is hard to use the collocation method to solve the unsaturated soil water movement equation directly. In this paper, a nonmesh Hermite collocation method with radial basis functions was proposed to solve the nonlinear unsaturated soil water movement equation with the Neumann boundary condition. By preprocessing the nonlinear diffusion term and using the Hermite radial basis function to deal with the Neumann boundary, the phenomenon that the collocation method cannot be used directly is avoided. The numerical results of unsaturated soil moisture movement with Neumann boundary conditions on the regular and nonregular regions show that the new method improved the accuracy significantly, which can be used to solve the low precision problem for the traditional collocation method when simulating the Neumann boundary condition problem. Moreover, the effectiveness and reliability of the algorithm are proved by the onedimensional and twodimensional engineering problem of soil water infiltration in arid area. It can be applied to engineering problems.
1. Introduction
Unsaturated soil moisture movement [1] referring to the water in the soil is not filled with all the pores of the water movement, and it is an important fluid movement form of porous media. It has important practical significance in extensive engineering sciences such as runoff calculation, groundwater resource evaluation, soil science, hydrology, farmland irrigation and drainage, and soil salinization prevention and management [2]. At present, the soil water movement equation mostly refers to the Richards’ equation, which is established based on the BuckinghamDarcy law and continuous equation. The specific expression with water content as dependent variable is shown as follows:
In recent years, differential method and finite element method have been used to solve this problem; for example, Lei and Yang [3] used the Ritz finite element method to solve the onedimensional unsaturated soil water flow problem, used finite difference method to establish a numerical model for the onedimensional flow of unsaturated soil water, and applied it to practical problems [4] successfully; Ma et al. [5] used ADI alternating implicit difference dispersion and GaussSeidel iteration to solve the soil water movement equation; Li et al. [6] proposed a radial basis collocation method, combined with the radial basis function and collocation method, for the onedimensional soil water movement equation; Su et al. [7] solved the twodimensional unsaturated soil water movement equation by the radial basis collocation with difference method; Hu et al. [8] solved the boundary value problem by weighted radial basis point method. Many scholars [9–20] have solved a series of the numerical solution problems based on the radial basis point method.
Compared with the finite element method, the meshless method [21, 22] has the following main advantages: no grid, preprocessing, and simple adaptive analysis. It is easy to expand the lowdimensional problem to the highdimensional problem and solve the problem that is difficult to solve by the traditional numerical methods. Therefore, it was a hot research at home and abroad in last twenty years and considered as a very promising numerical analysis method. Radial basis function meshless method [23] has two main forms: collocation method (strong form), Galerkin method (weak form). The radial basis function collocation method (RBCM) is a meshless method for solving the numerical solution of partial differential equations. The shape function formed by the interpolation of the radial basis functions (RBFs) satisfies the Kronecker Delta function property, so that the essential boundary condition can be applied very well [9]. For the high order boundary value problem, the Neumann boundary condition affects the numerical calculation result directly. However, the research results in [10–20] showed that the accuracy obtained by using direct collocation scheme is a bit poor especially on boundary. Therefore, in order to satisfy the Neumann boundary condition accurately, the essential boundary condition contains the description of the field function and the description of the derivative of the field function.
Approximation of RBFs possesses optimal error estimation and convergence. However, the accuracy of the derivative of the interpolating function is usually very poor on the boundary when collocation method is used. In order to enhance the accuracy of the derivative of the approximation functions, Chu et al. [24] used Hermite radial basis point interpolation method for the nonuniform material gradient plane plate vibration; a meshless method was developed based on the collocation method and ridge basis function interpolation by Wang et al. [25]; Xing [9] constructed the governing equations of laminated plates with large deflection bending problem, approximated these field variables with RBF and HRBF, and discreted the control equations with least square collocation method; Wang et al. [26] proposed an improved meshless method which neither needs the computation of integrals nor requires a partition of the region and its boundary, and this method is applied to elliptic equations for examining its appropriateness; Parzlivand and Shahrezaee [27] mentioned a numerical technique which is a combination of collocation method, radial basis functions, the operational matrix of derivative for radial basis functions, and the new computational technique. Then they got the solution of a parabolic partial differential equation with a timedependent coefficient subject to an extra measurement; Krowiak [28] proposed a Hermite radial basis function differential quadrature method for higher order equations and so on [29–32].
The traditional radial basis point method has a large error at the Neumann boundary condition, so it cannot be used to deal with the Neumann boundary condition well. However, due to the existence of partial guidance in the Hermite radial basis function, it can reduce the error better. In this paper, the Hermite radial basis function combined with the collocation method was used to solve the unsaturated soil water movement equation with Neumann boundary condition. The Hermite radial basis function interpolation not only is satisfied with the Kronecker Delta function, but also can improve the precision of the interpolation function significantly. In addition, due to the nonlinear diffusion term of the secondorder partial derivative term in the soil water equation, the chain rule was used to preprocess the nonlinear term of (1) to avoid the phenomenon that the collocation method cannot be used directly. The remainder of this paper is organized as follows: In Section 2, some preliminaries were provided regarding the Hermite radial basis collocation and nonlinear term processing methods. In Section 3, we introduced the Hermite radial basis collocation scheme for the unsaturated soil water movement equation with Neumann boundary condition. In Section 4, the numerical examples were used to examine the appropriateness and efficiency of the method proposed in this paper by comparing it with some traditional methods and applying it to the real background engineering problem. Conclusions are given at the end of the paper.
2. Hermite Radial Basis Collocation Method (HRBCM)
2.1. Radial Basis Function
RBF [33–37] is also called a distance basis function, which is a function of the distance from the node to the node . The RBFs, which have the advantages of having a simple form, being independent of space dimension and isotropy, and being suitable for applying to the numerical approximation. In general, the RBFs can be divided into two categories. One is compacted supported radial basis function, and the other one is the global radial basis function. In contrast to the compacted supported function, the global radial basis function has better accuracy in the interpolation calculation, but the computational complexity is larger. It is a function that depends on the radial coordinate :where is the Euclidean distance from node to node . Typical RBFs include multiquadrics function: , inverse multiquadrics function: , Gaussians function: , TPS (thinplate splines): ,
where is a constant () and is an integer.
2.2. Hermite Radial Basis Function Approximation
Using the radial basis collocation method (RBCM), the derivative of the approximation function will produce a large error at the boundary of the region. Thus, RBCM cannot be applied to solve the differential equation with Neumann boundary condition, but HRBCM can be used to solve this problem well [22]. Hermite interpolation of the field quantity function in the local support region near the point is approximated by the RBFs and its derivative at the Neumann boundary which can be expressed as [22]wherewhere is the number of nodes in the local support domain and is the number of nodes on the Neumann boundary; let ; represents the total number of nodes. is the interior points, is the boundary point, is the RBFs, is the undetermined coefficient of the corresponding RBFs, and is the number of undetermined coefficients corresponding to the derivatives of the RBFs on the Neumann boundary. and denoted the direction cosine in and directions on the Neumann boundary, respectively.
We can obtain algebraic equations by interpolating at the kth point in the local support domain:where , .
We also can be obtain algebraic equation by interpolating at the mth point on the Neumann boundaries:
As shown the Figure 1, “0” is a collocation point, “1”, “2,” and “3” are the Neumann boundary points, “4” and “5” are the Dirichlet boundary points, and the “hollow” points in the circle are the other collocation points used to construct the approximating function at the “0” point.
2.3. HermiteType Collocation
Many practical problems can be solved by the solution of partial differential equations with the initial and boundary conditions—a fixed solution problem. Consider a boundary value problem as follows:where is the differential operator defining the governing differential equations to be satisfied in the open domain Ω, is the differential operator defining the boundary condition on boundary is the problem unknown, is the source term, and is the term associated with boundary conditions. For the unknown , the approximation represented by HRBFs can be expressed aswhere is the radical basis function centered at source point and is the expansion coefficient. is the number of source points in the domain, and is the number of source points on the boundary.
Based on the main principle of the collocation method, the expressions can be obtained by putting (8) into (7) as follows:when and (9) is the overdetermined system to be solved using the least squares method. When , (9) can be summarized as where is the coefficients matrix formed by the operators and under the action, is the vector formed by the unknown quantity, and is the source term and the term associated with boundary conditions formed by the vector.
The basic principle of the collocation method [9, 22, 25] is to make the nodes satisfied with the differential equation in the domain and the nodes on the boundary satisfied with the boundary condition. Then the nonlinear equations about the initialboundary value problem of differential equation are obtained.
3. Hermite Radial Basis Collocation Scheme
3.1. Discretization Scheme
Considering the twodimensional unsaturated vertical infiltration, the Richards’ equation and initialboundary conditions are shown as follows:where is the soil water content, is the initial water content, and are the unsaturated soil water diffusivity and conductivity, is the root water absorption term, is the time, is the horizontal distance, is the underground depth (positive upward), is the moisture content of the surface due to wet condition remain unchanged, is the Dirichlet boundary, and is the Neumann boundary.
Step 1 (preprocessing the nonlinear term). Because of the nonlinear term, the collocation method cannot be used directly. It needs to preprocess the nonlinear term, and the nonlinear diffusion term is considered as a composite function and resolved by the chain rule as follows:where and .
Step 2 (constructing the approximation function). By taking the Gaussian function as an example, the approximation function at the point can be described aswhere and represents the distance between and . Thus, the approximate solution of the Richards’ equation can be expressed by Hermite radial basis function interpolation as follows:where is the number of boundary points, is the number of inner points, and at the nodes numbered .
Step 3 (discrete the Richards’ equation). The forward difference and collocation methods are used for the time discretization and space discretization, respectively. The discretization scheme for the Richards’ equation can be obtained as follows:Collating (15), where and .
Step 4 (discrete the initial and boundary condition). Initial ConditionBoundary Condition
(1) Dirichlet Boundarywhere represents the distance between and .
(2) Neumann Boundarywhere represents the distance between and . and are the cosine of the Neumann boundary in the and directions.
Step 5 (constructing nonlinear equations). Combining (16), (17), (18), and (19), the matrix equation can be written as follows: wherewherewhere and .
3.2. Linearization Scheme
The iterative method is used to solve systems of the nonlinear equations. The common method of solving nonlinear equations are direct iteration method, NewtonRaphson method, modified NewtonRaphson method, incremental method, arc length method, and so on. A predictor correction method is adopted in this paper. The specific steps are as follows:
Assuming , are nonlinear term and moment, let , ; use , , as the forecast value of , , ; get K(), D(θ_{i}^{n}), K(); make it as the forecast value of K(), D(), K(); the solution is solved by the first iteration; then get K(), D(), K() separately, and make it as the correction value of K(), D(), K() and as the next forecast value of . Repeat the above steps until the difference between the two iterations before and after the end of the solution is less than the specified allowable error, which iswhere is the number of iterations and is the allowable relative error, which is determined by the exact requirement of the calculation; for example, . If it satisfied (23), we can replace with .
4. Results and Discussion
Example 1 (onedimensional unsaturated soil water movement equation). where , , exact solution , the root water uptake item , and the initial value and boundary value , are determined by exact solution. Gaussian function was selected as the radial function. In this paper, we discuss the error of norm: The comparative result of calculation accuracy about different algorithm is shown in Figures 2, 3 and Tables 1, 2.

 
FDM: finite difference method, RBCM: radial basis collocation method. 
(a) Number of nodes is 21; time step is 0.01
(b) Number of nodes is 51; time step is 0.01
(c) Number of nodes is 51; time step is 0.01
(d) Number of nodes is 51; time step is 0.01
(a) Number of nodes is 21; time step is 0.01
(b) Number of nodes is 51; time step is 0.01
(c) Number of nodes is 51; time step is 0.01
(d) Number of nodes is 51; time step is 0.01
It can be seen from Figures 2, 3 and Tables 1, 2 that the algorithm of HRBCM has the highest calculation accuracy among the other methods for the Neumann boundary problem. It is obvious that the parameter will increase with increasing the space nodes when we choose the different time steps, space distributions, and parameters . From Table 1, we can see that the HRBCM error is better. For the same time steps 0.01 and the number of spatial nodes 11, 21, and 51, the HRBCM error, RBCM error, and FDM error range from 2.5677 × 10^{−4} to 3.1904 × 10^{−4}, from 3.5145 × 10^{−3} to 2.1202 × 10^{−2}, and from 1.0018 × 10^{−2} to 1.1933 × 10^{−2} when the computation time is 0.5, respectively. The RBCM and the FDM have larger error at the Neumann boundary condition, and the HRBCM reduces the error at the Neumann boundary condition (see Figures 3(b), 3(c), and 3(d)). Therefore, this method is relatively superior, and the error accuracy is higher compared with other algorithms.
Example 2 (simulation of the soil water infiltration process).
Two most common infiltration models used to calculate the unsaturated conductivity and diffusivity in practice are BrooksCorey (BC) model and van Genuchten (VG) model [1, 38, 39]. According to the BC model and VG model, we can obtain the expressions of nonlinear terms (27) and (28).
BrooksCorey Model
Van Genuchten Modelwhere in BC model, in VG model, is an index parameter representing the soil’s pore size distribution, in BC model, in VG model, and is bubbling pressure (cm).
The HYDRUS1D 4.xx [38, 39] software was used to solve (26), and the relevant data pertaining to the infiltration of onedimensional soil under a fixed water head was simulated. Table 3 shows the values of the physical parameters for two soils. Figure 4 shows the simulated results of the onedimensional vertical (Figures 4(b) and 4(d)) and horizontal (Figures 4(a) and 4(c)) infiltration for the two soils at different times.

(a) The horizontal infiltration
(b) The vertical infiltration
(c) The horizontal infiltration
(d) The vertical infiltration
Figures 4(a), 4(b), 4(c), and 4(d) are the infiltration processes of the BC model and VG model in the horizontal and vertical directions, respectively. According to the norm error, we can get the norm error about different models in different wetting front. The norm error for the vertical infiltration was 3.8396% and 5.0855% when the wetting front depth is 23.5 cm and 48 cm, and it was 2.7664% and 1.1824% for the horizontal infiltration when wetting front depth is 20 cm and 40 cm for silt soil, respectively. Similarly, the norm error for the vertical infiltration was 4.5637% and 1.6325% when the wetting front depth is 20 cm and 40 cm, and it was 2.9743% and 2.1441% for the horizontal infiltration when the wetting front depth is 25 cm and 35 cm for sandy loam soil, respectively.
From Figures 2, 3, and 4, we can see that the new method proposed in this paper adapts to the numerical examples, and it is feasible and effective in simulating the actual background problem within the allowable range of error. So the new method proposed in this paper can be applied to solve practical onedimensional infiltration problems in engineering.
Example 3 (twodimensional unsaturated soil water movement equation in regular region). where , , exact solution , root water absorption item , initial value , and boundary value , are determined by exact solution. The radial basis function chooses the Gaussian function , and the error was estimated using the norm. Figure 4 gives the comparison result of numerical solution and exact solution in time about different algorithms. Table 3 gives the error accuracy of this paper in and time when the time step is 0.01.
The main advantage of the nonmesh method is that less points were used as the collocation points to calculate the numerical solutions at all nodes in the whole region. It can be seen from Figure 5 that the numerical solution and the exact solution fit are better and the error precision is well high. With the increase of the spatial nodes, the parameter increases.
(a) Number of nodes is 441 in this method
(b) Number of nodes is 2551 in this method
(c) Comparison of FDM and RBCM
Table 4 shows the optimal numerical error. It can be seen that the HRBCM error is 9.0137 × 10^{−6} when the collocation points number is 16 and the calculated time is 1. Thus, the numerical solutions at the more nodes were simulated by the 16th collocation points, and the numerical errors were shown in Table 5. For the same time steps 0.01 and the number of spatial nodes 16, 36, 121, and 441, the HRBCM error, RBCM error, and FDM error range from 8.1384 × 10^{−6} to 9.4099 × 10^{−6}, from 6.1375 × 10^{−5} to 7.3059 × 10^{−4}, and from 3.0860 × 10^{−3} to 3.6695 × 10^{−2} when the computation time is 0.5, respectively. Therefore, the method of this paper is relatively superior, and the error accuracy is higher compared with other algorithm. The optimal numerical error was observed by choosing the optimal collocation points, the appropriate time step, and parameter .

