3.1.3. Solid Phase Mathematical Model

Solid particles were simulated using a Lagrangian approach and were treated as if their volume fraction were low compared to that of the continuous phase. Equation (15) was derived from the force balance on the Lagrangian reference frame.

$$\frac{\partial (v\_{\mathcal{P}})\_{i}}{\partial t} = \frac{18\mu}{\rho\_{\mathcal{P}}d\_{\mathcal{P}}^{2}} \frac{\mathbb{C}\_{D}\mathrm{Re}}{24} + \left(\frac{\rho}{\rho\_{\mathcal{P}}}\right)(v\_{\mathcal{P}})\_{i}\frac{\partial u\_{i}}{\partial x\_{i}} + \left(1 + \frac{\rho}{\rho\_{\mathcal{P}}}\right)g - \frac{1}{2}\frac{\rho}{\rho\_{\mathcal{P}}}\frac{\partial}{\partial t}(u\_{i} - v\_{\mathcal{P}i}) + F\_{z} \tag{15}$$

where:


The term *Fz* activates force terms in situations in which multiple reference frames are used and frame or mesh rotation is activated. The drag coefficient was estimated using the equations for the nonspherical drag law shown in Equation (16).

$$C\_D = \frac{24}{Re\_{sph}}(1 + b\_1 Re\_{sph}^{b\_2}) + \frac{b\_3 Re\_{sph}}{b\_4 + Re\_{sph}}\tag{16}$$

$$b\_1 = \exp(2.3288 + 6.4581\phi + 2.4486\phi^2) \tag{17}$$

$$b\_2 = 0.0964 + 0.5565\phi \tag{18}$$

$$b\_3 = \exp(4.905 - 13.8944\phi + 18.4222\phi^2 - 10.2599\phi^3) \tag{19}$$

$$b\_4 = \exp(1.4681 + 12.2584\phi - 20.7322\phi^2 + 15.8855\phi^3) \tag{20}$$

where:


Particle dispersion caused by turbulent flows can be estimated using the stochastic tracking model [24,25]. This method, also known as the discrete random walk model, takes into account the effect of turbulent velocity fluctuations on the trajectories of particles. Instantaneous fluid velocity, as shown in Equation (21), was used to integrate the particle trajectory equations along their path to predict the turbulent dispersion of particles. Random velocity fluctuations *u* were determined through Equation (22), where *ζ* is a normally distributed random number. Particle diffusivity was estimated using Equation (23), where the integral time scale as defined in Equation (24) describes the time the particle remains in turbulent motion along a path *ds*.

$$
\mu = \overline{\mu} + \mu'(t) \tag{21}
$$

$$
\mu' = \zeta \sqrt{2k/3} \tag{22}
$$

$$D(t) = \overline{u\_i' u\_j'} T \tag{23}$$

$$T = \int\_{o}^{\infty} \frac{v\_p'(t)v\_p'(t+s)}{\overline{v\_p'^2}} ds \tag{24}$$

3.1.4. Erosion Model

The erosion model developed by Oka, Okamura, and Yoshida [26,27] was used to determine the erosion for the present case. This model is one of the most frequently used to determine the erosion in CFD analyses where solid particles suspended in a liquid medium are present [28–30]. The equation developed by Oka is the following:

$$E(\mathfrak{a}) = \mathfrak{g}(\mathfrak{a})E\_{\mathfrak{A}0} \tag{25}$$

where:


and:

$$\mathbf{g}(\mathfrak{a}) = (\sin \mathfrak{a})^{n1} (1 + Hv(1 - \sin \mathfrak{a}))^{n2} \tag{26}$$

$$E\_{90} = K(aHv)^{k\_1 b} \left(\frac{v}{v'}\right)^{k\_2} \left(\frac{D}{D'}\right)^{k\_3} \tag{27}$$

$$k\_2 = 2.3(Hv)^{0.038} \tag{28}$$

$$2.1, n2 = 2.3(Hv)^{0.038} \tag{29}$$

where *Hv* is the material Vickers hardness in GPa, *k*<sup>2</sup> is a velocity exponent, *k*<sup>3</sup> is a diameter exponent, and constants *n*1 and *n*2 are model exponents used to calculate the impact angle influence on the erosion rate. *D* and *v* are the reference diameter and velocity, respectively. The calibrated values of these parameters for the present case, in which sand particles and stainless steel were considered, are the ones presented in Table 1.

**Table 1.** Oka model parameters.


#### *3.2. Geometry and Conditions*

Details of the general specifications of the turbine are presented in Table 2. In order to reduce the computational effort while using a more precise mesh for the numerical analysis, only one period of the turbine was simulated, which is comprises one stay vane, one guide vane, one runner blade, and an outlet domain representing the draft tube.

**Table 2.** Turbine specifications.


The computational domain is presented in Figure 3. The geometry was obtained performing a 3D scanning of the turbine. The obtained profiles were reconstructed using ANSYS BladeGen to obtain a smoother geometry optimized for mesh construction using TurboGrid.

#### *3.3. Operating Conditions*

Table 3 shows the operating conditions that were used to perform the simulations. These conditions were translated to mass flow inlet and pressure outlet boundary conditions. The atmospheric pressure from the region was used to define the outlet pressure.



**Figure 3.** Computational domain of San Francisco's Francis turbine.

The characteristics of the sediments used for this study are shown in Table 4. Particle mass flow rate was calculated as a function of the volumetric water flow rate of each case using the average particle concentration.

**Table 4.** Sediment characteristics.


#### *3.4. Mesh*

The mesh was generated using the Turbogrid module, which employs a high-fidelity hexahedric structured mesh with a uniform distribution. *y*<sup>+</sup> values were calculated for stay vanes, guide vanes, and runner blades through the following equation:

$$y^{+} = \frac{u\_{\!\!\! }y}{\nu} \tag{30}$$

where *uτ* is the friction velocity, *y* is the distance to the nearest wall, and *ν* is the kinematic viscosity. The first cell height of each domain was calculated to correctly compute the boundary layer with the selected turbulence model, as shown in Figure 3, by applying a locally refined region near the domain walls. The obtained *y*<sup>+</sup> values ranged between 45 and 125, which are considered appropriate for the selected turbulence model [31]. The quality of the mesh was also evaluated using the orthogonal quality model. The orthogonal quality of a cell was estimated as follows.

$$\min\left(\frac{\vec{A}\_i \vec{C}\_i}{|A\_i||\mathcal{C}\_i|}, \frac{\vec{A}\_i \vec{f}\_i}{|A\_i||f\_i|}\right) \tag{31}$$

where:


• *Ci* = vector from the centroid of the cell to the centroid of the adjacent cell.

A minimum orthogonal quality of 0.269 was obtained in 0.0000002% of mesh elements. The average orthogonal quality in all meshes was 0.93. The walls of the fully structured mesh are shown in Figure 4

Additionally, mesh independence studies were performed considering the pressure drop between the domain inlet and outlet. The structured mesh distribution was modified in all domain directions. Three different mesh resolutions were analyzed for each of the four individual operating conditions evaluated in the study. The results shown in Figure 5 indicate that the difference between the calculations of the fine and medium mesh resolutions is negligible.

**Figure 5.** Mesh independence analysis.
