Investigating soil fluidisation using a coupled TFM-KTGF model

$CD$

drag function

$do$

crack width

$dp$

mean particle diameter

$eqq$

restitution coefficient for phase $q$

$ess$

coefficient of restitution

$g→$

gravitational acceleration

$g0$

radial distribution function

$g0,qq$

radial distribution function for phase $q$

$g0,ss$

radial distribution function for solid phase

$hf$

height of fluidised cavity

$hs$

soil bed thickness

$hw$

water head

$I¯¯$

unit tensor

$Ksl$

drag coefficient

L

length of domain modelled

$Ps$

solid pressure

$p$

fluid pressure

Q

leakage rate

q*

normalised leakage rate

$Rep$

particle Reynolds number

$Rpq→$

interaction force between phases

$uq→$

velocity of phase q

$V$

total volume

$Vl$

volume of liquid phase

$Vs$

volume of solid phase

v

velocity of leaking fluid

$vpeak$

peak velocity of solid phase

$αl$

volume fraction of liquid phase

$αs$

volume fraction of solid phase

$αs,max$

maximum packing limit

β

interphase momentum exchange coefficient

$λq$

granular bulk viscosity of phase q

$λs$

bulk viscosity

$Θq$

granular temperature of phase q

$Θs$

granular temperature of solid phase

$μs$

total particle shear viscosity

$μs,kin$

kinetic viscosity

$μq$

granular viscosity of phase q

$μq,col$

collisional viscosity

σ_o’

initial effective stress

$ρq$

density of phase $q$

$ρs$

density of solid particles

$τ¯¯$

Reynolds stress tensor

$ϕ$

angle of internal friction

Water distribution pipelines are vulnerable elements in any urban infrastructure. Leakage of these pipes can result in fluidisation of the surrounding soils, flooding and a terrible disaster. There are several incidents that show the severity of the above impact. For example, the 2005 Helsingborg pipeline network failure softened the ground, resulting in the collapse of large acid storage tanks (Cui Beng, 2012). Kuliczkowska’s (2016) work analysed over 100 cases of road collapse in Poland due to pipe defects. These defects destabilised their surrounding soils, leading to pavement distress. Chai et al. (2018) analysed an accidental collapse above a tunnel excavation. It was observed that the collapse occurred due to the formation of high water pressure channels in poor soils. A water pipeline burst in 2023 fluidised the backfill soils to form a large cavity in Baltimore, which resulted in extensive traffic disruption along its roadway (Shiau et al., 2021). Failure of water supply pipelines is the major anthropogenic cause of either sinkhole formation or road collapse in urban environment (Chen et al., 2023; Dave and Juneja, 2023a,b; Wang et al., 2023). Bursting of pipe results in abrupt failure of the system under high-speed fluid jet. They form internal cavities, which reach up to the ground surface (Schulz et al., 2021; Van Zyl et al., 2013). A small but continuous leakage is often difficult to detect but it also weakens the surrounding soil, causing pipeline instability (Cui Beng, 2012).

Fluidised cavities are high-velocity zones containing vortex of rapidly rotating particles (Niven and Khalili, 1998; Philippe and Badiane, 2013; Van Zyl et al., 2013). Granular soils are more susceptible to fluidisation (Cui Beng, 2012; Jiang et al., 2023). Soil particles in non-plastic granular beds have the ability to reorganise and form smooth elliptical cavities. They transform into fluid due to the drag of the flowing liquid. Low permeable clays are normally avoided as bedding material, but when used, they resist fluidisation and instead, undergo hydraulic fracture due to preferential flow paths. The factors that affect soil fluidisation are pipe pressure and size of the defect, bed thickness, shape and size of the soil particles, and the ground water table level. In some of the previous investigations, leakage rate ranged between 15 and 250 L/min, and particle size varied between 0.25 and 1 mm, with up to 0.4 m thick soil bed (Akrami et al., 2023; Alsaydalani, 2010; Bailey and Van Zyl, 2015; Cui Beng, 2012; He et al., 2017; Jiang et al., 2023; Tang et al., 2017). These studies focused on assessing the variation in excess pore pressures within the soil bed (Alsaydalani, 2010). The critical leakage rate that initiates the formation of a fluidised cavity was still not well understood (He et al., 2017; Tang et al., 2017).

This study uses continuum-based two-fluid method (TFM), coupled with the kinetic theory of granular flow (KTGF), to model soil fluidisation under monotonic leakage. TFM or Euler–Euler method considers the solid and liquid phases as interpenetrating continua. Therefore, a sharp interface does not exist as in volume of fluid modelling technique which is used to track and capture the interface between two or more immiscible fluids. Instead, the interface is implicitly represented by the spatial variation of the phase volume fractions, where momentum exchange terms govern the interaction between the two phases. The volume fractions are assigned for both the phases, which is the ratio of volume of a particular phase to total volume. The sum of volume fraction of both the phase is 1.

Volume fraction of solid phase, $αs$

where V_s is volume of solid phase and V is the total volume. Similarly, $Vl$ and $αl$ are volume and volume fraction of liquid phase, respectively.

TFM uses macroscopic conservation equations that are valid throughout the flow domain. Since the concept of particles has disappeared completely in such a modelling, the effect of particle–particle interactions can only be included indirectly through KTGF, which expresses the pressure and the solids stress tensor as a function of the granular temperature, which is a measure of randomness of the particles (Gidaspow, 1994). Thus, combining KTGF with TFM considers solid phase as a fluid-like continuum, but applies constitutive relationships for viscosity and pressure, and therefore can provide certain key parameters like particle pressure and viscosity, and transport coefficients (Gidaspow, 1994). TFM-KTGF approach is computationally cost effective when the interaction within and between the phases plays a significant role in determining the hydrodynamics of the system. This method is suitable for simulating large-scale problems, requires fewer ad hoc adjustments and is computationally cost effective compared to other Lagrangian methods (Crowe, 1991; Haq et al., 2022; Ibrahim and Meguid, 2021).

Coupling TFM with the kinetic theory of granular flow is an effective promising framework for modelling multiphase flow problems at engineering scale (Ibrahim and Meguid, 2021). As the discrete character of the dispersed phase is lost due to the averaging procedure, appropriate closure models for drag law, granular temperature and coefficient of restitution can be applied for accurate solution (Moliner et al., 2019). It is effective in controlling the boundary conditions, and one has the ability to extract data at any step of the transient analysis. This method’s use in simulating fluidised bed in reactors is well established, which makes this tool suitable for simulating leakages induced cavity formation (Abdulrahman et al., 2022; Moliner et al., 2019). Tang et al. (2017) used TFM-KTGF framework to determine the effect of particle size and bed thickness on the expansion of granular bed. The increase in particle size was shown to significantly reduce the soil-bed expansion. Ibrahim and Meguid (2021) used TFM to model a polydisperse assembly of granular particles. They identified the fluidised cavity zone, in which the volume of solid fraction, α_s, reached zero. The effects of viscous and friction models were assessed, providing valuable insight into the use of continuum approach.

In two-fluid model, the continuity and momentum equations are solved for each phase. The equation of continuity equation for each phase q is (Euler 1757):

where q can be fluid, l or solid, s and α_s = 1 − α_l. ρ_q and $uq→$ are the density and velocity of phase q, respectively.

This equation states that the rate of change of mass of phase q in a control volume equals the net mass flux through its boundaries. Since there is no interphase mass transfer occurring in fluid–solid flows, the right-hand side is zero.

The equation for conservation of momentum equation for fluid phase is

where p is fluid pressure, $g→$ is the gravitational acceleration, $τ¯¯$ is the Reynolds stress tensor, $Rpq→$ is the interaction force between phases, $ρl$ is the density of liquid phase, $ul→$ and $us→$ are velocity of liquid and solid phase, respectively.

here $I¯¯$ is the unit tensor, μ_q is granular viscosity and $λq$ is granular bulk viscosity. μ_q is a summation three: collisional, kinetic and frictional viscosities. $λq$ represents the resistance the granular particles against compression or expansion.

The momentum equation for solid phase is similar to fluid phase with an additional solid pressure term, P_s as

The momentum equation contains the pressure gradient term $(−αs∇p)$ which accounts for the effect of the common pressure field shared by both liquid and solid phases; it drives the flow from regions of high to low pressure. The stress divergence term $(∇.τl¯¯)$ represents the internal stresses within each phase. For the liquid phase, this mainly includes viscous stresses, while for the solid phase it also includes additional stresses arising from particle–particle collisions and frictional contacts, as described by KTGF. The body force term $(αgρgg→)$ denotes the gravitational force acting on each phase, proportional to its local volume fraction and density. Finally, the interphase momentum exchange term $β(us→ − ul⇀)$ captures the interaction between the liquid and solid phases through drag forces, and sometimes through lift, virtual mass or turbulent dispersion effects. This term transfers momentum between phases, ensuring that the motion of one phase influences the other, and is central to coupling the two sets of equations in the TFM framework.

For modelling leakage-induced soil erosion, the continuity equations track how water and soil volume fractions evolve in space and time. Momentum equations describe how fluid pressure gradients and drag move the soil grains into the pipe. KTGF closures ensure the model captures dense to dilute transitions as soil is eroded.

Figure 1 represents the schematic diagram of the scenario considered for the present study. The model used in this study is a simplified representation of real field conditions. Since the two-fluid approach can handle only two phases, all regions outside the sand bed were approximated as water, resulting in a simplified domain compared to actual site conditions. A finite volume based software (Ansys, Fluent 2021R1) was used for Eulerian multi-fluid analysis. The problem was modelled as plane strain to simulate a longitudinal crack of width, d_o, located at the pipe’s crown. A system of saturated sand beds of uniform gradation and mean particle size, d_p, and bed height, h_s, were used. The particles were assumed to be spherical. The effect of peak internal friction angle, ϕ, was considered indirectly, using the closure of frictional viscosity. The different cases used in the parametric study are shown in Table 1. Frictional viscosity model was applied in all the cases except for runs 3, 6 and 11 in the table. Figure 2 shows the two-dimensional geometry and a typical mesh used in this study. The figure indicates length, L, of the numerical domain, soil bed thickness, h_s, crack width, d_o, and water head, h_w. The domain, except for the solid bed, was modelled as water. The number of nodes and elements in the mesh were dependent on h_s and d_o. A velocity inlet boundary condition was applied at the crack, as both the inflow velocity magnitude and direction were known. This condition is suitable for compressible and incompressible flows. The leakage rate through the crack was calculated through the continuity equation which expresses volumetric flow rate as the product of cross-sectional area and average flow velocity, v. The side walls were assigned a no-slip boundary condition, assuming zero tangential velocity along the fluid–wall interface. The top boundary was defined as a pressure outlet with atmospheric pressure. This condition is appropriate when the outlet static pressure is known and is particularly robust in situations where flow reversal or recirculation may occur near the outlet. The mesh was refined in the crack region and along the soil–water interface, while a coarser mesh was used farther away near the walls. The behaviour of particle-to-particle collision was defined using a restitution coefficient, e_ss, to represent inelastic collision of adjacent particles. Radial distribution function, g₀, provided a strong control over the particle volume fraction using the maximum packing limit, α_s,max, to achieve a precise simulation of the flow behaviour. Particle phase stresses were computed using the total particle shear viscosity, $μs$⁠, and bulk viscosity, $λs$⁠.

Properties of LBS-B sand was used to model the granular phase, as this soil has extensively been studied in the literature and its properties are well established (Alsaydalani, 2010; Bolton, 1986). The sand had d_p of 0.9 mm, specific gravity of 2.65 and ϕ of 35°. Its e_ss was 0.9. Chen and Yan (2021) suggested that α_s,initial be above 0.5 for the dense sand beds. Therefore, the initial solid volume fraction, α_s,initial, was assigned as 0.65 for all the cases in Table 1. This corresponds to initial porosity of 0.35 of the bed. The value of α_s,max was calculated as 0.684 using the minimum possible void ratio of 0.46 for LBS-B (Alsaydalani, 2010). These values simulate dense sand bed of 78% relative density. Leakage rates of 2 to 30 L/min in the table show the possible field behaviour of internal pipe pressures of 20–30 kPa.

Standard k–ε turbulent viscous model was used for all the cases. In KTGF-based two-fluid modelling, the granular-phase viscosity is the sum of kinetic, collisional and frictional components. These terms are derived from KTGF and are not meant to artificially force the phase into a static state. A static granular phase is not represented by assigning an extremely large viscosity. Instead, the granular shear viscosity naturally becomes very high when the solid volume fraction approaches the maximum packing limit, which strongly suppresses deformation. Thus, the phase behaves quasi-static due to the constitutive model itself, not because a manually imposed large viscosity is used.

Kinetic viscosity represents the momentum transport caused by the random motion of particles. It represents mean motion which is captured by solid velocity $us$ and random fluctuating motion which is quantified by granular temperature $Θ$⁠. Kinetic viscosity, μ_s,kin was calculated using Syamlal O’Brien’s model (Syamlal et al., 1993), which is

where Θ_q is the granular temperature, e_qq restitution coefficient for phase q, and g_0,qq is the radial distribution function for phase q (solid or fluid) which represents the repulsive function between molecules.

The collisional viscosity, $μq,col$⁠, is a viscosity contribution due to collisions between particles is taken from the kinetic theory of granular flow (Lun et al., 1984)

To model dense flow under low shear stress, that is, when stresses are generated by the friction between the particles, Schaeffer’s (1987) frictional viscosity model was applied. The frictional viscosity is the contribution of the friction between particles to the total shear viscosity. Kinetic theory assumes binary and instantaneous collisions between particles, but at high solids concentration, it fails to adequately describe the granular flow. This necessitated the use of a suitable frictional viscosity model to predict the hydrodynamics of a dense granular system. When the solid volume fraction, α_s, gets close to the maximum packing limit, α_s,max, the particles get very close to each other. Schaeffer’s model used the peak friction angle obtained from Mohr–Coulomb yield criterion. This frictional viscosity, μ_s,fr, model was implemented using Equation 10

where P_s is the solid pressure and equals the sum of the kinetic pressure due to the motion and inelastic collision between the particles using a restitution coefficient, and I_2D is the second invariant of deviatoric stress tensor.

Solid pressure, P_s, by Lun et al. (1984) used in the study is

where $αs$ is the solid volume fraction, $ρs$ is density of solid particles, $ess$ is the coefficient of restitution, $g0,ss$ is the radial distribution function and $ϴs$ is the granular temperature.

Gidaspow et al.’s (1992) drag model was used to simulate the interphase momentum transfer between the liquid and solid phases, where the interphase momentum exchange coefficient, β, is

where drag function, C_D, in Equation 12 is

where Re_p is particle Reynolds number.

The Gidaspow drag model is a hybrid drag correlation that blends two well-known drag models depending on the local solid volume fraction, α_s. This correlation is a combination of the works of Ergun (1952) and Wen and Yu (1966). Ergun equation is valid for dense particle assemblies, while Wen–Yu correlation is valid for dilute suspensions. Gidaspow et al. (1992) proposed using a switch at α_s = 0.2 (i.e. α_l = 0.8 for fluid), because the physics of drag are very different in dense and dilute conditions. It has been shown that this drag model is applicable for modelling dense solid beds (Busch and Johansen, 2020).

The programme was run until steady state, that is, when the fluidised cavity volume stopped increasing. This happened when the zone of solid volume fraction was less than the initial solid volume fraction, that is, when α_s < α_{s, initial} (Ibrahim and Meguid, 2021). For pressure–velocity coupling, phase-coupled SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm was adopted (ANSYS Inc., 2021). The maximum number of iterations were 35, although only 5–10 iterations were sufficient to achieve convergence. The simulations were conducted using a 32 GB RAM system with Intel Core I7-2.9 GHz processor.

Figures 3(a) and 3(b) show the profile of solid volume fraction and its phase velocity with bed height for laminar model and Reynolds averaged Navier–Stokes (RANS k-ε) model, respectively. In this sensitivity analysis, the unsteady laminar and turbulent RANS k-ε formulations for mass and momentum were solved separately for fluid and solid phases. Solid phase was modelled with KTGF that describes the fluctuation and collision between the particles. The momentum equations for the fluid and solid phases were obtained from modified Navier–Stokes equations using the terms for inter-phase momentum transfer. Erosion hydrodynamics were unlikely to be sensitive to the choice of flow models, as the solid volume and the solid velocity distributions are similar for laminar and RANS k-ε models in this problem (Abdulrahman et al., 2022). Therefore, only RANS k-ε model was used in all subsequent fluidisation studies.

Figures 4(a) and 4(b) show the effect of frictional viscosity model on the solid volume fraction and vertical velocity of solid phase, respectively, for 2–30 L/min flow. At high leakage rates of 30 L/min, cavity height was over-predicted by 20% in the absence of frictional viscosity model. The trend of solid-phase velocity remained unaffected for runs 3, 6 and 11 but achieved higher velocities compared to the runs 2, 5 and 10, when the frictional viscosity model was activated. The influence of frictional viscosity was found to diminish at low leakage rates. Frictional viscosity contributed to the friction between particles to the total shear viscosity. Kinetic theory assumes binary and instantaneous collisions between particles, but at high solids concentration, it fails to adequately describe the granular flow. This necessitated the use of a suitable frictional viscosity model to predict the hydrodynamics of a dense granular system. Frictional viscosity was therefore used in dense regimes as α_s was greater than 0.5, to determine the momentum exchange (Busch and Johansen, 2020).

Figure 5 shows the contours of solid velocity vectors captured at different time intervals. An upward water jet was formed to create a high velocity central region, which was surrounded by a slower circulating region. As shown in the figure, soil bed’s active region has a taper shape, with its cross-sectional area increasing from bottom to the top (Alsaydalani, 2010). The solid phase was entrained as a group to move up in the high-velocity central region at the start of fluidisation. The liquid phase forms eddies, which intensified the mixing process. The bed height increased as the solid phase move up with the flow until the flow reached a critical height. At this stage, the velocity of the liquid phase was reduced. The solid phase then took the shape of an umbrella, only to move in the opposite direction. The velocity vectors were then forced to enter the slow circulating region, to return back to the crack. The solid phase interacts and collides in its downward journey. The solid velocity and the solid concentration are symmetrical around the crack width in the figure. As also observed in Philippe and Badiane’s (2013) experiments, clockwise motion occurs on the one side of the water jet, while a counterclockwise motion occurs on the other side. These cycles continued for 20 s, after which a turbulent convection roll developed, allowing the fluidised cavity to collapse.

The shape of the vertical velocity profile of the solid phase is that of a tapered fluidised bed, which can be seen from Figure 6(a). The solid velocity profile shows a rapidly increasing starting zone to reach a plateau, only to rapidly decrease thereafter. The solid phase accelerates under high fluid velocity in the increasing zone near the crack. The contours gradually taper within the upper spout region. This shape remains similar irrespective of the use of different computational fluid mechanics tools (Moliner et al., 2019).

Figure 6(b) shows the variation of the solid volume fraction measured along the vertical plane above the crack. Here, α_s is zero at the crack location and achieves α_s,max of 0.684 above the cavity. The latter helps to compact the sand above the cavity. α_s is zero at and above the solid–water interface.

Figure 7 compares the cavity height to the leakage rate. Both the axes in the figure are normalised. Normalised leakage rate, q*, is represented by Equation 15 (Montellà et al., 2016)

where Q is in mm²/s and σ_o′ is the initial effective stress in Pa. The value of 0.001 in Equation 15 indicates the dynamic viscosity of water in Pa.s at 20°C. The figure shows that h_f/h_s is independent of q^* beyond a certain threshold. Indeed, the cavity height obtained in this study falls within the range of Alsaydalani’s (2010) and Schulz et al.’s (2021) experimental data. The data occupies a reasonably thick independent band between 0.1q^* and 1q^*, which further strengthens the reliability of the method used in this study. Figure 8 shows the effect of flowrate on the cavity development. It is a visual representation of fluidised cavity, where the x-axis shows the corresponding flowrate values and the grey colour region represents cavity.

Figures 9(a) and 9(b) show the effect of d_o and h_s on the height, volume and solid phase velocity of the cavity for leakage rate ranging from 0.5 to 30 L/min leakage through a 0.3 mm wide crack. TFM-KTGF framework uses particle diameter to calculate the drag coefficient, terminal velocity and inter-particle collisions. Any predefined solid volume fraction will account for the number of particles in the system (ANSYS Inc., 2021). Figure 9(a) shows that the cavity height reached the bed thickness for the bed of small particles. The bed was fluidised to form chimney-shaped blow-out failure when d_p was within 0.1–0.3 mm. The line showing blowout failure separates the region of blowout and hidden cavity. The cavity remained within the bed when d_p was greater than 0.3 mm. With the increase in d_p from 0.3 to 2 mm, the cavity height reduced drastically. Its shows that large particles sustained the cavity, even within the bed of 150 mm thickness for high leakage rates of 30 L/min. The cavity width, on the other hand, increases with increase in d_p. Van Zyl et al. (2013) also observed a similar cavity when Q ranged from 12 to 14 L/min in a thin bed of 1 mm dia. glass beads. They called this as the ‘hidden cavity’. It was found that the excess water pressure was barely 2.5 kPa even when the pipe pressure was over 235 kPa (Van Zyl et al., 2013). He et al. (2017) too did not find any blowout in 400 mm thick sand bed using d_p of 0.25 mm and Q of 20 L/min. Large particles seem to rapidly transfer the momentum, resulting in a faster dissipation of the energy. This resulted in small jet penetration and a small cavity. The cavity width did seem to increase with the increase in particle size, which shows wider spread of the water leakage. Figure 9(b) shows the velocity profile of the solid phase along the vertical axis. With the increase in d_p, the peak solid phase velocity decreases which results in reduction in momentum transfer. This seems to be the reason for the decrease in cavity volume.

Previous studies have shown the insignificant effect of bed height on fluidised cavity (Hong et al., 1997; Zhong and Zhang, 2005). Tang et al. (2017) assessed the reduction in the expansion of bed by 5% when the bed height was doubled. He et al. (2017) found that the onset cavity flowrate remains independent of the soil bed thickness. It is generally suggested that pipelines be installed at shallow depths to allow quick flow of water to reach the surface so that leaks can be found and repaired quickly (Alsaydalani, 2010).

Figures 10(a)–10(d) show the effect of d_o on the cavity dimensions for d_p of 0.9 mm and h_s of 150 mm for different Q. The range of crack width was chosen from the previous experimental studies of longitudinal pipe failures (Alsaydalani, 2010; Shao et al., 2019). The value of h_f reduces with increase in d_o. Narrow cracks produce high velocity fluid jet due to their small cross-sectional area, causing the cavity to grow higher compared to the results obtained for wide cracks. The cavity volume remained more or less unaffected due to increase in d_o. A few studies found that cavity diameter remains independent of the crack width when Q was sufficiently high to induce chimney fluidisation (Philippe and Badiane, 2013; Zoueshtiagh and Merlen, 2007), which is also evident form the cavity width results as shown in Figure 10(a), for d_o ranging from 1 to 5 mm, chimney regime was observed where the cavity width remained constant. As shown in Figure 10(c), the peak velocity decreased with increase in d_o. Figure 10(d) shows the effect of crack width on cavity shape when Q equals 12 L/min, for h_s as 150 mm. The cavity becomes wider with increase in crack width. Up to d_o of 3 mm not a significant change was observed; however, above d_o of 5 mm the significant change in height and width was observed.

Figures 11(a) and 11(b) show effect of friction angle on cavity hight using Schaeffer’s (1987) frictional viscosity model. The dimensions of the cavity decreased with the increase in the friction angle. There was little change in cavity height for ϕ greater than 40°. The increase in ϕ seems to affect the mechanical resistance and also strengthens the upper portion of the cavity (Philippe and Badiane, 2013). A dense system therefore has more resistance to movement and interparticle sliding. The increased resistance causes the jet to quickly lose energy, resulting in a shallow cavity height. However, higher friction angles can lead to more turbulent and disrupted flow, while lower friction angles allow for more stable and uniform jet penetration. Lower friction angles promote more uniform jet penetration and fluidisation. TFM-KTGF models dense systems using frictional stress. The results seem to agree with Tang et al.’s (2017) observation that the friction angle has little role on the onset of bed erosion. The critical velocity increased by 11% when the friction angle increased from 25° to 45° (Tang et al., 2017). The peak velocity of solid phase remained unaffected with change in ϕ as shown in Figure 10(b). As the value of ϕ was changed from 35° to 45°, sinϕ only increased from 0.57 to 0.70. As shown in Equation 10, the dominant factor is P_s and local dilation state due to which the velocity field hardly changed.

This study investigates the dynamics of upward water jets in sand beds. A computational fluid dynamics numerical tool was used to investigate the flow patterns and bed erosion processes for hidden cavity and burst. The numerical results show the capability of Eulerian coupled with kinetic theory of granular flows to approximately capture the phenomenon of soil fluidisation. The model was able to reproduce the phenomenology of soil fluidisation reported in several previous experimental observations. The limitation of the model did not permit to capture particle tracking and interparticle interaction. It was still able to quantify the cavity size, soil velocity vectors and volume fraction in the granular bed.

The simulation predicts well the changing trend of particle velocity along the bed axis, where the particles rapidly accelerated due to the high fluid velocity near the crack. This is followed by a constant velocity with a slow acceleration and then decelerate. The presence of bedding material surrounding a leaking pipe caused a significant change in the leakage rate.

Among the various parameters assessed, fluidisation seems to be greatly influenced by the particle size and crack width and less to the bed thickness and friction angle. Under the action of constant leakage rate by increasing the particle size from 0.1 to 2 mm, there was a significant 60%–90% reduction in the cavity height. Small crack produced high cavity; however, the volume of the cavity was found unaffected. The height of cavity reduces with increase in the crack width. For the leakage rate forming the chimney regime, the cavity width remained the same for different crack width. The bed thickness has minimal effect on the cavity dimensions. The effect of peak friction angle was incorporated via the closure of frictional viscosity which shows an effect on cavity dimensions as well as peak velocity of solid phase to be insignificant. It would therefore be prudent to use large granular particles and not small granular particles around underground pipes to control fluidisation.

This work was financially supported by PMRF and IIT Bombay. The software was provided by Computer Centre of the institute.