Analytical benchmark solutions for nanotube flows
with variable viscosity
I. V. Makeev, I. V. Blinova, I. Yu. Popov
ITMO University, 49, Kronverkskiy, St. Petersburg, 197101, Russia
popov1955@gmail.com
PACS 68.90.+g, 05.60.-k DOI 10.17586/2220-8054-2015-6-5-672-679
Three-dimensional Stokes equations with variable viscosity in cylindrical coordinates are considered. This case is natural for flow through a nanotube in biological applications. We obtain exact particular solutions -a benchmark for numerical approache.
Keywords: Stokes flow, variable viscosity, multigrid methods, benchmark solutions. Received: 10 October 2015
1. Introduction
Flow through nanotubes are intensively investigated, particularly due to its biological applications, e.g., as a transport channel for some compounds into the cell (see, e.g., [1]). Although there is no general approach in nanohydrodynamics [2-6], the Stokes approximation (with variable viscosity) is available in many cases. We consider the axisymmetric case, which is natural for the nanotube flow. Biologists need a reliable computational approach to describe and predict the mass transport through a nanotube "to" and "from" the cell. There are several computational schemes for solution of the Stokes and continuity equations, but in the case of strongly varying viscosity, these schemes are plagued by difficulties. One would like to have an instrument for testing of these algorithms. One convenient method is to suggest benchmark solutions for some particular cases [7-9]. In this paper, such a benchmark solution is constructed for the axisymmetric case. It is interesting to note that the same mathematical problem appears in geophysics [10, 11].
This paper deals with the Stokes and continuity equations for the case of variable viscosity and density having the following form:
(V- a) = -pG , (1)
V(pv) = 0, (2)
where v is velocity, n is a dynamic viscosity, a is the total stress tensor, p is a pressure, G is a gravitational force. We consider equations (1), (2) in cylindrical coordinates (r, <^,z) and construct a solution for the case when the variables depend only on the radius r.
2. Benchmark solutions
We construct particular solutions of the system of Stokes and continuity equations for specific density and viscosity distributions: n = n(r), P = p(r). Let vr = vr(r), vv = v^(r), vz = vz(r), P = P(r), n = n(r), P = p(r), G = G(r). Then equation (1) simplifies considerably:
2n1 vr + 2n'vr + 2nvr' - 2n-jvr - P = -pGr (3)
v'v'y—1 n'v^ + nv;; + ^n^- n^vv = —pg^ (4)
n1 + n'vZ + nvZ' = -PGz. (5) r
The continuity equation takes the form:
p1 Vr + p'v- + pvr = 0. (6) r
In this case, we obtain the following solutions of equations (1), (2):
Vr = rp, v^ = cif (r) + C2r + Ci(r)f (r) + C2(r)r,
r r2
Vz = - / ¿2(/ ripGzdri + Ci)dr2 + c2, 1 ' i
(7)
where
[' 1 1
P(r) = (pGr + + + 2nvr - vr)dr, (8)
r
[11 1 /(r) = expH (— + — --)dr2),
1 r2 nr2 /2 A dn + C
1 J r/r 1 1
C'(r) = / //) dr, C2(r) = - / //
Formulas (7), (8) gives us the solution of equations (1), (2). Derivation of the solution is presented in Appendix A in detail.
3. Multigrid method
In this work, we derive a procedure of multigrid method for solving the Stokes equations with variable viscosity in cylindrical coordinates. The analogous algorithm for the Cartesian coordinates was described in [10]. We derive a similar scheme for cylindrical coordinates.
As usual, the multigrid method algorithm contains three steps: 1) smoothing operation 2) restriction operation 3) prolongation operation. Cylindrical coordinates are orthogonal coordinates, thus the implementation of the prolongation and restriction operations in our method is not different from that in the case of the Cartesian coordinates. As for the smoothing operation, it differs. This procedure is described in Appendix B in detail.
The scheme for algorithm testing is as follows. Let us consider some particular analytical solutions (7), (8): vr = — i, vv = r2, vz = — pGz r5 — i r3, P (r) = pGr r + 1.2 in the domain 1 ^ r ^ 2, 0 ^ 0 ^ 1,0 ^ z ^ 1 (p = const, Gr = 10, Gv = 0, Gz = 10, n = r-3). We calculate the values for velocity and pressure given by our analytical solution and take these values as the boundary conditions for the numerical algorithm. The deviation of the numerical solution values from the analytical solution is related to the error of the multigrid scheme. The dependence of the relative error on the grid step for the multigrid scheme is shown in Figure.1. Positive curve slopes indicate a convergence for the algorithm.
Fig. 1. Error norm via the grid resolution in logarithmic scale: blue lines-pressure, red lines- vr , black lines - vv , green lines - vz; solid lines - Li-error, dashed lines- L^ -error, solid dotted lines- L2-error
Appendix A. Derivation of the Stokes equations solution
Integration of (6) gives us:
vr = —. (9)
rp
Substitution of (9) in (3) gives us the expression for pressure:
f 1 1
P(r) = (pGr+2rç-v; + v; + 2^v;' - vr)dr. (10)
J r r 2
Consider equation (4). Simple transformation gives us:
vy + v^( n + r) - r n + ) = T" (11)
We seek a solution of (11) in the form vy(r) = vy0 (r) + vy1(r), where vy0(r) is particular solution of (11) and vy1(r) is the general solution of the corresponding homogeneous equation:
,W' 1s ( 1 n 1s / s
vy + vy(n + ^ - v( 1 n + "2) = °. (12)
Let us make a replacement:
v'
y
v^
(13)
in (12). It transforms into the first-order equation:
z' + z(^ + !) + z2 - 1 'n- - -1- = 0. (14)
n r r n r 2
If one makes the substitution z = 1 y, where y = y(r) then (14) transforms into the Riccati equation:
y' + y2- - n' - - = 0. (15)
nr r
z
Function y = n is a particular solution of (15). Keeping this in mind, we seek the solution in the form y = u + n. For function u = u(r) one obtains the Bernoulli equation:
u' = —21 u — — u2. (16)
r nr
Note that u(r) = 0 is a solution. We seek a non-trivial solution of (16) in the form u = q , where q = q(r). The equation simplifies:
1
q + — q2 = 0,
giving us:
Coming back to (12), we get:
After integration, one obtains:
q - r i_
r3 1
M3 dr
J nr3
vi -!(,, + i
v^ nr r2 / "Ts dr
r
1 1 1
ln(K|)- / (- + ZT3 ^-)dr2 + C21
1 r2 nr3 / nTs dri + C11
3 (
1 nr3
Solution u(r) = 0 of (16) leads to the following solution of (12):
v^ - cr
As a result, we obtain the general solution of (12) in the following form:
v^i(r) = cif (r) + C2r,
r
where /(r) = exp^ (+ -Jg -F2-i-)dr2).
i -tg dri+C
1 nrg
We seek a particular solution of equation (11) by the Lagrange method in the form v^o(r) = Ci (r)/(r) + C2(r)r.
Here Ci(r), C2(r) satisfy the system of equation:
Ci / (r) + C2 r = 0
Ci /'(r) + C2 = ^ .
^ir^p^mg^ Ci(r) = / dr, C2(r) = — / n(/-%dr
As a result, we come to the following solution of equation (4):
Mr) = ci / (r) + C2r + Ci(r)/(r) + C2 (r)r, (17)
we let vz = w(r), equation (8) takes the form:
w' + (n' + 1)w = — pGz (18)
n r n
If one makes the substitution w = wnr, where w = w(r), then (11) transforms into equation:
w' 1 pg
w — =--Gz
nr n
1
Integration gives us:
w =--( / ripGzdri + ci),
nr J i
As a result, we get the following expression for the velocity component:
r T2
Vz = - / -( ripGzdri + Ci)dr2 + C2.
J nr2 J ii
(19)
Appendix B. Gauss-Seidel smoother
Smoothing operation can be implemented on the basis of Gauss-Seidel iterations. The respective iterative pressure and velocity update schemes for a regularly spaced grid can be derived:
T-) . * Tycontinuity ^continuity /nn\
P(i,j,l) + nn(i,j,l)^Rijl Vrelaxation (20)
rynew P(i,j,l)
a ^continuity _ ^continuity vr dvr 1 dvp dvz
ij = Rij - 7 - TT - - ,
öz
r(i,j,l)
V
<t>(i,j,l)
z(i,j,l)
A Rr-Stokes I ^Ri,j,l vStokes
vr(i,j,l) + ^ Vrelaxation
Cvr (i,j,l)
a r,^-Stokes + ^Ui,j,l vStokes
vp(i,j,l) + ^ V relaxation
Cvv(i,j,l)
A rz—Stokes + ^Ri,j,l vStokes
vz(i,j,l) + ^ Vrelaxation
Cvz (i,j,l)
(21) (22)
(23)
(24)
where CmX ^SaitAon are relaxation parameters. Cvr (i,j,l) , Cvv(i,j,l) A* (i,j,l) is the coefficients at vr(i,j,l)vv(i,j,l), vz(itjtl) in Stokes equations (1). For models with variable viscosity, the respective residuals in Eqs.(22)-(24) become:
AR
r—Stokes i,j,l
AR
sp—Stokes i,j,l
r Stokes
Rp—Stokes Ri,j,l
dr
Trr TC
pp
1 ÖT,
rp
r
ÖT,
rp
2t,
rp
1 0Tc
pp
a rz—Stokes i,j,l
R
z— Stokes i,j,
dr ÖTr.
dr
r
Trz
r
r
1 dx
pz
r
dTrz + dP öz ör '
dTpz + 1 dP öz rd^'
dTzz + dP
öz öz
(25)
(26) (27)
where t is a deviatoric stress tensor.
B.1 Discretization of the continuity equation
Fig. 2 shows an elementary volume (cell) of a 3D staggered grid that can be used for discretization. Using the stencil with six velocity nodes around a cell (Fig.2), equation (21) takes the form:
AR
continuity i,j,l
R
continuity i,j,l
r j+Ar/2
1_vV(i + 1,j + 1,l + 1)—vV(i,j + 1,l + 1)
rj +Ar/2
Ap
'r(i+1,j + 1,l + 1)+vr(i+1,j,l+1) 2
z(i+1,j+1,l+j)—vz(i + 1,j + 1,l)
Az
r(i+1,j+1,l+1)—vr(i+1,j,l + 1)
Ar
r
new
■
i
■
Fig. 2. Indexing of different variables for a 3D staggered grid
Fig. 3. Stencil of a 3D staggered grid used for the discretization of the r-Stokes equation with variable viscosity
B.2 Discretization of Stokes equations
An example of stencil for the discretization of the r-Stokes equation is shown in Fig.3. For terms in equations (22), (25), we obtain the following discretization:
1. Pressure derivative:
dP = P(i-i,j,i-i) — P(i-i,j-i,i-i)
dr Ar
Discretization for terms which contain deviatoric stress tensor components:
2.
/drrrs , o 1 f vr(i,j+1,Z) — vr(i,j,Z) vr(i,j,z) — vr(ij-1,0s
)i>j>1 = (dr(2n ))i'j'1 = 2Ar ^(i-i'^-i)-Ar--n(i-i,j-i,i-i)-Ar-)
3.
(t ) = (2ndvr) = (n + n ) vr(i'j+i'1) — vr(ij-U)
(Trr )i'j'i = (2n_dr )i'j,i = (n«(i-ij'i-i) + nn(i-i'j-i'i-i))-2A7-
4.
(rw) j = (2n( i ^ + ^ )) j = j (nn(i-i>j>i-i) + n^-ij-u-i) )( ^ ( —
— v^.o+W-1.^ )+ vr(i;j.0)
5.
( dTry ) =( (n( i dvr I dv^ _ ))) = (n , ,( i vr(i+1.j.1)-vr(i.j.1) + vy(i.j+1.Q-'Mi.j.O
( )i.j>1 ((n(r + dr r )))i.j.1 A^ (nr^(i,i,1-i)(r^ A^ + Ar
vy(i.j + 1.Q +vy(i.j.Q ) _ n ( vr(i.j.Q-vr(i-1.j.O + vy(i-1.j + 1.i)-vy(i-1.j.Q__i vy(i-1.j + 1.i)+vy(i-1.j.i) ))
rj 2 ) nr^(i-i,i,1-i)( rj. A^ + Ar rj 2 ))
6.
( dTrz ) =( (n( dVz I dVr ))) = A. n , J vz(i.j+1.1)-vz(i.j.1) + vr(i.j.1+1)-vr(i.j.Q )_
( dz ( dz(n( dr + dz )))«,j,1 Az //rz(t-i,j,i)( Ar + Az )
— nrz(i-i,j,i-i)(-Ar--1 Az-)
Discretization for coefficient CVr(i)j)1) in equation (16) takes the form:
C _ _2 nn(i-1.j.i-1)+nn(i-1.j-1.i-1) nn(i-1.j.i-1)+nn(i-1.
CVr(i,j,i) = 2 Ar2 r2
J
i nry(i.j.i-1)+nry(i-1.j.i-1) nrz(i-1.j.i)+nrz(i-1.j.i-1)
2 A^2 Az2
Discretization for the ^-Stokes and z-Stokes equations can be constructed and indexed analogously.
Acknowledgements
This work was partially financially supported by the Government of the Russian Federation (grant 074-U01), by Ministry of Science and Education of the Russian Federation (GOSZADANIE 2014/190, Projects No 14.Z50.31.0031 and No. 1.754.2014/K)), by grant MK-5001.2015.1 of the President of the Russian Federation. The authors thank Prof. P. Tackley and Prof. T. Gerya for interesting discussions.
References
[1] Zhang H., Xu Sh., Jeffries G.D.M., Orwar O., Jesorka A. Artificial nanotube connections and transport of molecular cargo between mammalian cells. Nano Comm. Net., 2013, 4(4), P. 197-204.
[2] Li T.-D., Gao J., Szoszkeiwicz R., Landman U., Riedo E. Structured and viscous water in subnanometr gaps. Phys. Rev. B, 2007, 75, P. 115415.
[3] Popov I.Yu., Chivilikhin S.A., Rodygina O.A., Gusarov V.V. Crystallite model for flow in nanotube caused by wall soliton. Nanosystems: Phys. Chem. Math., 2014, 5(3), P. 400-404.
[4] Popov I. Statistical derivation of modified hydrodynamic equations for nanotube flows. Phys. Scripta, 2011, 83, P. 045601/1-5.
[5] Mashl R.J., Joseph S., Aluru N.R., Jacobsson E. Anomalously immobilized water: A new water phase induced by confinement in nanotubes. Nano Letters, 2003, 3(5), P. 589-592.
[6] Kolesnikova S.A., Loonga C.K., de Souzaa N.R., Burnhamb C.J., Moravskyc A.P. Anomalously soft dynamics of water in carbon nanotubes, Physica B, 2006, 385-386, P. 272-274.
[7] Popov I., Makeev I. A benchmark solution for 2-D Stokes flow over cavity, Z. Angew. Math.Phys., 2014, 65, P. 339-348.
[8] Popov A.I., Gerya T., Lobanov I.S., Popov I.Yu. Benchmark solutions for nanoflows. Nanosystems: Phys. Chem. Math., 2014, 5(2), P. 391-399.
[9] Popov A.I., Gerya T., Lobanov I.S., Popov I.Yu. On the Stokes flow computation algorithm based on Woodbury formula. Nanosystems: Phys. Chem. Math. 2015, 6(1), P. 140-145.
[10] Gerya T. Introduction to Numerical Geodynamic Modelling. Cambridge: Cambridge University Press, 2010.