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

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

Том 11, часть 1, Специальный выпуск, 2006

ИСПОЛЬЗОВАНИЕ МЕТОДА КОНЕЧНЫХ

ОБЪЕМОВ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ГИПЕРЗВУКОВОГО МАГНИТОГАЗОДИНАМИЧЕСКОГО ТЕЧЕНИЯ

А.Б. Еремкин, Т.А. Коротаева, A.A. Маслов,

В.М. Фомин, A.n. Шашкин Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, Россия e-mail: korta@itam.nsc.ru, maslov@itam.nsc.ru fomin@itam.nsc.ru, shap@itam.nsc.ru

The problem of a three-dimensional supersonic gas flow around bodies embedded in an electromagnetic field imposes a number of requirements on the mathematical formulation and the method of solution. The computational (finite-volume) method proposed within the framework of three-dimensional Euler equations is rather simple for numerical implementation. A numerical scheme of the second order of approximation is used; the scheme is explicit in time and implicit in space. This scheme features a weak non-monotony and complete conservatism. Results of the numerical simulation of MHD-flow around a finite width wedge are presented as an example of application of the described numerical method.

Введение

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

Метод конечных объемов обладает свойством консервативности. При записи уравнений в консервативной форме и использовании дискретных преобразований, сохраняющих массу и другие величины, возможно получить решение, удовлетворяющее дивергентной форме исходных уравнений. Как показали Лаке и Вендрофф [1], решения уравнений, записанных в дивергентной форме, автоматически удовлетворяют условиям Рэнкина — Гнио-нио на любом скачке, который может возникнуть в потоке. Наиболее распространенным типом таких разрывов являются ударные волны. Следовательно, решение дискретных

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

уравнений автоматически улавливает поведение ударных волн, их интенсивность и скорость распространения в нестационарных течениях. Это позволяет использовать метод сквозного счета, т. е. проводить численное решение без предварительного выделения особенностей. Основная трудность применения методов сквозного счета состоит в получении резких профилей изменения параметров при переходе через скачок [2]. В частности, интерес к схемам, которые обеспечивают описание скачков уплотнения и других разрывов газодинамических параметров при отсутствии паразитных колебаний решения, обусловил широкое распространение схем минимизации полной вариации - TVD [2].

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

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

Разностная схема второго порядка как по времени, так и по пространственной переменной является неявной по пространству и явной по времени. В качестве примера использования метода конечных объемов в настоящей работе приведены результаты численного моделирования процессов, имеющих место при магнитогазодинамическом (МГД) воздействии на поток перед воздухозаборником [4].

Возникший в последнее время интерес к МГЛ-взаимодействиям в сверх- и гиперзвуковых потоках связан с возможностью использования МГД-систем для управления работой гиперзвуковых воздухозаборников и процессами в тракте сверхзвуковых прямоточных двигателей [5]. Этот способ управления потоком на входе в двигатель привлекателен, так как не требует изменения геометрии, в отличие от газодинамических способов регулирования воздухозаборников. Степень МГД-воздействия на поток зависит от значений параметров магнитного взаимодействия, Холла и коэффициента нагрузки, которые, в свою очередь, зависят от параметров потока, способа предыонизации и конкретной реализации МГД-системы.

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

В качестве простейшей модели воздухозаборника рассматривается клин конечной ширины. Электронный пучок расположен перед клином. Ток, возникающий в электронном пучке, направлен по нормали к набегающему на клин потоку. Предполагается, что пучок достаточно мощный, таким образом, он пересекает расчетную область насквозь, распределение плотности в поперечном сечении электронного пучка подчиняется распределению Гаусса. Два электрода расположены там же, где электронный пучок. Расположение электродов возможно в двух вариантах: вдоль пучка либо поперек, в плоскости симметрии клина. Магнитное поле направлено по нормали к плоскости симметрии клина, вдоль электронного пучка. Принята декартова система координат: ось X направлена вдоль набегающего потока, Y — то нормали к потоку и вдоль пучка, ось Z образует правую тройку. Схема задачи показана на рис. 1.

Рис. 1. Схема обтекания клина ионизированным воздухом (жирными линиями показаны электроды).

Исходными уравнениями являются трехмерные уравнения Эйлера. Система уравнений, описывающая движение повязкого, нестационарного, потоплопроводпого газа, дополнена 38 кинетическими уравнениями для 21-комнопептпой смеси, чтобы учесть физико-химические процессы, которые могут иметь место в потоке газа как результат воздействия на него электронного пучка. Учтены реакции пяти типов: химические реакции и реакции диссоциации; реакции ионизации при соударениях нейтральных частиц; реакции ионизации электронным ударом; реакции присоединения электрона и озона; возбуждение колебательных степеней свободы, электронным ударом. Список реакций и константы прямых и обратных реакций можно найти в |6|.

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

дЕ* ду

2п

Дге4 2

ту2

1п

Е*

где Е* — энергия электронного пучка; Zn — полное число электронов в единице объема;

ту2

0 — потенциал ионизации водорода;

— энергия электрона; е — его заряд. Затра-

2

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

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

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

Запишем исходную систему уравнений в терминах метода, который используется для получения численного решения. Пусть Ст(-) есть аналог оператора Гамильтона (V):

Ст (■) = V, (1)

где — направленный элемент поверхности элементарного объема V,

grad(ф) « Стф = Ст(ф), &у(Е) « Ст¥ = Ст ■ (Г), (w ■ ёга(1)Е « (w ■ Ст)¥.

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

Уравнение неразрывности для смеси в целом:

^ + Ст(р™) = 0. (2)

Уравнение сохранения количества движения:

д w 1

+ + -Ст(р) = Л х В. (3)

дЬ • р

Уравнение сохранения энергии:

дf I'2

-£ + СтШ-у, = +— + е1. (4)

Уравнение неразрывности для ¿-компоненты смеси:

Обобщенный закон Ома:

Л = а(Е + w х В) - О ((Л х В) + Стре)/1Б(6)

Здесь р, ю, р — плотность, скорость и давление газа соответственно; / — энтропийная функция, Ог = 1п(рг/Мгр); ^г = (Рг /Мгр), ^гг ~ раВНОВвСНОе ЗНаЧвНИв в реакции; NI — число реакций, в которых принимает участие ¿-й компонент; а — электрическая

ре В Л

тока; Е — напряженность электрического поля; аЕ — ток проводимости; а^ х В) — индуктивный ток; О — параметр Холла; I' — джоулево нагревание; £\ — энергия, полученная в результате химических реакций.

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

Для обезразмеривания приняты следующие параметры:

Ро = Р^, юо = = р^ю^,Ьо = Юо/Ьс,То =

/,0 = и}2, В0 = ЦеРоъи^, 30 = СГЮ0В0, Е0 = Ю0В0

(Ь — характерная длина). Здесь Ь — время; Т — температура; ре — магнитная проницаемость среды.

Уравнения решаются в пространстве физического потока (в декартовой системе координат). Обычный вид газодинамических граничных условий используется на всех границах. Условие непротекания принято на твердой стенке. На бесконечности параметры потока равны параметрам свободного потока. "Мягкие граничные условия" используются на внешних границах расчетной области. Граничные условия перед электронным пучком

Б

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

Зх = 0 = + О^ + Ву{Пи - гу),

для левой и правой границ

Ъ = 0 = % ~ + Ву(Пи! +

для нижней и верхней границ 0 —. На части границ, занятой электродами, р =

др ду

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

2. Численная реализация

Численное решение задачи проводится методом конечных объемов. Этот метод основан на представлении о том, что газодинамические параметры внутри элементарного объема являются среднеинтегральными и все изменения газодинамических параметров внутри ячейки определяются потоками через грани ячейки. Разбиение на элементарные объемы позволяет учитывать сложную конфигурацию области решения, не вводя криволинейную систему координат [11]. Рассматривается некоторая трехмерная область с границей, которая имеет, вообще говоря, изломы. Область разбита на ячейки, форма которых представлена в виде произвольных шестигранников. Для каждой расчетной ячейки уравнения газовой динамики записываются в интегральной форме (2)-(5). Обобщенная производная (1) может быть аппроксимирована соотношением

где к — номер грани г-й ячейки; — вектор площади к-й грани; — оператор осреднения по площадке. За номер ячейки, грани, ребра ячейки принимается номер узла с наименьшей суммой индексов.

Проводя дискретизацию в расчетной области, можно получить на множестве сеточных функций формулировку задачи (2)-(5):

д + С1 (¥д) = Н (д)

(7)

при условиях на границе Ф^(д) = 0 и начальных да иных дн = f (х1, х2, х3), где

( р \

Р и

V

ю

fs

V Р°г )

( рш + pvj + рюк \ (ри2 + p)i + рuvj + р^к рuvi + (рv2 + p)j + рvwk рию\. + рvwj + (рю2 + р)к

-Рfsui - Рfsvj - рЛюк \ -рСи - рСvj - рО^юк )

H(g)

0

(реЕ + Л х B)i (реЕ + Л х B)j (реЕ + Л х В)к

\

С^ т

\

ч 3 г

И

3 •

о

/

Подробный анализ используемой разностной схемы дан в [12], вде показано, что прпня-тая разностная схема имеет второй порядок аппроксимации и предназначена для расчета как стационарных, так и нестационарных трехмерных течений. Обсуждаются вопросы инвариантности, монотонности и консервативности схемы. Инвариантность и консервативность схемы — следствия принадлежности решения классу ограниченных функций (ВУ). Для усиления монотонности схемы используется экспоненциальный фильтр, устраняющий коротковолновые осцилляции в областях с большими градиентами параметров. Течение газа нестационарное, с заданным начальным полем параметров. Для получения стационарного решения применяется метод установления по времени. Алгоритм расчета химических процессов изложен в [13].

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

3. Расчеты обтекания клина потоком ионизированного газа в магнитном поле

Приведем основные параметры, которые позволяют оценить эффективность МГД-взаимо-действия.

Параметр магнитного взаимодействия (число Стюарта) N = оЕ21х/р^ю^ (1Х — характерный размер клина). Параметр магнитного взаимодействия определяет отношение магнитных сил к силам инерции.

д

р

Параметр Холла & = |В\а/(впе) (пе — количество электронов в единице объема). Этот параметр характеризует долю поперечного электрического ноля, которое возникает за счет отклонения электронов иод действием силы Лоренца. При больших значениях магнитной индукции в гиперзвуковом потоке плазмы роль эффекта Холла может быть значительной.

Коэффициент, нагрузки К = \Е|/(-ште|В|).

3.1. Влияние параметра магнитного взаимодействия

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

Полученные результаты представлены па рис. 2. Здесь показаны углы наклона головной ударной волны в плоскости симметрии в зависимости от значения параметра магнитного взаимодействия. Данные численного моделирования демонстрируют влияние электромагнитного ноля па положение головной ударной волны, генерируемой клипом. При N < 0.025 эффект магнитного взаимодействия слаб. Число Маха за областью магнитного взаимодействия уменьшается до М = 7.8 при N = 0.0028 и до М = 6.7 при N = 0.025. Со-

Параметры потока, электромагнитных нолей и геометрия тела

Обтекаемое тело

Симметри чный клин Ширина клина 0.1м Угол полураствора 15°

Параметры свободного потока

Число Маха М = 8 Высота полета 30 км Стати ческая температура 178.2 К

Параметры электронного пучка

Размеры поперечного сечения Напряжение электронной Плотность тока на выходе

(радиус) 0.05 м пушки 20 кВ 100 мА/см2

Параметры электромагнитного поля

Магнитная индукция 1. . 3 Тл Число Стюарта 0... 0.2 Коэффициент нагрузки К = 0

ответственно, угол наклона ударной волны увеличивается примерно на 1° при N = 0.0028 и примерно на 3° при N = 0.025. Положение головной ударной волны, полученное при N < 0.02, практически соответствует положению ударной волны при обтекании клина невозмущенным потоком. Если параметр магнитного взаимодействия растет, изменяется

N = 0.2

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

N

уп.иппения. В то же время наблюдается восстановление ударной волны вблизи кромок клина.

Рассчитанные в плоскости симметрии углы наклона скачка сопоставлены с данны-МП [14], полученными аналитически для одномерного случая и численно для плоского течения (рис. 2,6). Эти результаты согласуются между собой. Расхождение результатов тем более существенно, чем дальше расчетная точка от плоскости симметрии. Причины этого в пространственной форме области МГД-взаимодействия и влиянии концевых эффектов обтекания клипа конечной ширины. Исходя из этого можно заключить, что оценки угла наклона клина относительно малой ширины, полученные на основании одномерных или плоских аппроксимаций, могут быть недостаточно адекватными. Эффективность воздействия магнитного поля растет с увеличением параметра магнитного взаимодействия.

3.2. Электромагнитные силы и джоулево нагревание. Влияние параметра Холла

В работе предпринята попытка разделить джоулево нагревание и нетепловые эффекты для того, чтобы оценить, что приводит к перемещению ударной волны над клином. Изменение энергии или сил, которые могут быть результатом химических превращений, джо-улева нагревания и сил Лоренца и Ампера, учтено соответствующими членами в уравнениях движения и энергии. Таким образом, можно рассмотреть влияние тепловых и нетепловых эффектов раздельно, обнулив вклад соответствующих членов в уравнениях сохранения энергии или движения. Параметры, при которых проводились эти расчеты, приведены в таблице. Численные исследования проводились для условий, соответствующих высоте полета H = 47 км, так как было показано, что с уменьшением плотности эффективность МГД-воздействия возрастает, следовательно, роль тепловых и нетепловых эффектов может быть определена достаточно отчетливо. Выполнено три варианта расчетов.

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

2. С учет,ом, джоулева нагрева и, при, отсутствии действия электромагнитных сил. Этот случай показал, что при рассматриваемых условиях изменение наклона головного скачка за счет джоулева тепла невелико.

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

Рис. 3. Положение ударной волны.

Рис. 4. Течение газа в плоскости x = const при параметре Холла, не равном нулю.

Полученные результаты показаны на рис. 3. Кривая 1 — поверхность клина, кривая 2 — это положение ударной волны при обтекании клина однородным потоком, кривая 3 — положение ударной волны, когда обнулены члены, определяющие электромагнитные силы в уравнении движения. Положение ударной волны в магнитном поле (джоулево тепло обнулено) показано линией 4.

Электромагнитные силы приводят не только к изменению наклона и интенсивности ударных волн в гиперзвуковом потоке газа. При движении заряженного газа в магнитном поле возникают значительные силы Лоренца. На рис. 4 показано обтекание клина в одной из плоскостей x = const. Видно, что силы Лоренца и токи Холла приводят к смещению потока над поверхностью вправо и усиливают срыв потока за правым ребром клина. В потоке газа происходит активное образование вихрей.

Заключение

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

В работе представлены результаты численного моделирования влияния электромагнитного поля на обтекание клипа конечной ширины. Поток, обтекающий клин, предварительно ионизован электронным пучком. Результаты расчета показали, что в области электронного пучка рост температуры мал, так как толщина пучка мала и набегающий поток разрежен. Холодная плазма, образующаяся под воздействием электронного пучка, нагревается за счет торможения потока по мере движения над клином. Химические процессы не вносят какого-либо вклада в изменение положения головного скачка уплотнения. Наклон головного скачка над клином фактически определяется действием электромагнитных сил, во всяком случае, при тех условиях, при которых проводились расчеты (см. таблицу). Положение головной ударной волны может быть изменено при значениях параметра магнитного взаимодействия N > 0.01. При значениях магнитной индукции и плотности, соответствующих высоте полета 30 км, существенного МГД-влияния не наблюдается. Влияние токов Холла проявилось прежде всего в спрямлении линий плотности тока. Возникающие при движении заряженного газа в магнитном поле силы Лоренца и токи Холла приводят к несимметрии потока над клином и активному образованию вихрей.

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

[1] Lax P., Wendroff В. Systems of conservation laws // Commun. Pure and Appl. Math. 1960. Vol. 13, N 2. P. 217-237.

[2] Wang J.C.T., Widhopf G.F. A high-resolution TVD finite volume scheme for the Euler equations in conservation form // AIAA Pap. 1987. N 538. P. 1-17.

[3] Хемминг P.B. Численные методы: Пер. с англ. М.: Наука, 1968. 398 с.

[4] Fomin V.M., Maslov A.A, Malmuth N. et al. Numerical simulation of the MHD-effect upon flow around a finite-width wedge // ICMAR-2004 Proc. Pt I. 2004.

[5] Brichkin D.I., Kuranov A.L., Sheikin E.G. MHD-tecnology for scramjet control // AIAA Pap. 1998. N 98. P. 1642.

[6] MHD-control of a Flow Around a Wedge. Report ITAM ЕВ-2000. 2000. 62 p.

[7] mlgdal A.V. Qualitative Methods in Quantum Theory. M.: Fiz. Mat. Lit., 1975. 336 p.

[8] konovalov V.P., Son A.E. Degradation spectrums of electrons in gas // Khimiya Plasmy (Plasma Chemistry)/ Ed. by В.М. Smirnov. 1987. Vol. 14. P. 194-227.

[9] MHD-control of a Flow Around a Wedge. Report ITAM ЕВ-2001. 2001. 42 p.

[10] Shercliff J.A. A Textbook of Magnetohydrodynamics. N.Y.: Pergamon Press, 1965. 319 p.

[11] Chung T.J. Computational Fluid Dynamics. Cambridge: Cambridge Univ. Press, 2002.

[12] Fomin V.M., Maslov A.A., Korotaeva T.A., Shashkin A.P. Numerical simulation of a supersonic spatial nonuniform flow // CFD J. Special Issue. 2003. Vol. 12, N 2.

[13] Фомин B.M., Маслов A.A., Kopotaeba T.A. и др. Численные исследования обтекания клина конечной ширины низкотемпературной плазмой // Тр. Междунар. конф. RDAMM-2001. Новосибирск, 24-29 июня 2001. Т. 6, ч. 2. 2001 (спецвыпуск). С. 397-404.

[14] Malmuth N.D., Krivtsov V.M., Soloviev V.R. Quick, gridless estimations of MHD effects on hypersonic inlet ramp shocks //А1АА Pap. 2004. N 2004-0862.

Поступила в редакцию 12 апреля 2006 г.