A computational model for doctoring fluid films in gravure printing

The application of printing techniques to the production of electronics is a viable approach for producing low-cost, large-area, flexible electronics systems.^1–4 Gravure printing is a high-throughput, high-resolution, scalable printing technique that has garnered interest for its use in this industry.^5–11 Resolution of sub-5 μm features has been reported at a speed of 1 m/s, which is a significant improvement over other printing techniques viable for manufacturing of high performance devices, such as thin-film transistors.^12–15 Gravure printing is a contact printing technique with the pattern defined by recessed cells in a roller or a printing plate. It involves four steps: filling the cells with ink,¹⁶ wiping the excess ink from the roll surface, ideally only leaving ink inside the cells,¹⁷ transferring the ink to a substrate,¹⁸ and finally, ink drying on the substrate.^19,20 Of these steps, doctoring excess ink is arguably the most important factor limiting feature resolution desired of printed electronics. After wiping occurs, a residual ink layer, or residual fluid film, is typically left behind. Feature size sets an upper limit on the permissible thickness of this residual ink layer, as excessive ink undermines line-edge roughness. Thus, resolution of smaller features requires a thinner residual ink layer.

Residual ink layers not only affect pattern fidelity but also the printing performance and microstructure of subsequently printed active layers due to changes in surface roughness or surface energy. In addition, the thickness of the residual fluid film is generally not uniform across the width due to doctor blade defects that can lead to characteristic streaks aligned with the print direction and nonuniform blade loading. These streaks are often defects that can be several times larger than the features and can lead to short circuits that destroy electrical systems. The doctor blade defects that create these streaks may result during fabrication or from wear that occurs during printing. For very thin residual fluid films, blade wear is accelerated due the increased likelihood of blade-to-cylinder contact. A trade-off exists between minimizing undesirable residual films and reducing the rate of blade wear. Hence, fundamental understanding of the doctoring process is important in order to achieve optimal printing.

Kitsomboonloha and Subramanian¹⁷ studied the limits of the doctoring process with experiments and reduced-order mathematical models. Their setup consisted of a silicon printing plate, motorized linear stage, doctor blade, and pneumatic actuators. The actuators applied a force on the doctor blade directed toward the printing plate as it passed over the plate (Figure 1(a)). This configuration, so-called inverse direct gravure, had been chosen to allow the precise study of thin residual liquid films. It has been shown that results from such silicon printing plates carry over to metal-based rolls suitable for high-speed roll-to-roll printing.¹² Kitsomboonloha and Subramanian effectively produced sub-micron thick residual fluid films over a wide range of material properties and operating conditions including ink viscosity, loading force, and wiping speed. They reported that residual film thicknesses increase with dimensionless speed $U*$ in power-law fashion. Pranckh and Scriven's theoretical work²¹ on a blade coating system which is typically operated at higher speed than gravure printing using a more flexible blade, also predicted a similar trend between film thickness and elasticity number N_Es, a dimensionless parameter similar in form to the U^*, that is the film thickness increases as N_Es increases. However, they did not report power law scaling and their predicted film thickness is much thicker (in the order of 10 μm) than the sub-micron thicknesses obtained by Kitsomboonloha and Subramanian's experiments. Having compared their findings to tribology literature,^22,23 Kitsomboonloha and Subramanian¹⁷ hypothesized that the underlying physics of their system are governed by a flow/structure regime known as elastohydrodynamic lubrication.²⁴ This hypothesis was due to a power-law scaling between residual fluid film thickness and a dimensionless speed parameter U^* as predicted by Dowson and Higginson.²⁴ Our work is meant to verify the hypothesis, as well as elucidate any additional key physics governing the doctoring process performed in the experiments.

So-called lubrication flow can occur between two solid boundaries whose separation is significantly smaller than the other dimensions.^25,26 The Navier-Stokes system can be simplified, integrated in the gap direction, and rearranged into the reduced-order Reynolds equation. Elastohydrodynamic lubrication (EHL) stipulates that at sufficiently small gaps, the pressure generated is large enough to locally deform materials that are typically considered rigid, such as steel.²² At even smaller gaps, solid boundaries will come in contact and the fluid film will be negligible. In this case, solid-solid contact of the boundaries dominates in the lubrication layer resulting in boundary lubrication.²⁷ Such friction is undesirable for the doctoring process as it accelerates blade wear. Mixed lubrication occurs in the transition between EHL and boundary lubrication regimes where both the fluid film and solid-solid contact is of importance in the lubrication layer.²⁸ The transition between these lubrication regimes is typically illustrated in a Stribeck curve^29,30 (Figure 2) showing a coefficient of friction as a function of a film parameter. The film parameter is taken to be the ratio of the film thickness to the surface roughness of the boundaries, and as the film parameter is varied, a change of the trend of the coefficient of friction is observed demonstrating the transition between these lubrication regimes. It is perhaps noteworthy that traditional blade coating is typically operated in the far right of the EHL regime.

EHL regimes occur under exceedingly high fluid pressures, such that the change in pressure results in variations of liquid viscosity and/or elastic deformation of the solid.^22,23 Kitsomboonloha and Subramanian¹⁷ proposed that the doctoring process lies in the piezoviscous-elastic regime, or the “full EHL” regime. In this regime, high liquid pressure causes both local deformation of the elastic material and the variation of viscosity. Clearly, a predictive model will need to account for both local deformation of the doctor blade and the variation of viscosity with respect to pressure in this lubrication regime.

Mathematical models of EHL consider the flow of fluid between boundaries and local deformation of the solid in contact with the fluid.³¹ Typically, such models are not extended to the overall structural mechanics of the solid bodies not in contact with the fluid. In experiments performed by Kitsombolonloha and Subramanian, the blade was mounted to a holder and loaded with a downward force as shown in Figure 1(a). In this configuration, the blade's motion is not confined to the wetted region but will also extend throughout the blade body as the blade may vertically translate, as well as bend under the applied force. In the most comprehensive model presented in the Mathematical Model and Methods section, we include the effect of structural mechanics of the doctor blade body.

In this paper, we develop a comprehensive model for the doctoring process and discuss the approaches necessary to implement such a model in the context of the finite element method. The results of two simplified versions of the comprehensive model as well as the full model are compared with the experimental findings of Kitsomboonloha and Subramanian.¹⁷ Last, we discuss model outcomes, uncertainties, and future work.

A comprehensive model of the doctoring process must account for the coupling between fluid lubrication and the motion and deformation of the doctor blade. In the doctoring system, the disparate length scale between the thin-film fluid domain $(<1μm)$ and the bulk structure of the blade and blade-clamping mechanism (>1 cm) presents a computational challenge. Moreover, the fluid-structural interaction manifests at these two extremes: in elastohydrodynamic lubrication, deformation of the blade tip occurs at the fluid film scale in response to high liquid pressures, and at the blade-clamping mechanism scale, the blade motion and deformation occurs in response to the applied load (including gravity). This multiscale process is best approached with a model based on reduced-order formulations that take advantage of the large aspect ratios in the fluid film. In the following developments we do so principally with a thin-film lubrication formulation.

The Reynolds lubrication equation describes the variation in pressure through the gap and can be written as³² 

where h(x) is gap height between the printing plate and the wetted surface of the blade (Figures 1(c) and 1(d)), $∇II=I−nn$ is the surface gradient, p is liquid pressure, μ is the fluid viscosity, V_w is the printing speed, $t$ is surface tangent vector, I is the identity tensor, and n is the unit normal vector relative to the printing plate. Note that Equation (1) is written in the frame of reference of the blade, in which the plate is moving from left to right, as shown in Figure 1(a).

The printing plate is considered to be rigid and perfectly flat. The doctor blade is taken to be an elastic solid and its motion governed by conservation of momentum

where ρ is density. The blade is modeled as a linear, neo-Hookean solid with its stress tensor $σs$ taken as follows:

where G and λ are Lamé coefficients, and u is the displacement field.

The deformation of the blade and the flow of the fluid are coupled through fluid-structure interaction (FSI); correspondingly, the balance of stress is satisfied at the interface of the blade and the fluid

where $σf=Ip$ is the fluid stress and n is the unit normal at the interface. For the purposes of this study, we ignore fluid shear stresses. The FSI coupling also results in the change of gap height h(x) due to blade deformation

where $h0(x)$ is the initial, undeformed gap height between the plate and the blade. The minimum gap height, h_min, is taken to be the minimum of the deformed height function.

Boundary conditions are required at the inflow and outflow boundaries for the Reynolds equation (Eq. (1)). In this work, we assume the liquid pressure at those boundaries is atmospheric, ignoring any potential free surface effects, viz.,

where w is the width of the blade and p₀ is the atmospheric pressure. These boundary conditions are valid as long as we consider the gap heights to be infinity outside the width of the blade. As for the boundary conditions on Eq. (2), the motion of the blade is constrained at the top so that it only moves vertically, which is consistent with the experimental setup¹⁷ (Figure 1(a)). In addition, a downward loading force, F, is applied at the same location (Figure 1(b)).

The system of equations outlined above is solved using Goma 6.0,³³ a multiphysics, finite element software. The blade is modeled as a three-dimensional deforming structure with a two-dimensional lubrication model implemented in shell-elements over the blade tip surface.³² In this work, we constrain the motion of the blade to two-dimensional plane strain and disallow any liquid from entering or exiting the lubrication layer on the lateral edges. These constraints effectively reduce the dimension of the model to a two-dimensional deforming structure and a one-dimensional lubrication layer.

When EHL is considered, local deformation of the blade is on the scale of the thickness of the fluid film. The fluid films produced are sub-micron in thickness, and this requires the model resolution at nanometer to 10 nm scales in the lubrication layer. A typical finite element mesh used in our simulations is shown in Figure 1(b). Because the structural mechanics of the blade does not require the same nanometer resolution, the mesh is taken to be coarser in the structure of the blade body. Mesh refinement studies were performed, and doubling the mesh density resulted in a less than one percent change in the predicted residual fluid film thickness. Given the multiscale nature of this model, direct to steady state methods for numerical simulation have difficulty converging, and the simulations were typically performed in a transient manner to reach a steady state. Once steady state at a given set of parameters was achieved, continuation in a direct to steady state mode could be used.

As the printing plate passes by the doctor blade, it leaves behind a residual fluid film undergoing a transition from a combination of Couette and Poiseuille flow profiles at the blade exit into a thinner film with plug flow profile at a distance sufficiently far from the blade. The length of the transition zone can be estimated as the distance at which the boundary layer thickness grows into the final film thickness. The thickness is set by conservation of mass, in which volumetric flow rate of the liquid inside the gap between the blade and the plate must match that of the residual fluid film

Here, Q is the volumetric flow rate per unit width in the gap and δ is the residual fluid film thickness. With this relation, we need not include the flow transition region in our analysis in order to determine the residual fluid film thickness.

Physical parameters used in the simulations are taken from those used in experiments.¹⁷ The printing speed is varied from 0.01 to 1 m/s. The blade body length is taken to be $1300μm$⁠, angled $70°$ relative to the printing direction, and the blade body width is taken to be $w=50μm$⁠. In EHL, the viscosity-pressure characteristics of the fluid have been shown to have a significant impact on the predicted fluid film thicknesses.³¹ Kitsomboonloha and Subramanian¹⁷ hypothesized that the doctoring process lies in the piezoviscous elastic regime of lubrication.²³ The influence of viscosity-pressure characteristics is considered with a Barus model³⁴ for ink viscosity

where μ₀ is the viscosity at atmospheric pressure and α is the pressure-viscosity coefficient. μ₀ is varied between 15 and 1500 mPa s and $α=4.71GPa$ as determined by Kitsomboonloha and Subramanian.¹⁷ 

Johnson²³ defined a map of EHL regimes using dimensionless parameter groups. We use the dimensionless groups suggested by Dowson and Higginson;²⁴ a dimensionless speed parameter, U^*, is one such group and is defined as

Here, we take R = 25 μm is the blade's tip radius and $E′=220$ GPa is the elastic modulus of the blade. We focus our simulations on one value of the loading force, $F=21N$⁠, as this was the primary loading force used by Kitsomboonloha and Subramanian¹⁷ that produced a power law scaling for the residual fluid film thickness δ as a function of U^*. The power law took the form

where $b≈0.75$ was determined from a least-squares fit to experimental data of measured film thickness. In addition, the Dowson and Higginson²⁴ EHL results show a power law scaling of minimum gap height h_min as a function of U^*, and that exponent will be reported when a purely EHL model is considered.

The observed experimental dependence of the residual fluid film thickness on the dimensionless speed parameter U^* as mentioned is fit remarkably well with a power-law form (Figure 3(a)). In this section, the simulation results are compared with experimental results in terms of both the predicted residual fluid film thickness and the exponent b of the power law (Eq. (10)). Since the experiments performed by Kitsomboonloha and Subramanian¹⁷ measured residual fluid film thickness, we use the predicted residual fluid film thicknesses to compare directly with experimental results. We use predicted h_min values in order to compare directly with the results of Dowson and Higginson.²⁴ When blade deformation is considered, the predicted residual fluid film and minimum gap thicknesses are presented in Tables I and II, respectively.

In this model, the blade is assumed to be rigid with lubrication flow governing the residual fluid film thickness. The hydrodynamic lift generated in the gap must balance the applied force on the blade. Here, the viscosity is assumed to be constant, i.e., $μ=μ0$⁠.

At steady state, the one-dimensional version of Eq. (1) becomes

and an integral solution for p is

Applying the boundary conditions specified by Eq. (6) sets the value of the flow rate per unit width Q

By virtue of the relation given in (7), the residual fluid film thickness δ is

This model predicts that the residual fluid film thickness is independent of both ink viscosity and print speed, contrary to experimental observations; therefore, hydrodynamic lubrication theory alone is inadequate.

In this model, the wetted region of the blade is considered to be deformable local to the lubrication flow, but the blade body is assumed to be rigid. In other words, rigid body motion of the blade holder and blade deformation away from the wetted region is not considered, which is a similar approach to that of Dowson and Higginson.²⁴ We refer to this as our “purely EHL” model. The wetted region is assumed to be bounded by rounded edges taken to be parabolic arcs as shown in Figure 1(c). In this discussion, this shape will be referred to as a “rectangular tip.”

The predicted residual fluid film thicknesses are presented in Table I. A log-log plot of δ varying with U^* is presented in Figure 3(a), which corresponds to an exponent of b = 0.76 (Eq. (10)). The predicted final gap thickness is shown in Figure 3(b). Note that the deformation of the blade tip occurs on a scale of $10 nm.$

The minimum gap heights predicted are presented in Table II. The power law for h_min as a function of U^* is then

However, the values of δ are significantly larger than the experimental values. Due to the differences in δ, the role of overall deformation of the blade body is investigated next.

Both the blade body and wetted region are now taken to be deformable, i.e., the full mathematical model presented in the Mathematical Model and Methods section is solved. The tip is assumed to be rectangular as in Sec. III B (Figure 1(c)).

The predictions of residual fluid film thickness at varying U^* are shown in Figure 4(a). Fitting these to Eq. (10) results in b of 0.54. The predicted residual fluid film thicknesses δ are larger than the majority of the median reported values but are within range of the experimental measurements.

Deformation of the blade body is shown in Figure 5(a), magnified by 10 000-fold. The blade bends and deforms on a scale of less than 20 nm, which is much smaller than the blade length of 1300 μm and difficult to observe without such magnification. The predicted final gap thickness is shown in Figure 4(b). Compared to the predicted gap thickness when the blade body is considered rigid, Figure 3(b), there is a comparable amount of blade deformation in the wetted region on the order of 10 nm. These results suggest a strong dependence on blade tip shape.

Due to a lack of experimental observations of the blade tip during the wiping process, we can only speculate the extent of the wetted region. In Secs. III B and III C, both rounded edges of the blade are assumed to be submerged in the fluid as shown in Figure 6(a). Given that the blade is angled in experiments and the thin nature of the residual fluid films produced, it is likely the leading edge is submerged, but not necessarily the trailing edge (Figure 6(b)). In this section, we investigate the sensitivity of predicted film thickness with respect to the portion of the blade submerged in fluid.

The wetted region of the blade is taken to be a parabolic arc as shown in Figure 1(d) where the width is kept at $w=50μm$⁠, so that it is consistent with the rectangular tip case. The predicted residual fluid film thicknesses are plotted in Figure 7(a). When Equation (10) is fit to these results, an exponent b = 0.76 is predicted. The predicted layer thicknesses are within the range of experimental measurements.

The predicted final gap thickness is shown in Figure 7(b) and blade body deformation is shown in Figure 5(b). The blade body bends similar to the rectangular tip case (Figure 5(a)), but the wetted region shows little deformation.

The three computational models explored in the results section each revealed a key aspect to modeling the doctoring process in gravure printing. The analysis of lubrication theory without solid deformation, or hydrodynamic lubrication, predicts that the residual fluid film thickness δ only depends on gap height h(x) and not on dimensionless speed U^*. This result indicates that the printing speed can be taken arbitrarily high with no effect on the film thickness, and this is not what is observed in the experiments. Hence, hydrodynamic lubrication theory alone is inadequate, and blade deformation needs to be included in the analysis.

The analysis of lubrication theory with addition of the blade's local deformation, i.e., elastohydrodynamic lubrication (EHL), indicates that the doctoring step in gravure printing lies in the piezoviscous-elastic regime, that is, the generation of high liquid pressures results in elastic deformation as well as the variation of viscosity in the lubrication gap. The power law exponent for h_min as a function of U^* (Eq. (10)) predicted b = 0.77. Given that elastic deformation of the blade is observed, the EHL regime is either isoviscous-elastic or piezoviscous-elastic. The exponent b = 0.77 lies closer to the exponent of the piezoviscous-elastic regime, b = 0.7,²⁴ than the exponent of the isoviscous-elastic regime, b = 0.6.³⁵ This suggests the results of the EHL analysis lie in the piezoviscous-elastic regime, which is consistent with experimental conclusions. The value b = 0.7 is obtained with Hertzian contact assumptions, which is a theory for contact mechanics based on multiple simplifying assumptions.³⁶ The model presented here considers more general fluid-structure interaction (FSI) mechanics, and this is likely the reason for the differences in the exponent b. However, the predicted δ values relying solely on EHL do not show good agreement with the experimental results (Figure 3(a)). The values predicted are significantly larger, and we hypothesized the overall motion and deformation of the blade body contributes a significant effect on the residual fluid film thicknesses when considered in addition to the EHL model.

When the overall motion and deformation of the blade body is considered in addition to EHL, the model predictions agree with the experimental results at least to within experimental error. The values of δ predicted by the model are decreased by an order of magnitude compared to the purely EHL model, which places them within range of the experimental results. This suggests that the addition of the blade body motion and deformation is necessary for predicting sub-micron thick residual fluid films. The local deformation of the blade body is small, less than 20 nm over the 1300 μm length of the blade and is likely not the primary reason for this reduction in film thickness. The inclusion of the blade body allows for both translation of the body and the generation of a moment, which is more likely the cause for the reduction in film thickness. The exponential dependence δ on U^* deviates with that of the EHL model. Because the addition of the motion and deformation of the blade body creates additional structural response to the EHL model, this difference should be expected. However, the exponential dependence still deviates from that reported from the experiments.

Finally, when the wetted region is varied to a parabolic tip, the deviation from the experimentally determined exponential trend is rectified as shown in Figure 7(a). This indicates that there is both a need to simulate the blade structure and the fluid gap. When purely EHL is considered, the predicted residual fluid film thicknesses are an order of magnitude larger than the experiments. Inclusion of the blade structure improves the results, with predicted layer thicknesses on the right order with measured values, but the sensitivity to blade shape indicates the wetted region has a significant impact on the observed experimental trends. This adds further uncertainty to the model, as little is known of from the experiments of the blade tip interaction with the fluid.

The predicted residual fluid film thicknesses still exceed experimental values, but the experimental measurements contain uncertainty. The experimental quantification of the ink viscosity is a challenge at the high shear rates, here predicted to be on the order of $107s−1$⁠. Kitsomboonloha et al.¹⁴ measured the viscosity of the ink used here up to a shear rates on the order of $104s−1$ and found small but non-zero shear thinning at these shear rates. Thus, the effective ink viscosity in the wiping experiments was slightly smaller than the ink viscosity assumed to calculate the experimental value of U^*. This slight difference leads to a small under-estimate of the experimental residual film thickness for a given value of U^* consistent with our modeling results. Furthermore, the blade shape used in simulations is an idealized version of those used in experiments. The real shape of the blade tip is a more complex than a simple parabola. This not only impacts the blade shape that should be used for the simulations but also creates uncertainty in the blade tip radius R that should be employed to calculate U^*. The actual blades also have nonuniformly distributed manufacturing defects and surface roughness, as well as wear resulting from the blade contacting the cylinder or printing plate. Kitsomboonloha and Subramanian reported experimental values as the median of 28 measurements taken from the same print.¹⁷ This reduces the influence of large defects on the reported values since the tail of the thickness distribution is ignored; however, some statistical variability remains.

The effect of varying force was also explored in the context of the full model (results not shown). The residual fluid film thickness was found to decrease with increasing force and did not exhibit a power law scaling. Experimental results also show a decrease in residual fluid film thickness with no clear power law dependence. However, there were not enough experimental results to draw any definitive conclusions given the statistic variability of the experiments. The model requires more experimental data in order to further validate this aspect.

In the current modeling work, only deformation of the doctor blade was considered. The Young's modulus of the silicon printing plate is 70 GPa, which is on the order of that of steel. This suggests that local deformation of the printing plate should also be considered in the future. However, given that little deformation is observed in the parabolic tip case, Figure 7(b), deformation of the printing plate may have little effect on the predicted residual fluid film thicknesses.

A possible important factor that we did not consider in the model is van der Waals forces in the thin film known as disjoining pressure.³⁷ At residual film thicknesses on the order of 100 nm or less, it is likely these effects would play a significant role and may offer a possible explanation for the model's tendency to overpredict experimental results.

Reported root mean square measures of surface roughness were 614 nm for the doctor blades used by Kitsomboonloha and Subramanian.¹⁷ This is larger than the predicted minimum gap thicknesses (Table II) and would likely result in solid-solid contact, i.e., boundary or mixed lubrication (Figure 2), when surface roughness is considered. Our model has yet to be extended to include solid-solid contact effects due to the current form of Reynolds equation (1), as that would result in $h(x)→0$ and produce a singular pressure field. Boundary lubrication is especially significant in the regime of $U*<10−8$ where residual fluid film thickness is very thin, i.e., below 10 nm, and experimental data do not follow the same power law as for higher values of U^*. Whilst this regime exhibits the smallest possible residual fluid film thicknesses, doctor blade wear is accelerated making it less attractive for long print runs.

The model's sensitivity to blade shape offers an additional design space to explore. In all cases considered in this work, the predicted h_min was greater than δ at a given U^* as shown in Tables I and II. The model can be used to identify families of blade shapes that maximize the ratio of h_min to δ in an attempt to avoid boundary or mixed lubrication regimes.

Residual fluid films often vary in the direction perpendicular to the printing direction due to doctor blade defects. The suppression of these streaks is crucial to achieve high-yield circuits over large areas. This variation in the direction perpendicular to the printing direction can be captured by the current model, as it is already three-dimensional, but at a significantly higher computational cost. Our finite element mesh density is highly resolved to pick up nanometer-scale deformations associated with the EHL regime (Figure 1(b)), and this mesh density would have to be extended across the printing plate for the lubrication layer. Furthermore, the experiments were not designed to capture this variation of film thickness in relation to the location of doctor blade defects. More experimental data are needed to validate this aspect of the model.

In summary, we investigated three computational models for the doctoring process in gravure printed electronics. The model that incorporates lubrication flow and both deformation of the doctor blade body and wetted region is the most accurate and has been partially validated against experiments. The predicted film thickness is sensitive to the shape of the wetted region, which adds uncertainty to the model. Although we only restricted our analysis to two-dimensional effects, the model in its current form is already three-dimensional and can be used to further investigate other aspects of the doctoring process.

This work is based upon work supported primarily by the National Science Foundation under Cooperative Agreement No. EEC-1160494. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.