A computer-algebraic investigation of an oscillating boundary layer Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 20, № 1, 2015

A computer-algebraic investigation of an oscillating boundary layer

Karl G. Roesner Darmstadt, Germany

http://www.kgroesner.de, e-mail: kgroesner@t-online.de

For large Reynolds numbers analytical expressions are derived for the velocity components and the pressure in the boundary layer between two oscillating parallel planes. The linearization of the boundary layer equations leads to a coupled system of three parabolic differential equations with time-dependent or homogeneous boundary conditions. The solutions are derived using computer algebra software on a personal computer. The vortical character of the flow pattern — shown in experiments — resembles the famous Gortler vortices in the boundary layer along a curved wall. A criterion for vortex identification for this type of motion is at least fulfilled for half the period of the oscillating planes.

Keywords: boundary layer, time-periodic flow, pattern formation, Gortler vortices, parabolic differential equations, vortex identification, computer algebra.

Introduction

The following investigation is based on a computer-algebraic analysis of the time-dependent boundary layer equations, describing the oscillating flow in a fluid layer between two parallel plates induced by the torsional motion of the boundaries. The analytical solution of the linearized boundary layer equations results from the solution of three parabolic differential equations describing the initial-boundary value problems for the velocity components. The solution which defines the velocity components can be split into three terms, describing the transient, oscillating and stationary flow in radial and axial direction. The stationary flow is responsible for the mass transport in the boundary layer. Experiments show a formation of streamwise vortices in the boundary layer. Once the cylinder is set into oscillation, in the vicinity of the axis of rotation — close to top and bottom boundary — vortices are created, which propagate radially outward. This phenomenon resembles the well known Gortler vortices along concave walls. In the present case, the curved streamlines in the stagnation region play the role of a concave rigid wall, which causes Gortler vortices along the curved boundary. In the present case the vortex character of the flow can be analysed using the symmetric and skew symmetric part of the vector gradient of the Navier — Stokes equations. It can be shown that a local pressure minimum exists due to vortical motion [1]. In Fig. 1 the vortical structure in the vicinity of the axis of an oscillating flat cylinder results from the self organisation of the fluid under the torsional motion around the axis. The initially separated radial vortices interact while being pulled outward radially under the action of the steady radial flow enforced by the centrifugal force. This leads at last to a pattern of entangling interacting vortices.

© ict sb ras, 2015

Fig. 1. Streamwise vortices in an oscillating boundary layer

1. Basic equations and boundary conditions

In a cylindrical coordinate system (r,p,z), the Navier — Stokes equations in the following form describe the cylindrical flow of a fluid with circumferential symmetry character. This means that the functions u,v,w,p and all their partial derivatives with respect to p are not depending on the independent variable p.

d2v 1 dv

du du du v2 1 dp

— + u— + w----= —— + v

ot or dz r Q or

dv dv dv uv

— + u— + w— +--= v

at dr dz r

dw dw dw 1 dp

— + u— + w— = —— + v dt dr dz q dz

d 2u 1 du u d 2u

--1------1--

dr2 r dr r2 dz2

+ ---o +

d 2v

dr2 r dr r2 dz2

d2w 1 dw d 2w dr2 r dr dz2

(1) (2) (3)

u,v,w are the velocity components in r-, p- and z-direction, p is the pressure and v is the kinematic viscosity of the fluid. In addition, the continuity equation reads like

1 d (ru) dw

---—- +--= 0.

r dr dz

(4)

We assume that the motion of the two parallel planes, oscillating with the frequency n around the common axis, is limited to small angular amplitudes u ^ 1, that the following estimation holds:

v

Q(t) = u sin nt ^ 1.

The circumferential velocity component of boundary points in the planes z = 0 and z = d, is therefore prescribed by the function

v(r, 0, t) = v(r, d, t) = r cos nt, := un,

where d is the vertical distance of the parallel plates. This boundary condition prescribes a harmonic oscillation of the two planes parallel to each other with the frequency n and the velocity amplitude rQ0 around the z-axis.

For the radial and axial velocity components u(r,z,t) and w(z,t) we get homogeneous boundary conditions at z = 0 and z = d :

u(r, 0, t) = u(r, d, t) = 0, w(0, t) = w(d, t) = 0.

2. Nondimensional equations

The variables u(r,z,t),v(r,z,t),w(z,t) andp(z, r, t) are transformed into dimensionless form introducing the distance d of the two planes as length scale for the coordinates z and r and the frequency n for the time t.

y :=

-,, r

r -

nt.

The nondimensional velocity components and the pressure are defined as follows:

u :=

u(r, z, t) Q0d

v :=

v(r, z, t) Q0d

Introducing the parameter

u

w := n

w(z, t) Q0d

p :=

p(r, z, t) qQ0 d2 '

which represents the ratio of the amplitude of the angular velocity and the frequency of oscillation of the planes, into the dimensionless system of differential equations leads to the following system of equations:

1 du _3u dU v2

--T + u~ + w~ — —

u dt dr dy r

1 dV _3V _ dV uv

--r + u— + w— + —

u dt dr dy r

1 dw _ dw _ dw

---T- + u— + w —

u dt dr dy

1 d(ru) + dw

r dr- dy

dp

dr d2Qc\

d2 u 1 du

u

dr2

r dr

d2u dy2

d2tto

d 2v

v

- dp + dy d2Q0

1 dv

• dr d 2w dr2

v

. . . d2v dr2 r dr r2 dy2_

1 dw d2w r dr dy2

= 0.

(5)

(6)

(7)

(8)

The dimensionless boundary conditions are

u(r, o,t) = u(r, 1,t) = 0, v(r, 0,t) = v(r, 1,t) = r cos t, w(o,t) = w(1, t) = 0.

v

v

3. The ansatz functions

Introducing the ansatz functions [2] , [3], F(y,t), g(y), and P(y,t), the velocity components U, V, w and p are defined in the following way:

u(r,y, t)

v(r, y, t) w(y,t)

p(r,y,t)

This choice of ansatz functions is similar to the analytical approach of the flow near a rotating disk which was analysed the first time by Th. v. Karman.

= ru Fy(y,t),

= r9(y,~t),

= —2u F(y,t),

- 1 2 -= P (y,t) + - r2K (t)

The functions u,v depend linearly on r. w is independent of r according to the fact, that the boundary layer thickness 8 of the flow around the stagnation point at the plate is constant, not depending on the radius r. p contains the socalled centrifugal term, modified by the factor K(t), which takes into account the time dependent influence of the centrifugal force. Introducing the ansatz functions into the system (5)-(8) — omitting the bar on the dimensionless quantities and denoting the dimensionless time by t, we arrive at the following system of partial differential equations for the fluid under consideration with the Reynolds number

nd2

R :=

> 0,

(9)

Fyr + K - 2FFyy] - g2 9t + 2^2 [Fyg - Fgy]

Ft — luJ2 F Fy

-K (T) + RFyyy,

gyy,

R

1 P+i F 21 y + rf yy'

The boundary conditions for the ansatz functions are

F(0, t) = F(1,t) = 0, Fy (0, t ) = Fy (1,t ) = 0 g(0,T ) = g(1,T )

cos T.

(10) (11) (12)

The ansatz functions fulfill the continuity equation identically. In addition to these boundary conditions the functions have to show the following symmetry and antisymmetry character, respectively:

Fy (y,T )

F (y, t ) g(y,T )

P (y,T )

Fy(1 - y,T),

-F(1 - y,T),

g(1 - y,T ), P(1 - y,T).

4. The boundary layer approximation

Based on the assumption that is of order 0.1 s-1 = 5.7° s-1, we assume

u2 < 1, (13)

which can be realised experimentally for a moderate n of about 50 s-1. As a consequence of this assumption the differential equations are linearized and define the boundary layer equations for the flow in the vicinity of a fast oscillating plane:

Fy

yT

gT Ft

RFyyy + g2 (y,T) - K (t)

1

R

gyy,

— F + - P rf yy + 2P y.

(14)

(15)

v

1

These equations can be solved analytically. In (14) we introduce a new function

$(y,T ):= Fy (y,T) (17)

and get an inhomogeneous parabolic differential equation with homogeneous boundary conditions of the form

= R$yy + g2 (y,T) - K(t). (18)

The solution of this equation determines the radial component

u(r, y, T) = ru $(y, T). Integrating equation (14) with respect to y, we get

Ft = RFyy +J g2 (y,T) dy - yK (t) + G)(t), (19)

the solution of this equation determines the axial velocity component

w(y,T) = -2uF (y,T).

If we compare (16) with (19), the pressure function P(y,T) can be found by a twofold quadrature:

y / n

P(y,T) = 2j U g2 (Z,T) dZ I dn - y2K(t)+ yC0(T) + Ci(t). (20)

As only the pressure gradient is physically relevant, we can assume C1(t) = 0.

The three parabolic differential equations for the velocity components and the relation for the pressure function are solved successively in the following order:

_ 1

gT = Rgyy,

y / n

P(y,T) = 2 g2 (Z,T) dZ I dn - y2K (t)+ yC0(T),

1

R

Ft = RRFyy + y g2 (y,T) dy - yK (t) .

5. The circumferential velocity component

A solution of (15) which obeys the time dependent boundary conditions is defined by

g(y, T) := G1 (y) cos T - g2 (y) sin т, (21)

G1 (y) : = eXy (711 cos Xy - Y12 sin Ay) + e-Xy (721 cos Xy + Y22 sin Xy), G2 (y) : = eXy (712 cos Xy + 711 sin Xy) + e-Xy (722 cos Xy - Y21 sin Xy).

= r$yy + g2 (y,T) - K (t)

i l i l i l i l 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

-1.0

Fig. 2. ViM,!) >T = kl,k = 0,..., 8 rw 8

The coefficients Yik, i,k = 1, 2 are summarized in Appendix I. As initial condition we define g(y, 0) := Gi(y).

Writing (21) in a slightly different form we get

g(У, T) = eyX (Yii cos (y\ — t) - Yi2 (sin (yA + T))) + +e-yX (Y21 cos (y\ + t) + 722 sin (y\ — t)).

This can be interpreted as a superposition of two damped waves running from each of the boundaries into the interior of the fluid layer with the phase velocity

_ 1 _ l~2

vphase ± T" ±\ ~n .

A V R

The circumferential velocity distribution in the fluid layer is therefore

v(r, y,T) = r u (G1 (y) cos t — G2 (y) sin t).

(22)

Figure 2 shows a sequence of velocity distributions in the fluid layer for nine time steps. For a chosen Reynolds number of R = 500, the oscillating walls influence a fluid layer of approximately 20 % of the total thickness of the fluid.

Taking into account the initial condition for g(y, t), we have to add the transient part to the solution (21). The transient solution of the homogeneous differential equation is given by

„__ 2 2

tr 4 — n n T ■

gh = Ane R sin nny.

n=l

The coefficients An are listed in the Appendix II. After a transient time of t = 20 the maximal amplitude of ghr is approximately of order 10-4. Therefore, for further calculations only the time periodic solution (22) will be used.

6. The pressure function

Performing the twofold quadrature in (20) results in the following solution

P(y, t) = 2 (n (y) cos2 t + n2 (y) sin2 t - n (y) sin 2t) - y2K (t) + yCo(T).

The functions n, n2, n3 are listed explicitly in the Appendix III. From the symmetry character of P(y,T) one derives the two expressions for K (t) and C0 (t)

K (t) = (711721 - 712722) cos 2t, (23)

Co (t) = 2 (711722 + 712721) sin 2t. (24)

Finally, the pressure function P(y, t) is given by

P(y, T) = 2 (n (y) cos2 T + n (y) sin2 T - n (y) sin 2t) --y2 (711721 - 712722) cos 2t + 2y (711722 + 712721) sin 2t. (25)

7. The radial velocity component

Equation (18) is an inhomogeneous parabolic differential equation with homogeneous boundary conditions and the initial condition $(y, 0) = 0. The solution of this initial-boundary value problem is represented by an expansion with respect to the eigenfunctions sin (nny) , n E N, of the homogeneous equation. In a first step the inhomogeneity has to be developed into a series of the eigenfunctions which leads to coefficients U(k),V(k),W(k), which are listed in the Appendix IV.

g2 (y, t) - K (t) =

= 2 j (U(k) cos2 t + V(k) sin2 t - W(k) sin 2t) sin (2k + 1) ny +

k=0

+2 (711721 - 712722) j 4RR+ k2 e-RT sin (n (2k + 1) y) -k=0

2 (711721 - 712722) j R 2R ^ + Kr 2T sin (n (2k + 1) y). (26)

k=0

As mentioned at the beginning, the solution of equation (18) can be split into three separate parts which will be denoted by the abbreviations: tr for transient, os for oscillation, and st for stationary flow. The indices i and h stand for inhomogeneous or homogeneous, respectively.

The transient solution

~u(y,T )tr = ru

(X

2 j

k=0

( -RU(k) (2R2 + k2) - 2V(k)R3 + W(k)Rk2 \

k (k2 + 4R2) Rk (1 - cos n (2k + 1))

(4R2 + k2) (2k + 1) n

Y : = 711721 - 712722, k : = n2 (2k + 1)2 .

+

k

— TZT

e R sin (2k + 1) ny, (27)

0.15

o.io

0.05

o.oo

0.0 0.2 0.4 0.6 0.8 1.0 y

uf (r, y, T)

Fig. 3.

rw

radial component

The oscillating solution

—= rw

x

T

k=0

R

U(k)K cos 2t + 2R sin2r-2 -V(k)(K cos 2t + 2R sin2r - k) -

K +4R \ -2W(k)K cos 2t + 2Rsin2r

, (2Rsin 2t + k cos 2t) (1 - cosn (2k + 1))'

V

-27R-

(4R2 + k2) (2k + 1) n

sin (2k + 1) ny. (28)

The stationary solution

(X

T.

k=0

R (U (k) + V (k)) - 2

k

sin(n (2k + 1) y).

(29)

As (29) is a fast convergent series, it is sufficient to take into account only 20 terms. The approximate solution

20

v-^R(U (k) + V (k)) - 2 . . , J] ( ( ) „ ( ))-sin(n (2k + 1) y)

k=0

k

is shown in Fig. 3, where R = 500.

This velocity profile indicates the boundary layer character of the flow by a strong velocity gradient close to the walls.

The stationary flow in radial direction is responsible for the formation of streamwise vortices within the boundary layer. The transient part of the solution decays very fast, and the oscillating behavior is not responsible for the radial material transport, which is observed in the experiment. Therefore only the stationary term plays the major role in the analysis of the vortex formation.

8. The axial velocity component

The component w(y,T) can be calculated by solving the partial differential equation (19)

Ft = + g2 (y, T) dy - yK (t) + C0(t)

0.006

0.004 "

^A/0.5 0.6 0.7 0.8 0.9 1.0 y

-0.006

Fig. 4. -

wf (y, t)

w

axial component

with given functions K (t) and C0(t), defined by formulas (23), (24). The inhomogeneity

g2 (y,T) dy - yK (t) + Co(t)

must be developed into a series with respect to the eigenfunctions sin (nny) and contributes to the solution as a driving term. The complete solution is a superposition of the three parts: The transient solution

—1 wf (y,T ) =

u

2 E

k=1

( BR (2R2 + n4n4) + 2B2R3 + 2Q3R2n2n2

n2 (2k)2 (n4 (2k)4 + 4R2)

\

(7n2 (2k)2 + 8^R) R cos 2kn

V

4kn (4R2 + n4 (2k)

2 2 e-"r (2k) T sin 2kny.

/

The oscillating solution

(

1

--wi

u

\y,T) = 2 E

k=i

(Bi - 62)

r2 := (711722 - 712721) ,

R (n2 (2k)2 cos 2t + 2R sin2T) \

cos

r4 (OU\4

2 (n4 (2k)4 + 4R2)

+

V

^ 2R cos 2t + n2 (2k)2 sin2T ^

+B3R--—-+ 64

3 4R2 + n4 (2k)4 4

/

The stationary solution 1

—w: u

st

E R

k=1

,Bi + 62 + 2RK (t) n2(2k)2 n3 (2k)3

cos 2kn ) sin (2kny).

(30)

sin 2kny. (31)

(32)

The coefficients 01, 02,03 are summarized in the Appendix V. For the Reynolds number R = 500 the graph of the axial velocity component is given in Fig. 4.

4

9. Vortex identification

From several existing criteria, how to define a vortex in the flow field, the one which was stated by the authors Jeong and Hussain is the most convincing one and was used in the present investigation. By a decomposition of the vector gradient of the Navier — Stokes equations into a symmetric part S and a skew symmetric part A, the eigenvalues ^1,^2,^3 of the symmetric matrix

(

S2 + A2

, ,2 772 2

w Fy - g

0

0

2 2 2 w Fy - g

r (FyFyy^ + ggy) (gFyy -

r (Fy Fyy^2+ggy) ^

" 2 rw (gFyy - Fygy)

4w2F2

(33)

are calculated:

(y,T)

g2 5w2F2 1 -— +--y + -

- 2 2 2

A

g2 (g2 + r2g^) + g2 (r2F2y + 6F2) +

+r2F2g2 (yy

+-4F2 (r2F2y + 9F2)

+

^2 (y,T)

(y,T)

2

g2 + 5-2F2 - "2 +

2

A

g2 (g2 + r2g2) + g2 (r2F2y + 6 F2) +

+r2F2g2 (

+-4F2 (r2F2y + 9F2)

+

^2F2 - g2

(34)

According to w2 ^ 1 we can neglect all terms of order w2 and higher. Therefore we get

(y,T)

^2 (y,T) ^3 (y,T)

1

1

2 1

-ög2 + ög\/g2 + r2g2 > 0,

- 2 g2 - 2 g\/g2 + r2g2 < 0,

-g2 < 0, 0 < t < n/2.

(35)

(36)

j

(37)

(38)

(39)

Summing up the three eigenvalues we find

3

^ (y, t) < 0

i=1

which fulfills together with (37)-(39) the criterion for the existence of vortical structures in the flow field. The strong velocity gradient in the shear layer close to the oscillating plates is the region where vortices can be created and grow.

2

2

2

Acknowledgement

The manuscript was written using the software Scientific Workplace, Version 5.5, and most of the algebraic calculations were performed by the built-in MAPLE software.

References

[1] Jeong, J. and Hussain, F. On the identification of a vortex //J. Fluid Mech. 1995. Vol. 285. P. 69-94.

[2] v. Karman, Th. Laminare und turbulente Reibung // ZAMM. 1921. Vol. 1. P. 233-252.

[3] Rosenblat, S. Flow between torsionally oscillating disks //J. Fluid Mech. 1960. Vol. 8. P. 388-399.

Appendix I

Yii :

Yi2 :

Y21 :

Y22 :

mo -mo -

X •=* R

m0 (e-2A — cos 2X) + 2mi sinh X cos X

2 (cosh 2X — cos 2X) m0 cos X sin X — mi cosh X sin X cosh 2X — cos 2X '

m0 (e2A — cos 2X) — 2mi sinh X cos X 4 (sinh2 X + sin2 X) ,

mi cosh X sin X — m0 cos X sin X 2 (sinh2 X + sin2 X) ,

1 and mi = 1 • co-oscillation, 1 and mi = — 1 • counter-oscillation.

Appendix II

An

nn

(n2n2 — ex (R sin A + n2n2 cos A) cos nn) + \

Y11 R2 + n4n4 nn

+Y12^-4—j (R + ex (n2n2 sin A — R cos A) cos nn) +

R_n + n4n4 v '

—-—— (n2n2 + e-x (R sin A — n2n2 cos A) cos nn) +

R2 + n^n '

+Y22—ö-r~r [R — e-x (n2n2 sin A + R cos A) cos nn) ,

R2 + n4n4 v y /

+Y21

Appendix III

n1(y) : n2(y) :

Yl1^1 + Yl22^2 — Y11Y12 + Y21^4 + Y22^5 + YnY22^6 + + 2YHY21^7 + (Y11Y22 — Y12Y21) — 2Y12Y22^9 +

+ 2y2 (Y11Y21 — Y12Y22) ,

Yl22^1 + Yl21^2 + Y12Y11^3 + Y22^4 + Y21^5 — Y22Y21^6 + + 2Y12Y22^7 + (Y11Y22 — Y12Y21) ^8 — 2Y11 Y21^9 —

— öy2 (Y11Y21 — Y12Y22),

1

2 1

n3(y) : = 2

Ö (Y21 — Y12) ^3 + " (Y222 — Y21) ^6 + Y11Y12^10 +

2

+Y21Y22^11 + (Y11Y22 + Y12Y21) y,

2

:= :=

eV^Ry U + sin

8R

^2 :=

eV2Ry (2 - sin

8R

^3 := -

eV2Ry cos

4R ,

e-V2Ry (2 - sin V2Ry) e-V2Ry (2 + sin ^2Ry)

8R

cos

4R ,

:= -

8R

sin

2R

=V2Ry

sin

10 : =

4R

:= :=

e-V2Ry sin

e-V2Ry cos ^2Ry

4R ,

cos

4R ,

11

4R

Appendix IV

U (k) :

( 7n51 (k) + Y22q2 (k) - 71171253 (k) + \ +72154 (k) + 72255 (k) + 72172236 (k) + +271172137 (k) + (711722 - 712721) q8 (k) -

V

-271272259 (k)

/

V (k)

( Y^ (k) + 7n52 (k) + 71271153 (k) + \ +72254 (k) + 7I155 (k) - 72272156 (k) + +271272257 (k) + (711722 - 712721) 58 (k) --271172159 (k)

W (k) :

^ 711712510 (k) + 1 (7n - 7l22) 53 (k) + ^

21

+ (711722 + 721) 57 (k) + ^722 (721 - 712) 58 (k) + + (712721 + 722711) 59 (k) + 721722511 (k) + + 2 (7222 - 721) 56 (k)

K

n2 (2k + 1)2 ,

(

2 (8R2 + kR + k2) - eV2R

51 (k) = (2k + 1)n-

V

k2 (cos V2R + 1) + ^

+2Rk (cos + 2 sin V2R + +8R2 (2 + sin V2R)

2 (32R3 + k (16R2 + 2kr + k2))

2R (8R - k) - eV2R

52 (k)

(2k + 1)n-

(

\

K2 ^cos v/2R - 1J + +2Rk (2 sin + cos v^R) + +8R2 (sin - 2)

2 (32R3 + k (16R2 + 2kr + k2))

e^ ( 8R cos V2R — k sin V2R ) cos n(2k + 1) — 4R qa (k) = n(2k + 1)-^--J

16R2 + k2

e-V2n

q4 (k) = (2k + 1)n-

( k2 + cos + ^

+8R2 ^2 — sin V2R — — 2Rk (2 sin V2R — cos V2R)

+ 2 (8R2 + RK + K2)

2(32R3 + 16R2k + 2Rk2 + n6 (2k + 1)6)

^ (k2 (1 — cos V2R + 8R2 f 2 + sin V2R +\ e-V2R [ V J V ; I+2R (8R — K)

\ +2RK (2si^v/2R — co^V2R)

q5 (k) = (2k + 1)n-^-(-^-J--S--

y w v 2(32R3 + 16R2k + 2Rk2 + n6 (2k + 1)6)

e-V2R L sin V2R + 4R cos V2R + 2R q6 (k) = 2n(2k + 1)-^--'

qi (k)

qs(k) qo(k)

qio (k) = (2k + 1)n

16R2 + k2

^2R + k (cos V2R + 1)) — 2 (R — k) 2n(2k + 1) (k — 2R) ,

(1 — V2R + 2n(2k + 1)) sin

4R —2k ,

4R + k (cos — 1) 2n(2k + 1) (k — 2R) ,

k (1 + e/2R cos V2R + 4Re^2R sin

16R2 + k2 ,

(1 + e-V2R cos V2R — 4Re-V2R sin

k ( 1 + e v cos

qii (k) = (2k + 1)n

(16R2 + k2)

Appendix V

®i = YnTi + 722T2 — 7ii7i2Ta + y!^ + Y22 T5 + 72! 722T3+

+ 27ii72iTi + (7ii722 — 7i272i) TS — 27i2722Tg, @2 = 7i22Ti + 72iT2 + 7i27iiT3 + 7I2 T4 + 7|i T5 — Y22Y2iT6 + +2Yi2Y22Ti + (yiiY22 — Yi2Y2i) Ts — 2YiiY2iT9,

@3 = 1 (Yii — Y12) T3 + 2 (Y22 — Yii) T6 + YiiYi2Tio + Y2iY22Tii + YiiY22 + Yi2Y2i,

(2R (Ry sin 2t — Co(t)n2 (2k)2) + R (n2 (2k)2 K(t) + 8^R cos 2t)) cos 2kn

O4 - -(-^-

4kn (4R2 + n4 (2k)

2M3 + M1 + M2 _ 2M3 - M1 - M2 _ M2 - M1

T 1 — -, -, T 2 — -. -, T 3 — -. —,

^v/2R ^v^R ^v^R

= M5 - M4 - 2M9 T = M4 - M5 - 2M9 T = _ M4 + M5

4 = 4—2R , 5 = 4—2R , 6 = 2—2R ,

1 1 M8 M7 1

t7 = 1 m6 + —m7, t8 = -—ML, t9 = —M^ +1M6,

7 2 6 2—2R 7, 8 —2R 9 2—2R 2 6,

M1 + M2 M5 - M4

T10 = —"—, T11

^v/2R ' ^v/2R '

n2k2 (eV2R cos —2R - 1) + ReV2R sin —2R

cos

M1 = - kn

M4

2 (R2 + n4k4) R (1 - eV2R cos —2R + n^e^ sin —2R 2 (R2 + n4k4) ,

eV2R (—2R sin 2kn - 2kn) + 2kn = 4n2k2 + 2R ,

n2k2 (1 - 16n3k3e-V2R cos —2R + Re-V2R sin —2R 4 (R2 + n4k4) ,

n2k2e-V2R sin —2R + R (e-VR cos —2R - 1)

M2 = - kn

M3

M5 = - kn

2 (R2 + n4k4)

1 kn sin —2R (cos —- ^ _kn ^ - e-V2R)

M6 = - M7 = R - 2n2k2 , M8 = R - 2n2k2 , M9 = R + 2n2k2

Received for publication 21 January 2015