Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 0,940

MODELLING OF PARAMETERS OF THE PLASMA PISTON AT THE FRONT SHOCK WAVE OF THE COAXIAL DEVICE

Vasileva O.V. 1
1 National Research Tomsk Polytechnic University
Work is devoted to the solution of the one-dimensional equation of a hydraulic gas dynamics for the coaxial device - magneto plasma accelerator by means of Lax-Wendroff modified algorithm with an optimum choice of parameter of regularization – artificial viscosity. Replacement of the differential equations in private derivatives is made by final differences. Optimum parameter of regularization – artificial viscosity in the environment of MathCAD is added, using the exact known decision – a task Soda. The developed algorithm of calculation of thermodynamic parameters in a braking point is approved. On the basis of the offered algorithm in the environment of MathCAD thermodynamic parameters of a shock wave in front of the plasma piston are calculated at its departure from the coaxial magneto plasma accelerator. When modeling overwhelming unstable high-frequency fluctuations that allows to narrow area of heterogeneity are considered and to allocate only smooth decisions. Results of calculation of gas dynamic parameters in a point of braking coincide with literary data.
fluid dynamics.
artificial viscosity
shock wave
mathematical modeling
Magneto plasma accelerator

Введение

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

Методика

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

. (1)

Здесь ρ – плотность газа; p – давление газа; u – скорость распространения газа; ε – внутренняя энергия газа; t,x – время и координата.

Те же уравнения (1), записанные в векторной форме, удобной для численной реализации:

, (2)

где s – вектор консервативных переменных; f(s) – вектор потока.

Описание алгоритма:

1. Замена дифференциальных уравнений в частных производных конечными разностями (2).

2. Добавление оптимального параметра регуляризации – искусственной вязкости в среде MathCAD.

3. Выбор оптимальной искусственной вязкости, используя точное известное решение (задача Сода), рис. 1, б.

4. Апробация разработанного алгоритма расчета термодинамических параметров в точке торможения (аналитически), табл.

5. Расчет динамики изменения термодинамических параметров перед плазменным поршнем при вылете из ускорителя, рис. 5.

Для численного решения системы уравнений использовался модифицированный алгоритм Лакса–Уэндроффа [1], который заключается в том, что уравнения в частных производных заменяются конечными разностями. В конечных разностях появляется неустойчивость, т.е. влияние высокочастотных шумов из-за наличия сильных ударных волн. К данному алгоритму нами была добавлена искусственная вязкость. В режиме on-line был подобран оптимальный параметр регуляризации (искусственной вязкости) для обеспечения регуляризации решения и подавления его шумовой составляющей. Формируется массив значений для каждого слоя, используя среду MathCAD:

,

где LW – функция Лакса–Уэндроффа; N=200, M=200 – число точек разбиения пространственного и временного интервалов соответственно, h, τ – шаг по пространству и времени соответственно, γ=5/3 – показатель политропы, µ – искусственная вязкость, s, f – вспомогательные функции для формирования массива значений.

Величина искусственной вязкости определялась из невязки-рассогласования (рис. 1, а). Точное известное решение задачи Сода сравнивалось с нашим алгоритмом, если рассогласования составляли не более 10 %, то коэффициент регуляризации (искусственная вязкость) считался оптимальным, рис. 1, б.

Экспериментальная часть

Результат расчета динамики изменения термодинамических параметров в относительных единицах приведен на рис. 2. До столкновения волн величины давления, плотности и температуры в средах были одинаковы, а скорости волн одинаковые по величине, но разные по направлению (знаку). На рисунке приведен момент времени после столкновения двух ударных волн при числе Маха, равном 1,5. Осью симметрии рисунков является фронт плазменного поршня.

Проверка работы алгоритма была проведена на расчете критических параметров – давления, плотности и температуры торможения pТ, ρТ, TТ при заданных начальных данных невозмущенного газа p0=1/γ, ρ0=1, T0=1 (в относительных единицах) и заданном числе Маха M0, рис. 4, а. Параметры торможения вычислялись по известным соотношениям [2, 8]. Давление, плотность и температура из невозмущенной среды, через скачок уплотнения, пересчитываются через ударную адиабату Гюгонио-Ренкина (ударная волна):

, ,

.

а  

б  

Рис. 1. Восстановленные термодинамические параметры ударной волны и волны разряжения: а) с неоптимальной вязкостью µ=10–7; б) с оптимальной вязкостью µ=10–6


                 а                                б                                     в                             г

Рис. 2. Случай двух ударных волн: а) плотность ρ; б) давление p; в) скорость v; г) энергия E

аб

Рис. 3. Параметры ударной волны: а) критические параметры; б) распределение плотности на фронте волны при различных числах Маха M0

В нашем случае эти формулы вырождаются в простые удобные для инженерных расчетов соотношения: , , .

Для расчета параметров торможения нужно использовать адиабату Пуассона (волна разрежения): , , , здесь M1 – числа Маха после ударной волны, определяется через число Маха в невозмущенной среде выражением:. Для инженерных расчетов получаем:

, . (3)

Начальные данные приведем к относительным единицам p0=1/γ=0,6, ρ0=1, T0=1.

Результаты сравнения методов сведем в таблицу.

Таблица. Результаты сравнения инженерного (И) и программного (П) расчетов

Число

Маха M0

Расчет

1,5

1,536

2,280

1,714

2,172

1,495

1,750

И

1,542

2,287

1,717

2,175

1,497

1,752

П

2

2,850

3,807

2,286

2,719

2,078

2,333

И

2,914

3,893

2,306

2,744

2,106

2,364

П

3

6,600

8,204

3,000

3,418

3,667

4,000

И

6,557

8,151

2,995

3,413

3,649

3,980

П

5

18,600

22,300

3,571

3,982

8,680

9,333

И

18,522

22,206

3,570

3,980

8,647

9,298

П

10

74,850

88,397

3,883

4,291

32,123

34,333

И

75,050

88,643

3,884

4,291

32,210

34,426

П

Приведем расстояние между ударной волной и торцом поршня S, используя формулу Лунева. Выпишем коэффициент в выражении (3):

, .

Как видно из табличных данных, программа дает хорошие результаты. Графическое изображение скачков уплотнения приведены на рис. 3, б. Торец поршня находится в точке x=0,5 [7]. С помощью численной схемы Лакса–Уэндроффа и введенной искусственной вязкости рассчитаем давление, плотность и температуру среды непосредственно перед поршнем по значениям скорости ударной волны, полученным из эксперимента. Экспериментальные значения координаты ударной волны L(t) (субстанции) и скорости v(t) приведены на рис. 4, а. На этом же рисунке приведены значения L(t) аппроксимированные сплайнами и сглаженные с помощью фильтра «скользящее среднее» и кривая скорости v(t), полученная взятием производной от сплайновой кривой, рис. 4, б.

Рис. 4. Значения координаты L(t) и скорости ударной волны v(t) соответственно: а) экспериментальные значения; б) сглаженные значения

При расчете термодинамических параметров невозмущенная среда считалась одноатомным газом: постоянная политропы – , давление – p1=105 Па, плотность воздуха – ρ1=1,2 кг/м3, температура воздуха t=15 °C, T1=273,15+15=288,15 K, скорость звука в невозмущенной среде – c=340 м/с. На графиках приведена динамика изменения термодинамических параметров скорости ударной волны, давления, плотности и температуры, непосредственно перед плазменным поршнем, рис. 5.

Рис. 5. Динамика изменения термодинамических параметров: а) скорость ударной волны; б) давление; в) плотность; г) температура

Прокомментируем полученный результат. При вылете плазмы из электрода-ствола со сверхзвуковой скоростью при температуре 1300…3000 К плазма представляет собой газ, распространяющийся в виде струи, которая называется недорасширенной струей. Теория позывает, что при вылете газа из сопла Лаваля [1] со сверхзвуковой скоростью первоначально газ ускоряется, а затем замедляется. Этот процесс периодически повторяется, но уже с меньшей интенсивностью. Что и наблюдается на графике (рис. 4, рис. 5, а).

Результаты

На основе модифицированного алгоритма Лакса–Уэндроффа с введением искусственной вязкости плазмы проведено моделирование параметров плазменного поршня на фронте ударной волны коаксиального магнитоплазменного устройства.

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

Рецензенты:

Канев Ф.Ю., д.ф.-м.н., профессор кафедры ЭСиЭ ЭНИН ФГБОУ ВПО «НИ ТПУ», Национальный исследовательский Томский политехнический университет, г. Томск.

Усов Ю.П., д.т.н., профессор кафедры ЭСиЭ ЭНИН ФГБОУ ВПО «НИ ТПУ», Национальный исследовательский Томский политехнический университет, г. Томск.