Математическое моделирование детонационной волны в системе координат, связанной с лидирующим скачком тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Лопато Александр Игоревич

• Лопато Александр Игоревич
• кандидат наук
• 2017

Оглавление диссертации кандидат наук Лопато Александр Игоревич

ВВЕДЕНИЕ

ГЛАВА 1. ОБЗОР ОСНОВНЫХ РЕЗУЛЬТАТОВ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ОДНОМЕРНОЙ ПУЛЬСИРУЮЩЕЙ ВОЛНЫ ГАЗОВОЙ ДЕТОНАЦИИ

1.1. Линейный анализ устойчивости

1.2. Модель идеальной детонации

1.3. Описание детонационной волны в рамках двухстадийной модели кинетики

1.4. Описание детонационной волны в рамках одностадийной модели кинетики

1.5. Исследования волны детонации в системе координат, связанной с

лидирующей волной

ГЛАВА 2. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВОЛНЫ ДЕТОНАЦИИ В ЛАБОРАТОРНОЙ СИСТЕМЕ КООРДИНАТ

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

2.2. Математическая модель

2.3. Вычислительный алгоритм

2.4. ЕКО-реконструкция

2.5. Тестирование вычислительного алгоритма

2.5.1. Решение линейного уравнения переноса

2.5.2. Задачи Римана

2.5.3. Тест Шу

2.6. Результаты вычислительных экспериментов

2.6.1. Распространение детонационной волны в модельной ацетилено-воздушной смеси

2.6.2. Распространение детонационной волны в модельной водородно-воздушной смеси

2.6.3. Общий случай распространения устойчивой, слабо-неустойчивой,

нерегулярной и сильно неустойчивой детонационной волны

ГЛАВА 3. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВОЛНЫ ДЕТОНАЦИИ В СИСТЕМЕ КООРДИНАТ, СВЯЗАННОЙ С ЛИДИРУЮЩИМ СКАЧКОМ

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

3.2. Математическая модель

3.3. Вычислительный алгоритм

3.4. Решение Зельдовича-Неймана-Деринга

3.5. Тестирование вычислительного алгоритма

3.6. Результаты вычислительных экспериментов

ГЛАВА 4. СРАВНИТЕЛЬНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ ДВУХ

ПОСТАНОВОК

ГЛАВА 5. АНАЛИЗ ЭФФЕКТИВНОСТИ РАСПАРАЛЛЕЛИВАНИЯ

ВЫЧИСЛИТЕЛЬНОГО АЛГОРИТМА

ЗАКЛЮЧЕНИЕ

СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ

СПИСОК ЛИТЕРАТУРЫ

Введение диссертации (часть автореферата) на тему «Математическое моделирование детонационной волны в системе координат, связанной с лидирующим скачком»

ВВЕДЕНИЕ

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

Как известно из многочисленных экспериментальных и численных исследований, распространение ДВ в пространстве характеризуется формированием сложного нелинейного колебательного процесса, см. терминологию, например, в [1, 2], проявляющегося наличием:

- пульсаций параметров за фронтом ДВ в одномерных численных исследованиях, первые исследования подобного рода - [3, 4];

- поперечных волн сжатия, взаимодействующих с ЛВ, в двумерных расчетах и экспериментах по распространению ДВ в узких зазорах, основополагающие экспериментальные исследования - [5], численные - [6, 7];

- поперечной волны, двигающейся по спирали - спина, в трехмерном случае, основополагающие экспериментальные исследования - [1], численные - [8].

Неустойчивость фронта ДВ связана с сильной взаимосвязью между скоростью протекания химических реакций и интенсивностью ЛВ. Действительно, скорость протекания химических реакций со зависит от температуры газа непосредственно за фронтом ЛВ - температуры фон Неймана (ФН) [9] по закону

Аррениуса:

Г Е Л

со ~ ехр--.

Безразмерная энергия активации смеси Б/ () является, таким образом, основным определяющим параметром в теоретических исследованиях устойчивости распространения ДВ.

Актуальность темы исследования и степень ее разработанности. Множество численных исследований посвящено изучению механизмов распространения пульсирующей волны газовой детонации. Как правило, авторы используют постановку со стационарной ДВ, распространяющейся против набегающего потока, и методы сквозного счета, см., например, ссылки в [10]. Подробный анализ данных исследований представлен в Главе 1. Без сомнения, на сегодняшний день использование методов сквозного счета является безальтернативным при решении крупномасштабных многомерных проблем распространения волн горения в каналах и трубах.

Вместе с тем, как показано в [11] с использованием детонационных аналогий, разрешение зоны протекания химических реакций при распространении пульсирующей ДВ даже в одномерной постановке является на практике трудной задачей, поскольку масштаб зоны реакции меняется на порядки в процессе пульсаций в связи с возникновением локальных взрывов за фронтом ЛВ. Методы сквозного счета даже номинально высокого порядка аппроксимации (ПА) «размазывают» ЛВ, что может приводить к качественному изменению характера распространения ДВ. Кроме того, при использовании методов сквозного счета существует принципиальное ограничение на возможность определения мгновенной скорости ДВ с точностью выше, чем Ах/ Аt, где Ах - размер расчетной ячейки, А1 - шаг интегрирования по времени.

Следовательно, при изучении вопросов устойчивости ДВ и механизмов ее распространения существенным преимуществом обладают постановки, частично или полностью свободные от перечисленных выше недостатков. Оригинальным подходом, в котором решаются перечисленные выше проблемы, является переход в систему координат, связанную с ЛВ [12, 13]. При этом ЛВ становится

одной из неподвижных границ расчетной области, однако, приходится дополнительно интегрировать уравнение для эволюции скорости ЛВ.

Актуальность диссертационного исследования обусловлена рядом не до конца объясненных и понятных особенностей результатов вычислительных экспериментов по распространению волн детонации. К числу подобных особенностей можно отнести возможное затухание сильно неустойчивых ДВ при длительном распространении в одно- и двумерных постановках [14, 15]. Актуальность разработки высокоточной методики расчета пульсирующей ДВ связана также с рядом задач теории детонации, которые до настоящего времени не имеют окончательного решения, включая классификацию режимов распространения пульсирующей ДВ с использованием детальных моделей кинетики химических реакций, а также пульсирующей ДВ с потерями.

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

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

2. Разработка комплекса программ для численного исследования длительного распространения пульсирующей ДВ в одномерной постановке в ЛСК и в СКФ.

3. Проведение параметрических численных исследований характера пульсаций устойчивой, слабо неустойчивой, нерегулярной и сильно неустойчи-

вой ДВ при длительном распространении в одномерной постановке в ЛСК и в СКФ.

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

Теоретическая и практическая значимость результатов работы обусловлена результатами сравнительного анализа характера пульсаций ДВ при ее длительном распространении при расчетах в ЛСК и СКФ, что дает представление о масштабе неустранимой погрешности методов сквозного счета при решении подобного класса задач. Стоит отметить, что полученные результаты имеют отношение в большей степени к теории, чем к практике применения явления детонации.

Методология и методы исследования. В работе использовались методы математического моделирования и проведения вычислительных экспериментов на многопроцессорных ЭВМ, а именно: метод расщепления по физическим процессам, конечно-объемный метод Куранта-Изаксона-Рис расчета потоков, метод декомпозиции расчетной области.

На защиту выносятся следующие положения:

1. Вычислительный алгоритм второго порядка аппроксимации для исследования распространения пульсирующей ДВ в СКФ.

2. Комплекс программ для численного исследования длительного распространения пульсирующей ДВ в одномерной постановке в лабораторной системе координат и в системе координат, связанной с фронтом лидирующей волны. Программный комплекс адаптирован для проведения расчетов на современных многопроцессорных ЭВМ.

3. Сравнительный анализ характера пульсаций определяющих параметров при длительном распространении устойчивой, слабо неустойчивой, нерегулярной и сильно неустойчивой ДВ в ЛСК и СКФ.

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

Основные результаты диссертационной работы докладывались и обсуждались на следующих конференциях: 40-ая международная молодежная научная конференция «Гагаринские чтения» (Москва, 7 - 11 апреля 2014 г.); 9-ый Международный коллоквиум по пульсирующей и непрерывной детонации (Пушкин, 19 - 23 мая 2014 г.); Международная молодежная конференция «Современные проблемы прикладной математики и информатики» (Дубна, 25 - 29 августа 2014 г.); Ежегодные научные конференции отдела Горения и взрыва Института химической физики им. Н.Н. Семенова РАН (Москва, 11 - 13 февраля 2015 г., 10 - 12 февраля 2016 г.); 25-ый Международный коллоквиум по динамике взрывов и реагирующих систем (Лидс, Великобритания, 2 - 7 августа 2015 г.); XI Всероссийский съезд по фундаментальным проблемам теоретической и прикладной механики (Казань, 20 - 24 августа 2015 г.); XXII Международная конференция «Нелинейные задачи теории гидродинамической устойчивости и турбулентность» (Звенигород, 14 - 21 февраля 2016 г.); 10-ый Международный коллоквиум по пульсирующей и непрерывной детонации (Санкт-Петербург, 4 - 8 июля 2016 г.); Квазилинейные уравнения, обратные задачи и их приложения (Долгопрудный, 12 - 15 сентября 2016 г.); 59-ая Научная конференция МФТИ (Долгопрудный, 21 - 26 ноября 2016 г.).

Публикации автора по теме диссертационного исследования. Основные результаты опубликованы в 15 работах, включая 5 статьи в изданиях из Перечня ВАК ведущих рецензируемых научных журналов и изданий, в которых должны быть опубликованы основные результаты диссертации на соискание степеней кандидата наук (1 - 5), 1 статью в зарубежном издании, входящем в систему Web of Science (5), 3 публикации в сборниках статей (6 - 8), 6 публикаций в трудах конференций (9 - 14), 1 учебно-методическое пособие (15):

1. Лопато, А.И., Уткин, П.С. Математическое моделирование пульсирующей волны детонации с использованием ENO-схем различных порядков аппроксимации // Компьютерные исследования и моделирование. - 2014. - Т. 6, № 5. - С. 643 - 653.

2. Лопато, А.И., Уткин, П.С. Особенности расчета детонационной волны с использованием схем различных порядков аппроксимации // Математическое моделирование. - 2015. - Т. 27, № 7. - С. 75 - 79.

3. Лопато, А.И., Уткин, П.С. О двух подходах к математическому моделированию детонационной волны // Математическое моделирование. - 2016. - Т. 28, № 2. - С. 133 - 145.

4. Лопато, А.И., Уткин, П.С. Детальное математическое моделирование пульсирующей детонационной волны в системе координат, связанной с лидирующим скачком // Журнал вычислительной математики и математической физики. - 2016. - Т. 56, № 5. - С. 856 - 868.

5. Lopato, A.I., Utkin, P.S. Towards second-order algorithm for the pulsating detonation wave modeling in the shock-attached frame // Combustion Science and Technology. - 2016. - V. 188, No. 11 - 12. - P. 1844 - 1856.

6. Lopato, A.I., Utkin, P.S. Mathematical modeling of pulsating detonation wave propagation using monotone numerical methods of different approximation orders // Transient Combustion and Detonation Phenomena: Fundamentals and Applications. - Moscow: Torus Press, 2014. - P. 261 - 268.

7. Лопато, А.И., Уткин, П.С. Исследование пульсирующей волны детонации методами сквозного счета и в системе координат, связанной с лидирующей волной // Горение и взрыв. - 2015. - Т. 8, № 1. - С. 145 - 150.

8. Lopato, A.I., Utkin, P.S. Mathematical modeling of weakly unstable pulsating detonation wave propagation: laboratory versus shock-attached frame // Progress in Detonation Physics. - Moscow: Torus Press, 2016. - P. 100 - 105.

9. Лопато, А.И., Уткин, П.С. Математическое моделирование детонационной волны с использованием ENO-схем различных порядков аппроксимации // Научные труды Международной научной конференции XL Гагаринские чтения. Москва, 7 - 11 апреля 2014 г. - Т. 5. - М.: МАТИ, 2014. - С. 137 - 139.

10. Lopato, A.I., Utkin, P.S. Mathematical modeling of pulsating detonation wave using non-oscillating schemes of different approximation orders // Intern. Conf. Modern Probl. of Appl. Math. & Com. Science (MPAMCS'2014): Book of Abstracts. Dubna, August 25 - 29, 2014. - Dubna: JINR. - P. 164 - 168.

11. Lopato, A.I., Utkin, P.S. Towards second-order algorithm for the pulsating detonation wave modeling in the shock-attached frame // Proc. 25th Intern. Colloq. on the Dyn. of Expl. and React. Syst., 2 - 7 August 2015, Leeds, UK. - Paper 46.

12. Лопато, А.И. Математическое моделирование пульсирующей волны детонации в системе координат, связанной с лидирующей волной // Сб. докл. XI Всеросс. съезда по фунд. пробл. теор. и прикл. механики, 20 - 24 августа 2015 г., Казань. - Казань: Изд-во Каз. ун-та, 2015. - С. 2369 - 2372.

13. Лопато, А.И., Уткин, П.С. Исследование устойчивости распространения волны газовой детонации в двух постановках // Матер. XXII Межд. конф. «Нелин. зад. теор. гидрод. устойч. и турбул.». 14 - 21 февраля 2016 г., Звенигород. - М.: Изд-во МГУ, 2016. - Электр. изд. - C. 111 - 117.

14. Лопато, А.И., Уткин, П.С. Математическое моделирование течений с ударными и детонационными волнами с использованием высокопроизводительных вычислений // Труды 59-ой научной конференции МФТИ. Управление

и прикладная математика. - М.: МФТИ, 2016. - 2 с. URL: http://conf59.mipt.ru/ static/reports_pdf/3060.pdf.

15. Уткин, П.С., Лопато, А.И. Математическое моделирование детонационных волн: уч.-метод. пособие. - М.: МФТИ, 2017. - 24 с.

Личный вклад автора в публикациях с соавторами. Все работы, кроме одной, выполнены в соавторстве с научным руководителем. Личный вклад автора - разработка вычислительного алгоритма для расчета ДВ в СКФ, разработка соответствующего комплекса программ, проведение вычислительных экспериментов и анализ результатов. Основной вклад научного руководителя -в постановке задачи и интерпретации некоторых полученных результатов.

Глава 1. Обзор основных результатов математического моделирования распространения одномерной пульсирующей волны газовой детонации.

1.1. Линейный анализ устойчивости

Обзор результатов численных исследований нелинейной динамики детонационных волн уместно начать с основополагающей теоретической работы Эрпенбека (Erpenbeck) [16], в которой впервые был осуществлен линейный анализ устойчивости стационарного решения Зельдовича-Неймана-Деринга (ЗНД) [9] одномерных уравнений Эйлера, дополненных одностадийной моделью кинетики химических реакций. Анализ сводится к изучению поведения решения системы уравнений, линеаризованной вблизи решения ЗНД, при наложении малых возмущений. Отметим основные этапы данного анализа.

Математическая модель системы Эйлера для ЛСК записывается в следующем виде:

du df

— + — = f,

dt dx

(1.1)

р р v " 0 "

р v 4» pv2 + p 0

u = г , f = (p+г )v , s = -p Qcf

р Z _ p vZ ffOf

гp2s=p

, p = p—T, cf = -AZ exp /f

R T

г-* /V _£» ^

Здесь и далее и - вектор консервативных переменных, I - вектор потоков, 8 -

/V /у

вектор источников, х - время и координата в ЛСК соответственно, р -плотность газовой смеси, V - скорость смеси в ЛСК, е - полная энергия единицы объема смеси, Z - массовая доля реагирующего компонента смеси, р -

давление смеси, () - тепловой эффект химической реакции смеси, со - скорость

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

/V /V

ты смеси, R - универсальная газовая постоянная, // - молярная масса смеси, T

/V /V

- температура смеси, A - предэкспоненциальный множитель, E - энергия активации смеси.

После составления определяющей системы уравнений проводится процедура приведения ее к безразмерному виду. В качестве масштабов обезразмери-вания Эрпенбеком в [16] выбираются параметры газа перед скачком. Будем отмечать их индексом "0". Масштабом длины служит половина «длины полупревращения» (подробнее о ней в следующей главе). В дальнейшем, символы без верхней крышки считаются безразмерными.

Для исследования линеаризованной системы перейдем в СКФ:

% = x - Dt -у( t), т = t.

Здесь % - координата в СКФ, D - скорость ДВ, у - локус фронта волны. Линеаризация уравнений предполагает запись решения V = (р V p 2) определяющей системы, как суммы решения стационарной задачи V*(%) и возмущений:

V (%,т) = у(%) + V '(£т). (1.2)

Подстановка (1.2) в определяющую систему (1.1) приводит после математических преобразований к системе уравнений в частных производных для определения V'(%,т). Для дальнейшего исследования производится преобразование

Лапласа для величины у(т):

¥(т) = (2П)-1 \

и вектора V'(%,т). В результате преобразований анализ сводится к решению

системы обыкновенных дифференциальных уравнений. В соответствии с теорией дифференциальных уравнений, решение системы определяется комбина-

13

цией ее собственных векторов, умноженных на экспоненциальные множители, содержащие собственные числа. Выполнение условия ограниченности V'(%,т)

позволяет определить вид функции С(п). Это дробная функция комплексной переменной ц. Дальнейший анализ сводится к нахождению полюсов полученной дробной функции комплексной переменной ^(п), что эквивалентно нахождению корней её знаменателя V(ц): V(ц) = 0. По построению, функция С(п) является, в соответствии с преобразованием Лапласа, изображением у(тт. Поэтому, в случае положительной действительной части корня ц: Яе(ц)> 0 возмущение V'(%,т) неограниченно возрастает при ц ^да .

Рис. 1. Граница устойчивости решения линеаризованной системы (1.1) в плоскости ((, Q) при фиксированном значении показателя адиабаты у = 1.2 и энергии активации Е = 50 [16]. Цифры в области неустойчивости решения указывают на число корней V (ц) = 0 с положительной действительной частью.

14

Если же все корни имеют отрицательную действительную часть Яе (п)< 0, то

наложенные возмущения со временем подавляются, и решение носит устойчивый характер. Исследование значений полюсов при различных параметрах химической смеси позволяет определить границу устойчивости решения линеаризованной системы. Так, фиксируя показатель адиабаты у = 1.2 и энергию активации Е = 50, при варьировании теплового эффекта Q, степени пересжатия

/ = (О / ОС] )2 пересжатой ДВ, которая поддерживается движением левой границы области (подвижный поршень) со скоростью О, а Ос; - скорость Чепме-на-Жуге (ЧЖ), Эрпенбеком получена в [16] представленная на Рис. 1 граница устойчивости линеаризованной системы в плоскости (/, Q).

Таким образом, преобразование Лапласа является важным математическим инструментом при исследованиях Эрпенбеком линеаризованной системы. Используя его, Эрпенбеку в 1960-х удалось получить неизвестные раннее результаты. Дальнейшее развитие математического аппарата позволило прийти к новому подходу для проведения линейного анализа, область применимости которого может быть шире при варьировании энергии активации модельной смеси, как отмечено в [17].

В [18] авторами был проведен линейный анализ одномерной системы Эйлера с кинетикой Аррениуса при использовании метода наложения нормальных мод возмущения

V (4,т) = у*(^) + у'(4,т), у'(£,г) = v'(¿f)exp (ат), ц = ц' ехр (ат),

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

должно удовлетворять условию отсутствия в конце зоны реакции возмущений, идущих в сторону фронта ДВ - «radiation condition». Выполнение данного условия приводит к соотношению, накладываемому на компоненты вектора решения:

H (а,у, E, Q, f, u') = 0.

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

H (а(у, E, Q, f, u')) = 0

задает линию устойчивости («neutral stability curve») решения определяющей системы уравнений. Стоит, однако, заметить, что уравнение может иметь несколько корней а или частот. Каждая частота определяет свою границу устойчивости, поэтому для определения конечной области устойчивости задачи нужно знать расположение крайней границы устойчивости («neutral stability boundary»).

10-1 ..............-■■■■...................1.........

0 10 20 30 40 50

Activation energy, Е

Рис. 2. Границы устойчивости, полученные в ходе линейного анализа, при разных частотах а для диаграммы (E, Q) [18], f = 1. Граница, за которой ДВ становится неустойчивой, отмечена цифрой 1.

Так, при фиксированном показателе адиабаты у = 1.2, степени пересжатия / = 1 можно исследовать решение на устойчивость при варьировании теплового эффекта химической реакции Q и энергии активации Е. Результаты данного

исследования, представленные в [18], имеют вид, изображенный на Рис. 2. На Рис. 2 также видно, что при энергии активации смеси Е < 14 решение независимо от эффекта реакции Q всегда устойчиво. Также заметим, что на Рис. 2 отмечено 3 линии устойчивости. При выборе (Е, Q) справа от границы неустойчивости V' решение неустойчиво. Причем, характер неустойчивости - число мод - определяется числом линий устойчивости слева от выбранной точки.

1.2. Модель идеальной детонации

Несмотря на определенные недостатки данного инструмента исследования, линейный анализ дает возможность определить границы области устойчивости распространения ДВ. Данную информацию часто используют, в качестве ориентира при исследованиях. Однако, для получения более подробных результатов моделирования ДВ, в том числе отвечающих натурным экспериментам, необходимы другие модели ДВ и методы их исследования.

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

где и = р 1 - удельный объем, Q - удельная теплота химической реакции, «0» - параметры перед скачком, «Н» - параметры за скачком. Данное уравнение определяют кривую Гюгонио (см. [19]) в плоскости (р,и), на которой лежат

все возможные состояния вещества за скачком. Кривая Гюгонио состоит из двух ветвей - часть адиабаты МС, точки которой представляют собой параметры продуктов детонации (сверхзвуковая скорость) и кривой БЕ, соответствующей дефлаграционному горению (дозвуковая скорость относительно исходного вещества) (см. Рис. 3).

Рис. 3. (р,и)-диаграмма процессов детонации и горения: ветвь НС соответствует пересжатой детонации, точка Н - детонации ЧЖ, НМ - недосжатой детонации, БЕ - дефлаграционному горению [19].

Вполне определенной скорость детонации оказывается в случае установившейся детонации, которая соответствует точке касания Н при начальном состоянии газовой смеси (р0,и0). В таком случае, скорость продуктов превращения относительно скачка Б - V равна местной скорости звука с. Если V + с > Б, то дето-

национная волна называется пересжатой. Установившаяся детонация, соответствующая точке касания Н, была названа детонацией Чепмена-Жуге (ЧЖ). Чепмен и Жуге предполагали, что реальный процесс самоподдерживающейся детонации отвечает именно тем параметрам (рН ,иН), которые соответствуют

точке касания Н. В ряде случаев это оказалось верно, но не во всех случаях.

Использование данной модели позволило получить ряд полезных результатов. Так, исследование задачи о точечном взрыве при данной модели ДВ было дано в работах Коробейникова, Левина, Бишимова, Маркова. Например, в работе [20] представлено приближенное решение задачи о точечном взрыве. Представлена оценка расстояния, на котором происходит переход пересжатой ДВ в волну ЧЖ. В приближении, что весь газ собирается в узкий слой за волной, записываются основные законы сохранения для слоя - массы, импульса, энергии. Путем математических преобразований с заменами переменных удается свести систему уравнений к дифференциальному уравнению для скорости ДВ О относительно радиуса волны (координаты фронта). Однако используемое в работе приближение перестает работать при приближении значения искомой скорости к скорости волны ЧЖ.

Для другого исследования [21], где рассматривалась плоская и цилиндрическая ДВ, отметим, что распределение параметров за цилиндрической волной на больших временах расчета практически совпадает с распределением в задаче о распространении детонации от центра в режиме ЧЖ. Тогда как, в плоском случае переход к режиму ЧЖ не осуществляется даже на очень большом расстоянии, но имеется асимптотическое стремление волны к волне ЧЖ.

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

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

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

1.3. Описание детонационной волны в рамках двухстадийной модели кинетики

Список литературы диссертационного исследования кандидат наук Лопато Александр Игоревич, 2017 год

- -

-

-

-

1 III

-20 -15 -10 -5 0

Рис. 56. Рассчитанные профили ЗНД плотности, скорости, давления, доли реагента, отнесенные к параметрам ФН. Случай устойчивой ДВ.

В качестве примера распределения газодинамических величин в модели ЗНД, приведем профили ЗНД, полученные для числа ячеек N = 4000, энергии активации Е = 25, теплового эффекта 0 = 50, показателя адиабаты у = 1.2. На Рис. 56 представлены зависимости плотности, скорости, давления и доли реагента, отнесенные к рассчитанным по (2.1) значениям параметров фон-Неймана: Р^ = 8.738523, у^ = 6.030227, р^ = 42.062677 .

Результаты, представленные в данном разделе, опубликованы в [70].

3.5. Тестирование вычислительного алгоритма

В качестве тестовой задачи, демонстрирующей работоспособность описанной выше методики, рассмотрим задачу из [12] о взаимодействии двух УВ, одна из которых догоняет другую. Стоит также отметить, что данная задача качественно соответствует ситуации соударения внутренней ДВ, возникающей в результате локального взрыва, с ЛВ, имеющей место при нерегулярном пульсирующем режиме распространения ДВ. Длина расчетной области составляет Ь = 20. В начальный момент времени отрезок [-Ь /3; 0] отвечает области за первой УВ, распространяющейся со скоростью В1 = 6 (параметры параметры за фронтом первой УВ - р1 = 8.25, у = 5.27, р1 = 32.64), отрезок [-Ь;-Ь / 3] - области за второй УВ, двигающейся со скоростью В2 = 12 (параметры за фронтом второй УВ - Р2 = 44.29, у2 = 10.75, р2 = 336.45. Фоновые параметры перед первой УВ - Р0 = 1.0, у0 = 0.0, р0 = 1.0. Показатель адиабаты у = 1.2. Расчет проводился на сетке с числом ячеек N = 500.

Рис. 57 иллюстрирует зависимость скорости фронта ЛВ от времени. Аналитическое решение данной задачи показывает, что вторая УВ догоняет первую УВ к моменту времени 1.11. Столкновение волн при таких заданных параметрах задачи приводит к распаду разрыва с конфигурацией, содержащей УВ, КР и волну разрежения (см. Рис. 58).

D

13 12 11 10 9 8 7 6

1-гггг - 1 ■ . 1 ,],,,, - -

ПА 1

: i ПА 2 Тонное решение

: 1

1

:

: |

11

J_U_L. 1

1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25

T

Рис. 57. Тестовая задача о взаимодействии двух УВ. Динамика изменения скорости ЛВ.

Рис. 58. Тестовая задача о взаимодействии двух УВ. Профиль плотности в момент времени 2.5.

3.6. Результаты вычислительных экспериментов

Перейдем к результатам распространения пульсирующей ДВ с использованием разработанного и описанного выше инструментария. По аналогии с постановкой в ЛСК зафиксируем показатель адиабаты у = 1.2 и тепловой эффект 0 = 50, а значит, и температуру фон Неймана Тум, и будем варьировать величину энергии активации Е. Газодинамические величины в невозмущенной области составляют значения р0 = 1, у0 = 0, р0 = 1. Величина давления в пике ФН при таких параметрах смеси составляет рш = 42.062677. Скорость самоподдерживающейся детонации ЧЖ для данных параметров равна:

РСГ =у]г + 2( -1)0 + ^2( -1) * 6.809475.

Для всех расчетов раздела 3.6 число СБЬ выбиралось равным 0.3, количество «внутренних» шагов интегрирования при расчете химической кинетики -К = 50, а заданная точность для метода Ньютона - д = 10-6.

Устойчивая детонация

Рассмотрим случай распространения ДВ в модельной смеси с энергией активации Е = 25 . При длине расчетной области Ь = 20 и общем числе ячеек N = 4000 на характерный пространственный масштаб в одну единицу, соответствующий в размерных величинах "длине полупревращения", приходится = 200, что совпадает с соответствующей величиной в постановке в ЛСК.

Распределения величин за фронтом качественно представляют собой классические профили ДВ уже в начале расчета. Однако, количественно они отличаются от теоретических профилей. В частности, значения в пике ФН и точке ЧЖ количественно отличаются от теоретических значений в соответствующих точках (см. Таб. 10). Начиная со значения т = 500, картина меняется. Газодинамические величины за фронтом ДВ выходят на определенный уровень, и численные значения на этом уровне практически не отличаются от теоретических. В резуль-

тате, реализуется устойчивый режим детонации. Пространственное распределение давления, приведенное к давлению в пике ФН, представлено на Рис. 59.

Р Р V т В

Параметр ЧЖ в расчете 22.20 1.84 3.07 12.07 6.81

Теоретически рассчитанный параметр ЧЖ 21.53 1.79 3.01 12.03 6.81

Параметр ФН в расчете 42.08 8.73 6.03 4.82 -

Теоретически рассчитанный параметр ФН 42.06 8.74 6.03 4.81 -

Таб. 10. Полученные при численном расчете и теоретических исследованиях параметры ЧЖ и ФН для модельной смеси с параметрами у = 1.2, 0 = 50, Е = 25. Устойчивый режим детонации.

Р

VN

1.1

0-9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

| . м мм МЫ 1111 г \

; ;

:

; ;

;

; ;

: ;

:

, , ,' < ■ - .

20 -17.5 -15 -12.5 -1 0 -7.5 5 -2 .5 0

Рис. 59. Пространственное распределение давления за фронтом ДВ в момент времени т = 550 для постановки СКФ. Устойчивый режим детонации Е = 25 . N = 4000, Ь = 20, Ы1/2 = 200.

Рассматривая динамику изменения скорости ДВ, приведенную на Рис. 60, можно увидеть наличие осцилляций. Как видно, осцилляции имеют затухающий характер, и искомая скорость выходит на постоянную величину. Более того, скорость ДВ стремится к теоретическому значению скорости детонации ЧЖ. График зависимости представляет затухающую синусоиду.

Рис. 60. Динамика изменения скорости ЛВ для устойчивой детонации Е = 25 в постановке СКФ. N = 4000, Ь = 20, ИУ2 = 200. Красной линией отмечена кривая из [12], зеленой - расчет по описанной методике второго порядка, горизонтальной пунктирной линией отмечен уровень В'^огу.

Выбранная геометрия расчетной области и сетки позволяет сравнить полученные зависимости с результатами из [12]. Однако, в работе [12] рассматривается не давление за фронтом ДВ, а скорость фронта. Поэтому здесь также приведена скорость ДВ, а не давление за фронтом. После затухания осцилляций установившееся значение скорости в [12] В = 6.917 заметно отличается от теоретиче-

ского В'^01 = 6.809. В собственном же расчете, установившаяся скорость волны с точностью до третьего знака после запятой совпадает с теоретическим значением. Длительные расчеты до момента времени 5000 единиц подтверждают выход скорости волны на обозначенный уровень.

Для случая устойчивой детонации, подобно [13], возможно оценить «реальные», то есть реализующиеся в практических расчетах, свойства вычислительного алгоритма, а именно его ПА. Для расчета на сетке с числом ячеек N = 4000 скорость детонации выходит в момент времени около т = 500 единиц на значение Всас1 * 6.809970, и таким образом, погрешность по сравнению с теоретическим значением составляет:

А1

т-лса1с1 у^1Неогу

0.000495.

Также был проведен аналогичный расчет устойчивой детонации на вдвое более грубой сетке с числом ячеек N = 2000. Получившееся стационарное значение скорости в этом случае составило Вса2 * 6.811152, и значит:

А 2

са1с 2 ^еогу

0.001677.

В соответствии с классическим определением аппроксимации:

= С -А^,

тл са1с тл Леогу

где константа С не зависит от шага сетки А^ , а д - порядок аппроксимации

численного метода. Тогда для реализующегося порядка аппроксимации имеем оценку:

д = М^А) * 1.76.

1п2

Измеренный порядок аппроксимации, приблизительно равный двум, подтверждает, что предложенная реализованная численная методика позволяет получить результаты со вторым порядком точности. Стоит отметить, что подобный

порядок аппроксимации, основанный на анализе скорости ЛВ, является недостижимым для методов сквозного счета.

Слабо неустойчивая детонация

Увеличение энергии активации модельной смеси до значения Е = 26 приводит к изменению характера распространения ДВ. Геометрические характеристики оставим прежними - Ь = 20, N = 4000, ^/2 = 200. Как и в случае расчета при энергии Е = 25 начальный этап расчетов характеризуется наличием низкочастотных пульсаций типа синусоиды. Однако теперь пульсации не затухают, а напротив, растут. Причем рост амплитуды пульсации продолжается не неограниченно, а до определенного значения, по достижении которого рост прекращается. В результате, при Е = 26 наблюдается слабо неустойчивая детонация с выходом на стационарный предельный цикл. Причем выход на цикл начинается непосредственно с начала расчета и занимает примерно 100 единиц времени.

Удобно сравнить динамику изменения давления за фронтом ДВ с результатами из [12]. Численная методика из [12] позволяет получить аналогичный предельный цикл, параметры которого несколько отличаются от рассчитанных в данной работе - период пульсаций предельного цикла в данной работе чуть выше (Ат = 11.814 против 11.345 в [12]) и ближе к значению, полученному в

линейном анализе в [18] - Ат1теаг = 12.11. Такая разница в периоде пульсаций обусловлена, в том числе, использованием в нашем подходе более точного алгоритма интегрирования уравнения эволюции скорости ЛВ. Амплитуда же пульсаций, полученная по нашей методике, имеет здесь более высокое значение. Выход на предельный цикл в расчетах происходит заметно быстрее (т = 110 единиц времени против 175 единиц).

Результаты, представленные на Рис. 61, получены при сравнительно небольшой длине канала Ь = 20 . Были исследовано влияние длины расчетной области на результаты моделирования. Рассмотрим распространение ДВ в канале

большей длины - Ь = 600 . Размер ячейки расчетной сетки оставлен неизменным, поэтому число ячеек будет N = 120 000.

Рис. 61. Динамика изменения скорости ЛВ для слабо неустойчивой детонации Е = 26 в постановке СКФ. N = 4000, Ь = 20, ^/2 = 200. Красной линией отмечена кривая из [12], зеленой - расчет по описанной методике второго порядка.

Проведенный расчет динамики изменения давления за фронтом ДВ при длине канала Ь = 600 , результаты которого представлены на Рис. 62, дает результаты, мало отличающиеся от полученных при малой длине. В частности, период пульсаций давления составляет до третьего знака после запятой такую же величину, как и период пульсаций скорости при Ь = 20 - Ат = 11.814. Данный факт дает основание полагать, что дальнейшее увеличение длины канала при сохранении размера ячейки расчетной сетки слабо повлияет на результаты для выбранных параметров модельной смеси.

Рис. 62. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Слабо неустойчивый режим детонации Е = 26. N = 120 000, Ь = 600 , ^/2 = 200.

Рис. 63. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Слабо неустойчивый режим детонации Е = 26. Ь = 20, ^/2 = 200. Сплошной линией обозначен график при числе ячеек N = 4000, пунктиром - при N = 8000.

Было исследовано влияние сеточного разрешения на результаты моделирования. Помимо расчета распространения ДВ при энергии активации смеси Е = 26 с числом ячеек N = 4000 проведем также расчет с удвоенным числом ячеек - N = 8000. Сравнивая динамику изменения давления за фронтом для двух расчетов, представленных на Рис. 63, видим, что уменьшение размера

ячейки расчетной сетки от величины А^ = 5 -10 3 до А^ = 2.5 -10 3 приводит к тому, что выход давления на предельный цикл происходит позже (т = 120 единиц времени против 150 единиц). Период пульсаций также меняется при изменении размера ячеек, но меняется достаточно слабо - 11.812 при размере ячейки А^ = 5 -10-3 против 11.814 при размере А^ = 2.5 -10-3.

Рис. 64. Пространственное распределение давления за фронтом ДВ в момент времени 550 для постановки СКФ. Слабо неустойчивый режим детонации Е = 26. N = 120 000, Ь = 600, ^/2 = 200.

Таким образом, изменение числа ячеек в расчетной области при энергии активации Е = 26 оказывает более значительный эффект на результаты, чем изменение длины канала. Однако, основная разница в результатах касается начального промежутка времени, связанного с выходом газодинамических величин на предельный цикл. После выхода на цикл разница в результатах оказывается менее существенной.

Рис. 64 иллюстрирует пространственный профиль давления при расчете слабо неустойчивой ДВ в СКФ.

Нерегулярная детонация

Дальнейшее увеличение энергии активации приводит к более нерегулярному осуществлению реакции горения реагента. А именно, времена задержек и периоды активного сгорания теперь варьируются без какого-либо регулярного характера, что приводит к картине распространения нерегулярной пульсирующей ДВ и, соответственно, нерегулярному режиму детонации. Динамика изменения давления за фронтом ДВ для энергии активации Е = 28 при длине канала Ь = 20, общем числе ячеек N = 8000, и числе ячеек на характерный масштаб длины ^/2 = 400 представлена на Рис. 65. Картина подтверждает для достаточно долгого времени расчета отсутствие выхода газодинамических величин на какой-либо уровень, а также промежутков с регулярными взрывами за фронтом. Вместо этого представленная картина характеризуется наличием низкочастотных пульсаций различного размера, со значениями давления при взрывах от 1.1 рш до 1.7рш . Стоит отметить, что хотя размеры низкочастотных пульсаций варьируются в широких пределах, но промежутки времени между пульсациями слабо меняются. Таким образом, нет основания ожидать нарушение режима или затухания ДВ. Отметим также, что динамика пульсаций давления, представленная на Рис. 65, качественно хорошо соотносится с пульсациями скорости ДВ для модельной смеси с такими же физико-химическими параметрами из [12].

1.8 -1-1-|-1-1—1—|—1-

1.6

О 6 -1—1—1-1—1-!-1-1—1—1—1-

О 100 200 300 400 500

Г

Рис. 65. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Нерегулярный режим детонации Е = 28. N = 8000, Ь = 20, N1/2 = 400.

Рассмотрим теперь, как изменится картина при увеличении длины канала до Ь = 600 с сохранением размера ячейки расчетной сетки. Как видно из Рис. 66, отражающего динамику изменения давления за фронтом ДВ при энергии активации модельной смеси N = 28 и длине канала Ь = 600 , изменение длины канала не приводит к изменению режима распространения ДВ. По-прежнему реализуется именно нерегулярный режим низкочастотных пульсаций с размерами пульсаций, лежащих в таком же промежутке. Более того, качественное сравнение Рис. 65 и Рис. 66 показывает, что в начальный промежуток расчетов, в 50 единиц времени, истории давления за фронтом ДВ полностью совпадают.

Рис. 66. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Нерегулярный режим детонации Е = 28. N = 240 000, Ь = 600 , = 400.

Совпадение обусловлено тем фактом, что первичные возмущения за фронтом ДВ, возникающие в начальные моменты времени, еще не успели преодолеть расстояние, равное длине малого канала. После того, как возмущения достигли границы малого канала необходимо еще время, чтобы отраженные волны достигли фронта ДВ. Начиная с этого момента, начинается расхождение в историях давления за фронтом ДВ для малого и длинного каналов. Также как и для короткого канала Ь = 20, для длинного канала Ь = 600 характерны пульсации с почти одинаковыми промежутками времени между ними без каких-либо сбоев или напротив, учащений пульсаций. Таким образом, увеличение длины канала не привело к качественным изменениям в характере пульсаций газодинамических величин за фронтом ДВ.

г ■ . I . I . ■ ■ . I . . . ■ I ■ . I I I ■ I ■_

О 50 100 150 200 250

Т

Рис. 67. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Нерегулярный режим детонации Е = 28. Ь = 20, ЫУ2 = 400. Сплошной линией обозначен график при числе ячеек N = 8000, точками - при N = 16000.

На Рис. 67 представлены результаты моделирования при измельчении расчетной сетки.

Помимо динамики изменения величин за фронтом ДВ, исследуем также пространственный профиль давления в момент времени 550 единиц. Рассмотрим профиль давления за фронтом ДВ, представленный для энергии активации Е = 28 на Рис. 68. Отметим, что хотя пик ФН и зона реакции соответствуют качественной картине, представленной, например, на Рис. 32, профиль давления за точкой ЧЖ характеризуется здесь наличием множества структур, представляющих собой волны с резкими перепадами давления, с гребнями и ложбинами. В первую очередь, именно резкие перепады и ассиметрия структур отлича-

ют профили давления при энергии Е = 28 от профилей давления при Е = 26. Такая ассиметрия обусловлена значительным нарушением регулярности при осуществлении химической реакции. Однако, несмотря на эту нерегулярность химической реакции, параметры модельной смеси оказываются достаточными для поддержания нерегулярного режима распространения ДВ, что показывают длительные расчеты в 5000 единиц времени.

Р

Рш

0.9

0.8

0.7

0.6

0-5

0.4

0.3

0.2

0-1

_' 1 1 М 11 1 >1 1 14' I | 1 М их 1М1 1 1 .

-

-Г 1 . 1 1 , . , I .... .... .... -

00 -450 -400 -350 -300 -250 -200 -150 -100 -50 Э

Рис. 68. Пространственное распределение давления за фронтом ДВ в момент времени т = 550 для постановки СКФ. Нерегулярный режим детонации Е = 28. N = 240 000, Ь = 600, Ы1/2 = 400.

Сильно неустойчивая детонация

Пульсации параметров за фронтом ДВ при рассмотренных энергиях активации происходили достаточно часто, и, кроме того, с определенной периодичностью. Увеличение энергии активации до значения Е = 35 значительно меняет характер пульсаций - они становятся редкими, промежуток между соседними

пульсациями меняется в широком диапазоне, чего не происходило при исследовании смеси при рассмотренных низких энергиях активации. Рассмотрим распространение ДВ при следующей геометрии: длина расчетной области L = 600, общее число ячеек N = 160 000, на характерный пространственный масштаб приходится N1/2 = 266 ячеек. Представленная на Рис. 69 динамика давления за фронтом ДВ демонстрирует относительно небольшое число пульсаций на промежуток времени 500 единиц, и составляет 15, в то время как для нерегулярного режима на Рис. 66 число пульсаций больше 40. Амплитуда низкочастотных пульсаций, которые соответствуют моментам взрывов, находится в диапазоне от Ap/pvN = 0.1 до 0.9 для нерегулярного режима, тогда как для сильно нерегулярного режима амплитуда пульсаций меняется в другом, более широком, диапазоне - от 0.3 и вплоть до 4. Таким образом, картина в целом для E = 35 существенно отличается от динамики при E = 28, и соответствует сильно неустойчивому режиму детонации.

Pfic

4.5 4 3.5 3 2.5

PvN 2

1.5 1

0.5

1 т 1 s S

1_~ L j

Ц \ чЛ Ч

\ \ к к J и и у

1 .....

100

200 300

Г

400

500

Рис. 69. Динамика изменения давления за фронтом ДВ в расчетной области в постановке СКФ. Сильно нерегулярный режим детонации Е = 35. N = 160 000, Ь = 600 , N,2 = 266.

Рис. 70. Полученные по методике повышенного порядка аппроксимации пространственные профили давления (а) и скорости реакции с (с) в различные моменты времени для модельной смеси при сильно неустойчивом режиме распространения ДВ. Е = 35, N = 160 000, Ь = 600 , Nl/2 = 266. Представленные в [12] аналогичные профили давления (Ь) и скорости реакции (ё) в различные моменты времени для модельной смеси с такими же физико-химическими параметрами. N = 16 000, Ь = 60 . Постановка СКФ.

Значительный интерес представляет также начальный промежуток расчета до времени т = 200 единиц. Пространственные распределения давления и скорости химической реакции с, полученные при использовании описанной методики, в некоторые моменты времени представлены на Рис. 70 (а) и (с). В

начале расчета - моменты т = 0 и 8 единиц времени, давление за фронтом ЛВ существенно падает и возникают разумные основания полагать, что ДВ не будет инициирована. Однако, вместо этого, за фронтом ЛВ возникает постепенно растущая область с газом, где происходит накопление несгоревшей модельной смеси - моменты т = 75 и 90 единиц времени. Другими словами, происходит рост зоны индукции или зоны задержки воспламенения. К моменту времени т = 90 единиц зона индукции достигает размера А^ = 35, и происходит воспламенение смеси с образованием внутренней ДВ.

1,6

1.4

1.2

Р

PvN

0.0

0.6

0.4

0.2

1 ТТТТ] |ТТТТ| |-птт| ггггг тп-г гтттт^ 1 т

-,,,, 1111 . 1,, МП 11111 ■ ■ ■ ■ INI

-500 -450 -400 -350 -300 -250 -200 -150 -100 -50 О

Рис. 71. Пространственное распределение давления за фронтом ДВ в момент времени т = 550 для постановки СКФ. Сильно неустойчивый режим детонации E = 35. N = 160 000, L = 600, N1/2 = 266.

Для демонстрации характера протекания химической реакции на Рис. 70 (b) и (d) приведены также профили скорости реакции. На этапе интенсивного протекания реакции давление и скорость за внутренней волной существенно

увеличиваются, и внутренняя волна приближается в ЛВ, уменьшая при этом зону индукции - моменты т = 98 и 100 единиц времени. После момента времени т = 104 единиц, внутренняя ДВ догоняет фронт ЛВ. Происходит соударение с ЛВ, которое проводит к инициированию самоподдерживающейся ДВ. Образовавшаяся в ходе соударения ДВ распространяется при наличии сильных пульсаций газодинамических величин за фронтом ДВ. Описанная схема периодически повторяется, что позволяет поддерживать режим ДВ, и не допускать ее затухания.

Рис. 71 иллюстрирует пространственное распределение давления за фронтом сильно неустойчивой ДВ.

Глава 4. Сравнительный анализ результатов двух постановок.

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

В данной главе рассмотрим и сравним полученные результаты для модельной смеси со следующими физико-химическими параметрами - показатель адиабаты у = 1.2, тепловой эффект Q = 50. Также как в разделах 2.6.3 и 3.6, будем рассматривать смеси при различных энергиях активации. При сопоставлении будем использовать расчетные сетки с одинаковой разрешающей способностью.

Зависимость давления за фронтом ДВ от времени можно рассматривать, как сигнал для подробного исследования с использованием математического аппарата. Так, динамику изменения давления можно исследовать количественно разложением в Фурье-спектр. А именно, можно использовать тригонометрический ряд Фурье, описанный подробно, например, в [71] и широко применяемый при исследованиях сигналов различной природы.

Пусть исследуемая непериодическая функция /) определена на [а, Ь].

Продолжим данную функцию периодически до Е ) - периодической функции

Е (?) = / ), а < ? = ? - к -А? < Ь, к . В качестве периода А ?, выбираем число,

равное длине исходного промежутка А? = Ь - а. Тогда Е) = /(?) V е[а, Ь].

Полученную таким образом периодическую функцию можно разложить в тригонометрический ряд Фурье [72]:

А п --

Е()= о +ХАкс°8(ы+Рк) Ы=2п^к, Ак=\1 ак+Ьк2, (рк=агСё(ак1Ьк), 2 к=1

a, = — í f (t)• cos kn dt, bk = — í f (t)• sin-——— dt. k At ^ At/2 k At J At/2

a a

Здесь Ak - амплитуда k-той гармоники, ( - фазовый угол k-той гармоники, cok - угловая частота k-той гармоники, vk - частота k-той гармоники. На основании полученных значений амплитуд и фазовых углов можно построить, так называемые, амплитудный и частотный спектры сигнала, где по оси ординат откладывают Ak и ( соответственно, а по оси абсцисс - частоты vk либо угловые частоты c k , в зависимости от того, что интересует исследователя.

В данной работе рассматриваются амплитудные спектры давлений, главным образом, при частотах 0 <vk < 1, поскольку более высокие значения частот отвечают высокочастотным пульсациям самоподдерживающейся ДВ, а основной задачей спектрального анализа в данной работе является выявление наличия или отсутствия периодичности (регулярности) низкочастотных пульсаций при распространении ДВ. В качестве функции f (t) рассматривается зависимость давления от времени, отнесенного к давлению в пике ФН. Отрезком [a, b] служит промежуток времени [1000,5000], на котором динамика изменения давления в обеих постановках уже вышла на определенный режим. Устойчивая детонация

Рассмотрим модельную смесь с энергией активации E = 25. Результаты расчета для канала длины L = 500 с числом ячеек N = 100 000 в постановке в ЛСК, а также канала длины L = 20 с числом ячеек N = 4000 постановки СКФ представлены на Рис. 72. В обеих постановках давление выходит на уровень теоретического значения давления в пике ФН. При этом в постановке ЛСК давление выходит на уровень после момента времени t = 200 , тогда как в постановке СКФ выход на уровень происходит после т = 500 единиц времени, что значительно позже. Помимо различий в моменте установления ДВ, существенные различия связаны с характером установления. В постановке ЛСК происхо-

дит мгновенное сгорание смеси в узкой области, с образованием сильно пересжатой ДВ, которая ослабляется с выходом газодинамических величин на стационарные значения. Динамика изменения пикового давления при этом демонстрирует полное отсутствие низкочастотных пульсаций, и наличие из-за особенностей постановки ЛСК высокочастотных пульсаций. В постановке СКФ ДВ инициирована еще до начала расчетов, в качестве начального распределения величин. Однако, распределение соответствует решению ЗНД идеальной стационарной ДВ, которое может отличаться от установившегося решения нестационарной ДВ.

1.5 1.45 1.4 1.35 1.3

1.25

Рреак | 2

АлГ 1.15

1.1 1.05 1

0.95 0.9

О 100

. I.... I .

200 300 400 500 600 700 300 900 1000

X

1.01 1.0075 1.005 1.0025

Р ^"ОШ

1

РтК

0.9975 0.995 0.9925 0.99 0.9875 0.985

1 и 1 1 I 1 1 | 1 1 1 1 ■ ■ ■ ■'■II

: 1, ;

:

11

1

;

: ;

:|, , , ми III! II ......... 1 1 1 1 ММ II м

0 100 200 300 400 500 600 700 800 900 1000

Г

Рис. 72. Истории изменения давления в разных постановках для модельной смеси с энергией активации Е = 25, размер ячейки расчетной сетки 5 -10-3. Слева - история пикового давления в постановке в ЛСК при длине канала Ь = 500. Справа - история давления за фронтом ДВ в СКФ при длине канала Ь = 20.

Процесс установления представляет собой переход от решения ЗНД к решению нестационарной системы уравнений математической модели. Как видно из Рис. 72, процесс перехода характеризуется затухающими низкочастотными пульсациями. Высокочастотные пульсации при этом отсутствуют при установлении

ДВ, и не появляются после ее установления. Таким образом, затухающие низ-

133

кочастотные пульсации в СКФ и их отсутствие в ЛСК свидетельствуют о том, что эти низкочастотные колебания являются лишь численным эффектом, связанным с постановкой СКФ.

Сравнение результатов, представленных на Рис. 72, происходит больше на качественном уровне, чем на количественном. Более детальное сравнение может быть проведено при использовании Фурье-анализа.

О I I I I I | М I I М II I И I I I Г I I I II | I I I I | М I I I I I I I I I I I I

10 г —I ' ' ч

лск

" - СКФ

О 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Рис. 73. Фурье спектр давления, отнесенного к давлению в пике ФН, для промежутка времени t,ге[1000,5000] модельной смеси с энергией активации

Е = 25, размер ячейки расчетной сетки Ах = А^ = 5 • 10-3. Красный цвет соответствует постановке в ЛСК. Зеленым цветом обозначен результат для постановки в СКФ.

Фурье-спектры давления, представленные на Рис. 73, дают возможность увидеть новые результаты для постановок. Так, Фурье анализ подтверждает, что в постановке СКФ произошло установление ДВ. Имеется всего одна гармо-

ника - основная гармоника с амплитудой А0 = 1.0002. Для постановки ЛСК

также имеется основная гармоника. Ее амплитуда составляет А0 = 0.9984, что

более заметно отличается от значения 1.0, чем значение основной гармоники СКФ. На совокупность множества гармоник малой амплитуды - порядка

Ак < 10-7, составляющих шум, приходится остальная часть амплитуды давления. Полученный шум сравнительно малой амплитуды, отвечает всем частотам, как малым, так и большим, и являются особенностью постановки в ЛСК.

Слабо неустойчивая детонация

Рассмотрим модельную смесь с энергией активации Е = 26. Результаты расчета для канала длиной Ь = 1500 и общим числом ячеек N = 300 000 в постановке в ЛСК, а также канала длиной Ь = 600 с числом ячеек N = 120 000 для постановки в СКФ представлены на Рис. 74. В обеих постановках давление выходит на предельный цикл, когда низкочастотные пульсации давления оказываются регулярными во времени. При этом в постановке в ЛСК давление выходит на цикл после момента времени t = 280 единиц, тогда как в постановке СКФ выход на уровень происходит после т = 110 единиц времени, что существенно раньше. Помимо различий в моменте выхода пульсаций на предельный цикл, существенные различия связаны с характером установления. В постановке ЛСК выход на цикл происходит после ослабления пересжатой волны, возникающей вследствие мгновенного сгорания реагирующей смеси. Динамика изменения пикового давления при этом демонстрирует наличие низкочастотных пульсаций, начиная с этапа ослабления пересжатой волны. В постановке СКФ величины за фронтом ДВ начинают выходить на предельный цикл уже на начальном этапе расчетов. Как видно из Рис. 74, процесс перехода к предельному циклу в обеих постановках характеризуется низкочастотными пульсациями растущей амплитуды. Таким образом, наличие регулярных низкочастотных пульсаций в обеих постановках дает основания считать, что эти пульсации являются не численным, а физическим эффектом. Значения периода пульсаций

циклов в обеих постановках оказываются близкими. Так, для постановки ЛСК значение периода составляет Аt = 11.853, а для СКФ - Ат = 11.814. Что касается высокочастотных пульсаций, то они сопровождают постановку ЛСК в течение всего времени расчета, тогда как в постановке СКФ он полностью отсутствуют. Поэтому разумно полагать, что высокочастотные пульсации обусловлены особенностями постановки ЛСК и численного метода.

Рис. 74. Истории изменения давления в разных постановках для модельной смеси с энергией активации Е = 26, размер ячейки расчетной сетки Ах = = 5 • 10-3. Слева - история пикового давления в постановке в ЛСК при длине канала Ь = 1500. Справа - история давления за фронтом ДВ при длине канала Ь = 600 в постановке в СКФ.

Проведем более детальное сравнение результатов с использованием Фурье-анализа. Фурье-спектры давления, представленные на Рис. 75, позволяют увидеть новые результаты для исследуемых постановок. Так, Фурье анализ показывает наличие для обеих постановок выделенных гармоник кратной частоты, которые в сумме образуют регулярные пульсации, схожие по форме с сину-соидными колебаниям. При этом шум начинает играть роль по сокрытию выделенных гармоник при частотах ук > 0.8, причем в обеих постановках. Таким

образом, с достаточно хорошей точностью можно представить зависимость давления от времени, как сумму сравнительно небольшого числа слагаемых (п ^ 10) ряда Фурье.

Рис. 75. Фурье спектр давления, отнесенного к давлению в пике ФН, для промежутка времени t,ге[1000,5000] модельной смеси с энергией активации

Е = 26, размере ячейки расчетной сетки Ах = А£ = 5 • 10-3. Красный цвет соответствует постановке в ЛСК. Зеленым цветом обозначен результат для постановки в СКФ.

Сравнивая конкретные значения частот в двух постановках, можно отметить, что частоты оказываются практически одинаковыми, несмотря на различия в начальных распределениях и механизмов выхода на предельный цикл. Более точные значения частот представлены в Таб. 11. Так, для частот в промежутке 0 <ук < 1 разница оказываются в 3-4 знаках после запятой. Наличие

таких отчетливых гармоник в длительных расчетах подтверждает, что реализу-

137

ется слабо неустойчивый режим с предельным циклом изменения параметров за фронтом ДВ. Частоты гармоник совпадают в обеих постановках, что позволяет сделать вывод о физической природе данных гармонических колебаний.

к УГ уГ

0 0 0

1 0.08442 0.08464

2 0.16905 0.16929

3 0.25347 0.2537

4 0.33789 0.33834

5 0.42252 0.42299

6 0.50694 0.50764

7 0.59136 0.59228

8 0.67578 0.67693

9 0.76041 0.76133

10 0.92947 0.93064

Таб. 11. Частоты выделенных гармоник Фурье-спектра давления, отнесенного к давлению в пике ФН для промежутка времени [1000,5000] модельной смеси с энергией активации Е = 26 в обеих постановках. Нерегулярная детонация

Рассмотрим модельную смесь с энергией активации Е = 28. Результаты расчета для канала длиной £ = 750 и общем числом ячеек N = 300 000 для постановки в ЛСК, а также канала длиной £ = 600 с числом ячеек N = 240 000 для постановки в СКФ представлены на Рис. 76. В обеих постановках давление не выходит какой-либо регулярный режим, а напротив, испытывает нерегуляр-

ные пульсации во времени. В обеих постановках промежутки между соседними пульсациями составляют значения из отрезка t ,те[10.5,13.5] единиц времени, а давление пульсаций, приведенное к давлению в пике ФН, меняется в диапазоне р/руМ е[0.75,1.68]. Отсутствуют повторяющиеся участки пульсаций, когда соседние пульсации отличаются друг от друга, но пучок, состоящий из нескольких последовательных пульсаций, может повторяться через определенные промежутки времени.

Рис. 76. Истории изменения давления в разных постановках для модельной смеси с энергией активации Е = 28, размер ячейки расчетной сетки Ах = А^ = 2.5 -10-3. Слева - история пикового давления в постановке ЛСК при длине канала Ь = 750 . Справа - история давления за фронтом ДВ при длине канала Ь = 600.

В постановке в ЛСК давление выходит на нерегулярный режим при ослаблении пересжатой ДВ после t = 30 единиц времени, тогда как в постановке СКФ выход на режим происходит с начала расчетов и, соответственно, с момента времени т = 0 . Характер первых пульсаций одинаковый в обеих постановках. Сначала происходит 3 пульсации, пиковые значения давления в которых монотонно увеличиваются, а затем для четвертой пульсации пиковое значение стано-

вится меньше, чем у предыдущей. Начиная с пятой пульсации, не только количественно, но и качественно становится сложно сопоставлять картины разных постановок для выделения общих закономерностей и механизмов режима. Так, на участке ^ те [300,400] в обеих постановках происходит несколько пульсаций, схожих между собой в плане размера, пиковых, минимальных значений. Причем, характеристики данных пульсаций оказываются почти одинаковыми для двух различных постановок. Также стоит отметить, что пульсации в постановке в ЛСК на участке t е [300,500] качественно схожи с пульсациями в постановке в СКФ на участке t е [320,530]. Участки отмечены пунктирными линиями на Рис. 76. В обеих постановках участки начинаются с пульсации, пиковое значение приведенного давления которой составляет примерно Р/PvN = 1.66 . Затем следует пульсация с низким значением пикового давления - р/Рш = 1125. После чего происходит несколько вышеупомянутых пульсаций схожих между собой размеров. Затем идут две пульсации небольшой величины пикового давления - примерно р/рш = 1.3 и 1.125, между которыми пульсация высокого значения - около р/рш = 1.665. Затем идет несколько пульсаций, пиковые значения давления которых поочередно увеличиваются и уменьшаются.

В целом, картины давления для двух постановок во многом схожи между собой. Отметим, что как и для результатов предыдущих режимов детонации, постановка в ЛСК характеризуется наличием не только низкочастотных пульсаций, но и высокочастотных. У постановки в СКФ отсутствуют высокочастотные пульсации. Таким образом, наличие нерегулярных низкочастотных пульсаций в обеих постановках дает основания полагать, что эти пульсации являются физическим эффектом.

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

высоких значений. Так, у частот vk е[0,1] диапазон амплитуд приведенного давления для гармоник шума оказывается достаточно широким Ak е [10-6,3 -10-2]. Предельный цикл у такого режима детонации, несмотря на 5 выделенных гармоник, отсутствует, о чем свидетельствует наличие сильного шума.

Что касается выделенных гармоник, то, сравнивая значения частот в двух постановках, можно отметить, что частоты оказываются практически одинаковыми, несмотря на различия в начальных распределениях и механизмов выхода на предельный цикл. Более точные значения частот представлены в Таб. 12. Так, разница в значениях частот оказывается в 4 знаке после запятой.

ю1 10° ю1

Ю"2

А ю-3

Ю"4 105

10 6 н——————————1

-I о7 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Ъ

Рис. 77. Фурье спектр давления, отнесенного к давлению в пике ФН, для промежутка времени t,t е [1000,5000] модельной смеси с энергией активации

E = 28, размер ячейки расчетной сетки Ах = А^ = 2.5 -10-3. Красный цвет соответствует постановке в ЛСК. Зеленым цветом обозначен результат для постановки в СКФ.

к УГ уГ

0 0 0

1 0.08365 0.08371

2 0.16709 0.16721

3 0.25074 0.25093

4 0.33418 0.33464

Таб. 12. Частоты выделенных гармоник Фурье-спектра давления, отнесенного к давлению в пике ФН для промежутка времени ¿,те [1000,5000] модельной смеси с энергией активации Е = 28 в обеих постановках.

Сильно неустойчивая детонация

Рассмотрим модельную смесь с энергией активации Е = 35. Результаты расчета для канала длиной £ = 3600 и общим числом ячеек N = 960 000 в постановке в ЛСК, а также канала длиной £ = 600 с числом ячеек N = 160 000 в постановке в СКФ представлены на Рис. 78 в виде зависимостей давления от времени. В обеих постановках давление не выходит какой-либо регулярный режим, а подвергается нерегулярным пульсациям во времени. В отличие от нерегулярного режима промежутки между соседними пульсациями меняются в гораздо большем диапазоне. Так, при одинаковых по времени промежутках, на одном промежутке может произойти 1 пульсация, а на другом 4 пульсации. Однако, здесь стоит отметить, что по вопросу неравномерности пульсаций в постановке СКФ реализуется существенно большая неравномерность. Картина пульсации давления в СКФ включает в себя как участки очень частых пульсаций, так и участки очень редких. Кроме того, частота пульсаций в СКФ также выше, чем в ЛСК. На промежутке t ,те [0,1000] единиц времени число пульсаций в СКФ выше, чем в ЛСК более, чем на 30%, а на промежутках, соответствующих дальнейшим моментам времени, это число выше. Помимо особенно-

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

Рис. 78. Истории изменения давления в разных постановках для модельной смеси с энергией активации Е = 35, размер ячейки расчетной сетки Ах = А^ = 6.25 • 10-4. Слева - история пикового давления в постановке ЛСК при длине канала £ = 3600. Справа - история давления за фронтом ДВ при длине канала £ = 600 .

В постановке в ЛСК давление выходит на сильно нерегулярный режим при ослаблении пересжатой ДВ после ? = 60 единиц времени, тогда как в постановке в СКФ резкий выход на режим происходит после т = 100 единиц времени. Реи-

нициирование происходит за счет образования внутренней ДВ, которая догоняет ослабевающую ЛВ и соударяется с ней, инициируя ДВ (см. раздел 3.6).

После столкновения волн, которое характеризуется пульсацией давления большого размера - максимальное значение пульсации давления, приведенного к давлению в пике ФН, составляет почти р/= 4.5 - давление пересжатой ДВ на протяжении т = 100 единиц времени падает. Затем начинается участок частых пульсаций, после которого наступает этап, когда частые пульсации сменяются редкими и наоборот. Подобная динамика сильно неустойчивого режима представлена также, например, в [12].

кг5 |------------

10® _ц ' ' ' 1 ' 1 ' ' ' • ' ' ' ' _I I I I I I I I I I I I I_I I I I_м I I I I

О 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Рис. 79. Фурье спектр давления, отнесенного к давлению в пике ФН, для промежутка времени t,те[1000,5000] для модельной смеси с энергией активации

Е = 35, размер ячейки расчетной сетки Ах = А^ = 6.25 • 10-4. Красный цвет соответствует постановке в ЛСК. Зеленым цветом обозначен результат для постановки в СКФ.

Далее, рассмотрим результаты проведенного Фурье-анализа для приведенного давления обеих постановок. Фурье-спектры давления, представленные на Рис. 79, демонстрируют отсутствие в спектрах обеих постановок каких-либо выделенных частот, кроме основной, что свидетельствует о полной нерегулярности пульсаций. Другими словами, зависимость давления от времени нельзя представить с хорошей степенью точности в виде суммы сравнительно небольшого числа слагаемых тригонометрического ряда Фурье.

Вместо нескольких выделенных частот, Фурье спектры для обеих постановок характеризуются наличием множества гармоник достаточно высокой амплитуды. При этом, начиная с частоты ук = 0.2, амплитуды гармоник в постановке в СКФ в среднем оказываются выше, чем в ЛСК. В результате, какой-либо регулярный характер изменения давления у такого режима детонации в обеих постановках отсутствует, о чем свидетельствует наличие множества гармоник схожей амплитуды при слабо меняющейся частоте.

Таким образом, сравнение результатов используемых в работе постановок в ЛСК и СКФ, которое проведено в данной главе, позволяет сделать ряд заключений. Во-первых, несмотря на разницу в механизмах инициирования ДВ для разных постановок и в картинах изменения давления в начальные промежутки времени, начиная с некоторого момента времени, порядка 500 единиц, картины пульсаций давления становятся в целом похожими в разных постановках для всех рассмотренных режимов детонации. Более того, для устойчивого и слабо неустойчивого режимов значения газодинамических величин на участках после 500 единиц времени и вовсе практически совпадают. Данный факт дает основания полагать, что полученные результаты отражают действительный характер точного решения используемых математических моделей.

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

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

В-третьих, хотя в данной главе, как и в предыдущих, было мало слов о разнице в быстродействии вычислительных методик используемых постановок ЛСК и СКФ, нельзя не отметить, что постановка в СКФ дает возможность получить результаты, совпадающие с результатами в постановке в ЛСК, как минимум, на порядок быстрее, что достигается за счет использования существенно более короткой расчетной области. Постановка в ЛСК требует для расчетов наличия канала длины, достаточной, по крайней мере, для формирования ДВ путем перехода ее из пересжатого состояния, вызванного резким сгоранием смеси в узкой области у левого конца канала, в самоподдерживающееся состояние.

Глава 5. Анализ эффективности распараллеливания вычислительного алгоритма.

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

В рамках распараллеливания вычислительных алгоритмов, использованных в данной работе, см. разделы 2.3 и 3.3, была проведена процедура распараллеливания методом декомпозиции расчетной области. Расчетная область с числом ячеек N разбивается поровну между процессорами вдоль координатного направления х . Другими словами, разделим расчетную область на р подобластей, в каждой из которых будет N = NP ячеек. В результате, каждый процесс k отвечает за эволюцию газодинамических величин - плотности, скорости, давления, доли реагента - в подобласти с Np ячейками.

После инициирования процедуры распараллеливания в программном коде один из процессов - «нулевой» - считывает информацию из файла parame-ters.dat с описанием геометрии, постановки задачи, а также численным методом ее решения. В соответствии с информацией входного файла parameters.dat, нулевой процесс создает файл initial_distribution.dat с пространственным распределением определяющих газодинамических величин. При использовании постановки в ЛСК пространственные распределения величин могут быть построены сразу, поскольку представляют собой профили ступенчатого вида. Применение же постановки в СКФ требует построения профилей ЗНД.

Заполненный необходимой информацией файл initial_distribution.dat открывается каждым из процессов. Процесс k считывает из файла нужный набор из Nстрок, преобразует значения этих неконсервативных величин в консервативные, и заполняет этими значениями двумерный массив размерности N х 4.

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

Межпроцессорный обмен обусловлен следующими факторами. Во-первых, шаг интегрирования по времени рассчитывается в данной работе динамически. Поэтому процесс к должен сначала рассчитать минимальный шаг At" для своей подобласти из N ячеек, а затем послать значение At" нулевому

процессу. Нулевой процесс должен получить к значений шагов - At", к = 0,1,...,p - чтобы затем выбрать минимальное: At" = minк(Аt"). После вычисления At" нулевой процесс рассылает остальным вычисленное значение шага по времени. Таким образом, обмен данными между процессами связан, в данном случае, с расчетом шага по времени.

Во-вторых, при использовании ENO-реконструкции на гранях ячеек нужно получить значения консервативных величин близлежащих ячеек. Для процесса к в получении значений нет сложностей, если рассматриваются грани ячеек внутри подобласти. В этом случае используются элементы заполненного массива размерности Np х 4. Однако, для вычисления по ENO-реконструкции

значений на гранях ячеек, которые располагаются вблизи границы подобластей, требуется получить значения консервативных величин в ячейках соседних процессов (к -1) и (к +1). Другими словами, необходим обмен между соседними процессами значениями в приграничных ячейках для осуществления ENO-реконструкции. В данной работе преимущественно используется ENO-реконструкция второго порядка. Поэтому процессу к необходимо получить 2 х 4 значений консервативных величин от процесса к -1 и 2 х 4 значений от процесса к +1. Также нужно послать 2 х 4 значений двух крайних левых ячеек процессу к -1 и 2 х 4 значений двух крайних правых ячеек процессу к +1. Такой обмен нужно осуществлять не только на каждом временном шаге, но и на каждой стадии многостадийного метода интегрирования по времени. В данной

работе используется преимущественно метод Рунге-Кутты второго порядка аппроксимации, поэтому при переходе с временного слоя п на слой п +1 процессу k нужно получить 2 х 2 х 4 значений от процесса k -1 и столько же значений от процесса k +1.

В-третьих, при использовании постановки в СКФ на каждом шаге по времени требуется вычислять значение скорости ДВ Dn+1. Поскольку численный метод вычисления Dn+1 основан на использовании ячеек вблизи правой границы всей расчетной области, то разумно использовать на данном этапе процесс к = p, который содержит значения этих ячеек.

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

При реализации вышеописанного распараллеленного алгоритма в коде, его отладке и проведении тестовых расчетов для задач настоящей работы, а также для проведения основных расчетов данной работы, результаты которых представлены в разделах 2.6 и 3.6, использовался суперкомпьютер MBC-100k Межведомственного суперкомпьютерного центра РАН.

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

Sp = ГГ, (3.14)

p

т.е. как отношение времени решения задачи на скалярной ЭВМ T1 к времени Tp выполнения параллельного алгоритма на p процессорных ядрах. Другой важной характеристикой является эффективность использования параллельного алгоритма при решении задачи. Эффективность вычисляется как:

EP = ^ = (3-15)

P PTp

и определяет среднюю долю времени выполнения параллельного алгоритма, в течение которого процессоры реально используются для решения задачи.

При реализации параллельного алгоритма в коде используется интерфейс передачи сообщений - MPI (Message Passing Interface), который является наиболее распространенной технологией программирования параллельных компьютеров с распределенной памятью в настоящее время. MPI-программа представляет собой множество параллельных взаимодействующих между собой процессов. В соответствии с концепцией MPI основным способом обмена информацией между процессами является посылка сообщений - наборов данных некоторого типа.

Для получения количественных оценок основных характеристик параллельных версий вычислительных алгоритмов - ускорения и эффективности - на MBC-100k, рассмотрим тестовую задачу об устойчивом распространении ДВ в канале длины L = 403.2, числом ячеек N = 80640. Расчеты данного теста проводились для p процессорных ядер, где p = 1, 8, 16, 30, 48, 96 и 192. На Рис. 80

и Рис. 81 приведены результаты в виде графиков зависимости ускорения и эффективности от числа процессорных ядер для обеих постановок - ЛСК и СКФ. Анализ приведенных графиков показывает, что в обеих постановках при использовании небольшого числа ядер (до 16 шт.) ускорение близко к линейному (идеальному), а эффективность такого распараллеливания находится на высоком уровне (выше 90%). Однако, при введении большего числа процессорных ядер в обеих постановках ускорение отклоняется от линейного, и эффектив-

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

Неформальным критерием целесообразности использования того или иного числа процессорных ядер является величина эффективности 50%. Учитывая данный критерий, можно сказать, что описанные выше алгоритмы оказывается целесообразным использовать при числе процессорных ядер до р = 300 шт.

Рис. 80. Зависимость ускорения от числа процессорных ядер при расчете 10 единиц времени тестовой задачи (~ 80 000 ячеек) в двух постановках, исследуемых в работе, для суперкомпьютера МВС-100к, используемого для проведения основных вычислительных экспериментов.

Число процессорных ядер

Рис. 81. Зависимость эффективности от числа процессорных ядер при расчете 10 единиц времени тестовой задачи (~ 80 000 ячеек) в двух постановках, исследуемых в работе, для суперкомпьютера МВС-100к, используемого для проведения основных вычислительных экспериментов.

ЗАКЛЮЧЕНИЕ

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

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

3. Проведено численное исследование распространения пульсирующей волны детонации для четырех режимов - устойчивого, слабо неустойчивого, нерегулярного и сильно неустойчивого в лабораторной системе координат и в системе координат, связанной с лидирующим скачком. Для каждого из режимов для обеих постановок: (1) представлены типичные пространственные профили давления; (2) описана динамика пульсаций пикового давления или давления за фронтом; (3) сделано заключение о характере пульсаций для расчета длительного распространения на момент времени 5000 единиц; (4) рассчитан Фурье-спектр пульсаций давления. Результат расчетов устойчивого режима демонстрируют выход детонационной волны на уровень Чепмена-Жуге в обеих постановках. В результате расчета слабо неустойчивого режима получен выход на стационарный предельный цикл, причем обе постановки согласуются между собой и с известным аналитическим решением. При расчете нерегулярного режима методом сквозного счета отмечена тенденция к занижению пикового давления при пульсациях, что может являться причиной срыва детонации в одномерных расчетах сильно неустойчивых режимов.

СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ

Латинские символы

А - предэкспоненциальный множитель; амплитуда гармоники Фурье спектра в

разделе Глава 4

с - скорость звука

D - скорость лидирующей волны

Е - энергия активации

е - полная энергия единицы объема газа

/ - степень пересжатия детонационной волны

Н - функция Хевисайда

к - удельная энтальпия

к - номер «внутреннего» шага по времени в разделе 2.3, номер гармоники ряда Фурье в Главе 4, номер процессора в Главе 5

К - константа скорости реакции в разделе 1.3, число шагов по времени при интегрировании уравнений кинетики в разделе 2.3

Ь - длина расчетной области в разделах 2.1, 3.1, дифференциальный оператор в разделах 2.3, 3.3

1 - длина области инициирования

т - масса; индекс для формулы численного интегрирования в разделе 3.4 N - число ячеек расчетной области р - давление в Главах 1 - 4, число процессоров в Главе 5 Q - тепловой эффект химической реакции

д - сеточная функция в разделе 2.4, порядок аппроксимации численного метода в разделах 2.5, 3.6 Я - универсальная газовая постоянная Т - температура V - скорость

2 - массовая доля реагента

Греческие символы

Y - показатель адиабаты

G - внутренняя энергия единицы объема газа X - доля периода индукции или реакции ß - молярная масса

V - порядок реакции в разделе 1.3; частота гармоники Фурье спектра в разделе Глава 4

р - плотность и - обратная плотность

ç - величина критерия достижения заданной точности метода с - скорость протекания химических реакций

Векторы

f - вектор потока в дифференциальной задаче постановки, связанной с лабораторной системой координат

F - вектор численного потока в разностной задаче постановки, связанной с лабораторной системой координат s - вектор источниковых членов

u - вектор консервативных переменных в дифференциальной задаче U - вектор консервативных переменных в разностной задаче v - вектор примитивных переменных в дифференциальной задаче

Нижние индексы

0 - начальные параметры перед фронтом волны или аналитическое решение Зельдовича-Неймана-Деринга

B - пик фон Неймана

CJ - параметры Чепмена-Жуге

H - точка Чепмена-Жуге

1 - индукция в разделе 1.3; номер расчетной ячейки в разделе 2.3 r - реакция

s - характерные величины для приведения к безразмерному виду vN - параметры фон Неймана

Верхние индексы