Вычислительная модель радиационно-конвективного теплообмена в неоднородных магнитогазодинамических течениях Текст научной статьи по специальности «Физика»

Вычислительные технологии Том 10, № 6, 2005

ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ РАДИАЦИОННО-КОНВЕКТИВНОГО ТЕПЛООБМЕНА В НЕОДНОРОДНЫХ МАГНИТОГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЯХ

Е. Н. Васильев, Д. А. Нестеров Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: [email protected]

Non-stationary three-dimensional computation model of radiation-convection heat transfer in magnetogasdynamics flows with large gradients of physical parameters is introduced. The model is based on solution of equations of compressible, heat-conducting, inviscid gas dynamics combined with the system of equations of electrodynamics and equations of radiative transport. Modeling of a process of the MHD interaction between conducting layer with electric current and the gas flow revealed that a vortex-type flow is developed in a channel. Evolution of the Rayleigh — Taylor hydrodynamic instability is the cause of a partition of the current layer into two or more parts.

Введение

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

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

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

1. Математическая модель

1.1. Система уравнений газовой динамики

Рассмотрим нестационарное течение невязкого теплопроводного излучающего газа в канале постоянного прямоугольного сечения. Моделирование процесса проводится для плотной газовой среды при давлении, близком к атмосферному, и выше, поэтому для описания применимо магнитогидродинамическое приближение. Кроме того, в сильноточных разрядах выполняется условие термодинамического равновесия [1], при котором степень ионизации, электропроводность и другие свойства плазмы могут быть однозначно выражены через термодинамические параметры среды, например температуру и давление. В этом случае течение описывается системой уравнений газовой динамики в эйлеровых координатах (1)-(4), дополненной уравнениями состояния газа (5):

дИ + дЕ + д¥ + дО _ 8 д£ дх ду дг

(1)

и

р ри рv

ри ри2 + р рvм

рv , Е _ рмv , Г _ рv2 + р

ри рии рvw

. Ег _ _ (Ег + р)и + Уж _ _ (Ег + p)v + Уу _

Е _ р

О

и2 + v2 + и2

ри 0

рии /х

ри^ , й _ /у

ри2 + р /г

(Ег + р)и + ^ _ ф

2

+ е , ф _ QJ - фд + /ли + /уV + /и;

-Лдг

дх'

Уу

-Лдт

-Лдт;

дг'

р _ р(р,е), Т _ Т(р,е)-

(2)

(3)

(4)

(5)

Здесь р — плотность газа; — компоненты вектора скорости газа V; р — давление;

Ег — полная энергия единицы объема газа; е — внутренняя энергия единицы массы газа; Т — температура газа; /ж, /у, /г — компоненты вектора объемной силы ^; QJ — объемная мощность джоулевой диссипации, фд — изменение энергии за счет радиационного теплообмена; Л — коэффициент теплопроводности газа. Уравнения состояния газа (5) могут задаваться как в виде аналитических выражений, так и в виде табличных зависимостей.

Начальные условия для системы уравнений (1)-(5) задаются в виде распределений температуры, давления и скорости:

Т (х,У,г,^)|г=о _ ТоОх У Р(х,У,г,^)|г=о _ ро(х,У,г) у ^ ^ ^=о _ Граничные условия определяют на входе в канал параметры сверхзвукового потока:

Т(х,у,г,^)|ж=0 _ р(х,у,г,^)|ж=0 _ р^У^^ ^ЗД^^о _

У

У

ж

г

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

dT(x, y, z, t)

0 dp(x,y,z,t)

x=L dx

dx

На стенках канала задаются условия непроницаемости.

0 д v(x,y,z,t)

x=L дх

= 0.

x=L

1.2. Уравнения радиационного теплообмена

Для расчета изменения энергии единичного объема газа за счет радиационного переноса излучения решаются следующие уравнения:

П • grad(Iv) = ки(I„p - Iv); (6)

oo

W = J dv j a IvdQ; (7)

0

Qr = div(W). (8)

Здесь Iv — интенсивность энергии излучения частоты v; Ivp = 2hv3 /[c2(exp{hv / (kT)} — 1)] — интенсивность равновесного излучения; kv (v, T,p) — коэффициент поглощения излучения; a — единичный вектор, определяющий направление излучения для угла dQ; W — вектор потока энергии излучения.

Граничные условия для уравнения переноса излучения (6) определяются падающим извне излучением: Iv|г = I0v(a, и)|г, где Г — граница канала; a — вектор, определяющий направление излучения.

1.3. Система уравнений электродинамики

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

Eind = ~dtr - "V - 2ПТ (9)

с величиной характерного электрического поля в разряде, возникающего при движении газа в магнитном поле, v х B = uBext. Здесь для расчета величины индуцированного магнитного поля Bind применена формула для бесконечного прямолинейного проводника. В оценке используются характерные для процесса величины: r — радиус столба разряда, т = 10-4 с — время перестройки структуры разряда за счет силового взаимодействия, u = 103 м/с — скорость движения разряда, Bext = 3 Тл — величина внешнего магнитного поля, I = 104 А — разрядный ток, ^0 = 4п • 10-7 Гн/м — магнитная постоянная. Таким образом, отношение сравниваемых величин много меньше единицы:

- ЩИ - ю-2 < 1, (10)

UBext TUBext

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

j = a(Ee + Eext); (11)

rot Ee = 0; (12)

div Ee = 4npe- (13)

Здесь Ee(x,y,z) — напряженность электрического поля, определяемая внешней нагрузкой; Eext(x, y, z) — напряженность индуцированного электрического поля; j(x, y, z) — поле плотности тока; a(x,y,z) — электропроводность газа; pe(x,y,z) — поле плотности электрических зарядов.

Значение Eext определяется исходя из величины постоянного внешнего магнитного поля B = (0,By, 0) и скорости газа v, при этом индуцированное магнитное поле в силу малости в данной модели не учитывается

Eext = v х B. (14)

Используя скалярную функцию потенциала электрического поля <^(x,y, z), удовлетворяющую уравнению

-Vp = Ee, (15)

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

div(aV^) = div(aEext). (16)

Решение данного уравнения позволяет найти распределение потенциала y, z), по которому находятся электрическое поле Ee из (15) и распределение плотности токов j(x,y, z) с помощью (11). Величины джоулевой диссипации и объемной силы Лоренца, действующей на электропроводный газ, рассчитываются следующим образом:

Qj = j2/a, f = j х B, (17)

при этом не учитывается электростатическая сила peEe, поскольку она пренебрежимо мала по сравнению с силой Лоренца. В качестве граничных условий уравнения (16) на электродах задаются значения потенциала. На непроводящих границах, где отсутствует нормальная составляющая плотности тока, приравнивается к нулю соответствующая производная потенциала.

2. Вычислительный алгоритм

Все переменные в уравнениях (11)—(17) определены в прямоугольной области х £ [0,Ь], у £ [0, Ш], г £ [0,Н], где Ь, Ш и Н — длина, ширина и высота канала соответственно. Для численного решения уравнений вводится равномерная прямоугольная сетка

, Ну, Кг) _ (х* _ (г - 1/2)Кж, г _ 1, 2 ...,

у = а - 1/2)ку }3 = 1, 2, Ыу = (к - 1/2)кг ,к = 1, 2,...,^} (18)

с шагами кх = Ь/Ых, ку = Ш/Ыу, кх = Н/Мг и количеством узлов Ыу и вдоль осей х, у и г соответственно. Для задания граничных условий в некоторых случаях вводятся дополнительные узлы, находящиеся за границей расчетной области, с номерами 0 и N + 1.

Численное решение системы уравнений (1)-(5) на каждом временном шаге представляет собой три этапа:

1) по значениям температуры и давления с текущего временного слоя в каждой точке сетки из решения уравнений (6)-(8) находится распределение Qя;

2) из решения электродинамических уравнений (9)-(15) находятся распределение плотности токов и напряженность электрического поля, по которым в каждой точке расчетной области вычисляются величины объемной электродинамической силы ^ и джоулевой диссипации QJ;

3) исходя из критерия Куранта определяется величина шага по времени и решаются газодинамические уравнения (1)-(5).

2.1. Решение уравнений радиационного теплообмена

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

(6). Это уравнение численно решается с помощью метода характеристик (£п-метод) [2, 3], основанного на определении набора характеристических направлений, для каждого из которых решается одномерное уравнение переноса. При этом каждому характеристическому направлению в зависимости от пространственной конфигурации задачи соответствует свой телесный угол П8, в пределах которого интенсивность излучения считается постоянной 13. Интеграл по телесному углу в уравнении (7) / а принимает вид ^ (/8Ъ8), где

Ъ8 = / а^^. Векторы Ъ8 зависят только от выбора характеристических направлений, по-

п3

этому их целесообразно вычислять перед началом расчета. Таким образом, уравнения (6),

(7) сводятся к соотношениям

d^

^Г = К(1„р - /„); (19)

ds

18 = ^2 W = 5] Ъ818, Ъ8 = / asdtt. (20)

- 8 о,

Здесь s — характеристическое направление, вдоль которого рассчитывается интенсивность излучения; и — интенсивность излучения группы частот п и интегральная интенсивность по всем частотам в направлении s; — телесный угол, соответствующий направлению s; а8 — единичный вектор, определяющий направление s; Аи — интервал частот для группы п.

В соотношениях (18) суммирование проводится по всем группам частот и по всем направлениям. Порядок задания направлений определяется размерностью задачи и соотношением шагов пространственной сетки кх, ку, кх. В одномерном случае рассматриваются направления, ориентированные к оси х под углом в, как правило, с равным шагом по величине сое в. Для двумерной постановки пространственная дискретизация показана на рис. 1. В плоскости ху выбираются восемь направлений, согласованных с пространственной сеткой. В плоскостях, проходящих через ось г и каждое из восьми направлений, заданы еще

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

По каждому из направлений в находятся интенсивности излучения с помощью уравнения (17), для этого используется разностная схема:

I = + (1{-1 - 1р-1) • ехр(-1) - - 1рг-1)(1 - ехр(-/))//, (21)

где 1 = Н8(к— + к)/2 — оптическая толщина; Н8 — расстояние между соседними узлами в направлении в; к — коэффициент поглощения в рассматриваемом узле.

1 Тгж'

/ ч / ч ж\¡1/1+1,1 X

/ - - / 1 / А/ / / ¿Г 7 / АУ Ш \У/ 7 /

у/ /Ч У+и+1

Направление луча ^

(направление счета)

Рис. 1. Пространственная дискретизация при решении уравнения переноса излучения в двухмерном случае.

Рис. 2. Пространственная дискретизация при решении уравнения переноса излучения в трехмерном случае.

1

)

4— —^

] 10 20 30 : 0 60 70 80 90 100 110 120 130 140 1Е 0 170 180 190

1 1 Узлы сетки

1

-Решение уравнения переноса — Объемный излучатель

Рис. 3. Распределение источников и стоков энергии, обусловленных излучением.

В случае большого количества групп по частотам решение уравнений (17), (18) требует значительных вычислительных затрат. Для их сокращения эти уравнения следует решать только в горячей области потока и ближайшей окрестности, где излучение оказывает реальное влияние на течение процесса. Для иллюстрации особенностей радиационного теплообмена на рис. 3 представлены результаты решения уравнения переноса излучения для плоского слоя воздуха при давлении 105 Па. Слой имеет в центре однородную область горячего излучающего газа толщиной 10 см и с температурой 104 К, а на его краях находится относительно холодный газ с температурой 2 • 103 К. Также на рис. 3 для сравнения приведено распределение Qя, полученное в приближении объемного излучателя по коэффициентам черноты плоского излучающего слоя толщиной 10 см и с температурой 104 К. Значения коэффициентов поглощения к{у,Т,р) для 560 спектральных групп в интервале 250 ... 140 000 см-1 и коэффициентов черноты е(8, Т,р) взяты из [4] и введены в программу в виде таблиц.

Из рис. 3 видно, что для радиационного теплопереноса характерны повышенная мощность энергопотерь на краях горячей области и интенсивное поглощение излучения холодным газом в ее ближайшей окрестности. В центральной части горячего слоя распределение Qя практически однородно и близко по величине к значению, полученному по упрощенной модели объемного излучателя. Результаты данного тестового расчета свидетельствуют о том, что для корректного расчета параметров неоднородных высокотемпературных газовых потоков необходимо решать уравнение переноса излучения.

2.2. Решение электродинамических уравнений

Для решения электродинамического уравнения (14) использован метод установления, когда из соответствующего нестационарного уравнения

= а1у(аУ^) - а1у(аЕех4) (22)

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

В качестве начального условия для уравнения (??) используется установившееся распределение потенциала на предыдущем временном шаге газодинамической задачи, а в начале моделирования задается равномерное распределение потенциала между электродами

tfo(x, y, z) = ¥>i(1 - z/H) + tf 2(z/H), (23)

где H — расстояние между электродами; tf и tf2 — значения потенциала на электродах.

Для магнитогазодинамических течений с Т-слоями характерны резкие градиенты температуры, давления, скорости, при этом величины а и Eext могут отличаться в соседних точках на порядки. Как показали тестовые расчеты, в таких условиях использование значений а и Eext в "полуцелых" точках сетки или применение центральных разностей приводят к значительному искажению получаемых величин токов в точках, близких к градиентам параметров газа. И, как следствие, ошибочно определяются силы, действующие на газ, и величина джоулевой диссипации. Это может приводить к сильному "нефизическому" нагреву отдельных точек и "развалу" численного решения.

Наиболее корректное численное решение уравнения (??) получено при использовании неявных локально-одномерных схем следующего вида:

= hx Mtf+1 - tf**) - *-1(¥>Г - tf-1)) - ^ (* (EXxt) - а-1 (EXxtLi) ; (24) = ц (a+1(tf+1 - tf**) - a(tf** - tf-1)) - hX (am (EXxt)m - а (ET),) • (25)

Здесь tf * и tf * * — распределение потенциала до и после действия разностного оператора соответственно. В выражениях (??), (??) для краткости записи индексы j и k опущены. Разностные аппроксимации, записанные для всех i, совместно с граничными условиями образуют систему алгебраических уравнений, которую решают методом прогонки. Обозначим процедуру определения tf * * для оси x для всех j и k по соотношениям (??) как LX-, по соотношениям (??) — как LX+, а по осям y и z — соответственно Ly±, LZ±. Выполнение шага тп происходит при последовательном выполнении процедур по отдельным направлениям tfra+1 = LX±Ly±LZ±tfn, причем в ходе расчета для достижения симметрии поочередно используются все 8 возможных комбинаций: Lex~Ly-Lez~, Lex+ Ley~Lez~, Lex+ Ley+Lez~ и т.д.

Напряженность электрического поля определяется следующим образом. Если в направлении оси x использовался разностный оператор LX-, то компонента вектора (EX)ij fc = (tfi,j,k - tfi-1,j,k)/hx, если использовался оператор LX+, то (EX)ijjk = (tfi+1,j,k - tfi,j,k)/hx. Аналогично находятся компоненты векторов Ey и EZ для направлений y и z соответственно.

Сходимость решения электродинамической задачи контролируется по разнице в величине z — составляющей суммарного тока J в разных сечениях xy электропроводящей области. Итерации прекращаются, когда максимальное расхождение в значениях Jz для всех сечений становится меньше заданного ez.

Для иллюстрации работы алгоритма на рис. 4 приведено поле плотности тока в расчетной области 20 х 20 х 20 узлов, в ее центре расположена зона в форме куба размером 10 х 10 х 10 узлов с электропроводностью а1 = 2800 Ом-1м-1 (T1 = 104 К); за ее пределами электропроводность а2 = 89 Ом-1м-1 (T2 = 6 • 103 К). Шаг сетки hx = hy = hz = 4 мм, разность потенциалов Atf = 60 В.

Рис. 4. Поле плотности токов в среднем сечении уг, серая область соответствует области повышенной электропроводности.

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

2.3. Решение газодинамических уравнений

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

И. = их(т )Ц (тЩ (тЩ (т Щ (т )Щт . (26)

Здесь т — шаг по времени, определяемый из критерия Куранта на каждом временном слое. В двумерной постановке последовательность (23) содержит операторы Ь9х(т) и (т), а в одномерной — только Ь9(т).

Разностный оператор определен по двухшаговой схеме типа "предиктор-корректор":

— т

иг.,к = Иг.,к — (Ег+1,*,к — Ег,*,к) + т; (27)

^.9

и**

г.,к

и* ТТ** т (|.Л ** рл** \ , ох

г,*,к + Иг.,к 7 (Ег,*,к Ег-1,*,к) + т Бг,*,к П.

(28)

Здесь и* * к — значение вектора до действия оператора; и**к — значение вектора после действия оператора; Б.* к — правая часть уравнения для локально-одномерной схемы вдоль оси х. Для операторов Ьу(т) и Ьх(т) выражение имеет аналогичный вид, за исключением того, что вместо Б9 используются Бу и соответственно. Расщепление в

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

0

fx 0 0

_ QjR + fxU _

, S

y —

0

0

fy

0

QjR + fy v

0 0 0

fz

QjR + fzw _

, Qjr = (Qj - Qr) /3. (29)

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

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

Течение содержит области с большими градиентами физических параметров (давления, температуры, скорости), обусловленных наличием Т-слоев и ударных волн, поэтому для устранения осцилляций и увеличения точности расчетов используется метод коррекции потоков FCT (Flux-Corrected Transport) [7]. Этот метод является дополнительным шагом, который необходимо выполнять после применения каждого одномерного оператора в последовательности (23).

Свойства газодинамического алгоритма были исследованы на тестовых задачах.

Течение возле переднего уступа. Рассматривается двумерная задача обтекания уступа сверхзвуковым (M = 3) потоком в канале единичной ширины и длины L = 3. Уступ высотой l = 0.2 расположен на расстоянии d = 0.6 от входа в канал. На входе в канал заданы параметры сверхзвукового потока идеального газа с показателем адиабаты y = 1.4, на выходе применяются "мягкие" граничные условия. Начальные условия заданы в безразмерных переменных: скорость u = 3, p =1, р =1.4. Для оценки погрешности вычислительного алгоритма результаты моделирования приведены для двух последовательных сеток, имеющих соответственно 120 х 40 и 240 х 80 интервалов. На рис. 5 показаны изолинии плотности для момента времени t = 4. Как видно из рисунка, алгоритм позволяет хорошо рассчитывать течения с ударными волнами, фронт которых распределяется на два разностных интервала, при этом наблюдается достаточно хорошее соответствие положения скачков в сверхзвуковом потоке на сетках с различной величиной пространственного шага.

Обтекание цилиндра. Рассматривалась двумерная задача обтекания цилиндра дозвуковым потоком с числом Маха = 0.45 и давлением p= 105 Па. В качестве газа рассматривался идеальный газ с показателем адиабаты y = 1.4, молярной массой ^ = 28.9 г/моль. На боковых стенках канала и цилиндре задавались граничные условия на непроницаемой стенке, на входе задавались параметры невозмущенного потока, на выходе — "мягкие" условия. Результаты расчетов также представлены для двух сеток: 85 х 150 узлов и 170 х 300 узлов с пространственным шагом 5 и 2.5 мм; диаметр цилиндра d = 80 мм (рис. 6, а, б).

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

x

z

S

S

б

Рис. 5. Изолинии плотности: 25 значений плотности для изолиний от 0.6 до 6. Количество узлов расчетной сетки: а — 120 х 40, б — 240 х 80.

Рис. 6. Вихревое течение, формирующееся при обтекании цилиндра: а, б — расчет, в — эксперимент.

T « 3.2 • 10-3 с, который соответствует числу Струхаля Sh = d/TU^ = 0.17. Такое течение, называемое вихревой дорожкой Кармана, наблюдается экспериментально (рис. 6, в) при значении числа Рейнольдса Re > 40 [8]. В экспериментах при подобных условиях Sh = 0.18.

Результаты теста свидетельствуют об адекватном описании обтекания цилиндра газовым потоком. Факт формирования вихревой структуры потока интересен тем, что получен при моделировании невязкого течения, основанном на решении уравнений Эйлера. О возможности моделирования турбулентных течений на основе решения уравнений Эйлера известно из работ О.М. Белоцерковского [9], в которых этот подход развивается на основе метода крупных частиц.

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

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

= N

pm M am,

i=1

где N — количество спектральных интервалов; M — количество участвующих в расчете компьютеров; am — коэффициент, отражающий производительность компьютера m. Частотный диапазон и распределения температуры и давления передаются вспомогательным расчетным машинам, которые после получения задания начинают расчет. Основной компьютер в это время решает уравнения электродинамики и уравнения переноса для своего диапазона. После окончания счета полученные решения со вспомогательных компьютеров передаются основному, на котором проводится интегрирование результата по всему спектру. Далее на основном компьютере решаются уравнения газовой динамики.

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

Обмен данными между компьютерами происходит по сети с помощью общераспространенного протокола TCP/IP. На вспомогательных машинах расчетная программа (сервер) функционирует под операционной системой Win 2000/XP в виде сервиса.

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

3. Модель теплообмена в канале МГД-генератора

Представленная вычислительная модель была применена для моделирования МГД-взаи-модействия токового слоя с потоком неэлектропроводного газа в генераторном режиме. Для исследования процессов теплообмена в МГД-канале на основе описанного выше алгоритма была разработана и реализована трехмерная нестационарная вычислительная модель. Схема моделируемого процесса представлена на рис. 7.

В канал прямоугольного сечения со сверхзвуковой скоростью поступает относительно холодный непроводящий газ. Верхняя и нижняя стенки канала представляют собой электроды, замкнутые на нагрузку, боковые стенки — диэлектрики. В движущемся газе за счет энергии внешнего источника проводится инициирование токового слоя, т. е. локальная область газа нагревается до температуры порядка 104 К, что приводит к термической ионизации газа. В результате взаимодействия движущегося электропроводного газа с поперечным магнитным полем через контур Т-слой — электроды — нагрузка течет ток, обеспечивающий нагрев Т-слоя и компенсирующий тепловые потери. На проводящую область газа действует тормозящая сила Лоренца. Таким образом, образуется самоподдерживающийся токовый слой, который, частично перекрывая канал, выступает в роли плазменного поршня и позволяет преобразовывать энергию движущегося газа в электрическую.

В данной вычислительной модели в качестве рабочего тела использовался воздух. Уравнения (4), (5) задавались в приближении идеального газа с показателем адиабаты 7 =1.3 и молярной массой ^ = 28.9 г/моль. Коэффициенты электро- и теплопроводности задавались в виде таблиц а (Т,р) и А(Т,р) [10]. Значения коэффициентов поглощения ку (V, Т, р) взяты из [4]. Величина индукции внешнего магнитного поля равнялась В = 3 Тл, сопротивление нагрузки 0.05 Ом.

В начальном состоянии газодинамическое течение, имеющее температуру 103 К, включает в себя разрядную область, заданную в виде правильного цилиндра с радиусом г = 3 см и температурой Т = 104 К. Остальные параметры течения имеют однородное распределение: давление р = 0.1 МПа, компоненты скорости и = 2 • 103 м/с, V = 0, w = 0.

Взаимодействие с неэлектропроводным потоком и внешним магнитным полем приводит к деформации разряда и его расширению в поперечном направлении при сохранении однородности распределения параметров в вертикальном направлении до момента времени £ ^ 0.1 мс. Обтекание разрядной области холодным газом сопровождается прогревом тонкого слоя газа, прилегающего к разряду, и уносом образовавшейся горячей пелены вниз по течению. Кроме того, с токового слоя обтекающим потоком срываются и более крупные области горячего газа, которые, попадая в спутный след, образуют там круп-

Рис. 7. Схема моделируемого процесса: 1 — сверхзуковой поток газа, 2 — Т-слой, 3 — электроды, 4 — нагрузка.

Рис. 8. Распределение температуры в среднем сечении МГД-канала в момент времени Ь = 0.3 мс (на шкале приведено соответствие оттенков серого цвета температуре в 103 К).

Рис. 9. Поверхности постоянной температуры (Т = 7000 К), соответствующие моментам времени: а — 0.3 мс, б — 0.4 мс, в — 0.6 мс, г — 1.0 мс (по осям графиков отображены номера расчетных точек).

Рис. 10. Векторное поле скорости в МГД-канале в системе отсчета, связанной с набегающим потоком (момент времени £ = 0.4 мс).

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

Для границ токового слоя характерен большой перепад значений температуры и плотности, при этом на верхней по потоку границе наличие многократного перепада плотности сочетается с действием силы газодинамического давления со стороны более тяжелого "холодного" газа. Сочетание этих физических факторов является условием для развития на этой границе неустойчивости Рэлея — Тейлора, действие которой проявляется в том, что "языки" холодного газа вклиниваются в область разряда, разделяя ее на две части (рис. 9, в). Более крупная разрядная область имеет большие значения температуры и электропроводности и сильнее тормозится в магнитном поле. Вторая, меньшая, область сносится ниже по течению и попадает в вихревой спутный след первого разряда и из-за мощных конвективных и радиационных потерь с течением времени остывает и теряет электропроводность. Оставшийся разряд восстанавливает свое первоначальное (до развития неустойчивости) поперечное сечение благодаря увеличившейся джоулевой диссипации и прогреву ближайшего к разряду слоя холодного газа механизмами излучения и теплопроводности (рис. 9, г). Далее в данном режиме процесс развития неустойчивости, разделения разрядной области на две с последующим остыванием одной из них периодически повторяется.

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

Заключение

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

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

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

[1] Райзер Ю.П. Физика газового разряда. М.: Наука, 1987. 592 с.

[2] Четверушкин Б.Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1985. 304 с.

[3] Адрианов В.Н. Сеточный метод исследования радиационного и сложного теплообмена // Изв. АН СССР. Энергетика и транспорт. 1988. № 2. С. 142-150.

[4] Авилова И.В., Биверман Л.М. и др. Оптические свойства горячего воздуха. М.: Наука, 1970. 320 с.

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

[6] Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. М.: Мир, 1990. 392 с.

[7] Book D.L., Boris J.P., Hain K. Flux-Corrected Transport II. Generalization of the method // J. of Comput. Phys. 1975. Vol. 18. P. 248-283.

[8] Van Dyke M. Album of fluid motion. Stanford, California (USA): Parabolic Press, 1982.

[9] Белоцерковский О.М., Опарин А.М., Чечеткин В.М. Турбулентность: новые подходы. М.: Наука, 2003. 286 с.

[10] Соколова И.А. Коэффициенты переноса воздуха в области температур от 3000 до 25000 К и давлений 0.1, 1, 10, 100 атм // ПМТФ. 1973. Т. 13, № 2. С. 80-90.

[11] Sebald N. Measurement of the temperature and flow fields of the magnetically stabilized cross-flow N2 arcs // Appl. Phys. 1980. Vol. 21. P. 221-236.

Поступила в редакцию 11 апреля 2005 г., в переработанном виде — 23 июня 2005 г.