Abstract
The purpose of the present work is to develop a quasi-3D numerical method that can be used to study the diffusion mechanism of grout injection in a rock fracture based on the collocated structured grid of the finite volume method (FVM). Considering the characteristics of fracture in geometry that the aperture is much less than its length and width, the Hele-Shaw model is introduced to deduce the z-derivatives of velocities u and v at walls, which is a function of the relevant average velocity and the fracture aperture. The traditional difference scheme for the diffusive term is partly substituted with the derived analytical expressions; hence a three-dimensional problem of grout flow in the parallel fracture can be transformed into a two-dimensional one that concerns fracture aperture. The new model is validated by the analytical solution and experimental data on three cases of grouting in the parallel-plate fracture. Compared with the results from ANSYS-Fluent software, the present model shows better agreement with the analytical solution for the distribution of pressure and velocity. Furthermore, the new model needs less grid unit, spends less time, but achieves greater accuracy. The complexity of the grout flow field in the rock fracture is reduced; thus the computational efficiency can be improved significantly.
1. Introduction
Grouting is widely used as the main means of antiseepage reinforcement of fractured rock mass in the fields of civil engineering, water conservancy and mining. The knowledge of grout diffusion mechanism in the rock fracture can significantly guide the design and construction of the grouting project. The rock mass is usually considered as a discontinuity separated by fractures, in which the grout flows, spreads, and finally seals the channel. For decades, many researchers suggested that the grout propagation in a fracture may have a close relationship with the fracture aperture, the grout properties, and the grouting procedure such as the grouting pressure and time [1]. The theoretical analysis, experimental investigation, and numerical simulation are all effective ways to carry on the correlational study.
In terms of the analytical method, the first discussion of Bingham fluid flow in the fracture channel emerged with the work of Wallner [2]. The grout transport was described as linear flow in a channel, and the penetration length was calculated by the one-dimensional channel flow model. Dai and Bird [3] proceeded to study the radial flow rate of a Bingham fluid between the parallel plates. Taking the rheological properties of the grouts into account, Hässler [4] used similar flow rate equations to estimate the flow of cementitious grout in fractures. Subjected to pressure gradients, the transient flow of Bingham viscoplastic grouts in parallel fractures was solved analytically by Amadei and Savage [5]. The results show that the extent of a rigid central layer in fractures depends on the pressure gradient, the grout properties, the fracture aperture, and surface roughness. After that, the analytic relations of grout spread to time period were presented systematically by El Tani [6], including channel flow between parallel plates and radial flow between parallel discs, nonlinear Newtonian fluids, and Bingham material. Considering the influence of inertia forces in the governing equation, a simplified model for predicting grout flow was proposed by Xiao et al. [7]. Zhang et al. [8] derived the analytical expressions of evaluating the propagation of cement and sodium silicate grout in a planar fracture, revealing that the rheology of the grout has a most significant influence on the penetration length. All the above-mentioned models share a similar basic assumption that the grout propagation in the fracture is one-dimensional or two-dimensional flow.
Meanwhile, some field tests are conducted to verify the analytical model [9–11]. However, the concealment of rock fractures causes great difficulty for site test, while lab experiment is more feasible to observe the grout flow in the fracture. In this regard, the parallel-plates fracture, formed by two stiff rectangular plates or discs, is the most widely used model to investigate the grout flow and propagation [12–15]. On this basis, the effects of the grout properties, the grouting procedure, the geometry of the fracture, and the underground environmental conditions such as underground water on the diffusion behavior of grout in the fracture were studied [16–18].
Compared to the analytical or experimental method, the numerical simulation has great advantages in predicting the grout propagation with wide applicability, more efficiency, and low cost. It is flexible to simulate the grouting under different conditions and analyze the influences of various factors on the grouting effect. Based on the work of Wallner [2], Hässler et al. [19, 20] made the first effort to simulate the grout flow. The fracture plane in a rock mass was described as an orthogonal network of one-dimensional channels. The study was improved by Eriksson et al. [21]. In that work, a fracture network with varying aperture, constructed from the one-dimensional rectangular channel elements, was modeled to calculate the grout spread and predict the grouting result. Based on Monte-Carlo method, the stochastic fracture networks were generated to predict the grout penetration [22, 23]. Similarly, a stochastic discrete feature network model incorporating boreholes and the fractures intersecting them was presented to evaluate the foundation grouting performance in fractured rock [24]. However, this approach is limited by the requirement of site-specific fracture information. Mohajerani et al. [25, 26] developed an efficient algorithm of explicit grout forehead pressure to simulate the grout propagation within the fractures. In brief, whether applying a structured or stochastic network, the aforementioned works mostly take one-dimensional channel model as the basic unit. This base model always assumes that the grout flows in a straight line and occupies the entire cross-section along the channel. As a matter of fact, at the initial stage of grouting into a large-size fracture, the grout diffuses in two even three dimensions before the flow front of the grout reaches the fracture boundary. Only when the grout is constrained by the boundary and develops fully, it can be approximated as one-dimensional flow.
In recent years, the quick-setting grout, such as cement and sodium silicate grout, has been widely applied to prevent the seepage in fractured rock mass [27–29]. The curing time of such grouts is short, generally within a few minutes even in tens of seconds in the water disaster prevention engineering. Therefore, it is very difficult for the grout to reach a one-dimensional flow state during a short period of time before curing in the large-size fracture. In this case, it may be more appropriate to conduct two-dimensional or three-dimensional numerical analysis for the grout propagation process based on the computational fluid dynamics theory.
The aim of this paper is to provide an available numerical method for the study of such grouting engineering. Based on the conventional FVM, a quasi-3D numerical model of grout injection in the rectangular parallel fracture is presented considering the fracture aperture. Because the fracture is length and width far greater than the aperture, the Hele-Shaw approach is employed to deduce the partial derivative of velocity with respect to z, which is related to the average velocity and the fracture aperture. In this way, the three-dimensional flow of grout in the fracture is converted to a two-dimensional one. The proposed model is solved by the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm and the moving interface of the grout tracked by the volume of fluid (VOF) method. The new method solves faster with higher precision than ANSYS-Fluent. The numerical results, including the diffusion radius, the pressure and velocity distribution, are compared with the analytical solution and experimental data to verify the validity of the new method.
2. Quasi-3D Numerical Model
2.1. Governing Equations
For analyzing the grout flow in a parallel fracture illustrated in Figure 1, the following assumptions are made:(i)The grout is an incompressible Newtonian fluid.(ii)Flow is laminar.(iii)Fracture walls are stiff.(iv)No slip occurs at the solid interface.
It is also assumed that the velocity in the z-direction is zero; accordingly, there is no pressure gradient. It can be expressed as
The geometric feature of the rock fracture is that the aperture value is much smaller than the length and width; therefore the Hele-Shaw model [30] can be used for analyzing the grout flow in the fracture, with the following flow equations:
Based on Stokes’ hypothesis, the Newtonian flow of the grout in the fracture can be expressed by the mass and momentum conservation equations. The mass conservation equation is given bywhere is the density and u is the velocity vector.
The momentum equations can be written as follows:where p is the pressure and is the viscosity. u, v, and are the flow velocities in the x-, y-, and z- directions, respectively. Similarly, , , and are three components of the constant source term.
If introducing a general variable ϕ, (4) and (5) can be expressed as the following form:where is the diffusion coefficient and ϕ stands for velocity u, v, or in momentum equations. Equation (6) is called the transport equation for variable ϕ. It clearly highlights the various transport processes: the rate of change term and the convective term on the left-hand side and the diffusive term and the source term, respectively, on the right-hand side.
2.2. Three-Dimensional Control Volume Element
The simulation domain is divided into a number of discrete control volumes (CVs) based on a Cartesian grid of equally spaced horizontal and vertical grid lines. For each control volume, the continuity and momentum equations are approximated using algebraic expressions involving the values of the unknown u, v, p at the center of that CV and at the centers of neighbouring CVs [31]. A typical three-dimensional control volume element is shown in Figure 2. Positioned mid-way between adjacent nodes, the six faces of the control volume are labelled n, , n, s, t, and b, which represent East, West, North, South, Top, and Bottom. The centroid of the element is identified by the nodal point P and its neighbours in a three-dimensional geometry, the nodes to the west, east, south, north, top, and bottom, are denoted by W, E, S, N, T, and B, respectively. The positive directions along the coordinate are also given.
2.3. Analytical Expressions for -Derivatives of Velocity
The two sides of (2) and (3) are integrated twice with respect to z, which derives the formulation as follows:with the boundary condition of u = v = 0 at the fracture wall. The z-derivatives of flow velocities at the boundary are expressed as
According to the cubic law for Newtonian flow in the fracture [32], the averaged velocities in the x- and y-directions are separately given as
Respectively substituting (11) and (12) into (9) and (10) yields the following equations:
If the general variable is introduced, (13)-(14) can be usefully written in the following form:
2.4. Discretization of the Transport Equations
The key step of FVM is the integration of governing equations over a control volume and over a time interval to yield discretized equations at its nodal point P. For the three-dimensional control volume defined in Figure 2, from t to t+Δt this gives
Taking the average velocities and as unknown, the volume integral of the rate of change term can be written as
The resulting volume integral in the convective term is converted into an integral over the entire bounding surfaces of the control volume by using Gauss’s divergence theorem. Due to the assumed condition , the convective term can be transformed into the following equations:Similarly, applying Gauss’s divergence theorem and the central differencing method as well as substituting (15)-(16) into the diffusive term, the volume integral in the diffusive term can be written as follows:The source term is approximated by means of a linear form:The fully implicit scheme is applied to conduct the time integral of Eq. (17). Then two variables F and D are defined to represent the convective mass flux per unit area and diffusion conductance at cell faces, which are calculated with the following formulae:
By substituting (18)-(22) into (17), the three-dimensional numerical calculation is transformed into two-dimensional calculation, and the corresponding discretized momentum equation can be rearranged aswherewithThe neighbour coefficients of (23) for the exponential scheme are shown as follows:
The discretized form of continuity equation for the control volume gives
Based on SIMPLE algorithm, the discretized continuity equation will be turned into a pressure correction equation. A variable, pressure correction , is firstly defined as the difference between the correct pressure and the guessed pressure. The pressure correction values of nodal point P and its neighbours are denoted as , , , , and , the subscript of which indicates the position of the node in the three-dimensional geometry.
The pressure correction equation is treated in a similar manner. Its discretized equation is written aswherewith the coefficients as follows:
From the above calculation process, it can be seen that the three-dimensional flow has been transformed into a two-dimensional problem by introducing the hypothesis of Hele-Shaw model and deducing the analytical expression for partial derivatives of velocity with respect to z. The present model taking the fracture aperture into account only needs to make the two-dimensional numerical calculation; therefore it is called the quasi-3D model.
2.5. VOF Model
The grout diffusion in the fracture is a dynamic process. The flow problem involves two immiscible phases separated by a sharp interface. One phase is the grout, and the other is usually either air or water. As the grout propagates, the interface between grout and air in a fracture gradually changes. The grout flow front can be positioned by the interface tracking technology such as MAC method, Level-Set method, and VOF method [33]. In this study, the VOF method is selected to calculate the multiphase flow. A fractional volume function, which is defined as the fractional volume of one fluid in each cell, is used to reconstruct and track the interface at each moment. The evolution of the grout position and shape depends on how the function changes with time.
The volume of one phase in each cell is calculated firstly by building a scalar function as follows:where and represent the area occupied by air and grout, respectively. The scalar function satisfies the following equation:where t is time. u and v are the velocities in x- and y-direction, respectively.
The computational domain is divided into multiple noncoincident cells. The volume of a cell is denoted as . The integral of λ(x,y,t) over each cell Ii,j is Fi,j, which can be written asThe function Fi,j is called the fractional volume function, which obeys the following convection equation:
Equation (34) is the VOF equation and its solving results can be divided into three cases: if F = 1, only the grout is present in the cell; if F = 0, the cell is full of air; if , the phase interface can be in the cell. According to the value of F, the corresponding interface can be constructed.
It is obvious that the key to tracking the moving interface is solving the VOF equation [34]. There are many methods at present, such as the Hirt-Nichols method, the Flair method, the Youngs method, and the FCT method. The Youngs method based on geometric principles has higher computational accuracy compared with other methods [33, 35]. The basic implementation steps of the Youngs method include (1) constructing the fluid interface according to the value of F in the cell; (2) advancing the moving interface (updating the VOF function) according to the current flow velocity. Repeatedly execute the above two steps in each time step to track the movement of the interface between grout and air.
2.6. Numerical Details
The overall numerical procedure for predicting the grout diffusion in the fracture is as follows:(1)Fractional volume function F, velocity vector u, and pressure p are initialized at all cells and the boundary conditions are defined.(2)Material properties at all cells are initialized according to the time t and fractional volume function F.(3)Momentum equations are solved to obtain the current velocity.(4)Pressure correction equation is solved to obtain pressure correction value.(5)Update current velocity and pressure.(6)Repeat steps (3)-(5) until the flow field converges.(7)Update fractional volume function using Youngs method.(8)Return to step (2). Time step is advanced and the numerical process proceeds to the next iteration.
3. Results and Discussion
3.1. Validation of Numerical Model
3.1.1. Case 1: h = 2 mm
In order to validate the proposed model, an example of grouting a parallel fracture with a constant rate was taken. Its analytical solution is listed in Appendix A.
The fracture is 1 m long and 1 m wide with an aperture of 2 mm, as shown in Figure 3. A no-slip boundary condition is used for upper and lower wall surfaces and the pressure outlet boundary for four sides. A grouting hole, 10 mm in radius, is located at the center of upper wall. It is assumed that the fracture is initially filled with air under ambient pressure of zero. The grout with a density of 1500 kg/m3 and viscosity of 0.04 Pa·s is pumped through the grouting hole at a constant flow rate of 37.7 mL/s. Air viscosity at normal temperature is 1.8×10−5 Pa·s and its density is 1.205kg/m3. Considering the geometric symmetry of the fracture model, the simulation was simplified by modeling only a quarter of the flow field, as shown in Figure 4. Case 1 was computed separately by the present quasi-3D numerical model and the commercial code ANSYS-Fluent. In the new model the element length in the z-direction equaled the aperture, while the fracture was decomposed into 2 layers and 6 layers in the direction perpendicular to the wall before the calculation is performed using Fluent. All the calculations were finished by the same computer with the following basic configuration:(i)CPU: Intel Core i7 7700(ii)Mainboard: ASUS B150-PLUS(iii)Memory: 16 GB (DDR4 2401MHz), 400
In the quasi-3D numerical model, the computational domain was meshed into 2500 structured hexahedral elements at the interval of 10 mm along both length and width. The edge length in the z-direction was 2 mm, and the time step was set to 0.1 s. The calculation kept going until the front flow of grout reached the boundary. The total computing time was about 0.5 hours.
Simultaneously, the grout diffusion behavior in the fracture was simulated by ANSYS-Fluent 17.0 under the same boundary condition. The simulated domain was divided into 500064 hexahedral cells with an edge length of 1 mm in the z-direction. Most of the meshes had an edge length between 1 mm and 3.5 mm. A variable time step size limited by global Courant number was adopted, the maximum time step was less than 0.006 s, and the calculation finished after 48 hours. Correspondingly, when the fracture aperture was divided into 6 layers, the total number of cells was 1500264. The minimum edge length is 0.3 mm, the time step was less than 0.003 s, and the total computing time was more than 432 hours.
Figures 5–7 display the advancement of the grout flow front as a function of time in the fracture. The diffusion front positions of grout at t = 10 s, 20 s, 30 s, and 40 s are represented based on three different models. The results show that the filling range of grout gradually increases in the fracture with the continuously grouting. Besides, the diffusion shape is always circular, which is consistent with the diffusion pattern described by the analytical model.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The numerical results and analytical solution of the grout diffusion radius are shown in Table 1. According to these data, the variation curves of grout diffusion radius with time are plotted in Figure 8. The average relative errors of the Quasi-3D model and 2 layers and 6 layers model are 0.41%, 0.37%, and 0.38%, respectively. The numerical predictions of the three models are all in good agreements with the analytical solution for diffusion radius variation with time.
The grout pressure at different positions calculated by the three models is presented in Table 2 at t = 20 s and Table 3 at t = 40 s. The corresponding average flow velocity is summarized in Tables 4 and 5. The contrastive pressure distribution over radius is plotted in Figure 9. The comparative study for the distribution of average flow velocity at different times is illustrated by Figure 10. It shows that the presented model agrees with the analytical solution better than other models in the distributions of the pressure and velocity. For the new model, the maximum relative errors of the pressure and velocity are 4.44% and 1.91% at t = 20 s, respectively, while 4.25% and 1.90% at t = 40 s, respectively. By contrast, 6 layers model calculated by ANSYS-Fluent has a lower accuracy. At t = 20 s, the relative errors of the pressure and velocity are up to 7.30% and 5.63% separately. At t = 40 s, they are 7.52% and 5.77%, respectively. Then the solution accuracy of 2 layers model calculated by ANSYS-Fluent is further reduced. The maximum relative error is more than 33%.
(a)
(b)
(a)
(b)
In general, when grouting at a constant rate, both the grout pressure and the velocity decrease with the increase of the radius, and the rate of change declines gently. Their maximum values are located in the grouting hole. Obviously, as the diffusion radius increases, the grouting pressure rises up. The trend is in accord with the analytical model.
Table 6 shows an overview of the quasi-3D model and 2 layers model and 6 layers model about the implementation details and the average relative errors of solving results for the diffusion radius, the pressure distribution, and the average velocity. It is obvious that the quasi-3D model takes less time but obtains a higher solving precision. The new model has more advantage in simulating the grout diffusion in the fracture than others.
3.1.2. Case 2: h = 0.1 mm
The second verification example is to calculate the diffusion characteristics of grout in a parallel-plate fracture under constant grouting pressure, and its analytical solution is listed in Appendix B.
In Case 2, the grout was pumped into a fracture of size 1 m × 1 m × 0.1 mm with a constant grouting pressure of 200 kPa. As Case 1 has mentioned, the fracture was still initially full of air under ambient pressure of zero. The density and viscosity of the grout are 1500kg/m3 and 0.02 Pa·s, respectively. Other conditions were identical with Case 1. Still only a quarter of the fluid region was calculated. For the quasi-3D model, the number of the grid unit is 2500.
Figure 11 shows the change of grout diffusion front position with time under the grouting pressure of 200 kPa. The interface between the grout and air gradually goes forward away from the grouting hole in a circular way as the injection keeps going, which is the same as the diffusion behavior obtained by the analytical model.
(a)
(b)
(c)
(d)
(e)
Table 7 presents the comparison of the pumping rate and grout diffusion radius between the simulation results and the analytical solution, and their contrast curves are plotted in Figures 12 and 13, respectively. It can be seen that the diffusion radius gradually increases over time while the diffusion rate descends slowly under constant grouting pressure. The results of the quasi-3D model agree well with the analytical solution. The relative error of the diffusion radius is less than 1.25%, and of the pumping rate is 1.6%.
Tables 8 and 9 show the contrast of the radial pressure and average velocity calculated by the quasi-3D model and the analytical model at t = 20 s and t = 50 s, respectively. The distributions of the pressure and average velocity at different moments are plotted in Figures 14 and 15, respectively. It demonstrates that the developed numerical model matches well with the analytical result. The maximum average relative error of the pressure is about 1.91%, and of the average velocity is less than 1.95%. Figure 15 indicates that the average flow velocity in sections declines over the diffusion radius, which is due to the decrease in the pressure gradient of the flow field and pumping rate. To summarize, at different times the results of the pressure and average velocity obtained by the numerical model are in good agreement with the theoretical solution.
(a)
(b)
(a)
(b)
3.2. Comparison with Experimental Results
To further verify the applicability of the quasi-3D model, an experiment of grout injection in a fracture was carried out. The main device consists of a transparent 12 mm thick plexiglass plate, a 20 mm thick steel plate and 1mm thick rubber gasket in the middle of both to form the aperture, as shown in Figure 16. The grout with a viscosity of 0.2 Pa·s was injected through a 1 cm diameter grouting pipe mounted on the bottom of steel plate at a pumping rate 5 mL/s. The grout diffusion process was taken photos by a digital camera.
(a)
(b)
The quasi-3D model, set the same condition as experimental, was used to simulate the grout diffusion process. Due to the symmetric structure of the fracture, a quarter of simulated domain in the upper right was selected to implement the numerical work, as shown in Figure 17. The domain was divided into 50 cells in the x-direction and 30 cells in the y-direction.
Figure 18 shows the grout diffusion region and morphology at t = 5 s, 10 s, 15 s, 20 s, 25 s, and 30 s in the experiment. The grout diffuses in a circular shape in the fracture, which is consistent with the theoretical model. The grout front advances forward, then the filling area becomes larger over the time. Figure 19 displays the comparison of the grout diffusion region and morphology between numerical and experimental results. The solid line represents the experimentally measured front position of the grout. It indicates that the results of both methods are basically in agreement. The grout diffusion radius was obtained by averaging the measured values in four orthogonal directions. Table 10 shows the grout diffusion radius at different times obtained from the analytical, experimental, and numerical methods. The average relative error of the numerical results versus the analytical solution is 0.64%. The average relative error between the experimental results and the analytical solution is 0.86%. To summarize, the results of the grout diffusion distance at different times calculated by the quasi-3D model are in good accordance with the experimental and analytical results, demonstrating that the presented model can accurately simulate the grouting diffusion process
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusion
Numerical simulation is always an important way to study the diffusion mechanism of grout in the rock fracture. Based on the traditional three-dimensional FVM, the aim of this paper is to develop a quasi-3D model to simulate the grout injection in the parallel fracture considering its geometric characteristics. The present model is validated by the analytical model and the experimental measurement. The main work and conclusion are summarized as follows:(i)The geometric feature of fracture is that the length and width are much larger than the aperture. In the presented model, the average velocity parallel to the fracture surface was used as the unknown variable to discretize the governing equations. The partial derivatives of the velocity with respect to z at fracture walls were deduced by introducing the Hele-Shaw model. The derived equations were proven to be a function of the aperture and the average velocity. By substituting the expressions into the discretized momentum equations, the three-dimensional diffusion of the grout in the parallel fracture can be transformed into a two-dimensional problem. Combined with the VOF method of tracking the grout front, a quasi-3D simulation method was established to predict the flow and diffusion behavior of grout in the fracture.(ii)The quasi-3D model was validated by two cases under different conditions. It was also compared with the experimental results. The calculation results, including the diffusion radius, shape, pressure, and average velocity, show a good agreement with the analytical and experimental results. The contrast of the present model with the simulation by ANSYS-Fluent indicates that the new model needs less grid unit and spends much less time when they both satisfy the same accuracy requirement.(iii)The proposed model provides an efficient way to simulate the grout injection in the fracture. The process of grout diffusion under different environmental conditions and the influences of different factors on the grout, especially the quick-setting grout, will be focused on in the future study.(iv)This method succeeds to simulate the grout diffusion behavior in a fracture. On this basis, the quasi-3D model may be further extended to a three-dimensional intersected fracture network in the next phase; thus it is possible to simulate the grout injection in a more complex environment [36, 37].
Appendix
A. Constant Flow
The analytical expressions for grout pressure distribution, diffusion radius, and average flow velocity in different sections are given as follows, respectively [38]:where t is time and R is the diffusion radius at time t. is the kinematic viscosity of grout. q is the pumping rate. h is the fracture aperture. p is grout pressure at radius r, p0 is the ambient pressure, and u is the radial average velocity.
B. Constant Pressure
When the grouting pressure is constant, the analytical solutions of the grout diffusion radius, pumping rate, pressure distribution, and the average velocity are shown below, respectively [38]:where t is time and R is the diffusion radius at time t. is the kinematic viscosity of grout. q is the pumping rate. h is the fracture aperture. p is grout pressure at radius r, P0 is the grouting pressure, p0 is the ambient pressure, and u is the radial average velocity. r0 is the radius of the grouting hole.
Symbols
| : | Face area |
| : | Coefficient of general variable |
| : | Constant term in the momentum discretized equation |
| : | Source term in the pressure correction equation |
| : | Diffusion conductance at cell faces |
| : | Coefficient of momentum interpolation equation |
| : | Convective mass flux per unit area |
| : | Fractional volume function |
| : | Fracture aperture |
| : | Grouting pressure |
| : | Grout pressure |
| : | Ambient pressure |
| : | Pressure correction value |
| : | Pumping rate |
| : | Diffusion radius |
| : | Radius |
| : | Radius of the grouting hole |
| : | Source term |
| : | Constant term in the linearised form of the source term |
| : | Source term of -momentum |
| : | Source term of -momentum |
| : | Source term of -momentum |
| : | Slope of the function of source term with respect to at node P |
| : | Time |
| u: | Velocity vector |
| : | Velocity in the -direction |
| : | Average velocity in the -direction |
| : | Velocity in the -direction |
| : | Average velocity in the -direction |
| : | Velocity in the -direction |
| : | Coordinate value |
| : | Diffusion coefficient |
| : | Time interval |
| : | Volume of a grid cell |
| : | Volume of a cell |
| : | Edge length of a grid cell |
| : | Node spacing in the -direction |
| : | Node spacing in the -direction |
| : | Node spacing in the -direction |
| : | Scalar function |
| : | Viscosity |
| : | Density |
| : | General variable in the traditional model |
| : | General variable at node P at new time level |
| : | General variable at node P at (old) time level |
| : | Area occupied by air |
| : | Area occupied by air. |
| : | Node to the bottom |
| : | Bottom face |
| : | Node to the east |
| : | East face |
| : | Node to the north |
| : | North face |
| : | Nodal point P |
| : | Between node P and node E |
| : | Between node P and node N |
| : | Node to the south |
| : | Between node S and node P |
| : | South face |
| : | Node to the top |
| : | Top face |
| : | Node to the west |
| : | Between node W and node P |
| : | West face. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest to this work.
Acknowledgments
This research was supported by National Key Research and Development Plan [Grant no. 2016YFC0802207], the Program for Science and Technology Innovation Talents in Universities of Henan Province [Grant no. 19HASTIT041], the National Natural Science Foundation of China [Grants nos. 51878624 and 51878622], Key Research Projects of Higher Education in Henan Province [Grant no. 18A580001], the Natural Science Foundation of Henan province [Grant no. 182300410116], the Program for Innovative Research Team (in Science and Technology) at the University of Henan Province [Grant no. 18IRTSTHN007], and Major Scientific, Technological Special Project in Henan [Grant no. 181100310400], and the Postdoctoral Research Sponsorship in Henan Province [Grant no. 19030026].