A Computational Fluid Dynamics Investigation of Pneumatic Atomization, Aerosol Transport, and Deposition in Aerosol Jet Printing Process

Aerosol jet printing (AJP) is a droplet-based direct-write additive manufacturing process, which has emerged as a high-resolution method for the fabrication of a broad spectrum of electronics, e.g., interconnects [1–3], optical waveguides [4,5], optoelectronic devices [6–8], organic light emitting diodes [8], quantum dot devices [9,10], fine-pitch electronics [11], sensors [12,13], and transistors [14–16]. AJP has a large stand-off distance, allowing for fine material deposition and device fabrication on nonplaner surfaces (with feature sizes as small as 10 μm) [17–19]. It also accommodates a wide range of ink viscosity (approximately 1–2500 cP). This implies that novel and advanced ink materials, e.g., graphene oxide, PEDOT:PSS, carbon nanotube, and noble metal nanoparticles, can be (ultrasonically or pneumatically) nebulized and then additively deposited on a (flexible or rigid) surface toward device fabrication. The underlying physics behind ultrasonic atomization is based on high-frequency vibration, while in pneumatic atomization, a high-pressure flow of an inert gas results in droplet formation; these two mechanisms of aerosol generation will be delineated in Sec. 2. Figure 1 exemplifies flexible electronics, fabricated by the authors using the AJP process. Please note that while AJP has been reliably utilized for the fabrication of conformal electronics or printed 2-2.5D structures, complete three-dimensional (3D)-fabrication is still a looming challenge.

Despite the aforementioned unique advantages and the broad range of applications engendered by AJP, the process is highly unstable, complex, and prone to gradual drifts in machine behavior, leading to suboptimal (and sometimes unsatisfactory) material deposition and device fabrication [21,22]. Formation of overspray (and satellite particles), deposition of lines with discontinuous edges, excess generation of aerosols, and inconsistent flow collimation are examples of adverse phenomena, which cause a deviation from optimal regimes in AJP process. To a great extent, changes in material formulation and properties are the leading cause of the AJP process complexity. The fluctuation of material properties stems from, for example, ink temperature instability as well as solvent evaporation, leading to ink predrying and/or particle accumulation within the deposition head [22,23]. Consequently, implementation of physics-driven in situ process monitoring and control in the AJP process is inevitably a need. To realize this objective, the authors have already introduced an image-based platform for near real-time process monitoring as well as closed-loop control in AJP, as delineated in Refs. [21], [22], and [24–29]. The AJP system (machine) was integrated with temporal and image-based sensors, allowing for near real-time data acquisition from the process. In addition, a wide range of digital image/signal processing and control algorithms were developed for the characterization of two-dimensional (2D) and 3D line morphology traits. Although the established platform has significantly enabled AJP process monitoring and control, the underlying aerodynamics and physical phenomena behind AJP have remained unexplored. This gap is addressed in this study, which is a direct continuation of the authors' previous works [20–22], dedicated to development of a 3D compressible, turbulent multiphase flow computational fluid dynamics (CFD) model; the aim is to investigate the underlying physical phenomena behind aerosol generation, transport, and ultimately deposition in the AJP process.

There are several research works, which explore various aspects of AJP process monitoring as well as material deposition.

Wang et al. [30] developed a dual-material AJP process to fabricate tunable nanocomposites, utilizing both ultrasonic and pneumatic aerosol generation methods. The process was based on two independently controlled streams of materials, mixed prior to aerodynamic collimation (in the deposition head). The dual-material process allowed for deposition of nanocomposite structures as well as real-time control of the properties of the deposited structures. In a research work by Efimov et al., the AJP process was integrated with a spark discharge generator with the aim to deposit dry silver nanoparticles (15–100 nm in diameter) on heated silicon [31] and glass [32] substrates. As a result, this process does not require postdeposition sintering. While presenting novel and original AJP-based fabrication methods, these works did not explore the physics behind particle generation, collimation, and deposition on a surface.

Gu et al. [33] presented a new process monitoring method, based on inkwells, for on-demand deposition rate measurements. Chang [34] similarly presented a process control framework, using digital image processing techniques in addition to a vibration-based wavelet method to monitor the atomization process in AJP. The method was capable of not only estimating a stable fabrication time but also monitoring the morphology of AJ-printed structures. Using central composite design, they investigated the effects of AJP process parameters on printed line morphology. Overall, using the control framework, a 50% improvement in completion rate as well as a 20% reduction in process variation was obtained. However, no physical models were presented in their study in support of the experimental observations.

In a research work by Lombardi et al. [24], an image-based method was introduced for closed-loop control of the AJP process, including two types of controllers, i.e., proportional (P) as well as proportional-integral (PI). The aim was to control line width as a function of print speed (PS). The presented control method led to an enhancement of both linewidth and electrical properties of printed electronics. In a similar work, Thompson et al. [35] developed a velocity-regulated path planning algorithm to control the AJP process. Similarly, Smith et al. [36] presented a method for controlling the quality characteristics of AJ-fabricated flexible electronics. Sun et al. [37] developed a statistical method (processing microscopic images) for quality modeling of AJ-printed electronics. Using ensemble learning, a machine learning approach, Li et al. [38] introduced a method for AJP process monitoring and control. Verheecke et al. [39] as well as Mahajan et al. [40] demonstrated offline characterization and optimization of the AJP process for the fabrication of high-resolution silver interconnects. No physical models were discussed in all these works. Examination of the fundamental mechanisms and causes of process drifts could aid in more efficient characterization, optimization, and ultimately control of the AJP process.

Chen et al. [41] studied the effect of droplet size on overspray, an important characteristic of line morphology. In this regard, they developed a CFD model to investigate the principles of overspray formation as a function of the droplet size distribution and ShGFR.

Feng [42] studied the deformation of a sessile drop under an impinging jet. It was demonstrated that above a critical value, when viscous forces become dominant (compared to surface tension forces), unsteady droplet deformation occurs. Feng [43] also studied the dynamics of microdroplet impingement on a surface in regard to postdeposition phenomena of spreading, receding-relaxation, and wetting equilibrium. The results showed that when the influence of viscosity is negligible, periodic oscillation and aperiodic creeping-to-capillary equilibrium are discernable. In addition, contact angles as well as droplet viscosity were identified as two critical factors, influencing droplet bouncing after impingement remarkably. Feng, furthermore, investigated [44] the nature of negative pressure distribution during pneumatic aerosol generation with respect to geometric parameters, using a CFD model.

Discrete phase in a microcapillary was studied in a research work by Schulz et al. [45]. They showed that as aerosol flow velocity increases, the Saffman force becomes significant, leading to radially inward movement of the particles within the capillary. Similarly, Akhatov et al. [46,47] demonstrated that the Saffman force, drag force, as well as particle inertia are the main factors in the modeling of aerosol flow through long microcapillaries. In another research work, Akhatov et al. [48] also presented solutions for the problem of gas flow in a slowly converging microcapillary.

Chen et al. [41] studied the effect of droplet size on overspray with the aid of a CFD model. The objective was to explain the fundamental principles governing overspray formation in AJP as a function of sheath and aerosol gas flow rates as well as droplet size distribution. In addition, the influence of particle impact velocity and aerosol flow misalignment (stemming from an internal clearance in the deposition head assembly) on the formation of overstay were investigated. They exhibited that the aerosol flow collimation is enforced by Saffman lift force, acting more significantly on larger droplets, implying that smaller droplets contribute to the formation of overspray. They emphasized the existence of a tradeoff between future size and overspray, and thus the need for optimization of the atomization process to obtain an optimum particle size distribution.

Overall, the aforementioned studies have explored various aspects of the underlying physics behind AJP. However, they are based on 2D-CFD models and/or consider a small portion of the process. In addition, a few works [44,49] have been dedicated to pneumatic nebulization in the AJP process. This work is to address these gaps by forwarding a holistic 3D-CFD model that includes not only the 3D geometry of the deposition head as well as the pneumatic atomizer completely, but also all fundamental physical models behind aerosol nebulization and deposition, such as turbulence, compressible flow, and multiphase flow.

The goal of this research work is to establish a physics-based computational platform for prediction of aerosol flow regimes (as a function of AJP process parameters) and ultimately, physics-driven control of AJP process. In pursuit of this goal, the objective of the work is to forward a multiphysics CFD model to investigate the underlying phenomena behind aerosol generation, transport, and deposition in AJP process. The aim is to explain the causal aerodynamic phenomena behind experimental observations, using a compressible, turbulent multiphase CFD model developed for the flow of aerosols: (i) during aerosolization (atomization), (ii) within the deposition head, and (iii) after deposition on a free surface. This work tests the hypothesis that the complex aerodynamic states of AJP process is captured, using a multiphysics CFD model.

The rest of the paper is organized as follows: In Sec. 2, the AJP process along with its aerodynamic components (including an ultrasonic/pneumatic atomizer, virtual impactor, and deposition head) is discussed in addition to a review of the sensor-instrumented experimental setup used in this work. The development of a 3D-CFD model is delineated in Sec. 3, consisting of geometry modeling, meshing, and setting up the governing equations. Experimental results, model validation, as well as simulation results of aerosol generation, transport, and deposition are discussed in Sec. 4. Conclusions and future work are presented in Sec. 5.

The AJP process is centered on atomization/nebulization of an ink, which can be implemented pneumatically or ultrasonically. Figure 2 illustrates both the modes of atomization. In the pneumatic configuration, a high-pressure flow of an inert gas is injected into the ink reservoir via a head (as will be demonstrated later in Fig. 4 in Sec. 2.2) where due to the Venturi effect, the ink is drawn upward through the capillary and then sheared by the jet; as a result, a nonuniform mist/cloud of aerosols is generated. Therefore, an aerodynamic separator, called the virtual impactor (VI), is utilized to filter out droplets with low linear momentum using an exhaust flow and deliver a uniform stream of aerosols for deposition. In the ultrasonic mode of atomization (UA), a transducer is used for direct aerosol/droplet formation. The temperature of the UA process is kept constant using a chilled water system to ensure uniform ink viscosity and mitigate unavoidable process drifts during the deposition of material. A bubbler filled with ink solvent(s) is used to deliver a saturated gas flow for atomization and thus to compensate for solvent evaporation. The ultrasonically generated aerosols are more uniform (and to some extent finer) than those generated pneumatically [23,50,51]; this precludes the need for a separator such as the VI in the UA mode. The mean diameters of pneumatically and ultrasonically generated aerosols are approximately < 1 μm and < 0.2 μm, respectively [51]. Relatively monodisperse aerosols are generated using the ultrasonic mode of atomization, although it has limited ability to nebulize viscous materials [50]. A pneumatic nebulizer is typically used when the liquid viscosity is in the range of 20–1000 cP, and an ultrasonic atomizer is used for low viscosity fluids (1–20 cP) [51]. The uniform stream of aerosols flows toward the deposition head where a surrounding flow, called the sheath gas flow (ShGF), is introduced to collimate/focus the central aerosol flow to a narrow beam, deposited on free surface. The deposited material eventually experiences the postdeposition phenomena of spreading, receding, relaxation, and coalescence [43].

All printing experiments, discussed in the work, were carried out using an Optomec AJ-300 aerosol jet printer (Albuquerque, NM). Shown in Fig. 3, the experimental setup was further equipped with a 5 MP charge-coupled device (CCD) camera (Edmund Optics, Grasshopper, GS3-U3-50S5C-C, Barrington, NJ), to acquire online images from the printed structures. The camera was mounted on a 2.5$X$-10$X$ variable magnification lens (Edmund Optics, VZM 1000i, Barrington, NJ), illuminated by an LED ring light (AmScope, Irvine, CA). This light source has a maximum brightness of 30,000–40,000 Lux and a color temperature of 6000 K.

Demonstrated in Fig. 4, the pneumatic atomizer is composed of the following components: (i) a head with an internal capillary and an orifice, generating a multiphase jet of aerosols (based on the Venturing effect); (ii) a long stem, guiding the atomization gas flow toward the atomization head; (iii) a jar/reservoir, containing ink; and (iv) a cap, which ensures airtight operation and allows for ink temperature measurements via insertion of a thermocouple.

Figure 5 exhibits the main complements of the virtual impactor, including: (i) an impactor, which is basically a converging nozzle; (ii) a collector, which separates aerosols based on their linear momentum as a function of the exhaust gas flow rate (we note that for the central aerosol flow, not diverted toward the exhaust, this component acts as a diffuser, resulting in a decrease in the flow velocity); (iii) a stem, guiding the aerosol flow toward the impactor; and (iv) a retaining nut and housing. The outlet stream of aerosols (more uniform in particle size) flows toward the deposition head via a polymer tube, known as the mist tube. Accumulation of aerosols within the tube (resulting in droplet formation and ultimately nozzle clog) is one of the main causes of anomalies in the AJP process. Consequently, purging/cleaning the mist tube regularly is a need to maintain continuous AJP operation.

The deposition head is the last aerodynamic element in the AJP process, aiding in collimated deposition of aerosols. As demonstrated in Fig. 6, it consists of: (i) an upper plenum chamber (UPC); (ii) a lower plenum chamber (LPC); (iii) a combination chamber (CC); (iv) a long converging nozzle/microcapillary; and (v) a shell together with a retaining cap. The UPC and LPC guide the sheath gas flow toward the CC where it concentrically surrounds the aerosol flow. The collimated aerosol flow then passes through the nozzle (made up of ceramic), gaining momentum for deposition on a free surface. The exit diameter of the nozzle is an influential factor on the width and final morphology of the deposited material.

In this section, we forward a 3D CFD model to understand the underlying dynamics behind aerosol transport and deposition in AJP, and simulate formation of line topology under various process conditions. Building a CFD model involves: (1) geometry modeling; (2) meshing; (3) setting up the governing equations; (4) specifying the boundary conditions as well as material properties; and (5) selecting numerical schemes to solve a system of coupled, linearized algebraic equations relating cell center values (having discretized the governing partial differential equations (PDEs) of the problem). We note that the discretization of the governing equations is based on finite volume method. The simulation results are verified with the physical conditions of the problem, and then validated with offline experimental measurements. The modeling and simulations presented in this study were carried out in the ansys-fluent environment (V. 16.2). The objective is to develop the numerical solution to a compressible and turbulent multiphase flow problem: (i) in the deposition head (as an internal flow problem); (ii) during deposition on a rigid and movable substrate (as an external flow problem). In general, there are two types of errors, namely, discretization and linearization. The former reduces with refining the mesh, while the latter decreases by iteration. In our simulations, a linearization tolerance of $10−6$ was set for all conservation imbalances. This ensures obtaining an approximate yet reliable solution to the mathematical model (consisting of the governing equations and boundary conditions).

The 3D geometries of the deposition head (together with its surrounding environment) as well as the pneumatic atomizer were modeled based on accurate dimensions obtained from X-ray computed tomography (CT) imaging [22] and a patented design [52] (as will be illustrated in Figs. 8 and 9 later in Sec. 3.2.4). The geometry of the deposition head is composed of the following components: (1) ShGF inlet; (2) upper plenum chamber; (3) lower plenum chamber; (4) combination chamber; (5) nozzle; and (6) a relatively wide cylindrical volume (embedded beneath the nozzle to freely observe the effects of aerosol deposition).

The geometry of the pneumatic atomizer consists of: (1) atomization gas flow (AGF) port; (2) VI port; (3) jar/reservoir and cap (atomization environment); (4) gas flow stem; and (5) atomization head. Carried out using a high-resolution (submicron) X-ray system (GE Phoenix Nanome|X, Niskayuna, NY), the X-ray imaging allowed for not only observation and accurate measurement of the dimensions of the internal structures (of the assemblies and their components), but also understanding the flow of material within the assemblies (which aided in defining the boundary conditions of the CFD model, e.g., the walls and internal surfaces).

The entire volume of the geometries is then meshed using a mixture of smooth and soft quadrilateral elements compiled and refined with respect to both curvature and proximity. To obtain an accurate solution near the walls, three layers of inflation were considered on all surfaces (with a geometric growth rate of 1.2). This is to ensure that equilibrium turbulent boundary layers are covered with a sufficient number of cells [53]. In general, while generating inflation layers, to some extent, aids in obtaining a more accurate solution near a surface, it may lead to poor aspect ratio and thus deterioration of mesh element quality. In fact, there is a tradeoff between the number of inflation layers and mesh aspect ratio.

The term “soft” refers to the flexible behavior of a mesh element in terms of size or the number of divisions where the element mesh size may be affected by geometry features, e.g., curvature and proximity. In contrast, the “hard” behavior of a mesh element strictly preserves the element's size or shape. The term “smooth” refers to the geometric quality of a mesh element. To obtain a smooth mesh, the location and arrangement of nodes are iteratively changed with respect to their neighboring nodes or elements [53].

In Eq. (5), $u$⁠, $v$⁠, and $w$ are the $x$-, $y$-, and $z$-components of the velocity vector, $v$⁠, respectively.

where $k$ is the fluid turbulence kinetic energy [$J/kg$]; $μt$ is the turbulent/eddy viscosity [$kg/(s.m)$]; $Prk$ and $Prε$ are the turbulent Prandtl numbers for the turbulence kinetic energy and dissipation rate, respectively (Prandtl number is a dimensionless number, representing the ratio of viscous diffusion to thermal diffusion); $ε$ is the turbulence dissipation rate [$J/(kg.s)$]; $εD$ is the rate of dissipation resulting from the fluctuating dilatation of compressible turbulence [$J/(kg.s)$]; $Sk$ and $Sε$ are the turbulence kinetic energy and dissipation rate source terms, respectively. Note that $PrE$ is the turbulent Prandtl number for fluid energy. For the realizable $k$-$ε$ turbulent model, $Prk$ = 0.82, $Prε$ = 1.0, and $PrE$ = 0.85 [53]. The rate of strain of fluid elements is defined as: $S=2SijSij$where $Sij=0.5(∇v+∇vT)$⁠. Furthermore, $C1ε=1.44$ and $C2ε=1.92$⁠, while $C3ε$ is evaluated on the basis of the local flow direction (relative to the gravitational vector) [53]. In Eqs. (6) and (7), starting from the left-hand side, the first term—also known as the unsteady term—represents the rate of change of turbulence kinetic energy (⁠$k$⁠) and dissipation rate (⁠$ε$⁠), respectively; similarly, the second and third terms reflect the convection and diffusion transport, correspondingly; the fourth term is the rate of destruction; and the other remaining terms indicate the rate of production as well as sources. Compared to other $k$–$ε$ models, the realizable model reflects the spectral energy transfer more properly; this can be implicitly realized by looking at the sixth term of Eq. (7), i.e., $ρ {max(0.43,SkεSkε+5)}Sε$⁠, where there is no dependency on the eddy viscosity (⁠$μt$⁠) and thus, on the generation of turbulence kinetic energy. Note that the destruction term in Eq. (7), i.e., $ρ C2εε2k+μρε$⁠, is not of singularity (because of $μρε$ in the denominator).

We note that the whole term inside the outermost parentheses (called $Fμ$ hereafter) is constant in the standard $k$–$ε$ viscous model, while here is a complex function of $Ωij¯$ as well as $ωk$—in addition to $k$ and $ε$ known as the turbulence fields—which are the mean rate of rotation tensor and angular velocity of the system of rotation (moving reference frame), respectively.

In Eqs. (9) and (10), $v$ is the mean fluid velocity at a wall-adjacent cell of interest [$m/s$]; $Cμ$ is the turbulent viscosity constant; $k$ is the fluid turbulence kinetic energy at the wall-adjacent cell of interest [$J/kg$]; $τw$ is the wall shear stress [$kg/(m.s2)$] defined based on Newton's law of viscosity as $τw=μ(∂u∂y)y=0$where $u$ and $y$ are the $x$-component of velocity and the distance from the wall [$m$], respectively. Similarly, $y*$ represents the nondimensionalized distance from the wall. $ρ$ is the fluid density [$kg/m3$], and $μ$ is the fluid dynamic viscosity [$mPa·s$]. The term on the left hand-side of Eq. (9), i.e., $vCμ0.25k0.5τw/ρ$⁠, represents nondimensionalized velocity as a lumped combination of the flow and turbulence parameters.

A coupled two-phase flow model (allowing for interphase exchanges of mass, momentum, and energy between the continuous and discrete phases) was established with the aim for: (i) steady/unsteady particle tracking in the system and (ii) simulation and characterization of particle deposition under various process conditions. All of our discrete-phase simulations and calculations of particle trajectories were performed, having obtained a fully converged continuous-phase flow field.

Illustrated in Fig. 8, the following boundary conditions were defined for the transport and deposition of aerosols in the AJP process, modeled as a compressible, turbulent multiphase flow problem.

1. Mass flow inlet was the type of boundary set for the sheath gas flow (ShGF), being in absolute frame of reference. Based upon experimental sensor data (directly obtained from the Optomec AJ printer's data logger), both mass flow rate and gauge pressure were set as a function of the ShGFR. In addition, two turbulence parameters, i.e., turbulent intensity and turbulent viscosity ratio, were set at 1.5% and 10, respectively. The thermal condition of the inlet ShGF imposed a temperature of 300K.

2. Mass flow inlet was similarly chosen for the AGF. The mass flow rate and gage pressure were set according to the experimental data. The same turbulence parameters and thermal condition as the first boundary conditions were used.

3. Stationery wall with no-slip condition using the standard roughness model (having a roughness constant of 0.5) was the type of boundary set for the main body of the deposition head, which also includes the combination chamber as well as the nozzle (all being at 300K). Furthermore, “reflect” condition was the type of boundary set for the discrete phase model (DPM), utilizing the generic erosion model with a friction coefficient of 0.2.

4. Pressure outlet was the boundary of choice for the atmospheric environment surrounding the jet of aerosols depositing on the substrate, with “escape” condition for the DPM. The turbulence component of the boundary was defined using backflow turbulent intensity and viscosity ratio, set at 1.5% and 10, respectively (similar to the inlet conditions).

5. Moving wall with no-slip condition (utilizing the standard roughness model) was chosen to reflect the substrate conditions. In addition, “wall-film” condition was the type of boundary set for the discrete phase where the particle–substrate interactions—including four regimes of sticking, rebounding, spreading, and splashing [53]—were estimated using Stanton-Rutland model.

Illustrated in Fig. 9, the following boundary types were defined in the CFD model, reflecting the governing conditions of pneumatic atomization in AJP:

1. Mass flow inlet was the type of boundary defined for the AGF. The first phase's flow rate (i.e., nitrogen) as well as the flow pressure was set based on experimental sensor data (from Optomec AJ's data logger). In addition, the second phase's flow rate (i.e., functional ink's) was set to zero at this boundary.

2. Pressure outlet was the type of boundary defined for the pneumatic atomizer's exit flow.

3. Stationery wall with no-slip condition was defined for: (i) the body of the atomizer, (ii) the stem, and (iii) the head (see Sec. 2.2).

4. A volume fraction of 1 (i.e., 100% liquid) was patched for the ink volume to enforce the presence of the ink in the atomizer prior to simulation.

Please note that in the case of the UA, the CFD model would be changed in a twofold manner: (i) change in the geometry (see Fig. 2) and (ii) change in the mathematical model (including the governing equations and the boundary conditions, delineated in Sec. 3.2). The underlying physics behind the UA process is based on high-frequency vibration and thus generation of capillary waves (where the amplitude is the driving mechanism for aerosol generation). The frequency of vibration in addition to surface tension and fluid viscosity are the governing variables in the UA process. Hence, not only would the fundamental transport equations—i.e., the Naiver–Stokes equations, the continuity equation, the energy equation, the equation of state, as well as the equations used for turbulence and discrete-phase modeling—be used, but also the physical equations capturing high-frequency vibration and capillary wave formation would be an integral part of the mathematical model.

Supplied in-house, both sheath and atomization gas flows are composed of a pure (4.8 grade) and dry stream of nitrogen, flowing at ambient temperature (21 °C). The electronic structures were deposited on a flexible substrate, made up of polyimide (UPILEX 125S, Ube Plastics, Tokyo, Japan). A silver nanoparticle ink (Paru Co., Seoul, South Korea), having a particle size of 80±10 nm, was used for printing four-point probe test structures (delineated in Sec. 4.2). The ink consists of the following components, mixed with a ratio of 5:1 by weight, respectively:

• MicroPEPG007MOP: This part is composed of silver nanoparticles (∼66 wt%) and 1-methoxy-2-propanol (MOP) as solvent, and has a bulk density and viscosity of 1.5–2.5 g/ml and 50 cP, respectively.

• MicroPEPG007EG: This part is composed of silver nanoparticles (∼33 wt%) and ethylene glycol (EG) as solvent, and has a bulk density and viscosity of 1.2–3.3 g/ml and 39 cP, respectively.

Detailed in Table 1, a single-factor factorial experiment was designed and conducted with the objective: (i) to investigate the influence of ShGFR on line morphology and (ii) to validate the CFD model. In this experiment, the ShGFR was randomly changed at four levels of 40, 60, 80, and 100 standard cubic centimeter per minute (sccm). Each treatment combination of the experimental design was repeated three times. Atomization gas flow rate (AGFR) and exhaust gas flow rate (EGFR) were set at 500  sccm and 450  sccm, respectively; as a result, an aerosol flow rate of 20  sccm was obtained. Utilizing the pneumatic atomizer and a 300 μm nozzle, silver nanoparticle-based devices (illustrated in Fig. 10) were printed on a polyimide substrate (nitrogen-dried, having rinsed with isopropyl alcohol and acetone) at a speed of 4 mm/s. Working distance/standoff, i.e., the distance between the nozzle and the platen, was set at 3 mm. Subsequent to printing, the devices were sintered in a convection oven (Binder, Inc., Bohemia, NY) at 200 °C for 1 h. A 5-min equilibration time was taken between each set point change to obtain a steady-state flow of material. Prior to printing, the ink was ultrasonically sonicated for approximately 2 min to ensure uniform dispersion of particles in the ink medium. Leading to an optimal process operability regime, the aforementioned process parameters were selected on the basis of a review of literature as well as the authors' prior experimental observations as delineated in Refs. [20–22].

From the fluid dynamics perspective, the nozzle size (diameter) directly affects: (i) the flow of aerosols in the deposition head, (ii) the width of the collimated jet of aerosols issuing from the nozzle exit, and ultimately, and (iii) the width of the deposited material (i.e., line width). Particle coalescence is the primary mechanism of line formation on a free surface, followed by spreading, receding, relaxation, and wetting equilibrium [43]. The impact of nozzle size on deposited line width was studied by Mahajan et al. [40]; it was observed that the line width declined with a decrease in the nozzle diameter (from 200 μm to 100 μm).

Using the high-resolution imaging system (delineated in Sec. 2.1), online images were acquired from the central trace/line of each printed device, as exemplified in Fig. 11. It is evident that as the ShGFR increases from 40  sccm to 100  sccm, the line width decreases. This implies that under steady aerosol generation and transfer, thicker lines are deposited as a result of the increase in the ShGFR. No significant amount of overspray was observed in the range of 40–100 sccm. This is an indication of stable and collimated material deposition.

Using the 3D-CFD model forwarded in this study, the collimated deposition of aerosols on a free surface was simulated at the four levels of the ShGFR, as demonstrated later in Sec. 4.5. The changes in the line width were contrasted against both online and offline experimental observations, plotted in Fig. 12. The online measurements were based on the images captured using the CCD camera (discussed in Sec. 2.1). Similarly, the offline measurements were based on high-resolution images acquired using an optical microscope (Carl Zeiss M1M Axio Imager, Oberkochen, Germany). The online, offline, and CFD characterization of line width was carried out based on an algorithm (developed in matlab and python), which would analyze an image or the spatial coordinates of a line and subsequently calculate the average line width on the basis of individual column-by-column (or pixel-by-pixel) line widths, having detected the line edges, as delineated in Ref. [22].

The trend estimated by the CFD model closely mimics that of the offline measurements. In other words, there is not statistically significant difference between the two trends (ρ = 0.93 where ρ is the Pearson correlation coefficient). We note that the trend captured by the online images is slightly different, stemming from the image-based measurements of the line width. The difference between the online and offline measurements (which are based on CCD camera and microscope images, respectively) arises from different image properties. The difference, to a great extent, stems from the less accurately defined line edges in the CCD images. In addition, the intrinsic time difference between the online and offline image acquisition leads to different physical properties due to, for example, solvent evaporation and drying, perceived as different line widths [22,60].

During the aerosol transport process, losses can occur due to, for example, aerosol settlement and/or diffusion, which affect the transport efficiency, the flow rate, and ultimately the scalability of the process. It was observed, based on the authors' experimental observations, that losses significantly occur in the transfer line, known as the mist tube (see Fig. 3), during the aerosol transport process. To minimize the adverse effects of aerosol settlement in the mist tube during the validation experiment and to ensure repeatability of the results, the mist tube was cleaned prior to the experiment. The combined material transport and deposition efficiency were monitored by in situ quantification of the morphology of the deposited material with the aid of image acquisition and processing (as illustrated in Figs. 11 and 12). The CFD model was established on the basis of the assumption that the material transport losses would be negligible for short material deposition simulations (typically less than 10 s) once steady-state conditions have been reached.

The CFD model is utilized to investigate the influence of the ShGFR on the flow of aerosols within the deposition head and on the collimated jet of aerosols issuing from the nozzle (discussed in Sec. 4.4) as well as after deposition of aerosols on a moving free surface (discussed in Sec. 4.5). In addition, the underlying phenomena behind pneumatic aerosol generation are explored in Sec. 4.6.

Using the CFD model, several simulations were carried out to understand the underlying phenomena behind experimental observations. They were focused on the flow of material: (i) within the print head and nozzle; (ii) during the deposition on a free surface; and (iii) after the deposition (where the deposited aerosols experience the phenomena of coalescence, spreading, receding, relaxation, and wetting equilibrium [43]). Figure 13 depicts the influence of the ShGFR on the aerosol flow streamlines as well as velocity profile within the combination chamber.

As the ShGFR increases, the flow of aerosols gets more collimated in the form of a narrow beam, which is further accelerated when passing through the nozzle. Based on the aerosol flow diameter at the inlet and outlet of the combination chamber, a dimensionless variable is defined, called the collimation ratio, representing the ratio of the inlet flow diameter to the outlet flow diameter. An increase in the collimation ratio is observed from approximately 4 (at the ShGFR of 40  sccm) to 5 (at the ShGFR of 100  sccm).

To assess the stability of the aerosol flow, Reynolds number (Re) was monitored at the combination chamber exit as well as at the nozzle exit, plotted in Fig. 14 as a function of the ShGFR. It turned out that the collimated flow of aerosols would remain laminar even at the ShGFR of 100  sccm. This implies that stable material transfer and deposition happen with respect to the set experimental conditions (discussed in Table 1).

A comparison of the jet of aerosols (exiting from the nozzle) at the four levels of the ShGFR is demonstrated in Fig. 15. The jet width decreases as the ShGFR increases. This is an indication of narrower and thicker deposition of the material on the substrate (as observed experimentally in Sec. 4.3). We note that the streamlines ultimately diverge near the surface due to the surface barrier and pressure accumulation. This phenomenon gets less prominent as the collimation power (exerted by the sheath flow) increases from the ShGFR of 40  sccm to the ShGFR of 100  sccm.

The impact of the aerosols on the substrate was investigated, utilizing the Weber number (⁠$We$⁠)—defined as the ratio of aerodynamic forces to surface tension forces—as illustrated in Fig. 16. We note that there exists a wall film due to the continuous deposition of the aerosols. Since the substrate is not heated, the cold-wall-with-film model (also known as cold-wetting [53]) can properly reflect the physics behind the aerosol impingement. Under the presence of the wall film, when the Weber number is low (i.e., ≤ 2, at the ShGFR of 40  sccm and 60  sccm), the impinging aerosols are reflected. On the other hand, at high Weber numbers (observed at the ShGFR of 80  sccm and 100  sccm), the aerosols stick to the substrate and contribute to the wall film formation.

The Taylor analogy breakup (TAB) model [53]—where the droplet surface tension forces, drag force, and viscosity forces are analogous to the restoring force, external force, and damping force of a spring mass system, respectively—was employed for postdeposition analysis of the aerosol oscillation and distortion. Visualized in Fig. 17, the normalized displacement represents the displacement of the aerosol equator from its undisturbed position (normalized with respect to the initial, undisturbed equator). Note that those aerosols, deposited along the center of the line, are distorted much more significantly than the edge aerosols, remaining almost undistorted upon impingement. Consequently, the cross-sectional profile of the line would be nonuniform, as implied from Fig. 11. The magnitude of the distortion increases with the ShGFR.

The line morphology is adversely affected as the PS increases from 1 to 4 mm/s as illustrated in Fig. 18. This stems from the spreading effect of the PS, which deprives the impinging aerosols of initial wall film formation, subsequent coalescence, and ultimately formation of a continuous line trace. In addition, tracking the Weber number of the aerosols—defined as the ratio of the inertia of a fluid to the surface tension—revealed that the aerosols contribute to the formation of the thin wall film (and thus the line trace) less significantly when the PS increases.

The objective of this section is to discuss the underlying physical phenomena behind pneumatic atomization (aerosol generation) in the AJP process.

As discussed in Sec. 4.2, the physics (including the boundary conditions) of the model was setup based on experimental data. The veracity of the model was corroborated by: (i) setting experimental boundary conditions and (ii) obtaining a converged solution at the end of each transient simulation iteration. Please note that obtaining real-time sensor data (which reflects the dynamics of the pneumatic atomizer), other than the presented data, is beyond the scope of this work. This requires sophisticated sensors, inserted within the atomizer, allowing for measuring the velocity and other properties of the atomization jet flow. Note that the sensor insertion should be implemented such that the main atomization stream remains undisturbed. Particle-size distribution can be captured utilizing image-based sensors.

Particle-size distribution is influenced by the jet velocity (fluid's inertia), ink viscosity, and ink surface tension. CFD simulation of the jet velocity allows for estimation of the particle-size distribution, which affects the morphology (and the functional properties) of AJ-deposited electronic devices. In addition, simulation of the velocity field of the atomization chamber—as a function of the orifice's and the capillary tube's diameter as well as of flow conditions, such as flow rate and pressure—would be critical from the atomizer design standpoint.

Having passed through the orifice, the AGF becomes completely turbulent, as implied from Fig. 19. This turbulent jet flow—having a Reynolds number (Re) of about 7000 at the orifice exit—breaks the ink, drawn upward based on the Venturi effect, into a nonuniform mist (cloud) of aerosols. The aerosols, subsequently, are carried toward the atomizer's outlet (and next, toward the virtual impactor) by the AGF. This flow experiences expansion twice: first at the orifice's exit (as a single phase flow, composed of N₂), and second at the atomization chamber's exit (as a two-phase flow of N₂ plus ink).

As exhibited in Fig. 20, the velocity field and the streamlines of both phases indicate the presence of not only a turbulent two-phase flow but also a vertex above the capillary tube (pertaining to the ink flow). This phenomenon emerges when the ink flow is sheared by the atomization jet, resulting in formation of aerosols. The streamlines of Phase II (i.e., ink) reflect randomness and thus turbulence in the atomization/nebulization process. In addition, they imply that a significant portion of the generated aerosols may be taken inward (back to the chamber) due to the occurrence of the vortex. It seems there is no ink flow above the atomization jet (acting as an aerodynamic barrier); however, due to the presence of turbulence, the aerosols are mixed with the AGF jet, and eventually carried toward the atomizer's outlet. Please note that the capillary tube's diameter (under constant atomization conditions) influences the momentum of ink delivery to the chamber (based on the Venturi effect), the resulting ink vortex, and ultimately the rate of aerosol generation.

Demonstrated in Fig. 21, the velocity vectors of Phase I reveal another phenomenon, i.e., the occurrence of circulation above the atomization jet. This is indicative of the fact that some of the fluid in the atomizer is taken inward (back to the atomization chamber) from both the upper and the lower edges. However, the lower intake flow is much less significant than the upper one due to the presence of the capillary tube and thus ink volume in that region. Please note that this phenomenon is in addition to the intake flow resulting from the vortex.

Furthermore, the orifice's diameter influences the dynamics of the issuing jet. For example, as implied from Fig. 21, increasing the diameter would result in less (narrower) circulation and less significant flow intake toward the atomization chamber. Also, under constant atomization conditions, this would lead to more significant turbulent flow (according to Reynolds number), affecting the particle-size distribution.

As shown in Fig. 22, the ink in the capillary tube is drawn upward due to the Venturi effect (as a result of reduction in the fluid pressure). There is a rapid transition in the flow regime from developing to fully developed. This can be attributed to the fact that the ink flow is a viscous flow (the ink viscosity is approximately 50 cP). In addition, this phenomenon is implicitly apparent from the graph illustrated in Fig. 22, where the ink flow velocity initially increases (indicative of the developing region) and then approximately levels off (indicative of the fully developed region; this region is mathematically described as $∂w∂z=0$⁠, where $w$ is the z-component of velocity vector, $V$⁠. Once the ink flow has reached the turbulent region (at the end of the capillary tube), the velocity starts rising again in a random fashion.

In Eq. (16), $ρ$ is the fluid density [$kg/m3$]; $v$ is the velocity vector [$m/s$]; $D$ is the characteristic length (i.e., diameter) [m]; and $μ$ is the fluid dynamic viscosity [$mPa·s$]. Please note that as the dynamic viscosity of a fluid increases, the length of the hydrodynamic entry region, $L$⁠, decreases.

Salary et al. [22] demonstrated that optimal material deposition regimes exist when the ShGFR is in the range of 40–100 sccm. The veracity of this observation was corroborated in this study where it was shown that when the ShGFR increases (from 40 to 100  sccm), more focused and stable material deposition occurs leading to dense deposition of aerosols with almost no formation of overspray (as implied from Figs. 11 and 16). Furthermore, pressure accumulation and nozzle clog at high ShGFRs (> 100 sccm) were reported in [22]. As illustrated in Fig. 23, while the flow pressure in the deposition head increases as a result of an increase in the ShGFR, the aerosol flow streamlines still remain intact, confirming the existence of an optimal process operability regime in the range of 40-100 sccm. In addition, the trend shown in Fig. 18 implies that electronic traces (lines) with optimal morphological characteristics are obtained when the PS is at 1 mm/s; this trend is consistent with the experimentally observed trend (discussed in Ref. [22]), which concluded that a PS in the range of 1–1.5 mm/s would lead to optimal line morphology.

In this work, a 3D compressible, turbulent, multiphase flow CFD model was forwarded to investigate the aerodynamics behind: (i) aerosol generation (aerosolization), (ii) aerosol transport (within the deposition head), and (iii) aerosol deposition (on a moving free surface) in the AJP process. A factorial experiment (consisting of AJ-deposition of silver nanoparticles on a polyimide substrate) was conducted to observe the influence of ShGFR on line morphology, and also to validate the accuracy of the CFD model. A good agreement between the CFD prediction of line morphology and the experimental measurements was observed. The following conclusions were deduced from the results of this study:

• As the ShGFR increases, the flow of aerosols gets more collimated in the form of a narrow beam, which is further accelerated when passing through the nozzle. As a result, narrower and thicker structures are deposited on a substrate.

• Under the presence of wall film, at low Weber numbers, the impinging aerosols are reflected, while at high Weber numbers, they stick to the substrate and contribute to the wall film formation.

• The edge aerosols remain almost undistorted upon impingement, while those deposited along the line center, are distorted much more remarkably.

• The atomization gas flow experiences expansion twice: (i) at the orifice's exit; and (ii) at the chamber's exit. The atomization gas flow becomes completely turbulent once it has passed through the orifice. This turbulent flow breaks the ink flow into a nonuniform mist of aerosols.

• It turned out that there would be a turbulent, two-phase flow as well as a vertex above the capillary tube in the atomization chamber. Due to the presence of the vortex, a significant portion of the generated aerosols is pulled back toward the atomization chamber. In addition, some of the fluid in the atomizer is, additionally, drawn toward the chamber due to the occurrence of circulation above the atomization jet.

• The ink flow in the capillary tube becomes fully developed swiftly. However, once the ink flow has reached the turbulent region, the velocity starts climbing randomly.

Nebulization in the AJP process results in formation of aerosols, typically having a diameter of 1–5 μm. In order to observe the distribution of aerosols, the flow domain should be meshed with micro-elements (cells), which make the CFD model computationally expensive. This challenge will be addressed as part of the authors' future work. Furthermore, the effects of compressibility (in terms of density and pressure changes) during atomization will be investigated. This CFD study will be expanded with incorporation of more critical AJP process parameters, such as nozzle size as well as ink properties (viscosity, solvent evaporation, etc.).

This material is based upon work supported, in part, by Air Force Research Laboratory (Funder ID: 100006602) under agreement number FA8650-15-2-5401. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Air Force Research Laboratory or the U.S. Government.

Prahalada Rao thanks the National Science Foundation (NSF) (Funder ID: 100000147) for funding his work under awards OIA-1929172, CMMI-1719388, CMMI-1920245, CMMI-1739696, and CMMI-1752069. Modeling, sensing, and analytics of additive manufacturing was the major aspect of CMMI-1752069 (Program Officer: Kevin Chou). Supplemental funding for CMMI-1752069 obtained through the NSF INTERN program (Program Officer: Prakash Balan) and CMMI Data Science Activities (Program Officer: Martha Dodson) is greatly appreciated.

In addition, the authors would like to sincerely acknowledge the Small Scale Systems Integration and Packaging (S³IP), the Center for Advanced Micro-electronics Manufacturing (CAMM), the Center for Flexible Hybrid Medical Device Manufacturing (FlexMed), and the Analytical and Diagnostics Lab (ADL) at State University of New York (SUNY) at Binghamton.

.

 
• $A$ =

  area

 
• $CDg$ =

  drag coefficient

 
• $Cv$ =

  heat capacity at constant volume

 
• $CVM$ =

  virtual mass factor

 
• $Cμ$ =

  turbulent viscosity constant

 
• $D$ =

  characteristic length

 
• $dD$ =

  particle diameter

 
• $FDg$ =

  drag force per unit particle mass

 
• $Fg$ =

  gravitational force per unit particle mass

 
• $FPG$ =

  pressure gradient force per unit particle mass

 
• $FSL$ =

  Saffman lift force per unit particle mass

 
• $FVM$ =

  virtual mass force per unit particle mass

 
• $g$ =

  gravitational acceleration vector

 
• $k$ =

  fluid turbulence kinetic energy

 
• $n$ =

  surface normal vector

 
• $p$ =

  fluid pressure

 
• $Prε$ =

  turbulence dissipation rate Prandtl number

 
• $Prk$ =

  turbulence kinetic energy Prandtl number

 
• $R$ =

  universal gas constant

 
• $Re$ =

  Reynolds number

 
• $S$ =

  rate of strain of fluid elements

 
• $SE$ =

  energy source term

 
• $SM$ =

  momentum source term

 
• $Sk$ =

  turbulence kinetic energy source term

 
• $Sε$ =

  turbulence dissipation rate source term

 
• $T$ =

  fluid temperature

 
• T =

  time

 
• $u$ =

  x-component of the fluid velocity vector

 
• $v$ =

  y-component of the fluid velocity vector

 
• $V$ =

  volume

 
• $v$ =

  velocity vector

 
• $v$ =

  mean velocity vector

 
• $vD$ =

  discrete phase (particle) velocity vector

 
• $w$ =

  z-component of the fluid velocity vector

 
• $We$ =

  Weber number

 
• $y*$ =

  nondimensionalized distance from a wall

Greek Symbols
 
• $ε$ =

  turbulence dissipation rate

 
• $εD$ =

  turbulence dilatation dissipation rate

 
• $κ$ =

  fluid thermal conductivity

 
• $μ$ =

  fluid dynamic viscosity

 
• $μt$ =

  turbulent/eddy viscosity

 
• $ρ$ =

  fluid density

 
• $ρD$ =

  discrete phase density

 
• $τw$ =

  wall shear stress

 
• $Φ$ =

  dissipation function

 
• $ωk$ =

  angular velocity of the system of rotation

 
• $Ωij¯$ =

  mean rate of rotation tensor

Aerosol Jet Printing Terms
 
• AGF =

  atomization gas flow

 
• AGFR =

  atomization gas flow rate

 
• CC =

  combination chamber

 
• CT =

  computed tomography

 
• DPM =

  discrete phase model

 
• EGF =

  exhaust gas flow

 
• EGFR =

  exhaust gas flow rate

 
• LPC =

  lower plenum chamber

 
• PA =

  pneumatic atomizer

 
• PS =

  print speed

 
• ShGF =

  sheath gas flow

 
• ShGFR =

  sheath gas flow rate

 
• UA =

  ultrasonic atomizer

 
• UPC =

  upper plenum chamber

 
• VI =

  virtual impactor