A coupled VOF/embedded boundary method to model two-phase flows on arbitrary solid surfaces

One of the fundamental problems in simulating the motion of sharp interfaces between immiscible fluids is a description of the transition that occurs when the interfaces merge and reconnect. It is well known that classical methods involving sharp interfaces fail to describe this type of phenomena. Following some previous work in this area, we suggest a physically motivated regularization of the Euler equations which allows topological transitions to occur smoothly. In this model, the sharp interface is replaced by a narrow transition layer across which the fluids may mix. The model describes a flow of a binary mixture, and the internal structure of the interface is determined by both diffusion and motion. An advantage of our regularization is that it automatically yields a continuous description of surface tension, which can play an important role in topological transitions. An additional scalar field is introduced to describe the concentration of one of the fluid components and the resulting system of equations couples the Euler (or Navier–Stokes) and the Cahn–Hilliard equations. The model takes into account weak non–locality (dispersion) associated with an internal length scale and localized dissipation due to mixing. The non–locality introduces a dimensional surface energy; dissipation is added to handle the loss of regularity of solutions to the sharp interface equations and to provide a mechanism for topological changes. In particular, we study a non–trivial limit when both components are incompressible, the pressure is kinematic but the velocity field is non–solenoidal (quasi–incompressibility). To demonstrate the effects of quasi–incompressibility, we analyse the linear stage of spinodal decomposition in one dimension. We show that when the densities of the fluids are not perfectly matched, the evolution of the concentration field causes fluid motion even if the fluids are inviscid. In the limit of infinitely thin and well–separated interfacial layers, an appropriately scaled quasi–incompressible Euler–Cahn–Hilliard system converges to the classical sharp interface model. In order to investigate the behaviour of the model outside the range of parameters where the sharp interface approximation is sufficient, we consider a simple example of a change of topology and show that the model permits the transition to occur without an associated singularity.

A coupled Volume of Fluid and Immersed Boundary Method for simulating 3D multiphase flows with contact line dynamics in complex geometries

When boundary conditions arising from the usual hydrodynamic assumptions are applied, analyses of dynamic wetting processes lead to a well-known nonintegrable stress singularity at the dynamic contact line, necessitating new ways to model this problem. In this paper, numerical simulations for a set of representative problems are used to explore the possibility of providing material boundary conditions for predictive models of inertialess moving contact line processes. The calculations reveal that up to Capillary number Ca=0.15, the velocity along an arc of radius 10Li (Li is an inner, microscopic length scale) from the dynamic contact line is independent of the macroscopic length scale a for a>103Li, and compares well to the leading order analytical “modulated-wedge” flow field [R. G. Cox, J. Fluid Mech. 168, 169 (1986)] for Capillary number Ca<0.1. Systematic deviations between the numerical and analytical velocity field occur for 0.1<Ca<0.15, caused by the inadequacy of the leading order analytical solution over this range of Ca. Meniscus shapes produced from calculations in a truncated domain, where the modulated-wedge velocity field [R. G. Cox, J. Fluid Mech. 168, 169 (1986)] is used as a boundary condition along an arc of radius R=10−2a from the dynamic contact line, agree well with those using two inner slip models for Ca<0.1, with a breakdown at higher Ca. Computations in a cylindrical geometry reveal the role of azimuthal curvature effects on velocity profiles in the vicinity of dynamic contact lines. These calculations show that over an appropriate range of Ca, the velocity field and the meniscus slope in a geometry-independent region can potentially serve as material boundary conditions for models of processes containing dynamic contact lines.

A comprehensive numerical and experimental investigation on micrometer-sized water droplet impact dynamics and evaporation on an unheated, flat, dry surface is conducted from the standpoint of spray-cooling technology. The axisymmetric time-dependent governing equations of continuity, momentum, energy, and species are solved. Surface tension, wall adhesion effect, gravitational body force, contact line dynamics, and evaporation are accounted for in the governing equations. The explicit volume of fluid (VOF) model with dynamic meshing and variable-time stepping in serial and parallel processors is used to capture the time-dependent liquid-gas interface motion throughout the computational domain. The numerical model includes temperature- and species-dependent thermodynamic and transport properties. The contact line dynamics and the evaporation rate are predicted using Blake's and Schrage's molecular kinetic models, respectively. An extensive grid independence study was conducted. Droplet impingement and evaporation data are acquired with a standard dispensing/imaging system and high-speed photography. The numerical results are compared with measurements reported in the literature for millimeter-size droplets and with current microdroplet experiments in terms of instantaneous droplet shape and temporal spread (R/D(0) or R/R(E)), flatness ratio (H/D(0)), and height (H/H(E)) profiles, as well as temporal volume (inverted A) profile. The Weber numbers (We) for impinging droplets vary from 1.4 to 35.2 at nearly constant Ohnesorge number (Oh) of approximately 0.025-0.029. Both numerical and experimental results show that there is air bubble entrapment due to impingement. Numerical results indicate that Blake's formulation provides better results than the static (SCA) and dynamic contact angle (DCA) approach in terms of temporal evolution of R/D(0) and H/D(0) (especially at the initial stages of spreading) and equilibrium flatness ratio (H(E)/D(0)). Blake's contact line dynamics is dependent on the wetting parameter (K(W)). Both numerical and experimental results suggest that at 4.5 < We < 11.0 the short-time dynamics of microdroplet impingement corresponds to a transition regime between two different spreading regimes (i.e., for We < or = 4.5, impingement is followed by spreading, then contact line pinning and then inertial oscillations, and for We > or = 11.0, impingement is followed by spreading, then recoiling, then contact line pinning and then inertial oscillations). Droplet evaporation can be satisfactorily modeled using the Schrage model, since it predicts both well-defined transient and quasi-steady evaporation stages. The model compares well with measurements in terms of flatness ratio (H/H(E)) before depinning occurs. Toroidal vortices are formed on the droplet surface in the gaseous phase due to buoyancy-induced Rayleigh-Taylor instability that enhances convection.

In two-phase flows, the interface intervening between the two fluid phases intersects the solid wall at the contact line. A classical problem in continuum fluid mechanics is the incompatibility between the moving contact line and the no-slip boundary condition, as the latter leads to a nonintegrable stress singularity. Recently, various diffuse-interface models have been proposed to explain the contact line motion using mechanisms missing from the sharp-interface treatments in fluid mechanics. In one-component two-phase (liquid-gas) systems, the contact line can move through the mass transport across the interface while in two-component (binary) fluids, the contact line can move through diffusive transport across the interface. While these mechanisms alone suffice to remove the stress singularity, the role of fluid slip at solid surface needs to be taken into account as well. In this paper, we apply the diffuse-interface modeling to the study of contact line motion in one-component liquid-gas systems, with the fluid slip fully taken into account. The dynamic van der Waals theory has been presented for one-component fluids, capable of describing the two-phase hydrodynamics involving the liquid-gas transition [A. Onuki, Phys. Rev. E 75, 036304 (2007)]. This theory assumes the local equilibrium condition at the solid surface for density and also the no-slip boundary condition for velocity. We use its hydrodynamic equations to describe the continuum hydrodynamics in the bulk region and derive the more general boundary conditions by introducing additional dissipative processes at the fluid-solid interface. The positive definiteness of entropy production rate is the guiding principle of our derivation. Numerical simulations based on a finite-difference algorithm have been carried out to investigate the dynamic effects of the newly derived boundary conditions, showing that the contact line can move through both phase transition and slip, with their relative contributions determined by a competition between the two coexisting mechanisms in terms of entropy production. At temperatures very close to the critical temperature, the phase transition is the dominant mechanism, for the liquid-gas interface is wide and the density ratio is close to 1. At low temperatures, the slip effect shows up as the slip length is gradually increased. The observed competition can be interpreted by the Onsager principle of minimum entropy production.

Deposition of colloidal material in evaporating droplets is important in many applications including DNA sequencing and medical diagnostic testing. When colloidal droplets evaporate, the majority of material is often deposited at the periphery of the resultant deposition in a coffee-ring pattern. Formation of this pattern is the result of contact line pinning and the interplay between evaporative and surface tension effects in the droplet. When the contact line is pinned and the evaporative flux in the droplet is highest at the periphery, a radially outward flow is generated to conserve mass that deposits particles in the fluid at the contact line. Evaporation at the contact line can also create a temperature gradient across the droplet. This gradient gives rise to a surface tension driven flow that can resuspend particles in the droplet. When the evaporative flow dominates, particles are deposited at the contact line in a coffee-ring pattern. The presence of the coffee-ring pattern is undesirable in many printing and medical diagnostic processes. Suppression of the coffee-ring effect has been achieved by addition of surfactant, enhancement of surface tension flow, surface modification, alteration of particle shape, and application of an electric field. Manipulation of the coffee-ring effect has been achieved through the application of both AC and DC electric fields. One result of the presence of this field is the electrowetting force at the contact line which acts to reduce the contact angle and increase contact area. Since this force acts at the contact line, it may disrupt typical contact line dynamics, including evaporative dynamics, which are responsible for the formation of the coffee-ring effect. This work will experimentally examine contact line dynamics of evaporating droplets in the presence of DC electric fields. Droplets of water will be desiccated on a device with a photolithographically patterned electrode covered with a thin layer of SU8-3005. Experimental cases with applied DC fields will be compared with unactuated control cases to examine changes in transient interface shape and contact diameter.

The stream function solution for the inner region Stokes flow, for a locally plane moving fluid interface near the triple point, is derived considering three different boundary conditions: the Navier slip boundary condition (NBC), the super-slip boundary condition and the generalized Navier boundary condition (GNBC). The NBC, incorporating a slip length parameter λ , is a well-known method for regularization in the context of the three-phase dynamic contact line problem. It is demonstrated that the velocity field solution under this boundary condition maintains a C 0 continuity at the contact line, resulting in a logarithmic divergence of the pressure at the contact line. By contrast, the super-slip boundary condition establishes a proportional relationship between the wall velocity and the normal derivative of the shear stress, leading to a C 1 velocity field. Furthermore, the GNBC, which introduces an uncompensated Young stress to drive the contact line, yields a C 2 velocity field. The dominant terms are explicitly derived, and the analytical approach presented here can be extended to other bi-harmonic problems as well.

Challenges and opportunities in modelling of proton exchange membrane fuel cells (PEMFC)

In aero-engines, it is important to predict the behavior of shear flows in the different parts such as bearing chambers or gearboxes. In bearing chambers, the thickness distribution of wavy films is well studied as two-phase flows are still very hard to predict depending on the case. Experimental studies remain very expensive to carry out and Computational Fluid Dynamics (CFD) still struggles with two-phase flow prediction especially when a sharp interface between the two phases must be modelled. CFD is used to predict the oil film thickness distribution and interface velocity at different engine operating conditions. Currently Reynold-Averaged Navier-Stokes (RANS) CFD uses a semi-empirical method of turbulence damping, which is inaccurate for wavy films and so impacts the modelling of bearing chambers and gearboxes. With the objective of improving RANS models from Large Eddy Simulation (LES) methods, the Volume of Fluid (VOF) and Euler-Euler methods for two-phase flow modelling are investigated in this study. The Volume of Fluid (VOF) approach assumes a single set of momentum equations for the two phases and volume fractions are 1 or 0 everywhere except in the interface region. An alternative to VOF, is the Euler-Euler method with interface sharpening for shear flows. This approach assumes one set of momentum equations per phase but a shared field of pressure. The VOF and Euler-Euler approaches are compared in this study using LES with the CFD code OpenFOAM v6. The case study is based on experimental work investigating stratified flow in a horizontal channel that will be further detailed in the paper. In this study, a simplified 3D periodic channel filled with two distinct phases: air and water is used. A flow regime is studied in which flows are fully developed and the water phase has a much smaller velocity than the air phase in order to obtain a shear flow. Numerical results are compared with experimental measurements from the literature. With OpenFOAM, the VOF solver used for the study is interFoam and the Euler-Euler solver used is reacting-MultiphaseEulerFoam. Velocity profiles, shear-stress profiles and kinetic energy profiles are compared with experimental measurements for the assessment of the two flow solvers. Maps of vorticity magnitude are also provided to support the comparisons between the Euler-Euler and the VOF approaches as well as an appropriate vortex identification method.

With increasing computer capabilities, numerical modeling of two-phase flows has developed significantly over the last few years. Although there are two main categories, namely 'one' fluid and 'two' fluid methods, the 'one' fluid methods are more commonly used for tracking or capturing the interface between two fluids. Level set (LS), volume-of-fluid (VOF), front tracking, marker-and-cell (MAC) and lattice Boltzmann (LB) methods are all 'one' fluid methods. It is clear that there is no perfect method; each method has advantages and disadvantages which make it more appropriate for one kind of problem than for others. For instance, a LS method will accurately compute the curvature and the normal to the interface, but tends to loss mass which is physically incorrect. On the other hand, a VOF method will conserve mass up to machine precision, but the computation of the curvature and normal to the interface is not as accurate. In order to minimize the disadvantages of these methods, several authors have used two or more methods together to model two-phase flows. This is the case for the CLSVOF (Couple Level Set Volume Of Fluid) method, where LS and VOF are coupled together in order to better capture the interface. In CLSVOF, the level set function is used to compute the interface curvature and normal to the interface, while the volume of fluid function is used to capture the interface. For two-phase flows in microchannels, surface tension forces play an important role in determining the dynamics of bubbles whereas gravitational forces are generally negligible. Also it is very important to consider the interaction between the boundaries and the fluids by prescribing or computing the correct contact angle between them. The commercial CFD code FLUENT allows the use of static constant angles, or the use of User Defined Functions (UDF) to compute the dynamic contact angles. It is inappropriate to use a static contact angle to model cases involving moving contact lines. For such cases a dynamic contact angle scheme should be implemented. In this study, FLUENT was used to model adiabatic and diabatic, time dependent two-phase flows. Since FLUENT already contains a VOF method, a LS method was implemented and coupled with VOF into FLUENT via UDFs. Furthermore, since the LS function, used to compute the surface tension force, ceases to be a signed distance to the interface even after one time step, a re-initialization equation was solved after each time step. This involved using a fifth order WENO (Weighted Essentially Non Oscillatory) scheme to discretize the space derivatives (otherwise oscillations of the interface occurred), and a first order Euler method for the time integration. In another part of the study, a 3D dynamic contact angle model based on volume fraction, interface reconstruction, and experimentally available advancing and receding static contact angles was also developed and implemented into FLUENT via UDFs. Several validations for the developed CLSVOF method and dynamic contact angle model are presented in this thesis, these includes a static bubble, a bubble rising in a stagnant liquid for Morton numbers ranging from 102 to 10-11, droplet deformation due to a vortex flow field, droplets spreading over a wall under the gravity effect and droplets sliding over a wall due to gravity. These validations demonstrated the high accuracy and the stability of our methods for modeling these phenomena. A heat and mass transfer model was also implemented into the commercial CFD code FLUENT for simulating of boiling (and condensation) heat transfer. Several simulations were presented with water and R134a as working fluids. The influence of the contact angle and the wall superheat was also studied.

Label-free learning of elliptic partial differential equation solvers with generalizability across boundary value problems

A partial differential equation was developed and programed to compute the cumulation and flow of the liquid phase of a gas-condensate fluid. It was found that such cumulation near the wellbore may attain saturations permitting liquid flow within periods of the order of 1 day. The radius at which this saturation obtains advances toward the external reservoir boundary with continued recovery, although it may move only 2 to 5 percent of that distance during the productive life of the well. This relatively high saturation can cause a serious diminution in delivery capacity. The pressure differential required for flow of gas through a formation with such saturation may be as much as three times greater than the corresponding differential in a formation where saturation is less than 2 percent. Because of the radial geometry, total reservoir pressure differential required for flow with a zone of relatively high liquid saturation about the well was indicated to be two to three times greater than without the saturation. Introduction A computing program developed by the authors to integrate the partial differential equation for transient radial flow of gas-condensate fluids required a prevalued function, (1) related to cumulation of condensate in pores of the reservoir as pressure declines. Methods were available for relating mobility to saturation, saturation to liquid-vapor volume ratio, liquid-vapor volume ratio to pressure and, hence, mobility to pressure to obtain W(PD). However, certain assumptions had to be made in advance of program execution concerning movement of the liquid phase and change in radial distribution of liquid saturation as recovery progressed. Because reservoir fluid mobility and related pressure gradients are markedly dependent on liquid saturation, it was desirable to eliminate the necessity for these assumptions and in this respect improve the reliability of a result obtained by integration. Experimental research was begun to determine the relationship of variables that control the reservoir flow of low liquid-vapor ratio fluids. With the porous medium partly saturated, permeability to a two-phase flowing fluid could not be measured because the viscosity was unknown. Viscosity could not be measured without knowing permeability of the liquid-containing medium; but mobility, the ratio of permeability to viscosity, was measurable. This property, and liquid-vapor volume ratio of the flowing phases, have values determined by nature of the fluid and porous medium and, in addition, were found to be pressure, velocity and saturation dependent. A program developed by Kniazeff and Naville for computing multicomponent, two-phase flow used polynomials to represent saturation dependence of permeabilities for gas and for liquid phases, and an effluent viscosity constant under reservoir conditions. The authors developed a program based on a partial differential equation of similar mass balance origin but designed to represent distribution of an unlimited number of components between gas and liquid phases utilizing the experimental information available as to variable dependence. Attention was given to liquid-vapor volume ratio of the flowing phases, a property which is essential as a liquid saturation boundary condition. A program for computing implicitly saturation increase at uniform intervals along the drainage radius was combined with an existing program for calculating pressures, and the two were used alternately in an iterative solution for both saturations and pressures. The function W(PD) needed for calculating pressures was replaced by the ratio of mobility to compressibility factor. SPEJ P. 125ˆ

Partially wetting drops sliding down an inclined plane develop a “corner singularity” at the rear, consisting of two dynamic contact lines that intersect. We analyze the three-dimensional flow in the vicinity of this singularity by exploring similarity solutions of the lubrication equations. These predict a self-similar structure of the velocity field, in which the fluid velocity does not depend on the distance to the corner tip; this is verified experimentally by particle image velocimetry. The paper then addresses the small-scale structure of the corner, at which the singularity is regularized by a nonzero radius of curvature R of the contact line. Deriving the lubrication equation up to the lowest order in 1∕R, we show that contact line curvature postpones the destabilization of receding contact lines to liquid deposition, and that 1∕R increases dramatically close to the “pearling” instability. The general scenario is thus that sliding drops avoid a forced wetting transition by forming a corner of two inclined contact lines, which is regularized by a rounded section of rapidly decreasing size.

The oscillation of liquid/gas free surface in a partially filled storage tank caused by an abrupt drop of gravity level is of critical importance for the fluid management in space. In the present study, the dynamic behavior of free surfaces in a model tank (tube) is numerically investigated using volume of fluid (VOF) method in the context of dynamic contact angle (DCA) model. It is concluded that the dynamic behavior of free surface could be captured pretty well using the selected DCA model, as shown by comparison with the results of Drop Tower Beijing experiment. The temporal evolution of free surface reproduces exactly the characteristics of damping oscillations. The detailed dynamic deflections of meniscus reveal crucial dependency between the oscillation frequency of free surface and the boundary condition in the contact line. The oscillation frequency increases when the range of the moving contact line transfers from the spherical-shaped part to the cylindrical part of the tank and maintains constant when the moving contact line remains always at the cylindrical part of the tank. Meanwhile, the oscillation amplitude decreases in line with the increase of oscillation frequency.

CHAPTER 11 - PARTIAL DIFFERENTIAL EQUATIONS AND BOUNDARY-VALUE PROBLEMS