Melt Electrospinning Writing Process Guided by a “Printability Number”

The operating principle for the well-known electrospinning (ES) technology is based on the electrohydrodynamics (EHD) phenomenon [1,2]. This allows the electrostatic drawing of continuous liquid jet fibers under a divergent electrical field which is applied between a positively charged spinneret and a grounded stationary or moving collector plate [3]. ES is possible not only for polymer solutions but also for the processing of polymer melts [4,5]. In the latter case, the process is entitled melt electrospinning (MES). Various MES system configurations (laser MES [6–9], melt differential ES [10,11], gas-assisted MES [12], and twin-screw extrusion-based MES with multinozzle spinnerets [13,14]) have been reported. These systems have led to the production of nonwoven fiber meshes made from a wide variety of thermoplastic polymers [15] and characterized by fiber diameters ranging from the submicron to micron scales [16].

It has recently been demonstrated that MES can be performed in a direct writing mode by mounting the grounded collector plate on a dynamic mechanical device (x–y–z translational stage and/or using a rotational collector). This novel process is designated as melt electrospinning writing (MEW) [17]. The enabling MEW process harnesses both the ability of the EHD methods for generating fibers with diameters that are less than 100 μm and the precise fiber placement accuracy of the toolpath-controlled additive manufacturing (AM) methods, as indicated in Fig. 1. The latter is achieved by processing highly viscous, low-conductive, and solvent-free polymer melts that experience large straight jet paths [18]. Thus, the process can be regarded as a hybrid version of the melt extrusion-based AM methods (fused deposition modeling [19] and precision extrusion deposition [20]) in tandem with the various EHD-based methods (“near-field” electrospinning or E-jet printing) inspired by solution ES technology [15,21–23].

The physical mechanisms governing the MEW process can be understood by dividing the process into three discrete regimes [24], as shown in Fig. 2(a): (a) the polymer melt supply regime composed of the syringe and the needle tip, (b) the free-flow spinline regime, and (c) the collector plate regime. Initially, the thermoplastic polymer is maintained in the supply regime under a uniform melt state. The polymer melt is then extruded under a prescribed volumetric flow rate with an applied voltage potential on the positively charged needle tip. When the polymer melt enters the free-flow regime and the coulombic electrostatic forces acting at the jet-ambient air interface overcome the surface tension (along with the viscous and elastic forces), a straight, thin jet is formed between the needle tip and the collector plate [25]. In the spinline regime, the polymer melt is stretched due to its viscoelastic nature and solidified when the ambient temperature conditions are maintained below the polymer glass transition temperature [18,26,27]. The solidified fiber is deposited on the grounded collector plate, which is mounted on an automated x–y stage. Thus, the fiber can be printed in a direct, toolpath-controlled mode as a function of the stage translation.

When the process is performed under a stationary collector (zero translational stage speed—(U_T)), the jet initially experiences bending instabilities due to the occurrence of compressive buckling forces developed when the jet impacts the grounded collector plate [26,28,29]. As the process evolves (Figs. 2(b-I)–2(b-III)), attractive and repulsive electrostatic forces lead to the formation of conical fiber structures (Fig. 2(b-III)) [30]. When the process is performed under a moving collector, there exists a critical stage speed (U_CR) above which the bending and whipping motion instabilities are counterbalanced by the mechanical stage forces, resulting in the printing of aligned fibers (Figs. 2(b-IV) and 2(b-V)) [17,31]. When the jet speed approximates the critical stage speed, the result is precise localization of a fiber deposition point directly below the spinneret (Fig. 2(b-IV)) [25]. Finally, when the translational stage reaches speeds greater than the critical stage speed, a trailing edge is formed due to the viscoelastic nature of the material, which yields a time lag between the electrospinning and direct writing modes for fiber placement on the collector plate (Fig. 2(b-V)) [17,25].

The important process parameters of the MEW process are shown schematically in Fig. 3 and consist of parameters associated with the material, the process, and the geometry of the system. Generally applied trial and error methods to identify the efficient printability range are hampered by the significant number and interactions between the parameters. To reduce the number of independent variables, dimensionless numbers can be derived from the nondimensionalization of the conservation and constitutive equations.

The nondimensional governing equations of the electrospinning process, i.e., the conservation equations for mass, and momentum, in conjunction with the electric field equation for viscoelastic fluids have been presented by Feng [32] for solution electrospinning. Zhmayev et al. [18] have modified Feng's approach for the case of polymer melts on the basis of their relatively low electrical conductivity (K < 10⁻¹⁰ S/m). The nondimensional governing equations are detailed in the Appendix, and the derived scaled dimensionless numbers along with the characteristic quantities and their physical meaning are depicted in Table 1.

The superior placement accuracy compared to the chaotic fiber deposition observed in solution ES and the ability to generate fibers with diameters less than 100 μm, along with the solvent-free nature of the MEW process, are key characteristics of the MEW process [24]. These characteristics have promoted the use of MEW as an ABM technique for the printing of 3D biomaterial-based porous constructs, termed as scaffolds, that can be populated with cells to guide cellular function [33–35]. For example, 3D poly(ε-caprolactone) (PCL) fibrous scaffolds with lattice architecture and fiber diameters ranging from 0.8 to 40 μm have been printed via MEW [17,26,36–48]. These studies have mainly targeted scaffold-guided tissue engineering (TE) applications, with the reported micron-scale scaffolds promoting cell viability and demonstrating cell-invasive and favorable mechanical properties. Thus, MEW has the potential to furnish a new paradigm in the ABM field, since it can fill the missing dimensional scale window, whose lower and upper bounds are currently defined by more traditional EHD-based and AM methods, respectively (see Fig. 1).

MEW process development is still in its infancy and a step during the manufacturing process, whereby tuning single parameters in turn via a trial and error mode, is required until the desired scaffold geometry is obtained. This major hurdle in the MEW process stems from the large number and complex interplay of the parameters that need to be tuned. This, along with the prevailing variability between the reported custom MEW systems, makes deterministic parametric studies, where the effect of single parameters is studied by keeping the rest of the independent parameters constant, time consuming and not useful for optimization toward achieving “efficient printability.” Efficient printability is defined as the printing of layered meshes with well-defined microarchitectures marked by consistent fiber diameter and pore shape.

In the present study, dimensional analysis is used to overcome the previously mentioned roadblocks toward efficient printability by identifying dimensionless groups based on measurable polymer properties and controllable process parameters that govern the nature of the involved physical forces. The product of these dimensionless groups is defined as the Printability Number. The number is correlated with existing scaled dimensionless numbers obtained from the governing equations and provides guidance for tuning the processing parameters of MEW toward the fabrication of high-fidelity meshes. The proposed number is validated using a custom-designed MEW system through monitoring of the charged polymer jet in the free-flow regime and characterization of the final mesh morphology. Printing of meshes with well-defined microarchitectures is governed by steady electrospinning conditions in the free-flow spinline regime by balancing the applied voltage and the volumetric flow rate, and operating under a critical translational stage speed.

Poly(ε-polycaprolactone) (PCL) is selected for MEW on the basis of its Food and Drug Administration approval for in vivo applications, biocompatibility, long-term biodegradability, and relatively low and wide melt processing temperature window (60–90 °C). PCL material specifications with an average molecular weight of 45,600 g/mol and polydispersity of 1.219 are obtained from Perstorp Ltd. of Warrington, UK (Capa6500).

PCL pellets are molded into 8 mm and 25 mm circular disks using aluminum shims between Teflon surfaces and a Carver press at 120 °C for subsequent rheological characterization. The advanced rheological extended system (ARES) of Rheometric Scientific (currently TA Instruments) is used in conjunction with stainless-steel parallel disk fixtures with 25 mm diameter for small-amplitude oscillatory shear (SAOS) and steady torsional flow experiments. The force-rebalance transducer of the rheometer is capable of measuring simultaneously both the normal force and the torque. The oven temperature of the rheometer is controlled within ± 0.1 °C. The rheological characterization experiments are carried out at 70 °C, 80 °C, and 90 °C and using a constant 1 mm gap.

A high-resolution heat-assisted MEW system configuration is established. The process design is guided by detailed characterization of the thermorheological processing properties of the biomaterial substrate [31,49] along with the fluid dynamics, heat transfer, and electrostatics multiphysics phenomena governing the process under investigation. The overall system configuration is analyzed based on three defined discrete process regimes.

The polymer melt supply regime is composed of a glass Luer-lock 5 ml syringe (Hamilton, Reno, NV) and a stainless-steel needle tip with a plastic hub (McMaster Carr, Elmhurst, IL) attached to it. The polymer melt is maintained in a uniform melt state using an industrial heat gun (Steinel, HG 2510 ESD). In addition, a programmable syringe pump (Harvard Apparatus, Holliston, MA) is mounted vertically and used to set the volumetric flow rate by adjusting the speed of the plunger within a syringe (flow accuracy within 0.25% and reproducibility within 0.05%). The temperature is monitored both at the syringe barrel and the capillary tip with an infrared FLIR thermal camera (PM 290, Inframetrics, Thermacam).

In the free-flow regime, a high-voltage source (Gamma High Voltage Research, Ormond Beach, FL) is used for the application of a voltage potential between the needle tip and a grounded electrically conductive collector. An aluminum collector is mounted on an x–y programmable stage (ASI Applied Scientific Instrumentation, Eugene, OR) that is sequentially mounted on a lab jack (Newport Corporation, Irvine, CA) (Figs. 4(b) and 5(a)). The distance between the tip and the collector plate is monitored using a vertical digital meter (Figs. 4(b) and 5(a)) and set manually using the lab jack's turning knob with a vertical positioning resolution of 0.5 mm. To compensate for ambient conditions that might affect the process, the overall system configuration is placed on an antivibrating optical table with the spinning apparatus contained within a plexiglass enclosure. Furthermore, the temperature and humidity values within the enclosure are monitored using a multimeter (Extech Instruments, Waltham, MA) equipped with a type K thermocouple.

The heating element is composed of an industrial heat gun (HG) with controllable air flow (Q_HG) (0.002–0.008 m³/s) and adjustable air temperature (T_HG) settings (49–649 °C). The heat gun is mounted at the entrance of a heating tunnel housed by a transparent chamber constructed out of poly (methyl methacrylate) (Figs. 4(a) and 5(a)). The syringe passes through the heating tunnel, and a small portion of the syringe needle tip reaches the interior of the chamber through an electrically conductive tape covering a circular opening created at the ceiling of the chamber (Fig. 4(b)). Heating insulation tapes are applied onto the back wall and the floor of the heating tunnel in order to minimize heat losses. The area of the circular opening covered by the tape is kept tightly sealed in order to avoid disturbances along the spinline regime from the hot stream air.

The surface of the syringe is heated due to heat transfer via forced convection generated by the heat gun, and the ambient temperature conditions along the spinline are governed by free convection through to the heated tape (Fig. 4(b)). The heat transfer conditions are calibrated so that the temperature at the surface of the syringe hosting the PCL melt is maintained as the desired temperature. For example, it is determined that for the air flow rate of Q_HG = 0.0017 m³/s and air temperature of T_HG = 132 °C, the temperature on the syringe surface (T_s) is set and maintained at 78 ± 1 °C (Fig. 5(b)). Thermal imaging using the FLIR camera confirms that the temperature at the surface of the syringe does not vary outside of the T_s ± 1 °C over the experimental time course. Thermocouple measurements along the spinline coordinate z (Fig. 5(c)), where z = 0 mm is considered a measurement point under the tip (T_t = 40 ± 5 °C) and z = 30 mm is considered a measurement point on the surface of the collector plate (T_c = 30 ± 5 °C), demonstrate the presence of an exponentially decaying temperature profile (Fig. 5(d)). Due to the high thermal conductivity of glass and the small volume of polymer melt hosted in the syringe barrel, it is assumed that the temperature of the polymer melt (T_o) becomes equal to T_s, and the system reaches thermal equilibrium after 1 h. The latter is also confirmed by measuring a stable spinline temperature profile regularly after the heat gun is set over the time course of 2 h. In this way, the presence of temperature gradients higher than 5 °C along the process regimes that may yield variations in the temperature-dependent polymer viscosity and thus in the flow field along the process regimes are avoided.

Although studies that have used heated air systems have reported that the temperature at the spinneret may be difficult to control accurately using this approach [25,50], the present study demonstrates that a heat-gun based system is capable of maintaining uniform heating within the material head and a spinline temperature profile, whose higher end can be set close to the onset crystallization temperature of PCL. This capability can offer an alternative way of printing aligned fibers with submicron diameter by tuning the spinline temperature so as to induce prolonged stretching, through delayed “in-flight” fiber solidification.

Prior to the printing studies, pure PCL pellets are loaded into a glass syringe (Hamilton). Then, the syringe is placed in a laboratory convective oven and heated for 24 h to remove any bubbles that may affect the process stability and downstream structural formability of the melt electrospun fibers. After assuring the homogeneity of the polymer melt, a needle tip at a prescribed nominal inner diameter (21 gauge—0.514 mm) is adapted onto the syringe. The syringe with the attached tip is then placed in the material head of the system, which is preheated at a temperature (T_surf = 77.8 °C) with the heating element. At least 1 h is given to the system prior to initiating the printing studies in order to reach thermal equilibrium.

Thin glass coverslips are taped on the grounded aluminum plate and used as collectors for all the printing tasks in this study. In this way, the structural formability in terms of diameter and quality of the MEW fibers could be characterized using bright field microscopy. An inverted motorized microscope (IX83, Olympus, Tokyo, Japan) along with image processing software (cellsens 2.11) is used to image and characterize all samples. The fiber diameter is measured directly from the acquired images at five different points along the length of each fiber for statistical significance, and an average fiber diameter along with its standard deviation is reported. The apparent pore size is determined as the average value of the circle diameters that could be fitted inside each scaffold pore. All measurements are done using a 20×objective lens with the magnification set at 12.6.

The identification of relevant dimensionless groups can be carried out using a classical dimensional analysis technique starting with process and system-specific independent parameters [51,52]. The following definitions are employed. “n” is the number of independent variables relevant to the process. “j′” is the number of base dimensions found in the n variables. “j” is the number of variables necessary to be considered simultaneously. “k” is the number of the independent Π terms that can be identified to describe the process and is equal to n − j (k = n − j).

The variables that are relevant to the printing problem are listed in Fig. 3. The total number of independent variables, n, is equal to 12. Table 2 enumerates these variables and their base dimensions, where M stands for mass (SI unit: kilogram), L stands for length (SI unit: meter), T for time (SI unit: second), Θ for temperature (SI unit: Kelvin), and A for electric current (SI unit: Ampere). This number of base dimensions is equal to 5 with j′ = {L, M, T, A, Θ}. Next, j is determined by assuming that j = j′ and scanning for j repeating variables which do not form a dimensionless product. The prescribed number of five independent variables leads to the following independent variables j={d, Q, V_p, T_t, γ}. Thus, the number of independent dimensionless Π terms that could be formed would be equal to k = n − j = 12 − 5 = 7. The following step consists of the Π_i, i = 1, 2., 7 term formation. Each term is formed by forming a power product of the j repeating variables with the additional variable.

The solution of the system (Eqs. (3)–(7)) and its subsequent substitution in Eq. (1) yields a dimensionless term Π₁ shown in Table 3. The same procedure is followed for the formation of the remaining Π_i terms shown in Table 3. Thus, the product combination of the Π_i dimensionless terms can lead to a single dimensionless number Π

The formulation and calculation of two separate terms, N₁ and N₂, in turn enable the investigation of printability when the process is performed under a stationary (U_T = 0) and a moving collector (U_T > 0), respectively. In the former case, the N₁ term is a function of the independent process parameters that govern the polymer melt jet formation in the free-flow regime. In the latter case, the N₂ term additionally accounts for the translational stage speed (U_T), a process variable that quantitatively affects the fiber topography on the receiving substrate. To be sure, the initial N₁ term is defined for the preliminary procedural step of identifying the equilibrium state conditions in the free-flow regime to ensure stable jet formation. In the absence of this preliminary step, the direct application of N₂ for a stationary collector would yield a trivial printability value of zero.

The presented set of nondimensionalized equations (Eqs. (A1)–(A7), see the Appendix) enables the identification of the important dimensionless groups that need to be tuned toward efficient printability. These dimensionless numbers are summarized in Table 1.

Initiation of the printing process requires: (a) droplet emergence, (b) successful Taylor cone formation, and (c) subsequent emergence of a charged jet, which is electrostatically drawn across the spinline coordinate in the free-flow regime. All phenomena are dependent on the relative importance of the forces applied at the polymer melt jet.

Downstream pulling forces such as the gravitational and the electrostatic Coulombic forces are related to the Bond (Bo) number and the electrostatic force parameter (E_p), respectively. Upstream resistive forces such as the viscous, the elastic, and the surface tension forces are related to the Reynolds (Re) number, the Deborah (De) number, and the Capillary (Ca) number. According to the electrospinning operating principle, Taylor cone formation occurs when the electrostatic forces overcome the capillary forces [53]. Jet initiation and the electrostatic drawing of the polymer melt jet are strongly dependent on the viscoelasticity of the polymer melt. If the gravitational forces, along with the electrostatic drawing forces caused by the accumulation of the charges at the jet-ambient air interface, overcome the viscous and elastic stresses that are applied to the polymer melt, jet initiation occurs. Thus, the proposed Printability Number should assume values within a domain defined by a set of independent material, process, and geometry-related parameters for which the printing process can be realized.

where $ηp(Tm)$ denotes the melting temperature dependency of the polymer viscosity, and the characteristic jet radius just outside the needle tip, $Ro$⁠, is assumed to be equal to the needle tip diameter.

The dynamic properties of the PCL melt are shown in Fig. 6(a). Material functions of the Giesekus model are used for nonlinear fitting of the experimental data and the determination of model-specific input parameters [55]. As shown, the values of the loss modulus, G′′, i.e., the energy dissipated as heat, are higher than the values of the storage modulus, G′, i.e., the energy stored as elastic energy, over a broad range of frequencies (see Fig. 6(a)). In the linear viscoelastic region, i.e., relatively small strains and strain rates as would be encountered at the relatively low flow rate conditions of the melt electrospinning writing process (<50 μL/h), the shear viscosity of the polymer melt can be considered to be Newtonian, i.e., the zero-shear viscosity, $ηo(Tm)$ (it is indicated in Fig. 6(b) that up to a shear rate of 10 s⁻¹ the shear viscosity of PCL is constant). In the linear viscoelastic region, the uniaxial extensional viscosity of the melt, i.e., the Trouton viscosity, is equal to three times the Newtonian (zero-shear) viscosity, $ηo(Tm)$ [56]

$NPR,1$ can be computed using Eq. (16) for the melting range of PCL (70 °C ≤ T_m ≤ 90 °C) and a prescribed set of typical process and material parameters. The values of the material parameters (summarized in Table 4) are either derived from literature or through fitting of the rheological data. In order to assure that N_PR assumes values within a valid domain, each range is determined based on previously reported studies where PCL has been successfully processed by way of MEW. Validation of the previously reported ranges is performed through preliminary experiments with the present MEW system. Thus, a range of volumetric flow rates (25 μL/h ≤ Q ≤ 50 μL/h) are applied for a 21 gauge needle tip diameter (D_t = 2·R_o), for collector distances (d) of 10 mm to 30 mm and a range of applied voltage potentials (10 kV ≤ V_p ≤ 15 kV).

The results are depicted in Fig. 7(a), where the normalized Printability Number $NPR,1*$ is plotted for the minimum and maximum volumetric flow rates (Q_min = 25 μL/h and Q_max = 50 μL/h) as a function of the normalized melting temperature of PCL (T*). The normalized $NPR,1*$ is obtained by dividing the computed N_PR,1 value with the N_PR,1 value that defines the lower end of the printability window bounded by the material's melting range for T_ref = 70 °C and Q_max = 50 μL/h. The temperature of the polymer melt inside the reservoir (T_o) is normalized with respect to the reference temperature (T_ref = 70 °C), i.e., T* = T_o/T_ref. In Fig. 7(a), T* = T_m/T_ref since T_o assumes the melt temperature value (T_m). The printability window is seen to depend significantly on the volumetric flow rate, with the smaller Q (25 μL/h) yielding significantly larger $NPR,1*$ values compared to that obtained at the larger Q (50 μL/h). This trend is consistent with recent phenomenological observations that reflect stable printing by way of MEW under low volumetric flow rates [26,59]. As T* increases within each printability window, $NPR,1*$ increases exponentially due to the Arrhenius temperature dependence of the polymer melt viscosity, implying that for higher melt temperature conditions, the material can be electrospun more efficiently. This relationship indicates that for prescribed D_t, Q, and V_p settings, melt temperature conditions approaching the higher end of the material's melting temperature range (90 °C for PCL) enable earlier droplet emergence compared to the melt temperature conditions that approach the lower end of the material's melting temperature range, due to an increased volumetric flow rate inside the needle tip.

The N_PR,1 formulation (Eq. (16)) implies that the electrical field strength (V_p/d) and the volumetric flow rate (Q) are the key independent parameters toward efficient printability (fiber mesh printing with consistent dimensional characteristics) provided that the melting and ambient conditions in the polymer melt supply regime and the temperature profile along the spinline in the free-flow regime are not significantly perturbed during each printing event. N_PR,1 scales as N_PR,1 ∼ 1/Q and N_PR,1 ∼ V_p/d. This validates the physical significance of the derived number that expresses the key combinatorial role of electrostatic, viscous, and inertia forces toward steady electrospinning conditions as previously demonstrated for solution-based electrospinning systems [54,60–62]. Furthermore, all of the dimensionless groups are a function of Q-dependent inertial terms (see Table 1). Thus, the functional relationship between N_PR,1 and each dimensionless number is computed for the prescribed Q range and three different V_p values spanning the V_p range. The results are plotted for $NPR,1*$ as a function of the Re, Ca, De, and E_p numbers (Figs. 7(b)–7(e)) revealing that upon prescribing the melting conditions, a unique printability window can be defined for each V_p setting.

The printability window defined by the V_p range and the Q range remains to be optimized with respect to efficient printability through in situ process monitoring and morphological characterization of the printed meshes as described in Sec. 3.3 via the experimentally collected data.

During the experimental procedure, the temperature at the surface of the glass syringe (T_s) is set to 78 °C. The values for the remaining independent process parameters used for the N_PR,1 computation are set to the following for the initiation of the melt electrospinning process: (a) Q is equal to 50 μL/h, (b) V_p is equal to 12.5 kV, and (c) d is equal to 20 mm.

To enable real-time monitoring of the process in the free-flow regime, a universal serial bus microscope camera is mounted on the open side of the enclosure box and focused on the needle tip. Digital photographs showing the stages of the printing initiation process over time are shown in Fig. 8. For this set of experiments, the prescribed experimental conditions reside within the printability window corresponding to a $NPR,1*$ equal to 1.34 (Fig. 7). The observation of the process during its early stages contributes to an understanding of the dynamics and the underlying physical mechanisms governing the process, as described below.

Initially, the polymer melt enters the free-flow regime owing to the electrostatic volume forces developed from the applied voltage potential at the tip and the mechanical force applied by the syringe pump (proportional to the pressure drop of the melt through the needle). After approximately 12 min, the jet is observed to exhibit an elongated shape (Fig. 8-II), indicating that the downstream electrostatic forces along with the gravitational forces are sufficiently large to overcome the upstream resistive forces (viscous, elastic, and surface tension forces). At the point where the downstream forces overcome the upstream resistive forces, the formation of the Taylor cone is observed within approximately an additional 3 min (Fig. 8-III). Immediately following the formation of the Taylor cone, a jet emerges and is electrostatically drawn between the needle tip and the collector (Fig. 8-IV).

By extending the time course for process monitoring in the free-flow regime, it is observed that for the initially prescribed combination of process parameters, excess material is generated due to the imbalance between the downstream pulling and upstream resistive forces. Such excess material that is not fully stretched periodically disturbs (every 5 min) the flow field in the spinline regime (Fig. 9-I). After the excess material enters the free-flow regime, the cone shape is initially transformed into an oblique shape, and the partially stretched melt driven by the downstream forces leads to the relocation of the Taylor cone, i.e., the cone is pushed away from the tip of the needle, and a cylindrical body of melt occupies the distance between the tip of the needle and the Taylor cone. To achieve steady-state conditions, i.e., steady-state formation of the Taylor cone at the tip of the needle and the emerging polymer melt jet, the critical process parameters, V_p and Q, need to be optimized. Such optimization specifies the relevant Printability Number $NPR,1$ at which efficient printability is achieved. This optimization step thus aims to eliminate the perturbations observed under nonequilibrium conditions. Upon optimization, an equilibrium state, i.e., state at which the downstream pulling and upstream resistive forces are balanced, is achieved.

The tuning of the critical independent parameters in order to achieve equilibrium conditions is carried out in a stepwise manner. As a first step (Fig. 9-III), the collector is moved closer to the needle tip (d = 15 mm) to increase the electrical field strength. When the tip to collector distance d ≤ 10 mm, arching occurs due to excess ionized air molecules and dry ambient conditions (humidity < 25%). At such relatively small distances (d), the arching phenomenon becomes more pronounced for applied voltages that are ≥ 15 kV. By reducing the distance (d), a higher electrical field intensity facilitates stretching of the excess material collected at the tip. However, solely reducing the distance (d) is not sufficient to eliminate the periodic perturbations. To eliminate the perturbations and achieve equilibrium conditions, a reduction in the volumetric flow rate (Q) is also required. This is suggested by the relative importance of the Q-dependent inertial forces with respect to the V_p-dependent electrostatic forces, as guided by the N_PR,1 formulation denoted in Eq. (16). The decrease of the volumetric flow rate to Q = 25 μl/h results in the formation of a Taylor cone directly below the needle tip (Fig. 9-IV). However, chaotic jet movement occurs close to the collector plate, and stable jet cannot be achieved (Fig. 10-IV-inset). In order to eliminate the instabilities and establish equilibrium state conditions, the applied voltage potential is decreased to 11.5 kV, yielding stable cone–jet formation for a period of 30 min after which the printing process could start.

The observation that the optimized Printability Number (⁠$NPR,1*$  = 2.55) in Fig. 9-IV is smaller than the nonoptimized Printability Number (⁠$NPR,1*$  = 2.78) in Fig. 9-III raises the question of self-consistency of the proposed Printability Number. This finite difference in magnitude is attributed to a phenomenon that has previously been observed in electrohydrodynamic (EHD) cone–jets with highly viscous conductive polymer solutions [63,64]. Using voltage–flow rate (V_p–Q) operating diagrams, it has been shown that for a single Q value, steady cone–jets can be achieved for V_p values lying within a range of 1 kV. Similarly, in the present study, Q is maintained at 25 μL/h in Figs. 9-IV and 9-V , and the voltage is decreased from 12.5 kV to 11.5 kV. The decrement in V_p is aimed at eliminating the chaotic jet movement close to the collector plate (see Fig. 9-IV) that is likely caused by repulsive forces between the in-flight fiber and previously charged deposited material on the collector rather than periodic flow disturbances that are preliminarily eliminated. Therefore, the small difference between the $NPR,1*$ values does not affect the self-consistent scale of the Printability Number since it is related to the process physics, as previous work has demonstrated. Thus, by tuning the $NPR,1*$ number from an initial value of 1.34 (Fig. 9-I) to 2.55 (Fig. 9-V), the observed perturbations that may affect the downstream structural formability of the printed meshes and preclude efficient printability can be eliminated. This is validated hereafter by reproducibly printing layered meshes of woven and nonwoven topographies at both the optimum and nonoptimum printability settings.

A normalized Printability Number, $NPR,2*$⁠, is obtained by dividing the computed N_PR,2 value (based on Eq. (17)) with the N_PR,1 value (based on Eq. (16)) that defines the lower end of the printability window bounded by the material's melting range for T_ref = 70 °C and Q_max = 50 μL/h. The normalized Printability Number $NPR,2*$ is computed for each fiber pattern (Fig. 10) where optimum printability is achieved when U_T is tuned to its critical value (U_T = U_CR).

Layered meshes with woven and nonwoven architectures are fabricated using various $NPR,2*$ settings. Woven meshes with “0–90 deg” and “0–45–135–90 deg” pore architectures are fabricated using optimized and nonoptimized $NPR,2*$ settings. When $NPR,2*$ is not optimized, irregular structures are observed. This is shown in Fig. 11(a), which is obtained at an $NPR,2*$  = 31.9, where U_T = 25 mm/s < U_CR. When $NPR,2*$ is increased to 57.63 by independently tuning the stage speed (U_T = 85 mm/s ≥ U_CR) while neglecting equilibrium conditions in the free-flow regime, aligned structures with variable average fiber diameters (D_f = 27 ± 14 μm) are observed, as shown in Fig. 11(b). On the other hand, precise printing of mesh architectures composed of well-aligned fibers with uniform average diameters (D_f = 23 ± 3.7 μm) can be produced for an optimal Printability Number of $NPR,2*$  = 106. The produced fibers at this optimal printability setting are shown in Figs. 11(c) and 11(d) bearing the hallmarks of equilibrium state conditions in tandem with appropriate tuning of U_T at its critical value.

The proposed Printability Number N_PR assumes the form of N_PR,1 (see Eq. (12)) for a stationary collector (U_T = 0 m/s) and the form of N_PR,2 (see Eq. (17)) for a moving collector (U_T > 0 m/s). N_PR,1 is proposed in order to guide the user to achieve equilibrium state conditions in the free-flow regime. N_PR,2 is proposed in order to guide the user to produce well-defined microarchitectures (efficient printability). The latter Printability Number requires the process to operate at both equilibrium state conditions as well as at critical stage speed conditions. The ensuing discussion highlights the science and significance, along with the potential downstream biological applications enabled by proposed Printability Number N_PR,2 where the subscript 2 is omitted to denote the proposed Printability Number as N_PR. The values of N_PR are reported in its normalized form as $NPR*$⁠.

The science of the proposed number is based on the dynamics of the electrohydrodynamic phenomenon that governs the stability of the cone–jet formation across the spinline regime and its subsequent drawing on a collector using mechanical forces. The independent variables are tunable process and material parameters, whose magnitude defines the magnitude of key dimensionless numbers. These numbers in turn express the relative importance of all the forces acting on the polymer melt during the different process regimes. A discussion regarding the robustness of the Printability Number may raise questions as to the effectiveness of the proposed number under various settings for the most important independent parameters. For example, if V_p and Q are varied for the same $NPR*$ number, it is important to consider whether the same printing behavior will be observed. Based on the methodology followed for the optimization of the proposed Printability Number, the same printing behavior will be observed only if the alternate V_p and Q setting promotes equilibrium state conditions in the spinline of the free-flow regime (see Sec. 3.3). The optimum $NPR*$ equal to 106 suggests that for a new setting of V_p = 11 kV and Q = 15 μL/h, the same printing behavior should be observed by setting U_T = 51 mm/s. Indeed, multilayered meshes (ten layers, 100 fibers/layer) with consistent fiber diameter equal to 10 μm are printed, as shown in the micrographs of Fig. 12. This demonstrates the effectiveness of the $NPR*$ number directly enabling the printing of high-fidelity meshes without having to determine the U_CR at which aligned fibers are obtained. However, it is important to mention that the derived Printability Number neglects the charge distribution phenomena during MEW and particularly the effect of residual charge entrapped inside the printed polymer fibers. These phenomena warrant systematic investigation, since it has been shown to compromise the interfiber distance resolution of the process as well as the fidelity of the printed constructs upon subsequent layering [26,37]. This phenomenon is depicted with red insets at the indicated areas of the printed mesh shown in Fig. 12. Thus, the proposed $NPR*$ can allow efficient printability at an interfiber distance sufficiently large so that such charge interferences during fiber deposition are eliminated. A modified version of the MEW system that allows in situ charge monitoring is currently in development for the design of experiments that correlate entrapped charge with the interfiber distance resolution. The empirical observations on residual charge-dependent fiber deposition would introduce added complexity and predictive power to a modified N_PR number. The proposed method may be used in combination with numerical model prediction of MEW fiber diameters in the context of a semi-empirical dimensional analysis approach [54,60]. This approach could be used to obtain a single correlation of fiber diameter and final scaffold morphology with universal applicability across different polymer melt systems, since the proposed Printability Number formulation accounts for the Arrhenius temperature dependence of key rheological material properties. As the material library of MEW expands beyond PCL, future demonstration of printed structures by way of MEW processing of alternative polymer melts should be accompanied by detailed rheological characterization studies and in situ monitoring of the process dynamics in combination with the charge distribution along the process regimes.

The significance of the proposed number lies in the intuitive nature of the framework to guide and model the quality of a complex multiphysical process outcome. While equation-based models have been useful in predicting the fiber diameter profile of electrohydrodynamic polymer solution jets in steady-state, they cannot provide information related to the quality of the final process outcomes. In the context of melt electrospinning writing, this challenge is augmented since the characteristic quality of the outcome with respect to the uniformity of the geometrical features, herein defined as “printability,” is governed by the coupled dynamics of electrospinning and 3D printing. As the experimental results demonstrate, the proposed number captures both key phenomena, thereby enabling the fabrication of high-quality 3D microscale fibrous meshes while averting the need for trial and error methods to produce a prescribed construct. This represents an important step toward the development and future commercialization of melt electrospinning writing as a reliable high-resolution additive biomanufacturing process capable of fabricating hierarchical tissue constructs with minimal geometric deviations between the intended design and the produced manufacturing outcome [65,66].

Currently, the proposed number can be regarded as a metric for an emerging additive manufacturing process capable of producing mesh structures with cellular-relevant dimensional features (1–50 μm). Specifically, it addresses the key manufacturing challenge of reproducibility in producing 3D biomaterial scaffolds that can serve as physiologically relevant 3D cell culture models amenable for advanced tissue engineering applications and fundamental cell biology studies. Compared to random nanofiber meshes, the effect of dimensional features within the aforementioned range on cellular function has not been systematically explored. The reproducible fabrication of structurally homogeneous micron-scale meshes via MEW warrants optimized tuning of numerous process and material variables. To this end, the proposed reductionist approach serves as a tool whereby different Printability Numbers map to specific cellular responses.

The precise printing of fibrous mesh structures via melt electrospinning writing poses a significant manufacturing challenge. The system dynamics for the coupled electrospinning-3D printing process is poorly understood, resulting in the application of labor intensive trial and error methods. To methodically address these dynamics, a Printability Number $NPR*$ is defined on the basis of dimensional analysis and correlated with the dimensionless parameters associated with the conservation equations. A tuning procedure is established to determine an optimum Printability Number. At this setting, the downstream pulling and upstream resistive forces are shown to be balanced in combination with the tuning of the translational stage speed at its critical value in order to achieve a steady equilibrium printing state. The procedure is tested using a custom melt electrospinning writing system that employs a heat-gun based thermal management system. It is demonstrated that fibrous meshes with well-defined pore architectures can be printed under optimum $NPR*$ settings. Future work will demonstrate the biological relevance of the proposed Printability Number by testing cellular responses on scaffolds printed at different Printability Numbers and incorporating the charge distribution effect with advanced thermal and electrostatics control along the different process regimes.

We would like to thank Biruk Gebre (Research Engineer, Design & Manufacturing Institute, Stevens Institute of Technology), Milan Simonovic (Laboratory Supervisor, Stevens Institute of Technology), and Marshall Reid (Technician, Stevens Institute of Technology) for their valuable advice and help through the design process of the melt electrospinning system. The work presented in this paper was supported by the National Science Foundation under Award No. CMMI-MME-1554150. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the National Science Foundation.

$a$ =

mobility factor

$d$ =

distance between the spinneret and the collector (m)

$E$ =

electric field—z component (V/m)

$En$ =

electric field—normal on the jet surface (V/m)

$Et$ =

electric field—tangential on the jet surface (V/m)

$I$ =

jet current (A)

$K$ =

electrical conductivity of the polymer melt (S/m)

$p$ =

pressure (N/m²)

$Q$ =

volumetric flow rate (m³/s)

$R$ =

jet radius (m)

$rn$ =

viscosity ratio

$Ro$ =

jet radius at the origin (right outside of the spinneret)—radius of the spinneret (m)

$Tc$ =

temperature right above the collector plate (z = d mm) (°C)

$Tm$ =

melting temperature of PCL (°C)

$To$ =

temperature of the polymer melt inside the melt reservoir (°C)

$Ts$ =

temperature at the surface of the melt reservoir (°C)

$Tt$ =

temperature right under the tip (z = 0 mm) (°C)

$v$ =

jet velocity (m/s)

$Vp$ =

applied potential difference between the spinneret and the collector (V)

$z$ =

spinline coordinate (mm)

$γ$ =

surface tension of the polymer melt (N/m)

$ε$ =

absolute permittivity of the polymer melt (F/m)

$εo$ =

vacuum permittivity (F/m)

$εr$ =

relative permittivity of polymer melt

$ηo$ =

viscosity of the polymer at zero shear rate (Pa·s)

$ηp$ =

viscosity of the solution due to the polymer melt (Pa·s)

$λ$ =

relaxation time (s)

$ρ$ =

density of the polymer melt (kg/m³)

$σ$ =

surface charge density (C/m²)

$τp,rr$ =

normal stress—radial direction component (N/m²)

$τp,zz$ =

normal stress—axial direction component (N/m²)

A thin filament approximation is used, and by focusing on a small part of the melt electrospun stable jet region, a one-dimensional momentum balance is made by considering the various forces affecting the jet profile. The jet is subjected to: (a) Coulombic electrostatic, viscous, elastic, surface tension, and gravitational forces. Assuming axisymmetry along the path from the tip of the spinneret up to the surface of the collector (at distance, d) and using the characteristic quantities defined in Table 1, the dynamics of the melt electrospun jet can be modeled using the following system of nondimensional equations, where $R$ is the jet radius divided by the characteristic jet radius $Ro$ just outside of the needle tip, $v$ is the jet velocity divided by the characteristic velocity $vo$⁠, $R$ is the jet radius, and the prime indicates derivatives with respect to the spinline coordinate z:

1. (1)
   Conservation of mass—Continuity:

   $R2v=1$

   (A1)
2. (2)
   Conservation of momentum:

   $Revv′=Bo+3(1−rn)(R2v′)′R2+Tp′R2+CaR′R2+Ep(σσ′+βEE′+2EσR)$

   (A2)

   where Re, Bo, Ca, and E_p are defined in Table 1.
3. (3)
   Conservation of charge:

   $σ=R$

   (A3)
4. (4)
   Electric field:

   $Et=1(1+2z−z2/χ)1+(R′)2$

   (A4)

The viscous polymer stress τ_p denotes the elastic nature of the material due to normal stresses that arise during its deformation, and the strain rate tensor $γ˙$ is given by the sum of the velocity gradient and its reciprocal [55]. The input parameters of the Giesekus model (Eq. (1)) that are determined by fitting the experimental raw data on the basis of the corresponding material functions for each type of tested viscometric flow, are the following: $np$ represents the polymer viscosity, $λ$ the relaxation time, and $α$ the mobility factor, which is a parameter related to the anisotropic Brownian motion and/or hydrodynamic drag on the constituent polymer molecules [55].