Benchmark solutions for nanoflows Текст научной статьи по специальности «Физика»

BENCHMARK SOLUTIONS FOR NANOFLOWS

A.I. Popov1, I. S. Lobanov1, I.Yu. Popov1, T. V. Gerya2 1 ITMO University, Kronverkskiy 49, 197101, St. Petersburg, Russia 2 Institute of Geophysics, Department of Earth Sciences, Swiss Federal Institute of Technology Zurich (ETH), 5 Sonnegstrasse, CH-8092 Zurich, Switzerland

popov239@gmail.com, taras.gerya@erdw.ethz.ch PACS 47.10.ad, 47.11.Fg

Essential viscosity variation creates additional difficulties for numerical investigation of flows through nan-

otubes and nanochannels. Benchmark solutions of the Stokes and continuity equations with variable viscosity

are suggested. This is useful for testing of numerical algorithms applied to this problem.

Keywords: nanotube, Stokes flow, benchmark solution.

Received: 03 May 2014

Revised: 11 May 2014

1. Introduction

Flows through nanotubes and other nanostructures have many interesting peculiarities. One of them is the viscosity variation (see, e.g., [1], [2], [3]). Flows in nano-channels are influenced by local heterogeneity of molecular structure of the liquid if its size is compared with the channel width. A hypothesis about the existence of locally-ordered structures in liquid was put forward in [4]. Investigations of fluid flows in nano-sized domains show that it is strongly influenced by local ordering of nano-sized scale. Experiments [5], [6] show that the effective viscosity of water in nanochannel with hydrophilic walls is essentially greater than the corresponding macroscopic value. Experimental and theoretical investigations of water state in carbon nanotube [7], [8] show that there is an ice-like envelope with liquid water inside in the nanotube. Increasing of effective fluid viscosity via channel diameter was marked in [9] for channels of a few micrometers diameters. Thus, experiments confirm high viscosity variations for the flow in a nanotube, which creates computational problems. Namely, the convergence of the numerical algorithms in the case of strongly varying viscosity is not good, and, moreover, is not guaranteed ( [10], [11]). Correspondingly, one need an instrument to choose an appropriate numerical scheme. One can make a choice by using of benchmark solutions (see, e.g., [12], [13], [14])

In present work, we suggest methods for algorithm checking. The scheme of the algorithm testing is as follows. Consider a rectangular domain. Calculate the values of the benchmark solutions at the rectangle's boundary. Take these values as the boundary conditions. Due to the uniqueness theorem the solution of the boundary problem should coincide with our benchmark solution. So, we obtain a solution of the specific boundary problem. Note that we derived the solution analytically. Next, we solve the same boundary problem by a numerical algorithm, then we compare results and estimate the quality of the numerical algorithm.

We have found exact analytical solutions of the Stokes and continuity equations in the two-dimensional case for linearly varying viscosity. These solutions are convenient to use as benchmarks for numerical algorithm testing. The efficiency of the approach was

demonstrated on a numerical algorithm for calculations of the Stokes flow with varying viscosity.

2. Formulation of Stokes and continuity equations with variable viscosity

Consider the plane flow. 2D Stokes equations for the case of varying viscosity has the form:

2n

d 2vx dx2

+ 2

dn dvx

+ n

d2 vx

+n

d2u

y + dn dvx

+

d 2vy

' dx2

+n

dx dx dy2 dydx dy dy

dn dvy dP g

dy dx dx x'

d2vx dn dvx dn dvy dn dvy

+ T-, 7-, + T-, 7-, + 2T-, 7-, +

(1)

dydx dx dy dx dx dy dy

2r]d2vy _ dp = _pG

dy2 dy y'

(2)

+

dvy,

0.

(3)

dv,

dx ' dy

Here (vx,vy) is the flow velocity, n = n(x,y) is the viscosity, P is the pressure, p is the density, (Gx, Gy) is the gravitational force. Note that (3) is the continuity equation. Let us change the variables vx,vy,P in such a way that:

dvx 1 dux dvx 1 dux

n dx 1 duv

dx

dvy =__„

dx n dx ' 1 dP _ dP n dx dx}

dy

dvy =

dy 1 dP n dy

n dy 1 duy n dy _ dP dy .

(4)

(5)

(6)

The correctness conditions for such replacement are as follows:

d , 1 dux -( x

d , 1 dux -( x

d , 1 duy.

d , 1 duy.

dy n dx dx n dy dy n dx dx n dy

dPi\- dLi dL\

dy dx dx dy '

These conditions lead to the following correlations:

drq dux dn dux drq duy drq duy

dy dx

dx dy '

dn dP

dy dx

dy dx

dn dP

dx dy

dx dy

All conditions give one the same characteristic equation:

dn , dn ,

—dx + — dy dx dy

0.

Evidently, n(x,v) = C is an integral of the equation. Hence, the solutions of our equations, which predetermine the correctness of replacement suggested above, are

$(n), uy = ^(n), P = P(n).

ux

After replacement, the Stokes equations (1), (2) and the continuity condition (3) transform to the following form:

202ux + 02ux + d2Uy _ ^dp _ _ dx2 dy2 dy dx dx x'

d2uy + + 2d!uy _ ndi_ _pG , (8)

dx2 dydx dy2 dy y'

^ + ^ _ 0. (9)

dx Sy

Inserting the expressions for ux,uy into (7), (8), (9), one obtains the following equations:

2$'^ + 2$" f ^ 2 + $/ ^ + $// (^ 2 + d2x d2y \dy/

*/ 92n + *„ dv^V _ npi _ _pG (10)

y x y x - x - x

2*^ + 2* (^V + */ d!l + * ( ^V +

2y y 2x x

^ + $" drP- _ VP' IT _ _pGy, (11)

y x y x y

$/dr + *dr _0. (12)

x y

3. Exponentially varying viscosity

Let us construct the second benchmark solution. Next, we assume that the viscosity is the exponential function of the Cartesian coordinates:

n _ cexp (ax + by). (13)

General consideration up to (10), (11), (12) is the same as earlier. By inserting (13) into (10), (11), (12) and taking into account that:

dn dn

— _ an, TT _ bn, x y

one obtains the following system of equations:

(2a2 + b2)($"n2 + $/n) + ab(*"n2 + *n) _ aP'q2 _ _pGX) ab($"n2 + $/n) + (a2 + 2b2)(*"n2 + */n) _ bP'q2 _ _pGy,

a$/ + b* _ 0.

Using the last relation, we exclude * from the first two equations:

(a2 + b2)($"n2 + $/n) _ ap/n2 _ _pGx:,

a3 + ab2 ($"n2 + $/n) _ bP/n2 _ _pGy,

a

*/ _ _ a $/. (14)

b

One can see that we obtain a linear algebraic system with respect to ($"n2 + $'n) and P1. The solution is as follows:

P = ^, (15)

n2

$"n2 + $'n = bf (n). (16)

Remark. It is interesting that these formulas contain the same functions f (n), fi(n).

Equation (16) is a well-known Euler ordinary differential equation. One can get its solution for arbitrary function f:

r n f (ni)

ux = $(n) = b log( — )-dni + bci logn + c2.

i ni ni

Taking into account relation (14), one obtains u

y

uy

r n f(ni)

^(n) = -a log(—)-)dni - aci log n + C3.

i ni ni

(17)

(18)

Taking into account (4), (5), one obtains vx, vy:

r dni f1 f M 1 + +

b / _T dn2--bci- + bci + C2

'i m Ji n2 n

/V

dn2

f (n2) r dni 1+ +

- —5— bci- + bci + c2.

n2 JV2 n2 n

Vx

Hence, we get the expression for vx and analogously, for vy:

, f \ f (n2) n - n2 , 1 . , .

M dn2---bci- + bci + C2,

Ji n2 nn2 n f \ f (n2) n - n2 , 1 .

dn2---+ aci- + C3.

1 n2 nn2 n

Vy = - a

(19)

(20)

As for the pressure, we obtain it from (15) by taking into account (6):

p ind fi(ni) +

P = dni—2--+ C4.

1 n2

Hence,

P

f1(n1) dni—7--+ c4.

(21)

1n1

For a simple particular case (constant gravitational term), when f (n) = A = const, fi(n) Ai = const one has:

b(A + ci) bA log n

+ C2,

Vx =

n

n

a(A + ci) aA log n

Vy =--1---+ C3,

nn P = Ai log n + c4 — bi, A more complicated case is when the density is a linear function of the viscosity, p = Pin + fc, i.e.

f (n) = ain + a2, fi(n) = bin + b2,

where constants a1, a2, b1, b2 are the same as in the previous section. It is simple to evaluate integrals in (19), (20), (21). In such a way, one obtains:

b(ai — a2 — ci) , log n . ~

Vx = bai log n +---ba2--+ C2,

n

n

V

x

a(ai _ a2 _ ci) log n , ~

vy _ _aa1 log n---+ aa2--+ c3,

nn P _ bin + b2 logn + P4, where c2 _ c2 + bc1 + ba2 _ ba1, c3 _ c3 + aa1 _ aa2 _ ac1, c4 _ c4 _ b1.

4. Example problems and numerical convergence tests

The scheme of algorithm testing is as follows: initially, we have obtained particular solutions of the Stokes and continuity equations for the exponential type of viscosity variation. Let us choose a domain, e.g., a rectangle in 2D case. We calculate values for velocity and pressure given by our analytical solution and take these values as the boundary conditions. Then, due to the uniqueness theorem, the solution of the boundary problem in the domain should coincide with our analytical solution. Let us compute the solution of the boundary problem by a numerical method. Comparison of the result with the exact analytical solution shows the quality of the numerical algorithm.

4.1. Exponentially varying viscosity

Consider a simple example of such flow in a rectangle 0 ^ x ^ xsize, 0 ^ y ^ ysize. We assume that n _ ax + by + c. We will mark the exact solution obtained in Section 2 as vx,a, vy,a, Pa. It is the solution of the boundary problem in the rectangle Q with the following conditions at the boundary dQ _ {x _ 0, x _ xsize,y _ 0, y _ ysize} :

vy |dn _ vy,a, vx |dQ _ vx,a.

Let us compute the velocity and pressure using the finite-difference scheme. The corresponding solution is marked as vx n, vy n,

Pn. The deviation of these values from the exact solution (vx,n _ vx,a,vy,n _ vya,Pn _ Pa) is related with the error of the numerical scheme. We calculate the relative errors of three types: L^,L1,L2 for different viscosity contrasts , i.e. different values of the coefficients a, b. We test the program Stokes2D-variable-viscosity1 from [10]. The results are presented in Fig. 1-6. Namely, figures 1-3 correspond to low viscosity contrast, Figures 4-6 — to high viscosity contrast. Particularly, Fig. 1 and Fig. 4 show pressure and velocity components distributions. Fig. 2 and Fig. 5 characterize the viscosity and the density distributions. Fig. 3 and Fig. 6 contain plots of relative errors via the grid resolutions in logarithmic scale. The viscosity contrast, i.e. the values of the coefficients in the expression for the viscosity, is determined by the given values of the viscosity at three rectangle corners. The value of the viscosity at the initial rectangle corner is 1, n2,n3 are the values of the viscosity at two adjacent corners. For all figures, "n" means "numerical solution", "a" means "analytical solution" (benchmark).

For the case of exponentially varying viscosity, we made calculations for the following system parameters:

C _ n1,a _ (log(n3) _ log(n1))/x«ze,

b _ (log(n2) _ log(nl))/Уsize, n _ C exp (ax + by),

p _ n + ^2

xsize _ ysize _ 1 ,

Gx _ 10, Gy _ 10, n1 _ 1,^1 _ _ 3 x 103,

Fig. 1. Distribution of vx, vy and P; 2D case, exponentially varying viscosity, low viscosity contrast (n2 = n3 = 5).

Fig. 2. Distribution of viscosity n and density p; 2D case, exponentially varying viscosity, low viscosity contrast (n2 = n3 = 5).

One can see that there is rather high accuracy for the numerical approach. Figures 1-3 corresponds to the case of low viscosity contrast, figures 4-6 — to the case of high viscosity contrast. We observe the conventional situation — L^-error is the largest among the considered errors norms, and Li-error and L2-error are similar. The calculations show that one has good convergence of the numerical scheme for small viscosity contrast, but it is not so for high viscosity contrast (compare Fig. 3 and Fig. 6).

5. Conclusion

Numerical analysis of geophysical flows presents many difficulties. It is related with complex dependence of material parameters on spatial coordinates. Different schemes of numerical calculations are suggested. To establish the quality of suggested approach it is possible to compare the results of different numerical methods. More reliable examination of the approach is given by the comparison with the exact solution of the problem, similar to the considered one. For this purpose, one needs such a benchmark solution. In the present paper, we suggest a benchmark solution for the Stokes equation coupled with the continuity equation where the viscosity is exponentially dependent upon the spatial Cartesian coordinates. Comparison of the numerical result with this exact solution allows us to determine the order of convergence, the quality of discretization, etc.

Fig. 3. Logarithm of the relative error via logarithm of the grid step; 2D case, exponentially varying viscosity, low viscosity contrast (n2 = n3 = 5); blue line - pressure, green - vx, black - vy; line - Li-error, dashed line - Lœ-error, line with dots - L2-error.

Fig. 4. Distribution of vx, vy and P; 2D case, exponentially varying viscosity, high viscosity contrast (n2 _ n3 _ 100).

Fig. 5. Distribution of viscosity n and density p; 2D case, exponentially varying viscosity, high viscosity contrast (n2 _ n3 _ 100).

error norms convergence

■1 r

■1.5

■5-1-1-'-1-'-L

-1.7 -1 B -1.5 -1.4 -1.3 -1.2 -1.1

Fig. 6. Logarithm of the relative error via logarithm of the grid step; 2D case, exponentially varying viscosity, high viscosity contrast (n2 = = 100); blue line - pressure, green - vx, black - vy; line - Li-error, dashed line - L^-error, line with dots - L2-error.

Acknowledgement

This work was partially financially supported by the Government of the Russian Federation (grant 074-U01), by the Ministry of Science and Education of the Russian Federation (GOSZADANIE 2014/190, Project 14.Z50.31.0031), by grants of the President of Russia (state contracts 14.124.13.2045-MK and 14.124.13.1493-MK).

References

[1] Li, D. (Ed.) Encyclopedia of Microfluidics and Nanofluidics. Springer, New York (2008).

[2] Popov I.Yu. Statistical derivation of modified hydrodynamic equations for nanotube flows, Phys. Scripta., 83, P. 045601/1-3 (2011).

[3] I.Yu. Popov, S.A. Chivilikhin, V.V. Gusarov. Flows in nanostructures: hybrid classical-quantum models. Nanosystems: Phys., Chem., Math., 3(1), P. 7-26 (2012).

[4] Frenkel Ya. I. New development of the theory of liquid state. Appl. Uspekhi Fiz. Nauk., 25. P. 1-18 (1941).

[5] Antognozzi M., Humphris A. D. L., Miles M. J. Observation of molecular layering in a confined water film and study of the layers viscoelastic properties. Appl. Phys. Lett., 78, P. 300-302 (2001).

[6] Li T.-D., Gao J., Szoszkeiwicz R., Landmann U., Riedo E. Structured and viscous water in subnanometer gaps. Phys. Rev. B, 75, P. 115415-1115415-6 (2007).

[7] Kolesnikov A.I., Zanotti J.-M., Loong Ch.-K., and Thiyagarajan P.Anomalously Soft Dynamics of Water in a Nanotube: A Revelation of Nanoscale Confinement. Phys. Lett., 93, P. 035503-1035503-4 (2004).

[8] Kolesnikova A.I., C.-K. Loonga C.K., de Souzaa N.R., Burnhamb C.J. Moravskyc A.P. Anomalously soft dynamics of water in carbon nanotubes. Physica B, 385386, P. 272-274 (2006).

[9] Deryagin B. V., Ovcharenko F. D., Churaev V. N. (Eds.) Water in Disperse Systems. Moskow, Nauka (1989).

[10] Gerya T. Introduction to Numerical Geodynamic Modelling, Cambridge University Press, Cambridge (2010).

[11] Gerya T.V., May D.A., Duretz T. An adaptive staggered grid finite difference method for modeling geodynamic Stokes flows with strongly variable viscosity. Geochem., Geophys., Geosyst. 14, DOI: 10.1002/ggge.20078 (2013).

[12] Blankenbach B., Busse F., Christensen U. et al, A benchmark comparison for mantle convection codes. Geophysical Journal International 98, P. 23-38 (1989).

[13] Gerya T.V., Yuen D.A., Characteristic-based marker-in-cell method with conservative finite-difference schemes for modeling geological flows with strongly variable transport properties. Phys. Earth Planetary Interiors, 140, P. 293-318 (2003).

[14] Gerya T.V., Yuen D.A. Robust characteristics method for modeling multiphase visco-elasto-plastic thermo-mechanical problems. Phys. Earth Planetary Interiors, 163, P. 83-105 (2007).