3.1. Physical model

This study considers liquid flows through a semi-infinite capillary tube with flow rate $Q$ . A meniscus forms at the bottom of the tube. The meniscus deforms into a so-called Taylor cone under the action of the electric field $\boldsymbol{E}$ between the capillary tube (with an electric potential $V=V_{0}$ called printing voltage) and a substrate (with an electric potential $V=0 \mathrm{V}$ ). The Taylor cone emits a jet which will break up to produce droplets. Considering the symmetry of the physical system, a two-dimensional axisymmetric model can be established with coordinates $\{r,z\}$ , as shown in figure 2(a). The computational domain contains a stainless-steel capillary tube nozzle with an outer diameter of $2R_{o}$ , an inner diameter of $2R_{i}$ and a length of $L$ . An infinite flat plate is perpendicular to the tube and spaced at $z=H$ . The width of the computational domain other than the infinite elemental domain is $W$ . The role of the infinite elemental domain is equivalent to a region with infinite width, which is used to simulate the infinity of the capillary tube and the substrate. It has to be mentioned that the infinite element domain is only functional when solving the electric field.

Figure 2.

(a) Definition sketch for the EHD printing. (b) The mesh (left-hand side) of the computational domain and the phase field variable’s distribution (right-hand side) in the printing process.

Because EHD printing is a typical problem of two-phase flow, there are liquid phase $\Omega _{i}$ and gas phase $\Omega _{a}$ in the computational domain. These two phases are shown with different colours in figure 2(a). The density, dynamic viscosity, permittivity and conductivity of liquid phase and gas phase are $\{\rho _{i}, \mu _{i}, \varepsilon _{i}, \sigma _{i}\}$ and $\{\rho _{a}, \mu _{a}, \varepsilon _{a}, \sigma _{a}\}$ , respectively. The surface tension of the gas–liquid interface is $\gamma$ . It is assumed that the two phases are viscous incompressible Newtonian fluids with constant temperature, and the electrolytes in the liquid phase are fully dissociated in this study.

3.2. Governing equations and boundary conditions

Electrohydrodynamic printing is a typical process coupling flow and electric fields. First, for the dynamic behaviour of fluids, the axisymmetric Navier–Stokes equations for velocity $\boldsymbol{v}=(u, w)$ and pressure $p$ can be written as

(3.1) \begin{align} \nabla \cdot \boldsymbol{v}=0 , \\[-24pt] \nonumber \end{align}

(3.2) \begin{align} \rho \left[\frac{\partial \boldsymbol{v}}{\partial t}+\left(\boldsymbol{v}\cdot \nabla \right)\boldsymbol{v}\right]=-\nabla p+\nabla \cdot \left(\mu \nabla \boldsymbol{v}\right)+\rho \boldsymbol{g}+\boldsymbol{F}_{\gamma }+\boldsymbol{F}_{E} , \\[4pt] \nonumber \end{align}

where $\rho$ , $\mu$ , $\boldsymbol{g}$ and t denote the fluid’s density, fluid’s viscosity, gravitational acceleration pointing towards the z-direction and time. Variables $u$ and $w$ represent the components of velocity $\boldsymbol{v}$ in the r-direction and z-direction. Here $\boldsymbol{F}_{\gamma }$ and $\boldsymbol{F}_{E}$ represent the surface tension and electric volumetric forces (Saville Reference Saville1997; Jacqmin Reference Jacqmin1999), which can be expressed as

(3.3) \begin{align} \boldsymbol{F}_{\gamma }=G\nabla \phi , \\[-24pt] \nonumber \end{align}

(3.4) \begin{align} \boldsymbol{F}_{E}=\nabla \cdot \left[\varepsilon _{0}\varepsilon _{r}\left(\boldsymbol{EE}-\frac{1}{2}\left| \boldsymbol{E}\right| ^{2}\boldsymbol{I}\right)\right]=c\boldsymbol{E}-\frac{1}{2}\varepsilon _{0}\left| \boldsymbol{E}\right| ^{{2}}\nabla \varepsilon _{r} . \\[4pt] \nonumber \end{align}

The new symbols in (3.3) are used to model the two-phase interface, which will be introduced below. In (3.4), $\varepsilon _{r}$ represents the relative permittivity of the fluid, $c$ denotes the density of the space charges. The electric field intensity can be calculated by $\boldsymbol{E}=-\nabla V$ , where the variable V is the electric potential. The PNP equation is introduced to describe the conservation of different ionic species. For the k-species (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018; López-Herrera, Herrada & Gañán-Calvo Reference López-Herrera, Herrada and Gañán-Calvo2023; Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2023),

(3.5) \begin{align} \nabla \cdot \left(\varepsilon _{0}\varepsilon _{r}\nabla V\right)=-\sum _{k}ez_{k}n_{k} , \\[-28pt] \nonumber \end{align}

(3.6) \begin{align} \frac{\partial n_{k}}{\partial t}+\nabla \cdot \left(n_{k}\boldsymbol{v}+\omega _{k}ez_{k}n_{k}\boldsymbol{E}-\omega _{k}k_{B}T\nabla n_{k}\right)=r_{k} , \\[6pt] \nonumber \end{align}

where (3.5) is called the Poisson equation and (3.6) is called the Nernst–Planck equation. In these equations, $e$ is the elementary charge, $z_{k}$ is the valence of the k-species, $n_{k}$ is the concentration which is the number of that particles per unit volume, $\omega _{k}$ is its mobility, $k_{B}$ is the Boltzmann constant, $T$ is the absolute temperature and $r_{k}$ is the rate of production by chemical reactions. Three terms in the bracket of (3.6) represent the ion convection due to fluid flow, the ion migration due to electric field and the ion transportation caused by diffusion. Multiply (3.6) by $ez_{k}$ and sum up, a charge transport equation can be written as

(3.7) \begin{align} \frac{\partial c}{\partial t}+\nabla \cdot \left(c\boldsymbol{v}+\sigma \boldsymbol{E}-\sum _{k}ez_{k}\omega _{k}k_{B}T\nabla n_{k}\right)=\sum _{k}e{z_{k}}r_{k} \end{align}

for the space charge density $c=\sum _{k}ez_{k}n_{k}=-\nabla \cdot (\varepsilon _{0}\varepsilon _{r}\nabla V)$ and fluid’s conductivity $\sigma =\sum _{k}\omega _{k}e^{2}{z_{k}}^{2}n_{k}$ (Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012; López-Herrera et al. Reference López-Herrera, Herrada and Gañán-Calvo2023). As the thermal diffusion is neglected as compared with electromigration and the electrolyte in the liquid is assumed to be fully dissociated, it has $\sum _{k}ez_{k}\omega _{k}k_{B}T\nabla n_{k}\rightarrow 0$ and $r_{k}=0$ (Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012; Schnitzer & Yariv Reference Schnitzer and Yariv2015). Therefore, (3.7) can be expressed as

(3.8) \begin{align} \frac{\partial c}{\partial t}+\nabla \cdot \left(c\boldsymbol{v}+\sigma \boldsymbol{E}\right)=0 . \end{align}

The two-phase interface can be captured using the Cahn–Hillard equation (Cahn & Hilliard Reference Cahn and Hilliard1958),

(3.9) \begin{align} \frac{\partial \phi }{\partial t}+\boldsymbol{v}\cdot \nabla \phi =\nabla \cdot \xi \nabla G ,\end{align}

where the fluid phase is described by a phase field variable $\phi$ which varies from $-1$ to 1,

(3.10) \begin{align} \phi = \left\{\begin{array}{c@{\quad}l} 1, & \text{in liquid phase}, \\[4pt]-1\sim 1, & \text{across the gas-liquid interface},\\[4pt]-1, & \text{in air phase}, \end{array}\right. \end{align}

as illustrated in figure 2(b). Constant $\xi$ is the mobility, and $G$ is the chemical potential,

(3.11) \begin{align} G=\lambda \left[-\nabla ^{2}\phi +\frac{\phi \left(\phi ^{2}-1\right)}{\delta ^{2}}\right], \end{align}

where $\delta$ is the two-phase interface thickness’s length scale, $\lambda$ is the mixing energy density which is related with $\delta$ and surface tension $\gamma$ (Yue et al. Reference Yue, Feng, Liu and Shen2004),

(3.12) \begin{align} \lambda =\frac{3\delta \gamma }{\sqrt{8}} . \end{align}

The boundary conditions are summarized in table 1. The form ‘number-number’ in table 1 denotes the straight line between those two numbers, which are the points on the boundary marked in figure 2(b). ‘Type’ in table 1 refers to the boundary condition set in COMSOL Multiphysics software. $\boldsymbol{e}_{r}$ and $\boldsymbol{e}_{z}$ are the unit vectors in the r-direction and z-direction, respectively. Because the infinite elemental domain can only be applied to the electric field, the computational boundaries of the flow field do not contain 1–2, 4–1, 4–5, 7–8, 12–7, 12–13 or 13–11. As for the electric field, boundaries 5–2, 10–8 and 10–11 are inside the computational domain and do not belong to the computational boundary. The fluid velocity, pressure, electric potential and space charge density within the computational domain are initialized to zero.

Table 1.

Boundary conditions of the computational domain.

3.3. Physical properties

The fluid’s physical properties are set as functions of $\phi$ . The density and viscosity are defined using the weighted arithmetic mean (WAM) (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera, Popinet & Herrada Reference López-Herrera, Popinet and Herrada2011; Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2023),

(3.13) \begin{align} \rho =\rho _{a}\frac{1-\phi }{2}+\rho _{i}\frac{1+\phi }{2}=\rho _{a}+\left(\rho _{i}-\rho _{a}\right)\frac{1+\phi }{2} , \\[-24pt] \nonumber \end{align}

(3.14) \begin{align} \mu =\mu _{a}\frac{1-\phi }{2}+\mu _{i}\frac{1+\phi }{2}=\mu _{a}+\left(\mu _{i}-\mu _{a}\right)\frac{1+\phi }{2} , \\[-9pt] \nonumber \end{align}

where $(1+\phi )/2$ is the volume fraction of the liquid phase. However, if the conductivity and the permittivity are set using WAM, a non-realistic charge diffusion phenomenon will happen (Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2023; Suo et al. Reference Suo, Zhang, Huang, Wang, Jia, Yang, Zhang, Li, Tu and Song2023; Lu et al. Reference Lu, Zhong, Chong, Yang, Yue, Zhang and Li2024).

To prevent this phenomenon, a transition function $f(\phi )$ for $\sigma$ and fluid’s permittivity $\varepsilon$ are proposed,

(3.15) \begin{align} \sigma =\sigma _{a}+\left(\sigma _{i}-\sigma _{a}\right)f\left(\phi \right) , \\[-24pt] \nonumber \end{align}

(3.16) \begin{align} \varepsilon =\varepsilon _{0}\varepsilon _{r}=\varepsilon _{a}+\left(\varepsilon _{i}-\varepsilon _{a}\right)f\left(\phi \right) , \\[4pt] \nonumber \end{align}

where $f(\phi )$ is a piecewise function with continuous first-order derivative to ensure the continuity of $\nabla \varepsilon _{r}$ in the whole computational domain,

(3.17) \begin{align} f\left(\phi \right)= \left \{ \begin{array}{c@{\quad}l} 0, & \dfrac{1+\phi }{2}\leq \phi _{s}-\delta _{s},\\[3pt]f_{1}\left(\phi \right), & \phi _{s}-\delta _{s}\lt \dfrac{1+\phi }{2}\lt \phi _{s}+\delta _{s},\\[3pt] 1, & \dfrac{1+\phi }{2}\geq \phi _{s}+\delta _{s}, \end{array} \right. \end{align}

where $0\lt \phi _{s}\lt 1$ is the value of liquid’s volume fraction with $f_{1}(2\phi _{s}-1) = 0.5$ , and $\delta _{s}$ is half the range of the liquid’s volume fraction when $0 \lt f(\phi ) \lt 1$ . The expression of $f_{1}(\phi )$ is defined as

(3.18) \begin{align} f_{1}\left(\phi \right)=\frac{1}{2}+\frac{15}{16}\times \frac{\dfrac{1+\phi }{2}-\phi _{s}}{\delta _{s}}-\frac{10}{16} \left(\frac{\dfrac{1+\phi }{2}-\phi _{s}}{\delta _{s}}\right)^{3}+\frac{3}{16} \left(\frac{\dfrac{1+\phi }{2}-\phi _{s}}{\delta _{s}}\right)^{5} . \end{align}

The function presented here allows the conductivity and permittivity to vary smoothly across the two-phase interface. Moreover, unlike the functions suggested by other researchers (Huh & Wirz Reference Huh and Wirz2022; Suo et al. Reference Suo, Zhang, Huang, Wang, Jia, Yang, Zhang, Li, Tu and Song2023), the range of variation of fluid’s properties with $\phi$ can be fully customized by adjusting the values of $\phi _{s}$ and $\delta _{s}$ when applying the function proposed in this study. As a result, the fluid’s $\sigma$ and $\varepsilon$ can change rapidly to the value of gas’s $\sigma$ and $\varepsilon$ when the volume fraction of the liquid phase starts to decrease. This makes sure the charge will not transport outside the gas–liquid interface. With the usage of $f(\phi )$ , the leakage of the charge is well suppressed.

3.4. Non-dimensionalization

In order to improve the universality of the research, the physical quantities can be non-dimensionalized. Here, the length scale is $l_{c}=R_{i}$ , the time scale is capillary time $t_{c}=\sqrt{\rho _{i}R_{i}^{3}/\gamma }$ , the pressure scale is capillary pressure $p_{c}=\gamma /R_{i}$ , the electric field scale is $E_{0}=2V_{0}/[R_{i}ln(4H/R_{i})]$ considering Marginean, Nemes & Vertes (Reference Marginean, Nemes and Vertes2006), and the space charge density scale is $\varepsilon _{0}E_{0}/R_{i}$ . All physical quantities are made dimensionless as follows:

(3.19) \begin{align} \begin{array}{c} \!\!\!\left(\hat{r},\hat{z}\right)=\dfrac{\left(r,z\right)}{R_{i}},\ \hat{t}=\dfrac{t}{\sqrt{\rho _{i}R_{i}^{3}/\gamma }},\ \hat{\boldsymbol{v}}=\left(\hat{u},\ \hat{w}\right)=\dfrac{\left(u,w\right)}{R_{i}/\sqrt{\rho _{i}R_{i}^{3}/\gamma }},\ \hat{\rho }=\dfrac{\rho }{\rho _{i}},\ \hat{\mu }=\dfrac{\mu }{\mu _{i}},\\[4pt]\!\!\qquad\quad {\varepsilon_{r} }=\dfrac{\varepsilon }{\varepsilon _{0}} , \hat{\sigma }=\dfrac{\sigma }{\sigma _{i}},\ \hat{p}=\dfrac{p}{\gamma /R_{i}},\ \hat{q}=\dfrac{q}{\varepsilon _{0}E_{0}R_{i}^{2}},\ \hat{c}=\dfrac{c}{\varepsilon _{0}E_{0}/R_{i}},\ \hat{Q}=\dfrac{Q}{R_{i}^{3}/\sqrt{\rho _{i}R_{i}^{3}/\gamma }} ,\\[4pt]\!\!\! \hat{V}=\dfrac{V}{E_{0}R_{i}},\ \hat{\boldsymbol{E}}=\dfrac{\boldsymbol{E}}{E_{0}},\ \hat{\delta }=\dfrac{\delta }{R_{i}},\ \hat{\lambda }=\dfrac{\lambda }{\gamma R_{i}},\ \hat{\xi }=\dfrac{\xi }{\sqrt{R_{i}^{3}/\left(\rho _{i}\gamma \right)}},\ \hat{G}=\dfrac{G}{\gamma /R_{i}} , \end{array} \end{align}

where $\{\hat{\rho },\hat{\mu },{\varepsilon_{r}},\hat{\sigma }\}$ are the physical properties of fluid. The dimensionless forms of the physical quantities are denoted with hats. Moreover, five dimensionless numbers are introduced,

(3.20) \begin{align} \begin{array}{c} Oh=\dfrac{\mu _{i}}{\sqrt{\rho _{i}\gamma R_{i}}},We=\dfrac{\rho _{i}v_{flow}^{2}R_{i}}{\gamma },Ca_{e}=\dfrac{\varepsilon _{0}E_{0}^{2}R_{i}}{\gamma },\\[3pt] \alpha =\dfrac{t_{e}}{t_{c}}=\dfrac{\varepsilon _{i}}{\sigma _{i}\sqrt{\rho _{i}R_{i}^{3}/\gamma }},Bo=\dfrac{\rho _{i}gR_{i}^{2}}{\gamma }, \end{array} \end{align}

where $v_{flow}=Q/\pi R_{i}^{2}$ , $t_{e}=\varepsilon _{i}/\sigma _{i}$ , $g\approx 9.8\, \mathrm{m}\, \mathrm{s}^{-2}$ , are the average velocity of the liquid flowing into the capillary, the charge relaxation time and the value of the gravitational acceleration. Here $Oh$ is known as the Ohnesorge number which is used to measure the interrelationship of viscous forces with inertial forces and surface tension force. The Weber number $We$ measures the interrelationship of inertial force with surface tension force. Here $Ca_{e}$ is the electric capillary number which measures the interrelationship between electric field force and capillary force. The charge relaxation parameter $\alpha$ is defined as the ratio of the charge relaxation time to the capillary time. The Bond number $Bo$ is used to measure the interrelationship between gravitational force and surface tension force. Because the gravity is constant, Bo is set to be unchanged in this work. Finally, the dimensionless form of the governing equations can be written as

(3.21) \begin{align} \begin{array}{c} \hat{\nabla }\cdot \hat{\boldsymbol{v}}=0 \end{array} , \\[-24pt] \nonumber \end{align}

(3.22) \begin{align} \hat{\rho }\left[\frac{\partial \hat{\boldsymbol{v}}}{\partial \hat{t}}+\left(\hat{\boldsymbol{v}}\cdot \hat{\nabla }\right)\hat{\boldsymbol{v}}\right] & = -\hat{\nabla }\hat{p}+Oh\hat{\nabla }\cdot \left(\hat{\mu }\hat{\nabla }\hat{\boldsymbol{v}}\right)+Bo\boldsymbol{e}_{z}+\hat{G}\hat{\nabla }\phi \nonumber \\ & \quad + Ca_{e}\left(\hat{c}\hat{\boldsymbol{E}}-\frac{1}{2}\left| \hat{\boldsymbol{E}}\right| ^{2}\hat{\nabla }\varepsilon _{r}\right), \\[-24pt] \nonumber \end{align}

(3.23) \begin{align} \hat{\nabla }\cdot \left(\varepsilon _{r}\hat{\nabla }\hat{V}\right)=-\hat{c} , \\[-24pt] \nonumber \end{align}

(3.24) \begin{align} \frac{\partial \hat{c}}{\partial \hat{t}}+\hat{\nabla }\cdot \left(\hat{c}\hat{\boldsymbol{v}}\right)+\frac{\varepsilon _{i}}{\varepsilon _{0}\alpha }\hat{\nabla }\cdot \left(\hat{\sigma }\hat{\boldsymbol{E}}\right)=0 , \\[-24pt] \nonumber \end{align}

(3.25) \begin{align} \frac{\partial \phi }{\partial \hat{t}}+\hat{\boldsymbol{v}}\cdot \hat{\nabla }\phi =\hat{\nabla }\cdot \left(\hat{\xi }\hat{\nabla }\hat{G}\right) , \\[-12pt] \nonumber \end{align}

where $\hat{\nabla }=(\partial /\partial \hat{r},\partial /\partial \hat{z})$ .

3.5. Numerical method

The governing equations (3.21)–(3.25) are solved in COMSOL Multiphysics by using a finite-element method after being converted to weak form. To obtain the weak forms of the equations, test functions are introduced. After applying the Gaussian theorem and incorporating the boundary conditions, the weak forms of (3.21) and (3.22) can be expressed as

(3.26) \begin{align} \hspace{-60pt}-\int _{\Omega }\hat{\nabla }\Gamma _{\hat{p}}\cdot \hat{\boldsymbol{v}}\mathrm{d}\Omega +\int _{S}\Gamma _{\hat{p}}\left(\hat{\boldsymbol{v}}\cdot \boldsymbol{n}\right){\rm d}S=0 , \\[-24pt] \nonumber \end{align}

(3.27) \begin{align} &-\int _{\Omega }\boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}\,{\boldsymbol\cdot}\,\hat{\rho}\!\left[\frac{\partial \hat{\boldsymbol{v}}}{\partial \hat{t}}+(\hat{\boldsymbol{v}}\,\boldsymbol\cdot \hat{\nabla })\hat{\boldsymbol{v}}\!\right]{\rm d}\Omega +\int _{\Omega }\hat{p}\hat{\nabla }\,\boldsymbol\cdot \boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}{\rm d}\Omega \nonumber \\& +Bo\int _{\Omega }\boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}} \,\boldsymbol\cdot\, \boldsymbol {e_z}{\rm d}\Omega+\int _{\Omega } \! \boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}\,\boldsymbol\cdot\,(\hat{G}\hat{\nabla }\phi) {\rm d}\Omega \nonumber \\ & +Ca_{e} \!\int _{\Omega }\boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}\,\boldsymbol\cdot\, \left(\hat{c}\hat{\boldsymbol{E}}-\frac{1}{2}\left| \hat{\boldsymbol{E}}\right| ^{2}\hat{\nabla }\varepsilon _{r}\right)\!{\rm d}\Omega \nonumber \\ & -Oh\int _{\Omega }\hat{\nabla }\boldsymbol\Gamma _{\hat{\boldsymbol{v}}}\;:\; \left(\hat{\mu }\hat{\nabla }\hat{\boldsymbol{v}}\right){\rm d}\Omega \nonumber \\ & -\int _{S}\hat{p} \left(\boldsymbol {n}\,\boldsymbol\cdot\boldsymbol{\Gamma}_{\hat{\boldsymbol{v}}}\right){\rm d}S + Oh\int _{S}\boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}\,\boldsymbol\cdot\,[\boldsymbol{n}\,\boldsymbol\cdot\,(\hat{\mu }\hat{\nabla }\hat{\boldsymbol{v}})]{\rm d}S=0, \\[1pt] \nonumber \end{align}

where $\Omega$ and $S$ denote the whole computational domain and the boundary of $\Omega$ . n is the unit normal vector of S. Here ${\Gamma} _{\hat{p}}$ and $\boldsymbol {\Gamma} _{\hat{\boldsymbol{v}}}$ are the test functions of the pressure and the velocity, respectively. Similarly, the weak form of the Poisson equation for the electric field can be derived as

(3.28) \begin{align} -\int _{\Omega }\hat{\nabla }\Gamma _{\hat{V}}\cdot \left(\varepsilon _{r}\hat{\nabla }\hat{V}\right)\mathrm{d}\Omega +\int _{S}\Gamma _{\hat{V}}\left[\left(\varepsilon _{r}\hat{\nabla }\hat{V}\right)\cdot \boldsymbol{n}\right]{\rm d}S+\int _{\Omega }\Gamma _{\hat{V}}\hat{c}{\rm d}\Omega =0, \end{align}

where $\Gamma _{\hat{V}}$ represents the test function of the electric potential. The equation for charge transportation can be written as

(3.29) \begin{align} &\int _{\Omega }\Gamma _{\hat{c}}\frac{\partial \hat{c}}{\partial \hat{t}}\mathrm{d}\Omega -\int _{\Omega }\hat{\nabla }\Gamma _{\hat{c}}\cdot \left(\hat{c}\hat{\boldsymbol{v}}\right)\mathrm{d}\Omega -\frac{\varepsilon _{i}}{\varepsilon _{0}\alpha }\int _{\Omega }\hat{\nabla }\Gamma _{\hat{c}}\cdot \left(\hat{\sigma }\hat{\boldsymbol{E}}\right)\mathrm{d}\Omega \nonumber \\ & \quad +\int _{S}\Gamma _{\hat{c}}\left(\hat{c}\hat{\boldsymbol{v}}\right)\cdot \boldsymbol{n}{\rm d}S+\frac{\varepsilon _{i}}{\varepsilon _{0}\alpha }\int _{S}\Gamma _{\hat{c}}\left(\hat{\sigma }\hat{\boldsymbol{E}}\right)\cdot \boldsymbol{n}{\rm d}S=0 ,\end{align}

where $\Gamma _{\hat{c}}$ denotes the test function of the space charge density. For the equation of the phase field, its weak form reads

(3.30) \begin{align} \int _{\Omega } \Gamma _{\hat{\phi }}\frac{\partial \phi }{\partial \hat{t}}\mathrm{d}\Omega +\int _{\Omega }\Gamma _{\hat{\phi }}\left(\hat{\boldsymbol{v}}\cdot \hat{\nabla }\phi \right)\mathrm{d}\Omega +\int _{\Omega }\hat{\nabla }\Gamma _{\hat{\phi }}\cdot \left(\hat{\xi }\hat{\nabla }\hat{G}\right){\rm d}\Omega -\int _{S}\Gamma _{\hat{\phi }}\left(\hat{\xi }\hat{\nabla }\hat{G}\cdot \boldsymbol{n}\right){\rm d}S=0, \end{align}

where $\Gamma _{\hat{\phi }}$ is the test function for the phase field variable.

By combining the above equations with the given boundary and initial conditions, the mathematical description of this study is completed. A free triangular finite-element mesh is used to spatially discretize the computational domain. The mesh is refined along the path of the gas–liquid interface, as shown in figure 2(b). The variables $\hat{\boldsymbol{v}}$ , $\hat{p}$ , $\hat{c}$ , $\phi$ are discretized linearly, and $\hat{V}$ is discretized using the quadratic-order element.

The time discretization adopts a first-order or second-order implicit backward differentiation method with adaptive time stepping. All variables are solved using a segregated node with three substeps. First, $\hat{c}$ and $\hat{V}$ are solved using the constant Newton method with a damping factor of 0.8 and a multifrontal massively parallel sparse direct solver (MUMPS). Second, $\hat{\boldsymbol{v}}$ and $\hat{p}$ are solved using the constant Newton method with a damping factor of 0.8 and a parallel sparse direct solver (PARDISO). Finally, $\phi$ is solved using the constant Newton method with a damping factor of 0.8 and PARDISO.

4.1. Parameter independence analysis

The size of the computational domain is determined as $L=5R_{i}$ and $W=5R_{i}$ . It is proved that as the size increases, the space electric field approximately remains unchanged. Other parameters that need to be confirmed are $\delta$ , $\xi$ and the minimum grid size $\delta _{h}$ , which have a significant impact on the conservation of the liquid phase’s mass (Santra et al. Reference Santra, Mandal and Chakraborty2021). Here $\Delta M$ is used to measure the loss of the liquid, which is defined as the ratio of the liquid’s mass at the time $\hat{t}$ to the mass at the initial moment. As shown in figure 3(a), $\Delta M$ is nearly 99.9 % when $\delta _{h}=R_{i}/80$ at the time $\hat{t}=100$ . This is more than enough for the simulation in this study. Therefore, $\delta _{h}$ is fixed to be $R_{i}/80$ in our numerical settings. Here $\delta =2\delta _{h}$ and $\chi =0.1\, \mathrm{m}\,\mathrm{s}\, \mathrm{kg^{-1}}$ are chosen for $\Delta M$ has almost the maximum value at the time $\hat{t}=100$ seen in figure 3(b), where $\chi =\xi /\delta ^{2}$ .

Figure 3.

(a) The changes in $\Delta M$ over time with different minimum grid size $\delta _{h}$ when $\delta =2\delta_h$ and $\chi =0.1\, \rm {m\,s\,kg}^{-1}$ . (b) The relationship between $\Delta M$ and $\chi$ with different $\delta$ at the time $\hat{t}=100$ when $\delta_h = R_i /80$ .

In order to effectively prevent unreal charge leakage, the conductivity and permittivity have to decrease rapidly with the decrease of $\phi$ . Therefore, $\delta _{s}$ is set as $1-\phi _{s}$ . In this way, $f(\phi )$ will decrease immediately as the liquid’s volume fraction decreases and its first-order derivative keeps continuous. The value of $\phi _{s}$ is determined by considering the change of a suspended droplet’s deformation D with relative conductivity $\sigma _{r}$ under the action of an electric field. This relationship, proposed by Taylor et al. (Reference Taylor, McEwan and de Jong1966) and Ajayi (Reference Ajayi1978), was used to validate the accuracy of the numerical results by many scholars (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera et al. Reference López-Herrera, Popinet and Herrada2011; Das, Dalal & Tomar Reference Das, Dalal and Tomar2021; Chen et al. Reference Chen, Wang, An, Yin and Fang2024). Here $D=(A-B)/(A + B)$ is the deformation parameter of the droplet, A and B are the lengths of the droplet along and perpendicular to the electric field. Here $\sigma _{r}$ is the ratio of the droplet’s conductivity to the surrounding medium’s conductivity. Here, the same setting is adopted as Chen et al. (Reference Chen, Wang, An, Yin and Fang2024). As shown in figure 4, the numerical results agree well with the first-order theory by Taylor and the second-order theory by Ajayi. When $\phi _{s}=0.85$ , the charge leakage is non-ignorable. When $\phi _{s}=0.95$ , the change of $f_{1}(\phi )$ is too fast to get a smooth solution. As a result, the values of $\phi _{s}$ and $\delta _{s}$ are selected as 0.9 and 0.1 in our simulation.

Figure 4.

The variation of droplet’s deformation D with relative conductivity $\sigma _{r}$ , where the colourful dots denote results using the numerical method proposed in this study with $\phi _{s}=0.85, 0.9, 0.95, 0.99$ . The solid and dash line represent the theoretical solution by Taylor et al. (Reference Taylor, McEwan and de Jong1966) and Ajayi (Reference Ajayi1978), respectively.

4.2. Comparison with the experiments

The accuracy of the numerical printing process is also verified by comparison with the experimental results of Hijano et al. (Reference Hijano, Loscertales, Ibáñez and Higuera2015). The experimental platform we set up has high precision to generate droplets with small size, high frequency and high speed. However, the printing process is clearer when the droplets are large and slow. Therefore, the experimental results with a large droplet by Hijano et al. (Reference Hijano, Loscertales, Ibáñez and Higuera2015) are used instead to validate the accuracy of the numerical printing process. In the experiment, the microdripping emission was observed with flow rate Q of 2.1 ml h^–1, printing voltage $V_{0}=3.39\, \mathrm{kV}$ , printing height $H=9\,\text{mm}$ , nozzle’s length $L=20\, \mathrm{mm}$ , nozzle’s outer diameter 500 μm and nozzle’s inner diameter 200 μm. The experimental liquid is a mixture of water, ethanol and ethylene glycol with 1.7 mg g^–1 hydrochloric acid. The properties of the experimental liquid are surface tension $\gamma =41.35\pm 0.42\, \mathrm{mN}\, \mathrm{m^{-1}}$ , conductivity $\sigma _{i}=4.13\pm 0.66\, \mathrm{mS}\, \mathrm{cm^{-1}}$ , viscosity $\mu _{i}=2.28\pm 0.17\, \mathrm{mPa} \cdot \mathrm{s}$ and density $\rho _{i}= 1065\, \mathrm{kg}\, \mathrm{m}^{-3}$ . The relative permittivity is taken as the weighted mean value of the mixture’s components with 63.94. The comparison between the simulation and experiment with the same parameters is presented in figure 5. Numerical results show a good agreement with the experimental observation and obtain more tiny droplets or thin jet below the Taylor cone. The simulation with these tiny droplets shows its higher accuracy than the experiment limited by optical resolution. In the following context, more details about continuous EHD printing are discussed based on the numerical calculation. Moreover, the accuracy of the numerical model’s printing frequency is validated by simulating the experimental conditions in figure 1(b), which can be found in Appendix A for details.

Figure 5.

The comparison of the printing process’s temporal evolution between (a) experimental results and (b) the simulations. Reprinted with permission from Hijano et al. (Reference Hijano, Loscertales, Ibáñez and Higuera2015). Copyright 2015 American Physical Society.

Figure 6.

(a) The profiles of the Taylor cone at $\hat{t}=3.25, 5.25, 7.25, 9.25$ . (b) The changes of the contact angle θ _c with time $\hat{t}$ . The case here has the parameters with Oh = 0.05, $\alpha=0.5$ , Ca _e = 4.5, We = 0.002 and H = 15R _i .