CC BY
51
2014 Математика и механика № 6(32)
УДК 621.45, 519.63
Л.Л. Миньков, Э.Р. Шрагер, А.Е. Кирюшкин
О ДВУХ ПОДХОДАХ К МОДЕЛИРОВАНИЮ ГРАНИЦЫ ГАЗОПРИХОДА1
Обсуждаются два подхода к моделированию поверхности газоприхода на основе решения одномерных уравнений газовой динамики: с помощью граничных условий и с помощью правых частей в уравнениях с привлечением дельта-функции Дирака. Проводится сравнение точных и численных решений одномерных стационарных уравнений газовой динамики и численных решений нестационарных уравнений газовой динамики. Рассматриваются случаи, когда газоприход задается постоянным и изменяющимся во времени.
Ключевые слова: граница газоприхода, численное моделирование, газовая динамика, дельта-функция Дирака.
Для моделирования горящей поверхности топлива или границы, с которой осуществляется газоприход, при газодинамическом подходе, как правило, используют граничные условия, задавая расход и полную энтальпию продуктов горения топлива [1, 2]. В работе [3] показано, что задание значений массового расхода и полной энтальпии постоянными на границе области позволяет получить стационарное решение задачи о течении газа в каналах с помощью решения нестационарных уравнений газовой динамики методом установления.
Такой способ моделирования поверхности, с которой осуществляется приход газа, требует, в общем случае, привлечения теории характеристик для определения значений параметров на границе области [3, 4], что связано с определенными трудностями при использовании разностных схем для решения соответствующих задач.
С другой стороны, известно, что при решении задач внутренней баллистики РДТТ в одномерном или квазиодномерном приближениях газо- и энергоприход с боковой поверхности горящего топлива моделируется через правые части в уравнениях сохранения [1, 2, 5].
В работе [6] было показано, что решение однородного дифференциального уравнения в частных производных с неоднородным граничным условием первого рода совпадает с решением неоднородного дифференциального уравнения в частных производных с соответствующим однородным граничным условием. Причем неоднородность в дифференциальном уравнении содержит в качестве множителя дельта-функцию Дирака 5(х).
В данной работе проводится сравнение двух подходов к моделированию границы газоприхода, а именно, через граничные условия и через правые части в уравнениях газовой динамики с привлечением дельта-функции Дирака.
Стационарная задача. Постоянный газоприход
Подход 1. Рассмотрим следующую модельную стационарную задачу. В канал постоянного сечения длиной I через левую границу х = 0 поступает газ с известным массовым расходом О0 и известной полной энтальпией Н0, а через правую
1 Работа выполнена при финансовой поддержке Минобрнауки России в рамках выполнения государственного задания № 10.1329.2014/К.
границу х = I газ вытекает в область постоянного давления рн. Требуется определить параметры газа в случае установившегося течения в трубе. Постановка этой краевой задачи имеет вид
( к р и
0 <*</: ^ = 0; = Ри1 ^Р + ^ = 0;
ёх
Сх
ёх
(1)
к р и
к = 0: ри = О0; + = Н0; (2)
к -1 р 2
х = I: р = Рн . (3)
Разрешая систему дифференциальных уравнений (1) относительно производ-
С р Си Ср
ных —, —, —, получим, что они равны нулю при условии, что скорость газа
Сх Сх Сх
не равна нулю или не совпадает со скоростью звука, и Ф^[кр/р . Отсюда следует, что плотность, скорость и давление газа не будут изменяться вдоль трубы. Тогда давление в канале будет равно давлению рн и краевые условия (2) можно переписать в виде
ри = 00;
к Рн
+—= Н0
к -1 р 2
разрешая которые относительно р и и , получим
к Рн
к -1 О0
1 + 2 Н 0
( к -1 00 ^
-1
(4)
О,
Р = Р1 = -
(5)
Подход 2. Систему однородных уравнений (1) с граничными условиями (2), (3) согласно [6] запишем в виде системы следующих неоднородных уравнений:
0 < х < I:
Сри= О05(х); С(ри + Р) = О0и15(х);
Сх
Сх
, , к р и
а ри I-— +—
\к -1 р 2
Сх
= О0Н 08(х);
(6) (7)
с однородным левым граничным условием
х = 0: и = 0;
и правым граничным условием (3):
х = I: р = Рн . (8)
Физический смысл задачи (6) - (8) следующий. В начале трубы постоянного сечения, левый конец которой закрыт «твердой» поверхностью непроницаемой
и
для газа, а правый открыт, образуется газ с некоторыми известными импульсом и энтальпией.
Интегрируя систему (6) от 0 до I и воспользовавшись условиями (8), получим следующую систему алгебраических уравнений:
^ 2 ^ к Рн и тт
ри = Со; ри + Рн = ОоЩ + Ро ; -;—;-— +— = Н0 .
(9)
к -1 р
Первое и третье уравнения этой системы эквивалентны граничным условиям (2). Тогда, согласно (4), и = и1 и из второго уравнения системы (9) следует, что
р0 = рн. Поэтому решения задач (1)-(3) и (6)-(8) совпадают.
Стационарная задача. Газоприход зависит от давления.
Подход 1. Рассмотрим случай, когда газоприход на левой границе зависит от давления по закону т0ру, где т0 и V - константы, а энтальпия имеет постоянное значение Н0. Задание граничных условий в таком виде соответствует модели поверхности горения, используемой во внутренней баллистике РДТТ [1, 2]. Тогда краевые условия на левой границе (2) принимают следующий вид:
х = 0:
2
V к р и
Ри = т0 Р ; т-т- + -Т = Н0. к -1 р 2
Разрешив (3) и (10) относительно р и и , получим
к р,
¡-V
к -1 т0
1 + 2Н 0
(к-1 л2 к рн
-1
(10)
(11)
р = р2 ^ "
т0 рн
(12)
Подход 2. Вместо системы однородных уравнений (1) с граничными условиями (10) и (3) рассмотрим систему следующих неоднородных уравнений:
0 < х < I:
Л ри ) а (ри 2 + р ) V х( ) ■ = т0р о(х); —--- = т0р и2о(х);
Лх
Лх
, , к р и
а ри I-— +--
\к -1 р 2
Лх
= т0 РV Н 08(х)
(13)
с граничными условиями (7) и (8).
Интегрируя систему (13) от 0 до I и воспользовавшись условиями (7) и (8), получим систему нелинейных алгебраических уравнений
v 2 v к рн и2 ТТ
ри = т0Р(>; ри + р„ = т0Р0и2 + Р(>; -¡—^ +^т = Н0
к -1 р 2
(14)
которую разрешим относительно трех неизвестных величин: р, и , р0. Введем безразмерную переменную х = р0/рн и два безразмерных параметра:
и
2
и = и2т0!у , к = Но (т0!р1и
Тогда система уравнений (14) сводится к одному нелинейному уравнению относительно х:
к 1 (_ х -1А 1 (_ х -1" ,
--1 и +-|+ —I и +-| = к .
к -1 ху ) 21 х"
(15)
Решение этого уравнения для 0 <у< 1, 1 < к < 2 и к < 1 имеет вид х = 1. Поскольку при х = 1 уравнение (15) сводится к уравнению для нахождения и2 (11), то р0 = рн и и = и2. Следовательно, в стационарной постановке моделирование течения газа в канале с горящим левым торцом можно проводить как с помощью краевой задачи (1), (10), так и с помощью краевой задачи (13), (7), (8). Результат решения в обоих случаях будет один и тот же.
Численное решение стационарной задачи методом установления
Рассмотрим задачу о втекании газа в канал постоянного сечения через левую границу с известным газоприходом т0 ру и полной энтальпией Н0. Через правую границу газ вытекает в область постоянного давления рн. В начальный момент времени газ покоится.
Данная задача описывается с помощью следующей системы нестационарных уравнений газовой динамики, граничных и начальных условий:
Подход 1.
ф + фи = 0_ д г д х
дри д(ри2 + р2
д г
д х
= 0:
ф
1 р и | - ( к р и2
-— + —| дри I-— + —
к -1 Р 2 ) + у к -1 р 2 у = 0
д г
д х
Начальные условия:
Граничные условия:
г = 0: р = Рн: и = 0; р = рн .
х = 0: ри = т0 р
к р-+- = Н0:
Подход 2.
к -1 р
Ф , дри V*/ ч
— + = т0р б(х): 5 г д х 0 v 7
х = I:
р = рн
дри д(ри 2 + р 2
д г
д х
= т0руиЪ(х2:
(16)
(17)
(18)
(19)
(20)
(21) (22)
p u2 | _ | k p u2
5Р1т~7 + V I дРи\^~л i
k -1 p 2 I I k -1 p 2 , v
y+-^---^ = m0 pv H0S(x). (23)
dt dx Начальные условия:
t = 0: p = pH; u = 0; p = pH. (24)
Граничные условия:
x = 0: u = 0; x = l: p = pH. (25)
Задачи (16) - (20) и (21) - (25) решаются численно по разностной схеме первого порядка точности на основе метода конечных объемов [4], причем потоки на границах ячеек разностной сетки определяются по методу Ван-Лира [7]. Шаг по времени определяется из условия: т = Khxjmax (|ui -ci |, |uz- +ci |), число Куранта
K задается равным 0,5. Дельта-функция Дирака аппроксимируется обратным размером приграничной ячейки 1/hx . Стационарное решение этих задач находится
методом установления при значениях параметров равными: m0 = 10-2; H0 = 102; v = 0,5; pIJ = 1; рн = 1; k = 1,4; l = 1.
Точное стационарное решение этой задачи дается формулами (11) и (12), по которым получаются следующие значения плотности, давления и скорости: p = 0,03501428; p = 1,0000000; u =0,2855978. Такой же результат получается при решении задачи с использованием подхода № 1, (16) - (20).
Расчеты, выполненные на разностной сетке, состоящей из 20 ячеек, показали, что традиционная реализация левого граничного условия (твердая стенка) с помощью решения задачи Римана в подходе №2, (21) - (25), приводит к тому, что значения плотности, давления и скорости в приграничной ячейке могут значительно отличаться от точных: p = 0,0365043; p = 1,0525237; u = 0,1418703, в то время как во всех остальных ячейках решение близко к точному: p = 0,03501503; p = 1,0000000; u = 0,2929958. При этом измельчение разностной сетки не влияет на численное решение. Точного решения удается добиться, если давление на левой границе полагать равным давлению в центре приграничной ячейки на нижнем временном слое.
Численное решение нестационарной задачи
Для сравнения двух подходов в моделировании границы газоприхода в случае нестационарного процесса рассмотрим задачу о втекании газа в канал постоянного сечения через левую границу с постоянной полной энтальпией H0 и газоприходом, изменяющимся во времени по закону G = G0 (1 - exp(-at)), что соответствует возрастанию расхода от нулевого значения до значения G0 при t ^го. Через правую границу газ вытекает в область постоянного давления рн. В начальный момент времени газ покоится.
Данная задача описывается системой уравнений, граничными и начальными условиями (16) - (20) в случае подхода №1 и системой уравнений, граничными и начальными условиями (21) - (25) в случае подхода № 2, где вместо m0 pv следует записывать G0 (1 - exp(-at)). Параметр a отвечает за скорость нарастания газоприхода на левой границе области.
Для получения численного решения задачи будем использовать метод, изложенный в предыдущем разделе.
Расчеты проводились при следующих значениях параметров к = 1.4, 00 = 10-2, Н0 = 102, рн = 1, рн = 1 на трех разностных сетках: грубой (число ячеек N = 20), средней (Ы = 200) и на мелкой (Ы = 2000).
На рис. 1 и 2 показаны изменения соответственно давления и скорости в первой ячейке разностной сетки в зависимости от времени для значения а = 1. Видно заметное различие в результатах расчетов, полученных по двум способам моделирования границы газоприхода на грубой разностной сетке (рис. 1, а и 2, а), которое исчезает по мере измельчения сетки (рис. 1, б, в и рис. 2, б, в). Моделирование границы газоприхода с помощью правых частей системы уравнений (21) -(23) дает завышенные значения давления и заниженные значения скорости в первой ячейке по сравнению с моделированием границы газоприхода через граничные условия, хотя качественное поведение параметров от времени одинаковое. Различие в скоростях, вычисленных по двум подходам на грубой разностной сетке, не превышает 15 %, а различие в давлениях не превышает 3 %. Эти различия уменьшаются при измельчении разностной сетки и при увеличении времени физического процесса, что говорит о том, что моделирование границы газоприхода с помощью системы уравнений (21) - (23) и граничных условий (25) выполнено корректно.
Рис. 1. Изменение давления от времени в первой ячейке разностной сетки, а = 1; а - N = 20, б - N = 200, в - N = 2000; кр. 1 - подход № 1, кр. 2 - подход № 2
Рис. 2. Изменение скорости от времени в первой ячейке разностной сетки, а = 1: а - N = 20, б - N = 200, в - N = 2000; кр. 1 - подход № 1, кр. 2 - подход № 2
На рис. 3 и рис. 4 показаны изменения давления и скорости во времени в первой ячейке для двух различных подходов, когда скорость нарастания газоприхода более интенсивная, а = 10. Увеличение скорости нарастания газоприхода несколько увеличивает различие в поведении параметров от времени в первой ячейке, даже на качественном уровне (рис. 3, а и рис. 4, а), которое также исчезает при измельчении разностной сетки, (рис. 3, б, в и рис. 4, б, в).
Р-1,21,110,90,8-
г , _1
у" -- 2
ч /
а
Р
1,21,1 1
0,9 0,8
Г -- 2
/ \
\
\
б V
Р
1,2' 1,1 1
0,9 0,8
г \ - - 2
\
\
\ и
в 1 ■ ■
0
1
0
1
0
1
2
Рис. 3. Изменение давления от времени в первой ячейке разностной сетки, а = 10; а - N = 20, б - N = 200, в - N = 2000; кр. 1 - подход №1, кр. 2 - подход № 2
0,30,20,1 0-
. — 1
/V //
Г— —
1 а
0,30,20,10
--2
Г
б
0,3 0,20,1 0
0
1
0
1
0
1
2
Рис. 4. Изменение скорости от времени в первой ячейке разностной сетки, а = 10; а - N = 20, б - N = 200, в - N = 2000; кр. 1 - подход № 1, кр. 2 - подход № 2
■ - - 2
/
1 в
Р-1,31,21,1 1
1 о-о 2
\
1 I I I а ^ 1
0,350,30,250,20,15-
е-о-о 1 2 /
/
У
V у
б
Р-1
0,80,60,40,2
»11 1 2 /
{
/
/ в
0 0,2 0,4 0,6 0,8 х
0 0,2 0,4 0,6 0,8 х
0 0,2 0,4 0,6 0,8 х
Рис. 5. Зависимость параметров от пространственной координаты на момент времени Г = 1: а - давление; б - скорость; в - плотность; кр. 1 - подход № 1, кр. 2 - подход № 2; N = 20; а = 10
Г
Г
г
и
и
и
Г
г
г
и
Для сравнения результатов расчета в других точках области на рис. 5 представлены распределения давления, скорости и плотности вдоль координаты x на момент времени t = 1, полученные на грубой разностной сетке. Как видно, два различных подхода к моделированию границы газоприхода дают результаты, хорошо согласующиеся между собой во всей области решения, за исключением первой расчетной ячейки, где имеют место расхождения, как уже отмечалось, для давления и скорости (рис. 5, а и б).
Интересно отметить, что использование в правой части в уравнении (22) переменного значения скорости u, а не постоянного значения u7, как это было сделано в (6), не повлияло существенным образом на результаты вычислений.
Заключение
Рассмотрены два подхода к моделированию границы газоприхода на основе решения одномерных уравнений газовой динамики, в одном из которых газоприход реализуется с помощью граничных условий, а в другом - с помощью правых частей уравнений газовой динамики с привлечением дельта-функции Дирака.
Показано, что точные решения стационарных задач, полученные в двух подходах, совпадают друг с другом для случая постоянного газоприхода и газоприхода, задаваемого по закону m0 pv. Показано, что численное решение стационарной задачи, полученное методом установления, в подходе, реализующем границу газоприхода через правые части уравнений, совпадает с аналитическим решением во всей области, если давление на левой грани приграничной ячейки задавать равным давлению в центре этой ячейки.
На основе численного моделирования границы газоприхода с изменяющимся во времени расходом газа, проведенного по разностной схеме первого порядка точности, показано, что численные решения, полученные по двум рассмотренным подходам, сходятся друг к другу при измельчении разностной сетки.
ЛИТЕРАТУРА
1. Райзберг Б.А., Ерохин Б.Г., Самсонов К.П. Основы теории рабочих процессов в ракетных системах на твердом топливе. М.: Машиностроение, 1972. 383 с.
2. Ерохин Б.Т., Липанов А.М. Нестационарные и квазистационарные режимы работы РДТТ. М.: Машиностроение, 1977. 200 с.
3. Рынков А.Д. Математическое моделирование газодинамических процессов в каналах и соплах. Новосибирск: Наука, 1988. 222 с.
4. Годунов С.К. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с.
5. Ерохин Б.Т. Теория внутрикамерных процессов и проектирование РДТТ: учебник для высших учебных заведений. М.: Машиностроение, 1991. 560 с.
6. Бутковский А.Г. Характеристики систем с распределенными параметрами. М.: Наука, 1979. 224 с.
7. Van Leer B. Flux-vector splitting for the euler equations // Lecture Notes in Physics. 1982. Vol. 170. P. 507-512.
Статья поступила 03.11.2014 г.
Minkov L.L., Shrager E.R., Kiryushkin A.E. ON TWO APPROACHES TO MODELLING THE INFLOW BOUNDARY
Two approaches to modeling inflow boundary through solving one-dimensional equations of gas dynamics by means of boundary conditions and right-hand sides in the equations involving the Dirac delta function are discussed. Comparisons of numerical solutions of one-dimensional
102
n.n. MuHbHOB, Э.Р. Wparep, A.E. KupmwHUH
stationary equations of gas dynamics with the exact ones, as well as between numerical solutions of unsteady gas dynamics equations are performed for two approaches. Three kinds of boundary conditions are considered, namely, a constant inflow, pressure-dependent inflow, and inflow varying in time. It is shown that the numerical solutions obtained based on the finite-difference scheme of first order accuracy by the two considered approaches converge to each other with the mesh refinement. The numerical solution for the steady state problem coincides with the analytical one if the pressure at the boundary cell face is set equal to the pressure in the center of the boundary cell.
Keywords: inflow boundary, numerical modeling, gas-dynamics, Dirac delta-function.
MINKOVLeonidLeonidovich (Doctor of Physics and Mathematics, Assoc. Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
SHRAGER ErnstRafailovich (Doctor of Physics and Mathematics, Assoc. Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
KIRYUSHKIN Aleksandr Evgenievich (student, Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
REFERENCES
1. Rayzberg B.A., Erokhin B.G., Samsonov K.P. Osnovy teorii rabochikh protsessov v raketnykh sistemakh na tverdom toplive. Moskow, Mashinostroenie Publ., 1972. 383 p. (in Russian)
2. Erokhin B.T., Lipanov A.M. Nestatsionarnye i kvazistatsionarnye rezhimy raboty RDTT. Moskow, Mashinostroenie Publ., 1977. 200 p. (in Russian)
3. Rychkov A.D. Matematicheskoe modelirovanie gazodinamicheskikh protsessov v kanalakh i soplakh. Novosibirsk, Nauka Publ., 1988. 222 p. (in Russian)
4. Godunov S.K. Chislennoe reshenie mnogomernykh zadach gazovoy dinamiki. Moskow, Nauka Publ., 1976. 400 p. (in Russian)
5. Erokhin B.T. Teoriya vnutrikamernykh protsessov i proektirovanie RDTT: uchebnik dlya vysshikh uchebnykh zavedeniy. Moskow, Mashinostroenie Publ., 1991. 560 p. (in Russian)
6. Butkovskiy A.G. Kharakteristiki sistem s raspredelennymi parametrami. Moskow, Nauka Publ., 1979. 224 p. (in Russian)
7. Van Leer B. Flux-vector splitting for the euler equations. Lecture Notes in Physics, 1982, vol. 170, pp. 507-512.
