Thermal free-surface immersed-boundary lattice Boltzmann method for free surface flows with a liquid-solid phase transition

. This paper reports on the progress of the liquid-solid phase transition modeling of water in open channel flow by using the lattice Boltzmann method with the immersed boundary modification. The phase transition in a fluid flow has a moving interface between the liquid and solid state, which leads complicated treatments in existing numerical models. By applying the immersed boundary modification in the lattice Boltzmann method and the non-iterative enthalpy approach for the separation of the states, the moving boundary of the melting or solidification front is solved without any difficulty. The ice bed and the submerged ice cover under dynamic flow conditions is exercised to demonstrate the model performance. The model is extremely suitable in the formulation in terms of its simple and compact framework extendable to any dimensions.


Introduction
The ice and water mixture is the one of the common objects in cold regions. Ice acts like a solid in water flow and the situation will become even more difficult when thermodynamic is involved in the system. This complicated physics of ice and water is the motivation of this study. The liquid-solid phase changes are often referred as the Stefan problem [1], which has complexity having moving boundary showing a liquid-solid interface in time. Solving moving body or boundary brings many difficulties for conventional methods [2]. State-of-Art methods for phase change are limited to enclosed region and are based on fixed grid-based approach and front tracking approach. However, these approaches require to solve additional sets of equation and boundary condition on the moving boundary interface [3,4]. Recently, particle-based method appears to provide solutions and visual simulations of phase change in free surface flow and performance of this approach will cost in large scale simulations [5]. The phase change of water in free surface flow is important physical phenomena in engineering and nature, however, no particular computational effort is presented in the literature to the best of authors knowledge.
In this paper, we present a numerical procedure for liquid-solid phase change of water by using the Lattice Boltzmann method (LBM) modified by the immersed boundary term. The numerical procedure is composed as follows. Two distribution functions approach was used to account flow field and temperature field on fixed grids in two dimensions. The D2Q9 lattice arrangement was used for both field of distribution functions (refer to [6] for details). The enthalpy-based non iterative method was used for phase change of fluid counting on the local enthalpy updates [7]. For open channel fluid flow, free surface formulation of single-phase LBM is employed [8]. Integration of the immersed boundary modification and enthalpy-based non iterative technique in the free surface flow is straightforward.

Lattice Boltzmann method for free surface flows with phase transitions
For fluid flow with free surface representation, a discretized Boltzmann equation with the Bhatnagar-Gross-Krook collision operator [9] and modification of the immersed moving boundary [10] is to be solved by collision and streaming steps as: where is the discrete unit velocity in the direction, is the dimensionless relaxation time with respect to the lattice viscosity and is adjusted with the sub-grid scale turbulent model and is the parameter given by in which , is the liquid fraction value of the cell, which takes a value between 0 and 1. Liquid fraction values of 0 and 1 represent ice and water, respectively. In Eq. (2), the total relaxation can be used instead of the relaxation time . The immersed boundary modification can be used for not only dynamic separation of solid (ice) and liquid (water) phases, but also for a moving body (moving ice) condition in a Thermal free-surface immersed-boundary Lattice Boltzmann method 3 fluid flow. Therefore, an additional collision term is for cells partially or fully covered by a solid, i.e., ice cell, is given as where is the velocity of the moving solid [10]. The macroscopic variables, namely density ρ and velocity , can be computed by the orders of the moments of the distribution functions as Free surface of the liquid is tracked by volume fraction values defined by the ratio of cell mass to cell density. The cell mass is only evaluated on interface cells by distribution functions and is a unit for every liquid cell as water or ice. Solution of free surface in LBM is known as a single-phase model, so that no physical variables will be defined for empty cell which represents gas state as air. Influence of gas state realized by a free surface boundary condition which assumes that the fluid has a much lower kinematic viscosity than the gas state and is expressed in terms of the following distribution function on a free surface interface node: where is the gas density implicitly acting as air pressure onto the free surface. In the modelling of heat transport with phase change, the temperature field is considered to be an essential variable and can be calculated by the following thermal lattice Boltzmann equation with latent heat of fusion [11]: where , is the distribution function for the temperature field, 3 1/2 is the dimensionless relaxation time with respect to the thermal diffusivity , is the dimensionless latent heat of fusion, and is the specific heat capacity of water or ice. The equilibrium distribution function of the temperature field can be given as and the macroscopic temperature can be converted from dimensionless temperature as follows: 4 Ayurzana Badarch, Tokuzo Hosoyamada . 8 After the dimensionless temperature evolution, the local enthalpy [12], obtained by , ∆ , can be used to linearly interpolate the liquid fraction, and the liquid fraction defines the liquid (water) and solid (ice) phases in the domain. The phase, i.e., ice or water, is assigned to F cells because ice, which acts like a solid, will become water after melting. Moreover, if ice interacts with the air (G cells), boundary cells between the ice and air must be IF cells because IF cells have a certain water content. At the interface between ice and water, takes a value of between 0 and 1. The interaction between the free surface flow module and the heat transport with the phase change module is such that the temperature difference produces a buoyance force in the flow field, and the flow field affected by the buoyance force forms a convective temperature field in the domain. Although the buoyance force is negligible in a turbulent flow, it must be included in the computation. The lattice viscosity is related to the lattice thermal diffusivity of a fluid as / , where is the Prandtl number, so that the relation between the computational modules is maintained.

Numerical simulations and discussions
The thermal free-surface LBM is implemented as computation code and validation can be found in [6] and [13]. The validated code was then extended to include a phase transition without the latent heat source term and tested on the melting of the ice-slab and the melting of the ice cube in an ambient air [14]. With the following numerical tests, the application of the proposed model is elaborated.

Melting of an ice bed
In the early stage of model development, ice bed melting by flow over a weir was simulated in [6] for ice that was ready to melt, i.e., the ice temperature was set to 0°C, and the latent heat was ignored. In the present study, we included a latent heat source term in heat transport and extended the temperature range. The initial and boundary condition is indicated by the problem geometry, as shown in Fig. 1a. Initially, temperatures of -30°C, 20°C, and 30°C were set for ice, air, and water, respectively. The inlet and outlet boundary of the flow field was imposed with a velocity boundary condition [15], whereas the wall and the surface of the weir were modeled as slip walls. Water at a temperature of 30°C was supplied to the inlet, where the thermal boundary was given by the Dirichlet boundary condition [16]. The secondorder extrapolation boundary condition was imposed on the outlet boundary for the heat transport module [17]. The other walls and the surface of the weir were assumed to be adiabatic. Weir flow by the free surface LBM was carefully investigated and validated by [8], and we herein used a weir flow with a Froude number of 0.13. In order to examine grid independence, we considered two grid resolutions, namely, ℎ 60 and ℎ 80, where h is the grid number used for the weir height, as shown in Fig. 1a. The relaxation times for the grid resolutions for the flow field were chosen as = 0.526 and 0.534, respectively, which is adjusted to by the sub-grid scale model [18]. The relaxation times for the heat transfer module was determined by using the thermal diffusivity of water, which can be connected to the lattice viscosity by the Prandtl number. The ratio of the remaining ice area to the initial ice area was measured and is shown with respect to melting time in Fig. 2. Depending on the parameters, both grid resolution and the selected relaxation times, the total times of melting differed by approximately 0.6 min, and in the case of ℎ 60, the ice lasted 3.8 min. Based on these considerations, choosing appropriate parameters based on their physical relations is more important than grid resolution. The melting rates for these two cases have similar melting rate tendencies, as shown in Fig. 2, but after approximately 3.3 min melting rates are changed due to the low melting rate of the ice located directly behind the weir. Between the nappe entrance and the weir, where the flow is partially circulated, the convective heat transfer between the water and the ice was small due to the low velocity in this region.
Heat transfer between water and ice can clearly be explained in terms of the Nusselt number. The local Nusselt number at a point on the melting front is defined as [19]: where is the grid number for the length of the ice, and is the grid number for the depth of the water above the ice. The local Nusselt number and ice depth are plotted with respect to the elapsed time in Fig. 3. The local Nusselt number at 0.7 m fluctuates with a higher frequency and amplitude than that at 1.0 m, and the tendency in both cases is to increase with time until decreasing suddenly at the location at which the ice depth decreases. Generally, this tendency is due to the fact that heat increases in ice and water near the interface, and the sudden drop is due to the disappearance of the ice. As the frequency increases, ice at approximately 0.7 m from the origin quickly melted because the convective heat Thermal free-surface immersed-boundary Lattice Boltzmann method 7 transport at this position is high. The heat transfer coefficient expressed in terms of in Fig. 4 was approximately ℎ / 1,090.8 WK -1 m -2 , where is the characteristic depth of water on the ice. The temperature field and ice/water phase with a free surface at three different times are plotted in Fig. 5, where the interaction of the flow structure and the thermal behavior of phases are shown clearly.

Melting of a submerged ice cover
A submerged ice cover melting simulation was conducted in order to determine the applicability of the numerical model to real field problems, in which ice in a river or a reservoir is mixed with free surface flows. In this case, the fixed position of the ice will help in the visualization of the freezing and melting of water. Therefore, the ice does not break up and move freely. The motion of ice in the model was studied, where we analyzed floating ice on the water surface [20] in short time.
As shown in Fig. 1b, the condition and geometry of the submerged ice cover are such that heat is absorbed by the bottom of the ice and the water is expected to freeze downward due to water being trapped under the ice cover, where natural convection flow may dominate. Excluding the outlet velocity condition, all of the parameters and given conditions were the same as in the ice bed melting simulation.   6 shows the melting rate of ice for the two simulation cases, ℎ 60 and 80, and a discrepancy between the time courses appears at around 1.0 min in the simulation. However, a comparison between those two cases of the general melting shape of ice cover revealed no significant differences. The differences of the curves show that the finer grid might provide more detailed dynamics of ice cover melting and freezing in shape. Fig.7 shows the general melting process of the ice cover in the simulation. The ice cover absorbs heat from water flowing over it, and the cooled water flows out of the outlet boundary. At the same time, beneath the ice cover, heat can only be transported by conduction, as shown in Fig. 7a, until a natural convection flow forms in the water region beneath the ice cover, because the water in this region is trapped by the circulating water flow near the outlet boundary, as shown in Fig. 7d. The circulation flow near the outlet boundary transports heat and momentum into the closed region beneath the ice cover. Another circulation flow was observed where the water flows over the weir. This circulation convects heat to the ice behind the weir. As water flows over the ice cover, melting occurs and the ice is gradually eroded by the overtopping Thermal free-surface immersed-boundary Lattice Boltzmann method 9 flow (Fig. 7b). The melting rate of the upper surface of the ice cover was higher than in other parts of the melting ice. This situation continues until a natural convection flow forms in the water region beneath the ice cover. As expected, freezing occurred on the bottom surface of the ice cover due to convection. However, an opening eventually formed near the back face of the weir through which water could flow into the region of water trapped by the circulation flow, as shown in Fig. 7e. The flow then surrounds the ice cover and melts the ice from all sides. The outlet circulation flow was an active heat transporter and melted and sharpened the tail of the ice cover, as shown in Figs. 7a and 7e. Due to erosion, the ice cover was split into two pieces by the overtopping flow after approximately 1.6 min, which is indicated by the recurve shape of the lines in Fig. 6. The piece near the outlet boundary quickly melted because it was surrounded by an active flow field. The small piece of ice remaining can be seen in Fig. 7f, and its effect on the temperature field can be seen in Fig. 7c. The piece of ice near the weir eventually extended downward to the bottom boundary, as shown in Fig. 7f. Fig.8 shows the temperature and velocity profiles at 0.7 m and 1.0 m at various times. Fig.8 shows that the velocity and temperature profiles have the same tendency due to the massive amounts of heat transport by convection in turbulent flow. The vertical temperature gradient can be high where the flow velocity is high, as indicated near the upper surface of the ice cover in Figs. 8a and 8b. Unrealistic velocity and temperature decreases in the middle of the profiles appear at 1.0 m in Fig. 8c because of the piece of ice remaining in the flow, which cannot move with the water flow. As such, this piece of ice influenced the velocity field as an obstacle, resulting in a cooler temperature distribution. Thermal free-surface immersed-boundary Lattice Boltzmann method 11

Conclusion
We have tested our numerical procedure on the ice bed and submerged cover melting in open channel flow. Validation of numerical procedure on liquid-solid phase changes in closed region can be found in our previous publications. The numerical procedure does not require any additional equation or condition to treat moving boundary of melting or solidification fronts. Exploiting the liquid fraction value is useful for multiphysics simulation involving solid and liquid motion in the system. In that concept, we have adopted the immersed boundary modification for the LBM to separate liquid and solid states in assistance of liquid fraction values. The immersed boundary modification provides stable dynamic interactions between liquid and solid states. The performance of the numerical model on the two example tests was promising. Nevertheless, the accuracy of the numerical procedure needs to be validated against experimental study in the open channel flow with ice. In future studies, ice motion and ice-ice interaction need to be considered in order to study ice jam or ice water mixture in the river. The large simulation involving the large time and space scale can be solved in parallel computational architectures, in which current model is perfectly suitable.