# Unsplit Youngs Method for Tracking the Moving Interface of Expansible Grout

Authors: Xiaolong Li; Tiantian Zhang; Meimei Hao; Bei Zhang; Wenyi Yao; Zhenzhou Shen
Journal: Complexity
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–3]. With the development of modern macromolecular science, polymer grouting material with the characteristic of self-expansion and its high-pressure injection technology have more extensive application in recent years [4–6]. When the two-component polymer grout is injected into the gap or fracture, it rapidly expands and solidifies due to the chemical reaction. The excellent properties can be utilized to fill the gap, compact the soil, strengthen the rock mass, and seal the flow channels.The diffusion of polymer grout mainly depends on the secondary pressurization caused by self-expansion force, which is produced by the chemical reaction. This is different from the nonexpansive grout which is driven by the grouting pressure [7, 8]. The expansion property of the polymer is influenced by environmental pressure, the extent of reaction, and other factors. Therefore, 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]. The 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. Therefore, 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. The 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. Then, the interface between grout and air is advanced based on the fractional volume. The 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.


## 2. 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–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. The fractional volume of one fluid in each cell is expressed by the function of the volume of fluid. The 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 Ω. The area where fluid A locates is Ω1, and the area where fluid B locates is Ω2. A scalar function is defined as follows:(1)λx,y,t=1,x,y∈Ω1,0,x,y∈Ω2.When the flow field consists of two immiscible fluids, the scalar functionλ (x, y, t) follows continuity equation:(2)∂λ∂t+u∂λ∂x+v∂λ∂y=0,where t is time, u is the velocity in X, and v is the velocity in Y.The calculation area is divided into multiple noncoincident cells,ΔVij is the volume of a cell, and the integral of λx,y,t on one cell can be written as follows:(3)Fij=1ΔVij∫Iijλx,y,tdV.F is the fractional fluid volume and follows the continuity equation:(4)∂F∂t+u∂F∂x+v∂F∂y=0.The above equation is the VOF equation. The volume fractionF is solved by equation (4). There 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.


## 3. Youngs Method

The key to tracking the moving interface is to solve the VOF equation [21]. There 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].The 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.


### 3.1. Construction of the Fluid Moving Interface

The 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)
The cell is full of fluidA, F = 1(2)
The cell is an interface unit, 0 <F < 1(3)
The cell is filled with the fluidB, F = 0Case 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(5)nijx=Fi+1,j+1+2Fi+1,j+Fi+1,j−1−Fi−1,j+1−2Fi−1,j−Fi−1,j−1δx.nijy=Fi+1,j+1+2Fi,j+1+Fi−1,j+1−Fi+1,j−1−2Fi,j−1−Fi−1,j−1δy.Then, the angleβ that the interface makes with the x-axis is determined according to the normal vector:(6)tanβ=−nxny.The angleβ is normalized to ω:(7)tanω=tanβδxδy,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.Figure 1
Schematic of the fluid cell.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. The four types of interfaces are shown in Figure 2. The discriminant method is as follows.ifω < π/4ifF ≤ 0.5 tanωcase Ielse ifF ≤ 1.0 − 0.5 tanωcase IIelsecase IVelseifF ≤ 0.5 cotωcase Ielse ifF ≤ 1.0 − 0.5 cotωcase IIIelsecase IVFigure 2
Four types of interface.
(a)(b)(c)(d)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 (st, sb, sl, sr) can be calculated according to the angle ω and the VOF function. The 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:(8)st=F−0.5cotωsr=1sb=F+0.5tanωsl=0.After obtaining the side fractions of each cell, the straight-line segment that represents the fluid interface in the cell can be plotted.


### 3.2. Advancement of the Moving Interface (Update VOF Function)

According to the current velocity, the fluxes through sides of the cell to adjacent cells in a time step are calculated. The VOF function of each cell is updated by accumulating the inflow and the outflow. The moving interface is determined and advances according to the updated VOF function.The 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.


#### 3.2.1. 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 ft, fb, fl, and fr, respectively, as shown in Figure 3.Figure 3
Schematic of fluid transport calculation for split Youngs algorithm.Taking case III (in Figure2) as an example, according to the geometric law, the fluxes can be calculated as follows.(1)
Top face:ifut > 0ft=utδtstδx+0.5utδxcotβ(2)
Bottom face:ifub < 0fb=absubδtsbδx−0.5absubδtcotβ(3)
Left side:iful < 0if abs (ul)δt ≤ (1 − sb)δxfl=0elseifabsulδt≤1−stδxfl=0.5absulδt−1−sbδx2tanβelsefl=absulδtδy−1−cδxδy(4)
Right side:ifur > 0ifurδt≤stδxfr=urδtδyelseifurδt≤sbδxfr=urδtδy−0.5urδt−stδx2tanβelsefr=Fδxδywhereut, ub, ur, and ul represent the interface velocity in the top, bottom, left, and right, respectively. t is the time step. The fluxes of cases I, II, and IV in four directions can be calculated in the same way.(2) Update VOF Function. The 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 Ft 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 ft, fb, fl, and fr, respectively. The volume fraction in the cell at time t + δt is(9)Ft+δti,j=Fti,j−ft+fb+fl+frδxδy.The fractions of fluid volume in the four adjacent cells at timet + δt are, respectively, corrected as follows:(10)Ft+δti,j+1=Fti,j+1+ftδxδyFt+δti,j−1=Fti,j−1+fbδxδyFt+δti+1,j=Fti+1,j+frδxδyFt+δti−1,j=Fti−1,j+flδxδy.Conducting the same calculation for each cell can update the fractional volume of the flow field.The 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. The fluid flow towards the oblique cell is ignored, which may lead to some deviations from the real.


#### 3.2.2. Unsplit Youngs Algorithm

In order to improve tracking accuracy, an unsplit operator algorithm was introduced to advance the interface. The 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 Figure4. The fluxes in eight directions are denoted as Ct, Cb, Cl, Cr, Ccor1, Ccor2, and Ccor3, Ccor4, 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.Figure 4
Schematic of fluid transport calculation for unsplit Youngs algorithm.Hu et al. [16] assumed the unit fluid form is a triangle in the cell containing the interface. Five cases were listed when Ut > 0 and Ur > 0. The relationship between the above eight fluxes was described and ft, fb, fl, and fr 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. The 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. The implementation is shown as follows.(1) Fluid Fluxes in the Orthogonal Directions. The fluxes from the central cell to four adjacent cells were supposed as Ct, Cb, Cl, and Cr. Figure 5 shows the pattern of fluid transport from the central cell (Figure 5(a)) to the upper right cell (Figure 5(b)), the upper left cell (Figure 5(c)), the lower left cell (Figure 5(d)), and the lower right cell (Figure 5(e)) from tn to tn+1.Figure 5
Schematic of fluid transport from the central celli, j to the four oblique cells: (a) the position at time tn, (b) the position at time tn for case 1, (c) the position at time tn+1 for case 2, (d) the position at time tn+1 for case 3, and (e) the position at time tn+1 for case 4.
(a)(b)(c)(d)(e)According to Figures5(b)–5(e), the equations can be expressed as follows:(11)Ct=ft−Ccor1,Cr=fr−Ccor1,(12)Ct=ft−Ccor2,Cl=fl−Ccor2,(13)Cl=fl−Ccor3,Cb=fb−Ccor3,(14)Cr=fr−Ccor4,Cb=fb−Ccor4.It indicates that the fluxes in four oblique directions need to be calculated first.(2) Fluid Fluxes in the Oblique Directions. Taking the case of fluid flow from the central cell to the upper right cell (Figure 5(b)), the oblique flux Ccor1 was derived. The fluid could be in nine locations after the oblique movement, as shown in Figure 6. “+” indicates the grid line, the dotted polygon is the position before the fluid moves, solid polygon indicates the current position, and the shaded portion shows the fluid flowing into the upper right cell. The shaded area is the required flux Ccor1.Figure 6
Nine cases in which unit fluid may appear in calculatingCcor1 for case III.
(a)(b)(c)(d)(e)(f)(g)(h)(i)In Figure7, lm is the height of the polygon fluid beyond the intersected horizontal cell line in case III, lb is the length of the horizontal cell line intersected with the fluid polygon, and θ is the angle of the fluid polygon bevel inclination.Ccor1 is calculated as follows:ifft>0andfr>0then(excluding Figure6 (1), (3), (6), (8), (9))lm=utδtlb=lmcotθ+stδxiflb>urδt (two cases as Figure 6 (2), (5))Ccor1=fr−urδtδy−utδtelse (two cases as Figure6 (4), (7))Ccor1=ftFigure 7
Parameters’ definition for calculatingCcor1.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)), Ccor2, Ccor3, and Ccor4, can be written as follows.ForCcor2,ifft>0andfl>0thenlm=utδtlb=lmcotθ+stδxiflb+ulδt>δxifulδt≥δxCcor2=ftelseCcor2=fl−0.5δy−utδtsbδx+lb+2ulδt−2δxForCcor3,iffb>0andfl>0thenlm=δy−ubδtlb=lmcotθ+stδxifulδt<δxiflb>δx−ulδtCcor3=fb−ubδtδx−ulδtelseCcor3=flelseCcor3=fbForCcor4,iffb>0andfr>0thenlm=δy−ubδtlb=lbcotθ+stδxifurδt<δxiflb≥urδtCcor4=ubδturδtelseCcor4=fr−1/2lslb+stδxelseCcor4=fbSimilarly, the flux expressions of case I, case II, and case IV in four oblique directions can be solved.(3) Update the VOF Function. After the fluid fluxes from the central cell (i, j) to eight adjacent cells during a time δt are calculated, the VOF function can be updated. Given fractions of the fluid volume in the central cell and the surrounding eight adjacent cells at time t are Ft (i, j), Ft–(i, j + 1), Ft (i, j − 1), Ft (i − 1, j), Ft (i + 1, j), Ft (i + 1, j + 1), Ft (i − 1, j + 1), Ft (i − 1, j − 1), and Ft (i + 1, j − 1), respectively, at time t + δt, they are(15)Ft+δti,j=Fti,j−Ct+Cb+Cl+Cr+Ccor1+Ccor2+Ccor3+Ccor4δxδy,(16)Ft+δti+1,j=Fti+1,j+CrδxδyFt+δti,j+1=Fti,j+1+CtδxδyFt+δti−1,j=Fti−1,j+ClδxδyFt+δti,j−1=Fti,j−1+CbδxδyFt+δti+1,j+1=Fti+1,j+1+Ccor1δxδyFt+δti−1,j+1=Fti−1,j+1+Ccor2δxδyFt+δti−1,j−1=Fti−1,j−1+Ccor3δxδyFt+δti+1,j−1=Fti+1,j−1+Ccor4δxδy.Repeat the calculation for each cell; then, the update for the VOF function in the whole field can be completed.


## 3.1. Construction of the Fluid Moving Interface

The 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)
The cell is full of fluidA, F = 1(2)
The cell is an interface unit, 0 <F < 1(3)
The cell is filled with the fluidB, F = 0Case 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(5)nijx=Fi+1,j+1+2Fi+1,j+Fi+1,j−1−Fi−1,j+1−2Fi−1,j−Fi−1,j−1δx.nijy=Fi+1,j+1+2Fi,j+1+Fi−1,j+1−Fi+1,j−1−2Fi,j−1−Fi−1,j−1δy.Then, the angleβ that the interface makes with the x-axis is determined according to the normal vector:(6)tanβ=−nxny.The angleβ is normalized to ω:(7)tanω=tanβδxδy,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.Figure 1
Schematic of the fluid cell.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. The four types of interfaces are shown in Figure 2. The discriminant method is as follows.ifω < π/4ifF ≤ 0.5 tanωcase Ielse ifF ≤ 1.0 − 0.5 tanωcase IIelsecase IVelseifF ≤ 0.5 cotωcase Ielse ifF ≤ 1.0 − 0.5 cotωcase IIIelsecase IVFigure 2
Four types of interface.
(a)(b)(c)(d)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 (st, sb, sl, sr) can be calculated according to the angle ω and the VOF function. The 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:(8)st=F−0.5cotωsr=1sb=F+0.5tanωsl=0.After obtaining the side fractions of each cell, the straight-line segment that represents the fluid interface in the cell can be plotted.


## 3.2. Advancement of the Moving Interface (Update VOF Function)

According to the current velocity, the fluxes through sides of the cell to adjacent cells in a time step are calculated. The VOF function of each cell is updated by accumulating the inflow and the outflow. The moving interface is determined and advances according to the updated VOF function.The 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.


### 3.2.1. 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 ft, fb, fl, and fr, respectively, as shown in Figure 3.Figure 3
Schematic of fluid transport calculation for split Youngs algorithm.Taking case III (in Figure2) as an example, according to the geometric law, the fluxes can be calculated as follows.(1)
Top face:ifut > 0ft=utδtstδx+0.5utδxcotβ(2)
Bottom face:ifub < 0fb=absubδtsbδx−0.5absubδtcotβ(3)
Left side:iful < 0if abs (ul)δt ≤ (1 − sb)δxfl=0elseifabsulδt≤1−stδxfl=0.5absulδt−1−sbδx2tanβelsefl=absulδtδy−1−cδxδy(4)
Right side:ifur > 0ifurδt≤stδxfr=urδtδyelseifurδt≤sbδxfr=urδtδy−0.5urδt−stδx2tanβelsefr=Fδxδywhereut, ub, ur, and ul represent the interface velocity in the top, bottom, left, and right, respectively. t is the time step. The fluxes of cases I, II, and IV in four directions can be calculated in the same way.(2) Update VOF Function. The 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 Ft 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 ft, fb, fl, and fr, respectively. The volume fraction in the cell at time t + δt is(9)Ft+δti,j=Fti,j−ft+fb+fl+frδxδy.The fractions of fluid volume in the four adjacent cells at timet + δt are, respectively, corrected as follows:(10)Ft+δti,j+1=Fti,j+1+ftδxδyFt+δti,j−1=Fti,j−1+fbδxδyFt+δti+1,j=Fti+1,j+frδxδyFt+δti−1,j=Fti−1,j+flδxδy.Conducting the same calculation for each cell can update the fractional volume of the flow field.The 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. The fluid flow towards the oblique cell is ignored, which may lead to some deviations from the real.


### 3.2.2. Unsplit Youngs Algorithm

In order to improve tracking accuracy, an unsplit operator algorithm was introduced to advance the interface. The 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 Figure4. The fluxes in eight directions are denoted as Ct, Cb, Cl, Cr, Ccor1, Ccor2, and Ccor3, Ccor4, 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.Figure 4
Schematic of fluid transport calculation for unsplit Youngs algorithm.Hu et al. [16] assumed the unit fluid form is a triangle in the cell containing the interface. Five cases were listed when Ut > 0 and Ur > 0. The relationship between the above eight fluxes was described and ft, fb, fl, and fr 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. The 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. The implementation is shown as follows.(1) Fluid Fluxes in the Orthogonal Directions. The fluxes from the central cell to four adjacent cells were supposed as Ct, Cb, Cl, and Cr. Figure 5 shows the pattern of fluid transport from the central cell (Figure 5(a)) to the upper right cell (Figure 5(b)), the upper left cell (Figure 5(c)), the lower left cell (Figure 5(d)), and the lower right cell (Figure 5(e)) from tn to tn+1.Figure 5
Schematic of fluid transport from the central celli, j to the four oblique cells: (a) the position at time tn, (b) the position at time tn for case 1, (c) the position at time tn+1 for case 2, (d) the position at time tn+1 for case 3, and (e) the position at time tn+1 for case 4.
(a)(b)(c)(d)(e)According to Figures5(b)–5(e), the equations can be expressed as follows:(11)Ct=ft−Ccor1,Cr=fr−Ccor1,(12)Ct=ft−Ccor2,Cl=fl−Ccor2,(13)Cl=fl−Ccor3,Cb=fb−Ccor3,(14)Cr=fr−Ccor4,Cb=fb−Ccor4.It indicates that the fluxes in four oblique directions need to be calculated first.(2) Fluid Fluxes in the Oblique Directions. Taking the case of fluid flow from the central cell to the upper right cell (Figure 5(b)), the oblique flux Ccor1 was derived. The fluid could be in nine locations after the oblique movement, as shown in Figure 6. “+” indicates the grid line, the dotted polygon is the position before the fluid moves, solid polygon indicates the current position, and the shaded portion shows the fluid flowing into the upper right cell. The shaded area is the required flux Ccor1.Figure 6
Nine cases in which unit fluid may appear in calculatingCcor1 for case III.
(a)(b)(c)(d)(e)(f)(g)(h)(i)In Figure7, lm is the height of the polygon fluid beyond the intersected horizontal cell line in case III, lb is the length of the horizontal cell line intersected with the fluid polygon, and θ is the angle of the fluid polygon bevel inclination.Ccor1 is calculated as follows:ifft>0andfr>0then(excluding Figure6 (1), (3), (6), (8), (9))lm=utδtlb=lmcotθ+stδxiflb>urδt (two cases as Figure 6 (2), (5))Ccor1=fr−urδtδy−utδtelse (two cases as Figure6 (4), (7))Ccor1=ftFigure 7
Parameters’ definition for calculatingCcor1.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)), Ccor2, Ccor3, and Ccor4, can be written as follows.ForCcor2,ifft>0andfl>0thenlm=utδtlb=lmcotθ+stδxiflb+ulδt>δxifulδt≥δxCcor2=ftelseCcor2=fl−0.5δy−utδtsbδx+lb+2ulδt−2δxForCcor3,iffb>0andfl>0thenlm=δy−ubδtlb=lmcotθ+stδxifulδt<δxiflb>δx−ulδtCcor3=fb−ubδtδx−ulδtelseCcor3=flelseCcor3=fbForCcor4,iffb>0andfr>0thenlm=δy−ubδtlb=lbcotθ+stδxifurδt<δxiflb≥urδtCcor4=ubδturδtelseCcor4=fr−1/2lslb+stδxelseCcor4=fbSimilarly, the flux expressions of case I, case II, and case IV in four oblique directions can be solved.(3) Update the VOF Function. After the fluid fluxes from the central cell (i, j) to eight adjacent cells during a time δt are calculated, the VOF function can be updated. Given fractions of the fluid volume in the central cell and the surrounding eight adjacent cells at time t are Ft (i, j), Ft–(i, j + 1), Ft (i, j − 1), Ft (i − 1, j), Ft (i + 1, j), Ft (i + 1, j + 1), Ft (i − 1, j + 1), Ft (i − 1, j − 1), and Ft (i + 1, j − 1), respectively, at time t + δt, they are(15)Ft+δti,j=Fti,j−Ct+Cb+Cl+Cr+Ccor1+Ccor2+Ccor3+Ccor4δxδy,(16)Ft+δti+1,j=Fti+1,j+CrδxδyFt+δti,j+1=Fti,j+1+CtδxδyFt+δti−1,j=Fti−1,j+ClδxδyFt+δti,j−1=Fti,j−1+CbδxδyFt+δti+1,j+1=Fti+1,j+1+Ccor1δxδyFt+δti−1,j+1=Fti−1,j+1+Ccor2δxδyFt+δti−1,j−1=Fti−1,j−1+Ccor3δxδyFt+δti+1,j−1=Fti+1,j−1+Ccor4δxδy.Repeat the calculation for each cell; then, the update for the VOF function in the whole field can be completed.


## 3.2.1. 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 ft, fb, fl, and fr, respectively, as shown in Figure 3.Figure 3
Schematic of fluid transport calculation for split Youngs algorithm.Taking case III (in Figure2) as an example, according to the geometric law, the fluxes can be calculated as follows.(1)
Top face:ifut > 0ft=utδtstδx+0.5utδxcotβ(2)
Bottom face:ifub < 0fb=absubδtsbδx−0.5absubδtcotβ(3)
Left side:iful < 0if abs (ul)δt ≤ (1 − sb)δxfl=0elseifabsulδt≤1−stδxfl=0.5absulδt−1−sbδx2tanβelsefl=absulδtδy−1−cδxδy(4)
Right side:ifur > 0ifurδt≤stδxfr=urδtδyelseifurδt≤sbδxfr=urδtδy−0.5urδt−stδx2tanβelsefr=Fδxδywhereut, ub, ur, and ul represent the interface velocity in the top, bottom, left, and right, respectively. t is the time step. The fluxes of cases I, II, and IV in four directions can be calculated in the same way.(2) Update VOF Function. The 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 Ft 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 ft, fb, fl, and fr, respectively. The volume fraction in the cell at time t + δt is(9)Ft+δti,j=Fti,j−ft+fb+fl+frδxδy.The fractions of fluid volume in the four adjacent cells at timet + δt are, respectively, corrected as follows:(10)Ft+δti,j+1=Fti,j+1+ftδxδyFt+δti,j−1=Fti,j−1+fbδxδyFt+δti+1,j=Fti+1,j+frδxδyFt+δti−1,j=Fti−1,j+flδxδy.Conducting the same calculation for each cell can update the fractional volume of the flow field.The 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. The fluid flow towards the oblique cell is ignored, which may lead to some deviations from the real.


## 3.2.2. Unsplit Youngs Algorithm

In order to improve tracking accuracy, an unsplit operator algorithm was introduced to advance the interface. The 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 Figure4. The fluxes in eight directions are denoted as Ct, Cb, Cl, Cr, Ccor1, Ccor2, and Ccor3, Ccor4, 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.Figure 4
Schematic of fluid transport calculation for unsplit Youngs algorithm.Hu et al. [16] assumed the unit fluid form is a triangle in the cell containing the interface. Five cases were listed when Ut > 0 and Ur > 0. The relationship between the above eight fluxes was described and ft, fb, fl, and fr 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. The 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. The implementation is shown as follows.(1) Fluid Fluxes in the Orthogonal Directions. The fluxes from the central cell to four adjacent cells were supposed as Ct, Cb, Cl, and Cr. Figure 5 shows the pattern of fluid transport from the central cell (Figure 5(a)) to the upper right cell (Figure 5(b)), the upper left cell (Figure 5(c)), the lower left cell (Figure 5(d)), and the lower right cell (Figure 5(e)) from tn to tn+1.Figure 5
Schematic of fluid transport from the central celli, j to the four oblique cells: (a) the position at time tn, (b) the position at time tn for case 1, (c) the position at time tn+1 for case 2, (d) the position at time tn+1 for case 3, and (e) the position at time tn+1 for case 4.
(a)(b)(c)(d)(e)According to Figures5(b)–5(e), the equations can be expressed as follows:(11)Ct=ft−Ccor1,Cr=fr−Ccor1,(12)Ct=ft−Ccor2,Cl=fl−Ccor2,(13)Cl=fl−Ccor3,Cb=fb−Ccor3,(14)Cr=fr−Ccor4,Cb=fb−Ccor4.It indicates that the fluxes in four oblique directions need to be calculated first.(2) Fluid Fluxes in the Oblique Directions. Taking the case of fluid flow from the central cell to the upper right cell (Figure 5(b)), the oblique flux Ccor1 was derived. The fluid could be in nine locations after the oblique movement, as shown in Figure 6. “+” indicates the grid line, the dotted polygon is the position before the fluid moves, solid polygon indicates the current position, and the shaded portion shows the fluid flowing into the upper right cell. The shaded area is the required flux Ccor1.Figure 6
Nine cases in which unit fluid may appear in calculatingCcor1 for case III.
(a)(b)(c)(d)(e)(f)(g)(h)(i)In Figure7, lm is the height of the polygon fluid beyond the intersected horizontal cell line in case III, lb is the length of the horizontal cell line intersected with the fluid polygon, and θ is the angle of the fluid polygon bevel inclination.Ccor1 is calculated as follows:ifft>0andfr>0then(excluding Figure6 (1), (3), (6), (8), (9))lm=utδtlb=lmcotθ+stδxiflb>urδt (two cases as Figure 6 (2), (5))Ccor1=fr−urδtδy−utδtelse (two cases as Figure6 (4), (7))Ccor1=ftFigure 7
Parameters’ definition for calculatingCcor1.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)), Ccor2, Ccor3, and Ccor4, can be written as follows.ForCcor2,ifft>0andfl>0thenlm=utδtlb=lmcotθ+stδxiflb+ulδt>δxifulδt≥δxCcor2=ftelseCcor2=fl−0.5δy−utδtsbδx+lb+2ulδt−2δxForCcor3,iffb>0andfl>0thenlm=δy−ubδtlb=lmcotθ+stδxifulδt<δxiflb>δx−ulδtCcor3=fb−ubδtδx−ulδtelseCcor3=flelseCcor3=fbForCcor4,iffb>0andfr>0thenlm=δy−ubδtlb=lbcotθ+stδxifurδt<δxiflb≥urδtCcor4=ubδturδtelseCcor4=fr−1/2lslb+stδxelseCcor4=fbSimilarly, the flux expressions of case I, case II, and case IV in four oblique directions can be solved.(3) Update the VOF Function. After the fluid fluxes from the central cell (i, j) to eight adjacent cells during a time δt are calculated, the VOF function can be updated. Given fractions of the fluid volume in the central cell and the surrounding eight adjacent cells at time t are Ft (i, j), Ft–(i, j + 1), Ft (i, j − 1), Ft (i − 1, j), Ft (i + 1, j), Ft (i + 1, j + 1), Ft (i − 1, j + 1), Ft (i − 1, j − 1), and Ft (i + 1, j − 1), respectively, at time t + δt, they are(15)Ft+δti,j=Fti,j−Ct+Cb+Cl+Cr+Ccor1+Ccor2+Ccor3+Ccor4δxδy,(16)Ft+δti+1,j=Fti+1,j+CrδxδyFt+δti,j+1=Fti,j+1+CtδxδyFt+δti−1,j=Fti−1,j+ClδxδyFt+δti,j−1=Fti,j−1+CbδxδyFt+δti+1,j+1=Fti+1,j+1+Ccor1δxδyFt+δti−1,j+1=Fti−1,j+1+Ccor2δxδyFt+δti−1,j−1=Fti−1,j−1+Ccor3δxδyFt+δti+1,j−1=Fti+1,j−1+Ccor4δxδy.Repeat the calculation for each cell; then, the update for the VOF function in the whole field can be completed.


## 4. Case Studies

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


### 4.1. Calculation Model

The 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 Figure8. The 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 R0. According to Mitani and Hamada [23], the grout density changes with time as below:(17)ρt=ρ01+αt,where ρ0 is the initial grout density with a value of 1100 kg/m3, α is the expansion rate, t is the time, and ρt is the grout density at time t. The grout viscosity remains unchanged at 1.0 pa·s before gelation.Figure 8
Fracture model.Considering the fracture model is geometrical symmetry, the solving could be simplified by only calculating a quarter of the simulated region, as shown in Figure9. Therefore, there are two symmetry boundaries and the other boundary conditions are pressure outlet.Figure 9
Calculation region.The grouting quantity is represented by the initial diffusion radius of expansible grout in a fracture. The grouts with different components have different expansion rates. The input parameters are shown in Table1. Two types of expansible grout with different initial radius were used to validate the unsplit Youngs algorithm.Table 1
Input material parameters.
MaterialsDensity (kg/m3)Viscosity (Pa·s)Expansion rate (s−1)Initial diffusion radius (m)Air1.2051.8 × 10−5——Grout1100 (initial)1.00.10.250.10.550.20.250.20.55The analytical solution of the diffusion distance of the expansible grout in a fracture is expressed as follows [24]:(18)Rt=R01+αt,where R0 is the initial diffusion radius and Rt is the diffusion radius at time t.


### 4.2. Numerical Calculation

The diffusion behavior of the expansible grout can be described by one continuity equation and three momentum conservation equations. The continuity equation is given as(19)∂ρtF∂t+∇⋅ρtFU=0,where F is the volume fraction of the expansible grout and U is the velocity vector.The two phases of grout and air share one single set of momentum equations, which are written as follows:(20)∂ρU∂t+∇⋅ρUU=−∇p+∇⋅μ∇U+∇UT+ρg,where g is the gravity acceleration, ρ and μ are the volume fraction weighted average of the density and viscosity, respectively, and p is the pressure.The 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). The 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 interfaceReturn to step (2), update the time, and the numerical process proceeds to the next iteration. The numerical calculation was programmed by FORTRAN. The procedure is shown in Figure10.Figure 10
Numerical procedure.


### 4.3. Comparison and Discussion

Since different grouting quantities and expansion rates lead to different diffusion rates, in Cases 1–4, the grout takes different times to reach the boundary. Considering the limitation of the simulation domain, the calculation time is set as 240 s, 100 s, 55 s, and 24 s for Cases 1–4.Figures11–14 are cloud charts of the diffusion process of the expansible grout in Case 1 (t = 80 s, 160 s, 240 s), Case 2 (t = 20 s, 60 s, 100 s), Case 3 (t = 15 s, 35 s, 55 s), and Case 4 (t = 8 s, 16 s, 24 s), respectively. It shows that the grout diffusion rate is quite different under different conditions of the initial filling range and the grout expansion rates. 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. The grout diffusion distance calculated by the unsplit Youngs algorithm is in good agreement with the analytical solution. The 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. The diffusion ranges calculated by the split algorithm at different times lag behind the theoretical results. The grout contour calculated by the split algorithm always lies inside the analytical solution, and the deviation increases gradually with time.Figure 11
Numerical results of grout spread, shape, and range of Case 1. (a) Numerical results of the split algorithm att = 80 s, 160 s, and 240 s. (b) Numerical results of unsplit algorithm at t = 80 s, 160 s, and 240 s.
(a)(b)Figure 12
Numerical results of grout spread, shape, and range of Case 2. (a) Numerical results of split algorithm att = 20 s, 60 s, and 100 s (b) Numerical results of unsplit algorithm at t = 20 s, 60 s, and 100 s.
(a)(b)Figure 13
Numerical results of grout spread, shape, and range of Case 3. (a) Numerical results of the split algorithm att = 15 s, 35 s, and 55 s. (b) Numerical results of unsplit algorithm at t = 15 s, 35 s, and 55 s.
(a)(b)Figure 14
Numerical results of grout spread, shape, and range of Case 4. (a) Numerical results of the split algorithm att = 8 s, 16 s, and 24 s. (b) Numerical results of unsplit algorithm at t = 8 s, 16 s, and 24 s.
(a)(b)The average diffusion radius at different times was calculated for Cases 1–4, shown in Tables2–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. The maximum absolute errors are 10.71 cm, 7.99 cm, 4.05 cm, and 4.13 cm, respectively.Table 2
Calculation results of grout spread radius at different times in case 1.
Time (s)Analytical (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)200.24490.22330.02178.850.23970.00522.14400.33170.29740.034310.340.32510.00651.97600.40000.35470.045311.330.39230.00771.91800.45830.40390.054411.860.45030.00801.751000.50990.44790.062012.160.50160.00831.631200.55680.48740.069412.460.54770.00911.63140060000.52400.076012.670.59090.00911.511600.64030.55830.082012.810.63120.00911.421800.67820.59060.087612.920.66880.00951.402000.71410.62080.093413.080.70450.00961.352200.74830.64990.098413.150.73910.00931.242400.78100.67400.107113.710.77180.00921.18Table 3
Calculation results of grout spread radius at different times in case 2.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)100.25500.23730.01766.900.24890.00602.37200.34640.31740.02908.370.33910.00732.10300.41830.38030.03819.100.40990.00842.01400.47960.43340.04629.640.47040.00911.91500.53390.48080.05319.950.52430.00961.80600.58310.52390.059210.150.57310.01001.72700.62850.56340.065110.350.61820.01031.64800.67080.60050.070310.490.66030.01051.57900.71060.63510.075510.630.69990.01071.511000.74830.66840.079910.680.73730.01101.47Table 4
Calculation results of grout spread radius at different times in case 3.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)50.30000.29270.00732.440.29640.00361.21100.37420.36120.01303.470.37070.00350.93150.43590.41870.01723.940.43180.00410.94200.48990.46910.02084.240.48600.00390.79250.53850.51400.02454.550.53450.00400.74300.58310.55530.02784.760.57980.00330.57350.62450.59400.03054.890.62100.00350.56400.66330.63020.03325.000.66020.00320.48450.70000.66400.03605.150.69660.00340.48500.73480.69640.03855.230.73140.00350.48550.76810.72760.04055.270.76460.00350.46Table 5
Calculation results of grout spread radius at different times in Case 4.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Numerical results (m)Absolute error (%)20.28980.28310.00672.330.28740.00250.8540.35780.34600.01173.280.35560.00210.6060.41470.39890.01593.830.41210.00260.6380.46480.44500.01974.240.46340.00130.28100.50990.48660.02334.570.50880.00110.21120.55140.52520.02614.740.55030.00110.20140.58990.56080.02914.940.58730.00260.44160.62610.59440.03175.070.62360.00250.40180.66030.62580.03455.220.65840.00190.29200.69280.65580.03705.340.69100.00180.26220.72390.68490.03895.380.72250.00140.20240.75370.71230.04135.480.75240.00120.17Figure 15
Comparison of grout spread radius at different times: (a)R0 = 0.1 m, α = 0.25; (b) R0 = 0.1 m, α = 0.55; (c) R0 = 0.2 m, α = 0.25; (d) R0 = 0.2 m, α = 0.55.
(a)(b)(c)(d)Figure16 shows the relative error between the numerical method and analytical solution with respect to the diffusion radius under the four cases. The 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 diffusion radius, while the relative error of the unsplit algorithm declines slightly.Figure 16
Relative error with respect to different grout spread radii: (a)R0 = 0.1 m, α = 0.25; (b) R0 = 0.1 m, α = 0.55; (c) R0 = 0.2 m, α = 0.25; (d) R0 = 0.2 m, α = 0.55.
(a)(b)(c)(d)The 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 direction. This reduces the calculation error effectively and improves the accuracy of tracking the moving interface significantly.It should be noted that the present algorithm is still complex. There are four types of interface when it locates in the cell. For each type, fluid transport could appear in eight directions. Each fluid flux from the central cell to four oblique direction cells could have nine cases. In addition, the proposed method can only be used to solve two-dimensional interface tracking. The three-dimensional Youngs method has to take into account the fluid flux in thez-direction, which needs to reprogram algorithm. This 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.


## 4.1. Calculation Model

The 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 Figure8. The 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 R0. According to Mitani and Hamada [23], the grout density changes with time as below:(17)ρt=ρ01+αt,where ρ0 is the initial grout density with a value of 1100 kg/m3, α is the expansion rate, t is the time, and ρt is the grout density at time t. The grout viscosity remains unchanged at 1.0 pa·s before gelation.Figure 8
Fracture model.Considering the fracture model is geometrical symmetry, the solving could be simplified by only calculating a quarter of the simulated region, as shown in Figure9. Therefore, there are two symmetry boundaries and the other boundary conditions are pressure outlet.Figure 9
Calculation region.The grouting quantity is represented by the initial diffusion radius of expansible grout in a fracture. The grouts with different components have different expansion rates. The input parameters are shown in Table1. Two types of expansible grout with different initial radius were used to validate the unsplit Youngs algorithm.Table 1
Input material parameters.
MaterialsDensity (kg/m3)Viscosity (Pa·s)Expansion rate (s−1)Initial diffusion radius (m)Air1.2051.8 × 10−5——Grout1100 (initial)1.00.10.250.10.550.20.250.20.55The analytical solution of the diffusion distance of the expansible grout in a fracture is expressed as follows [24]:(18)Rt=R01+αt,where R0 is the initial diffusion radius and Rt is the diffusion radius at time t.


## 4.2. Numerical Calculation

The diffusion behavior of the expansible grout can be described by one continuity equation and three momentum conservation equations. The continuity equation is given as(19)∂ρtF∂t+∇⋅ρtFU=0,where F is the volume fraction of the expansible grout and U is the velocity vector.The two phases of grout and air share one single set of momentum equations, which are written as follows:(20)∂ρU∂t+∇⋅ρUU=−∇p+∇⋅μ∇U+∇UT+ρg,where g is the gravity acceleration, ρ and μ are the volume fraction weighted average of the density and viscosity, respectively, and p is the pressure.The 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). The 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 interfaceReturn to step (2), update the time, and the numerical process proceeds to the next iteration. The numerical calculation was programmed by FORTRAN. The procedure is shown in Figure10.Figure 10
Numerical procedure.


## 4.3. Comparison and Discussion

Since different grouting quantities and expansion rates lead to different diffusion rates, in Cases 1–4, the grout takes different times to reach the boundary. Considering the limitation of the simulation domain, the calculation time is set as 240 s, 100 s, 55 s, and 24 s for Cases 1–4.Figures11–14 are cloud charts of the diffusion process of the expansible grout in Case 1 (t = 80 s, 160 s, 240 s), Case 2 (t = 20 s, 60 s, 100 s), Case 3 (t = 15 s, 35 s, 55 s), and Case 4 (t = 8 s, 16 s, 24 s), respectively. It shows that the grout diffusion rate is quite different under different conditions of the initial filling range and the grout expansion rates. 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. The grout diffusion distance calculated by the unsplit Youngs algorithm is in good agreement with the analytical solution. The 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. The diffusion ranges calculated by the split algorithm at different times lag behind the theoretical results. The grout contour calculated by the split algorithm always lies inside the analytical solution, and the deviation increases gradually with time.Figure 11
Numerical results of grout spread, shape, and range of Case 1. (a) Numerical results of the split algorithm att = 80 s, 160 s, and 240 s. (b) Numerical results of unsplit algorithm at t = 80 s, 160 s, and 240 s.
(a)(b)Figure 12
Numerical results of grout spread, shape, and range of Case 2. (a) Numerical results of split algorithm att = 20 s, 60 s, and 100 s (b) Numerical results of unsplit algorithm at t = 20 s, 60 s, and 100 s.
(a)(b)Figure 13
Numerical results of grout spread, shape, and range of Case 3. (a) Numerical results of the split algorithm att = 15 s, 35 s, and 55 s. (b) Numerical results of unsplit algorithm at t = 15 s, 35 s, and 55 s.
(a)(b)Figure 14
Numerical results of grout spread, shape, and range of Case 4. (a) Numerical results of the split algorithm att = 8 s, 16 s, and 24 s. (b) Numerical results of unsplit algorithm at t = 8 s, 16 s, and 24 s.
(a)(b)The average diffusion radius at different times was calculated for Cases 1–4, shown in Tables2–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. The maximum absolute errors are 10.71 cm, 7.99 cm, 4.05 cm, and 4.13 cm, respectively.Table 2
Calculation results of grout spread radius at different times in case 1.
Time (s)Analytical (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)200.24490.22330.02178.850.23970.00522.14400.33170.29740.034310.340.32510.00651.97600.40000.35470.045311.330.39230.00771.91800.45830.40390.054411.860.45030.00801.751000.50990.44790.062012.160.50160.00831.631200.55680.48740.069412.460.54770.00911.63140060000.52400.076012.670.59090.00911.511600.64030.55830.082012.810.63120.00911.421800.67820.59060.087612.920.66880.00951.402000.71410.62080.093413.080.70450.00961.352200.74830.64990.098413.150.73910.00931.242400.78100.67400.107113.710.77180.00921.18Table 3
Calculation results of grout spread radius at different times in case 2.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)100.25500.23730.01766.900.24890.00602.37200.34640.31740.02908.370.33910.00732.10300.41830.38030.03819.100.40990.00842.01400.47960.43340.04629.640.47040.00911.91500.53390.48080.05319.950.52430.00961.80600.58310.52390.059210.150.57310.01001.72700.62850.56340.065110.350.61820.01031.64800.67080.60050.070310.490.66030.01051.57900.71060.63510.075510.630.69990.01071.511000.74830.66840.079910.680.73730.01101.47Table 4
Calculation results of grout spread radius at different times in case 3.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Absolute error (m)Relative error (%)50.30000.29270.00732.440.29640.00361.21100.37420.36120.01303.470.37070.00350.93150.43590.41870.01723.940.43180.00410.94200.48990.46910.02084.240.48600.00390.79250.53850.51400.02454.550.53450.00400.74300.58310.55530.02784.760.57980.00330.57350.62450.59400.03054.890.62100.00350.56400.66330.63020.03325.000.66020.00320.48450.70000.66400.03605.150.69660.00340.48500.73480.69640.03855.230.73140.00350.48550.76810.72760.04055.270.76460.00350.46Table 5
Calculation results of grout spread radius at different times in Case 4.
Time (s)Analytical results (m)Split algorithmUnsplit algorithmNumerical results (m)Absolute error (m)Relative error (%)Numerical results (m)Numerical results (m)Absolute error (%)20.28980.28310.00672.330.28740.00250.8540.35780.34600.01173.280.35560.00210.6060.41470.39890.01593.830.41210.00260.6380.46480.44500.01974.240.46340.00130.28100.50990.48660.02334.570.50880.00110.21120.55140.52520.02614.740.55030.00110.20140.58990.56080.02914.940.58730.00260.44160.62610.59440.03175.070.62360.00250.40180.66030.62580.03455.220.65840.00190.29200.69280.65580.03705.340.69100.00180.26220.72390.68490.03895.380.72250.00140.20240.75370.71230.04135.480.75240.00120.17Figure 15
Comparison of grout spread radius at different times: (a)R0 = 0.1 m, α = 0.25; (b) R0 = 0.1 m, α = 0.55; (c) R0 = 0.2 m, α = 0.25; (d) R0 = 0.2 m, α = 0.55.
(a)(b)(c)(d)Figure16 shows the relative error between the numerical method and analytical solution with respect to the diffusion radius under the four cases. The 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 diffusion radius, while the relative error of the unsplit algorithm declines slightly.Figure 16
Relative error with respect to different grout spread radii: (a)R0 = 0.1 m, α = 0.25; (b) R0 = 0.1 m, α = 0.55; (c) R0 = 0.2 m, α = 0.25; (d) R0 = 0.2 m, α = 0.55.
(a)(b)(c)(d)The 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 direction. This reduces the calculation error effectively and improves the accuracy of tracking the moving interface significantly.It should be noted that the present algorithm is still complex. There are four types of interface when it locates in the cell. For each type, fluid transport could appear in eight directions. Each fluid flux from the central cell to four oblique direction cells could have nine cases. In addition, the proposed method can only be used to solve two-dimensional interface tracking. The three-dimensional Youngs method has to take into account the fluid flux in thez-direction, which needs to reprogram algorithm. This 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.


## 5. 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. The diffusion process of expansible grout in a parallel-plate fracture is studied by the proposed algorithm and split Youngs algorithm. The 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.This 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.



---

*Source: 2014959-2020-10-29.xml*