Proceedings of ITC/ISHW2007

Critical transition model of edge shear flow formation
Luis GARCIA1) , José M. DELGADO1) , Benjamin A. CARRERAS2) , Iván CALVO3) ,M. Angeles PEDROSA3) and Carlos HIDALGO3)
1)

Universidad Carlos III, 28911 Leganés, Madrid, Spain
2)
BACV Solutions Inc., Oak Ridge, TN 37830, USA
3)
Asociación EURATOM-CIEMAT, 28040 Madrid, Spain
October 15, 2007
Recently, the experimental results for the emergence of the plasma shear flow layer in TJ-II have been
explained as a second-order phase transition like process by using a simple model of envelope equations forthe fluctuation level, the averaged poloidal velocity shear, and the pressure gradient [Phys. Plasmas 13, 122509(2006)]. Here, we extend this model by incorporating radial coupling. The model is applied to the study ofthe turbulence-shear flow interaction when the energy flux is low. Transition dynamics and their concomitantthresholds are examined within the context of this model. The effect of an external torque has also been
considered. In particular, we analyze the damping rate of the shear flow once the external torque has been
removed.
Keywords: Anomalous transport, Shear flow, Phase transition

1Introduction

The importance of shear flows in magnetically confinedplasmas is widely recognized. This in part due to the rolethat these flows play in improving confinement by turbulence suppression [1, 2]. A current problem in plasma
physics is to understand how shear flows are created andhow they interact with turbulence. This is a complex problem that can be approached at different levels, going fromdirect numerical simulations of the plasma to the development of simplified reduced models. In the present paperwe follow the latter approach.
The spontaneous formation of shear flows with nontrivial radial structure is studied. Due to the lack of sufficiently high spatial and temporal resolution, it is difficult
to formulate an experimental test that will definitively select the dominant mechanisms responsible for the transition. However recent TJ-II experimental results offer thepossibility of improved diagnostics with increased radialresolution during and after the transition [3]. Motivated bythis, we investigate in this paper the radial structure of theturbulence and poloidal flow in the context of phase transition models.
When the input power is low, we have a transition
from L-mode to a different regime. In this regime, the fluctuation level decreases or increases at a slower rate than theinput flux, and shear flow is spontaneously developed. We
focus our attention on this transition in our starting model.
The model we use for our study is a modified versionof the fluctuation-flow model with radial structure [6]. The
fluctuation-flow model consists of three coupled partial
differential equations of the reaction-diffusion type withauthor’s e-mail: lgarcia@fis.uc3m.es

nonlinear diffusivities for the averaged poloidal velocity
shear, the envelope of the turbulence fluctuations level, andthe pressure. We include the dependence of the coefficientson the pressure gradient [7]. We have derived an expression for the Reynolds stress term in the averaged poloidal
flow equation which conserves angular momentum.

2

Transport model

The model is formulated in terms of three
E1/2the averD fields:, whereaged turbulence fluctuation level E ≡ (ñ/n0 )2n0 is the equilibrium plasma density and ñ is the fluctuation
density; the poloidal flow shear, σ ≡ ∂ hVθ i /∂r, where hVθ idenotes poloidal and toroidal average over a magnetic flux
surface; and the averaged pressure, p.
∂E=N 2/3 E − N −1/2 E 2 − N −1/3 σ2 E
∂t
"
#
∂E
∂
(D1 E + D0 )+∂x∂x∂σ∂2 
= − σ − α3 2 N −4/3 E 2 σ∂t∂x" ∂2 σ #∂2 −5/3 2− 2 D2 N
E + D3∂x∂x2"#∂p∂∂p(D1 E + D0 )=S +,∂t∂x
∂x

(1)

(2)

(3)

where N ≡ |∂p/∂x|. This is a one-dimensional model inwhich quantities are assumed to depend only on a normalized radial coordinate, x. The equations are written interms of dimensionless variables, and x goes from 0 to 1,
which corresponds to the shear layer region at the plasmaedge.

Proceedings of ITC/ISHW2007
The first term in the equation for E represents a pressure driven, linear growth of the fluctuations. The second
term on the right-hand side models the saturation of turbulence in absence of shear flow, and the third term models
the turbulence suppression by shear flow. Here we assume
low energy flux so we can neglect diamagnetic effects. The
first term on the right-hand side of equation (2) models
the poloidal flow shear damping due to magnetic pumping, while the other two terms are contributions from theReynolds stress. The second term on the right-hand sideof equation (2) is a negative viscosity term and represents
the generation of poloidal flow by Reynolds stress, andthe third term is a hyperdiffusivity and damps the poloidalflow.
We include the dependence of the coefficients on thepressure gradient as in Ref. [7]. Because of the bad magnetic field line curvature at the stellarator edge, we assume
that the basic instability underlying the turbulence at theedge of TJ-II is the resistive interchange mode. From thedependence of the coefficients on the linear growth rate, γ,and mode width, Wk , we can derive the dependence of thecoefficients on the pressure gradient. We assume that allother dependencies are weak and we take the coefficientsto depend only on N. The expression for Reynolds stresscontributions is derived by a quasi-linear calculation similar to he one proposed in Ref. [8] and has two terms. The
one responsible for the generation of flow is a negative viscosity term and the other is a hyper-viscosity term. The
Reynolds stress contributions to the poloidal shear flow
equation conserve angular momentum.
The terms involving spatial derivatives are diffusion
terms. In equations (1) and (3) these terms have the standard Fick’s Law structure, with D0 representing the collisional diffusion and D1 E representing a renormalized
turbulent diffusion. The structure of the negative diffusion operator in equation (2) is due to the fact that is anequation for the shear, the derivative of the momentum.
Hyper-diffusivity, with fourth order in the derivative of theshear, keeps a relation two in the orders of source/diffusion
derivatives.
We assume that the energy source term S (x) is zero inthis layer and that the system is driven by an energy flux Γ0from the core, which determines the boundary condition at
x = 0 according to:∂pΓ0 = − (D1 E + D0 ),∂x 0p(1) = 0(4)For the other equations we use zero derivative boundary
conditions.

Apart from the trivial solution, E = σ = 0, Equations(1) and (2) have one fixed point solution with σ = 0, E =E0 , and N = N0 . From Eq. (1),E0 = N07/6(6)

By substitution of Eq. (6) in the relation (5), we obtainthe (constant) density gradient for the fixed point solution,
D1 N013/6 + D0 N0 − Γ0 = 0.

(7)

For solutions close to the fixed point is possible to derive a simplified description of the system using a multiple
scale perturbation analysis. As a first step, we introduce asmall parameter δ representing the size of the perturbation,and consider the following expansion
σ = δσ1 ,E = E0 + δ2 E2 ,N = N0 + δ2 N2(8)By substituting the expansion in Eq. (2), and taking
into account Eq. (6), we get at first order ∂4 σ1∂σ1∂2 σ1 = −σ1 −α3 N0 2 − D2 N02/3 + D3(9)
∂t
∂x
∂x4
We try as solutions modes like σ1 = σ10 cos(kπx)
which satisfy the boundary conditions. To study their stability properties, we consider a temporal and spatial dependence likeσ(x, t) = δσ10 eγt cos(kπx)(10)

Then, we get from Eq. (9),
γ = −1 + α3 N0 (kπ)2 − D2 N02/3 + D3 (kπ)4(11)This means that the range of possible unstable modes
is given by the relationk− < k < k+ ,(12)where(k± π)2 =α3 N0 ±q(α3 N0 )2 − 4 D2 N02/3 + D3


2 D2 N02/3 + D3(13)To have instability, there should be an integer k between k− and k+ . This gives us a critical value for N, Nc ,and, consequently, a threshold value for the flux, Γc .
By applying Eq. (13) to different k-values, we can obtain the threshold for Γc . We have done a scan in D2 withk going from 1 to 7 (Fig. 1), the rest of parameters being

Analytic and numerical solutions

D0 = 10−3 , D1 = 10−2 , D3 = 10−6 , α3 = 0.0175

For stationary solutions, Eq. (3) is readily integrated, giving a relation between E(x) and N(x),We look now for stationary solutions with a given kvalue, that is, with the poloidal flow shear at lowest order given by σ(x) = σ10 cos(kπx). From the stability
properties of Eq. (2), the stationary solution is the fixed3(D1 E + D0 ) N = Γ0(5)Proceedings of ITC/ISHW2007
101

3.5

k=7

k=6

100

k=5

k=4

k=3

3

k=2
k=1

2.5-1

10Γc

σ10

2
1.5

10-2

110-3
0.5

10-4

Analy
Num

010-5

10-4

10-3

1

1.5

D2Fig. 1 Flux threshold vs. D2 for different k-values. The rest ofparameters are given in the text.

point solution (7) for values of the flux below Γc . Since
σ2 = σ210 [1 + cos(2kπx)]/2, we get from Eqs. (1) and(5) that at first order, E(x) = E s + E21 cos(2kπx), andN(x) = N s + N21 cos(2kπx) for values of the flux above Γc .
By substituting these expansions in Eqs. (2) and (5), andtaking into account that D0  D1 E, and D3  D2 N −5/3 E 2for the parameters we are using, we get the following expressions,!2
Γ0 5,(14)E s = EcΓc!3Γc Γ0 5Ns =,(15)D1 Ec Γc
!3!− 1!37 Γc Γ0 5
Γc 6 Γ0 10 26σ10 = 2 − Ec (16)D1 Ec ΓcD1Γc 
where− 53Ec = (kπ)!2 !− 132 − 103Γc 5 Γα3 − D2 (kπ)2 c
D1D1(17)

This gives a good approximation to the numerical results for Γ0 close to the threshold. Since the experimental
profiles of the shear flow and fluctuations close the shear
layer have few oscillations, we will concentrate in valuesof the parameter space such as the most unstable modes
are k = 1 or 2. In Fig. 2 we compare the analytical andnumerical results for D2 = 1.5 × 10−3 and k = 1. The flux
threshold for these parameters is 1.358. The analytical results are obtained from Eq. (16), and the numerical results
are obtained by advancing numerically Eqs. (1) to (3) until
a stationary solution is reached.

4

External torqueTo study the effect of an external torque, we add a term
τ = τ0 cos(kπx), to the r.h.s. of Eq. (2), so the Equation for

22.5

Γ03

3.5

4

4.5

Fig. 2 Comparison of analytical (solid line) and numerical(dots) values of the shear flow for D2 = 1.5 × 10−3 andk = 1.

the shear flow is now∂σ∂2 = τ0 cos(kπx) − σ − α3 2 N −4/3 E 2 σ∂t∂x" ∂2 σ #∂2 −5/3 2− 2 D2 N
E + D3
∂x
∂x2

(18)

In the rest of the calculations of this paper, we apply
the torque during a time t = 10, and then we remove thetorque to analyze the decay of the shear flow. The evolution of the integral of σ2 for different values of τ0 whenΓ0 = 1 is shown in Fig. 3. For this scan, D2 = 1.5 × 10−3 ,and k = 1, so we are below the threshold flux in absenceof external torque (subcritical regime). In most of thecases the shear flow has two decay scales and the change
between them is more pronounced as τ0 increases. The
square root of the integral of σ2 decays like eγ1 t just after removing the external torque (first decay region), andlike eγ2 t at larger times (second decay region). The first
decay rate is easily understood from Eq. (18). As we remove the torque, the instantaneous exponential decay rate
will be γ1 = −τ0 /σ1 , where σ1 corresponds to the stationary state with external torque τ0 . The second exponential
decay rate is very similar for all the cases, with γ2 -values
between −0.16 and −0.18.
The evolution of the integral of σ2 when Γ0 = 1.35
is shown in Fig. 4. This value of Γ0 is very close to thecritical value, Γc = 1.35794. For this scan, the evolutionof the integral of σ2 when we suppress the external torqueis no longer exponential, and can be fitted to the algeabricexpression C/(1 + t/T ), where C is a constant (value ofthe integral when we suppress the external torque), and wetake as origin of t the time when we remove the externaltorque. This is due to the fact that we are approaching thecritical point.
We have tried also scans in τ0 with values of Γ0 above
the threshold (supercritical regime), in particular Γ0 = 1.37
and 1.5. For these cases, the evolution of the integral of
σ2 when we suppress the external torque can be fitted toProceedings of ITC/ISHW2007
101102
τ0 = 0.5100

τ0 = 1τ0 = 2101

10-1

τ0 = 310-2

τ0 = 4
τ0 = 5100

10-310-410

Γ0 = 1.510-6

Γ0 = 0.8

Γ0 = 2

Γ0 = 1Γ0 = 2.5

Γ0 = 1.2

Γ0 = 32410051015

Γ0 = 1.37

10-710-2

Γ0 = 1.35

Γ0 = 0.4
Γ0 = 0.6-5-1

Γ0 = 0.20206810tt

Fig. 3 Evolution of the integral of σ2 for different values of τ0when Γ0 = 1. The external torque is suppressed at t = 10.

Fig. 5 Evolution of the integral of σ2 for different values of Γ0when τ0 = 0.5. The origin of t is taken at the time thatthe external torque is suppressed.

100
τ0 = 2

τ0 = 0.2

τ0 = 3

τ0 = 0.3

τ0 = 4

τ0 = 0.5

τ0 = 5

τ0 = 0.7

τ0 = 7

τ0 = 1τ0 = 10

10

32.5

T

τ0 = 0.1

2

1
1.5

1

0.1

00

2

468

10

0.5

1

1.5

2

2.5

3

3.5

Γ0t

Fig. 4 Evolution of the integral of σ2 for different values of τ0when Γ0 = 1.35. The origin of t is taken at the time thatthe external torque is suppressed.

C/(1 + t/T ) + C0 , where C0 is a constant with a value close
to the integral of σ2 for the stationary state in absence ofexternal torque.
Finally, the results for a scan in Γ0 are shown in Fig. 5.
In the y-axis, we represent the integral of σ2 when we suppress the external torque minus the value for the stationarystate in absence of external torque. The slower decay corresponds to the flux value closer to Γc , Γ0 = 1.35. The
decay time tends to infinity as Γ0 tends to Γc , which corresponds to the singularity of the transition point. The dataof the evolution of the squared root of the integral of σ2 foreach Γ0 -value is fitted to C exp(−t/T ) + C0 , as it is done inbiasing experiments in TJ-II [9]. The results for the decay
time T are shown in Fig. 6. The slower decays correspondto flux values closer to Γc . The damping is close to the viscous damping (one in our units) only when the flux is far
above the threshold. The results are similar when the valueof the external torque is changed.

Fig. 6 Results of the exponential fit of the decay of the squaredroot of the integral of σ2 for different values of Γ0 whenτ0 = 2.

[1] K. H. Burrell, Phys. Plasmas 4, 1499 (1997).
[2] P. W. Terry, Rev. Mod. Phys. 72, 109 (2000).
[3] B. Gonçalves, C. Hidalgo, M. A. Pedrosa, R. O. Orozco, E.
Sánchez and C. Silva, Phys. Rev. Lett. 96, 145001 (2006).
[4] Y. Koide, M. Kikuchi, M. Mori et al., Phys. Rev. Lett. 72,3662 (1993).
[5] K. H. Burrell, M. E. Austin, D. P. Brennan et al., Phys. Plasmas 8, 2153 (2001).
[6] D. del-Castillo-Negrete, B. A. Carreras. Physics of Plasmas
9, 118 (2002).
[7] B. A. Carreras, L. Garcia, M. A. Pedrosa, C. Hidalgo.
Physics of Plasmas 13, 122509 (2006).
[8] P. H. Diamond, S.-I. Itoh, K. Itoh, T. S. Hahm, Plasma Phys.
Controlled Fusion 47, R35 (2005).
[9] M.A. Pedrosa et al., Plasma Phys. Controlled Fusion (2007).