Unsplit Youngs Method for Tracking the Moving Interface of Expansible Grout

Accurately tracking the moving interface is one of the most critical challenges in numerically modeling grout diffusion in fractures. The traditional Youngs method does not consider fluid transport to the oblique adjacent cell, which leads to repeated calculation and error accumulation. Based on the interface reconstruction technology, this paper introduced the idea of the unsplit algorithm to capture the grout-air interface. The unsplit interface advancement formulas were derived considering the fluid fluxes from the central cell to the oblique adjacent cells in the structured cell system. The numerical method based on the unsplit algorithm was validated by examples of predicting the propagation of expansible grout in a fracture. Results show that the diffusion pattern and grout front position calculated by the proposed method are in better agreement with the analytical solution compared with the traditional split algorithm. The new method can effectively improve the grout-air interface tracking accuracy.


Introduction
Grouting has been widely used to reinforce and control the seepage of rock mass in the field of infrastructure construction, such as traffic and water conservancy [1][2][3]. With the development of modern macromolecular science, polymer grouting material with the characteristic of selfexpansion and its high-pressure injection technology have more extensive application in recent years [4][5][6]. When the two-component polymer grout is injected into the gap or fracture, it rapidly expands and solidifies due to the chemical reaction. e excellent properties can be utilized to fill the gap, compact the soil, strengthen the rock mass, and seal the flow channels. e diffusion of polymer grout mainly depends on the secondary pressurization caused by self-expansion force, which is produced by the chemical reaction. is is different from the nonexpansive grout which is driven by the grouting pressure [7,8]. e expansion property of the polymer is influenced by environmental pressure, the extent of reaction, and other factors. erefore, the grouting mechanism of polymer grout is more complex compared with the traditional nonexpansive grout. Although great progress has been made in grouting equipment, technology, and engineering practices, the study on the grouting mechanism still lags significantly behind it due to the concealment of grouting construction and the particularity of polymer material.
Numerical simulation is an important way to study the diffusion mechanism of grout in rock fracture. Hassler et al. established a two-dimensional fracture network model to predict the grout diffusion process [9]. Eriksson et al. further studied the impact of filtration and fracture aperture on the diffusion behavior of Bingham grout in the fracture network [10]. Based on the Monte Carlo method, the two-dimensional or three-dimensional discrete fracture network was generated randomly to investigate the grouting mechanism [11]. e abovementioned studies construct fracture networks in a deterministic or random way based on the one-dimensional channel model, and the grout flow is assumed as a monophasic flow. However, the grout may not be a one-dimensional flow in a large-size fracture. erefore, the replacement between grout and air (or water), which reflects the movement of the grout front, should be paid attention.
In recent years, considering the grout diffusion as multiphase flow, the study on the refined simulation of grout diffusion behavior has attracted much attention. Liu et al. simulated the change of grout front with time and pressure distribution in a fracture by COMSOL software [12]. Similarly, Zhao [13] explored the diffusion performance of cement grout in the three-dimensional cross fracture. Accurate interface tracking is crucial for the prediction of the grout front. Li [14,15] simulated expansible grout spread in a two-dimensional flow field based on the split Youngs algorithm. However, the split algorithm only calculated the fluid fluxes from the central cell to the top, bottom, left, and right cells, which do not fit the actual fluid flow.
In order to capture the moving interface more accurately, the unsplit algorithm was developed [16,17]. Depending on the adopted coordinate system and unit fluid form, the pattern and ways of implementation can be different. e new algorithm considers the additional fluid transport to the oblique adjacent cell and avoids the problem of the repeated calculation of the fluid flux in the traditional Youngs method.
In this paper, an unsplit Youngs method is used to track the moving interface of expansible polymer grouting a fracture. Assuming that the unit fluid is a trapezoid in the cell containing the interface, the expressions of fluid flux from the central cell to the oblique adjacent cells are derived under the structured cell system. en, the interface between grout and air is advanced based on the fractional volume. e diffusion process of expansible grout in a fracture is studied by the unsplit algorithm. Compared with the analytical solution and the traditional split Youngs algorithm, the unsplit Youngs method significantly improves the accuracy of tracking the moving interface.

Volume of Fluid Method
Tracking the moving interface is required when the flow of two-phase or multiphase fluid is modeled. Currently, the interface tracking method includes the marker and cell (MAC) method, Level-Set method, and Volume of Fluid (VOF) method [18][19][20], among which the VOF method is the most widely used. It can guarantee the sharp, nonoscillatory interface that is desirable in immiscible multifluid flow simulations and essential in free-surface problems on stationary meshes. e fractional volume of one fluid in each cell is expressed by the function of the volume of fluid. e interface is reconstructed by the change of the volume of fluid with time; thus, the position, shape, and change of the moving interface can be tracked.
In a two-dimensional system, the calculation area is assumed as Ω. e area where fluid A locates is Ω 1 , and the area where fluid B locates is Ω 2 . A scalar function is defined as follows: When the flow field consists of two immiscible fluids, the scalar function λ (x, y, t) follows continuity equation: where t is time, u is the velocity in X, and v is the velocity in Y. e calculation area is divided into multiple noncoincident cells, ΔV ij is the volume of a cell, and the integral of λ(x, y, t) on one cell can be written as follows: F is the fractional fluid volume and follows the continuity equation: e above equation is the VOF equation. e volume fraction F is solved by equation (4). ere are three cases: when F � 1, the fluid occupies the entire cell; when F � 0, the cell is empty; and when 0 < F < 1, there is two-phase fluid in the cell. Based on the value of F, the corresponding interface can be constructed.

Youngs Method
e key to tracking the moving interface is to solve the VOF equation [21]. ere are many methods at present, such as the Hirt-Nichols method, the Flair method, the Youngs method, and the FCT method. However, the Youngs method based on geometric principle has higher accuracy compared with other methods [19,20]. e basic implementation steps of the Youngs method include (i) Constructing the fluid interface according to the fractional volume in the cell (ii) Advancing the moving interface by updating the VOF function according to the current flow velocity (iii) Repeatedly executing the above two steps in each time step.

Construction of the Fluid Moving
Interface. e Youngs method constructs the approximate interface with a sloped line segment in a cell base on the geometric principle. As mentioned above, there are three cases in the calculated cell, namely, (1) e cell is full of fluid A, F � 1 (2) e cell is an interface unit, 0 < F < 1 (3) e cell is filled with the fluid B, F � 0 2 Complexity Case 1 and Case 3 do not reconstruct the interface. For Case 2, the normal vector of the interface in the cell is calculated firstly, as shown in en, the angle β that the interface makes with the x-axis is determined according to the normal vector: e angle β is normalized to ω: where δx is the length of a fluid cell in the x-direction and δy is the length of the cell in the y-direction, as shown in Figure 1.
According to the angle ω and VOF function in the cell, there are 16 types of interface combinations in the cell, which can be simplified into four cases by symmetry and inversion. e four types of interfaces are shown in Figure 2. e discriminant method is as follows.
Once the case has been determined, the four side fractions of the cell, which is the fraction of the top, bottom, left, and right sides of the cell that lie within the fluid (s t , s b , s l , s r ) can be calculated according to the angle ω and the VOF function. e subscripts t, b, l, and r represent up, down, left, and right, respectively, as shown in Figure 2(d) below.
Taking case III as an example, the side fractions of the cell can be calculated as follows: After obtaining the side fractions of each cell, the straight-line segment that represents the fluid interface in the cell can be plotted. Function). According to the current velocity, the fluxes through sides of the cell to adjacent cells in a time step are calculated.

Advancement of the Moving Interface (Update VOF
e VOF function of each cell is updated by accumulating the inflow and the outflow. e moving interface is determined and advances according to the updated VOF function. e Youngs method was proposed in 1982, but very few details of flux calculation were given in the original paper [22]. For this reason, Rudman [20] deduced the relevant formulas by the direction split method. From then on, the Youngs method began to draw attention and be used widely.

Split Youngs Algorithm
(1) Calculate Fluxes. Considering the horizontal and vertical flow of the fluid through the interface of the cell, the fluxes from the central cell to the adjacent ones, including the top, bottom, left, and right, are calculated based on the split Youngs algorithm, represented by f t , f b , f l , and f r , respectively, as shown in Figure 3.
Taking case III (in Figure 2) as an example, according to the geometric law, the fluxes can be calculated as follows.
(1) Top face: 2 tan β else f r � Fδxδy where u t , u b , u r , and u l represent the interface velocity in the top, bottom, left, and right, respectively. t is the time step. e fluxes of cases I, II, and IV in four directions can be calculated in the same way.

Complexity
(2) Update VOF Function. e fractional volume in the cell can be updated by calculating the fluxes from the central cell to adjacent ones during the time δt.
Given the volume fraction of the fluid in the cell (i, j) is F t at time t, after the time δt, the increments of fractional volume in the cells (i, j + 1), (i, j − 1), (i − 1, j), (i + 1, j) caused by the fluid motion are f t , f b , f l , and f r , respectively. e volume fraction in the cell at time t + δt is e fractions of fluid volume in the four adjacent cells at time t + δt are, respectively, corrected as follows: Conducting the same calculation for each cell can update the fractional volume of the flow field. e computing process of fluid fluxes shows that the split Youngs algorithm just includes fluxes to four adjacent cells (top, bottom, left, and right) from the central cell. e fluid flow towards the oblique cell is ignored, which may lead to some deviations from the real.

Unsplit Youngs Algorithm.
In order to improve tracking accuracy, an unsplit operator algorithm was introduced to advance the interface. e method not only calculates the fluxes from the central cell to the top, bottom, left, and right cells but also considers the transport to four oblique direction cells including upper right, upper left, lower left, and lower right, as shown in Figure 4. e fluxes in eight directions are denoted as C t , C b , C l , C r , C cor1 , C cor2 , and C cor3 , C cor4 , and the subscripts cor1, cor2, cor3, and cor4 represent the upper right corner, the upper left corner, the lower left corner, and the lower right corner, respectively.
Hu et al. [16] assumed the unit fluid form is a triangle in the cell containing the interface. Five cases were listed when U t > 0 and U r > 0. e relationship between the above eight fluxes was described and f t , f b , f l , and f r were calculated. However, the analysis of the pattern of fluid flow is less comprehensive, and the flux calculation for different cases has not been explained in detail. Liu et al. [17] studied the unsplit transport method in the Eulerian method. e expressions of fluid flux were given in the cylindrical coordinate system. We further analyze the unsplit transport pattern and derive the formulas of interface advancement based on trapezoid unit fluid (see Figure 2(c)) in rectangular coordinates. e implementation is shown as follows.
(1) Fluid Fluxes in the Orthogonal Directions. e fluxes from the central cell to four adjacent cells were supposed as C t , C b , C l , and C r . Figure 5 shows the pattern of fluid transport from the central cell ( Figure 5 It indicates that the fluxes in four oblique directions need to be calculated first.  Figure 6 (2), (5)) else (two cases as Figure 6 (4), (7)) In the same way, the flux expressions of case III during the time δt from the central cell to the upper left ( Figure 5(c)), lower left ( Figure 5(d)), and lower-right corner ( Figure 5(e)), C cor2 , C cor3 , and C cor4 , can be written as follows.
For C cor2 , Similarly, the flux expressions of case I, case II, and case IV in four oblique directions can be solved.   Complexity j + 1), F t (i − 1, j − 1), and F t (i + 1, j − 1), respectively, at time t + δt, they are Repeat the calculation for each cell; then, the update for the VOF function in the whole field can be completed.

Case Studies
e expansive diffusion process of expansible grout in a fracture was tracked to validate the precision of the unsplit Youngs algorithm. e results were compared with the analytical solution and the traditional split Youngs method.

Calculation Model.
e large-size and well-developed fracture can be modeled as a parallel-plate fracture. Four side edges of the fracture are set as the pressure outlet, as shown in Figure 8. e length and width of the fracture are both 1.6 m and the fracture aperture is 3 mm. A grouting hole is located in the center of the model, through which a certain amount of expansible grout is injected into the fracture. Assuming that the fracture is initially filled with air, the injected grout will partially displace it expansible and occupy a circle area with the radius R 0 . According to Mitani and Hamada [23], the grout density changes with time as below: where ρ 0 is the initial grout density with a value of 1100 kg/ m 3 , α is the expansion rate, t is the time, and ρ t is the grout density at time t. e grout viscosity remains unchanged at 1.0 pa·s before gelation. Considering the fracture model is geometrical symmetry, the solving could be simplified by only calculating a quarter of the simulated region, as shown in Figure 9. erefore, there are two symmetry boundaries and the other boundary conditions are pressure outlet. e grouting quantity is represented by the initial diffusion radius of expansible grout in a fracture. e grouts with different components have different expansion rates. e input parameters are shown in Table 1. Two types of expansible grout with different initial radius were used to validate the unsplit Youngs algorithm. e analytical solution of the diffusion distance of the expansible grout in a fracture is expressed as follows [24]: where R 0 is the initial diffusion radius and R t is the diffusion radius at time t.

Numerical Calculation.
e diffusion behavior of the expansible grout can be described by one continuity equation and three momentum conservation equations. e continuity equation is given as 8 Complexity where F is the volume fraction of the expansible grout and U is the velocity vector.
e two phases of grout and air share one single set of momentum equations, which are written as follows:  Complexity where g is the gravity acceleration, ρ and μ are the volume fraction weighted average of the density and viscosity, respectively, and p is the pressure. e solution domain was meshed into 6400 hexahedral elements with an edge length of 0.01 m. Based on the finite volume method, the governing equations were solved iteratively by the Semi-Implicit Method for Pressure Linked Equations algorithm (SIMPLE). e procedure is listed below: (1) Define the boundary conditions and initialize the fractional volume of each cell, velocity vector, and pressure (2) Determine the material properties of expansible grout according to the time and the VOF function (3) Solve the discretized momentum equation to obtain the current velocity (4) Solve the pressure correction equation to obtain the pressure correction value (5) Update the current velocity and pressure (6) Repeat Steps (3)-(5) until the flow field converges (7) Use the unsplit Youngs algorithm to update the VOF function and advance the interface Return to step (2), update the time, and the numerical process proceeds to the next iteration. e numerical calculation was programmed by FORTRAN. e procedure is shown in Figure 10 However, centered on the grouting hole, the grout is circular diffusion and spreads symmetrically towards the surroundings, which is consistent with the theoretical diffusion pattern. e grout diffusion distance calculated by the unsplit Youngs algorithm is in good agreement with the analytical solution. e contour of the grout front calculated by the two methods almost completely coincides. However, the calculation result from the split algorithm brings a significant deviation from the analytic solution. e diffusion ranges calculated by the split algorithm at different times lag behind the theoretical results. e grout contour calculated by the split algorithm always lies inside the analytical solution, and the deviation increases gradually with time. e average diffusion radius at different times was calculated for Cases 1-4, shown in Tables 2-5, and plotted in Figures 15(a)-15(d), respectively. Obviously, the diffusion radius calculated by the unsplit Youngs algorithm always agrees well with the analytical solution with the increasing time, and their contours of the grout front almost coincide. From Case 1 to Case 4, the maximum absolute errors are 0.92 cm, 1.1 cm, 0.41 cm, and 0.26 cm, respectively. By contrast, the diffusion radius calculated by the split algorithm is always smaller than the analytical solution, and the deviation gradually increases with time.
e maximum absolute errors are 10.71 cm, 7.99 cm, 4.05 cm, and 4.13 cm, respectively. Figure 16 shows the relative error between the numerical method and analytical solution with respect to the diffusion radius under the four cases. e maximum relative errors calculated by the unsplit algorithm are 2.14% and 2.37% and 1.21% and 0.85%, respectively, while they are 13.71%, 10.68%, 5.27%, and 5.48% calculated by the split algorithm, respectively, which are 6.4 times, 4.5 times, 4.36 times, and 6.45 times, respectively, higher than the unsplit algorithm. Besides, it can be observed that the relative error of the split algorithm is gradually increasing with the increase of the Complexity diffusion radius, while the relative error of the unsplit algorithm declines slightly. e above comparison of relative errors shows that the numerical results of the split Youngs algorithm deviate from the analytical solution more and more as the diffusion distance increases. Because the split Youngs algorithm does not consider fluxes in the diagonal direction, the arithmetic errors will be accumulated. By comparison, the unsplit Youngs algorithm better guarantees the conservation of fluid mass in the flow by adding the fluxes in the diagonal  dimensional Youngs method has to take into account the fluid flux in the z-direction, which needs to reprogram algorithm. is paper is our preliminary result based on the idea of the unsplit interface. We will further develop the unsplit algorithm of the moving interface to make it more widely applicable.

Conclusion
How to track the moving interface accurately is crucial to predict the grout diffusion behavior. Considering the fluxes from a central cell to oblique adjacent ones, unsplit interface advancement formulas were deduced under the structured cell system. Based on the traditional split Youngs algorithm, the unsplit algorithm improves the accuracy of tracking the moving interface by introducing the basic principle of the direction unsplit method. e diffusion process of expansible grout in a parallel-plate fracture is studied by the proposed algorithm and split Youngs algorithm. e results show that the diffusion distance calculated by the unsplit Youngs algorithm agrees better with the analytical solution under different cases. Compared with the traditional split Youngs algorithm, the new method can reduce the calculation error effectively and improve the tracking accuracy of the grout front significantly. is study can be used to predict the grout front in a fracture, so the arrangement of grouting hole can be designed reasonably. It has been used in the repair and maintenance of transportation infrastructure, e.g., repair for cavity beneath cement concrete pavement. On this basis, we will try to apply the unsplit Youngs algorithm to the tracking of the grout moving interface in fracture network in the future; thus, it is likely to capture the grout diffusion process in more complex geology.

Data Availability
All data in our study are available from the corresponding author by request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.  16 Complexity