key: cord-0045739-m7iegp6k authors: Churuksaeva, Vladislava; Starchenko, Alexander title: Numerical Modeling of the Two-Phase Flow of Water with Ice in the Tom River date: 2020-06-15 journal: Computational Science - ICCS 2020 DOI: 10.1007/978-3-030-50417-5_16 sha: 161b63e3b74d2f4aef031fee11530d818048df6c doc_id: 45739 cord_uid: m7iegp6k A new mathematical model and a numerical method were constructed for numerical investigation of a two-phase turbulent flow in an open channel. Solid particles with a density close to that of water were considered a continuous phase with effective properties. This new model is based on a continuum-mechanics approach, a hydrostatic assumption, and equations averaged by the flow depth. Turbulent closure of the equations was done with a two-parameter k − ε turbulence model modified by Pourahmadi and Humphrey to account for the influence of the particles on the turbulent structure of the flow. The new numerical method is based on partial elimination algorithm for computing areas of the two-phase flow free of ice particles and uses semi-implicit approximation in time. The influence of the dynamic parameters of the dispersed phase on the structure of the flow was also investigated by computing several scenarios of the flow in an open channel with a 90-degree bend. Applications of the approach to the modeling of riverside flooding due to sudden increase in the river depth after a release of an ice jam illustrate the capabilities of the model. 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. 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 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 finiteelement Eulerian approach and Lagrangian model of ice particles; the threedimensional 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 twodimensional 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-todate 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 q 0 i ¼ 910 kg/m 3 is less than the density of water q 0 l ¼ 1000 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] @h 0 @t þ r Á h 0w For the dispersed phase of light particles (ice) h 00 @w i @t Here g is the gravitational acceleration;w ¼ w 1 ; w 2 ð Þis the velocity of phase; q 0 l is the density of water. z b ¼ z b ðx; yÞ is the function that describes the bathymetry; and h is the of the layer of ice particles. c f ¼ 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 i Þ, whereã l is the depthaveraged volume fraction of water (0\ã l 1). are the components of the viscous and turbulent stress tensor for the liquid/ice particles; m 0 l;i is the viscosity of liquid/dispersed phase; m t l;i is the eddy viscosity of the flow (liquid/ice); d 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. The force term in the momentum equation for the dispersed phase is F i ¼F A þF l þF VM þF C , which is the sum of the following forces [13] : where c D ¼ max 24 Re p 1 þ 0:15R 0:687 and d i is the characteristic diameter of the ice particles; Re p ¼w l Àw i The model is closed with the high-Reynolds k − e turbulence model for depthaveraged 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 k À e model [17] . At the initial moment t = 0 the following conditions are used: 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] . 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 0 ; h 00 ; k; e are defined. 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]: 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) U 0 e and W 0 e combine approximations of convective, diffusive, and source terms, b 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. Þ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 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 e > 0 such that if the depth becomes less than e, 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 e. 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. 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 m 3 /s. These parameters define turbulent flow with Reynolds number Re h = 77170 defined by the depth of the flow. Spherical polypropylene particles (q = 900 kg/m 3 ) with the diameter d i = 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) . 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. 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 U 0 = 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. Computations led to the following conclusions: 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 (d i = 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. 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. With decrease in f i from f i = 1 to f i = 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] . 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. 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. Twophase 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. Comparative testing of numerical models of river ice jams Progress in studies on ice accumulation in river bends Rivers ice jams: theory, case studies, and applications Wave-propagation in ice-covered channels Dynamic transport of river ice Numerical modeling of ice jams SPH simulation of river ice dynamics Shokotsu river ice jam formation Mathematical modeling of river ice processes Two-phase flow numerical simulation of a bend-type ice sluice in the diversion water channel of powerhouse Hydrodynamic model of a river breakup for studying ice jams Modeling floating objects at river structures Dynamics of Multiphase Media Numerical modeling of the two-phase flow of a liquid with particles in an open channel. Tomsk State Univ Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions Numerical investigation of the aerodynamics in the circulating fluidized bed installations A depth-averaged mathematical model for the near field of side discharges into open channel flow Modelling solid-fluid turbulent flows with application to predicting erosive wear The numerical computation of turbulent flows Towards the ultimate conservative difference scheme. Part II. Monotonicity and conservation combined in a second order scheme High resolution schemes for hyperbolic conservation laws Depth averaged modelling of turbulent shallow water flow with wet-dry fronts A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography Bend ice jams: laboratory observations Mathematical modeling of a river stream based on a shallow water approach Characteristics of flow around open channel 90°b ends with vanes Acknowledgments. The numerical method was developed and the computations were performed by V. Churuksaeva, and this work was funded by Russian Foundation for Basic Research under research project N 18-31-00386.The problem was stated and the mathematical model was formulated by A. Starchenko, and this work was supported by the Ministry of Science and Higher Education of the Russian Federation under project of Regional Scientific and Educational Mathematical Center of Tomsk State University.The authors wish to thank Jean Kollantai of Tomsk State University for her assistance with the style of the paper.