Electrohydrodynamic jet (e-jet) printing is a microscale additive manufacturing technique used to print microscale constructs, including next-generation biological and optical sensors. Despite the many advantages to e-jet over competing microscale additive manufacturing techniques, there do not exist validated models of build material drop formation in e-jet, relegating process design and control to be heuristic and ad hoc. This work provides a model to map deposited drop volume to final spread topography and validates this model over the drop volume range of 0.68–13.4 pL. The model couples a spherical cap volume conservation law to a molecular kinetic relationship for contact line velocity and assumes an initial contact angle of 180 deg to predict the drop shape dynamics of dynamic contact angle and dynamic base radius. For validation, the spreading of e-jet-printed drops of a viscous adhesive is captured by high-speed microscopy. Our model is validated to have a relative error less than 3% in dynamic contact angle and 1% in dynamic base radius.

## Introduction

Electrohydrodynamic jet (e-jet, although some works use EHD-jet) printing is a microscale additive manufacturing technique with applications in fabricating new, low-cost, customized, printed electronics, and biological and optical sensors [1–3]. Similar to electrospinning in its actuation and inkjet printing in its application, e-jet printing yields drops of diameter 0.05*–*50 *μ*m, making it a high-resolution technique with orders of magnitude smaller drop volume than inkjet printing [4]. Demonstrated e-jet build materials include photocurable polymer precursors, solutions of polymers or biomolecules, and colloidal suspensions of conductors or semiconductors [1,5–13].

E-jet-printed topographies are formed by a sequential process of drop ejection, spreading, and coalescence. The printing input signal can be related to drop volume and synchronized with stage motion to eject drops of expected volume at given positions and times [14]. To e-jet print a desired microscale three-dimensional topography, the operator specifies a sequence of drops to print, each with its own volume, position, and ejection time. The resulting topography is formed by drop spreading and coalescence processes, which the operator cannot modulate midprint. Lacking models of drop spreading and coalescence, the current practice for achieving a desired printed topography is a heuristic tuning of the inputs. The tuning process can be expensive for a multimaterial, multilayer topography, with no guarantee of success. Furthermore, a change to the desired topography, materials, or environment may require retuning of the printing inputs. Hence, the development of drop spreading and coalescence models for e-jet printing will facilitate faster tuning and provide the foundation for alternative solutions through advanced control techniques.

Few authors have modeled drop spreading for e-jet printing, and they have not progressed beyond isolated drop systems. One work predicted final drop shape in response to a constant printing voltage [15]. Another work used an overall energy balance model with a hydrodynamic treatment of viscous loss to predict dynamic contact angle for known drop volume but lacked experimental validation using e-jet-printed drops [16].

The wetting dynamics literature contains detailed empirically derived and first-principles models of drop spreading for significantly larger volumes than e-jet-printed drops. These models leverage hydrodynamic, molecular kinetic, and combined approaches [17]. Remaining gaps in extending these models to the e-jet printing process include: (1) the development of new methodologies to identify experimentally driven model parameters and select appropriate initial conditions for the models; and (2) experimental validation of these models through well-designed e-jet printing experiments.

The specific contributions of this work to the field of e-jet-printed drop spreading include: (1) derivation of a model to predict dynamic contact angle from drop volume; (2) high-speed (20 kHz) measurement of picoliter drop spreading dynamics; (3) demonstration that a molecular kinetic relationship can fit contact line velocity to obtain generalized spreading parameters; and (4) the first validation of a model to predict dynamic contact angle.

This work is organized as follows: Section 2 describes the e-jet printing process. Section 3 develops a molecular kinetic theory model of drop spreading for e-jet printing. Section 4 describes an experimental method for capturing e-jet-printed drop spreading data and provides a model validation through the comparison of experimental and simulated drop spreading dynamics under various e-jet printing conditions. Finally, conclusions are made, and future directions are discussed in Sec. 5.

## Electrohydrodynamic Jet Printing Process

The actuation signal for e-jet printing is a large, typically 100–2000 V, voltage established on a metal-coated glass capillary nozzle suspended above an electrically grounded substrate by a small standoff distance, typically 10–200 *μ*m. The resulting electric field draws the meniscus into a *Taylor cone*, which ultimately breaks into a jet of smaller diameter than the nozzle orifice. Build material flows through this jet to the substrate. The jet ejection can be localized in time by applying to the nozzle a voltage pulse from low voltage *V _{l}* to high voltage

*V*for a pulse width

_{h}*T*[14] and further localized in space by halting stage motion during the voltage pulse, giving drop-on-demand printing. Figure 1 shows the setup, voltage pulse, jet, and drop spreading for e-jet printing with a nozzle of inner diameter 30

_{p}*μ*m, giving a picoliter-scale drop volume. Common nozzle inner diameters for high-resolution e-jet printing are 0.1

*–*5

*μ*m, which yield femtoliter drop volumes. However, imaging the spreading of femtoliter drops is optically limited.

## Dynamical Model of Drop Spreading

The modeled system is a nonvolatile liquid drop spreading on a flat, impermeable surface in the presence of an air atmosphere. The model's assumptions are:

*The drop is transported to the substrate as a discrete volume*$\Omega $

*, which remains constant during spreading. Volume conservation is expressed in the following equation:*

Assumption 2. *The drop shape is a spherical cap at all times, with*$\Omega $*related to dynamic base radius R(t) and dynamic contact angle*$\theta (t)$*by Eq.**(2)**[**18 **,**19 **]. R(t) and*$\theta (t)$*are shown in Fig.*2

Assumption 3. *The initial condition is a sphere in point-contact with the substrate, so*$\theta (0)=180\u2009deg$*, with t = 0 at the instant of the jet retraction from the substrate. This conservative assumption of the largest possible contact angle as the initial condition avoids the need for an accurate prediction of the initial contact angle.*

*R*(

*t*) as follows. Differentiating both sides of Eq. (2) with respect to time, applying Eq. (1), and isolating $(d\theta /dt)$ gives the following equation:

where *κ* is the equilibrium frequency of molecular displacements in one direction (either the wetting or the unwetting direction), *λ* is the average molecular displacement length, *θ _{e}* is the equilibrium contact angle, $\gamma lv$ is the surface tension at the liquid–vapor interface,

*k*is the Boltzmann constant, and

*T*is the absolute temperature [21]. The equilibrium contact angle is the contact angle attained by a drop at thermodynamic equilibrium [22].

*θ*. Model parameters are

_{e}*κ*,

*λ*,

*γ*, and

_{lv}*T*

The model prediction for *R*(*t*) is given by inserting the solution $\theta (t)$ into Eq. (4).

## Experimental Validation

To validate the model, high-speed images of drop spreading are captured and processed to measure $\theta (t)$, *R*(*t*), $\Omega $, and *θ _{e}*. One build material, nozzle, and substrate are used, with printing parameters $(Vh,Tp)$ varied to give a range of $\Omega $. A fitting of $\theta (t)$ and the derived $(dR/dt)$ give model parameters

*κ*and

*λ*. Next, simulations of $\theta *(t)$ and $R*(t)$ are calculated, with $*$ denoting simulated values. The error between the measured $\theta (t)$ and

*R*(

*t*) and their simulated values is calculated to quantify model accuracy.

### Printing Method.

The printer used in this work is a custom-built, high-resolution e-jet printer housed at the University of Michigan and shown in Fig. 3. This printer holds the nozzle fixed and moves the substrate with precision motion stages (PlanarDL-200XY, Aerotech, Pittsburgh, PA), with ± 0.1 *μ*m bidirectional position repeatability. To image with high-frame rate (20 kHz) and spatial resolution (0.65 *μ*m/pixel at the object), a high-speed camera (Phantom V9.0, Vision Research, Wayne, NJ) is used with 20× magnification (see Table 1). A high-power white light emitting diode spotlight (SL162, Advanced Illumination, Rochester, VT) and custom glass light guide are used for backlighting.

Component | Mag. | Manufacturer |
---|---|---|

EO M Plan HR $\u221e$ objective | 10× | Edmund Optics |

InfiniTube FM-200 tube lens | 1× | Infinity Photo-Optical |

DL doubler tube | 2× | Infinity Photo-Optical |

Component | Mag. | Manufacturer |
---|---|---|

EO M Plan HR $\u221e$ objective | 10× | Edmund Optics |

InfiniTube FM-200 tube lens | 1× | Infinity Photo-Optical |

DL doubler tube | 2× | Infinity Photo-Optical |

Pulsed e-jet printing [14], shown in Fig. 1, is used. The build material used in the validation is the UV-curable polyurethane prepolymer Norland Optical Adhesive (NOA) 81 (Norland Products, Cranbury, NJ). The viscosity of NOA 81 is measured as 0.453 Pa·s on an ARES-LS rheometer (TA Instruments, New Castle, DE). NOA 81 is used because of its low volatility and its large viscosity. Low volatility is required for Assumption 1. Increased viscosity slows spreading, which improves the capture of spreading dynamics by high-speed camera. The substrate is a polished silicon wafer with p-type boron doping, orientation 〈100〉, and resistivity 0.005–0.020 Ohm·cm. Stage motion is stopped during and after each pulse to image drop spreading.

Printing parameters and their corresponding drop volumes are selected to satisfy constraints imposed by high-speed digital microscopy. Drops are large enough for profile measurement, yet not so large that they spread out of the camera frame. Furthermore, the desired e-jet behavior is also limited to a single jet ejection per drop. The printing parameters are given in Table 2. The complete range of *V _{h}* and

*T*tested is shown as shaded bins in Fig. 4. For each pair of input parameters ($Vh,\u2009Tp$), 20 drops are printed and recorded for video analysis. The bottom left corner of Fig. 4, where zero drops are printed, contains low $Vh$ and low $Tp$ that are insufficient for e-jet printing. The upper right quadrant of Fig. 4 is not used because high $Vh$ and high $Tp$ inputs cause multidrop ejections that violate Assumption 1 and large drops that spread out of the image frame.

_{p}Parameter | Value | Unit |
---|---|---|

Nozzle inner diameter | 30 | μm |

Standoff height | 150 | μm |

Low voltage, $Vl$ | 525 | V |

High voltage, $Vh$, incr. of 50 | 1100–1600 | V |

Pulse width, $Tp$, incr. of 0.1 | 1–3 | ms |

Pulse period, T_{d} | $Tp+91$ | ms |

Parameter | Value | Unit |
---|---|---|

Nozzle inner diameter | 30 | μm |

Standoff height | 150 | μm |

Low voltage, $Vl$ | 525 | V |

High voltage, $Vh$, incr. of 50 | 1100–1600 | V |

Pulse width, $Tp$, incr. of 0.1 | 1–3 | ms |

Pulse period, T_{d} | $Tp+91$ | ms |

### Data Processing.

The high-speed video data are processed frame by frame with constant image processing filter settings for each of the 20 drops of each of the 45 ($Vh,Tp$) pairs. A Gaussian filter with a standard deviation of one pixel is applied to the 8-bit grayscale image, and then a binary threshold is applied. A noise filter removes connected regions smaller than a threshold area. Next, connected regions touching the image border are removed. From the single remaining region of each frame, any internal holes due to bright glare are filled, giving only the image region occupied by the drop and its reflection. The vertical coordinate of the region centroid is assumed to be the axis of reflection. Using the least squares method of Ref. [23], a circle is fit to the perimeter points above the axis of reflection. Then, the equation of the circle and the axis of reflection are used to calculate $\theta j(t)$ and $Rj(t)$ for drop *j* at frame time *t*. Figure 2 shows the identified $R1(t)$ and $\theta 1(t)$ at two times. Time is initialized by setting *t* = 0 at the first frame after jet retraction.

*t*= 0–50 ms. The temporal mean of $\theta \xaf(t)$ from $t=$ 37.5–50 ms is denoted $\theta \xafe$. Likewise, the temporal mean of $R\xaf(t)$ from $t=$ 37.5–50 ms is denoted $R\xafe$. $\Omega \xaf$ is calculated according to Assumption 2 as

### Parameter Identification.

The molecular kinetic relationship Eq. (5) is fitted to $(dR\xaf/dt)$ and $\theta \xaf(t)$ in Fig. 5, giving spreading parameters $\kappa =1.91\xd7105\u2009Hz$ and $\lambda =1.08\xd710\u22129\u2009m$. Model parameters *T* = 298 K, and $\gamma lv=39\u2009mN\u22c5m\u22121$, measured by the pendant drop method on a Ramé-Hart Model 260 tensiometer (Ramé-Hart Instrument Co., Succasunna, NJ), are used for the fit. The derivative $(dR\xaf/dt)$ is calculated as a centered difference upon a moving average of $R\xaf(t)$. Curve fitting is performed using the matlab (The MathWorks Inc., Natick, MA) Statistics and Machine Learning Toolbox fitnlm function for nonlinear least squares regression with a robust fitting option using bisquare weights [24]. A robust fit is used because the data are concentrated around the origin, where the spreading is low-speed, and the dynamic contact angle is near its steady-state value. There is a poor fit for $(cos\u2009\theta \xafe\u2212cos\u2009\theta \xaf(t))>0.5$, where less than 0.5% of the data are located. Furthermore, the extreme $(cos\u2009\theta \xafe\u2212cos\u2009\theta \xaf(t))$ data represent the initial spreading of the drop after the jet retracts. This initial spreading experiences inertial and shape change effects not captured in the model.

### Simulation.

Simulations of Eq. (6) are performed in matlab using ode45 [25] to produce one $\theta *(t)$ for each ($Vh,\u2009Tp$) pair. The same globally identified *κ* and *λ* from Fig. 5 are used in each simulation. According to Assumption 3, the simulation initial condition is $\theta *(0)=180\u2009deg$. Each simulation of $\theta *(t)$ uses the measured $\Omega \xaf$ and $\theta \xafe$ unique to its ($Vh,\u2009Tp$) pair as model inputs $\Omega $ and *θ _{e}*. An $R*(t)$ is calculated for each $\theta *(t)$ using Eq. (4).

### Results.

Three representative comparisons between the simulated and experimentally captured time evolution of *R*(*t*) are provided in Fig. 6. The first frame of measurement is the first frame after the jet has separated from the drop, so $R\xaf(0)>0$ and $\theta \xaf<180\u2009deg$. The average measured initial contact angle is $\theta \xaf(0)=72\u2009deg\u2009\xb1\u200915\u2009deg$, which is significantly smaller than the simulation initial condition $\theta *(0)=180\u2009deg$ specified by Assumption 3. The conservatism of Assumption 3 gives a prediction of drop spreading dynamics without requiring an accurate prediction of the initial contact angle. Despite the initial error, the simulation quickly evolves to track the measured values. The relative errors, *e*, of the simulation variables $\theta *$ and $R*$ are defined using root mean square error in Eqs. (12) and (13), where *i* is the first frame included in the relative error computation, $\Delta t$ is the time-step, and *N* is the final frame number

Both the simulation and measurements use $\Delta t=50\u2009\mu s$ and *N* = 1000. Relative errors are plotted in Fig. 7 for each ($Vh,\u2009Tp$) pair. Despite the significant deviation between initial measured data and the simulation's initial condition, the simulation tracks the measured data well, as seen in the inset plot of average relative error for a varying number of excluded initial frames *i* in Fig. 7. Using only one excluded initial frame, the simulated $\theta *(t)$ tracks the measured values with an average relative error below 3%, while the simulated $R*(t)$ tracks the measured values with an average relative error below 1%, as shown in Table 3.

## Conclusion

A dynamical model of drop spreading for e-jet printing is formulated using a molecular kinetic relationship coupled to volume conservation and a spherical cap drop shape assumption. The model's molecular kinetic relationship parameters *κ* and *λ* are fitted to high-speed microscopy measurements. The model successfully predicts e-jet-printed drop shape evolution over a 20-fold range of drop volume, 0.68–13.4 pL. Despite Assumption 3 of a $\theta (0)=180\u2009deg$ initial condition, which is not supported by observations, the dynamical model demonstrates high-tracking correlation with the experimentally captured time evolution of $\theta (t)$ and *R*(*t*) as evidenced by low relative errors. It should be noted that the simulation uses the measured $\Omega \xaf$ and $\theta \xafe$ as inputs, which drive the simulation to the measured steady state.

For application of this model to conditions outside of the experimental validation of this work, it is instructive to re-examine model assumptions. The model requires the process maps from inputs ($Vh,\u2009Tp$) to drop properties ($\Omega $,*θ _{e}*), which may be collected by microscope imaging. Model application also requires a prediction of the time of jet retraction, which can be obtained by high-speed imaging or by nanoampere current measurement for sufficiently conductive build materials [15]. For conditions of varying substrate wetting,

*θ*describes the spreading limit. For complete wetting $\theta e=0$, and a liquid drop spreads to a molecular film on the substrate, so the spherical cap shape Assumption 2 does not hold. However, $\theta e=0$ is undesirable for high-resolution e-jet applications, which favor partial wetting, defined by $\theta e>0$. Another possible cause of deviation from Assumption 2 is hydrostatic pressure within a large drop of size approaching its capillary length, which is 2 mm for the build material in this work. As for a lower limit for applicability of the model, molecular dynamics simulations support the molecular kinetic theory of drop spreading down to the nanoscale [26]. Experimental validation of the model with smaller e-jet-printed drops remains an interesting challenge. Furthermore, the model assumes drop volume conservation in Assumption 1. For some applications involving volatile build materials, the spreading dynamics are significantly faster than evaporation dynamics so that they can be separated into two phases. Assumption 1 is applicable to the spreading phase of a separable system. Furthermore, this work provides a baseline model from which to develop future models of e-jet drop spreading under complex conditions such as simultaneous evaporation, absorption, chemical reaction, or drop coalescence.

_{e}This work provides new insights into the physics that govern the spreading of e-jet-printed drops. High-speed images provide the first demonstration that e-jet-printed drop spreading fits a molecular kinetic relationship between contact line velocity and dynamic contact angle. Furthermore, this is the first validation of a dynamical model of drop spreading for e-jet-printed drops. Such a model is key for understanding the coalescence or noncoalescence of a drop with neighboring topography. Direct applications of this model include: (1) predicting the position, timing, and dynamic behavior of drop contact for systems with one or more spreading drops and known stationary topography; (2) facilitating process planning for drop-on-demand e-jet printing of structures in which drop contact is desired at certain locations and forbidden at other locations, such as closely spaced lines; and (3) providing a baseline model to be extended for a more advanced understanding of e-jet-printed drop coalescence. Furthermore, the wetting and spreading formalism of this model can be used to design e-jet drop volume, position, and time sequences that yield the desired topography with a tolerance for perturbations to build material and substrate wetting properties. Finally, new build material and substrate combinations can be characterized in the manner of this work to find their molecular kinetic spreading properties *κ* and *λ*, so that this model can predict the spreading of these materials with minimal experimentation cost.

## Acknowledgment

Maxwell Wu and Monica Piñon assisted in data collection. Patrick Sammons and Zhi Wang assisted in editing.

## Funding Data

Division of Civil, Mechanical and Manufacturing Innovation, National Science Foundation (Grant Nos. 1434660 and 1434693).

NSF Graduate Research Fellowship Program, Division of Graduate Education (Grant No. 1256260).