Three-dimensional cellular automaton simulation of coupled hydrogen porosity and microstructure during solidification of ternary aluminum alloys

Hydrogen-induced porosity formed during solidification of aluminum-based alloys has been a major issue adversely affecting the performance of solidification products such as castings, welds or additively manufactured components. A three-dimensional cellular automaton model was developed, for the first time, to predict the formation and evolution of hydrogen porosity coupled with grain growth during solidification of a ternary Al-7wt.%Si-0.3wt.%Mg alloy. The simulation results fully describe the concurrent nucleation and evolution of both alloy grains and hydrogen porosity, yielding the morphology of multiple grains as well as the porosity size and distribution. This model, successfully validated by X-ray micro-tomographic measurements and optical microscopy of a wedge die casting, provides a critical tool for minimizing/controlling porosity formation in solidification products.

growth was considered in the model, and porosity evolution with growth of multi-grains over time was simulated as well. Solute and hydrogen distributions were obtained based on diffusion during solidification. Porosity morphology was then visualized from the simulated results, and the porosity size and distribution were computed. The effects of initial hydrogen concentrations and cooling rates on porosity evolution were investigated. Finally, the simulation results were validated by comparing with the experimental results from a wedge die casting.

Results and Discussion
Grain nucleation/growth as well as porosity nucleation/growth were simulated using the CA model established in this paper. The calculation domain is comprised of a 200 × 200 × 20 mesh cube with a uniform cubic mesh size 5 μm. Symmetrical boundary conditions are used on all the six surfaces of the calculation domain. The temperature is assumed to be homogeneous in the entire domain with a constant cooling rate 50 K/s, and the initial temperature is the liquidus temperature. The porosity nucleation parameters used are as follows: the initial hydrogen concentration in the melt is set as 0.5 mL/100 g Al, the maximum pore nucleation density = × N 1 10 . As solidification proceeds, the melt temperature decreases, nucleation and growth of grains occurs, and hydrogen pores start to form and grow within the calculation domain. Figure 2 shows simulated results of both grains and porosity at different time steps, allowing the examination of the grain and porosity morphologies, solute concentrations (Si and Mg), as well as hydrogen concentrations at different time steps.
It can be seen from Fig. 2(a1-a4,b1-b4) that the growth of grains and porosity are related and their structures are intertwined. The formation of porosity hinders the growth of grains, and vice versa. As the gas pores are surrounded by solidifying grains, they become enclosed and are unable to escape from the liquid, resulting in porosity. Figure 2(c1-c4,d1-d4,e1-e4) show that the hydrogen distribution is significantly different from the Si and Mg distributions. As the temperature decreases, the solute concentration in the remaining liquid is always found to increase due to solute rejection from solid to liquid during liquid-solid transformation. Furthermore, the nucleation and growth of pores do not affect the distributions or the diffusion coefficients of Si and Mg in Al. When the temperature reaches the eutectic temperature, the Si concentration in the remaining liquid increases to the eutectic concentration of 12.6 wt.% as shown in Fig. 2(c4). Conversely, the hydrogen concentration first increases and then decreases in the latter stages of solidification. During the liquid-solid transformation the temperature of the melt will decrease and is accompanied by grain growth and an initial increase in the hydrogen concentration. However, when the hydrogen concentration in the liquid is increased to a value larger than the hydrogen saturation limit and the nucleation criterion is fulfilled, the nucleation of hydrogen gas pore occurs. Stable hydrogen pores are then able to grow and act as sinks where hydrogen is absorbed within the pore, resulting in a decreased hydrogen concentration in the remaining liquid. In addition, the hydrogen diffusivity is much larger than solute diffusivity, which results in a hydrogen concentration that is less accumulated in the liquid near S/L interface compared to the Si and Mg concentrations as shown in Fig. 2(c1,d1,e1).
The corresponding solid fraction and porosity percentage with time were recorded quantitatively as shown in Fig. 3. During the solidification process, Fig. 3 shows that the solid fraction of primary Al increases first, and then converges to a constant value. It can be seen that the solid fraction increases after the beginning of solidification, and the increasing rate of the solid fraction increases from zero, and then decreases to zero. This phenomenon can be described by the nucleation and growth of grains and the effects of undercooling and cooling rate. Since the existence of Si particles, Mg 2 Si, and eutectic solidification were ignored in the simulation, the solid fraction of primary Al could not increase when the Si concentration in the remaining liquid reaches eutectic concentration of Si (12.6 wt.% in this case) and all remaining liquid is considered as eutectic. Such an assumption explains in the flat stage of the solid fraction that appears at the end of solidification.
Also as shown in Fig. 3, the percentage of porosity increases from zero to a constant value first, and is followed by an exponential increase. The solid fraction increases earlier than the percentage of porosity, which is expected due to the larger hydrogen solubility limit within the liquid. In the initial stage of the solidification, the porosity percentage remains zero, which means there is an incubation time for porosity nucleation. As solidification continues and the liquid-solid transformation increases the hydrogen concentration in the liquid, the nucleation of the pores will occur as long as the local hydrogen concentration is larger than the local hydrogen saturation and the nucleation criterion is satisfied. It should be noted that after the pore nucleation, the hydrogen solubility changes greatly. The small initial radius of the nucleated pore leads to a large internal pressure, and results in a large hydrogen solubility. On the other hand, the nucleation of the pores will absorb hydrogen around them. Thus, pores are not able to grow immediately after nucleation as the local hydrogen concentration is lower than the local hydrogen solubility. The incubation time which appears in pore nucleation and in pore growth has been previously reported 27 . As solidification proceeds, the local hydrogen concentration in the melt continues increasing. When the local hydrogen concentration near the pore is higher than local hydrogen solubility limit again, the pore will start to grow, resulting in higher percentage of porosity. It can be seen that the porosity percentage rapidly increases, and the initially increased rate of solid fraction slows down as the temperature is lower than 870 K (597 °C). Figure 4 shows the simulated morphologies of porosity and hydrogen concentration fields when the temperature is decreased to the eutectic temperature point with various cooling rates of 50 K/s, 10 K/s and 5 K/s. It can be clearly observed that as the cooling rate is decreased, the total number of individual pores decreases, yet the size of individual pores increases. Similarly, the lower cooling rate is also shown to increase the grain size as shown in Fig. 4(d,e). Conversely, it is known that a higher cooling rate results in increased nucleation of grains, faster grain growth velocity, and thus a finer grain size. The increased number of solidified grains hinders the diffusion of hydrogen in the remaining liquid and the growth of the nucleated pores. Due to the effect of the interfacial tension at Gas/Liquid interface, the nucleated pores should grow spherically. With the growth of grains and www.nature.com/scientificreports www.nature.com/scientificreports/ pores, the growth of pores is blocked, leading to irregular shapes of the pores. The hindered growth leads to an increase in the total number of individual pores as the melt has less time for hydrogen diffusion and pore growth, thus, smaller size and increased number of porosity. Less time for hydrogen diffusion also has a consequence on the hydrogen concentration in the remaining liquid, as more hydrogen will be present in the remaining liquid www.nature.com/scientificreports www.nature.com/scientificreports/ compared to lower cooling rates. The larger the pore radius, which is associated with lower cooling rates, the lower the local hydrogen saturation. Therefore, more hydrogen molecules will contribute to the growth of pores, which results in the lower hydrogen concentration within the melt at the end of solidification. On the other hand, several pores have large radii, while the other pores are very small as shown in Fig. 4(c). It should be noted that the internal pressure P G of a large pore is lower than that of a small pore, which leads to a lower solubility of a large pore than that of a small pore. The lower the internal pressure and solubility, the larger the gas volume increment. Following the nucleation of pores, several pores grow larger than the others, and then they will continue growing and grow faster than the other pores. It can be found that the hydrogen concentration in the remaining liquid as shown in Fig. 4(f) is around 0.45 mL/100 g Al, which is even lower than the initial hydrogen concentration of 0.5 mL/100 g Al. The lower levels of hydrogen can be explained by porosity growth. As the local hydrogen saturation decreases, the pores continue to act as sinks and absorb hydrogen atoms, causing a decrease in the hydrogen concentration of the remaining liquid. Figure 5(a) plots the percentage of porosity vs. temperature with different cooling rates shown in Fig. 4. It can be seen that each of the percentage of porosity curves have an incubation time for the nucleation of pores and for the growth of pores as analyzed above. With the decrease of cooling rate, the final percentage of porosity at the eutectic temperature increases. This is due to that the higher cooling rate limits the time for the nucleated pores to grow. In the higher cooling rate situation where hydrogen atoms contribute less in the growth of pores, they remain in the liquid and lead to the higher hydrogen concentration. On the other hand, with a lower cooling rate, pores begin to nucleate at a lower temperature, and the nucleation period is shorter. This occurs because with a lower cooling rate, the hydrogen accumulated at the S/L interface will have more time to diffuse away at a slightly elevated temperature. It takes more time for the local hydrogen concentration to increase to a value larger than the local hydrogen saturation, and thus pores start to nucleate at a lower temperature. Furthermore, at a lower cooling rate, the hydrogen concentration in almost all of the remaining liquid is higher than the hydrogen saturation before the pore begins to nucleate, which means it will take less time for the hydrogen concentration to increase beyond the hydrogen saturation again after the nucleation of pores. Figure 5(b) shows the porosity percentage vs. temperature with different initial hydrogen concentrations of 0.4, 0.5, and 0.6 mL/100 g Al. It is evident that the final porosity percentage at the eutectic temperature increases with the increase of initial hydrogen concentration. With a lower initial hydrogen concentration, pores start to nucleate at a lower temperature. This is explained by hydrogen diffusion, as with less hydrogen it is able to diffuse away from the S/L interface relatively quickly and maintain a hydrogen concentration below the critical value. When the initial hydrogen concentration is increased, even though hydrogen is diffusing away from the interface, some hydrogen remains and accumulates to a value which satisfies the critical hydrogen concentration point for nucleation. Thus, higher levels hydrogen will nucleate quicker and at higher temperatures when all other conditions held constant. All percentages of porosity resulting from different initial hydrogen concentrations continue to increase as the temperature decreases, even at the end of the solidification. The continued rise in porosity is due to the growth of the pores and the decrease of temperature. This phenomenon leads to the observation that the hydrogen concentration is always higher than the hydrogen saturation, and the hydrogen atoms absorbed by the pores continue contributing to the growth of porosity even following solidification. With the decrease of initial hydrogen concentration, both the total number of pores and the size of pores will decreases. Regarding the hydrogen concentration field, the final hydrogen concentration in the melt is lower when the initial hydrogen concentration is higher. The reason is that the higher initial hydrogen concentration leads to a larger pore radius, which results in a lower local hydrogen saturation, and thus a lower final hydrogen concentration in the remaining melt.
Porosity in casting products can be reduced by controlling the gas quantity in the melt (such as degassing), or optimizing the solidification conditions (such as increasing cooling rates). Since cooling rates vary considerably www.nature.com/scientificreports www.nature.com/scientificreports/ at different locations of an individual casting due to varying section thickness/cooling lines etc., it will undoubtedly result in varying porosity size, morphology, and formation within the casting. Such a porosity distribution will lead to location-specific mechanical properties which are often unaccounted for in casting design and  www.nature.com/scientificreports www.nature.com/scientificreports/ manufacturing. Figure 6 shows the comparison of the porosity morphology between X-ray micro computed tomography (microCT) measurements and CA simulated results. As shown in X-ray microCT results, with the highest cooling rate in Fig. 6(a), the pores are constrained between the dendrites. The growth of porosity is blocked by the dendrite morphology, and vice versa. With the decrease of cooling rate, the porosity appears to have a larger size and a more rounded morphology. Porosity with transparent dendrites at the end of solidification can be clearly observed in Fig. 6(b,e,h). The pores appear to be located between dendrites. As the cooling rate decreased, the porosity size increased as expected, while the total number of individual pores decreases. This result is in good agreement with the analyses of Atwood et al. 3 . Since the pores surrounded by dendrites are not able to escape from the remaining liquid, they form porosity in the final casting. It shows that CA simulation can be used to visualize the porosity morphology by simultaneously simulating both porosity and dendrites, and the simulated results agree well with the experimental results.
In the optical microstructure, both dendrites and porosity can be clearly observed. Several pores were identified and examined as shown in Fig. 7. It also can be seen that the dendrite morphologies in the simulated results match well with the experiment results as shown in the white triangle of Fig. 7(c,d). With decreasing cooling rate, the dendrites coarsen and are associated with larger pores. However, the pores did not appear as round morphologies. This is attributed to the 2-D cross-section can only show a slice of complete 3-D information 30 . www.nature.com/scientificreports www.nature.com/scientificreports/ Nevertheless, it shows that the 2-D simulated results present a good agreement with the optical microstructures, both the dendrite morphology and porosity morphology.
It is clear that the cooling rate affects the final percentage of porosity as analyzed above. However, it should be noted that the porosity size is not uniform at a slow cooling rate. Several pores grow to have large sizes, while others are not able to grow as large as seen in Fig. 4(c). With a larger radius, the local hydrogen saturation is lowered in the melt as discussed above, and results in the pore absorbing more hydrogen and growing to a larger size. Meanwhile, pores with small radius are not able to grow. Smaller pores have higher local hydrogen saturations, and the hydrogen atoms in the liquid melt are absorbed by the larger pores when there is enough time for hydrogen diffusion at a slow cooling rate. Accordingly, the competitive growth among the nucleated pores becomes more obvious with the decrease in cooling rate, resulting in an uneven distribution of pore radius (Fig. 4(c)) and an increased maximum porosity radius (Fig. 7). Similar phenomenon was noticed by Karagadde et al. 23 who built a model to simulate the hydrogen bubble evolution in micro-scale domain, and studied the effect of cooling rate on average pore radius. They reported that, with increasing cooling rate, the average final pore radius decreases. In the present study the pore radius varied greatly in both the experimental and simulated results. Figure 8 clearly shows that increasing cooling rate www.nature.com/scientificreports www.nature.com/scientificreports/ decreases the percentage of porosity present within the casting. Both the simulation and experimental results in this study were compared to the experimental results collected by Gao et al. 26 and Sigworth et al. 31 , with excellent agreement. Gao et al. 26 calculated the percent porosity of Al-4.5wt.%Cu alloy with 0.3 mL/100 g Al and a cooling rate of 1.0 K/s by building a theoretical model, where the green triangle in Fig. 8 shows the porosity percentage is 1%. Sigworth et al. 31 studied the mechanisms of porosity formation during solidification, and showed the values of the porosity percentage during solidification of A356 alloy with 0.3 mL/100 g Al, which used the same material and initial hydrogen content as in this study. It should be pointed out that the porosity percentage from the simulation is slightly lower than the experimental results. This discrepancy can also be found in porosity morphology shown in Figs 6 and 7, which is attributed to micro-shrinkage. Micro-shrinkage occurs at the end of solidification due to the lack of liquid feeding and increases the total percentage of final porosity. It should also be noted that several assumptions are adopted in the present modeling and simulation which is different from the real non-equilibrium solidification. The solidification process is assumed to be equilibrium, and only hydrogen gas porosity is considered in the present model. As shown in Fig. 3, approximately 48% of the melt has not solidified at the end of simulation, which is supposed to solidify as eutectic. Further eutectic solidification is not considered in the model since the formation of gas porosity has mostly completed before that. Although shrinkage porosity is not considered in the present porosity model, the evolution of hydrogen porosity (which is dominant in most solidification products of aluminum alloys) can be predicted, and the morphology and distribution of hydrogen porosity can be quantitatively obtained. Validating the present model offers an excellent opportunity to predict location-specific microstructure and mechanical properties of final castings. By understanding where and how pores will evolve in a cast structure, it will help optimize the design of castings and reduce scrap rates in casting production. Therefore, the validated porosity model will provide an indispensable part from location-specific microstructure to location-specific mechanical properties in ICME design and manufacturing of lightweight Al castings.
In summary, a 3-D CA model was built to simulate the formation and evolution of hydrogen porosity coupled with grain growth during solidification process of Al-7wt.%Si-0.3wt.%Mg ternary alloy. The concurrent nucleation and growth of both grains and pores were examined, and the redistribution and diffusion of solute and hydrogen were obtained using the present model. By applying the 3-D CA model, multi-grain growth with porosity evolution can be simulated. From simulated results, porosity morphology and distribution can be visualized qualitatively, and porosity size and percentage can be computed quantitatively. It is found that, with lower cooling rates, the final porosity percentage and the maximum porosity size increases, while the total number of individual pores decreases. In addition, for higher cooling rates, porosity nucleates at higher temperatures and starts to grow at lower temperatures. With higher initial hydrogen concentrations, the final percentage of porosity increases, and both the number and the size of porosity also increase. In addition, wedge die casting experiments were conducted to generate samples with different cooling rates, and X-ray microCT measurements were performed on the samples to obtain quantitative porosity information. The percentage and the distribution of porosity of simulated results agree well with the experimental results. This model can serve as an indispensable part in connecting location-specific process conditions and location-specific microstructure models for ICME of cast products. www.nature.com/scientificreports www.nature.com/scientificreports/ Methods A 3-D cellular automaton model is developed to investigate the nucleation and growth of both hydrogen gas pores and grains. By applying the model, the nucleation and growth of hydrogen gas pores and grains during solidification of ternary aluminum alloys can be simulated in 3-D. Additionally, the solute and hydrogen redistribution can be calculated and the size, percentage and distribution of porosity can be obtained. Further, the effect of cooling rate and initial hydrogen concentration on the porosity size and morphology can be evaluated. In CA modeling, a discrete number of cells are given states including: (1) gas cell (f G = 1), (2) solid cell (f S = 1), (3) liquid cell (f G + f S = 0), (4) the solid/liquid (S/L) interface cell (0 < f S < 1, and f G = 0), (5) the gas/liquid (G/L) interface cell (0 < f G < 1, and f S = 0), (6) the gas/solid (G/S) interface cell (f G + f S = 1), and (7) the gas/liquid/solid (G/L/S) interface cell (0 < f G + f S < 1), where f G and f S are the gas fraction and solid fraction of the cell, respectively. The state of each cell is then updated or transformed based on a set of neighborhood conditions as the solidification process continues. More detailed explanation of the CA process can be found in the authors' previous publications 32,33 . Assumptions adopted in this paper include: (1) only hydrogen pores are considered; (2) the effect of fluid flow is not considered; (3) shrinkage porosity due to solidification shrinkage and inadequate feeding during solidification is ignored; and (4) the eutectic solidification is not considered. As solidification proceeds, grains begin to nucleate and grow. Solute (Si and Mg) is rejected from the solid, and accumulates at the S/L interface due to the phase equilibrium. Since hydrogen is far less soluble in solid Al than in liquid Al, the excess atomic hydrogen is rejected from the solid into the remaining liquid phase. With the increase of hydrogen concentration in the liquid, molecular hydrogen pores form in the gaseous state when the atomic hydrogen reaches the solubility limit in the liquid. Based on the local hydrogen concentration levels and the diffusion rate, pores are able to grow or shrink. As the temperature decreases, thermodynamically stable grains and the hydrogen pores will grow. Firstly, the nucleation and the growth of solid grains is primarily controlled by solute and temperature distributions. The continuous nucleation distribution proposed by Rappaz 34,35 was used. Due to the fact that the rejected solute during solidification increases the solute concentration 36 , the liquidus temperature decreases since the equilibrium partition coefficient is lower than 1, which leads to a decrease in local undercooling 32,[37][38][39] . The local undercooling was adopted to accurately describe the location specific nucleation phenomena. Local equilibrium at the S/L interface and solute diffusion can be expressed as: where k i is the solute partition coefficient of the solute element i, ⁎ C i S and ⁎ C i L are the solute compositions at the S/L interface in solid and liquid respectively, C i E is the solute concentrations in phase E (superscript E represents S in solid or L in liquid), and D is the diffusion coefficient 33 . Based on the lever rule of both solute elements Si and Mg, the solid fraction and the solid fraction increment can be calculated 39,40 .
Secondly, gaseous pore nucleation only occurs when the gas dissolves in the liquid (C H L ) and exceeds the critical supersaturation 27  is the internal pressure of a gas pore, P 0 is the standard atmospheric pressure, γ is the surface tension of the G/L interface, r p is the pore radius, and C Si L and C Mg L are the solute concentrations of Si and Mg in melt Al, respectively. The supersaturated hydrogen originates from liquid-solid transformation 40,41 . Furthermore, any hydrogen that is located within the liquid state is able to be absorbed by adjacent pores (as shown in Fig. 1(c)). The diffusion of hydrogen can be expressed as: , where C and D are the hydrogen concentrations and the diffusion coefficient respectively, and k H is the partition coefficient of hydrogen 40,41 . Then the gas volume increment of the pore, ∆V, can be related to the hydrogen absorbed at the interface cell, V G , in one time step interval based on the ideal gas law: , where N G means the amount of hydrogen, R is the gas constant, Σ includes all G/L interface cells of the pore, V cell is the volume of a cell, and ρ is the density of liquid Al melt. The volume change of the G/L interface cell (i, j, k) of the pore A at any given increment can be evaluated as: is geometrical factor of the G/L interface cell (i, j, k). The geometrical factor G G for 2-D system proposed by Zhu et al. 27 is extended to 3-D system as follows: . The time step for hydrogen diffusion, ∆t H , is deter- where D H L is the diffusion coefficient of hydrogen in melt Al. Due to the fact that the diffusion coefficient of hydrogen in the melt Al is much larger than that of solute (Si or Mg), the time step used in hydrogen diffusion is much smaller than that of solute diffusion. To increase the computational efficiency, two different time steps are used in the simulations. The diffusion of hydrogen is calculated for = ∆ ∆ t t Nt / C H times using ∆t H . ∆t C is used in the calculation of the diffusion of solute Si, and the evolution of grains and gas pores. In this CA model, the anisotropy 36,42-45 from the regular coordinate system in a 3-D system was employed as detailed in previous research 46 . Based on the definitions above, porosity evolution is coupled with grain growth. Ternary Al-7wt.%Si-0.3wt.%Mg was used as an approximate composition in the above modeling and simulations. The material parameters in the literatures 14 www.nature.com/scientificreports www.nature.com/scientificreports/ To validate the CA model, a V-shaped wedge casting experiment was performed. Primary A356 ingots were melted at 973 K (700 °C) in graphite crucibles coated in boron nitride using an electric resistance furnace, and the surface dross was removed prior to casting. K-type thermocouples with 0.2 mm diameter were placed in three different heights along the wedge casting, corresponding to different thickness, i.e. different cooling rates. Thermocouple wires were shielded in alumina tubes, while the thermocouple junctions were exposed directly to the melt for faster response and improved accuracy. The cooling curves were collected using a computer data acquisition system (NI 9219 from National Instruments) with a sampling rate of 50 readings per second. The instantaneous cooling rates at different locations calculated directly at the liquid temperature were 64.8, 10.7, and 2.5 K/s, respectively.
Prior to casting, the initial hydrogen concentration was measured using the Straube-Pfeiffer vacuum solidification test 48 or equivalently the reduced pressure test (RPT) 49 . In this study, the RPT was carried out on three specimens, followed by density measurements to obtain the average density of the melt. For each test, approximately 100 g of A356 alloy was first melted in a muffle furnace at 973 K (700 °C). Then the sample melt was poured into a steel cup and solidified at a pressure of 100 mmHg for 20 minutes. The average density was calculated based on the Archimedes Principle after the dry weight and wet weight of each sample was measured. The hydrogen content can be obtained as C H 0 = 0.3 mL/100 g Al. After wedge casting experiments, X-ray microCT measurements were performed to characterize the porosity in the casting samples.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.