Физика плазмы, 2019, T. 45, № 12, стр. 1072-1087

1. ВВЕДЕНИЕ

Процессы переноса частиц и энергии остаются одной из ключевых проблем в исследованиях по магнитному удержанию плазмы, направленных на создание эффективного термоядерного реактора. Поперечные потоки тепла и частиц, наблюдаемые в токамаках и других системах магнитного удержания высокотемпературной плазмы, как правило, значительно превышают соответствующие неоклассические потоки и, потому, называются аномальными и ассоциируются с тем или иным видом турбулентности плазмы. К настоящему времени выполнено большое число работ, посвященных анализу турбулентности плазмы и построению транспортных моделей и кодов для описания и моделирования аномальных процессов переноса в замагниченной плазме. Сравнительный анализ таких работ приведен, в частности, в обзорах [1–5].

Многие известные в настоящее время транспортные коды позволяют описывать глобальные нестационарные транспортные процессы в различных областях плазмы токамака с учетом реальной формы магнитных поверхностей. При этом, благодаря использованию потоковых координат, эти коды, по-существу, остаются одномерными по радиальной потоковой координате, а вся специфика магнитной конфигурации учитывается набором метрических коэффициентов. К этому классу кодов относятся такие широко известные коды как ASTRA [6], TRANSP [7], TOPICS-IB [8], CRONOS [9], JINTRAC [10] и другие.

Однако, для описания потоков частиц и тепла в таких кодах, как правило, используется диффузионное приближение с локальными эмпирическими транспортными коэффициентами. В то же время, значительное число экспериментов, выполненных на токамаках (см., например, [11–13]) и на других системах магнитного удержания высокотемпературной плазмы (стеллараторах [14, 15], тандемных ловушках [16, 17] и т.п.), показывают, что диффузионное приближение с локальными транспортными коэффициентами не позволяет достаточно адекватно описать аномальные процессы переноса, вызываемые низкочастотной (НЧ) турбулентностью. Для изучения НЧ-турбулентности и результирующих аномальных транспортных процессов в замагниченной плазме разумно использовать прямое компьютерное моделирование нелинейной динамики плазмы. При этом важно не просто развивать коды, адекватно описывающие турбулентность, но и включать рассчитываемые с их помощью турбулентные потоки частиц и тепла в самосогласованное описание глобальных транспортных процессов, как, например, это предполагается делать в работах [18, 19].

В отличие от большого числа работ, в которых моделирование НЧ-турбулентности основано на гирокинетическом подходе (см., например, обзоры [1, 2]), в наших работах для описания нелинейной динамики плазмы используется значительно более простая адиабатически-редуцированная динамическая модель (ARD-модель), которая требует существенно более скромных компьютерных ресурсов и допускает моделирование транспортных процессов на временах, сравнимых или бóльших энергетического времени удержания плазмы. Первоначально, в работах [20, 21] была построена ARD-модель турбулентной конвекции плазмы в цилиндрическом квазиравновесном плазменном шнуре, удерживаемом в чисто полоидальном магнитном поле. Устойчивость плазмы в такой системе (“пинч без продольного поля”) была исследована Б.Б. Кадомцевым в работе [22], где было обнаружено существование равновесного состояния плазмы с гранично-устойчивым (MS) профилем давления, спадающим к периферии. В развитие концепции MS-равновесий, В.Д. Шафрановым была предложена модель MS-равновесного удержания плазмы в “компактном торе” (или в “конфигурациях с обращенным полем” (FRC)). Эта идея и расчеты MS-равновесных состояний в FRC были представлены в работе [23]. Более позднее численное моделирование эволюции турбулентной плазмы в рамках цилиндрической ARD-модели [20, 21, 24] также было поддержано В.Д. Шафрановым и выявило ряд процессов самоорганизации плазмы, таких как формирование самосогласованных профилей давления, L-H переходы и т.д., которые ранее наблюдались в экспериментах на токамаках.

Совокупность полученных результатов создала предпосылки для применения упомянутой выше простой турбулентной модели к расчетам недиффузионного транспорта в основной горячей области токамака. С этой целью, исходная ARD-модель была модифицирована и адаптирована к цилиндрической модели токамака, и на ее основе был разработан код CONTRA-C, предназначенный для моделирования НЧ-турбулентности и транспортных процессов в токамаках с большим аспектным отношением $A = R{\text{/}}a$, где R и a обозначают, соответственно, большой и малый радиусы плазмы, и почти круглым полоидальным сечением магнитных поверхностей. Моделирование проводилось для условий ряда реальных экспериментальных разрядов на токамаке Т-10 с аспектным отношением $A = 5$. Описание модифицированной ARD-модели и полученные результаты представлены в работах [25–28]. Расчеты продемонстрировали тенденцию к формированию и поддержанию самосогласованных профилей давления плазмы, близких к наблюдаемым в экспериментах, быстрый отклик на изменение мощности и условий нагрева плазмы, быстрый отклик на выбросы энергии из области $q < 1$, вызываемые пилообразными колебаниями, и ряд других специфических свойств, весьма похожих на те, что наблюдаются в экспериментах на установках токамак. Далее, было проведено более детальное изучение процессов формирования и поддержания самосогласованных профилей давления плазмы в рамках нашей относительно простой теоретической модели НЧ турбулентной конвекции плазмы в токамаке и сравнение результатов моделирования с экспериментальными результатами, полученными в различных разрядах на токамаке Т-10. Результаты этих исследований представлены в работах [29, 30].

Отметим, что, наряду с кодом CONTRA-C, на основе ARD-модели был также разработан турбулентный блок CONTRA-A, который работает совместно с известным транспортным кодом A-STRA и позволяет рассчитывать турбулентную конвекцию плазмы и результирующие транспортные процессы с учетом тороидальных эффектов в токамаках с произвольной формой полоидального сечения шнура. При этом адаптированная версия кода ASTRA выступает в качестве транспортной оболочки и служит для расчета эволюции квазиравновесных параметров плазмы, силового равновесия, метрических коэффициентов и источников тепла и частиц. Описание физической модели турбулентного блока CO-NTRA-A и модифицированной версии кода ASTRA дано в работе [31].

В любой транспортной модели наряду с уравнениями, описывающими процессы переноса частиц, энергии, а в некоторых моделях и импульса, важную роль играют граничные условия к этим уравнениям. В случае тороидальных систем (токамак, стелларатор и т.п.) внутренние граничные условия для транспортных уравнений должны соответствовать отсутствию потоков частиц и тепла на магнитную ось. Что касается внешних граничных условий, то это непростой и принципиально важный вопрос. Как показали результаты наших предшествующих работ, выбор внешних граничных условий существенно влияет на энергетическое время удержания плазмы в проводимых расчетах. Данная работа посвящена анализу связи внешних граничных условий с характеристиками удержания плазмы в SOL и более детальному изучению влияния этих граничных условий на потоки тепла из основной плазмы и на самосогласованную эволюцию турбулентной плазмы в целом. Проведено сравнение результатов моделирования эволюции плазмы в переходных режимах, вызванных изменением мощности дополнительного нагрева, с результатами конкретных экспериментов. С помощью кода CONTRA-C, хорошо зарекомендовавшего себя при выполнении работ [25–30], и с использованием специально подобранных внешних граничных условий, проведено численное моделирование эволюции турбулентной плазмы для трех реальных разрядов с включением ЭЦР-нагрева плазмы в токамаке Т-10. Результаты, полученные для разрядов с различной мощностью ЭЦР-нагрева, продемонстрировали, что предложенные нелинейные граничные условия приводят к вполне разумной и согласующейся с результатами экспериментов зависимости энергетического времени удержания плазмы от мощности ЭЦР-нагрева.

В разделе 2 кратко представлена используемая ARD-модель плазмы токамака и вводятся внешние нелинейные граничные условия, основанные на упрощенном скейлинге удержания плазмы в SOL. В разделе 3 обсуждаются результаты компьютерного моделирования эволюции турбулентной плазмы для трех разрядов в Т-10 с разными уровнями вводимой ЭЦР-мощности и с учетом различных процессов, влияющих на время удержания плазмы. Раздел 4 кратко суммирует основные результаты.

2. РЕДУЦИРОВАННАЯ ДИНАМИЧЕСКАЯ МОДЕЛЬ ТУРБУЛЕНТНОЙ КОНВЕКЦИИ ПЛАЗМЫ ТОКАМАКА И УСОВЕРШЕНСТВОВАННЫЕ ВНЕШНИЕ ГРАНИЧНЫЕ УСЛОВИЯ

Как отмечалось во Введении, для описания процессов переноса тепла и частиц в токамаке мы используем относительно простую ARD-модель турбулентной конвекции плазмы, основанную на системе нелинейных адиабатически-редуцированных уравнений магнито-гидродинамического типа. Основные принципы используемой модели были развиты в работах [20–30] сначала применительно к системам с чисто полоидальным магнитным полем, а затем модель была модифицирована и адаптирована к цилиндрической модели токамака, и на ее основе был разработан код CONTRA-C, предназначенный для моделирования НЧ-турбулентности и транспортных процессов в токамаках с большим аспектным отношением и почти круглым полоидальным сечением магнитных поверхностей, к числу которых относится и токамак Т-10. Результаты моделирования для ряда реальных экспериментальных разрядов в Т‑10 представлены в работах [25–30]. В данной работе мы также ориентируемся на условия экспериментов на токамаке Т-10. Поэтому для моделирования снова использовался хорошо зарекомендовавший себя код CONTRA-C.

Достаточно подробное описание и обсуждение физической модели, заложенной в код CONTRA-C, содержится в работах [26, 28, 30]. Поэтому здесь лишь кратко напомним основные положения этой модели. В цилиндрической модели токамака рассматривается плазменный столб, удерживаемый в радиально-неоднородном магнитном поле, образующим цилиндрические магнитные поверхности. Предполагается, что плазменный столб имеет длину $2\pi R$ и отождествленные торцы. Это позволяет ввести эффективный “тороидальный”, угол $\varphi = z{\text{/}}R$, где z – координата, направленная вдоль оси плазменного столба, и считать все физические величины, включая флуктуации, периодическими функциями φ. Кроме этого, предполагается, что плазма самосогласованно поддерживается вблизи турбулентно-релаксированного (TR) состояния, которое гранично-устойчиво относительно идеальной перестановочной моды, инициируемой градиентом давления плазмы. Как и в работах [22, 23], TR-состояние соответствует пространственной однородности функции $S = p{{U}^{\gamma }}$, где $p = n({{T}_{e}} + {{T}_{i}})$ – полное давление плазмы, n – плотность частиц плазмы, ${{T}_{{e,i}}}$ – температуры электронов и ионов, $U(\psi ) = dV(\psi ){\text{/}}d\psi $ = $\oint {d{{l}_{p}}{\text{/}}{{B}_{p}}} $ – удельный объем силовой трубки, ψ – полоидальный магнитный поток. Сжимаемость плазмы характеризуется показателем адиабаты γ. При моделировании токамаков разумно выбрать эффективное значение $\gamma = 2$. В этом случае результаты моделирования лучше согласуются со многими экспериментами на токамаках, в частности, описанными в работах [34–37], а также с концепцией “канонических профилей” [38–41], применяемой для интерпретации транспортных процессов в токамаках. Полоидальное магнитное поле можно представить в виде ${{{\mathbf{B}}}_{p}} = [\nabla \psi \times \nabla \varphi ]{\text{/}}2\pi $. Функция S является однозначной функцией энтропии плазмы, заключенной в удельном объеме U. Поэтому рассматриваемое TR-состояние можно назвать состоянием турбулентного равнораспределения (equipartition) энтропийной функции S (см. [42, 43]).

Процедура адиабатического редуцирования уравнений движения базируется на предположении, что относительная величина флуктуаций давлений электронов и ионов в процессе развития турбулентной конвекции остается малой как ${{\varepsilon }^{2}}$, где параметр адиабатичности $\varepsilon \ll 1$ связан с обратным числом Пекле следующим соотношением: ${{\varepsilon }^{3}} \sim {{\chi }_{{bg}}}{\text{/}}a{{c}_{s}}$. Здесь ${{\chi }_{{bg}}}$ – фоновая (ионная неоклассическая) локальная температуропроводность, ${{c}_{s}}$ – скорость звука. Наличие малого параметра адиабатичности ε позволяет ограничиться рассмотрением НЧ-конвекции плазмы, описываемой адиабатическим полем скоростей ${{{\mathbf{v}}}_{a}} \sim \varepsilon {{c}_{s}}$, которое имеет вид, задаваемый формулой (1) работы [30]. Такое адиабатическое поле скоростей не возбуждает “быстрые” (высокочастотные) устойчивые магнитозвуковые и альфвеновские волны, а также выравнивает плотность и давления электронов и ионов по объему силовой трубки, формируя вытянутые вдоль магнитного поля движущиеся плазменные филаменты. Величина и пространственно-временная структура ${{{\mathbf{v}}}_{a}}$ определяются 2D-функцией ϕ, которая имеет смысл электрического потенциала и зависит от эффективного “тороидального” угла φ и эффективного (потокового) малого радиуса $\rho = $ $ = \sqrt {\Psi (\psi ,t){\text{/}}\pi {{B}_{0}}(t)} $, введенного в описании транспортного кода ASTRA [6], где Ψ – тороидальный магнитный поток, а ${{B}_{0}}(t)$ – вакуумное тороидальное магнитное поле в центре камеры токамака. При переходе от полоидального потока ψ, используемого в качестве радиальной потоковой координаты в более ранних работах [20–30], к эффективному малому радиусу ρ в уравнениях для флуктуаций появляется коэффициент $h(\rho ,t) = $ $ = 2\pi \rho {\text{/}}({{\partial }_{\rho }}\psi )$ = $q(\rho ,t){\text{/}}{{B}_{0}}$, где $q(\rho ,t)$ – традиционный коэффициент запаса устойчивости токамака, и удельный объем слоя толщиной $d\rho $: $V{\kern 1pt} ' \equiv {{\partial }_{\rho }}V = 2\pi \rho U{\text{/}}h$.

Адиабатически редуцированное уравнение движения плазмы в ARD-модели представлено уравнением (6) работы [30] и, по существу, является обобщением уравнения Эйлера для двумерных несжимаемых течений на случай наличия специфической сжимаемости замагниченной плазмы, обусловленной радиальной неоднородностью удельного объема U, а также наличию силы, связанной с градиентом давления плазмы и вызывающей радиальную конвекцию плазмы. Это уравнение имеет вид локального закона сохранения для величины $w$, названной в работах [26, 28, 30] динамической завихренностью и связанной с электрическим потенциалом следующим соотношением:

где ${{m}_{i}}$ – масса иона, c – скорость света.

Редуцированные уравнения для давлений электронов и ионов, а также плотности плазмы записываются в терминах более адекватных переменных: энтропийных функций ${{S}_{{(e,i)}}} = {{p}_{{(e,i)}}}{{U}^{2}}$ и числа частиц в объеме силовой трубки $D = nU$, введенных вместо давлений ${{p}_{{(e,i)}}}$ и плотности частиц n соответственно. Эти переменные позволяют явно учесть инвариантные свойства исходных МГД-уравнений. Обе эти функции состоят из средних по поверхности медленно меняющихся функций ${{\bar {S}}_{{(e,i)}}}(t,\rho )$, $\bar {D}(t,\rho )$ и малых флуктуирующих компонент ${{S}_{{(e,i)f}}}(t,\rho ,\varphi )$, ${{D}_{f}}(t,\rho ,\varphi )$:

При этом предполагается, что функции ${{\bar {S}}_{{(e,i)}}}(t,\rho )$ близки к MS-состоянию: $\left| {\nabla {{{\bar {S}}}_{{(e,i)}}}} \right| \sim {{\varepsilon }^{2}}{{\bar {S}}_{{(e,i)}}}{\text{/}}a$. Динамическая завихренность и электрический потенциал также представляются в виде сумм средних и флуктуирующих компонент: $w = \bar {w}(t,\rho ) + {{w}_{f}}(t,\rho ,\varphi )$ и $\bar {\phi }(t,\rho ) + {{\phi }_{f}}(t,\rho ,\varphi )$.

Механизм аномального переноса тепла основан на конкуренции, когда нагрев плазмы, фоновая (неоклассическая) теплопроводность и внешние граничные условия искажают начальный профиль давления, делая его слабонеустойчивым, а неустойчивость возбуждает и поддерживает квазидвумерную нелинейную конвекцию, которая стремится восстановить TR-состояние с MS-профилем давления и приводит к аномальному недиффузионному переносу тепла и плотности плазмы. Строго говоря, следует рассматривать отдельные уравнения переноса тепла как для электронной, так и для ионной компонент плазмы, связанные между собой лишь слагаемыми, описывающими столкновительный теплообмен между компонентами. Однако основной задачей данной работы является исследование влияния внешних граничных условий на интегральное энергетическое время удержания плазмы ${{\tau }_{E}}$. Поэтому, для упрощения расчетов, мы проводим моделирование, считая, что отношение температур ${{T}_{i}}{\text{/}}{{T}_{e}}$ задано и фиксировано. Это не оказывает существенного влияния на величину ${{\tau }_{E}}$, однако позволяет ограничиться суммарным уравнением для объединенной энтропийной функции $S = {{S}_{e}} + {{S}_{i}}$.

В соответствии с работами [26, 28, 30], уравнение переноса плотности получается из уравнения непрерывности путем перехода к потоковым координатам, интегрирования по удельному объему силовой трубки и усреднения по магнитной поверхности:

Здесь ${{\bar {Q}}_{n}}$ – поверхностно-усредненная плотность источника частиц, а полный поток частиц $\Gamma $ представлен в виде суммы “фонового” потока ${{\Gamma }_{{bg}}}(\rho ,t)$ и усредненного по магнитной поверхности турбулентно-конвективного потока ${{\Gamma }_{{turb}}}(\rho ,t)$, который имеет вид

где флуктуации ${{D}_{f}}$ рассчитываются с помощью уравнений (12) работы [30], а флуктуации потенциала восстанавливаются по функции w путем решения уравнения (1). В качестве “фонового” потока ${{\Gamma }_{{bg}}}(\rho ,t)$ может рассматриваться неоклассический диффузионный поток, а также какой-либо иной поток диффузионного типа, полученный из модельных соображений.

В отличие от уравнения непрерывности для плотности, в исходных уравнениях, описывающих перенос тепла и флуктуации давления, конвективные члены, содержащие ${{{\mathbf{v}}}_{a}}$ и задающие аномальный характер процессов переноса тепла, нельзя представить в форме полной дивергенции потоков. Известно, что полностью дивергентный вид имело бы уравнение, соответствующее закону сохранения полной энергии, которое представляет собой сумму тепловых энергий электронов и ионов, кинетической энергии плазмы и энергии магнитного поля. Однако использовать такое уравнение для полной энергии в численных расчетах довольно сложно и неудобно. В отсутствии турбулентной конвекции, кинетическая энергия, связанная со скоростью “фонового” потока частиц, пренебрежимо мала, что позволяет в традиционной версии кода ASTRA записать отдельные уравнения для переноса тепла как в электронной, так и в ионной компонентах плазмы в дивергентной форме. В присутствии турбулентной конвекции, кинетической энергией, связанной с ${{{\mathbf{v}}}_{a}}$ и ее изменениями на временах масштаба времени жизни плазмы, пренебречь нельзя. В этом случае, в соответствии с результатами предшествующих работ [26–30], вместо уравнений для переноса тепла, предпочтительнее использовать уравнения переноса энтропийных функций ${{S}_{{(e,i)}}}$, в которых доминирующий член, описывающий турбулентную конвекцию, имеет дивергентный вид. В данной работе также использовалось уравнение для переноса объединенной энтропийной функции в форме

где ${{q}^{{turb}}}$ – суммарный турбулентно-конвективный поток тепла, выраженный через турбулентно-конвективный поток энтропии и определяемый выражением ${{q}^{{bg}}} = q_{e}^{{bg}} + q_{i}^{{bg}}$ – суммарный “фоновый” поток тепла c единым для ионов и электронов коэффициентом теплопроводности, примерно соответствующим неоклассической ионной теплопроводности, а $\bar {P}$ – поверхностно-усредненная полная плотность мощности источников нагрева, включающая мощность электрон-циклотронного резонансного нагрева (ЭЦР нагрева) ${{P}_{{ECR}}}$, мощность омического нагрева (OH) ${{P}_{{OH}}}$, мощность вязкого нагрева вследствие диссипации кинетической энергии турбулентных течений ${{P}_{{visc}}}$, определяемую выражением (16) работы [30], и потери на излучение ${{P}_{{rad}}}$: $\bar {P} = {{P}_{{ECR}}} + {{P}_{{OH}}} + {{P}_{{visc}}} - {{P}_{{rad}}}$. Для расчета флуктуаций ${{S}_{f}}$ в выражении (6) используется уравнение (13) работы [30]. Источник вязкого тепловыделения ${{P}_{{visc}}}$ возвращает кинетическую энергию флуктуаций, почерпнутую из тепловой энергии плазмы в результате работы сил давления, обратно в тепловую энергию плазмы, поддерживая закон сохранения полной энергии.

Перейдем теперь к обсуждению граничных условий. Как отмечалось во Введении, в тороидальных системах внутренние граничные условия для транспортных уравнений должны соответствовать отсутствию потоков частиц и тепла на магнитную ось. Это достаточно легко реализуется при учете того, что и турбулентные флуктуации также должны обращаться в нуль на магнитной оси. Наиболее корректными внешними граничными условиями для уравнений переноса (3) и (5) в основной горячей области были бы условия непрерывности потоков частиц и тепла на границе этой области со “слоем обдирки” (SOL). Однако, строго говоря, такой подход требует решения самосогласованной транспортной задачи в SOL, что представляет собой самостоятельную весьма сложную задачу. Поэтому, при построении внешних граничных условий для уравнений переноса, можно попытаться ограничиться некоторыми модельными предположениями относительно удержания плазмы в SOL.

В рамках традиционного подхода, основанного на диффузионной модели процессов переноса с заданными транспортными коэффициентами, используемого, в частности, в упомянутых выше кодах ASTRA, TRANSP, TOPICS-IB, CRONOS, JINTRAC, часто ограничиваются простейшими граничными условиями 1-го рода, соответствующими обращению температур электронов и ионов (а часто и плотности) в нуль на границе с областью SOL. В транспортной модели, основанной на уравнениях диффузионного типа, использование таких граничных условий оправдано тем, что температуры электронов и ионов в области SOL значительно ниже, чем в основной плазме, и учет конечной величины этих температур лишь незначительно влияет на величину потоков тепла из основной плазмы в SOL, которые в этом случае определяются коэффициентами теплопроводности и градиентами температур в основной плазме. Однако, при наличии в основной плазме доминирующих турбулентно-конвективных потоков тепла, такие упрощенные граничные условия становятся принципиально некорректными, поскольку величины этих доминирующих потоков могут изменяться в несколько раз при почти неизменных градиентах температур (например при L-H-переходах, сбросах энергии вследствие пилообразных колебаний и других явлениях) и существенно зависят от близости профиля суммарного давления плазмы к TR-состоянию, при котором температуры электронов и ионов имеют конечную величину на границе с SOL. В этом случае, при построении граничных условий, действительно необходимо учитывать непрерывность тепловых потоков на границе основной горячей плазмы с SOL, которую следует дополнить некоторой, возможно упрощенной, но физически оправданной моделью удержания плазмы в SOL.

Достаточно естественно предположить, что потери частиц и энергии из области SOL пропорциональны средним значениям плотностей частиц и тепловой энергии в SOL соответственно, а они, в свою очередь, пропорциональны плотностям частиц и тепловой энергии на границе основной горячей области и SOL (при $\rho = {{\rho }_{L}}$). Тогда условия непрерывности полных потоков частиц и тепла на внешней границе расчетной области могут быть записаны в форме обобщенных граничных условий 3-го рода

где коэффициенты ${{\nu }_{n}}$ и ${{\nu }_{E}}$, по своему физическому смыслу, обратно пропорциональны эффективным временам удержания частиц и энергии в SOL. В зависимости от модели эти коэффициенты могут быть выбраны в форме констант, заданных функций времени, а также зависеть от значений $\bar {n}$ и ${{\bar {T}}_{{(e,i)}}}$ на границе, в соответствии с тем или иным скэйлингом потерь частиц и энергии в SOL. Абсолютные значения коэффициентов ${{\nu }_{n}}$ и ${{\nu }_{E}}$ подбираются так, чтобы интегральные характеристики энергосодержания, плотности и энергетического времени удержания плазмы ${{\tau }_{E}}$ в проводимых расчетах соответствовали условиям базовых режимов в моделируемых экспериментах.

В предшествующих работах [20–30] мы предполагали, что все флуктуации обращаются в нуль как на магнитной оси (при $\rho = 0$), так и на внешней границе основной области плазмы (при $\rho = {{\rho }_{L}}$), что приводит к обращению в нуль на обеих границах и турбулентно-конвективных потоков ${{\Gamma }_{{turb}}}$ и ${{q}^{{turb}}}$. Условие обращения флуктуаций в нуль на внешней границе не столь очевидно и мотивируется тем, что пространственные структуры флуктуаций в основной области плазмы, где имеются замкнутые магнитные поверхности, и в области SOL, где силовые линии магнитного поля пересекают материальные поверхности, должны существенно различаться. При отсутствии явной связи между флуктуациями в основной области плазмы и в области SOL следует ожидать значительного снижения уровня флуктуаций на граничной поверхности $\rho = {{\rho }_{L}}$. Как показали результаты моделирования эволюции турбулентной плазмы, полученные в предшествующих работах, предположение об обращении флуктуаций в нуль на внешней границе оказывается вполне разумным и не нарушает самосогласованность задачи в целом. Поэтому в данной работе мы также используем нулевые граничные условия для флуктуаций на внешней границе и, соответственно, имеем ${{\left. {{{\Gamma }_{{turb}}}} \right|}_{{\rho = {{\rho }_{L}}}}} = 0$ и ${{\left. {{{q}^{{turb}}}} \right|}_{{\rho = {{\rho }_{L}}}}} = 0$ на границе основной плазмы с SOL.

Кроме того, принимая во внимание, что в моделируемых экспериментах плотность плазмы поддерживалась примерно постоянной на протяжении всего разряда, мы дополнительно упростили задачу. С этой целью предполагается, что источник частиц ${{\bar {Q}}_{n}}$ и поток $\Gamma (t,\rho )$ в уравнении (3) равны нулю на протяжении всего процесса. Тогда уравнение переноса частиц (3) и уравнение (12) работы [30] для флуктуаций плотности будут удовлетворяться при $\bar {D}(t,\rho ) = {\text{const}}$ и нулевых флуктуациях ${{D}_{f}}$. Поэтому, в рамках указанных предположений, уравнение (3) и уравнения (12) работы [30] можно полностью исключить из рассмотрения, и, без нарушения самосогласованности задачи, полагать, что в остальных уравнениях $\bar {D}(t,\rho ) = {\text{const}}$ и $\Gamma (t,\rho ) = 0$.

В предшествующих работах моделирование проводилось при граничном условии (8) с ${{\nu }_{E}} = {\text{const}}$. В этом случае успешно моделировались такие процессы, как формирование самосогласованных профилей давления плазмы [25–28], быстрый отклик на радиальное перераспределение вводимой мощности нагрева плазмы и быстрый отклик на выбросы энергии из области $q < 1$, вызываемые пилообразными колебаниями [28–30], а также ряд других процессов. Однако при моделировании переходных режимов с существенным изменением полной мощности нагрева, например, перехода от стадии омического нагрева к стадии с мощным ЭЦР-нагревом, ${{\tau }_{E}}$ в расчетах менялось относительно слабо. Это не согласуется со снижением (или деградацией) ${{\tau }_{E}}$ при увеличении мощности нагрева, которое наблюдается в экспериментах и учитывается в различных скейлингах для ${{\tau }_{E}}$. В целях усовершенствования расчетной методики и получения лучшего согласия результатов расчетов с экспериментами, в работах [32, 33] было предложено использовать обобщенные граничные условия 3-го рода, учитывающие специально подобранную нелинейную зависимость потоков тепла через внешнюю границу из основной плазмы в SOL от температур электронов и ионов на границе с SOL. Другими словами, предполагается, что ${{\nu }_{E}}$ зависит от температур электронов и ионов на границе. Там же было отмечено, что нелинейную зависимость потоков тепла через внешнюю границу от температур электронов и ионов на границе основной плазмы с SOL можно попытаться связать с закономерностями удержания плазмы в SOL, выраженными в форме некоторых степенных скейлингов.

К вопросу о разумном выборе нелинейной зависимости потока тепла из основной плазмы в SOL от температур электронов и ионов можно подойти двумя путями. Первый путь базируется на накопленном экспериментальном опыте, который дает ряд скейлингов, связывающих важнейшие интегральные характеристики удержания плазмы с параметрами плазмы и вводимой мощностью. Принято считать, что в стационарных режимах имеет место степенная зависимость энергетического времени удержания плазмы от вводимой мощности нагрева

где ${{Q}_{E}} = \int {({{P}_{{OH}}} + {{P}_{{ECR}}})V{\kern 1pt} 'd\rho } $ – полная вводимая мощность, а угловые скобки означают усреднение по объему плазмы. В частности, во многих случаях ожидаемую величину ${{\tau }_{E}}$ оценивают по так называемому многомашинному скейлингу стационарной Н-моды в ИТЭР ($\mathop {{{\tau }_{E}}}\nolimits^H (ITER - 98(y,2))$) с показателем $\alpha = 0.69$. В соответствии с (9), для средней по объему плотности тепловой энергии в стационарных режимах можно написать следующий скейлинг $\left\langle {n({{T}_{e}} + {{T}_{i}})} \right\rangle \propto {{({{Q}_{E}})}^{{1 - \alpha }}}$ и представить скейлинг (9) в форме: ${{\tau }_{E}} \propto {{\left\langle {n({{T}_{e}} + {{T}_{i}})} \right\rangle }^{{ - \alpha /(1 - \alpha )}}}$.

Тенденция к поддержанию самосогласованных (“канонических”) профилей давления, наблюдавшаяся во многих экспериментах [34–37], а также в нашем предшествующем моделировании (см., например, [30]), дает основания полагать, что давление плазмы на внешней границе основной области плазмы ${{\left. {n({{T}_{e}} + {{T}_{i}})} \right|}_{{\rho = {{\rho }_{L}}}}}$ должно поддерживаться примерно пропорциональным объемно-усредненному давлению в этой области $\left\langle {n({{T}_{e}} + {{T}_{i}})} \right\rangle $. В связи с этим, можно полагать, что в стационарных режимах должны выполняться следующее соотношение: ${{\left. {n({{T}_{e}} + {{T}_{i}})} \right|}_{{\rho = {{\rho }_{L}}}}} \propto {{({{Q}_{E}})}^{{1 - \alpha }}}$. С другой стороны, при ${{\left. \Gamma \right|}_{{\rho = {{\rho }_{L}}}}} = 0$ и ${{\left. {{{q}^{{turb}}}} \right|}_{{\rho = {{\rho }_{L}}}}} = 0$, а также нулевом потоке кинетической энергии флуктуаций скорости через внешнюю границу, что обеспечивается граничным условием сохранения интегральной завихренности (см. [28–30]), в стационарных режимах должно выполняться соотношение ${{\left. {{{q}^{{tot}}}} \right|}_{{\rho = {{\rho }_{L}}}}} = {{Q}_{E}}$. Как правило, все варианты наших расчетов эволюции турбулентной плазмы начинаются с квазистационарной ОН-стадии с ${{Q}_{E}} = {{Q}_{{OH}}}$, после которой включаются различные режимы нагрева плазмы. В этом случае, в соответствии с представленными выше соотношениями, нелинейное граничное условие (8) можно представить в следующей форме:

Важно подчеркнуть, что внешнее граничное условие для потока тепла в форме (10) следует использовать после окончания квазистационарной ОН-стадии на протяжении всей последующей нестационарной эволюции плазмы с произвольной последовательностью включения источников дополнительного нагрева. При этом ожидается, что граничное условие в форме (10) будет автоматически обеспечивать выход на квазистационарный скейлинг (9) для ${{\tau }_{E}}$. Граничное условие (10) легко обобщается на случай раздельного рассмотрения процессов переноса тепла в электронной и ионной компонентах плазмы и принимает следующий вид для потоков тепла в электронном и ионном уравнении соответственно:

В данной работе, мы хотим не просто формально ввести некоторые модельные нелинейные граничные условия 3-го рода, а сделать дополнительный шаг и связать конкретный вид этих граничных условий с некоторой, возможно упрощенной, но физически оправданной моделью удержания плазмы в SOL. Поэтому мы выбираем альтернативный путь для нахождения зависимости потоков тепла из основной плазмы в SOL в граничном условии (8) от локальных значений плотности плазмы и температур электронов и ионов на границе с SOL. Этот путь связан с анализом и оценкой физических механизмов, ответственных за удержание плазмы в SOL.

Прежде всего, следует заметить, что полная тепловая энергия плазмы в SOL всегда значительно меньше тепловой энергии плазмы в основной области. Поэтому в общем тепловом балансе любыми изменениями тепловой энергии в SOL можно пренебречь и считать, что суммарный поток тепла из основной плазмы в SOL равен суммарным потерям тепла из SOL не только в стационарных, но и в переходных режимах. Конкретную зависимость коэффициента ${{\nu }_{E}}$ в граничном условии (8) от локальных значений температур на границе с SOL получим в предположении, что главным каналом тепловых потерь из SOL является поток тепла вдоль силовых линий магнитного поля на материальные поверхности (лимитер или диверторные пластины).

В большинстве экспериментов и, в частности, в моделируемых далее экспериментах на токамаке Т-10 длина свободного пробега электронов вдоль силовых линий в SOL $( \sim {\kern 1pt} {{v}_{{Te}}}{{\tau }_{{ee}}})$ значительно меньше характерной длины пробега вдоль силовых линий до материальной стенки (диверторных пластин или лимитера) ${{L}_{c}}$, которую обычно называют “длиной связи” (connection length). В этом случае главным механизмом потерь тепла из SOL является хорошо известная классическая продольная электронная теплопроводность, при которой основной поток тепла из SOL пропорционален следующей величине:

где n и ${{T}_{e}}$ – некоторые эффективные средние значения плотности и температуры электронов в SOL, а ${{\Delta }_{s}}$ – характерная толщина SOL. Длина связи ${{L}_{c}}$ является геометрическим параметром магнитной системы и, как правило, не зависит от ${{T}_{e}}$. Полагая, что ${{\Delta }_{s}}$ также не зависит ${{T}_{e}}$, и что величина $({{T}_{e}} + {{T}_{i}})$ на границе с SOL примерно пропорционально ${{T}_{e}}$ в SOL, получаем степенную зависимость для $q_{\parallel }^{{SOL}} \propto ({{T}_{e}} + {{T}_{i}})_{{\rho = {{\rho }_{L}}}}^{{7/2}}$ от температур на границе с SOL. Из сопоставления выражения (10) для ожидаемой величины полного потока тепла из основной плазмы в SOL с приведенной оценкой для $q_{\parallel }^{{SOL}}$, получим следующее соотношение для показателя α: $1{\text{/}}(1 - \alpha ) = 7{\text{/}}2$, что дает значение $\alpha = 5{\text{/}}7 \approx 0.71$, близкое к значению $\alpha = 0.69$ в скейлинге H(ITER-98(y,2)).

Полученная оценка величины показателя α конечно не единственная и может зависеть от различных уточнений и изменений модели удержания плазмы в SOL. Так, в ряде оценок квазистационарной толщины SOL, учитывающих возможность развития дрейфовых неустойчивостей, толщина SOL оказывается пропорциональной ионно-звуковому ларморовскому радиусу, что дает ${{\Delta }_{s}} \propto \sqrt {{{T}_{e}} + {{T}_{i}}} $. Нетрудно видеть, что такое изменение модели удержания плазмы в SOL приводит к несколько большей величине $\alpha = 3{\text{/}}4 = 0.75$. Полученные значения показателя $\alpha $ сохраняются и в случае раздельного рассмотрения процессов переноса тепла в электронной и ионной компонентах плазмы с использованием граничных условий (11).

3. МОДЕЛИРОВАНИЕ ЭВОЛЮЦИИ ПЛАЗМЫ В СЦЕНАРИЯХ С НЕСТАЦИОНАРНЫМ ДОПОЛНИТЕЛЬНЫМ НАГРЕВОМ В ЭКСПЕРИМЕНТАХ НА ТОКАМАКЕ Т-10

В данном разделе представлены результаты моделирования эволюции турбулентной плазмы в трех сценариях, соответствующих трем разрядам на токамаке Т-10, который имеет большой радиус плазмы $R = 150$ см и круглое поперечное сечение плазмы с малым радиусом a = 30 см, ограниченным диафрагмой. Вследствие большого аспектного отношения $A = R{\text{/}}a \approx 5$ и относительно небольшого шафрановского смещения магнитной оси относительно граничной магнитной поверхности в моделируемых разрядах в Т‑10, для эффективного малого радиуса на внешней границе справедливо соотношение ${{\rho }_{L}} \approx a$. Эксперименты проводились при тороидальном поле ${{B}_{T}} \approx 2.41$ T. Моделируемые разряды имели близкие величины среднехордовой плотности плазмы $\bar {n} \approx 1.8 \times {{10}^{{19}}}$ м^–3, а энергетическое время удержания плазмы в них на ОН-стадии было на уровне ${{\tau }_{{EOH}}} \sim 30$ мс. При этом мощности нестационарного ЭЦР-нагрева плазмы в этих разрядах были существенно разными.

Моделирование проводилось для трех разрядов на токамаке Т-10 со следующими параметрами:

– в разряде #61203 вводимая мощность ${{Q}_{E}}$ изменялась в 2.6 раза от ${{Q}_{{EOH}}} \approx 0.26$ МВт на ОН-стадии до ${{Q}_{E}} \approx 0.68$ МВт на стадии ЭЦР-нагрева ($ \approx {\kern 1pt} 0.52$ МВт ЭЦР-мощности плюс $ \approx {\kern 1pt} 0.16$ МВт джоулева нагрева), время жизни на ОН-стадии составляло ${{\tau }_{{EOH}}} \approx 32$ мс, мощность нецентрального ЭЦР-нагрева вводилась вблизи поверхности $\rho \approx 6.5$ см;

– в разряде #61200 вводимая мощность ${{Q}_{E}}$ изменялась в 4.6 раза от ${{Q}_{{EOH}}} \approx 0.26$ МВт на ОН-стадии до ${{Q}_{E}} \approx 1.20$ МВт на стадии ЭЦР-нагрева (1.10 МВт ЭЦР плюс 0.10 МВт ОН), время жизни на ОН-стадии составляло ${{\tau }_{{EOH}}} \approx 28$ мс, мощность нецентрального ЭЦР-нагрева вводилась вблизи поверхности $\rho \approx 6.5$ см;

– в разряде #61208 вводимая мощность ${{Q}_{E}}$ изменялась в 8.3 раза от ${{Q}_{{EOH}}} \approx 0.24$ МВт на ОН-стадии до ${{Q}_{E}} \approx 1.98$ МВт на стадии ЭЦР-нагрева (1.95 МВт ЭЦР плюс 0.03 МВт ОН), время жизни на ОН-стадии составляло ${{\tau }_{{EOH}}} \approx 29$ мс. ЭЦР-нагрев состоял из комбинации центрального нагрева 0.85 МВт вблизи $\rho \approx 1.5$ см и 1.10 МВт нецентрального ЭЦР-нагрева вблизи поверхности $\rho \approx 6.5$ см.

Описание моделируемых разрядов дано в работе [44]. Радиальные профили поглощаемой ЭЦР-мощности в эксперименте рассчитывались квазилинейным кодом OGRAY [45], а мощности джоулева нагрева на стационарных стадиях экспериментальных разрядов рассчитывались с помощью кода ASTRA [6].

Как отмечалось выше, компьютерное моделирование проводилось с использованием кода CONTRA-C (цилиндрическая модель токамака с круглым сечением плазмы и объединенным уравнением теплопереноса для электронов и ионов при фиксированном отношении ${{T}_{e}}{\text{/}}{{T}_{i}}$). Каждый из расчетных сценариев начинался со стадии омического нагрева длительностью ${{t}_{{OH}}} = 45$ мс, на которой происходило первоначальное развитие и установление режима самосогласованной турбулентной конвекции. Внешнее граничное условие для потока тепла бралось в форме (10). На начальной ОН-стадии величина суммарной температуры на внешней границе ${{\left. {{{{({{T}_{e}} + {{T}_{i}})}}_{{OH}}}} \right|}_{{\rho = {{\rho }_{L}}}}}$ в этом граничном условии подбиралась так, чтобы обеспечить стационарное поддержание параметров плазмы (включая среднюю плотность тепловой энергии и энергетическое время удержания плазмы ${{\tau }_{{EOH}}}$), соответствующими стационарной OH-стадии моделируемого разряда в Т-10. На последующих стадиях (при $t > {{t}_{{OH}}}$) полученная величина ${{\left. {{{{({{T}_{e}} + {{T}_{i}})}}_{{OH}}}} \right|}_{{\rho = {{\rho }_{L}}}}}$ сохранялась неизменной, а текущее значение ${{\left. {({{T}_{e}} + {{T}_{i}})} \right|}_{{\rho = {{\rho }_{L}}}}}$ в граничном условии (10) самосогласованно менялось в результате изменения вводимой мощности, что обеспечивало автоматическое изменение ${{\tau }_{E}}$ в процессе эволюции плазмы в соответствии с меняющимися параметрами плазмы.

При численном моделировании радиальные профили вводимой ЭЦР-мощности брались в форме гауссовых распределений с $\sigma = 1.5$ см для нецентрального ЭЦР-нагрева и $\sigma = 3.0$ см для центрального нагрева, что достаточно хорошо аппроксимировало профили, рассчитанные кодом OGRAY. Предполагалось, что вводимые мощности ЭЦР-нагрева выходят на номинальное значение через 8 мс после их включения. В данном моделировании мы не проводили детальных расчетов эволюции джоулева нагрева, а использовали данные о мощностях джоулева нагрева, полученные на стационарных стадиях соответствующих экспериментальных разрядов, полагая, что время изменения мощностей джоулева нагрева при включении ЭЦР-нагрева составляет 15 мс. Это время превышает время нарастания вводимой ЭЦР-мощности, а, с другой стороны, примерно соответствует величине ${{\tau }_{E}}$ на стадии ЭЦР-нагрева. Указанные предположения не влияют на конечную величину ${{\tau }_{E}}$ на квазистационарной стадии ЭЦР-нагрева, а лишь обеспечивают достаточно естественный ход эволюции ${{\tau }_{E}}$ в переходном процессе.

В реальных разрядах #61200 и #61203 на Т-10 время поддержания квазистационарной стадии ЭЦР-нагрева составляло почти 200 мс. Полная длительность в численном моделировании сценариев разрядов также составляла вполне макроскопическую величину 130 мс, которая была вполне достаточной для целей данной работы, поскольку значительно превышала характерные времена релаксации к квазистационарным состояниям. В разряде #61208 длительность ЭЦР-стадии с наибольшей вводимой мощностью была короче и составляла 138 мс. Причем сначала отключалась примерно половина полной мощности, а через 9 мс остальная часть. Это обстоятельство навело на мысль продемонстрировать самосогласованное изменение ${{\tau }_{E}}$ при двухступенчатом отключении ЭЦР-мощности и в численном моделировании. Однако мы сократили длительность ЭЦР-стадии с наибольшей вводимой мощностью до 45 мс, что вполне достаточно для демонстрации полной релаксации к квазистационарному состоянию на этой стадии, а интервал между последовательным отключением сначала 0.85 МВт центрального нагрева, а потом 1.10 МВт нецентрального нагрева увеличили до 34 мс, чтобы продемонстрировать релаксацию ${{\tau }_{E}}$ на промежуточной стадии отключения ЭЦР-нагрева.

Ниже будет представлена эволюция энергетического времени удержания плазмы ${{\tau }_{E}}(t)$ в двух определениях. Первое “стационарное” определение $\tau _{E}^{{st}}(t) = 3V\left\langle {n({{T}_{e}} + {{T}_{i}})} \right\rangle {\text{/}}2{{Q}_{E}}$ было введено выше формулой (9). Оно соответствует часто используемой при описании экспериментов характеристике удержания плазмы на той или иной стационарной стадии разряда. Однако, при описании переходных процессов в условиях нестационарного дополнительного нагрева плазмы более адекватной характеристикой удержания плазмы представляется величина $\tau _{E}^{{tr}}(t)$, задаваемая следующим соотношением нестационарного теплового баланса:

а введенная ранее величина $\tau _{E}^{{st}}(t)$ выступает лишь как вспомогательная формальная характеристика, позволяющая более наглядно проследить переход в квазистационарное состояние, в котором $\tau _{E}^{{tr}} \approx \tau _{E}^{{st}}$.

В предшествующих работах [32, 33] хорошее согласие с результатами экспериментов на токамаке Т-10 было получено при выполнении численного моделирования эволюции ${{\tau }_{E}}(t)$ с показателем $\alpha = 0.7$ во внешнем граничном условии в форме (10). Однако, позднее было обнаружено, что в указанной серии расчетов имела место некоторая некорректность, связанная с небольшой утечкой на магнитную ось кинетической энергии тороидального вращения, что влияло на общий энергобаланс. В расчетах, представляемых в данной работе, указанная некорректность была устранена, но после исправления оказалось, что на стадии ЭЦР-нагрева расчетное ${{\tau }_{E}}$ при $\alpha = 0.71$ стало заметно превышать значения ${{\tau }_{E}}$ в соответствующих экспериментальных разрядах. Поэтому было решено провести моделирование эволюции плазмы в обсуждавшихся выше сценариях с увеличенным показателем $\alpha = 0.75$, учитывающим увеличение толщины SOL с ростом температуры плазмы.

На рис. 1 представлены результаты моделирования эволюции энергетического времени удержания плазмы ${{\tau }_{E}}(t)$ при α = 0.75 для обсуждавшихся выше сценариев разрядов плазмы в Т-10. Черные сплошные кривые 1 соответствуют $\tau _{E}^{{tr}}(t)$, задаваемому соотношением (12). Серые сплошные кривые 2 представляют эволюцию $\tau _{E}^{{st}}(t)$. Штриховыми линиями с крупными штрихами показаны стационарные уровни ${{\tau }_{E}}$, рассчитанные на основе экспериментальных наблюдений в соответствующих разрядах на Т-10 с указанием погрешности измерений на стадии ЭЦР-нагрева. Штриховыми линиями с мелкими штрихами показаны стационарные уровни ${{\tau }_{E}}$, соответствующие скейлингу (9) для зависимости ${{\tau }_{E}}$ от ${{Q}_{E}}$ с показателем $\alpha = 0.69$.

Рисунок 1 показывает, что во всех трех сценариях величины $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$ к концу ОН-стадии достаточно надежно выходят на квазистационарный уровень, соответствующий уровню ${{\tau }_{E}}$ на ОН-стадии соответствующего разряда в Т-10. После включения ЭЦР-нагрева величина $\tau _{E}^{{st}}(t)$ демонстрирует глубокий провал, связанный с тем, что мощность нагрева ${{Q}_{E}}$ нарастает во времени значительно быстрее, чем полное энергосодержание плазмы. Как и следовало ожидать, величина $\tau _{E}^{{tr}}(t)$ при нагреве снижается более плавно и монотонно, чем формальная величина $\tau _{E}^{{st}}(t)$. Небольшой, но быстрый сброс $\tau _{E}^{{tr}}(t)$ на начальной стадии ЭЦР-нагрева связан с тем, что на этой стадии прирост тепловой энергии плазмы за счет нарастающей ЭЦР-мощности почти полностью компенсируется расходом тепловой энергии на увеличение кинетической энергии турбулентности. Последующий, относительно плоский участок на кривой эволюции $\tau _{E}^{{tr}}(t)$ связан с тем, что уже возросшие турбулентно-вихревые течения, за счет вязкости, начинают отдавать свою увеличенную кинетическую энергию назад в тепловую энергию плазмы. В упрощенном тепловом балансе (12) процесс передачи тепловой энергии в кинетическую энергию турбулентно-вихревых течений и обратно не учитывается, поскольку его практически невозможно рассчитать по экспериментальным данным. Однако его проявление в форме некоторой немонотонности кривой эволюции $\tau _{E}^{{tr}}(t)$, определяемой уравнением упрощенного теплового баланса (12), может наблюдаться не только в расчетах, но и при анализе экспериментальных данных.

По истечению времени, составляющего примерно $2.5{{\tau }_{E}}$ на стадии дополнительного нагрева, величины $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$ практически сравниваются и выходят на квазистационарный уровень, который, несмотря на увеличенное значение показателя $\alpha = 0.75$, во всех трех сценариях оказывается несколько выше экспериментального значения ${{\tau }_{E}}$ и выходит за пределы погрешности измерений.

Рисунок 1 также показывает, что при моделировании сценария, соответствующего разряду #61208, после отключением 0.85 МВт центрального ЭЦР-нагрева величины $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$, как и ожидалось, начинают автоматически расти и релаксируют к квазистационарному уровню, близкому к квазистационарному уровню ${{\tau }_{E}}$ в разряде #61200, имеющему аналогичный уровень мощности ЭЦР-нагрева ≈1.10 МВт. На стадии окончательного отключения ЭЦР-нагрева (при $t > 124$ мс) величины $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$ начинают эволюционировать к уровню ${{\tau }_{E}}$ на исходной ОН-стадии. Таким образом, внешнее граничное условие в форме (10) обеспечивает автоматическую подстройку величин $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$ под уровень вводимой мощности, как при включении, так и при отключении дополнительного нагрева.

Превышение расчетного ${{\tau }_{E}}$ над экспериментальным на ЭЦР-стадии разряда в проведенной серии расчетов может быть связано с тем, что, наряду с потоком тепла в SOL, не были учтены дополнительные потери энергии. С этой целью была проведена новая серия расчетов, в которой были учтены потери на излучение, играющие доминирующую роль среди ранее неучтенных дополнительных потерь энергии. Учет потерь на излучение сводится к тому, что в плотность мощности источников нагрева в уравнении переноса энтропийной функции (5) нужно включить плотность мощности потерь на излучения ${{P}_{{rad}}}$, а во внешнем граничном условии в форме (10) мощность нагрева на ОН-стадии ${{Q}_{{E{\text{OH}}}}}$ заменить на $({{Q}_{{E{\text{OH}}}}} - {{Q}_{{rad{\text{OH}}}}})$, где ${{Q}_{{rad{\text{OH}}}}} = \int {{{P}_{{rad{\text{OH}}}}}V{\kern 1pt} 'd\rho } $. На основании имеющихся экспериментальных данных была получена следующая оценка величины потерь на излучение на ОН-стадии ${{Q}_{{rad{\text{OH}}}}}$: примерно 85 кВт для разрядов #61200 и #61208 и 95 кВт в разряде #61203, которые мы и использовали в последующем моделировании. Оценка потерь на излучение на ЭЦР-стадии дала следующие значения для ${{Q}_{{rad{\text{OH}}}}}$: 140 кВт для разряда #61203, 150 кВт для разряда #61200 и 180 кВт для разряда #61208. Плотность мощности потерь на излучения ${{P}_{{rad}}}$ из основной плазмы предполагалась локализованной в слое толщиной примерно 5 см, граничащем с область SOL. Как выяснилось, учет потерь на излучение не только не снизил, а в итоге увеличил превышение расчетного ${{\tau }_{E}}$ над экспериментальным на ЭЦР-стадии разряда. Это связано с тем, что, несмотря на увеличение мощности потерь на излучение на ЭЦР-стадии, отношение мощности потерь на излучение к полной мощности нагрева при увеличении нагрева снижается.

Другая причина превышения расчетного ${{\tau }_{E}}$ над экспериментальным на ЭЦР-стадии разряда проиллюстрирована рис. 2, на котором представлены нормированные радиальные профили энтропийной функции $\bar {S}$ на ОН-стадии и на поздней стадии ЭЦР-нагрева, полученные в ходе моделирования импульса #61200. Штриховыми линиями показаны профили турбулентно-релаксированных состояний с теми же величинами объемно-усредненной энтропийной функции. Видно, что расчетные профили нормированных функций $\bar {S}$ на рис. 2, как и предполагалось, близки к TR-состояниям почти всюду, что соответствует ожидаемой самосогласованности профилей давления. Исключение демонстрирует тонкий слой вблизи внешней границы, где образуется повышенный градиент энтропийной функции, в результате чего, при нагреве плазмы, рост ${{\left. {\bar {S}} \right|}_{{\rho = {{\rho }_{L}}}}}$ (и, соответственно, ${{\left. {\bar {n}({{{\bar {T}}}_{e}} + \bar {T})} \right|}_{{\rho = {{\rho }_{L}}}}}$) на внешней границе основной плазмы существенно отстает от роста объемно-усредненной энтропийной функции. Большой градиент $\bar {S}$ вблизи внешней границы является следствием предположения, что все флуктуации, а, соответственно, и турбулентный поток тепла, задаваемый выражением (6), обращаются в нуль при $\rho = {{\rho }_{L}}$. Этот результат можно интерпретировать как появление мощного внешнего транспортного барьера (ETB) на ЭЦР-стадии наших расчетных сценариев. В действительности в моделируемой серии экспериментальных разрядов сколько-нибудь заметного ETB не наблюдалось. Как показывает эксперимент, все флуктуации хотя и значительно снижаются при $\rho = {{\rho }_{L}}$, но не обращаются в нуль, сохраняя уменьшенную, но конечную величину турбулентного потока тепла, задаваемого выражением (6). Чтобы учесть это обстоятельство в рамках предлагаемой в данной работе концепции, достаточно в нелинейное граничное условие для потока тепла из основной области в SOL в форме (10) подставлять температуры, рассчитанные не в точке $\rho = {{\rho }_{L}}$, а в точке $\rho = ({{\rho }_{L}} - \delta )$, слегка смещенной внутрь основной плазмы, где турбулентные флуктуации уже значительно снижены, но турбулентный поток тепла еще не обратился в нуль. Тогда граничное условие (10) с учетом потерь на излучение и небольшой коррекции в выборе “граничной точки” примет вид

На рис. 3, для сравнения, приведена эволюция ${{\tau }_{E}}(t)$, рассчитанная в трех вариантах сценария моделирования разряда #61208 в токамаке Т-10 при $\alpha = 0.75$. Исходный сценарий соответствует внешнему граничному условию в форме (10) и демонстрирует превышение расчетного ${{\tau }_{E}}$ над экспериментальным на ЭЦР-стадии разряда. Второй сценарий выполнен с дополнительным учетом потерь на излучение. В нем превышение расчетного ${{\tau }_{E}}$ над экспериментальным оказывается еще больше, чем в первом сценарии. В обоих этих сценариях превышение расчетного ${{\tau }_{E}}$ выходит за пределы погрешности измерений. Сценарий #61208-3 выполнен с использованием модифицированного внешнего граничного условия в форме (13), как с учетом потерь на излучение, так и с коррекцией “граничной точки”. В данных расчетах была выбрана величина смещения “граничной точки” $\delta = 5$ мм, что составляло примерно 1.5% от величины малого радиуса плазмы. Видно, что в этом сценарии превышения расчетного ${{\tau }_{E}}$ над экспериментальным нет.

На рис. 4 представлены результаты моделирования эволюции ${{\tau }_{E}}(t)$, выполненные с использованием модифицированного внешнего граничного условия в форме (13) при $\alpha = 0.75$ для всех трех обсуждавшихся выше разрядов в Т-10. Видно, что этом случае, при всех уровнях вводимой ЭЦР-мощности, различие между расчетными и экспериментальными величины ${{\tau }_{E}}$ на поздней ЭЦР-стадии существенно меньше, чем при моделировании этих же разрядов с использованием внешнего граничного условия в форме (10) (см. рис. 1) и не выходят за пределы погрешности экспериментальных измерений. Таким образом, модифицированное граничное условие в форме (13) дает хорошую основу для численного моделирования (в том числе, и прогнозирования) процессов эволюции турбулентной плазмы в токамаках в различных режимах с нестационарным дополнительным нагревом плазмы. Обобщение граничного условия в форме (13) на случай расчетов с раздельными уравнениями переноса тепла в электронной и ионной компонентах выполняется достаточно тривиально, как простая модификация граничных условий (11).

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

В целях дальнейшего исследования физической природы аномальных транспортных процессов в токамаках и совершенствования методики самосогласованного компьютерного моделирования низкочастотной турбулентной конвекции плазмы и результирующих аномальных процессов переноса тепла проведено численное моделирование эволюции турбулентной плазмы в основной горячей области токамака. Моделирование проводилось в нестационарных режимах удержания плазмы, включающих переход от квазистационарной стадии омического нагрева к стадии мощного дополнительного ЭЦР-нагрева.

В разделе 2 отмечено, что в присутствии развитой турбулентной конвекции, демонстрирующей тенденцию к поддержанию самосогласованных профилей давления плазмы, в отличие от традиционной диффузионной модели процессов переноса, тепловые потоки через внешнюю границу основной области плазмы в область SOL могут существенно зависеть от механизма тепловых потерь плазмы из самой области SOL. Там же предложен и проанализирован специальный класс обобщенных внешних граничных условий 3-го рода вида (8), которые, вследствие относительно малой собственной теплоемкости области SOL, позволяют связать тепловые потоки через внешнюю границу основной области плазмы с тепловыми потерями плазмы из области SOL и выразить эту связь в форме нелинейной зависимости коэффициента ${{\nu }_{E}}$ в граничном условии (8) от локальных величин температур и плотности на границе основной области плазмы и области SOL. Принимая во внимание, что в моделируемых экспериментах на токамаке Т-10 плотность плазмы поддерживалась примерно постоянной во времени на протяжении всего разряда, задача была дополнительно упрощена путем полного исключения процессов переноса частиц, выполненного без нарушения самосогласованности задачи, и использования заданного и независящего от времени профиля плотности плазмы. Это позволило акцентировать основное внимание в данной работе на вопросах переноса тепла.

Конкретная зависимость коэффициента ${{\nu }_{E}}$ в граничном условии (8) от локальных значений температур на границе с SOL получена в данной работе в предположении, что главным каналом тепловых потерь из SOL является поток тепла вдоль силовых линий магнитного поля на материальные поверхности (лимитер или диверторные пластины) и что этот поток задается механизмом классической продольной электронной теплопроводности. В этом случае плотность потока тепла из SOL пропорциональна $T_{e}^{{7/2}}$, что приводит к степенной зависимости коэффициента ${{\nu }_{E}} \propto T_{e}^{{5/2}}$. В разделе 2 также показано, что, в силу тенденции к поддержанию самосогласованных профилей давления плазмы, такая степенная зависимость коэффициента ${{\nu }_{E}}$ довольно хорошо коррелирует со степенной зависимостью ${{\tau }_{E}} \propto {{({{Q}_{E}})}^{{ - \alpha }}}$ с α = $ = 5{\text{/}}7 \approx 0.71$, которое близко к так называемому многомашинному скейлингу стационарной Н-моды в ИТЭР ($\tau _{E}^{H}({\text{ITER}} - 98({\text{y}},2))$) с показателем $\alpha = 0.69$. Различные изменения и уточнения модели удержания плазмы в SOL могут приводить к изменению оценки величины показателя α. В частности, ряд оценок квазистационарной толщины SOL, учитывающих возможность развития дрейфовых неустойчивостей, приводит к толщине SOL, пропорциональной ионно-звуковому ларморовскому радиусу, что увеличивает показатель степени в оценке величины полного продольного потока тепла из SOL и делает его пропорциональным $T_{e}^{4}$. Такое изменение модели удержания плазмы в SOL приводит к величине $\alpha = 3{\text{/}}4 = 0.75$.

Для тестирования и анализа свойств предложенных граничных условий, дополнивших существующий компьютерный код CONTRA-C, было проведено численное моделирование эволюции энергетического времени удержания плазмы ${{\tau }_{E}}(t)$ в сценариях с переходными режимами, вызванными ЭЦР-нагревом плазмы. Эти сценарии соответствовали трем реальным экспериментам на токамаке Т-10, проведенным с существенно различным уровнем вводимой ЭЦР-мощности. Выполненное моделирование продемонстрировало, что внешнее граничное условие в форме (10) обеспечивает автоматическую подстройку величин $\tau _{E}^{{tr}}(t)$ и $\tau _{E}^{{st}}(t)$, соответствующих двум возможным определениям энергетического времени удержания плазмы, под уровень вводимой мощности как при включении, так и при отключении дополнительного нагрева. Показано, что для описании переходных процессов, вызываемых включением и выключением дополнительного нагрева плазмы, более адекватной характеристикой удержания плазмы представляется величина $\tau _{E}^{{tr}}(t)$, определяемая соотношением (12).

В отличие от предшествующих расчетов, представленных в работах [32, 33], в данной серии расчетов моделирование эволюции плазмы проводилось при граничном условии (10) с увеличенным показателем $\alpha = 0.75$, который учитывал увеличение толщины SOL с ростом температуры плазмы. Результаты моделирования, однако, показали, что и в этом случае граничное условие в форме (10) приводит к небольшому превышению расчетного ${{\tau }_{E}}$ над экспериментальным на квазистационарной стадии ЭЦР-нагрева.

В ходе анализа причин превышения расчетного ${{\tau }_{E}}$ над экспериментальным было проведено моделирование возможного влияния потерь на излучение, как доминирующего механизма дополнительных потерь энергии из основной плазмы, на эволюцию и величину расчетного ${{\tau }_{E}}$. Плотность мощности потерь на излучения ${{P}_{{rad}}}$ предполагалась локализованной в относительно тонком слое, граничащем с область SOL. Как выяснилось, учет потерь на излучение не только не снизил, а в итоге увеличил превышение расчетного ${{\tau }_{E}}$ над экспериментальным на ЭЦР-стадии разряда.

Далее было показано, что более существенной причиной некоторого отличия ${{\tau }_{E}}$ в наших расчетах от эксперимента связано с появлением на ЭЦР-стадии наших расчетных сценариев довольно значительного ETB, которого нет в моделируемой серии экспериментальных разрядов. Появление такого фиктивного ETB в нашем моделировании связано с предположением об обращении флуктуаций и, соответственно, турбулентных потоков в нуль на границе с SOL. Поскольку такие граничные условия для флуктуаций оказывают очень локализованное влияние на профиль давления плазмы, введение небольшой коррекции в выборе “граничной точки”, позволило устранить влияние фиктивного ETB в расчетных сценариях.

Применение модифицированного внешнего граничного условия (13), учитывающего потери на излучение и коррекцию в выборе “граничной точки”, позволило существенно снизить различия между расчетными и экспериментальными величинами ${{\tau }_{E}}$ так, что на поздней ЭЦР-стадии эти различия перестали выходить за пределы погрешности экспериментальных измерений. Таким образом, модифицированное граничное условие в форме (13) дает хорошую основу для численного моделирования (в том числе, и прогнозирования) процессов эволюции турбулентной плазмы в токамаках в различных режимах с включением дополнительного нагрева плазмы.

Полученные результаты показывают, что предложенный подход, в принципе, позволяет более адекватно моделировать самосогласованную эволюцию турбулентной плазмы на макроскопических временах при включении и выключении различных комбинаций дополнительного нагрева. Существенно, что при использовании такого подхода вместо решения весьма сложной самосогласованной транспортной задачи в SOL, достаточно ограничиться эмпирическим подбором численного коэффициента во внешнем граничном условии (8) на начальной квазистационарной стадии с омическим нагревом.

Проведенный анализ может быть весьма полезен как для анализа существующих экспериментов, так и для планирования новых экспериментов в новых, еще не освоенных режимах.

Работа выполнена при частичной поддержке Госкорпорации “Росатом”.