Drop-on-Demand Inkjet Printhead Performance Enhancement by Dynamic Lumped Element Modeling for Printable Electronics Fabrication

The major challenge in printable electronics fabrication is the print resolution and accuracy. In this paper, the dynamic lumped element model (DLEM) is proposed to directly simulate an inkjet-printed nanosilver droplet formation process and used for predictively controlling jetting characteristics. The static lumped element model (LEM) previously developed by the authors is extended to dynamic model with time-varying equivalent circuits to characterize nonlinear behaviors of piezoelectric printhead. The model is then used to investigate how performance of the piezoelectric ceramic actuator influences jetting characteristics of nanosilver ink. Finally, the proposed DLEM is applied to predict the printing quality using nanosilver ink. Experimental results show that, compared to other analytic models, the proposed DLEM has a simpler structure with the sufficient simulation and prediction accuracy.


Introduction
The ability of inkjet technology to deposit materials with diverse physical and chemical properties on a substrate has made it an exciting technology for mechanical engineering, life sciences, and electronics industry [1].Low cost metal coating and rapid prototyping are popular applications of inkjet technology in mechanical engineering [2].In the life science field, it has been used for printing DNA structures and making artificial skin and other important organs by jetting live cells [3,4].In the past decade, electronics fabrication using the drop-on-demand (DoD) piezoelectric (PZT) inkjet printhead to print conductive ink on flexible substrates is an enabling technology that will provide desired high-volume, low-cost production of flexible electronic devices, ranging from radio frequency identification (RFID) tags, solar cells, and organic light emitting diodes (OLED) to packaging flexible electronic devices such as healthcare equipment [5][6][7].
A typical DoD piezoelectric inkjet printhead comprises a number of ink channels that are arranged in parallel [8].Each ink channel is connected with a piezoelectric actuator.When applying a standard actuation voltage pulse, the actuator can generate pressure oscillations inside the ink channel to push the ink droplet out of the nozzle.Nowadays the developments of DoD printing are moving towards higher productivity and quality, which require adjustable small droplet sizes fired at high jetting frequencies.Generally, the print quality delivered by an inkjet printhead depends on jetting characteristics, that is, the jetting direction, the droplet velocity, and the droplet volume.However, in order to meet the challenging requirements of printable electronics fabrication, the jetting characteristics of the conductive ink have to be tightly controlled for higher inkjet performance [9].The performance improvement, which is associated with the design and operation of DoD printhead, is severely hampered by several operational issues, such as residual vibrations and cross talk [10].

Mathematical Problems in Engineering
Researchers have developed many approaches to solve these issues.The most typical approaches include using analytical or numerical techniques to model the DoD piezoelectric inkjet printhead [12,13].Based on solving the nonlinear Navier-Stokes equations, numerical modeling of droplet formation process has been focused on the liquid filament evolution with known pressure input boundary conditions, assumed to be either known or an input correction for certain specific droplet generator.Numerical models are often programmed in CFD packages using complex meshes and solving techniques.This approach provides visualized information about the acoustic pressure wave travelling inside the channel [14].However, numerical models are very complex and therefore computationally expensive.In the last decade, analytical modeling approach is gradually used to describe the ink channel dynamics, although the accuracy of analytical models is lower than numerical models [15].Based on several assumptions and simplifications, analytical models can provide simple and timesaving ways to modelbased control droplet generator with a sufficient accuracy.In this category, lumped element modeling [16] and feedforward two-port network modeling [17] technique are two typical analytic printhead models to analyze the acoustic-mechanical behaviors of droplet generator.
The lumped element modeling (LEM) approach introduces an equivalent electric circuit including capacitors, resistors, and inductors in series or parallel to describe the dynamics of the ink channel [19].The LEM results explain the driving force difference between different operating conditions.Lumped element modeling technique has been a suitable method to simulate the operating mechanism of droplet generator due to its outstanding advantages such as moderate model complexity and high simulation speed.However, the accuracy of these classical LEM models is low for the improved performance requirements of printable electronics fabrication [20].The two-port model is proposed to describe the ink channel dynamics utilizing the concept of bilaterally coupled systems [21].Each part of droplet generator is, respectively, modeled as a two-port subsystem and finally cascaded into the whole system.Compared with LEM technique, two-port network owns the advantage of satisfactory simulation accuracy and the disadvantages of more complex structures and inevitably low simulation speed.
This work is aimed at building a dynamic LEM (DLEM) model to characterize nonlinear behaviors of piezoelectric inkjet printhead.In the proposed DLEM, the linear equivalent circuit of piezoceramic in original LEM is replaced into nonlinear time-varying capacitance, inductance, and resistor in series.Compared to other analytic models such as twoport network, the DLEM has a simpler structure with the sufficient simulation accuracy.Additionally, the DLEM is used to search the optimal combinations of high-frequency driving waveform to effectively minimize the residual vibration for nanosilver ink.With comparisons of simulation and experimental results, the significant merits of the proposed methods lie in the following aspects.
(1) By applying dynamic nonlinear structure, the simulation accuracy of LEM is remarkably improved.That is, the simulation curves of droplet volume and velocity can well overlap the experimental measurements.
(2) The obvious overshoot observed in the outflow curves can sufficiently indicate the break-off phenomenon, which divides the droplet formation process into the meniscus protruding, filament breaking-off, and droplet falling-down stages.
(3) The proposed DLEM can effectively predicate optimal combinations of high-frequency driving waveform with high jetting characteristics.
The rest of this paper is organized as follows.In Section 2, a review of original lumped element modeling droplet generator is present, and three equivalent circuits are given to characterize the nonlinear behaviors of piezoelectric printhead.Next, a schematic of experimental method is introduced in Section 3.After that, the simulation results of droplet volume and velocity decide the most suitable model in Section 4. The prediction under restrictive condition is built in Section 5. Finally, discussing and concluding remarks are collected in Section 6.

Original Lumped Element
Model.LEM provides a simple method to estimate the low-frequency response of printhead, which assumes that the characteristic length scales of the governing physical phenomena are larger than geometric dimension.For example, the pressure wave length (in the order of 100 mm) of KM series IJ Head of Konica, which is determined by the sonic speed of the fluid (about 1400 m/s) and the highest operation frequency of the droplet generator (12.8 kHz), is much larger than nozzle diameter (27 𝜇m).So this assumption for LEM application is well satisfied.
A classical LEM is given to simulate jetting characteristics of PZT printhead by Gallas, as shown in Figure 1.The equivalent circuit model is constructed by the energy storage elements and ideal dissipative terminal.In this electroacoustic system, pressure and voltage are independent variables, while current and volumetric flow rate are dependent variables.Model structure shows that the energy converts from electrical energy to mechanical energy, then to fluidic/acoustic energy, and finally to kinetic energy, as described in Figure 1(b).The droplet generator structure can be characterized by equivalent acoustic mass (representing stored kinetic energy) and acoustic compliance (representing stored potential energy), in which the corresponding equivalent circuit models are supported by various fluid mechanisms, as represented in Figures 1(a) and 1(c).Furthermore, the piezoceramic model is constructed based on the electrofluidic/acoustic theory [22].The neck model is built on the velocity profile function [23].The nozzle model is constructed according to the end correction for an open tube theory [24].
As the equivalent circuit model shown in Figure 1, excitation voltage  ac is applied to piezoelectric ceramic to create mechanical deformation.Coupling coefficient   represents a conversion from mechanical domain to acoustic domain. eb is the blocked electrical capacitance of the piezoelectric material.In the acoustic domain,  aD and  aC represent the acoustic compliance of piezoceramic and channel. aD ,  aN , and  aO are the acoustic resistance due to structural damping, neck tapering, and fluid flowing out of nozzle, respectively. aD ,  aN , and  aRad represent the acoustic mass of piezoceramic, neck, and nozzle in proper order, respectively.
The dimensions of droplet generator and physical properties of nanosilver ink provided by the LEED-PV Corporation are listed in Table 1, and calculation formulas of LEM model parameters are provided in Table 2.
According to above lumped element model, a transfer function describing the response relationship between the flow rate  out and input AC voltage  ac is obtained as follows: where The frequency response of droplet generator driven by the sinusoidal signal, as shown in Figure 2, reveals that natural frequency (around 48.6 kHz) is higher than operating frequency (below 12.8 kHz).This means that the meniscus can obtain sufficient energy to overcome viscosity and surface tension and then to generate consistent droplets.
As the above simulation models,  aD ,  aD , and  aD are regarded as static parameters based on the approximate linear behaviors of piezoelectric ceramic with a fix path of P-E hysteresis loop.However, in fact the piezoelectric ceramic is a nonlinear material with multiple electric field/strain curves (parallel and perpendicular to the field), as shown in Figure 3.It has been proved that the above original static model has insufficient ability to simulate the piezoelectric ceramic dynamic response in time-domain below linear voltage [11].
Our research finds that the maximum error between the experimental measurements and simulation results can reach eight percent.Additionally, when the break-off phenomenon , where  0 is volume of the cavity Piezoelectric , where  is the elastic modulus, V is Poisson's ratio, and ℎ is the thickness , where  is the experimentally determined damping factor [18] Nozzle where  1 is the Bessel function of the first kind and  = / 0 , where  2 is electroacoustic coupling factor

Improved LEM with Nonlinear Equivalent Circuits.
In order to resolve above problems, in this section, three equivalent circuits with different characters are incorporated into the classical LEM to characterize the nonlinear effect of the piezoelectric ceramics.
(1) Proportional Model.In order to describe the nonlinear relationship between electric field and strain, a proportional model (illustrated in Figure 4) as a resonance was defined in the IEEE standards publication [25].This analytical model is derived based on the constitutive relations of piezoelectricity and the following assumptions: (a) Euler-Bernouli beam assumption [26]; (b) negligible external excitation from air damping; (c) proportional damping (i.e., the strain rate damping and viscous damping are assumed to be proportional to the bending stiffness and mass per length of the beam); (d) uniform electric field through the piezoelectric thickness.Then the circuit element values of the proposed lumped parameter model are formulated in (1), where T is the adjustment coefficient: (2) Multiple-Resonator Model.Yang and Tang [27] give the multiple-resonator model to characterize piezoelectric ceramic in piezoelectric energy harvesters, as shown in Figure 5. His research finds that the network of piezoelectric ceramic driving circuit consists of infinite branches, and each is composed of simplified equivalent circuit with an inductor, a capacitor, and a resistor.When the voltage signals from several branches drive the piezoelectric ceramic simultaneously, the higher order vibration modes are excited.It is realized that working at higher order modes causes the nonlinear effect of piezoelectric ceramic.To characterize the different higher order modes, a circuit of multiple resonators is constructed to match the external branches.In this literature, to model the droplet generator, nearly a double circuit may meet the simulation accuracy.Then the circuit element values of the proposed lumped parameter model are formulated in (3), where  and N are both the adjustment coefficients: (3) Dynamic Model.Based on the works in [28,29], an equivalent electrical circuit with nonlinear resonance oscillations in series, namely, the dynamic model, is given to analyze power BAW (bulk acoustic wave) resonators, as shown in Figure 6.Three major nonlinear effects are considered in this equivalent circuit: the amplitude-frequency effect (A-F effect) [30], the intermodulation effect [31], and the bias-frequency effect (B-F effect) [32].The dynamic values of resistance and capacitance in the motional branch are dependent on the motional current in the branch.Then derivation formulas are given in (4), where A 1 , A 2 , B 1 , and B 2 are the adjustment coefficients: Then, the improved LEM structure can be illustrated as in Figure 7.

Experimental Method
3.1.Experimental Setup.In this experiment, the investigated printhead is KM series IJ Head of Konica and the used ink is nanosilver ink provided by the LEED-PV Corporation, and the according parameters are listed in Table 1 in Section 2. A schematic overview of the experimental setup to observe the nanosilver droplet formation process is shown in Figure 8.The visualization system consists of a cold source of halogens (XAO LG150) and a CCD camera (SVS ECO424), which can synchronize to the driving signal.The volume and velocity of droplet and the profile of meniscus are measured through the images captured by CCD with exposure time of 1 us.The high-frequency driving waveform is programmed by arbitrary wave generator (RIGO DG2014A).In order to repeatedly generate uniform droplets, ink supply unit connects air pressure unit, which offers suitable negative force to prevent the ink leaking from the printhead.

Edge Detection and Image Analysis as a Qualitative
Measure.The leading edge detection is to determine the nozzle position and meniscus position.A straight line is generated along the nozzle in the region of interest (ROI) with CCD camera image analysis.By setting a suitable threshold gradient of brightness, the meniscus location can be detected.The satellite droplets and filament can be detected through farther image processing such as cropping, filtering, and threshold recontrasting (shown as in Figure 9).In image analysis system, the droplet volume is calculated by which hypothesizes that liquid jet is axis-symmetric along nozzle centerline.And the droplet velocity is calculated based on the position change of the moving-forward leading edge in two consecutive images with specified sampling interval.Note that the information is lost if the meniscus is hidden inside the nozzle, so edge detection only measures the protruding meniscus.

High-Frequency Driving Waveform as a System
Excitation.The high-frequency driving waveform is composed of positive and negative periods including three pulse edges and two dwell times, as shown in Figure 10.The first rising edge expands the channel to get the ink meniscus into the channel.
Then the falling edge works to eject the droplet from the nozzle.This period of waveform design is based on the half period of channel resonance.Thus, the waveform design issue is focused on calculating the optimum value of dwell time  1 =  0 /, where  0 is the length of channel and  is the sonic speed of the used ink.The falling edge of the second negative voltage generates larger pressure in channel to accelerate the droplet with higher speed than conventional single trapezoidal waveform.The second rising edge, which sets the voltage back to 0, works to restrain the meniscus vibration [33].
Theoretically, the suitable slops of the rising edge and falling edge can restrain the meniscus vibration [34].However, in practical applications, these two values are too large to control accurately.
As mentioned above, theoretical optimum value for dwell time  1 is set as 5.2 us.The corresponding dwell time  2 is recommended as 10.4 us (twice of  1 ) and the voltage magnitude is set as 12 V according to the instruction of Konica printhead.In this paper, the driving waveform with these parameters is called the "standard" waveform [35].

Simulation Result
4.1.Coefficient Identification.The coefficients of original LEM model can be determined experimentally.However, the three improved models in this work are somewhat complicated, and their coefficients reveal some important features of these equivalent circuits.For this reason, this section is dedicated to identify the adjustment coefficients in these corresponding circuits.
In this experiment, a group of experimental measurements with printhead driven by the standard high-frequency waveform are chosen as the basic reference.The coefficients of each proposed improved LEM model should be adjusted to minimize the deviations between the simulation results and the selected experimental measurements shown in Figure 11.The nonlinear least squares method [36], using the Levenberg-Margquardt and Gauss-Newton Algorithms, is applied to find the optimal coefficients set for each proposed equivalent circuit, which makes the simulation curves overlap the experimental measurements as far as possible.
The resulting coefficients are listed in Table 3.

Droplet Volume Simulation.
In many applications of inkjet printing, the volume of droplet is a critical restrictive condition.Droplets with controllable volume affect not only the printing resolution but also the depositional thickness.Therefore, how to accurately control droplet volume is a key problem to be solved by the proposed method.According to experimental measurements in Figure 11, the comparison between the original LEM model and the proposed three improved models is presented in Figure 12.In this figure, the solid line is the outflow simulation curves based on the LEM technique and the star line is the actual measured volume values during the droplet formation process.It should be noted that, in Figures 12(a)-12(d), the first three points of actual measurement curve are default to zero because no fluid flowing outside the nozzle can be captured by our image system.
As the actual experimental measurements indicated in Figure 11, the droplet formation process is divided into three stages: the meniscus protruding stage (i.e.,  1st = 15∼30 us), the filament breaking-off stage (i.e.,  2nd = 30∼35 us), and the droplet falling-down stage (i.e.,  3rd = 35∼60 us).That   is, ideally, the outflow curve simulated by analytic models should include overshoot phase, which can clearly demonstrate that the part of filament is inhaled back into nozzle or splits into satellite droplets by viscosity and surface tension in the filament break-off stage.However, according to the simulation result of the original LEM in Figure 12(a), there is no obvious decreasing after reaching the maximum in the outflow curve, and the outflow curve of the original model cannot well overlap the actual measurement curve.
In Figures 12(b) and 12(d), the overshoot phase can be observed clearly.This performance improvement can be straightforwardly explained by the fact that poles and zeros of the transfer functions are manipulated to speed up the dynamic response speed of the proportional and dynamic models.From Figure 12, both the multiple-resonator and dynamic models can obtain precise simulation results.That is, more well-designed coefficients are used to characterize the nonlinear behaviors of piezoelectric printhead and the adjustable ranges are extended to a satisfied degree of both models.Additionally, it should be noted that the outflow curve of dynamic model has not only the overshoot phase but also high overlap ratio to the experimental measurement.Moreover, the simulation curves reveal that the meniscus moves inward before the droplet appears.

Flow and Droplet Velocity Simulation.
The velocity of droplet is the other important aspect to evaluate the printing effect.Droplets with a high speed can remarkably improve the printing efficiency and the quality of impact with substrate.Therefore, this experiment gives the velocity simulation using the proposed LEM models on the same actual experimental measurement in Section 4.1.
The comparisons between the flow velocity curves simulated by the original LEM and three improved models are illustrated in Figure 13.The excitation is the standard driving waveform.For the actual measured velocity curves, the first two points are also default to zero because no fluid flowing outside the nozzle can be captured by our image system.
It should be noted that all the values of actual measured velocity are positive.However, from the flow velocity curves in Figure 13, the values of flow velocity in nozzle can be positive or negative in different phases.This positive-tonegative velocity variation indicates the changes of the flow direction that leads to the break-off phenomenon.Hence, in order to correctly simulate the velocity in leading edge, we propose a new calculation formula which combines both of outflow curve and flow velocity curve: where  out is the outflow value and  is the flow velocity. 1 is the moment when the sum of  out changes from negative to positive in outflow curve and  2 after  1 is the moment when flow velocity changes to negative in nozzle.
From Figure 13, it can be obviously observed that, in the beginning of droplet formation process, the multipleresonator and dynamic models exhibit better performance in simulating the droplet velocity rather than the original model.That is, the average velocity curves obtained by these two models are almost coincident to the actual measured velocity curves.
Based on the above analysis and comparisons, the dynamic model based LEM has the most accuracy in simulating the droplet volume and velocity and can carry sufficient information to distinguish the break-off phenomenon occuring in the nanosilver droplet formation processs.Therefore, among these three improved models, the dynamic lumped element model (DLEM) is the most suitable improved equivalent structure for simulating jetting characteristics.

Prediction
For industry applications in printable electronics fabrication, the printhead must work at a certain status to meet many restrictive conditions.Among them, the most important  problem is choosing the appropriate combinations of the driving waveform parameters for the used conductive inks.However, an exhausting manual selective process inevitably wastes a lot of time.Therefore, the computer-aided methods are urgently needed for searching these appropriate combinations efficiently and robustly.In this section, the proposed DLEM is applied to predict jetting characteristics for nanosilver ink.Step 1.According to droplet generator parameters and properties of fluid, calculate most of the parameters in DLEM.
Step 2. According to experimental measurements of droplet volume driven by a standard waveform, determine the rest of the adjustment coefficients and damping coefficients.
Step 3. Based on the constructed DLEM, predicate the values of droplet volume and velocity.
Step 4. According to the need of applications and restrictive conditions, find out several groups of the appropriate parameters of driving waveform from predictive data lists of volume and velocity.
Step 5.After making a subtle adjustment in practice, decide the most appropriate parameters of driving waveform.

Prediction Case 1: Standard Waveform with Different
Dwell Time.As mentioned above, the DLEM model is chosen to simulate jetting characteristics for the nanosilver ink.Here, an illustrative example is given to display the effect of predicting the typical combinations of waveform parameters with different dwell times ( 1 = 3∼8 us,  2 = 2 1 ) and the same amplitude (12 V).The properties of nanosilver ink provided by the LEED-PV Corporation have been listed in Table 1.
In Figure 14, the predicted outflow and average velocity curves with various dwell times are depicted.We can observe from Figure 14(a) that, before the droplets rush out of the nozzle, the depth of the inhaled meniscus is proportional to dwell time.From Figure 14(b), the time of the droplet appearing is postponed with the increase of dwell time.That is, both two simulation phenomena consist with the theoretical analysis results in [37].
The comparisons between actual measured and predicted jetting characteristics with various dwell times are given in Figures 15(a 4, the deviation of the predicted jetting characteristics is within acceptable range.

Prediction Case 2: Driving Waveform under Restrictive
Condition.In this experiment, an illustrative example under restrictive conditions is presented to demonstrate the feasibility and effectiveness of the proposed method for predicting jetting characteristics.The dimensions of droplet generator and properties of the nanosilver ink are listed in Table 1.In order to reach a high printing resolution and good quality of droplet impact with the substrate, the restrictive conditions involve that the droplet volume is smaller than 14 pL and the droplet velocity is larger than 5 m/s.Based on the prediction procedure, the predictive values of drop volume/velocity with various dwell times  1 and  2 are listed in Table 5.According to the restrictive conditions, three combinations of parameters are found out:  1 = 5us/ 2 = 13us,  1 = 5.5 us/ 2 = 14us, and  1 = 5 us/ 2 = 7.5∼8 us.After a small range search around predictive values in the actual test, the combination of  1 = 4.9 us/ 2 = 7.7 us is finally chosen.The dynamic effect of jetting characteristics driven by this combination is shown in Figure 16.From this figure, the jetting characteristics satisfy the specified restrictive conditions quite well.Meanwhile, almost no satellite droplet emerges after jetting the main droplet.
The 3D plot of predictive jetting characteristics with various dwell times  1 and  2 is shown in Figure 17.With the increase of dwell time  1 , peak positions of volume and velocity have a tiny shift in 9.5∼11.5 us.With the same dwell time  2 , there is also a parabolic profile in predictive volume and velocity.

Conclusion
The lumped element modeling based on linear model of piezoelectric ceramic has been successfully used to simulate piezoelectric droplet generator.However, the nonlinear model of piezoelectric ceramic is necessary to diminish deviation between the simulated and measured results.In order to characterize the nonlinear behavior of piezoelectric It should be noted that the proposed nonlinear model is focused on relationship between the droplet volume/velocity and the dwell time of driving waveform, while the voltage amplitude is also a decisive factor to droplet formation process.This should be discussed in our future work.

Figure 1 :
Figure 1: Schematic overview of lumped-element modeling for PZT printhead: (a) droplet generator structure, (b) model structure, and (c) equivalent circuit model.

Figure 2 :
Figure 2: Frequency response of lumped element model.

1 Figure 4 :
Figure 4: Equivalent electrical circuit model of PZT nonlinear effect with proportional coefficient.

Figure 5 :
Figure 5: Equivalent electrical circuit model of PZT nonlinear effect with multiple parallel resonances.

Figure 6 :
Figure 6: Equivalent electrical circuit model of PZT nonlinear effect with dynamic coefficients.

Figure 11 :
Figure 11: Pictures of nanosilver droplet formation process driven by standard waveform.

Figure 12 :
Figure 12: Outflow change of measurements and simulations of (a) original model, (b) proportional coefficient model, (c) multiple resonators model, and (d) dynamic coefficient model.

Figure 14 :
Figure 14: Predictive droplet volume and velocity with different and various dwell times  2 = 2 1 : (a) simulation curves of outflow change and (b) simulation curves of droplet average velocity.

Figure 15 :
Figure 15: Predictive droplet volume and velocity compared with measured volume and velocity.
) and 15(b), which both illustrate a parabolic change rule with the increase of dwell time.It should be noted that, at the dwell time  1 = 5.6 us, near the theoretical optimum value   1 = 5.2 us, the droplet volume and velocity simultaneously reach the maximum.From the performance of prediction accuracy listed in Table

Figure 17 :
Figure 17: Predictive droplet volume and velocity with various dwell times  1 and  2 : (a) predictive droplet volume and (b) predictive droplet velocity.

Table 3 :
Coefficients for three lumped element models.

Table 4 :
Performance of prediction accuracy.In this work, the droplet volume/velocity curves simulated by dynamic model can carry sufficient information to distinguish nanosilver droplet formation process clearly.Although the complexity of system increases, the outflow curves can obviously reflect the filament breakingoff stage that occurs, which is coincident with sampled images.Such nonlinear LEM model is also successful in predicting jetting characteristics under restrictive conditions.Comparing the droplet volume/velocity curves with different dwell time of high-frequency driving waveform, the linear relationship between dwell time and meniscus appearing moment is discovered.Meanwhile, a parabolic change rule of volume/velocity with the increase of dwell time is in good agreement with present theoretical analysis.
Defaults are because the combinations of parameters in driving waveform are out of the operation range.