A Comparative Study of Multiphase Lattice Boltzmann Methods for Bubble-Dendrite Interaction during Solidiﬁcation of Alloys

: This paper presents a comparative study between the pseudopotential Shan-Chen model and the phase ﬁeld multiphase lattice Boltzmann method for simulating bubble dynamics during dendritic solidiﬁcation of binary alloys. The Shan-Chen method is an efﬁcient lattice Boltzmann multiphase method despite having some limitations, including the generation of large spurious currents. The phase ﬁeld model solves the Cahn-Hilliard equation in addition to the Navier-Stokes equation to track the interface between phases. The phase ﬁeld method is more accurate than the Shan-Chen model for simulation of ﬂuids with a high-density ratio since it generates an acceptable small spurious current, though at the expense of higher computational costs. For the simulations in this article, the multiphase lattice Boltzmann model was coupled with the cellular automata and ﬁnite difference methods to solve temperature and concentration ﬁelds. The simulated results were presented and compared regarding the ability of each model to simulate phenomena at a microscale resolution, such as Marangoni convection, the magnitude of spurious current, and the computational costs. It is shown that although Shan-Chen methods can replicate some qualitative features of bubble-dendrite interaction, the generated spurious current is unacceptably large, particularly for practical values of the density ratio between ﬂuid and gas phases. This occurs even after implementation of several enhancements to the original Shan-Chen method. This serious limitation makes the Shan-Chen models unsuitable to simulate ﬂuid ﬂow phenomena, such as Marangoni convection, because the large spurious currents mask completely the physical ﬂow.


Introduction
Formation of micro defects during solidification processes has gained the attention of many researchers since they affect the mechanical properties of solidified components [1][2][3]. Numerous attempts have been made to minimize porosity in castings [4][5][6]. For example, understanding bubble formation and motion during dendritic solidification would help to improve the quality of cast products [7].
In-situ observation provides useful information about the bubble formation and morphology during solidification. Several methods are available to observe porosity during solidification, such as optical methods [8,9] and micro focus X-ray imaging [10][11][12]. Generally, in-situ experiments suffer from limitations regarding tracking the bubble motion, transparency of alloys, and the thickness of sample [13].
Numerical simulations offer an alternative to investigate the mechanism behind bubble formation and their interactions during solidification of metallic alloys. Karagadde et al. [14] proposed a hybrid level set enthalpy-based model to simulate the evolution of hydrogen bubbles during the solidification of aluminum alloys in 2D. They estimated the final pore shape based on the pore radius and cooling rate. Tiedje et al. [15] simulated the evolution of porosity for Al-Si. They identified three zones in the domain. In zone one, a small-sized porosity was observed. In zone two, which they called the transition zone, the porosity was elongated. In zone three or the central zone, the pores are rounded. Du et al. [16] simulated the dendritic growth in pure aluminum with pre-existing pores using the phase field method. They presented the effects of pressure on gas bubble formation, nucleation, and evolution. They modeled and discussed the interaction between bubbles and the solid-liquid interface. They showed that the presence of gas bubbles in solid-liquid interfaces could change the dendritic growth significantly.
One of the challenges in modeling bubble dynamics during dendritic growth is how to simulate multiphase flow. Bubble dynamics simulations are often very computationally demanding. One reason for the cost is the difficulty regarding tracking the interface between different phases. Conventional computational fluid dynamics (CFD) models to solve multiphase flow are divided into two main categories: Interface tracking methods and interface capturing methods.
In interface tracking methods, the location of the interface is marked with a separate grid or set of surface meshes from an initial condition and calculated explicitly during the simulation. The boundary integral methods [17], the arbitrary Lagrangian-Eulerian (ALE) method [18], and all front tracking methods [19][20][21][22][23][24] belong to this category.
In interface capturing methods, instead of tracking a sharp interface, a function represents the thin, but non-zero, interface between phases. This leads to the fluid properties changing continuously from one fluid to another across the interface. The volume of fluid (VOF), phase field, and level set methods belong to this category [25][26][27].
As an alternative to more conventional methods, the lattice Boltzmann (LB) method is successfully used in many different multiphase problems [28,29]. The LB method, with its roots in kinetic theory, is based on a solution of the Boltzmann kinetic equation for a group of imaginary particles in a discretized domain, which can recover the Navier-Stokes equation. Therefore, LB retrieves macroscopic properties of the fluid while making use of the microscopic method. Due to its simplicity, local structure, ease in dealing with complex geometries, explicit nature, and suitability to use in distributed memory architectures, it has been applied to simulate various physical phenomena [13,28,30,31].
In the color gradient model, each phase is assigned a color, and for each color, a distribution function is introduced to calculate the interaction between different phases. The interface is determined using the color gradient. In the free energy model, two particle distribution functions are utilized. One distribution function calculates an order parameter, which defines each phase, and the other distribution function is responsible for the predicted velocity without considering a pressure gradient. To find the exact pressure and velocity, another Poisson-type equation must be solved in each iteration. In the phase field method, in addition to the fluid flow equation, another equation is solved to determine the interface location. The fluid flow distribution function for the local density and momentum is transformed to the mean density and momentum to improve stability, and then a cohesion force is introduced in the mean field flow. The He-Chen-Zhang (HCZ) model is a variation of the phase field-based multiphase method, in which two sets of distribution functions are used to recover the Navier-Stokes (NS) and Cahn-Hilliard (CH) equations [29,36,37,[39][40][41].
The Shan-Chen model and all of its variations are among the most used models because the interface is not introduced as a boundary condition and no extra equation is required to track the interface between different phases [29]. However, the original Shan-Chen model has several issues, including, but not limited to: Large spurious current, limitations regarding simulating fluids with a high density and viscosity ratio, and avoidance of coalescences of all bubbles [28]. Chen et al. [28] classified the ways to improve the original Shan-Chen model. The improvements consist of the following: Using a realistic equation of state [42], using an interaction force with a higher order of isotropy [43], and introducing a mid-range repulsive force [44].
The phase field LB model has been successfully used to simulate many problems of multiphase flow [45][46][47][48]. Most of the researchers have solved the CH equation by another set of LB distribution functions. However, solving the CH equation by the LB technique has drawbacks that may compromise the results. First, selecting the relaxation parameter for the CH equation is not apparent. Many articles use the same fluid flow relaxation time for the CH equation. This selection is not physically acceptable since different values of relaxation time may lead to different results. Next, the LB model cannot entirely recover the CH equation, but the recovered CH equation always has higher order terms. It is also hard in the LB model to discretize the convection term in the CH equation by upwind schemes since the streaming process implies a central discretization. Finally, it is hard to use higher order explicit time stepping for the LB method since LB models usually use second-order explicit time stepping, which makes the time step size significantly small for the CH equation [30,41,49]. In this article, the CH equation is solved by a weighted essentially non-oscillatory (WENO) scheme for the convection term and a third order Runge-Kutta time stepping scheme for the time marching.
There are few articles in which the cellular automaton (CA) and LB methods were coupled to simulate bubble dynamics during dendrite growth; two of them are Refs. [7,50]. In these cases, the model predicted gas-liquid-solid interactions correctly and the evolution of bubbles during solidification. However, in both articles, the original Shan-Chen model was used, hence suffering from all the problems mentioned earlier. Most importantly, they did not discuss the high spurious current generated by the Shan-Chen model. This spurious current can compromise the results since it is hard to distinguish between real and spurious velocity. Another problem is that their studies were limited to a low density ratio of 10. Therefore, it is still a long way until current Shan-Chen models can simulate phenomena in a microscale resolution with real physical properties.
The purpose of this paper is to compare Shan-Chen-based and phase-field-based LB models to conduct a feasibility study of each model for simulating phenomena in the microscale, such as Marangoni convection. First, various enhancements available for the Shan-Chen LB scheme for modeling multiphase flow were applied to build an improved version of the Shan-Chen based model to study bubble dynamics during dendritic growth. These enhancements have been utilized separately before, but to the best of the authors' knowledge, they have never been combined for use in simulations. By applying these improvements, the enhanced model can handle a larger density ratio while maintaining a low spurious current. The phase field LB model was also developed to compare between the results. An interface-tracking CA algorithm and an LB model for transport equations were used to simulate and track the dendrite growth, while the finite difference method was applied to solve the heat transfer equation for all models. After validating the Shan-Chen based and phase field models against the published literature, the different variations of multiphase models developed in this article were utilized to simulate bubble dynamics during dendritic solidification of a binary alloy.

Cellular Automaton Model for Dendritic Growth
The CA scheme was used to track the solid-liquid interface, as explained in [13,[51][52][53]. In this model, dendritic growth is driven by the difference between the local interface equilibrium solute concentration, C * l , and local actual solute concentration, C l : where D is the solute diffusivity and ∆f s is the increase in the solid fraction in each node located at the interface. The last term in Equation (2) represents the solute rejected to the interface due to the solute partitioning taking place between the solid and liquid. ∆f s depends on the actual liquid concentration, C l , and the local interface equilibrium solute concentration, C * l , and can be obtained by: where k is the partitioning coefficient. The growth angle, χ, is calculated as: Based on the solid fraction, fs, the interface curvature of a cell is obtained by:

Single Phase Lattice Boltzmann Equation
In LB, the fluid characteristics are described by sets of imaginary particles by propagation and relaxation in the d-dimensional lattice. For instance, in the D 2 Q 9 lattice, which was used in this work, each node at position x has a distribution function, f α (x, t), in any of the nine discrete directions, e α . Each node in each iteration not only passes information to itself, but also communicates with its eight neighboring nodes. The LB equation with the Bhatnagar-Gross-Krook (BGK) collision model in each iteration is written as: where F α (x, t) is the force exerted on each node. τ υ = 1 ω υ is the relaxation time and is related to the local kinematic viscosity (in LB units) as ν LB = c 2 s (τ v −0 .5). For D 2 Q 9 lattice, the discrete direction, e α , is defined as:  The equilibrium distribution function in the D 2 Q 9 model is defined as [54]: where ω α is a weight coefficient defined as:

The Original Shan-Chen Model for Multiphase Flow
The original Shan-Chen model [34] introduces a cohesion force term, F(x, t), responsible for phase separation. This cohesion force satisfies the non-ideal equation of state. Based on the pressure difference, the phase separation between different phases takes place: The parameter, ψ, is called the effective mass and is related to the equation of state. Here, G 1 controls the force exerted at each node by surrounding particles, with a positive (negative) value leading to a repulsive (attractive) force between particles. In this model, the phase separation occurs when G 1 is higher than a critical value. Using this scheme, only eight neighbor nodes are considered for calculating the cohesion force in the D 2 Q 9 lattice.
Shan and Chen [34] proposed the following equation for the effective mass: Other researchers [7] have suggested different equations for effective mass as: , or even as ψ(ρ) = ρ 0 . In this article, Equation (11) was used for calculating effective mass, where ψ 0 and ρ 0 are constant. The pressure term in this model has an extra term compared to the ideal gas equation of state. The pressure term is defined as: To model the contact angle, the scheme developed by Benzi et al. [55] was employed. They suggested that adhesion force can be implemented in a similar way to cohesion force as follows: Here, s(x + e α ∆t, t) is an indication parameter and has the value of 0 or 1 for fluid or solid nodes, respectively. Ψ(ρ w ) is the effective mass at the wall; by changing ψ(ρ w ), different contact angles can be achieved.
Depending on the node type (solid or liquid), Equation (10) or (13) was used to simulate the interaction between nodes.
The fluid velocity is calculated as: (14) As mentioned in the introduction, the original Shan-Chen model generates a high spurious current that can mask other important phenomena. Another issue with this model is the unrealistic rate of Appl. Sci. 2019, 9, 57 6 of 24 coalescence between bubbles (drops) that ultimately results in a single bubble (drop) remaining in the whole domain. These drawbacks motivated us to look for a more robust Shan-Chen model for dendritic solidification simulations.

Realistic Equation of State
Based on Yuan and Schaefer [42], to model a higher density ratio and reduce the spurious current at the same time, a different equation of state (EOS) is utilized. Here, we use the Carnahan-Starling (C-S) EOS since it is stable, easy to implement, and generates a lower spurious current compared to other EOS [42].
No matter what EOS is used, the effective mass is expressed as: The C-S EOS is expressed as: where a = 0.4963(RT cr ) 2 p cr , b = 0.1873RT cr /p cr . P cr and T cr are the critical pressure and temperature, respectively. Without losing the generality of the model, we assume a = 1, b = 4, R = 1 and by reducing the temperature in this equation, a higher density ratio is achieved [28,29].

Force with a Higher Order of Isotropy (E8 Force Scheme) and Middle-Range Repulsion Force
In the original Shan-Chen method, to compute the cohesion force, only eight neighbor nodes are considered. However, the cohesion force can include any number of neighbor nodes. By communicating only with the eight neighboring nodes in the D2Q9 lattice, the highest isotropy order that can be achieved is four (E4 force scheme). By considering the second layer (24 neighbor nodes), tensors of the eighth order can be produced (E8 force scheme) [28,43].
So, for the E8 force scheme, instead of using Equation (10), the force term in the D 2 Q 9 lattice is expressed by: where G 2 is a negative coefficient representing the attractive force, and ω γ is the weighting coefficient. The direction of e γ is shown in Figure 1 and values for the weighting coefficients for the E8 force scheme are provided in Equation (18).
Chibbaro et al. [44] introduced a mid-range repulsion force between fluid particles. In the original Shan-Chen model, the attractive force causes the phase separation. Without adding any repulsive force, after some iterations, only one big bubble (drop) remains in the domain.
Basically, in Ref. [44], the authors introduced G1 and G2 parameters to model cohesion force. The G1 is negative and represents the attractive force of the first belt of neighboring nodes, while G2 is positive and represents the repulsive force of the second belt of neighboring nodes. The negative or positive signs here are, respectively, corresponding to attractive and repulsive forces between fluids.
Therefore, Equation (10) is modified as: The pressure and effective mass for both cases are defined as: and: Chibbaro et al. [44] introduced a mid-range repulsion force between fluid particles. In the original Shan-Chen model, the attractive force causes the phase separation. Without adding any repulsive force, after some iterations, only one big bubble (drop) remains in the domain.
Basically, in Ref. [44], the authors introduced G 1 and G 2 parameters to model cohesion force. The G 1 is negative and represents the attractive force of the first belt of neighboring nodes, while G 2 is positive and represents the repulsive force of the second belt of neighboring nodes. The negative or positive signs here are, respectively, corresponding to attractive and repulsive forces between fluids.
Therefore, Equation (10) is modified as: The pressure and effective mass for both cases are defined as: and: Appl. Sci. 2019, 9, 57 8 of 24 Note that in the case of the E8 force scheme without mid-range repulsion, G 1 = 0 and G 2 < 0. The weighting coefficient and the direction definitions are identical to Equation (18) and Figure 1. By introducing this midrange repulsive force, phenomena, such as spray formation, and the soft-glassy system can be modeled, which is impossible to model by short-range attraction force alone as it is mentioned in Ref [44].
In this article, the realistic EOS and E8 force scheme with midrange repulsion force were implemented to improve the Shan-Chen model for bubble-dendrite interactions during alloy solidification. To be succinct, this model will be referred to as "the enhanced model".

Phase Field Lattice Boltzmann Method
In this section, the phase field LB method is presented. The fluid flow is simulated by the model of Shao et al. [56], which uses the mean density of two phases in the distribution function of the LB method. The interface between phases is captured by solving the CH equation.

Solving the Cahn-Hillard Equation for Interface Capturing with the WENO Scheme
In the phase field model, a local order parameter, ∅, is defined to differentiate between two phases, i.e., ∅ l and ∅ g , referring to the liquid and gas phase, respectively. The local order parameter controls the kinetics and evolution of each phase throughout the simulation. In the CH equation, the motion of the interface is defined as: This equation contains terms related to convection and diffusion of the interface. λ is a diffusion coefficient called mobility. µ φ is the chemical potential, which is defined as the derivative of the free energy (Ψ φ ) with respect to the order parameter as µ ∅ = Ψ φ (∅) − κ∇ 2 ∅. If the free energy takes the double well form as Ψ φ = β∅ 2 (∅ − 1) 2 , then the chemical potential becomes 4β∅(∅ − 0.5)(∅ − 1) − κ∇ 2 ∅, where both β and κ are related to the surface tension, σ, and interface width, W, as β = 12σ W(ρ l −ρ g ) 4 and κ = 1.5 Wσ (ρ l −ρ g ) 2 . For multiphase flow with large density ratios between phases, the order parameter distribution across the interface may show a high gradient, which causes instability in the simulation. To overcome this problem, it is preferred to use an upstream condition for the advection term. In this article, the upwind WENO scheme with third order Runge Kutta total variation diminishing (TVD) [57] is implemented to discretize the convection and temporal term, respectively. Equation (22) can then be rewritten as: The term, u∅, in Equation (22) is labeled as l for simplification. The discretization of the first term in the right-hand side of Equation (23) can be achieved in a variety of ways. In this paper, we consider the Lax-Friedrichs flux splitting, which uses three stencils, formed by five points [58]. In this method, the convection term is calculated by flux terms, which are determined by: where m is selected as max|(u∅) | to make the scheme stable. Since the discretization in the any coordinate direction follows the same procedure, the discretization in the x direction is going to be explained. By splitting the flux term, l, into positive and negative terms, we have: are approximated by five points as: And: The stencil weights can be calculated as: where: The parameter, IS k , is defined as: And: In this way, l can be calculated in the same way.
Next, we focus on the diffusion term in Equation (23). By considering the constant mobility, the diffusion term can be rewritten as: The second order derivative in the diffusion term can be expressed as: For marching in time, the third order TDV Runge-Kutta scheme is applied as: where ∅ t is ∅ value at time t (∆t is the time interval), and ∅ (1) and ∅ (2) are the transitional values in each time step. With a stable spatial and temporal discretization of the CH equation, the evolution of the interface for multiphase flow with a large density ratio can be achieved [41].

Phase Field LB for the Flow Field
Here, we briefly introduce the flow part of the model presented by Zheng et al. [59]. In this model, the LB equation for fluid flow is written as: where the chemical potential, µ ∅ , was defined before. Γ α (u) is given as: The −φ∇µ ∅ term is the interfacial force between phases where φ is the local order parameter. The equilibrium distribution function can be expressed as: In Equation (37), the parameters, ρ 0 and ρ, represent the mean and local density at each point. The mean density is initialized as (ρ l +ρ g ) 2 where l and g refer to the liquid and gas phase. The macroscopic properties can be obtained by: The relationship between the local density and local relaxation parameter at each point with the local order parameter are defined as: It is known that through the Chapman-Enskog expansion analysis, the following macroscopic equation can be recovered [56]: Now that every method in this article introduced, the flow-chart of simulation is illustrated in Figure 2. The multiphase solver can be Shan-Chen based or the phase field model.

Results
In this section, after validating the enhanced model by a phase separation problem in a periodic domain, and the phase field model with the Rayleigh instability problem, the simulations of dendritebubble interaction are presented with the original Shan-Chen model, the enhanced model, and the phase field model.

Rayleigh Instability
The phase field LB model was validated by simulating the Rayleigh instability problem with different density ratios. The Rayleigh instability is defined when a denser fluid with density, ρ l , is placed over a less dense fluid with density, ρ g , in the domain. Here, the domain size was set as L*4L, where L is the number of nodes in the x-direction. The dimensionless At number, which represents the density ratio between phases, was defined as At = (ρ l −ρ g ) (ρ l +ρ g ) , the characteristic velocity was assigned as U = √gL, and the Reynold number was Re = UL υ . Periodic boundary conditions were imposed on the left and right boundaries while the wall boundary condition with bounce back was applied to the bottom and top boundaries. With 10% initial perturbation, the initial interface shape was y = 2L + 0.1Lcos ( 2πx L ) [41].
For the first case, the domain size was selected as 128 × 512 and the Reynolds number was Re = 256. The At number was set as At = 0.5, which corresponds to a density ratio of three. A characteristic velocity of U = 0.04 in the LB unit was adopted, which makes the kinematic viscosity, υ = 0.02. All the results are reported in nondimensional time, which was normalized by √ L g . Figure 3 shows the location of the interface during the simulation. The results are in good agreement with the benchmark solution reported in Ref. [46].

Results
In this section, after validating the enhanced model by a phase separation problem in a periodic domain, and the phase field model with the Rayleigh instability problem, the simulations of dendrite-bubble interaction are presented with the original Shan-Chen model, the enhanced model, and the phase field model.

Rayleigh Instability
The phase field LB model was validated by simulating the Rayleigh instability problem with different density ratios. The Rayleigh instability is defined when a denser fluid with density, ρ l , is placed over a less dense fluid with density, ρ g , in the domain. Here, the domain size was set as L*4L, where L is the number of nodes in the x-direction. The dimensionless At number, which represents the density ratio between phases, was defined as At = (ρ l −ρ g ) (ρ l +ρ g ) , the characteristic velocity was assigned as U = gL, and the Reynold number was Re = UL υ . Periodic boundary conditions were imposed on the left and right boundaries while the wall boundary condition with bounce back was applied to the bottom and top boundaries. With 10% initial perturbation, the initial interface shape was y = 2L + 0.1Lcos 2πx L [41]. For the first case, the domain size was selected as 128 × 512 and the Reynolds number was Re = 256. The At number was set as At = 0.5, which corresponds to a density ratio of three. A characteristic velocity of U = 0.04 in the LB unit was adopted, which makes the kinematic viscosity, υ = 0.02. All the results are reported in nondimensional time, which was normalized by L g . Figure 3 shows the location of the interface during the simulation. The results are in good agreement with the benchmark solution reported in Ref. [46]. To validate our results quantitively, the spike tip and bubble front locations are shown in Figure  4. These locations are in excellent agreement with the results reported by He et al. [37].  Figure 3 and comparison with He et al. [37] results. Axis variables are non-dimensional.
The model was also validated for the same problem, but with a density ratio of 1000, which corresponds to At = 0.999. The Reynolds number was kept at Re = 256. The evolution of the interface is shown in Figure 5, which is consistent with the results of Figure 6 in Ref [41]. The results also show the dependency of the interface shape with the density ratio. For a small density ratio, the tendency of mixing between the two fluids was large. While for a larger density ratio, a more coherent structure between the two fluids was observed. To validate our results quantitively, the spike tip and bubble front locations are shown in Figure 4. These locations are in excellent agreement with the results reported by He et al. [37]. To validate our results quantitively, the spike tip and bubble front locations are shown in Figure  4. These locations are in excellent agreement with the results reported by He et al. [37].  Figure 3 and comparison with He et al. [37] results. Axis variables are non-dimensional.
The model was also validated for the same problem, but with a density ratio of 1000, which corresponds to At = 0.999. The Reynolds number was kept at Re = 256. The evolution of the interface is shown in Figure 5, which is consistent with the results of Figure 6 in Ref [41]. The results also show the dependency of the interface shape with the density ratio. For a small density ratio, the tendency of mixing between the two fluids was large. While for a larger density ratio, a more coherent structure between the two fluids was observed.  Figure 3 and comparison with He et al. [37] results. Axis variables are non-dimensional.
The model was also validated for the same problem, but with a density ratio of 1000, which corresponds to At = 0.999. The Reynolds number was kept at Re = 256. The evolution of the interface is shown in Figure 5, which is consistent with the results of Figure 6 in Ref [41]. The results also show the dependency of the interface shape with the density ratio. For a small density ratio, the tendency of mixing between the two fluids was large. While for a larger density ratio, a more coherent structure between the two fluids was observed.

The Phase Separation Problem with the Enhanced Model
The phase separation between the liquid and vapor phases of an initial mixture of both phases was studied in this section. The domain was initially at rest with an average density plus a random variation in in all nodes. All boundary conditions were periodic. The system was unstable based on EOS, and phase separation occurred. The final shape of bubbles or drops was circular since the free energy of the domain tends to minimize, and a circle has the minimum surface area compared to other shapes.

The Phase Separation Problem with the Enhanced Model
The phase separation between the liquid and vapor phases of an initial mixture of both phases was studied in this section. The domain was initially at rest with an average density plus a random variation in in all nodes. All boundary conditions were periodic. The system was unstable based on EOS, and phase separation occurred. The final shape of bubbles or drops was circular since the free energy of the domain tends to minimize, and a circle has the minimum surface area compared to other shapes.  Table 1 at T = 0.9T cr . (a) case 1, (b) case 2, (c) case 3, (d) case 4, (e) case 5, (f) case 6. (blue: gas, red: liquid).

The Phase Separation Problem with the Enhanced Model
The phase separation between the liquid and vapor phases of an initial mixture of both phases was studied in this section. The domain was initially at rest with an average density plus a random variation in in all nodes. All boundary conditions were periodic. The system was unstable based on EOS, and phase separation occurred. The final shape of bubbles or drops was circular since the free energy of the domain tends to minimize, and a circle has the minimum surface area compared to other shapes.
The problem was simulated by the enhanced model in six different cases. In each case, the initial density of the mixture was different, and a random variation in the order of 0.01 was added to the density. The temperature was selected as T = 0.9T cr which corresponds to a density ratio of 6.
Based on the value of G 1 and G 2 in Table 1, different scenarios of phase separation occurred as shown in Figure 6. For larger G 2 , corresponding to a bigger repulsive force, the coalescence of the drops was prevented in some of the cases. These results are in agreement with Ref. [44], which simulated the same problem. One difference with the simulations of Ref. [44] is that they used the original Shan-Chen EOS while the more realistic C-S EOS was used in this work. However, since the density ratio was almost the same in this study as the one in Ref. [44] (six in this study and five in Ref. [44]), the results can be compared directly. The results show that the midrange-repulsive force can be used to prevent the coalescences of drops, but not bubbles. In the case of bubbles, this method only delays the merging process, but, eventually, just one bubble remains in the domain. For case (a) and (b), one bubble was present in the domain. However, for other cases, the formation of drops was observed. The formation of bubble or drops was related to the initial density. If the initial density is larger than a critical value, after some iterations, only one bubble will result in the whole domain.
As mentioned before, the enhanced model reduced the spurious current, but even this reduced amount of spurious current was too large for modeling phenomena, such as Marangoni convection, in this scale. The maximum amount of spurious current was observed in the interface between the bubbles and liquid, which is the place where the Marangoni force is exerted. The magnitude of the spurious current is presented in Table 2, showing approximately one order of magnitude reduction of the spurious current in the enhanced Shan-Chen model with respect to the original one, but still a large value at 5 × 10 −3 . On the other hand, when this same problem was simulated with the phase field lattice Boltzmann method, the generated spurious current was in the order of 10 −6 , which is acceptable for most cases of buoyancy and Marangoni flows [31]. The total computational time of the simulation was also a critical parameter. The simulation time was non-dimensionalized with the time needed to solve the original Shan-Chen model. In Table 2, it is observed that the model enhancements make it about 250% more expensive in CPU time than the original Shan-Chen model.

Columnar Dendritic Growth with Original Shan-Chen Model
The growth of five columnar dendrites in a rectangular domain was simulated using the original Shan-Chen model. The simulation was performed for Al-3wt% Cu binary alloy with initial undercooling of ∆T = 2.0 K. All the external boundary conditions were considered as stationary walls for the velocity field and insulated for the temperature field except for the left wall where a temperature gradient of 1200 K.m −1 was imposed. The boundary condition for the concentration field was assumed as insulated in all directions. The preferential crystallographic orientation was 0 • with respect to the x-axis (horizontal). The domain size was 91.2 µm × 72 µm, discretized with 304 × 240 cells for solving fluid flow and solute transport equations, and 38 × 30 cells for solving the energy equation. Five columnar dendrites and 10 bubbles with a radius between 4 µm to 6 µm were initially placed in the domain (10). The contact angle was selected as 150 • . Other parameters are shown in Table 3. Since the focus of this section is to compare the ability of the original and enhanced Shan-Chen models to simulate dendrite-bubble interaction, the gas component rejection into the liquid phase during solidification was not considered in these simulations. The initial condition is depicted in Figure 7. was assumed as insulated in all directions. The preferential crystallographic orientation was 0° with respect to the x-axis (horizontal). The domain size was 91.2 μm × 72 μm, discretized with 304 × 240 cells for solving fluid flow and solute transport equations, and 38 × 30 cells for solving the energy equation. Five columnar dendrites and 10 bubbles with a radius between 4 μm to 6 μm were initially placed in the domain (10). The contact angle was selected as 150°. Other parameters are shown in Table 3. Since the focus of this section is to compare the ability of the original and enhanced Shan-Chen models to simulate dendrite-bubble interaction, the gas component rejection into the liquid phase during solidification was not considered in these simulations. The initial condition is depicted in Figure 7. The interaction between dendrites and bubbles is shown in Figure 8. During the simulation, some bubbles tended to dissolve in the fluid or merge to other bubbles. This process made some bubbles bigger while eliminating the rest. The reason for this phenomenon is that the pressure of the bulk fluid was greater than the corresponding saturated value for the smaller bubbles. Therefore, these bubbles condensed. On the other hand, the pressure of the bulk fluid was lower than the  Figure 8. During the simulation, some bubbles tended to dissolve in the fluid or merge to other bubbles. This process made some bubbles bigger while eliminating the rest. The reason for this phenomenon is that the pressure of the bulk fluid was greater than the corresponding saturated value for the smaller bubbles. Therefore, these bubbles condensed. On the other hand, the pressure of the bulk fluid was lower than the saturated value of larger bubbles, causing the larger bubbles to grow. This is related to the EOS and the fact that no repulsive force was present in the domain. Therefore, during the simulations, some bubbles became smaller or vanished, while others became larger. The rate of bubble merging obtained with the original Shan-Chen model is not realistic as evidenced by experimental results. Coalescence of bubbles occurred in the experiments, but in a different time scale compared to what is presented in Figure 8. Based on the simulation results, all the small bubbles in the middle of the domain coalesced after 0.005 s while in the experiment, this phenomenon happened in the order of 10 s [60,61]. This model also generates a high spurious current, especially in the interface. The spurious current was in the order of 5 mm/s, which is very high and can completely mask the actual physical flow.

Dendritic Growth with the Enhanced Shan-Chen Model
The same problem of the previous section was modeled with the enhanced version of the Shan-Chen model. Since the enhanced model allows the simulation of higher density ratios, a density ratio of 40 was used for the results shown in Figure 9. The enhanced model also has the ability to control the coalescence of bubbles. The values used for the G parameters were G1 = −1.4 and G2 = 1.
(a) (b) The rate of bubble merging obtained with the original Shan-Chen model is not realistic as evidenced by experimental results. Coalescence of bubbles occurred in the experiments, but in a different time scale compared to what is presented in Figure 8. Based on the simulation results, all the small bubbles in the middle of the domain coalesced after 0.005 s while in the experiment, this phenomenon happened in the order of 10 s [60,61]. This model also generates a high spurious current, especially in the interface. The spurious current was in the order of 5 mm/s, which is very high and can completely mask the actual physical flow.

Dendritic Growth with the Enhanced Shan-Chen Model
The same problem of the previous section was modeled with the enhanced version of the Shan-Chen model. Since the enhanced model allows the simulation of higher density ratios, a density ratio of 40 was used for the results shown in Figure 9. The enhanced model also has the ability to control the coalescence of bubbles. The values used for the G parameters were G 1 = −1.4 and G 2 = 1.

Dendritic Growth with the Enhanced Shan-Chen Model
The same problem of the previous section was modeled with the enhanced version of the Shan-Chen model. Since the enhanced model allows the simulation of higher density ratios, a density ratio of 40 was used for the results shown in Figure 9. The enhanced model also has the ability to control the coalescence of bubbles. The values used for the G parameters were G1 = −1.4 and G2 = 1. The effect of the contact angle can be observed through the shape of bubbles that were in contact with dendrites. The magnitude of the spurious current in this simulation was about 3.5 mm/s, similar to the one obtained with the original model, but for a density ratio 10 times smaller. As observed in the original Shan-Chen model, dissolution of small bubbles and growth of large bubbles was observed here as well. However, the rate of merging could be controlled by the parameters, G1 and G2, giving the model the ability to reproduce physical results. The enhanced model clearly improved the results by reducing the spurious current and avoiding all bubbles coalescing together. The spurious current acted as a barrier for simulating higher density ratios. The enhanced model handled higher density ratios than the original model with a similar magnitude of spurious current. However, the spurious current was still too large in both models. Since the velocity and concentration fields were coupled, the high spurious current affected the concentration field and, consequently, the morphology of the dendrites and reliability of the results.

Dendritic Growth with the Phase Field Model
The same dendrite growth problem, as defined in Section 7.4, was solved with the phase field LB model. For comparison purposes, the density ratio was selected as 40 (ρ l = 1 and ρ g = 0.025) even though the phase field model can simulate much higher density ratios. The results are shown in Figure 10 The effect of the contact angle can be observed through the shape of bubbles that were in contact with dendrites. The magnitude of the spurious current in this simulation was about 3.5 mm/s, similar to the one obtained with the original model, but for a density ratio 10 times smaller. As observed in the original Shan-Chen model, dissolution of small bubbles and growth of large bubbles was observed here as well. However, the rate of merging could be controlled by the parameters, G 1 and G 2 , giving the model the ability to reproduce physical results. The enhanced model clearly improved the results by reducing the spurious current and avoiding all bubbles coalescing together. The spurious current acted as a barrier for simulating higher density ratios. The enhanced model handled higher density ratios than the original model with a similar magnitude of spurious current. However, the spurious current was still too large in both models. Since the velocity and concentration fields were coupled, the high spurious current affected the concentration field and, consequently, the morphology of the dendrites and reliability of the results.

Dendritic Growth with the Phase Field Model
The same dendrite growth problem, as defined in Section 7.4, was solved with the phase field LB model. For comparison purposes, the density ratio was selected as 40 (ρ l = 1 and ρ g = 0.025) even though the phase field model can simulate much higher density ratios. The results are shown in Figure 10, for the same time values of Figure 8.
Unlike the Shan-Chen models, the dissolution and reappearance of bubbles related to the EOS were not observed in this model. Therefore, the size of the bubbles remained almost constant during the simulation. A small movement was observed for bubbles beyond the dendrite tips while the bubbles between dendrites were trapped and did not move.

Dendritic Growth with the Phase Field Model
The same dendrite growth problem, as defined in Section 7.4, was solved with the phase field LB model. For comparison purposes, the density ratio was selected as 40 (ρ l = 1 and ρ g = 0.025) even though the phase field model can simulate much higher density ratios. The results are shown in Figure 10, for the same time values of Figure 8.   Table 4 summarizes the magnitude of the spurious current, the simulated density ratio, and the computational time. The original Shan-Chen model could only handle the density ratio of four while it diverged for the higher density radios. The phase field method was capable of simulating much higher density ratios, but since the enhanced Shan-Chen model results diverged for the higher density ratios, the selected density ratio was 40 for both models. It can be observed that the original Shan-Chen model was the fastest model, but it produced an unacceptable large spurious current, even for a low density ratio. The phase field model was the most precise model, though at the expense of a significantly higher computational cost. To compare the magnitude of the spurious current among all models, the velocity profile at t = 0.0013 s was plotted along the centerline of the domain in both the x and y-direction as shown in Figure 11.
The velocity component was formed by the physical velocity and the spurious current. The reason for the existance of any physical velocity was related to the pressure difference throughout the domain, which was modeled in Shan-Chen models by the EOS and through Equation (38) in the phase field method. From the figure, it is observed that the magnitude of the spurious current was not constant in the whole section. The maximum spurious current took place near the bubble/fluid interface in all the models i.e., in the bottom of the domain in Figure 11a,c. As expected, the phase field method generated the least spurious current, in the order of 5 × 10 −5 m/s. To compare the magnitude of the spurious current among all models, the velocity profile at t = 0.0013 s was plotted along the centerline of the domain in both the x and y-direction as shown in Figure 11. Figure 11. Velocity magnitude along the centerline in the width and height direction (a) X*=1/2 Shan-Chen model, (b) Y*=1/2 Shan-Chen model, (c) X*=1/2 phase field, (d) Y*=1/2 phase filed, at t = 0.0013 s.
The velocity component was formed by the physical velocity and the spurious current. The reason for the existance of any physical velocity was related to the pressure difference throughout the domain, which was modeled in Shan-Chen models by the EOS and through Equation (38) in the phase field method. From the figure, it is observed that the magnitude of the spurious current was

Comparison with PFMI Experiments
To have a better understanding of the influence of the spurious current in the simulation, the pore formation and mobility investigation (PFMI) experiment performed at the International Space Station is presented here. In the experiment, 10 mm width quartz tubes of 30 cm length were filled with SCN-0.24wt% H 2 0 alloy, and bubbles of nitrogen were injected into the sample to ensure that porosity would be present during the microgravity experiments. During the solidification, a clockwise circulation of small bubbles was observed in front of the tip of dendrites as seen in Figure 12 [61].
Since the experiment took place in the space station, the only possible explanation for fluid flow is Marangoni convection. The average bubble velocity due to Marangoni convection in the experiment was in the order of 0.2 mm/s. In enhanced Shan-Chen simulation, the magnitude of spurious current for the density ratio of 40 was around 3.5 mm/s. This implies that the Shan-Chen models generate a high spurious current, which makes them useless to simulate Marangoni convection. Therefore, a phase field model should be utilized for this situation.
However, the Shan-Chen model could successfully simulate elongated porosity based on the in situ X-ray tomography result of Ref. [62]. The experiment refers to directional columnar dendrites of Al-30wt%Cu alloy contained between two plates as shown in Figure 13. The dendrites were growing downward at an angle of about 60 • from the horizontal. The contrast of the image was proportional to the atomic number of the elements (i.e., low Cu concentration appear white in the image). The eutectic line was a nearly straight thin line, shown at about the middle of the Figure 13 images. During the solidification, the shape of the left bottom bubble changed from circular to worm-shape as it can be seen in Figure 13 from the left to right images. The elongated shape can be seen for both left and right bubbles. Since the model parameters were not the same as the experiment, it is hard to compare the results directly. However, the Shan-Chen models can predict that the shape of the bubble changes from a circular to an elongated elliptical shape due to adhesive force between the bubbles and dendrite.
Since the experiment took place in the space station, the only possible explanation for fluid flow is Marangoni convection. The average bubble velocity due to Marangoni convection in the experiment was in the order of 0.2 mm/s. In enhanced Shan-Chen simulation, the magnitude of spurious current for the density ratio of 40 was around 3.5 mm/s. This implies that the Shan-Chen models generate a high spurious current, which makes them useless to simulate Marangoni convection. Therefore, a phase field model should be utilized for this situation. However, the Shan-Chen model could successfully simulate elongated porosity based on the in situ X-ray tomography result of Ref. [62]. The experiment refers to directional columnar dendrites of Al-30wt%Cu alloy contained between two plates as shown in Figure 13. The dendrites were growing downward at an angle of about 60° from the horizontal. The contrast of the image was proportional to the atomic number of the elements (i.e., low Cu concentration appear white in the image). The eutectic line was a nearly straight thin line, shown at about the middle of the Figure 13 images. During the solidification, the shape of the left bottom bubble changed from circular to worm-shape as it can be seen in Figure 13 from the left to right images. The elongated shape can be seen for both left and right bubbles. Since the model parameters were not the same as the experiment, it is hard to compare the results directly. However, the Shan-Chen models can predict that the shape of the bubble changes from a circular to an elongated elliptical shape due to adhesive force between the bubbles and dendrite. Mathiesen, Sintef, [62]).
In conclusion, the Shan-Chen models were capable of modeling phenomena observed in experiments, such as the narrow and long worm-like bubbles in the interdendritic regions. However, they were unable to simulate phenomena that produce characteristic velocity in the order of 1 mm/s or lower since the spurious current masks the physical flow.

Conclusions
This paper has compared the original Shan-Chen method, the enhanced Shan-Chen model, and the phase field lattice Boltzmann method for the problem of bubble-dendrite interaction during solidification. The original Shan-Chen model has already been studied for this problem, but the results were limited to small density ratios between the fluid and gas phases, which do not correspond to real physical properties. In this article, to attempt the simulation of higher density ratios and reduce the magnitude of the spurious current, several individual enhancements available for the original Shan-Chen model were, for the first time, used in a combined form to develop an enhanced version of the original model. It was found that the enhancements, including the use of a realistic equation of state, a high order isotropy force scheme, and a mid-range repulsion force, were able to improve the accuracy and capability of the original Shan-Chen method. Even with implementation of all these enhancements, the spurious currents were still large, which makes the model unsuitable to study the Marangoni effect in its current form and the enhanced versions implemented in this manuscript. A more reliable and accurate alternative is the phase field lattice Boltzmann method. Although computationally more expensive than the Shan-Chen model, it is able to produce accurate results not affected by spurious currents. In conclusion, the Shan-Chen models were capable of modeling phenomena observed in experiments, such as the narrow and long worm-like bubbles in the interdendritic regions. However, they were unable to simulate phenomena that produce characteristic velocity in the order of 1 mm/s or lower since the spurious current masks the physical flow.

Conclusions
This paper has compared the original Shan-Chen method, the enhanced Shan-Chen model, and the phase field lattice Boltzmann method for the problem of bubble-dendrite interaction during solidification. The original Shan-Chen model has already been studied for this problem, but the results were limited to small density ratios between the fluid and gas phases, which do not correspond to real physical properties. In this article, to attempt the simulation of higher density ratios and reduce the magnitude of the spurious current, several individual enhancements available for the original Shan-Chen model were, for the first time, used in a combined form to develop an enhanced version of the original model. It was found that the enhancements, including the use of a realistic equation of state, a high order isotropy force scheme, and a mid-range repulsion force, were able to improve the accuracy and capability of the original Shan-Chen method. Even with implementation