CC BY
95
Численное моделирование течений детонации газовых смесей методом конечного объема
С. Н. Мдртюшов Государственная дума РФ, Москва, Россия e-mail: martyush@mail.ru
Я. Г. МАРТЮШОВА Московский авиационный институт, Россия e-mail: martyush@rambler.ru
System of equations for an ideal non viscous gas and the kinetic equation in an integral form were used for investigation of an axially symmetric detonation in gas mixes. For numerical simulation of this flows the TVD Harten's scheme of second order of accuracy was used. Calculation grids where constructed on the basis of the Thompson algorithm. Two problems of engineering importance were calculated, namely, flows arising during an ignition of the pulsing detonation engine and the detonation process arising under the containment of a nuclear reactor. Chapman-Juguet wave propagation was calculated as a test problem.
Введение
Система уравнений идеального газа и кинетических уравнений в интегральной форме для осесимметрических течений может быть представлена в следующем виде [1, 2]:
4 / QdV + / nFdS + Ф = 0; (1)
dt J J
V S
Q = (p, m, pe, рв, pa); F = (m, (mm)/p + PI, m(e + P)/p, (mB)), B = (в,а); Ф = (0, 0, 0, 0, pw/з,pwa); P = pRT, e = RT/(y — 1) + V2/2 + pq; da — 1
Wa = =-= —kiP exp(—Ei/RT);
dt Tind
в 1 Г>2п2 f—E2n m 2 / E2 + q
We = ~T = —k2P в exp —— — (1 — в) exp(—
dt 1\ЯТУ,/1УЯТ
где р,Р,е, т,Т, Я — плотность, давление, удельная энергия, вектор количества движения, температура и газовая постоянная соответственно. Параметры а, в характеризуют продвижение реакции: в периоде индукции в = 1, 0 < а < 1, на стадии экзотермической реакции а = 0, 0 < в < 1; 7 — отношение удельных теплоемкостей — постоянная, равная 1.4, для разбавленных смесей кислорода с водородом; д — удельный тепловой эффект химической реакции; к\,к2, Е\, Е2 — константы реакции.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
1. Разностный алгоритм
Поверхностный интеграл (1) для контрольного объема в виде пространственного шестигранника в форме клина (в осесимметрических координатах — четырехугольник) может быть представлен следующим образом:
4
nFdS T-1HSj + SkE; j=i
где H = (pu,pu2 + P, рvu,рu + Pu, рви, pau)T, E = P(0, 0,1, 0,0, 0)T, Sk — площадь боковой грани шестигранника, T — матрица перехода от цилиндрических координат к криволинейным, n = (nx, ny), т = (тх, ту) — векторы нормали и касательной к грани Sj, u = nV, v = тV. Вектор потока на грани H вычислялся по формуле:
H = /нг + Hr + ^ Sj |Aj |(
V j=i
где Aj, eAj — собственные значения и правые собственные векторы матрицы dH/dQ, вычисленные на грани ячейки с использованием газодинамических величин, полученных осреднением Ройе [3]:
Pj+1/2 = y/PjPj+1, K = v/pj/(^JJj + VЛ, (2)
A, A,A, /3,a = Ku, v, h, в, aj + (1 - K)u,v, h, A, aj+1, h = yP/(y - 1)р + V2/2, V = (u, v).
Вектор a = LAQ — произведение вектора неизвестных системы (1) в дельта-форме на матрицу левых собственных векторов, компоненты 0j вектора a имеют вид
ai = -Ар + AApvA2, § = (y - 1)/2А2, (3)
=1 - + a2Ар + (y - 1)(AApu + AApv - Аре + Ар/)/A2, O = §(A2 + a2 - А/ААр/2 + A/2 - A$Aрu - §(AAрv - Аре + Ар/), 04 = 03 + (ААр - Ар^/А, 05 = - А(0з + 04) + Ар/, 06 = Ара - 0(03 + 04),
где, например, Aрu = р^^^-^ - рjuj.
Правые собственные векторы матрицы дH/dQ Aj вычисляются по формулам
А = (0, 0, А, А2, 0,0)T, А = (1,а,а,а2, 0, 0)T, (4)
ез = (1, А + A, A, H + АА, /3,0)T, е4 = (1, А - A, A, H - АА, A, 0)T,
е5 = (0,0,0,1,1,0)T, А = (0,0,0,0,0,1)T, а1)2,5,б = А, а3 = А + А, а4 = А - А.
Численный алгоритм, разработанный на основе метода конечного объема с использованием TVD-схемы [4] второго порядка точности для системы уравнений идеального невязкого газа, использован для численной аппроксимации системы уравнений (1). Оператор шага по времени разностной аппроксимации системы (1) расщепляется на симметричную последовательность операторов шага в направлении [5]:
дп+2 = L(2At)Qn, L(2At) = L¿(Aí)Lj(At)Lj(Aí)L¿(Aí),
Li(—t) = I — — [Fi+l/2Si+l/2 — ^-1/23^1/2 + йкЕ],
А/2 = 1/2[_Ро + А + Я?1/2Ф а/2 ]. (5)
Обозначим компоненты ФО/2 как (ФО/2)к, тогда: ( ФО/2)к = 1/2Ф(ак/2)(^ок + д\) - Q(ak1/2 + 1Ч/2Н/2, дк = (1+ ике*)Ит11ег(акц/2,ак-1/2),
к = Г 1/2Ф(а1/2)(дк - дк)/а\/2, а\/2 = 0
71/2 = \ 0, а*/2 = 0,
Ф(г) = Q(z) - Лг2, ак/2 = Ь^СИ - Ио), ек = |ак/2 - а-1/2|/(|ак/2| + |а-1/2|), ^ =
где € [0, 2] — параметры искусственного сжатия, которые могут принимать различные значения для каждого характеристического переменного. Обычно применялись значения ик = 2, однако для более точных расчетов в течениях, содержащих зоны разрежения, все эти параметры полагались нулевыми. Параметр Л определяет точность схемы по времени: при Л = 1 — схемы второго порядка точности по времени, при Л = 0 — первого порядка. В этом случае схема используется для расчета стационарных течений методом установления.
Использовались операторы — ограничители двух типов: оператор, предложенный Хартеном:
Нт^ег(ж,у) = ттто^я, у) = sign(ж) тах[0, тт(|ж|,у sign(x)], (6)
и оператор, предложенный Рое [3] (БирегЬее):
Нш^ег(ж, у) = ттто^ттто^2ж, у), ттто^я, 2у)]. (7)
Расчетные сетки строились по алгоритму [6], подробное описание реализации дано в [5].
2. Результаты расчетов
Для расчетов в уравнениях (1) был проведен переход к безразмерным переменным, значения масштабных коэффициентов выбирались в соответствии с [2]. После перехода к безразмерным переменным значения констант в (1) принимают значения, соответствующие стехиометрической смеси водород—кислород—азот: 20, 10 и 70 % соответственно: д = 2, Е1 = 1.7, Е2 = Е 1/4.9, К1 = 3.6, К2 = К/2.
2.1. Расчет цикла работы пульсирующего детонационного двигателя
Рассчитывалось течение, соответствующее такту работы пульсирующего детонационного двигателя [7, 8]. Схема двигателя взята из [8] и показана на рис. 1. Газовая смесь, содержащая горючий компонент (водород или ацетилен), подготавливается в реакторе (рис. 1, а, 1). Подготовка заключается в подогреве и разложении молекул горючей смеси на химически активные составляющие. В реактор смесь поступает из топливного смесителя (рис. 1, а, 4), сжатый воздух к которому поступает от компрессора. Из реактора
Рис. 1. Схема устройства пульсационного детонационного двигателя: а — газогенератор с тяговыми модулями; б — блок резонатора; 1 — реактор; 2 — резонатор; 3 — кольцевое сопло; 4 — топливный смеситель
подготовленная газовая смесь поступает в блок резонатора через кольцевое сопло (на рис. 1, а, 3) в виде кольцевой сверхзвуковой струи. Конструкция самого импульсного детонационного двигателя ( рис. 1, б) состоит из резонатора в форме полусферы, в который поступает топливная смесь в виде кольцевой сверхзвуковой струи. По замыслу авторов [8], выходящая из кольцевого сопла сверхзвуковая струя должна фокусироваться в окрестности центра сферического резонатора, порождая ударную волну, распространяющуюся к стенкам камеры-резонатора. Отраженная от стенок сферы ударная волна в свою очередь фокусируется около центра сферы, порождая детонацию воздушно-водородной смеси. Возникшая детонационная волна, которая опять распространяется к стенкам резонатора, отражается от него и образует отраженную ударную волну. Эта волна, двигаясь направо по продуктам сгорания, разрывает струйную завесу, образуемую струей из кольцевого сопла. Прореагировавшая газовая смесь выбрасывается из полости резонатора со сверхзвуковой скоростью, одновременно начинается новый цикл работы двигателя.
Расчетная сетка (9600 узлов), на которой рассчитывалось течение в резонаторе, изображена на рис. 2, 11. Выход кольцевого сопла — А, край резонатора — В. Параметры, определяющие течение, сводятся (не считая концентрации горючей смеси) к отношению давлений и температур в реакторе и внешнем пространстве и отношению выходного и критического сечений кольцевого сопла. Начальные условия газодинамических параметров в резонаторе соответствуют атмосферным значениям. На стенке резонатора (верхняя граница расчетной области до точки В на рис. 2, 11) задаются условия непротекания: Уп = 0, д//дп = 0, / = Р, V, а, в, Р на выходе — известные условия вытекания [5], на оси симметрии — условия симметрии, в выходном сечении кольцевого сопла А (на рис. 2, 11) газодинамические параметры определяются как решение задачи о сопле Лаваля [9] по значениям давления и температуры в реакторе и отношению выходного и критического сечений сопла. На рис. 2, 1-10 изображены линии уровня параметра в (черным цветом), характеризующего выгорание горючей смеси, и изолинии плотности газовой смеси (светло-серый цвет). В приведенном варианте расчета были выбраны параметры: Рг/Р0 = 4, Тг/Т0 = 3, Sk/Sexit = 0.25. На рис. 2, 1-3 видны выход
Рис. 2. Последовательные стадии колебательного течения в пульсационном двигателе
струи горючей смеси и образование на ее границе вихревых зон несгоревшей смеси. На рис. 2, 4 наблюдаются инициация детонации в вихревых зонах на границах выходящей струи и формирование вихревой зоны повышенного давления в донной части резонатора. На рис. 2, 5 показаны полное выгорание горючей смеси в струе и первый выброс продуктов сгорания сквозь кольцевую газовую струю. На рис. 2, 6 кольцевая струя формируется вновь, выбрасываемая из донной части резонатора струя продуктов горения утончается. На рис. 2, 7 вновь происходит запирание продуктов сгорания в донной части резонатора. На рис. 2, 8 кольцевая струя вновь выгорает и образуются вихри продуктов сгорания. На рис. 2, 9 и 10 наблюдается повторный выброс струи продуктов сгорания. Таким образом, приведенные картины изолиний плотности и массовой концентрации показывают, что при выборе параметров течения, близких к указанным в [8], наблюдается пульсационный характер течения, однако детонация горючей смеси происходит не в окрестности центра сферы, а в вихревых структурах выходящей кольцевой струи. Такой механизм возникновения детонации представляется более реалистичным, так как существенно нестационарный и пульсационный характер течения позволяет предполагать отсутствие в течении такой степени симметрии, которая привела бы к точной фокусировке ударной волны в центре сферы.
2.2. Расчет горения и детонации водородсодержащей газовой смеси под оболочкой реакторного зала
Проведено численное моделирование различных режимов горения и детонации смеси воздух—водород под защитной оболочкой (контейнментом) реакторного зала АЭС для реакторов типа ВВР. После аварии на Три-Майл Айленд, где имело место горение водорода, и аварии на Чернобыльской АЭС для обеспечения безопасности работы АЭС уделяется внимание даже маловероятным ситуациям, в том числе связанным с проблемой выделения водорода (хотя существуют достаточно надежные методы борьбы с накоплением водорода под оболочкой реакторного зала). На рис. 3 приведена схема помещений реактора из [10]. Расчетная область представляет собой цилиндр (реакторная шахта), сопряженный с полусферой (купол контейнмента), без служебных помещений. Режимы самопроизвольного возгорания и детонации водорода моделировались частью расчетной области, в которой в нулевой момент времени происходит моментальное выгорание водорода с повышением температуры (горение) и возрастанием нормальной к границе раздела компоненте скорости до значений, соответствующих волне Чепмена—Жуге (детонация). Целью расчета являлась оценка пиковых значений давления на стенках защитной оболочки для прогноза ее разрушения. В соответствии с этим оценивалось превышение давления вдоль купола контейнмента в различные моменты времени. Проектный запас прочности конструкции предполагает сохранение целостности конструкции до превышения значения давления в 30-40 раз [10]. На рис. 4, 7 изображена расчетная сетка (9600 узлов), правая часть которой соответствует куполу контейнмента, а левая — шахте реактора. Были рассчитаны варианты возникновения и распространения детонационной волны в смеси водород—воздух, близкой к стехио-метрической со значениями констант (8) для системы уравнений (1), дополняющие результаты [10]. Начальные условия возникавшего очага детонации задавались точными решениями плоской волны Чепмена—Жуге, распространявшейся в трубе с закрытым
Рис. 3. Схема реакторных помещений
Рис. 4. Инициация детонации в верхней части купола
концом, и сферической волны Чепмена—Жуге [9]. Визуализация течения производится, как и для первой задачи, при помощи линий уровня параметра, соответствующего массовой концентрации водорода (черные) и изолиний плотности (светло-серые).
Инициация детонации в верхней части купола контейнмента. Расчет соответствует периоду времени, достаточно удаленному от момента аварии, когда выделявшийся при нагревании элементов реактора водород проник из шахты реактора в купол и распространился во всей защитной оболочке (купол и шахта реактора). На рис. 4, 1-5 даны картины изолиний плотности и массовой концентрации в последовательные моменты времени для инициации детонации в верхней части купола контейнмента. На рис. 4, 6 показано распределение давления вдоль купола от его верхней точки (правая нижняя точка расчетной области), соответствующей нулю оси абсцисс графика, к нижней точке свода, соответствующей максимальному значению абсциссы графика. Для сравнения на графике приведено значение давления в непрореагировавшем газе — горизонталь Р0. Максимальное превышение значений давления примерно в 10-12 раз наблюдается в высшей точке купола, что соответствует результатам, полученным в [10].
Инициация детонации в середине купола контейнмента при заполнении горючей смесью верхней части оболочки. Этот вариант расчета соответствует периоду времени, когда водород проник из шахты реактора в купол и успел заполнить только верхнюю часть купола. На рис. 5, 1-4 даны картины изолиний плотности и массовой концентрации в последовательные моменты времени, когда горючая
Давление вдоль купола
Рис. 5. Инициация детонации в середине купола контейнмента при заполнении горючей смесью верхней части оболочки
смесь заполнила верхнюю часть купола контейнмента. На границе раздела воздух— горючая смесь заданы условия равенства давлений. На рис. 5, 5 показано распределение давление вдоль купола от его центра (правая нижняя точка расчетной области) — нуль оси абсцисс, к нижней точке свода — максимальное значение абсциссы графика. Для сравнения на графике приведено значение давления в непрореагировавшем газе — горизонталь Р0. Максимальное превышение значений давления примерно в 8 раз наблюдается в высшей точке купола, уменьшение значения пика давления обусловлено большим расстоянием от центра инициации и ослаблением детонационной волны при распространении в шахту реактора, не содержащую горючей смеси.
Инициация детонации в середине шахты реактора при заполнении горючей смесью всего объема защитной оболочки. Этот вариант расчета соответствует моменту времени, когда водород проник из шахты реактора в купол и заполнил всю оболочку, однако инициация произошла не в верхней части, а в шахте реактора. На рис. 6, 1-9 даны картины изолиний плотности и массовой концентрации в последовательные моменты времени. Образующаяся в этом случае сложная структура взаимодействия ударных и детонационных волн становится причиной частичного гашения пиковых нагрузок. На рис. 6, 10 приведено распределение давления вдоль купола в различные моменты времени, наблюдаемое превышение давления — в 8-9 раз.
Инициация детонации во всей шахте реактора при заполнении горючей смесью шахты реактора. Этот вариант расчета соответствует начальному периоду аварии, когда водород еще не проник из шахты реактора в купол. На рис. 7, 1-6 приводятся картины изолиний плотности и массовой концентрации в последовательные моменты времени. На рис. 7, 7 показано распределение давления вдоль купола в различные моменты времени, цифрами обозначены графики давления вдоль купола
для картин изолиний 4-6. При этом варианте расчета превышение давления составляет 30 раз и для такого сценария разрушение купола возможно. В целом структура течения соответствует задаче о выходе детонационной волны из трубы с узким сечением в более широкую полость [9]. Наблюдается ослабление детонационной волны и опережение ее ударной волной. Взаимодействие отраженной от купола ударной волны и падающей детонационной ведет к сложной неустойчивой картине взаимодействия разрывов, отраженной на рис. 7, 5 и 6. Численное моделирование такого течения проведено, в частности, в [11], и там же даны сходные картины изолиний параметров течения.
з. Обсуждение результатов
Модель двухстадийной химической реакции является весьма грубым приближением
и, по-видимому, плохо описывает химическую структуру детонации, однако в сочетании с использованием разностных схем повышенного порядка точности вполне успешно применяется для моделирования течений со сложной структурой разрывов [11, 12] и
Рис. 7. Инициация детонации во всей шахте реактора при заполнении горючей смесью шахты реактора
достаточно широко используется для моделирования качественной структуры течения. Применение в численном моделировании модели, включающей систему уравнений химических реакций для большого числа компонентов смеси, по-видимому, лучше описывает химические процессы горения и детонации, однако в сочетании со схемой Годунова [13] или схемой крупных частиц [10] не может дать хорошего разрешения сложных структур разрывов течений. Приведенные результаты расчетов показывают работоспособность методики для широкого диапазона параметров идеального реагирующего газа. Разрешающая способность используемой схемы Хартена позволяет наблюдать четкую структуру образующихся детонационных и ударных волн и оценивать величину пиковых значений газодинамических величин.
Список литературы
[1] Левин В.А., Марков В.В. Возникновение детонации при концентрированном подводе энергии // Физика горения и взрыва. 1975. Т. 11, № 4. С. 623-633.
[2] Таки С., Фудзивара Т. Численный анализ двумерных нестационарных детонационных волн // Ракетная техника и космонавтика. 1975. Т. 11, № 1. С. 93-98.
[3] Roe P.L. Approximate Riemann solvers, parameter vectors and difference scheme //J. Comp. Phys. 1981. Vol. 43. P. 357-372.
[4] Yee H.C., Warming R.E., Harten A. Implicit total variation diminishing (TVD) schemes for steady-state calculation //J. Comp. Phys. 1985. Vol. 57. P. 327-361.
[5] МАртюшов С.Н. Расчет двух нестационарных задач дифракции явным алгоритмом второго порядка точности // Вычисл. технологии. 1996. Т. 1, № 4. С. 82-89.
[6] Thompson J.F., Warsi Z.U.A., Mastin C.W. Numerical Grid Generation. N.Y.: North Holland, 1985. 306 p.
[7] Левин В.А., Нечаев Ю.И. Тарасов А.И. Новый подход к организации рабочего процесса пульсирующих детонационных двигателей // Химическая физика. 2001. Т. 20, № 6. С. 90-98.
[8] Марчуков Е., Нечаев Ю., Полев А., Тарасов А. Пульсирующие детонационные двигатели // Двигатель. 2003. Т. 1. С. 14-17.
[9] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. М., 2001. Т. 6. Гидродинамика. 736 c.
[10] Гальвурт В.А., Иванов М.Ф., Минеев В.Н. и др. Воздействие взрыва водорода на защитную оболочку реакторного зала АЭС // Математическое моделирование. 2002. Т. 14, № 1. С. 73-86.
[11] Pantov E.G., Fisher M., Kratzel T. Decoupling and recouhling of detonation waves associated with sudden expansion // Shock Waves. 1996. Vol. 6. P. 131-137.
[12] Jones D.A., Kemister G., Sichel M., Oran E.S. The influence of cellular structure on detonation transmission // Shock Waves. 1996. Vol. 6. P. 119-130.
[13] Левин В.А., Марков В.В., Осинкин С.Ф. Инициирование детонации в водородно-воздушной смеси взрывом сферического заряда ТНТ // Физика горения и взрыва. 1995. Т. 31, № 2. С. 91-95.
Поступила в редакцию 13 января 2005 г., в переработанном виде —15 июня 2006 г.
