Numerical simulation of electrified jets: An application to electrospinning D. Borzacchiello, S. Vermiglio, F. Chinesta, S. Nabat, and K. Lafdi Citation: AIP Conference Proceedings 1769, 050004 (2016); View online: https://doi.org/10.1063/1.4963432 View Table of Contents: http://aip.scitation.org/toc/apc/1769/1 Published by the American Institute of Physics Articles you may be interested in Electrospinning and electrically forced jets. I. Stability theory Physics of Fluids 13, 2201 (2001); 10.1063/1.1383791 Taylor cone and jetting from liquid droplets in electrospinning of nanofibers Journal of Applied Physics 90, 4836 (2001); 10.1063/1.1408260 Bending instability of electrically charged liquid jets of polymer solutions in electrospinning Journal of Applied Physics 87, 4531 (2000); 10.1063/1.373532 The stretching of an electrified non-Newtonian jet: A model for electrospinning Physics of Fluids 14, 3912 (2002); 10.1063/1.1510664 Electrospinning and electrically forced jets. II. Applications Physics of Fluids 13, 2221 (2001); 10.1063/1.1384013 Numerical Simulation of Electriﬁed Jets: An Application to Electrospinning D. Borzacchiello1,a) , S. Vermiglio1 , F. Chinesta1 , S. Nabat2 and K. Lafdi2 2 1 Ecole Centrale de Nantes, GeM UMR 6183, 1 rue de la Noe, F-44300 Nantes, France Chemical and Materials Engineering, University of Dayton,300 College Park, Dayton, Ohio, USA a) Corresponding author: domenico.borzacchiello@ec-nantes.fr Abstract. This paper concerns the numerical simulation of electriﬁed jets with application to the electrospinning process for the fabrication of ﬁbers with controllable size, diameter, and cross section shape. Most numerical models used to simulate electrospinning rely on the Upper Convected Maxwell model (UCM) which is ﬁt to model polymer melts. However, in most electrospinning processes the ﬂuid is a polymer solution with a Newtonian solvent that evaporates after the ﬁber is deposited on the collector. In this work we propose to describe the ﬂuid rheology using Giesekus model, which predicts the properties of polymer solutions more accurately, and show the impact of the rheological model on the prediction of the ﬁber radius and size. INTRODUCTION The physics of electrically driven jets is extensively investigated in the seminal work of Taylor [1]. Experimental evidence shows that rectilinear jets are unstable under small bending disturbances developing a so called whipping instability. The reason for the observed instability is attributed to the repulsive Coulombic forces between the electrical charges in the ﬂuid and can be explained by Earnshaw’s theorem [2]. In an ideally rectilinear electriﬁed jet, the net force exerted on a single charge by the remaining ﬂuid acts along the jet axial direction. The jet is constantly elongated due to the repulsive interaction between the charges. However, if the axial symmetry is perturbed due to small amplitude disturbances, the net electrostatic force has a lateral component. If Coulomb interactions are strong enough to counteract the stabilizing eﬀect of the surface tension, the initially small disturbances are magniﬁed and the instability develops. A more detailed description of the instability mechanism can be found in [3]. The main interest in this phenomenon is the possibility to use whipping instabilities in electriﬁed jets to produce a variety of continuous polymeric ﬁbers by the process of electrospinning. Fibers are formed by the creation and elongation of an electriﬁed ﬂuid jet. The advantage is that, by using an adequate voltage diﬀerence of tens of kV, the jet can be elongated to degree that the resulting ﬁbers have a thinner diameter (from nanometer to micrometer) than ﬁbers produced with other spinning techniques. The basic apparatus, schematically depicted in ﬁgure 1, consists of a high voltage power supply, a spinneret and a grounded collecting plate. A electrically charged ﬂuid, often a polymer solution, is accelerated trough a straight capillary or a tapering cone and is cast towards the collector of opposite polarity. Diﬀerent types of collectors can be employed, including conductive papers, conductive cloths, wire meshes, parallel or grided bars, rotating rods and rotating wheels. Other setup conﬁgurations are reported in [4]. As the free jet is stretched by the electrical ﬁeld, the elongation results into a diameter reduction, therefore the ﬁnal diameter can be predicted and controlled by adjusting the distance traveled by the jet, the voltage source, the spinneret outlet section diameter and the ﬂow rate. Additional constituents can be included in the ﬂuid like chemical reagents, other polymers or dispersed ﬁllers, therefore allowing for a ﬁne tuning not only of the ﬁber geometry but also of the mechanical, thermal and electrical properties. Electrospun ﬁbers are used in various ﬁeld of application like tissue engineering, ﬁltration, protective clothing, drug delivery, energy generation and cosmetics [4, 5, 6]. Several mathematical models have been proposed to describe bending instabilities in electriﬁed jets in an attempt to predict the ﬁnal ﬁbers geometry through numerical simulation and select optimal process parameters. The physics of electriﬁed jets is governed by the following equations : ESAFORM 2016 AIP Conf. Proc. 1769, 050004-1–050004-6; doi: 10.1063/1.4963432 Published by AIP Publishing. 978-0-7354-1427-3/$30.00 050004-1 HC End of capillary Electriﬁed jet Collector R V FIGURE 1. Schematics of the electrospinning process. • • • • Conservation of mass Conservation of electric charge Balance of linear momentum Evolution of the stress tensor Numerical diﬃculties associated to the solution of these equations are mainly due to nonlinear coupling, moving meshes and the fact that the characteristic ﬁlament length is by several orders greater than the diameter, which would either require extremely stretched elements or very ﬁne meshes. Therefore, some simplifying assumptions are in order to eﬃciently treat the problem. Several simpliﬁed models have been proposed for both rectilinear [7, 8, 6, 9, 10] and whipping jets [11, 12, 13] for both linear stabilities studies and direct numerical simulations. In particular, the models that deals with direct simulation of the bending instabilities, like [11] are classically based on the hypothesis of one dimensional elongational ﬂow along the jet axis, while the variations of the ﬂow ﬁeld in the cross section are neglected. Further simpliﬁcations are made considering the jet as formed by discrete rectilinear cylinders with constant radius with both mass and electric charge lumped at the two ends. In this way the continuous equations are reduced to a set of discrete equations for beads connected through spring-dashpot elements modeling the viscoelastic behavior of the ﬂuid. In the majority of works, viscoelasticity is taken into account with an Upper Convected Maxwell model (UCM). Polymer solutions can be modeled more accurately with Oldroyd-B model which considers the total viscoelastic stress as the sum of a contribution coming from a Maxwellian ﬂuid and a contribution coming from a Newtonian solvent. In fast deformation with strong elongational rate, the elongational viscosity of both Maxwell’s and Oldroyd-B models tends to inﬁnity. The Giesekus model is a generalization of Oldroyd-B model based on the concept of deformation dependent mobility [14]. This allows to have a bounded elongational viscosity as the elongational rate increases. In this paper we consider the bending instability of an electriﬁed jet of a Giesekus ﬂuid. First, the governing equations are reviewed and the characteristic numbers arising from dimensional analysis are presented. The adopted numerical method is explained in detail and simulation results for electriﬁed jets using both UCM and Giesekus model are reported. GOVERNING EQUATIONS With reference to ﬁgure 2, the electriﬁed jet is modeled as a chain of beads connected trough rectilinear elements, as in [11]. Each element is considered cylindrical with constant radius ai and length li . For incompressible ﬂuids, the conservation of mass can be expressed by saying that the volume of the element is constant, hence πa2i li = πa20 l0 050004-2 (1) 600 mi+1 , ei+1 500 G, α 400 μ s li μp 300 200 z 100 x y mi , ei 0 15 10 5 0 -5 -10 -15 -15 -10 -5 0 5 10 15 FIGURE 2. Discrete model for the electriﬁed jet. Leftmost panel: the polymer solution is modeled, according to the Giesekus model, using a nonlinear spring and dashpot element in parallel with a second dashpot to account for the Newtonian solvent viscosity. Rightmost panel: Example result of a jet instability simulation using the proposed model. in which the a0 and l0 are the radius and length of the initial pending drop. Hence, given the length li = (xi+1 − xi )2 + (yi+1 − yi )2 + (zi+1 − zi )2 , the element radius is readily available: ai = a20 l0 , li (2) (3) The position of the beads ri = xi i + yi j + zi k evolves according to the equation dri = vi , dt (4) where vi is the equation of the i-th bead. The equilibrium of forces on each bead reads as: mi dvi = FCi + FEi + FVi + FTi , dt where FCi = i j ei e j κ ri − r j |ri − r j |3 (5) (6) is the net of the Coulomb forces exerted by all beads on bead i, is the electrostatic force, V FEi = − k h (7) ⎡ 2 ⎤ a2 σi−1 ⎢⎢ a σi ⎥⎥ (ri − ri−1 )⎥⎥⎦ FVi = π ⎢⎢⎣ i (ri+1 − ri ) − i−1 li li−1 (8) 050004-3 is the viscoelastic force and FTi = − πγ (ai + ai−1 )2 Ki 1/2 ri · (i + j) 4 xi2 + y2i (9) is the surface tension. In the previous expressions mi and ei are respectively the mass and electric charge of the i-th bead, κ is the Coulomb constant, V is the intensity of the voltage source, h is the distance between the spinneret tip and the collector plate, γ is the surface tension, Ki is the local curvature of the jet and σ is the viscoelastic stress modeled with the Giesekus model [14]. μ s dli + τi li dt dτi αλ 2 μ p dli λ = −τi − τ + dt μp i li dt σi = (10) (11) The viscoelastic stress is split into a Newtonian contribution due to a solvent of viscosity μ s and an extra stress τ due to the polymer. This is characterized by relaxation time λ viscosity μ p . The coeﬃcient α takes into account the non isotropic mobility of the molecular chains and varies between 0 and 1. Note that for α = 0 the Oldroyd-B model is recovered and for α = 0 and η s = 0 the model is equivalent to a Maxwellian ﬂuid. Dimensional analysis Dimensionless groups can be formed on the basis of the following reference quantities: μ = μs + μ p , t∗ = λ , L ∗ = a0 v∗ = , a0 λ , σ∗ = μ . λ The problem can be therefore described in terms of the following dimensionless numbers: Re = m a0 λμ , V= eVλ2 ha0 m , Q= e2 λ2 a20 m , A= γπλ2 m , β= μs , μ where Re is the Reynolds number, V, Q and A express the ratio of the electrostatic, Coulombic, and surface tension forces with respect to the inertial forces, and β is the solvent fraction. NUMERICAL SOLUTION The governing equations are integrated using an explicit ﬁrst order Euler scheme. Using a temporal step Δt, the time is discretized as t = nΔt with n = 0, 1, . . . , Nt and Nt being the number of time steps. The superscript n is used to indicate that a certain quantity is evaluated at time tn . Given the beads positions rni and velocities vni at time tn , ﬁrst the beads position are updated using the following equation: The velocities are then updated rn+1 = rni + Δtvni . i (12) vn+1 = vni + mΔt FCi + FEi + FVi + FTi n i (13) In which all the forces are evaluated using the formulas 6,7,8 and 9. Finally the stress is updated: lin+1 − lin αλ n 2 n n τ + 2 = τ − Δt τ − τn+1 i i i μp i lin+1 + lnI (14) μ s lin+1 − lin + τn+1 . i lin+1 Δt (15) σn+1 = i The time stepping is started with a single segment of length l0 and diameter a0 , that is a pending drop. Every time that the length of the pending drop is doubled the ﬁrst element is split in two and a new bead with mass m and 050004-4 1 UCM Giesekus = 0.4 = 0.1 0.9 0.8 a/a 0 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 200 400 600 800 1000 1200 1400 1600 1800 s/a 0 FIGURE 3. Jet radius a as a function of the curvilinear coordinate s along the jet centerline: comparison between Upper Convected Maxwell and Giesekus models. The two ﬂuids have the same zero shear viscosity. charge e is added. In this way the bead chain grows throughout the simulation. The position of the newly added bead is perturbed by adding a small lateral displacement of the form xi = 10−3 l0 sin(ωt) yi = 10−3 l0 cos(ωt) (16) (17) This is needed for the onset of the instability. This procedure is iterated until the tip of the jet reaches the bottom. The evolution of the jet radius can be computed in a post-processing phase, once the numerical integration is over. In this paper simulations are run using reference values for the dimensionless groups from [11]. More speciﬁcally: Re = 0.083 , V=2 , Q = 12 , A = 0.9 , h/a0 = 629 , l0 /a0 = 1 , ωλ = 100 . In addition a Giesekus ﬂuid with identical zero shear viscosity and relaxation time is considered. For this the solvent fraction is chosen as β = 0.4 and the mobility α = 0.1. Figure 3 reports the evolution of the dimensionless jet radius as a function of the curvilinear coordinate along the jet centerline. It is evident the the jet of Giesekus ﬂuid is far more stretched than the one of Maxwellian ﬂuid. This is because, as mentioned before, the elongational viscosity of the UCM model grows indeﬁnitely as the elongational rate increases, whereas it remains bounded in Giesekus model. Therefore a higher stretching is required in case of the Giesekus ﬂuid in order to balance the Coulombic and electrostatic forces. As a consequence, the magnitude of the lateral swinging is ampliﬁed and the radius reduced in the case of Giesekus model. CONCLUSIONS In the present paper we modiﬁed the multi-bead model presented in [11], to include a more realistic constitutive law for the evolution of the stress tensor. The Giesekus model is known to be able to predict the behavior of polymer solutions more accurately than the UCM model. The results presented show how the diﬀerent rheological model signiﬁcantly aﬀects the radius of the jet, the magnitude of the lateral swings due to bending instabilities and ultimately the jet length. REFERENCES [1] G. I. Taylor, Proc. R. Soc. Lond. A: Math. Phys. Sci. 313, 453–475 (1969). 050004-5 [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] J. Jeans, The Mathematical Theory of Electricity and Magnetism (Cambridge University Press, Cambridge, 1958). D. H. Reneker and A. L. Yarin, Polymer 49, 2387 – 2425 (2008). Z. Huang, Y.-Z. Zhang, M. Kotaki, and S. Ramakrishna, Comp. Sci.Tech. 63, 2223 – 2253 (2003). N. Bhardwaj and S. C. Kundu, Biotechnology Advances 28, 325 – 347 (2010). K. Garg and G. L. Bowlin, Biomicroﬂuidics 5, p. 13403 (2011). D. A. Saville, Annu. Rev. Fluid. Mech. 29, 27–64 (1997). J. I. Ramos, Int. J. Numer. Meth. Fluids 23, 221–239 (1996). J. J. Feng, Phys. Fluid. 14, p. 3912 (2002). Q. G. Y. Wan and N. Pan, Int. J. Nonlin. Sci. Num. 5, p. 5 (2004). H. F. D. H. Reneker, A. L. Yarin and S. Koombhongse, J. Appl. Phys. 87, p. 4531 (2000). L. Xu, Y. Wu, and Y. Nawaz, Computers & Mathematics with Applications 61, 2116 – 2119 (2011). G. C. R. M. M. Hohman, M. Shin and M. P. Brenner, Phys. Fluid. 13, p. 2201 (2001). H. Giesekus, J. Non-Newtonian Fluid Mech. 11, 69–109 (1982). 050004-6

1/--страниц