Физика плазмы, 2021, T. 47, № 8, стр. 684-692

1. ВВЕДЕНИЕ

Среди многочисленных исследований в области термоядерного синтеза наиболее актуальной задачей является поиск ключевых параметров плазмы токамака, изменение которых позволяет перейти в режим улучшенного удержания, реализуя так называемый L–H-переход. Несмотря на то, L–H-переход впервые наблюдался [1] более 30-ти лет назад, окончательного теоретического объяснения он пока не получил. Знание операционной области параметров плазмы для эффективного доступа в зону Н-моды играет важную роль при проектировании токамака ITER. Учитывая сложную нелинейную природу указанного явления, изучение физических процессов его сопровождающих проводят с помощью анализа результатов численного моделирования турбулентной динамики плазмы на основе анализа решений двухжидкостных редуцированных МГД-уравнений Брагинского. Как правило, рассматривается узкая пристеночная область токамака, разделенная на две зоны: основную плазму (с замкнутыми силовыми линиями) и зону SOL-плазмы.

Так, например, результаты численного моделирования проведенные в указанных зонах с учетом нагрева ионов показали [2–4], что эволюция турбулентности проходит в двух режимах. В первом режиме доминирует крупномасштабная турбулентная конвекция и величина скорости ЕхВ‑потока сравнительно небольшой величины. Во втором режиме, при превышении порога в величине ионного нагрева, интенсивность турбулентности увеличивается, одновременно возрастает скорость перехода турбулентной энергии флуктуаций в кинетическую энергию среднего течения (ZF) ЕхВ-потока. Как результат, наступает динамическая бифуркация, которая приводит к падению уровня флуктуаций и уменьшению турбулентных радиальных потоков частиц и тепла. Указано, что ключевую роль в этом процессе играет турбулентная сила напряжений Рейнольдса и геодезическая работа сил давления ионов ${{p}_{i}} \cdot {\text{div}}{{{\mathbf{V}}}_{E}}$. Кроме того отмечено, для получения в расчетах L–H-перехода необходимо использовать модель обобщенного вихря, включающего помимо электрического дрейфа также скорость диамагнитного дрейфа ионов [3]. Определенные черты рассмотренного механизма появления L‒H-перехода наблюдались при исследовании решений нелинейных одномерных уравнений хищник–жертва [4, 5]. Показано, что ZF играют роль “хищника”, растущего за счет поглощения кинетической энергии (“жертва”) флуктуаций и, в конечном счете, инициирующих L–H-переход.

В работе [6] изучалась динамика L–H-перехода и образования транспортного барьера за счет нагрева ионов в рамках двухполевой {ϕ, p_i} 3D модели уравнений Брагинского. Показано, что в процессе генерации скорости полоидального вращения за счет турбулентной силы Рейнольдса и силы неоклассического трения образуется положительная обратная связь между скоростью ЕхВ-дрейфа и градиентом ионного давления, которая при превышении порога нагрева приводит к появлению низкочастотных квазипериодических колебаний флуктуаций (LCO limit-cycle oscillations, I-фаза) с последующим переходом в режим улучшенного удержания.

В работах [7, 8] исследовалась динамика плазмы в токамаке при наличии полоидально-неоднородных источников в уравнениях для плотности, импульса и энергии ионной компоненты. Показано, что при соответствующей полоидальной ассиметрии этих источников происходит эффективное увеличение геодезической силы, связанной с неоднородностью тороидального магнитного поля (сила Стрингера–Винзора F_SW [9]), под действием которой растет полоидальная скорость ионов. В работе [10] рассмотренный случай полоидально-неоднородной зависимости коэффициента диффузии $D = {{D}_{0}} \cdot [1 + \delta \cdot \cos (\theta )]$ в уравнении для плотности обнаружил появление бифуркации в величине полоидальной скорости, что можно трактовать как переход из L-ежима в H.

В данной работе изучается механизм появления L–H-перехода при наличии полоидально-неоднородного нагрева ионов в тороидальном пристеночном слое токамака. Оказалось, что в условиях такого нагрева основную роль в генерации полоидальной скорости вращения играют не турбулентная сила напряжений Рейнольдса F_RE, а геодезическая сила F_SW и неоклассическая сила F_NEO ~ (T_i0)^5/2 возникающая за счет продольной вязкости, величина которой существенно растет в процессе нагрева ионов. Показано, что при достижении определенного порога в величине нагрева ионов, за счет роста средних скоростей диамагнитного вращения ионов $cp_{i}^{'}{\text{/}}enB$ и шировой скорости дрейфа ${{cE} \mathord{\left/ {\vphantom {{cE} B}} \right. \kern-0em} B}$ происходит подавление турбулентных флуктуаций и, как следствие, наблюдается L–H-переход.

2. ОСНОВНЫЕ УРАВНЕНИЯ

В работе используется 3-х полевая $\left\{ {\phi ,{{p}_{e}},{{p}_{i}}} \right\}$ система редуцированных двухжидкостных нелинейных МГД-уравнений Брагинского, описывающих поведение плазмы в условиях сильной столкновительности в тороидальном пристеночном слое токамака. Основные уравнения приведены в работах [11–14]. Расчеты показали, что в исследуемой низкотемпературной области токамака введенных переменных достаточно для адекватного описания взаимодействия дрейфовой турбулентности с зональными потоками и изучения L–H-перехода при нагреве ионов.

Указанную систему удобно привести к безразмерному виду, используя следующее преобразование переменных: $\frac{n}{{{{n}_{*}}}}$, $\frac{p}{{{{p}_{*}}}}$, $\frac{\phi }{{{{\phi }_{*}}}}$, $\frac{t}{{{{t}_{*}}}}$, $\frac{r}{{{{L}_{*}}}}$, где ${{n}_{*}} = {{n}_{{13}}}\frac{{{{L}_{*}}}}{d}$, ${{p}_{*}} = {{n}_{*}}{{T}_{*}}\frac{{{{L}_{*}}}}{d}$, ${{\phi }_{*}} = \frac{{{{B}_{0}}}}{c}\frac{{L_{*}^{2}}}{{{{t}_{*}}}}$, ${{t}_{*}} = \gamma _{B}^{{ - 1}} = C_{S}^{{ - 1}}\sqrt {\frac{{{{R}_{0}}d}}{2}} $ обратный инкремент ${{\gamma }_{B}}$ идеальной баллонной моды.

Для коэффициентов поперечного переноса по соображениям численной устойчивости были выбраны величины ${\text{[}}{{\nu }_{{ \bot ,}}},{{\chi }_{{e,i}}}_{ \bot }] \cong 0.01{\kern 1pt} - {\kern 1pt} 0.04{{D}_{{BOHM}}}$.

После процедуры нормализации МГД-уравнения преобразуются к виду

Величина плотности n₀(r) считается заданной функцией. Значения величин ${{Z}_{{{\text{WB}}}}},$ ${{Z}_{{eB}}},$ ${{Z}_{{iB}}}$ указаны ниже. В работе использована тороидальная система координат $({\text{r,}}\theta {\text{,}}\varsigma )$, в котором равновесное магнитное поле имеет вид ${\mathbf{B}} = \frac{{{{B}_{0}}{{R}_{0}}}}{R}{{{\mathbf{e}}}_{\varsigma }} + \frac{{\varepsilon {{B}_{0}}}}{{q(r)}}{{{\mathbf{e}}}_{\theta }}$, $R = {{R}_{0}} + r \cdot \cos \theta $_, где r – малый и R – большой радиусы токамака, q – коэффициент запаса устойчивости. Функции ${{\Psi }_{j}}$ в уравнениях (1a)–(1f) включают эффекты, связанные с кривизной магнитного поля токамака, где оператор кривизны C( f ) определяется как $С(f{\text{)}} = - \frac{{{{R}_{0}}{{B}_{0}}}}{2}{\text{rot}}\frac{b}{B} \cdot \nabla f$ = = ${\text{sin}}\theta \frac{{\partial {\text{f}}}}{{\partial {\text{r}}}} + \cos \theta \frac{{\partial f}}{{r\partial \theta }}$.

В системе (1а)–(1f) использованы следующие обозначения $\frac{D}{{Dt}} = \frac{\partial }{{\partial t}} + \left\{ \phi \right\}$, $\left\{ {A,B} \right\} = \frac{{\partial A}}{{\partial r}}\frac{1}{r}\frac{{\partial B}}{{\partial \theta }} - $ $ - \;\frac{{\partial B}}{{\partial r}}\frac{1}{r}\frac{{\partial A}}{{\partial \theta }}$, ${{\nabla }_{{||}}} = \frac{\partial }{{\partial \theta }}$, $\alpha = \frac{{{{\rho }_{S}} \cdot {{C}_{S}} \cdot {{t}_{*}}}}{{d \cdot {{L}_{*}}}}$, $\xi = \alpha \frac{d}{{{{L}_{*}}}}$, $\sigma = \frac{{g \cdot \alpha }}{\nu }$.

Нормализованные коэффициенты продольной диссипации для электронов и ионов имеют вид

Безразмерное выражение для продольной вязкости ${{\Pi }_{{||}}}$ [10] принимает вид:

Расчеты проводились в прямоугольной пристеночной области r₀ < r < a, –π < θ < +π, d = a – r₀ – ширина расчетного слоя. Значения нормировочных величин для плотности, электронной температуры и магнитного поля принимались равными ${{n}_{{13}}} = {{10}^{{13}}}$ см^–3, ${{T}_{*}} = 100$ эВ, ${{B}_{0}} = $ = 2 Тл.

Были выбраны размеры пристеночной области r₀ < r < a, r₀ = 24 см, a = 30 см при учете SOL-слоя. Незамкнутые силовые линии магнитного поля (SOL-слой) занимали область r_SOL < r < a, r_SOL = 29 см. Таким образом, ширина расчетной области равнялась d = a – r₀ = 6 см, а ширина SOL-слоя d_SOL = a – r_SOL = 1 см. Диссипация в SOL-слое связана с продольным движением ионов и электронов вдоль незамкнутых силовых линий магнитного поля (sheath current),которые соединяют поверхности пластин дивертора. Указанный ток описывают источниковые члены ${{\Lambda }_{j}}$ в области SOL-слоя и имеют вид [14, 15]:

Так как вне области SOL ${{\Lambda }_{j}}$ отсутствуют, в работе эти члены умножались на ступенчатую функцию вида:

Для численного расчета турбулентной динамики преобразуем систему (1a–1f), используя подход, развитый в работе [16], основанный на введение “медленных” и “быстрых” полевых переменных. В этом случае в тороидальной системе координат {r, θ, ς} полевые переменные могут быть представлены в виде суммы “медленной” аксиальносимметричной (не зависящей от угла ς) переменной f_A(r, θ, t) и “быстрой” баллонной переменной ${{f}_{B}}(r,\theta ,\varsigma ,t)$ имеющей зависимость от тороидального угла ς. Рассмотрение существенно упрощается, если изучать только одну тороидальную баллонную моду вблизи резонанса m = nq. Тогда для полевых переменных описывающих плазменные поля справедливо представление

После замены переменных $r \to r{\kern 1pt} '$, $\theta \to \theta {\kern 1pt} '$, $\varsigma \to \lambda $ и усреднения по углу λ полная система МГД‑уравнений разбивается на две взаимодействующие подсистемы уравнений для аксиальносимметричных переменных f_AX и баллонных f_S,C. При использовании представления плазменных полей в виде разложения (2) появляется возможность экономии временного ресурса при численных расчетах, за счет перехода от решения пространственно трехмерной {r, θ, ς} задачи к двумерной {r, θ}. При этом сохраняя физически важную трехмерную природу турбулентных флуктуаций. Кроме того, появляется возможность рассмотрения наиболее неустойчивых баллонных мод с $m \gg 1$.

Используемое разложение (2) полевых переменных на аксиально-симметричные и баллонные моды, позволяет разбить исходную систему МГД-уравнений (1a)–(1f) на две взаимодействующие между собой подсистемы мод. Система эволюционных уравнений для баллонных мод ${{f}_{{S,C}}}(r,\theta ,t)$, полученная из исходной системы (1a)–(1f) усреднением по углу λ в приближении $m \gg 1$ приведена в Приложении. В работе численное моделирование проводилось на основе совместного решения системы (1а)–(1f) для аксиально-симметричных мод и системы (Π.1)–(Π.3) для баллонных мод. Описание численной схемы, используемой для решения вышеприведенных уравнений подробно рассмотрена в [13].

3. ПОЛОИДАЛЬНОЕ ВРАЩЕНИЕ И ЭЛЕКТРИЧЕСКОЕ ПОЛЕ

Для удобства рассмотрения, выделим среднюю по полоидальному углу величину в каждой из рассматриваемых полевых переменных

Усредняя уравнение (1a) по угловым переменным получим эволюционное уравнение для средней величины вихря ${{w}_{0}}(r,t) = \left\langle {w(r,\theta ,t)} \right\rangle $, где по определению $w(r,\theta ,t) = {\text{div}}[{{{\mathbf{V}}}_{ \bot }} \times {\mathbf{b}}]$, ${{{\mathbf{V}}}_{ \bot }} = {{c[{\mathbf{E}} \times {\mathbf{b}}]} \mathord{\left/ {\vphantom {{c[{\mathbf{E}} \times {\mathbf{b}}]} B}} \right. \kern-0em} B} + $ $ + {{c[\nabla {{p}_{i}} \times {\mathbf{b}}]} \mathord{\left/ {\vphantom {{c[\nabla {{p}_{i}} \times {\mathbf{b}}]} {e{{n}_{0}}B}}} \right. \kern-0em} {e{{n}_{0}}B}}$:

где ${{\Pi }_{{RE}}} = {{\Pi }_{{AX}}} + {{\Pi }_{B}}$ – сумма турбулентных потоков импульса, связанных с аксиально-симметричными модами ${{\Pi }_{{{\text{AX}}}}} = \left\langle {{{{\tilde {V}}}_{{Er}}} \cdot \tilde {w}} \right\rangle $ и баллонными – ${{\Pi }_{B}} = \left\langle {m{\text{/}}2r \cdot ({{\phi }_{C}}{{w}_{S}} - {{\phi }_{S}}{{w}_{C}})} \right\rangle $.

Величина ${{w}_{0}}(r,t)$ связана с усредненными скоростями электрического дрейфа ${{V}_{E}}(r,t)$ и диамагнитного дрейфа ионов ${{V}_{{Di}}}(r,t)$ соотношениями:

где ${{V}_{0}}(r,t)$ – скорость полоидального вращения ионов. Далее с учетом малой ширины слоя $d \ll a$ считаем $\frac{1}{r}\frac{d}{{dr}}r(...) \approx \frac{d}{{dr}}(...{\text{)}}$.

Интегрирование (3) по радиусу r приводит к уравнению для полоидальной скорости ${{V}_{0}}(r,t)$

Здесь ${{F}_{{RE}}} = - ({{\Pi }_{{AX}}} + {{\Pi }_{B}})$ – турбулентная сила напряжений Рейнольдса, ${{F}_{{SW}}} = - {{\left\langle {({{p}_{e}} + {{p}_{i}}) \cdot \sin \theta } \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {({{p}_{e}} + {{p}_{i}}) \cdot \sin \theta } \right\rangle } {{{n}_{0}}}}} \right. \kern-0em} {{{n}_{0}}}}$ – геодезическая сила (сила Стрингера–Винзора SW), ${{F}_{{NEO}}} = - {{\nu }_{{||}}} \cdot ({{V}_{0}} - {{V}_{{NEO}}}){\text{/}}{{n}_{0}}$ – неоклассическая сила за счет продольной вязкости, ${{F}_{{SOL}}} = $ $ = \int_{{{r}_{0}}}^r {\left\langle {{{\Lambda }_{w}}} \right\rangle dr} $ – сила за счет продольного тока в диверторе, ${{\nu }_{ \bot }}$ – поперечная вязкость ионов.

Интегрируя равенство ${{\partial {{\phi }_{0}}} \mathord{\left/ {\vphantom {{\partial {{\phi }_{0}}} \partial }} \right. \kern-0em} \partial }r = {{V}_{0}} - {{V}_{{Di}}}$ по радиусу находим величину (безразмерного) электростатического потенциала

Граничное условие выбираем в виде ${{\phi }_{0}}(r = a) = $ $ = \xi \cdot {{T}_{{e0}}} \cdot {{\Lambda }_{e}}$.

Важно отметить, что электрическое поле в токамаке возникает вследствие электроиндукции E ~ VB/c при движении плазмы в магнитном поле и разделения зарядов за счет силы Лоренца и градиента давления, а не за счет появления объемного заряда из-за неамбиполярных потерь заряженных частиц ($\left\langle {{{j}_{r}}} \right\rangle = 0$). Хотя последний эффект может при определенных условиях возникать вблизи стенок камеры токамака. Из расчетов следует, что в обычном омическом режиме силы F_RE > 0 и F_SW < 0 направлены в разные стороны, компенсируя друг друга. Однако полной компенсации не происходит, и суммарная сила, как правило, положительна, что приводит к полоидальному вращению ионов V₀ > 0 в сторону диамагнитного дрейфа электронов.

4. ЧИСЛЕННЫЕ РАСЧЕТЫ ДИНАМИКИ L–H-ПЕРЕХОДА

Численные расчеты проводились в пристеночной области r₀ = 24 cм < r < a = 30 cм при следующих параметрах токамака Т-10: R = 150 cм, d = 6 cм, B = 2.1T, q = 4.5, n₁ = n₀(r₀) = 0.3, n₂ = = n₀(а) = 0.2 (плотность в единицах 10¹³ см³), T_be(r₀) = 60 эВ, T_bi(r₀) = 30 эВ, T_we(a) = T_wi(a) = = 10 эВ. Рабочий газ-дейтерий m_i = 2, Z_eff = 1.6. Профиль плотности выбран в виде: n₀(r) = n₁ + + (n₂ – n₁)*(1 – x²)² , 0 < x = (r – r₀)/d < 1. Величина номера полоидальной гармоники баллонной моды m = 30 выбиралась из условия максимальности величины турбулентного теплового потока электронов. В турбулентном режиме рассматриваемые полевые переменные осциллируют во времени, поэтому при построении графиков использовалась формула усреднения:

В работе источник нагрева ионов S_pi(r, θ) был выбран в виде следующей зависящей от координат x, θ, t функции: S_pi(r, θ, t) = s₀f(t)(1 – x²)²g(θ), x = (r – r₀)/d. Считалось, что нагрев включается в момент времени t = t₀. Для функции f(t) использовались формулы в виде; f(t) = 0 при 0 < t < t₀, и f(t) = 1 – exp[–(t – t₀)/Δt] при t > t₀. Полоидально-неоднородная зависимость источника нагрева учитывалась с помощью функции g(θ) = exp[–y²], где y = (θ – θ₀)/Δ. В расчетах Δ = π/3.

Из эволюционного уравнения (1c) для p_i(r, θ, t) нетрудно получить уравнение для синус компоненты $\left\langle {{{p}_{i}}\sin \theta } \right\rangle (r,t)$, являющейся геодезической акустической (ГА) модой давления ионов. Очевидно, что при достаточно высоких уровнях S_pi знак величины$\left\langle {{{p}_{i}}\sin \theta } \right\rangle (r,t)$, определяется знаком величины $\left\langle {{{S}_{{pi}}}\sin \theta } \right\rangle $ – синус компоненты источника нагрева

Ясно, что в результате нагрева указанной моды будет происходить одновременное воздействие на величину силы ${{F}_{{SW}}}\sim \left\langle {{{p}_{i}}\sin \theta } \right\rangle $, а значит и на величину полоидальной скорости V₀. Нетрудно убедиться, что при расположении источника тепла в верхней полуплоскости камеры токамака (в области θ₀ ≈ +π/2) для данной гауссообразной функции g(θ) величина $\left\langle {{{S}_{{pi}}}\sin \theta } \right\rangle $ принимает положительное значение. И наоборот, при расположении этого источника в нижней полуплоскости (в области θ₀ ≈ –π/2) величина $\left\langle {{{S}_{{pi}}}\sin \theta } \right\rangle $ отрицательна. Соответственно в первом случае, как следует из уравнения (7), при достаточно высоком уровне нагрева величина ГА моды ионов величина $\left\langle {{{p}_{i}}\sin \theta } \right\rangle $ становится положительной (F_SW < 0), во втором – отрицательной (F_SW > 0). Рисунок 1 показывает поведение величин $\{ \left\langle {{{p}_{{e,i}}} \cdot \sin \theta } \right\rangle \} $, где $\{ A\} = \frac{1}{{{{\Delta }_{r}}}}\int_{{{r}_{0}}}^a {A(r)} rdr$, ${{\Delta }_{r}} = ({{a}^{2}} - r_{0}^{2}){\text{/}}2$, в течение времени импульса нагрева для трех различных случаев: отсутствия нагрева (s₀ = 0) и нагрев (s₀ = 4.0) для θ₀ = –π/2 и θ₀ = +π/2. Видно, что при θ₀ = = +π/2 усредненная по слою величина ГА моды давления ионов $\{ \left\langle {{{p}_{i}} \cdot \sin \theta } \right\rangle \} $ при выбранном уровне нагрева s₀ = 4.0 принимает положительное значение, а при θ₀ = –π/2 отрицательное, соответственно определяя знак силы F_SW.

Рис. 1.

Временная зависимость средних величин ГА мод давления электронов $\left\{ {\left\langle {{{p}_{e}}\sin \theta } \right\rangle } \right\}$ и ионов $\left\{ {\left\langle {{{p}_{i}}\sin \theta } \right\rangle } \right\}$ для случаев a) s₀ = 0, б) s₀ = 4, θ₀ = –π/2, в) s₀ = 4, θ₀ = +π/2.

На рис. 2 изображены радиальная зависимость основных сил F_RE, F_SW и F_NEO, баланс которых, определяет величину полоидальной скорости V₀ для указанных выше трех различных случаев. Видно, что в случае отсутствия нагрева s₀ = 0 достижение квазистационарного состояния турбулентных колебаний V₀ осуществляется за счет баланса сил F_RE > 0 и F_SW < 0, сила F_NEO – мала и в балансе практически не участвует. При нагреве (s₀ = 4.0) в двух других случаях θ₀ = –π/2 и θ₀ = = +π/2 геодезические силы F_SW существенно возрастают и согласно вышесказанному относительно знака $\left\langle {{{S}_{{pi}}}\sin \theta } \right\rangle $, имеют противоположные знаки. Сильная зависимость неоклассической силы от температуры ионов F_NEO ~ (T_i0)^5/2 приводит в этих двух случаях (нагрева ионов) к ее существенному росту. В первом случае θ₀ = –π/2, V₀ > 0, F_SW > 0 при этом сила F_NEO < 0. Во втором случае θ₀ = +π/2 V₀ < 0 F_SW < 0 и сила F_NEO > 0. В обоих случаях сила F_NEO исполняет функцию торможения полоидального вращения. Таким образом, в каждом из рассмотренных случаев различного расположения источника нагрева S(r, θ, t) геодезические силы и силы продольной вязкости различным образом участвуют в балансе, компенсируя друг друга из-за противоположных знаков. При этом в этих случаях сила напряжений Рейнольдса знака не меняет F_RE > 0 и также участвует в балансе, оставаясь, однако, значительно меньше двух других по величине. Можно сделать вывод о том, что основную роль в генерации полоидальной скорости вращения (в рассматриваемых случаях полоидально-неоднородного нагрева) играет не турбулентная сила напряжений Рейнольдса F_RE как в случае полоидально-однородного нагрева [2–4], а силы F_SW и F_NEO.

Рис. 2.

Радиальные профили сил F_SW, F_RE, F_NEO для случаев a) s₀ = 0, б) s₀ = 4, θ₀ = –π/2, в) s₀ = 4, θ₀ = +π/2.

Рисунок 3 демонстрирует радиальную зависимость скоростей V₀(r), V_Di(r), V_E(r) на ширине расчетного слоя. Из расчетов следует, что величина V_E(r) средней шировой скорости ЕхВ-дрейфа в случае нагрева ионов существенно возрастает по сравнению со случаем отсутствия нагрева. Следствием этого при нагреве происходит подавление турбулентных флуктуаций и наблюдается переход плазмы в режим улучшенного удержания, именуемым L–H-переходом. Расчеты показывают, что при превышении определенного порога в величине нагрева резко уменьшаются амплитуды тепловых колебаний баллонных мод p_es и p_ec. Эффект L–H-перехода представлен на рис. 4, на котором изображено временное поведение тепловых потоков электронов для случаев (s₀ = 4.0) θ₀ = –π/2 и θ₀ = +π/2. Видно, что после включения импульса (t₀ = 250 мкс) в течение некоторого времени колебания продолжаются с прежней амплитудой, демонстрируя так называемую I-фазу. Указанные колебания отчетливо наблюдаются в экспериментах, предшествуя L–H-переходу. Из рис. 4 видно, что только через время t ~ 300 мксек для θ₀ = –π/2 и t ~ 450 мксек для θ₀ = +π/2 колебательный режим I-фазы исчезает и происходит переход в режим улучшенного удержания с характерным резким уменьшением теплового потока электронов $\left\{ {\left\langle {{{Q}_{e}}} \right\rangle } \right\}$. Очевидно, что при одной и той же величине теплового импульса L–H-переход в случае θ₀ = −π/2 реализуется эффективней. Это связано с тем, что в случае θ₀ = –π/2 полоидальная скорость на промежутке слоя положительна V₀ > 0, по сравнению с θ₀ = +π/2, когда V₀ < 0 на большем промежутке слоя (см. рис. 3). Как результат, в соответствии с равенством ${{V}_{E}} = {{V}_{0}} - {{V}_{{Di}}}$, ${{V}_{{Di}}} = \tau \cdot {{d{{p}_{{i0}}}} \mathord{\left/ {\vphantom {{d{{p}_{{i0}}}} {dr}}} \right. \kern-0em} {dr}} < 0$, средняя скорость V_E ~ ~ ЕхВ‑дрейфа в первом случае больше, чем во втором, что приводит к более эффективному подавлению турбулентных флуктуаций. Отметим, что важную роль в получении L–H-перехода играет использование модели обобщенного вихря, включающего помимо электрического дрейфа также скорость диамагнитного дрейфа ионов ${\text{w}} = \nabla _{ \bot }^{{\text{2}}}\phi + \tau \cdot \nabla _{ \bot }^{{\text{2}}}{{p}_{i}}$, $\tau = {\alpha \mathord{\left/ {\vphantom {\alpha {{{n}_{0}}}}} \right. \kern-0em} {{{n}_{0}}}}$ [3]. При α = 0 L–H-переход не наблюдался. На рис. 5 изображена временная зависимость величины электронного потока в случае изменения величины параметра α в момент времени t = 800 мксек. Видно, что уменьшение величины α в 3 раза приводит к обратному переходу из H-режима в L-режим. В заключение заметим, что рассмотрение явления L‒H-перехода в рамках более адекватной 4-х полевой модели $\left\{ {\phi ,\;n{\text{,}}\;{{p}_{e}},\;{{p}_{i}}} \right\}$ c учетом турбулентных флуктуаций плотности практически не меняет рассмотренных выше результатов.

Рис. 3.

Радиальные профили скоростей V₀, V_Di, V_E для случаев a) s₀ = 0, б) s₀ = 4, θ₀ = –π/2, в) s₀ = 4, θ₀ = +π/2.

Рис. 4.

Временная зависимость средних величин теплового потока электронов для случаев a) s₀ = 4, θ₀ = –π/2, б) s₀ = 4, θ₀ = +π/2.

Рис. 5.

Временная зависимость средней величины теплового потока электронов для случая уменьшения параметра α → α/3 (s₀ = 4, θ₀ = –π/2).

5. ЗАКЛЮЧЕНИЕ

В работе с помощью численного моделирования решений редуцированных двухжидкостных уравнений Брагинского изучено воздействие полоидально-неоднородного нагрева ионов на турбулентную динамику пристеночной плазмы токамака. Показано, что при достижении определенного порога в величине нагрева ионов происходит динамическая бифуркация в решениях рассматриваемой нелинейной МГД-системы, что приводит к росту величины скорости течения – среднего ЕхВ-дрейфа. При этом реализуется сценарий “хищник-жертва”, в котором указанное течение играет роль “хищника”, поглощающего “жертву” – энергию турбулентных флуктуаций. В результате происходит подавление турбулентных флуктуаций и, как следствие, наблюдается L‒H‑переход. Из расчетов следует, что в условиях несимметричного нагрева основную роль в генерации полоидальной скорости вращения играют не турбулентная сила напряжений Рейнольдса F_RE, а геодезическая сила F_SW и неоклассическая сила, возникающая за счет продольной вязкости F_NEO ~ (T_i0)^5/2, величина которой существенно растет в процессе нагрева ионов. Следует отметить, что сила F_SW ~ (p_emn + p_imn) n = 0, m = 1 по сути является ГА-модой суммарного давления. Таким образом можно утверждать, что указанная мода играет ключевую роль в рассматриваемом явлении L–H-перехода. Проведенные расчеты для различной локализации источника тепла (в верхней или нижней полуплоскости камеры токамака), показали, что найденное различие приводит к различию в эффективности подавления турбулентного потока тепла электронов и временной длины предшествующих L–H-переходу квазипериодических колебаний (I-фаза).