Universal scaling laws of keyhole stability and porosity in 3D printing of metals

Metal three-dimensional (3D) printing includes a vast number of operation and material parameters with complex dependencies, which significantly complicates process optimization, materials development, and real-time monitoring and control. We leverage ultrahigh-speed synchrotron X-ray imaging and high-fidelity multiphysics modeling to identify simple yet universal scaling laws for keyhole stability and porosity in metal 3D printing. The laws apply broadly and remain accurate for different materials, processing conditions, and printing machines. We define a dimensionless number, the Keyhole number, to predict aspect ratio of a keyhole and the morphological transition from stable at low Keyhole number to chaotic at high Keyhole number. Furthermore, we discover inherent correlation between keyhole stability and porosity formation in metal 3D printing. By reducing the dimensions of the formulation of these challenging problems, the compact scaling laws will aid process optimization and defect elimination during metal 3D printing, and potentially lead to a quantitative predictive framework.


Supplementary Method 1 Multiphysics modeling of keyhole-mode laser melting
A multiphysics model is developed to capture the complex physical phenomena in keyholemode laser melting, including multiple reflections of the laser beam, heat and mass transfer, Marangoni flow, vaporization-induced recoil pressure, and transport of metal vapor plume. In the computational implementation, the computational domain, including a gas domain and a substrate domain, is compared with the x-ray experimental configuration in Supplementary Fig.  3. The volume fractions of gas and substrate phases are defined as and , respectively, and they sum to unity (1) The volume fraction equations for gas and substrate phases can be written as where is the time, is the velocity, ̇ is the mass transfer rate due to vaporization, and denotes the density of the gas phase and denotes the density of the substrate phase, such that the total density at each point = + . The unit normal vector and curvature of the interface can be calculated from the volume faction: where is the gravitational acceleration, mush is the mushy zone constant, is a small number to prevent division by zero, and liq is the liquid fraction that can be calculated from the temperature (1). The interfacial forces include the surface tension force , the Marangoni force , and the recoil force due to vaporization : where is the surface tension coefficient, is the temperature, d d is the temperature coefficient of surface tension defined as the partial derivative of the surface tension coefficient with respect to temperature, recoil is the vaporization-induced recoil pressure, and atm is the atmospheric pressure.
To approximate the mass transfer rate ̇ and the recoil pressure recoil across the Knudsen layer, the Hertz-Langmuir relation (2) is typically used. Since the Hertz-Langmuir relation is only valid at high vaporization intensities (when temperature is much higher than boiling point or in vacuum), the ̇ is calculated by bridging the vaporization regimes with a smoothed thirdorder polynomial: where 1 , 1 , 1 , and 1 are fitting coefficients, is the retro-diffusion coefficient, is the molar mass of the vaporized species, is the gas constant, and the temperature thresholds and represent the low and high vaporization intensity regimes. The temperature-dependent saturated vapor pressure sat is calculated with the Clausius-Clapeyron law (3) as where is the latent heat of vaporization, and is the boiling point at atmospheric pressure.
To consider the effects of the atmospheric pressure, the recoil pressure can be expressed as where 2 , 2 , 2 , and 2 are fitting coefficients. Energy conservation equation can be written as where is the enthalpy of the material, is the thermal conductivity, ℎ ref is the reference enthalpy with respective to the reference temperature ref , is the heat capacity, and is the latent heat of melting. The volumetric source term represents the energy source applied on the interface between substrate and gas , which includes radiative source rad , laser source laser , and evaporative source evp : = rad + laser + evp (17) Those sources can be expressed as where SB is the Stefan-Boltzmann constant, is the material emissivity, ∞ is the ambient temperature. To account for the multiple reflection absorptions of laser beam, the absorbed energy flux laser is calculated using a ray tracing method (4,5): where represents the laser energy flux, represents the unit vector of the beam, represents the angle between the laser beam and normal of the keyhole interface , is the unit normal vector of the interface . The subscript 0 denotes the incident beam and denotes the th reflections. The Fresnel absorption coefficient Fr is applied, and 0 is a coefficient related to the types of lasers and materials, is the laser power, 0 is the laser spot radius, is the scan speed, is the radial coordinate, and is the vertical z-coordinate. More detailed calculation of the ray tracing method and its validation are provided elsewhere (6). An illustrative result of the ray tracing is shown in Supplementary Fig. 4. A conservation equation for metal vapor concentration in the gas phase is coupled with the momentum conservation equation: where represents the mass fraction of substrate species (one major component in the substrate is considered in this study) and is the mass diffusion coefficient. The gas density is defined using the ideal gas law for an incompressible flow: (25) where op is the operating pressure in the experimental chamber, is the molar mass of the gas (e.g., argon), and is the molar mass of the substrate (e.g., aluminum) species. External boundaries of computational domain are assumed to be adiabatic since they are sufficiently far from the heat source and the processing time is sufficiently short. A no slip condition for momentum equations and a zero diffusive flux for the species equation are set on the external boundaries except the top boundary, which is set as a fluid outlet with ambient pressure and zero mass fraction of the substrate species.
The governing equations are solved by the finite volume method using the non-iterative PISO scheme within ANSYS FLUENT 2020 R1 (7) using user-defined functions (UDFs). The Second Order Implicit Scheme is used for the transient formulation. The Least Squared Cell-Based scheme is used to compute gradients, PRESTO is used to compute pressures, and Second Order Upwind is used for momentum spatial discretization. An explicit Volume of Fluid (VOF) solver is applied with the CICSAM discretization scheme (8). The energy equation is discretized using the Power Law scheme, and the species equation is discretized using the Second Order Upwind scheme. An automatic local grid refinement technique is used. A uniform hexahedral mesh with an edge length of 8 μm is used initially, and then two levels of local refinement (mesh edge length down to 2 μm) is specified when the temperature associated with the region is higher than the solidus temperature of the material. A variable time step ranging from 1×10 -10 s to 1×10 -8 s is used such that the global Courant number is smaller than one.
Thirty simulations with different laser power and scan speed are conducted (ten for the Al6061 substate and twenty for the Ti-6Al-4V substrate). The thermophysical properties are provided in Supplementary Table 1, and the computational parameters are given in  Supplementary Table 2. Quantitative comparisons of keyhole aspect ratio and melt pool size between x-ray experiments and multiphysics simulations are presented in Supplementary Figs. 5 and 6.

Supplementary Method 2 Energy balance calculation and approximation
The energy balance in laser melting or additive manufacturing can be expressed as laser = reflect + convect + radiate + evaporiate + spatter + conduct (26) where laser is the total laser energy deposited during the process, reflect is the reflected energy, convect , is the convection energy losses, radiate is the radiation energy loss and spatter is the spattering energy loss. The portion of energy transferred within the substrate due to conduction is denoted by conduct .
A power balance can be obtained by deriving the energy balance equation with respect to time, such that a transient power balance or averaged power balance during a period of time can be analyzed: laser = reflect + convect + radiate + evaporate + spatter + conduct (27) where laser is equal to laser power . Similarly, the laser power consists of power losses due to reflection reflect , convection convect , radiation radiate , evaporation evaporate , spattering spatter , and transferred power due to conduction conduct .
Those powers can be approximated based on the multiphysics model: where the subscript denotes the gas phase, is the interface between substrate and gas, and ′ is a flat surface at the top of the substrate. A schematic of those two interfaces is shown in Supplementary Fig. 7. Other parameters have been described previously. These integrals are computed at each time step after the keyhole depression completes a rapid growth. Mean and standard deviation of the integrals during the quasi-steady state are recorded. The power loss due to spattering spatter cannot be accurately calculated from the current model. Alternatively, the average power loss due to spattering ̅ spatter can be approximated as where spatter is the volume of spattered droplets that can be approximated from the x-ray images by assuming that the separated droplets are spherical (an illustrative result is shown in Supplementary Fig. 8). The subscript denotes the liquid phase, and 0 is the observation time of x-ray imaging. Vaporization temperature is used in the approximation because most of the spattered droplets are ejected from the keyhole depression region where the temperature is near the . Experimentally, average power due to conduction can be measured by attached thermocouples (6). Those experimental data are used to validate the predicted data from the model.

Supplementary Method 3 Mathematical derivation of scaling parameters
To correlate the keyhole size and aspect ratio with process parameters and material properties in a compact form, scaling laws are derived from normalized governing equations and boundary conditions with appropriate assumptions and simplifications. We begin with a heat conduction problem in a semi-infinite region without considering the powder layer, Marangoni flow, evaporation and spattering. Previous analysis of the energy balance has shown that energy loss due to evaporation and spattering can be ignored as comparison with heat conduction. We consider a well-developed keyhole or transition regime where the laser energy is high enough to melt the powder around the laser, and thus the effect of particle morphology on the keyhole is neglectable. Ultrahigh speed x-ray observation (9) and qualitative arguments (10) support this assumption. The effect of the Marangoni flow will be analyzed later. We neglect the latent heat of melting because it is much smaller than the energy required to heat material to the melting point and consequently only affects the temperature distribution near the mushy zone. We also assume temperature-independent thermophysical properties at the melting point in the derivation, a previous study indicates that the inclusion of the temperature-dependent thermophysical properties does not change the keyhole depth qualitatively (11).
We consider a laser beam scanning a substrate with the scan speed along direction . Direction is the transverse coordinate, and is the normal into the substrate surface. In a reference frame that moves with the laser beam at a fixed relative location, the temperature field is governed by = − • (35) where is the density of the solid, is the heat capacity, is the temperature, is the velocity of the reference frame, is the speed of the frame which is equal to scan speed of the heat source, and is the unit vector in the x direction. We substitute Equation 35 into Equation 34 and consider a steady state: where is the thermal diffusivity that is defined as = .
The boundary condition considering a Gaussian heat source can be expressed as where is the effective laser absorptivity, is the laser power, 0 is the laser spot radius, 0 is the preheat temperature, and is defined as = √ 2 + 2 + 2 .
The goal of normalization is to acquire an equivalent but compact representation governed by a minimal set of dimensionless parameters. We define the dimensionless groups by dividing the temperature and spatial variables by their natural scaling factors, which are combinations of the dimensional parameters in this problem, as We are also interested in the size of a specific isotherm, for example liquidus temperature , so we define * = The maximum depth of the melting isotherm * is determined by the relation * = ( * , * , * , ) Thus, based on the above dimensional analysis the keyhole depth normalized by thermal diffusion length is a universal function of the normalized enthalpy enth * , normalized diffusion length * , and the absorptivity . We found a linear relationship between * and the ratio of enth * and * as * ∝ enth * * where Rayleigh number Ra is associated with buoyancy-driven flow, Marangoni number Ma characterizes the Marangoni flow, and Prandtl number Pr affects the morphology of the melt pool driven by Marangoni flow (12). The Weber number We, Bond number Bo and Capillary number Ca affect the shape of the melt pool free surface (12). Supplementary Table 3 lists the values of those dimensionless numbers for three substrate materials investigated in this study: Al6061, Ti6Al4V, and SS316. Those dimensionless numbers are roughly on the same order of magnitude, implying that fluid flow is similar for different materials under the same process conditions. Thus, ignoring the parameters associated with fluid flow in the melt pool ought not significantly affect the scaling for keyhole depth or aspect ratio. Based on a procedure similar to those discussed above, scaling laws for the front angle and inlet length of keyhole can also be obtained as shown in Supplementary Fig. 9 (Data S1). The keyhole front angle is approximately proportional to the keyhole aspect ratio, and thus it correlates well with the Keyhole number. We also conclude that the normalized keyhole inlet length * = 0 is dominated by the normalized diffusion length * .