Keywords

1 Introduction

Dispersed two-phase flows occur in such areas of environmental modeling as the atmosphere (clouds, airborne particulate matter), sediment transport in rivers, and river flow with ice as well as in technological facilities (coolants in cooling systems, fuel combustion). The majority of the flows mentioned above are turbulent. In spite of geophysical flows involving sediment transport in rivers being among the first flows observed from the mechanical point of view, river flow with ice is much less researched than multiphase flows in technological facilities and flows with sediment transport [1, 2]. Nevertheless several hydrodynamic approaches to modeling river ice processes have been developed. One of the most important applications of these models is predicting riverside flooding caused by ice jams.

Static Models

Classical theoretical research on surface ice jams formed on a river during the breakup is based on the assumption that an ice jam is a one-dimensional static formation of ice floes with constant porosity. River flow velocity and characteristics of ice cover are considered to be constant across the river.

This approach allows evaluating the thickness of the existing ice jam but cannot indicate the initiation of an ice jam in a particular place and the conditions in which a jam could form, because they cannot correctly describe transport of ice and ice interaction with a flow. Static models are applicable when there is much data about the object of modeling such as ice thickness, density, friction between ice particles in the jam, and some other characteristics of ice cover. Wide review and comparison of the static models based on test cases is found in [1, 3].

Dynamic Models

Dynamic models of ice processes in rivers include a hydrodynamic model of the flow and a model of moving particles. Because multiphase flow is a very complex phenomenon, dynamic models began to develop from one-dimensional spatial approximation [4,5,6]. One-dimensional models, however, have limited application, because the transport of ice particles is a significantly two-dimensional process due to the wall friction, bathymetry, and mutual influence of the flow and the ice jam.

An advanced approach to modeling river flow with ice involves a multidimensional hydrodynamic model of the flow and a model to describe the movement of solid particles. Some of the multidimensional dynamic models that are precise in describing a flow with ice particles are that of Shen et al. [7,8,9], which is built within the finite-element Eulerian approach and Lagrangian model of ice particles; the three-dimensional unsteady Eulerian two-phase model of the flow with ice particles suggested by Wang et al. [10]; the model based on two-dimensional inviscid shallow water equations and Lagrangian model for ice particles by Shlychkov et al. [11]; and the two-dimensional hydrodynamic code with DEM for ice particles by Stockstill et al. [12].

Mathematical models that are part of much current hydrological software are quite robust in mathematical description of physical processes and require high-quality input data about a hydrological object, especially its bathymetry. Therefore developing mathematical models for hydrodynamic investigation of small rivers without precise bathymetric data about them is of particular interest.

Developing effective numerical algorithms for solving hydrodynamic equations are also of particular importance, especially in cases of river breakup, floods, and flows with ice particles, because of the complexity of these problems.

In this work, both a new mathematical model and a numerical method for computing two-phase flow of water with light particles densely placed at the water surface are presented. This approach meets all the requirements that exist for an advanced, up-to-date model.

2 Mathematical Model

Two-phase isothermal flow of the mixture of water and light particles in an open channel (river bed) is considered. Thermal exchange between phases is not considered because the temperatures of the water and the environment are close and change slightly during the period modeled. Because the considered density of ice \( \uprho_{i}^{0} = 910\,{\text{kg/m}}^{3} \) is less than the density of water \( \uprho_{l}^{0} = 1000\,{\text{kg/m}}^{3} , \) ice particles are considered to be packed in the upper layer of water and their concentration remains constant at the inlet of the channel (or the section of the river). Interactions between particles (collisions and friction) are also accounted for. Horizontal dimensions of the area of modeling are presumed to be much greater than the water depth. The size of ice particles is significantly less than the characteristic linear dimension of the channel (river bed).

Hydrostatic balance was applied because of the significant difference in the scale of the process in vertical and horizontal, and therefore all the terms in the equation for vertical velocity are negligible except ones that describe pressure and the force of gravity.

The mathematical model of the process described is based on continuum mechanics approach [13]. Light particles densely placed on the water surface are regarded as a continuum with effective properties.

Equations that describe the flow of the liquid phase (water) are [14]

$$ \frac{{\partial h^{\prime}}}{\partial t} + \nabla \cdot \left( {h^{\prime}{\vec{\mathbf{w}}}_{l} } \right) = 0; $$
(1)
$$ h^{\prime}\left( {\frac{{\partial {\vec{\mathbf{w}}}_{l} }}{\partial t} + \left( {{\vec{\mathbf{w}}}_{l} \cdot \nabla } \right){\vec{\mathbf{w}}}_{l} } \right) = - gh^{\prime}\uppsi \nabla \left( {z_{b} + h} \right) + \nabla \cdot h^{\prime}{\varvec{\uptau}}_{{\mathbf{l}}} \left( {\vec{w}_{l} } \right) - c_{f} \left| {{\vec{\mathbf{w}}}_{l} } \right|{\vec{\mathbf{w}}}_{l} + {\vec{\mathbf{F}}}_{{\mathbf{l}}} . $$
(2)

For the dispersed phase of light particles (ice)

$$ \frac{{\partial h^{\prime\prime}}}{\partial t} + \nabla \cdot \left( {h^{\prime\prime}{\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right) = 0; $$
(3)
$$ h^{\prime\prime}\left( {\frac{{\partial {\vec{\mathbf{w}}}_{{\mathbf{i}}} }}{\partial t} + \left( {{\vec{\mathbf{w}}}_{{\mathbf{i}}} \cdot \nabla } \right){\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right) = - \frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}gh^{\prime\prime}\uppsi \nabla \left( {z_{b} + h} \right) + \nabla \cdot h^{\prime\prime}{\varvec{\uptau}}_{{\mathbf{i}}} \left( {\vec{w}_{i} } \right) - c_{f}^{i} \left| {{\vec{\mathbf{w}}}_{i} } \right|{\vec{\mathbf{w}}}_{i} + {\vec{\mathbf{F}}}_{i} . $$
(4)

Here g is the gravitational acceleration; \( {\vec{\mathbf{w}}} = \left( {w_{1} ,\,w_{2} } \right) \) is the velocity of phase; \( \uprho_{l}^{0} \) is the density of water. \( z_{b} = z_{b} (x,\,y) \) is the function that describes the bathymetry; and \( h \) is the water depth. \( \int\limits_{{h - h_{i} + z_{b} }}^{{h + z_{b} }} {\upalpha_{i} dz = \bar{\upalpha }_{i} h_{i} = h^{\prime\prime}} , \) hʹ = h − hʺ where hi is the characteristic depth of the layer of ice particles. \( c_{f} = \frac{{gn^{2} }}{{h^{0.333} }} \) is the bed friction of the liquid phase, and n > 0 is the Manning coefficient.

The bed friction of the dispersed phase in shallow regions is computed as \( c_{f}^{i} \left| {{\vec{\mathbf{w}}}_{i} } \right|{\vec{\mathbf{w}}}_{i} \), where \( c_{f}^{i} = 0.0025\left| {{\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right|^{ - 1} \) [15]. \( \uppsi = \tilde{\upalpha }_{l} + \left( {\frac{{\uprho_{i}^{0} }}{{\uprho_{l}^{0} }} - 1} \right)\left( {1 - \tilde{\upalpha }_{l} } \right) \), where \( \tilde{\upalpha }_{l} \) is the depth-averaged volume fraction of water (\( 0 < \tilde{\upalpha }_{l} \le 1 \)).

$$ \uptau_{kj} (\vec{w}) = \left( {\upnu_{{}}^{0} + \upnu_{{}}^{t} } \right)\left[ {\left( {\frac{{\partial w_{k} }}{{\partial x_{j} }} + \frac{{\partial w_{j} }}{{\partial x_{k} }}} \right) - \frac{2}{3}\updelta_{kj} div\,\vec{w}_{{}} } \right] - \frac{2}{3}\updelta_{kj} k,(j,k = 1,2) $$

are the components of the viscous and turbulent stress tensor for the liquid/ice particles; \( \upnu_{l,i}^{0} \) is the viscosity of liquid/dispersed phase; \( \upnu_{l,i}^{t} \) is the eddy viscosity of the flow (liquid/ice); \( \updelta_{kj} \) is the Kronecker delta; and k is the turbulent kinetic energy. Viscous stress tensor in the dispersed phase appears due to collisions of particles and friction between them.

2.1 Force Terms in the Momentum Equations for Phases

The force term in the momentum equation for the dispersed phase is \( {\vec{\mathbf{F}}}_{{\mathbf{i}}} = {\vec{\mathbf{F}}}_{{\mathbf{A}}} + {\vec{\mathbf{F}}}_{{\varvec{\upmu}}} + {\vec{\mathbf{F}}}_{{{\mathbf{VM}}}} + {\vec{\mathbf{F}}}_{{\mathbf{C}}} \), which is the sum of the following forces [13]:

  • buoyancy

    $$ {\vec{\mathbf{F}}}_{{\mathbf{A}}} = h^{\prime\prime}\frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\left( {\frac{{D_{l} {\vec{\mathbf{w}}}_{{\mathbf{l}}} }}{{D_{l} t}} - {\vec{\mathbf{g}}}} \right);\,\frac{{D_{l} {\vec{\mathbf{w}}}_{{\mathbf{l}}} }}{Dt} = \frac{{\partial {\vec{\mathbf{w}}}_{{\mathbf{l}}} }}{\partial t} + \left( {{\vec{\mathbf{w}}}_{{\mathbf{l}}} \cdot \nabla } \right){\vec{\mathbf{w}}}_{{\mathbf{l}}} ; $$

  • fluid drag force [16]

    $$ {\vec{\mathbf{F}}}_{{\varvec{\upmu}}} = \left\{ \begin{aligned}& \left[ {150\frac{{h^{\prime\prime}\tilde{\upalpha }_{i} \left( {1 - \tilde{\upalpha }_{l} } \right)\upmu_{l} }}{{\tilde{\upalpha }_{l} \uprho_{i}^{0} \left( {d_{i} f_{i} } \right)^{2} }} + 1.75\frac{{h^{\prime\prime}\uprho_{l}^{0} \left| {{\vec{\mathbf{w}}}_{{\mathbf{l}}} - {\vec{\mathbf{w}}}_{i} } \right|}}{{\uprho_{i}^{0} d_{i} f_{i} }}} \right]\left( {{\vec{\mathbf{w}}}_{{\mathbf{l}}} - {\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right),\,\,\tilde{\upalpha }_{l} \le 0.8; \\ &\frac{3}{4}\frac{{h^{\prime\prime}\uprho_{l}^{0} c_{D} \tilde{\upalpha }_{l}^{ - 1.65} }}{{\uprho_{i}^{0} d_{i} f_{i} }}\left| {{\vec{\mathbf{w}}}_{{\mathbf{l}}} - {\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right|\left( {{\vec{\mathbf{w}}}_{{\mathbf{l}}} - {\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right),\,\,\tilde{\upalpha }_{l} > 0.8; \hfill \\ \end{aligned} \right. $$

    where \( c_{D} = \hbox{max} \left[ {\frac{24}{{\text{Re}_{p} }}\left( {1 + 0.15R_{p}^{0.687} } \right),\,0.44} \right] \) is the dimensionless drag coefficient, and di is the characteristic diameter of the ice particles; \( \text{Re}_{p} = \frac{{\left| {{\vec{\mathbf{w}}}_{{\mathbf{l}}} - {\vec{\mathbf{w}}}_{{\mathbf{i}}} } \right|d_{i} \uprho_{l}^{0} \upalpha_{l} }}{{\upmu_{l}^{0} }}; \)

  • virtual mass

    \( {\vec{\mathbf{F}}}_{{{\mathbf{VM}}}} = c_{VM} h^{\prime\prime}\frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\left( {\frac{{D_{l} {\vec{\mathbf{w}}}_{{\mathbf{l}}} }}{Dt} - \frac{{D_{i} {\vec{\mathbf{w}}}_{{\mathbf{i}}} }}{Dt}} \right) \), CVM ≈ 0.5 for spherical particles;

  • coriolis force \( {\vec{\mathbf{F}}}_{{\mathbf{C}}} = h^{\prime\prime}\left( {{\vec{\mathbf{w}}}_{{\mathbf{i}}} \times {\vec{\varvec{\upomega }}}} \right) \).

    The model is closed with the high-Reynolds k − ε turbulence model for depth-averaged equations [17] with a modification suggested by Pourahmadi and Humphrey [18] to account for the influence of particles on the flow. Constants of the model are the same as in the standard high-Reynolds \( \bar{k} - \bar{\upvarepsilon } \) model [17].

Initial and Boundary Conditions

At the initial moment t = 0 the following conditions are used:

  • hʺ = hice, hice is a known value;

  • \( u_{i} = u_{i0} ,\,\,v_{i} = v_{i0} ;\,h^{\prime} = h - h_{ice} ;\,u_{l} = u_{l0} ,\,\,v_{l} = v_{l0} . \)

At the inlet of the channel, parameters of the liquid phase are considered to be known; at the outlet normal derivatives of the parameters are set to zero. At the solid boundaries shear stresses are considered both for liquid and for particles. Both normal and tangential velocities on the wall are set to zero. In the liquid phase, near the wall the stresses and turbulent characteristics of the liquid phase were set by Launder-Spalding wall functions [19].

3 Numerical Method

The area of the flow is contained in the rectangle covered with structured mesh. The equations of the model are discretized with the finite volume method on staggered mesh (Fig. 1), that is, finite volumes for velocity components are half-cell shifted from the node P, which is the center of the cell where scalar values \( h^{\prime},\,\,h^{\prime\prime},\,\,\bar{k},\,\,\bar{\upvarepsilon } \) are defined.

Fig. 1.
figure 1

Mesh stencil. Uppercase letters are for centers of the finite volumes, lowercase letters are for midpoints of their edges.

Full size image

All the terms in the Eqs. (1)–(4) were approximated explicitly in time except for the drag force in the momentum Eqs. (2) and (4) which was approximated implicitly.

Convective fluxes in the momentum equations and advection-diffusion equations are approximated the MUSCL scheme [20].

To reach the third order accuracy in space in regions where the functions are monotonic, variables on the edges of the finite volumes are computed with linear interpolation from the values in the centers of the volumes (Fig. 1) [21]:

$$ \begin{aligned} \Phi_{e} = \Phi_{P} + \frac{\Delta x}{2}\upsigma_{e}^{ + } ,\,\,\bar{u}_{e} > 0; \hfill \\ \Phi_{e} = \Phi_{E} - \frac{\Delta x}{2}\upsigma_{e}^{ - } ,\,\,\bar{u}_{e} \le 0; \hfill \\ \end{aligned} $$

where \( \upsigma_{e}^{ + } = \frac{{\updelta_{e} }}{\Delta x}\upphi \left( {\uptheta_{e} } \right),\,\,\upsigma_{e}^{ - } = \frac{{\updelta_{e} }}{\Delta x}\upphi \left( {\uptheta_{ee}^{ - 1} } \right) \) are the slopes, limited by the function \( \upphi \left( {\uptheta_{e} } \right);\,\,\updelta_{e} = \Phi_{E} - \Phi_{P} \). Slope limiter \( \upphi \left( \uptheta \right) \) was chosen to satisfy the sufficient condition of the Harten’s theorem [22] \( \left( {0 \le \frac{{\upphi \left( {\uptheta_{e} } \right)}}{{\updelta_{e} }} \le 2,\,\,0 \le \upphi \left( {\uptheta_{e} } \right) \le 2} \right) \), which means that reconstructed value is within the values used for reconstruction. In this work the function \( \upphi \left( \uptheta \right) = \hbox{max} \left[ {0,\,\hbox{min} \left( {2\uptheta ,\frac{2 + \uptheta }{3},2} \right)} \right] \) [20] was used as a limiter. Here \( \,\uptheta_{e} = \frac{{\updelta_{e} }}{{\updelta_{w} }}\,\left( {\updelta_{w} \ne 0} \right), \) and \( \updelta_{e} \) could be approximated with left, right, or central finite difference.

First order upwind scheme was used to approximate convective terms in the equations of the turbulence model.

Central finite difference scheme

$$ \left[ {gh\frac{{\partial (z_{b} + h)}}{\partial x}} \right]_{e} \approx \left\{ {\begin{array}{*{20}l} {gh_{e} \frac{{h_{E} - h_{P} + z_{bE} - z_{bP} }}{\Delta x},\,\,h_{E} \,{\text{and}}\,h_{P} > \upvarepsilon_{wd} > 0\,({\text{if the cells with}}} \hfill \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {\text{centers}}\,{\text{in E and P are wet)}},} \hfill \\ {0,\,\,h_{E} \,{\text{or}}\,h_{P} < \upvarepsilon_{wd} \, ( {\text{if one of the cells is dry);}}} \hfill \\ {\upvarepsilon_{wd} \,{\text{is the small positive value}}} \hfill \\ \end{array} } \right. $$

was used to approximate the source terms, which represent the influence of the bed slope in momentum equations.

The diffusive terms are discretized with the second order central-difference scheme.

In order to reduce a limitation on the time step for momentum equations, terms that express the dynamic interaction between phases (friction) was approximated with an implicit scheme and the partial elimination algorithm that allows solving equations in the areas with no dispersed phase (hʺ = 0). A brief description of this approach is below.

Consider Eqs. (2) and (4) for the longitudinal velocity component for the inner node e of the mesh shown at Fig. 1.

$$ \frac{{h_{e}^{'0} u_{le} - h_{e}^{'0} u_{le}^{0} }}{\Delta t}\Delta x\Delta y = \Phi_{e}^{0} + \upbeta^{0} h_{e}^{''\,0} \left( {u_{ie} - u_{le} } \right); $$
(5)
$$ \frac{{h_{e}^{''0} u_{i\,e} - h_{e}^{''0} u_{i\,e}^{0} }}{\Delta t}\Delta x\Delta y = \Psi_{e}^{0} + \frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\upbeta^{0} h_{e}^{''0} \left( {u_{l\,e} - u_{i\,e} } \right). $$
(6)

\( \Phi_{e}^{0} \,\,{\text{and}}\,\,\Psi_{e}^{0} \) combine approximations of convective, diffusive, and source terms, β0 is the coefficient of difference between velocities of the phases in the expression for the drag force. Upper index ‘0’ is for the discrete values from the previous time step. Here and below in the section the over-bar for averaging is omitted.

In the right side of the discrete Eqs. (5)–(6), terms \( \upbeta^{0} h^{\prime\prime}_{e} \left( {u_{i\,e} - u_{l\,e} } \right) \) and \( \frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\upbeta^{0} h_{e}^{''0} \left( {u_{le} - u_{ie} } \right) \) describe dynamic interactions between phases and contain the differences between velocities of phases. Explicit approximation of these terms leads to a more strict convergence condition than the Courant–Friedrichs–Lewy condition

$$ \uptau < \frac{0.5\,\Delta x\,\Delta y}{{\hbox{max} \left| {u_{e} } \right|\Delta y + \hbox{max} \left| {v_{n} } \right|\Delta x + \sqrt {gh} \left( {\Delta x + \Delta y} \right)}} $$

in the case of flow with particles with large inertia.

In the case of hʺ = 0 Eq. (6) becomes an identity. To avoid this, rewrite Eqs. (5)–(6):

$$ \left( {h_{e}^{'0} \Delta x\Delta y + \Delta t\upbeta^{0} h_{e}^{''0} } \right)u_{le} - \Delta t\upbeta^{0} h_{e}^{''0} u_{ie} = \Delta t\Phi_{e}^{0} + h_{e}^{'0} u_{le}^{0} \Delta x\Delta y; $$
(7)
$$ \left( {\Delta x\Delta y + \frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\Delta t\upbeta } \right)u_{ie} - \frac{{\uprho_{l}^{0} }}{{\uprho_{i}^{0} }}\Delta t\upbeta u_{le} = {{\Delta t\Psi_{e}^{0} } \mathord{\left/ {\vphantom {{\Delta t\Psi_{e}^{0} } {h_{e}^{''0} }}} \right. \kern-0pt} {h_{e}^{''0} }} + u_{ie}^{0} \Delta x\Delta y. $$
(8)

In (8), the term \( {\raise0.7ex\hbox{${\Psi_{e}^{0} }$} \!\mathord{\left/ {\vphantom {{\Psi_{e}^{0} } {h_{e}^{''0} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${h_{e}^{''0} }$}} \) \( \Psi_{e}^{0} \) is equated to zero if \( h_{e}^{''0} < \upvarepsilon \), where ε is an infinitely small positive value. The determinant of the SLAE (7)–(8) is nonzero and the system is being solved by Cramer’s rule. The SLAE for vl n, vi n was solved in the same way.

The wet-dry boundary treatment for an unsteady flow is of particular complexity because the solution becomes unstable due to very little water depth in the boundary cell [23, 24]. One of the simplest approaches to wet-dry boundary treatment is choosing a small positive value ε > 0 such that if the depth becomes less than ε, the cell is considered dry and excluded from computations. In the case of the two-phase flow considered in this article, hʹ, which is the depth of the liquid phase, was compared to ε. When the cell became dry, the depth of the dispersed phase hʺ was also equated to zero.

Several test cases of unsteady flow in open channels were computed with the model and the method presented above. The results were compared to experimental data to evaluate the model.

4 Results and Discussion

4.1 Two-Phase Flow in a 180° Bend Flume

Unsteady turbulent flow in a 180° bend flume with polypropylene particles to model ice was computed. This test case was investigated both experimentally and numerically in [25]. The depth of the flow at the inlet of the channel was set to 0.45 m, and the volumetric flow rate was 0.16 m3/s. These parameters define turbulent flow with Reynolds number Reh = 77170 defined by the depth of the flow. Spherical polypropylene particles (ρ = 900 kg/m3) with the diameter di = 0.005 m were used to model ice floes. Particles were tossed with constant rate into the steady flow at the end of the first straight section of the channel. A wire-mesh screen was placed after the bend to initiate ice jam development and support its downstream end.

Results of the computations made were compared with those of Urroz and Ettema [25] (Fig. 2).

Fig. 2.
figure 2

Jam-head profiles around the bend (a) experiment [25], (b) computations

Full size image

Moving along the flume, particles tend to accumulate near the left wall because of the centrifugal force, which also causes the increase in the velocity towards the left wall. After crossing the longitudinal axis of the channel, the flow tends to move closer to the right bank (Fig. 2).

This showed that the mathematical model and the computational method proposed accurately predicted both the velocity field and the distribution of the particles in the channel: they showed the increase in ice jam thickness both downstream and toward the inner bank of the bend as it was observed in experiments and gave a jam-head shape similar to those observed.

In [25] it was also noted that a similar shape of the front of ice particles was observed on the Iowa River during the breakup.

4.2 Two-Phase Flow in Open Channel with a 90° Bend

Two-phase flow in the channel with a 90° bend was computed. The first section of the channel was 5.555 m long and 0.86 m wide with a flat bottom, and the second section was 4.43 m long and 0.72 wide with a flat bottom. At the end of the first section there was a step with a change in the bed elevation of 0.013 m. One-phase flow in this channel was investigated in detail in [23, 26]. In the present article, the two-phase flow of the water with ice particles in the described channel is discussed. At the inlet, the longitudinal velocity of the flow was U0 = 0.2 m/s, the initial depth of the flow was h = 0.175 m, and the depth of the dispersed phase was hʺ = 0.04 m. Parameters of the dispersed phase are shown in Table 1. Structured mesh of 254 × 208 nodes was used in computations. Time step was chosen automatically from the Courant–Friedrichs–Lewy condition.

Table 1. Parameters of the test cases: 1 basic case, 2–3 parametric computations
Full size table

Computations led to the following conclusions:

  1. 1.

    Velocities of the smaller and the larger particles did not differ significantly near the bend. A recirculation area formed behind the bend near the right wall, and recirculating flow was more intensive for the larger particles (di = 0.1 m).

    In the corner of the channel, velocities of the phases were also different, but for cases (1) and (2) the particles did not accumulate in the corner because the particles had little inertia and wall friction to resist the dragging force of the water. It was also found that the flow with smaller particles had less turbulent kinetic energy behind the bend.

  2. 2.

    The change in bottom elevation influenced the velocities and the depth of the dispersed phase layer more than the size of the particles. Contours of the velocity magnitudes of both phases were not smooth near the step on the bottom because particles with little inertia in the same way as the liquid phase reacted to the change in the flow conditions (Fig. 3). At the same time, the distribution of hʺshowed that with the increase in size of particles (and consequently their inertia) hʺ increases more smoothly near the step and again near the bend.

    Fig. 3.
    figure 3

    Liquid phase velocity field (a), dispersed phase velocity field (b), depth of the dispersed phase layer (c) for the basic variant (1), and computation with smaller particles (2)

    Full size image

Evaluation of the Influence of the Shape of Particles

With decrease in fi from fi = 1 to fi = 0.16, resistance of the dispersed phase to the dragging force of the liquid phase increased (Fig. 4). The increase in the depth hʺ in the corner of the channel was more efficient for spherical particles. In the case of the flow with spherical particles, the intensity of the recirculating flow was greater and the contours of hʺ near the step were smoother. Computed distribution of the free surface in the bend was in accordance with experimental data [27].

Fig. 4.
figure 4

Liquid phase velocity field (a), dispersed phase velocity field (b), and the depth of the dispersed phase (c) for the basic variant (1) and computation with spherical particles (3)

Full size image

4.3 Modeling the Floods of the Tom River Near Tomsk, Russia

The abrupt increase in the flow rate of the river at the inlet of the section studied was computed. Such a phenomenon occurs when a sudden failure of an ice jam upstream releases large quantities of water and ice and often causes severe damage in the floodplain. As a precondition, flow with a longitudinal velocity of 1.25 m/s and free surface elevation of 69.8 m above sea level was computed for 20 h. After flow stabilization, the depth at the inlet of the section studied was increased by 2 m. The structured mesh of 221 × 180 nodes was used in computations.

With the increase in volume of fluid moving from the inlet, the depth of the flow in the section studied increased and water flooded the areas of the floodplain with the lowest heights: islands near 12,000 m and lowlands on the left bank of the river between 18,000 and 20,000 m (Fig. 5). The depth of the lake on the left bank of the river (near 6,000) increased and the lake connected to the river. Flooding of these areas is usually observed during the spring breakup.

Fig. 5.
figure 5

Free surface level before increasing water level, 9.61 h after increasing water level, and 16.181 h after increasing the water level

Full size image

5 Conclusions

A new mathematical model that is based on mechanics of interacting, interpenetrating continua was developed and tested. It is one of the first Eulerian hydrodynamic models of the two-phase flow of water with ice particles within the shallow water approach that takes into account interaction between phases, interactions between ice particles, interaction of both phases with the river bed, and the turbulence of the flow.

A computational algorithm based on finite volume method, semi-implicit discretization in time, and original partial elimination technique was developed to avoid uncertainty in the areas with no dispersed phase. The novelty of the computational method is in combining the technique for detecting the moving boundary of the river, which is of particular importance in modeling floods, with the partial elimination technique that ensures the correctness of computations throughout the domain in cases of areas with no ice particles.

The approach presented in this work was applied to numerical investigation of the Tom River during the breakup in spring near the city of Tomsk (Russia), where ice jams often form and cause localized flooding in the area behind the blockage. Two-phase river flow modeling showed that the approach developed correctly modeled changes in hydrodynamic characteristics of the river during the breakup and gave accurate results for the cases of flooding of the floodplain that involve change in the boundary of the river.

The investigation of the two-phase flow in the channel with 90° bend showed that the dispersed phase is more influenced by the bottom relief than by sharp change in flow direction. It was also shown that the presence of the ice particles in the flow increased the nonuniformity in the free surface level.