Shan–Chen model is a numerical scheme to simulate multiphase fluid flows using lattice Boltzmann approach. The original Shan–Chen model suffers from inability to accurately predict behavior of air bubbles interacting in a non-aqueous fluid. In the present study, we extended the Shan–Chen model to take the effect of the attraction–repulsion barriers among bubbles in to account. The proposed model corrects the interaction and coalescence criterion of the original Shan–Chen scheme in order to have a more accurate simulation of bubbles morphology in a metal foam. The model is based on forming a thin film (narrow channel) between merging bubbles during growth. Rupturing of the film occurs when an oscillation in velocity and pressure arises inside the channel followed by merging of the bubbles. Comparing numerical results obtained from proposed model with metallography images for aluminum A356 demonstrated a good consistency in mean bubble size and bubbles distribution.

Demanding for advanced materials are increasing rapidly via new technologies. Closed cell metal foams gained a lot of interest as one of the major branches of advanced materials due to their unique physical and mechanical properties, including high specific strength and compressibility along with good energy absorption capability [1–4]. Despite the advantages, the employment of metal foams in industrial applications is limited due to the inhomogeneity of the structure which results in the deviation of the mechanical properties of the foams from what predicted by the scaling relations. This is mainly due to the morphological defects such as missing or wavy distortions of the cell walls and non-uniform shape and size of the cells which results in poor reproducibility of foam structures [5]. In metals, unlike ionic liquids, the formation mechanism of metal foam has not yet fully understood [5].

Bubble stability is the primitive challenge in understanding the mechanism of metal foam formation. A variety of studies and researches have been performed by scientists in order to investigate and analyze the parameters affecting bubble stabilization [4–14]. Most of the investigations have focused on formation of single bubble in ionic liquid environment, especially water, with no impurity [8,10–14]. However, even in the purest condition, molten metal contains dozens of different impurities. Another important issue which has been neglected during various studies is the multi-bubble nature of metal foam formation process, which mostly appears in bubbles interactions. In addition, due to the presence of metallic bond in metal melt, there is no ionic or polar attraction and repulsion forces, which causes a different behavior of the liquid–gas interface in metal melts in comparison with aqueous solutions [4–7].

In order to have a computational study on metal foam formation process, a basic understanding of bubble stability conditions in the presence of particles is required. Therefore, a computational model based on known dynamics of bubbles and improving it using the computational–experimental approach has to be built, which is accomplished by adding some constraints that are focused on the boundary of each bubble according to available theories and verify the selected ones using an experiment [5].

Most of the studies in the field of bubble dynamics investigated single bubble dynamics in aqueous solutions or water based liquids and only a few are conducted to study multi bubble dynamics and foam formation process.

Chine and Monno [9] developed an axial symmetry model to simulate the behavior of a single gas bubble expansion, embedded in a viscous fluid using finite element method. Ghosh and Das [15] conducted a numerical investigation of single bubble dynamics using lattice Boltzmann method. Their model contains a rising gas bubble inside a tube filled with liquid. Chahine [10] studied the dynamics of clouds of bubbles via both analytical technique and numerical simulation using 3D boundary element method (BEM). He also investigated influence flow field to interaction of bubbles [11]. Ida [12] conducted a mathematical modeling using a nonlinear multi-bubble model from pressure pulses perspective. Bremond et al. [13] performed an outstanding computational investigation on bubble interactions and validated it using a novel experimental procedure based on Rayleigh–Plesset equation. Kim et al. presented an immersed boundary method (IBM) to simulate and predict the 2D structure of a dry foam using a set of thin boundaries to partition the gas into discrete cells or bubbles [16].

The lattice Boltzmann (LB) simulation has been used extensively to simulate the kinetic effects on bubbles [17]. Advantage of LB approach lies in the fact that there are no global systems of equations which have to be solved. Besides, boundaries in simulation domain do not have an effective impact on the computation time. These features are essential for foam formation simulation regarding the complex internal structure of foams. Anderl et al. [18] used lattice Boltzmann to simulate flow in a simplified single phase model as an alternative to a liquid–gas two phase model to analyze bubble interaction in protein foams in order to determine bubble coalescence conditions. Leung et al. [19] studied bubbles nucleation, growth, stability conditions and interaction in plastic foaming process. Beugre et al. [20] developed a 3D lattice Boltzmann code to simulate fluid flow in metal foam. Pressure drop was the criteria used to compare the obtained results with experimental measurements. Körner [5] conducted a thorough research using lattice Boltzmann method to simulate the growth and formation of an aluminum foam using a single phase model. Diop et al. simulated solidification process of metal foams using lattice Boltzmann method [21].

In this study we constructed a hybrid experimental–computational model to predict the structure of a two-dimensional closed-cell aluminum foam. Most of the simulations conducted earlier in this field are based on single phase models, while a two-phase computational model is adopted in the present study. To this end, a modified and improved Shan–Chen model is developed for multiphase simulations. Besides, a novel boundary condition is utilized to simulate particle network effect on the interaction of two bubbles. Moreover, some additional phenomena such as random nucleation, growth, coalescence, and aging of bubbles are also developed in this model.

2Numerical modeling2.1Lattice Boltzmann approachThe LB method has shown to be suitable for foam formation problems [22]. Random micro bubbles in a virtual medium nucleate and interact within a set of rules. If correct physics is applied in the simulation, spontaneous hydrodynamic behavior can be expected. It can be said that LB method is a mesoscopic approach that is between macroscopic CFD approaches and microscopic molecular dynamics. Many multiphase models exist that use the LB method, such as immiscible lattice Boltzmann (ILBM) [23] Shan–Chen model [24,25], free energy model [26,27], chromodynamic model [28] and HSD model [29].

LB method models fluid dynamics by evaluating particle distribution function f (x, v, t) at each lattice point, where f is the probability of finding a moving fluid particle with velocity v at point x and time t. By knowing f, one could get the values of density and momentum. The distribution function used in LB method, fi, is a discretized form of the main continues function.

Shan–Chen model is based on incorporating long-range attractive forces (F) between distribution functions. In the original Shan–Chen model, the interacting force is approximated using the following equation [24,25]:

where b is the number of nearest sites with equal distance, D is the dimension of the space and g is proportional to the interaction strength. The function ψ is so-called pseudopotential and is a function of time and location.By analogy with classical mechanics, the potential of the force can be introduced as:

where coefficient G controls the strength of the attraction. Since the gradient terms in Eq. (2) are in small compared to the leading terms (the characteristic length of the interface is longer than the lattice spacing, as in all diffuse-interface methods), the Eq. (2) can be approximated (ρ is the density, cs is the speed of sound in solids, and p is the pressure for the equation of state for the LBE):By a suitable choice of the pseudo-potential ψ (x), this equation can describe the separation of phases. One simple and usual choice can be ψ=ρ. By using this pseudo-potential function, the momentum flux tensor resembles diffuse interface method. However, the choice of ψ=ρ (x, t) is not the best in terms of stability. When ψ equals ρ, it becomes larger for larger ρ. Thus, the attractive potential contains a malfunctioning loop: the larger density ρ leads to a larger ψ which causes larger gradients and instabilities. ψ=ρ is good for small gas–liquid density ratios.

In the case of aluminum liquid and hydrogen gas (two-phase system), density ratio is considerably higher. Therefore, to handle the pseudopotential function ψ for larger ρ while preserving its ratio for smaller ρ, the following choice of the pseudopotential is used:

which is for small ρ equals to ψ (x)=ρ and for large densities, ψ (x)=1. This choice of the pseudopotential allows separation of gas and liquid in larger density ratios (if not more than 60–70) [30].The critical value, when separation occurs, can be calculated from the thermodynamic theory by these two equations [31]:

By substituting Eq. (4) in Eq. (3) and for D2Q9 lattice (cs=1):

And from Eq. (5) and Eq. (6) one would get:

Solving these equations lead to Gcritical=−4 and ρcritical=ln2. This means that if the system is initialized with the liquid density more than ln2 and the gas density less than ln2 in simulations with G≤−4, the result is stable and separation will occur.

2.2Bubble nucleation and growthNucleation occurs at random points in the domain. After nucleation, bubble growth will begin. For growth modeling in each time step, gas will be added by a virtual blowing agent (proportional to blowing agent weight and gas production rate) in each bubble, and then the bubble volume will be increased due to the pressure balance [32]. Gas blowing or growth will continue until all virtual blowing agent gases are added to the domain.

Number of bubble nucleation sites depends on the initial amount and size of the blowing agents. This number which is based on experimental results, is initially inserted into the main procedure. A time random subroutine is used to determine nucleation site positions in a 2D domain. These locations called domain gas points are in fact virtual blowing agents which in lattice domain are defined as gas nuclei and hydrogen density resulting from decomposition of the blowing agents, is accordingly calculated for these points. Other lattice points are set as liquid and give aluminum density. After this step, all numbers and parameter are changed to dimensionless parameters by open source OpenLB code [33]. Hydrogen gas release rate is calculated in each time step and added to each lattice point by pressure increment.

where dp/dt is pressure increase rate for each lattice point that refers to gas, A is a constant which depends on gas behavior (in case of ideal gas A=RT/V where R is the gas constant, T is the temperature and V is the volume), dn/dt is the gas release rate in mol/s and N is the population of lattice gas points. Thus, because of pressure increase, pressure expansion is defined by multiphase code that has been developed in this work.Present code has been developed base on OpenLB open source code. The Shan–Chen model was incorporated and some modifications were added to the core structure of this algorithm. The modifications were performed to handle the interaction of bubbles in melt liquid.

2.3Bubbles interaction and modification of Shan–Chen modelInteraction of bubbles in pure liquid, without suspended solid, modeling techniques is completely different from liquids containing floating particles (e.g. SiC particles in molten aluminum). In case of two moving bubbles in pure aqueous liquid, if they are closed to each other, common behavior is increment in their surface curvature, that leads to merging and coalescence phenomena from the contact tip (see Fig. 1). However, as shown in Fig. 2 for molten metals, due to existence of a lot of solid particles (impurities and inclusions), interaction between bubbles have a different behavior. Experimental observations during aluminum foam production show a particles network see Fig. 2 between bubbles, which is often called particle network [7]. This network at interface acts as a mechanical barrier and prevents further cell wall thinning. Therefore, main mechanism of foam stabilization between bubbles is due to particle confinement (see Fig. 2).

To cover this phenomenon, some simple conditions are defined. First of all, it is assumed that each bubble interacts with liquid domain only, i.e. any numerical or logical conflicts between the bubbles are neglected. Secondly, each bubble possesses an interaction zone in the liquid phase as a result of its dynamic and velocity vectors. When these domains reach each other, the attracting force between the bubbles begins its performance. As mentioned earlier, a barrier of oxide-networks is formed between these domains that prevents bubbles’ coalescence. This effect could be modeled as an imaginary pressure in thin walls. By a simple condition the oxide-network effect can be simulated in the LBM code. This condition states in order to calculate corresponding Shan–Chen force at each lattice point when the interaction zones of bubbles collide, one would use the nearest bubble to that point and the effects of the rest of near bubbles are ignored, because their effect is practically neutralized by the oxide-network. This statement yields acceptable results in the final simulation.

The computational procedure is conducted by two separate lattices, one for the melt and the other for the gas. This separation requires the utilization of the pseudopotential function, described in Eq. (4) and at the end of each time step, the lattices are coupled.

Next step is the detection of the lattice points having the material between the melt and the gas (according to their velocity and density) and computing a new velocity for these points:

Now the interaction potential could be calculated at each lattice point of the phases:

The final stage is to compute the velocities according to the calculated interaction potential and the external forces. But in the modified model, a correction is applied on the computed values. New lattices are created for each bubble, and contribution of each lattice (i.e. each bubble) to the velocity of the desired point in liquid lattice is computed:where n is the number of bubbles in the domain. And the final velocities could be calculated:When the interaction condition is applied to simulate the particle network, the bubbles would never cross the particle network barrier and no coalescence occurs in the main domain. To solve this problem, another condition should be defined. If this condition is defined based on experimental data, it would dictate if the distance between two bubbles reach a critical number, the particle network wall condition would be removed, which leads the normal procedure of bubble dynamics to proceed. In mathematical and computational declarations, various variables could be processed in each time step in order to determine when to remove the wall condition. In this study, the pressure field and its second derivative in normal direction are the chosen criteria to develop a condition to model removal of thin film and bubble coalescence. This condition is removing the wall when the second derivative of pressure equals zero:

By this removal of interaction condition, bubbles rapidly merge and their dynamic effect on liquid domain could be observed.

3Experimental modelTo study the structure of a porous metal foam, formgrip method is chosen to produce aluminum foam. This method is a combination of powder metallurgy and melting method. At first, precursor of A356 with TiH2 is produced by melting. Then the precursor containing bubble nuclei is placed in a mold and is heated in furnace. Near the liquidus temperature, the precursor suddenly begins to blow. Bubble growing continues until they reach each other and before coarsening stage, the part is solidified. Then it is prepared for metallography images from cross-section of A356 foam after etching process. Finally, the experimental metallography images will be compared with bubbles images obtained from a simulation code that is developed in the present study. Also accuracy of the present code is evaluated as quantitative.

For producing the precursor, around 350g aluminum 356 alloy is melted in the furnace, around 1.5wt.% blowing agent is mixed with aluminum powder with fraction of 0.5. For improving blowing agent wetting property in aluminum alloy melt and better uniform distribution, at 700°C the blowing agent mixture is added to the melt. In the next step, the furnace temperature is set to 600°C. When the temperature reaches 620°C, the mixer rpm is set to 1500 and stirs for 1–2min. In this step, more stirring causes more gas release. The resulting melt is rapidly casted into a metal mold.

The produced solid is called precursor. The precursor is cut for foaming process according to the mold size (a cylinder with a radius of 40mm and 80mm height). Then the mold and the precursor are heated in furnace to blow. In the present study, the effect of blowing temperature on the stability of foam for 675, 725 and 775°C is investigated. For each temperature, 2–3 samples are produced to estimate the optimized foam processing time for producing foam with minimum density and stable cell structure. Schematic of foam producing steps is shown in Fig. 3. Detailed experimental parameters are listed in Table 1.

4ResultsIn present study, modified Shan–Chen model in LB method is used to simulate the behavior of two in line bubbles in a foam like situation to validate the modified code. Simulation conditions are similar to Table 2. In Fig. 4, the result of the pressure and velocity fields is demonstrated by the plot of the values across a vertical line, before and after the instability caused by the perturbation field of the interacting bubbles.

Simulation conditions of small domain of metal foam by LB method and modified Shan–Chen model.

OpenLB | ||
---|---|---|

Domain dimension (lattice parameter) | 750×500 | |

Fluid density (g/cm3) | 2.7 | |

Gas density (g/cm3) | 0.089 | |

Initial condition | Melt | V=0mm/s P=1atm μ675=1.20×10−3Ns/m2 μ725=1.10×10−3Ns/m2 μ775=1.02×10−3Ns/m2 |

Bubble | V=0mm/s P=1atm μ675=1.82×10−5Ns/m2 μ725=1.87×10−5Ns/m2 μ775=1.93×10−5Ns/m2 | |

Boundary condition | Mirror |

Plot of pressure and velocity across a vertical line in the middle of two bubbles (all of the quantities, including length, pressure, and velocity, are dimensionless in these plots). Left: before instability, right: after instability. Red line indicates the pressure and green line is vertical velocity.

The simulation results of bubble growth in a small domain of closed-cell aluminum foam structure with 6 primary bubbles in the mirror boundary conditions are shown in Fig. 5 and the simulation data are listed in Table 2. In Fig. 6 magnified picture of aluminum foam cellular structure is shown to evaluate the accuracy of this study's code in regards to the detection of cell wall thinning stage in bubble coarsening (see top right of Fig. 6). Small domain of Fig. 6, as mentioned before is simulated by mirror boundary condition which dictates that if any melt comes out of a wall it has to come from an opposite one. Thus, by repeating the simulated mirror domain, domain of a few millimeters can be transformed, by a good approximation, to a several centimeters domain. The result of such repetitions is shown in Fig. 7. Blue and red color of phase contours have been changed to gray scale to resemble the color of aluminum foam cross sections. In this case, it could be assumed that the speed of solidification is fast enough to freeze the molten aluminum foam into solid state as the cell structure maintains its molten state. Thus Fig. 7 can be a part of a solid aluminum foam structure. In other words, Fig. 7 can be considered as a simulation metallography picture of porous aluminum foam structure. Fig. 8 shows the experimental metallography pictures of aluminum A356 foam that produced by Formgrip method in 675, 725 and 775°C, for comparison with the simulated metallographic pictures. The results of this comparison are shown in Figs. 9 and 10.

Imaginary metallographic structure of porous aluminum foam that created from mirror repeating of Fig. 6.

In lattice Boltzmann method, the macroscopic properties of the domain of interest could be predicted by solving the LB equations at mesoscopic scale. In this investigation, this numerical method is used and the Shan–Chen scheme for multiphase modeling is modified to present a code to predict the correct behavior of metal foams. The results of conventional codes are factual for every liquid without any solid impurity, but this behavior is not acceptable for metals melts.

To validate the results of the developed code, the pressure and velocity values across a vertical line is shown in Fig. 4. It is observed that for a specific thickness an instability will occur, which is dependent on the size of the bubbles. This instability could be seen in the pressure profile (right image of Fig. 4). Thus, in the beginning of the instability, the second derivative of pressure (∇2P) profile in the lattice point would be zero, which indicates the initialization of bubbles wall rupture. Then it is possible and more comfortable to determine the wall rupture by using the pressure profile and pressure second derivative instead of detecting the critical thickness. By applying this criterion to the domain of foam formation, the merge condition of the gas bubbles would be detected based on the second derivative of pressure field at each lattice points, which leads to an improved modeling of metal foam formation process.

Fig. 7 shows simulation results of a cross-section of A356 foam based on the modified Shan- Chen method at present code. In this image, a small portion of the domain has been repeated by using periodic boundary condition. Fig. 5 shows foaming stage for small domain of metal foam by time. Bubbles with random distribution began to grow after nucleation. Bigger bubbles have a higher growth rate than the others. Each bubble has an affected zone, and when these zones reach each other, bubbles interaction will begin. Simulation of bubbles’ growth will continue until first cell reaches the wall rupture criterion as like Fig. 6. After this time (best foaming time), foaming simulation process would suddenly enter the aging step and bubbles coalescence, which leads to a sever drainage in metal foam structure. These simulation metallography images could be used for predicting of experimental metallography images (as like Fig. 8) and optimum foaming time. These comparisons are represented in Fig. 9 for 675°C and Fig. 10 for 725 and 775°C, respectively. In addition to visual results, quantitative results are extracted from Fig. 9 simulation data in order to validate. These extracted results are processed by using MATLAB image processor and are illustrated in Fig. 11 and Table 3. These data show a minor error between the simulation and the experiment results. Therefore, present code based on modified Shan–Chen model and LB method could simulate and predict foamy structures for aluminum A356 foams during isothermal foaming process.

Bubble size distribution of simulation and experimental samples of Fig. 9.

Comparison of extracted data from simulation and experimental metallographic cellular structures of Fig. 9.

Simulation | Experimental | Error | |
---|---|---|---|

Bubble percentage | 44.3% | 45.5% | 2.64% |

Foam density (g/cm3) | 1.50 | 1.47 | 2.04% |

Mean bubble size (mm) | 3.03 | 3.05 | 0.66% |

In this investigation, the lattice Boltzmann method (LBM) was utilized for understanding of foaming process by simulation of different stages of the Aluminum A356 foam production process at micro and meso scales. Therefore, to predict the structure of metal foam during foaming process, a model is established to simulate the dynamics of bubbles interaction based on development of Shan–Chen model. The presented model can consider the effect of the attraction–repulsion barriers among bubbles into the A356 aluminum foam liquid due to the solid particles network. In order to validate the presented model, results were compared by both some other reference codes and experimental data. Comparison of cellular structure obtained from the experimental route (experimental metallography) and the numerical code (simulation metallography) shows a good consistency. The presented code is also capable of simulating and presenting virtual metallography images for all aluminum alloys foams. Therefore, this software can be used for controlling and predicting density of foams combined with uniform distribution of bubbles at the metal foams.

Conflicts of interestThe authors declare no conflicts of interest.