CC BY
56
УДК 621.384.647.001.5
МОДЕЛИРОВАНИЕ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ И ОЦЕНКА ТЕРМОДИНАМИЧЕСКИХ ПАРАМЕТРОВ УДАРНОЙ ВОЛНЫ ПЛАЗМЕННОГО ГАЗА КОАКСИАЛЬНОГО МАГНИТОПЛАЗМЕННОГО УСКОРИТЕЛЯ
А.А. Сивков, Ю.Н. Исаев, О.В. Васильева, А.М. Купцов
Томский политехнический университет E-mail: vasileva.o.v@mail.ru
Работа посвящена решению одномерного уравнения гидрогазодинамики для коаксиального магнитоплазменного ускорителя с помощью модифицированного алгоритма Лакса~Уэндроффа с оптимальным выбором параметра регуляризации - искусственной вязкости. На основе предложенного алгоритма в среде MathCAD рассчитаны термодинамические параметры ударной волны перед плазменным поршнем при ее вылете из коаксиального магнитоплазменного ускорителя.
Ключевые слова:
Магнитоплазменный ускоритель, математическое моделирование, ударная волна, искусственная вязкость, гидрогазодинамика. Key words:
Magneto plasma accelerator, mathematical modeling, shock wave, artificial viscosity, fluid dynamics.
При вылете плазменной субстанции из ствола коаксиального магнитоплазменного ускорителя перед ее головной частью образуется отсоединенная ударная волна [1]. Произведем оценку термодинамических параметров за ударной волной. Для этого примем некоторые упрощения - субстанцию условно будем называть затупленным телом или поршнем, расчет будем производить в одномерном случае. Если перейти к системе координат связанной с поршнем, то невозмущенный газ-воздух будет двигаться на поршень со скоростью поршня. При моделировании движения газовой волны на твердую преграду данную модель можно представить как движение двух одинаковых волн на встречу друг другу. Приведем схему этой задачи называемой нестационарным газодинамическим разрывом. Схематически картину течения можно изобразить в виде изображений, представленных в табл. 1 [2-4]. Конфигурация содержат контактный разрыв, на котором имеет место разрыв плотности р, а давление р и скорость V остаются непрерывными, табл. 1. Для решения одномерной нестационарной газодинамической задачи будем использовать уравнения сохранения массы, импульса и энергии, записанные в дивергентной форме:
Таблица 1. Конфигурации нестационарного газодинамического разрыва
Схематическое изображение Комментарии
Случай наблюдается при столкновении двух масс газа, движущихся навстречу друг другу.
График характеристики волн: УВ - ударная волна; КР - контактный разрыв (пунктирная линия - скачек плотности р). На границе раздела - контактном разрыве не может быть скачка давления р.
На границе раздела - контактном разрыве не может быть скачка скорости V
Ф і д(up) = 0
дї дx
д(up) і дjpu21 p) = о
дї дx
д(p(е + pu2/2))
(1)
д(ш(pe і pu2 /21 p))
д? дх
Здесь р - плотность газа; р - давление газа; и -скорость распространения газа; є - внутренняя энергия газа; /, х - время и координата.
Те же уравнения (1), записанные в векторной форме, удобной для численной реализации:
( п \
s =
p up
p(e і pu2 / 2)
—
f (s) =
дs
——------------1 “
дї дx
up s2
pu21 p = s2 / s і p
Ku(pe і pu2 /21 p) j кs2(s31p)/ s j
f (s)
= 0,
(2)
где s - вектор консервативных переменных; f(s) -вектор потока.
Описание алгоритма для расчета термодинамических параметров ударной волны коаксиального магнитоплазменного ускорителя в программно-интегрированной среде MathCAD:
1. Замена дифференциальных уравнений в частных производных конечными разностями (2).
2. Добавление оптимального параметра регуляризации - искусственной вязкости в среде MathCAD.
3. Выбор оптимальной искусственной вязкости, используя точное известное решение (задача Сода).
4. Апробация разработанного алгоритма расчета термодинамических параметров в точке торможения (аналитически), табл. 2.
5. Расчет динамики изменения термодинамических параметров перед плазменным поршнем при вылете из ускорителя.
Для численного решения системы уравнений использовался модифицированный алгоритм Лак-са-Уэндроффа [5], который заключается в том, что уравнения в частных производных заменяются конечными разностями. В конечных разностях появляется неустойчивость в виде высокочастотных шумов из-за наличия сильных ударных волн. К данному алгоритму нами была добавлена искусственная вязкость. В режиме online был подобран оптимальный параметр регуляризации (искусственной вязкости) для обеспечения регуляризации решения и подавления его шумовой составляющей.
Формируется массив значений для каждого слоя, используя среду MathCAD:
LW(р,p, u, N, h,T, M, у, ц):= si ^ s(p, p, u, N, у) for i e 0..M
f i ^ f (si, Y)
S0,5 ^ F0,5(S1, f 1 h,T) f0,5 ^ f (S0,5 , Y)
for к e 0..N -1
(s2 > s2 к ,0 ,0
s2 k ,i sik ,i T — h
Vs2к,2 J Vsik,2 J
f>A-f0,5 k-
i f0,5k
- f0,5 k
Vf0,5k+l,2 f0,5k,2 J
TЦ ~2h
( ' 0 ^ (
if к — i < 0, 0
0
V V J V
sL
— 2 sik,0 + sik-i,
Vi
si — 2 si + si
k+i,i 2 sik,i + k—i,i
sL
— 2 sik ,2 + sik —i,
( vi ^ \
Sin s2k ,0
S 2k S 3,
s2k ,i s2
s2 'I s2 N ,0 (si ^ siN ,0
s2 s2 N ,i siN ,i
s2N,2 J VsiN,2 J
SIn,i л (s2 ^ s 2 N ,0
s 2n ,i s 2 s 2 N ,i
S 3n ,i j s2 V s2 N,2 J
( f lN,0 — f lN—i,0 ^ f lN ,i — f lN—i,i f lN,2 — f lN—i,2 J
51 ^ 52
( 51 ^
5 2 5 3
где 1Ж - функция Лакса-Уэндроффа; N=200; М=200 - число точек разбиения пространственного и временного интервалов соответственно; к, т -шаг по пространству и времени соответственно; 7=5/3 - показатель политропы; /л - искусственная вязкость; 5, / - вспомогательные функции для формирования массива значений.
Величина искусственной вязкости определялась из невязки-рассогласования (рис. 1). Точное известное решение задачи Сода сравнивалось с нашим алгоритмом, если рассогласования составляли не более 10 %, то коэффициент регуляризации (искусственная вязкость) считался оптимальным, рис. 2.
Результат расчета динамики изменения термодинамических параметров в относительных единицах приведен на рис. 3. До столкновения волн величины давления, плотности и температуры в средах были одинаковы, а скорости волн одинаковые по величине, но разные по направлению (знаку). На рисунке приведен момент времени после столкновения двух ударных волн при числе Маха, равном 1,5. Осью симметрии рисунков является фронт плазменного поршня.
Проверка работы алгоритма была проведена на расчете критических параметров - давления, плотности и температуры торможения рТ, рТ, ТТ при заданных начальных данных невозмущенного газа р0=1/у, р0=1, Т0=1 (в относительных единицах) и заданном числе Маха М0, рис. 4, а. Параметры торможения вычислялись по известным соотношениям [6, 7]. Давление, плотность и температура из невозмущенной среды через скачок уплотнения пересчитываются через ударную адиабату Гюго-нио-Ренкина (ударная волна):
2y , ,2 Y — i
pi = p01 7 M0 —
Y + i
Y + i
T = T Ti T0
2 J
pi p0
Y — i
Y +i
(Y — i) + 2/M0
Y + i
2Y
i
Y+i M0 J{y—i
Рис. 1. Восстановленные термодинамические параметры ударной волны и волны разряжения с неоптимальной вязкостью Л=10-
Акг/д:
Р, МПа
О 0,5 1 1,5 2 х,м
Рис. 2. Восстановленные термодинамические параметры ударной волны и волны разряжения с оптимальной вязкостью }л=Ю~6
а б в
Рис. 3. Случай двух ударных волн: а) плотность р; б) давление р; в) скорость у; г) энергия Е
Ро
Vo
ро
х,м
а
^Мп=10
Рис. 4. Параметры ударной волны: а) критические параметры; б) распределение плотности на фронте волны при различных числах Маха М0
В нашем случае эти формулы вырождаются в простые и удобные для инженерных расчетов соотношения:
1 2 4
рі - —1) ро, рт - "
1 + -
М2
р0,
тт - — т 16
1
1 -М2 I(3 + М0)то-
Для расчета параметров торможения нужно использовать адиабату Пуассона (волна разрежения):
Рт - I 1 + М1
2 7-1 >-1
7-1
рт- I 1 + М1
2 7-1
т„ — 1+м2 у-1
Р1>
?1-
здесь Д - числа Маха после ударной волны, определяется через число Маха в невозмущенной среде выражением:
Мх -,
2 + (у- 1)Мо2 >Мо2 - (у-1)'
Для инженерных расчетов получаем:
М1 - I3 + М
5М02 -1
(
Рт =
1+МТ
3
(
1 -2
1+-3-1 р.
т =
±т
1+-г. і т.
(3)
Начальные данные приведем к относительным единицам p0=1/7=0,6, р0=1, Г0=1.
Результаты сравнения методов сведем в табл. 2.
Таблица 2. Результаты сравнения инженерного (И) и программного (П) расчетов
Число Маха М0 Р Рт Р1 Рт т Тт Расчет
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):
— -1/К -
Р (
_4__________ К -1
л ’ КР л
1 +
-о2 у
1+-
М
о о
5 -^Кр (1 + о,6Кр).
Таблица 3. Значения расстояния между ударной волной и торцом поршня Б и коэффициентом плотности торможения Кр в зависимости от числа Маха М0
М0 1,5 2 3 5 10
Кр 0,583 0,437 0,333 0,280 0,258
5 1,031 0,835 0,693 0,618 0,586
Как видно из табличных данных (табл. 3) программа дает хорошие результаты. Графическое изображение скачков уплотнения приведено на рис. 4, б. Торец поршня находится в точке х=0,5 [8].
С помощью численной схемы Лакса-Уэндроф-фа и введенной искусственной вязкости рассчитаем давление, плотность и температуру среды непосредственно перед поршнем по значениям скорости ударной волны, полученным из эксперимента.
Экспериментальные значения координаты ударной волны Щ) (субстанции) приведены на рис. 5, а. На этом же рисунке приведены значения ДО, аппроксимированные сплайнами и сглаженные с помощью фильтра «скользящее среднее», рис. 5, б.
1.(1) ,СМ'
100-
У,
*
6 1*^
20
40
60
80
І, МКС
Рис. 5. Значения координаты ударной волны 1(1): а) экспериментальные значения для субстанции от времени; б) сглаженные
Рис. 6. Значения скорости ударной волны у^): а) полученные разностной аппроксимацией; б) сглаженные
На следующем рис. 6, а, приведены значения скорости v(t), полученные по экспериментальным данным (ломаная кривая) с помощью аппроксимации производной центральными разностями. На этом же рисунке приведена кривая скорости v(t), полученная взятием производной от сплайно-вой кривой, рис. 6, б.
При расчете термодинамических параметров невозмущенная среда считалась одноатомным газом: постоянная политропы - 7=5/3; давление -p1=105 Па; плотность воздуха - р;=1,2 кг/м3; температура воздуха г?=15 °С; ^=273,15+15=288,15 К; скорость звука в невозмущенной среде -с=340 м/с. Результаты расчета сведем в табл. 4.
На графиках приведена динамика изменения термодинамических параметров скорости ударной волны, давления, плотности и температуры непосредственно перед плазменным поршнем, рис. 7.
Прокомментируем полученный результат. При вылете плазмы из электрода-ствола со сверхзвуковой скоростью при температуре 1300...3000 К плазма представляет собой газ, распространяющийся в виде струи, которая называется недорасширен-ной струей. Теория позывает, что при вылете газа из сопла Лаваля [5] со сверхзвуковой скоростью первоначально газ ускоряется, а затем замедляется. Это процесс периодически повторяется, но уже с меньшей интенсивностью. Что и наблюдается на графике (рис. 6, 7, а).
Выводы
На основе модифицированного алгоритма Лак-са-Уэндроффа с введением искусственной вязкости плазмы проведено моделирование газодинамических процессов и оценка термодинамических параметров ударной волны плазменного газа.
Таблица 4. Значения термодинамических параметров невозмущенной газовой среды
Время, мкс Скорость УВ, м/с Число Маха, М0 Плотность, кг/м3 Давление, МПа Температура, К
1,35 1251,017 3,371 4,334 1,715 1379,953
4,05 1281,213 3,453 4,371 1,796 1433,296
6,75 1377,217 3,712 4,477 2,068 1611,342
9,45 1514,445 4,081 4,603 2,491 1888,169
12,15 1688,459 4,550 4,727 3,086 2276,987
14,85 1863,405 5,022 4,824 3,749 2710,476
17,55 1981,377 5,339 4,878 4,233 3026,899
20,25 1972,403 5,316 4,874 4,195 3002,147
22,95 1820,368 4,906 4,803 3,579 2599,877
25,65 1589,745 4,284 4,661 2,740 2051,233
28,35 1372,093 3,698 4,472 2,053 1601,514
31,05 1230,526 3,316 4,308 1,660 1344,479
33,75 1150,249 3,099 4,195 1,457 1211,150
36,45 1090,288 2,938 4,100 1,313 1117,429
39,15 1015,346 2,736 3,966 1,145 1007,344
41,85 934,749 2,519 3,800 0,978 897,697
44,55 875,289 2,359 3,661 0,863 822,618
47,25 856,273 2,308 3,614 0,828 799,646
49,95 868,344 2,340 3,644 0,850 814,169
52,65 880,228 2,372 3,673 0,873 828,667
55,35 865,884 2,333 3,638 0,846 811,193
58,05 820,842 2,212 3,520 0,765 758,193
60,75 764,325 2,059 3,359 0,670 695,694
63,45 718,564 1,936 3,217 0,598 648,354
66,15 690,637 1,861 3,125 0,556 620,900
68,85 675,264 1,819 3,073 0,534 606,252
71,55 665,917 1,794 3,040 0,521 597,506
74,25 669,298 1,804 3,052 0,526 600,655
76,95 685,696 1,848 3,109 0,549 616,156
79,65 699,027 1,884 3,153 0,569 629,034
82,35 678,689 1,829 3,085 0,539 609,486
85,05 615,025 1,657 2,856 0,452 552,028
87,75 544,339 1,467 2,579 0,366 494,858
90,45 499,798 1,347 2,395 0,317 462,414
I, МКС
/, МКС
р, кг/м‘
Г, К
5- •
4 ‘
З1 ■ \
2"
1"
3200- ’
240(Г
1600- -
800
О 10 20 30 40 50 60 70 80 90
в К мкс
О 10 20 30 40 50 60 70 80 90
г 1, мкс
Рис. 7. Динамика изменения термодинамических параметров: а) скорость ударной волны; б) давление; в) плотность; г) температура
При моделировании учтены подавляющие неустойчивые высокочастотные колебания, что позволяет сузить область неоднородности и выделить только гладкие решения. Результаты расчета газо-
динамических параметров в точке торможения совпадают с литературными данными.
Работа выполнена при финансовой поддержке гранта РФФИ проект № 09-08-01110.
СПИСОК ЛИТЕРАТУРЫ
1. Сивков А.А., Сайгаш А.Я., Пак А.А., Евдокимов А.А. Прямое получение нанодисперсных порошков и композиций в гипер-скоростной струе электроразрядной плазмы // Нанотехника. -2009. - № 2 (18). - С. 38-44.
2. Морозов А.И. Введение в плазмодинамику. - М.: Физматлит, 2008. - 613 с.
3. Колесников П.М. Электродинамическое ускорение плазмы. -М.: Атомиздат, 1971. - 388 с.
4. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика Т. 8: Электродинамика сплошных сред. - М.: Наука, 1992. - 664 с.
5. Зализняк В.Е. Основы вычислительной физики. Ч. 1. Введение в конечно-разностные методы. - М.: Техносфера, 2008. - 224 с.
6. Черняк В.Г, Суетин П.Е. Механика сплошных сред. - М.: Физматлит, 2006. - 352 с.
7. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. - М.: Физмат-лит, 2008. - 656 с.
8. Трубников Б.А. Теория плазмы. - М.: Энергоатомиздат, 1996. -464 с.
Поступила 18.05.2011 г.
УДК 621.762.4.04.016.2
ОПРЕДЕЛЕНИЕ ОПТИМАЛЬНЫХ РЕЖИМОВ ИЗГОТОВЛЕНИЯ ВЫСОКОПЛОТНОЙ КЕРАМИКИ ИЗ ПОРОШКА КАРБИДА БОРА МЕТОДОМ СПЕКАНИЯ В ПЛАЗМЕ ИСКРОВОГО РАЗРЯДА
О.Л. Хасанов, Э.С. Двилис, А.О. Хасанов, З.Г. Бикбаева,
В.В. Полисадова, В.М. Соколов, А.А. Качаев, Я.В. Валова
Томский политехнический университет E-mail: aokhasanov@tpu.ru
Приведены экспериментальные результаты по спеканию в плазме искрового разряда керамики из промышленного порошка карбида бора на установке Spark Plasma Sintering System SPS-515S. Проведён подбор режимов спекания: температуры, времени идавления прессования. Показано, что при оптимальном технологическом режиме (температуре спекания 1900...1950 °С, давлении прессования 45...90 МПа, времени спекания 10 мин.) микротвёрдость керамики достигает значения Нм=35,45...36,50ГПа при трещиностойкости Kt=4,22...5,62 МПа-м’/2иотносительной плотности ротн=9в,4...98,8 %. SPS-спекание способствует понижению температуры и времени спекания керамики из порошка карбида бора по сравнению с методом горячего прессования и формирует изотропную зеренную структуру с хорошо сформированными межзеренными границами.
Ключевые слова:
Керамика, карбид бора, искровое плазменное спекание, микротвердость, трещиностойкость.
Key words:
Ceramic, boron carbide, spark plasma sintering, microhardness, crack growth resistance.
В последние годы в России и за рубежом проводятся интенсивные работы в области создания противопульной бронезащиты на основе керамических материалов. Благодаря меньшему весу керамическая бронезащита может применяться в средствах индивидуальной защиты и защиты транс -портных средств, где есть ограничения по весу и одновременно требуется высокая пулестойкость. Наиболее предпочтительным материалом является карбид бора, сочетающий высокую твердость (третий по твердости материал после алмаза и кубического нитрида бора) с низким удельным весом (2,52 г/см3). Температура плавления карбида бора 2350 °С; карбид бора обладает высокой износостойкостью и химической инертностью. Благодаря этим свойствам карбид бора находит применение в технике, хотя относительно высокая стоимость ограничивает его применение [1, 2].
Высокие физико-технические характеристики и эксплуатационные качества керамики могут быть реализованы в материале с тонкой однородной структурой, хорошо сформированными меж-зёренными границами и плотностью, близкой к теоретической. В ряде работ [2, 3] отмечалось, что формирование плотной изотропной мелкокристаллической структуры в процессе спекания приводит к созданию керамики с высокими значениями твердости и прочности, что определяет ресурс службы керамических изделий.
Одним из предпочтительных способов получения высокоплотных керамик из порошка карбида бора является способ горячего прессования [3, 4]. Однако этот способ обладает рядом существенных недостатков: высокие температуры спекания (2250. 2280 °С), длительность спекания больше 20 мин [1, 3]. Достоинствами метода являются от-
