Continuum-Based Approach to Model Particulate Soil–Water Interaction: Model Validation and Insight into Internal Erosion

Abstract

Resolving the interaction between soil and water is critical to understanding a wide range of geotechnical applications. In cases when hydrodynamic forces are dominant and soil fluidization is expected, it is necessary to account for the microscale interactions between soil and water. Some of the existing models such as coupled Computational Fluid Dynamics–Discrete Element Method (CFD-DEM) can capture microscale interactions quite accurately. However, it is often computationally expensive and cannot be easily applied at a scale that would aid the design process. Contrastingly, continuum-based models such as the Two-Fluid Model (TFM) can be a computationally feasible and scalable alternative. In this study, we explored the potential of the TFM to simulate granular soil–water interactions. The model was validated by simulating the internal fluidization of a sand bed due to an upward water jet. Analogous to leakage from a pressurized pipe, the simulation was compared with the available experimental data to evaluate the model performance. The numerical results showed decent agreement with the experimental data in terms of excess pore water pressure, fluidization patterns, and physical deformations in violent flow regimes. Moreover, detailed soil characteristics such as particle size distribution could be implemented, which was previously considered a shortcoming of the model. Overall, the model’s performance indicates that TFM is a viable tool for the simulation of particulate soil–water mixtures.

1. Introduction

The interaction between soil and water is fundamental to understanding the mechanical behavior of soils. The state of stresses in the soil depends, to a great extent, on the interaction between the soil particles and the water found within pores (e.g., seepage and consolidation). In civil and geotechnical engineering, the relationships governing water–soil interaction in static and quasi-static cases are well established. A key assumption adopted in the majority of these relations is that no significant deformations occur in the soil matrix such that Darcy flow is valid. However, in more dynamic cases, when hydrodynamic forces become dominant (e.g., debris flow and internal fluidization), the soil integrity is most likely to be compromised and can no longer be characterized by conventional parameters such as hydraulic conductivity or permeability. In such cases, the need to carry out the analysis at the microscale level emerges, where accounting for the particulate nature of soil as well as the interaction forces between soil and water becomes necessary.

One of the examples of these dynamic interactions, which is the main focus of this study, is the internal fluidization and erosion in the soil around leaking water mains. Water leakage under pressure can lead to the erosion of the surrounding soil and decrease in effective stresses. This, in turn, can result in soil softening threatening the water main and the nearby structures [1]. Characterizing and predicting internal fluidization of soils can be particularly challenging considering the complex underlying physics that govern fluidization [2]. These physics involve the interaction mechanisms, momentum transfer, and force balance between soil particles and water, or other fluids [3]. Micromechanical analysis has been widely adopted in fluidized beds, where fluidization is induced by gas or liquid injection through an opening [4,5,6]. In ground engineering and earth structure, with a less controlled environment than that of fluidized beds, the analysis becomes more complex as the properties of soils such as cohesion and particle size variability come to play a major role [7,8]. Furthermore, if cyclic pressure is to be considered (e.g., hydraulic transient), more complex interactions are to be expected such as soil infiltration into the pipe accompanied by wash-out of soil particles in the vicinity of leakage.

Analyzing internal fluidization requires access to information related to the micromechanics of the system and resolving the particle–fluid interaction. Several attempts have been made to investigate these interactions experimentally [9,10,11], numerically [2,12], and using simplified analytical models [13,14]. Given the limited flexibility of experimental and analytical solutions in obtaining microscale information and dealing with a wide range of initial and boundary conditions, numerical analysis is mostly adopted. The two most common approaches are pure continuum analysis and coupled discrete-continuum analysis. In pure continuum analysis, the solid and fluid phases are treated as a single continuum with collective properties that account for the presence of both phases (i.e., a mixture). Alternatively, the two phases are interpreted as inter-penetrating continua each representing a single phase such as the case of the Two-Fluid Model (TFM) [15,16,17]. On the other hand, in coupled discrete-continuum analysis, the solid phase receives a discrete Lagrangian treatment while the fluid phase is considered as a continuum. Examples for this type of analysis involve the coupled Discrete Element Model–Computational Fluid Dynamics (DEM-CFD) [18] and Discrete Element Model–Lattice Boltzmann Model (DEM-LBM) [19].

Although utilizing discrete methods such as DEM-CFD allows for a relatively more accurate representation of the soil–water coupling compared to continuum-based approaches, it comes at a high computational cost [20]. On the contrary, coupled continuum modeling such as the TFM is much more computationally feasible with a code execution time of approximately 7% of that of CFD-DEM [21]. Thus, it allows for analyzing larger and more realistic systems and producing results that can be incorporated into conventional analysis needed for design and risk assessment. The model’s computational feasibility arises not only from its small execution time compared to discrete models but also from the flexibility in mesh refining that allows for directing the computational load to certain areas in the domain. This is particularly useful when dealing with large-scale systems where the mesh refinement might be required only locally. The application of the TFM, however, has thus far been limited to a few cases in simulating sediment transport and riverbed morphology [22,23]. Rigorous testing for the model’s potential to capture the basic characteristics of soil–water interactions is still needed.

In this study, we aimed to explore the potential of the continuum-based two-fluid model to simulate granular soil–water interaction in a variety of flow and soil conditions. Model validation was first performed, in order to confirm the model performance, by comparing the calculated results with experimental data related to the internal erosion of granular soil due to an upward water jet. The ability of the model to simulate internal erosion is then discussed. Finally, the response of the sand bed, including the evolution of excess pore pressure, fluidized zones, and particle size distribution, are investigated.

2. Numerical Analysis

The granular soil–water mixture is represented as interpenetrating continua of fluid and solids, each characterized by volume fractions of ${\epsilon }_f$ and ${\epsilon }_s$ for fluid and solids, respectively. For variations in soil particle sizes and subsequent interactions with the fluid phase to be accounted for, the solid phase can be further divided into a number (n) of sub-phases, such that

$${\epsilon }_f+{\displaystyle \sum _{i=1}^n{\epsilon }_{si}}=1$$

(1)

The governing equations for this mixture consist of the volume-average continuity and momentum equations proposed by Anderson and Jackson [15]. The volume-averaging of the governing equations for the solid phase results in a fluid-like resemblance of the solid phase, and hence the model is also known as the Two-Fluid Model (TFM). When there is the assumption that there are no mass flux sources or sinks, the locally averaged continuity equations are expressed as

$$\frac{\partial }{\partial t}\left({\epsilon }_f{\rho }_f\right)+\nabla \cdot \left({\epsilon }_f{\rho }_f{\mathrm{v}}_f\right)=0$$

(2)

$$\frac{\partial }{\partial t}\left({\epsilon }_{si}{\rho }_{si}\right)+\nabla \cdot \left({\epsilon }_{si}{\rho }_{si}{\mathrm{v}}_{si}\right)=0$$

(3)

where ${\epsilon }_f$ is the volume fraction of the fluid phase (e.g., porosity of saturated soil); ${\rho }_f$ and ${\rho }_{si}$ are the density of the fluid and ith solid phases, respectively; ${\epsilon }_{si}$ is the volume fraction of the ith solid phase; and ${\mathrm{v}}_f$ and ${\mathrm{v}}_{si}$ are the velocities of the fluid phase and the ith solid phase, respectively. For the fluid phase, the conservation of linear momentum appears in a similar form to that of a single-phase flow averaged by the fluid volume fraction. However, the outstanding difference is the inclusion of the particle–fluid interaction forces that account for the momentum transfer between the solid phase(s) and the fluid phase. The locally averaged momentum equation for the fluid phase is given as

$$\frac{\partial }{\partial t}\left({\epsilon }_f{\rho }_f{\mathrm{v}}_f\right)+\nabla \cdot \left({\epsilon }_f{\rho }_f{\mathrm{v}}_f{\mathrm{v}}_f\right)=-\nabla p_f+\nabla \cdot {\tau }_f+{\epsilon }_f{\rho }_f\mathrm{g}-{\displaystyle \sum _{i=1}^n{f}_{fs}}$$

(4)

where $p_f$ is the fluid pressure, ${\tau }_f$ is the shear stress tensor of the fluid, $\mathrm{g}$ is the gravitational acceleration, and $f_{fs}$ is the particle–fluid interaction force. The interaction forces can include a wide range of different forces depending on the scope of simulation [20]. For the interaction between soil and water in geotechnical engineering applications, buoyancy and drag forces are considered sufficient to describe the interaction forces. The particle–fluid interaction forces are expressed as

$$f_{fs}=-{\epsilon }_{si}\nabla p_f-F_d\left({\mathrm{v}}_{si}-{\mathrm{v}}_f\right)$$

(5)

where $F_d$ is the momentum transfer coefficient. The first term on the right-hand side of Equation (5) denotes the buoyancy force, while the second term is the drag force. Different expressions have been reported for the momentum transfer coefficient [24,25]. Here, we adopt the drag expression proposed by Syamlal and Obrien [26]:

$$F_d=\frac{3{\epsilon }_{si}{\epsilon }_f{\rho }_f}{4v_{rs}^2{d}_{pi}}{C}_d\left(\frac{{\mathrm{R}}_{es}}{v_{rs}}\right)\left|{\mathrm{v}}_{si}-{\mathrm{v}}_f\right|$$

(6)

where $d_{pi}$ is the average particle diameter of the ith solid phase, ${\mathrm{R}}_{es}$ is the particle Reynold’s number of the ith solid phase, $v_{rs}$ is the terminal velocity of the ith solid phase as given by Garside and Aldibouni [27], and $C_d$ is the drag coefficient [28]:

$$C_d={\left(0.63+\frac{4.8}{\sqrt{{\mathrm{R}}_{es}}}\right)}^2$$

(7)

The fluid shear stress tensor, under the assumption of a Newtonian fluid, is given as

$${\tau }_f=2{\epsilon }_f{\mu }_f{\overline{\overline{D}}}_f+{\epsilon }_f{\lambda }_f\mathrm{tr}\left({\overline{\overline{D}}}_f\right)\mathrm{I}$$

(8)

where ${\mu }_f$ and ${\lambda }_f$ are the fluid’s shear and bulk viscosity, respectively; $\mathrm{I}$ is the identity matrix; and ${\overline{\overline{D}}}_f$ is the strain rate tensor, such that

$${\overline{\overline{D}}}_f=\frac{1}{2}\left(\nabla {\mathrm{v}}_f+{(\nabla {\mathrm{v}}_f)}^{\mathrm{T}}\right)$$

(9)

Similar to the description of fluid motion, the momentum equation of the solid phase is given as

$$\frac{\partial }{\partial t}\left({\epsilon }_s{\rho }_s{\mathrm{v}}_s\right)+\nabla \cdot \left({\epsilon }_s{\rho }_s{\mathrm{v}}_s{\mathrm{v}}_s\right)=-\nabla p_s+\nabla \cdot {\tau }_s+{\epsilon }_s{\rho }_s\mathrm{g}+{\displaystyle \sum _{i=1}^n{f}_{fs}}$$

(10)

where $p_s$ is the equivalent solid pressure resulting from particle contact and collision, and ${\tau }_s$ is the shear stress tensor of the solid phase. It is important to note that the last term on the right-hand side of Equation (10) is essentially equal to that in Equation (4) in magnitude and opposite in direction such that Newton’s third law of motion is not violated.

The fluid-like representation of the solid phase in Equation (10) requires equivalent values for the pressure field, shear viscosity, and bulk viscosity similar to those of the fluid phase. However, obtaining these equivalent values is not as straightforward considering the inherently discrete (particulate) nature of the solid phase. Therefore, closures are needed to express the behavior of the solid phase (e.g., contact pressure and collision) in terms of continuous pressure and viscosity. One of the most adopted approaches to achieve this is the Kinetic Theory of Granular Flow (KTGF) [17,29,30]. Before incorporating the KTFG to obtain the needed closures, we first need to identify two main granular flow regimes: (i) plastic or frictional flow and (ii) viscous flow (Figure 1). In the plastic flow regime, particle contact forces are dominant, which is suitable for capturing the behavior of static or quasi-static systems where no significant deformations take place. In the viscous flow regime, solid particles are considered to be in a state of fluidization, that is, particles are no longer in contact, and their movement is only affected by collisions and interaction between the particles and the surrounding fluid.

In contrast to discrete element modeling, the particle contact cannot be identified explicitly due to the averaging technique adopted in continuum modeling. Thus, there needs to be a porosity threshold that separates the frictional regime from the viscous regime. This threshold is set to be slightly less than the initial packing volume fraction, ${\epsilon }_s^{*}$, at which the granular assembly is assumed to be in a static condition. For volume fractions greater than or equal to ${\epsilon }_s^{*}$, the plastic regime is assumed to be in effect (marked with the superscript $p$), otherwise the viscous solid regime is valid (marked with the superscript $v$). Under these conditions, the shear stress tensor of the solid phase is given as

$${\tau }_{si}^v=2{\mu }_{si}^v{\overline{\overline{D}}}_{si}+{\lambda }_{si}^v\mathrm{tr}\left({\overline{\overline{D}}}_{si}\right)\mathrm{I}$$

(11)

$${\tau }_{s1}^p=2{\mu }_{s1}^p{\overline{\overline{D}}}_{s1}$$

(12)

where ${\mu }_{si}^v$ and ${\lambda }_{si}^v$ are the shear and bulk viscosities of the ith solid phase in a viscous flow regime, respectively; ${\overline{\overline{D}}}_{si}=1/2\left[\nabla {\mathrm{v}}_{si}+{(\nabla {\mathrm{v}}_{si})}^{\mathrm{T}}\right]$ is the strain rate tensor of the ith solid phase; and ${\mu }_{s1}^p$ is the shear viscosity of the first solid phase in plastic flow regime. A summary of the expressions for the pressure field, and shear and bulk viscosities is given in

Table 1. The expressions for solid pressure, shear viscosity, and bulk viscosity of the solid phase after Jenike [31] and Schaeffer [32].

Parameter	Viscous Regime	Plastic Regime
$p_s$	$p_{si}^v=K_{1i}{\epsilon }_{si}^2{\mathsf{\Theta }}_i$	$p_{si}^p=A{\epsilon }_{si}{\left({\epsilon }_f^{*}-{\epsilon }_f\right)}^n$
${\mu }_s$	${\mu }_{si}^v=K_{3i}{\epsilon }_{si}\sqrt{{\mathsf{\Theta }}_i}$	${\mu }_{s1}^p=\frac{A{\left({\epsilon }_f^{*}-{\epsilon }_f\right)}^n\sin \left(\varphi \right)}{2\sqrt{I_2{}_{\left({\overline{\overline{D}}}_s\right)}}}$
${\lambda }_s$	${\lambda }_{si}^v=K_{2i}{\epsilon }_{si}\sqrt{{\mathsf{\Theta }}_i}$	N/A

, following the formulation of Jenike [31] and Schaeffer [32]. In Table 1, $\mathsf{\Theta }$ denotes the granular temperature, $A={10}^{25}$, $n=10$, $\varphi$ is the angle of internal friction of the granular assembly; and $I_2{}_{\left({\overline{\overline{D}}}_s\right)}$ is the second invariant of the strain rate tensor. The evolution of the granular temperature is given as [17]

$$\frac{3}{2}\frac{\partial }{\partial t}{\displaystyle \sum _{i=1}^n{\epsilon }_{si}{\rho }_{si}{\mathsf{\Theta }}_i}+\frac{3}{2}\nabla \cdot {\displaystyle \sum _{i=1}^n{\epsilon }_{si}{\rho }_{si}{\mathsf{\Theta }}_i{\mathrm{v}}_{si}}={\displaystyle \sum _{i=1}^n\left[\left(-p_s+{\tau }_s\right):{\mathrm{v}}_{si}-\nabla \cdot q_{\mathsf{\Theta }}-{\gamma }_{\mathsf{\Theta }i}+{\varphi }_{fi}+{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^n{\varphi }_{ij}}\right]}$$

(13)

where ${\gamma }_{\mathsf{\Theta }i}$ is the rate of granular energy dissipation due to particle collision, $q_{\mathsf{\Theta }}$ is granular energy flux, ${\varphi }_{fi}$ is the rate of energy transfer between the ith solid phase and the fluid phase, and ${\varphi }_{ij}$ is the rate of energy transfer between the ith and jth solid phases. In this research, the open-source code [33] was used to carry out the numerical simulation, incorporating the governing equations above. As the code contains a large set of sub-models and numerical schemes for calculating the solid friction and drag forces, we here adhered to the Syamlal–O’Brien drag model and Shaeffer solid friction model [32]. More information related to the sensitivity of the results to the sub-models can be found in MFiX documentation [33]. Additional details on the theoretical background and the numerical implementation can be found in Syamlal et al. [34].

3. Model Validation

The model validation was performed by simulating submerged sand bed subjected to upward water jet. The calculated results were compared with experimental data for the same setup.

Simulation of Internal Soil Erosion Around Leaking Pipes

Approaching the problem of internal erosion around a leaking source (Figure 2) is particularly challenging since the compromised soil is not visible in contrast to backward erosion resulting from high exit gradients [9]. From a theoretical point of view, the problem poses further challenges to accurately determine the velocities and pressure in the vicinity of an orifice with exit blockage. Moreover, it is numerically and theoretically challenging to capture the transition between initial seepage flows and potential fluidization that may occur at high inlet pressure. A few studies have attempted to numerically tackle the problem of internal soil fluidization due to upward water jet using various methods such as the coupled Lattice Boltzmann Method–DEM (LBM-DEM) [2,35] or continuum-based models [12], or using simplified analytical techniques [13]. A few issues, however, still need to be addressed in order to be able to practically incorporate these models in real-life problems. For example, in LBM-DEM simulations, relatively small systems are often adopted in order to avoid the high computational cost. For continuum-based and semi-analytical models, the flexibility of the model to capture flow transitions with realistic boundary conditions and soil characteristics (e.g., particle size distribution) remains questionable.

Alsaydalani and Clayton [11] carried out an experimental investigation of the internal fluidization of sand due to water injection from a rectangular slot, which can be fairly represented by a two-dimensional simulation. A similar investigation was carried out by Van Zyl, Alsaydalani, Clayton, Bird, and Dennis [10] using a circular orifice opening at the inlet instead of a rectangular slot. In this study, we limited our numerical investigation to a two-dimensional case; however, a three-dimensional simulation should be doable using the same governing equations and numerical tools reported in this study.

The experimental setup proposed by Alsaydalani and Clayton [11], which was adopted for the simulation in this study, consists of a box filled with submerged silica sand (sand bed) subjected to water injection at the bottom, similar to that shown in Figure 2. Water is injected at controlled pressures provided by a small pump through the rectangular slot opening, while flowrates are calculated from the water volume collected from the overflow at the top. Probes to measure the excess pore water pressure are placed along the centerline of the box and connected to sight-tubes to measure the development of excess pore water pressure throughout the experiment. Snapshots of the deformation of the sand bed are taken throughout the experiment to monitor the surface heave, the extent of the fluidized zone, and the patterns of sand boiling at high flowrates. For the numerical simulation considered in this study, the height of sand inside the box was set to 0.3 m with a width of 0.6 m. Sand was considered to be monodispersed with an average particle size of 0.9 mm and a slot opening of 0.62 mm.

The simulation was carried out using MFiX platform [33], which allows for the governing equations to be solved using implicit Euler scheme for temporal discretization. For convective acceleration terms, first-order upwind scheme with superbee flux limiter was used with a maximum continuity and momentum residual at convergence of 0.0001. For the numerical solution of the discretized equations, we used Biconjugate Gradient Stabilized solver (BICGSTAB). The simulation was carried out using a uniform mesh of size of 0.005 m × 0.005 m and a time step of 0.00001 s to ensure numerical stability of the simulation. However, for the computational mesh, the mesh was refined in the x-direction in the vicinity of the opening to a width of 0.2 mm to capture the response around the orifice. The total number of computational cells used in the simulation was 19,440 cells. This discretization was compatible with the mesh dimensions recommended by Tang, Chan, and Zhu [12] for the same simulation setting where further mesh refinement was found to not affect the results significantly. A summary of the simulation parameters is provided in

Table 2. Summary of the experimental parameters [11].

Parameter	Value
Bed width (m)	0.6
Orifice opening (mm)	0.62
Bed height (m)	0.3
Fluid density (kg/m3)	1000
Solid density (kg/m3)	2680
Initial porosity [-]	0.35
Applied upstream pressure (kPa)	10, 27, 60, and 190
Time step (s)	0.00001
Mesh cell size (m)	0.005 × 0.005
Temporal discretization	Implicit Euler
Particle friction angle (degrees)	30
Coefficient of restitution [-]	0.6
Solid friction model	Schaeffer (1987)
Drag model	Syamlal–O’Brien (1988)

.

Initially, the system was set up with a fluid volume fraction of 0.38, which was higher than the desired porosity prescribed in the experiment ${\epsilon }_f=0.35$, in order to allow for interparticle forces to properly develop as solid packing occurred. Afterward, the solid particles were left to settle down to the prescribed volume fraction of 0.65, such that the total height of the silica sand submerged underwater was 0.3 m. After the particles settled such that interparticle forces were properly developed and no initial motion was detected in the sand bed, water was injected at constant upstream pressure applied at the slot opening. Here, we chose to impose pressure boundary conditions instead of the velocity boundary condition adopted by Tang, Chan, and Zhu [12]. This was not to only align with the experimental procedure but also to provide a practical aspect to the simulation parameters. In the case of a local pressurized leakage, there would be many uncertainties related to estimating the leakage velocity, such as the size of the rupture and the soil conditions in the vicinity of the opening. On the other hand, the upstream pressure at the leakage condition could be easily determined from the pipe hydraulics. The boundary pressure values considered here were 10, 27, 60, and 190 kPa.

4. Results and Discussion

The results of the numerical analysis, including excess pore pressure, flowrate, and the onset of sand boiling, are summarized below.

4.1. Excess Pore Water Pressure and Fluidized Zone

The excess pore water pressure, represented by the difference between the calculated dynamic pressure and the initial hydrostatic pressure, are shown in Figure 3. The results show a good agreement with the experimentally reported values in terms of the magnitudes and distribution pattern of excess pore water pressure. At low upstream pressure (e.g., 10 kPa), pressure build-up was observed in the vicinity of the location of injection. This build-up of pressure continued until it was large enough to mobilize the sand particles, creating a small fluidized zone around the injection slot. Following fluidization, a drop in the pressure was observed in the fluidized zone, which indicated pressure relief accompanied by a movement of the peak of excess pore water pressure to the top of the fluidized zone where the soil was still intact. This was consistent with the porosity distribution shown in Figure 4 as the porosity peaked to 0.9 near the location of injection. It was also observed that the peak porosity value was not located right after the inlet but rather shifted upwards before it declined again and dropped sharply at the end of the fluidized zone. This can be interpreted as particles rearranging under gravity (i.e., soil self-healing) near the boundaries of the mobilized zone.

The upstream pressure required to initiate fluidization reported by Alsaydalani and Clayton [11] is approximately 100 kPa at a distance of 10 mm from the inlet and is characterized by the occurrence of a pressure drop at 53 mm from the inlet. In the current analysis, however, the occurrence of the first pressure drop was observed at upstream pressure of 27 kPa and a distance 11.3 mm from the inlet (Figure 5). Although this pressure seemed to be much lower than that measured in the experiment, it was noted that the experimental excess pore water pressure data were measured using probes, of which, the lowest was placed at a height of 10 mm and 53 mm from the inlet. At the inlet location, however, pressure measurements were not available, as placing a probe at the inlet will block the flow. Therefore, experimental detection of pressure drop initiation could only be obtained by comparing the first two probes, while pressure variations near the inlet were not accounted for. It was also observed that, above the location of the first probe, the simulation results seemed to agree quite decently with the experimental data (Figure 4). This suggests that fluidization had already started at an upstream pressure of 27 kPa but was not captured due to limitations in pressure measuring tools, where the pressure drop was reported exactly at the same distance of the second measurement probe. This was confirmed by the good agreement observed between the simulation and the experiment at upstream pressure of 190 kPa, at which the pressure drop occurred above 53 mm from the inlet location. A more important implication of the results was that the fluidization can happen within a limited domain around the location of the leakage. Although this was hardly observable at a relatively short distance from the leakage source, the increase in excess pore water pressure affected the entire sand bed, implying a significant reduction in effective stresses and shear strength of the soil.

4.2. Flowrate

Flowrate was calculated for each upstream pressure boundary value upon achieving a steady-state flow similar to that followed in the experiment. It was obtained by multiplying the inlet velocity by the width of the slot as for two-dimensional simulation. The results for flowrate compared to the experimental results are shown in Figure 6. It can be observed that the calculated flowrates deviated from the experimentally measured values by 29%, 10%, and 1.5% for upstream pressures of 10, 27, and 60 kPa, respectively. This discrepancy was most likely caused by the different approaches adopted in calculating the flowrates in both the model and the experiment. For the experimental results, the flowrate was obtained by measuring the volume of water passing through the bed over 5 min at steady-state flow condition. This essentially accounted for water flow through the entire sand bed during the 5 min period as opposed to the calculated values of inlet velocity multiplied by the slot width. Other reasons might arise from velocity averaging in the TFM and the two-dimensional simulation instead of a more accurate three-dimensional representation. Overall, the results indicated that the model can reasonably reproduce results comparable to those measured in the experiment. Interestingly, both the experimental data and numerical simulation indicated a relatively high increase in flowrate beyond 27 kPa, which, again, suggests that fluidization occurred at a much lower pressure than experimentally reported.

4.3. Sand Boiling at High-Inlet Velocities

Another important aspect of the model’s performance that was tested here was its ability to capture high deformations and fragmentations. An example of this is sand surface boiling as a water jet of high velocity penetrates through the sand bed. The experimental results reported by Alsaydalani [9] included a qualitative description of surface boiling in a shallow sand bed. The experimental setup and simulation parameters were the same as those used previously, except for the sand bed height and slot opening, which were reduced to 150 mm and 0.33 mm, respectively. In this simulation case, a velocity boundary condition was imposed at the inlet to reproduce the flowrate boundaries of 700 L/h, 900 L/h, and 1200 L/h. For each case, the simulation results were displayed at the instant of water jet penetration out of the sand surface.

The results in Figure 7 show the physical deformation of the sand bed compared to the photographs taken for the corresponding inflow values in the experiment [9]. It can be seen that the model could capture the basic characteristics of the bed deformation to a good extent. At inlet flowrates of 700 L/h and 900 L/h, the central plunging zone and the soil heave on both sides could be well-replicated by the model. The same behavior was also observed in the case of inflow of 1200 L/h; however, the model was not able to replicate the sand fragmentation around the central plunge. This shortcoming stemmed from the numerical diffusion involved with the discretization of the convective acceleration terms.

4.4. The Effect of Polydispersity

Polydispersity is an essential aspect that should be considered in conducting microscale analysis of granular soils. It is difficult to implement the variations in particle size distribution of soil when continuum-based approaches are adopted. This is mainly because of the volume averaging techniques that account for the effect of particle size implicitly as contributing factors to the interparticle and particle–fluid interaction forces. However, when the governing equations are written in terms of multiple solid phases as presented in this study, accounting for polydispersity becomes possible through assigning various particle sizes to the division of solid sub-phases. The sum of volume fractions of the fluid phase and the solid phases should sum up to 1 as per Equation (1).

The numerical results presented thus far in this study only consider mono-dispersed granular soils. Although the effect of particle size variations has been reported both experimentally and numerically [11,12], uniform particle size was always used for each simulation case. This uniform particle size is characterized by an average value or more precisely the 50th percentile of particle sizes (D₅₀). The outstanding question is how representative this uniform value is of a granular assembly, especially with different assemblies that share the same D₅₀. To answer this question, we synthesized three different particle size distributions (PSD1, PSD2, and PSD3) such that they all had the same D₅₀ (Figure 8). The particle size distributions were designed to represent three different soil types. The distribution PSD1 represented a well-graded soil with particle diameters ranging from 0.5 to 3 mm. The distribution PSD2 was also well-graded; however, it had a narrower range of particle sizes to mimic fine sand. The same particle size range of PSD1 was used for PSD3, but instead it was characterized by two particle sizes of 0.5 and 3 mm to represent gap-graded soil. Three cases were conducted using the same system described earlier for the three particle size distributions with a slot opening of 3 mm. A summary of the synthetized particle size distributions is provided in

Table 3. The diameter and volume fraction allocated for particle size distributions PSD1, PSD2, and PSD3 for overall porosity of 0.4.

	PSD1		PSD2		PSD3
	Diameter (mm)	Volume Fraction	Diameter (mm)	Volume Fraction	Diameter (mm)	Volume Fraction
solid 1	0.5	0.153	0.85	0.153	0.5	0.288
solid 2	0.9	0.147	0.9	0.147	0.9	0.012
solid 3	1.2	0.153	1	0.153	1.2	0.012
solid 4	3	0.147	1.5	0.147	3	0.288

. At the inlet boundary, a water jet of a constant velocity of 3 m/s was set. The inlet velocity was set to be relatively high such that water jet penetration through the entire bed was achieved.

The evolution of porosity and velocity distributions across the sand bed is shown in Figure 9 and Figure 10. Fundamental differences can be seen in the physical deformation of the bed and the velocity field. For PSD2, where well-graded fine particles existed, fast propagation of the water jet was observed with a narrow fluidization zone and minimal lateral spread. For PSD1 and PSD2, where larger particles appeared to be more dominant, slower propagation of the water jet was observed with a wider lateral spread of fluidization. The width of the fluidized zone was maximum in PSD3, which represented gap-graded soils with the largest portions of both fine and coarse particles of all three distributions. Although it was evident that the three distributions displayed different behavior, it was not clear what was exactly the underlying parameter. Tang, Chan, and Zhu [12] attribute the slow jet propagation when large particles are dominant to an increase in soil permeability that helps to dissipate the excess pore water pressure. This might be valid for mono-dispersed systems; however, a quick look at the evolution of pressure field in Figure 11 shows no major differences between the three particle size distributions. One possible reason for the differences in behavior could have been the wash-out of finer particles as they were easier to mobilize. This aligned the lateral spread observations under nearly the same pressure distribution, especially that a full cavity was not developed, rather a region of higher porosity. Another interesting observation was seen in the velocity field that was the development of a trapezoidal-like wedge that characterized the mobilized zone. This is consistent with the assumption of Cui, Li, Chan, and Chapman [2] in their analysis for the critical inlet velocity. However, the breadth of the wedge was found to depend also on the particle size distribution, and therefore a modifier should be implemented to account for such variation if a single-value particle size is to be used.

5. Conclusions

Continuum-based modeling is seldom applied to coupled microscale simulation of coupled soil–water mixtures. In this study, we presented a model validation as well as a series of numerical simulations to test different aspects of the model’s performance and limitations. The following conclusions can be identified:

• The numerical results presented show the model’s capability to capture the basic flow characteristics of soil–water mixtures as well as transitions in flow regime.
• The model could successfully reproduce the evolution of internal soil fluidization due to local leakage.
• More violent flows were also decently captured. However, a shortcoming is the model’s inability to capture soil fragmentation accurately due to numerical diffusion.
• Polydispersity could be successfully implemented, and the simulation results suggest that using a mono-dispersed assembly is not sufficient to describe the soil–water properly.
• Further testing and validation of the model is still required to ensure its robustness in dealing with geotechnical applications.
• Overall, the model shows good agreement with theoretical solutions and experimental data and can be a viable tool for investigating particulate soil–water flows in geotechnical applications.