A mathematical formulation of the boundary integral equations for a compressible stokes flow

A mathematical formulation of the boundary integral equations for a compressible stokes flow

Francisco Ricardo Cunha^I; Aldo João de Sousa^I; Michael Loewenberg^II

^IUniversidade de Brasília, Departamento de Engenharia Mecânica-FT Campus Universitário, 70910-900 Brasília, DF, Brasil

^IIYale University, Department of Chemical Engineering New Haven, CT, 06520-8286, USA

E-mail: [e-mail]/ [e-mail]/ [e-mail]

ABSTRACT

A general boundary integral formulation for compressible Stokes flows is theoretically described within the framework of hydrodynamic potentials. The integral equation is implemented numerically to the study of drop expansion in compressible viscous flows. Marker point positions on the drop interface are involved by using the boundary integral method for calculation of fluid velocity. Surface discretization is adaptive to the instantaneous drops shapes. The interplay between viscous and surface tension and its influence on the evolving emulsion microstructure during its expansion is fundamental to the science and technology of foam processing. In this article the method is applied for 3D simulations of emulsion densification that involves an uniform expansion of a viscous fluid containing spherical drops on a body centered cubic lattice (BCC).

Mathematical subject classification: 76N99, 76N15, 45K05.

Key words: boundary integral, compressible drops, emulsion, densification.

1 Introduction

The dynamics of interface deformation in low Reynolds number flow is of interest in a wide variety of fields including chemical and petroleum engineering, solid-earth geophysics, hydrology and biology. Typical applications span an immense range of length scales from microns to hundreds of kilometers: biological studies of cell deformation; chemical engineering studies of flotation, coating flows and the dynamics of thin films. More recently, attention has been focused on problems of foam processing where a dense phase of viscous drops are dispersed throughout a second fluid, and in particular in the generation of a dense emulsion and the determination of foam rheology.

In the low Reynolds number limit the motion is governed by the Stokes and continuity equations. Although time-dependence does not appear explicity in Stokes equations, it is consistent to study time-dependent interface deformation. This quasi-static assumption requires that the time scale for the diffusion of vorticity

²/n, where n is the kinematic viscosity and characteristic length of the flow, is much less than a typical time scale for a drop deforms significantly, m/G. G is the surface tension. Physically, the absence of the temporal derivative in the equation of the motion does not necessarily imply that the flow is steady, but merely reflects the fact that the forces exerted on fluid parcels are in a state of dynamic equilibrium as a result of the rapid diffusion of momentum (or vorticity). Consequently the instantaneous structure of the flow depends upon the current boundary configuration and boundary conditions.

The boundary integral method relates velocities at points within the fluid to the velocity and stress on the bounding surfaces. It is an ideal method for studying free-boundary problems [1]- [2]. Advantages of the technique include the reduction of the problem dimensionality, the direct calculation of the interfacial velocity, the ability to track large surface deformations, and the potential for easily incorporating interfacial tension. The boundary integral formulation for incompressible Stokes flow was theoretically introduced by Ladyzhenskaya [3] within the framework of hydrodynamic potentials. This integral equation method was first implemented numerically by Rallison and Acrivos [4]. In recent years the number of applications including simulations of dilute and concentrated emulsion has increased enormously [5]-[9].

This work proceeds by considering the extension of the boundary integral formulation for a compressible Stokes flow. We provide some of the details of the basic integral equations, along with the original modifications necessary for the study of dense emulsions. We use an adaptive mesh restructuring algorithm in order to match the instantaneous surface configuration. For deformable drops, mesh restructuring is based on a curvature-rule, and is fully independent from the history of deformation i.e., velocity calculations [10]. A gap-rule is also implemented when near contact motion of interacting drops needs to be accurately resolved. In the last part of the article, we outline the numerical techniques typically applied in densification of emulsions.

The mathematical formulation, governing equations, boundary conditions and assumptions are discussed in §2. The reciprocal theorem and the integral representation for a compressible newtonian fluid are proposed in §3. In §4 we describe the numerical implementations and apply the method for generating configurations of dense emulsions and foams. Concluding remarks are made in §5.

2 Balance equations and boundary conditions

The theoretical formulation discussed in the present article will be applied to investigate densification of emulsions formed by drops of viscosity lm and radius a (undeformed shape at time t = 0) immersed in a second immiscible fluid of viscosity m with an externally imposed velocity field u^¥ (see figure (1)). In the following analysis, it is assumed that the Reynolds numbers for the flows inside and outside the drop are both extremely small, Re = ru^¥a/m, Re = ru^¥a/lm << 1. It will be helpful to keep in mind that l denotes the viscosity ratio between the internal and the external flow, and when l = 0 or ¥ the particle becomes a frictionless bubble or a rigid body, respectively.

2.1 Governing equations

In the regime of low Reynolds number, compressible fluid motions are governed by the Stokes and continuity equations

V denotes a closed region of fluid bounded by a surface S, u is the Eulerian velocity field, p is the pressure, b is an external body force per unit of mass, m and r are the fluid viscosity and density, respectively, and D = V^-1DV/Dt is the rate of expansion of the flow.

2.2 Boundary conditions

The boundary conditions on a drop interface S with surface tension G require a continuous velocity across the interface and a balance between the net surface traction and surface tension forces that express the discontinuity in the interfacial surface forces [2], [11]. Mathematically, these conditions are expressed as

where u¢ denotes the flow inside the drop. For an active interface free of surface viscosity, surface elasticity and surface module of bending and dilatation, the traction jump D = [[n · s]] constitutive equation is written as [12]

The notation [[ ]] denotes a jump in flow quantities, G is the interfacial tension acting between the drop and fluid phases, n is the unit normal vector to S, s is the Eulerian stress tensor, D = n · s is the surface traction, (I – nn) · Ñ denotes the gradient operator Ñ^s tangent to the interface, I is the identity tensor and Ñ^s · n denotes the mean curvature of the interface . In general, it can be expressed as the sum of the inverse principal radii of curvature = ().

The force balance at the interface (3) might include both a normal traction jump (n · D)n and a tangential stress jump (I – nn) · D due to gradients of surface tension G along the interface. The dynamic effect on interfacial tension gradients are known as Marangoni effects. These variations here are associated with the presence of surfactants in the fluid. In a real flow system, we must often expect G to vary from point to point on the interface, and it is important to consider how gradients of G may influence the deformation of the drop. To describe Marangoni stress effects additional convection-diffusion equation is necessary for determining the surfactant distribution along the interface [6].

A kinematic constraint relates changes in the interface position to the local velocity. Thus interface evolution of drops are described with a Lagrangian representation

3 Integral equations

The boundary integral formulation developed in this article provides a powerful method for computing compressible Stokes flow by solving integral equations for functions that are defined over the boundaries. The important benefits of this extended approach is the ability to track deformation and swelling of drops in very dense emulsions or foam.

3.1 Reciprocal theorem for a Newtonian fluid

In this section we derive the Lorentz reciprocal theorem for the general case of compressible Newtonian fluids. The calculations comprise an extension of the demonstration presented in [13]. The reciprocal identity find extensive applications in the study of Stokes Flows. The major strength of the reciprocal identity is that it allows us to obtain information about a flow without having to solve the equations of the motion explicitly, but simply by using information about another flow.

Consider a closed region of fluid V bounded by a surface S. Now consider two unrelated compressible flows of two different Newtonian fluids with densities r and r¢ and viscosities m and m¢, and stress fields s and s¢, respectively:

Flow 1: u, s; (r, m). The equations for conservation of mass and momentum in terms of the material derivative, D/Dt, for the flow 1 are respectively:

Here, locally, u is the velocity, s is the stress field and b is the external body force per unit of mass. The Newtonian constitutive equation for a compressible flow is given by Batchelor [14]

where I is the identity, E = (Ñu + Ñ^Tu) is the rate of strain tensor and Ñ^Tu denotes the transpose of the tensor Ñu.

Flow 2: u¢, s¢; (r¢, m¢). Similarly, the conservation and constitutive equations for the flow 2 are respectively

where E¢ = (Ñu¢ + Ñ^Tu¢) is the rate of deformation for the fluid 2. First consider the tensorial operation

but for a Newtonian compressible flow

D denotes the rate of expansion of the flow [14]. Thus (9) may be written as

If the same steps are applied to s¢ : E; it must reduce in an analogous fashion to

One require E : E¢ = E¢ : E so that

Now, substituting (13) into (12), ones obtain

As a consequence of the symmetry of the stress tensor s : Ñu¢ = s^T : Ñu¢ = s : Ñ^Tu¢. Then, one may write that

In this step we make the dot product of Cauchy's equation (5) by u¢ in order to define the last term in the RHS of (15). Therefore after substituting back the result of this operation into (15), it gives

By reversing the role of the primed and unprimed variables, it is also possible to obtain

Now, substituting (16) and (17) into (14), multiplying the resultant equation by m and make few manipulations, we obtain an expression for the generalized Lorentz reciprocal theorem for a Newtonian compressible flow (reciprocal identity):

For low Reynolds number flow (i.e. Stokes flows) the reciprocal identity (18) takes the simpler form

3.2 Integral representation for a compressible Stokes flow

Consider the particular flow of interest with velocity u and stress tensor s. The known flow is the one due to a point force with strength h, and located at a point x_o. Suppose that the inertia of both fluids has a negligible influence on the motion of the fluid elements, and by convenience takes m = m¢ and r = r¢. Flow 1 and flow 2 for this particular situation are described as following.

Flow 1: u, s. The equations for conservation of mass and momentum in terms of the material derivative, D/Dt, for the flow 1 are respectively:

where rb = B is the body force by unit of volume and D = Ñ · u is the flow rate of expansion.

Flow 2: Fundamental solution of the Stokes flow; u¢, s¢. The fundamental solution for Stokes equations correspond to the velocity and stress fields at a point x produced by a point force h located at x_o:

with | u¢ | 0 and | s¢ | ¥ as | x | ¥. The solution of such equations may be derived, for example, using Fourier transforms [2]

where, the stokeslet, G and the stresslet T are defined by the following expressions:

The above functions are the kernels or the free-space Green's functions that maps the force h at x_o to the fields at x in an unbounded three-dimensional domain. Here = x – x_o, and r = | |. Physically u = G() · h expresses the velocity field due to a concentrated point force hd(x – x_o) placed at the point x_o, and may be seen as the flow produced by the slow settling motion of a small particle. T_ijk is the stress tensor associate with the Green's function G_ij and s_ik(x) = T_ijkh_j is a fundamental solution of the Stokes produced by the hydrodynamic dipole D · Ñd(x – x_o). T_ijk = T_kji as required by symmetry of the stress tensor s. It is straightforward to show that the Reciprocal theorem for the present case, equation (19), takes the form

Now, considering the body force exerted on the flow (u, s) the gravity force B = Ñ(rg · x), substituting the expressions of the point-force solution into (25) and discarding the arbitrary constant h ones obtain

Note that we have used for the first term on the RHS of equation (26), the incompressibility of the singular solution Ñ · G = 0, so that G · Ñ(rg · x) = Ñ · [G(rg · x)]. The above equation is valid everywhere except at the singular point x_o.

Consider a material volume of fluid V bounded by the singly or multiply connected surface S (see figure 2 a,^b) in order to evaluate the integration of equation (26). There are two situations to be considered next.

Point force outside V. For this case d(x – x_o) = 0 inside V and thus after integrating equation (26) the integral representation of the Reciprocal theorem takes the form

The volume integrals in equation (27) can be converted to the surface integrals over S, by using the divergence theorem obtaining

where n is the unit outward normal to the surface S. Equation (28) is the integral representation of the flow if the singularity is outside V. It will be shown that the integral equation (28) is a useful identity for developing new integral equations in terms of jump conditions on the interface.

Point force inside V. In order to formally determine the integral representation for the situation which d(x – x_o) ¹ 0 inside V, we must integrate again the equation (26). Applying the divergence theorem and d-distribution property, one obtains

Equation (29) is the integral representation for compressible Stokes in terms of four boundary distributions involving the Greens's functions G, the stresslet T and the potential source 1/r. The first distribution on the RHS of (29) is termed the single-layer potential, the second distribution is termed double layer potential, whereas the last new term is a potential source distribution due to the compressibility of the flow with a constant rate of expansion.

3.3 Integral representation in terms of the traction jump

External flow representation. Using the reciprocal identity (28) for the internal flow u¢ (inside the particle) at a point x_o that is located exterior to the particle, one obtain

Now, applying equation (29) for the external flow subject to an ambient flow, u^¥, and combining the result with equation (30), the integral representation is obtained as a function of the traction jump D,

Internal flow representation. We repeat the above procedure for the internal flow. Hence, the integral representation of the internal flow is obtained when equation (29) is applied,

Using the reciprocal identity (26) for the external flow u¢ (outside the particle) at a point x_o that is located in the interior of the particle, one obtain after dividing the full equation by l, that

The integral representation of the internal flow as a function of the jump condition is obtained by combining (32) and (33)

3.3.1 Integral representation for the interface.

The integral representation for the flow solution at the interface is found by applying the jump condition (1/2) [u(x_o) + lu¢(x_o)] to the equations (31) and (34) (see figure 1). For the limit of x_o going to the interface, u(x_o) = u¢(x_o) (continuity of velocity) and the traction discontinuity D is given by the equation (3). Under these conditions only the integral representation for the fluid-fluid interface S need to be considered, hence

where A = 1/4pm and B = (l – 1). It should be noted that when the viscosity ratio l = 1.0 and Dr = 0 the double layer and the single layer integrals (related to the buoyancy force) vanish and the flow is expressed merely in terms of a single-layer potential with known density force D and the source potential integral. This means that the same fluid is occupying all space, but with a membrane of points force and sources provide by the singularities at the positions of the interface.

4 Application

In order to generate realistic high-volume-fraction microstructures, one propose to simulate the centrifugation process by which high-density emulsions are produced from low-density materials [15]. The extraction of the continuous-phase fluid during this densification process is equivalent to a distributed sink of continuous-phase fluid. Thus, we will simulate this densification process by describing the evolution in a system with a uniformly-distributed sink of continuous-phase fluid.

4.1 Interacting drops

All simulations rely on the boundary integral method. Periodic boundary conditions are enforced through the use of periodic Greens functions. These are obtained by Ewald summation [16] using accurate computationally-efficient tabulation of the nonsingular background contribution [7]. The appropriate formulation was derived in § 3. The resulting boundary integral formulation are now capable of describing the dynamics of dense emulsions for uniform compressibility. Accordingly, the evolution of M neutrally buoyant deformable drops is described by time-integrating the fluid velocity u(x_o) on a set of interfacial marker points x₀ on each drop surface. In the present application it is considered the case in which there is a rapid equilibrium of insoluble surfactants (incompressible surfactants). Therefore, Marangoni stresses and adsorption-desorption of surfactants can be ignored.

All quantities below are made dimensionless using the (volume-averaged) drop size a and the relaxation rate G/ma. The relevant physical parameters that describe the simulated system simulated are: l, f and the compressibility parameter, Ca_o = amD/G (ratio of viscous to surface tension stress). Ca_o is the appropriate capillary number for the densification process. In the absence of an imposed flow (i.e. u^¥(x_o) = 0), the dimensionless fluid velocity is governed by the second-kind integral equation on the interfaces S_m (m = 1, ¼, M) of all simulated drops.

where

G^P and T^P are respectively the periodic stokeslet and stresslet defined as in reference [17], F^P is the periodically-replicated r^–1 potential and f(f) = –[1 + l(f^–1 – 1)].

The volume-averaged stress tensor from the dispersed phase S is obtained from an integral of the traction jump and fluid velocity over the drop interfaces [14],

Equation (38) is the contribution of the dispersed phase to the macroscopic stress of the emulsion due to the dipole stresslet that each drop torque free generates in the flow. The (volume-averaged) non-equilibrium osmotic pressure during densification is tr(S) – 2f, where tr(S) is the trace of the stress tensor, and 2f is the contribution from the capillary pressure of spherical drops. When densification stops, the drop shapes relax, and the stress relaxes to the equilibrium osmotic pressure which depends only on the drop shapes.

4.2 Numerical results

Evolution of a drop surface S was simulated by means of a surface discretization with initial number of marker points N_o [10], [17]. Marker points are convected with the fluid velocity. During the simulation mesh restructuring is performed on S. After each time step, first marker points are added/subtracted on S as required by condition (39); then global mesh equilibration and reconnection are performed. The surface discretization is equilibrated as a dynamical system of springs. An equilibrium configuration is found by direct numerical simulation of evolution of the system of springs using a second order Runge-Kutta scheme to preserve accurate description of the interface for evaluation of the normal vector n at each integration step. The iterative process is stable and converges quickly to an equilibrium configuration for a smooth distribution of traction. The normal vector and curvature were calculated by the local surface-fitting algorithm of [18]. The fluid velocity on the drop interfaces are obtained by an iterative solution of (36) using the GMRES algorithm (a generalization of the conjugate gradient method to non-symmetric matrices) to achieve convergence for the closely-spaced interface configurations that characterize dense emulsions. Once fluid velocity is known, positions of marker points are evolved by a second Runge-Kutta scheme. An appropriate time step that is proportional to the shortest edge length is set in order to ensure stability.

The adaptive surface triangulation algorithm, described in [10], has been extended to construct efficient simulations of dense systems. Accordingly, a new marker point density function was defined that resolves the minimum local length scale everywhere on the drop interfaces. For the proposed problem the minimum local length scale may depend on the local curvature or local film thickness h. Thus, the marker point density function [10] should be generalized to:

where R₁, R₂ are the local principal radii of curvature, and C₁ are O(1) constants whose precise value is unimportant. The proposed marker point density function (39) resolves the rim of dimpled regions where the film thickness varies rapidly and the lubrication length scale in regions where the film thickness is slowly-varying (R = min[R₁, R₂]). Only rim regions, not flat regions, require high resolution.

An inspection of the result in figure 3, which was obtained using the density function (39), illustrates this point: high resolution is needed on the plateau borders and junctions, not on the flat films between drops. This observation is important for ensuring the feasibility of accurately resolving the dynamics of the closely-spaced interfaces that characterize the dense systems that we have explored.

The surface integrations in Equations (36), involve singular integrand. Trapezoid-rule integration with singularity subtraction and near-singularity subtraction for closely-spaced interfaces of drops (if surface tension gradients are not present) [7] can be used to accurately evaluate the integrals. Equation (36) has eigensolutions that cause unphysical changes in the dispersed-phase volume at small viscosity ratios, corrupt numerical solutions at large viscosity ratios, and slow the iterative convergence. These effects were eliminated by implementing Wielandt eigenvalue deflation described in [2] in order to purge the solutions corresponding to l = 0 and l = ¥. The resulting surface integration algorithm is economical and O(1/N) accurate, consistent with the triangulated representation of the drop interface.

A wide range of densification processes can be simulated through the appropriate f dependence (i.e., time-dependence) of D. It is also possible to explore microstructural manipulation through an applied shear during densification. Here, we consider constant compressibility, quiescent flow conditions and the case l = 1.0 because the numerical implementation is simpler. Densification with prescribed non-equilibrium osmotic pressure has been also explored since it corresponds most closely to the experimental procedure that has been developed [15]. Under these conditions, densification proceeds until the applied osmotic pressure is balanced by the equilibrium osmotic pressure of the emulsion. The Kelvin-cell and Weaire-Phelan microstructure, depicted in figures (6) and ⁽⁵⁾, were obtained using an algorithm based on the compressible formulation discussed in § 3.

Spherical drops on a BBC lattice make first contact with eight nearest neighbors and form precursors of hexagonal faces when f < f_c = p3^1/2/8 = 0.68 (maximum packing for BCC). The results for Ca_o = 0.5 contained in figure (6) show the evolution of drop shape from spheres to Kelvin cells. Once the emulsion expansion has stopped, the drop shape will continue to relax toward equilibrium. This process is illustrated in figure (6); the emulsion expand and then relaxes for twenty time units. We have no experimental results with which compare our numerical predictions.

5 Conclusion

The formulation of boundary integral equations for compressible Stokes flow has been discussed. The approach was applied for generating dense compressed emulsion structure in viscous flows with periodic boundary conditions. We have made substantial progress on the problem of emulsion densification. The results demonstrate the feasibility of simulating high-volume-fraction systems. A study of densification may have interesting materials processing applications that will be pursued in a future study.

6 Acknowledgment

The authors are grateful to the CAPES-Brazil, CNPq-Brazil and CT-Petro-Finep for their generous support of this work. We would like to thanks Dr. Vittorio Cristini and Dr. Jerzy Blawzdziewicz for helpful discussions on their algorithm for adaptive mesh restructuring of three-dimensional drop surface.

Received: 25/IX/01.

#539/01.

• [1] G.K. Youngren and A. Acrivos, Stokes flow past a particle of arbitrary shape: A numerical method of solution, J. Fluid Mech., 69 (1975), pp. 377-403.
• [2] C. Pozrikidis, Boundary integral and singularity methods for linearized viscous flow, Cambridge University Press, Cambridge, (1992).
• [3] O.A. Ladyzheskaya, The mathematical theory of viscous incompressible flow, Gordon & Breach, New York, (1969).
• [4] J.M. Rallison and A. Acrivos, A numerical study of the deformation and burst of a viscous drop in an extensional flow, J. Fluid Mech., 89 (1978), pp. 191-200.
• [5] J.M. Rallison, The deformation of small viscous drops and bubbles in shear flows, Ann. Rev. Fluid Mech., 16 (1984), pp. 45-66.
• [6] H.A. Stone, Dynamics of drop deformation and breakup in viscous fluids, Ann. Rev. Fluid Mech., 26 (1994), pp. 65-102.
• [7] M. Loewenberg and E.J. Hinch, Numerical simulations of a concentrated emulsion in shear flow, J. Fluid Mech., 321 (1996), pp. 395-419.
• [8] F.R. Cunha, M. Neimer, J. Blawzdziewicz and M. Loewenberg, Rheology of an Emulsion in Shear Flow, In: Proceedings of American Institute of Chemical Engineers, Annual Meeting, Fundamental Research in Fluid Mechanics: Particulate & Multiphase Flow II., paper 124j., October (1999).
• [9] F.R. Cunha and M. Loewenberg, Emulsion Drops in Oscillatory Shear Flow In: 52nd American Physical Society, Annual Meeting of the Division of Fluid Dynamics, paper GK. 08, session: suspension I, November (1999).
• [10] V. Cristini, J. Blawzdziewicz and M. Loewenberg, Drop breakup in three-dimensional viscous flows, Phys. Fluids, 10 (1998), pp. 1781-1783.
• [11] G.L. Leal, Laminar flow and convective transport processes, Butterworth-Heinemann, Boston, (1992).
• [12] L.E. Scriven, Dynamics of fluid interface, Chem. Eng. Science, 12 (1960), pp. 98-108.
• [13] J. Happel and H. Brenner, Low Reynolds number hydrodynamics, Prentice Hall, (1965).
• [14] G.K. Batchelor, An introduction to Fluid Dynamics, Cambridge University Press, Cambridge, (1967).
• [15]T.G. Mason, M.D. Lacasse, S.G. Grest, D. Levine, J. Bibette and D.A. Weitz, Osmotic pressure and viscoelasticity shear moduli of concentrated emulsions, Physical Review E, 56 (1997), pp. 3150.
• [16] C.W.J. Beenakker, Ewald sum of Rotne-Prager tensor, J. Chem. Phys., 85 (1986), pp. 1581-1582.
• [17] M.P. Almeida and F.R. Cunha, Numerical simulation of non-colloidal deformable particles at low Reynolds numbers, In: Proceedings, 7th Brazilian Congress of Engineering and Thermal Sciences, 1 (1998), pp. 963-970.
• [18] A.Z. Zinchenko, M.A. Rother and R.H. Davis, A novel boundary-integral algorithm for viscous interaction of deformable drops, Phys. Fluids, 9 (1997), pp. 1493-1511.

Publication Dates

• Publication in this collection
  19 July 2004
• Date of issue
  2003

History

• Received
  25 Sept 2001