A thermal fluid dynamics framework applied to multi-component substrates experiencing fusion and vaporisation state transitions

The fluid dynamics of multi-component alloy systems subjected to high energy density sources of heat largely determines the local composition, microstructure, and material properties. In this work a multi-component thermal fluid dynamics framework is presented for the prediction of alloy system development due to melting, vaporisation, condensation and solidification phenomena. A volume dilation term is introduced into the continuity equation to account for the density jump between liquid and vapour species, conserving mass through vaporisation and condensation state changes. Mass diffusion, surface tension, the temperature dependence of surface tension, buoyancy terms and latent heat effects are incorporated. The framework is applied to describe binary vapour collapse into a heterogeneous binary liquid, and a high energy density power beam joining application; where a rigorous mathematical description of preferential element evaporation is presented. Modelling the dynamics of multi-component alloy systems under extreme conditions and predicting the resulting compositions is a complex task. Here, the authors develop a theoretical framework to predict vaporisation and condensation in multi-component flows during electron beam welding and demonstrate how these processes lead to chemical heterogenisation.

H igh energy density manufacturing processes, such as power beam welding and additive manufacturing, utilise concentrated sources of heat, often a laser or electron beam, to initiate localised melting 1 . With sufficient power density, vaporisation of the metallic substrate occurs, permitting the heat source to penetrate further into the substrate material and generate a highly dynamic thermo-capillary structure of molten and vaporised material [2][3][4] . The thermo-capillary, often referred to as a 'keyhole', is maintained by a combination of recoil pressure, as the material vaporises normal to the liquid-vapour interface, and surface tension forces, against body forces such as gravity [4][5][6] . The complex heat and mass transfer dynamics in the region of the thermo-capillary govern its stability, but also determine the heat and mass transfer through the various state transitions and, ultimately, the final configuration of the material in its solid state [7][8][9][10][11][12] . The complexity of the problem is amplified when multiple component species are present, such as in alloy systems or in the joining of chemically dissimilar materials, where mixing between the materials in the fluid state governs the final composition. With large enough power density of an applied heat source, multi-component alloy systems may experience preferential elemental evaporation due to the different vaporisation temperatures of the constituent elements 13 . Once the heat source is extinguished, the collapse of the thermo-capillary and condensation of the various vaporised metallic components is a potential source of significant porosity development in the joint. Understanding the highly coupled and highly transient thermal, compositional, energy, and momentum fields generated in metallic substrates during these high energy density processes is key to optimising such processes, and material systems, for nextgeneration applications, such as nuclear fusion power generation components.
The high temperatures, large magnitude momentum fields, and generally highly transient nature of the thermo-capillary formation, maintenance and collapse make it difficult to observe the complex flow behaviour within the keyhole experimentally 1,11,14 . Therefore, mathematical modelling of these high energy density processes offers unique opportunities to understand the thermocapillary behaviour and the key criteria for keyhole stability. Many mathematical frameworks exist for simulating the thermal fields induced in components due to the application of heat. Traditional approaches neglect the fluid dynamics of the molten and vaporised regions entirely, instead using phenomenological heat source descriptions to represent the complex physics involved in the fluid region and obtain thermal predictions by solving a simple heat conduction problem either numerically or analytically [15][16][17][18] . However, as the fluid dynamics governs the heat and mass transfer, and therefore the thermodynamic evolution of the system, these phenomenological approaches are limited in their ability to be predictive as they do not appropriately capture the physics of the processes. The most advanced descriptions of high energy density state change processes, which aim to capture the complex heat and mass transfer, are thermal fluid dynamics approaches 2,4,8,9,19,20 . Various mathematical approaches have been utilised in the literature to study state changes between condensed and dispersed/vaporised configurations of matter in fluid dynamics frameworks, including: Lagrangian tracking of individual droplets and bubbles [21][22][23][24] , two-fluid Eulerian averaging of dispersed features [25][26][27][28] , and interface capturing and tracking techniques [29][30][31] . Current state of the art approaches used to simulate the thermal fluid dynamics of high energy density processes remain limited. In particular, the vaporisation and subsequent condensation of the metallic substrate are neglected, meaning that the complex fluid-vapour dynamics that maintain the stability of the thermo-capillary (the subsequent collapse of which largely determines the level of porosity accommodated during processing) are not captured appropriately 32,33 . Current thermal fluid dynamics frameworks also only consider the substrate as a single component. Therefore chemical segregation, and compositional heterogenisation are neglected and dissimilar metal systems cannot be considered. Furthermore, preferential elemental evaporation in multi-component alloy systems could never be simulated and understood with current frameworks 6,7,33 . No mathematical framework currently describes a multicomponent metallic substrate experiencing melting/solidification and vaporisation/condensation state changes, with appropriate volume dilation and constriction due to the vaporisation and condensation, respectively. In this work, we present a mathematical framework that fully describes such behaviour, in addition to component diffusion in the liquid and vapour states. The present effort focuses on the volume-of-fluid interface capturing approach, solving a number of fully coupled component advection-diffusion equations corresponding to the number of phases present in the domain, as well as momentum and energy equations for the whole domain.
This thermal multi-component fluid approach is presented and applied to two physically significant scenarios pertinent to advanced manufacturing processes. The model considers the effect of surface tension, temperature-dependent surface tension (often referred to as Marangoni forces), buoyancy, and recoil velocity on the formation of the thermo-capillary. A Gaussian representation of the power beam heat source is utilised in the energy equation. Volume changes due to vaporisation and condensation are also captured mathematically, as is the chemical diffusion between species. It is shown that the vaporisation and condensation state transitions in multi-component systems lead to chemical heterogenisation, primarily due to the varying vaporisation temperatures of the constituent elements in the material. A mathematical framework is presented that is capable of predicting preferential elemental evaporation, as is observed during high energy density advanced manufacturing processes.

Results
An accurate description of multi-component substrates experiencing fusion and vaporisation state changes allows the investigation of numerous physical phenomena that play key roles in determining the evolution of an alloy substrate, when subjected to intense sources of heat. The purpose of this work is to introduce such a framework, applying the governing equations to technologically relevant scenarios, in order to glean the fundamental mechanisms involved in multi-component state change. The evolution of chemical species in general is not only governed by diffusion mechanisms, their transport is strongly determined by fluid dynamics. Consequently, the composition fields will be governed by an advection-diffusion equation for each phase. Source terms are added to the advection-diffusion equations describing transformations between the condensed and vapour states. A global energy equation is utilised, determining the evolution of temperature in the domain, with a global momentum equation, that is solved to determine the velocity field in the domain. For N components in the flow, the N coupled advection-diffusion equations are solved using a finite volume approach, in a volume of fluid framework. An alternative to the volume of fluid approach for tracking the complex interfaces between liquid and gaseous states is the level set method 34,35 . The volume of fluid approach has the advantage of mass conservation, while with the level set method, the surface tension force can often be calculated more accurately 19,36 . In this work, the volume of fluid method is adopted as the primary concern is predicting the chemical heterogenisation during multi-component vaporisation and condensation. In order to ensure that the surface tension force is accurately represented in the model, it is ensured that the mesh is sufficiently refined such that the interface falls over at least three computational cells. The following multicomponent transport equation is postulated, and used to describe the evolution of the components, α k , in the flow mixture as a function of velocity vector field, U, and temperature, T where δ kl = 1, if l is the vapour counterpart to phase k, and δ kl = 0 otherwise. Here, the local phase volume fraction, α k , refers to a component of the flow, of which there are N components. In this work, diffusion between phase pairs α k ; α l f gis accounted for by a constant diffusivity D kl , with summation over all possible phase interactions, in this way heterogeneous diffusivities may be considered. A mathematical description of the diffusive fluxes between multiple chemical components, and the influence of chemistry on chemical mobilities, can be accounted for through the Onsager linear theory. D kl may be considered as a simplification of this theory for an ideal fluid, where the summation of the volumes of the components is constant (where no condensation or vaporisation state changes are occurring). It is further assumed that solid state diffusion is negligible relative to diffusion in the liquid and vapour states. The fluid fraction, ϵ 1 , and the vapour fraction, ϵ 2 , are governed by the temperature, which is determined by an energy transport equation. The mass transfer rates, _ m α evap and _ m α cond , for a given component α k refer to the vaporisation and condensation rates, respectively, depending on the local temperature of a given component in the flow. The volume dilation rate, v Á , is obtained from the mass transfer rates of all the components. In this work, it is assumed that a density change occurs during vaporisation and condensation for a component in the flow; although it is assumed that the density is not a function of temperature, apart from a buoyancy term, in any other way. The third term on the left hand side of the transport equation is an interface compression term that operates at the liquid/vapour interface. The evolution of the velocity field in the domain is governed by the Navier-Stokes equations, and a continuity equation with ∇ ⋅ U = v Á . Additional terms are added to the traditional free-flow Navier-Stokes equations governing the momentum dampening due to solidification, and the temperature-dependent surface tension phenomenon known as the Marangoni convection. The latent heat effects due to the fusion and vaporisation state changes are incorporated into the energy transport equation. Individually, the governing equations are nontrivial to solve, and hence the highly coupled nature of the advection-diffusion, energy transport, and momentum equations creates a uniquely challenging problem. The presented framework, describing multi-component metallic flows experiencing fusion and vaporisation state changes, is applied to the well-posed problem of a binary vapour bubble condensing and rising through a heterogeneous binary liquid, before being applied to a scenario representative of power beam welding manufacturing processes. The vapour bubble collapse scenario is both technologically relevant and also provides validation of the system of partial differential equations that form the framework, demonstrating that the mass of the system is conserved through the condensation state transformation, a key finding. It is not clear that other approaches to similar problems appropriately conserve mass. The power beam welding of 316L (designated by the American Iron and Steel Institute) stainless steel is then investigated, where the substrate material is treated as a mixture of the four most abundant metallic phases in a typical 316L steel for Nuclear power generation applications. In this scenario, the mechanism by which chemical heterogenisation may occur through the vaporisation and condensation state changes, is described. The thermo-physical properties used in all of the simulations are detailed in Table 1.
The surface tension, σ, temperature-dependent surface tension coefficient, dσ dT , mass diffusivity, D, and interface sharpening parameter, C α , between the phases utilised in this work are shown in Table 2. It should be noted, that surface tension parameters, namely the surface tension and temperature dependence of this surface tension, are implemented as a weighted average of the phases present in the interface of the mixture, as described in Eq. (5). It is believed that this is a good approximation, given that this approach will reproduce the pure substance behaviours at the extreme, and the interpolation between these behaviours in the mixture states. During welding, certain surface active elements, such as sulfur can drastically influence the flow dynamics, to a Table 1 Thermo-physical properties, including base S.I. units, used in the thermal fluid dynamics simulations: density, ρ, kinematic viscosity, ν, thermal conductivity in the solid phase, k solid , thermal conductivity in the fluid phase, k fluid , specific heat in the solid state, c p solid , specific heat in the fluid state, c p fluid , melting temperature, T melt , vapourisation temperature, T vap , latent heat of fusion, L fus , latent heat of vapourisation, L vap , molar mass, m, thermal expansion coefficient, β.
The thermal conductivity's for the molten state, kfluid, of Ni, Mo, and Cr are estimated, as little data for these values is present in the literature 60-63 . Table 2 Interface and diffusion coefficients between the phases used in the thermal fluid dynamics simulations: surface tension, σ, temperature dependence of surface tension, dσ dT , diffusivity coefficients, D, and interface compression values, greater degree than other species. In this work, these elements are neglected due to their low concentration in the substrate considered.
Binary vapour bubble collapse in heterogeneous binary molten liquid. Initially, the behaviour of a binary metallic vapour phase, entrained within a heterogeneous binary metallic liquid phase was investigated. The advantages of investigating such a scenario are twofold: the scenario is interesting within its own right, and meaningful conclusions can be drawn on the behaviour of multi-component vapour phase behaviour during thermo-capillary collapse; while a secondary benefit of such a system is that the problem is well posed, with clearly defined final conditions for the volume fractions of the components in the mixture at the end of the simulation. Therefore this scenario can be used as a numerical validation case for the mathematical framework.
In this scenario, the collapse of a nickel/chromium (Ni/Cr) vapour bubble, in an iron/molybdenum (Fe/Mo) liquid mixture was investigated. A cuboid vapour bubble, comprising 50% Ni vapour and 50% Cr vapour, was initialised in liquid metal, comprising 50% Fe and 50% Mo, under an argon atmosphere. Figure 1 shows the initial geometry and domain for the binary vapour collapse scenario. Figure 2 shows the time evolution of the domain as the binary bubble rises against gravity, and collapses, through the liquid metal column until all of the vapour components are transformed into their respective liquid states.
Following the start of the simulation, the Ni and Cr vapour phases begin condensing rapidly into their liquid phases, as the temperature at the vapour/liquid interface falls below their respective vaporisation temperatures. As the vapour bubble rises through the liquid column, the surface tension between the vapour and liquid states minimises the surface area of the bubble. As the density of the Ni and Cr vapour phases is three orders of magnitude lower than that of their liquid phases, a volume dilation occurs, conserving the mass of the system. The effect of the volume change induced during the condensation can be seen in Fig. 2; the total volume of the metallic phases reduces and Ar enters the top boundary to compensate (permitted by an open container boundary condition for U, p, and all α k phases). Figure 3 shows the volume fraction evolution of the components through the simulation.
As can be seen in Fig. 3, during the initial stages of the bubble rise; Ni condensation dominates over Cr condensation, as the vaporisation temperature of the Ni is higher than that of Cr. The condensation rate of Ni vapour to Ni is greater that that of Cr The red region in the domain contains 50% Cr vapour and 50% Ni vapour, the green region contains 100% Ar, and the remainder of the domain contains a homogeneous mixture of 50% Fe and 50% Mo. The temperature of the liquid metal, and vapour regions, at t = 0 s is 2400 and 3000 K, respectively The domain was uniformly discretised into 512,000 cells with a grid spacing of 2.5 × 10 −4 m. vapour. As the bubble rises and cools further, the rate of condensation for both vapour phases increases, largely due to the ever increasing surface area to volume ratio of the bubble. Eventually, the vapour bubble fully condenses, after finally colliding with the boundary as shown in the included video of this simulation, and all Cr vapour and Ni vapour is fully transformed into Cr and Ni, respectively. At this point the divergence of the velocity field becomes zero.
As can be deduced from Fig. 3, with knowledge of the vapour and liquid densities, the mass of all components in the domain is conserved through the volume constriction due to the condensation, and the state transition is fully captured. It is pertinent to look at the distribution of the elements present in the vapour bubble, as the simulation proceeds. Figure 4 shows the volume fractions of the Ni and Cr vapour at t = 0.2 s and the Ni and Cr liquid fractions at t = 0.435 s, respectively.
At t = 0.2 s almost all of the Ni has condensed out of the rising vapour bubble, and the remaining vapour is almost 100% Cr as shown in Fig. 4a, b. This difference in the vaporisation temperatures of the components within the vapour bubble generates a considerable amount of chemical heterogeneity. Figure 4c, d shows the volume fractions of Ni and Cr at 0.435 s, following complete collapse of the vapour bubble. It can be seen that as the Ni condensation initially dominates, there is a greater volume fraction of Ni at the base of the fluid column relative to Cr. If the simulation were allowed to proceed, the advective and diffusive terms in the multi-component transport equation would homogenise these distributions. However should a solidification front trail the condensation field, one could predict that this heterogeneity would be frozen into the domain.
Multi-component metallic power beam welding. The presented framework is now applied to a technologically and societally important advanced manufacturing scenario intended to be representative of power beam welding application, such as laser welding or reduced pressure electron beam welding, used to join metallic components. An area where understanding the complex heat and mass transfer in welding applications is particularly important, is in nuclear power generation. As such, in this work, a stainless steel substrate is chosen for investigation during the power beam scenarios. The multi-component description of the substrate permits alloy systems to be properly investigated; solving the coupled sets of component transport equations means that chemical mixing and preferential vaporisation of the lighter components can be modelled. In this scenario, a domain is initialised with a 10 mm thick metallic substrate, with a chemical composition intended to be representative of nuclear grade stainless steel, namely 2.5% Mo, 16% Cr, 12% Ni, and 69.5% Fe. Figure 5 shows the simulation domain used in this power beam welding scenario. The temperature of the entire domain is initialised to 323 K, therefore, the substrate is in its solidified state. As stated in the 'Methods' section, the positiondependent properties are obtained from a linear interpolation of the phase fractions present. Given the composition stated for the 316L substrate, such an interpolation would yield a thermal conductivity  of~85 kg m s −3 K −1 for the metallic substrate. It is known that the thermal conductivity of 316L is~15 kg m s −3 K −1 . To remedy this discrepancy between the known conductivity and interpolated conductivity, in this scenario the thermal conductivity of the metallic phase is set to 15 kg m s −3 K −1 such that the linear interpolation of the phase fractions reproduces the known thermal conductivity of the substrate. This adjustment of the conductivity is required in the absence of a properly coupled thermodynamic database. Note that the other physical properties do not suffer the same degree of deviation between the linearly interpolated estimations from their components and the known values for 316L. The density and thermal conductivity of the Ar phase was also reduced to 0.1 kg m −3 and 0.0017 kg m s −3 K −1 , respectively, to reflect the lower density, partial vacuum conditions usually employed during these power beam applications.
Starting at t = 0 s, a high energy density heat source, as described by Eq. (16) with radius a hs = 2 × 10 −4 m, is applied at the interface between the liquid (or initially solid) metallic substrate and the atmospheric/vapour phases. The heat source is stationary in the x-z plane. This surface Gaussian heat source initially causes melting, and subsequently vaporisation of the metallic substrate; generating a thermo-capillary through the substrate, and permitting the heat source to propagate through this thermo-capillary and effectively 'dig' through the domain; the heat source location is updated every time step, ensuring the heat is always deposited at the interface between the condensate and vaporised/shielding gas phases. As the thermo-capillary develops, progressively more of the substrate material vaporises. The volumetric change of the component elements into their vapour states (all vapour phases are assumed to have densities three orders of magnitude lower than their condensed counterparts in this work, as shown in Table 1) causes extremely large magnitude velocity fields inside the thermo-capillary, as the vaporised material escapes the substrate, of hundreds of metres per second. The recoil effect of the vaporising material has a stabilising effect on the molten material on the keyhole walls, and aids in maintaining the keyhole structure. In this simulation, the thermo-capillary fully penetrated the substrate material in 1 × 10 −2 s. Figure 6 shows the time evolution, at a cut plane of z = 0 m, where the Ar phase is represented as transparent.
The thermo-capillary that is formed has various forces acting within it that either stabilise or destabilise. Recoil velocity of the vaporising material helps to maintain the metallic liquid on the thermo-capillary walls with the help of surface tension, which acts to minimise the surface area of the liquid interface, and Marangoni forces (due to the temperature dependence of the surface tension coefficient); while gravitational body forces act on the fluid and exert a pull vertically downward towards the base of the thermo-capillary. Once the heat source is terminated, following complete penetration of the domain, at t = 0.01 s, the velocity field dramatically reduces in magnitude. Without the recoil effects of the vaporising material maintaining the liquid metal on the thermo-capillary walls, the molten material flows down into the base of what was the thermocapillary. In Fig. 6d, the vapour plume is still condensing, depositing latent heat of vaporisation into the substrate. As the vapour phases condense their volume decreases, as per the assumed density difference between the liquid and vapour states; this leaves a considerable void above the base of the weld region. This void is often observed at the termination of electron beam welds.
In the presented framework, the vaporisation and condensation rates of individual given chemical components in the substrate are treated independently through the coupled advection-diffusion equations for each component (Eq. (1)) with the vaporisation and condensation rates, a strong function of the components super-heating and under-cooling about their respective vaporisation temperatures. This permits the effect of chemical heterogenisation, due to the different state transition temperatures, and rates, of the constituent metallic components, to be captured implicitly in the system of partial differential equations. While the goal of this research is not to fully predict this heterogenisation in power beam applications, it is an interesting output from the framework. Figure 7 shows the volume fraction of the Cr and Fe phases at the base of the substrate, following complete vapour collapse and solidification at t = 0.026 s.
Further investigations into the preferential elemental evaporation are required, but the presented framework is a suitable vehicle by which this behaviour may be predicted numerically.

Discussion
A mathematical framework that describes fusion and vaporisation state changes in multi-component metallic substrates, permits the investigation of a wide variety of important phenomena with unprecedented fidelity. Alternative frameworks in the literature employed to investigate these types of phenomena have considered single-component metallic systems, and have largely neglected the volumetric changes associated with vaporisation and condensation. In this work, utilising a multi-component transport description, it has been possible to incorporate any number of metallic phases; permitting the simulation of dissimilar metal welding and additive manufacturing scenarios. It is now possible to consider the vaporisation and condensation for each species explicitly. The treatment of the volume dilation, v Á , means that the recoil at the fluid/vapour interface, which determines the thermo-capillary stability is implicitly captured in terms of vaporisation and condensation events. This framework is a step change in the understanding, and prediction of multi-component vaporisation and condensation phenomena, and is crucial to understand the complex dynamics governing advanced manufacturing processes, indeed any process influenced by vaporisation and condensation events. The proposed framework captures all of the dominating physics in multi-component flows experiencing state changes; with such a model, the most important input becomes accurate knowledge of the physical properties. Application of the model with incorrect physical properties will cause the predictions of the framework to deviate from reality. A volume of fluid approach has been adopted in this work over a level set approach for the interface tracking. While the volume of fluid approach is preferred over level set for mass conservation purposes, it is known that the level set approach has the potential to capture the surface tension effects of the flow better. It is possible to couple the volume of fluid and level set techniques, and this will be investigated further in the future. This coupling between the volume of fluid and level set techniques may also permit more accurate descriptions of the dynamic Knudsen layer that exists between the vaporising and condensing phases. The physical properties of the flow are determined from the properties of the pure substances through a simple weighted linear interpolation. For the majority of properties, this gives very good estimates. For example, the density of the 316L substrate used in the power beam welding application is recorded in the literature as 7970 kg m −3 . The linear interpolation used in this work, for the 316L substrate with 2.5% Mo, 16% Cr, 12% Ni, and 69.5% Fe yields a density of 7948 kg m −3 , a difference of~0.2%; this is clearly a very good first-order approximation. However certain properties do not agree as well when found by interpolating the pure substance values, for example, thermal conductivity. In the power beam welding scenario presented in this work, the thermal conductivity of the metallic phases is adjusted such that the thermal conductivity of the substrate matches the known thermal conductivity of 316L. In the future the presented thermal fluid dynamics framework will be coupled to a thermodynamic database for properties that cannot be accurately found through simple interpolations of the pure substance values. The framework assumes that the mass diffusivities between pairs of components in the simulation domain is constant, although it may take different values for different phase pairs. In reality, the diffusivity of elemental species is a strong function of temperature, and diffusion also occurs in the solid state at elevated temperatures.
In this work, the rate of change between a condensed and vaporised phase, i.e., Ni and Ni vapour, _ m α evap and _ m α cond was found by integrating the Clausius-Clapeyron equation around the vaporisation temperature. Other studies have used an empirical rate-parameter-based approach with the similar form for condensation, where r is a rate constant 37 . The form of the rate change equations does not have a significant effect on the overall behaviour of the system, provided that the rate is sufficiently large that the transition occurs rapidly; suggesting that the dominating factor governing the state change rate in the presented framework is the heterogeneous vaporisation temperatures of the components present in the flow.
Considering the well-bounded problem of a binary vapour bubble rising in a heterogeneous binary liquid, it is shown that the different condensation rates, and temperatures, of the species present induce heterogeneity in the final species distribution. Initially the Ni species, with a higher vapourisation temperature, prefferentially condenses from the Ni-Cr binary vapour. As the vapour bubble continues to cool and therefore decrease in volume due to the transition of the vapour species to their condensed counterparts, eventually the bubble is almost entirely formed of  Cr vapour. As the volume of the bubble decreases the condensation rates of the species increase, due to the increased surface area to volume ratio of the rising bubble, enabling more effective heat transfer at the interface. It is shown that the presented framework properly conserves the mass of all components through the state transitions. Next, the framework is applied to the industrially important scenario of power beam welding of a 10 mm thick nuclear grade steel. Both the predicted velocity fields, and the time taken for complete penetration of the domain are in good agreement with values measured in the literature 34,[38][39][40][41][42] . During vaporisation, the three orders of magnitude increase in volume of the metallic phases creates a large velocity normal to the fluid interface (captured through the v Á term in the continuity equation); this high velocity ejects the vapour plume, simultaneously creating a recoil effect on the liquid metal that stabilises the thermo-capillary against collapse. The predicted vapour plume ejection velocity is in good agreement with experimental and numerical studies in the literature 38,42,43 . The maximum velocity recorded in the present simulation was 124 m s −1 . Similarly, the time taken for the heat source to penetrate the entirety of the 10 mm domain is in good agreement with penetration times measured experimentally using x-ray transmission methods [39][40][41] . The components of the substrate do not vaporise homogeneously. As the temperature at the heat source location increases, the elements with the lowest vaporisation temperatures will tend to evaporate preferentially. The latent heat absorbed during this state transition has to be overcome before the temperature can increase further and vaporise the element with the next highest vaporisation temperature. This pattern repeats, with sufficient energy, until all components vaporise. With insufficient energy, some components will remain in their condensed state, and these tend to be the elements that are most difficult to vaporise due to their respective latent heat, molar mass, and vaporisation temperatures. It is shown that Cr is depleted in the region of the power beam weld in the presented case, which is in agreement with experimental observations 44 .
In summary, a complete thermal fluid dynamics description of multi-component alloy systems experiencing fusion and vaporisation state changes is presented, where the vaporisation and/or condensation of any number of chemical species, may be considered. Incorporated into the framework are the volume dilation and constriction due to vaporisation and condensation, respectively, the surface tension and temperature-dependent surface tension between phases, species diffusion in the fluid phases, buoyancy, latent heat effects, and a representative heat source description for power beam type advanced manufacturing processes.

Methods
Simulation of multi-component vaporisation and condensation phenomena utilising the volume of fluid approach requires the solution of the governing momentum, energy, and phase equations; these equations are strongly coupled to one another. In the presented model, the continuity equation is modified with a source term, v Á , representing the volumetric dilation/constriction per unit volume due to the vaporisation/condensation phase change, respectively, as shown in Eq.
(2), where U is the velocity.
Mathematical framework. It is assumed in this work that the densities of the liquid and vapour states is constant, and therefore the flow is approximated as incompressible in the liquid and vapour states, when no vaporisation or condensation is occurring with v Á ¼ 0, and therefore ∇ ⋅ U = 0. The regions of the domain where v Á is nonzero are where condensation and/or vaporisation state change is occurring 37,45 . The momentum equation is formulated as shown in Eq.
(3) with additional terms introduced for physical phenomena, such as the freezing of the flow field during solidification, and buoyancy effects through a Boussinsq approximation. The momentum equation is formulated accounting for varying fluid density and viscosity 36,37,45 .
The density and dynamic viscosity are given by ρ and μ, respectively. g is the body force due to gravity (modified by buoyancy term β according to the Boussinsq approximation). A viscous stress, τ, assuming laminar behaviour is considered in the present work; the form of which is shown in Eq. (4).
Surface forces are accommodated in the F s term; these include the normal surface tension, and the temperature-dependent surface tension (Marangoni) forces. Note that the recoil pressure is implicitly accounted for in the v Á volume dilation term.
where σ kl ¼ σ 0 kl þ dσ kl dT T À T m ð Þis the surface tension coefficient between phases k and l corrected for the relative temperature dependence. A Darcy source term, S m , is used in Eq. (3) to capture the momentum dampening in the 'mushy' zone and solidified regions according to the fluid fraction, ϵ 1 , using the Carman-Kozeny equation 46,47 , where K 1 and K 2 are the permeability coefficients 7 . This momentum porosity function ensures that the velocity remains negligible in the regionsm where the substrate is in a solid state. ϵ 1 and ϵ 2 represent the fluid fraction and vapour fraction, respectively. For non-isothermal state change, the fluid and vapour fractions are described as linear functions between the melting and vaporisation bounds, respectively.
where T s and T l are the lower and upper bounds for melting, respectively, and are calculated based on the phase fractions of the metallic phases present. In this work, a simplified approach is used to estimate the solidus and liquidus temperatures for the multi-component flow. This involves linearly interpolating the melting temperatures for the pure substances present in any given region, weighted by the fraction of phase present, and applying a temperature range of 20 K with the solidus temperature taken as the lower bound, and the liquidus as the upper bound of this range. In this manner, the composition-dependent solidus and liquidus may be approximated to first-order accuracy without the need to couple to a thermodynamic database. HðT À T α vap Þ is the Heaviside step function about the vaporisation temperature for the particular phase, α k , T α vap . For a fully vaporised flow ϵ 2 = 1, while for a fully condensed flow ϵ 2 = 0. The phase equation is given by Eq. (10), and incorporates diffusion between phase pairs and phase change between the condensed and vaporised chemical species states, for example, between Fe liquid and Fe vapour, and vice versa 48 .
In Eq. (10), α k represents the fraction of phase k. An additional compression velocity term is applied at the interface between the immiscible components. No geometric interface reconstruction or tracking is performed; rather, a compressive velocity field is applied in the vicinity of the interface. The value for the interface 'compression velocity' U c is given by U c ¼ min C α kl jUj; max jUj ð Þ ∇α j∇αj . The magnitude of the velocity is used since dispersion of the interface (which is being counteracted by the compression velocity) can only occur as fast as the magnitude of the local velocity, in the worst case. The coefficient C α kl is the primary means of controlling the interfacial compression between phases k and l. C α kl has a maximum physically permissible magnitude of 1 (refs. 35,49 ). For a given phase pair interaction α k ; α l f g, C α kl ¼ 0 imposes no interface compression. In the presented work, interface compression is performed between the dispersed/vapour and condensed/liquid phases, similar to the model proposed by Wardle et al. 49 .
Diffusion is included in the framework, with different diffusion constants permitted between different phase pairs. D kl is the diffusion constant between phase pairs in the interaction α k ; α l f g. The effective diffusion constant D for phase α k in a given region is found by considering all the phase interactions through a summation; in the same manner, as the compression term, these double summation operations are omitted for clarity and readability. In Eq. (10), _ m α evap and _ m α cond are the phase mass transfer rates between phases k and l due to vaporisation/condensation state changes, i.e., the rate at which a condensed phase, say liquid Fe, vaporises into its dispersed phase, say Fe vapour, or vice versa. The volume dilation/constriction term v Á is then calculated from these m Á rate terms. v Á has units of inverse time, v Á is negative for condensation and positive fore vaporisation, as shown by Eqs. (11)- (13). The vaporisation/condensation rates, _ m α evap and _ m α cond , are obtained from the Clausius-Clapeyron equation and the Herts-Knudsen relation, in order to correlate a surface mass transfer flux to a volumetric flux over a representative length scale above the interface, b d (refs. 50,51 ). _ m α evap and _ m α cond are therefore obtained from a mechanistic model, with a physical basis. Strictly speaking, the derivation of the mass transfer equations assumes the transfer flux is integrated above a flat surface, this assumption is relaxed in the current study, as the interface is not necessarily flat; however, the sensible mechanistic form of the equations provides confidence in their applicability. The transfer rate terms are functions of the difference between the local temperature and the vaporisation temperature of phase α k and α l , where α k is the condensed (liquid or solid) state of a mixture component and α l is the vaporised state of that same chemical component 37 . Note that in this framework the vaporisation temperature of α k is equal to that of α l , that is to say the phases transform at the same temperature, T α vap , during vaporisation and condensation.
where R is the gas constant. M α and L α vap are the molar mass and latent heat of vaporisation of the component, respectively, where ρ k and ρ l are the density of the condensed and vapour states. It is known that the Knudsen layer, between the liquid and vapour phases during evaporation and condensation is a highly dynamic layer with significant shear stress [52][53][54] . In this framework, the description of this layer is accounted for in the mass transfer rates described in Eqs. (11) and (12), which consider the Hertz-Knudsen formula in their derivation. At the resolution considered in this framework, where the domain is treated as a continuum, it is believed that this description is adequate. A simple energy transport equation is utilised, that incorporates the latent heat effects of fusion and vaporisation. Viscous dissipation and pressure transport terms are assumed negligible compared to the large convective and conductive heat fluxes, typically encountered during high energy density advanced manufacturing processes, as well as the large latent heat fluxes during to the state changes.
where c p and k are the specific heat and thermal conductivity of the multicomponent mixture, respectively. S h accommodates the heat contribution from the fusion and vaporisation state changes L fus is the latent heat of the melting and fusion transition 6,55-58 . q represents the rate of heating term due to the incident power beam heat source. In this work, a surface Gaussian distribution, q, is used where the heat is deposited at the interface between the dispersed (vaporised) and condensed phases.
where a hs is the radius of the distribution, Q is the input power. b g is the offset of the heat source reference frame in the x-axis. In this coordinate system, the heat source distribution is applied normal to the interface of condensed/dispersed phases in the y-direction. In power beam welding scenarios, the high energy density source of heat which causes vaporisation state changes, either laser or electron beam, is applied at the interface between the condensed and vapour/atmospheric phases. At this interface, the photons of laser light are absorbed, or the electrons transfer their kinetic energy to the substrate, causing intense local heating. It is important then, if one is to try to capture the important physical phenomena caused by such a heat source, that a representative mathematical description of such a surface distribution of heat is incorporated. The Gaussian surface distribution of heat assumed in this work, is coupled with a computationally efficient heat source application algorithm, to ensure that as the surface evolves, the heat is applied appropriately at the condensed state boundary. This involves a firstorder ray tracing approach, where a collimated Gaussian distribution is applied at the interface, at the first cell in the depth direction. In this approach, it is assumed that the energy reflected following the initial deposition of heat is negligible.
Numerical Implementation. The thermodynamics and transport properties are computed in each computational cell using the volumetric phase fraction as a weighting factor. The linear average is used for density and harmonic interpolation for the conductivity's. Specific heat capacity and kinematic viscosity are linearly weighted by mass 59 . As well as being functions of phase fraction, the specific heat and thermal conductivity's are also functions of fluid fraction in the condensed state. In this manner, the thermal properties of the alloy system can be considered for the solid, liquid, and vapour phases.The governing equations are discretized using the OpenFOAM® library following the finite volume approach. The equations are evaluated in an iterative manner, correcting the velocity, and pressure fields until convergence is reached, using a PISO (pressure implicit with splitting of operators) approach for multiple iterations at each time step, combined with a semi-implicit method for the pressure linked equation (SIMPLE), known as the PIMPLE algorithm. During the formation of the pressure-Poisson equation, the v Á term must be accounted for, so the pressure equation in the framework becomes where 1

A D
are the diagonal entries of the momentum matrix equation and ∇ ⋅ ϕ is the numerical equivalent of the divergence of the velocity field for a cell. The volume dilation term, and multi-phase transport are handled explicitly, requiring extremely small Courant numbers. The multidimensional limiter for explicit solution approach is utilised for the solution of the multi-component transport equation. Figure 8 shows a schematic diagram of teh computational framework.

Data availability
The authors declare that the data supporting the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code used to generate the results presented in this manuscript is not widely available. However, any parties wishing to implement, wholly or partly, the framework presented in this manuscript are invited to contact the authors, who will make every effort to assist with implementation of the presented equations, in any discretisation scheme, where reasonably possible. Fig. 8 Schematic diagram for the solution procedure. In the presented framework, an explicit solution procedure is used for teh phase species transport equations. This is fully coupled to the momentum and energy transport equations.