Multi-phase-field simulation of microstructure evolution in metallic foams

This paper represents a model for microstructure formation in metallic foams based on the multi-phase-field approach. The model allows to naturally account for the effect of additives which prevent two gas bubbles from coalescence. By applying a non-merging criterion to the phase fields and at the same time raising the free energy penalty associated with additives, it is possible to completely prevent coalescence of bubbles in the time window of interest and thus focus on the formation of a closed porous microstructure. On the other hand, using a modification of this criterion along with lower free energy barriers we investigate with this model initiation of coalescence and the evolution of open structures. The method is validated and used to simulate foam structure formation both in two and three dimensions.

The paper is organized as follows. "Simulation method" section presents the new method and its important ingredients. The maturity of the proposed approach to simulate the formation of a complex foam microstructure is presented in "Results and discussion" section. A summary compiles the most important findings of this work.

Simulation method
The present simulation methodology combines ideas from multiphase flows [20][21][22] with the well-established multiphase-field method for microstructure evolution [23][24][25] . While the presence of fluid flow, capillarity, and wetting phenomena such as bubble-bubble coalescence and triple-phase contact angles in the formation of foams motivate the need for concepts in the field of multiphase flows, involving the multi-phase-field (MPF) method may, at least at first sight, appear less obvious. Therefore, we first introduce the MPF method and highlight the advantage of using such an approach. Dynamical equations, which dictate the time evolution of the system, are given after this subsection.
A multi-phase-field model. As mentioned above, the present approach builds upon the well-established multi-phase-field method. We use a version of this approach, which has first been proposed by Steinbach and coworkers 23,24 . For this purpose, we consider each individual bubble as a separate 'phase' which occupies a certain region of space. The associated phase field function, φ α (x) , then takes the value of 1 if the point x is completely occupied by the bubble with index α ( α = 1, 2, . . . , N − 1 ; the index N being reserved for the surrounding liquid). Similarly, φ α (x) = 0 in the opposite case of complete absence of bubble α at point x . Of more interest is, of course, the situation, where the point x lies at the interface of bubble α with other ones, in which case 0 < φ α (x) < 1 (Fig. 1). We adopt here the interpretation that φ α (x) is the fraction of volume element dV at x which is filled by the phase α . From this, it immediately follows that N α=1 φ α (x) = 1 , where the index α = N accounts for the presence of ambient liquid.
One of the main purposes of setting up a free energy functional in the present study is that it allows to derive, in a systematic way, an expression for the pressure tensor 27 . As will be seen below, the divergence of pressure tensor, ∇ · P , plays a central role in updating the fluid velocity field. Exploring the translational invariance of L and accounting for the constraint that the sum of all phase fields is conserved at any point in space, one obtains (for a derivation see Supplementary Appendix-A) Since all forces and thus ∇ · P vanish at equilibrium, the right hand side of Eq. (2) can be understood as the sum of forces, which arise due to deviations from equilibrium state. By considering changes of the functional integral, Eq. (1), under small deviations from thermodynamic equilibrium, one obtains, Derivatives on the right hand side of Eq. (3) are to be understood as partial derivatives in ordinary sense. Under equilibrium conditions, the left hand side of Eq. (3) vanishes and one obtains the well-known Euler-Lagrange equation.
As the next necessary step, the free energy density, L , must be specified. A standard choice here is a squaregradient model as in the Ginzburg-Landau theory of phase transitions 28 . When extended to the case of multiple phases, it can be written as 19,24,29 where σ αβ is the interface energy between the phases α and β and η denotes the width of the interface. |φ α φ β | is known as the double obstacle (DO) potential 30 and the square bracket gives the contribution to pressure arising from the pair of phases α and β , with h being an interpolation function. It can be shown that h does not affect the equilibrium state of a planar interface 19 . For a sphere in static equilibrium with its surrounding medium, on the other hand, the function h enters the force balance condition (see Supplementary Appendix-A.2).
It is important to note that the pressures p α and p β in Eq. (4) must be evaluated via the equation of states (EOS) corresponding to the phases α and β , respectively. In the present study, the same ideal gas EOS is used for all bubbles ( α = 1, 2, . . . , N − 1 ). For the surrounding liquid ( α = N ), we use the well-known van der Waals equation of state, where c s,G is the sound speed in the gas phase and a, b, and c are the thermodynamic constants, characterizing the van der Waals liquid. It is noteworthy that a and c have the dimension of velocity square but b is a unreachable threshold density, at which fluid pressure diverges.
Density is determined via, with m α (t) being the total mass of the α-th phase at time t, given by m α (t) = t + m α (0) . Here, m α (0) is the initial mass within the bubble and is a constant mass release rate (see also "Results and discussion section"). The use of density from Eq. (6) in evaluating pressure is justified because we consider processes which are slow compared to the speed of sound. Spatial variations of density are thus assumed to homogenize instantaneously both within the bubbles and in the surrounding liquid. The fact that we do not need to survey the process of sound propagation allows the choice of a coarser grid and a larger discretization time step and thus provides a major enhancement of computational efficiency.
Dynamical equations. In the present model, microstructure evolution is governed by the dynamics of phase fields, φ α , with α ∈ {1, 2, . . . , N} . The rule to update the φ-fields accounts on the one hand for the dynamics of fluid and on the other hand for thermodynamic driving forces. Following 31,32 , one writes, where u is the fluid velocity, δF/δφ α stands for functional or variational derivative as given by Eq. (3) and M αβ is the interface mobility. Importantly, Ñ is the number of phase fields present at the point x and must be distinguished from the total number of phase fields N. The fluid velocity field, u(x) , is tracked everywhere in space by solving Navier-Stokes (NS) equations, where µ is dynamic viscosity, ρ is the fluid density and f ext is the external force per unit volume. In regions of space where only a single phase, say α , is present, µ = µ α and ρ = ρ α . In general, however, these quantities are defined via phase-averages, ρ(x) = N α=1 φ α (x)ρ α and µ(x) = N α=1 φ α (x)µ α . Even though it may not be apparent at first glance, Eq. (8) also contains interface effects on the right hand side. This information is encoded in the divergence of pressure tensor, ∇ · P , which accounts both for the variation of hydrostatic pressure with density via the EOS and for interface forces by considering curvature terms 33 , see Eq. (A.14).
Scientific Reports | (2020) 10:19987 | https://doi.org/10.1038/s41598-020-76766-z www.nature.com/scientificreports/ Non-coalescing case: the contact angle. The property that, within multi-phase-field method, one can assign a different phase index to each individual bubble provides a natural way to completely prevent the coalescence process. By doing so, each bubble can be treated as an independent entity with its physical properties. Using this feature, it is easy to consider the effect of additives at least on a qualitative level. For this purpose, we remark that when two bubbles α and β come closer than a threshold distance, additive molecules which surround them repel each other. Here, for the sake of simplicity, we mimic this effect by a free energy penalty per unit area, σ αβ . By tuning this parameter, one can adjust the strength of additive-induced repulsion, thereby also influencing the so-called contact-angle. This is illustrated in Fig. 2 for the case of two bubbles in contact. At static equilibrium, it follows from the force balance that contact angles θ α and θ β obey (Fig. 2a), The idea of an interface energy penalty between two domains occupied by the same phases but having a sort of "mismatch" at the contact zone has been used in 34 for the case of solid state precipitation in Ni-base superalloys. Here, we generalize and apply this idea for the first time to liquid-gas systems.
The above equations (9) and (10) for contact angle deserve some comments. First, it can be easily verified that Eq. (9) reduces to the well-known Young-equation 35 for the three-phase contact angle on a flat solid. To see this, it is sufficient to identify the phase fields α and β with vapor and solid phases, respectively, and use the fact that θ β = 0 for a planar substrate. It is noteworthy that, in this case, Eq. (10) is not valid since the normal projection of σ αN is not compensated by that of σ βN but by elastic forces of the solid body, which resist deformation along the perpendicular direction. Second, Eq. (9) makes sense only if σ αβ ≤ σ αN + σ βN (recall that cosine function cannot exceed unity). This condition is always satisfied for σ αβ = 0 . In this case, and recalling that all interface energies are non-negative, Eq. (9) is solved by θ α = θ β = π/2 (Fig. 2b). However, if the free energy arising from additives at the gas-gas interface exceeds the sum of two gas-liquid interfaces energies, a melt layer forms between the bubbles, thereby reducing the total interface energy (Fig. 2c). This means that, from the perspective of bubble-bubble contact, a complete dewetting process takes place. The corresponding non-wetting condition can be expressed as with q being a non-negative empirical parameter. To see the physical meaning of q, let us consider that by some fluctuation an amount of liquid penetrates into the contact zone between two adjacent bubbles and separates them, thus generating two new liquid-gas interfaces (Fig. 2a-c). If q = 0 , the total free energy will not change by this process, so that a fluctuation in the reverse direction can occur with equal probability and can restores the previous situation. A positive value of q changes this balance in favor of dewetting. As will be shown below, it ensures the existence of a dewetting-force or, equivalently, a free energy barrier against bubble-bubble coalescence.
(9) σ αN cos(θ α ) + σ βN cos(θ β ) = σ αβ , Figure 2. A schematic view of two bubbles in contact, embedded in a melt pool. The presence of additives at the interface between two adjacent gas-bubbles leads to a free energy penalty per unit area, which is represented by σ αβ . (a) At equilibrium, force balance at three phase contact line (the point C in this 2D projection) must be satisfied. When projected onto the σ αβ -line, one obtains for the relation between specific surface free energies, σ αβ = σ αN cos(θ α ) + σ βN cos(θ β ) . Balance of forces along the direction normal to the αβ-interface leads to σ αN sin(θ α ) = σ βN sin(θ β ) . In panel (b), the energy penalty associated with additives vanishes ( σ αβ = 0 ). It then follows from force balance that θ α = π/2 and θ β = π/2 . This means that the melt-gas and gas-gas interfaces must be perpendicular to each other. The horizontal dashed lines serves to highlight the fact that, within the multi-phase-field framework, bubbles α and β remain distinguishable even though their interface energy vanishes. Panel (c) shows what happens if σ αβ > σ αN + σ βN . In this case, bubbles go apart since the two newly formed liquid-gas interfaces have both together a lower free energy than that associated with the agglomeration of additives between two interfaces.
A model which completely avoids coalescence is an idealization which allows to focus on a time window and parameter range, where bubbles rearrange and change their shapes in order to accommodate with boundary conditions, while at the same time keeping their individual character as a physical entity. This way, formation of closed structures containing a disconnected set of bubbles can be investigated. Such structures do usually have good mechanical properties, while at the same time being lighter than the corresponding bulk metallic solid.
Condition for coalescence. If one is interested in a study of open structures, a model that allows for coalescence is needed. An important application of open porous structures is heterogeneous catalysis, where it is desirable that the reactants be able to enter and exit the porous body in order to come in contact with the entire catalytic surface. Other applications use open structures to reduce weight while at the same time optimizing mechanical properties.
Thanks to the flexibility of the multi-phase-field method, it is easy to include the coalescence phenomenon into the model, while at the same time having a safe control over its rate. This latter aspect is important since, as already mentioned above, the coalescence process usually proceeds quite fast, making it difficult to influence the foam structure by tuning the process parameters.
Following a similar approach, which was used in a study of superalloys 34 , we introduce a simple criterion, for the initiation of coalescence. The basic idea is that two bubbles will coalesce if (i) they come sufficiently close so that their distance falls below a certain threshold and (ii) if the force pushing their respective liquid-gas interfaces towards each other, p film , is sufficiently high to overcome the free energy barriers (often referred to as disjoining pressure 12 ), which tends to keep the bubbles apart.
A natural choice for the threshold distance for coalescence is the interface width η . To estimate the disjoining pressure, which tries to push apart adjacent bubbles, we consider the slicing of an additive-rich gas-gas interface (which occurs at the contact zone between two adjacent bubbles) into two liquid-gas interfaces. This process can be regarded as completed when the separation distance reaches the interface thickness, η . The change in free energy during this process is dF = (σ αN + σ βN − σ αβ )A = −q(σ αN + σ βN )A , where we used Eq. (11). A is the interface area considered in this problem. The volume created during this process is dV = Aη . The (disjoining) pressure, which is responsible for this process can now be obtained from the standard thermodynamic relation, � disj = −dF/dV . This gives, It is important to note that the empirical parameter q provides the possibility to freely adjust the disjoining pressure. Of course, in order for disj to be effective, it must be of the same order of magnitude as the free energy densities involved in the problem. Thus, reasonable values of q will be of the order of unity.
In order to obtain a closed expression for coalescence criterion, we estimate the driving force, p film , which pushes two neighboring gas-liquid interfaces, say α and β , towards each other. For this purpose, we first evaluate the net force per unit area on each of these interfaces. For the α th bubble, this force is given by �p αN = p α − (p N + κ α σ LG ) , where p α is the inner bubble pressure, p N is that in the surrounding liquid and κ α is the curvature of the gas-liquid interface on the side close to the α th bubble. Note that, as expected, this force vanishes at static equilibrium, where the Young-Laplace law holds. This is in line with the idea that thermodynamic equilibrium is the state of matter where all forces are in balance (resulting in zero net force) and all macroscopic variables are time-independent. Similarly, �p βN = p β − (p N + κ β σ LG ) is the net force per unit area on the β th gas-liquid interface on the side facing the bubble α . A careful analysis now reveals that the net force, which pushes the two bubbles towards each other is given by Putting all this together, coalescence takes place within the present model if where d is the distance between adjacent liquid-gas interfaces. If the condition (15) is satisfied, the region of space occupied by the two bubbles is identified as a single bubble. Technically, this is achieved by assigning the region of space occupied by φ β to φ α (obviously, the reverse assignment α → β is equally valid). The redundant phasefield (in this example φ β ) is then removed from the list of phase fields. Moreover, since the interface associated with the contact area of the phases α and β disappears, the corresponding specific free energy, σ αβ , has no effect anymore and is thus removed from the list of parameters. Finally, depending on details of implementation, a re-indexing or other schemes can be applied to optimize memory usage.

Scientific Reports
| (2020) 10:19987 | https://doi.org/10.1038/s41598-020-76766-z www.nature.com/scientificreports/ As seen from the above description, in this work, we account for additive effects on foam stability at a qualitative level only. It would be interesting to consider in future studying additives explicitly. A force approach here would be to include additives as small solid particles. However, due to the large size difference between bubbles and particles, this would be computationally very demanding. A more computationally applicable choice would be to consider the time evolution of the concentration field associated with additive particles. This would provide access to local variations of the interface energy and would thus open the way to study a broader range of situations including Marangoni-type effects. This completes the model section. Next, we will show that the model is indeed capable of both completely hindering or conditionally allowing the coalescence between neighboring bubbles (Fig. 4). Then, we will apply the thus developed computational tool to study evolution of microstructure in two and three dimensions.

Results and discussion
As a necessary step, we perform a number of Benchmark simulations to validate the approach proposed above. Aiming to mimic the case of an open melt pool, the density (and consequently pressure) of the melt pool is kept constant. This corresponds to the situation that liquid can leave the simulation domain upon bubble growth. Interestingly, since the sum of all phase fields at any point in space is constrained to unity (see the first paragraph in "A multi-phase-field model" section), the spatial domain occupied by φ N (liquid phase) reduces automatically as the total gas volume increases. The insertion of mass into the bubbles is stopped when the volume occupied by the gas phase exceeds a certain threshold.
The current model is designed to study generic trends which arise from a variation of different parameters such as gas release rate, liquid-gas density ratio, gas-liquid interface energy, and viscosity. Therefore, reduced parameters are used throughout this study. However, for liquid-gas density ratio, we chose a relatively large value of 10,000, which is close to density ratios in typical metallic foams. This already solves one of the challenging issues with regard to physical parameters, namely how to efficiently implement a high liquid-gas density ratio.
Unless otherwise stated, the grid spacing is x = 0.01 and the time step is t = 10 −5 . The unit of mass is set to one. All quantities given below are expressed in these reduced units. This means that to obtain the numerical value of a length, it must be multiplied with x , and velocities have a non-written factor of �x/�t . With this convention, the size of simulation domain is L x = L y = 800 in 2D and L x = L y = L z = 300 in 3D. The numerical interface thickness is η = 6 both in 2D and 3D simulations. The speed of sound is set to c s = 0.12 and constants of the van der Waals fluid, Eq. (5), are set to a = 6.4 , b = 6000 , and c = 0.75 × 10 −6 . For simplicity, we set viscosity of gas and liquid phases identical, µ = µ Gas = µ Liq = 1 . The initial radius of all bubble-nuclei is set to R(t = 0) = 6 both in 2D and 3D simulations. A 'newborn' gas bubble thus has a radius equal to the interface thickness. Gas densities within these bubbles are randomly assigned from the intervals ρ Gas (t = 0) ∈ [2.6 − 2.9] in 2D and ρ Gas (t = 0) ∈ [3.24 − 3.6] in 3D. As mentioned in "Simulation method" section, the amount of gas inside each bubble is then increased linearly with a constant rate of = 0.32 × 10 −4 in 2D and = 0.16 × 10 −4 in 3D. This process is stopped after bubbles grow sufficiently to come into close contact and occupy the entire simulation cell. The density of liquid phase is kept constant through the simulation, ρ Liq = 5200.
As a test of physical consistency, it is instructive to use the above parameters and estimate the speed of sound in the liquid phase. Using the van der Waals EOS in Eq. (5), we obtain www.nature.com/scientificreports/ Inserting the above values for a, b, c and ρ Liq , one obtains c 2 s,Liq ≈ 0.06 and thus c s,Liq ≈ 0.245 , which is roughly twice the speed of sound in the gas phase. It is also noteworthy that the parameter c has essentially no effect on the value obtained for c s,Liq . This is a consequence of the fact that the last term in Eq. (5) becomes important only at moderate densities relevant for condensation processes and has hardly an effect for phase behavior of a dense melt pool considered in the present study.
The liquid-gas surface free energy is set to σ LG = 20.5 in 2D and 10.25 in 3D. For bubble-bubble interface energy, we use σ SS = 2σ LG (1 + q) , which is a special case of Eq. (11). It must be emphasized here that the q-parameter influences the nucleation of coalescence in our model via two closely related mechanisms. On the one hand, a larger q leads to a higher coalescence barrier via disjoining pressure, Eq. (13). On the other hand, it increases σ SS and thus makes the approaching motion of two bubbles towards each other energetically unfavorable.  (15). In panel (a) bubbles remain separated throughout the entire process of growth and deformation because, in addition to choosing a high σ SS , the phase fields representing the two bubbles are not allowed to merge into a single one. In (b) and (c), however, this possibility is introduced via the coalescence condition, inequality (15). The only difference between (b) and (c) is the use of different initial bubble-densities. Otherwise, identical simulation parameters are employed. In (b) coalescence condition is satisfied at a time of t ≈ 225 t and the bubbles merge into a single one (see also panel d). In (c), the initial densities within the bubbles are lower than in (b) so that driving forces always remain below the coalescence-threshold. (d) Variation of the p film with time for the two cases (b) and (c). The horizontal lines marks the threshold pressure, disj which must be overcome for initiation of coalescence. In all the cases shown, the initial bubble radius is R(t = 0) = 30 = 5η. The color code represents the density field, dark blue for the liquid phase and light blue for the gas phase. From top to bottom, each row corresponds to t = 0 , t = 0.9 × 10 4 t , t = 1.5 × 10 4 t , and t = 10 5 t , respectively. The system size is 800 × 800 in the units of grid spacing, the interface width is set to η = 6 �x , and ρ Liq /ρ Gas ≈ 10, 000.
Scientific Reports | (2020) 10:19987 | https://doi.org/10.1038/s41598-020-76766-z www.nature.com/scientificreports/ Benchmark tests. The first benchmark test deals with the static equilibrium between two non-coalescing bubbles. As discussed above, in such a situation, the contact angle satisfies Eq. (12). A validation of this relation is provided in Fig. 3a, where simulation results for contact angle are plotted versus σ SS /2σ LG . For each ratio of σ SS /2σ LG , two sets of simulations are performed using different bubble sizes. It is visible from this plot that larger bubbles reproduce the analytic result more closely. This is related to the fact that determination of angle between curved lines is more accurate if the radius of curvature is larger. In the opposite limit, one would have a systematic bias towards larger angles as curves go more quickly apart with distance from the crossing point. In agreement with this interpretation, it is seen from Fig. 3a that angles obtained for smaller bubbles lie systematically above those for larger R. In Fig. 3b-d, the equilibrium configurations of three different contact angles ( θ/2 = 30 , 45, and 60 • ) for larger radius R = 45 are shown. When two bubbles come into close contact, depending on forces which act upon them, they may undergo coalescence or remain separate. In the present model, the possibility for this bifurcation in dynamic behavior, which plays a key role in structure evolution in multi-bubble systems, is controlled via inequality (15). This issue is addressed in Fig. 4, where two bubbles grow with time until they meet. In the first case shown (panel (a)), the non-wetting condition, Eq. (11) is used with q = 1 . Apparently, the system minimizes its free energy by keeping the bubbles apart. Starting with the same configuration as in panel (a), simulations are repeated using the coalescence condition, inequality (15). Two very similar cases that differ only in the initial gas density are considered. In one of these cases, the two bubbles merge and form a bigger one (Fig. 4b). In the other case, which started with a lower gas density, bubbles remain separated during the entire simulation time (Fig. 4c). This different behavior can be rationalized by a survey of the force, p film , which drives the coalescence process. As seen in Fig. 4d, in the coalescing case, there is a time where p film exceeds the disjoining pressure. In the setup corresponding to Fig. 4c, however, the initial gas density inside bubbles is low and the driving force for coalescence remains below the threshold during the whole simulation.
In the instant of coalescence in Fig. 4b, one of the two phase-fields is deleted and the entire gas domain is assigned to the other one. This is shown in Fig. 5a, where the green lines represent the phase-field profiles of the bubbles at time t = 0 and the black line corresponds to the final state of the single bubble, which is the product of coalescence. All the profiles are plotted along the horizontal line passing through the center of coalescing bubbles. Density profiles before and after coalescence are also shown in Fig. 5b.
The results discussed above clearly show that the present model has the capability to adequately account for static equilibrium and dynamic behavior of bubbles. Most importantly, it provides a physically motivated model to study structure formation in foams by completely preventing bubble-coalescence. At the same time, the model also contains a barrier-controlled criterion for coalescence, which can be used to tune the rate of coalescence and the resulting microstructure. Next, we show the results of many-bubble simulations both in two and three dimensions which underline the maturity of our model in dealing with complex structures.
Many-bubble simulations. We employ the model to simulate formation of foam microstructures for different disjoining pressures. This is achieved by changing the value of q in Eqs. (11) and (13). Since disjoining pressure determines the coalescence barrier, Eq. (15), the number of bubbles which merge will vary with q. This, in turn, will affect the number of pores and their size distribution, leading to variable foam structures. Here, we performed a set of simulations using three different values of q, q 1 = 0.25 , q 2 = 0.24 , and q 3 = 0.23 . Except for the value of q, all the other parameters and conditions were identical in the three simulations reported below. In all the cases shown, bubble nuclei grow due to the increase of mass until they come into contact, deform and rearrange, Fig. 6a-c. The process continues until the system is filled with bubbles separated by the liquid films (see the last row in Fig. 6a-c). For the case of q = 0.25 , coalescence is hindered most effectively and the final structure contains a homogeneous distribution of pores, Fig. 6a. For the smaller value of q = 0.24 (Fig. 6b), coalescence is less suppressed and, consequently, the number of pores is decreased. In this case, distribution of www.nature.com/scientificreports/ bubble size is broadened since larger pores are produced due to coalescence. This results in a less homogeneous foam structure, as compared to the case of larger disjoining pressure (cf. panels a and b in Fig. 6). For even lower disj ( q = 0.23 ), still more bubbles merge into each other and this creates a foam including pores with a very large range of size, Fig. 6c. These results suggest that homogeneity of pore structure is deteriorated by lowering the coalescence barrier. Moreover, as highlighted in Fig. 7, this loss of homogeneity is accompanied by a broader bubble size distribution.
It is noteworthy that foam structures obtained in our simulations (last row in Fig. 6) are qualitatively similar to the experimental results for the aluminum foam (see Fig. 7 of 18 ). Moreover, simulations reproduce characteristic Figure 8. Result of 3D simulation for microstructure evolution of a foam representing bubbles in (a) and the liquid film around them in (b). Each raw, beginning from the top, corresponds to t = 0 , 1.5 × 10 4 , and 5 × 10 4 t . The system size is 300 × 300 × 300 in the units of numerical resolution. The interface width is set to η = 6 �x . The initial density of each bubble is assigned from the range of ρ α (t = 0) = 3.24 − 3.6 and the density of the liquid is constant throughout the simulation, ρ l . Thus, the density ratio is around ρ l /ρ g ≈ 10, 000.
Scientific Reports | (2020) 10:19987 | https://doi.org/10.1038/s41598-020-76766-z www.nature.com/scientificreports/ foam structures such as a homogeneous distribution of polygonal bubbles (Fig. 6a) and a non-homogeneous round bubbles with different sizes and curvy liquid films (Fig. 6c). As the last example, we provide results on a 3D foam-microstructure simulated with the present model. To obtain a homogeneous structure and based on the knowledge gained from the above 2D simulations, we choose a very high free energy barrier for coalescence by setting q = 1 in Eqs. (11) and (13). The result of this simulation is shown in Fig. 8. The left column corresponds to bubbles/pores which evolve as a result of mass addition until they fill the simulation box. During the growth process, as the distance between two adjacent bubbles falls below η , the energy barrier prevents their coalescence. Thus, bubbles deform by growing further until the structure reaches a semi-stable configuration (last row in Fig. 8).

Conclusion
In this paper, we represent a 3D multi-phase-field model for simulation of microstructure evolution in metallic foams. First, a formulation of the model is presented which allows to completely avoid coalescence between bubbles. The model is then extended by a rule to allow for coalescence if (i) the liquid film between the bubble becomes thinner than a threshold and (ii) the force, which pushes the bubbles towards each other, surpasses the disjoining pressure. This approach is validated through benchmark tests and is shown to reproduce the analytically expected results. The thus established model is used to study structure evolution in two and threedimensional cases. The present model is designed such that it allows to study the effects arising from a variation of the effective materials parameters (e.g., liquid-gas interface energy, equation of state in the gas phase) and process conditions (e.g., gas increase rate, the total released gas) on the overall structure of the foam as it evolves with time. It also allows for a study of the effects arising from different nucleation mechanisms on the evolution of foam microstructure. The model can be also easily extended to investigate dynamic evolution of the additive field and its effect of the local interface properties. This way, Marangoni-type effects will become accessible to the model.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.