Modelling Time–Dependent Flow through Railway Ballast

1. Introduction

The understanding of water flow behavior through porous media is crucial in many applications to minimize damage caused by seepage [1,2,3,4,5]. In this context, the present study focuses on water flow through railway ballast, which is commonly modeled as a porous media. Of particular importance is the scouring of railway ballast due to flood waters, which can occur when the upstream height exceeds the railway ballast’s height [6,7]. Heterogeneous ballast gradations can also cause erosion for the ballast, even if the upstream height is less than the railway ballast’s height [8]. Scouring of railway ballast can lead to catastrophic outcomes, as discussed in [9], particularly if the ballast becomes too weak to support the weight of a train.

However, analyzing the problem of porous media flow is much more challenging when the free surface needs to be calculated. The problem becomes even more complicated when considering the time-dependent motion of the problem. Nevertheless, solving the time-dependent problem, which involves determining the free surface position of water flow, can be used to provide the steady-state solution (which is commonly calculated iteratively) [10]. In our present study, we also use the method developed to solve for the steady-state unknown free surface.

Nonetheless, the primary challenge is to solve the unsteady problem, where the actual time is considered. This problem represents realistic flow through porous media and the time-dependent change in the position of the free surface. We must also consider the specific yield $S_y$ when investigating the motion of the free surface for unsteady flow, whereas this is not necessary for determining the steady-state solution [10]. The specific yield $S_y$ models the drainage capacity of the porous material [11]. In the present work, we particularly focus on validating our calculations by comparing them with results in the literature, particularly those in [12,13] for three cases in which they solved for the steady-state flow through a rectangular dam. For the unsteady calculations, we use the experimental results calculated in [14] and the theoretical results of [15].

The aim of this work is to develop numerical methods to determine the unsteady flow through railway ballast with a moving free boundary. We consider realistic ballast models that represent the effect of ballast fouling, which is regarded as critical to ballast flooding. We also thoroughly and carefully validate our present study by comparing it with previous studies [12,13,14,15] and our previous study [9], which solved the steady free boundary problem. We use the Finite Element Method (FEM) for all our calculations.

2. Methods and Numerical Modelling of Flow through a Ballast

The region under consideration is depicted in Figure 1 and can be defined by the domain $D:\left[\right(x,z\left)\right|0\le x\le a,0\le z\le b]$. The saturated part of the region is denoted by $\Omega \left(t\right)\subset D$, which is a function of time t, where $\Omega \left(t\right):\left[\right(x,z\left)\right|0\le x\le a,0\le z\le \phi (x,t)\le b]$. Here, $\phi (x,t)$ represents the water level inside the rectangular ballast. The free surface within the porous media is denoted by ${\Gamma }_5\left(t\right)$, while ${\Gamma }_1\left(t\right)$ and ${\Gamma }_3\left(t\right)$ represent the upstream and downstream regions, respectively, with a fluid boundary. The seepage face is denoted by ${\Gamma }_4\left(t\right)$. The water heights are represented by $H_1\left(t\right)$ and $H_2\left(t\right)$, while the impermeable boundary is represented by ${\Gamma }_2$, where there is no flux.

The governing equation is derived from the divergence of the seepage velocity and is given by

$$\nabla ·K_s\nabla \varphi =0,$$

(1)

where $K_s$ and $\varphi$ are the hydraulic conductivity and the hydraulic head, respectively [16]. We use the following boundary conditions [12]

$$\left\{\begin{array}{cc}\varphi (x,z,t)=H_1\left(t\right)\hfill & \mathrm{on}{\Gamma }_1\left(t\right),\forall t\in [0,T],\hfill \\ \varphi (x,z,t)=H_2\left(t\right)\hfill & \mathrm{on}{\Gamma }_3\left(t\right),\forall t\in [0,T],\hfill \\ \varphi (x,z,t)=z\hfill & \mathrm{on}{\Gamma }_4\left(t\right)\cup {\Gamma }_5\left(t\right),\forall t\in [0,T],\hfill \\ \frac{\partial \varphi }{\partial v}=0\hfill & \mathrm{on}{\Gamma }_2,\forall t\in [0,T].\hfill \end{array}\right.$$

(2)

The final condition is that for the movement of the free surface

$$\frac{\partial \phi }{\partial t}=[(\frac{K_s}{S_y}\nabla \varphi )·\mathbf{n}],$$

(3)

where $S_y$ is the specific yield. We will discuss this term shortly. We also need to specify the initial value of the free surface $\phi (x,0)={\phi }_0\left(x\right)$

The solution of Eq.(1) with boundary conditions(2) is obtained using the finite element method, as described in [16]. Time marching is employed to solve the boundary equation and track the movement of the free boundary over time.

$$\phi (x,t+\delta t)=\phi (x,t)+\delta t[(\frac{K_s}{S_y}\nabla \varphi )·\mathbf{n}],$$

(4)

[12]. The movement of the free boundary function $\phi$ with respect to time is shown in Figure 2.

In order to investigate the effect of variable inflow levels on the unsteady-state flow, it is necessary to define the evolution of the water level with respect to time. The inclusion of the specific yield $S_y$ is important for our railway ballast problem [11], although not all previous studies have considered this term. Caution must be exercised when comparing our results with those reported in the literature. The importance of $S_y$ is discussed in [10]. The method for obtaining the free surface in Eq. 4 by dividing hydraulic conductivity by the specific yield was developed by France et al. [17] and Mills [18]. This method does not have any impact on the numerical difficulties encountered. Unless explicitly stated otherwise, we will assume $S_y$ to be unity.

3. Results for Rectangular Regions

We will commence by examining the three cases that have been previously explored by [12,13], in all of which the inflow height is constant. We will compute these three cases and conduct a thorough comparison with the aforementioned studies. The specific parameters for each case are outlined in Table 1.

The results of Case 1 are presented in Figure 3, where the free surface elevations obtained by our method are compared with those obtained by Di Nucci [13], who employed an approximate method to solve the same problem. A slight discrepancy between the free surface elevations is observed, with Di Nucci’s results yielding higher elevations than ours. However, such a difference is not surprising given that the two methods solve slightly different problems. Conversely, the results of Chakib and Nachaoui [12], depicted in Figure 4, exhibit discrepancies with our outcomes shown in Figure 3. Nonetheless, our results match closely with those of Di Nucci, as previously mentioned. We also note that our final free surface elevations align with those of Chakib and Nachaou, as depicted in Figure 5. We thus infer that the errors in Chakib and Nachaou’s outcomes stem from their time steps, while the consistency of our results with Di Nucci’s supports the accuracy of our method.

The results obtained for Case 1 enable a comparison of the free surface evolution, as shown in Figure 5. The figure includes the free surface results from this study, as well as those from previous investigations, including the steady-state results reported in [16,19]. Figure 6 displays the free surface convergences at $\phi (a,z)$ over the time interval $t\in [0,2]$, as well as the steady-state free surface. The variation of $\phi (a,t)$ is relatively constant for $0.5\le t$ in both this study and Di Nucci’s study, while Chakib and Nachaoui’s study shows a slightly different pattern, with the variation of $\phi$ being constant for $1\le t$, as evident from Figure 6. The solution converges to the desired accuracy when $t=2$, as illustrated in Figure 5 and Figure 6.

The free surface elevations for Cases 2 and 3 are presented in Figure 7 and Figure 8, respectively. The steady-state free surface solutions for these cases are computed through the iterative algorithm proposed in [9]. Our findings indicate that the results and conclusions for Cases 2 and 3 are comparable to those of Case 1.

3.1. Unsteady-State Flow Through a Rectangular Dam

In this study, we compare our results with the findings reported in [15], concerning the water flow through a homogeneous rectangular dam with a hydraulic conductivity of 1 m/s. The dam dimensions are 500 m in width and 12 m in height. In this scenario, the inflow water level rises abruptly from an equilibrium value of 7 m to 12 m at the first time step. The time steps used in this case are half an hour. This case simulates a sudden flow in saturated railway ballast.

The figure in Figure 9 displays the solutions for the free surface of a rectangular dam with a specific yield $S_y$ of 0.1 at different times. The dam has a width of 500 m and a height of 12 m with hydraulic conductivity equal to 1 m/s. Additionally, we included the analytical solutions of Ye et al. [15] to compare them with the solutions from our method. The free-surface elevations were computed at $t=0s$, $t=5$ days, $t=20$ days, and $t=50$ days using a time step of half an hour. Our results exhibit excellent agreement with those of Zuyang et al. [15]. The calculations for various time steps are presented in Figure 10. This figure demonstrates that both time step values yield the same outcomes for the free-surface position. Selecting a higher value for $\delta t$ during the calculation may cause a numerical instability error. Therefore, it is essential to choose an appropriate time step value to ensure good convergence. This issue is comprehensively discussed in Appendix A.

3.2. Unsteady-State Flow Through a Trapezoidal Dam

The porous medium of railway ballast can be accurately modeled as trapezoidal in shape. In this study, we aim to analyze the unsteady-state flow through such a trapezoidal porous region. To achieve this, we apply the same boundary conditions as presented in eq.2 but adapted to the trapezoidal shape, as depicted in Figure 11. The scenario under investigation involves a linear increase in the upstream variation, and the goal is to determine the free-surface position at different times. The obtained results will be carefully compared with previous studies in the field.

Desai et al. [14] conducted a comprehensive investigation of unsteady-state flow through various types of dams with differing hydraulic conductivities. Their study provides an opportunity to verify the validity of eq.4 using diverse scenarios for the upstream variation. In one of their experiments, a trapezoidal dam was subjected to a linear increase in upstream head from 0 to 15 cm, as depicted in Figure 12. Analysis of the resulting increase in water level indicates a rate of increase of approximately $1.28$ cm/min, according to [15]. Thus, the water level at 7 and 9 minutes was approximately $8.96$ cm and $11.52$ cm, respectively. This type of upstream variation mimics the flow of water through unsaturated railway ballast.

The experiment conducted by Desai et al. [14] involved a trapezoidal dam of length 34 cm and height 26.2 cm with a hydraulic conductivity of 0.1 cm/s and a specific yield of 0.4. They recorded the free-surface elevations at 7 and 9 minutes. In this study, we aim to compare our predictions for the free-surface elevation of the saturated region of the dam with their experimental observations. Initially, we assume a straight line at a height of 3.84 cm, representing the water level at 3 minutes. Using equation 4 with $\delta t=1s$, we determine the free surface. Our results, along with those of Desai et al. and Zuyang et al., are shown in Figure 13. The calculation of the free surface for this case is presented in Appendix B.

The results obtained from our analysis closely correspond with the experimental observations and the solutions proposed by Zuyang et al., as illustrated in Figure 13. While Desai et al.’s solution exhibited limitations after 9 min due to various factors that are discussed in [14], our method used a time step of $\delta t=1s$, and we assumed a downstream height of $H_2\left(t\right)=0.05$ cm as the boundary condition ${\Gamma }_3\left(t\right)$. It should be noted that the free-surface positions in our study are marginally lower downstream because of the maintenance of a downstream height of 0.05 cm. Consequently, we have established that our code is operational, and we have successfully completed the validation process.

4. Water Flow Through Realistic Railway Ballast

4.1. Steady-State Flow for Realistic Railway Ballast

Here, we investigate the time-dependent water flow through railway ballast, taking into account the realistic dimensions of railway ballast, as depicted in Figure 14. Initially, we solve the steady-state flow problem and compare our results with the findings of [9]. We consider the porosity of the ballast to be affected by fouling, as illustrated in Figure 15, following the work of [9,20]. Subsequently, we analyze the free surface positions for various models, including clean ballast, Model 1, Model 2, and Model 3 (refer to Figure 15), as they evolve over time. The results are presented in Figure 17.

The red dashed lines presented in Figure 17 depict the free surface positions calculated using the iterative algorithm, as proposed in [9]. The time step, $\delta t$, used in all models was set to $0.05s$. Interestingly, each model resulted in different times for the final position of the free surface. For instance, the free surface of the clean ballast was achieved after $86.4s$, while for Models 1 and 2, it was achieved after 7 minutes. Our findings indicated that the free-surface positions for the first three models closely corresponded to those reported in [9]. However, in the last model, convergence was slower, and the iterative algorithm did not converge even after 22 minutes.

Figure 16. Fouling ratio or Void Contaminant Index (VCI) and the equivalent hydraulic conductivity $K_s$ [9].

Figure 16. Fouling ratio or Void Contaminant Index (VCI) and the equivalent hydraulic conductivity $K_s$ [9].

Figure 17. Free surface of steady-state flow with respect to time for different models.

Figure 17. Free surface of steady-state flow with respect to time for different models.

The convergence rate of the time-dependent free surface in comparison to the steady-state free surface is significantly influenced by the hydraulic conductivity value, as demonstrated in Figure 17. The free surface of the clean ballast model converged at $86.4s$ due to its homogeneous material with high conductivity. In contrast, Model 3 is nonhomogeneous and complex, with regions of low conductivity, resulting in a much slower convergence rate. Therefore, we suggest utilizing the iterative algorithm proposed in [9] for nonhomogeneous hydraulic conductivity, but only if the steady-state solution is required.

4.2. Unsteady-State Flow for Realistic Railway Ballast

In this section, we model a sudden flood in railway ballast by considering a sudden rise in water level upstream. We incorporate the specific yield $S_y$ to compute the free surface at various times. Specific yield represents the capacity of railway ballast to drain water. For realistic railway ballast, we assume a specific yield of 0.2, as coarse sand and gravel typically have specific yields of up to 0.2 [11]. We acknowledge that specific yield may vary with contamination, but we keep it constant at 0.2 for all models. We investigate the sudden flow for clean ballast and Model 1 in Figure 15, where the effect of fouling ratio is similar to that depicted in Figure 16 for Model 1. We employ eq. 4 to determine the position of the free surface with respect to time. The hydraulic conductivity of the clean ballast is 0.3 m/s [9].

Figure 18 shows the results of our investigation on sudden flood scenarios for clean ballast and Model 1, where the sudden upstream change occurs from a stable water level height of 0.005 m to the height of the railway ballast of 0.3 m at the first time step of $\delta t_1$, while the downstream is held constant at 0.005 m. We compared the free surface of both cases at specific time intervals of 0.1, 1, 5, 10, and 90 s. The impact of fouling in Model 1 greatly affected the free surface compared to the clean ballast case, in which the fouling ratio was considered. This effect is supported by previous studies on railway ballast scouring, which primarily occurs when the upstream height exceeds the height of the ballast [6,7]. We also observed that the ballast fouling in Model 1 led to a slower flow of water for a significant time compared to clean ballast, which may cause overflow in a continuous flooding situation.

5. Conclusions

The solution of time-dependent flow with a variable free surface in fluid flow through porous media is a powerful technique that enables the calculation of steady- and unsteady-state flows. In this study, we have investigated both flow conditions and compared the results with experimental data. The steady-state flow results showed good agreement with the iterative algorithm proposed in [9], except in regions of highly variable hydraulic conductivity. For unsteady-state flow, we examined two inflow situations: linear and sudden increases. Specifically, we analyzed the latter case as it pertains to sudden floods that can occur in railway ballast. We demonstrated how fouling impacts the rate of water flow through ballast, as illustrated in Figure 18. Furthermore, we observed that the convergence of the free surface over time depends on the value of the specific yield, $S_y$, for the unsteady state, and that the time steps, $\delta t$, need to be chosen with care to obtain smooth curves for the free surface.