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 medium. 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 of the ballast, even if the upstream height is less than the railway ballast’s height [
8]. The 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 the realistic flow through porous media and the time-dependent change in the position of the free surface. We must also consider the specific yield
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
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 obtained from three cases in [
12,
13], 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 Modeling of Flow through Ballast
The region under consideration is depicted in
Figure 1 and can be defined by the domain
. The saturated part of the region is denoted by
, which is a function of time
t, where
. Here,
represents the water level inside the rectangular ballast. The free surface within the porous medium is denoted by
, while
and
represent the upstream and downstream regions, respectively, with a fluid boundary. The seepage face is denoted by
. The water heights are represented by
and
, while the impermeable boundary is represented by
, where there is no flux.
The governing equation is derived from the divergence of the seepage velocity and is given by
where
and
are the hydraulic conductivity and the hydraulic head, respectively [
16]. We use the following boundary conditions [
12]
The final condition is that for the movement of the free surface,
where
is the specific yield. We will discuss this term shortly. We also need to specify the initial value of the free surface
The solution of Equation (
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:
in Ref. [
12]. The movement of the free boundary function
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
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
is discussed in [
10]. The method for obtaining the free surface in Equation (
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
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 which the inflow height is constant for all three. 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 with 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, as 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
over the time interval
, as well as the steady-state free surface. The variation in
is relatively constant for
in both this study and Di Nucci’s study, while Chakib and Nachaoui’s study shows a slightly different pattern, with the variation in
being constant for
, as evident from
Figure 6. The solution converges to the desired accuracy when
, 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 flow of water 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 through saturated railway ballast.
Figure 9 displays the solutions for the free surface of a rectangular dam with a specific yield
of 0.1 at different times. The dam has a width of 500 m and a height of 12 m, with a 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
,
days,
days, and
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
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 Equation (
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 Equation (
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
cm/min, according to [
15]. Thus, the water levels at 7 and 9 min were approximately
cm and
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 34 cm in length and 26.2 cm in height, 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 min. 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, and minimizing Equation (
4) with
, 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
, and we assumed a downstream height of
cm as the boundary condition
. 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.
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 these 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 with 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,
, for the unsteady state and that the time steps,
, need to be chosen with care to obtain smooth curves for the free surface.