{'fraction_non_alphanumeric': 0.03577962285392626, 'fraction_numerical': 0.014037433155080214, 'mean_word_length': 4.480046269519954, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Boltzmann Models of Charged Multiphase Flows for Electrospinning Simulations

We present a Shan-Chen model for charged leaky dielectric multiphase fluids in the context of electrospinning simulations. The role of non-linear rheology on the dynamics of electrified jets is considered by exploiting the Carreau model for pseudoplastic fluids. We report exploratory simulations of charged droplets at rest and under a constant electric field, and we provide results for charged jet formation under electrospinning-like conditions.


I. INTRODUCTION
The dynamics of charged jets presents a major interest, both as an outstanding problem in non-equilibrium thermodynamics, as well as for its numerous applications in science and engineering [1][2][3]. In particular, the recent years have witnessed a surge of interest towards the manufacturing of electrospun nanofibers, mostly on account of their prospective applications, such as tissue engineering, air and water filtration, optoelectronics, drug delivery and regenerative medicine [3][4][5]. As a consequence, several experimental studies have focused on the characterization and production of one-dimensional elongated nanostructures [5][6][7][8][9].
Electrospun nanofibers are typically produced at laboratory scale via the uniaxial stretching of a jet, which is ejected at a nozzle from an electrified charged polymer solution. The charged jet elongates under the effect of an external electrostatic field applied between the spinneret and a conductive collector and eventually undergoes electromechanical (e.g., whipping) instabilities due to various sources of disturbance, such as mechanical vibrations at the spinneret, hydrodynamic friction with the surrounding fluid and others [10]. While such instabilities can be detrimental in some respect, making an accurate position of individual fibers on target substrates very hard, in other experiments they are sought for, since they result in thinner cross sections, hence finer electrospun fibres, as they hit the collector [4].
This follows from a plain argument of mass conservation: whipping instabilities generate longer jets, hence thinner cross sections [11].
The computational modelling of the electrospinning process is based on two main families of techniques: particle methods and Lagrangian fluid methods. The former is based on the representation of the polymer jet as a discrete collection of discrete particles (beads) connected via elastic springs with frictional coupling (dissipative dashpots) and interacting via long-range Coulomb electrostatics [10,[12][13][14]. The latter, on the other hand, describe the jet as a continuum media, obeying the Navier-Stokes equations for a charged fluid in Lagrangian form [15][16][17]. Both methods are grid-free, hence well suited to describe abrupt changes of the jet morphology without taxing the grid resolution, as it is the case for Eulerian grid methods.
In this respect, grid-based methods, such as Lattice Boltzmann (LB), are not expected to offer a competitive alternative to the two aforementioned class of methods. Nevertheless, owing to its efficiency, especially on parallel computers, and its flexibility towards the in-clusion of physical effects beyond single-phase hydrodynamics, it appears worth exploring the possibility of using LB also in the framework of electrified fluids and jets. For instance, in the last decade significant improvements in LB methods for modelling microfluidic flows containing electrostatic interactions have been achieved [18,19], opening new applications of LB methods in electrohydrodynamic problems [20][21][22][23][24]. In particular, LB methods were successfully employed to simulate deformations and breakup of conductive vapor bubbles, bubble deformation due to electrostriction, dynamics of drops with different electric permittivity. All these investigations usually exploit the approach originally introduced by A.
L. Kupershtokh and D. A. Medvedev [19], where dielectric liquids are assumed with zero free charge density, so that the charge carriers are essentially locally bounded to the material [25]. Within this assumption, charge carriers are explicitly modeled by a convective transport equation solved by a second LB solver, taking into account the rates of ionization and recombination of charge carriers fluctuating around a local value (distribution) of equilibrium.
In the 1960s, G. I. Taylor provided several considerations for dealing with electrified fluid in a series of papers [26][27][28]. In particular, Taylor discovered that a moving charged fluids cannot be considered either as a perfect dielectric or as a perfect conductor. Instead, the fluid acts as a "leaky dielectric liquid", where a non-zero free charge is mainly accumulated on the interface between the charged liquid and the gaseous phase [29]. As a consequence, the charge produces electric stresses different from those observed in perfect conductors or perfect dielectrics. Indeed, in the last cases the charge induces a stress which is perpendicular to the interface, altering the interface shape to balance the extra stress. In the electrospinning process, the non-zero electrostatic field tangent to the liquid interface produces a non-zero tangential stress on the interface which is balanced from the viscous force [16].
The present work introduces a LB method for simulating charged leaky dielectric liquids, which is of main interest for modeling the electrospinning process. In this context, the largest part of the charge is modeled to lie along the interface between the liquid and the gaseous phase in similarity with previous works [30,31]. Further, the present LB approach is generalized to the case of non-Newtonian flows with shear-thinning viscosity in order to account the rheological properties of electrospun jets. In this context, we adopt an entropic approach [32,33] in order to preserve locally the second principle (H-theorem) also in presence of sharp changes in the fluid viscosity and structure.
The paper is organised as follows. In section II we present the basic features of the LB extension to the case of charged multiphase fluids. In sections III, we present results for the case of charged multiphase fluids at rest, and we report on preliminary results for charged multiphase jets under conditions related to electrospinning experiments.

II. MODEL
We consider a single species, charged fluid as composed of point-like particles and neglect correlations stemming from excluded volume interactions. Following Boltzmann's description, the state of the fluid is determined by the distribution function f p ( r, t) being the probability of finding at time t the fluid at position r and moving with discrete velocity c p , with p = 1, b given b the number of lattice directions. Here, the velocities c p are also viewed as vectors connecting a lattice site r to its lattice neighbors. The LB equation reads where the product αβ plays the role of a collision frequency, S p is a source term (see below), and f eq p is the continuum Maxwell-Boltzmann distribution computed at density ρ and velocity u [32]. The macroscopic variables are given by the density ρ = p f p and the fluid velocity u = 1/ρ p c p f . In the following, we refer to lattice units where the mesh spacing and timestep ∆t are conveniently set to unity. Also, we adopt the so-called D2Q9 scheme, composed by 8 discrete speeds (connecting first and second lattice neighbors) and one extra null vectors accounting for particles at rest. In this scheme, Here, the f eq p are chosen as a second-order Mach-number expansion where the w p are weights equal to 4/9 for the rest particles, 1/9 and 1/36 respectively for the smallest and largest velocities c p , and c s is the speed of sound that in lattice units is equal to 1/ √ 3. At the same time, we consider a unit fluid molecular mass, so that the thermal energy is equal to k B T = c 2 s with k B the Boltzmann constant and T the temperature. Following the approach of Refs. [32][33][34], the factor β in Eq. 1 depends on the kinematic viscosity ν by the relation while α is the largest value of the over-relaxation parameter so that the local entropy reduction can be avoided, ensuring the H-theorem. In particular, α is computed as the root of the scalar nonlinear equation [33,35] where H denotes the Boltzmann's entropy function, defined in discrete form [36] as In Eq. 1, the source term S p takes into account the global effect of all the internal and external forces F . This is assessed by the exact difference method proposed by Kupershtokh et al. [37], which reads where ∆ u = F /ρ. In bulk conditions, the LB method is intrinsically second-order accurate in space and time, and, in order to ensure the same accuracy in presence of forces, the local velocity is taken at half time step Here, the total body force F = F int + F el includes the inter-particle force F int and the electric force F el . The electric force acting on the boundary point r between a gas and a fluid with the local non-uniform permittivity ε( r) in an electric field E reads [25,29] where q is the local free charge carried on the fluid. In the last Eq., the vacuum permittivity ε 0 was assumed equal to 1 as in the Gaussian centimetre-gram-second (cgs) unit system, so that the charge in lattice units is length 3/2 mass 1/2 time −1 in similarity with the statcoulomb definition (note that Coulomb's constant is also 1). For the sake of convenience, we report in Appendix the units conversion Table in cgs dimensions from lattice units for several physical quantities shown in the following. As in Ref. [19], we consider a fluid with permittivity ε = 1 + ρ/ρ 0 with ρ 0 an arbitrary constant (in the following taken for simplicity equal to 1) so that it is ρ(∂ε/∂ρ) = ε − 1. As a consequence, Eq. 8 reduces to In the following, we assume that the magnetic induction effects can be neglected so ∇ × E = 0, and the system follows the Gauss law ∇ · (ε E) = q. Since E = −∇φ with φ the electric potential, the Poisson equation div(ε( r)∇φ) = −q( r) can be solved at each lattice node r, given the boundary conditions of the system and the local charge q( r) at the node (specified below). In particular, we determine the electric potential by solving numerically the two-dimensional Poisson equation by means of a SOR (Successive Over-Relaxation) algorithm and the Gauss-Seidel method [38]. Note that the Poisson equation includes the non-uniformity of the permittivity ε( r), and it is solved on-the-fly during the simulation. Hence, the electric force F el = q(−∇φ) is added into the LB scheme by Eq. 6.
Since we are modeling a leaky dielectric fluid, we assume that the free charge in the system is mainly distributed over the liquid-gaseous interface. Further, in similarity to previous electrospinning models [16,17], the relaxation time of free charge in the system is assumed to be irrelevant. In other words, the free charge in bulk liquid relaxes to the liquid interface in a smaller time than any other characteristic time in the system [39]. This is a well-established assumption of a leaky dielectric fluid (for further details see Ref. [29]). The liquid charge in the point r is given as which is the sum of a surface charge q s and a small bulk term q b . The bulk term q b is taken as where Q b denotes the total charge in the bulk, the denominator acts to keep constant the charge due to the charge conservation principle, and θ(ρ; ρ 0 ) denotes a smoothed version of the Heaviside step function switching from zero to one at ρ 0 (equal to 1 in all the following simulations) in order to select only the liquid phase. The term q s is modeled as a proportional to the absolute density gradient where Q s denotes the total charge over the surface, and the denominator ensures the charge conservation principle as in the previous case. This approach is usually referred to as the constant surface charge model, and it was already adopted in Refs [30,31] as a strategy to simplify the charge transport and distribution on the droplet interface. Nonetheless, the constant surface charge model fails in describing a distributed charge on the drop interface whenever the charge density is high, since the curvature surface alters the local charge density [16]. In order to address the issue, we assume that the curvature biases the surface charge density as in a conductive liquid, following the power-law introduced by I. W. McAllister [40], which states Here, K denotes the mean curvature K = ∇ · n with the local interface normal n = ∇ρ( r)/|∇ρ( r)| [41], while q s,max is the maximum surface charge at the maximum curvature K max chosen as a reference value for the system under investigation. It is worth to emphasize that treating a leaky dielectric as a conductive liquid is a simplification already made by several authors (e.g., G. Taylor [26], A. Yarin et al. [42] , etc.). For the sake of simplicity, we take in the following the maximum curvature K max equal to K d value, defined as the curvature doubling the local surface charge density. Thus, we rewrite Eq. 10 as The total charge of the system is conserved and equal to Q = Q b + Q s .
In addition, the fluid is subjected to an internal thermodynamic force F int promoting a phase separation. The phase separation force is accounted for by means of the Shan-Chen method [43]. We construct the local force as with the sum p∈f luid running over lattice nodes where the fluid is allowed, that is, not belonging to the wall, and p∈wall runs over nodes belonging to the wall. G and G w are fluidfluid and fluid-wall interaction strengths, respectively. In Eq 15, ψ is an effective number density, which is taken for simplicity ψ(ρ) = ρ 0 [1 − exp(−ρ/ρ 0 )], being ρ 0 an arbitrary constant [44] (in the following assumed equal to 1).

A. Extension to Non-Newtonian flows
In the electrospinning process, the rheological behavior of polymeric liquid with shearrate-dependent viscosity is expected to play a significant role on jet dynamics. As a consequence, we now generalize the present model to the case of non-Newtonian flows, in similarity with the approach reported in Refs [45][46][47]. The shear-rateγ is a functional of the density distribution function f . In particular, the strain tensor Γ η,δ reads [45,48] where is the the stress tensor with η and δ running over the spatial dimensions. Note that τ in Eq. 16 is defined as the inverse of the product αβ, where α was computed by Eq. 4, and β depends by Eq. 3 on the kinematic viscosity ν. We now rewrite Eq. 16 aṡ whereγ and Π are computed as matrix 2-normγ = 2||Γ|| 2 and Π = ||Π|| 2 of the shear and stress tensor [45], respectively. Note that in the last Eq. we exploit a constitutive relation between the kinematic viscosity ν and the shear-rateγ, so that τ = τ (γ). As a consequence, γ is computed as the root of the scalar nonlinear Eq. 18.
We should now consider the general trend observed in electrospun polymeric filaments [49]. As main features, we highlight that a polymeric spinning solution at low shear rate behaves as a quasi-Newtonian fluid with viscosity ν 0 , since the initial condition can be recovered, while at high shear rate a non-reversible disentanglement is present. In particular, it is possible to identify a retardation time λ at which the shear-thinning starts, which is equal to the inverse value of the shear rate at that instant. At very high shear rate, a quasi-Newtonian behavior is again observed as soon as the alignment of the polymer chains is extremely high (almost complete). The last region is characterized by a final viscosity value (infinite viscosity ν ∞ ), which is lower than ν 0 .
In the present investigation, we exploit the Carreau model [50], which is able to describe all the mentioned rheological properties. The Carreau model states that where n is the flow index (n < 1 for a pseudoplastic fluid). Obtainedγ by resolving Eq.
18 and ν(γ) by Eq. 19, and assuming a slow variation of ν(γ) over the time ∆t, the local parameter β is finally estimated by Eq. 3. Note that a validation of a similar implementation in LB scheme of the Carreau model was given in Ref. [45].

A. Charged drop
In order to assess the properties of our implementation for charged multiphase systems, we have initially run a set of simulations modeling a charged leaky dielectric fluid system obeying the Shan-Chen equation of state [44]. In order to assess the static behavior, we take a two-dimensional periodic mesh made of 320 x 320 nodes, and prepare the system by Since we aim to model a leaky dielectric fluid, the ratio Q s /Q b is taken equal to 10, so that the largest part of the charge lies over the surface. The total charge Q = Q s + Q b was set equal to 2.13, and q( r) computed by Eq. 14. Whenever the Poisson equation is solved, a uniform negative charge is added to obtain a system with net charge zero. Hence, the electric force F el = q(−∇φ) is added into the LB scheme by Eq. 6, where q is computed by Eq. 14 with K d equal to 1. The liquid is Newtonian with kinematic viscosity ν = 1/6. The stationary configuration of the described system is obtained after 1000 time steps.
Hence, we inspect the electric field (see Fig. 1) at rest conditions in order to analyze the balance of forces acting at the interface, including Shan-Chen pressure and capillary and electrostatic forces, the latter pointing normal to the interface (see panel b of Fig. 1). Here, we observe that the largest part of the electric field is located over the liquid surface where the charge distribution is higher. In particular, at the boundary of the drop the magnitude of the electric field | E| is equal to 3.5 · 10 −3 .
It is of interest to estimate the various forces which concur to provide a stable configuration of the droplet. The mechanical balance reads as follows: where p L and p V are the liquid and vapour pressure, respectively, p cap = σ/R is the capillary pressure (given the surface tension σ and the drop radius R), and p el is the repulsive electrostatic pressure. The latter can be estimated by standard considerations in electrostatics, namely: where E s is the electric field at the surface.
In actual numbers with σ = 5.8·10 −2 , this ratio is equal to 1/50. This shows that electrostatic forces act as a small perturbation on top of the neutral multiphase physics.
Next, we investigate the effects of a uniform external electric field E ext of magnitude pointing along the x axis. Using the previous configuration at the equilibrium as starting point of our simulation, we set E ext at two different values equal to 0.1 and 0.5. For each one of the two cases, we report a snapshot of the fluid density ρ taken as soon as the liquid drop touches the point of coordinates (280,160) in lattice units. The set in Fig. 2 highlights the significant motion of the charged drop towards right in accordance with the direction of the electric field. In the figure, we note that a sizable change in the drop shape is present only for the case E ext = 0.5. In order to elucidate this effect, it is instructive to assess the strength of the electrostatic field in units of capillary forces, namely: In actual numbers, this ratio is equal to 0.5 and 3 for the case at lower and higher E ext , respectively. This shows that the electric force magnitude is sufficiently large to provide an alteration of its shape only in the case at higher E ext . In particular, the shape shows an elongation towards the direction of the electric field E ext , which results from the effect of the curvature on the surface charge.
In Fig. 3, we report the alteration of charge density due to the mean curvature term K of Eq. 14. The alteration is estimated as δ q = q − q K=0 , where q K=0 is computed with Eq. 14 with K = 0 everywhere, corresponding to a constant surface charge model without curvature effect correction. Here, we note an accumulation of charge on the rightest part of the drop, where the mean curvature K shows a maximum value equal to 5.17 · 10 −2 , which corresponds to a charge accumulation δ q equal to 1 · 10 −3 , in the following denoted δ + q . The accumulation of charge is counterbalanced by a negative charge δ − q distributed over the almost straight surface part of the drop (just behind the rightest protrusion in Fig. 3).
Both partial charges δ + q and δ − q favor the presence of a protrusion in the drop shape. In order to analyze this effect, we report in Fig. 4 the mean curvature computed at the same time t = 1250δt for two simulations, both at E ext = 0.5 differing for the inclusion of the curvature effects in the constant surface charge model of Eq. 14. Even though the circular shape of the drop is deformed by the external electric field in both cases, the curvature effects increases the protrusion on the drop shape (see Fig. 4 panel b). Further, the charge differences provide a shift in the electric force acting on the drop surface, the effect of which accumulates in time, so that the alteration in the drop shape increases in time.

B. Charged multiphase jet in electrospinning setup
We set up a system modeling the electrospinning process, containing a charged Shan-Chen fluid. The system is a mesh made of 320 x 320 nodes (see Similarly, we set outlet nodes on the right side (at x = 320) where the jet will impinge under the effect of the external electric field. Such electric field is chosen to mimic the potential difference that is normally applied between the nozzle and a conductive collector in the real electrospinning setup [51]. The computational setup is quite sensitive to the choice of    The initial density of the two phases are 2.0 and 0.16 for the liquid and gaseous phase, respectively. The initial configuration consists of the liquid phase filling the inner space of nozzle with a liquid drop just outside the needle (see Fig. 5 panel a). All remaining fluid nodes are initialized to gaseous density.
Both the Shan-Chen constants for the fluid-fluid G and fluid-wall G w interaction are set to −5. As in previous section, for the resolution of the Poisson equation a uniform negative charge is added to the system in order to counterbalance the positive charge and obtain a system with net charge zero. Further, we impose Dirichlet boundary condition in the following form: we impose the electric potential φ l = 0 on the left side (x = 0), while the electric potential on the right side (x = 321) was set equal to φ r = −32.2, providing a background electric field E back = (φ r − φ l )/322, which is imposed between the two opposite sides (left-right) of the system. On the upper and bottom sides, the Dirichlet boundary conditions are set equal to the φ u (x) = φ b (x) = x(φ r − φ l )/322. Note that the last condition is equivalent to impose an electric field E of magnitude 0.1 oriented along the x-axis.
Since in a typical electrospinning setup the liquid phase is always connected to a generator addressing a charge, it is reasonable to assume that, whenever the stretching of the liquid jet increases the jet interface, extra charge rapidly reaches the liquid boundary in order preserve the value of charge surface density. As a consequence, the charge conservation condition can not be applied (the liquid jet is not isolated). Instead, we assume the conservation condition of the surface charge density value for the same mean curvature, so that Eq. 14 is rewritten as where we have adopted a similar condition also for the bulk charge term, being ξ b and ξ s two proportionality constants, in the following taken equal to 1 · 10 −4 and 6 · 10 −2 , respectively.
Note that the two proportionality constants were tuned in order to obtain a mean ratio  The internal electric field computed by the Poisson solver (see Section 2) is computed on-the-fly during the simulation. In panel b of Fig. 5 we report the electric field magnitude | E| for the initial configuration. Here, we note a maximum value in | E| close to the drop interface, which is due to the higher surface charge density. Further, a lower value of | E| is observed in the nozzle as a consequence of the larger dielectric constant ε 3 in the liquid phase (versus ε 1 in the gaseous phase).
We now report in Fig. 6 several snapshots taken over the time evolution of the system labeled case b. In all the cases, we observe the formation of a liquid charged jet, which is  [10,52]. After the jet touches the collector, the jet shape fluctuates around a mean profile, providing a stationary regime. In particular, at this stage the jet shows an hyperbolic profile (see Fig. 8 panel b) which appears to be in qualitative agreement with the characteristic shape of the jet experimentally observed close to the injecting nozzle by the Rafailovich and Zussman groups [53] (see Fig. 8 panel a) and in consistency with previous theoretical results on the jet conical shape [11,[54][55][56][57][58].
In order to characterize the stationary regime, we report in panel c of Fig. 8 the magnitude of the velocity field, and the line integral convolution (LIC) visualization technique [59], highlighting the fine details of the flow field. As expected, we observe a higher value of u(x, z) along the jet towards the collector. In particular, we analyze the profile of the velocity in a jet cross section along what is generally observed in the experimental process.
We investigate the effect of pseudoplastic rheology on the stress tensor Π. In panel a of Fig. 9 we report the mean value of the stress Π measured along the central axis of the jet y = 160, and averaged over a time interval of 15000 steps in the stationary regime for the three cases under investigation. Here, we note a decreasing trend of the the stress by increasing the pseudoplasticity of the fluid (decreasing the flow index n). Nonetheless, we observe a small shift in the stress magnitude. This is essentially due to the low value of retardation time λ in the Carreau model adopted in our simulation, which provides a small decrease in the kinematic viscosity ν. In order to clarify the this point, we report in panel b of Fig. 9 the mean value of the kinematic viscosity < ν > again assessed along the central axis of the jet y = 160, and averaged over the stationary regime time, where we observe a small decrease of < ν > along the stretching direction as a function of the pseudoplastic behavior in the fluid. These results look promising, and we plan to investigate systematically the effect of the rheological parameters on the jet dynamics in a future work.

IV. SUMMARY
Summarizing, we developed a Shan-Chen model for charged leaky dielectric fluids mainly aimed at modeling the electrospinning process. The curvature effects on the charge surface were included in our theoretical treatment, and we generalized the model to non-Newtonian flows in order to account for the peculiar rheological behavior. Different scenarios were investigated to test the model. We initially investigated the effect of strong electric fields on the droplet shape evolution. We also probed the jet formation under electrospinninglike conditions, obtaining a good agreement with both experimental results and previous theoretical works present in literature. At first glance, the pseudoplastic behavior alters the jet dynamics, although a more systematic investigation requires an extensive test of the rheological parameters. Work along these lines is currently underway.
The preliminary applications of the presented LB model look promising, although more systematic numerical investigations, as well as theoretical analysis need to be undertaken.
Nonetheless, the actual effort can be regarded as a significant forward step to extend the applicability of the LB method to the context of electrospinning systems, providing a useful computational tool in completion of the others presently available in literature. ("COPMAT").  The conversion to 3-d dimension units is obtained by multiplying the 2-d square area ∆l 2 for an extra node ∆l along the z-axis. Note that ∆m denotes the mass unit equal to 5 · 10 −10 g (water density), ∆l the length unit equal to 10 −3 cm, and time step ∆t is the time unit equal to 8.3·10 −7 s.

Appendix: Units Conversion
The last number was tuned in order to set the kinematic viscosity ν = 1/6 ∆l 2 ∆t −1 = 0.2 cm 2 s −1 in similarity with the viscosity value adopted in the Lagrangian model of Ref. [12].