Lattice Boltzmann Modeling for Melting/Solidification Processes

The phenomena of melting and solidification are associated with many practical applications, such as metal processing, castings, environmental engineering, thermal energy storage system in space station and many more. In these processes, matter is subject to a change of phase and consequently, a boundary separating two different phases evolves and moves within the matter. Mathematical modeling of such 'moving boundary problems' are always a challenging task because of the dynamic evolution of the phase separating boundary, complex boundary conditions as well as varying thermophysical properties. Many macroscopic mathematical modeling strategies for the solidification/melting problems can be found in the contemporary literatures. An excellent review in this regard can be found in Hu & Argyropoulos (1996). Early efforts in melting/solidification modeling initiated with a moving/deforming grid approach (Rubinsky & Cravahlo, 1981; Voller & Cross, 1981; Voller & Cross, 1983; Weaver & Viskanta, 1986; Askar, 1987), in which independent conservation equations for each phase need to be initially formulated, and are to be subsequently coupled with appropriate boundary conditions at the inter-phase interfaces. However, such multiple region solutions require the existence of discrete interfaces between the respective phases. In fact, a major difficulty with regard to their implementation is associated with tracking of the phase interfaces (which are generally unknown functions of space and time). The need for moving numerical grids and/or coordinate mapping procedures complicates the application of this technique further, and generally, simplifying assumptions regarding the geometric regularity of the interfaces are made. Additionally, a serious limitation exists for modeling phase change behavior of multicomponent systems, since, unlike pure substances; such systems do not exhibit a sharp interface between solid and liquid phases, in a macroscopic sense. The phase-change behavior of such systems depends on many factors including the phase-change environment, composition, and thermodynamic descriptions of specific phase transformations. Moreover, solidification occurs over extended temperature ranges and solid formation often occurs as a permeable crystalline-like matrix which coexists with the liquid phase. In such cases, it would be virtually impossible to track a morphologically complex zone in a macroscopic framework, using any moving grid technique. In contrast, in fixed-grid mathematical models of phase change (Comini et al., 1974; Morgan et al., 1978; Roose & Storrer, 1984; Dalhuijsen & Segal, 1986; Pham, 1986; Dhatt et al., 1989; Comini et al., 1990; Voller et al., 1990), transport equations for individual phases are volume-averaged to

come up with equivalent single-phase conservation equations that are valid over the entire domain, irrespective of the constituent phases locally present. A separate equation for evolution of liquid fraction is solved in conjunction with the above set of conservation equations, which implicitly specifies and updates the interfacial locations with respect to space and time. It can be noted at this point that although simulation strategies mentioned as above have become somewhat standardized over the past few decades, solution of phase-change problems over multiple length scales still poses serious challenges, primarily because of the disparate and coupled length scales characterizing the entire sequence of transport processes. To overcome such difficulties, phase-field models of dendritic solidification have been developed and proposed by several researchers (Mikheev & Chernov, 1991;Kim et al., 1999;Harrowell & Oxtoby, 1987;Khachaturyan, 1996;Beckermann et al., 1999;Tong et al., 2001). Advantage of the phase field models lies in the fact that computational difficulties associated with front tracking are eliminated by introducing an auxiliary order parameter (the so-called phase field) that couples with the evolution of the thermal field. Dynamics of the phase field are designed to follow the evolving solidification front, thereby eliminating the necessity of any explicit front tracking. Recently, the multiscale mesoscopic lattice Boltzmann (LB) (Kendon et al., 2001;Sankaranarayanan et al., 2002;Barrios et al., 2005) method has emerged to offer huge potentials for solving complex thermofluidic problems involving morphological development of complicated phase boundaries such as the problem of phase separation of two immiscible fluids . Such a method, typically, considers volume elements of fluid comprising of a collection of particles that are represented by characteristic particle velocity distribution functions defined at discrete grid points. The rules governing the collisions and subsequent relaxations are designed such that the time-averaged motion of fluid particles becomes consistent with that predicted by the Navier-Stokes equation. Further advancements in LB modeling of fluid flow enabled the research community to explore more complicated problems addressing flow through porous medium and a few generic cases of multi-phase flow (Gunstensen et al., 1991;Shan & Chen, 1993;Ferreol & Rothman, 1995). In this context, it can be mentioned here that a distinct advantage of the LB method for modeling solid-liquid phase transitions, in comparison to a classical continuum based formulation, lies in the fact that the LB method is fundamentally based on microscopic particle models and mesoscopic kinetic equations, which means that micro and meso-scale physics of phase transitions can elegantly be incorporated. Another important advantage is that it does not require an immediate explicit calculation of fluid pressure, leading to time-efficient computational simulations. Further, LB models are inherently parallelizable, which renders their suitability to address phase change processes over largescale computational domains. The LB approaches proposed so far for modeling solid-liquid phase transition problems can broadly be categorized into two major groups, viz. (a) phase field based methods following the Ginzburg-Landau theory and (b) enthalpy based methods. De Fabritiis et al. (1998) developed a thermal LB model for such problems by employing two types of quasiparticles for solid and liquid phases, respectively.  proposed a simple reaction LB model with enhanced collisions, using a single type of quasiparticle and a phase field approach. Further work proceeded along similar lines (Miller & Schroder, 2001;Miller, 2001;Miller & Succi, 2002;Miller et al., 2004;Rasin et al., 2005;Medvedev & Kassner, 2005), with the phase-field model acting as a pivotal basis for determining the evolution of respective phase fractions. It needs to be emphasized that one of the major problems in implementing the phase field based methods in the context of solid-liquid phase transition problems is the requirement of limitlingly finer grid spacing for resolving the interfacial region to reproduce the dynamics of the sharp interface equations. Consequently, adaptive mesh refinement strategies involving computationally involved data structures are required for problems subjected to small undercooling. On the other hand, enthalpy-based models have been extensively used to solve complex solidification problems over macroscopic and mesoscopic length scales. An extended LB methodology, in conjunction with an enthalpy formulation for treatment of solid-liquid phase change aspects in case of diffusion dominated problems, was first introduced by Jiaung et al. (2001). Subsequently, taking the computational advantage of the enthalpy-based technique, Chatterjee & Chakraborty (2005, Chakraborty & Chatterjee (2007) and Chatterjee (2009Chatterjee ( , 2010 proposed a series of LB models primarily applicable to a wide range of melting-solidification problems. It was started with an enthalpy based LB model for diffusion dominated phase transition problems, followed by a hybrid LB method for generalized convection-diffusion transport processes pertinent to melting/solidification problems. In the diffusion models, the temperature field was obtained from an evolution equation of a single particle density distribution function (DF), whereas in the convection-diffusion models the thermal field is described by a novel enthalpy density DF through a kinetic equation based on the total enthalpy of the phase changing system or alternatively from an evolution equation of temperature. The analysis of solidification in a semitransparent material using the enthalpy based LB method was performed by Raj et al. (2006). The radiative component of the energy equation in the LB formulation was computed using the discrete transfer method in their model. Recently, Huber et al. (2008) developed a multiple DF LB model for coupled thermal convection and pure-substance melting, where the two DFs were interrelated through the buoyancy term and the equilibrium DF of the temperature kinetic equation. In this chapter, we describe a straightforward technique for simulating solid-liquid phase transition, by coupling a passive scalar based thermal LB model with a fixed-grid enthalpyporosity approach (Brent et al., 1988) that is consistent with the microscopic solvability theory. The macroscopic density and velocity fields are simulated using a single particle density DF through a kinetic equation, while the macroscopic temperature field is obtained from a separate temperature DF through another kinetic equation (Chatterjee, 2010). The phase change aspect is numerically handled by the enthalpy-porosity technique with an adapted enthalpy-updating scheme. The source terms originating out of the physical situation are incorporated into the respective kinetic equations by the most formal technique following the extended Boltzmann equation. Test cases for one and two-dimensional solidification problems are presented and compared with the analytical and available numerical solutions. Finally, simulation results for a popular solid-liquid phase change problem, such as the Bridgman crystal growth in a square crucible are also shown to establish the capability of the model.

Model formulation
In this section we first present the generalized convection-diffusion macroscopic conservation equations governing the transport processes occurring during phase change, followed by the corresponding lattice Boltzmann formulation.

Macroscopic conservation equations
The equivalent single phase volume-averaged macro-scale continuity, momentum and energy conservation equations for a nonparticipating phase changing system, assuming a laminar, incompressible and Newtonian flow, can be presented as: where u , p and T denote the macroscopic velocity, pressure and temperature,  SK u , where K represents the permeability tensor. Components of the tensor K depend on the specific morphology of the phase changing domain, for which any appropriate formulation for flow through a porous medium can be effectively invoked. A common approach followed in the literature is to adopt the celebrated Darcy model (or some of its variants) for flow through a porous medium, in association with the Cozeny-Karman equation (Voller & Prakash, 1987) , where  is the dynamic viscosity, where T s and T l are the solidus and liquidus temperatures respectively. The above formulation effectively ensures that in the phase changing cells, the porous medium resistance term dominates over the transient, convective and diffusive effects manifested by molecular interaction mechanisms, thereby forcing the velocity field to imitate the Cozeny-Karman law. On the other hand, in totally solid elements   0 l f  , the high porous medium resistance forces any velocity predictions effectively to zero. In a fully liquid element   1 l f  , however, this term has no consequence, and the usual form of the Navier Stokes equation can be retrieved. The latent-heat evolution is accounted for by introducing a source term in the macroscopic energy conservation equation (final term on the right-hand side of Eq. (3)) as,

The lattice Boltzmann model
A statistical description of a fluid system can be made in terms of a particle density DF, which satisfies the continuous Boltzmann equation with a single-relaxation-time BGK (Bhatnagar-Gross-Krook) model (Bhatnagar et al., 1954) as: is a single particle density DF from which the macroscopic properties of the fluid can be obtained, ξ is the microscopic velocity, f  is the relaxation time, F and   ,, are the external forcing parameter and the Maxwell-Boltzmann type equilibrium DF, given by, where R is the gas constant and D is the dimensionality. The macroscopic variables are obtained by taking (microscopic velocity) moments of the density DF f as: However, it is a well known fact that the temperature field obtained from the second moment of the DF f yields a fixed Prandtl number, implying that the thermal diffusivity cannot be adjusted independent of the kinematic viscosity (He et al., 1998), which restricts its applicability to a limited class of problems only.
To this end, we define a new temperature DF   ,, g t x ξ following the passive scalar approach of He et al. (1998), which obeys a kinetic equation of the form: where g  is the relaxation time,

 
,, eq gt x ξ is the is the Maxwell-Boltzmann type equilibrium DF, given by, Here II II I I    accounts for the viscous heating ( I  ), compression work ( II  ) and kinetic energy ( III  ) contributions given as (Shi et al., 2004): It should be mentioned that in the incompressible limit (which is the limit of small Mach number, 0 Ma  ), II  and III  become negligible and hence I   . In Eq. (7), Q is the www.intechopen.com energy source term that features out of the physical situation within the participating fluid media. The macroscopic temperature can be obtained from the temperature DF g as: Since two separate kinetic equations with corresponding DFs are used to describe the flow and thermal fields respectively, the kinematic viscosity and thermal diffusivity can be independently adjusted, which makes the model suitable for varying Prandtl number flows. Eqs. (5) and (7) can now be represented by the following generic discrete-velocity form: for the respective kinetic equations and   1, ib  stands for the b base vectors of the underlying lattice type. Eq. (10) is subsequently integrated along its characteristic using the second order trapezoidal rule (He et al., 1998) to yield the discrete evolution equation: where t  denotes the time step. The forcing parameter and the discrete equilibrium DFs can be constructed as: Finally, the energy sources can be formulated as: The weights i w and the discrete velocities i ξ correspond to the D2Q9 configuration (He et al., 1998) In order to avoid implicitness of Eq. (11), we further introduce (He et al., 1998) 22 The macroscopic flow and thermal quantities are obtained from i f and i g as: www.intechopen.com

Enthalpy update
For accurate prediction of the liquid fraction, the latent enthalpy content of each computational cell needs to be updated according to the predicted macroscopic value of temperature each iteration within a time step. For that purpose, an enthalpy updating scheme in accordance with the formulation of Brent et al. (1988) is used, which is of the , where n is the iteration level characterizing the updation stage, h s is the sensible enthalpy of the concerned cell, h  is the latent heat contained by the cell, and  is a suitable relaxation factor to smoothen convergence. In the above formulation, Fh   is an appropriate mathematical function, which needs to be constituted in consistency with microscopic phase change considerations. A detailed guideline in this regard can be found in Chakraborty & Dutta (2001

Boundary conditions
A no-slip hydrodynamic boundary condition and both Neumann and Dirichlet type thermal boundary conditions at the walls can be used. The non-equilibrium extrapolation method (Guo et al., 2002), which has a good numerical accuracy and stability, can be adopted to implement the above mentioned boundary conditions in the LB framework. According to this method, the non-equilibrium part of the DF at a boundary node can be well approximated by the same at the nearest neighboring node in the fluid region along the discrete velocity. As an example, if b x represents a boundary node and f x its nearest neighboring fluid node, then x can be given as: where the equilibrium part of the DF is determined by imposing the macroscopic boundary conditions through the auxiliary density   , velocity  u or temperature T  . For example, if with ff x as the next neighboring fluid node of b x in the same direction, and xx .  u and T  are specified according to the given boundary conditions for u and T . A Neumann boundary condition can be implemented by transferring it to the Dirichlet type boundary condition by using a conventional second order finite difference approximation to obtain the boundary temperatures (Shu et al., 2002), in an iterated manner. Regarding interface conditions, it is apparent that the solid/liquid interface in phase change problems www.intechopen.com acts as a wall, and the same needs to be treated appropriately. However, according to the enthalpy-porosity formulation, one does not need to track the interface separately and impose hydrodynamic or thermal boundary conditions on the same, since the interface comes out as a natural outcome of the solution procedure itself.

Numerical scheme
The simulation starts with the prescribed initial values of the temperature   The relaxation parameters should lie in the range 0.5 1, , i i fg    such that positive distribution functions can be obtained close to the local equilibrium, thereby ensuring nonlinear stability of the numerical scheme (Higuera et al., 1989). It should be emphasized here that a rigorous, exact, theoretical analysis of nonlinear stability of the scheme is impossible, for it would amount to solving the lattice Boltzmann equations (LBE) itself. However, a number of general guiding criteria prove fairly useful. One of these criteria is the conservativeness of the scheme. The streaming operators in the LBE are perfectly conservative and the collision operators are also conservative. This makes the method an exactly conserving numerical scheme which automatically protect against numerical blowups in the actual simulation (Chatterjee, 2009

Case studies
We present here some case studies such as 1-D and 2-D solidification/melting problems for which analytical solutions are available and some other benchmark problems in melting and solidification.

1-D directional solidification
A one-dimensional (1-D) directional solidification problem is solved for which analytical solution is available. The schematic of the problem is shown in Fig. 2. Initially, the material is kept in a molten state at a temperature T i (= 1) higher than the melting point T m (= 0.5). Heat is removed from the left at a temperature T 0 , which is scaled to be zero. The onedimensional infinite domain is simulated by a finite domain (considering a domain extent of 4). The analytical solutions for the interface position   t  , the solid (T s ) and liquid (T l ) temperatures are given by (Voller, 1997;Palle & Dantzig, 1996):  Fig. 3   (a, b). An excellent agreement is found between the present simulation and the analytical solution which in turn demonstrates the effectiveness of the proposed method. www.intechopen.com

2-D solidification problem
A two-dimensional (2-D) solidification problem for which analytical (Rathjen & Jiji, 1971) and numerical (LB) (Jiaung et al., 2001;Lin & Chen, 1997) solutions are available in the literature is now presented. Fig. 4 shows the schematic diagram of the problem with the boundary conditions. The material is kept initially at a uniform temperature T i which is higher than or equal to the melting temperature T m. The left (x = 0) and bottom (y = 0) boundaries are lowered to some fixed temperature   0 m TT  and consequently, solidification begins from these surfaces and proceeds into the material. Setting the scaled temperatures 0.3 i T  , 0 1 T  and 0 m T  as considered in Jiaung et al. (2001) and Rathjen & Jiji (1971) and assuming constant material properties, we obtain the LB simulation results following the proposed methodology. Fig. 5a

Melting of pure gallium
Melting of pure gallium in a rectangular cavity is a standard benchmark problem for validation of phase change modeling strategies, since reliable experiments in this regard (particularly, flow visualization and temperature measurements) have been welldocumented in the literature (Gau & Viskanta, 1986). Brent et al. (1988) solved this problem numerically with a first order finite volume scheme, coupled with an enthalpy-porosity approach, and observed an unicellular flow pattern, in consistency with experimental findings reported in Gau & Viskanta (1986), whereas Dantzig (1989) obtained a multicellular flow pattern, by employing a second order finite element enthalpy-porosity model. Miller et al. (2001), again, obtained a multicellular flow patterns while simulating the above problem, Fig. 4. Schematic of the two-dimensional solidification problem (Chatterjee, 2010) www.intechopen.com for the twodimensional solidification problem (Chatterjee, 2010) by employing a LB model in conjunction with the phase field method. In all the above cases, nature of the flow field was observed to be extremely sensitive to problem data employed for numerical simulations. Here, simulation results (Chakraborty & Chatterjee, 2007) are shown with the same set of physical and geometrical parameters, as adopted in Brent et al. (1988). The study essentially examines a two-dimensional melting of pure gallium in a rectangular cavity, initially kept at its melting temperature, with the top and bottom walls maintained as insulated. Melting initiates from the left wall with a small thermal disturbance, and continues to propagate towards the right. The characteristic physical parameters are as follows: Prandtl number (Pr) = 0.0216, Stefan number (St) = 0.039 and Raleigh number (Ra) = 6 × 10 5 . Numerical simulations are performed with a (56 × 40) uniform grid system, keeping the aspect ratio 1.4 in a 9 speed square lattice (D2Q9) over 6 × 10 5 time steps (corresponding to 1 min of physical time). The results show excellent agreements with the findings of Brent et al. (1988). For a visual appreciation of flow behavior during the melting process, Fig. 6 is plotted, which shows the streamlines and melt front location at time instants of 6, 10 and 19 min, respectively. The melting front remains virtually planar at initial times, as the natural convection field begins to develop. Subsequently, the natural convection intensifies enough to have a pronounced influence on overall energy transport in front of the heated wall. Morphology of the melt front is subsequently dictated by the fact that fluid rising at the heated wall travels across the cavity and impinges on the upper section of the solid front, thereby resulting in this area to melt back beyond the mean position of the front. After 19 min, the shape of the melting front is governed primarily by advection. Overall, a nice agreement can be seen between numerically obtained melt front positions reported in a benchmark study executed by Brent et al. (1988) and the present simulation. Slight discrepancies between the computed results Fig. 6. Melting of pure gallium in a rectangular cavity (Chakraborty & Chatterjee, 2007) (both in benchmark numerical work reported earlier and the present computations) and observed experimental findings (Gau & Viskanta, 1986) can be attributed to threedimensional effects in experimental apparatus to determine front locations, experimental uncertainties and variations in thermo-fluid properties. However, from a comparison of the calculated and experimental (Gau & Viskanta, 1986) melt fronts at different times (refer to Fig. 7), it is found that both the qualitative behavior and actual morphology of the experimental melt fronts are realistically manifested in the present numerical simulation.

Bridgman crystal growth
Results are presented for simulation of transport processes in a macroscopic solidification problem such as the Bridgman crystal growth in a square crucible (Chatterjee, 2010). The Bridgman crystal growth is a popular process for growing compound semiconductor crystals and this problem has been solved extensively as a benchmark problem. The typical problem domain along with the boundary condition is shown schematically in Fig. 8. Initially, the material is kept in a molten state at a temperature T i (= 1) higher than the melting point T m . Since initially there is no thermal gradient, consequently, there is no convection. At t = 0 + , the left, right and the bottom walls are set to the temperature T 0 , which is scaled to be zero, while the top wall is assumed to be insulated. This will lead to a new phase formation (solidification) at the walls with simultaneous melt convection. The characteristic physical parameters (arbitrary choice) for the problem are the Prandtl number  (Gau & Viskanta, 1986) results (dotted line) and continuum based numerical simulation (Brent et al., 1988) predictions (solid line) (Chakraborty & Chatterjee, 2007)  for the Bridgman crystal growth process (Chatterjee, 2010) the characteristic dimension of the simulation domain. Numerical simulations are performed on a (80 × 80) uniform grid systems with an aspect ratio of 1, in a 9 speed square lattice (D2Q9) over 6 × 10 5 time steps corresponding to 1 min of physical time. For a visual appreciation of the overall evolution of the transport quantities in this case, Fig. 9 is plotted, which shows the representative flow pattern and isotherms at a normalized time instant of 0.05 t  . The interval between the contour lines is 0.05 units (dimensionless). Larger isotherm spacing is observed in the melt which is a consequence of the heat of fusion released from the melt as well as a subsequent convection effect. The isotherms are normal to the top surface since the top surface is an adiabatic wall. Two counter rotating symmetric cells are observed in the flow pattern which is consistent with the flow physics. The melt convection will become weaker as the solidification progresses since there is very little space for convection. Also the thermal gradient will become small at this juncture. The calculation continues until the melt completely disappears and the temperature of the entire domain eventually reaches T 0 . In order to demonstrate the capability of the proposed method in capturing the interfacial region without further grid refinement as normally required for the phase field based method or any other adaptive methods, Fig. 10 is plotted in which the comparison of the isotherm obtained from the present simulation for the Bridgman crystal growth and from an adaptive finite volume method (Lan et al., 2002) is shown. Virtually there is no deviation of Fig. 10. Comparison of isotherm from the present calculation (solid lines) and from an adaptive finite volume method (dashed lines) (Lan et al., 2002) (Chatterjee, 2010) www.intechopen.com the calculated isotherm form that obtained from the adaptive finite volume method (Lan et al., 2002) has been observed. This proves that the present method is quiet capable of capturing the interfacial region without further grid resolution.

Crystal growth during solidification
In this section, the problem of crystal growth during solidification of an undercooled melt is discussed (Chatterjee & Chakraborty, 2006). Special care is taken to model the effects of curvature undercooling, anisotropy of surface energy at the interface and the influence of thermal noise, borrowing principles from cellular automaton based dendritic growth models (Sasikumar & Sreenivasan, 1994;Sasikumar & Jacob, 1996), in the framework of a generalized enthalpy updating scheme adopted here. Numerical experiments are performed to study the effect of melt convection on equiaxed dendrite growth. Since flow due to natural convection (present in a macroscopic domain) can be simulated as a forced flow over microscopic scales, a uniform flow is introduced through one side of the computational domain, and its effect on dendrite growth morphology is investigated. Computations are carried out in a square domain (50 × 50 uniform grid-system) containing initially a seed crystal at the center, while the remaining portion of the domain is filled with a supercooled melt. The physical parameters come from the following normalization of length (W) and time (τ) units: , where δ is the interfacial length scale (typically O (10 -9 m)), μ k is a kinetic coefficient (typically O (10 -1 m/s.K)) and g  is the Gibbs-Thompson coefficient (typically O (10 -7 m.K)). Exact values of the above parameters have been taken from Beckermann et al. (1999). The degree of undercooling corresponds to 0.515 K. Fig. 11 demonstrates the computed evolution of dendritic arms under the above conditions. In absence of fluid flow, the dendrite arms grow in an identical manner  (Chatterjee & Chakraborty, 2006) ( Fig. 11a), simply because of isotropic heat extraction through all four boundaries. Fig. 11(b) illustrates the effect of convection on the above dendritic growth. In the upstream side (top), convection opposes heat diffusion, which subsequently reduces the thermal boundary layer thickness and increases local temperature gradients, eventually, leading to a faster growth of the upper dendritic arm. Evolution of the downstream arm (bottom), on the other hand, is relatively retarded, for identical reasons. For a more comprehensive validation of the quantitative capabilities of the present LB model to simulate dendritic growth in presence of fluid flow, results predicted by the present model are compared with those reported in Beckermann et al. (1999), and a visual appreciation of the same is depicted in Fig. 12. It is revealed from Fig. 12 that the solid fraction contours and isotherms based on the present model match excellently with the dendritic envelopes depicted in Beckermann et al. (1999). These results, further, indicate excellent convergence properties of the present LB based method, over a wide range of Reynolds and Prandtl numbers.  (Chatterjee & Chakraborty, 2006) www.intechopen.com

Summary
This chapter briefly summarizes the development of a passive scalar based thermal LB model to simulate the transport processes during melting/freezing of pure substances. The model incorporates the macroscopic phase changing aspects in an elegant and straightforward manner into the LB equations. Although the model is developed for twodimensional phase change problems, it can be easily extended to three-dimension. These features make the model attractive for simulating generalized convection-diffusion melting/solidification problems. Because of its inherent simplicity in implementation, stability, accuracy, as well as its parallel nature, the proposed method might be a potentially powerful tool for solving complex phase change problems in physics and engineering, characterized by complicated interfacial topologies. Compared with the phase field based LB models, the present scheme is much simpler to implement, since extremely refined meshes are not required here to resolve a minimum length scale over the interfacial regions. Although a finer mesh would definitely results in a better-resolved interface and a more accurate capturing of gradients of field variables, the mesh size for the present model merely plays the role of a synthetic microscope to visualize topological features of the interface morphology.

Acknowledgment
The author gratefully acknowledge Dr. The constant evolution of the calculation capacity of the modern computers implies in a permanent effort to adjust the existing numerical codes, or to create new codes following new points of view, aiming to adequately simulate fluid flows and the related transport of physical properties. Additionally, the continuous improving of laboratory devices and equipment, which allow to record and measure fluid flows with a higher degree of details, induces to elaborate specific experiments, in order to shed light in unsolved aspects of the phenomena related to these flows. This volume presents conclusions about different aspects of calculated and observed flows, discussing the tools used in the analyses. It contains eighteen chapters, organized in four sections: 1) Smoothed Spheres, 2) Models and Codes in Fluid Dynamics, 3) Complex Hydraulic Engineering Applications, 4) Hydrodynamics and Heat/Mass Transfer. The chapters present results directed to the optimization of the methods and tools of Hydrodynamics.