4.1. Computation of steady-state solutions in the absence of an electric field

In the absence of the electric field, using that $(x_i')^2+(z_i')^2=1$ and $x_i' x_i''+z_i' z_i''=0$ for $i=1,2$ and introducing $c_0\equiv - p_0/\gamma _{lg}$ , we can rewrite the equations for the interfacial shapes as a system of first-order ordinary differential equations for the following unknowns:

(4.17) \begin{align} \boldsymbol{X}&=\left (X_1,X_2,X_3,X_4 \right )=\left (x_1,x_1',z_1,z_1' \right )\!, \\[-12pt]\nonumber \end{align}

(4.18) \begin{align} \boldsymbol{Y}&=\left (Y_1,Y_2,Y_3,Y_4 \right )=\left (x_2,x_2',z_2,z_2' \right )\!. \end{align}

The system of equations takes the form

(4.19a,b) \begin{align} \displaystyle X_1' &= X_2, \hspace {3.5cm} \displaystyle Y_1' = Y_2, \\[-12pt]\nonumber \end{align}

(4.20a,b) \begin{align} \displaystyle X_2' &= -\frac {\rho g}{\gamma _{lg}} X_3 X_4 - c_0 X_4, \hspace {1.0cm} \displaystyle Y_2' = -\frac {\rho g}{\gamma _{lg}} Y_3 Y_4 - c_0 Y_4, \\[-12pt]\nonumber \end{align}

(4.21a,b) \begin{align} \displaystyle X_3' &= X_4, \hspace {3.5cm} \displaystyle Y_3' = Y_4, \\[-12pt]\nonumber \end{align}

(4.22a,b) \begin{align} \displaystyle X_4' &= \frac {\rho g}{\gamma _{lg}} X_2X_3 + c_0 X_2, \hspace {1.3cm} \displaystyle Y_4' = \frac {\rho g}{\gamma _{lg}} Y_2Y_3 + c_0 Y_2. \end{align}

The boundary conditions for the new system are

(4.23a,b) \begin{align} X_1(0)&=0, \hspace {3.3cm} Y_1(0)=0, \\[-12pt]\nonumber \end{align}

(4.24a,b) \begin{align} X_4(0)&=0, \hspace {3.3cm} Y_4(0)=0, \\[-12pt]\nonumber \end{align}

(4.25a,b) \begin{align} X_1(l_1)&=d+r\left (1-\cos {\alpha _1} \right )\!, \hspace {0.6cm} Y_1(l_2)=d+r\left (1-\cos {\alpha _2} \right )\!, \\[-12pt]\nonumber \end{align}

(4.26a,b) \begin{align} X_3(l_1) &= r \sin {\alpha _1}, \hspace {2.2cm} Y_3(l_2) = -r \sin {\alpha _2}, \\[-12pt]\nonumber \end{align}

(4.27a,b) \begin{align} X_4(l_1) &= \cot {(\alpha _1+\theta _0)}X_2(l_1), \hspace {0.4cm} Y_4(l_2) = -\cot {(\alpha _2+\theta _0)} Y_2(l_2). \end{align}

The constraints for the arclengths $l_1$ and $l_2$ take the following form:

(4.28a,b) \begin{align} \int _0^{l_1} \big(X_2^2+X_4^2\big)^{1/2} \text{d}s =l_1, \quad \int _0^{l_2} \big(Y_2^2+Y_4^2\big)^{1/2} \text{d}s =l_2. \end{align}

We also have the area constraint (4.14), where $A_1$ and $A_2$ take the following form:

(4.29) \begin{align} A_1 &= \int _0^{l_1}X_2(s)X_3(s) \text{d}s-\left ( \frac {r^2 \alpha _1}{2}-\frac {1}{2}r^2 \cos {\alpha _1}\sin {\alpha _1}\right )\!, \\[-12pt]\nonumber \end{align}

(4.30) \begin{align} A_2 &= -\int _0^{l_2}Y_2(s)Y_3(s) \text{d}s-\left ( \frac {r^2 \alpha _2}{2}-\frac {1}{2}r^2 \cos {\alpha _2}\sin {\alpha _2}\right )\!. \end{align}

We note that Cooray et al. (Reference Cooray, Huppert and Neufeld2016) obtained analytical solutions of the YL equations in terms of elliptic integrals of the first and second kinds to describe liquid bridges between horizontal cylinders. Here, however, we obtain the shapes of liquid bridges by adopting the following numerical approach to solve problem (4.19)–(4.28), (4.14): we apply a Newton iteration procedure to determine $c_0$ , $\alpha _1$ , $\alpha _2$ , $l_1$ and $l_2$ . Given initial estimates for these variables, we solve the boundary-value problem (4.19a )–(4.22a ), (4.23a )–(4.26a ) for $\boldsymbol{X}$ and the boundary-value problem (4.19b )–(4.22b ), (4.23b )–(4.26b ) for $\boldsymbol{Y}$ . Then, a further Newton procedure step is applied to correct $c_0$ , $\alpha _1$ , $\alpha _2$ , $l_1$ and $l_2$ , in order to satisfy the conditions (4.27a , b ), (4.28a , b ) and (4.14). The iterations are continued until convergence is achieved. We compared with the analytical solutions of Cooray et al. (Reference Cooray, Huppert and Neufeld2016) in order to validate the accuracy of our numerical implementation.

4.2. Computation of steady-state solutions in the presence of an electric field

In the presence of the electric field, there appear additional non-local terms $\mathcal{E}_1$ and $\mathcal{E}_2$ in the modified YL (or electrified YL, referred to as eYL) equations, (4.6) and (4.7), resulting from the jumps in the Maxwell stress tensor across the liquid–air interfaces. To compute these contributions, we need to solve Laplace’s equations (3.8a , b ) subject to appropriate boundary conditions. This can be done using a boundary-element method. First, using a standard methodology (see e.g. Pozrikidis Reference Pozrikidis2002), we reformulate the problem for the electric field in the form of boundary-integral equations. Due to symmetry with respect to the $z$ -axis, it is sufficient to consider only the right half of the $(x,z)$ -plane, taking into account that the electric potentials must vanish on the $z$ -axis. We denote the halves of the upper and lower liquid–air interfaces $\varSigma _1$ and $\varSigma _2$ belonging to the right-half of the $(x,z)$ -plane by $\sigma _1$ and $\sigma _2$ , respectively, and we denote the wetted and dry parts of the right cylinder by $c_1$ and $c_2$ , respectively.

The boundary-integral equation for the liquid region at a point $\boldsymbol{x}_0=(x_0,z_0)$ lying on $\sigma _1\cup \sigma _2\cup c_1$ takes the form

(4.31) \begin{eqnarray} \displaystyle \frac {1}{2} \varphi _l(\boldsymbol{x}_0)&=&\displaystyle \int _{\sigma _1\cup \sigma _2} \varphi _l(\boldsymbol{x}) \left [\boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot }\boldsymbol{\nabla }\mathcal{G}\left (\boldsymbol{x},\boldsymbol{x}_0\right ) \right ] \; \text{d} l(\boldsymbol{x})+ \frac {V}{2} \int _{c_1} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot }\boldsymbol{\nabla }\mathcal{G}\left (\boldsymbol{x},\boldsymbol{x}_0\right ) \; \text{d} l(\boldsymbol{x})\nonumber \\ &&+\displaystyle \int _{\sigma _1\cup \sigma _2\cup c_1} \mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0) \boldsymbol{E}_{ln}(\boldsymbol{x}) \; \text{d} l(\boldsymbol{x}), \end{eqnarray}

where $\boldsymbol{x}=(x,z)$ is a variable point on the respective boundaries, $\boldsymbol{n}(\boldsymbol{x})$ is the unit normal vector to the respective boundaries pointing into the liquid and $\boldsymbol{E}_{ln}=-\boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _l$ is the normal component of the electric field in the liquid at the respective boundaries. Also, $\mathcal{G} (\boldsymbol{x},\boldsymbol{x}_0 )$ is the Green’s function for Laplace’s equation that vanishes at $x=0$ , and it is given by

(4.32) \begin{equation} \mathcal{G} (\boldsymbol{x},\boldsymbol{x}_0)=-\frac {1}{2\pi } \ln |\boldsymbol{x}-\boldsymbol{x}_0|+\frac {1}{2\pi } \ln |\boldsymbol{x}-\boldsymbol{x}_0^{im}|, \end{equation}

where $\boldsymbol{x}_0^{im}=(-x_0,z_0)$ is the image of the point $\boldsymbol{x}_0$ with respect to the symmetry $z$ -axis.

The boundary-integral equation for the air region at a point $\boldsymbol{x}_0=(x_0,z_0)$ lying on $\sigma _1\cup \sigma _2\cup c_2$ takes the form

(4.33) \begin{eqnarray} \displaystyle \frac {1}{2}\varphi _g(\boldsymbol{x}_0)&=& -\displaystyle \int _{\sigma _1\cup \sigma _2} \varphi _g(\boldsymbol{x}) \left [{\boldsymbol{n}}(\boldsymbol{x}) \boldsymbol{\cdot }\boldsymbol{\nabla }\mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0)\right ] \; \text{d} l (\boldsymbol{x}) +\frac {V}{2} \int _{c_2} \boldsymbol{n}(\boldsymbol{x}) \boldsymbol{\cdot }\boldsymbol{\nabla }\mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0) \; \text{d}l (\boldsymbol{x})\nonumber \\ && - \displaystyle \int _{\sigma _1 \cup \sigma _2} \mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0) \boldsymbol{E}_{gn} \; \text{d} l(\boldsymbol{x}) + \int _{c_2} \mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0) \boldsymbol{E}_{gn}(\boldsymbol{x}) \;\text{d}l (\boldsymbol{x}), \end{eqnarray}

where $\boldsymbol{n}(\boldsymbol{x})$ points into the liquid if $\boldsymbol{x}\in \sigma _1 \cup \sigma _2$ , as before, and $\boldsymbol{n}(\boldsymbol{x})$ points into the air if $\boldsymbol{x}\in c_1$ . Also, $\boldsymbol{E}_{gn}=-\boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\nabla }\varphi _g$ is the normal component of the electric field in the air at the respective boundaries. Note that to obtain these boundary-integral equations we took the advantage of the fact the electric potentials vanish on the $z$ -axis and tend to zero at infinity, i.e. when $(x^2+z^2)\rightarrow \infty$ , and we also used that $\mathcal{G} (\boldsymbol{x},\boldsymbol{x}_0)=0$ when $\boldsymbol{x}$ is on the $z$ -axis.

Due to the boundary conditions, $\varphi _g$ and $\boldsymbol{E}_{gn}$ at $\sigma _1 \cup \sigma _2$ can be expressed in terms of $\varphi _l$ and $\boldsymbol{E}_{ln}$ as $\varphi _g=\varphi _l$ and $\boldsymbol{E}_{gn}=\tilde {\epsilon }_l \boldsymbol{E}_{ln}$ . Hence, we are left with solving the boundary-integral equations for the following unknowns: the electric potential in the liquid, $\varphi _l$ , at $\sigma _1 \cup \sigma _2$ , the normal component of the electric field in the liquid, $\boldsymbol{E}_{ln}$ , at $\sigma _1 \cup \sigma _2 \cup c_1$ and the normal component of the electric field in the air, $\boldsymbol{E}_{gn}$ , at the cylinder $c_2$ . Knowledge of the solution furnishes the electric-field contribution to the eYL equations, (4.6) and (4.7). In these equations, the tangential components of the electric field at the liquid–air interfaces, i.e. $\boldsymbol{E}_{lt}$ at $\sigma _1$ and $\sigma _2$ , are needed. These can be computed by numerically differentiating the electric potential $\varphi _l$ on $\sigma _1$ and $\sigma _2$ along the tangential directions.

More specifically, given guesses for the liquid–air interface shapes, the electric-field problem is solved numerically using a boundary-element method, in which $\sigma _1$ , $\sigma _2$ , $c_1$ and $c_2$ are discretised into a fixed number of straight elements connected by a sequence of nodes. The unknown variables are approximated with constant values over each element. The integral equations (4.31) and (4.33) are enforced at the midpoints of the elements, with the integrals of $\mathcal{G}(\boldsymbol{x},\boldsymbol{x}_0)$ , $\mathcal{G}_x(\boldsymbol{x},\boldsymbol{x}_0)$ and $\mathcal{G}_y(\boldsymbol{x},\boldsymbol{x}_0)$ computed using a Gauss–Legendre quadrature. This leads to a system of linear equations for the values of $\varphi _l$ and $\boldsymbol{E}_{ln}$ on elements of the top and bottom liquid–air interfaces, and the values of $\boldsymbol{E}_{ln}$ and $\boldsymbol{E}_{gn}$ on the elements of the wetted and dry cylinder parts, respectively, which is solved using Gaussian elimination.

The a priori unknown locations of the liquid–air interfaces $\sigma _1$ and $\sigma _2$ (along with the unknown lengths $l_1$ and $l_2$ of $\sigma _1$ and $\sigma _2$ , respectively, and the unknown pressure $p_0$ ) can then be computed numerically by solving a system of nonlinear equations resulting from a discretisation of the eYL equations, (4.6) and (4.7), using finite differences for the derivatives with respect to $s_1$ and $s_2$ and employing boundary conditions (4.8)–(4.12) and discretised versions of the integral constraints (4.13) and (4.14), and with the electric contributions computed using the boundary-element method described above. The nonlinear system is solved using a Newton iteration procedure, with the electric-field contributions recomputed at each iteration step using the boundary-element method. Numerical convergence tests were performed and it was found that discretising $\sigma _1$ and $\sigma _2$ into $40$ elements each, and, using $40$ elements for the wetted part of the cylinder is sufficient to obtain numerically accurate solutions. The number of the elements for the dry cylinder part is chosen to be $\lceil 40\, (2 \pi -\alpha _1-\alpha _2)/(\alpha _1+\alpha _2) \rceil$ , so that the length of the element on the dry part is approximately the same as the length of the element on the wetted part.

Figure 6.

Bifurcation diagrams (a,c) and solution profiles (b,d) for the non-electrified case with $B_e=0$ (a,b) and the electrified case with $B_e=12$ (c,d). See table 2 for the values of all the dimensionless parameters. The bifurcation diagrams show the dependence of the bridge thickness in the middle, $H_b$ , on the area of the bridge cross-section, $A$ . The solid and dashed lines correspond to stable and unstable solutions, respectively. Solution profiles are shown for different values of $A$ . The blue solid lines correspond to $A=1$ (see the blue filled circles in the bifurcation diagrams). The green dotted lines correspond to the turning points in the bifurcation diagrams (indicated by the green filled diamonds), for which $A=10.3$ and $A=12.1$ for the non-electrified and electrified cases, respectively. The red dashed lines correspond to the points of pinching of the lower liquid–air interface (indicated by the red empty circles in the bifurcation diagrams), for which $A=5.6$ and $A=7.3$ for the non-electrified and electrified cases, respectively.

4.3. Numerical steady-state results

We now present steady-state solutions for liquid bridges in both non-electrified and electrified configurations. In particular, we examine how the bridge shape varies with cross-sectional area. To characterise the solutions, we use the dimensionless bridge thickness in the middle, $H_b=[z_1(0)-z_2(0)]/r$ , and construct bifurcation diagrams (using the pseudo-arclength continuation method) showing how $H_b$ depends on the dimensionless cross-sectional area $A=A_0/r^2$ . We compare the non-electrified and electrified cases in figure 6. For these results, the dimensionless half-distance between the rods is $\tilde d =0.5$ , the equilibrium contact angle is $\theta _0=\pi /12$ and the Bond number is $B=0.5$ . Additionally, for the electrified case, we assume that the permittivity ratio is $\tilde {\epsilon }_l=3$ and that the electric Bond number is $B_e=12$ ; see table 2 for all the values of the dimensionless parameters.

Table 2.

Dimensionless parameters for the steady-state numerical simulations presented in figure 6.

Table 2. Long description

A table with five rows and three columns listing dimensionless parameters for steady-state numerical simulations. The columns are labeled Parameter, Symbol, and Value. The parameters listed are Half-distance between the rods with symbol d-tilde and value 0.5, Equilibrium contact angle with symbol theta-naught and value pi/12, Bond number with symbol B and value 0.5, Electric Bond number with symbol B-e and value 12, and Permittivity ratio of the liquid with symbol e-tilde-l and value 3.

Panels (a) and (c) of figure 6 represent bifurcation diagrams showing the dependence of the bridge thickness $H_b$ on the area of the bridge cross-section $A$ for the non-electrified and electrified cases, respectively. The black solid lines indicate stable solutions up to the turning points (represented by the green filled diamonds). For these points $A=10.3$ and $A=12.1$ , for the non-electrified and electrified case, respectively. These points correspond to the maximal cross-sectional bridge areas (or to the points of the maximum trapping capacity, as defined by Cooray et al. Reference Cooray, Huppert and Neufeld2016) for the respective parameter values. The black dashed lines indicate unstable solutions. They start at the points corresponding to the maximal cross-sectional bridge areas and continue backwards up to the points where the solution profile for the lower liquid–air interface pinches. These points are represented by the red empty circles. For these points $A=5.6$ and $A=7.3$ , for the non-electrified and electrified case, respectively. The branches of unstable solutions can be continued beyond these point to smaller $A$ values, but the profiles become self-intersecting and, therefore, not physically relevant. Thus, we terminate the bifurcation diagrams at these points.

Panels (b) and (d) of figure 6 show some solution profiles of the liquid bridge for the non-electrified and electrified cases, respectively. The blue solid lines correspond to $A=1$ . The respective points in the bifurcation diagrams are indicated by the blue filled circles. The green dotted lines correspond to the maximal cross-sectional bridge areas (see the green filled diamonds in the bifurcation diagrams). Finally, the red dashed lines correspond to the bridges at the point of pinching. As has been already mentioned, the respective points in the bifurcation diagrams are indicated by the red empty circles.

The results indicate that electrified liquid bridges can hold more liquid, i.e. the maximum trapping capacity is increased in the presence of an electric field. Also, it can be observed that, for relatively small bridge cross-sectional areas, the effect of the electric field is to pull the liquid bridge upwards and to make the liquid–air interfaces flatter (compare, for example, the blue solid lines for $A=1$ in panels (b) and (d) of figure 6).

4.4. Comparison with experiments

We now compare experimental observations with steady-state solutions of the eYL model for electrified liquid bridges. Specifically, we compare the 2-D cross-sections of 3-D liquid bridges at their midpoints (obtained from the experimental photos) with the corresponding eYL solutions. Note that in the experiments calibration measurements are first performed: the rod diameter is measured in order to accurately extract 2-D rod and liquid bridge cross-sections from the experimental photos. These measured lengths and curve shapes can then be used to determine the cross-sectional area of the liquid bridge and the other relevant dimensionless parameters required for the model.

Table 3.

Physical properties of the $500\,\mathrm{cSt}$ silicone oil used for the experimental results presented in figures 7, 8 and for the DNS results presented in figure 14.

Table 3. Long description

A table with four rows and four columns detailing the physical properties of silicone oil. The columns are labeled Parameter, Symbol, and Value. The parameters listed are Density, Surface tension, Kinematics viscosity, and Permittivity ratio. The values provided are 970 kg m^-3 for Density, 0.0212 N m^-1 for Surface tension, 500 cSt for Kinematics viscosity, and 3 for Permittivity ratio.

Table 4.

Geometrical and dimensionless parameters for the results presented in figures 7, 8 and 14.

Table 4. Long description

The table presents geometrical and dimensionless parameters for non-electrified and electrified scenarios at 600 V, 1000 V, and 4000 V. It has 6 rows and 5 columns. The columns are labeled as Parameter [unit], Non-electrified, 600 V, 1000 V, and 4000 V. The rows are labeled with parameters: r [mm], d [mm], A [1], θ0 [Rad], B [1], and Be [1]. Row 1: r [mm], Non-electrified, 1; 600 V, 1; 1000 V, 1; 4000 V, 1. Row 2: d [mm], Non-electrified, 0.53; 600 V, 0.53; 1000 V, 0.53; 4000 V, 0.53. Row 3: A [1], Non-electrified, 4.07; 600 V, 2.79; 1000 V, 1.90; 4000 V, 1.10. Row 4: θ0 [Rad], Non-electrified, 0.40; 600 V, 0.40; 1000 V, 0.40; 4000 V, 0.40. Row 5: B [1], Non-electrified, 0.45; 600 V, 0.45; 1000 V, 0.45; 4000 V, 0.45. Row 6: Be [1], Non-electrified, 0; 600 V, 0.48; 1000 V, 1.33; 4000 V, 21.25.

Figure 7 shows a comparison between experimentally extracted interface locations and the corresponding eYL solutions for the 2-D cross-sections of $500\,\mathrm{cSt}$ silicone-oil bridges between copper rods. The relevant properties of the oil are given in table 3, and the corresponding dimensionless and geometrical parameters are listed in table 4. Figure 7(a) corresponds to the non-electrified case, while figure 7(b–d) corresponds to cases where the potential differences between the rods are 600, 1000 and $4000\,\mathrm{V}$ , respectively. When the electric field is applied, the liquid bridge moves upward, causing both the upper and lower interfaces to flatten. At the same time, the cross-sectional area decreases due to the stretching of the liquid bridge along the rods. This is a 3-D effect, the detailed analysis of which is left for future investigation. Nonetheless, we observe that the eYL model captures the influence of the electric field very well, showing good agreement with the experimental results. Figure 8 further illustrates this good agreement, where we overlay the eYL model prediction onto a photograph of a liquid bridge for the case of a $600\,\mathrm{V}$ potential difference.

Figure 7.

Experimental interface locations for a $500\,\mathrm{cSt}$ silicone-oil bridge between copper rods compared with the eYL model solutions for different electric potential differences: (a) non-electrified, (b) $600\, \mathrm{V}$ , (c) $1000\, \mathrm{V}$ and (d) $4000\, \mathrm{V}$ . The corresponding geometrical and dimensionless parameters are shown in table 4. The physical parameters for the $500\,\mathrm{cSt}$ silicone oil are given in table 3.

Figure 8.

The solution of the eYL model overlaid onto a photograph of a $500\,\mathrm{cSt}$ silicone-oil bridge between copper rods with a $600\,\mathrm{V}$ applied electric potential difference. The various relevant oil and rod parameters are given in tables 4 and 3.

Figure 9.

The solution of the eYL model overlaid onto a photograph of a $5\,\mathrm{cSt}$ silicone-oil bridge between platinised titanium rods with a $3000\,\mathrm{V}$ applied electric potential difference.

In figure 9, we show the eYL model prediction overlaid on a photograph of a somewhat different situation: that of a silicone oil with much lower viscosity ( $5\, \mathrm{cSt}$ ) between a pair of platinised titanium electrodes, with a potential difference of $3000\,\mathrm{V}$ applied between them. The geometrical and dimensionless parameters for this case are: $r=1\,\mathrm{mm}$ , $d_0 = 0.5411\,\mathrm{mm}$ , $\theta _0=20 ^{\circ }$ , $A=1.0311$ , $B=0.45$ and $B_e=11.95$ . We again observe that the eYL model is in very good agreement with the experimental result.

5.1. Reduced-order model based on the Onsager variational principle

The dynamics of the system in the overdamped regime (i.e. assuming negligible inertia) can be calculated by taking the minimum of the Rayleighian according to the Onsager variational principle, see e.g. Doi (Reference Doi2013), Lopes et al. (Reference Lopes, Thiele and Hazel2017), Wu et al. (Reference Wu, Di, Man and Doi2019). The Rayleighian is

(5.1) \begin{equation} \mathcal{R}= \dot {F}+\varPhi , \end{equation}

where $\dot {F}\equiv \text{d} F/\text{d}t$ is the time derivative of the free energy and $\varPhi$ is the dissipation function. Note that, more generally, the dissipation function and the Rayleighian are both functionals, but in the present case we simplify so that they are functions of just a few state variables, which we introduce below.

We consider the non-electrified case, for which the free energy, $F$ , is given by

(5.2) \begin{equation} F=\mathcal{F}+G-\lambda ({\mathcal A}-A_0). \end{equation}

Here, $\mathcal{F}$ is the interfacial free energy, $G$ is the gravitational potential energy and $\mathcal A$ represents the area of the cross-section of the bridge (which depends on the state of the system), with $\lambda$ denoting a Lagrange multiplier, which is identified with the pressure $p_0$ in (4.6) and (4.7). The last term in (5.2) is introduced to ensure conservation of volume, so that ${\mathcal A}=A_0$ , where $A_0$ is the specified cross-sectional area of the bridge. Note that we do not incorporate the effects of any electric fields into the Onsager model. It would be at this stage, via the free energy, that such contributions could be included, but we defer doing this to future work.

The interfacial free energy is given by

(5.3) \begin{equation} \mathcal{F}= \, 2 \,l_1 \gamma _{lg}+2 \,l_2 \gamma _{lg} + 2\,l_3\, (\gamma _{sl}-\gamma _{sg}), \end{equation}

which is composed of several terms each consisting of an interfacial tension multiplied by the corresponding length of that interface. We recall that $l_1$ and $l_2$ are the half-lengths of the upper and lower liquid–air interfaces, respectively (see figure 10), and $l_3$ is the full length of the wetted part of the right cylinder. The factor $2$ in each of the terms in (5.3) is due to the symmetry assumption, introduced earlier. Also, $\gamma _{lg}$ , $\gamma _{sl}$ and $\gamma _{sg}$ are the liquid–air, solid–liquid and solid–air surface tensions, respectively.

Figure 10.

Two-dimensional sketch of the liquid bridge under the assumption of the circular interfaces indicating various parameters for this geometry.

The gravitational potential energy is

(5.4) \begin{equation} G={\mathcal A} \, \rho \, g\, \bar {z}, \end{equation}

where $\rho$ is the uniform density of the liquid and $g$ is the gravitational acceleration. Furthermore, $\bar {z}$ is the $z$ coordinate of the centre of mass of the liquid bridge, which can be calculated as $ \bar {z}=\int _{\mathcal{W}_l} z \, \rho \, \text{d}A/ {\mathcal A}\, \rho$ , where $\mathcal{W}_l$ denotes the cross-section of the liquid bridge. Therefore, the potential energy is

(5.5) \begin{equation} G= {\mathcal A} \, \rho \, g \, \frac {\int _{\mathcal{W}_l} z \, \rho \, \text{d}A}{{\mathcal A}\,\rho } = \rho \, g\int _{\mathcal{W}_l} z \, \text{d}A. \end{equation}

Note that minimisation of the free energy, $F$ , given by (5.2) with $\mathcal{F}$ and $G$ defined by (5.3) and (5.5), respectively, leads to the equations for the interfacial shapes (4.6) and (4.7) (with $\mathcal{E}_1=0$ and $\mathcal{E}_2=0$ , as we consider the non-electrified case here). In addition, the transversality condition at each of the three-phase contact lines simply implies that the equilibrium contact angle that the liquid makes with the cylinder surface, $\theta _0$ , satisfies Young’s equation

(5.6) \begin{equation} \cos \theta _0 = \frac {(\gamma _{sg}-\gamma _{sl})}{\gamma _{lg}}. \end{equation}

To obtain a reduced-order model for the dynamics, we approximate the shapes of the interfaces by arcs of circles. Note that, in the absence of gravity, the interfaces are indeed arcs of circles, which we can show analytically. Hence, we expect that this assumption remains appropriate for small liquid bridges. Under this assumption, the state of the system is fully determined by four angles: the angles to the horizontal $\alpha _1$ and $\alpha _2$ that determine the points of contact between the right cylinder and the upper and lower interfaces, respectively, and the contact angles $\theta _1$ and $\theta _2$ between the cylinder and the top and bottom interfaces, respectively – see figure 10. So under this assumption, the state variable is

(5.8) \begin{equation} \boldsymbol{\xi } = ( \alpha _1, \alpha _2, \theta _1, \theta _2 )^T. \quad \end{equation}

Computing the arclengths $l_i$ , $i=1,2,3$ , and using Young’s equation (5.6), leads to the following expression for the interfacial free energy:

(5.9) \begin{align} \begin{split} \mathcal{F}(\boldsymbol{\xi })= & \,2 \left (\pi /2-\alpha _1-\theta _1\right ) r_{1} \gamma _{lg} + 2 \left (\pi /2-\alpha _2-\theta _2 \right ) r_{2} \gamma _{lg} \\ &-2 \gamma _{lg} (\cos {\theta _0}) r \alpha _1 -2 \gamma _{lg} (\cos {\theta _0}) r \alpha _2. \end{split} \end{align}

Here, $r_{1}$ and $r_{2}$ are the radii of the upper and lower circular interfaces, given by

(5.10a,b) \begin{align} r_{1}=\frac {r+d-r\cos \alpha _1}{\cos \left ( \theta _1+\alpha _1\right )}, \quad \quad \quad r_{2}=\frac {r+d-r\cos \alpha _2}{\cos \left ( \theta _2+\alpha _2\right )} . \end{align}

The gravitational potential energy under these assumptions takes the form

(5.11) \begin{align} \begin{split} G(\boldsymbol{\xi })&=\rho g \left [ \vphantom{-b_2\left (d_2\sqrt {r_{2}^2-d_2^2}-r_{2}^2\arctan {\frac {d_2\sqrt {r_{2}^2-d_2^2}} {d_2^2-r_{2}^2}}\right ) } \big(b_1^2+r_{1}^2+d^2+2rd \big)d_1-\left (d+r\right )d_1^2 \right. \\ &-b_1\left (d_1\sqrt {r_{1}^2-d_1^2}-r_{1}^2\arctan {\frac {d_1\sqrt {r_{1}^2-d_1^2}}{d_1^2-r_{1}^2}}\right ) \\&-\big (b_{2}^2+r_{2}^2+d^2+2rd\big )d_2+\left (d+r\right )d_2^2 \\& \left. -b_2\left (d_2\sqrt {r_{2}^2-d_2^2}-r_{2}^2\arctan {\frac {d_2\sqrt {r_{2}^2-d_2^2}} {d_2^2-r_{2}^2}}\right ) \right]\!, \end{split} \end{align}

where $d_1$ and $d_2$ are the $x$ coordinates of the three-phase contact points at the right cylinder (on the upper and lower interfaces, respectively) given by

(5.12a,b) \begin{align} d_1=d+r-r \cos {\alpha _1}, \quad d_2=d+r-r \cos {\alpha _2}. \end{align}

Also, $b_1$ and $b_2$ are the $z$ coordinates of the centres of the circles determining the upper and lower liquid–air interfaces, respectively (see figure 10), given by

(5.13a,b) \begin{align} b_1=r_{1} \cos {\beta _1}+r \sin {\alpha _1}, \quad b_2=-r_{2} \cos {\beta _2}-r \sin {\alpha _2}. \end{align}

With the assumption of circular interfaces, the area is fully determined by the four angles (considering $r$ and $d$ as parameters), and it can be shown that

(5.14) \begin{align} \begin{split} {\mathcal A}(\boldsymbol{\xi }) = &\,2 r_{1} r \sin \beta _1 \sin \alpha _1-r^2 \left ( \alpha _1-\sin \alpha _1\cos \alpha _1\right ) \\ &-r_{1}^2\left ( \beta _1-\sin \beta _1\cos \beta _1\right )+2r_{2}r\sin \beta _2\sin \alpha _2 \\ &-r^2\left (\alpha _2-\sin \alpha _2\cos \alpha _2 \right )-r_{2}^2 \left ( \beta _2-\sin \beta _2\cos \beta _2\right )\!, \end{split} \end{align}

where $\beta _1=\pi /2-\theta _1-\alpha _1$ and $\beta _2=\pi /2-\theta _2-\alpha _2$ .

For the energy dissipation function, as usual we use a quadratic approximation (Doi Reference Doi2013)

(5.15) \begin{equation} \varPhi (\dot {\boldsymbol{\xi }})=\frac {1}{2}\, \dot {\boldsymbol{\xi }} \boldsymbol{\cdot }\boldsymbol{D} \boldsymbol{\cdot }\dot {\boldsymbol{\xi }}, \end{equation}

where $\boldsymbol{D}\equiv (d_{ij})_{1\leqslant i,j \leqslant 4}$ is a symmetric positive definite friction coefficient or dissipation matrix.

In order to compute the state evolution equation, we minimise the Rayleighian function with respect to the rates $\dot {\alpha _1}$ , $\dot {\alpha _2}$ , $\dot {\theta _1}$ and $\dot {\theta _2}$ . Using that the rate of change of the free energy is

(5.16) \begin{equation} \dot {F}(\boldsymbol{\xi })=\frac {\partial F}{\partial \alpha _1} \dot {\alpha _1}+ \frac {\partial F}{\partial \alpha _2} \dot {\alpha _2}+ \frac {\partial F}{\partial \theta _1} \dot {\theta _1}+ \frac {\partial F}{\partial \theta _2} \dot {\theta _2}, \end{equation}

the dynamics is given by the following system of coupled equations:

(5.17a) \begin{align} \frac {\partial \mathcal{R}}{\partial \dot {\alpha _1}} &= \frac {\partial F}{\partial \alpha _1}+ d_{11} \dot {\alpha _1}+ d_{12} \dot {\alpha _2} + d_{13} \dot {\theta _1}+ d_{14} \dot {\theta _2} = 0, \\[-12pt]\nonumber \end{align}

(5.17b) \begin{align} \frac {\partial \mathcal{R}}{\partial \dot {\alpha _2}} &= \frac {\partial F}{\partial \alpha _2}+ d_{12} \dot {\alpha _1}+ d_{22} \dot {\alpha _2} + d_{23} \dot {\theta _1}+ d_{24} \dot {\theta _2} = 0, \\[-12pt]\nonumber \end{align}

(5.17c) \begin{align} \frac {\partial \mathcal{R}}{\partial \dot {\theta _1}} &= \frac {\partial F}{\partial \theta _1}+ d_{13} \dot {\alpha _1}+ d_{23} \dot {\alpha _2} + d_{33} \dot {\theta _1}+ d_{34} \dot {\theta _2} = 0, \\[-12pt]\nonumber \end{align}

(5.17d) \begin{align} \frac {\partial \mathcal{R}}{\partial \dot {\theta _2}} &= \frac {\partial F}{\partial \theta _2}+ d_{14} \dot {\alpha _1}+ d_{24} \dot {\alpha _2} + d_{34} \dot {\theta _1}+ d_{44} \dot {\theta _2} = 0. \end{align}

Here, we have exploited the fact that the dissipation matrix $\boldsymbol{D}$ is symmetric.

Using that $(\boldsymbol{\nabla }A)^T \dot {\boldsymbol{\xi }}=0$ (which follows from volume conservation), the dynamics equations can be written as

(5.18) \begin{equation} \dot {\boldsymbol{\xi }}= -\boldsymbol{D}^{-1}(\boldsymbol{\nabla }(\mathcal{F}+G)-\lambda \boldsymbol{\nabla }{\mathcal A}), \end{equation}

where

(5.19) \begin{equation} \boldsymbol{\nabla }=\left (\frac {\partial }{\partial \alpha _1}, \frac {\partial }{\partial \alpha _2}, \frac {\partial }{\partial \theta _1}, \frac {\partial }{\partial \theta _2} \right )^T \!, \end{equation}

and the Lagrange multiplier $\lambda$ is given by

(5.20) \begin{equation} \lambda = \frac {(\boldsymbol{\nabla }{\mathcal A})^T \boldsymbol{D}^{-1} \boldsymbol{\nabla }(\mathcal{F}+G)} { (\boldsymbol{\nabla }{\mathcal A})^T \boldsymbol{D}^{-1} \boldsymbol{\nabla }{\mathcal A}}. \end{equation}

Figure 11.

(a) The bridge thickness in the middle $H_b$ , varying with the area of the bridge cross-section $A$ , for the same parameter values as in figure 6. The black solid and dashed lines indicate stable and unstable solutions, respectively, of the full governing equations, while the red dash-dotted line indicates steady-state solutions of the Onsager model. The blue filled squares correspond to $A=4$ and $A=6$ , and the respective solutions are shown in panels (b) and (c). The green filled diamond corresponds to the maximal trapping capacity $A=10.3$ and the respective solution is shown in panel (d). As in figure 6, the red empty circle corresponds to the unstable solution where the profile for the lower liquid–air interface pinches.

5.2. Steady-state validation of the Onsager model

Before discussing the dynamics predicted by the reduced-order Onsager model, we first compare steady-state solutions of this model with those of the full governing equations for the non-electrified case (i.e. with solutions of the boundary-value problem described in § 4.1). We consider the same values of the dimensionless parameters as in table 2, but with $B_e=0$ (i.e. without electric field). In figure 11(a), we plot the dimensionless liquid bridge thickness $H_b$ versus the dimensionless cross-sectional area $A$ . As before (see figure 6 a), the black solid and dashed lines correspond to stable and unstable solutions, respectively, of the full governing equations. The green filled diamond indicates the turning point at $A=10.3$ , which is the maximal cross-sectional bridge area. The red empty circle corresponds to the point where the solution profile for the lower liquid–air interface pinches. The red dash-dotted line corresponds to steady-state solutions of the reduced-order Onsager model. Unlike for the full model, this line continues without a turning point beyond the maximal trapping capacity $A=10.3$ . However, we see very good agreement between the results from the Onsager model and from the full governing equations for smaller values of $A$ . In fact, the red dash-dotted line is nearly on top of the black solid line for $A\lesssim 4$ . The steady-state solutions for $A=4$ are shown in figure 11(b). For larger values of $A$ , discrepancies between the solutions increase. An example is shown in figure 11(c) for $A=6$ . There are visible differences between the bottom liquid–air interfaces for the full and reduced-order models, while the agreement remains very good between the top liquid–air interfaces for the two models. The differences between the solutions for the two models further increase as $A$ increases. Figure 11(d) compares the solutions at the maximal trapping capacity for the full governing equations, i.e. at $A=10.3$ . The lower interfaces differ significantly and we now also observe noticeable disagreement between the top interfaces too.

We conclude that the model provides good predictions for moderate values of $A$ ( $A\lesssim 4$ ), but becomes less accurate for larger $A$ , where the assumption of circular interface shapes is less well satisfied due to stronger gravity-induced deformations of the interfaces, particularly of the lower one.

Figure 12.

The dynamics of the liquid–air interfaces obtained using the reduced-order Onsager model. The values of the dimensionless parameters are the same as in figure 6 (but with $B_e=0$ ), and the dimensionless cross-sectional area of the bridge is $A=3$ . Panel (a) shows locations of the interfaces at different times. The initial interfaces are shown with the blue dashed lines. The near-equilibrium state at $t=100$ is shown with the red thick solid lines. The green thin solid lines indicate intermediate liquid–air interface positions at $t=0.05$ , $0.1$ , $0.2$ , $0.4$ and $0.8$ . Panel (b) shows time evolution of the energy $\mathcal{F}+G$ . Panel (c) shows time evolution of the upper and the lower contact angles $\theta _1$ and $\theta _2$ and the angles $\alpha _1$ and $\alpha _2$ that determine the upper and lower three-phase contact points. The equilibrium contact angle $\theta _0$ is also indicated. The angles are given in radians. Panel (d) shows time evolution of the bridge thickness $H_b$ .

Figure 13.

The dynamics of the liquid-air interfaces obtained using the reduced-order Onsager model, as a counterpart to the case presented in figure 12, but here with the cross-section area of the bridge $A=4$ .

5.3. Dynamics predicted by the Onsager model

We now consider the evolution of liquid bridges using the reduced-order Onsager model. Examples are shown in figure 12 (for $A=3$ ) and in figure 13 (for $A=4$ ), where the same parameter values as in figure 6(a,b) are used. For these simulations, we set the dissipation matrix, $\boldsymbol{D}$ , to be equal to the identity matrix, with the overall prefactor (a mobility coefficient) being absorbed into the time scale. This coefficient can be approximated by matching with either the experimental results or the results of DNS. We return to this below. The initial locations of the liquid–air interfaces are indicated by the blue dashed lines in panels (a) of figures 12 and 13 and these correspond to the steady-state liquid bridges in the absence of gravity for $A=3$ and $A=4$ , respectively. The red lines correspond to the liquid–air interfaces at time $t=100$ , at which the liquid bridges have practically reached the final equilibrium states. Intermediate liquid–air interface positions at $t=0.05$ , $0.1$ , $0.2$ , $0.4$ and $0.8$ are shown with the green lines in panels (a).

Panels (b) of figures 12 and 13 show the evolution of the sum of the interfacial energy $\mathcal{F}$ and the gravitational potential energy $G$ . They confirm that $\mathcal{F}+G$ monotonically decreases, as expected. Panels (c) of figures 12 and 13 show the evolution of the state variables for the Onsager model, i.e. the evolution of the upper and the lower contact angles $\theta _1$ and $\theta _2$ and the angles $\alpha _1$ and $\alpha _2$ that determine the upper and lower three-phase contact points. The equilibrium contact angle $\theta _0$ is also shown as the horizontal line, and it has been used as the initial value for both $\theta _1$ and $\theta _2$ . It can be observed that both for $A=3$ and for $A=4$ the upper contact angle, $\theta _1$ , initially decreases to a minimum value and then monotonically increases and converges to the values $0.2737$ and $0.2738$ for $A=3$ and $A=4$ , respectively, which are close to the equilibrium contact angle $\theta _0=\pi /12\approx 0.2618$ . The lower contact angle, $\theta _2$ , initially monotonically increases to a maximum value and then converges to the values $0.321$ and $0.421$ for $A=3$ and $A=4$ , respectively, which differ significantly from $\theta _0$ . The fact that the contact angles $\theta _1$ and $\theta _2$ for the equilibrium state of the Onsager model differ from $\theta _0$ can be explained as follows: Although the assumption of circular shapes of the liquid–air interfaces proves to be good for moderately small values of $A$ , it is, strictly speaking, not valid in the presence of gravity. It turns out that the minimum state for the reduced free-energy function is obtained for $\theta _1$ and $\theta _2$ different from $\theta _0$ . However, this state still gives good overall agreement for the liquid–air interfaces with the steady-state solution profiles for the full governing equations for a wide range of $A$ values (for $A\lesssim 4$ , as already mentioned). Note, however, that, while for $A=3$ the upper contact angle remains positive for all $t$ (nearly reaching zero at its minimum point), for $A=4$ the upper contact angle becomes negative for $t\in (0.301,11.248)$ . This is of course non-physical, and indicates some failure of the model. We conclude that although the reduced-order model gives good prediction of steady-state solutions for $A\lesssim 4$ , it is only valid for $A\lesssim 3$ , as far as the dynamics is concerned. To improve the reduced-order model, it must be extended by having more degrees of freedom in the shapes of the liquid–air interfaces. This is left as a topic for future investigation. Finally, we note that panels (d) of figures 12 and 13 show the evolution of the bridge thickness $H_b$ , for $A=3$ and $A=4$ , respectively. In both cases, $H_b$ initially decreases monotonically to its minimum value before converging to the equilibrium state.

5.4. Comparison with direct numerical simulations

To complement our benchmark experimental results in § 2 and further validate the reduced-order Onsager model, which enables a highly efficient exploration of the parameter space, we perform DNS for the non-electrified case within the open-source VOF Basilisk framework (Popinet Reference Popinet2009, Reference Popinet2015). The DNS also allow the dynamics to be examined in regimes where the Onsager model begins to break down, and provide access to details that are difficult to resolve experimentally. The details of the DNS set-up are discussed in Appendix A.

A fully coupled electrohydrodynamic formulation is not included in the DNS. While such simulations have been performed in Gerris/Basilisk and other frameworks (see e.g. 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; Lin et al. Reference Lin, Skjetne and Carlson2012; Esmaeeli & Behjatian Reference Esmaeeli and Behjatian2020), they typically involve simpler geometries. In the present configuration, the presence of curved embedded solids and moving contact lines makes a fully coupled formulation significantly more involved. We therefore restrict our attention to the non-electrified case, leaving the DNS implementation of the fully coupled multi-physics problem for future work.

Figure 14.

Direct numerical simulation snapshots of the liquid bridge evolution (see also associated supplementary movies) for the non-electrified case of a $500\,\mathrm{cSt}$ silicone-oil bridge between copper rods, with relevant parameters given in tables 3 and 4 (see also corresponding experiments in § 4.4). The evolution is shown at dimensionless times $(a)\ t=0.0$ , $(b)\ t=2.5$ , $(c)\ t=5.0$ and $(d)\ t=20.0$ , from the initially prescribed condition matching the equivalent reduced-order model set-up towards the final stable configuration, which is compared with the results from the various other approaches in figure 15.

A typical numerical evolution is illustrated in figure 14, which shows four snapshots of the interface shape of the liquid bridge within the computational domain at the dimensionless times $t=0.0$ , 2.5, 5.0 and 20.0. These capture the key stages of the evolution within the full simulation time scale of $50.0$ time units. Provided the viscosity is sufficiently large, the interfaces move monotonically. The initial condition used here is identical to that used for the reduced-order Onsager model. Almost immediately in its early time evolution, the bridge shape breaks its initial horizontal symmetry, subsequently falling under gravity towards a final state that is depicted in figure 15. This figure also shows the final state obtained from all the other techniques used in our study. When starting with an identical initial condition as the model, we find the lower part of the bridge stops slightly short (by approximately 5 $\%$ of the solid cylinder arclength) of the equivalent results obtained from the alternative techniques. By contrast, choosing instead an initial condition that already matches the final state at the lower contact point location (with an area-matched shape), leads to very good agreement across the entirety of the solution. The discrepancy between the observed final states arises because the contact lines become pinned at slightly different locations due to the effective roughness introduced by discretising curved solid boundaries on a structured grid. While the influence of this numerically introduced roughness can be reduced by increasing the resolution, practical limitations arise due to the computational cost of DNS in this configuration (see also the discussion in Appendix A).

Figure 15.

Comparisons of theoretical approaches with experimental data for the non-electrified case of a $500\,\mathrm{cSt}$ silicone-oil bridge between copper rods. The corresponding geometrical and dimensionless parameters are shown in table 4. The physical parameters for the $500\,\mathrm{cSt}$ silicone oil are given in table 3.

As briefly touched upon above in our discussion of the Onsager model, the viscosity (or equivalently the dissipation matrix $\boldsymbol{D}$ ) sets the overall time scale for the bridge equilibration dynamics. However, when the viscosity is sufficiently low, there is a substantial change in the dynamics. Namely, as the viscosity is decreased, inertial effects become significant, and the system undergoes a transition from overdamped monotonic relaxation dynamics to underdamped oscillatory behaviour that is not captured by the Onsager model. Figure 16 shows the evolution of the position of the lowest point of the bridge $y$ for four different values of the dynamic viscosity, going down to two orders of magnitude lower compared with the experimental base case. We see that as the dynamic viscosity is decreased, the dynamics of $y$ over time begins to exhibit an overshoot and small-amplitude oscillations as it relaxes toward equilibrium, while for even lower viscosities, a large overshoot is observed, accompanied by oscillatory decay towards equilibrium, corresponding to an underdamped dynamics. For the lowest viscosity value, secondary higher-frequency oscillations can also be observed superposed upon the larger-amplitude lower-frequency main oscillations. We defer the understanding of this transition from an overdamped to an underdamped dynamics to a future study.

Figure 16.

Position $y$ of the mid-point of the lower interface over time, for a liquid bridge relaxing towards equilibrium, for several values of the dynamic viscosity, as given in the legend. These are all for a non-electrified $500\,\mathrm{cSt}$ silicone-oil bridge between copper rods. The largest dynamic viscosity value corresponds to the case described previously in figure 15. With decreasing dynamic viscosity, we see that the evolution of $y$ undergoes a transition from monotonic to underdamped oscillatory behaviour.

Figure 17.

Panel (a) shows the variation over time of the bridge thickness $H_b$ , for the non-electrified case of castor oil suspended between copper rods previously presented in figure 3, comparing the DNS result with that from the Onsager model. They both have the same initial condition and equilibrium contact angle $\theta _0 = 40^{\circ }$ . Panel (b) summarises the final liquid bridge thickness results for a range of values of $\theta _0$ , comparing the Onsager model results with those from the DNS (VOF) approach, which we see starts to fail for $\theta _0\lesssim 30^{\circ }$ (indicated using a dashed vertical line).

We must emphasise that the DNS has no fitting parameters, and as such it is a valuable calibration tool for determining the overall magnitude of the dissipation matrix $\boldsymbol{D}$ in the Onsager model. In figure 17(a), we plot the evolution of $H_b$ over time, as obtained from the DNS together with the result from the Onsager model. For the latter, in order to go from dimensionless time to the real time, we are free to adjust the overall magnitude of the dissipation matrix. Doing this, we observe good qualitative agreement. Note that the DNS dynamics exhibits a delayed onset of the decrease in $H_b$ , followed by an indication of stick–slip-type motion. This behaviour is related to how the contact line is treated numerically. In particular, representing curved boundaries on a structured grid introduces a small effective roughness of the solid surfaces, which affects the contact line dynamics. Although this can be reduced by increasing the resolution, it may still affect the dynamics at low velocities. We also note that accurate numerical treatment of moving contact lines in VOF-type methods remains an active area of research (Tavares et al. Reference Tavares, Josserand, Limare, Lopez-Herrera and Popinet2024; Chen et al. Reference Chen, Han, Pan, Fuster and Zaleski2025; Huang et al. Reference Huang, Han, Zhang and Ni2025), even for simpler geometrical configurations. In figure 17(b), we show how the final bridge thicknesses compare between the DNS and the Onsager model, for varying equilibrium contact angle $\theta _0$ . We observe that for the imposed equilibrium contact angle values in the range $30^{\circ } \lt \theta _0 \lt 50^{\circ }$ , there is remarkably good agreement between the two approaches. The DNS approach with embedded boundaries is known to have limits of applicability for very low or very large values of the contact angle. For our problem, we have identified the lower limit to be approximately $30^{\circ }$ across a range of numerical experiments.

The comparisons in this section show that the Onsager model provides an accurate description of the dynamics in the overdamped regime, corresponding to sufficiently large liquid viscosity. Recall also that, for equilibrium configurations, the Onsager model is also in good agreement with the YL model, provided that the liquid volume is small enough (see § 5.2). Taken together, these cross-validations establish the validity of the Onsager model in the overdamped regime for moderate liquid volumes. Extending the analysis to include electric-field effects in the dynamics is left for future work.