- Email: [email protected]
Low Reynolds number deformation of viscous drops in a bounded flow region under surface tension
MATHEMATICAL COMPUTER MODRLLJNG PERGAMON
Mathematical
and Computer
Modelling 31 (2000)
99-118 www.elsevier.nl/locate/mcm
Low Reynolds Number Deformation of Viscous Drops in a Bounded Flow Region under Surface Tension A. R. M.
PRIMO*
AND
L.
C.
WROBEL
Brunel University, Department of Mechanical Engineering Uxbridge, Middlesex UB8 3PH, United Kingdom H. POWER Wessex Institute of Technology, Ashurst Lodge Ashurst, Southampton SO40 7AA, United Kingdom Abstract-This paper presents a novel formulation of the completed indirect boundary element method to study the deformation of viscous drops on a slow viscous flow in a bounded region subject to surface tension effects. The theoretical background to the technique is explained in detail, including mathematical proofs of existence and uniqueness of solutions of the second-kind Fkedholm integral equations generated. The computer implementation of the formulation is simple and efficient. Numerical results of several tests are included. @ 2000 Elsevier Science Ltd. All rights reserved.
Keywords-Boundary
element method, Viscous drops, Stokes Aow, Surface tension.
1. INTRODUCTION The practical importance of studying the motion of drops and bubbles is due to their common occurrence in many industrial and biological systems, as well as in a number of technological processes. Mechanical engineers have studied droplet behaviour in connection with combustion operations, and bubbles in electromachining and boiling. Chemical and metallurgical engineers rely on drops and bubbles for such operations as distillation, absorption, flotation, and spray drying. The investigation of the mechanics of drops and bubbles has a long record in fluid dynamics, and still involves a substantial portion of pure and applied research. A general review of the subject can be found in [l]. Rallison and Acrivos [2] developed a numerical solution for the low Reynolds number deformation of a viscous drop, using Green’s integral representation formulae for the fluids inside and outside the drop. The external fluid region was infinite in extent, and the viscosity ratio between the two fluids was defined by X = p~/pr, with ~1 and ~2 the dynamic viscosities of the external fluid and the drop, respectively. A second-kind F’redholm integral equation was derived for the unknown surface velocity at the interface, for given external fluid velocity and interfacial surface tension. The interface deformation was found from a kinematic condition relating the rate of displacement of the interface to the normal velocity component. Although Rallison and Acrivos [2] only conjectured that their integral equation did not have any eigensolution for 0 < X < 00, since *On leave from Department
of Mechanical
Engineering,
Universidade
Federal de Pernambuco,
08957177/00/$ - see front matter @ 2009 Elsevier Science Ltd. Ail rights reserved. PII: SO895-7177(99)00218-S
Typeset
Recife, Brazil. by A&-‘B$
100
A. R. M. PRIMO et al.
their numerical solution encountered no difficulties for all values of X tested in this range, this fact was later proved analytically by (31. The main objective of the present work is to extend Rallison and Acrivos’ [2] formulation for the case when the external fluid region is finite and its exterior boundary is subject to known surface tension. For this problem, the cases when X ‘= 0 and X = 00 have a direct application for viscous sintering with air bubbles and solid inclusions, respectively; the cases when 0 < X < co could be seen as viscous sintering problems with fluid impurities. Integral equation formulations for viscous sintering using the direct boundary element method have been developed by [4-61. A more recent formulation for the shrinking of air bubbles in a sintering flow (i.e., X = 0) was proposed by [7]. The present approach is based on an indirect integral equation formulation for the secondkind boundary-value
problem in terms of only one surface potential.
Recasting the Stokes flow
equations in integral form produces a system of Fredholm integral equations of the second kind for the unknown surface density which is shown to possess eigenvalues corresponding to rigid body motions, restricted to the exterior surface. The eigenvalues are removed by the addition, to the corresponding integral equation, of a term which is linearly proportional to the rigid body motions. Theoretical proofs are included that the addition of this new term does not alter the solution of the problem, and that although the original integral equation and the completed one are identical, the former does not have a unique solution while the latter does. The computational implementation of the completed indirect BEM formulation is simple and efficient. Results of several tests are included to verify the theoretical and numerical procedures.
2. MATHEMATICAL
FORMULATION
A compound drop is considered consisting of a fluid drop with viscosity ~2 completely coated by another drop of viscosity ~1 (see Figure 1). The present analysis can be directly extended to a finite number of internal drops of different viscosities without any difficulty. Let 7~ be the interfacial surface tension between fluids 2 and 1, and yr between fluid 1 and the exterior environment. The governing equations for the fluid velocity ??and pressure p are the Stokes equations auln =O a 9 dXi
(1) where at&'" J+aXj
(
U$=-pmbij+j.Lm
dU7 dXi
)
,
(2)
with m = 1 for a point x E Rr and m = 2 when x E 02. The flow fields have to satisfy the following boundary conditions at every point 5 E I’r: ui'j (Ebj(E)
=
%K(t)ni(5)9
(3)
and the following matching conditions at every point c E I’s: u;(5) = G(C)
(4)
and (5) In the above set of equations, Ir is the exterior surface and l?s the interface between fluids 1 and 2, both assumed to be Lyapunov surfaces.
Low Reynolds Number Deformation
101
Figure 1. A viscous drop embedded in a finite viscous region.
The motion of the drops surface cannot be arbitrary since at each instant in time the force and torque yielded by the surface tractions have to satisfy the following compatibility conditions:
and aFjnjqf dS = 6,
s l-2 where 7,
1 = 1,2,3,
are the three linearly-independent ‘Pf =
417
cp (x) =
(22,
-3
(7)
rigid body vectors, taken as
for 1 = 1,2, -21).
Equation (7) follows from the fact that the flow field inside lY2 is a regular interior Stokes flow field yielding zero force and torque upon its surface container, while equation (6) is due to the balance of force and torque between the surfaces l?r and I’z, taking into account that the flow between these two surfaces has to be a regular Stokes flow field. Adding equations (6) and (7), using the boundary condition (3) and the matching condition (5), leads finally to the solvability condition
(8) where the normal vector ?? has been defined outwardly to the domain bounded by the corresponding surface boundary, i.e., I’1 and l?s. In order to guarantee continuity of the velocity field across the surface l?s, equation (4), the velocity field inside the domain fir U Cl2 will be defined by the following single-layer potential:
t&(x) = -
s rlur2
for every 3: E 521 or &. is given by
u{(z,Y)+j(Y)dSv,
The fundamental solution of the Stokes system of equations (Stokeslet)
A. R. M. PRIMO et al.
102
where T = )z - Yf. For convenience, the dynamic viscosity of each fluid will be defined in the corresponding pressure field. Therefore, at each fluid domain
(11) with
(12) and 2 = 1,2. Using the above representation for the flow field in each fluid region and letting a point z E Or approach a point < E I’r, taking the boundary condition (3) into consideration, the following equation is obtained:
where the kernel is given by Kji(Y4
= -;
1 (Yi - 6) (Yj - Ed(Yk - 5d
nh(E)
r4
Similarly, letting a point z E Ql approach a point < E l?z, an expression for the surface traction at such surface is obtained in the form
Letting a point x E R2 approach a point t E l?2
Subtracting equation (15) from equation (14), taking into account the matching condition (5), the following expression is obtained for a point < E l?2: -r2f4lh(t)
-
(111 -
cl2) [I
ra
Kji (Y, Wj
= -
(I4 ;
(Y) dSu +
p2) I&(<)
J rl
Kji(y, O@j(y)
dsv
1 .
Equations (13),(16) g ive a system of second-kind Fredholm integral equations for the unknown surface density 3 over I’i U r2. For convenience, these equations will be normalized with respect to the viscosity ~1. In this way, the equations can be rewritten as n49W)
=
;Q~(E)
-
E)Qj (Y) d& y Jrl Kj~(Y,Wj(Y)dSg +Jr2 Kji(Yv I
(17)
when < E I?r, and -T2fi(S)%(S)
=
--
(l; 3 o&g (18)
Kji (Y,t)cpjb) dsy +
J rl
when < E I’z, where 57 = TL/~I, 1 = 1,2, and A = p2/p1,
Kji(Y, Wj (Y>dsg
1
7
LOW hynolds
Number Deformation
103
According to F’redholm’s alternative, to analyse the uniqueness of solution of the above system of integral equations, the following homogeneous system has to be considered: &W
-
WY,
E)@;(Y) &,
+
J r2
K&I,
S)@;(Y) dS,
I
= 0,
(19)
for every < E rr, and
for every 6 E l?s. The above homogeneous system of integral equations implies that the flow field defined by
t&x) =
s r1
&,
Y>@~"(Y)ds,
(21)
UPa
and P’(Z) = w
I rlur2
q%c, y)@;(y) dS,,
(22)
where 1 = 1,2 when z E RI or z f Rs, respectively, and 3 is a nontrivial solution of the homogeneous system (19), (20), is a regular interior Stokes flow in Rr U &, having zero surface traction at the surface Fr and continuous surface traction across the surface l?z. Therefore, according to Green’s first identity for the Stokes equation in each flow region,
and (24) From the above two equations, it follows that
where the continuity of the surface traction across the surface I?2 has been used. Taking into account the zero surface traction at the surface I’r
Therefore, the flow field defined by equations (21),(22) has to behave as a rigid body motion -i_- - cpn, n = 1,2,3 (where continuity of the velocity field across the 1_ - 21 in s1r U 02: i.e., u surface r2 has been used), with an associated zero pressure all over sir U 522. This is due to the zero traction condition at I’1 and continuity of the surface traction across rs, since such flow field is the only solution of the Stokes flow
compatible with the given homogeneous conditions.
A. R. M. PRIMO et al.
104
From the continuity property of the velocity field defined by equation (21) across its densitycarrying surfaces, it follows that the limiting value of such velocity field when a point of the domain exterior to I’r, z E R,, tends to a point 5 E I’r, is also a rigid body motion, i.e., u&, The surface traction
Here, oij (uf,#)(,)
n = 1,2,3,
= CPXSL
for every < E l?r.
at I’r due to such external flow field is given by
is the limiting value of the stress oij(Z,p’)
rl,
(25)
when a point z tends to a
point < E coming from R,, with the auxiliary exterior domain filled with the same fluid of domain Sir. From equations (19),(26), it follows that Uij
(
f
u
,p
1
>
celnj
=
for every 5 E
-@p(<),
rl.
(27)
J?rom the uniqueness of solution of the exterior first boundary-value problem for the twodimensional Stokes system of equations satisfying the corresponding conditions at infinity and according to Chang and Finn’s [8] interpretation of the Stokes’ paradox, it follows that the vector density $
at the surface l?r has to be given by any of the eigensolutions of ;W)
-s,,
for every 5 E
Kji(Y, E)Q$ (Y) dSu = 9,
rl,
(28)
or any linear combination of them. The above homogeneous equation has precisely three independent linear solutions, which are generally unknown, since its corresponding adjoint homogeneous equation only admits the rigid body motions as nontrivial solutions. To prove the above, we consider the following single-layer potential:
Q(x) = J l-1
&xv Y)@;(Y) ds,
defining a regular Stokes flow in the interior or exterior of rr filled with the same fluid, with $ being any nontrivial solution of equation (28). The above singlelayer linearly independent flow fields in the domain interior to J?r, i.e., T(x)
1 = 1,2,3
= J(z),
and
potential
defines three
p(x) = 0,
(39)
according to the zero surface traction condition (28), yielding zero internal stresses, i.e.,
Uij
(T(X))= -i S, (Yi
-
xi) (Yj - Xj) (Y/C - Xk)
Q%(Y)ds,
r4
1
= 0,
(31)
for every x inside l?r. By the continuity pro erty of the velocity field due to a single-layer potential, it follows that the limiting value of 3 (x) when a point z E Q2, tends to a point < E I’r is also a rigid body motion, i.e., W)(e)
On the other hand, the surface traction oij(q)nj point z crosses l?r, given by aij
nj
($1
-
for every [ E rr.
1 = 1,2,3,
= 4(E),
U+,j
(32)
corresponding to (29) experiences a jump as a
(e)
nj
=
K&i,
1 = 1,2,3,
Low Reynolds Number Deformation
105
where the limiting value of the surface traction when a point belonging to the interior of l?l tends to a point on the surface I’l, i.e., rrij( V)(~nj, Therefore,
from the above equation,
is given by the left-hand side of equation
the surface traction
(28).
exerted upon I’1 by the exterior flow
field defined by the velocity (29) with its corresponding pressure is given by flij
(e)
nj =
-$,
According to Chang and Finn’s [8] interpretation
1 = 1,2,3. of the Stokes’ paradox, the above exterior
flow field representing the flow due to the motion of the surface rl as a rigid body, and having a logarithmic
behaviour at infinity, always yields a nonzero total force, i.e.,
I’$(Y)(P~(Y)Q, #
for 1 = 1,2,3
0,
and
n = 1,2.
(34)
* l-1
Comparing equations (25) and (27) with equations (32) and (33), taking into account the logarithmic behaviour of the flow at infinity and the uniqueness of solution for such unbounded flow field, it fol%ws that the vector density $ at every point < E I’1 has to be equal to the eigenfunctions Q’, 1 = 1,2,3, or any linear combination of them. It is easy to see that such eigenfunctions are the only nontrivial solutions of the homogeneous system (19),(20), since it follows from equation (31) that the second integral in equation (20) has to be identically zero, given that the curve I’2 is completely embedded in I’l, i.e.,
.I
K,,(Y,
W;(Y) dSy = --
nj (0
l-1
7r
s
(Yi -
l-1
Ei)(YJ - 0) (Yk - Ek)Q;(y) r4
=
dS
Y
o
1
for every < E r2. Therefore, equation (20) has been reduced to the homogeneous form of the adjoint of the secondkind integral equation proposed by Rallison and Acrivos [2] to solve the problem of viscous drop deformation, which is known to admit only a trivial solution when 0 < X < 00, i.e., G(e) E 0 for every < E l?z (see [3]). The cases X = 0 and X = co, for which the integral equation does not possess a unique solution, correspond to the problems of a gas bubble and a solid particle, respectively. In this way, equation (19) reduces to (28) showing that 8-i are the only eigensolutions of the homogeneous system (19),(20). To remove the above eigenvalues, which are yly linearly proportional
to the rigid body vector
restricted to the contour r‘l, a term which is
cp’, 1 = 1,2,3, is added to equation (17), i.e.,
where
It is easy to prove that the addition of this new term to equation (17) does not alter the nature of the problem since for any admittable nonhomogeneous term, i.e., satisfying the compatibility condition (6), the integral fll in equation (35) will end up being zero. A similar approach has been proposed by Botte and Power [9] in order to complete the deficient range of the doublelayer integral operator when the solution of the interior first-kind boundary-value problem for the Stokes system of equation is represented in terms of a double-layer potential alone. In such a case, the original second-kind Fredholm integral equation, coming from the jump property of the double-layer potential, is completed by adding a term linearly proportional to the surface normal vector. The coefficient of proportionality is identically equal to zero when the nonhomogeneous
A. R. M. PRIMO et al.
106
term, i.e., given surface velocity, satisfies the corresponding compatibility condition, i.e., the noflux condition.
Although for any admittable given surface velocity the original integral equation
and the completed one are identical, the former does not have a unique solution since its corresponding homogeneous equation has a nontrivial solution, whilst the completed one possesses a unique solution since its corresponding homogeneous equation admits only a trivial solution. Multiplying the expanded form of equation (17) by the rigid body vector, qP,
at a point on
the surface I’1 and integrating the resulting expression with respect to 5, i.e., evaluating the total force and torque upon the surface J?r, gives the following expression:
where the order of integration in the double integrals has been interchanged and the first two terms of the original expression are cancelled since
Kji(Y, J)$(S)
1 dSc dSg = 2
J r1
@j(~)Ipjp(~) d%
according to the adjoint of equation (28). Relation (36) can be reduced to -
Qj(Y)(pjp(Y>dSy +Pl J Jrl Jl-2 (P~P(Y)v;(Y)~% =
rl TI~Y)~~(Y)Y$‘(Y)~%
(37)
by using the property that a double-layer potential with a rigid body vector as density yields an interior flow field that behaves like a rigid body motion, i.e.,
J Kji(Y, t)(PP(t)dsc= rl
for every y inside rr.
'P;(Y),
For the following analysis, it is important to take notice that such double-layer identically equal to zero outside its density-carrying surface. Now, multiplying equation (18) by the rigid body vector qp
potential is
at a point of the surface I?2 and
integrating the resulting expression with respect to < gives -
J @i(t)dK)d% = r2
J r2
W4~)ni(O,d’(~) dsc,
(38)
where
for every y outside I?2 as noticed above. Therefore, from (37), (38), and the solvability condition (8), it follows that P41
‘P%$P~(<) d% = 9.
(39)
The above linear algebraic system for pr, &, ,Ba only admits the trivial solution, because the determinant of equation (39) has the term
J cpldds, r1
n,l = 1,2,3,
Low Reynolds Number Deformation
107
as element in its lth row and nth column, and is thus the Gram determinant for the vector functions 2,
2,
and 2
with a nonvanishing value, on account of the linear independence of 2,
k = 1,2,3. Having proved that the additional term in equation (17) does not alter the nature of the original problem, it can now be shown that the corresponding expanded homogeneous system, i.e., when the term
&)
s,, q(Y)Pg(Y) dS?J
is added to equation (19), only admits the trivial solution. From the above analysis it also follows that
J r1
Q;(Y)cP:(Y)
dsy = 0,
since the solvability condition is automatically
for n = 1,2,3,
satisfied when Ti = Tz = 0.
Thus, the new system always reduces to the original homogeneous system, with nontrivial solution $ = -7 ‘p , 1 = 1,2,3, when E E l?i and zero otherwise. In this way, the orthogonal relations (40) remove the eigensolutions according to the nonzero total force condition (34), in order to comply with the Stokes’ paradox. Therefore, it has been established that the expanded form of the nonhomogeneous system (17),(18) h as a unique continuous solution for any physical admittable data when 0 < X < 03. The above procedure can be directly expanded to the case of n different internal drops by expressing the density-carrying surface of the single-layer potential (9) in terms of the union of all the surfaces involved and using the corresponding pressure field. In this case, the solvability of the corresponding system of integral equations will be restricted to the values of the parameters O
The system of
(41)
and
TzCi -R;(l -- (1 - A) (1 + A)
+x)
[I rZ
Kji(Y, C)@j(Y) d% +
= -&(C)
J
Kji(y,
tP)j(y)
d&,
r1
1 ,
(42) t. E
I-2.
The unique solution of the above system is
(43) and (44) where the fact that every rigid body motion is orthogonal to the normal vector was used. This is a consequence of the rigid body motion eigensolution of the adjoint equation of equation (28), i.e., &(Y&#Y)
ds, = 0,
c E rl
and
1 = 1,2,3,
A. R. M. PRIMOet al.
108
and of the fact that the only nontrivial solution of the above homogeneous equation is the normal vector, i.e., L
Jr1 integrating over l?i, and using the adjoint equation
Multiplying the above equation by J, of (28), it follows that
n for 1 = 1,2,3.
Since the density vector is proportional to the normal vector at each of the surfaces, it follows that the flow velocity defined by the single-layer potential is equal to zero, as expected from this stable configuration. It is important to observe that the above solution also holds if the circles are nonconcentric. Similar conditions are encountered when n circular drops are embedded in a large circular drop. In all these cases, the resulting flow velocity is the state of rest.
3. NUMERICAL 3.1. Indirect
Boundary
Element
STMULATIONS
Method
The simulations carried out by using the Indirect Boundary Element Method (IBEM) were based on the expanded form of equations (17) and (18). I n a more generalized form, these equations can be written as
and -TsK(E)ni(E)
= -(I
+ %j(r)@j(E)
- (1 - A)
sr
Kji(g,Wj(Y)d%,
5cr2,
(46)
with r = rl ur2. The term C’ij corresponds to the characteristic term (see [lo, p. 166]), given as &j, Cij(6)
=
c E Q,
Cij(J),
e E r, 5 E RC
{ 0,
(complement).
For a smooth boundary, c&J) = &; otherwise, c~ is given by cij(J) = %I+
&
[J (cpd- J (cp2)l ,
(47)
where cpi and cps are the angles between the xi-axis and the tangents tl and t2 at a point 5, and A(p = cpi - (~2 (see [lo, p. 1131). I is the identity matrix and J(v) is equal to (- sin 291 + sin 2~2) J(cp) = [ (cos 292 - cos 2qi)
(cos 292 -- cos 2Pi) (+ sin 2pi - sin 29s)
1’
The first term in equations (45),(46) comes from the given boundary conditions, expressing the tractions acting on the respective boundary represented by the vectors x, WE(E) = Tlfi(f)%(C),
c E rl,
k = 1,2
Low Reynolds
Number Deformation
109
and Wzi(E)
When equations (45),(46)
=
t E r2.
T24l)%(l),
are discretized for any boundary point <, they can be presented as
and
-T24W(J)
= -0
+ 4Cij(wq5)
-
(1 - X) 5
[I
d=l
Kj.(y,E)mi(y)dSV]
,
I E
r2,
(49)
I-A
with A4 and N the number of elements in Ii and I, respectively. The density functions Qj are unknown in the above equations. In order to evaluate the integrals involved, it was assumed that the density functions are approximated by a quadratic polynomial
where Ni , N2, Ns are the standard Lagrangian interpolation functions. Substituting equation (50) into equations (48) and (49) gives
P:(YWYW,~
t
E
(51)
h,
and
After performing the integrations over each element, the terms multiplying the same nodal value 6$ are assembled, producing the coefficients
where the last term of the above equation is only evaluated for points < E Ii. When h,&rd is calculated, the last node of each element also receives contributions of the integrals calculated for the first node of the next element. Equations (51) or (52) applied to a collocation point < can be written in matrix form as
(53)
where the indice P refers to the total number of boundary nodes and fit, by assembling the coefficients htj,Ed.
is the matrix formed
A. R. M. PRIMO et al.
110
In a more compact form, the above equation can be expressed as P H
c
=
(54
wt.
x=1
where
. HFx = HEX1
for x # E
and Hg = I&
+q.
When equation (54) is collocated at each node, it generates a 2N x 2N system of equations, which can be expressed in matrix form as H+=W.
(55)
Gauss elimination was utilized to solve the above system of equations where the density 9 is obtained. In order to find the surface velocities, equation (9) is discretized over the boundary (56) After applying the above equation to all boundary nodes, following the same procedures used to evaluate the integrals for the density function 9, a matrix L is formed which will be multiplied by the density function in order to provide the velocity at the boundary nodes u=L*. 3.2.
Calculation
of Surface
(57)
Tension
The surface tension driving the deformation of the fluid mass is expressed as products of the curvature by the normal, as given in boundary condition (3) and matching condition (5). The curvature and the normal vector were calculated at each nodal point by fitting a fourth-order Lagrangian polynomial. The point where the curvature/normal is calculated is considered the central point of the polynomial. Four more nodal points were taken for the fitting, two at each side of the central node. The normal was calculated by the expression (53) and the curvature was calculated by
(59) where
z:(o) =
$(Lrt-
z:'(O)= ;(-x; 3.3. Time
+16x;
- 304
+16x;
-x5).
Displacement
The deformation of the boundary was obtained by taking into account the Lagrangian representation of the fluid, giving the path followed by the particle AZ = 2(i’(r), E E I, (60) At where 7?(r) is the computed velocity at every boundary node and At is a fixed interval of time. A simple Euler scheme was applied to give the displacement of the fluid particle as the product of the time interval by the computed velocity. Although simple, this scheme proved to be stable for all cases studied.
Low Fkynolds Number Deformation
4. RESULTS
AND
111
DISCUSSION
This part of the present work simulates the deformation of systems formed by viscous drops embedded in a finite region of a fluid of different viscosity. The integral formulation developed in this paper was numerically solved through the Indirect Boundary Element Method (IBEM). Quadratic elements were used in all simulations. The physical parameters were made dimensionless, following [4]. A characteristic velocity 2cC, a characteristic pressure p,, and a characteristic time t, were defined as
where y is the surface tension coefficient and b the dynamic viscosity of the external fluid, and 1 is a characteristic length. Viscosity ratios in the range from X = 0.1 to X = 10.0 were chosen to generate the present results, what means that viscosities ten times lower (X = 0.1) to viscosities ten times higher (X = 10.0) than the external medium are taken into account. The case where the internal drop and the external medium have the same viscosity (X = 1.0) also received attention. This last case is apparently of lesser importance but it will prove to be interesting, especially when the internal drops are noncircular. At first sight, it appears that cases where the interior drops have the same viscosity as the exterior medium should not exhibit any deformation or displacement. However, the interior drop surface acts like a membrane, therefore susceptible to deformation. When not indicated otherwise, all the tests were carried out for external boundaries of elliptical (u = 0.75; b = 1.3) or circular form (r = 1.3). The interior viscous drops were circular discs (T = 0.2) or ellipses (o = 0.2, b = 0.4). For the numerical modelling, the external boundary was divided into 32 elements while the internal boundaries were divided into 16 elements each. Almost all the figures were produced by running the program from time t = 0.0 to t = 5.0 in steps of At = 0.1. The simplest example of deformation of multiple viscous forms is the case of a circular viscous disc containing a concentric circular viscous drop. This configuration has an exact solution provided in the theoretical modelling of the present integral formulation, leading to zero flow velocity. The numerical simulations through the IBEM confirmed the theoretical predictions, showing no deformations or displacements during the simulation for viscosity ratios in the range x = 0.1 to x = 10.0. Tables l-4 compare numerical results for the density function 3 with those predicted by equations (43),(44) for TI = Tz = 1. The results at the external boundary do not depend on X; the numerical values in Table 1 are for X = 10.0, although very similar results were obtained for all values of X tested. Tables 2-4 depict results at the internal surface for X = 0.1, X = 1.0, and X = 10.0, respectively. Very good agreement was achieved in all simulations. In the next series of tests, an ellipse was chosen as exterior domain. The interior domain continues to be a circular disc. Figure 2 describes the case when the external ellipse and the internal circle are concentric, with X = 10.0. The surface tension acting on the external boundary causes it to deform until the stable form of a circle is achieved. The interior circle does not move but it receives some influence of the deformation of the ellipse. When the circular drop is nonconcentric, it moves in different directions according to its initial position inside the domain. The direction of movement is dictated by the shape of the external region. Two positions were chosen to demonstrate the directions of movement inside the domain. Figure 3 shows a circle placed along the major axis of the ellipse. The circle does not deform but is forced to move towards the centre of the system. It is interesting to note that although moving while the ellipse is being transformed into a circle, the circular fluid disc stops moving before reaching the origin of the external region. In this figure, the viscosity ratio is x = 1.0.
A. R. M. PRIMO et al.
112
Table 1. Comparison between analytical (equation (43)) and numerical results for the density function along the external surface (X = 10.0).
r
Coordinates Xl
x2
T
-I-
Numerical @l
I
*pa
Ana @l
1 *z
ical
1.3000
0.0000
0.7692
0.0000
0.7692
0.0006
1.2937
0.1274
0.7655
0.0754
0.7655
0.0754
1.2750
0.2536
0.7544
0.1502
0.7544
0.1501
1.2440
0.3774
0.7361
0.2233
0.7361
0.2233
1.2010
0.4975
0.7107
0.2945
0.7106
0.2943
1.1465
0.6128
0.6784
0.3626
0.6784
0.3626
1.0809
0.7222
0.6397
0.4275
0.6396
0.4274
1.0049
0.8247
0.5946
0.4880
0.5946
0.4880
0.9192
0.9192
0.5440
0.5440
0.5439
0.5439
Table 2. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 0.1).
I
Coordinates Xl
Numerical a1
22
0.2000
0.0000
Analytical @2
92
56.8193
0.0000
o.oooo
-11.1000
0.1962
-0.0390
55.7275
-11.0849
0.1848
-0.0765
52.4578
-21.6566
-21.7730
0.1663
-0.1111
47.2435
-31.5671
-31.6207
0.1414
-0.1414
40.0825
-40.0825
-40.2446
Table 3. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 1.0).
Coordinates Xl
0.2000
Numerical
Analytical Ql
x2 o.oooo
5.0004
0.1962
-0.0390
4.9043
0.0000
-0.9755
5.0000
0.0000
4.9050
-0.9750
0.1848
-0.0765
4.6166
-1.9059
4.6200
-1.9125
0.1663
-0.1111
4.1577
-2.7781
4.1575
-2.7775
0.1414
-0.1414
3.5275
-3.5275
3.5350
-3.5350
Table 4. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 10.0).
r
1 Xl
@2
0.2000
0.0000
-0.1923
0.0000
0.1962
0.0375
-0.1886
0.0375
0.1848
0.0733
-0.1777
0.0735
0.1663
0.1068
-0.1599
0.1068
0.1414
0.1357
-0.1360
0.1360
Low Reynolds Number Deformation
113
Figure 2. A circular viscous disc concentric with the exterior elliptical viscous disc (X = 10.0).
Figure 3. A circular viscous disc along the major axis of the elliptical external region, initially placed near the boundary (X = 1.0).
Figure 4 outlines the deformation of an elliptical disc containing a circle placed on the first quadrant of the ellipse, with viscosity ratio X = 1.5. The circle is forced to move to the outside direction, following the enlargement of the minor axis of the exterior elliptical region. Figure 5 depicts an elliptical viscous disc containing six circular drops of radius r = 0.1 with X = 8.5 for each drop. The drops can be seen to be taking different directions until the motion ceases when the outer boundary becomes circular. In the case of a circular external boundary containing elliptical drops, the external boundary, although circular, will receive some effects of the deformation of the internal drops. Figure 6 presents a circular viscous region containing a concentrical elliptical drop with viscosity ratio X = 0.5. The stability of the exterior circle is slightly influenced by the deformations of the interior drop during the process. However, it recovers its initial shape when the internal drop
114
A. R. M. PRIMOet al.
Figure 4. A circular viscous disc placed in the first quadrant of the exterior elliptical viscoue dii (X = 1.5).
Figure 5. An elliptical viscous disc containing six internal circular drops (X = 8.5). achieves the final circular stable form. When the internal drop (X = 2.0) is placed adjacent to the external boundary (Figure 7) the instabilities are more visible alongside the boundary near the drop. It is interesting to note that the deformation of the ellipse is followed by a little displacement of its centre. Other positions inside the circle showed the same behaviour.
The following results concern the deformation of an elliptical viscous region containing ellip tical viscous drops. Here, more intense deformations and displacements take place, except for concentric ellipses. Starting from the simplest, case, Figure 8 displays two concentric elliptical discs (X = 0.1). As expected, although the inner ellipse achieves the circular form first, its stability is affected by the deformations of the outer ellipse, ss occurred in Figure 2. Figure 9 shows a set of four elliptical interior drops inside the exterior elliptical domain (X = 7.0). It is interesting to notice that the interior discs not only deform to a stable circle but also move while deforming.
Low Reynolds Number Deformation
Figure 6. A circular viscous region containing a concentrical elliptical drop (X = 0.5).
Figure 7. A circular viscous region containing a nonconcentrical elliptical drop (X = 2.0).
Figure 8. An external viscous elliptical region containing a concentric elliptical drop (X = 0.1).
115
116
A. R. M. PRIMO et al.
Figure 9. Four elliptical drops inside an external elliptical boundary (X = 7.0).
Figure 10. A circular viscous region containing two nonconcentric elliptical discs of different viscosity ratios. The drop on the left-hand side has X = 7.5 and X = 2.5 for the other drop.
The present boundary integral formulation also allows to consider drops of different viscosity ratios. Figure 10 shows the case of a circular drop containing two viscous elliptical discs of different viscosities, placed along the axis 22 = 0. The drop on the left-hand side has a viscosity ratio X = 7.5, while the drop on the right-hand side has a viscosity ratio X = 2.5. The drop on the right, lighter than the other, achieves its steady state faster. Finally, nine elliptical drops (a = 0.25, b = 0.10) of different viscosity ratios are placed inside a circular exterior domain (Figure 11). The central drop has X = 10.0 and the other drops have X = 1.0 or X = 5.0. The deformation of the drops is more evident than their motion. It is also clear that the rate of displacement and deformation depends on X. Figure 12 shows an exterior elliptical region (a = 2.0, b = 1.3) containing interior elliptical drops (a = 0.25, b = 0.10). The lighter drops achieve the stable circular form first. In this figure, the central drop has X = 5.0 and the others have X = 2.5 or X = 7.5. Again, because of the strong deformation of the exterior boundary, more intense motion of the drops is seen when compared to the previous figure. Figure 13 outlines the displacements of nine circular drops (T = 0.1) inside an elliptical domain (o = 2.0, b = 1.3). The central drop has X = 1.5 and the others have X = 6.5 or X = 9.5. A similar behaviour to Figure 5 is seen with lighter drops reaching the steady state faster than heavier drops.
Low EteynoIdsNumber Deformation
117
Figure 11. Nine elliptical drops of different viscosity ratios inside a circular external domain. The central drop has X = 10.0, while the other drops have X = 1.0 or x = 5.0.
Figure 12. An exterior elliptical region containing nine elliptical viscous drops of different viscosity ratios. The central drop has X = 7.5, while the other drops have x = 5.0 or x = 10.0.
Figure 13. Nine circuIar drops of different viscosity ratios inside an eliipticsl region. The centrai drop has X = 1.5, while the other drops have X = 6.5 or X = 9.5.
5. CONCLUSIONS This paper has presented a new integral equation formulation, in terms of only one surface potential, which provides unique solutions to the deformation of viscous drops in a bounded viscous flow region subject to surface tension. The formulation has been v&dated by application to several problems involving drops of arbitrary shapes and viscosity ratios. Numerical results obtained by the boundary element method displayed the expected behaviour, with convergence to circular steady forms.
118
A. R. M. PRIMO et al.
It has been shown that the formulation is of simple computational implementation, with the added term for completion of the system of integral equations requiring only extra regular integrations which are computed simultaneously with the standard integrations of the indirect boundary element scheme. Therefore, the additional computing time for the completed formulation is negligible.
REFERENCES 1. R. Clift, J.R. Grace and M.E. Weber, Bubbles, Drops and Particles, Academic Press, New York, (1978). 2. J.M. Rallison and A. Acrivos, A numerical study of the deformation and burst of a viscous drop in an extensional flow, Journal of Fluid Mechanics 89, 191-200, (1978). 3. H. Power, On the Rallison and Acrivos solution for the deformation and burst of viscous drops in an extensional flow, Journal of Fttid Mechanics 185, 547-550, (1987). 4. H.K. Kuiken, Viscous sintering: The surface-tension-driven flow of a liquid form under the influence of curvature gradients at its surface, Journal of Fluid Mechanics 214, 503-515, (1990). 5. G.A.L. van de Vorst, R.M.M. Mattheij and H.K. Kuiken, A boundary element solution for two-dimensional viscous sintering, Journal of Computational Physics 100, 50-63, (1992). 6. G.A.L. van de Vorst and R.M.M. Mattheij, Numerical analysis of a 2-D viscous sintering problem with non-smooth boundaries, Computing 49, 239-263, (1992). 7. A.R.M. Primo, L.C. Wrobel and H. Power, An indirect boundary element method for slow viscous Row in a bounded region containing air bubbles, Journal of Engineering Mathematics (to appear). 8. I.D. Chang and R. Finn, On the solution of a class of equations occurring in continuum mechanics, with application to the Stokes paradox, Arch. Rational Mech. Anal. 7, 388-401, (1961). 9. V. Botte and H. Power, A second kind integral equation formulation for three-dimensional interior flows at low Reynolds number, Boundary Elements Communications 6, 163-166, (1995). 10. F. Hartmann, Introduction to Boundary Elements, Springer-Verlag, Berlin, (1989).
Mathematical
and Computer
Modelling 31 (2000)
99-118 www.elsevier.nl/locate/mcm
Low Reynolds Number Deformation of Viscous Drops in a Bounded Flow Region under Surface Tension A. R. M.
PRIMO*
AND
L.
C.
WROBEL
Brunel University, Department of Mechanical Engineering Uxbridge, Middlesex UB8 3PH, United Kingdom H. POWER Wessex Institute of Technology, Ashurst Lodge Ashurst, Southampton SO40 7AA, United Kingdom Abstract-This paper presents a novel formulation of the completed indirect boundary element method to study the deformation of viscous drops on a slow viscous flow in a bounded region subject to surface tension effects. The theoretical background to the technique is explained in detail, including mathematical proofs of existence and uniqueness of solutions of the second-kind Fkedholm integral equations generated. The computer implementation of the formulation is simple and efficient. Numerical results of several tests are included. @ 2000 Elsevier Science Ltd. All rights reserved.
Keywords-Boundary
element method, Viscous drops, Stokes Aow, Surface tension.
1. INTRODUCTION The practical importance of studying the motion of drops and bubbles is due to their common occurrence in many industrial and biological systems, as well as in a number of technological processes. Mechanical engineers have studied droplet behaviour in connection with combustion operations, and bubbles in electromachining and boiling. Chemical and metallurgical engineers rely on drops and bubbles for such operations as distillation, absorption, flotation, and spray drying. The investigation of the mechanics of drops and bubbles has a long record in fluid dynamics, and still involves a substantial portion of pure and applied research. A general review of the subject can be found in [l]. Rallison and Acrivos [2] developed a numerical solution for the low Reynolds number deformation of a viscous drop, using Green’s integral representation formulae for the fluids inside and outside the drop. The external fluid region was infinite in extent, and the viscosity ratio between the two fluids was defined by X = p~/pr, with ~1 and ~2 the dynamic viscosities of the external fluid and the drop, respectively. A second-kind F’redholm integral equation was derived for the unknown surface velocity at the interface, for given external fluid velocity and interfacial surface tension. The interface deformation was found from a kinematic condition relating the rate of displacement of the interface to the normal velocity component. Although Rallison and Acrivos [2] only conjectured that their integral equation did not have any eigensolution for 0 < X < 00, since *On leave from Department
of Mechanical
Engineering,
Universidade
Federal de Pernambuco,
08957177/00/$ - see front matter @ 2009 Elsevier Science Ltd. Ail rights reserved. PII: SO895-7177(99)00218-S
Typeset
Recife, Brazil. by A&-‘B$
100
A. R. M. PRIMO et al.
their numerical solution encountered no difficulties for all values of X tested in this range, this fact was later proved analytically by (31. The main objective of the present work is to extend Rallison and Acrivos’ [2] formulation for the case when the external fluid region is finite and its exterior boundary is subject to known surface tension. For this problem, the cases when X ‘= 0 and X = 00 have a direct application for viscous sintering with air bubbles and solid inclusions, respectively; the cases when 0 < X < co could be seen as viscous sintering problems with fluid impurities. Integral equation formulations for viscous sintering using the direct boundary element method have been developed by [4-61. A more recent formulation for the shrinking of air bubbles in a sintering flow (i.e., X = 0) was proposed by [7]. The present approach is based on an indirect integral equation formulation for the secondkind boundary-value
problem in terms of only one surface potential.
Recasting the Stokes flow
equations in integral form produces a system of Fredholm integral equations of the second kind for the unknown surface density which is shown to possess eigenvalues corresponding to rigid body motions, restricted to the exterior surface. The eigenvalues are removed by the addition, to the corresponding integral equation, of a term which is linearly proportional to the rigid body motions. Theoretical proofs are included that the addition of this new term does not alter the solution of the problem, and that although the original integral equation and the completed one are identical, the former does not have a unique solution while the latter does. The computational implementation of the completed indirect BEM formulation is simple and efficient. Results of several tests are included to verify the theoretical and numerical procedures.
2. MATHEMATICAL
FORMULATION
A compound drop is considered consisting of a fluid drop with viscosity ~2 completely coated by another drop of viscosity ~1 (see Figure 1). The present analysis can be directly extended to a finite number of internal drops of different viscosities without any difficulty. Let 7~ be the interfacial surface tension between fluids 2 and 1, and yr between fluid 1 and the exterior environment. The governing equations for the fluid velocity ??and pressure p are the Stokes equations auln =O a 9 dXi
(1) where at&'" J+aXj
(
U$=-pmbij+j.Lm
dU7 dXi
)
,
(2)
with m = 1 for a point x E Rr and m = 2 when x E 02. The flow fields have to satisfy the following boundary conditions at every point 5 E I’r: ui'j (Ebj(E)
=
%K(t)ni(5)9
(3)
and the following matching conditions at every point c E I’s: u;(5) = G(C)
(4)
and (5) In the above set of equations, Ir is the exterior surface and l?s the interface between fluids 1 and 2, both assumed to be Lyapunov surfaces.
Low Reynolds Number Deformation
101
Figure 1. A viscous drop embedded in a finite viscous region.
The motion of the drops surface cannot be arbitrary since at each instant in time the force and torque yielded by the surface tractions have to satisfy the following compatibility conditions:
and aFjnjqf dS = 6,
s l-2 where 7,
1 = 1,2,3,
are the three linearly-independent ‘Pf =
417
cp (x) =
(22,
-3
(7)
rigid body vectors, taken as
for 1 = 1,2, -21).
Equation (7) follows from the fact that the flow field inside lY2 is a regular interior Stokes flow field yielding zero force and torque upon its surface container, while equation (6) is due to the balance of force and torque between the surfaces l?r and I’z, taking into account that the flow between these two surfaces has to be a regular Stokes flow field. Adding equations (6) and (7), using the boundary condition (3) and the matching condition (5), leads finally to the solvability condition
(8) where the normal vector ?? has been defined outwardly to the domain bounded by the corresponding surface boundary, i.e., I’1 and l?s. In order to guarantee continuity of the velocity field across the surface l?s, equation (4), the velocity field inside the domain fir U Cl2 will be defined by the following single-layer potential:
t&(x) = -
s rlur2
for every 3: E 521 or &. is given by
u{(z,Y)+j(Y)dSv,
The fundamental solution of the Stokes system of equations (Stokeslet)
A. R. M. PRIMO et al.
102
where T = )z - Yf. For convenience, the dynamic viscosity of each fluid will be defined in the corresponding pressure field. Therefore, at each fluid domain
(11) with
(12) and 2 = 1,2. Using the above representation for the flow field in each fluid region and letting a point z E Or approach a point < E I’r, taking the boundary condition (3) into consideration, the following equation is obtained:
where the kernel is given by Kji(Y4
= -;
1 (Yi - 6) (Yj - Ed(Yk - 5d
nh(E)
r4
Similarly, letting a point z E Ql approach a point < E l?z, an expression for the surface traction at such surface is obtained in the form
Letting a point x E R2 approach a point t E l?2
Subtracting equation (15) from equation (14), taking into account the matching condition (5), the following expression is obtained for a point < E l?2: -r2f4lh(t)
-
(111 -
cl2) [I
ra
Kji (Y, Wj
= -
(I4 ;
(Y) dSu +
p2) I&(<)
J rl
Kji(y, O@j(y)
dsv
1 .
Equations (13),(16) g ive a system of second-kind Fredholm integral equations for the unknown surface density 3 over I’i U r2. For convenience, these equations will be normalized with respect to the viscosity ~1. In this way, the equations can be rewritten as n49W)
=
;Q~(E)
-
E)Qj (Y) d& y Jrl Kj~(Y,Wj(Y)dSg +Jr2 Kji(Yv I
(17)
when < E I?r, and -T2fi(S)%(S)
=
--
(l; 3 o&g (18)
Kji (Y,t)cpjb) dsy +
J rl
when < E I’z, where 57 = TL/~I, 1 = 1,2, and A = p2/p1,
Kji(Y, Wj (Y>dsg
1
7
LOW hynolds
Number Deformation
103
According to F’redholm’s alternative, to analyse the uniqueness of solution of the above system of integral equations, the following homogeneous system has to be considered: &W
-
WY,
E)@;(Y) &,
+
J r2
K&I,
S)@;(Y) dS,
I
= 0,
(19)
for every < E rr, and
for every 6 E l?s. The above homogeneous system of integral equations implies that the flow field defined by
t&x) =
s r1
&,
Y>@~"(Y)ds,
(21)
UPa
and P’(Z) = w
I rlur2
q%c, y)@;(y) dS,,
(22)
where 1 = 1,2 when z E RI or z f Rs, respectively, and 3 is a nontrivial solution of the homogeneous system (19), (20), is a regular interior Stokes flow in Rr U &, having zero surface traction at the surface Fr and continuous surface traction across the surface l?z. Therefore, according to Green’s first identity for the Stokes equation in each flow region,
and (24) From the above two equations, it follows that
where the continuity of the surface traction across the surface I?2 has been used. Taking into account the zero surface traction at the surface I’r
Therefore, the flow field defined by equations (21),(22) has to behave as a rigid body motion -i_- - cpn, n = 1,2,3 (where continuity of the velocity field across the 1_ - 21 in s1r U 02: i.e., u surface r2 has been used), with an associated zero pressure all over sir U 522. This is due to the zero traction condition at I’1 and continuity of the surface traction across rs, since such flow field is the only solution of the Stokes flow
compatible with the given homogeneous conditions.
A. R. M. PRIMO et al.
104
From the continuity property of the velocity field defined by equation (21) across its densitycarrying surfaces, it follows that the limiting value of such velocity field when a point of the domain exterior to I’r, z E R,, tends to a point 5 E I’r, is also a rigid body motion, i.e., u&, The surface traction
Here, oij (uf,#)(,)
n = 1,2,3,
= CPXSL
for every < E l?r.
at I’r due to such external flow field is given by
is the limiting value of the stress oij(Z,p’)
rl,
(25)
when a point z tends to a
point < E coming from R,, with the auxiliary exterior domain filled with the same fluid of domain Sir. From equations (19),(26), it follows that Uij
(
f
u
,p
1
>
celnj
=
for every 5 E
-@p(<),
rl.
(27)
J?rom the uniqueness of solution of the exterior first boundary-value problem for the twodimensional Stokes system of equations satisfying the corresponding conditions at infinity and according to Chang and Finn’s [8] interpretation of the Stokes’ paradox, it follows that the vector density $
at the surface l?r has to be given by any of the eigensolutions of ;W)
-s,,
for every 5 E
Kji(Y, E)Q$ (Y) dSu = 9,
rl,
(28)
or any linear combination of them. The above homogeneous equation has precisely three independent linear solutions, which are generally unknown, since its corresponding adjoint homogeneous equation only admits the rigid body motions as nontrivial solutions. To prove the above, we consider the following single-layer potential:
Q(x) = J l-1
&xv Y)@;(Y) ds,
defining a regular Stokes flow in the interior or exterior of rr filled with the same fluid, with $ being any nontrivial solution of equation (28). The above singlelayer linearly independent flow fields in the domain interior to J?r, i.e., T(x)
1 = 1,2,3
= J(z),
and
potential
defines three
p(x) = 0,
(39)
according to the zero surface traction condition (28), yielding zero internal stresses, i.e.,
Uij
(T(X))= -i S, (Yi
-
xi) (Yj - Xj) (Y/C - Xk)
Q%(Y)ds,
r4
1
= 0,
(31)
for every x inside l?r. By the continuity pro erty of the velocity field due to a single-layer potential, it follows that the limiting value of 3 (x) when a point z E Q2, tends to a point < E I’r is also a rigid body motion, i.e., W)(e)
On the other hand, the surface traction oij(q)nj point z crosses l?r, given by aij
nj
($1
-
for every [ E rr.
1 = 1,2,3,
= 4(E),
U+,j
(32)
corresponding to (29) experiences a jump as a
(e)
nj
=
K&i,
1 = 1,2,3,
Low Reynolds Number Deformation
105
where the limiting value of the surface traction when a point belonging to the interior of l?l tends to a point on the surface I’l, i.e., rrij( V)(~nj, Therefore,
from the above equation,
is given by the left-hand side of equation
the surface traction
(28).
exerted upon I’1 by the exterior flow
field defined by the velocity (29) with its corresponding pressure is given by flij
(e)
nj =
-$,
According to Chang and Finn’s [8] interpretation
1 = 1,2,3. of the Stokes’ paradox, the above exterior
flow field representing the flow due to the motion of the surface rl as a rigid body, and having a logarithmic
behaviour at infinity, always yields a nonzero total force, i.e.,
I’$(Y)(P~(Y)Q, #
for 1 = 1,2,3
0,
and
n = 1,2.
(34)
* l-1
Comparing equations (25) and (27) with equations (32) and (33), taking into account the logarithmic behaviour of the flow at infinity and the uniqueness of solution for such unbounded flow field, it fol%ws that the vector density $ at every point < E I’1 has to be equal to the eigenfunctions Q’, 1 = 1,2,3, or any linear combination of them. It is easy to see that such eigenfunctions are the only nontrivial solutions of the homogeneous system (19),(20), since it follows from equation (31) that the second integral in equation (20) has to be identically zero, given that the curve I’2 is completely embedded in I’l, i.e.,
.I
K,,(Y,
W;(Y) dSy = --
nj (0
l-1
7r
s
(Yi -
l-1
Ei)(YJ - 0) (Yk - Ek)Q;(y) r4
=
dS
Y
o
1
for every < E r2. Therefore, equation (20) has been reduced to the homogeneous form of the adjoint of the secondkind integral equation proposed by Rallison and Acrivos [2] to solve the problem of viscous drop deformation, which is known to admit only a trivial solution when 0 < X < 00, i.e., G(e) E 0 for every < E l?z (see [3]). The cases X = 0 and X = co, for which the integral equation does not possess a unique solution, correspond to the problems of a gas bubble and a solid particle, respectively. In this way, equation (19) reduces to (28) showing that 8-i are the only eigensolutions of the homogeneous system (19),(20). To remove the above eigenvalues, which are yly linearly proportional
to the rigid body vector
restricted to the contour r‘l, a term which is
cp’, 1 = 1,2,3, is added to equation (17), i.e.,
where
It is easy to prove that the addition of this new term to equation (17) does not alter the nature of the problem since for any admittable nonhomogeneous term, i.e., satisfying the compatibility condition (6), the integral fll in equation (35) will end up being zero. A similar approach has been proposed by Botte and Power [9] in order to complete the deficient range of the doublelayer integral operator when the solution of the interior first-kind boundary-value problem for the Stokes system of equation is represented in terms of a double-layer potential alone. In such a case, the original second-kind Fredholm integral equation, coming from the jump property of the double-layer potential, is completed by adding a term linearly proportional to the surface normal vector. The coefficient of proportionality is identically equal to zero when the nonhomogeneous
A. R. M. PRIMO et al.
106
term, i.e., given surface velocity, satisfies the corresponding compatibility condition, i.e., the noflux condition.
Although for any admittable given surface velocity the original integral equation
and the completed one are identical, the former does not have a unique solution since its corresponding homogeneous equation has a nontrivial solution, whilst the completed one possesses a unique solution since its corresponding homogeneous equation admits only a trivial solution. Multiplying the expanded form of equation (17) by the rigid body vector, qP,
at a point on
the surface I’1 and integrating the resulting expression with respect to 5, i.e., evaluating the total force and torque upon the surface J?r, gives the following expression:
where the order of integration in the double integrals has been interchanged and the first two terms of the original expression are cancelled since
Kji(Y, J)$(S)
1 dSc dSg = 2
J r1
@j(~)Ipjp(~) d%
according to the adjoint of equation (28). Relation (36) can be reduced to -
Qj(Y)(pjp(Y>dSy +Pl J Jrl Jl-2 (P~P(Y)v;(Y)~% =
rl TI~Y)~~(Y)Y$‘(Y)~%
(37)
by using the property that a double-layer potential with a rigid body vector as density yields an interior flow field that behaves like a rigid body motion, i.e.,
J Kji(Y, t)(PP(t)dsc= rl
for every y inside rr.
'P;(Y),
For the following analysis, it is important to take notice that such double-layer identically equal to zero outside its density-carrying surface. Now, multiplying equation (18) by the rigid body vector qp
potential is
at a point of the surface I?2 and
integrating the resulting expression with respect to < gives -
J @i(t)dK)d% = r2
J r2
W4~)ni(O,d’(~) dsc,
(38)
where
for every y outside I?2 as noticed above. Therefore, from (37), (38), and the solvability condition (8), it follows that P41
‘P%$P~(<) d% = 9.
(39)
The above linear algebraic system for pr, &, ,Ba only admits the trivial solution, because the determinant of equation (39) has the term
J cpldds, r1
n,l = 1,2,3,
Low Reynolds Number Deformation
107
as element in its lth row and nth column, and is thus the Gram determinant for the vector functions 2,
2,
and 2
with a nonvanishing value, on account of the linear independence of 2,
k = 1,2,3. Having proved that the additional term in equation (17) does not alter the nature of the original problem, it can now be shown that the corresponding expanded homogeneous system, i.e., when the term
&)
s,, q(Y)Pg(Y) dS?J
is added to equation (19), only admits the trivial solution. From the above analysis it also follows that
J r1
Q;(Y)cP:(Y)
dsy = 0,
since the solvability condition is automatically
for n = 1,2,3,
satisfied when Ti = Tz = 0.
Thus, the new system always reduces to the original homogeneous system, with nontrivial solution $ = -7 ‘p , 1 = 1,2,3, when E E l?i and zero otherwise. In this way, the orthogonal relations (40) remove the eigensolutions according to the nonzero total force condition (34), in order to comply with the Stokes’ paradox. Therefore, it has been established that the expanded form of the nonhomogeneous system (17),(18) h as a unique continuous solution for any physical admittable data when 0 < X < 03. The above procedure can be directly expanded to the case of n different internal drops by expressing the density-carrying surface of the single-layer potential (9) in terms of the union of all the surfaces involved and using the corresponding pressure field. In this case, the solvability of the corresponding system of integral equations will be restricted to the values of the parameters O
The system of
(41)
and
TzCi -R;(l -- (1 - A) (1 + A)
+x)
[I rZ
Kji(Y, C)@j(Y) d% +
= -&(C)
J
Kji(y,
tP)j(y)
d&,
r1
1 ,
(42) t. E
I-2.
The unique solution of the above system is
(43) and (44) where the fact that every rigid body motion is orthogonal to the normal vector was used. This is a consequence of the rigid body motion eigensolution of the adjoint equation of equation (28), i.e., &(Y&#Y)
ds, = 0,
c E rl
and
1 = 1,2,3,
A. R. M. PRIMOet al.
108
and of the fact that the only nontrivial solution of the above homogeneous equation is the normal vector, i.e., L
Jr1 integrating over l?i, and using the adjoint equation
Multiplying the above equation by J, of (28), it follows that
n for 1 = 1,2,3.
Since the density vector is proportional to the normal vector at each of the surfaces, it follows that the flow velocity defined by the single-layer potential is equal to zero, as expected from this stable configuration. It is important to observe that the above solution also holds if the circles are nonconcentric. Similar conditions are encountered when n circular drops are embedded in a large circular drop. In all these cases, the resulting flow velocity is the state of rest.
3. NUMERICAL 3.1. Indirect
Boundary
Element
STMULATIONS
Method
The simulations carried out by using the Indirect Boundary Element Method (IBEM) were based on the expanded form of equations (17) and (18). I n a more generalized form, these equations can be written as
and -TsK(E)ni(E)
= -(I
+ %j(r)@j(E)
- (1 - A)
sr
Kji(g,Wj(Y)d%,
5cr2,
(46)
with r = rl ur2. The term C’ij corresponds to the characteristic term (see [lo, p. 166]), given as &j, Cij(6)
=
c E Q,
Cij(J),
e E r, 5 E RC
{ 0,
(complement).
For a smooth boundary, c&J) = &; otherwise, c~ is given by cij(J) = %I+
&
[J (cpd- J (cp2)l ,
(47)
where cpi and cps are the angles between the xi-axis and the tangents tl and t2 at a point 5, and A(p = cpi - (~2 (see [lo, p. 1131). I is the identity matrix and J(v) is equal to (- sin 291 + sin 2~2) J(cp) = [ (cos 292 - cos 2qi)
(cos 292 -- cos 2Pi) (+ sin 2pi - sin 29s)
1’
The first term in equations (45),(46) comes from the given boundary conditions, expressing the tractions acting on the respective boundary represented by the vectors x, WE(E) = Tlfi(f)%(C),
c E rl,
k = 1,2
Low Reynolds
Number Deformation
109
and Wzi(E)
When equations (45),(46)
=
t E r2.
T24l)%(l),
are discretized for any boundary point <, they can be presented as
and
-T24W(J)
= -0
+ 4Cij(wq5)
-
(1 - X) 5
[I
d=l
Kj.(y,E)mi(y)dSV]
,
I E
r2,
(49)
I-A
with A4 and N the number of elements in Ii and I, respectively. The density functions Qj are unknown in the above equations. In order to evaluate the integrals involved, it was assumed that the density functions are approximated by a quadratic polynomial
where Ni , N2, Ns are the standard Lagrangian interpolation functions. Substituting equation (50) into equations (48) and (49) gives
P:(YWYW,~
t
E
(51)
h,
and
After performing the integrations over each element, the terms multiplying the same nodal value 6$ are assembled, producing the coefficients
where the last term of the above equation is only evaluated for points < E Ii. When h,&rd is calculated, the last node of each element also receives contributions of the integrals calculated for the first node of the next element. Equations (51) or (52) applied to a collocation point < can be written in matrix form as
(53)
where the indice P refers to the total number of boundary nodes and fit, by assembling the coefficients htj,Ed.
is the matrix formed
A. R. M. PRIMO et al.
110
In a more compact form, the above equation can be expressed as P H
c
=
(54
wt.
x=1
where
. HFx = HEX1
for x # E
and Hg = I&
+q.
When equation (54) is collocated at each node, it generates a 2N x 2N system of equations, which can be expressed in matrix form as H+=W.
(55)
Gauss elimination was utilized to solve the above system of equations where the density 9 is obtained. In order to find the surface velocities, equation (9) is discretized over the boundary (56) After applying the above equation to all boundary nodes, following the same procedures used to evaluate the integrals for the density function 9, a matrix L is formed which will be multiplied by the density function in order to provide the velocity at the boundary nodes u=L*. 3.2.
Calculation
of Surface
(57)
Tension
The surface tension driving the deformation of the fluid mass is expressed as products of the curvature by the normal, as given in boundary condition (3) and matching condition (5). The curvature and the normal vector were calculated at each nodal point by fitting a fourth-order Lagrangian polynomial. The point where the curvature/normal is calculated is considered the central point of the polynomial. Four more nodal points were taken for the fitting, two at each side of the central node. The normal was calculated by the expression (53) and the curvature was calculated by
(59) where
z:(o) =
$(Lrt-
z:'(O)= ;(-x; 3.3. Time
+16x;
- 304
+16x;
-x5).
Displacement
The deformation of the boundary was obtained by taking into account the Lagrangian representation of the fluid, giving the path followed by the particle AZ = 2(i’(r), E E I, (60) At where 7?(r) is the computed velocity at every boundary node and At is a fixed interval of time. A simple Euler scheme was applied to give the displacement of the fluid particle as the product of the time interval by the computed velocity. Although simple, this scheme proved to be stable for all cases studied.
Low Fkynolds Number Deformation
4. RESULTS
AND
111
DISCUSSION
This part of the present work simulates the deformation of systems formed by viscous drops embedded in a finite region of a fluid of different viscosity. The integral formulation developed in this paper was numerically solved through the Indirect Boundary Element Method (IBEM). Quadratic elements were used in all simulations. The physical parameters were made dimensionless, following [4]. A characteristic velocity 2cC, a characteristic pressure p,, and a characteristic time t, were defined as
where y is the surface tension coefficient and b the dynamic viscosity of the external fluid, and 1 is a characteristic length. Viscosity ratios in the range from X = 0.1 to X = 10.0 were chosen to generate the present results, what means that viscosities ten times lower (X = 0.1) to viscosities ten times higher (X = 10.0) than the external medium are taken into account. The case where the internal drop and the external medium have the same viscosity (X = 1.0) also received attention. This last case is apparently of lesser importance but it will prove to be interesting, especially when the internal drops are noncircular. At first sight, it appears that cases where the interior drops have the same viscosity as the exterior medium should not exhibit any deformation or displacement. However, the interior drop surface acts like a membrane, therefore susceptible to deformation. When not indicated otherwise, all the tests were carried out for external boundaries of elliptical (u = 0.75; b = 1.3) or circular form (r = 1.3). The interior viscous drops were circular discs (T = 0.2) or ellipses (o = 0.2, b = 0.4). For the numerical modelling, the external boundary was divided into 32 elements while the internal boundaries were divided into 16 elements each. Almost all the figures were produced by running the program from time t = 0.0 to t = 5.0 in steps of At = 0.1. The simplest example of deformation of multiple viscous forms is the case of a circular viscous disc containing a concentric circular viscous drop. This configuration has an exact solution provided in the theoretical modelling of the present integral formulation, leading to zero flow velocity. The numerical simulations through the IBEM confirmed the theoretical predictions, showing no deformations or displacements during the simulation for viscosity ratios in the range x = 0.1 to x = 10.0. Tables l-4 compare numerical results for the density function 3 with those predicted by equations (43),(44) for TI = Tz = 1. The results at the external boundary do not depend on X; the numerical values in Table 1 are for X = 10.0, although very similar results were obtained for all values of X tested. Tables 2-4 depict results at the internal surface for X = 0.1, X = 1.0, and X = 10.0, respectively. Very good agreement was achieved in all simulations. In the next series of tests, an ellipse was chosen as exterior domain. The interior domain continues to be a circular disc. Figure 2 describes the case when the external ellipse and the internal circle are concentric, with X = 10.0. The surface tension acting on the external boundary causes it to deform until the stable form of a circle is achieved. The interior circle does not move but it receives some influence of the deformation of the ellipse. When the circular drop is nonconcentric, it moves in different directions according to its initial position inside the domain. The direction of movement is dictated by the shape of the external region. Two positions were chosen to demonstrate the directions of movement inside the domain. Figure 3 shows a circle placed along the major axis of the ellipse. The circle does not deform but is forced to move towards the centre of the system. It is interesting to note that although moving while the ellipse is being transformed into a circle, the circular fluid disc stops moving before reaching the origin of the external region. In this figure, the viscosity ratio is x = 1.0.
A. R. M. PRIMO et al.
112
Table 1. Comparison between analytical (equation (43)) and numerical results for the density function along the external surface (X = 10.0).
r
Coordinates Xl
x2
T
-I-
Numerical @l
I
*pa
Ana @l
1 *z
ical
1.3000
0.0000
0.7692
0.0000
0.7692
0.0006
1.2937
0.1274
0.7655
0.0754
0.7655
0.0754
1.2750
0.2536
0.7544
0.1502
0.7544
0.1501
1.2440
0.3774
0.7361
0.2233
0.7361
0.2233
1.2010
0.4975
0.7107
0.2945
0.7106
0.2943
1.1465
0.6128
0.6784
0.3626
0.6784
0.3626
1.0809
0.7222
0.6397
0.4275
0.6396
0.4274
1.0049
0.8247
0.5946
0.4880
0.5946
0.4880
0.9192
0.9192
0.5440
0.5440
0.5439
0.5439
Table 2. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 0.1).
I
Coordinates Xl
Numerical a1
22
0.2000
0.0000
Analytical @2
92
56.8193
0.0000
o.oooo
-11.1000
0.1962
-0.0390
55.7275
-11.0849
0.1848
-0.0765
52.4578
-21.6566
-21.7730
0.1663
-0.1111
47.2435
-31.5671
-31.6207
0.1414
-0.1414
40.0825
-40.0825
-40.2446
Table 3. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 1.0).
Coordinates Xl
0.2000
Numerical
Analytical Ql
x2 o.oooo
5.0004
0.1962
-0.0390
4.9043
0.0000
-0.9755
5.0000
0.0000
4.9050
-0.9750
0.1848
-0.0765
4.6166
-1.9059
4.6200
-1.9125
0.1663
-0.1111
4.1577
-2.7781
4.1575
-2.7775
0.1414
-0.1414
3.5275
-3.5275
3.5350
-3.5350
Table 4. Comparison between analytical (equation (44)) and numerical results for the density function along the internal surface (X = 10.0).
r
1 Xl
@2
0.2000
0.0000
-0.1923
0.0000
0.1962
0.0375
-0.1886
0.0375
0.1848
0.0733
-0.1777
0.0735
0.1663
0.1068
-0.1599
0.1068
0.1414
0.1357
-0.1360
0.1360
Low Reynolds Number Deformation
113
Figure 2. A circular viscous disc concentric with the exterior elliptical viscous disc (X = 10.0).
Figure 3. A circular viscous disc along the major axis of the elliptical external region, initially placed near the boundary (X = 1.0).
Figure 4 outlines the deformation of an elliptical disc containing a circle placed on the first quadrant of the ellipse, with viscosity ratio X = 1.5. The circle is forced to move to the outside direction, following the enlargement of the minor axis of the exterior elliptical region. Figure 5 depicts an elliptical viscous disc containing six circular drops of radius r = 0.1 with X = 8.5 for each drop. The drops can be seen to be taking different directions until the motion ceases when the outer boundary becomes circular. In the case of a circular external boundary containing elliptical drops, the external boundary, although circular, will receive some effects of the deformation of the internal drops. Figure 6 presents a circular viscous region containing a concentrical elliptical drop with viscosity ratio X = 0.5. The stability of the exterior circle is slightly influenced by the deformations of the interior drop during the process. However, it recovers its initial shape when the internal drop
114
A. R. M. PRIMOet al.
Figure 4. A circular viscous disc placed in the first quadrant of the exterior elliptical viscoue dii (X = 1.5).
Figure 5. An elliptical viscous disc containing six internal circular drops (X = 8.5). achieves the final circular stable form. When the internal drop (X = 2.0) is placed adjacent to the external boundary (Figure 7) the instabilities are more visible alongside the boundary near the drop. It is interesting to note that the deformation of the ellipse is followed by a little displacement of its centre. Other positions inside the circle showed the same behaviour.
The following results concern the deformation of an elliptical viscous region containing ellip tical viscous drops. Here, more intense deformations and displacements take place, except for concentric ellipses. Starting from the simplest, case, Figure 8 displays two concentric elliptical discs (X = 0.1). As expected, although the inner ellipse achieves the circular form first, its stability is affected by the deformations of the outer ellipse, ss occurred in Figure 2. Figure 9 shows a set of four elliptical interior drops inside the exterior elliptical domain (X = 7.0). It is interesting to notice that the interior discs not only deform to a stable circle but also move while deforming.
Low Reynolds Number Deformation
Figure 6. A circular viscous region containing a concentrical elliptical drop (X = 0.5).
Figure 7. A circular viscous region containing a nonconcentrical elliptical drop (X = 2.0).
Figure 8. An external viscous elliptical region containing a concentric elliptical drop (X = 0.1).
115
116
A. R. M. PRIMO et al.
Figure 9. Four elliptical drops inside an external elliptical boundary (X = 7.0).
Figure 10. A circular viscous region containing two nonconcentric elliptical discs of different viscosity ratios. The drop on the left-hand side has X = 7.5 and X = 2.5 for the other drop.
The present boundary integral formulation also allows to consider drops of different viscosity ratios. Figure 10 shows the case of a circular drop containing two viscous elliptical discs of different viscosities, placed along the axis 22 = 0. The drop on the left-hand side has a viscosity ratio X = 7.5, while the drop on the right-hand side has a viscosity ratio X = 2.5. The drop on the right, lighter than the other, achieves its steady state faster. Finally, nine elliptical drops (a = 0.25, b = 0.10) of different viscosity ratios are placed inside a circular exterior domain (Figure 11). The central drop has X = 10.0 and the other drops have X = 1.0 or X = 5.0. The deformation of the drops is more evident than their motion. It is also clear that the rate of displacement and deformation depends on X. Figure 12 shows an exterior elliptical region (a = 2.0, b = 1.3) containing interior elliptical drops (a = 0.25, b = 0.10). The lighter drops achieve the stable circular form first. In this figure, the central drop has X = 5.0 and the others have X = 2.5 or X = 7.5. Again, because of the strong deformation of the exterior boundary, more intense motion of the drops is seen when compared to the previous figure. Figure 13 outlines the displacements of nine circular drops (T = 0.1) inside an elliptical domain (o = 2.0, b = 1.3). The central drop has X = 1.5 and the others have X = 6.5 or X = 9.5. A similar behaviour to Figure 5 is seen with lighter drops reaching the steady state faster than heavier drops.
Low EteynoIdsNumber Deformation
117
Figure 11. Nine elliptical drops of different viscosity ratios inside a circular external domain. The central drop has X = 10.0, while the other drops have X = 1.0 or x = 5.0.
Figure 12. An exterior elliptical region containing nine elliptical viscous drops of different viscosity ratios. The central drop has X = 7.5, while the other drops have x = 5.0 or x = 10.0.
Figure 13. Nine circuIar drops of different viscosity ratios inside an eliipticsl region. The centrai drop has X = 1.5, while the other drops have X = 6.5 or X = 9.5.
5. CONCLUSIONS This paper has presented a new integral equation formulation, in terms of only one surface potential, which provides unique solutions to the deformation of viscous drops in a bounded viscous flow region subject to surface tension. The formulation has been v&dated by application to several problems involving drops of arbitrary shapes and viscosity ratios. Numerical results obtained by the boundary element method displayed the expected behaviour, with convergence to circular steady forms.
118
A. R. M. PRIMO et al.
It has been shown that the formulation is of simple computational implementation, with the added term for completion of the system of integral equations requiring only extra regular integrations which are computed simultaneously with the standard integrations of the indirect boundary element scheme. Therefore, the additional computing time for the completed formulation is negligible.
REFERENCES 1. R. Clift, J.R. Grace and M.E. Weber, Bubbles, Drops and Particles, Academic Press, New York, (1978). 2. J.M. Rallison and A. Acrivos, A numerical study of the deformation and burst of a viscous drop in an extensional flow, Journal of Fluid Mechanics 89, 191-200, (1978). 3. H. Power, On the Rallison and Acrivos solution for the deformation and burst of viscous drops in an extensional flow, Journal of Fttid Mechanics 185, 547-550, (1987). 4. H.K. Kuiken, Viscous sintering: The surface-tension-driven flow of a liquid form under the influence of curvature gradients at its surface, Journal of Fluid Mechanics 214, 503-515, (1990). 5. G.A.L. van de Vorst, R.M.M. Mattheij and H.K. Kuiken, A boundary element solution for two-dimensional viscous sintering, Journal of Computational Physics 100, 50-63, (1992). 6. G.A.L. van de Vorst and R.M.M. Mattheij, Numerical analysis of a 2-D viscous sintering problem with non-smooth boundaries, Computing 49, 239-263, (1992). 7. A.R.M. Primo, L.C. Wrobel and H. Power, An indirect boundary element method for slow viscous Row in a bounded region containing air bubbles, Journal of Engineering Mathematics (to appear). 8. I.D. Chang and R. Finn, On the solution of a class of equations occurring in continuum mechanics, with application to the Stokes paradox, Arch. Rational Mech. Anal. 7, 388-401, (1961). 9. V. Botte and H. Power, A second kind integral equation formulation for three-dimensional interior flows at low Reynolds number, Boundary Elements Communications 6, 163-166, (1995). 10. F. Hartmann, Introduction to Boundary Elements, Springer-Verlag, Berlin, (1989).