Whipping of electrified visco-capillary jets in airflows

TL;DR

This paper introduces a novel model framework using a Cosserat rod description to analyze the whipping instability in electrified visco-capillary jets, enabling efficient computation of stable stationary solutions considering multiple physical effects.

Contribution

It proposes a new modeling approach for jet whipping instability based on a Cosserat rod framework with a continuation-collocation numerical method.

Findings

The model captures the qualitative behavior of jet whipping.

Numerical simulations demonstrate the method's efficiency and robustness.

The approach allows exploration of high-dimensional parameter spaces.

Abstract

An electrified visco-capillary jet shows different dynamic behavior, such as cone forming, breakage into droplets, whipping and coiling, depending on the considered parameter regime. The whipping instability that is of fundamental importance for electrospinning has been approached by means of stability analysis in previous papers. In this work we alternatively propose a model framework in which the instability can be computed straightforwardly as the stable stationary solution of an asymptotic Cosserat rod description. For this purpose, we adopt a procedure by Ribe (Proc. Roy. Soc. Lond. A, 2004) describing the jet dynamics with respect to a frame rotating with the a priori unknown whipping frequency that itself becomes part of the solution. The rod model allows for stretching, bending and torsion, taking into account inertia, viscosity, surface tension, electric field and air drag. For the resulting parametric boundary value problem of ordinary differential equations we present a continuation-collocation method. On top of an implicit Runge-Kutta scheme of fifth order, our developed continuation procedure makes the efficient and robust simulation and navigation through a high-dimensional parameter space possible. Despite the simplicity of the employed electric force model the numerical results are convincing, the whipping effect is qualitatively well characterized.

Figures

Tables

Table 1. Table 3.1. Performance of continuation procedure. Computing time of each Part A-D and the number of respective continuation steps to reach the reference parameter tuple p r ​ e ​ f superscript p 𝑟 𝑒 𝑓 \mathrm{p}^{ref} (cf. Fig.  3.1 ).

Part	Computing Time [s]	Percentage	Continuation Steps
A	173.64	8.6%	61
B	660.14	32.6%	136
C	190.96	9.4%	25
D	1000.04	49.4%	112
total	2024.78	100%	334

Equations

$${\partial }_s\mathbf{r}={\mathbf{d}}_3,{\partial }_s{\mathbf{d}}_{\mathbf{i}}=\mathbf{\kappa }\times {\mathbf{d}}_{\mathbf{i}},$$

$${\partial }_{t}A+{\partial }_s(uA)$$

$$\rho {\partial }_t(A\mathbf{v})+\rho {\partial }_s(uA\mathbf{v})$$

$$\rho {\partial }_t(\mathbf{J}\cdot \mathbf{\omega })+\rho {\partial }_s(u\mathbf{J}\cdot \mathbf{\omega })$$

$$A=\pi a^2,$$

$$\mathbf{n}\cdot {\mathbf{d}}_3=3\mu A{\partial }_{s}u,$$

$${\mathbf{f}}_{ca}=\pi \gamma {\partial }_s(a{\mathbf{d}}_3),$$

$$\mathbf{F}(\mathbf{\tau },\mathbf{w})=w_{\nu }{r}_{\nu }(w_{\nu })\mathbf{\nu }+w_{\tau }{r}_{\tau }(w_{\nu })\mathbf{\tau }$$

$$r_{\nu }(w_{\nu })$$

$$r_{\tau }(w_{\nu })$$

$$r_{\nu }^S=\frac{4\pi }{\log (4/\delta )}-\frac{\pi }{{\log }^2(4/\delta )},r_{\tau }^S=\frac{2\pi }{\log (4/\delta )}+\frac{\pi /2}{{\log }^2(4/\delta )}$$

$$p_{\nu ,0}$$

$$p_{\tau ,0}$$

$$q_{k,0}$$

$${\partial }_t{\mathbf{d}}_{\mathbf{i}}^s=(\mathbf{\omega }-u\mathbf{\kappa }+\Omega {\mathbf{d}}_3^s)\times {\mathbf{d}}_{\mathbf{i}}^s,i=1,2,3.$$

$$\mathbf{y}=\sum _{i=1}^3{y}_i{\mathbf{d}}_{\mathbf{i}}^s=\sum _{i=1}^3{\breve{y}}_i{\mathbf{a}}_{\mathbf{i}}\in {\mathbb{E}}^3$$

$$\mathbf{M}(t)=\cos (\Omega t){\mathbf{P}}_0(0,0)+\sin (\Omega t){\mathbf{a}}_3\times {\mathbf{P}}_0(0,0)+{\mathbf{a}}_3\otimes {\mathbf{a}}_3$$

$$\mathbf{r}(0)=0,{\mathbf{d}}_{\mathbf{i}}^s(0)={\mathbf{a}}_{\mathbf{i}},\mathbf{\kappa }(0)=0,u(0)=U,\mathbf{n}(L)=0,\mathbf{m}(L)=0.$$

$${\breve{r}}_1(L)=0,{\breve{r}}_3(L)=H.$$

$$L^{-1}\mathsf{R}\cdot {\partial }_s\breve{\mathsf{r}}$$

$$L^{-1}{\partial }_s\mathsf{R}$$

$$L^{-1}{\partial }_s\kappa$$

$$L^{-1}{\partial }_{s}u$$

$$L^{-1}{\partial }_s\mathsf{n}$$

$$+\mathrm{Re}\frac{{\Omega }^2}{u}\mathsf{R}\cdot \left({\mathsf{e}}_3\times ({\mathsf{e}}_3\times \breve{\mathsf{r}})\right)-{\mathsf{f}}_{ca}-{\mathsf{f}}_{el}-{\mathsf{f}}_{air},$$

$$L^{-1}{\partial }_s\mathsf{m}$$

$$\displaystyle\quad-\frac{\mathrm{Re}}{4}\frac{\Omega}{u}\mathsf{P_{2}}\cdot\bigg{(}\frac{1}{3}\mathsf{R}\cdot\mathsf{e_{3}}n_{3}-\frac{1}{3}\mathsf{e_{3}}n_{3}+\left(\mathsf{\kappa}-\frac{\Omega}{u}\mathsf{e_{3}}\right)\times\mathsf{R}\cdot\mathsf{e_{3}}\bigg{)}$$

$$-\frac{\mathrm{Re}}{4}\left(\frac{1}{u^2}{\mathsf{P}}_2\cdot (u\kappa -\Omega {\mathsf{e}}_3+\Omega \mathsf{R}\cdot {\mathsf{e}}_3)\right)\times \left(u\kappa -\Omega {\mathsf{e}}_3+\Omega \mathsf{R}\cdot {\mathsf{e}}_3\right)$$

$${\mathsf{f}}_{ca}$$

$${\mathsf{f}}_{el}$$

$${\mathsf{f}}_{air}$$

$$\breve{\mathsf{r}}(0)$$

$${\breve{r}}_1(1)$$

$$\mathrm{Re}$$

$$\mathrm{M}$$

$$\displaystyle\mathsf{R}(\mathsf{q})=\left(\begin{array}[]{ccc}q_{1}^{2}-q_{2}^{2}-q_{3}^{2}+q_{0}^{2}&2(q_{1}q_{2}-q_{0}q_{3})&2(q_{1}q_{3}+q_{0}q_{2})\\ 2(q_{1}q_{2}+q_{0}q_{3})&-q_{1}^{2}+q_{2}^{2}-q_{3}^{2}+q_{0}^{2}&2(q_{2}q_{3}-q_{0}q_{1})\\ 2(q_{1}q_{3}-q_{0}q_{2})&2(q_{2}q_{3}+q_{0}q_{1})&-q_{1}^{2}-q_{2}^{2}+q_{3}^{2}+q_{0}^{2}\end{array}\right),$$

$$L^{-1}{\partial }_s\mathsf{q}=\mathcal{A}(\kappa )\cdot \mathsf{q},\mathsf{q}(0)=(1,0,0,0),$$

$$\frac{\mathrm{d}}{\mathrm{d}s}\mathsf{y}=\mathsf{f}(\mathsf{y}),\mathsf{g}(\mathsf{y}(0),\mathsf{y}(1))=0$$

$${\mathsf{y}}_{i+1}-{\mathsf{y}}_i-h_{i+1}\sum _{j=1}^4{b}_j{\mathsf{k}}_j=0,\mathsf{g}({\mathsf{y}}_0,{\mathsf{y}}_N)=0$$

$$\displaystyle\text{with }\mathsf{k}_{j}=\mathsf{f}\bigg{(}\mathsf{y}_{i}+h_{i+1}\sum\limits_{k=1}^{4}a_{jk}\mathsf{k}_{k}\bigg{)}$$

Full text

Keywords. electrified jets, lateral instabilities, electrospinning, viscous Cosserat rod model, parametric boundary value problem, homotopy method

AMS-Classification. 34B08, 65Lxx, 76-XX

1. Introduction

The interaction of an intense electrical field with the interface between a conducting liquid and a dielectric medium has been known to exist since Gilbert [7] reported in 1600 the formation of a conical meniscus when an electrified piece of ambar was brought close enough to a water drop. The deformation of the interface is caused by the force that the electric field exerts on the net surface charge induced by the field itself. This phenomenon is at the base of modern devices for the production of micro- and nanostructures of interest in several technological fields [18, 14]. As schematized in Fig. 1.1 (left), these devices consist essentially of a high-voltage power supply, a metallic needle (spinneret) and a grounded collector (counter-electrode). The metallic needle is connected to a syringe pump through which a conducting liquid can be fed at a constant and controllable rate. When a high voltage is applied, the electric field induces an electric current in the liquid that accumulates electric charge at the surface and causes an electric force that elongates the pendant drop at the needle’s exit in the direction of the field. It is observed that if the field strength is below a certain threshold value the balance of electrostatic and surface tension forces gives rise to a motionless conical shape commonly known as the Taylor cone [28]. However, above the threshold, the large electrostatic forces concentrated near the cone tip overcome the surface tension stresses and force the ejection of an electrified liquid jet from the cone tip. For certain values of the applied voltage and imposed liquid flow rate, the jet emanating from the cone tip is stationary and breaks into spherical droplets at some distance downstream due to axisymmetric Rayleigh-Plateau (varicose) instabilities corrected to account for the presence of surface charge. This so-called cone-jet mode forms the basis of the electrospray technique [5, 6] for generating small monodisperse drops with great applications in fine coatings, synthesis of powders, micro- and nanocapsules, etc. However, non-symmetric perturbation modes can also grow due to the net charge carried by the jet. Indeed, if a small portion of the charged jet moves slightly off axis, the charge distributed along the rest of the jet will push that portion farther away from the axis according to Earnshaw’s theorem, thus leading to a lateral instability known as whipping or bending instability. If the growth rate associated to this whipping instability is larger than that associated to varicose jet break-up – as may happen, for example, for sufficiently high values of the applied voltage or of the liquid viscosity –, the off-axis movement of the jet becomes the most significant aspect of its evolution, see Fig. 1.1 (middle). The whipping mode manifests itself in the form of fast and violent slashes which give rise to very large tensile stresses and to a dramatic jet thinning. This is of fundamental importance in the electrospinning process [3, 11, 15, 21, 27, 32], where micro- or nanofibers of a polymeric fluid are produced by solidification of the jet issuing from the Taylor cone before it breaks up into droplets. The reduction of the jet diameter, that is typically several orders of magnitude makes the electrospinning technique very competitive with other existing ones (such as phase separation or self-assembly) and subject of research. Another mode, the so-called coiling mode [13], is observed when the ground electrode is located sufficiently close to the needle’s exit so that the liquid jet reaches the plate before being set into chaotic motion by the whipping instability, see Fig. 1.1 (right). The situation is then the same as if a thin stream viscous fluid such as honey is poured onto a surface from a certain height [22, 23, 24]. Rather than approaching the surface vertically the jet builds on it a helical structure which resembles a pile of coiled rope.

In the last decades intensive work has been spent on experimental studies of electrospinning, whose complexity makes empirical determination of the effects of the parameters very difficult [27, 15]. Theoretical models were developed to predict the behavior of the charged jet, see the trend-setting papers [21, 10, 4]. The models were set up in terms of balance equations on basis of slender-body descriptions for the jet, including the electric field and surface tension, differing in the constitutive laws, such as linear and nonlinear, Newtonian and upper-convected Maxwell, models for viscous and viscoelastic jet behavior, etc. Also effects like evaporation and solidification have been taken into account [32]. It turned out [31] that the final jet radius is mainly determined by the volumetric charge density, the distance from the nozzle to the collector, the orifice radius, the relaxation time and the viscosity. The external field and mutual electric interaction of multiple charged jets influence the jet path and evolution [29]. The whipping instability, which causes the jet thinning, was documented and analyzed in [21, 32] by applying perturbation frequencies. On top of a linear stability analysis operating diagrams for different polymer solvents were developed in [10, 11]. The instability was also explored in temporal stability analysis in [16, 17]. Numerical simulations are given in [33]. However, it must be realized that the nature of the whipping regime makes it very difficult to unravel its detailed structure. Recent experiments show that there are some circumstances which greatly enhance the parametric range for which the bending leads to a stable steady-state helicoidal structure with a constant opening angle as, for example, when the conducting liquid is surrounded by a dielectric bath [25] or by another coflowing liquid [8] (see Fig. 1.1, middle). The liquid bath makes the frequency of the jet oscillations several orders of magnitude lower than those found in typical electrospinning experiments in air [11, 20]. In the phenomenological study on a glycerine jet in a hexane bath [25] the whipping instability (frequency, amplitude, wavelength) is characterized in terms of flow rate, applied voltage and electric conductivity.

The observation of the stable steady-state helical structure motivates the idea of this paper. We propose a model framework in which the whipping instability can be computed straightforwardly as the stable stationary solution of an asymptotic Cosserat rod description. For this purpose we formulate the jet’s whipping as a stationary process by introducing a frame that rotates with the a priori unknown whipping frequency, so the frequency becomes part of the solution. A similar transformation has been used in the investigation of viscous rope coiling [22]. In the whipping instability this results in a parametric boundary value problem of ordinary differential equations which describes the jet behavior in terms of 19 variables in dependence on viscosity, surface tension, external electric field, self-repulsion of the induced charges and the resistance of the surrounding fluid. The used model for the electric and capillary forces is taken from [32], the material and geometrical models for the incompressible viscous rod come from the work [22] on viscous rope coiling and [2, 1] on rotational spinning. For the resistance of the jet surrounding airflow we employ a drag model developed in [19]. For the numerical treatment of the problems we propose a continuation algorithm that allows the efficient and robust simulation and automatic navigation through a high-dimensional parameter space. The underlying collocation is performed with a Lobatto IIIa formula (implicit Runge-Kutta scheme of fifth order), and the resulting nonlinear system is solved with a Newton method. We explore our framework in a parameter study. Although the model is physically comparatively simple in the electric force model, the numerical results are very convincing. They show qualitatively well the characteristic jet behavior in terms of whipping frequency, elongation and throwing range. The presented model framework is able to give a quantitative explanation for the observed strong jet thinning. Indeed, the analysis of the periodic equilibrium states in combination with a temporal global stability analysis reveal the helical structure in the whipping instability in more details as it is done in previous works. In addition, it enables the derivation of an analytical solution for the ’longtime’ behavior of the jet.

The paper is structured as follows. Starting with a short introduction into the viscous Cosserat theory, we derive the stationary model framework for the jet’s whipping in Sec. 2. In Sec. 3 we present a continuation-collocation method for the numerical solution of the parametric boundary value problem and demonstrate the computational efficiency of the proposed algorithm. Numerical results are shown and discussed in Sec. 4 with respect to parameter studies and a detailed investigation of the whipping effect. In particular, an analytical solution for the jet’s ’longtime’ behavior is stated. Moreover, we include a temporal stability analysis of the underlying transient Cosserat rod model in App. A. Appendix B gives details to apparent boundary layers in the numerical solutions.

2. Viscous Cosserat Rod Model

2.1. Electrified visco-capillary jet in airflow

As a jet is a slender long object, its dynamics can be reduced to an one-dimensional description by averaging the underlying balance laws over its cross-sections. In the special Cosserat rod theory there are two constitutive elements: a curve $\mathbf{r}:\mathcal{D}\rightarrow\mathbb{E}^{3}$ specifying the jet position (e.g. midline) and an orthonormal director triad $\{\mathbf{d_{1}},\mathbf{d_{2}},\mathbf{d_{3}}\}:\mathcal{D}\rightarrow\mathbb{E}^{3}$ that is attached to the curve and characterizes the orientation of the cross-sections in the three-dimensional Euclidean space $\mathbb{E}^{3}$. The Euclidian space $\mathbb{E}^{3}$ can be identified with $\mathbb{R}^{3}$ by choosing a basis. We particularly use the domain $\mathcal{D}=\{(s,t)\in(\mathbb{R}_{0}^{+})^{2}\}$ where $t$ is the time and $s$ is the arc length parameter to impose an Eulerian (spatial) description. In the following we consider an electrified visco-capillary jet with circular-shaped cross-sections surrounded by an airflow. To describe its dynamics we proceed from the incompressible viscous Cosserat rod model of [2, 1] that was derived for jets in rotational spinning processes on the basis of the work [22] on viscous rope coiling. The rod model consists of balances for mass (cross-section), linear and angular momentum and allows for stretching, bending and torsion. The electric and capillary forces are included according to [32], whereas the air drag model is taken from [19]. The resulting rod system is given by the following four kinematic and three dynamic equations

[TABLE]

that are supplemented with a relation between cross-sections and moments of inertia, with viscous material laws and with outer forces

[TABLE]

with the scaling tensors $\mathbf{P}_{k}=\mathbf{d_{1}}\otimes\mathbf{d_{1}}+\mathbf{d_{2}}\otimes\mathbf{d_{2}}+k\,\mathbf{d_{3}}\otimes\mathbf{d_{3}}$, $k\in\mathbb{R}$.

The unknowns of the system (2.1) are the jet curve $\mathbf{r}$, triad $\{\mathbf{d_{1}},\mathbf{d_{2}},\mathbf{d_{3}}\}$, curvature $\boldsymbol{\kappa}$, cross-section $A$, linear $\mathbf{v}$ and angular $\boldsymbol{\omega}$ velocities as well as the convective speed $u$ and the normal contact force components $\mathbf{n}\cdot\mathbf{d_{i}}$, $i=1,2$. The curve and the triad are coupled by $\partial_{s}\mathbf{r}=\mathbf{d_{3}}$. In the chosen Eulerian description this relation for the jet tangent contains the arc length parameterization and the generalized Kirchhoff constraint that allows for stretching and prevents shearing as jet deformation. In this context the convective speed (parameter speed) $u$ can be viewed as the Lagrange multiplier (unknown) to the arc length parametrization. The further kinematic equations relate the jet curve and the triad to the curvature $\boldsymbol{\kappa}$ and the linear and angular velocities $\mathbf{v}$, $\boldsymbol{\omega}$. The mass density $\rho$ is considered to be constant and the cross-section $A$ to be circular-shaped of radius $a$. The coupling of the angular momentum line density with the moment of inertia $\mathbf{J}$ preserves the jet’s incompressibility. The tangential contact force $\mathbf{n}\cdot\mathbf{d_{3}}$ and the couple $\mathbf{m}$ are specified by a linear material law in the strain rate variables with dynamic viscosity $\mu$, whereas the normal force components are the Lagrange multipliers to the generalized Kirchhoff constraint. The acting outer line force densities consist of the capillary line force density $\mathbf{f}_{ca}$ with surface tension coefficient $\gamma$, the electric line force density $\mathbf{f}_{el}$ and the aerodynamic line force density $\mathbf{f}_{air}$. The last is based here on a stationary non-moving flow situation where the density $\rho_{\star}$ and the dynamic viscosity $\mu_{\star}$ of the air are considered to be constant.

The electric forces are assumed to be split into two parts, the effects due to the external unperturbed electric field $\mathbf{E}$ and the ones due to the Coulomb interactions of the induced charges on the jet. The self-repulsion is modeled by help of a local interaction approximation [32] in terms of the surface charge density $\sigma$, the permittivity $\varepsilon_{p}$ and the distance between nozzle and counter-electrode (device height) $H$ as typical (cut off) length. The external electric field depends on the applied voltage $\Phi$ and the device geometry. We assume it to be constant, i.e., $\mathbf{E}=E\mathbf{a_{3}}$ with $E=\Phi/H$ and $\|\mathbf{a_{3}}\|=1$, see Fig. 1.1 (left) for the device-specific direction $\mathbf{a_{3}}$. Nevertheless, the Cosserat rod model (2.1) is still not closed, because the surface charge density $\sigma$ is in no situation a constant parameter, but an unknown in the problem. Instead of $\sigma$, we consider the electric current $I$ resulting from convection and conduction, i.e., $I=2\pi a\sigma u+\pi a^{2}\lambda\mathbf{E}\cdot\mathbf{d_{3}}$ with the jet’s conductivity $\lambda$. In the intended transition to stationarity the current $I$ as well as the flow rate $Q=Au$ become constant. Whereas the flow rate is prescribed at the nozzle, the current is still unknown but can be measured in experiments. In literature different phenomenological relations of the form $I\sim Q^{m}\Phi^{n}$, $m,n\in\mathbb{R}$ are documented, see e.g. [6, 11, 27, 30, 25]. The powers vary with the material properties and device geometry, the proportionality constant crucially depends on the experimental conditions.

The dimensionless air drag function

[TABLE]

is expressed in terms of the tangential $w_{\tau}=\mathbf{w}\cdot\boldsymbol{\tau}$ and normal velocity components $w_{\nu}=\sqrt{\mathbf{w}\cdot\mathbf{w}-w_{\tau}^{2}}$ with normal vector $\boldsymbol{\nu}=(\mathbf{w}-w_{\tau}\boldsymbol{\tau})/w_{\nu}$. We particularly use the regularized air resistance coefficients $r_{\nu}$, $r_{\tau}$ given in [19]

[TABLE]

with transition points w_{0}=2\exp\big{(}2.0022-{4\pi}/{r_{\nu}^{S}}\big{)}, $w_{1}=0.1$, $w_{2}=100$, the Stokes limits

[TABLE]

the function $S(w_{\nu})=2.0022-\log(w_{\nu})$ and the regularization parameter $\delta=3.5\cdot 10^{-2}$. The other parameters $p_{k,j}$ and $q_{k,j}$ ($k\in\{\nu,\tau\}$, $j\in\{0,1,2,3\}$) ensure smoothness and are

[TABLE]

[TABLE]

2.2. Jet’s whipping

In the whipping regime the electrified visco-capillary jet forms a helical structure, see for example the experiments in [11, 32, 25] and Fig. 1.1 (middle). In previous works stability analysis has been performed using classical perturbation theory [11, 21, 32, 17]. We also find a strong connection between the whipping instability observed in experiments and the unstable solutions of our electrospinning model (2.1), for details we refer to the temporal stability analysis in App. A. However, our main idea in this paper is different. To explore the instability numerically, we formulate the jet’s whipping as the stationary solution of the Cosserat rod model (2.1). For this purpose we consider a spun jet of certain – a priori unknown – length $L$ with stress-free end. At the nozzle it is straight. Because of the electric field $\mathbf{E}=E\mathbf{a_{3}}$, we have a fixed predominant direction in the device, i.e., $\mathbf{a_{3}}=\mathbf{d_{3}}(0,t)$ is the jet tangent at the nozzle $s=0$ for all times $t$. We introduce a time-dependent outer basis $\{\mathbf{a_{1}}(t),\mathbf{a_{2}}(t),\mathbf{a_{3}}\}$, $\partial_{t}\mathbf{a_{i}}=\boldsymbol{\Omega}\times\mathbf{a_{i}}$ that rotates with the jet’s – a priori unknown – whipping frequency $\Omega$, $\boldsymbol{\Omega}=\Omega\mathbf{a_{3}}$, $\Omega\in\mathbb{R}$ (see Fig. 1.1, left). In addition, the director triad for the jet is modified by incorporating the respective spin in order to get stationary boundary conditions [22], it becomes

[TABLE]

A representation in these director and outer bases eliminates the time-dependencies and yields a stationary set-up, but it obviously introduces fictitious body forces and couples, such as Coriolis, centrifugal and spin-associated ones, due to inertia in the model equations. The director and outer bases are related by the tensor-valued rotation $\mathbf{R}$, i.e., $\mathbf{R}=\mathbf{a_{i}}\otimes\mathbf{d}_{\mathbf{i}}^{s}$. For any quantity we use the following coordinate terminology:

[TABLE]

with $\mathsf{y}=(y_{1},y_{2},y_{3})\in\mathbb{R}^{3}$ and $\mathsf{\breve{y}}=(\breve{y}_{1},\breve{y}_{2},\breve{y}_{3})\in\mathbb{R}^{3}$ where $\mathsf{y}=\mathsf{R}\cdot\mathsf{\breve{y}}$ and $\mathsf{R}=(R_{ij})=(\mathbf{d}_{\mathbf{i}}^{s}\cdot\mathbf{a_{j}})\in SO(3)$. A similar transformation into a stationary set-up has been performed in the investigation of viscous rope coiling [22]. We emphasize that with this approach we do not artificially insert whipping into the problem. Because the whipping frequency belongs to the solution of the problem, it can also be $\Omega=0$ in certain parameter settings. However, the ansatz for the rotation is only valid if $\mathbf{d_{3}}(0,t)$ is aligned with $\mathbf{a_{3}}$ as shown with mathematical arguments in [26].

Remark 1** (Stationarity).**

The periodic rotation of the system around the symmetry axis $\mathbf{a_{3}}$ allows alternatively also the following approach to obtain stationarity. We can express any scalar $y$ and vector-valued $\mathbf{y}$ variables as $y(s,t)=y(s,0)$ and $\mathbf{y}(s,t)=\mathbf{M}(t)\cdot\mathbf{y}^{\circ}(s)$, where

[TABLE]

represents the rotation tensor with respect to the jet’s whipping frequency $\Omega$ and the fixed reference triad at the nozzle $\mathbf{d}_{\mathbf{i}}^{\circ}(0)$ at $t=0$ with $\mathbf{d}_{\mathbf{3}}^{\circ}(0)=\mathbf{a_{3}}$, here $\mathbf{P}_{0}(0,0)=\mathbf{d}_{\mathbf{1}}^{\circ}(0)\otimes\mathbf{d}_{\mathbf{1}}^{\circ}(0)+\mathbf{d}_{\mathbf{2}}^{\circ}(0)\otimes\mathbf{d}_{\mathbf{2}}^{\circ}(0)$. Hence, we get $\partial_{t}y=0$ and $\partial_{t}\mathbf{y}(s,0)=\Omega\mathbf{a_{3}}\times\mathbf{y}^{\circ}(s)$ and, in particular, $\mathbf{v}^{\circ}-u^{\circ}\mathbf{d}_{\mathbf{3}}^{\circ}=\Omega\mathbf{a_{3}}\times\mathbf{r}^{\circ}$ and $\boldsymbol{\omega}^{\circ}-u^{\circ}\boldsymbol{\kappa}^{\circ}=\Omega(\mathbf{a_{3}}-\mathbf{d}_{3}^{\circ})$ for the linear and angular velocities. In this consideration the relevant two frames are the director triad $\mathbf{d}_{\mathbf{i}}^{\circ}(s)$ and the reference triad $\mathbf{d}_{\mathbf{i}}^{\circ}(0)$.

The derivation of the stationary model equations is analogue to the one for rotational spinning, see [2] for details. In the stationary set-up the flow rate is constant, i.e., $Q=Au=const$, such that $A$ can be replaced by $1/u$ in the equations. In addition, the linear and angular velocities can be expressed in terms of the other variables. Incorporating the viscous material laws leads to a boundary value problem of ordinary differential equations for jet curve $\mathbf{r}$, triad (rotational group) $\mathbf{R}$, curvature $\boldsymbol{\kappa}$, convective speed $u$, contact forces $\mathbf{n}$ and couples $\mathbf{m}$. For the spun jet that leaves the nozzle vertically and straight-lined we have the following geometric and kinematic boundary conditions at the nozzle $(s=0)$ as well as stress-free conditions at the end $(s=L)$

[TABLE]

To determine the further two unknowns of the problem that are the jet length $L$ and the whipping frequency $\Omega$, we impose two additional geometric conditions on the curve’s end point. We prescribe its phase and its height in the device geometry

[TABLE]

In spite of the posed condition on the phase, the choice of the outer basis $\{\mathbf{a_{1}}(t),\mathbf{a_{2}}(t),\mathbf{a_{3}}\}$ still leaves one further degree of freedom, e.g. the sign of $\breve{r}_{2}(L)$. Moreover, the sign of $\Omega$ is free. Thus, four stationary solutions are similarly possible in this model framework: $\Omega\geq 0$ or $\Omega\leq 0$ as well as $\breve{r}_{2}(L)\geq 0$ or $\breve{r}_{2}(L)\leq 0$. The form of the solutions is invariant, the resulting jets only differ in the rotational direction of the whipping and in the position of the end point with respect to the $\mathbf{a_{2}}$-halfspace. We aim for the case: $\Omega\geq 0$ and $\breve{r}_{2}(L)\geq 0$. This implies a positive rotation in the outer frame, a positive Rossby number that characterize the relation between the inertial and rotational forces as well as an end point in the positive $\mathbf{a_{2}}$-halfspace.

For the numerical treatment it is convenient to deal with the dimensionless model equations that are stated in the director basis. Thus, we make the system dimensionless using the three problem-relevant lengths (jet length $L$, distance between nozzle and counter-electrode $H$, nozzle diameter $D$) and the jet speed at the nozzle $U$. We introduce the dimensionless quantities as $\tilde{\mathsf{y}}(\tilde{s})=\mathsf{y}(\bar{s}\tilde{s})/\bar{y}$. The reference values are $\bar{s}=L$, $\bar{r}=H$, $\bar{\kappa}=1/H$, $\bar{u}=U$, $\bar{n}=\pi\mu UD^{2}/(4H)$ and $\bar{m}=\pi\mu UD^{4}/(16H^{2})$. The dimensionless jet length is scaled with the device height, the dimensionless whipping frequency corresponds to the inverse Rossby number, i.e., $\tilde{L}=L/H$ and $\tilde{\Omega}=\Omega H/U$. To keep the notation simple we suppress the label $\tilde{~{}}$ and also use the symbols $L$ and $\Omega$ for the dimensionless unknowns in the following. The dimensionless stationary model for the jet’s whipping regime is then given by

[TABLE]

with the capillary, electric and air resistance forces

[TABLE]

and the boundary conditions

[TABLE]

where $\mathsf{P}_{k}=\mathrm{diag}(1,1,k)$, $k\in\mathbb{R}$ and $\{\mathsf{e_{1}},\mathsf{e_{2}},\mathsf{e_{3}}\}$ is the canonical basis in $\mathbb{R}^{3}$. Moreover we write the air drag function $\mathbf{F}$ with respect to the spin-associated director triad using $\mathsf{F}=(F_{1},F_{2},F_{3})$ with \sum_{i=1}^{3}F_{i}(\tau,\mathsf{w})\mathbf{d}_{\mathbf{i}}^{\mathbf{s}}=\textbf{F}\big{(}\sum_{i=1}^{3}\tau_{i}\mathbf{d}_{\mathbf{i}}^{\mathbf{s}},\sum_{i=1}^{3}w_{i}\mathbf{d}_{\mathbf{i}}^{\mathbf{s}}\big{)}. Then, the jet’s whipping model (2.2) is characterized by eight dimensionless parameters

[TABLE]

that are the Reynolds number $\mathrm{Re}$ as ratio between inertial and viscous forces, the surface-tension associated number $\Gamma$ (scaled inverse Capillary number), the potential-current-associated numbers $\Xi$, $\Theta$ and $\Lambda$, the slenderness ratio $\epsilon$ between nozzle diameter and device height as well as the air drag-associated numbers $\mathrm{M}$ and $\mathrm{Re}_{\star}$. The last is obviously also a Reynolds number but with respect to the airflow quantities.

The conductive effects are represented by $\Lambda$. This characteristic number shows the limitation of our comparatively simple electric force model that assumes a constant external electric field and the current as a parameter. It must hold $\Lambda/4<1$. Otherwise the acting electric forces at $s=0$ point into the nozzle which is not physical. In the following we neglect the term, i.e., $\Lambda=0$, and hence overestimate the electric forces by assuming pure convection.

Remark 2** (Quaternions for rotation).**

To parameterize the rotation $\mathsf{R}\in SO(3)$ we use unit quaternions, i.e., $\mathsf{q}=(q_{0},q_{1},q_{2},q_{3})\in\mathbb{R}^{4}$ with $\|\mathsf{q}\|=1$. Consider

[TABLE]

the equations $L^{-1}\partial_{s}\mathsf{R}=-\mathsf{\kappa}\times\mathsf{R}$ with $\mathsf{R}(0)=\mathsf{P}_{1}$ are replaced in (2.2) by

[TABLE]

3. Continuation-Collocation Method

The derived Cosserat rod model for the jet’s whipping (2.2)-(2.7) is a parametric boundary value problem of ordinary differential equations with 19 variables. The numerical challenge lies in solving the problem for arbitrary parameter settings, which requires suitable initial guesses of the respective solutions. The difficulty of finding a suitable guess was already addressed in [22] for viscous rope coiling and treated there manually. We present here a continuation-collocation method that makes the efficient and robust simulation and automatic navigation through a high-dimensional parameter space possible.

3.1. Collocation scheme

To solve a boundary value problem of the form

[TABLE]

we use a four-stage Lobatto IIIa formula as collocation scheme [9]. It is an implicit Runge-Kutta method. The collocation polynomial provides an once continuously differentiable solution that is fifth-order accurate uniformly in $s\in[0,1]$. Mesh selection and error control are based on a scaled residual and the true error of the continuous solution [12]. Thus, we have

[TABLE]

with collocation points $0=s_{0}<s_{1}<...<s_{N}=1$, mesh size $h_{i}=s_{i}-s_{i-1}$, the abbreviation $\mathsf{y}_{i}=\mathsf{y}(s_{i})$ and the coefficients

[TABLE]

The resulting nonlinear system of $N+1$ equations for $(\mathsf{y}_{i})_{i=0,...,N}$ is solved using a Newton method with analytically prescribed Jacobian. This is a classical approach that is provided in the software MATLAB by the routine bvp5c.m (see www.mathworks.com). Its applicability depends on the convergence of the Newton method that is crucially determined by the initial guess. We aim for adapting the initial guess iteratively by means of a continuation method, solving a sequence of slightly varying boundary value problems.

3.2. Continuation procedure

In the continuation method we embed the boundary value problem of interest (A.2) into a family of problems by introducing a continuation parameter tuple $\mathsf{c}\in[0,1]^{n}$

[TABLE]

Here, $\mathsf{\underline{1}}$ denotes the $n$-dimensional tupel of ones. The functions $\mathsf{f}_{0}$, $\mathsf{g}_{0}$ are chosen in such a way that for $\mathsf{c}=\mathsf{0}$ an analytical solution is known. Given this starting solution, we seek for a sequence of parameter tuples $\mathsf{0}=\mathsf{c}_{0},\mathsf{c}_{1},\ldots,\mathsf{c}_{m}=\mathsf{\underline{1}}$ such that the solution to the respective predecessor boundary value problem provides a good initial guess for the successor. The solution associated to $\mathsf{c}=\mathsf{\underline{1}}$ finally belongs to the original system. By help of the continuation parameters certain terms in the ordinary differential equations can be first excluded, then included. Also the boundary conditions can be varied. The choice of the continuation path decides about failure or success because there are not always existing solutions and several meaningful ways. The core of a robust continuation procedure are the step size control and the choice of a continuation path to navigate through a high-dimensional parameter space.

To explain the used procedures, we consider at first an one-dimensional parameter space $c\in[0,1]$. Proceeding from an initial continuation step size $\Delta c_{0}$, a boundary value problem is always solved twice by using one full step and two half steps. If the full step requires more Newton iterations than both half steps together or $k_{1}$-times more collocation points than the second half step, the continuation step is reduced by a factor $k_{2}$, otherwise it is increased by $k_{2}$ for the further computation. If the Newton method fails, the step size is reduced by a factor $k_{3}$ and the computation is repeated. Certainly, it can happen that no solutions exist for $c>c_{crit}\geq 0$, thus the algorithm for the adaptive step size control has a stopping criterion that is based on a minimal step size. In particular, we use $k_{1}=0.1$, $k_{2}=1.5$, $k_{3}=10$, $\Delta c_{0}=10^{-1}$, $\Delta c_{min}=10^{-14}$.

In a high-dimensional parameter space $\mathsf{c}\in[0,1]^{n}$, $n\gg 1$ there are numerous possible but also impossible continuation paths. We consider an equidistant grid on the parameter space with size $\delta c$ in every dimension and apply a tree search algorithm where the grid points (parameter tuples) are the nodes in the tree. The traversal of the tree is described by a recursion. Because from a given node there is more than one possible next node when changing one parameter, we use a heuristic to order the nodes. We apply a depth-first search to explore a continuation path as far as possible before backtracking when the solution to a boundary value problem cannot be found. The step from one node to the next node is performed by solving the associated boundary value problem either directly or by help of the adaptive step size control described above. To accelerate the algorithm we use a heuristic that prefers the diagonal navigation through the parameter space and register all visited nodes on a forbidden list. As grid size we take $\delta c=10^{-1}$, implying $10^{n}$ nodes in total.

3.3. Model-dependent continuation strategy

In our application we consider a generalized rod model as basis for the continuation. We extend the Cosserat rod equations (2.2)-(2.7) by introducing continuation parameters for the viscous, capillary, electric and air forces $c_{\mathrm{Re}},c_{\Gamma},c_{\Xi},c_{\Theta},c_{\epsilon},c_{\mathrm{M}}\in[0,1]$. This means we replace the Reynolds number $\mathrm{Re}$ by the term $c_{\mathrm{Re}}\mathrm{Re}+(1-c_{\mathrm{Re}})\mathrm{Re}_{0}$, analogously for $\Gamma$, $\Xi$, $\Theta$, $\epsilon$ and $\mathrm{M}$ with starting values $\mathrm{Re}_{0}$, $\Gamma_{0}$, $\Xi_{0}$, $\Theta_{0}$, $\epsilon_{0}$, $\mathrm{M}_{0}$. This system is supplemented with the originally posed boundary conditions

[TABLE]

and the following – slightly modified – ones

[TABLE]

Consequently, we embed the jet model into a family of boundary value problems in a high-dimensional parameter space, $\mathsf{c}=(c_{1},...,c_{7},c_{\mathrm{Re}},c_{\Gamma},c_{\Xi},c_{\Theta},c_{\epsilon},c_{\mathrm{M}})\in[0,1]^{13}$. The choice of the continuation parameters allows the variation of the boundary conditions (jet end with lay-down or without stresses) as well as the exclusion and inclusion of the viscous, capillary, electric and air drag effects that mainly dominate the jet dynamics. We summarize the dimensionless model parameters in the tuple

[TABLE]

The remaining model parameter $\mathrm{Re}_{\star}$ is always directly set as desired and not considered in the continuation.

The starting solution for the continuation that belongs to $\mathsf{c}=\mathsf{0}$ is taken from the study on viscous rope coiling [22], i.e.,

[TABLE]

It describes a (non-coiling) jet having the form of a quarter circle in the absence of inertia, surface tension and outer forces, $\mathrm{Re}_{0}=\Gamma_{0}=\Xi_{0}=\mathrm{M}_{0}=0$. The initialization is consistent to our aim $\Omega\geq 0$ and $\breve{r}_{2}(1)\geq 0$. The corresponding parameter tuple is denoted by $\mathrm{p}_{0}=(\mathrm{Re}_{0},\Gamma_{0},\Xi_{0},\Theta_{0},\epsilon_{0},\mathrm{M}_{0})$ with

[TABLE]

using intermediate parameters $\bar{\Theta}$, $\bar{\epsilon}$. In addition we introduce the tuple $\mathrm{\bar{p}}$ of intermediate parameters that allow for moderate physical effects,

[TABLE]

In particular we choose $\bar{\mathrm{p}}=(0.2,4,2000,400,0.1,0)$. The starting solution (3.2) satisfies the desired stress-free condition $\mathsf{n}(1)=\mathsf{m}(1)=\mathsf{0}$ that belongs to $c_{7}=1$. However, the change of each continuation parameter has different effects on the form of the solution and hence on finding appropriate initial guesses in the continuation. To navigate through the high-dimensional parameter space from $\mathsf{c}=\mathsf{0}$ to $\mathsf{c}=\mathsf{\underline{1}}$ we follow therefore a strategy that consists of four parts:

1. (A)

from $\mathsf{c}=\mathsf{0}$ to $\mathsf{c}^{A}=(1,1,1,1,1,1,0,0,0,0,0,0,0)$:

By changing the boundary conditions we obtain the solution associated to viscous rope coiling with dimensionless model parameters $\mathrm{p}_{0}$ and vanishing linear and angular velocities at the jet end, $c_{i}=1$ for $i=1,...,6$. 2. (B)

from $\mathsf{c}^{A}$ to $\mathsf{c}^{B}=(1,1,1,1,1,1,0,c_{\mathrm{\bar{Re}}},c_{\bar{\Gamma}},c_{\bar{\Xi}},0,0,c_{\mathrm{\bar{M}}})$:

By increasing the parameters $c_{\mathrm{Re}}$, $c_{\Gamma}$, $c_{\Xi}$, $c_{\mathrm{M}}$ to intermediate values we incorporate moderate forming of the viscous, capillary, electric and air drag effects and achieve a solution corresponding to the parameter tuple $\mathrm{\bar{p}}$. 3. (C)

from $\mathsf{c}^{B}$ to $\mathsf{c}^{C}=(1,1,1,1,1,1,1,c_{\mathrm{\bar{Re}}},c_{\bar{\Gamma}},c_{\bar{\Xi}},0,0,c_{\mathrm{\bar{M}}})$:

By changing $c_{7}$ we switch from lay-down to stress-free boundary conditions. 4. (D)

from $\mathsf{c}^{C}$ to $\mathsf{c}=\mathsf{\underline{1}}$:

By the final increase of the parameters $c_{\mathrm{Re}}$, $c_{\Gamma}$, $c_{\Xi}$, $c_{\Theta}$, $c_{\epsilon}$, $c_{\mathrm{M}}$ we incorporate all inertial, capillary, electric and aerodynamic effects as desired.

As mentioned we use a heuristic that seeks for the diagonal continuation path through the parameter space. In Part A we use the tree search algorithm on a grid with $10^{6}$ nodes and compute each step directly without refinement. In Part B it turns out that balancing the physical effects strictly leads to the diagonal path such that the parameter space can be reduced to $c_{\mathrm{Re}}=c_{\Gamma}=c_{\Xi}=c_{\mathrm{M}}$. The corresponding one-dimensional search is performed with the adaptive step size control, the same holds for Part C. In Part D we also use the one-dimensional search by increasing every parameter to the target value successively. A good heuristic is to start with the parameter $c_{\Theta}$ and continue increasing the remaining parameters in the fixed order $c_{\Gamma}$, $c_{\Xi}$, $c_{\mathrm{Re}}$, $c_{\mathrm{M}}$, $c_{\epsilon}$.

Remark 3** (Selective preconditioning).**

Due to boundary layers occurring in the solution to the variables $\mathsf{n}$ and $\mathsf{m}$ at the nozzle, which we analyze later, the system of ordinary differential equations (2.2) is badly scaled with respect to these quantities. Therefore our numerical approach includes scaling of these variables using the component-wise maximum norm. This means we solve the system (2.2) with variables $n_{i}$ and $m_{i}$ ($i=1,2,3$) that are scaled equal to one. This procedure equals a selective preconditioning.

3.4. Performance of continuation algorithm

The forthcoming numerical simulations are performed on an Intel Xeon X5675 CPU (6 cores, 12 threads) and 98 GBytes of RAM. The collocation-continuation algorithm is implemented in MATLAB version R2015b and the solver bvp5c.m is mainly used with the default values. The only modification is that we update the Jacobian every Newton step.

As an example for the performance study of the continuation algorithm we choose a parameter tuple inducing moderate emergence of all involved effects in the numerical solution. In particular we set

[TABLE]

Figure 3.1 shows the jet curve during the continuation procedure. We initialize with the prescribed solution (3.2) forming a quarter circle without any inertia, surface tension and outer forces. After changing the boundary conditions in continuation Part A the jet shows a coiling onto the lay-down plane. Incorporating viscous, capillary and outer forces in Part B causes the onset of the whipping mode due to the Coulomb interactions, whereas the external electric field straightens the jet. After Part C we clearly see a stress-free jet end. Finally, the inclusion of all effects in their full magnitude gives rise to the whipping of the jet after Part D. Considering the computational performance of the algorithm (Tab. 3.1) we see that the preparatory Part A takes a CPU time of 174 sec. Each step (including rejection) requires about 2-4 sec. In view of the huge tree with $10^{6}$ nodes and the fact that the continuation path is by no means the diagonal, this is a spectacularly good performance. Part C is in general not very sensitive: approx. 191 sec. CPU time, the minimal step size is of the order $\min\Delta c\sim\mathcal{O}(10^{-3})$. Part B and D contain the actual continuation and need the highest effort (together 82% of the total computation time). These parts take several minutes, $\min\Delta c\sim\mathcal{O}(10^{-8})$. However, with only around hundred steps in total each (instead of $10^{8}$ for the respective fixed step size), it clearly stresses the efficiency of the presented adaptive step size control.

4. Results

Employing our model framework we perform a parameter study. Special attention is paid to the numerical simulation of the whipping of very thin jets and of jets exposed to high aerodynamic drag, as it is typical for electrospinning applications.

4.1. Influence of dimensionless model parameters

Proceeding from the reference solution to $\mathrm{p}^{ref}$ we investigate the jet characteristics under variation of the single model parameters. We particularly study the jet curve, the elongation of the jet end and the throwing range, see Figs. 4.1, 4.2 and 4.3 respectively. Thereby the elongation is given by the Euclidean norm of the jet’s tangent in material (Lagrangian) description at the jet end, which corresponds to the convective speed $u(1)$. The throwing range characterizes the width of the envelope enclosing the three-dimensional jet curve and is defined as the distance from the jet’s end to the vertical axis, i.e., $\breve{r}_{2}(1)$.

As we see, an increase of the Reynolds number $\mathrm{Re}$ induces a higher number of turning points coming along with higher elongation and higher throwing range, which can be explained as follows: Increasing $\mathrm{Re}$ equals decreasing the viscous forces under constant inertia. This reduces the thinning of the jet and thus leads to a larger jet surface. As a result we expect a higher surface charge, which comes along with higher Coulomb repulsion inducing a more intense bending profile including higher elongations and throwing ranges. In contrast an increase of the inverse capillary number $\Gamma$ reduces the elongation and the throwing range of the jet including a complete narrowing of the jet’s envelope. This takes place due to increasing surface tension forces, which form a counteracting force to the Coulomb repulsion. Furthermore we see that increasing the parameter $\Xi$ induces a more uniform bending due to the stabilizing effect of the outer electric field coming along with a decrease of the Coulomb repulsion. This stabilization procedure is indicated by the formation of a maximal throwing range in the middle of the tested parameter region. Furthermore the elongation reaches an asymptotic bound for large $\Xi$. Varying the parameter $\Theta$ simply means increasing or decreasing the effect of Coulomb repulsion. The effect of stronger Coulomb repulsion is indicated by a higher number of turning points, larger elongation and wider throwing range. Decreasing the typical length scale associated parameter $\epsilon$ leads to a intense increase of turning points. This is clear because following the definition a decrease of this parameter equals an increase of the distance from the nozzle to the ground collector while holding the diameter of the nozzle constant. Therefore the bending jet has more space in $\mathbf{a}_{\mathbf{3}}$-direction to propagate and thus forms more turning points in total. This higher number of turning points leads to a higher elongation of the jet’s end, which reaches a maximum at $\epsilon=2.6\cdot 10^{-2}$ approximately. The throwing range forms a maximum for values in the middle of the chosen parameter region around $\epsilon=5.8\cdot 10^{-2}$. Because an increase of $\mathrm{Re}_{\star}$ simply means a higher air drag acting on the jet and the analogon for decrease, we examine the effect of the surrounding air by varying the parameter $\mathrm{M}$ only. Increasing $\mathrm{M}$ narrows the envelope of the jet considerably and leads to a extensive increase of number of turning points. Considering the curve for the highest value of $\mathrm{M}$ the jet seems to form a helix after a short transient effect. A higher air drag acts as higher counteracting force to the Coulomb repulsion and thus narrows the envelope, which becomes cylindrical as soon as air drag and Coulomb repulsion balance. This effect is also indicated by the rapid drop of elongation and throw range for large air drag parameters until they reach an asymptotic bound.

4.2. Distinct whipping

In electrospinning we are faced with parameter settings involving high whipping frequencies (high numbers of turning points in the jet curve) and large elongations at the jet end. Our stationary model and the proposed numerical scheme are able to handle such situations. However, the required resolution is often computationally very demanding. Hence we make use of the numerical results to formulate an analytical solution for the jet’s ’longtime’ behavior.

Motivated from the parameter study in the previous subsection, we consider a setting with a highly developed whipping, namely $\mathrm{Re}=3.5$, $\Gamma=250$, $\Xi=10^{5}$, $\Theta=170$ and $\mathrm{Re}_{\star}=1.5$. Concerning the remaining two dimensionless model parameters $\epsilon$ and $\mathrm{M}$ we point out: On the one hand our air model neglects any coupling between the jet and the surrounding air, that means we overestimate air drag forces in general. On the other hand as seen previously a high air drag-associated number $\mathrm{M}$ and a small slenderness ratio $\epsilon$ severally induce large numbers of turning points of the jet curve, which can easily require a high numerical resolution and thus lead to an unpractical number of discretization points in our collocation method. Thus we focus our considerations on two examples:

[TABLE]

In Example 1 any air drag-associated effects are neglected, see Fig. 4.4 (top). The jet curve bends strongly and its envelope can nearly be described by a circular right cone. The dimensionless convective velocity $u$ indicating the jet’s elongation describes a rapid growth directly at the nozzle with lower increase afterwards. At the jet end we reach approximately $u(1)=3.4\cdot 10^{4}$ (compared to $u(0)=1$). In Example 2 the additional air drag causes a considerable increase of the number of turning points under decrease of the throwing range, see Fig. 4.4 (bottom). Due to the additional resistance caused by the surrounding air the convective speed remains in more moderate regions compared to the first example. In the components of curvature $\kappa$, contact force $\mathsf{n}$ and couple $\mathsf{m}$ a boundary layer arises directly at the nozzle (Figure 4.5). This boundary layer is not a numerical artefact but an essential characteristic of our solution (see also App. B).

Moreover, the contact force $\mathsf{n}$ and couple $\mathsf{m}$ vanish beyond the boundary layer. This serves as motivation to consider the rod model in terms of vanishing contact forces and couples. Because the resulting solution is expected to describe the jet’s ’longtime’ behavior, we neglect all initial conditions and also leave the jet length $L$ and whipping frequency $\Omega$ arbitrary. Surprisingly, fixing $\mathsf{n}=\mathsf{m}=0$ in our rod model (2.2) leads to an analytically solvable system. A solution thereof is

[TABLE]

with 13 parameters $\mathsf{x}=(a_{0},a_{1},b_{0},b_{1},c_{0},c_{1},d_{3},k_{1},k_{2},k_{3},z,\Omega,L)\in\mathbb{R}^{13}$ being implicitly given by the highly nonlinear system of equations

[TABLE]

with air drag associated normal velocity components

[TABLE]

With the particular choice

[TABLE]

the system (4.2) is approximately fulfilled. The corresponding solution (LABEL:eq:analyticSol) forms a helical jet curve that is approached by the numerical computed jet curve of Example 2, see Fig. 4.6. So in the regime of distinct whipping we are able to give an analytical description of the jet’s ’longtime’ behavior.

Summing up, despite the simplicity of the employed electric force model where the current is considered as a constant parameter and conduction is neglected, the numerical results are convincing and characterize qualitatively well the known whipping effects. However, for quantitative comparisons with experiments we lack respective measurements. Experimental studies in literature investigated intensively the whipping characteristics in dependence of various parameter settings, but the underlying potential-current relations of the used set-ups are mostly not quantitatively documented. An exception is [10, 27] where experiments with a solvent of polyethylene oxide in water were performed. But due to the high conductivity of the solvent, we can not benefit from these results. A combined experimental and numerical study is left to future research.

5. Conclusion

In this paper we showed that the whipping instability observed in experiments with electrified visco-capillary jets can be directly approached with the stable stationary solutions of a Cosserat rod model. Our proposed continuation-collocation method for the respective boundary value problems assures efficient and accurate simulations including the automatic navigation through high-dimensional parameter spaces. We qualitatively verified the jet characteristics (whipping frequency, elongation, throwing range) to changes of the model parameters and showed the computability of parameter settings involving drastic jet thinning. Even analytical solutions for the ’longtime’ behavior of the whipping could be derived. In view of quantitative predictions a combined experimental and numerical study is aimed at in future.

Appendix A Global Linear Stability Analysis

The whipping instability observed in experiments was investigated by means of classical perturbation theory in e.g. [11, 21, 32, 17]. In this context, we point out that the unstable solutions of our transient rod model (2.1) are directly related to the whipping instability as we can conclude from the following global stability analysis.

Proceeding from the time-dependent rod system (2.1) in which the geometric model, the material laws and the capillary force are inserted we replace the kinematic equations for the time derivatives of curve and triad by the compatibility conditions [1], yielding

[TABLE]

with the electric and aerodynamic forces

[TABLE]

We supplement the system with the spinning-associated boundary conditions, i.e., inflow at nozzle ($s=0$) and stress-free at the jet end ($s=L$), keeping in mind that the jet length $L$ is an unknown and requires an additional condition,

[TABLE]

In particular we consider here a stationary outer basis $\{\mathbf{a_{1}^{\circ}},\mathbf{a_{2}^{\circ}},\mathbf{a_{3}^{\circ}}\}$ with predominant spinning direction $\mathbf{a_{3}^{\circ}}=\mathbf{a_{3}}$. To preserve the time-dependencies we formulate (A.1) with respect to the director basis and the stationary outer basis using the coordinate terminology $\mathbf{y}=\sum_{i=1}^{3}\hat{y}_{i}\mathbf{d_{i}}=\sum_{i=1}^{3}\breve{y}_{i}^{\circ}\mathbf{a_{i}}^{\circ}$ with $\hat{\mathsf{y}}=(\hat{y}_{1},\hat{y}_{2},\hat{y}_{3})\in\mathbb{R}^{3}$ and $\breve{\mathsf{y}}^{\circ}=(\breve{y}_{1}^{\circ},\breve{y}_{2}^{\circ},\breve{y}_{3}^{\circ})\in\mathbb{R}^{3}$. The corresponding tensor-valued rotation is denoted by $\hat{\mathbf{R}}=\mathbf{a_{i}}^{\circ}\otimes\mathbf{d_{i}}$ with associated matrix $\hat{\mathsf{R}}$. To the dimensionless quantities $\mathsf{\tilde{y}}(\tilde{s},\tilde{t})=\mathsf{y}(\bar{s}\tilde{s},\bar{t}\tilde{t})/\bar{y}$ of Sec. 2 we introduce the additional reference values $\bar{v}=U$, $\bar{a}=D$, $\bar{\omega}=U/H$ and $\bar{t}=L/U$. Suppressing the label $\tilde{~{}}$, the respective dimensionless system reads

[TABLE]

with the outer forces

[TABLE]

and boundary conditions

[TABLE]

Obviously, the boundary value problem has the general conservation form

[TABLE]

For the temporal stability analysis we assume the solution form

[TABLE]

with $\lambda\in\mathbb{C}$ and $0<\varepsilon\ll 1$. By means of this ansatz and Taylor expansion around $\hat{\mathsf{y}}_{0}$, (A.2) can be split into two systems, i.e.,

[TABLE]

for the stationary solution $\hat{\mathsf{y}}_{0}$ as well as

[TABLE]

for the transient correction $\hat{\mathsf{y}}_{1}$ and $\lambda$, where $\mathsf{H}(s)=\partial_{\hat{\mathsf{y}}}{\mathsf{h}}(\hat{\mathsf{y}}_{0}(s))$, $\mathsf{J}(s)=\partial_{\hat{\mathsf{y}}}{\mathsf{j}}(\hat{\mathsf{y}}_{0}(s))$, $\mathsf{K}(s)=\partial_{\hat{\mathsf{y}}}{\mathsf{k}}(\hat{\mathsf{y}}_{0}(s))$, $\mathsf{A}=\partial_{1}{\mathsf{g}}(\hat{\mathsf{y}}_{0}(0),\hat{\mathsf{y}}_{0}(1))$ and $\mathsf{B}=\partial_{2}{\mathsf{g}}(\hat{\mathsf{y}}_{0}(0),\hat{\mathsf{y}}_{0}(1))$. Thereby, the sign of the real part of $\lambda$ decides about stability or instability of $\hat{\mathsf{y}}_{0}$ for $t\rightarrow\infty$.

A solution of (A.3) is obviously a stationary jet that forms a straight line from the nozzle towards the collector, i.e.,

[TABLE]

with $u_{0}$, $\hat{n}_{0,3}$ being prescribed by the boundary value problem

[TABLE]

We determine $u_{0}$ and $\hat{\mathsf{n}}_{0,3}$ with the continuation-collocation method of Sec. 3 and initial values $u_{0}=1$, $\hat{\mathsf{n}}_{0,3}=0$. In order to solve (A.4) we consider an equidistant grid of mesh size $h$ with grid points $s_{i}$, $i\in\{0,...,N\}$. Integrating the differential equation over the intervals $[s_{i-1},s_{i}]$ we apply the fundamental theorem of calculus and the trapezoidal quadrature rule respectively

[TABLE]

with $\hat{\mathsf{y}}_{1,i}=\hat{\mathsf{y}}_{1}(s_{i})$, $\mathsf{H}_{i}=\mathsf{H}(s_{i})$, $\mathsf{J}_{i}=\mathsf{J}(s_{i})$, and $\mathsf{K}_{i}=\mathsf{K}(s_{i})$. After inclusion of the boundary condition we obtain a generalized eigenvalue problem

[TABLE]

with

[TABLE]

and $N_{v}$ the number of rod unknowns. We compute the solution using the MATLAB routine eig.m. Since $Z$ is singular it is convenient to consider the inverse problem $\tilde{\lambda}\mathsf{W}\cdot\hat{\mathsf{Y}}=\mathsf{Z}\cdot\hat{\mathsf{Y}}$, $\tilde{\lambda}=1/\lambda$ transforming infinite eigenvalues to zero.

The observed whipping effect turns out to be a instability in the stated mathematical sense. Figure A.1 (left) shows the eigenvalue spectrum to the model parameters of Example 2, belonging to a mesh with $N=200$. Obviously there exist (inverse) eigenvalues $\tilde{\lambda}$ with positive real part. The corresponding solution $\hat{\mathsf{y}}_{\varepsilon}$ moves away from the stationary solution $\hat{\mathsf{y}}_{0}$ for infinite times, which is hence unstable. The (discrete) eigenfunction $\breve{r}_{1}^{\circ}$ related to the eigenvalue $\tilde{\lambda}$ with largest real part reveals the formation of the characteristic whipping, see Fig. A.1 (right). We employ the transition from stable to unstable solutions in dependence on the electric force related parameters $\Xi$ and $\Theta$. For the computation of the phase diagram we choose discrete test parameters $\Xi_{i}=1000i$, $i=0,...,100$ and determine the corresponding smallest $\Theta_{i}$, for which a solution of the eigenvalue problem with $N=200$ has an eigenvalue with positive real part. In Fig. A.2 the lower region bounded by the data points and the horizontal axis characterizes the stable solutions, the upper region above the data points corresponds to the unstable solutions. As we have seen in Sec. 4 an increasing outer electric field stabilizes the jet, hence for larger values of $\Xi$ the transition from stable to unstable solutions takes place for larger values of $\Theta$.

Appendix B Boundary Layers

For distinct whipping boundary layers arise in the variables $\kappa$, $\mathsf{n}$ and $\mathsf{m}$ at the nozzle, as seen in Fig. 4.5 for Example 2. To verify that the layers are model-based and that our numerical scheme resolves the effect correctly, we re-compute the results in material (Lagrangian) description.

For the transformation of our model (2.2) from Eulerian to Lagrangian description, we apply the same concepts as in [1]: We introduce a general bijective mapping $\Phi(\cdot,t):(\sigma_{a}(t),\sigma_{b}(t))\rightarrow(s_{a}(t),s_{b}(t))$, $\sigma\mapsto\Phi(\sigma,t)$ that satisfies the relations

[TABLE]

where $t_{in}(\sigma)$ prescribes the time of the material point $\sigma$ entering the steady flow domain $\Phi(\sigma,t)\in[0,L]$. Thus, $\Phi$ depends only on the run time $\zeta(\sigma)=t-t_{in}(\sigma)$ setting $\Phi(\sigma,t)=\hat{\Phi}(t-t_{in}(\sigma))=\hat{\Phi}(\zeta(\sigma))$. To ensure that all physical meanings are preserved in the material description we use the concept of type-$n$-fields, $n\in\mathbb{Z}$, i.e. to an arbitrary field $f$ in Eulerian description the associated Lagrangian field $\tilde{f}$ is related according to

[TABLE]

Thereby, our rod variables $\breve{\mathsf{r}}$, $\mathsf{R}$, $\mathsf{v}$, $\mathsf{\omega}$, $\mathsf{n}$, $\mathsf{m}$, $J$, $A$ are treated as type-[math]-fields and $\mathsf{\kappa}$, $\mathsf{f}$, $j$ type-$1$-fields. The exclusive dependence on the run time $\zeta$ is also valid for an arbitrary Lagrangian field, $\tilde{f}(\sigma,t)=\hat{f}(t-t_{in}(\sigma))=\hat{f}(\zeta(\sigma))$. Using $t_{in}(\sigma)=-\sigma/\tilde{u}(0)=-\sigma/U$ with $U=1$ in the dimensionless formulation we have $\partial_{t}\tilde{f}(\sigma,t)=\partial_{t}\hat{f}(t+\sigma)=\partial_{\zeta}\hat{f}(\zeta)=\partial_{\sigma}\hat{f}(t+\sigma)=\partial_{\sigma}\tilde{f}(\sigma,t)$ and therefore

[TABLE]

for all $\sigma,t$. Applying these concepts, regarding the elongation $\hat{e}(\zeta(\sigma))=\tilde{e}(\sigma,t)=j(\sigma,t)$ and dropping the label $\hat{~{}}$, the model for the jet’s whipping in the Lagrangian description reads

[TABLE]

with the outer forces

[TABLE]

and boundary conditions

[TABLE]

The dimensionless end time is $T=t_{end}U/H$ with $t_{end}$ addressing the unknown time of the material point reaching the collector at height $H$.

The Lagrangian description allows to zoom into the region at the nozzle, where the boundary layers are present for distinct whipping set-ups (see Fig. B.1 (left) for the parameters of Example 2 and $\epsilon=0.17$). Mapping the Lagrangian parameters $\zeta_{i}$ used in the simulation as grid points onto the associated Eulerian arc length parameters $s_{i}=\hat{\Phi}(\zeta_{i})$ by means of the nonlinear transformation $\hat{\Phi}$, we see that both descriptions resolve the boundary layer in a similar good manner, Fig. B.1 (right). The reason why we generally prefer to use the Eulerian description is the desired consideration of high jet thinning which goes along with the development of boundary layers at the jet end in the Lagrangian description. The transformation from spatial to material description shifts layers from the nozzle region to the jet end.

Acknowledgments

J. Rivero-Rodríguez and M. Pérez-Saborid would like to thank Profs. A. Fernández-Nieves, and J. Guerrero-Millán for many helpful discussions. This work has been supported by the Ministry of Science and Innovation of Spain (project DPI 2010-20450-C03-02) and the German Research Foundation (DFG, project 251706852, MA 4526/2-1, WE 2003/4-1).