Nanoparticle Tracing during Laser Powder Bed Fusion of Oxide Dispersion Strengthened Steels

The control of nanoparticle agglomeration during the fabrication of oxide dispersion strengthened steels is a key factor in maximizing their mechanical and high temperature reinforcement properties. However, the characterization of the nanoparticle evolution during processing represents a challenge due to the lack of experimental methodologies that allow in situ evaluation during laser powder bed fusion (LPBF) of nanoparticle-additivated steel powders. To address this problem, a simulation scheme is proposed to trace the drift and the interactions of the nanoparticles in the melt pool by joint heat-melt-microstructure–coupled phase-field simulation with nanoparticle kinematics. Van der Waals attraction and electrostatic repulsion with screened-Coulomb potential are explicitly employed to model the interactions with assumptions made based on reported experimental evidence. Numerical simulations have been conducted for LPBF of oxide nanoparticle-additivated PM2000 powder considering various factors, including the nanoparticle composition and size distribution. The obtained results provide a statistical and graphical demonstration of the temporal and spatial variations of the traced nanoparticles, showing ∼55% of the nanoparticles within the generated grains, and a smaller fraction of ∼30% in the pores, ∼13% on the surface, and ∼2% on the grain boundaries. To prove the methodology and compare it with experimental observations, the simulations are performed for LPBF of a 0.005 wt % yttrium oxide nanoparticle-additivated PM2000 powder and the final degree of nanoparticle agglomeration and distribution are analyzed with respect to a series of geometric and material parameters.


Introduction
Powder-based laser additive manufacturing techniques such as Laser Powder Bed Fusion (LPBF) [1] or Direct Energy Deposition (DED) [2] have been recently established as methods that allow the strengthening of metal alloys by modification of the microstructure [3]. Often, the strengthening is achieved by introducing lattice-matched nanoparticles within the surrounding matrix [1]. Another methodology used for strengthening is the introduction of exogenic dispersoids into the metallic matrix leading to the retardation of dislocation movements. Dislocations interact with the impenetrable dispersoids by the formation of a dislocation loop in between neighboring dispersoids only allowing dislocation loops with equilibrium diameter below the dispersoids interspacing to bypass the obstacle. This effect, known as the Orowan mechanism, describes that a fine dispersion results in efficient hardening [4]. Oxide dispersion strengthened (ODS) steels make use of this mechanism. The improved mechanical properties at high temperatures of ODS alloys can be related to the nature of the steel matrix as well as the composition, size, and distribution of the dispersoids (i.e., nanoparticles). Introduced nanoparticles in ODS steels are typically composed of yttrium-based oxides, exhibiting low solubility in the steel matrix and having a low potential for coarsening by Ostwald ripening [5]. In combination with the refining agent Titanium, nanoparticles with the chemical composition Y 2 Ti 2 O 7 [6], Y 2 TiO 5 [6] are formed [7] with typical diameters of 2-3 nm in the 14YWT alloy [8]. The addition of aluminum enables the formation of various Y-Al-O-based compounds (YAlO 3 , Y 3 Al 5 O 12 ) [9] leading to coarser dispersoids (∼ 15 nm) and a reduced number density in the PM2000 alloy [8], resulting in increased ductility at the cost of lower strength.
The main fabrication route for ODS steels is the powder metallurgy route [10]. However, this expensive and complex fabrication route often leads to a low fracture toughness in the fabricated part and is inflexible in the part design [11]. As an alternative, additive manufacturing techniques are capable of producing ODS materials offering high solidification rates in combination with strong Marangoni forces within the melt pool, potentially leading to a homogeneous distribution of the introduced nanoparticles [12]. The distribution of these nanoparticles is highly influenced by the additive manufacturing technique employed and the parameters selected for processing due to differences in the melt pool dynamics. However, recent studies on this topic show that it is difficult to achieve the optimum nanoparticle size according to the Orowan mechanism due to segregation or agglomeration of the nanoparticles, which in turn was found to deteriorate the mechanical properties of the part [13]. Approaches to optimize the process such as evaluating the influence of the powder characteristics [14] and process parameters [15], or alternative additivation routes such as light mixing [16], improve the dispersion of the nanoparticles; however, they still lack control over the nanoparticle size.
The size and dispersion of the nanoparticles in the ODS steels are influenced during the steps that the initial powder undergoes until the fabrication of the ODS steel. First, the nanoparticles are supported on the steel microparticles. Ball milling of the yttriumbased nanoparticles with the steel powder is the most common approach to achieve it [17,18]. While widely employed, the control of the final nanoparticle size and degree of dispersion by this methodology is limited, and the steel microparticles size and shape can also be affected [19]. To address this drawback, alternative supporting procedures have been proposed [20] such as resonant acoustic mixing [21], solid-liquid reaction [22], or colloidal dielectrophoretic deposition [23,24]. Once the nanoparticle-additivated powder is obtained, the processing technique employed to generate the ODS steel samples and the experimental parameters selected determine the evolution of the nanoparticles, their final size, and dispersion in the ODS steel [12,25]. Even techniques with a similar working principle to LPBF like direct energy deposition (DED)-both laser additive manufacturing techniques-lead to differences in the nanoparticle agglomeration during processing due to the higher cooling rates of the melt pool achieved in LPBF that favor the preservation of the nanoparticle dispersion [12]. Even though these experimental observations provide an insight into the nanoparticle behaviors, further investigations would be required to completely understand the undergoing nanoparticle capture (nanoparticles trapped interior the microstructure), enrichment (local concentration increase of nanoparticles), and agglomeration processes [26][27][28][29]. Since the nanoparticles are not accessible during processing to perform in situ measurements, the combination of simulations with the experimental characterization of the nanoparticle-additivated powder and the generated ODS steel [13,30] represents the best approach to understand and control the LPBF processing of ODS steels.
Intuitively, the simulation of the LPBF processing of ODS steels and the underlying behaviors of nanoparticles, including their drift in the melt pool, captured during re-solidification, enrichment, and agglomeration, requires the proper modeling of the physical phenomena during the LPBF process and interactions among nanoparticles. It already remains a great challenge to model the underlying phenomena due to their sophisticated and interactive nature, covering a broad range of time and length scales. Notably, the thermal/mass transfers and material transformations (including melting, solidification and evaporation) dominate on a length scale of hundreds to thousands micrometers over a few dozen milliseconds. However, nanoparticle attraction and repulsion take place in a different spatial and temporal scale, around several tenths of nanometers and within microseconds. Finally, the manufacturing of an ODS steel sample is explained in a completely different scale of centimeters and takes hours or even days [31]. Resultant morphologies also reveal themselves in a multiscale fashion [32]. In this sense, the existent simulation schemes for the LPBF process are more or less subjected to strong simplifications and segregated modeling schemes, considering either only selected aspects or with the thermal history taken as input from separate numerical approaches. Those schemes usually feature a computational fluid dynamics method to simulate the thermal-fluid coupled spatial-transient evolution, including the arbitrary Lagrangian-Eulerian [33][34][35][36] and Lattice Boltzmann [37][38][39], incorporated with another method to simulate the accompanied microstructure evolution (mostly the polycrystalline re-solidification), such as the cellular automata [40][41][42][43][44] and phase-field [45]. Recently, we proposed a new phase-field model considering coupled processes among heat transfer, melt flow dynamics, and microstructure evolution (noted as "heat-melt-microstructure-coupled processes"), which shows the possibility to simulate the LPBF process in a unified and thermodynamic consistent route and the ability to recapitulate various experimental observations via simulation, such as high-gradient temperature field, tilted columnar grains, and lack-of-fusion pores with irregular shape [46].
On the other hand, very few works have been conducted regarding nanoparticle drift and interactions in the melt pool, which is, however, the central aspect to investigate the behaviors such as enrichment and agglomeration. Xu et al. [29,47] aimed to convey fundamental understandings regarding nanoparticle capture during re-solidification in order to obtain a nanocomposite-dispersed metal bulk, bringing the consideration of the Van der Waals (abbreviated as VDW hereinafter) effect, the Brownian effect, and thermodynamic analysis into the modeling of interactions in the nanocomposite-melt dispersion system. There are also experimental works revealing the VDW [48] and electrostatic interactions [49][50][51] among inorganic nanocomposites in the liquid metal. Nevertheless, those works disregard the influence of the driving effect from the melt flow dynamics when describing the nanocomposite-liquid/molten metal dispersion as well as the inter-particle interactions. In this sense, a simulation scheme combining the coupled phenomena (especially the melt pool dynamic and microstructure evolution) and the nanoparticle behaviors is still to be developed for the investigation of LPBF processing of ODS steels.
Joining our heat-melt-microstructure (HMM) coupled phase-field model with nanoparticle kinematics, we present and apply in this work a simulation scheme for tracing the nanoparticle drift and interactions in the melt pool during the LPBF process of the ODS steels. This proposed scheme is aimed to demonstrate the chronological and spatial variation of multiple traced nanoparticles with respect to various factors (such as chemical composition and size distribution), which enables graphical and statistical analysis on the nanoparticle migration and further effects, such as nanoparticle capture, enrichment, and agglomeration. The nanoparticle compositions used in the simulation are chosen based on our experimental studies on the additive manufacturing of ODS steels, including the use of Y 2 O 3 [12,23] and YIG (yttrium iron garnet, Y 3 Fe 5 O 12 ) [23] nanoparticles. Influences from different types of nanoparticle size distribution, including the monomodal, normal, and log-normal distributions fitted from the experimental measurement [23], are also investigated and discussed in this work.

Non-Isothermal Phase-Field Model for Stable LPBF Processes
To properly simulate the microstructure of a polycrystalline material, manufactured under a stable LPBF process, a conserved order parameter (OP) ρ is employed to represent the substance and atmosphere/pores; a set of non-conserved OPs φ S and φ L is employed to represent the solid and liquid phases, respectively; and a series of non-conserved OPs {η α } are employed to represent the orientation distribution among the powders/grains. Two constraints should be applied on those OPs to properly recapture the reality (Figure 1a), i.e., substance constraint (1 − ρ) + φ L + φ S = 1 to restrict the existence of liquid as well as solid within the substance only, and the polycrystal constraint (1 − φ S ) + ∑ α η α = 1 to restrict polycrystalline orientations inside the solid substance. Notice that the substance constraint always satisfies when φ L = ρ − φ S . Therefore, only one φ to represent the solid substance and (ρ − φ) to represent the liquid substance is sufficient. The other constraint can be fulfilled via some practical numerical method, e.g., the Lagrange multiplier or the penalty method.
We consider a stable LPBF process (Figure 1a), i.e., without significant vaporization and resultant phenomena, such as scattering and keyholing. In this sense, the heat/mass transfer through various diffusion paths (e.g., surface, grain interior, and grain boundary) and melt flow, the dynamics of melt flow as well as the bubbles, the melting-solidification of the grains, and some inter-coupling effects, such as thermocapillary (Marangoni effect) and thermophoresis (Soret effect), would take the major role determining the resultant microstructure of the manufactured sample. The nonlinear kinetic system in simulating the stable LPBF processes is then adapted from our HMM-coupled phase-field model [46]. It is worth noting that the HMM-coupled phase-field model, derived under the thermodynamicconsistent framework, can be regarded as the combination of Navier-Stokes-Cahn-Hilliard (NSCH) and Navier-Stokes-Allen-Cahn (NSAC) systems with inter-coupling effects integrated. In this work, however, mentioned inter-coupling effects are tentatively dropped under the consideration of the computational stability, consumption, and complicity control in variants. These effects would be explicitly covered in our separate (e.g., the thermophoresis in [52]) or upcoming works. The adapted nonlinear kinetic system for the velocity field of the melt u, the temperature field T, and fields of OPs ρ, φ and {η α }, eventually presents as with dimensionless quantities, namely the Reynolds number (Re), the Froude number (Fr), the Stefan number (Ste), the Péclet numbers for thermal (Pe T ) and mass (Pe ρ ) transfer, and Allen-Cahn numbers for melting-solidification (Ac φ ) and grain growth (Ac η ). Here, g stands for the unit vector of gravitation direction. These quantities are employed not only to parameterize the nonlinear system but also to characterize the ratio between their corresponding physical processes to the chosen rate (by default, the characteristic rate of the fluid). Detailed parameterization of these quantities will be explained in Section 3.3.  The explicit formulation of non-isothermal free energy F is formulated as with where T M is the melting temperature, L is the latent heat of the material, and c r is the relative volumetric specific heat. Therefore, the variational derivatives of the free energy appearing in Equation (2) yield δF δρ Model parameters A, B, G, C pt , C cf , D pt , D cf , H pt , H cf , as well as the gradient constants (κ ρ , κ φ , and κ η ), are obtained from given diffusive interface width and the experimentally measured interface energies. Coefficient ξ is employed to favor the determination of model parameters by fitting the experimental results (sufficiently summarized in Ref. [46] and supplementary information of Ref. [53]). Φ L is the interpolation function indicating the spatial landscape of the liquid/melt (see Section 3.3).
Finally, the thermal effect is equivalently treated as an internal heat source term q v moving with the scan velocity v l in which P l is the laser beam power reaching the surface of the powder bed, β is the attenuation coefficient. Φ ss is the interpolation function for the substance.
x is an arbitrary point on the projected plane of the laser beam on the powder bed surface, while z is an arbitrary depth from the plane. (v l t, z v ) is the moving center of the beam following the morphology of the surface. ζ is the characteristic penetration depth and normally takes the value of the powder bed thickness. Notice that parameter Π is utilized to adjust the concentration of the deposited power inside the circular beam spot with nominal radius R l , e.g., as suggested by the ISO standard [54], Π = 2, indicating 86.5% of the concentrated power within the spot. P l , R l , and the mode of the scan velocity v l = |v l | (scan speed) are thereby regarded as the major processing parameters of the laser scan in this work.

Nanoparticle Kinematics
The coupled evolution among polycrystalline microstructure, melt flow dynamics, and temperature transfer is calculated on the mesoscale (0.1-100 µm) using the phase-field model presented above. Since the size of nanoparticles is smaller by several orders, the possible influence of nanoparticles on these mesoscale effects is ignored. As an important part of nanoparticle kinematics, the drift effect of melt flow is inherited from the phase-field simulations. Moreover, nanoparticles can interact with each other by different mechanisms. Based on former research for dispersed non-metallic particles [29,47,[55][56][57] neglecting the rotation, the kinematic equation for a dispersed rigid nanoparticle (labeled as i) with a volume V i , the density i , and the translation velocity v i in the melt, in which the r.h.s. terms are the driving force densities due to the gravity, the melt flow, the VDW, and the electrostatic interaction, respectively ( Figure 1b). Considering the boundary of the nanoparticle as ∂V i , the force density due to meltflow driven is formulated as f M = 1 V i ∂V i σ · n dS, in which the Cauchy stress tensor of the fluid reads σ = −pI + 1 Re ∇ 2 u. This is widely employed in simulating nanoparticle dispersion, including lattice-Boltzmann [57,58] and finite element method using static [59] or fictitious domain [60,61]. Assuming an infinitesimal size of the nanoparticle compared to the characteristic length scale of the melt, i.e., V i → 0, f M can be thus calculated as follows in this work according to the definition of the divergence where x i is the center point of the nanoparticle, and Du/Dt| x i is the undisturbed time difference of the melt velocity at point x i . The VDW and electrostatic force densities, depending on respective inter-particle potential U A and U E , can be formulated as where d C is the inter-particle center distance, Figure 1b. ∇ d represents the gradient operator in the inter-particle direction. Derived by Hamaker [62], the inter-particle potential of VDW attraction between nanoparticles i and j with corresponding radius r i and r j is where A p:m is referred to as the Hamaker coefficient of the nanoparticle in a medium, which is dependent on the permittivity of involved materials and the intervening medium. According to Dzyaloshinskii-Lifshitz-Pitaevskii interpretation [63,64], in which nanoparticles are treated as continuous media (rather than atomic structure as [62]), and inter-particle forces are derived in terms of permittivities and refractive indices, the non-retarded (instantaneous) Hamaker coefficient for VDW interactions between two nanoparticles of the same material through a medium (denoted as A p:m ) is estimated as with the permittivities ε p and ε m , and the refractive indices n p and n m of the nanoparticle and the medium, respectively. k B andh are Boltzmann constant and reduced Planck constant, while ω e is the electronic absorption frequency, ranging in 3 ∼ 5 × 10 15 Hz. It is worth noting that Equation (8a) is not employable due to the difficulties in the practical measurement of permittivity and refractive index of the molten metal. As an alternative, A p:m is calculated from the ones for nanoparticle (A m:v ) and medium (A p:v ) that obtained separately in the vacuum, according to the following combining relation [48,65] where A p:v can be obtained from experimentally measured optical data [56], while A m:v for dispersed metallic particle/droplet in the vacuum is derived from Dzyaloshinskii-Lifshitz-Pitaevskii theory by taking the estimated frequency-dependent permittivity for metals as ε m (ω) ≈ 1 − ω 2 e /ω 2 with two critical conditions, i.e., ε m (0) = ∞ and ε m (∞) = 1, which is valid for plasma and metal. Then, A m:v eventually yields Obviously, A p:m obtained from either Equations (8a) and (8b) is always positive with a typical value ranging in 10 −19 ∼ 10 −20 J, demonstrating an ever attractive VDW forces between two nanoparticles of the same material through a melt (taking negative value as attraction and positive one as repulsion for all nanoparticles if without further interpretation).
Unlike in the aqueous solution, electrostatic repulsion in the molten metal has rarely been studied. Although there are works [29,47] assuming negligible electrostatic interaction due to the strong screening effect of background free electrons in the melt, there is reported evidence of charged surface on the oxide nanoparticle in the liquid metal [49,50]. Due to accompanied electron defects (like quasi-free electrons or holes) in the oxide, electron flow occurs across the interface when the oxide nanoparticles are in contact with the liquid metal, resulting in the surface potential ψ i and corresponding profiles at the interface [49,51], as shown in Figure 1c. Combining the above viewpoints from existing studies, the following assumptions are made in the work to model the electrostatic interaction between nanoparticles through melt:

1.
Electrons in the molten metal behaves as a free electron gas with a density n e , receiving only the contributions from valence electrons of all metallic elements in the melt.

2.
Charged surface exists on the oxide nanoparticles in the melt, even though its effect is weak or absent due to the strong screening effect. The value of this surface charge ψ i is assumed to be equal to that obtained from a neutral aqueous dispersion.
Then, the screened-Coulombic potential between two nanoparticles, applied for electrostatic interaction in the electron-screened system (like plasma), is adopted from [66], i.e.
where the interaction coefficient Z is analogous to the Hamaker coefficient, reading as which also presents a positive value for nanoparticles of the same material in the melt, demonstrating an ever repulsive electrostatic interaction. λ P is the plasma Debye length, which depends on the free electron density n e , i.e., where z l , c l , and m l are the valence electron number, the atom fraction and the atom mass of the element l, respectively. N A is the Avogadro constant, ε 0 is the vacuum permittivity, and e 0 is the elementary charge. Taking the formulation of the interaction potentials in Equations (7) and (9), the VDW and electrostatic force densities are formulated as According to its physical meaning, λ P characterizes the strength of the screening effect due to the electrostatic interaction received by one nanoparticle, which is reflected by the decaying length scale of the f E , as shown in Figure 1d. It is calculated that for the alloy such as PM2000, λ P is at the scale of 10 −3 nm, indicating a very short-range electrostatic repulsion. In addition, f E trends to the ordinary Coulombic potential f Coul E with a finite value when two nanoparticles are getting in touch, while f A tends to infinity. This demonstrates that the attractive VDW interaction dominates the resultant effect between two nanoparticles, which would destabilize the dispersion system from a perspective of colloidal science [29,48,65] and lead to agglomeration of the nanoparticles. This is against the preference of forming a homogeneous oxide dispersion during the LPBF process. Therefore, the validation of the proposed model relies on further experimental insights into the interactive behavior of the nanoparticle in the molten/liquid metal, which would be covered in our future works.

Numerical Scheme and Implementation
Assuming negligible counter-effect of the nanoparticles on the melt flow, the simulation is designed in a subsequent scheme: an HMM-coupled phase-field simulation for the LPBF process and a subsequent nanoparticle kinematics simulation with information received from the phase-field simulation, as shown in Figure 2. The information includes the initial condition (IC) of the powder bed as a set V(X I , R I ) with the center X I and the radius R I of the powder I to create IC (centers and radii) for the additivated nanoparticles, and nodal values (ρ J , φ J , {η J } and f J M of nodal J) of every time-step to provide the melt-flow driving force as well as the on-site phase information. Due to the infinitesimal-volume assumption, the drifted nanoparticles are represented and traced by corresponding center position x i . The trajectories of the nanoparticles are then calculated numerically by discretizing the kinematic model outlined in Section 2.2 in the backward differences fashion: Note here the time difference δt should be by default no larger than the time interval ∆t of the phase-field simulation to ensure the accuracy of the tracing. This tracing program repeats along with the phase-field simulations till reaching the stop criteria, i.e., the end time of simulation or interrupted due to non-converge situation during finite element method calculation. In addition to the trajectories, another important piece of information is the relative position of a nanoparticle in mesoscopic microstructure, for instance, grain interior, grain boundary region, pore, or surface. The position of a nanoparticle can change during its drift and is thus described by a chronological variable, which is termed here as "position indicator".
The phase-field model is numerically implemented via the finite element method within the program NIsoS developed by authors based on MOOSE framework (Idaho National Laboratory, ID, USA) [67]. Four-node quadrilateral Lagrangian elements are chosen to mesh the geometry. The Cahn-Hilliard equation is solved in a split way [68,69].
A transient solver with preconditioned Jacobian-Free Newton-Krylov (PJFNK) method and second-order backward Euler algorithm has been employed to solve the non-isothermal phase-field problems. Adaptive meshing and time-stepping schemes are used to reduce the computation costs. In addition, to stabilize the calculation of NSCH and NSAC systems, streamline-upwind Petrov-Galerkin and pressure-stabilized Petrov-Galerkin methods are introduced associated with the weak forms of the Navier-Stokes equations [70]. For now, the subsequent nanoparticle tracing program is coded by Python (ver. 3.7.10), which is independent of NIsoS. It is planned to integrate the tracing program as a subroutine of the phase-field simulation program in the future.

Simulation Setup
As a preliminary step, we apply first our nanoparticle tracking scheme for a 2D phase/field simulation of LPBF following our former work [46], which can recapture certain characteristic powder bed features, e.g., particles with multiple sizes and various pores due to the particle packing. The simulation domain has a size of 500 × 100 µm. Particles inside the domain are generated with the random close packing procedure. Due to the uncertainty of the initial grain structure of a single particle, we simply treat each particle as a monocrystal with a unique random orientation following the reported simulation works [46,71]. With the help of the minimum coloring algorithm and grain tracking algorithm [53,72], five η α are sufficient to uniquely represent all the particles/grains for these simulations. The zero Neumann boundary condition (BC) for ρ and zero Dirichlet BC for u, representing a close BC for the system, are applied on the boundary Γ = ∂Ω of the whole simulation domain Ω, wheren is the normal vector of the boundary Γ, and 0 is the null vector. The heat convective BC allows heat dissipation as heat convection with the atmosphere on the top boundary Γ where h is the convective coefficient and T E is the environment temperature. The Dirichlet BC on temperature with a fixed pre-heating temperature T P is applied on the rest of the boundary Γ (Γ = Γ ∪ Γ ) to emulate the contact with a semi-infinite heat reservoir (e.g., the substrate) which helps to restrain the melt pool size for better demonstration in this work.
In addition, the pinning pressure BC on the top-left corner C (C ∈ Γ ) is applied on the hydrodynamic pressure p as p| C = 0 (16) to avoid the difficulties associated with the non-trivial nullspace of the operator prespecified in the PJFNK solver as suggested in Ref. [70].

Parameters and Properties
The dimensionless quantities are employed in the nonlinear kinetic system for HMMcoupled phase-field simulations as explained in Section 2.1. Following the conventions of NSCH and NSAC system, they are formulated explicitly as wherev,l andT are the characteristic velocity, length and temperature of the melt flow, while is characteristic width of the diffusive interface corresponding to the given characteristic surface tensionγ. ν is dynamic viscosity, k is the thermal conductivity, and M is the isotropic diffusivity. L φ and L η are the isotropic Allen-Cahn mobility of the solid-liquid interface and grain boundaries. b is the magnitude of the resultant body forces acting on the melt flow. For convenience, we use a set of simpler reference quantities to re-define those characteristic quantities by substituting the following relations in Equation (17a): Notice that the C Notice here that material properties c r , k, , ν, should be phase-dependent and thereby formulated in a direct interpolated fashion as where c p (·) , k (·) , (·) , and ν (·) are respectively the specific enthalpy, thermal conductivity, density, and dynamic viscosity of the corresponding phase. Similarly, the effective value of mobility L η and M through possible paths are adopted correspondingly from the selfdiffusivities D eff (·) and grain boundary mobility G eff gb [46,53,73,74], i.e., while the quantity of L φ is tentatively given as 20/tC T M pt , which makes resultant 1/Ac φ sufficiently larger than 1/Ac η to emulate a relatively faster melting-solidifying process than grain growth [46]. This is due to the lack of a quantitative description of migration mobility of the liquid. The subscript represents the quantities of the corresponding phases, e.g., 'ss' as the substance, 'at' as the atmosphere/pore, 'sf' as the surface, 'gb' as the grain boundary, 'S' as the solid and 'L' as the liquid/melt. Then, the interpolating functions Φ ss , Φ at , Φ gb , Φ sf , Φ S , and Φ L can be simply formulated as Note that the constraints on the OPs should be also applied on the interpolation functions, i.e., 1 = Φ ss + Φ at and Φ ss = Φ L + Φ S .

Results and Discussion
Simulation results are presented here for a single scan LPBF nanoparticle-decorated PM2000 alloy in an argon atmosphere. The composition of the PM2000 adopted in this work (presented in Table 1 following our former experimental investigations [12,23]) presents a ferritic structure in the high-temperature range (from 1000 K to the melting temperature) without solid-state phase transition, according to the Fe-Cr-Al ternary phase diagram reported in Ref. [75]. Therefore, the material properties for ferritic PM2000 alloy displayed in Table 2 are employed for the simulations. The reference length scale selected for the simulation isl = 1 µm, and the time scalet = 1 µs. The LPBF parameters are selected according to the LPBF experiments in [12]. The characteristic radius of the beam is set as R l = 80 µm, the penetration depth ζ is defined by the thickness of the powder bed, and the laser power and scanning speed are P l = 160 W and v l = 800 mm/s, respectively. The attenuation coefficient of the laser is β = 0.65. The environment temperature T E as well as the pre-heating temperature T P are both set to 353 K. The single scan LPBF of the monomodal Y 2 O 3 (15 nm in radius) nanoparticle-additivated PM2000 powder is first discussed. Chronological microstructure evolution along with the traced nanoparticles and their trajectory are depicted in Figure 3. To provide a further insight into the steel and nanoparticle evolution during processing, three PM2000 particles from different locations with radii of 17-19 µm are chosen as the "parent particle" and marked as 'A', 'B', and 'C'. Fifty equispaced nanoparticles are placed on their surfaces in order to simulate a perfectly homogeneous initial nanoparticle dispersion. Parameters for nanoparticle interactions, including the Hamaker coefficient and the surface potential can be found in Table 3. As the HMM-coupled phase-field simulations [46] show, the microstructure evolution of LPBF-processed PM2000 features multiple phenomena, including melt flow convection as well as accompanied behaviors of pores/bubbles, partially melted particles, re-solidification, and sintering of particles/grains, resulting in columnar grains and trapped irregular-shape pores due to lack of fusion. Those phenomena are also highlighted in Figure 3a 1 -a 4 . The focus here is on how the nanoparticles drift and eventually migrate under the influence of those effects. Initially, all the nanoparticles are located onto their corresponding parent particles before the laser interaction and generation of the melt pool (Figure 3a 1 ). Once fully immersed by the melt, the nanoparticle drift is firstly driven by the melt flow, i.e., there is a driving force along the tangent direction of the streamline (Figure 3a 2 ,b 1 -d 1 . This effect is illustrated in Figure 1b 1 -b 4 ). Due to the sufficiently large spacing, other interaction contributions among nanoparticles are negligible. After the initial drift due to the melt front propagation, the evolution of the nanoparticles is associated with different physical phenomena depending on their location and surroundings. Nanoparticles around a pore/bubble, which are identified as the "bubble-carried" ones, would follow the floating, deforming, or even splitting of bubbles (Figures 3a 2 . It is also possible that the melt flow brings multiple nanoparticles into a narrow region at the time, which is already within the range of nanoparticle interactions (specifically the VDW attraction since the electrostatic repulsion has an even shorter range around 10 −3 nm). In this sense, the trajectory of certain nanoparticles would be redirected abruptly. Meanwhile, nanoparticles near the bottom of the melt pool present very little migration comparing to others. The reason can be the very short immersing time (thus less driven by melt-flow) or the partial melting of the corresponding parent particles. Next, when the laser front scans away, the local temperature drops and, once below T M , the resolidification occurs, forming the tail of the melt pool. Once the re-solidifying front goes through the migrating nanoparticles, they are immediately captured and become either interior or grain boundary (GB) nano-inclusions, while the ones carried by the bubbles stay as they are and become pore-trapped nanoparticles. Nanoparticles can also be found on the surface driven by either melt flow or emerging bubbles (Figure 3a 3 ). Figure 3d 2 -d 4 present an additional case where nanoparticles, initially located on the surface, migrate in accordance with the deformation of the surface morphology. During this process, the nanoparticles with a relatively higher speed (driven by the melt flow) might enter the melt pool and become the in-grain nano-inclusions, while others remain on the surface. Finally, Figure 3a 4 presents traced nanoparticles with all sorts of position indicators (denoted by colors) and, notably, the locally enriched nano-inclusions and several ones with overlapped trajectories, implying potential agglomeration effects. One can readily tell from Figure 3a 2 ,a 3 that such enrichment majorly can be attributed to the melt-flow driving force, where multiple nanoparticles follow similar trajectories governed by the transient streamlines. However, it is worth noting that the overlap of the trace markers (indicating the center locations of the nanoparticles rather than the sizes) do not sufficiently reflect nanoparticle agglomeration effects, which should be explicitly determined by an adjacency test on the real scale of the nanoparticles' size. This will be discussed in the following demonstrations.    [83] ν at ∼ 7.53 × 10 −5 (J s)/m 3 [84] * Temperature-dependent data form [78] and scaled based on the value at T M from [77]. † Data from ferrite. ‡ Linearly interpolated from temperature-dependent measurements on Fe-Cr melt with 21 at% Cr.  Figure 4a,b presents the interactive forces vs. surface distance between two nanoparticles with different radii and composition, which are Y 2 O 3 and YIG, their properties and related coefficients are listed in Table 3. Notice here the Hamaker coefficients in the melt for both Y 2 O 3 and YIG are calculated according to Equation (8b) using the coefficients in vacuum for both compositions and the melt, where A m:v = 37 × 10 −20 J is calculated by Equation (8c). It can be observed that, for the same surface distance, increasing nanoparticle radius leads to a higher VDW attraction as well as electrostatic repulsion. Comparing nanoparticles of the same size but different composition, YIG nanoparticles present less VDW attraction and more electrostatic repulsion compared to Y 2 O 3 ones. Nevertheless, the interaction range of the electrostatic repulsion is limited to the 10 −3 nm scale, and this interaction highly decreases in the nanoscale due to strong screening effects from background free electrons (reflected by very small λ P ), where the VDW attraction still takes the major role. This fact explains the instability of homogeneously dispersed nanoparticles on the 10 nm scale, as explained in Section 2.2. In the case of VDW attraction, it presents a considerable decay when d S ≥ 2 nm, hence its influence over the inter-nanoparticle attractions is low, especially when the nanoparticles are sparsely distributed. Therefore, melt-flow-driven effects are the dominant mechanism, and the effects from nanoparticles' size (radius) and density would be significant.
In Figure 4c Figure 4c 9 -c 12 , on the other hand, indicate the potential interaction between particles. In particular, the latter pair presents a merged trajectory after one abrupt redirection, implying a potential agglomeration due to a shortrange VDW attraction. A similar effect is also spotted for the case with relatively higher density (i.e., YIG) by comparing again the pair A1-A2 between Figures 4d 1 ,d 2 and pair C41-C42 between Figures 4d 5 ,d 6 , demonstrating the potential enhanced interaction in an increased size for nanoparticles of identical composition/using YIG compared to Y 2 O 3 for nanoparticles of identical size. Apart from these effects, a more complicated pattern variation in the overall trajectories change shall be discussed with respect to the change in nanoparticle size/composition. Unfortunately, information for current simulations (noting the nanoparticle are sparsely distributed) remains insufficient to deduct such a pattern change, which will be sufficiently covered in our upcoming works.  with varying radius and composition on different parent particle, i.e., for Y 2 O 3 nanoparticles on parent particle A with varying radius of (c 1 ) 10 nm, (c 2 ) 20 nm, (c 3 ) 30 nm, (c 4 ) 40 nm; on parent particle B with radius of (c 5 ) 10 nm, (c 6 ) 20 nm, (c 7 ) 30 nm, (c 8 ) 40 nm; on parent particle C with radius of (c 9 ) 10 nm, (c 10 ) 20 nm, (c 11 ) 30 nm, (c 12 ) 40 nm; and for nanoparticles with monomodal radius of 15 nm on parent particle A with composition (d 1 ) Y 2 O 3 , (d 2 ) YIG; on parent particle B with composition (d 3 ) Y 2 O 3 , (d 4 ) YIG; on parent particle C with composition (d 5 ) Y 2 O 3 , (d 6 ) YIG. it demonstrates that the pattern of trajectories changes due to varying size/composition. By comparing the selected trajectories (labeled uniformly), an enhanced nanoparticle floating is observed when increasing the size of nanoparticles.
A statistic investigation of the nanoparticle evolution during processing is conducted to deduct features such as enrichment and agglomeration. The simulations shown in Figure 5 are performed for PM2000 powder with nine selected parent particles decorated with 0.005 wt % Y 2 O 3 nanoparticles due to the limited computational efficiency, as presented in the inset of Figure 5a. Notice here that the weight percentage is calculated adapting the 2D scenario, i.e., wt % = In this sense, nanoparticles of N p = 1240 are traced simultaneously. Three types of nanoparticle size distributions are investigated: the normal and the log-normal fittings of the large nanoparticle fraction measured in [23] and a monomodal obtained from the fitted mean radiusr i = 13 nm of the distribution. To reduce the contamination to the statistics of the drift destinations, alloy powders located lower than the melt pool depth are not decorated with nanoparticles (inset in Figure 5a). Figure 5b presents the statistics of the drift destinations of the Y 2 O 3 nanoparticles tested with different size distributions. Employing the adjacency test, in which the center distances (d C ) between nanoparticles are individ-ually compared with the corresponding sum of mutual radii (r i + r j ), the agglomerated nanoparticles can be identified. Following this criterion, the nanoparticles are classified as agglomerated/non-agglomerated, showing the corresponding fraction in Figure 5b. The results show that ∼ 55% of the nanoparticles are captured inside the grains after the LPBF process. Comparing different size distributions, monomodal and log-normal types have almost the same amount of in-grain nanoparticles (nano-inclusions), while the normal type shares the same agglomeration ratio as the log-normal type. However, only ∼ 2% of the nanoparticles end up in the grain boundaries. Within the nanoparticles in the grain boundaries, the log-normal distribution presents more agglomerated nanoparticles, reaching 46% (0.6% out of in-total 1.3%). In addition, there are nanoparticles that end up in the pores, ∼ 30%, and the surface, ∼13%. It is worth noting that these statistical results with a high fraction of in-grain nano-inclusions and a very low fraction of GB-captured ones may be attributed to the lack of driving effects from the migration of various interfaces, such as melt-grain interface and grain boundaries, which might force some of the nanoparticles to translate towards the migrating directions, pushing them away from the grain interior and eventually capturing them after the fusion of the powder bed.  In addition to the results provided in Figure 5b, the nanoparticle evolution during processing is graphically displayed in Figure 5c. This methodology is not only useful to predict the nanoparticle agglomeration but can be also extended to study the nanoparticle capture, enrichment, and agglomeration. However, even the evaluation of the agglomeration requires attentions on not only the trajectories overlapping but also the vast differences in spatial scales between the microstructure and the nanoparticles. For instance, Figure 5c presents the trajectories and positions of nanoparticles (log-normal type) dispersed in the LPBF-processed PM2000 polycrystalline matrix, in which two scopes are taken and magnified to visualize the nanoparticle scale. If only the trajectories are evaluated, one might prematurely conclude that Scope 2 has profound agglomerations comparing to Scope 1. However, after a magnified look, agglomerated nanoparticles only appear in Scope 1, while Scope 2 only contains locally enriched isolated nanoparticles.
The issue of identifying the nanoparticle agglomeration/enrichment should be also addressed in the experimental observation, yet the large spatial scale difference makes it difficult to characterize both the steel microstructure and nanoparticle size and dispersion in a single measurement. An example of an experimental SEM-EDX measurement from an LPBF processed sample of PM2000 steel decorated with 0.08 wt % Y 2 O 3 is presented in Figure 5d. Four inclusions (labeled as I1-I4) that can be initially thought of as nanoparticle agglomerations during processing can be observed [12]. However, to evaluate whether the inclusions are pores or agglomerated nanoparticles, an elemental mapping with techniques such as EDX is required. The EDX maps show that, while no Y or O content is found in I3, these elements are present in I1, I2, and I4, and so it can be concluded that they are nanoparticle agglomerated structures. Furthermore, there is solely a small fraction of Y content found on the surrounding of I3 while O content is fully presented in the interior, implying agglomerated Y 2 O 3 nanoparticles on the surface of another oxide inclusion that might be attributed to impurities. Another interesting point is that the diameter of some nanoparticle agglomerations spotted in Scope 1 of Figure 5c are less than 1 µm, consisting of merely countable nanoparticles, while the ones (esp. I4 in Figure 5d) experimentally observed are larger, and even much larger if compared with the original size of the Y 2 O 3 nanoparticles-even though there are some (e.g., I2 in Figure 5d) showing consistency in scale. A possible reason is the relatively low mass fraction of the additivated nanoparticles in the tracing simulation (0.005%) compared to the experimental one (0.08%), which increases the probability of the nanoparticle interactions leading to agglomeration during processing. In future steps, the nanoparticle concentration employed in the simulation will be increased to match the experimental conditions.

Conclusions
In this work, we proposed a simulation scheme joining the heat-melt-microstructurecoupled phase-field model and the nanoparticle kinematics to trace nanoparticle during the LPBF process of the ODS steels, which is experimentally inaccessible. Simulations on stable LPBF single scan of a ferritic PM2000 nanoparticle-additivated powder bed were conducted for factors such as the nanoparticle composition and size distribution. The following conclusions can be drawn from this combined numerical and experimental study: Simulation results provide the chronological location and located phases of the traced nanoparticle. This helps to depict nanoparticle drift associated with the evolution of local melt flow as well as the morphology, such as migrating nanoparticles driven by melt-flow or carried by floating/deforming bubbles, or stationary ones in the melt pool bottom area. Events such as nanoparticle capture (by grain/pore/grain boundary), enrichment, and potential agglomeration can be also visualized via trajectories and position indicator.

2.
The drift and interactions of nanoparticles with different sizes and compositions (Y 2 O 3 and YIG) are analyzed. By comparing the trajectories and positions of selected nanoparticles (or nanoparticle pairs) among cases, some preliminary discussions can be conducted regarding the influences of the nanoparticle size and compositions.
An enhanced nanoparticle floating is observed when the size is increased for nanoparticles of identical composition or when using Y 2 O 3 for nanoparticles of identical size. In addition, an enhanced nanoparticle interaction is observed when increasing the size of nanoparticles of identical composition or using YIG for nanoparticles of identical size. Note that the above conclusions are made under the condition that nanoparticles are sparsely decorated (i.e., nanoparticles are sufficiently spaced). In this scenario, the melt-flow-driven effect is expected to dominate the process. 3.
LPBF simulations of a PM2000 powder bed are conducted, in which nine parent particles are additivated with 0.005% Y 2 O 3 nanoparticles. Three size distributions are evaluated i.e., monomodal, normal, and log-normal distributions. The results show that ∼ 55% of the nanoparticles are eventually captured by a grain, while merely ∼ 2% ones end up in the grain boundaries. Although the differences of nanoparticle location for the different size distributions are small, the monomodal case presents a relatively higher agglomeration ratio in grain-captured nanoparticles (nano-inclusions), while the log-normal type shows a higher agglomeration ratio in GB-captured nanoparticles.

4.
By visualizing the traced nanoparticles on a nanometric scale, nanoparticle agglomeration and enrichment spotted in the simulation are distinguished. Comparisons between the simulations and experimental results show promising similarities, proving the potential of the simulation methodology to optimize the LPBF process parameters in order to reduce agglomeration effects and maximize the material reinforcement achieved in ODS steels.
The proposed scheme and models should be further extended in the near future for different aspects, e.g., implementation of nanoparticle tracing code in a computationallyefficient way to enable decoration with a larger mass fraction of the nanoparticles. Further effects such as the driving forces from migrating interfaces, and interactions between nanoparticles and the surface of the parent particles (as ideally the semi-infinity large interfaces) should also be considered.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Symbols
The following symbols are used in this manuscript: ρ Conserved order parameter indicating the substance and pores/atmosphere φ S , φ L , φ Non-conserved order parameters indicating the solid (φ S ) and liquid (φ L ) phase. Due to the substance constraint φ L = ρ − φ S , only one φ is sufficient to represent the solid substance and (ρ − φ) to represent liquid substance {η α } Non-conserved order parameters indicating grain orientations u Velocity field of melt flow p Hydrodynamic pressure T Temperature field T M , T E , T P Melting temperature, environmental temperature, and pre-heating temperature, respectively f loc , f grad , f ht Local, gradient and thermal terms of Helmholtz free energy density A, B, G Model parameters with no dimension C pt , D pt , H pt Model parameters with dimension of free energy density C cf , D cf , H cf Model parameters with dimension of entropy density κ ρ , κ φ , κ η Gradient constants of corresponding order parameters Φ ss , Φ at , Φ sf , Φ gb , Φ S , Φ L Interpolation functions with subscript 'ss' standing for the substance, 'at' for the atmosphere/pore, 'sf' for the surface, 'gb' for the grain boundary, 'S' for the solid and 'L' for the liquid/melt P l , v l , v l Laser power, scan velocity and speed, respectively. v l = |v l | R l Nominal laser radius Relative volumetric specific heat L Latent heat during melting/solidification G eff gb Grain boundary mobility D ss , D at , D sf , D gb Self-diffusivity coefficients with subscripts indicating the corresponding phases M Isotropic Cahn-Hilliard mobility L φ , L η Isotropic Allen-Cahn mobilities b Magnitude of the resultant body forces acting on the melt flow Re, Ste, Pe T , Pe ρ , Ac φ , Ac η Reynold number, Froude number, Stefan number, Péclet numbers for thermal and mass transfer, Allen-Cahn numbers for melting-solidification and grain growth, respectively V i , i , v i Volume, density and translation velocity of the i-th dispersed rigid nanoparticle in the melt respectively f M , f A , f E Driving force densities due to the melt flow, the VDW and the electrostatic interaction, respectively σ Cauchy stress tensor of the melt d C , d S Inter-particle center distance and surface distance, respectively U A , U E Inter-particle potential of VDW and the electrostatic interactions, respectively A p:v , A m:v , A p:m Hamaker coefficients with subscript 'p:v' standing for nanoparticle in the vacuum, 'm:v' for medium/melt in the vacuum, and 'p:m' for nanoparticle in the medium/melt Z Screened-Coulombic interaction coefficient ψ i Nanoparticle surface potential ε p , ε m , ε 0 Permittivities of the nanoparticle, medium and vacuum, respectively N p Total amount of the nanoparticle, n p , n m Refractive indices of the nanoparticle and the medium, respectively λ P Plasma Debye length n e , ω e Free electron density and electronic absorption frequency of the melt z l , c l , m l Valence electron number, atom fraction and the atom mass of the l-th element in the alloy, respectively k B ,h, e 0 , N A Boltzmann constant, reduced Planck constant, elementary charge, and Avogadro constant, respectively V(X I , R I ) Powder bed as a set, containing the center coordinates (X I ) and the radius (R I ) of I-th powder ∇, ∇ d Gradient operators. Subscript 'd' represents the gradient along the inter-particle direction