Elsevier

Journal of Computational Physics

Volume 306, 1 February 2016, Pages 137-147
Journal of Computational Physics

A numerical method to study the dynamics of capillary fluid systems

https://doi.org/10.1016/j.jcp.2015.11.048Get rights and content

Abstract

We propose a numerical approach to study both the nonlinear dynamics and linear stability of capillary fluid systems. In the nonlinear analysis, the time-dependent fluid region is mapped onto a fixed numerical domain through a coordinate transformation. The hydrodynamic equations are spatially discretized with the Chebyshev spectral collocation technique, while an implicit time advancement is performed using second-order backward finite differences. The resulting algebraic equations are solved with the iterative Newton–Raphson technique. The most novel aspect of the method is the fact that the elements of the Jacobian of the discretized system of equations are symbolic functions calculated before running the simulation. These functions are evaluated numerically in the Newton–Raphson iterations to find the solution at each time step, which reduces considerably the computing time. Besides, this numerical procedure can be easily adapted to solve the eigenvalue problem which determines the linear global modes of the capillary system. Therefore, both the nonlinear dynamics and the linear stability analysis can be conducted with essentially the same algorithm. We validate this numerical approach by studying the dynamics of a liquid bridge close to its minimum volume stability limit. The results are virtually the same as those obtained with other methods. The proposed approach proves to be much more computationally efficient than those other methods. Finally, we show the versatility of the method by calculating the linear global modes of a gravitational jet.

Introduction

Capillary flows are not only ubiquitous in nature [1], [2] but also of fundamental importance in varied technological fields ranging from chemical and industrial engineering, pharmacy, biotechnology, food industry, etc. In these flows, the coupling between the surface tension and other dynamical effects, such as fluid inertia, viscosity [3], thermal [4] and solute [5] gradients, or viscoelasticity [6], complicates the system behavior, and gives rise to a number of interesting phenomena. Among them, one can mention the pinching and coalescence of interfaces [7], splashing of drops on solid surfaces [8], Marangoni convection due to temperature [4] or surfactant concentration [5] gradients, elasto-capillary dynamics [6], etc. Because surface tension is the driving force in all these cases, the interface evolution must be calculated very accurately. The development of accurate numerical methods for analyzing capillary flows is of great relevance at both the fundamental and technological levels.

The numerical methods devised to deal with the interface evolution are typically categorized into two major classes: interface-tracking [9] and interface-capturing [10], [11]. In an interface-tracking technique, the mesh follows the interface by advecting with the flow a discrete set of points distributed along the interface. On the contrary, interface-capturing methods make use of fixed spatial domains, and introduce an auxiliary function to determine the interface location. Specifically, this location corresponds to the topological region where that function takes a certain value. The temporal evolution of the auxiliary function is calculated from an advection equation on the Eulerian grid, and the surface tension can be incorporated into the Navier–Stokes equations through the continuum surface force approximation [12]. The Volume-of-Fluid (VoF) method [10] is an interface-capturing technique where the auxiliary function is a scalar field, which makes the interface location suffer from numerical diffusion. To improve the simulation accuracy, use is made of adaptive mesh refinement [13], which accumulates the grid points in the vicinity of the interface. An excellent review of numerical methods for modeling multiphase capillary flows can be found in Ref. [14].

The interface evolution can also be accounted for by introducing boundary fitted coordinates [15], [16]. In this case, one defines a time-dependent coordinate transformation such that the physical space occupied by the fluids is mapped onto a fixed computational domain. Typically, the boundaries are aligned with a constant coordinate line, which simplifies the treatment of boundary conditions. The interface position is explicitly characterized with a contour function whose temporal evolution is calculated in the course of the simulation. In this approach, one avoids the numerical diffusion of the interface and takes no approximation to compute the surface tension force, which makes the method very appropriate for capillarity-dominated flows. However, the coordinate transformation renders the hydrodynamic equations much more complicated, and the simulation cannot surpass topological changes (e.g., the pinching of the interface), where the coordinate transformation cannot be applied.

The methods for the spatial discretization of the hydrodynamic equations are traditionally classified into three main categories: finite difference, finite volume, and finite element schemes [17]. In the latter case, one writes the flow variables as a sum of certain “basis functions” multiplied by coefficients whose values are calculated to satisfy the hydrodynamic equations. Spectral methods [18], [19] are related to the finite element approach. The essential difference between them is the fact that the basis functions in the finite element method vanish except on small subdomains, while they are different from zero over the whole domain in the spectral techniques. Chebyshev polynomials [20] are frequently used in the direction normal to the interface as the basis functions of the spectral method. The Chebyshev spectral collocation technique [21] provides very accurate results when the gradients of the hydrodynamic quantities increase significantly in the vicinity of the interface, because it naturally accumulates the grid points in that region. It can also be applied to the other two axes of the problem, or combined with other discretization methods, such as the (forward/central/backward) finite differences.

The analysis of a fluid-dynamical problem frequently involves two fundamental aspects: (i) the linear stability of perturbations around the base flow, and (ii) the nonlinear evolution of the system from given initial conditions. In the linear stability analysis, one first calculates the base flow from the steady nonlinear hydrodynamic equations, and then the corresponding eigenmodes of the linearized problem [22]. The first step entails the resolution of the same equations as those of the nonlinear evolution problem, but dropping the time derivatives. On the contrary, the realization of the second step requires the numerical integration of a completely different set of (linear) equations. The derivation of such equations may become a lengthy and tedious symbolic calculation, especially when boundary fitted coordinates are used. In addition, their resolution implies the implementation of a numerical algorithm considerably different from that of the nonlinear problem.

Capillary flows have been successfully resolved in varied situations by making use of boundary fitting coordinates and the Chebyshev spectral collocation technique. We can mention, for instance, the nonlinear dynamics of axisymmetric liquid bridges [23], [24], gravitational capillary jets [25], the cone-jet mode of electrospray [26], and the tip streaming formation in an electrified pendant drop [27]. That numerical strategy has also been applied to the linear stability analysis of liquid bridges [24], rivulets [28], coflowing streams [29], shear boundary layers [30], and compound jets [31]. In most of the cases, the numerical results were compared with experiments, showing remarkable agreement. It must be pointed out that two considerably different algorithms were used for the nonlinear analysis and the calculation of the linear modes in the above mentioned problems. The set of (nonlinear) algebraic equations were solved with the Matlab subroutine Fsolve to conduct the nonlinear analysis. To calculate the linear modes, the linearized hydrodynamic equations were derived analytically in the first place, and then the corresponding discrete eigenvalue problem was solved with the Matlab subroutine Eigs.

In this paper, we propose a numerical method to solve 2D (axisymmetric) unsteady capillary flows. The method allows one to calculate both the nonlinear evolution of the system and its linear stability with essentially the same algorithm. The Jacobian of the discretized system of equations are symbolic functions calculated before running the simulation and evaluated numerically in the Newton–Raphson iterations. This novel strategy reduces considerably the computing time. The flexibility and computational efficiency of the proposed method constitute two important advantages over the previous approach, and are obtained without losing the high accuracy proved by the Chebyshev spectral collocation technique in the analysis of capillary systems [23], [24], [25], [26], [27], [28], [29], [30], [31]. As applied to the dynamics of a liquid bridge, our method provides virtually the same results as those given by the two procedures specifically designed to deal with the nonlinear and linear stability problems. In the nonlinear case, the computing time is much shorter than that of the reference technique. We also show the capabilities of the proposed method by calculating the global linear modes of an axisymmetric capillary jet driven by the gravitational force (a problem which has been solved only with the 1D approximation [32]).

The paper is organized as follows. We describe the numerical method in Sections 2 and 3. The algorithm is validated in Section 4 by analyzing the dynamics of a liquid bridge at the minimum volume stability limit. This section also allows us to show how the general procedure is applied to a specific problem. The versatility of the method is illustrated in Section 5, where the linear global modes of a gravitational jet are calculated. The paper closes with some conclusions in Section 6.

Section snippets

Mapping and formulation

Consider a 2D (axisymmetric) unsteady flow with two immiscible phases. A fundamental element of our numerical procedure is the instantaneous mapping of the physical space (x,y) onto a numerical domain (ξ,η) which does not depend on time. For that purpose, one defines the coordinate transformation ξ=ξ(x,y;t) and η=η(x,y;t) so that the interface location corresponds to a fixed line in the numerical domain. The fluid-dynamical problem is re-written in terms of the transformed coordinates asF1(t,ξ,η

Numerical method. Linear stability analysis

The procedure described in Section 2 can also be applied to the calculation of the linear global modes with minor modifications. In this case, one assumes the temporal dependence U=U0+εδueiωt (ε1). Here, U represents any hydrodynamic quantity, U0 and δU stand for the base (steady) solution and the spatial dependence of the perturbation, respectively, while ω is the (complex) eigenfrequency. The goal now is to calculate both the eigenfrequencies ω and the corresponding eigenmodes δU as a

Validation of the numerical method

In this section, we apply the numerical method described in Sections 2 and 3 to study the dynamics of a liquid bridge close to its minimum volume stability limit [24]. We consider both the nonlinear dynamics and the linear stability analysis of this fluid configuration. It must be noted that both problems were studied with other numerical procedures in Ref. [24]. The comparison with experimental results showed remarkable agreement in all the cases analyzed. Therefore, those numerical

Global stability of a gravitational capillary jet

The stability analysis presented in Section 4.3 might be regarded as simple because the base solution is that of equilibrium. In this section, we briefly illustrate the capabilities of our numerical method by calculating the global modes of an axisymmetric capillary jet ejected by a feeding capillary and driven by the gravitational force. The global stability of this system has been examined theoretically in several occasions [38], [39], [40], [32], [41]. However, and to the best of our

Conclusions

In this paper, we proposed a general numerical approach to solve the Navier–Stokes equations and the corresponding boundary conditions for a 2D (axisymmetric) capillary system. The method entails two major advantages: (i) it reduces very significantly the computing time in applying the Newton–Raphson method to solve the nonlinear discrete equations; and (ii) it allows the calculation of both the base flow and its eigenmodes using essentially the same algorithm as that of the nonlinear regime.

Acknowledgements

Partial support from the Ministry of Science and Education, Junta de Extremadura, and Junta de Andalucía (Spain) through Grant Nos. DPI2013-46485, GR10047, and P08-TEP-04128, respectively, is gratefully acknowledged. We thank Javier Rivero for his helpful comments on the analytical capabilities of Matlab.

References (41)

  • H.A. Stone

    Dynamics of drop deformation and breakup in viscous fluids

    Annu. Rev. Fluid Mech.

    (1994)
  • M.F. Schatz et al.

    Experiments on thermocapillary instabilities

    Annu. Rev. Fluid Mech.

    (2001)
  • H.A. Stone

    A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface

    Phys. Fluids A

    (1990)
  • L.E. Rodd et al.

    Capillary break-up rheometry of low-viscosity elastic fluids

    Appl. Rheol.

    (2005)
  • J. Eggers

    Nonlinear dynamics and breakup of free-surface flows

    Rev. Mod. Phys.

    (1997)
  • A. Yarin

    Drop impact dynamics: splashing, spreading, receding, bouncing…

    Annu. Rev. Fluid Mech.

    (2006)
  • M. Worner

    Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications

    Microfluid. Nanofluid.

    (2012)
  • J.F. Thompson et al.

    Numerical Grid Generation. Foundations and Applications

    (1985)
  • J.H. Ferziger et al.

    Computational Methods for Fluid Dynamics

    (2013)
  • G. Karniadakis et al.

    Spectral/hp Element Methods for CFD

    (1999)
  • Cited by (72)

    • A numerical investigation of a liquid bridge solidifying with volume change

      2023, International Journal of Heat and Mass Transfer
    • Effects of geometry in the operation of coaxial electrosprays

      2023, Journal of Aerosol Science
      Citation Excerpt :

      A useful model is presented here to investigate the influence of tubes’ geometry on the coaxial electrospray. As in a previous work (López-Herrera et al., 2020), we use the Taylor–Melcher leaky dielectric model (LDM) and the numerical scheme of Herrada and Montanero (2016). We analyze in detail two effects for a given set of liquid properties: (i) the protruding and retracted position of the inner tube at the exit of the outer tube, and (ii) the geometry of the tubes’ edges (blunt or sharp).

    • Fully coupled interface-tracking model for axisymmetric ferrohydrodynamic flows

      2022, Applied Mathematical Modelling
      Citation Excerpt :

      The properties of the solution in the containers, that employ 10 nm particles, are reported in Table 1. The ferrohydrodynamic simulation framework here introduced extends the capillary model presented in Ref. [38] by considering the magnetic interaction. After verifying that the fluid-dynamic results in the absence of magnetic fields are consistent with previous works, the verification and validation processes focus on the magnetic modules.

    • An efficient, robust and high accuracy framework for direct numerical simulation of 2D and 2D axisymmetric immiscible flow with large property contrast

      2021, Computers and Fluids
      Citation Excerpt :

      The incompressible multiphase flow occurs widely in the chemical and physical processes within numerous industrial and scientific applications, including ship hydrodynamics [1], liquid jets [2,3], geodynamic phenomena [4,5], viscoelastic free surface flows [6,7], and immiscible fluid–fluid displacement in porous media [8–10]. As an indispensable complement to the theoretical study, numerical simulation provides a crucial tool to capture, understand, and control the physics involved in the incompressible multiphase problems [11–15]. Over the past decades, various numerical techniques have been developed in order to improve the accuracy and efficiency of the simulation of multiphase flow [16–22].

    View all citing articles on Scopus
    View full text