**4. RNG based k-ε erosion model equations for the HET**

The RNG *k* � *ε* model differs from the standard model by the special form of the transport equations, which contain the additional term *Rε*. These equations write

$$\frac{\partial k}{\partial t} + \frac{1}{r} \frac{\partial (r k u)}{\partial r} + \frac{\partial (k v)}{\partial \mathbf{z}} = \frac{1}{r} \frac{\partial}{\partial r} \left( a \mu\_t r \frac{\partial k}{\partial r} \right) + \frac{\partial}{\partial \mathbf{z}} \left( a \mu\_t \frac{\partial k}{\partial \mathbf{z}} \right) + a \mu\_t \frac{k}{r^2} + \mu\_t \mathbf{S}^2 - \varepsilon \tag{1}$$

$$\frac{\partial \boldsymbol{\varrho}}{\partial t} + \frac{\mathbf{1}}{r} \frac{\partial (\boldsymbol{r} \boldsymbol{u} \boldsymbol{u})}{\partial r} + \frac{\partial (\boldsymbol{\varepsilon} \boldsymbol{v})}{\partial \mathbf{z}} = \frac{\mathbf{1}}{r} \frac{\partial}{\partial r} \left( a \mu\_t r \frac{\partial \boldsymbol{\varepsilon}}{\partial r} \right) + \frac{\partial}{\partial \mathbf{z}} \left( a \mu\_t \frac{\partial \boldsymbol{\varepsilon}}{\partial \mathbf{z}} \right) + a \mu\_t \frac{\boldsymbol{k}}{r^2} + \mathbf{C}\_{1\varepsilon} \frac{\boldsymbol{\varepsilon}}{\rho \boldsymbol{k}} - \mathbf{C}\_{2\varepsilon} \frac{\boldsymbol{\varepsilon}^2}{\rho \boldsymbol{k}} - \frac{\mathbf{1}}{\rho} \mathbf{R}\_{\varepsilon} \tag{2}$$

With

$$R\_{\varepsilon} = \frac{C\_{\mu}\rho S^3 k^2 \varepsilon (\eta\_0 \varepsilon - \text{S}k)}{\eta\_0 \left(\varepsilon^3 + \beta \text{S}^2 k^3\right)}\tag{3}$$

$$\mathbf{S}^2 = 2\left(\frac{\partial \mathbf{u}}{\partial r}\right)^2 + 2\left(\frac{\partial \mathbf{v}}{\partial \mathbf{z}}\right)^2 + \left(\frac{\partial \mathbf{v}}{\partial r} + \frac{\partial \mathbf{u}}{\partial \mathbf{z}}\right)^2 \tag{4}$$

where u and v are average radial and axial flow velocities, r and z are axial and radial coordinates, t is time, *ε* the rate of dissipation of turbulent kinetic energy, *k* is the turbulent kinetic energy, *ρ* is the fluid density, *μ<sup>t</sup>* the total kinematics viscosity, *α* is the inverse effective Prandtl number for both *k* and *ε*, and *C*1*<sup>ε</sup>*, *C*2*<sup>ε</sup>*, *Cμ*, *η*<sup>0</sup> and *β* are constants.

To calculate the tangential shear stress distribution on the cylinder inner wall, boundary conditions must be introduced. When the shear stress is calculated using ANSYS Fluent software, the classical linear erosion law is used to estimate the erosion rate. This law states the estimated erosion rate, defined as the mass departure of particles due to erosion per unit time and per unit area, is expressed by the following formula: *ε*\_*er* ¼ *cer*ð Þ *τ* � *τcr* where *cer* and *τcr* are constant depending on the considering soil material. The rate *ε*\_*er* can be related to time variation of local radius by *ε*\_*er* ¼ *ρddR=dt* where ρ*<sup>d</sup>* represents the dry density of the soil and *R* the inner radius of the hole. The law of erosion states that the rate of erosion ε*er* is proportional to the shear stress, exceeding the critical shear *τcr* for which erosion begins [14].
