Numerical simulation of electrospinning process in commercial and in-house software PAK

The aim of this research was to investigate if it is possible to implicitly determine the homogeneity of the obtained electrospun fibers based on jet shape during electrospinning. Experiments were performed with 10 wt% PVA solution, and four variations in process parameters were investigated in order to examine their effect on fiber structure. Data obtained during experiments were used as input for computational simulation. The simulation results, both using commercial ANSYS and in-house software PAK, show good agreement with experiments in terms of outcome-no fiber differences in experiments were present when different voltage pairs were used, and similar jet shapes were obtained during simulations. Shapes of the electric field potential for all the used voltage pairs were very similar, due to the uniformity of the field, which is in agreement with the experiment, as no differences in fiber structure are observed in these cases. This confirms the hypothesis that based on jet shape during electrospinning, it is possible to implicitly determine the homogeneity of the obtained electrospun fibers. Differences that may occur between experiments and simulation can be a result of simplifications in simulations, influence of uniform and non-uniform electric field etc. This kind of two-phase simulation could be useful in reducing the trial-and-error approach and maintenance costs in electrospinning experiments.


Introduction
Electrospinning is widely known technique, which is used to produce micro-and nanofibers through the process of polymer elongation due to the electrostatic field (Spivak et al 2000, Ghaly 2015, Liverani and Boccaccini 2016. It is a multi-physics and multi-phase process that involves several phenomena-mass and heat diffusion and transfer, evaporation, electrohydrodynamics etc. (Xu et al 2014). Generally, the parameters that govern the whole electrospinning process can be divided into three main groups: solution parameters (polymer concentration, solvent system, conductivity, etc.), process parameters (flow rate, applied voltage, distance between needle tip and grounded collector and collector design), and the ambient parameters (temperature and relative humidity). These parameters control the uniformity, continuity and range of fiber diameter and fabrication of aligned fibers (Rajput 2012, Liverani andBoccaccini 2016).
Industrial applications of electrospinning are still relatively rare, mainly because of the unresolved problem of very low fiber throughput for existing devices Zussman 2004, Kowalewski et al 2005). Reports show that 70% of the electrospun polymer fibers are used for tissue engineering application such as tissue engineering scaffolds, wound dressing etc (Rafiei et al 2013, Ghaly 2015, Hu 2015. Many researches focused on applications of electrospun fibers in bone engineering in order to develop bone substitute or implant (Salgado et al 2004), cartilage tissue engineering (Hutmacher 2000), nerve tissue engineering (Groeber et al 2011) and skin tissue engineering (Ananta et al 2011).
Many experimental studies have been conducted to fully understand the electrospinning process and some mathematical models have been developed to describe the mechanism of the process (Sanders et al 2003, Wu et al 2011, Xu et al 2011. As far as modelling is concerned, Spivak et al (2000) modeled the straight electrified jet, Fridrikh et al (2003) analyzed the terminal diameter of a thinning electrified jet at the last stage of bending, while Shin et al (2001) developed the models of the dynamics of bending jets in electrospinning. He and Wan (2004) examined the allometric relationship between jet radius and the axial distance from spinnerette. Special attention is payed to the electric field distribution as it has great influence on the fiber diameter and behavior during electrospinning (Zheng et al 2013, Jiang et al 2016. They concluded that more uniform electric field distribution produces finer and more homogeneous fibers. Wang et al (2013) proposed electrical force formula after simulation of the distribution of electrostatic field in electrospinning using Finite Element Method (FEM). Weimin et al (2013) also used FEM method for simulation of electric field intensity and electrostatic force in electrospinning system.
The analysis of experimental data has been done as an approach towards developing the practical tools for finding optimal parameters in the electrospinning process (Kowalewski et al 2005). The numerical simulation of electrospinning process, which includes the shape of the jet travelling from the needle to the collector will be useful in terms of optimization of the electrospinning process, reducing the trial-and-error approach, having beneficial effect in decreasing the amount of chemicals, amount of waste management and equipment maintenance costs. To the best of authors' knowledge, all topics related to numerical simulations in literature include single-phase flow (simulated only one fluid-electrospinning polymer), and are not able to offer indepth insight into physical understanding of many complex electrospinning phenomena, which cannot be fully explained experimentally. Those approaches are mostly discrete modeling with partial differential equations (PDE) or finite element method of either electric field of Taylor cone alone, but no research has used computational approach that included two-phase flow (air and electrospinning polymer are the two fluid phases) and analyzed the whole jet in electrospinning. In that sense, this paper is the first to introduce the coupling between the electrospinning polymer and surrounding air in simulations. Therefore, this paper proposes the use of Finite Element Method (FEM) in examining the electrospinning jet trajectory. In addition, not only the commercial software is used, but also the governing equations are implemented in the in-house software PAK, which enables us to further expand the model adding more phenomena (e.g. evaporation, solidification etc) in future researches. The extension of implemented models would not be possible by using just commercial software. With all this said, this paper examines the two-phase flow in the process of electrospinning by computational approaches, in order to gain significant insight and a better understanding of the process.

Experimental setup
Electrospinning was performed using a commercially available EC-CLI device (IME Technologies, Geldrop, The Netherlands). The device is equipped with a climate control that allows the setting and control of the temperature and relative humidity (RH) inside the electrospinning chamber ( figure 1(A)). The device is also equipped with a gas shield module and the nitrogen flux was used in experiments PVA_2 to PVA_4. The distance between the nozzle and the collector was 14 cm and drum collector (with 9 cm diameter) was rotating at constant speed of 500 rpm. Needle diameter (G) used was 21 G, with inner diameter of 0.51 mm and outer diameter of 0.82 mm. Figure 1 with details B-D depicts needle with solvent drop during experiment ( figure 1(B)), gap between nozzle and drum collector before experiment (figure 1(C)) and fibers obtained after experiment ( figure 1(D)). These dimensions, along with needle dimensions provide data for geometrical model necessary for computational simulation. The duration of the electrospinning process was fixed for all samples at 5 min.
To evaluate the influence of different electrospinning parameters on fibers morphology, four experimental conditions of electrospinning a polyvinyl alcohol (PVA) dissolved in deionized water were performed. In particular, the solution concentration was fixed at 10 wt% of PVA (average molecular weight 72 000 Da, AppliChem GmbH, Germany), while voltage and presence of gas shield were of primary concern to investigate the influence on fibers morphology. Parameters used in each experiment are given in table 1.
Additionally, temperature and relative humidity (RH) were kept constant and were fixed at 25°C and 25% RH respectively. In experiments when the gas shield was included, nitrogen flux was 8 ml min −1 . Moving part connected to the needle had the moving distance from −30 mm to 50 mm with speed 10 mm s −1 .

Modelling the electrospinning process
Computational domain-geometry and mesh: Based on obtained dimensions, a simplified geometrical twodimensional model was created. At the beginning of simulation, nozzle was loaded with polymer solution, whilst the area between nozzle and drum collector is loaded only with air. Air domain was with a diameter D=90 mm and the tip-to-collection distance is 14 cm. The syringe needle of internal diameter is 0.51 mm.
Area of interest was stability and instability region, which was assumed to grow starting from nozzle up to the collector, and was assigned smaller mesh with triangular boundaries (figure 2). The problem is threedimensional; however, two-dimensional models have the advantage of simplicity and smaller number of elements and therefore smaller computational necessity. Two-dimensional mesh consisted of 70 896 elements and 71 150 nodes. Element quality analysis showed satisfactory results with average value of 0.887 with standard deviation of 0.076.
Calculations have been done both using in-house software PAK, built at the University of Kragujevac, where electric field equations were added for the purposes of this study, and commercial software ANSYS Academic Research Mechanical, Release 16.1, ANSYS, Inc.; Fluent with included add-on MHD module (Fluent 2012).
Models: In order to simulate the interaction between the air and the polymer in ANSYS, three models were adopted: (1) Magnetohydrodynamics (MHD) model-used to simulate the interaction between flow field and electric field (2) Turbulent model (k-ε)-used to describe turbulent flows in electric fields (3) Volume of Fluid (VOF) model-used to determine the interface between two fluids  Electromagnetic fields are generally described by Maxwell's equations given by: where q e is the electric charge, J is the current, E is the electric field, B is the magnetic induction, H is the magnetic field, D is the atomic displacement vector, and c is the velocity of light in a vacuum c 2.998925 0.00006 10 m sec . 8 = ( ) Introducing additional source terms to the fluid momentum equation and energy equation achieves the MHD coupling. Lorentz force is the additional source term in the fluid momentum equation, given by: Joule heating rate is the additional source term for the energy equation, given by: where s is the electric conductivity of the media. In general, electric field E  is given by: where j and A  are the scalar potential and the vector potential, respectively. For a static field and assuming b B , 0      Ohm's law can be adapted to the form: The boundary condition for the electric potential j is given by: for an insulating boundary, where n  is the unit vector normal to the boundary, and 0 j j = for a conducting boundary, where 0 j is the specified potential at the boundary.
The VOF method is based on fraction function , q a which is defined as the integral of fluid's characteristic function in the control volume (volume of a computational grid cell). When the cell is empty, there is no trace of fluid inside and 0; q a = when the cell is full, 1 q a = and when the interphase interface cuts the cell, then 0 1. Because PAK software is based on the Finite Element Method (FEM), mixed formulation for velocity, pressure and electric potential is: where r is the mass density kg m , j is the electric potential, u is the velocity, k is the constitutive constant and q is the heat flux vector W m .
which are coupled with Navier-Stokes equations. Additionally, an equation for the electrostatic force, which influences the flow field, that depends on the force in y axis F , B ( ) coefficient constants (β and 0 q ), should be incorporated. H is the interpolation matrix, and ρ is the density of the material.
The implemented algorithm in PAK software is a specific algorithm of jet movement using Gauss points of 2D axisymmetric elements.
Boundary conditions: For all walls, no-slip boundary conditions are used. As described previously, three models were used to describe the electrospinning process. Magnetohydrodynamics (MHD) model is initialized as electric potential model, with prescribed values on inlet and outlet. Volume of Fluid (VOF) model is initialized such that the needle is fully loaded with only polymer 1,

Electrospun fiber diameter calculation
The morphology of all obtained fibers was investigated using a scanning electron microscope (SEM, Auriga 0750, Zeiss, Germany) after sputtering the electrospun samples with gold using a sputter coater (Q150T, Quorum Technologies, Darmstadt, Germany). SEM analysis of the obtained electrospun fibers is reported in figure 3. Fiber average diameters (table 3)  These results show homogenous structure with average diameters that do not differ, despite the change in process parameters. When the gas shield module was used, there was no difference in fibers size and structure in neither of the samples. It should be also stated that different positive and negative charges (and same voltage difference) did not influence fibers structures (nor difference was present when 15 kV and 0 kV and 13 kV and −2 kV potentials were used).
In PVA_4, some thicker fibers and defects were noticed (figure 4), which can be explained by using higher voltage (20 kV instead of 15 kV). Overall, these defects were only occasional and did not influence general homogenous structure.

Simulation results
The fraction of each phase is shown in figure 5. The results show that the jet starts with formation of Taylor cone (figure 5(A)), which is further broken and the jet is elongated into stability region ( figure 5(B)). It can be seen that now the mixture of air and polymer takes up most of the surrounding area of the jet, while the solely PVA jet is of micro order of magnitude. Up to the last step, which corresponded to 0.1 s of the simulation, the jet is further elongated to the collector ( figure 5(C)). This means that the whole process is very fast which is in accordance to the literature. Previous researches studying the jet profile showed that the stability of the jet and the cone-like surface of the jet have significant effects on fiber qualities (Theron et al 2004). It is reported that bending instability develops due to the mutually repulsive forces resulting from the electric charges of the jets , Yarin et al 2001. However, bending instability was not always observed (Theron et al 2004). The reason for disappearance of the bending instability are higher electric field strengths and higher flow rates, resulting in a straight jet over the whole electrode-to-ground distance. Additionally, the bending instability appeared more readily in solutions with lower molecular weights (Theron et al 2004). Other factors such as electric conductivity, magnetic permeability, viscosity etc could also contribute to the shape of the jet. Since these parameters were not able to be experimentally determined for the used solution, used values were adopted from literature as stated in table 2. Another possible explanation for this kind of jet behavior lies in simplification of simulated behavior of the needle, which was modeled stationary during simulation, whereas in experiment the sub-setup with needle was moving along 80 mm long path with constant velocity. This kind of movement causes non-uniform electric filed with maximum value at the droplet/Taylor cone tip (Theron et al 2004). In contrast to the experiment, and due to high computational necessity for these simulations, stationary needle causes uniform field that contributes to having only stable, non-bending jet.  It can be seen that the electrospinning has two main stages-the initial stage where higher velocities are present due to the polymer passing through a small surface orifice and further down the shown path, where almost constant velocity is obtained ( figure 6). Velocities in the upper part, obtained by simulation in ANSYS, have the maximal value 2.67 m s −1 . These results are in correspondence with experimental literature results obtained in (Kowalewski et al 2005), which show the velocity around 2.5 m s −1 at 10 kV, in the upper parts.
Comparing the results of the velocity during the process obtained by the ANSYS and PAK software, we can see that maximal velocities, obtained in the upper parts of the jet, are 2.96 m s −1 . Here is only right half of the model shown, due to symmetry, where the y axis is the symmetry axis. Figure 7(A) shows the case at the beginning of the simulation, where Taylor cone is formed, figure 7(B) shows the 60th iteration, while the figure 7(C) shows the last iteration when the jet has reached the cathode.
Additionally, when non-optimal parameters (non-optimal in the sense that experiments gave no fibers collected on the collector or the fibers were very little collected) were used in the simulation, it was also shown that electrospinning jet was broken and back-flow was present during simulation, meaning that no jet reached the collector and no fibers would be found there. The experiments confirmed the simulation results, no fibers were collected on the collector using these conditions. One example of such a case is shown in figure 8. Figure 8(A) shows phase fraction and it is visible that there are breaks in the fiber (circled in red). This results in turbulent reversed flow, as such small elements are carried away by the electric field (figures 8(B) and (C)).
In comparison to the ANSYS results, the results obtained with PAK software also show backflow ( figure 10(A)), with a detail on the left side ( figure 9(B)), which is the case when no jet reaches the cathode and therefore no fibers are collected.
Electric potential contours after at the moment that the jet reached the collector show regular decrease of electric potential from the needle to the collector ( figure 10(A)). Since the positive voltage is applied at the inlet, it is expected that the field shape at the beginning of the air region is a half circle (figure 10(C)), which gradually decreases to the final destination and when it reaches grounded collector, it has a value of 0 kV. The shape of the field throughout the middle part of the area of interest and in the middle of the collector is affected by the jet, which causes disturbances in regular shape ( figure 10(B)).
In comparison to the used 15 kV on needle and grounded collector, two additional simulations included applied 13 kV on needle and −2 kV on collector, as well as 20 kV on needle and 0 kV on collector. As far as the consistency of the experimental and simulation results, it was observed that shapes of the electric field potential for all used voltage pairs were very similar, which is in agreement with the experiment, as no differences in fiber diameters are observed in any cases. These mentioned three voltage pairs (15 kV and 0 kV, 13 kV and −2 kV and 20 kV and 0 kV) cover three boundary cases in the electrical sense-first two differ only by upper and lower value, while potential (the difference) between the cathode and anode is the same (15 kV). This means that the electric field shape in theory should be the same and only the maximal and minimal value should be different. The simulation confirmed this, as shapes of the fields were the same. Third voltage pair was different from the first two, as the potential was 20 kV. The results in these three cases showed little difference again in the shape of the field, meaning the shape of the jet should be the same and in the experimental sense that would mean collected fibers of similar diameter. The experiment confirmed this hypothesis, as no fiber diameter differences were observed in any case of the voltage variations.
Results of electric filed distribution for all three cases are given in figures 11(A)-(C) respectively. There are no shape differences in the electric field with 15 kV and 0 kV pair, 13 kV and −2 kV pair and 20 kV and 0 kV pair, which is in agreement with the experiment, as no differences in fiber structures are observed in these cases.
It was also obtained that the diameter of the jet decreases with the fiber travelling towards collector, phase fraction reduces along the jet. In addition, only slight differences were observed between three simulated situations -maximum phase fraction in first case voltage pair was 0.99, while in the second and third case, maximum phase fraction was 0.97. Measurements on all three simulations showed 0.02 mm-0.05 mm diameter of the jet at about 2 mm from the orifice of the needle, which is in correspondence with available literature   (Theron et al 2004). All this indicates that based on jet shape during electrospinning, it is possible to implicitly determine the homogeneity of the obtained electrospun fibers.
Limitations of this study include assumed material properties of the PVA for viscosity, electric conductivity and magnetic permeability, as well as for the surface tension that were adopted from literature, as no  experimental values could be obtained at this point. It could be said that the limitation of the study is that only voltage pairs and gas shield were investigated. This is because we wanted to first validate the proposed model and see if the two-phase model could be used to implicitly determine if the shape of the electrospinning jet is a predictor of the quality of the fibers collected on the collector. Generally, no models like this existed previously in the literature, and the aim was to introduce and validate the proposed model, with the statement of limitations that not all the influencing parameters were investigated. However, we do plan to examine many more influencing parameters in future research. It could be discussed that some parameters like polymer concentration, conductivity or flow rate could be easily implemented in the system, as those are constants in the model, which will be for sure added in the future studies. Other parameters like distance between needle tip and grounded collector change the geometrical part of the model, which again was not of interest in the initial model, but will be observed in future. Investigations of some parameters like temperature and relative humidity was not possible at the moment as those parameters are at a higher level of complexity and no literature on simulations has included these phenomena in models yet. In future, we will take into account solvent evaporation and solidification of the jet, as well as potential disturbances caused by air movement inside the chamber, which play pivotal roles in determining the jet shape and indirectly internal fiber morphology. Future research will also consider the effect of gas shield in the form of higher electric field density as boundary condition that reduces surface charge density of the jet and increases stable jet length (ElectrospinTech 2018). In addition, high-speed digital imaging systems that are equipped with zoom lenses are reported to be capable of capturing both the stability an instability parts of the jet (Theron et al 2004, Kowalewski et al 2005. The same systems can be applied in order to compare the behaviour of the jet during experiments and during simulation, under same specific conditions, which will be the theme for future research.

Conclusions
Previous numerical studies have dealt only with one-phase flows or had too many simplifications when it comes to electrospinning process. Parameters used in experimental studies are obtained by collecting data from many experiments, which results in an increase in material consumption, and can lead to higher equipment maintenance costs in long term. This was all motivation for research, which would apply computational approaches, using both ANSYS-Fluent and in-house software PAK, to simulate the electrospinning jet from the syringe nozzle to the collector. The shape of the jet in stability and instability region can be a factor of indication whether the chosen parameters would result in fibers with good quality.
The results of the simulation with both ANSYS and PAK software show good agreement with experiments. In the sense of simulation, no difference was obtained in the jet shape using aforementioned three voltage pairs, and in the sense of corresponding experimental results, there was no statistically significant differences in fibers collected on the collector. Furthermore, with conditions that gave no fibers collected on the collectors in experiments, simulations also showed that jet did not reach the collector.This confirms previous researches that claimed it can be implicitly determined if the electrospinning fiber quality will be satisfying, just based on jet shape during electrospinning. Differences that may occur between experiments and simulation can be a result of simplifications in simulations; influence of uniform/non-uniform electric field, as well as adopted parameters that were used based on literature values. It can be concluded that the computational approaches are a powerful research mechanism of electrospinning process and can be helpful in choosing optimized parameters for electrospinning.