Numerical Study on Flow-Melt Characteristics of Ice Slurry in Horizontal Straight Pipe with a Local Large Heat Flux Segment

Abstract

Ice slurry is a high energy density coolant with excellent flow, phase change, and thermophysical properties. In order to investigate the recovery of ice slurry flowing through a local large heat flux segment, a 3D Eulerian-Eulerian model based on the granular kinetic theory considering the flow and melting phase change of ice slurry is developed. Sensitivity analysis of interphase forces is carried out. A comparison of the pressure drop, solid phase velocity, and heat transfer coefficient with empirical data is carried out, respectively. The calculated results are in good agreement with the experimental results, indicating that the numerical model could accurately describe the flow and melting characteristics. Thermophysical field distributions, the axial variation of ice volume fraction (IVF), recovery curve, the average heat transfer coefficient, as well as the re-uniformization length are obtained. After passing through local large heat flux segment, due to shear stress action, the IVF and the particle uniformity of the cross section have recovery characteristics. The gradient of the recovery curve decreases with increasing inlet IVF as well as with increasing Reynolds number. After the local large heat flux increases to a certain value, its effect on the recovery curve of the ice slurry is small. The re-uniformization length increases as the local large heat flux increases. The average heat transfer coefficient of local large heat flux segment increases due to damage of the boundary layer. These results can provide a theoretical basis for the design of ice slurry systems in practical application.

1. Introduction

Ice slurry, which is an appealing phase change cold storage material, consists of a mixture of carrier liquid and ice particles (which usually has a diameter of 0.1–1 mm) [1,2,3]. The solid and liquid two phase flow characteristics of ice slurry endow it with good fluidity, thermal, and transport properties. The substantial latent heat ice slurry (334 kJ/kg) [4] can effectively reduce the amount of pipe used in transit, lower the energy consumption of pumping, and reduce the size of the heat exchanger. It is used as both cold storage medium and heat exchange medium in various fields, being widely studied in recent years and presents a broad application prospect [5,6,7,8].

Ma et al. [9] provided a comprehensive review of the thermal fluid properties of different shell-less phase change slurry PCS from various aspects. Wang et al. [6,10,11] reviewed the characteristics of flow and heat transfer for ice slurry in cooling systems and studied the fluidity characteristic of ice slurry in horizontal, vertical, and 90° elbow pipes, respectively. It is easy to see that the characteristics of the ice slurry (additive concentration, additive type) as well as the structure of the flowing pipe and different thermal boundary conditions have a large influence on the heat transfer and flow characteristics of the ice slurry.

Flow and heat transfer experiments on ice slurry are commonly performed with round tubes, and scholars have conducted a series of experiments with different ice slurry. Grozdek et al. [12,13] used the 10.3% aqueous solution of ethanol while Lee et al. [14] experimentally studied the 6.5% ethanol-water solution ice slurry. Niezgoda-Zelasko [15,16] developed dimensionless relations based on experimental results to estimate local heat transfer coefficient values for different conditions. Kumano et al. [17] investigated the effect of ice slurry additive mass concentration on heat transfer and found that the effect was not significant. In addition to ice slurry transport in flat pipes, scholars have also conducted experimental studies on specially shaped pipes and heat exchanger coils. Mi et al. [18] conducted many experiments to study the performance of ice slurry flow through a plate heat exchanger with different mass concentrations. Ohira et al. [19] studied the effect of different shapes of pipes on pressure as well as heat transfer. The study of solid-liquid two-phase slurry is not limited to ice slurry. Experimental studies have also been conducted on cryogenic slurry. Li et al. [20] investigated experimentally and numerically the heat transfer characteristics of slush nitrogen and obtained a modified empirical correlation equation for the Nusselt number.

Researchers have also developed numerical models to simulate and analyze the flow characteristics of ice slurry, and the more commonly used model is the Euler-Euler model. Hu et al. [21] performed isothermal flow simulations of ice pigging using the Euler-Euler model. The numerical results indicated an uneven shear stress distribution along the flow. Gao et al. [22] employed an Euler-Euler approach considering interphase transfer mechanisms. It was found that the ice particle distribution as well as homogeneity is related to the velocity and the ice content. Cai et al. [23] conducted a numerical study on the heat transfer enhancement of ice slurry in horizontal pipes and found that pulsating flow can largely improve the heat transfer characteristics of ice slurry. Yadav et al. [24] investigated ice slurry flow in elliptical pipes under isothermal and non-isothermal conditions. The results showed that the melting of ice slurry in elliptical tubes was inhibited compared to that in circular pipes. Xu et al. [25] numerically simulated the flow of ice slurry in the inlet straight pipe of shell-and-tube heat exchanger. It was shown that the appearance of turbulence in the outlet cross-section and the degree of ice slurry stratification were related to the inlet velocity. Jin et al. [26,27] numerically investigated the flow characteristics of slush nitrogen. It was seen in some published experimental work that the pressure drop of slush nitrogen is larger than that of subcooled liquid nitrogen. Shi et al. [28,29] numerically investigated the characteristics of heat transfer of hydrate slurry in horizontal 90° bends and U-pipe.

It was found that the curved section leads to the appearance of secondary flows, vortices, and boundary layer separation. Simulation studies have also been performed using other numerical modeling methods, such as VOF models, single-fluid models, Boltzmann method, and other complex models. Li et al. [30] used VOF model to numerically study the heat transfer characteristics of ice slurry. Onokoko et al. [31] developed a 3D single phase model to study isothermal ice slurry flow in laminar and turbulent conditions. Onokoko et al. [32,33] proposed a simpler model for studying CFD models of isothermal ice slurry flow. Suzuki et al. [34] used the thermal immersed boundary–lattice Boltzmann method to simulate the interactions between particles and liquid in ice slurry flow as well as to study the Nusselt number correlation equation. Bordet et al. [35,36] used a model that captured more complex flow characteristics, particularly near-wall boundary layers and secondary flow, to study the ice slurry dynamical behavior. Languri et al. [37] quantified and discussed the heat transfer data of MPCM additives with different mass fractions in the base solution in detail, and obtained flow rate and heat transfer curve relationships.

The above is presented in an experimental study as well as a numerical study. The experimental part is classified by pipe structure into circular pipes as well as other complex pipes. The next part focuses on the numerical simulation studies, which are mainly introduced by numerical method classification, focusing on the studies conducted using the Euler-Euler model, enumerating the simulation studies, such as the VOF method and the single fluid method. A general review of the study of heat transfer characteristics of ice slurry is presented. However, there are few studies on the characteristics of flow-melt and recovery for ice slurry in the transportation pipe with local large heat flux segment, e.g., flowing through pumps, heat exchangers, heating elements, etc.

In the present paper, to further investigate the flow and melting characteristics of ice slurry as it flows through the pipe with local large heat flux segment, a 3D numerical model based on granular kinetic theory and considering flow and melting is proposed in Section 2. The numerical model developed is validated against the experimental data in Section 3. In Section 4, the results are given and discussed. The effects of various parameters on the thermophysical field distribution, axial IVF variation, average heat transfer coefficient, and the re-uniformization length are investigated in particular. The conclusion is provided eventually in Section 5.

2. Mathematical Model and Numerical Implementation

2.1. Physical Model

As ice slurry flows through the pipeline system, it may pass through pumps, flanges, heat exchanger, or other local heat leakage units. The horizontal straight pipe is divided into three sections, as shown in Figure 1. The inner diameter and length of the pipe are fixed to D = 16 mm, L = 1500 mm, respectively. The length of local large heat flux segment (below referred to as Section II) is fixed to L_heat = 100 mm, and the two low heat flux segments (below referred to as Section I and Section III, respectively) are each 700 mm long and symmetrically distributed on both sides of Section II.

2.2. Control Equations

The liquid and ice particles both are treated as interpenetrative continua in the Euler–Euler approach with interphase momentum and heat exchange. It is assumed that the flow of ice slurry is turbulent as well as incompressible, and that the ice particles are smooth, inelastic, and spherical.

• Continuity equation

Ice slurry is an incompressible fluid. The mass conservation equation satisfies the Navier–Stokes equation (N-S equation), and the expression of each phase equation is shown below:

$$\frac{\partial }{\partial t}({\alpha }_l{\rho }_l)+\nabla \cdot ({\alpha }_l{\rho }_l{\overrightarrow{\nu }}_l)={\dot{m}}_{sl}$$

(1)

$$\frac{\partial }{\partial t}({\alpha }_s{\rho }_s)+\nabla \cdot ({\alpha }_s{\rho }_s{\overrightarrow{\nu }}_s)={\dot{m}}_{ls}$$

(2)

$${\dot{m}}_{sl}=-{\dot{m}}_{ls}$$

(3)

where $\rho$ is the density, $\overrightarrow{\nu }$ is the local velocity, $\alpha$ is the local phase volume fraction, $\dot{m}$ is mass transfer rate caused by melting, and the subscript l represents liquid phase and the subscript s represents solid phase.

The sum of the volume fractions of each phase is one and the relationship is shown as follows:

$${\alpha }_l+{\alpha }_s=1$$

(4)

• Momentum conservation

The momentum equation of ice slurry can still be expressed by the N-S equation. The conservation of momentum equation for each phase is expressed as follows:

For liquid phase

$$\frac{\partial }{\partial t}({\alpha }_l{\rho }_l{\overrightarrow{\nu }}_l)+\nabla \cdot ({\alpha }_l{\rho }_l{\overrightarrow{\nu }}_l{\overrightarrow{\nu }}_l)=-{\alpha }_l\nabla P+\nabla \cdot {\tau }_l+{\alpha }_l{\rho }_{l}g+F_{sl}+({\dot{m}}_{sl}{\overrightarrow{\nu }}_s-{\dot{m}}_{ls}{\overrightarrow{\nu }}_l)$$

(5)

where $P$ is pressure, $g$ is gravity acceleration, and $F_{sl}$ is the interaction between each phase.

${\tau }_l$ is the pressure-strain tensor of the liquid phase, which expression is shown below:

$${\tau }_l={\alpha }_l{\mu }_l[\nabla {\overrightarrow{\nu }}_l+{(\nabla {\overrightarrow{\nu }}_l)}^T]-\frac{2{\alpha }_l{\mu }_l(\nabla \cdot {\overrightarrow{\nu }}_l)I}{3}$$

(6)

where ${\mu }_l$ is bulk viscosity, $({\dot{m}}_{sl}{\overrightarrow{\nu }}_s-{\dot{m}}_{ls}{\overrightarrow{\nu }}_l)$ is the momentum transfer introduced by the mass exchange due to solid phase melting.

Solid phase

$$\frac{\partial }{\partial t}({\alpha }_s{\rho }_s{\overrightarrow{\nu }}_s)+\nabla \cdot ({\alpha }_s{\rho }_s{\overrightarrow{\nu }}_s{\overrightarrow{\nu }}_s)=-{\alpha }_s\nabla P-\nabla P_s+\nabla \cdot {\tau }_s+{\alpha }_s{\rho }_{s}g-({\dot{m}}_{sl}{\overrightarrow{\nu }}_s-{\dot{m}}_{ls}{\overrightarrow{\nu }}_l)$$

(7)

where $P_s$ is the solid pressure derived from granular dynamics theory.

${\tau }_s$ is the solid phase stress, denoted as [38]

$${\tau }_s=(-P_s+{\zeta }_s\nabla \cdot {\overrightarrow{\nu }}_s)I+{\alpha }_s{\mu }_s[\nabla {\overrightarrow{\nu }}_s+{(\nabla {\overrightarrow{\nu }}_s)}^T]-\frac{2(\nabla \cdot {\overrightarrow{\nu }}_s)I}{3}$$

(8)

${\zeta }_s$ is the bulk viscosity of solid phase, which can be formulated as [39]

$${\zeta }_s=\frac{4}{3}{\alpha }_s{\rho }_s{d}_s{g}_0(1+e_{ss}){(\frac{{\theta }_s}{\pi })}^{\frac{1}{2}}$$

(9)

${\mu }_s$ is the solid shear viscosity which consists of three parts and can be represented as [40]

$${\mu }_s=\frac{4}{5}{\alpha }_s{\rho }_s{d}_s{g}_0(1+e_{ss}){(\frac{{\theta }_s}{\pi })}^{\frac{1}{2}}+\frac{10{\rho }_s{d}_s\sqrt{\pi {\theta }_s}}{96{\alpha }_s{g}_0(1+e_{ss})}{[1+\frac{4g_0(1+e_{ss}){\alpha }_s}{5}]}^2$$

(10)

where $g_0$ is the radial distribution function, $e_{ss}$ is the coefficient of restitution for collisions between particles, ${\theta }_s$ is the granular temperature, $d_s$ is the diameter of particles. $g_0$ modifies the probability of particle–particle interaction with the following expression [41].

$$g_0={[1-{(\frac{{\alpha }_s}{{\alpha }_{s,\max }})}^{\frac{1}{3}}]}^{-1}$$

(11)

• Energy conservation

The energy equation for each phase of the ice slurry is expressed as follows:

Solid phase

$$\begin{array}{l}\frac{\partial }{\partial t}({\alpha }_s{\rho }_s{H}_s)+\nabla \cdot ({\alpha }_s{\rho }_s{\overrightarrow{\nu }}_s{H}_s)=\nabla \cdot ({\lambda }_{eff,s}\nabla T_s)+{\tau }_s\cdot \nabla {\overrightarrow{\nu }}_s\hspace{1em}\\ -h_v(T_s-T_l)-({\dot{m}}_{sl}{H}_s-{\dot{m}}_{ls}{H}_l)\end{array}$$

(12)

Liquid phase

$$\begin{array}{l}\frac{\partial }{\partial t}({\alpha }_l{\rho }_l{H}_l)+\nabla \cdot ({\alpha }_l{\rho }_l{\overrightarrow{\nu }}_l{H}_l)=\nabla \cdot ({\lambda }_{eff,l}\nabla T_l)+{\tau }_l\cdot \nabla {\overrightarrow{\nu }}_l\hspace{1em}\hspace{1em}\\ -h_v(T_l-T_s)+({\dot{m}}_{sl}{H}_s-{\dot{m}}_{ls}{H}_l)\end{array}$$

(13)

where $H$ is the enthalpy, ${\lambda }_{eff}$ is the effective thermal conductivity, $h_v$ is the coefficient of volumetric interphase heat transfer, and $({\dot{m}}_{sl}{H}_s-{\dot{m}}_{ls}{H}_l)$ is the energy transfer between two phases introduced by mass transfer.

2.3. Turbulence Model

Since the $\mathrm{SST}\mathrm{k}-\mathsf{\omega }$ turbulence model (mixture) has higher accuracy and reliability in a wide range of flow fields, the conservation equations are shown as follows:

$$\frac{\partial }{\partial t}({\rho }_{m}k)+\frac{\partial }{\partial x_i}({\rho }_{m}k{\overrightarrow{\nu }}_m)=\frac{\partial }{\partial x_j}[({\mu }_m+\frac{{\mu }_m}{{\sigma }_k})\frac{\partial k}{\partial x_j}]+G_k-Y_k+S_k$$

(14)

$$\frac{\partial }{\partial t}({\rho }_m\omega )+\frac{\partial }{\partial x_i}({\rho }_m\omega {\overrightarrow{\nu }}_m)=\frac{\partial }{\partial x_j}[({\mu }_m+\frac{{\mu }_{t,m}}{{\sigma }_{\omega }})\frac{\partial k}{\partial x_j}]+G_{\omega }-Y_{\omega }+S_{\omega }$$

(15)

The constants used in equations are ${\sigma }_{k,1}=1.176$, ${\sigma }_{\omega ,1}=2.0$, ${\sigma }_{k,2}=1.0$, ${\sigma }_{\omega ,2}=1.168$, ${\alpha }_1=0.31$, ${\beta }_{i,1}=0.075$, ${\beta }_{i,2}=0.0828$.

2.4. Phase Interfacial Forces

The interaction between solid and liquid phase has a great impact on the solid phase kinetic. Therefore, interfacial effects need to be considered, mainly including drag force, lift force, and virtual mass force.

The expressions are given as follows.

• Drag force

Drag force plays an important role in two-phase flow [42], which is written as

$$F_{D,sl}=K_{sl}({\overrightarrow{\nu }}_s-{\overrightarrow{\nu }}_l)$$

(16)

$$F_{D,ls}=K_{ls}({\overrightarrow{\nu }}_l-{\overrightarrow{\nu }}_s)$$

(17)

where $K_{ls}(K_{sl})$ is the coefficient of interphase momentum exchange.

For ${\alpha }_s\le 0.2$ [38]

$$K_{sl}=\frac{3}{4}{C}_D\frac{{\alpha }_s{\alpha }_l{\rho }_l\left|{\overrightarrow{\nu }}_s-{\overrightarrow{\nu }}_l\right|}{d_s}{\alpha }_l^{-2.65}$$

(18)

For ${\alpha }_s>0.2$

$$K_{sl}=150\frac{{\alpha }_s^2{\mu }_l}{{\alpha }_l{d}_s^2}+1.75\frac{{\alpha }_s{\rho }_l\left|{\overrightarrow{\nu }}_s-{\overrightarrow{\nu }}_l\right|}{d_s}$$

(19)

$C_D$ is the drag force coefficient and ${\mathrm{Re}}_s$ is solid phase Reynolds number. The expressions are as follows

$$C_D=\frac{24}{{\alpha }_l{\mathrm{Re}}_s}\times [1+0.15{({\alpha }_l\times {\mathrm{Re}}_s)}^{0.687}]$$

(20)

$${\mathrm{Re}}_s=\frac{{\rho }_s\left|{\overrightarrow{\nu }}_s-{\overrightarrow{\nu }}_l\right|d_s}{{\mu }_s}$$

(21)

• Lift force

The formulation of lift force is given below [42]

$${\overrightarrow{F}}_{lift}=-C_l{\rho }_l{\alpha }_s({\overrightarrow{\nu }}_l-{\overrightarrow{\nu }}_s)\times (\nabla \times {\overrightarrow{\nu }}_l)$$

(22)

where $C_l$ takes a constant value of 0.2 [43].

2.5. Phase Change Model

Gunn [44] provides a heat transfer correlation for particles moving in liquid, as follows:

$$Nu=(7+10{\alpha }_l+5{\alpha }_l^2)(0.7{\mathrm{Re}}_s^{0.2}{\mathrm{Pr}}^{1/3}+1)+(-2.4{\alpha }_l+1.33+1.2{\alpha }_l^2)({\mathrm{Re}}_s^{0.7}{\mathrm{Pr}}^{1/3})$$

(23)

$$Nu=\frac{h_{sl}{d}_s}{{\lambda }_l}$$

(24)

$$h_{sl}=\frac{{\lambda }_l}{d_s}[(7+10{\alpha }_l+5{\alpha }_l^2)(0.7{\mathrm{Re}}_s^{0.2}{\mathrm{Pr}}^{1/3}+1)+(-2.4{\alpha }_l+1.33+1.2{\alpha }_l^2)({\mathrm{Re}}_s^{0.7}{\mathrm{Pr}}^{1/3})]$$

(25)

$\mathrm{Pr}$ is given as follows:

$$\mathrm{Pr}=\frac{{\mu }_{l}C{p}_l}{{\lambda }_l}$$

(26)

$h_v$ is given as below:

$$h_v=\frac{6{\alpha }_s{h}_{sl}}{d_s}$$

(27)

${\dot{m}}_{sl}$ is given as below:

$${\dot{m}}_{sl}=h_v\frac{T_l-T_s}{\Delta H}$$

(28)

The effective thermal conductivity of the liquid and solid phases in the near-wall region as well as in the main flow region depends on the local IVF and is expressed as follows [45].

• Near-wall region ($L/d_s\le 0.5$)

Energy transfer between particles and the wall could be considered based on point contact [46], as follows:

$${\lambda }_{eff,s}^w=3{\lambda }_l{\alpha }_{s}S\phi$$

(29)

$${\lambda }_{eff,l}^w={\lambda }_l{\alpha }_l$$

(30)

$$S=2L/d_s$$

(31)

$$\phi ={(\frac{A}{A-1})}^2[(1-\frac{1-S}{A})\ln (A)-\frac{A}{A-1}S]$$

(32)

where $L$ is the distance from wall, and ${\lambda }_{\mathrm{l}}$ is the real thermal conductivity of the liquid phase.

• Main flow region ($L/d_s>0.5$)

  $${\lambda }_{eff,s}^m={\lambda }_l\sqrt{{\alpha }_s}(A\beta +(1-\beta )K)$$

  (33)

  $${\lambda }_{eff,l}^m={\lambda }_l(1-\sqrt{{\alpha }_l})$$

  (34)

  $$K=\frac{1}{1-\frac{B}{A}}[\frac{A-1}{{(1-\frac{B}{A})}^2}\frac{B}{A}\ln \frac{A}{B}-\frac{B-1}{1-\frac{B}{A}}-0.5(1+B)]$$

  (35)

  $$A=\frac{{\lambda }_l}{{\lambda }_s}$$

  (36)

  $$B=1.25{(\frac{{\alpha }_s}{{\alpha }_l})}^{10/9}$$

  (37)

  $$\beta =7.26\times {10}^{-3}$$

  (38)

  where ${\lambda }_s$ is the real thermal conductivity of ice particles.

2.6. Boundary and Initial Conditions

The imposed boundary conditions have been vividly labeled, as shown in Figure 1. A uniform inlet velocity is specified for the inlet boundary as well as an inlet volume fraction. Turbulence intensity and hydraulic diameter are used. The turbulent intensity is taken as 3%, while hydraulic diameter is taken as 16 mm. Pressure outlet is used and specified as 0.1 MPa for the outlet boundary. No-slip condition and constant heat flux are adopted at the wall. In order to save calculation time, the initial conditions of the ice slurry in the flow field are specified as the same as the inlet conditions. The heat flux of the three sections is set to 2 kW/m², 40 kW/m², and 2 kW/m², respectively.

2.7. Numerical Solution Strategy

The three-dimensional numerical model of circular pipe with local large heat flux segment is performed using the commercial CFD solver (FLUENT-2020) to investigate the flow-melt characteristics of ice slurry. In the numerical calculation, the second-order upwind schemes are used to discretize all the governing equations. The equations are solved by phase-couple SIMPLE method. The UDF technology is applied to describe melting and thermal properties change of ice slurry. In the simulations, the time step is taken as 10⁻³ s and the numerical calculation is considered to converge when all the residuals are lower than 10⁻⁴ s and energy residual is lower than 10⁻⁸.

3. Numerical Model Validation

Ice particle distribution is more difficult to measure, while solid phase velocity as well as pressure drop are easier to obtain, so most of the previous experimental data revolve around the latter two. In this paper, the solid velocity and pressure drop are chosen as criteria for the validation of the flow model, together verifying the accuracy of the numerical model. Many scholars have made experimental measurements of ice slurry heat transfer properties, mostly under constant wall heat flux, and the heat transfer coefficient is used as the criterion for the validation of the phase change model in this paper.

3.1. Mesh Generation and Grid Independence

In order to improve the calculation accuracy of mesh and reduce the calculation cost, a meshing scheme of hexahedron grid + boundary layer for the horizontal straight pipe is used, as shown in Figure 2.

The number of grids has a significant impact on the accuracy and time-consumption of numerical calculation, so it is necessary to verify the independence of grids before numerical research. In this paper, seven grid schemes from 0.03 million to 1.6 million grids are compared in terms of solid velocity, IVF distribution, and pressure drop along the central axis, as shown in Figure 3. It can be seen that when the grid number reaches to 1.42 million, it has little effect on the calculation results, so the grid number of 1.42 million is adopted to study the following numerical calculation.

3.2. Sensitivity Validation

In this section, three interphase forces (virtual mass force, drag force, and lift force) are considered to test the sensitivity of the model. The drag force is set as imperative with a model of Gidaspow [38]. The tests are divided into three groups, considering simultaneously virtual mass force, drag force, and lift force, considering only virtual mass force and drag force, and considering only drag force and lift force. Figure 4 shows a comparison in terms of the distribution of solid velocity and IVF along the vertical cross-section of the outlet profile for the three groups. It is observed from Figure 4a that the selection of interphase forces has little effect on the solid velocity distribution. The distribution of solid velocity is more uniform about the center with the addition of lift force. It is seen from Figure 4b that the model considering drag force and lift force behaves the same as the model considering the three interphase forces simultaneously. Since the density difference between solid and liquid phases is not large, the virtual mass force can be neglected. The model considering only virtual mass forces and drag forces shows a different situation, with an increase in IVF in the upper and lower part of the pipe. IVF distribution is significantly influenced by lift force, which makes a smooth transition change in IVF. Therefore, drag forces and lift forces are considered in this numerical model.

3.3. Validation for Isothermal Ice Slurry Flow

Figure 5a shows a comparison of results between the numerical calculation of the ice slurry solid velocity profile with the experimental results of Vuarnoz [47] obtained by adopting the ultrasonic velocity profiler (UVP) to determine the velocity profile of ice particles. The solid phase velocity appears to be large in the center and small on both sides. As can be seen in Figure 5a, the numerical results are in good agreement with the experimental results, with errors within 10%.

Figure 5b indicates a comparison of results between the numerical calculation of ice slurry pressure drop with experimental results of Niezgoda-Zelasko [15] obtained from many experiments in a 0.016 m diameter horizontal pipe using 15.1% inlet IVF. It is seen from Figure 5b that the general trend of the solid velocity calculated agrees with the experiment values and shows a good linear regularity. The experimental data from the inlet velocity of 1–1.3 m/s shows that the data are not linearly varying, and it is conjectured that this may be caused by measurement errors. Since irregularly shaped ice particles are used in the experiment, consuming a large quantity of kinetic energy due to the collision of ice particles, the pressure drop data are smaller [21]. Neglecting the measurement error points, the remaining data are considered acceptable, mostly falling within 10% of the simulation data, so the numerical model can describe the ice slurry flow well.

3.4. Validation for Non-Isothermal (with Melting) Ice Slurry Flow

Figure 6 shows a comparison of results between the numerical calculation of the local heat transfer coefficient with experimental results of Niezgoda-Zelasko [16] obtain from many experiments in a 0.016 m diameter horizontal pipe using ethanol-water ice slurry with various inlet velocity of 0.5–2.5 m/s and 0.1 mm particle diameter. It can be seen from Figure 6 that the local heat transfer coefficient increases with the inlet velocity due to the increase in turbulence intensity. The calculated results match well with the experimental values, within 5% error from Figure 6.

The local heat transfer coefficient is represented by the equation below:

$$h_{lo,f}=\frac{\dot{q}}{T_w-T_{m,f}}$$

(39)

where $\dot{q}$ is the constant heat flux of the wall, $T_w$ is the temperature of the wall, and $T_{m,f}$ is the mean temperature of fluid on the center of cross-section.

4. Results and Discussion

The base case settings, where inlet IVF is 15%, inlet velocity is 1 m/s (Re_is = 1950 at the inlet) with q_I,III = 2 kW/m² and q_II = 40 kW/m², are presented and analyzed. The variation of IVF of sections I and II and the whole pipe is represented below as ΔIVF-I, ΔIVF-II, and ΔIVF-all respectively, which is defined as the inlet IVF minus the outlet IVF in the corresponding segment. The recovery characteristic of IVF after passing through the local large heat flux segment is defined as the outlet IVF minus the inlet IVF of Section III. The initial solution consists of 10.3% ethanol aqueous liquid, and the ice particle diameter is fixed at 0.1 mm.

4.1. Distribution Characteristics of Ice Slurry Thermophysical Field

Figure 7 illustrates the distribution characteristics of the physical fields for pressure field, IVF, temperature field, and velocity field at different x cross-sections of the pipe with Re = 1950 and IVF_in = 15%. For pressure field, the relative pressure is evenly distributed in the same cross section and gradually decreases along the flow direction. The pressure loss is attributed to the flow resistance caused by the fluid viscosity. Along the flow direction, collision friction between particles and between particles and wall can also lead to pressure loss.

For ice volume fraction field, it is seen from Figure 7 that the IVF in the mainstream for Section I is not influenced and remains the same as that of the inlet of the pipe. The ice particles absorb heat and melt mainly near the wall. Due to the density difference, the ice particles are suspended above the pipe. Therefore, the reduction of IVF for the lower side is slightly faster than that for the upper side. For Section II, due to the large local heat flux, the ice particles near the wall melt considerably and a liquid phase region appears near the wall along the flow direction. For Section III, under the effect of turbulent dispersion and buoyancy, the mainstream ice slurry is in a squeezing state for the upper near-wall liquid area due to a density difference, while it is in a flashing state for the lower near-wall liquid area. Due to the collision between the particles and the wall, the near-wall liquid phase region is disrupted, and after a certain recovery distance, the distribution of the IVF in the cross section becomes homogeneous again. Figure 8 shows the distribution of IVF along the flowing direction with local large heat flux value of 40 kW/m² at z = 0. It can be seen clearly that the IVF of upper region remains larger than that of the lower part owing to the two states of the mainstream ice slurry on the liquid area and density difference.

For temperature fields, it is seen that the near-wall ice slurry has an obvious temperature rise due to the large heat flux in Section II. The ice particles are easily melted, and the liquid phase is warmed up quickly due to the heat obtained. The temperature of the mainstream ice slurry hardly changes due to the weak conduction and most heat is used to melt the near wall ice particles.

For velocity fields, it is seen that the near-wall velocity approaches 0 m/s due to the no-slip boundary condition. The maximum value of the flow velocity and the proportion are gradually reduced along the flow direction due to the friction loss between the ice slurry and the wall and the momentum exchange between the particles in collision. High velocity zone gradually shrinks in a circular shape toward the center of the pipe. The velocity becomes consistent for a large area in the center of outlet of the pipe.

4.2. Effect of Inlet Velocity on the Melting Characteristics of Ice Slurry Flow

In order to investigate how the recovery of ice slurry is influenced by inlet velocity, five kinds of inlet velocity ($0.5\mathrm{m}/\mathrm{s},1\mathrm{m}/\mathrm{s},1.5\mathrm{m}/\mathrm{s},2\mathrm{m}/\mathrm{s},2.5\mathrm{m}/\mathrm{s}$) are applied. The corresponding Reynolds numbers varies between 975 and 4874, the inlet IVF is fixed to 15%, and the local large heat flux is set to 40 kW/m².

Figure 9a shows a comparison of the axial evolution of IVF for different inlet velocities. It can be seen that ΔIVF-I is strongly influenced by inlet velocity. The amount of melting increases with the decrease of Reynolds numbers. Taking Reynolds number 1950 as delineation, for the low Reynolds number region, the melting amount still increases more even in the low heat flux situation. For Section III, it is found that the IVF rises again along the flow direction which is the recovery phenomenon of IVF of ice slurry flowing after the local large heat flux segment. The gradient of recovery curve for different inlet velocity increases along with the decrease in velocity. For some high velocity conditions, the time for ice slurry flowing through the local large heat flux segment is reduced, the heat absorption is limited, the near-wall liquid phase region cannot be formed, and the introduced IVF variation is very small. Thus, the IVF distribution is less affected. The presence of local large heat flux does not cause much change or damage to the IVF field, so the recovery curve is flat. Under certain low velocity conditions, ice content of Section II is greatly influenced by large heat flux, which lead to the formation of a liquid phase region near the wall. When the velocity is low, the low heat flux also causes relatively obvious melting. The gradient of the recovery curve is still larger, as can be seen from Figure 9a, indicating that the degree of influence on the variation of IVF introduced by the local large heat flux far exceeds that caused by the melting in the low heat flux section.

Figure 9b shows the variation of IVF for each section with different Reynolds number. As can be seen from Figure 9b, ΔIVF-all under the studied conditions is largely determined by the local large heat flux segment. The mainstream ice slurry is less influenced by the heat flux and remains at a relatively high ice concentration, and the decrease in heat flux in section III also leads to a decrease in near-wall melting. Therefore, the IVF of cross-section increases and shows a recovery phenomenon of IVF along the flow direction. It is found that the recovery characteristics of IVF for Section III decrease gradually with increasing Reynolds number.

4.3. Effect of Inlet IVF on the Melting Characteristics of Ice Slurry Flow

To investigate how the recovery characteristics of the ice slurry are affected by the inlet IVF, five kinds of inlet IVF (5%, 10%, 15%, 20%, and 25%) are applied. The heat flux of Section II is set to 40 kW/m² and the inlet velocity is fixed to 1 m/s.

Figure 10a indicates the comparing of the axial evolution of IVF for different inlet velocities. It is seen that the amount of particle melting in Section II increases with the increase of inlet IVF. Maximum IVF variation is found in Section II. The variation curves of IVF of Section I and Section III are approximated horizontally in Figure 10a. Due to the limited contact heat transfer area and the point contact between the ice particles and the walls, the effect on heat transfer is insignificant even if the collision probability increases with increasing inlet IVF.

Figure 10b shows variation of IVF for each section with different inlet IVF. For section III, when the inlet IVF is low, the recovery curves of IVF are negative, indicating that the ice particles are continuously melting along the flow direction. It can be presumed that only if the inlet IVF exceeds a certain critical value will it have the recovery characteristic after a local large heat flux section. From Figure 10b, it can be seen that the critical value is approximately 8% for the present boundary conditions. The lower ice content in the mainstream region is not enough to offset the effects of the near-wall liquid phase region caused by the local large heat flux section. Thus, the IVF continues to decrease along the way. The reduction of IVF is 0.215% in Section I and 0.189% in Section III, which becomes smaller after the large heat flux section. This also reveals the recovery characteristics of volume fraction of ice slurry. As the inlet IVF increases, the recovery of the IVF becomes positive and increases accordingly. It is worth noting that ΔIVF-all remains almost constant. It can be seen that the inlet IVF is not an important factor for the amount of melting of the pipe compared to the inlet velocity. From the gradient of the corresponding pipe segment curve, it is clear that the inlet IVF has the same influence on ΔIVF-II and the recovery of IVF.

At the inlet IVF of 15%, it is observed that ΔIVF-all and ΔIVF-II are almost the same. ΔIVF-all is smaller or larger than ΔIVF-II when the inlet IVF is greater or less than 15%. Therefore, it is necessary to transport ice slurry with an ice content greater than 15% to minimize melting losses due to heat leakage through pumps and other connected components.

4.4. Effect of Local Heat Flux on the Melting Characteristics of Ice Slurry Flow

Figure 11 shows the IVF variation of ice slurry passing through each section with different values of the heat flux multiplier in section II, where the values of the heat flux multiplier in section II are 5, 10, 20, 30, 40, and 50 times the heat flux in section I, at an inlet velocity of 1 m/s and an inlet IVF of 15%, respectively. It can be seen that as the local large heat flux increases, the IVF decreases along the flow direction and the melting of ice particles near the walls increases. After the heat flux is greater than 20 times, the liquid phase zone is observed near the wall, which can be seen in Figure 8. From Figure 11a, it can also be observed that there is a significant IVF recovery along the Section III. As the local large heat flux increases, the gradient of the recovery curves first increases and then changes slightly. It is presumed that the gradient of the recovery curve is more dependent on the liquid phase region, and the gradient is essentially constant after the appearance of the liquid phase region (i.e., after the local large heat flux is greater than 20 times)

Figure 11b indicates the variation of IVF for each section with different local large heat flux. It can be seen more directly that the IVF recovery value increases first with increasing heat flux and then remains essentially constant. It can be deduced that increasing the heat flux has a negligible effect on IVF recovery when the local large heat flux exceeds a certain value. The ΔIVF-all increases linearly with the increasing multiplier.

Figure 12 shows the distribution of IVF for five local large heat flux multiples (5, 10, 20, 40, 50 times) at different x cross-sections. The inlet and outlet of the local large heat flux segment and three cross-sections of Section III are chosen to display the IVF. It can be seen from the legend of Figure 12 that the minimum value of IVF gradually increases along the flow direction in Section III. This indicates that the proportion of liquid phase in the near-wall region gradually decreases, and the cross-sectional IVF has a tendency to be re-uniform. From Figure 12c, under shear stress action as well as the effect of mixing, the near-wall liquid phase region is disrupted and the IVF in the near-wall region is increased. While the IVF around the mainstream decreases, the area of high ice content of the mainstream region shrinks toward the center. It is presumed from Figure 12d,e that the boundary layer is re-established after x = 1.267 m and the distribution of IVF tends to remain stable. Therefore, after the large heat flux, the IVF of the ice slurry demonstrates re-uniform distribution along the cross section.

Figure 13 shows the axial variation of IVF at six monitoring points on the same vertical line for six large heat flux multiples (5, 10, 20, 40, 50 times, respectively). In order to focus on the re-uniform of ice volume fraction distribution after passing through the local large heat flux segment, six near wall monitoring points are selected for comparison. The boundary layer is destroyed by the melting near the wall. A certain distance is needed to weaken the effect of inhomogeneous stratification due to melting, when the boundary layer is re-established and the distribution of IVF returns to a uniform state. Owing to slight difference in ice slurry in the mainstream, the near-wall region is applied to determine the re-uniformity of ice fraction distribution, and ice fraction distribution at x = 1.48 m is chosen as criterion. It is stipulated that a relative error within 2% is considered that the distribution of IVF is returned to uniformity. Relative error is defined as $({\alpha }_{x,y}-{\alpha }_{x=1.48,y})/{\alpha }_{x=1.48,y}\times 100\%$. In the case of different multiples, the axial position corresponding to the re-uniformization of 5 and 10 times is roughly x = 1.1 m, for 20 times is roughly x = 1.35 m, and for 30, 40 and 50 times are roughly x = 1.4 m. It is found that as the multiplicity of large heat flux increases, there-uniformization length becomes longer.

The firm line represents the upper data points, while the dashed line represents the lower data points. It is seen clearly that the IVF of the upper part is larger than that of the lower part due to the density difference. It is interesting to find that the minimum value of IVF does not appear at the outlet of the local large heat flux section but has a certain delay, as seen in Figure 13, basically at x = 0.85 m. It can be found that there is a relatively obvious recovery of the IVF in the x cross-section after the local large heat flux section, as well as a re-uniformity.

Figure 14 indicates the average heat transfer coefficient for each section at different multiples of local large heat flux. It can be noticed that in Section II, the average heat transfer coefficient increases with the increasing multiples of heat flux due to the near-wall melting. This is because the boundary is severely damaged, which strengthens the heat and mass transfer between the mainstream region and the near-wall region. The average heat transfer coefficient of Section III changes slightly for 5–20 multiples while decrease with the increase of heat flux for multiples bigger than 20. This is presumed to be caused by the thermophysical properties of the ice slurry as well as changes in the distribution of ice particles [32]. The viscosity as well as the effective thermal conductivity of the ice slurry increases with increasing ice concentration. In Section III, the temperature of the ice slurry increases, the IVF decreases, and the viscosity decreases, leading to a corresponding increase in the heat transfer coefficient. However, the effective thermal conductivity decreases with decreasing ice content, and the heat transfer coefficient should decrease accordingly. From Figure 14, it can be seen that when the heat flux is greater than 20 times, the heat transfer coefficient decreases. From this, it can be inferred that the effective thermal conductivity may bring more influence than the viscosity on average heat transfer coefficient.

5. Conclusions

The flow-melting characteristics are simulated for ice slurry in a straight horizontal pipe with a localized large heat flux section. The flow-melt characteristic of ice slurry under different boundary conditions is obtained. Axial recovery and re-uniformization of distribution of IVF after local large heat flux segment are investigated. The conclusions drawn are as follows.

(1)

A 3D Eulerian-Eulerian model based on the granular kinetics theory considering the $\mathrm{SST}\mathrm{k}-\mathsf{\omega }$ turbulence model and phase change is developed, which can accurately describe the heat transfer characteristics of ice slurry.

(2)

When the local large heat flux is 20 times higher than the low heat flux, a liquid phase region will appear near the wall.

(3)

After passing through the local large heat flux segment, the IVF and the particle uniformity of the cross section have recovery characteristics due to shear stress action. The degree of recovery is determined by the degree of disruption of the IVF field by the local large heat flux.

(4)

The re-uniformization length increases with the increasing local large heat flux. Interestingly, the minimum value of IVF does not appear at the outlet of the local large heat flux section, but with a certain delay.

(5)

The average heat transfer coefficient of the local large heat flux segment increases due to the boundary damage.