Анализ влияния эффекта Холла на структуру неравновесного плазменного слоя в канале МГД-генератора Текст научной статьи по специальности «Физика»

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

Том 12, № 4, 2007

АНАЛИЗ ВЛИЯНИЯ ЭФФЕКТА ХОЛЛА НА СТРУКТУРУ НЕРАВНОВЕСНОГО ПЛАЗМЕННОГО СЛОЯ В КАНАЛЕ МГД-ГЕНЕРАТОРА

И. А. Кузоватов, Т. А. МиловидовА

Красноярский государственный технический университет, Россия

e-mail: kuzovatov@yandex.ru

В. С. Славин Институт теплофизики СО РАН, Новосибирск, Россия

2D modelling of non-stationary interaction processes between non-equilibrium plasma layer, a magnetic field and a turbulent supersonic gas flow in a MHD generator channel is carried out. The influence of the Hall effect on plasma layer structure is analyzed.

Введение

Многолетний предыдущий этап исследования МГД-процессов в неоднородных газоплазменных потоках, в которых плазменные слои чередуются со слоями неэлектропроводного газа, был основан на численном моделировании одномерных процессов [1, 2]. При этом молчаливо предполагалось, что плазменный слой обладает пространственной устойчивостью, т. е. сохраняет поршневую структуру на протяжении всего пролетного времени в канале МГД-генератора. Важным итогом одномерного моделирования неравновесных плазменных слоев (П-слоев) явилось открытие эффекта "замороженной ионизации" [i], суть которого — в резком замедлении процесса рекомбинации в плазме чистого (без щелочной присадки) инертного газа при высоких значениях электронной температуры (> 5QQQK).

Одномерное численное моделирование генераторного процесса в рабочей части МГД-канала показало, что при магнитном поле до б Тл и частоте появления П-слоев, обеспечивающей одновременное присутствие в канале четырех—пяти плазменных слоев, достигается оптимальный режим, в котором степень преобразования энтальпии приближается к 4Q%, а адиабатическая эффективность составляет более 7Q% [2]. Такие характеристики соответствуют требованиям к эффективному МГД-генератору.

Однако, как было сказано выше, одномерная постановка задачи исключала возможность анализа влияния турбулентности, явления Холла и других пространственных эффектов на структуру плазменных слоев. Поскольку параметр Холла в П-слое имеет

(( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

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

1. Постановка задачи

Будем исследовать неоднородный газоплазменный поток с неравновесным плазменным слоем, движущийся в канале постоянного сечения. Расчетная область представляет собой прямоугольник, ограниченный сверху и снизу электродными стенками канала, а слева и справа — контактными границами потока, движущимися вместе с плазменным слоем. Выберем направление осей следующим образом: ось х — вдоль течения, ось у — перпендикулярно электродным стенкам (рис. 1).

Постоянное внешнее магнитное поле однородно заполняет весь объем МГД-канала. На электродных стенках канала задается постоянная разность потенциалов.

Математическая модель основывается на системе двумерных уравнений электродинамики, магнитной газодинамики и уравнениях на концентрацию и температуру электронов.

Закон Ома с учетом эффекта Холла имеет вид

j = а (Е + и х В) - в ( х В),

где а — скалярная электропроводность среды, и (и, и) — вектор скорости, в — параметр Холла.

Рис. 1. Расчетная область

Соответственно, компоненты вектора плотности тока, принимая во внимание связь между напряженностью и потенциалом, выражаемую уравнением Е = —grad т, можно записать в следующем виде:

(

а

1+ в2

а

V 1 + в2

— дт + в^т + В (ви + и)

дх ду

0ТГ + + В (и — ви)

х у

\

/

Теперь, используя соотношение div j = 0, получим эллиптическое диффузионно-конвективное уравнение на потенциал. Данное уравнение в дивергентной форме записи конвективных слагаемых примет вид

д ( а дт д

+ 7^

дх\ 1 + в2 дх ду\ 1 + в2

а

д

в т — тт-

а дт д

у 1 + в2 у х 1 + в2

а

в т

д д — — (аВ (ви + и)) + — (аВ (и — ви)). х у

(1.1)

Уравнения магнитной газодинамики могут быть представлены в векторной форме:

ди дБ дС дРи дС и

--1---1---1--- +---

дЬ дх ду дх ду

Б,

(1.2)

и

р ри ри

ри , Р = ри2 + Р , С = рии ри2 + Р , Б =

ри рии

е и(е + Р) и(е + Р)

0

—З*в

X (Те — Т) + ЗуиВ

Компоненты вектора вязких сил и связанных с ними тепловых потоков представляются как

Р-

0

Тх

— Тх

ху

Ях (итхх + иТху)

С -

0

Ту

ух

Т

уу

Яу — (иТух + иТуу)

Здесь используются стандартные обозначения основных газодинамических величин: р — плотность; Р — давление; Т — температура газа; Те — температура электронов; е — полная удельная энергия; ^, V — молекулярная и турбулентная вязкость; тензор напряжений Тг^ определяется как сумма тензора вязких напряжений т— и рейнольдсова тензора турбулентных напряжений тТ, т. е. тг^ = т— + тТ; аналогично поток тепла яг определяется как сумма молекулярного потока тепла и турбулентного потока тепла ; х — коэффициент энергообмена электронов и тяжелых частиц.

В настоящей работе для описания турбулентных характеристик течения использовалась однопараметрическая дифференциальная модель турбулентности ^-92, описывающая эффективную турбулентную вязкость [4]. Система уравнений газовой динамики (1.2) дополняется уравнением, описывающим эволюцию турбулентной вязкости:

д д д ~дЬ (т) + дх (рищ) + ду (Р^ —

-щ»+а э - +а щ)- ^ ™

где выражение (Ри - Пи) учитывает генерацию и диссипацию турбулентной вязкости [4].

Температуру электронов определим из уравнения энергобаланса для электронного газа:

^^ - X (Те - Т) - (/ + 3 къ Т^ (ие)3 = 0, (1.4)

где I — потенциал ионизации, къ — постоянная Больцмана, (Пе)3 — скорость изменения концентрации электронов, определяемая при решении уравнений ионизационно-реком-бинационной кинетики [5]. Уравнение на электронную температуру включает джоулеву диссипацию, передачу энергии в упругих столкновениях электронов с атомами и ионами, выделение энергии при рекомбинации и ее поглощение при ионизации.

В уравнении переноса концентрации электронов помимо амбиполярной диффузии учитывается также и турбулентный перенос:

д \ * д ( ^ ■ дпе

-7^'Пе + ТГ" (иПе) + 7Т" (™Пе) - "Т" \ Оа +--1

г\. ' ve i с\ \y-vi i i \ »ve ¡ о \ I ai /о

ot ox oy dx W Ote ) dx

-Í ((Da + ^) ^ 1 =(Ue), , (1.5)

dy W Ote) dy

где Da — коэффициент амбиполярной диффузии; ote = 1 — число Шмидта для концентрации электронов.

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

В качестве начальных условий задано стационарное турбулентное течение неэлектропроводного газа в канале постоянного сечения.

На стенках канала ставятся естественные для вязкого газа граничные условия — условия "прилипания" для скорости. Граничные условия на стенке дополняются условиями изотермичности, равенства нулю поперечного градиента давления и равенства нулю турбулентной вязкости:

T = Tw, ^ = 0, Vt = 0. dy

2. Алгоритм численного решения системы уравнений 2.1. Общая схема алгоритма

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

Для дискретизации уравнений в расчетной области Q = {0 < x < L, 0 < y < H} введем прямоугольную сетку

Üh = {(xi, yj) | xi = i • hx, yj = j • hy, 0 < i < Nx, 0 < j < Ny} .

Исследуемая система дифференциальных уравнений включает в себя уравнения разных типов, т. е. является системой составного типа. Дискретизация каждого из уравнений, входящих в систему (1. 1)—(1.5), производится наиболее подходящим способом. Об основных моментах этого процесса, а также о методах решения полученных сеточных уравнений речь пойдет в следующих пунктах.

Рассмотрим последовательность численного интегрирования уравнений системы (1.1)—(1. 5) в виде общего алгоритма.

1. Задаем начальные условия: П-слой задается в виде синусоидального возмущения концентрации электронов, на которое накладывается изобарическое возмущение температуры газа. Начальное распределение остальных газодинамических параметров задается однородным.

2. Решается блок электродинамических уравнений (1.1), (1.4), (1.5).

2.1. Решается уравнение (1.4) на электронную температуру Те, определяется распределение в расчетной области параметра Холла З, а также (пе)3.

2.2. Решается уравнение (1.1), по найденному распределению потенциала т вычисляем значение напряженности электрического поля Е и плотность тока j.

2.3. Из решения уравнения (1.5) находим распределение электронной плотности пе.

3. Решаются уравнения газовой динамики (1.2) и турбулентной вязкости (1.3).

Излагая общую схему алгоритма, необходимо упомянуть, что разномасштабность

характерного времени протекания электродинамических и газодинамических процессов требует при одном шаге по времени для расчета уравнений газодинамики выполнения многократного пересчета блока уравнений электродинамики.

2.2. Уравнение электродинамики

Введем обозначения:

/ 5

а

1 + в2

Ь

у 1 + в2 а

а З^

в

д д f = — (аВ (/и + и)) + — (аВ (и — /и)). х у

/

х 1 + в2

Тогда уравнение электродинамики (1.1) можно записать в следующем виде:

—div (а^т — Ьт) = f■

(2.1)

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

а

1 + з2

—дт + здг + В (З3и + и)

х у

0.

Уравнение (2.1) относится к числу диффузионно-конвективных уравнений с малым коэффициентом диффузии, численное решение которых при помощи традиционных (например, конечно-элементных) методов сопряжено со значительными трудностями. Центрально-разностная аппроксимация конвективного слагаемого приводит к услов-

Ь • к

но устойчивой разностной схеме при выполнении условия Ие^, < 1, где Ие^, = - —

2а

а

разностное число Рейнольдса. Выполнение данного условия накладывает жесткое ограничение на шаг сетки. Применение односторонних, направленных разностей ведет к появлению численной диффузии.

Один из способов решения данной проблемы — метод экспоненциальной подгонки, а одна из первых публикаций, где приводится строгое обоснование данного подхода, — работа [6]. Суть метода состоит в локальной перестройке разностного оператора; достаточно полное изложение различных способов применения экспоненциальной подгонки при построении разностных схем, доказательство абсолютной устойчивости методов решения и равномерной сходимости по малому параметру приведены в обзорной монографии [7]. Симметризованный способ записи разностных схем с экспоненциальной подгонкой обсуждался в [8]. В данной работе применение метода экспоненциальной подгонки осуществляется в симметризованном виде по каждой пространственной переменной в соответствии с одномерной схемой:

д / ь

дх V дх

- ехр

а ядя~

ЪяЬх

ая

- ехр

+ аьдь~

ЪьНх

аь

кх

Здесь введены следующие обозначения:

ая =

2 aiai+l +

ь=

2 а i аi-l а i + а^1

Ъ + Ъ

я

г+1

Ъь =

2 :

дя

к

-1

е-ь(х, дя

Ъяк

ая (1 - в-ьпк*/ап)'

-1

дь Кх у е-ь(х-хг,/а(1х

, дь =

ЪьКх

аь (еь^н*/а^ - 1)'

Аппроксимация производных по координате у двумерного диффузионно-конвективного уравнения производится аналогично.

Полученное разностное уравнение решается методом минимальных поправок [9].

2

н

х

0

0

2.3. Уравнения магнитной газодинамики

Уравнения (1.2) составляют систему дифференциальных уравнений в частных производных. Для решения этой системы был использован алгоритм, в основе которого лежит идея операторного расщепления исходной эволюционной задачи по физическим факторам [10], при этом распределение газодинамических параметров на новом временном слое определяется в два этапа.

На первом этапе решается подсистема уравнений Эйлера, в которой учитываются невязкие газодинамические потоки:

Ф-00 + + ^ = (2.2)

т дх ду

Здесь U — вектор газодинамических величин; т — шаг по времени, определяемый из условия Куранта:

min (hx, hy)

т <

max (Cij + |Uj |): ij

где с — местная скорость звука, |и| = л/и2 + т2.

Используя найденное распределение газодинамических параметров II1, на втором этапе решаем систему уравнений, в которую входят вязкие потоки:

12 -11 + дк (11) + да у (11) = о

т дх ду ' '

Подсистема Эйлера (2.2) решается с помощью явной ТУБ-схемы Хартена второго порядка точности [11], позволяющей получать квазимонотонное решение со вторым порядком аппроксимации во всей расчетной области и обеспечивающей хорошее разрешение контактных разрывов и ударных волн.

Для решения подсистемы уравнений газовой динамики, содержащей диффузионные слагаемые, запишем уравнения импульса и энергии в дивергентном виде:

^^ - а1у ((^ + у%9) уи) =1V ((^ + у%9) а^и), дЬ 3

де

— - а1у (кУе) = 0.

дЬ К '

Для построения разностной схемы применим проекционно-сеточный метод Галёрки-на с билинейными базисными функциями [12]. Решение уравнения представим в следующем виде:

lh _

(xi,yj )&w

uh _ E Ui,j (t) • фг,з (x,y),

СН = Е ег,з (Ь) • фг,з (х,у) ,

где

фг,з = фг(х) • Фз (У).

Для дискретизации диффузионного уравнения согласно методу конечных элементов умножаем уравнение на базисную функцию и интегрируем его по области . Для упрощения разностного уравнения используем метод концентрирующих операторов для билинейных базисных функций [13]. Разностная схема имеет семиточечный шаблон (рис. 2).

Полученная система разностных уравнений решается попеременно-треугольным методом:

^^- (ри)П + 1 (¿п+1/2из+1/2 + ци) = 0, (ри)П+1 +1 (ьп+1/2ига+1/2 + ьп+1ип+1) = о,

где Ь1 и Ь2 — нижняя и верхняя треугольные матрицы разностного оператора, аппроксимирующего диффузионные потоки.

Рис. 2. Разностный шаблон

Уравнение энергии, содержащее вязкие потоки, решалось по аналогичному алгоритму.

Уравнения, описывающие процесс переноса электронной плотности и турбулентной вязкости (1.3) и (1.5), имеют одинаковое математическое выражение и могут быть записаны в общем виде как нестационарное диффузионно-конвективное уравнение:

дв

— - (аУв - Ьв) = ¡,

где в обозначает V или пе соответственно, а(х,у) и Ь(х, у) — соответствующие диффузионный и конвективный коэффициенты, f — соответствующая правая часть.

Разностная аппроксимация последнего уравнения по пространственным переменным выполняется аналогично уравнению (2.1). Полученное в результате нестационарное разностное уравнение решается попеременно-треугольным методом.

3. Результаты численных расчетов

Численное моделирование двумерных процессов, возникающих в окрестности П-слоя из-за эффекта Холла, с учетом кинетики ионизационно-рекомбинационных процессов выполнялось в двух вариантах: в условиях узкого канала, когда толщина плазменного слоя 5 приблизительно равна расстоянию между электродами, т.е. И/5 = 1, и в условиях канала, близкого к дисковому, когда толщина плазменного слоя много меньше межэлектродного расстояния, в данном расчете оно принималось И/5 = 5. Остальные параметры режима МГД-генератора представлены в таблице.

Параметры генераторного режима

Параметры МГД-генератора Значение параметра

Рабочее тело — идеальный газ Неон

Температура торможения на входе Ts, К 2000

Давление торможения на входе Ps, МПа 1.0

Мах потока на входе Мо 1.5

Длина расчетной области Ь, м 0.4

Магнитное поле В, Тл 1

3.1. Вариант узкого МГД-канала (Н/8 < 1)

Ввиду сложности реального процесса в плазменном слое, что является результатом наложения многих эффектов (эффекта Холла, ионизационно-рекомбинационных и газодинамических эффектов), целесообразно было анализировать каждый из них отдельно, в отсутствие остальных. Кроме того, нас интересовали столь малые времена, что газодинамику можно было считать замороженной. Поэтому сначала рассматривалась постановка задачи без учета явлений газодинамики. Затем для уточнения роли газодинамических эффектов был рассмотрен вариант с учетом гидродинамики. Разность потенциалов между электродами составила 100 В, расстояние между ними — 10 см.

Из характера распределения эквипотенциальных линий, представленных на рис. 3, можно заключить, что эффект Холла приводит к формированию объемного заряда в диаметрально противоположных точках П-слоя. Эти области соответствуют концентрическим эквипотенциальным линиям в приэлектродных областях на рис. 3. ЭДС Холла обусловливает появление в приэлектродных областях отличной от нуля горизонтальной

Рис. 3. Распределение потенциала (р, В) в уз- Рис. 4. Распределение плотности тока в узком ком канале без учета газодинамики канале без учета газодинамики

Рис. 5. Распределение параметра Холла в уз- Рис. 6. Распределение потенциала в узком каком канале без учета газодинамики нале с учетом газодинамики

составляющей плотности тока — . В случае узкого канала эти приэлектродные зоны перекрываются, что приводит к наклону линий тока по всей их длине. Соответственно возникает наклон и всего плазменного слоя относительно стенок МГД-канала (рис. 4).

П-слой ограничен неэлектропроводным газом, что подавляет холловский ток в центральной части плазменного сгустка. Это приводит к выравниванию токового канала по всей области П-слоя, кроме приэлектродных областей. Здесь холловский ток может протекать вдоль электрода, что формирует в диаметрально противоположных приэлектродных областях П-слоя зоны концентрации линий тока, располагающихся косо-симметрично относительно ядра П-слоя. На рис. 5 зоны концентрации соответствуют областям пониженного значения параметра Холла. На рис. 6 представлены эквипотенциальные линии, полученные с учетом гидродинамических эффектов. Здесь симметричная картина распределения газодинамических параметров в пространстве между электродами, к поверхности которых происходит "прилипание" потока, накладывается на кососимметричную картину электродинамических характеристик. В результате общая картина распределения потенциала теряет симметрию.

3.2. Вариант широкого МГД-канала (Н/8 > 1)

Вариант широкого МГД-канала рассматривался без учета газовой динамики, разность потенциалов между электродами задавалась равной 500 В, расстояние между ними было 50 см. В случае, когда И/8 > 1, приэлектродные зоны с повышенным значением параметра Холла значительно удалены от центральной части П-слоя, в которой распределение линий тока имеет практически одномерный характер (рис. 7). Однако влияние приэлектродных областей с повышенной плотностью тока обусловливает сдвиг линий тока в зоне коммутации линий тока с ядром П-слоя. Это приводит к контракции разряда и постепенному стягиванию П-слоя в узкий токовый канал, т. е. к развитию перегревной неустойчивости, которая в данной модели канала постоянного сечения не подавляется охлаждением газа из-за расширения. Термин "контракция" (от лат. соПгасНо — сжатие, стягивание) распространен в физике газового разряда и означает уменьшение поперечного размера области, заполненной разрядным током.

0.3 -

.1е105

§

0.2

X, м

0.4

Рис. 7. Распределение плотности тока в широком канале без учета газодинамики

Таким образом, речь может идти об обнаружении нового физического явления — контракции неравновесного токового слоя в канале МГД-генератора под действием эффекта Холла. Было показано, что за время, значительно меньшее пролетного (0.1 мс), толщина П-слоя уменьшается более чем в 3 раза. Такая токовая пелена уже не может эффективно тормозить набегающий поток. Однако этот эффект нуждается в проверке в условиях расширяющегося канала.

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

Заключение

В заключение перечислим основные результаты данного исследования.

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

В "широком" канале с разведенными на значительное расстояние электродными стенками удается подавить эффекты Холла в ядре потока, но из-за приэлектродных явлений происходит процесс контракции разряда.

Эффект Холла в сочетании с гидродинамическими явлениями в неоднородном газоплазменном потоке вызывает нарушение симметрии распределения эквипотенциальных линий.

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

Список литературы

[1] СлАвин В.С. Неоднородный газоплазменный поток инертного газа в канале МГД-генератора jj ТВТ. 1998. Т. Зб, № 4. С. 647-654.

[2] СлАвин В.С., Данилов В.В., Кузоватов И.А. и др. Космические энергетические и транспортные системы, основанные на МГД-методе преобразования энергии // ТВТ. 2002. Т. 40, № 5. С. 810-825.

[3] Velikhov E.P. Hall instability of current carrying slightly ionized plasmas jj 1st Intern. Conf. of Magnetoplasmadynamic Electrical Energy Conversion. Newcastle-upon-Tyne, 1962.

[4] Гуляев А.Н. К созданию универсальной однопараметрической модели для турбулентной вязкости jj Изв. РАН. МЖГ. 1993. № 2. С. 69.

[5] Slavin V.S. Ionization kinetics in the nonuniform flow of noble gas carrying T-layers in the Faraday channel of MHD generator jj Proc. of 13th Intern. Conf. on MHD Power Generation and High Temperature Technologies. Beijing, China. 1999. Vol. 2. P. 539-550.

[6] Ильин А.М. Разностная схема для дифференциального уравнения с малым параметром при старшей производной j j Математические заметки. 1969. Т. 6, № 2. С. 237-248.

[7] Дулан Э. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

[8] Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск: Ин-т математики СО РАН, 2000. 345 с.

[9] Самарский А.А. Теория разностных схем. М.: Наука, 1989.

[10] Марчук Г.И. Методы расщепления. М.: Наука, 1988. 263 с.

[11] Harten A. High resolution schemes for hyperbolic conservation laws // J. of Comp. Physics. 1983. Vol. 49. P. 357-393.

[12] Флетчер К. Численные методы на основе метода Галеркина. М.: Мир, 1988. 352 с.

[13] Лаевский Ю.М. Проекционно-сеточные методы решения двумерных параболических уравнений. Новосибирск: ВЦ СО АН СССР, 1987. 138 с.

Поступила в редакцию 12 декабря 2006 г., в переработанном виде — 22 апреля 2007 г.