- Email: [email protected]
Simulation of the Taylor–Couette flow in a finite geometry by spectral element method
Applied Numerical Mathematics 33 (2000) 241–249
Simulation of the Taylor–Couette flow in a finite geometry by spectral element method E. Magère ∗ , M.O. Deville 1 EPFL, LMF/IMHEF/DGM, CH-1015 Lausanne, Switzerland
Abstract The purpose of the present numerical method is to simulate the flow in between two concentric cylinders of finite axial length. The numerical integration of the transient three-dimensional incompressible Navier–Stokes equations is performed by a Legendre spectral element method in space while the time marching scheme is fully explicit. By taking full advantage of the tensor product bases, it is possible to set up fast diagonalization techniques in order to speed up the solver. Numerical results are compared favorably with experimental data. 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Spectral element method; Direct numerical simulation; Flow transition
1. Introduction The Taylor–Couette flow occurring in between two coaxial differentially rotating cylinders is one of the fundamental problem in fluid mechanics. When only the inner cylinder is rotating and the outer one is fixed, a transition occurs from the laminar Couette flow to a periodic superposition of axisymmetric vortices, called the Taylor vortex flow, as the angular velocity is increased. Increasing further the velocity of the inner cylinder, the Taylor vortices are subject to a second transition and produce wavy vortices. However, when the cylinders are counter-rotating, the physical situation is less clear. Some experimental results are available [2], but further investigations need to be carried out. This paper aims to design a spectral element algorithm in order to handle this case. The paper is organized as follows. Section 2 describes the weak formulation of the Navier–Stokes (NS) equations. Section 3 presents the spatial discretization, while Section 4 treats the time marching scheme. In Section 5, the boundary conditions are set up to avoid the singularity between the inner cylinder and the top and bottom plates. The results are given in Section 6. ∗ Corresponding author. 1 E-mail: [email protected]
0168-9274/00/$20.00 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 8 9 - 6
242
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
2. Weak formulation of the basic equations The governing equations are the NS equations. The momentum and continuity equations are made dimensionless by introducing characteristic variables, the gap separating the cylinders, d, the viscous time d 2 /ν, where ν is the kinematic viscosity, the velocity of the outer cylinder of radius ro and angular velocity ωo , V = ro ωo , and the pressure µV /d, µ being the dynamic viscosity. The radius of the inner cylinder is ri and its angular velocity is ωi . Writing the basic equations in vector form, and taking the boundary and initial conditions into account, one is left with: ∂v + Re nl(v) = −∇p + ∇ 2 v ∂t
∀(r, t) ∈ D × [0, T ],
∇ ·v=0
(1) ∀(r, t) ∈ D × [0, T ], ∀(r, t) ∈ ∂D×]0, T ], v = vi ∀(r, t) ∈ D × {0}. The notation nl refers to the nonlinear term. The cylindrical coordinates, r = (r, θ, z), are defined in D =]ri , ro [×]0, 2π [×]0, L[. The symbol Re = dV /ν denotes R the Reynolds number. Introducing the scalar product of unit weight: hf, gi = D f g, we resort to a weak formulation of Eq. (1). We search the velocity solution of these equations in the Sobolev space [H1 (D)]3 = {Π v ∈ L2 (D), R∂Π v /∂r ∈ L2 (D)}3 and the pressure in the Lebesgue space with zero mean, L20 (D) = {Π p ∈ L2 (D), D Π p = 0}. v = vb
∂v ∂t
+ ∇p − s(v), Π
v
3
= 0 ∀Π v ∈ H1 (D) ,
∇ · v, Π p = 0
(2)
∀Π p ∈ L2 (D),
where s(v) = −Re nl(v) + ∇ 2 v, and Π v and Π p are high order polynomial test functions. 3. Spatial discretization The azimuthal direction being periodic, we first approximate the velocity v and the pressure p by Fourier series: v(r, t) =
NθX /2−1
ikθ
v k (r, z, t)e ,
p(r, t) =
k=−Nθ /2
NθX /2−1
pk (r, z, t)eikθ .
k=−Nθ /2 v
p
The test functions Π and Π are expressed as products of Fourier bases by continuous functions of [H1 (Ω)]3 and L2 (Ω), respectively, where Ω =]ri , ro [×]0, L[. We obtain a two-dimensionalSproblem for each of the Nθ Fourier modes. Ω is further decomposed in E rectangular elements: Ω = Ee=1 Ω e . Lagrange interpolation on the Gauss–Lobatto–Legendre (GLL) points is applied to the Fourier modes: v k (r, z, t)|Ω e =
Nz Nr X X
v eij k (t)φije (r, z),
i=0 j =0
where = hi (x)hj (y). (r(x), z(y)) represents the affine mapping from ]−1, 1[×]−1, 1[ to Ω e . The basis hi (respectively hj ) is the Lagrange interpolant of degree Nr (respectively Nz ) associated to the ith (respectively j th) GLL point of ]−1, 1[. In order to avoid spurious pressure modes we use the φije (r(x), z(y))
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
243
(PN /PN−2 ) formulation, but instead of taking the pressure at the Gauss–Legendre points and the velocity at the GLL points, we apply the modified version of Azaïez et al. [3] which consists in using the inner GLL points for the pressure. Each pressure Fourier mode is then expanded as follows: pk (r, z, t)|Ω e =
NX z −1 r −1 N X
pije k (t)ψije (r, z),
i=1 j =1
ei (x)h ej (y) and h ei (h ej ) is the Lagrange interpolant of degree Nr − 2 (Nz − 2) where ψije (r(x), z(y)) = h associated to the ith (j th) inner GLL point. We use a unique discrete inner product defined at the GLL points of each element:
hf, gid =
Nz Nr X E X X
fije gije rie ρi ρj ,
e=1 i=0 j =0
where rie = r(xi ), r belonging to Ω e and ρi and ρj represent the GLL weights. The pressure is therefore extrapolated on the interfaces of the elements to obtain the pressure gradient on the velocity grid. The discrete form of Eq. (2) is
∂v k
+ ∇pk − s k (v), φ ∂t h∇ · v k , ψid = 0
= 0 ∀φ ∈ (Xd )3 , d
(3)
∀ψ ∈ Yd .
In the previous equations, Xd = H1 (Ω) ∩ PN,E (Ω), Yd = L2 (Ω) ∩ PN−2,E (Ω), PN,E (Ω) = {φ, ∃e ∈ {1, . . . , E}, φ|Ω e ∈ PN (Ω e ), φ|Ω−Ω e = 0}. PN (Ω e ) is the space of polynomials of degree Nr in the radial direction and Nz in the axial direction defined on Ω e . The trial functions are chosen to be the same as the test functions. We introduce the following matrices locally to each Ω e . The mass matrices are defined as rB r
ij
= ρj r(xj )δij
e
ge , 2 e
Br r
= ij
ge 1 ρj δij r(xj ) 2
and
Bijz = ρj δij
he , 2
e
where g is the width of Ω and h , its height. The first-order derivative matrices are dijr =
∂hj (xi ), ∂x
dijz =
∂hj (yi ). ∂y
The second-order derivative matrices: Arij
Nr 2 X r r =− e r(xk )ρkr dk,i dk,j , g k=0
Azij
Nz 2 X =− e ρzd z d z . h k=0 k k,i k,j
The mass matrices incorporating the extrapolation from the pressure grid to the velocity grid are
e e ge fr ij = ρi h fj (xi ) g , fz ij = ρi h fj (yi ) h . B B , ij 2 2 2 The first-order derivative matrices going from the pressure space to the velocity space are
fr rB
fr = − G
fj (xi ) = ρi r(xi )h
2 r T fr fr d rB − B , ge
fz = − G
2 zT fz d B . he
244
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
We perform the direct stiffness [6] on these matrices. We use calligraphic letters for the resulting global matrices. The NS equations become: duk fr fz (rB r ) ⊗ B z + G ⊗ B pk − (rB r ) ⊗ B z skr = 0, dt (rB r ) ⊗ B z dvk + ik B fr ⊗ B fz pk − (rB r ) ⊗ B z s θ = 0, k dt r z dwk r z z fr fz (rB ) ⊗ B dt + (r B ) ⊗ G pk − (rB ) ⊗ B sk = 0, fr T ⊗ B fz T uk + ik B fr T ⊗ B fz T vk − (r B fr )T ⊗ G fz T wk = 0, −G
(4)
where ⊗ denotes the tensor product of matrices. The velocity components uk , vk and wk are defined on the grid Ωd made of all the GLL points in each of the elements. The pressure field pk is defined on the e d made of the inner GLL of all the elements. We now define grid Ω fr fz G ⊗B e fr ⊗ B fz G = ik B fr ⊗ G fz rB
and
rB r ⊗ B z B = rB r ⊗ B z . rB r ⊗ B z
In compact matrix form, the discrete equations are now
◦ dv k e k − Bs k = 0 in Ω d , + Gp dt T ed . Ge V k = 0 in Ω
B
(5) ◦
In (5), v k = (uk , vk , wk ) are the three components of the velocity defined in Ω d , the finite set that contains the grid points of Ωd , the boundary points excluded. The notation V k represents the three components of the velocity defined in Ωd .
4. Time discretization A fully explicit treatment of the source term, s k (V ) = −Re nlk (V ) + B −1 Av k , is chosen instead of the classical implicit linear viscous/explicit non linear decomposition. Here, A denotes the Laplacian operator. In order to be able to carry out long term integration over large time scales, the skew-symmetric form of the nonlinear term is chosen. In the general case of projection methods, a Helmholtz operator for the velocity has to be inverted. In cylindrical coordinates, this operator cannot be expressed in a separate form. Hence, the fast diagonalization technique [4,5] cannot be applied to invert this matrix. An iterative method is then needed. For a three-dimensional problem with one periodic direction, each matrix– vector multiplication involved in the algorithm requires 2N 4 operations (1 addition and 1 multiplication counts for one operation), with N 3 , the total number of grid points. Each of the iteration steps can be achieved with the order of 10 matrix–vector multiplications. In the case of the totally explicit choice, the pseudo-Laplacian operator for the pressure can be inverted with the fast diagonalization method and requires therefore 4N 4 operations. The solver in itself is more efficient in the fully explicit case than in the implicit/explicit one. However, the time step requirements are more stringent for the former. Indeed, when an explicit treatment of the discrete NS equations is undertaken, two stability constraints
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
245
have to be enforced on the time step 1t. The first condition is imposed on the diffusive part of the operator −2 −2 1t 6 CSD Re 1rmin + (r1θ)−2 min + 1zmin
−1
,
CSD < 0.4,
(6)
while the second condition comes from the nonlinear term (CFL condition)
umax vmax wmax −1 1t 6 CSN + + , CSN < 0.2. (7) 1rmin r1θmin 1zmin Usually, 1rmin = 1zmin and, in the case of the Taylor–Couette geometry, 1rmin (r1θ)min . Moreover, we have for a Legendre GLL grid: 1rmin ∼ 1/N 2 . Hence, the viscous constraint (6) can be written as: 1tmax ∝ Re/N 4 . Let us assume umax ∼ wmax ∼ 1, then we obtain from (7) 1tmax ∝ 1/N 2 . For a polynomial degree, N ' 10, and a Reynolds number, Re ' 100, the diffusion constraint is of the same order as the CFL constraint. Consequently, an explicit treatment of the linear viscous term is possible, for the Reynolds range we are interested in: between 100 and 1000. In practice, at a Reynolds number of 100, the viscous constraint imposes a time step about three times smaller. Dropping the indices k of the Fourier modes, the semi-discrete NS equations give rise to the set of ordinary differential equations: ◦ dv = −B −1 Gp e + s(V ) in Ω d , (8) dt T e ed . G V =0 in Ω We have chosen a second-order Runge–Kutta (RK2) scheme for the time discretization. It reads as follows: n+1 e n+1/2 + s n+1/2 , v = v n + 1t −B −1 Gp (9) GeT V n+1 = 0, with 1t 1t n+1/2 n n+1/2 n p =p t + , V =V t + and s n+1/2 = s V n+1/2 . 2 2 n+1/2 Using the fractional step method, we first obtain V from V n and then V n+1 from V n+1/2 . The first step can be decomposed in three stages: a prediction, where we obtain v ∗ = v n + (1t/2)s n which is not divergence-free, a projection on the divergence-free velocity space, where we get the pressure 2 eT v∗ n T −1 e −1 e p = G B G G , (10) n+1/2 vb 1t and, finally, a correction step, where we find the divergence-free velocity 1t −1 e n v n+1/2 = v ∗ − B Gp . 2 So that, n+1/2 v V n+1/2 = n+1/2 . vb The second and last step is also performed in three successive stages. We first obtain: v ∗ = v n + 1t s(V n+1/2 ), then: ∗ −1 1 T v p n+1/2 = Ge T B −1 Ge Ge , (11) v n+1 1t b
246
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
Fig. 1. ωi = 1, ωo = −2, ri = 7, d = 1, ε = 0.05d. e n+1/2 . So that, and, finally, we find: v n+1 = v ∗ − 1t B −1 Gp n+1 v V n+1 = . v n+1 b ◦
The projection equation is obtained from the combination of the momentum equation in Ω d with the ◦ boundary conditions on Ωd − Ω d . (
e n+1/2 v n+1 = v ∗ − 1t B −1 Gp v n+1 = v n+1 b
◦
in Ω d , ◦ on Ωd − Ω d .
(12)
We then apply the divergence operator, GeT , to this equation to recover (11). 5. Boundary conditions We impose Dirichlet boundary conditions on the velocity. The rotation velocity of the inner cylinder, ωi , is different from the rotation velocity common to the outer cylinder and the end plates, ωo , resulting in a singularity at (r, z) = (ri , 0) and (r, z) = (ri , L). We replace the rotation velocity of the end plates, ωo , by ω(r) on [ri , ri + ε], ε d, where ω(r) = ωC (r)fi (r) + ωo fo (r), with fi and fo two C 1 functions defined so that ω evolves smoothly from ωC to ωo on [ri , ri + ε] (see Fig. 1), ωC being the rotation velocity of the circular Couette flow defined between the radii ri and ri + ε: ωC (r) = a +
b , r2
a=
(ri + ε)2 ωo − ri2 ωi (ri + ε)2 − ri2
and
dω dωC (ri ) = (ri ), dr dr
ω(ri + ε) = ωo
b=
(ωi − ωo )ri2 (ri + ε)2 . (ri + ε)2 − ri2
We impose: ω(ri ) = ωC (ri ) = ωi , We choose:
π ri + ε − r π ri + ε − r fi,o = Ai,o sin + Bi,o cos . 2 ε 2 ε The four constraints determine the constants Ai,o and Bi,o .
and
dω (ri + ε) = 0. dr
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
247
Fig. 2. Stream traces in the meridian plane. Re = 165. Γ = L/d = 1.281.
Fig. 3. – experiments of Aitta et al. – simulations of Streett et al. 1 – ours.
6. Results To validate our code, we have reproduced the simulations of Streett et al. [7] in a short aspect ratio Taylor–Couette geometry (see Figs. 2 and 3). Our simulations compare well to theirs, and also to the experiments of Aitta et al. [1]. We observe a sub-critical transition from a symmetric to an asymmetric flow. They define an asymmetry parameter, RL
ψ = R L0 0
w(r, z) dz
|w(r, z)| dz
,
to track this transition. In Fig. 4, we also show the beginning of the instability in the counter-rotating case. The aspect ratio is Γ = 12, The inner Reynolds number, based on the velocity of the inner cylinder, is 400 and the outer Reynolds number, based on the outer cylinder velocity, is −612. The polynomial degree we use
248
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
Fig. 4. Stream traces and azimuthal velocity contour in the meridian plane.
is Nr = Nz = 12, the number of elements is E = 3 × 26. The number of Fourier modes is 32. The instability starts at mid-height and propagates towards the end plates. The region of instability, as given by the non-viscous theory, is confined near the inner cylinder. The flow forms spirals after a sufficient time.
7. Conclusions In this paper, we have detailed a numerical algorithm based on Fourier-spectral element discretization in space. The velocity nodes are the GLL nodes while the pressure nodes are the inner GLL points. This formulation avoids the spurious pressure modes. The explicit time discretization relies on a two-stage Runge–Kutta scheme. In order to speed up the pressure calculation, a fast diagonalization method is devised. The numerical method is well adapted to the problem at hand. The algorithm is efficient, the CPU time spent on a CRAY J90 per time step, per node is 1.3 · 10−5 s. The physics is resolved correctly in space and time. A good agreement is found between numerical results and experimental data. References [1] A. Aitta, G. Ahlers, D.S. Cannel, Tricritical phenomena in rotating Couette–Taylor flow, Phys. Rev. Lett. 54 (1985) 673–676. [2] C.D. Andereck, S.S. Liu, H.L. Swinney, Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech. 164 (1986) 155–183. [3] M. Azaïez, A. Fikri, G. Labrosse, A unique grid spectral solver of the nD Cartesian unsteady Stokes system. Illustrative numerical results, Finite Elements in Analysis and Design 16 (1994) 247–260. [4] D.B. Haidvogel, T. Zang, The accurate solution of Poisson’s equation by expansion in Chebyshev polynomials, J. Comput. Phys. 30 (1979) 167–180.
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
249
[5] R.E. Lynch, J.R. Rice, D.H. Thomas, Direct solution of partial difference equations by tensor product methods, Numer. Math. 6 (1964) 185–199. [6] E.M. Rønquist, Spectral element methods for the unsteady Navier–Stokes equations, Von Karman Institute Lecture Series, 1991. [7] C.L. Streett, M.Y. Hussaini, A numerical simulation of finite-length Taylor–Couette flow, Comput. Fluid Dynamics (1988) 663–675.
Simulation of the Taylor–Couette flow in a finite geometry by spectral element method E. Magère ∗ , M.O. Deville 1 EPFL, LMF/IMHEF/DGM, CH-1015 Lausanne, Switzerland
Abstract The purpose of the present numerical method is to simulate the flow in between two concentric cylinders of finite axial length. The numerical integration of the transient three-dimensional incompressible Navier–Stokes equations is performed by a Legendre spectral element method in space while the time marching scheme is fully explicit. By taking full advantage of the tensor product bases, it is possible to set up fast diagonalization techniques in order to speed up the solver. Numerical results are compared favorably with experimental data. 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Spectral element method; Direct numerical simulation; Flow transition
1. Introduction The Taylor–Couette flow occurring in between two coaxial differentially rotating cylinders is one of the fundamental problem in fluid mechanics. When only the inner cylinder is rotating and the outer one is fixed, a transition occurs from the laminar Couette flow to a periodic superposition of axisymmetric vortices, called the Taylor vortex flow, as the angular velocity is increased. Increasing further the velocity of the inner cylinder, the Taylor vortices are subject to a second transition and produce wavy vortices. However, when the cylinders are counter-rotating, the physical situation is less clear. Some experimental results are available [2], but further investigations need to be carried out. This paper aims to design a spectral element algorithm in order to handle this case. The paper is organized as follows. Section 2 describes the weak formulation of the Navier–Stokes (NS) equations. Section 3 presents the spatial discretization, while Section 4 treats the time marching scheme. In Section 5, the boundary conditions are set up to avoid the singularity between the inner cylinder and the top and bottom plates. The results are given in Section 6. ∗ Corresponding author. 1 E-mail: [email protected]
0168-9274/00/$20.00 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 8 9 - 6
242
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
2. Weak formulation of the basic equations The governing equations are the NS equations. The momentum and continuity equations are made dimensionless by introducing characteristic variables, the gap separating the cylinders, d, the viscous time d 2 /ν, where ν is the kinematic viscosity, the velocity of the outer cylinder of radius ro and angular velocity ωo , V = ro ωo , and the pressure µV /d, µ being the dynamic viscosity. The radius of the inner cylinder is ri and its angular velocity is ωi . Writing the basic equations in vector form, and taking the boundary and initial conditions into account, one is left with: ∂v + Re nl(v) = −∇p + ∇ 2 v ∂t
∀(r, t) ∈ D × [0, T ],
∇ ·v=0
(1) ∀(r, t) ∈ D × [0, T ], ∀(r, t) ∈ ∂D×]0, T ], v = vi ∀(r, t) ∈ D × {0}. The notation nl refers to the nonlinear term. The cylindrical coordinates, r = (r, θ, z), are defined in D =]ri , ro [×]0, 2π [×]0, L[. The symbol Re = dV /ν denotes R the Reynolds number. Introducing the scalar product of unit weight: hf, gi = D f g, we resort to a weak formulation of Eq. (1). We search the velocity solution of these equations in the Sobolev space [H1 (D)]3 = {Π v ∈ L2 (D), R∂Π v /∂r ∈ L2 (D)}3 and the pressure in the Lebesgue space with zero mean, L20 (D) = {Π p ∈ L2 (D), D Π p = 0}. v = vb
∂v ∂t
+ ∇p − s(v), Π
v
3
= 0 ∀Π v ∈ H1 (D) ,
∇ · v, Π p = 0
(2)
∀Π p ∈ L2 (D),
where s(v) = −Re nl(v) + ∇ 2 v, and Π v and Π p are high order polynomial test functions. 3. Spatial discretization The azimuthal direction being periodic, we first approximate the velocity v and the pressure p by Fourier series: v(r, t) =
NθX /2−1
ikθ
v k (r, z, t)e ,
p(r, t) =
k=−Nθ /2
NθX /2−1
pk (r, z, t)eikθ .
k=−Nθ /2 v
p
The test functions Π and Π are expressed as products of Fourier bases by continuous functions of [H1 (Ω)]3 and L2 (Ω), respectively, where Ω =]ri , ro [×]0, L[. We obtain a two-dimensionalSproblem for each of the Nθ Fourier modes. Ω is further decomposed in E rectangular elements: Ω = Ee=1 Ω e . Lagrange interpolation on the Gauss–Lobatto–Legendre (GLL) points is applied to the Fourier modes: v k (r, z, t)|Ω e =
Nz Nr X X
v eij k (t)φije (r, z),
i=0 j =0
where = hi (x)hj (y). (r(x), z(y)) represents the affine mapping from ]−1, 1[×]−1, 1[ to Ω e . The basis hi (respectively hj ) is the Lagrange interpolant of degree Nr (respectively Nz ) associated to the ith (respectively j th) GLL point of ]−1, 1[. In order to avoid spurious pressure modes we use the φije (r(x), z(y))
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
243
(PN /PN−2 ) formulation, but instead of taking the pressure at the Gauss–Legendre points and the velocity at the GLL points, we apply the modified version of Azaïez et al. [3] which consists in using the inner GLL points for the pressure. Each pressure Fourier mode is then expanded as follows: pk (r, z, t)|Ω e =
NX z −1 r −1 N X
pije k (t)ψije (r, z),
i=1 j =1
ei (x)h ej (y) and h ei (h ej ) is the Lagrange interpolant of degree Nr − 2 (Nz − 2) where ψije (r(x), z(y)) = h associated to the ith (j th) inner GLL point. We use a unique discrete inner product defined at the GLL points of each element:
hf, gid =
Nz Nr X E X X
fije gije rie ρi ρj ,
e=1 i=0 j =0
where rie = r(xi ), r belonging to Ω e and ρi and ρj represent the GLL weights. The pressure is therefore extrapolated on the interfaces of the elements to obtain the pressure gradient on the velocity grid. The discrete form of Eq. (2) is
∂v k
+ ∇pk − s k (v), φ ∂t h∇ · v k , ψid = 0
= 0 ∀φ ∈ (Xd )3 , d
(3)
∀ψ ∈ Yd .
In the previous equations, Xd = H1 (Ω) ∩ PN,E (Ω), Yd = L2 (Ω) ∩ PN−2,E (Ω), PN,E (Ω) = {φ, ∃e ∈ {1, . . . , E}, φ|Ω e ∈ PN (Ω e ), φ|Ω−Ω e = 0}. PN (Ω e ) is the space of polynomials of degree Nr in the radial direction and Nz in the axial direction defined on Ω e . The trial functions are chosen to be the same as the test functions. We introduce the following matrices locally to each Ω e . The mass matrices are defined as rB r
ij
= ρj r(xj )δij
e
ge , 2 e
Br r
= ij
ge 1 ρj δij r(xj ) 2
and
Bijz = ρj δij
he , 2
e
where g is the width of Ω and h , its height. The first-order derivative matrices are dijr =
∂hj (xi ), ∂x
dijz =
∂hj (yi ). ∂y
The second-order derivative matrices: Arij
Nr 2 X r r =− e r(xk )ρkr dk,i dk,j , g k=0
Azij
Nz 2 X =− e ρzd z d z . h k=0 k k,i k,j
The mass matrices incorporating the extrapolation from the pressure grid to the velocity grid are
e e ge fr ij = ρi h fj (xi ) g , fz ij = ρi h fj (yi ) h . B B , ij 2 2 2 The first-order derivative matrices going from the pressure space to the velocity space are
fr rB
fr = − G
fj (xi ) = ρi r(xi )h
2 r T fr fr d rB − B , ge
fz = − G
2 zT fz d B . he
244
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
We perform the direct stiffness [6] on these matrices. We use calligraphic letters for the resulting global matrices. The NS equations become: duk fr fz (rB r ) ⊗ B z + G ⊗ B pk − (rB r ) ⊗ B z skr = 0, dt (rB r ) ⊗ B z dvk + ik B fr ⊗ B fz pk − (rB r ) ⊗ B z s θ = 0, k dt r z dwk r z z fr fz (rB ) ⊗ B dt + (r B ) ⊗ G pk − (rB ) ⊗ B sk = 0, fr T ⊗ B fz T uk + ik B fr T ⊗ B fz T vk − (r B fr )T ⊗ G fz T wk = 0, −G
(4)
where ⊗ denotes the tensor product of matrices. The velocity components uk , vk and wk are defined on the grid Ωd made of all the GLL points in each of the elements. The pressure field pk is defined on the e d made of the inner GLL of all the elements. We now define grid Ω fr fz G ⊗B e fr ⊗ B fz G = ik B fr ⊗ G fz rB
and
rB r ⊗ B z B = rB r ⊗ B z . rB r ⊗ B z
In compact matrix form, the discrete equations are now
◦ dv k e k − Bs k = 0 in Ω d , + Gp dt T ed . Ge V k = 0 in Ω
B
(5) ◦
In (5), v k = (uk , vk , wk ) are the three components of the velocity defined in Ω d , the finite set that contains the grid points of Ωd , the boundary points excluded. The notation V k represents the three components of the velocity defined in Ωd .
4. Time discretization A fully explicit treatment of the source term, s k (V ) = −Re nlk (V ) + B −1 Av k , is chosen instead of the classical implicit linear viscous/explicit non linear decomposition. Here, A denotes the Laplacian operator. In order to be able to carry out long term integration over large time scales, the skew-symmetric form of the nonlinear term is chosen. In the general case of projection methods, a Helmholtz operator for the velocity has to be inverted. In cylindrical coordinates, this operator cannot be expressed in a separate form. Hence, the fast diagonalization technique [4,5] cannot be applied to invert this matrix. An iterative method is then needed. For a three-dimensional problem with one periodic direction, each matrix– vector multiplication involved in the algorithm requires 2N 4 operations (1 addition and 1 multiplication counts for one operation), with N 3 , the total number of grid points. Each of the iteration steps can be achieved with the order of 10 matrix–vector multiplications. In the case of the totally explicit choice, the pseudo-Laplacian operator for the pressure can be inverted with the fast diagonalization method and requires therefore 4N 4 operations. The solver in itself is more efficient in the fully explicit case than in the implicit/explicit one. However, the time step requirements are more stringent for the former. Indeed, when an explicit treatment of the discrete NS equations is undertaken, two stability constraints
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
245
have to be enforced on the time step 1t. The first condition is imposed on the diffusive part of the operator −2 −2 1t 6 CSD Re 1rmin + (r1θ)−2 min + 1zmin
−1
,
CSD < 0.4,
(6)
while the second condition comes from the nonlinear term (CFL condition)
umax vmax wmax −1 1t 6 CSN + + , CSN < 0.2. (7) 1rmin r1θmin 1zmin Usually, 1rmin = 1zmin and, in the case of the Taylor–Couette geometry, 1rmin (r1θ)min . Moreover, we have for a Legendre GLL grid: 1rmin ∼ 1/N 2 . Hence, the viscous constraint (6) can be written as: 1tmax ∝ Re/N 4 . Let us assume umax ∼ wmax ∼ 1, then we obtain from (7) 1tmax ∝ 1/N 2 . For a polynomial degree, N ' 10, and a Reynolds number, Re ' 100, the diffusion constraint is of the same order as the CFL constraint. Consequently, an explicit treatment of the linear viscous term is possible, for the Reynolds range we are interested in: between 100 and 1000. In practice, at a Reynolds number of 100, the viscous constraint imposes a time step about three times smaller. Dropping the indices k of the Fourier modes, the semi-discrete NS equations give rise to the set of ordinary differential equations: ◦ dv = −B −1 Gp e + s(V ) in Ω d , (8) dt T e ed . G V =0 in Ω We have chosen a second-order Runge–Kutta (RK2) scheme for the time discretization. It reads as follows: n+1 e n+1/2 + s n+1/2 , v = v n + 1t −B −1 Gp (9) GeT V n+1 = 0, with 1t 1t n+1/2 n n+1/2 n p =p t + , V =V t + and s n+1/2 = s V n+1/2 . 2 2 n+1/2 Using the fractional step method, we first obtain V from V n and then V n+1 from V n+1/2 . The first step can be decomposed in three stages: a prediction, where we obtain v ∗ = v n + (1t/2)s n which is not divergence-free, a projection on the divergence-free velocity space, where we get the pressure 2 eT v∗ n T −1 e −1 e p = G B G G , (10) n+1/2 vb 1t and, finally, a correction step, where we find the divergence-free velocity 1t −1 e n v n+1/2 = v ∗ − B Gp . 2 So that, n+1/2 v V n+1/2 = n+1/2 . vb The second and last step is also performed in three successive stages. We first obtain: v ∗ = v n + 1t s(V n+1/2 ), then: ∗ −1 1 T v p n+1/2 = Ge T B −1 Ge Ge , (11) v n+1 1t b
246
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
Fig. 1. ωi = 1, ωo = −2, ri = 7, d = 1, ε = 0.05d. e n+1/2 . So that, and, finally, we find: v n+1 = v ∗ − 1t B −1 Gp n+1 v V n+1 = . v n+1 b ◦
The projection equation is obtained from the combination of the momentum equation in Ω d with the ◦ boundary conditions on Ωd − Ω d . (
e n+1/2 v n+1 = v ∗ − 1t B −1 Gp v n+1 = v n+1 b
◦
in Ω d , ◦ on Ωd − Ω d .
(12)
We then apply the divergence operator, GeT , to this equation to recover (11). 5. Boundary conditions We impose Dirichlet boundary conditions on the velocity. The rotation velocity of the inner cylinder, ωi , is different from the rotation velocity common to the outer cylinder and the end plates, ωo , resulting in a singularity at (r, z) = (ri , 0) and (r, z) = (ri , L). We replace the rotation velocity of the end plates, ωo , by ω(r) on [ri , ri + ε], ε d, where ω(r) = ωC (r)fi (r) + ωo fo (r), with fi and fo two C 1 functions defined so that ω evolves smoothly from ωC to ωo on [ri , ri + ε] (see Fig. 1), ωC being the rotation velocity of the circular Couette flow defined between the radii ri and ri + ε: ωC (r) = a +
b , r2
a=
(ri + ε)2 ωo − ri2 ωi (ri + ε)2 − ri2
and
dω dωC (ri ) = (ri ), dr dr
ω(ri + ε) = ωo
b=
(ωi − ωo )ri2 (ri + ε)2 . (ri + ε)2 − ri2
We impose: ω(ri ) = ωC (ri ) = ωi , We choose:
π ri + ε − r π ri + ε − r fi,o = Ai,o sin + Bi,o cos . 2 ε 2 ε The four constraints determine the constants Ai,o and Bi,o .
and
dω (ri + ε) = 0. dr
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
247
Fig. 2. Stream traces in the meridian plane. Re = 165. Γ = L/d = 1.281.
Fig. 3. – experiments of Aitta et al. – simulations of Streett et al. 1 – ours.
6. Results To validate our code, we have reproduced the simulations of Streett et al. [7] in a short aspect ratio Taylor–Couette geometry (see Figs. 2 and 3). Our simulations compare well to theirs, and also to the experiments of Aitta et al. [1]. We observe a sub-critical transition from a symmetric to an asymmetric flow. They define an asymmetry parameter, RL
ψ = R L0 0
w(r, z) dz
|w(r, z)| dz
,
to track this transition. In Fig. 4, we also show the beginning of the instability in the counter-rotating case. The aspect ratio is Γ = 12, The inner Reynolds number, based on the velocity of the inner cylinder, is 400 and the outer Reynolds number, based on the outer cylinder velocity, is −612. The polynomial degree we use
248
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
Fig. 4. Stream traces and azimuthal velocity contour in the meridian plane.
is Nr = Nz = 12, the number of elements is E = 3 × 26. The number of Fourier modes is 32. The instability starts at mid-height and propagates towards the end plates. The region of instability, as given by the non-viscous theory, is confined near the inner cylinder. The flow forms spirals after a sufficient time.
7. Conclusions In this paper, we have detailed a numerical algorithm based on Fourier-spectral element discretization in space. The velocity nodes are the GLL nodes while the pressure nodes are the inner GLL points. This formulation avoids the spurious pressure modes. The explicit time discretization relies on a two-stage Runge–Kutta scheme. In order to speed up the pressure calculation, a fast diagonalization method is devised. The numerical method is well adapted to the problem at hand. The algorithm is efficient, the CPU time spent on a CRAY J90 per time step, per node is 1.3 · 10−5 s. The physics is resolved correctly in space and time. A good agreement is found between numerical results and experimental data. References [1] A. Aitta, G. Ahlers, D.S. Cannel, Tricritical phenomena in rotating Couette–Taylor flow, Phys. Rev. Lett. 54 (1985) 673–676. [2] C.D. Andereck, S.S. Liu, H.L. Swinney, Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech. 164 (1986) 155–183. [3] M. Azaïez, A. Fikri, G. Labrosse, A unique grid spectral solver of the nD Cartesian unsteady Stokes system. Illustrative numerical results, Finite Elements in Analysis and Design 16 (1994) 247–260. [4] D.B. Haidvogel, T. Zang, The accurate solution of Poisson’s equation by expansion in Chebyshev polynomials, J. Comput. Phys. 30 (1979) 167–180.
E. Magère, M.O. Deville / Applied Numerical Mathematics 33 (2000) 241–249
249
[5] R.E. Lynch, J.R. Rice, D.H. Thomas, Direct solution of partial difference equations by tensor product methods, Numer. Math. 6 (1964) 185–199. [6] E.M. Rønquist, Spectral element methods for the unsteady Navier–Stokes equations, Von Karman Institute Lecture Series, 1991. [7] C.L. Streett, M.Y. Hussaini, A numerical simulation of finite-length Taylor–Couette flow, Comput. Fluid Dynamics (1988) 663–675.