- Email: [email protected]
The Numerical Calculation of Resistance and Mobility Functions
293
A10. T H E N U M E R I C A L C A L C U L A T I O N OF R E S I S T A N C E AND M O B I L I T Y F U N C T I O N S
To conduct Stokesian dynamics simulations of non-dilute colloidal dispersions, we have to take into account hydrodynamic interactions among colloidal particles.
This means that the values
of resistance functions or mobility functions for arbitrary particle-particle separations are indispensable in advancing the time step in simulations.
Their tabulation as a function of
particle-particle separations is quite useful from the viewpoint of computation time. Appropriate values of the resistance and mobility functions can be obtained from an interpolation procedure with the tabulated data and are used in actual Stokesian dynamics simulations.
The data of the resistance and mobility functions for a system composed of two
equal spherical particles, which were numerically obtained, have already been shown in Sec. 5.3.
In this appendix, we outline the calculation method of the resistance and mobility
functions for a two rigid spherical particle system.
We concentrate our attention mainly on
the side of the methodology of the calculation method and not discuss in detail the mathematical aspect such as the derivation of each equation, since the mathematical manipulation is not straightforward and is beyond the objective of the present book. We use the coordinate system shown in Fig. A10.1.
The disturbance velocity at a point x in
Fig. A10.1 is expressed using Lamb's general solution as [1-3]
"='
~{
. . . . . .
'
n(2n - l ) r, pi',~_, - 2n(2n -1) r, Vp[~_,
(n+,) r , p,:,
_,,_
+ ,--, Vcl)~2' , + V • (rzz'_~-'_,)+ n(2n - 1) -
-n-
--
(A10.1)
(n-2) , Vp.....,},
2n(2n
1)
2
_[2)
in which v | is the ambient velocity field, r~=x-x~, r~= ]r~ ]for a=l,2, and the spherical harmonics p~), . . .cI)C~'l and x.c~) . .... ~are expanded as follows
p~)_, = ~ r - " -a ' p ' ( c o s O a ) { a ' a ' 6O0n m
+ a'~' sinmcp} , rnn
m=O
*~)_, = ~r-"-'p~,'(cosO~){b '~ '~' sin mop}, ~ 0,,)fi0m + b m,
(A10.2)
m =0
=
a
)Cm,, COSm~o.
m =0
Equation (A10.1) with Eq. (A10.2) satisfies the Stokes equation (4.20) and the equation of continuity (4.6).
If the unknown constants of am,,. - ~ ' b ,l~ ' , and Cm, -~) are determined properly,
the solution (A10.1) can satisfy the boundary condition on the particle surfaces.
For each
294
r~
rj
Z z
y/
-Y
jr
Figure A10.1 A coordinate system for two spherical particle interaction.
resistance function, only one particular value of m is actually required.
The final goal of the
present discussion is, therefore, to obtain algebraic equations for the unknown coefficients, whereby the resistance and mobility functions can be calculated numerically. The disturbance velocity at the particle surface, v,, is equal to the difference between the imposed ambient velocity field and the velocity of the particle motion, and is expressed for the boundary condition at the surface of particle 1 as
v~= ~-'~-'[V{r[PT(cosO,)(A,o6om+ A,,, sin m(p,)}
(A10 3)
/=1 m = 0
+v x
cos
, }],
in which PF' (cos0~) is the associated Legendre function, (0~, q)j) are the zenithal and azimuthal angles.
The appropriate choice of the coefficients At~ and Bern gives all surface velocities
associated with disturbance fields; for example, setting all coefficients zero except A~0=1 gives a translational velocity along the particle-particle axis [1,2]. The use of the cylindrical coordinate system (z,p,q)), shown in Fig A10.1, makes the mathematical manipulation more straightforward.
The solution for the disturbance velocity,
Eq. (A10.1), has to satisfy the boundary condition (A10.3).
This requirement concerning the
boundary condition leads to the following algebraic equations which must be satisfied by the unknown coefficients, a ~) mn
,
b~, ) and c I~" ,
mn
~--'[a~ (-S) ~-' (n+l) _,, p~ (n-2) _, 2 m ,,=, ~=, 2(2n - 1) r~ ~: (~:~ 1- 2n(2n - 1) r~ (1 - ~:~ )(P,, (~:~11'
b2'~,~f'(-S) ~-' l- (n+l)r~ ' a=l
-n-')
m
"~:~P,, (~:~)+r~
-n-2
2
(I-~:~)(P/
m
{l) (~:a))' }+mCmn
(A10.4)
295
n=l
u,..
x t~~ , x
2(2n
4=1
--
1) r4 sin 0~ P~ (~:4) + 2n(2n
(~)+msinO~P,
(~)
}}
+b':>.,
1) r4
r4
a=l
(n+l)sinO4P/~ ( ~ ) - ~ P , ; ' + ' ( ~ ) - r n s i n O ~ P / ~
(~)
c~,
r~
(~)
4=1
= A~., {/sin0,P+: (~:,)- . (~,p.,+l , (~,)+msin0,P+ (~, >)}- B,., p"+'.` ,..,rr-', (A10.5)
m ,,,u~.
Ps ( ~ ) ( s i n 0
t sin0~(P,~ ( ~ ) ) '
+c,..
4=1
4=-1
= A, mP+" (~,) / sin 0, + B,~ sin 0, (P~" (4:,))'.
(A10.6) in which ~=cos0a and S is the symmetry parameter; S=1 for problems with mirror symmetry and S= -1 for problems with mirror anti-symmetry [1,2].
Equation (A10.4) comes directly
from the z-component boundary condition, and Eq. (A10.6) from the {p-component boundary condition.
Equation (A10.5) is obtained by subtracting the q~-component equation from the p-
component equation.
. . . . . b~,,~,~ and % , can be determined by the algebraic If the constants a .~'~ _~4~
equations (A10.4) to (A10.6), then the solution for the velocity field for the two equal-particle system may be obtained. It is seen from Eq. (13.7) that the velocity field can be expanded as a function of forces, torques, stresslets, etc.
The resistance functions can. therefore, be correlated to Um.~4~, b,.,,(4~.
and c<,.=.) by comparing Eq. (A10.1) with the equation expressed by forces, torques, stresslets. etc.. such as Eq. (13.7).
The final results for the resistance functions are as follows:
296
X{42" = -~1 {a0~ (1,-1) - a0~ (1,1)},
1 (1,-1) + a0, (1,1)}, X~" = ~-{a0,
1
Yzf" =-~{a,,1 (1,-1) + a,, (1,1)}, Y~" - -~-{a,, (1,-1)- a, , (1,1)}, Y~" : -{c,, (1,-1) + c,, (1,1)}, Y~" : C,l (1,-1) - c,, (1,1),
1 X~* = - ~1 {a02(1,-1)+ ao2 (1,1)}, X~" = -~-{ao2 (1,-1)- %2 (1,1)},
(A10.7)
y~. = -~-{a,2 1 (1,-1) + a,2 (1,1)}, Y~"" 1 (1,1)}, 2 = -4{a,2(1,-1)-a,2
X,'#"
1 (2,1), + X;7 . = -f-dao2
1
y(~r.+ yff. = i-O a'2 (2,-1),
1 Z,#" + Z,#" = -{-~ a 22 (2,1),
in which the superscript (1) of a(l~, b~, and-(" C mn has been dropped for simplicity, and they are a function of l and S, expressed as a,,,(l,S), bm,(l,S), and c,,,(l,S).
The results shown in Eq.
(A10.7) have been obtained for the set of At,,=l and B~,,,=O. From the other set of Btm=l and A~m=0, the remaining resistance functions are derived as
1 (1,-1) + c0~(1,1)}, X c• = -~{c0~
X(2*= -~{Coz 1 (1,-1)- Co~(1,1)},
1
,.
=-~{c~(1,-1)+c~(1,1)}, y~zU. =
1
Y~2 =
--~{a,2(1,-1)+a,2(1,1)},
{c, , (1,-1) - c, . (1,1)},
(110.8)
1 Y~" = - { a , 2 (1,- 1) - a,z (1,1)}.
These equations of Eq. (A10.7) and (A10.8) complete the calculation of the resistance functions. To solve the algebraic equations of Eqs. (A10.4) to (A1 0.6), the summation of n has to be truncated at N terms.
We consider a special case of l=m=S=l, Arm=l, and B~=0, to understand
the way of solving these equations.
In this case, the 3N unknown coefficients, a~, a~2, ..., a~N,
b ~ , bl2 , ...,
b~x, c~, c~2.... , c~,v, have to be determined by the truncated version ofEqs. (A10.4) to
(A10.6).
Thus N points at the particle surface, at which the solution of the velocity field, Eq.
(A10.1), satisfies the boundary condition, are required.
If we take such points as 0~=krd(N-1)
(k=-0,1,...,N-1), we can get the 3N algebraic equations necessary for calculating the 3N unknown coefficients. Finally, it should be noted that more boundary points (that is, a larger value of N) are necessary as the two particles approach to a nearly touching situation; in other words, this means that the solutions do not converge rapidly as N increases, when the particles are almost
297
touching.
Once the resistance functions from these numerical procedures are obtained, the
mobility functions are calculated from the resistance functions using the relationships between the resistance and mobility functions which have already been given in Sec. 5.3.
A sample
FORTRAN program for calculating the resistance functions is shown in Appendix A11.7. References
1. S. Kim and R. T. Mifflin, The Resistance and Mobility Functions of Two Equal Spheres in Low-Reynolds-Number Flow, Phys. Fluids, 28 (1985) 2033 2. S. Kim and S. J. Karrila, "Microhydrodynamics," Butterworth-Heinemann, Stoneham, 1991 3. P. Ganatos, R. Pfeffer, and S. Weinbaum, A Numerical-Solution Technique for ThreeDimensional Stokes Flows with Application to the Motion of Strongly Interacting Spheres in a Plane, J. Fluid Mech., 84 (1978) 79
A10. T H E N U M E R I C A L C A L C U L A T I O N OF R E S I S T A N C E AND M O B I L I T Y F U N C T I O N S
To conduct Stokesian dynamics simulations of non-dilute colloidal dispersions, we have to take into account hydrodynamic interactions among colloidal particles.
This means that the values
of resistance functions or mobility functions for arbitrary particle-particle separations are indispensable in advancing the time step in simulations.
Their tabulation as a function of
particle-particle separations is quite useful from the viewpoint of computation time. Appropriate values of the resistance and mobility functions can be obtained from an interpolation procedure with the tabulated data and are used in actual Stokesian dynamics simulations.
The data of the resistance and mobility functions for a system composed of two
equal spherical particles, which were numerically obtained, have already been shown in Sec. 5.3.
In this appendix, we outline the calculation method of the resistance and mobility
functions for a two rigid spherical particle system.
We concentrate our attention mainly on
the side of the methodology of the calculation method and not discuss in detail the mathematical aspect such as the derivation of each equation, since the mathematical manipulation is not straightforward and is beyond the objective of the present book. We use the coordinate system shown in Fig. A10.1.
The disturbance velocity at a point x in
Fig. A10.1 is expressed using Lamb's general solution as [1-3]
"='
~{
. . . . . .
'
n(2n - l ) r, pi',~_, - 2n(2n -1) r, Vp[~_,
(n+,) r , p,:,
_,,_
+ ,--, Vcl)~2' , + V • (rzz'_~-'_,)+ n(2n - 1) -
-n-
--
(A10.1)
(n-2) , Vp.....,},
2n(2n
1)
2
_[2)
in which v | is the ambient velocity field, r~=x-x~, r~= ]r~ ]for a=l,2, and the spherical harmonics p~), . . .cI)C~'l and x.c~) . .... ~are expanded as follows
p~)_, = ~ r - " -a ' p ' ( c o s O a ) { a ' a ' 6O0n m
+ a'~' sinmcp} , rnn
m=O
*~)_, = ~r-"-'p~,'(cosO~){b '~ '~' sin mop}, ~ 0,,)fi0m + b m,
(A10.2)
m =0
=
a
)Cm,, COSm~o.
m =0
Equation (A10.1) with Eq. (A10.2) satisfies the Stokes equation (4.20) and the equation of continuity (4.6).
If the unknown constants of am,,. - ~ ' b ,l~ ' , and Cm, -~) are determined properly,
the solution (A10.1) can satisfy the boundary condition on the particle surfaces.
For each
294
r~
rj
Z z
y/
-Y
jr
Figure A10.1 A coordinate system for two spherical particle interaction.
resistance function, only one particular value of m is actually required.
The final goal of the
present discussion is, therefore, to obtain algebraic equations for the unknown coefficients, whereby the resistance and mobility functions can be calculated numerically. The disturbance velocity at the particle surface, v,, is equal to the difference between the imposed ambient velocity field and the velocity of the particle motion, and is expressed for the boundary condition at the surface of particle 1 as
v~= ~-'~-'[V{r[PT(cosO,)(A,o6om+ A,,, sin m(p,)}
(A10 3)
/=1 m = 0
+v x
cos
, }],
in which PF' (cos0~) is the associated Legendre function, (0~, q)j) are the zenithal and azimuthal angles.
The appropriate choice of the coefficients At~ and Bern gives all surface velocities
associated with disturbance fields; for example, setting all coefficients zero except A~0=1 gives a translational velocity along the particle-particle axis [1,2]. The use of the cylindrical coordinate system (z,p,q)), shown in Fig A10.1, makes the mathematical manipulation more straightforward.
The solution for the disturbance velocity,
Eq. (A10.1), has to satisfy the boundary condition (A10.3).
This requirement concerning the
boundary condition leads to the following algebraic equations which must be satisfied by the unknown coefficients, a ~) mn
,
b~, ) and c I~" ,
mn
~--'[a~ (-S) ~-' (n+l) _,, p~ (n-2) _, 2 m ,,=, ~=, 2(2n - 1) r~ ~: (~:~ 1- 2n(2n - 1) r~ (1 - ~:~ )(P,, (~:~11'
b2'~,~f'(-S) ~-' l- (n+l)r~ ' a=l
-n-')
m
"~:~P,, (~:~)+r~
-n-2
2
(I-~:~)(P/
m
{l) (~:a))' }+mCmn
(A10.4)
295
n=l
u,..
x t~~ , x
2(2n
4=1
--
1) r4 sin 0~ P~ (~:4) + 2n(2n
(~)+msinO~P,
(~)
}}
+b':>.,
1) r4
r4
a=l
(n+l)sinO4P/~ ( ~ ) - ~ P , ; ' + ' ( ~ ) - r n s i n O ~ P / ~
(~)
c~,
r~
(~)
4=1
= A~., {/sin0,P+: (~:,)- . (~,p.,+l , (~,)+msin0,P+ (~, >)}- B,., p"+'.` ,..,rr-', (A10.5)
m ,,,u~.
Ps ( ~ ) ( s i n 0
t sin0~(P,~ ( ~ ) ) '
+c,..
4=1
4=-1
= A, mP+" (~,) / sin 0, + B,~ sin 0, (P~" (4:,))'.
(A10.6) in which ~=cos0a and S is the symmetry parameter; S=1 for problems with mirror symmetry and S= -1 for problems with mirror anti-symmetry [1,2].
Equation (A10.4) comes directly
from the z-component boundary condition, and Eq. (A10.6) from the {p-component boundary condition.
Equation (A10.5) is obtained by subtracting the q~-component equation from the p-
component equation.
. . . . . b~,,~,~ and % , can be determined by the algebraic If the constants a .~'~ _~4~
equations (A10.4) to (A10.6), then the solution for the velocity field for the two equal-particle system may be obtained. It is seen from Eq. (13.7) that the velocity field can be expanded as a function of forces, torques, stresslets, etc.
The resistance functions can. therefore, be correlated to Um.~4~, b,.,,(4~.
and c<,.=.) by comparing Eq. (A10.1) with the equation expressed by forces, torques, stresslets. etc.. such as Eq. (13.7).
The final results for the resistance functions are as follows:
296
X{42" = -~1 {a0~ (1,-1) - a0~ (1,1)},
1 (1,-1) + a0, (1,1)}, X~" = ~-{a0,
1
Yzf" =-~{a,,1 (1,-1) + a,, (1,1)}, Y~" - -~-{a,, (1,-1)- a, , (1,1)}, Y~" : -{c,, (1,-1) + c,, (1,1)}, Y~" : C,l (1,-1) - c,, (1,1),
1 X~* = - ~1 {a02(1,-1)+ ao2 (1,1)}, X~" = -~-{ao2 (1,-1)- %2 (1,1)},
(A10.7)
y~. = -~-{a,2 1 (1,-1) + a,2 (1,1)}, Y~"" 1 (1,1)}, 2 = -4{a,2(1,-1)-a,2
X,'#"
1 (2,1), + X;7 . = -f-dao2
1
y(~r.+ yff. = i-O a'2 (2,-1),
1 Z,#" + Z,#" = -{-~ a 22 (2,1),
in which the superscript (1) of a(l~, b~, and-(" C mn has been dropped for simplicity, and they are a function of l and S, expressed as a,,,(l,S), bm,(l,S), and c,,,(l,S).
The results shown in Eq.
(A10.7) have been obtained for the set of At,,=l and B~,,,=O. From the other set of Btm=l and A~m=0, the remaining resistance functions are derived as
1 (1,-1) + c0~(1,1)}, X c• = -~{c0~
X(2*= -~{Coz 1 (1,-1)- Co~(1,1)},
1
,.
=-~{c~(1,-1)+c~(1,1)}, y~zU. =
1
Y~2 =
--~{a,2(1,-1)+a,2(1,1)},
{c, , (1,-1) - c, . (1,1)},
(110.8)
1 Y~" = - { a , 2 (1,- 1) - a,z (1,1)}.
These equations of Eq. (A10.7) and (A10.8) complete the calculation of the resistance functions. To solve the algebraic equations of Eqs. (A10.4) to (A1 0.6), the summation of n has to be truncated at N terms.
We consider a special case of l=m=S=l, Arm=l, and B~=0, to understand
the way of solving these equations.
In this case, the 3N unknown coefficients, a~, a~2, ..., a~N,
b ~ , bl2 , ...,
b~x, c~, c~2.... , c~,v, have to be determined by the truncated version ofEqs. (A10.4) to
(A10.6).
Thus N points at the particle surface, at which the solution of the velocity field, Eq.
(A10.1), satisfies the boundary condition, are required.
If we take such points as 0~=krd(N-1)
(k=-0,1,...,N-1), we can get the 3N algebraic equations necessary for calculating the 3N unknown coefficients. Finally, it should be noted that more boundary points (that is, a larger value of N) are necessary as the two particles approach to a nearly touching situation; in other words, this means that the solutions do not converge rapidly as N increases, when the particles are almost
297
touching.
Once the resistance functions from these numerical procedures are obtained, the
mobility functions are calculated from the resistance functions using the relationships between the resistance and mobility functions which have already been given in Sec. 5.3.
A sample
FORTRAN program for calculating the resistance functions is shown in Appendix A11.7. References
1. S. Kim and R. T. Mifflin, The Resistance and Mobility Functions of Two Equal Spheres in Low-Reynolds-Number Flow, Phys. Fluids, 28 (1985) 2033 2. S. Kim and S. J. Karrila, "Microhydrodynamics," Butterworth-Heinemann, Stoneham, 1991 3. P. Ganatos, R. Pfeffer, and S. Weinbaum, A Numerical-Solution Technique for ThreeDimensional Stokes Flows with Application to the Motion of Strongly Interacting Spheres in a Plane, J. Fluid Mech., 84 (1978) 79