Uniform Droplet Spraying of Magnesium Alloys: Modeling of Apollonian Fractal Structures on Micrograph Sections

Abstract

A variety of advanced manufacturing processes have been developed based on the concept of rapid solidification processing (RSP), such as uniform droplet spraying (UDS) for the additive manufacturing of metals and alloys. This article introduces a morphological simulation of fractal dendric structures deposited by UDS of magnesium (Mg) alloys on two-dimensional (2D) planar sections. The fractal structure evolution is modeled as Apollonian packs of generalized ellipsoidal domains growing out of nuclei and dendrite arm fragments. The model employs descriptions of the dynamic thermal field based on superposed Green’s/Rosenthal functions with source images for initial/boundary effects, along with alloy phase diagrams and the classical solidification theory for nucleation and fragmentation rates. The initiation of grains is followed by their free and constrained growth by adjacent domains, represented via potential fields of level-set methods, for the effective mapping of the solidified topology and its metrics (grain size and fractal dimension of densely packed domains). The model is validated by comparing modeling results against micrographs of three UDS-deposited Mg–Zn–Y alloys. The further evolution of this real-time computational model and its application as a process observer for feedback control in 3D printing, as well as for off-line material design and optimization, is discussed.

1. Introduction

A variety of advanced manufacturing processes have been developed based on the concept of rapid solidification processing (RSP), where the fast liquid-to-solid phase transformation (melt quenching) is utilized for the fabrication of parts with a unique microstructure and properties, as compared to the counterparts fabricated by conventional casting methods. Particularly, RSP-based metal additive manufacturing (AM) techniques, such as uniform droplet spraying (UDS) [1,2,3,4], laser powder bed fusion [5,6,7], and direct energy deposition [8,9,10], have attracted significant attention in recent years due to their exceptional design flexibility and rapid prototyping capability. The microstructure and properties of AM-fabricated metallic parts are closely connected to the equilibrium thermodynamics of solidification and/or non-equilibrium kinetics in RSP [11,12]. Therefore, the understanding of the process–microstructure–property relationship in these RSP-based AM processes is essential for the purpose of off-line material design, real-time process control, and property optimization.

In the rapid solidification of metals and alloys, the chaotic RSP evolution is often attributed to its nonlinear dynamics, in which random spatio-temporal variations in its initial and/or boundary conditions grow unstably to dominate the microstructural responses [13,14]. Chaos phenomena originate from the non-unique, probabilistic nature of impulsive, local nucleation by heterogeneous (or, rarely, homogeneous) mechanisms, as well as arm fragmentation of dendritic structures, out of which subsequent grain growth ensues. These effects are further accentuated by nonlinearities of the solidified fraction and composition with respect to the temperature drop in liquid–solid transformation areas of the alloy phase diagram, along with concentration and temperature dependencies of bulk conduction and diffusion, as well as solidification interface properties (such as the Marangoni effect). Away from atomic/molecular scales, descriptions of such chaotic mechanisms lack a characteristic length; i.e., they function similarly at dimensional scales in the millimeter–micron–submicron range, producing multi-scale self-similar (or fractal) structures upon solidification. As a measure of the scalar range within which such a morphological iteration is observed, the concept of fractal dimension (Hausdorff) [15] was introduced, bridging the gap between integer dimensionalities (0D/1D/2D/3D) to describe the wiggled interface complexity of space-filling forms.

The evolution of such fractal dendrites in RSP leads to the formation of a Brownian branch network. Fundamentally, the origin of this branching in dendritic crystallites is traced to Mullins–Sekerka instabilities [16] of the solid–liquid interface, akin to more familiar Hele–Shaw boundary patterns between two different viscous media [17]. On a topologically smooth boundary, the random appearance of a local protuberance of the equivalent radius $\varrho$, because of the surface tension γ effect, is accompanied by a pressure difference p = γ/$\varrho$ between the two domains. Upon volumetric growth dV, this protrusion yields a respective mechanical work dG = −p·dV, which reduces the Gibbs free energy G at the interface and therefore causes its spontaneous growth, resulting in the generation of a fractal Brownian network. To better understand such phenomena during RSP, dendric growth formulations or models were proposed in literature for specific RSP-based manufacturing processes. For instance, a free dendritic growth formulation was developed to describe the rapid solidification during UDS [18]. This is based on liquidus and solidus curve descriptions in a binary alloy phase diagram, incorporating droplet enthalpic losses and latent heat of fusion effects, along with solidification phase changes. This sharp-interface model describes the tip radius of a paraboloid branch morphology for a dendrite solid–liquid boundary, its unconstrained growth velocity kinetics, and the thermal and solutal profiles across its surface.

Moreover, the chaotic fracture and further growth of previously developed dendrite arms during RSP contribute another generative mechanism for fractal grain structures, in the form of Apollonian globular assemblies. Such fragmentation of dendrite arms is initiated through their local remelting by the latent heat release during recalescence, in combination with the local depression of the melting point due to the supersaturation of the crystal. The arm break-up is accomplished by liquid flow field forces, i.e., capillary action at the solid–liquid interface, resulting in Rayleigh-type instabilities of a liquid column, separating the arm fragment from the dendrite trunk. In literature, a fragmentation criterion has been established [19], whereby the requisite break-up time t_b predicted by the model should fall within a plateau duration t_p for such a fracture, i.e., t_b < t_p. Dendrite arm fragments, together with grain nuclei appearing on heterogeneous boundaries, and nucleant and container surfaces (or seldom in the bulk of the homogeneous melt), as described by the classical nucleation theory, are further grown freely into globular agglomerates. However, such growth is topologically confined by the expansion of adjacent grains, yielding densely packed domain morphologies. A framework for such constrained grain growth has recently been developed, with the eventual confined grain size connected with differential attributes of the dynamic temperature field during RSP [1,4]. Thus, arm fragmentation and nucleation-initiated constrained growth complements Brownian networks of branching dendrites with Apollonian packs of globular grains in a space-filing, hybrid fractal material structure [20,21].

However, the full 3D experimental imaging of such complicated fractal interfaces, necessary to corroborate theoretical solidification volume models, is technologically challenging. Instead, 2D planar sections of the fractal assembly displaying surface domains of the solidified grains are readily available as ordinary micrographs in the laboratory. To better understand the fractal dendric growth during RSP, in this paper, a 2D model is developed to simulate the Apollonian fractal structures on a specified flat section of RSP-fabricated magnesium (Mg) alloys. This real-time simulation employs level-set methods suited to planar formulations, and uses the thermomechanical process conditions to computationally predict solidified domain morphologies, along with their average grain size and fractal dimension. UDS deposition experiments of ternary alloys of the Mg–Zn–Y system were conducted, and the proposed 2D model is validated by experimental data (grain sizes and fractal dimensions). The potential applications of the developed model in the predictive design and optimization of Mg–Zn–Y alloy properties, and the in-process feedback for the UDS process control were discussed.

2. UDS Experiments

The RSP of metals and alloys is exemplified in this work by UDS of Mg–Zn–Y alloys. The experimental results are applied to validate and calibrate the theoretical model proposed in Section 3. In UDS, from an orifice of section area $w$ at the bottom of a heated crucible, the molten alloy is ejected under gas pressure at flow rate $\dot{V}$ as a jet, broken up via piezoelectric vibration at frequency F into a linear stream of uniform-sized spherical droplets of volume $V=\dot{V}/F$ (Figure 1a). Travel of the droplets in the UDS chamber at speed $U=\dot{V}/w$ cools them down to a desired thermal state at uniform superheated temperature T_o, which is controlled via the flight distance. At this initial state, droplets are quenched in an oil bath into mono-sized solidified spheres. Alternatively, the droplets impinge on a flat cooling substrate translated at speed v in a scanning motion, and coalesce into a deposit consisting of adjacent linear beads (Figure 1b). Single-bead deposits display an approximately half-elliptical cross section with half-axes A and B, aspect ratio ς = A/B, and area $W=\dot{V}/v$ = ½πAB. The solidified 3D structure of these materials consists of densely packed dendrite crystals intertwined with globular grains in a hybrid fractal topology.

The UDS experimental setting consists of an insulated stainless-steel crucible, a water-cooled inductive heater system, a piezoelectric vibration transducer controlled via a digital function generator, and a thermostatic controller with a K-type thermocouple. The system features an inert gas-pressurized cylindrical chamber, a stroboscopic camera monitor, an infrared pyrometer for the droplet stream, and a servo-driven deposition table. The latter carries a planar stainless-steel cooling substrate plate, elevated manually in the Z-direction by a scissor jack to control the droplet flight distance. This deposition system is moved by a motorized X–Y table driven by two micro-stepping motors through a motion control board in a laboratory PC workstation. The control software allows for the flexible design of the scanning pattern and speed of deposition, and includes a graphical simulator of the deposit geometry [20,21].

In the experimental tests, as-cast bars of three Mg alloys (Mg₉₇ZnY₂, Mg₈₈Zn₁₀Y₂, and Mg₇₆.₅Zn₂₀Y₃.₅) are diced under water cooling, ground, ultrasonically cleaned, and fed into the UDS crucible. The UDS setup is then sealed and repeatedly filled with ultra-high purity Ar gas, purging oxygen to pyrophoricity-safe levels. The thermostatic heating, piezoelectric transducer, stroboscopic monitor, and droplet deposition systems are properly initialized. Melt ejection is started via inert gas pressure in the crucible, and vibration maintains a steady droplet break-up and flight during the tests (Figure 1a) until the melt is depleted. The in-flight thermal state of the droplet stream is assessed by the IR pyrometry camera. The UDS process conditions in three tests are shown in

Table 1. Experimental UDS process conditions in laboratory tests [20].

UDS Process Condition	Test 1	Test 2	Test 3
Orifice Diameter [μm]	348	500	348
Pressure Difference [psi]	3	3	2
Piezo Frequency F [kHz]	2.39	1.03	2.39
Melt Temperature [K]	1023	973	973
Droplet Diameter D [μm]	700	1000	700
Droplet Flight Distance [m]	0.30	0.25	0.25
Droplet Landing Temperature To [K]	1012	967	965/917/844

, while the material properties of the three Mg alloys are summarized in

Table 2. Physical properties of Mg–Zn–Y alloys [22].

Alloy Properties	Mg97ZnY2	Mg88Zn10Y2	Mg76.5Zn20Y3.5
Liquidus temperature TL [K]	903	861	813
Solidus temperature TS [K]	713	697	680
Latent heat of fusion H [MJ/m3]	637	603	567
Liquid density $\varrho$L [kg/m3]	1680	2120	2570
Solid density $\varrho$S [kg/m3]	1850	2330	2820
Heat capacity c [kJ/kg K]	1.026	1.058	1.089
Thermal conductivity k [W/mK]	64	59	53
Thermal diffusivity α [10−6 m2/s]	35.4	25.0	18.1

[22]. After each test, the solidified deposits (Figure 1b) are separated from the substrate, sectioned transversely across the bead deposit direction under water cooling. The specimens are ultrasonically cleaned and encapsulated, polished, and etched, and optical micrographs are obtained by a metallurgical digital microscope. A material analysis of the sections was conducted by standard image processing software, through grayscale blob segmentation for edge detection by tessellated spatial occupancy methods [23].

3. Computational Modeling Framework

3.1. UDS-Deposited 2D Structure to Be Modeled

Figure 2 displays a small, rectangular, fixed section of the deposit on which the 2D structure is to be modeled and compared with experimental data. Upon deposition, this section is initially in a superheated liquid state at temperature T_o (right of Figure 2). With the cooling deflation of the temperature field T in the material, this section gradually solidifies as it is intercepted by the solidification front zone between moving isotherms $T_L^{*}$ and $T_S^{*}$ (middle of Figure 2), within which the solidification phenomena take place. These include grain nucleation, dendrite formation and arm fragmentation, free and constrained growth, and finally the solidification of inter-granular eutectic phases in residual interstices among solidified grains (left of Figure 2). This evolving structure is eventually frozen over the entire section as the solidification front moves past it. The morphology of this eventual structure is the final output of the model. This topological model does not explicitly determine the elemental concentration distribution in the section, along with its composition variations, e.g., by diffusion as it cools down to ambient temperature T_a.

In Figure 2, the 3D grain components growing freely in the solidification zone, such as dendrite trunks and branches, develop out of the initial nuclei or fragmented arms by a growth velocity which may be a directional but smooth polar function, and are thus assumed to be topologically quadratic (i.e., paraboloid, spheroidal, ellipsoidal, cylindrical, conical, toroidal, etc.). Upon the intersection with the modeled plane of Figure 2, all these components therefore initially yield 2D planar sections of a generalized elliptical form as individual domains. As their development becomes limited by adjacent domains, their constrained growth defines boundary segments of adjoining domains which are also shaped as elliptical arcs. Therefore, eventual domains on the solidified section exhibit generalized ellipsoidal shapes, with their boundaries delimited by and consisting of random elliptical arc segments, in self-similar patterns over a range of dimensional scales. These are tiled in space-filling arrangements, with the space among large domains filled by smaller ones and so on, i.e., in a generalized Apollonian globular pack. This topological argument is further elaborated in the constrained growth section below.

3.2. Thermal Modeling

Since the solidification phenomena are driven by temperature field dynamics, and because of the essentials-only nature of the model, simplifying assumptions are made for thermal modeling. Specifically, its initial composition and temperature state are assumed to be known and determined by the process inputs, including UDS droplet volume V, frequency F, and superheated temperature T_o. The boundary conditions, including convection in a cooling fluid, radiation to ambient space, and conduction to a solid substrate, are also assumed to be known. If quasi-stationary solidification prevails inside the material, with an insignificant internal melt flow field and convective transport to the solid phase, the heat transfer is dominated by Laplace linear conduction [1,4]:

$$\rho c\dot{T}=-k{\nabla }^{2}T+\dot{H},\mathrm{with}\overrightarrow{q}=-k\nabla T=\overrightarrow{n}\left[h_a\left(T-T_a\right)+\epsilon \sigma \left(T^4-T_a^4\right)\right]$$

(1)

where T is the local temperature. ρ, c, k, and α = k/ρc are the density, heat capacity, thermal conductivity, and thermal diffusivity of the material, respectively. For simplicity, these properties are assumed to be independent of the temperature, elemental concentration, as well as solid fraction in the solidification zone, in which the properties are assigned the average value of the solid phase properties at T_S and the liquid phase properties at T_L. In the solidification zone, H is the volumetric latent heat of fusion released during recalescence. At the surface of the material and in its outward normal direction n, the heat loss flux q may consist of convection to a fluid with heat transfer coefficient h_a, and radiation at emissivity ε and Stefan constant σ to an ambient temperature T_a.

Numerical solution of Equation (1) under the boundary conditions and initial temperature distribution T(x,y,z;0) at time t = 0 is typically achieved via off-line numerical techniques. For example, level-set methods [24] have been implemented to determine the coupled thermo-fluid fields in UDS droplets impinging and solidifying on a flat cooling substrate. For such simple UDS droplet conditions and for real-time computational efficiency, semi-analytical techniques have also been successfully employed [1], based on the linear time-varying superposition of Green’s functions. Similarly, in the deposition of UDS droplets over a translating flat substrate into a bead deposit of the semi-elliptical cross section, the temperature field was determined semi-analytically by the superposition of traveling Green’s (or Rosenthal) functions G [4]:

$$T\left(x,y,z;t\right)=T\left(x,y,z;0\right)+{\int }_0^t\eta \left(\tau \right)\dot{Q}\left(\tau \right)\mathrm{G}\left(r\left(\tau \right);t-\tau \right)d\tau \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}Q=\left[{\rho }_L{c}_L\left(T_o-T_L\right)+H\right]V$$

(2)

$$\mathrm{G}\left(r\left(\tau \right);t-\tau \right)=\frac{1}{2\pi kr}\exp \left[-\frac{v\left(r+X\right)}{2\alpha }\right]\mathrm{with}r\left(\tau \right)=\sqrt{{\left(x-X\left(\tau \right)\right)}^2+{\left(y-Y\right)}^2+{\left(z-Z\right)}^2}$$

(3)

The Rosenthal function G (Equation (3)) reflects the steady-state, Eulerian (moving with the source) thermal field generated by a unit power (1 W) point source at a location (X,Y,Z), translating at speed v in the X direction (X(τ) = X(0) + vτ), i.e., the temperatures at radial locations r from the source in a semi-infinite conductive medium (Figure 3a). In UDS, the source energy Q (Equation (2)) conveys the enthalpy of droplets of volume V due to their superheating at initial temperature T_o over the liquidus $T_L$, and latent heat of fusion H released. In Equation (2), the boundary conditions are enforced by the method of images, i.e., by adding fictitious mirror images of the heat source Q (Figure 3b), with respect to the conductive substrate interface (image Q′) and the convective and radiative deposit surface (Q″). The location and power of all images is reflected in Equation (2) by the modulation efficiency field η (Figure 3b). Figure 4 illustrates a typical steady-state temperature field and isotherm contours T over a moving half-ellipse section (Y,Z) of the deposit at the source location (X = 0) in UDS bead deposition.

Such Green’s-function-based superposition thermal modeling exhibits advantageous computation efficiency, because of its suitability for independent semi-analytical calculations locally at the section area of Figure 2 only, without the need for global temperature simulation over the entire material. For example, the solidification front zone in Figure 2 sweeps the section at velocity u, which is determined from the respective moving isotherms of the dynamic thermal field above by level-set methods as:

$$\dot{T}=-\overrightarrow{u}\nabla T\underset{}{\Rightarrow }u=-\dot{T}/(\overrightarrow{\mathrm{n}}\nabla \mathrm{T})=-\dot{T}/(\frac{\partial \mathrm{T}}{\partial \mathrm{n}})$$

(4)

where u depends on local cooling rate $\dot{T}$ and temperature gradient in direction n normal to the front.

3.3. Modeling of Nucleation and Fragmentation

In the solidification zone (Figure 2), the delimiting isotherm lines $T_L^{*}$ and $T_S^{*}$ can be determined from temperature distribution T, as intersections of the respective isotherm surfaces with the section plane. Because of non-equilibrium thermal conditions during RSP, these temperatures are depressed relative to the liquidus T_L and solidus T_S points of the phase diagram (Figure 5), due to supercooling ΔT = $T_L^{*}$ − T_L, comprising curvature, thermal, solutal, and kinetic components [20]. The pseudo-binary phase diagram of Figure 5 for various alloys of the Mg–Zn–Y system corresponds to a section of the ternary diagram between the Mg-rich corner and an intermetallic compound such as the I, X, or W phase as elaborated in Section 4. Note that, although the concentration profile across growing domains of Figure 2 as temperature T drops in the phase diagram (dashed line) can be determined (double curve in Figure 5), the model does not pursue such a composition simulation, as its predictions cannot be verified via plain micrographs alone. Figure 5 also indicates, via the lever rule at temperature $T_S^{*}$, the solid fraction and compositions of the a-Mg and e-pseudoeutectic phase with compound I as m_a/m_e = l₂/l₁ and m_Mg/m_I = L₂/L₁. This a-e pseudo-eutectic blend fills interstitial spaces among grains upon final solidification.

In the liquid phase of the solidification zone between $T_L^{*}$ and $T_S^{*}$ in Figure 2, according to the classical solidification theory, grain nucleation takes place when a critical nucleus size of radius r* is reached through a critical molar activation Gibbs free energy G*. This maximizes the balance of bulk latent heat release H versus surface energy γ developing at the solid–liquid interface, and, along with the liquid diffusion molar activation energy D, it determines the volumetric rate of nucleation $\dot{N}$. This, in turn, can be expressed as a probability density function (PDF) f(T) for nucleation over the solidification zone liquid area ΔS (Figure 2) and transformation time $∆t=-\left(T_L^{*}-T_S^{*}\right)/\dot{T}$, expressed as:

$$\dot{N}\left(T\right)={\dot{N}}_o\exp \left(-\frac{G^{*}+D}{RT}\right)\underset{}{\Rightarrow }f\left(S,t\right)=\zeta \dot{N}\left(T\right)\mathrm{with}\zeta ={\left({\int }_0^{∆t}{\iint }_{∆S}^{}\dot{N}\left(T\left(S,t\right)\right)dSdt\right)}^{-1}$$

(5)

$$G^{*}=\left(\frac{16\pi {\gamma }_{}^3{T}^2}{3H^2}\right)\frac{1}{{∆T}^2}\left[\frac{\left(2+\cos \theta \right){\left(1-\cos \theta \right)}^2}{4}\right]\mathrm{and}{r}^{*}=\left(\frac{2\gamma T}{H}\right)\frac{1}{∆T}$$

(6)

where ${\dot{N}}_o$ is a nucleation rate constant depending on the nucleant potency and population or extent, and ζ is a probability normalization factor. R is the universal gas constant, θ the wetting angle of the heterogeneous-nucleated solid–liquid interface with the nucleant surface, while the bracketed term in Equation (6) is absent of homogeneous nucleation. Figure 6 illustrates the profile of nucleation PDF f over a transverse section ($T_S^{*}$ to $T_L^{*}$) of the solidification zone (Figure 2), while f is assumed to be uniform in the direction of the isotherms.

In RSP, the pronounced cooling rate $\dot{T}$ is associated with high supercooling ΔT, with $∆T\approx \beta \dot{T}$ [1], with β as a proportionality constant. This leads to dendritic crystal growth, accompanied with fragmentation of dendrite arms because of supersaturation and capillary forces. According to fragmentation theory [19], this occurs when the arm dendrite break-up time t_b falls within a material-dependent critical plateau time t_p:

$$t_b\approx \frac{3R_t^3}{\alpha d_o}\left(1+\frac{m_l{c}_o\left(1-k_o\right)}{H/c}\frac{\delta }{\alpha }\right)<t_p$$

(7)

where R_t is the dendrite trunk radius (Figure 2), d_o is the arm capillary length, c_o is the composition of the binary alloy (Figure 5), k_o is the equilibrium partition coefficient, m_l is the equilibrium slope of the liquidus curve on the phase diagram (Figure 5), and δ is the solute diffusivity in the liquid phase. As with nucleation, dendrite fragmentation is driven by a thermal activation barrier similar to G* (Equation (5)) for a melting phase change involving latent heat H and solid–liquid surface tension effects γ (Equation (6)), along with solute diffusion similar to D (Equation (5)), to cause arm break-up at a critical dendrite size similar to r* (Equation (6)). Thus, due to its thermo-kinetic similarity to nucleation, fragmentation is assumed to initiate new grain growth out of fragmented arms analogous to that of nucleants above, i.e., governed by a fragmentation rate similar to $\dot{N}$ (Equation (5)), with physical parameters equivalent to G* + D and ${\dot{N}}_o$.

Consequently, nucleation and fragmentation are jointly considered in the model for initiating new grains in the liquid phase of the solidification zone (Figure 2). The PDF f sequentially determines the random centroid location (X,Ψ) of each initiator (nucleus or fragmented arm), also oriented randomly at azimuthal angle ϕ (Figure 7a) following a uniform PDF in the range [0, π). With no loss of generality as per prior discussions, each initiator is assumed to bear an elliptical shape, adjustable from the spheroidal nuclei to cylindrical dendrite arms through the values of its half-axes (a*,b*) (Figure 7b). The latter are randomly selected in the range (r*,R*), with R* as the maximum solidification zone thickness $R^{*}=\left(T_L^{*}-T_S^{*}\right)/\left(\frac{\partial \mathrm{T}}{\partial \mathrm{n}}\right)$, again following a uniform PDF.

3.4. Modeling of Domain Growth

Initiators are subsequently freely grown through liquid-phase solute diffusion, into generalized ellipsoidal domains. The growth velocity normal to the solid–liquid interface of the domain nυ is assumed to follow Arrhenius’ law, dependent on temperature.

$$\upsilon ={\upsilon }_o\exp \left(-\frac{\mathrm{D}}{RT}\right)\mathrm{i}.\mathrm{e}.\overrightarrow{n\upsilon }=\frac{d\overrightarrow{r}}{dt}=\overrightarrow{n}\upsilon \underset{}{\Rightarrow }\overrightarrow{r}\left(t\right)=\overrightarrow{r}\left(0\right)+{\int }_0^t\overrightarrow{n}\left(\tau \right)\upsilon \left(\tau \right)d\tau$$

(8)

At the domain boundary (x,y), described relatively to the centroid (X,Ψ) as (χ,ψ) or (r,φ) in polar co-ordinates (Figure 7a), the growth velocity can be derived by first-order Taylor series expansion as:

$$\begin{array}{c}\upsilon \left(x,y\right)=\upsilon \left({\rm X},\Psi \right)+\mathsf{\Delta }\upsilon \left(\chi ,\psi \right)\\ \mathsf{\Delta }\upsilon =\frac{\partial \upsilon }{\partial T}\nabla T\overrightarrow{r}=\frac{{\upsilon }_{o}D}{R{T}^2}\exp \left(-\frac{\mathrm{D}}{RT}\right)\left[\frac{\partial T}{\partial x}\chi +\frac{\partial T}{\partial y}\psi \right]=\frac{{\upsilon }_{o}Dr}{R{T}^2}\exp \left(-\frac{\mathrm{D}}{RT}\right)\left[\frac{\partial T}{\partial x}\cos (\varphi +\phi )+\frac{\partial T}{\partial y}\sin (\varphi +\phi )\right]\end{array}$$

(9)

For the free growth of such a single domain, its expanding boundary r(φ,t) (Equation (8)) can be illustrated by the level-set method as the intersection of the micrograph plane z(x,y) = 0 with a potential function field P(x,y;t), rising vertically in the z direction with time at speed ω, e.g., ω = ξu, with ξ as a rate constant (Figure 7b). The evolution of the domain boundary can thus be derived from the topology and motion of the potential field as:

$$\dot{P}=-\overrightarrow{\upsilon }\nabla P=\omega \underset{}{\Rightarrow }\upsilon =-\omega /\left(\frac{\partial P}{\partial n}\right)\mathrm{and}\frac{\partial P}{\partial n}=-\frac{\omega }{\upsilon }\underset{}{\Rightarrow }d\overrightarrow{r}=-\frac{\overrightarrow{\upsilon }}{\omega }dP$$

(10)

with the last expression used to gradually construct the potential field from its successive isolines.

When several domains i expand simultaneously in the solidification region, the growth of their boundary is confined by that of adjacent domains. Figure 8 illustrates such constrained growth of i = 3 domains. When the boundaries of two domains 2 and 3 initially touch each other tangentially at location r₂₃ at time t₂₃, the trajectories of their intersection points for t > t₂₃ define the locus of the joint interface. When such bilateral interface trajectories among three domains 1, 2, and 3 meet at location r₁₂₃ at time t₁₂₃, then they are terminated at that fixed triple point (Figure 8b):

$${\overrightarrow{r}}_{23}\left(t\right)\equiv {\overrightarrow{r}}_2\left(t\right)={\overrightarrow{r}}_3\left(t\right)\mathrm{for}{t}_{23}<t<t_{123}\mathrm{such}\mathrm{that}{\overrightarrow{r}}_{123}\equiv {\overrightarrow{r}}_1={\overrightarrow{r}}_2={\overrightarrow{r}}_3\mathrm{at}{t}_{123}$$

(11)

Since the intersections of domain boundaries r_i on the micrograph plane at time t correspond to the intersections of their potential fields P_i at the respective level P = P(0) −ωt, the composite potential field P of all evolving domains can be used to identify their joint interfaces and triple points. This is illustrated in Figure 8a, where all domain boundaries can be conveniently and simultaneously identified as the projections of intersections of potential fields on the section plane (or just the top view of the composite potential field P):

$$\overrightarrow{r},t:P_1\left(\overrightarrow{r},t\right)=P_2\left(\overrightarrow{r},t\right)=\dots =P\left(\overrightarrow{r},t\right)\mathrm{with}P\left(\overrightarrow{r},t\right)\equiv \underset{i}{\max }\left[P_i\left(\overrightarrow{r},t\right)\right]$$

(12)

Therefore, the level-set method provides an elegant tool for the facile mapping of complicated domain boundaries on the solidified section, along with their formation times during solidification through their levels on the composite potential field.

As illustrated in Figure 8, a planar domain solidified under constrained growth conditions by other adjacent domains maintains a generalized elliptical shape out of its free growth development, but its boundary consists of elliptical contact arcs with its neighboring domains. As can be proven by geometric inversion, the curvature 1/R of such an arc as the locus of intersection points between two arcs of signed curvatures 1/R₁ and 1/R₂, respectively, expanding at the same rate (dictated by temperature T at the contact, Equation (8)), is the mean of the two adjoining curvatures:

$$\frac{2}{R}=\frac{1}{R_1}\pm \frac{1}{R_2}\mathrm{and}R=\frac{2R_1{R}_2}{R_1\pm R_2}$$

(13)

where the radius R is the harmonic mean of the signed R₁ and R₂. This establishes the self-similarity of domain boundaries in a generalized Apollonian pack of the fractal grain section (Figure 2). Finally, as this solidification section is progressively swept by the $T_S^{*}$ isotherm, its areas not occupied by already developed domains (i.e., with P = 0) gradually solidify to the interstitial phase at the (pseudo-)eutectic composition e (Figure 5).

4. Computational Implementation

The setting of the computational simulation follows the theoretical framework outlined in the previous section. Initially, the UDS process conditions (Table 1) are used to calculate the geometry of a bead deposit, with aspect ratio ς of its cross section calibrated to match the experimental wetting conditions (e.g., ς = 1, 1.28, and 0.43 in Tests 1, 2, and 3, respectively, for alloy Mg₉₇ZnY₂). Along with the material properties (Table 2), this is employed to determine the initial thermal state of the UDS droplets in agreement with the infrared measurements, upon impinging on the deposit, along with the boundary heat transfer from the latter to the substrate and chamber ambient (Equation (1)). The dynamic thermal field in the deposit is modeled by a superposition of Green’s/Rosenthal functions (Equations (2) and (3), and Figure 3a), with the mirror-image efficiency distribution calibrated via the droplet thermal state pyrometrically (e.g., modulated to average <η> = 0.38, 0.35, and 0.34 in Tests 1, 2, and 3 for alloy Mg₉₇ZnY₂; Figure 3b). The resulting temperature distribution on the bead deposit section (Figure 4) is analyzed for its local differential attributes (gradients and cooling rates) and trajectory of its solidification front isotherms (Equation (4)) in the region of the examined section area (Figure 2), and in relation to the phase diagrams of the alloys (e.g., Figure 5).

Based on this dynamic thermal analysis, the model next implements the simulation of domain initiation in the solidification zone of the section area, while this is swept by the respective front (Figure 2), via the probability density function f (Equation (5)). Its activation energy exponent G* + D for nucleation (Equation (6)) and dendrite fragmentation (upon checking the condition in Equation (7)) is calibrated by matching the model-predicted number of initiated domains over the section area with that measured on the respective experimental micrographs as in [1,4]. The initiation rate parameter ${\dot{N}}_o$ is calibrated through the normalization factor ζ (e.g., to ζ = 0.12, 0.09, and 0.11 in Tests 1, 2, and 3 for alloy Mg₉₇ZnY₂) so that the total solidification initiation probability during the sweeping of the current unsolidified region areas equals unity. This is effected by random sampling for each initiator centroid location (X,Ψ) following PDF f over the solidification zone (Figure 6), for its orientation ϕ uniformly in [0, 2π), and for its initial half-axes (a* and b*) size uniformly in the range (r*, R*) as in the previous section. When such a random selected initiator is found to overlap previously solidified regions, repeat sampling is used for its replacement. Next, free domain growth (Equations (8) and (9)) is simulated with velocity factor υ_o and diffusion exponent D calibrated by matching the model-predicted average size of eventual domains with that on laboratory micrographs as in [1,4]. The level-set potential field P is constructed for each domain (Equation (10)) and elevated at rate constant ξ = 1, so as to assess boundary intersection conditions (Equation (11)) upon constrained domain growth, and to map joint interfaces of adjacent domains (Equation (12)) over the top view of the composite potential field (Figure 7 and Figure 8).

The simulation model is implemented on commercial programming software (Matlab^®), and it is capable of real-time computation efficiency relative to the experimental UDS process on modern PC hardware. Simulated sections are processed off-line for the apparent size DS of crystal domains of arbitrary geometry, from their average section area S. As a measure of the dimensional range of self-similarity in the resulting complex boundary patterns of generalized Apollonian packed domains (Equation (13)), the model determines the fractal dimension FD of the solidified sections (Figure 9). This is performed off-line by the tessellation of the domain interface patterns by orthogonal quadtree rasters with a progressively finer resolution, i.e., with pixel size (s₂ < s₁), and using box counting for the increasing number of pixels intercepting the boundary pattern (N₂ > N₁):

$$DS\equiv 2\sqrt[]{\frac{S}{\pi }},FD\equiv 1-\underset{s_2/s_1\to 0}{\lim }\frac{ln\left(N_2/N_1\right)}{ln\left(s_2/s_1\right)}$$

(14)

The same method is also used to determine the average domain size and fractal dimension of the boundary patterns on experimental micrographs.

5. Results, Discussion, and Conclusions

The structure of simulated sections, along with the apparent domain size and fractal dimension, is validated against those of the respective optical micrographs across the bead deposits of the three Mg–Zn–Y alloys at the three UDS test conditions (Table 1 and Table 2). The resulting average values and variation ranges of the attributes domain size and fractal dimension are comprehensively summarized in

Table 3. Experimental and simulated average and range of domain size and fractal dimension.

Material/Test	Test 1	Test 2	Test 3
Mg97ZnY2	Exp DS = 22.7 ± 8.1 μm Sim DS = 21.9 ± 5.7 μm	Exp DS = 19.1 ± 5.0 μm Sim DS = 19.7 ± 3.8 μm	Exp DS = 20.5 ± 5.8 μm Sim DS = 20.8 ± 4.4 μm
	Exp FD = 2.40 ± 0.11 Sim FD = 2.35 ± 0.09	Exp FD = 2.31 ± 0.08 Sim FD = 2.28 ± 0.08	Exp FD = 2.36 ± 0.09 Sim FD = 2.37 + 0.10
Mg88Zn10Y2	Exp DS = 21.4 ± 7.4 μm Sim DS = 20.5 ± 5.0 μm	Exp DS = 18.3 ± 5.3 μm Sim DS = 18.9 ± 3.6 μm	Exp DS = 19.5 ± 6.2 μm Sim DS = 20.1 ± 4.0 μm
	Exp FD = 2.46 ± 0.13 Sim FD = 2.42 ± 0.11	Exp FD = 2.37 ± 0.09 Sim FD = 2.35 ± 0.08	Exp FD = 2.40 ± 0.12 Sim FD = 2.44 ± 0.11
Mg76.5Zn20Y3.5	Exp DS = 15.3 ± 9.5 μm Sim DS = 14.6 ± 6.3 μm	Exp DS = 12.7 ± 6.5 μm Sim DS = 13.0 ± 5.0 μm	Exp DS = 13.1 ± 7.4 μm Sim DS = 11.3 ± 6.1 μm
	Exp FD = 2.53 ± 0.16 Sim FD = 2.53 ± 0.13	Exp FD = 2.40 ± 0.13 Sim FD = 2.43 ± 0.12	Exp FD = 2.45 ± 0.14 Sim FD = 2.51 ± 0.13

, while Figure 10 and Figure 11 compare typical experimental microstructures at various magnifications with the respective computational model predictions. In particular, Figure 10 illustrates representative micrographs and modeled sections for Mg₉₇ZnY₂ alloy deposits under Test 1, 2, and 3 conditions, while Figure 11 shows these results for all Mg₉₇ZnY₂, Mg₈₈Zn₁₀Y₂, and Mg₇₆.₅Zn₂₀Y₃.₅ alloy deposits under Test 3 UDS conditions. In the simulated sections, bright areas show crystalline grain domains (a-Mg phase), while dark regions correspond to the intergranular, pseudo-eutectic (e-phase), eventually solidified interstices.

In general, the sections in Figure 10 and Figure 11 indicate deposits with full density. This is attributed to sufficiently superheated droplets at T_o, with the liquid e-phase filling gaps among previously deposited domains. Such pseudo-eutectics in intergranular areas of the Mg₉₇ZnY₂ alloy have been found to contain the I-phase (Mg₃Zn₆Y) and X-phase (Mg₁₂ZnY), while, for the Mg₈₈Zn₁₀Y₂ and Mg₇₆.₅Zn₂₀Y₃.₅ alloys, they have the I- and W-phases (Mg₃Zn₃Y₂) [20]. All sections display finer microstructures (DS = 10–60 μm) compared to the coarse-grained micrographs of the original cast ingots [20], attributed to the RSP conditions of UDS with higher cooling rates. Such RSP structures promise improved mechanical properties via grain boundary strengthening. In Table 3, it appears that highly superheated conditions of Test 1 increase the solidified domain size, while the lower superheating in Test 2 has the opposite effect, with respect to Test 3. Additionally, lower values of thermal conductivity for highly alloyed materials (Mg₈₈Zn₁₀Y₂ and Mg₇₆.₅Zn₂₀Y₃.₅) decrease cooling rates, favoring nucleation events, and drastically reducing the solidified domain size relative to Mg₉₇ZnY₂. In general, in Table 3, the model-predicted average domain size DS matches the laboratory measurements within ±14% over the various testing conditions of alloys. Further improvement can be obtained by more specific, individual calibration of the model parameters (η,ζ) for each combination of alloy–test condition in UDS.

The predominantly globular domain assemblies (generalized Apollonian packs) in Figure 10 and Figure 11 appear to map 2D planar sections of 3D crystalline dendrites (Brownian branch networks). Such Apollonian domain agglomerates are dominant in lower solute concentration alloys, i.e., Mg₉₇ZnY₂, and Mg₈₈Zn₁₀Y₂ (Figure 10 and Figure 11a,b), in which the growth of Mg-rich a-phase grains prevails. On the other hand, the higher-solute-content Mg₇₆.₅Zn₂₀Y₃.₅ alloys (Figure 11c) with significant e-phase interstitial portions blend domain globules with Brownian branches in hybrid fractal dendritic patterns. In all UDS-processed Mg alloys, domains are preferentially grown out of dendrite arm fragments rather than heterogeneous nuclei, because of the prevalent RSP fragmentation conditions due to pronounced thermal gradients. In Table 3, the small superheated UDS droplets in Test 1 increase the solidified fractal dimension, while the larger droplets in Test 2 have the opposite effect relative to Test 3. In addition, the lower thermal diffusivity of highly alloyed materials (Mg₈₈Zn₁₀Y₂ and Mg₇₆.₅Zn₂₀Y₃.₅) increase temperature gradients, favoring dendrite fragmentation, and markedly increasing the boundary fractal dimension with respect to Mg₉₇ZnY₂. The simulation-estimated fractal dimensions match the experimental values within ±3% in Table 3. Therefore, the morphological simulation model of RSP appears to adequately reflect the fundamental solidification distribution dynamics.

However, a lack of consideration to the compositional field evolution at its current development stage prevents the model from predicting certain interesting features of the experimental micrographs in Figure 10 and Figure 11. In particular, inter-pass regions in deposit sections between overlapping beads appear, especially at a higher superheating T_o of the droplets (e.g., Figure 10a), because of bead surface oxidation in the UDS chamber. Nevertheless, because of the impulsive periodic impact of the UDS droplets on the deposit, such apparently discontinuous or porous oxide films still permit epitaxial domain growth across inter-pass surfaces (e.g., Figure 10c and Figure 11a). This enhances local thermal conductivity and therefore cooling rates in these regions, reducing the local domain size of deposited beads (e.g., Figure 10a) overlaid on previous ones, acting as cooling substrates.

In general, finer domains and higher fractal dimensions are observed in areas of higher cooling rates, promoting the effective nucleation and fragmentation of dendrite arms. Moreover, such phenomena are intensified during solidification at lower droplet superheating T_o (Figure 10b,c), because of slower growth rates and more confined domain expansion. High cooling rates at low solidification temperatures also promote more cohesive boundaries among adjacent domains, i.e., with less extensive intergranular phase areas (e.g., Figure 10c and Figure 11a) than during slower solidification at higher temperatures. The latter conditions yield higher solute concentrations at the end of solidification, promoting (pseudo-)eutectic phase areas in intergranular regions. Internally. in the grain domains, occasional striations (Figure 10) appear. particularly near domain edges rather than at their center, indicating the precipitation of a second phase on habit planes. This is because of higher solidification temperatures towards the end of domain growth than at its initiation, presumably due to the enhanced effects of local latent heat released during recalescence. This condition favors dendritic arm fragmentation, and causes serrations of domain boundaries through growth at preferential crystallographic directions (e.g., Figure 10a,b). These could prevent intergranular slip and creep, particularly at elevated operation temperatures, improving the mechanical behavior of Mg alloys.

Based on the results reported in this paper, future research in progress includes a full parametric study of UDS-processed alloys of the Mg–Zn–Y system, intended towards a comprehensive database for the specific calibration of the model to improve the accuracy of its predicted features in comparison to laboratory measurements. Further, the off-line coding of thermal boundary and solidification rate conditions into look-up tables for on-line model reference is presently pursued, to enable real-time computation efficiency of the simulation including in-process computation of apparent domain size and fractal dimension. Moreover, a model of the nucleation rate explicitly increasing with subcooling until recalescence is currently being implemented. Finally, the incorporation of phase and elemental concentration data on the potential fields during gradual domain solidification is currently examined in the model, for comparison with local X-ray diffraction spectra of experimental sections and the study of material composition effects previously identified. The results of these studies will be documented under separate articles.

Considering the potential applications of the developed model, in its current implementation, the simulation offers a usable modeling tool for the direct process–structure relation of RSP binary alloys, along with its inverse structure–process design of UDS deposition via (trial-and-error) sensitivity analysis [25,26,27]. In combination with established material structure–properties models, the simulation can be useful for the direct predictive estimation of new alloy properties and applications by computational means, as well as inverse material design and optimization. The real-time computation performance of the model supports its in-process use as an RSP observer, e.g., for the feedback control of a UDS system to insure the desirable deposited alloy structure during additive manufacturing processes. Since real-time microstructure analysis of materials by non-destructive and non-invasive methods is technically impractical, the model can be run as an RSP observer in parallel to the actual UDS process, in order to provide real-time estimates of structural features, such as domain size and fractal dimension, as surrogate feedback to the closed-loop controller during 3D printing [28,29,30,31,32]. An adaptation scheme can also be designed to improve the simulation fidelity to the RSP process by the real-time identification and adjustment of its calibration parameters. Research towards adaptive feedback control of the UDS process in additive manufacturing is also currently underway.