Численное моделирование трехмерных течений с волнами детонации на многопроцессорной вычислительной технике

УПРАВЛЕНИЕ, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И ИНФОРМАТИКА

УДК 534.222

И. В. Семенов, И. Ф. Ахмедьянов, П. С. Уткин, А. Ю. Лебедева

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ ТЕЧЕНИЙ С ВОЛНАМИ ДЕТОНАЦИИ НА МНОГОПРОЦЕССОРНОЙ ВЫЧИСЛИТЕЛЬНОЙ ТЕХНИКЕ

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

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

Как известно, реальная детонация в газах сопровождается образованием сложной неодномерной и нестационарной структуры течения за ее передним фронтом [1]. В частности, с 1926 года известно явление спиновой детонации, суть которого заключается в том, что в трубах вблизи детонационных пределов наиболее яркое свечение фронта волны сосредотачивается у стенки - в «голове» спина, вращающейся по окружности одновременно с поступательным движением фронта [2]. Спиновый режим впервые наблюдался С. Кемпбеллом, Д. В. Вудхедом и А. С. Финчем как периодические неоднородности на фоторазвертках самосвечения детонации в смесях окиси углерода с кислородом. Другим примером являются экспериментально обнаруженные в 1957-1958 гг. ячеистые струк-

Контактная информация: semenov@icad.org.ru Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии-2010». Работа выполнена в рамках Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» (конкурс НК-100П, контракт № П-359).

туры детонационного фронта вдали от пределов, обусловленные наличием сильных неоднородностей - поперечных волн в зоне за лидирующим скачком [2].

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

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

Исследование проводится в рамках идеологии вычислительного эксперимента [3] на многопроцессорных ЭВМ с распределенной памятью с использованием разработанного авторами программного комплекса для численного исследования динамики потоков реагирующих многокомпонентных сред в многомерных областях. В [4] изложены модели, численный метод и ал-

горитм распараллеливания для решения двумерных осесимметричных задач, связанных с инициированием и распространением волн газовой детонации. Данная работа служит дальнейшим развитием и расширением исследования [4] для трехмерного случая.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Исследование проводится с использованием системы уравнений, описывающей трехмерные нестационарные течения невязкой сжимаемой многокомпонентной реагирующей газовой смеси, которая в декартовой системе координат (х, у, z) имеет следующий вид:

Эц ЭГ ЭГ ЭГ,

— + — + — + — = 8,

Эt Эх Эу Эz

N

e = I ^e, 1=1 P

здесь ґ — время; их, иу и иг — составляющие скорости; р, р и Е — плотность, давление и полная удельная энергия газовой смеси соответственно; рь — плотность и скорость изменения плот-

ности /-го компонента смеси в результате химических реакций; д — тепловой эффект химических реакций. Полная удельная энергия газовой смеси определяется как:

U2 + U2 + U2

E = x y z 2

+ e,

где е — внутренняя удельная энергия газовой смеси.

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

N

1 = 1

e, = cv,T, i = I- N,

где Т - температура; Я — универсальная газовая постоянная; е, еУ1 и ц — удельная внутренняя энергия, удельная теплоемкость при постоянном объеме и молярная масса /-го компонента смеси; N — число компонентов смеси.

Химические реакции моделируются одностадийной кинетикой горения пропана [5]:

C3H8 + 5O2 + 5 • 3.762N2

3CO2 + 4H2O + 5 • 3.762N2 + Q,

где Q = 46,6-10 Дж/кг - теплота сгорания, рассчитанная на единицу массы топлива.

Таким образом, число компонентов N газовой смеси равно 5. Им приписываются индексы в соответствии со схемой: С3Н8 (/ = 1), 02 (/ = 2), N2 (/ = 3), СО2 (/ = 4), Н2О (/ = 5).

Скорость изменения плотности пропана ^ определяется как:

( „ „ А

w =т1<Щ‘ =m

kr •Pirn m2

k = —7 • 108 p-°-2264 exp (— E*/RT) E* = 190.3 1 03 Дж/моль,

моль • с

mole

где ^>1 - скорость изменения его молярной

концентрации, р - давление в атмосферах, Т -температура в градусах Кельвина. Скорости изменения плотностей остальных компонентов смеси определяются через м™010 и стехиометрические коэффициенты в реакции:

^2 = 5ц2^5ПО10<“2 = 0<

ш4 = -Зц4йГ10,й5 = -4ц5<010.

Тепловой эффект реакции определяется как:

& = -0“1-

2. ЧИСЛЕННЫЙ МЕТОД

Для решения задачи (1)-(3), (5)-(7) используется метод расщепления по физическим процессам. На шаге по времени п в текущей расчетной ячейке / сначала решается система уравнений газовой динамики (1) без учета члена 8 в правой части. Затем полученное решение

газ п+1

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

Для дискретизации системы уравнений газовой динамики по пространственным переменным используется метод конечных объемов (МКО). Основными преимуществами класса подобных схем являются их вычислительная

3

м

робастность, консервативность, а также возможность работы на неструктурированных сетках [6]. Используется так называемая центрированная схема МКО, когда конечные объемы совпадают с ячейками расчетной сетки. Второй порядок точности на гладких решениях достигается путем сочетания использования кусочнолинейной аппроксимации величин внутри ячейки с двухшаговым пересчетом по времени, который также называется предиктор-корректор [7]. Потоки вычисляются на основе точного решения задачи Римана о распаде произвольного разрыва методом типа Г одунова [8]. При расчете градиента вектора-решения чП в текущей расчетной ячейке / на текущем шаге по времени п для экстраполяции значений газодинамических параметров из центра ячейки на грани применяется метод наименьших квадратов. Будем обозначать вектор экстраполированных величин на грань о как чП^. Шаг интегрирования по времени Дt определяется динамически в процессе вычислений для выполнения условия устойчивости [8].

Этап предиктора может быть записан в следующем виде:

чп+1 = ап — — V Ёп  5.

тт/ О Т/7" /,& /,& ^

2V/ а

где V/ - объем текущей ячейки, суммирование ведется по всем граням о ячейки, 5/о - площадь соответствующей грани, - поток через

грань о:

= ТХ, ^ = Г №,*), Qn,s = TX,

ii,S’

где То - матрица преобразования от лабораторной системы координат к локальной системе координат, один орт которой ориентирован по внешней нормали, а другие лежат в плоскости грани о.

Затем следует этап корректора:

Аґ ,

где

qr,n+1 = qn ——I Fn

Y .• s

Fn rfi—1 е n+1 е n+1 P ( n +1 \

i,s = Ts *i,s , *i,s = f1 (Qi,s ),

а вектор ф п++1 - решение задачи Римана о распаде произвольного разрыва, соответствующее

,в*'п + 1 ,в*'п + 1

начальным данным ч и ч ] а по разные стороны от разрыва, определяемого гранью о. Здесь ч Та обозначает вектор экстраполированных величин на грань о, рассчитанный по зна-

чениям ч, с шага предиктора и градиенту Vчn, а индекс ] соответствует параметрам в ячейке, граничащей по грани о с рассматриваемой ячейкой /.

На следующем этапе производится учет источниковых членов в правой части (1), обусловленных химическими реакциями:

;/{| ( / I

“г V / ____£1 + + + 1 (+ \ __^газ,п+1

dt

здесь tn — текущее время. В качестве начального условия для системы обыкновенных дифференциальных уравнений (8) берется решение, полученное на предыдущем газодинамическом этапе. Для решения системы ОДУ используется многошаговый неявный метод Гира. Построенное решение, взятое в момент времени tn+i, будет решением всей задачи на следующем шаге по времени:

n+1 _ хим ( j \

Ч/ _ Ч/ Vn+1 ).

3. РАСПАРАЛЛЕЛИВАНИЕ РАСЧЕТНОГО АЛГОРИТМА И СЕТКИ

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

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

сплошной среды, таких как ANSYS, для осуществления декомпозиции расчетной области. Оценки эффективности разбиения сетки с помощью METIS по сравнению с аналогами (JOSTLE [10], PARTY [11], SCOTCH [12]) оказываются намного лучше.

На рис. 1 приведен пример используемой расчетной О-сетки и ее разбиения с помощью METIS на K = 10 частей для трубы с параболическим сужением и коническим расширением, процесс инициирования детонации в которой будет рассмотрен далее. Рассмотрим количественные характеристики качества разбиения сетки с помощью METIS. Размер тестовой сетки, представленной на рис. 1, составляет 2 929 360 ячеек, общее число граней - 8 734 740. В результате разбиения отклонение числа ячеек во всех частях от равномерного 2 929 360 / K не превышает 2%. При этом общее количество граней, по которым граничат ячейки из разных частей, составляет 98 693. Оценим минимальное число граней, по которым будут граничить ячейки из различных частей, в случае одномерной декомпозиции и использования структурированной сетки. Число ячеек в поперечном сечении составляет примерно 10 000, значит, количество граничных ячеек при разбиении на 10 частей - около 90 000, так что фактическое число граничных ячеек отличается от оценочного минимального для одномерной декомпозиции менее чем на 10%.

Рис. 1. Пример расчетной О-сетки и ее разбиения с помощью METIS на 10 частей

Для расчетов используются от 200 до 600 процессорных ядер СК МВС-100к Межведомственного суперкомпьютерного центра РАН или СК «СКИФ» МГУ «Чебышев». Характерное пространственное разрешение процессов составляет 0,3 мм, временное - 5-10 нс.

На рис. 2 представлена зависимость ускорения от числа используемых процессорных ядер при решении тестовой задачи на сетке с числом ячеек порядка 9 млн. Стоит отметить, что при использовании 190 процессорных ядер ускорение составляет 94, а эффективность, соответственно, около 50%, что является приемлемым

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

80

60

0)

S

X

CD

I 40

20

50 100 150

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

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

4. СПОСОБЫ ИНИЦИИРОВАНИЯ ДЕТОНАЦИИ

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

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

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

смеси достаточно мощную УВ, способную вызвать практически мгновенное воспламенение смеси [13]. Мощная УВ, инициирующая интенсивные экзотермические реакции, возникает при сильном плоском, цилиндрическом или сферическом взрыве в реагирующей газовой смеси [14]. Такой способ инициирования детонации принято называть прямым инициированием или инициированием с помощью концентрированного подвода энергии. ДВ, образующаяся в результате взрыва, является пересжатой. Численные и экспериментальные исследования показывают, что пересжатая ДВ переходит к само-поддерживающейся детонации Чепмена-Жуге (ЧЖ) на конечном расстоянии от места взрыва, составляющем несколько десятков поперечных размеров детонационных ячеек. Это проявляется не только в уменьшении скорости ДВ, но и в развивающейся неустойчивости фронта и появлении системы поперечных волн, типичных для самоподдерживающейся детонации в газе [1]. Важным параметром является критическая энергия взрыва, ниже которой самопод-держивающийся детонационный режим не формируется [14]. Существование критической энергии связано с конечной скоростью протекания химических реакций в сжатой и нагретой ударной волной смеси. Эти реакции включают в себя процессы диссоциации с образованием активных центров, протекающие во время так называемого периода индукции без выделения тепла и даже с поглощением тепловой энергии среды, и последующие экзотермические реакции рекомбинации, при которых образуются продукты детонации. Прямое экспериментальное определение критических условий инициирования и распространения детонации могло бы ответить на массу вопросов, возникающих у проектировщиков силовых установок и инженеров по технике безопасности, но это непозволительно дорогое решение проблемы. В этом заключается одна из причин, по которой для предсказания поведения горючих смесей в конкретных условиях, представляющих интерес, исследователи концентрируют усилия, в основном, на численном моделировании и полуэмпи-рических методах [15].

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

передается теплопроводностью, нагревающей несгоревший газ, и диффузией, поставляющей в него химически активные частицы из горящего слоя. Явление перехода горения в детонацию (ПГД) заключается в том, что при распространении пламени в длинной трубе скорость его постепенно увеличивается, и на некотором расстоянии от места зажигания возникает ДВ. Преддетонационное расстояние зависит не только от состояния исходной смеси, но и от гидродинамических условий, при которых происходит распространение пламени, от диаметра трубы, состояния стенок (гладкие, шероховатые) и т. д. [13]. Как отмечается многими исследователями [16], механизм ПГД - одна из основных нерешенных проблем современной теории горения и детонации. Связано это, прежде всего, с тем, что качественное предсказание ПГД в реагирующих газах - очень сложная научная задача, которая требует знаний о динамике взаимодействия пламен, ударных волн, пограничных слоев, турбулентности и т. д. При этом с точки зрения практических приложений, нежелательный, неконтролируемый ПГД имеет огромный разрушительный потенциал. Классический механизм ПГД в прямой трубе включает несколько стадий, а именно: 1) вынужденное зажигание смеси с образованием ламинарного пламени, 2) прогрессирующее увеличение скорости горения вследствие проявления неустойчивости и последующей турбулизации течения перед фронтом пламени, 3) образование и усиление УВ перед ускоряющимся фронтом пламени и 4) самовоспламенение ударно-сжатой смеси в области между УВ и фронтом пламени (образование «горячих точек» и «взрыв во взрыве»), приводящее к образованию пересжатой ДВ, а затем и 5) самоподдерживающейся детонации ЧЖ [17].

Другой подход к инициированию детонации заключается в том, чтобы на конечном этапе ПГД, когда перед пламенем уже сформировалась УВ, воздействовать на эту УВ внешними факторами и перевести ее в детонационную, не ожидая, пока сработает механизм образования «горячих точек». В [18] для этих целей было предложено использовать электрические разряды принудительного зажигания, распределенные по длине трубы. Тщательная синхронизация разрядов с моментами прохождения УВ обеспечивает инициирование детонации в трубах с гладкими стенками на коротких расстояниях без дополнительной турбулизации потока. При этом суммарная энергия, выделенная в разрядах, существенно ниже, чем энергия прямого

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

5. ИНИЦИИРОВАНИЕ ДЕТОНАЦИИ В ТРУБЕ С ПАРАБОЛИЧЕСКИМ СУЖЕНИЕМ И КОНИЧЕСКИМ РАСШИРЕНИЕМ

В [20] на основе вычислительного эксперимента в осесимметричной постановке с помощью программного комплекса [4] был выявлен механизм инициирования детонации в трубе с параболическим сужением и коническим расширением, а также обнаружена «оптимальная» форма параболического сужения, обеспечивающая инициирование детонации для минимального числа Маха инициирующей УВ 2,65 для блокировки трубы 0.75. Под блокировкой трубы понимается БЯ =1 - (й / П)2, где П -диаметр максимального по площади сечения трубы, а й - минимального. В качестве топлива использовалась стехиометрическая пропановоздушная смесь. В дальнейшем результаты этих вычислительных экспериментов были подтверждены в натурных опытах [21]. Для «оптимальной» формы параболического сужения из [20] и угла раствора конуса расширительной секции 20° было продемонстрировано, что существует критическое значение числа Маха инициирующей УВ 2,85, выше которого наблюдается переход ударной волны в детонационную, а ниже - нет. Для оптимизации процесса инициирования в экспериментальной работе [22] была существенно увеличена длина расширяющейся конической секции. При этом критическое число Маха инициирующей УВ снизилось до величины 2,0.

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

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

условиях (см. рис. 3). Профиль стенки в секции 2 характеризуется квадратичной зависимостью ^(г), которая строится однозначно так, чтобы получались: а) заданная величина блокировки трубы БЯ; б) заданное значение угла наклона профилированного элемента ф; в) фокус полученной параболы лежал на оси симметрии трубы. Для исследования была выбрана «оптимальная» форма сужения ф = 45°, БЯ = 0,75 [20]. Предполагается, что движение в трубе возникает в результате вхождения в нее инициирующей УВ с числом Маха М и нулевым градиентом параметров непосредственно за ее фронтом.

1 г *Г и г(П\ 2* 3

1— _._ц * С) "ч! ось симметрии

а = 10°{

£| = 4 см стенка 1, = 7 см

Рис. 3. Схема профилированной трубы с параболическим сужением и коническим расширением

В проведенных расчетах используются гек-саэдрические сетки с числом ячеек до 20 млн.

Результаты трехмерного численного исследования в целом подтверждают основные этапы и особенности механизма инициирования детонации в трубе с параболическим сужением и коническим расширением, выявленные в расчетах в осесимметричной постановке в [20]. Рис. 4 иллюстрирует результаты трехмерного численного исследования для случая ф = 45°, БЯ = 0,75, М = 2,8. Для трехмерной визуализации описываемых процессов используются изоповерхности плотности пропана, а также поля давления и температуры в продольных сечениях трубы и на указанных изоповерхностях. Изозначения в каждом конкретном случае выбираются для наиболее наглядной визуализации описываемой стадии инициирования. Время отсчитывается от момента вхождения в трубу инициирующей УВ.

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

расширения (см. рис. 4, а). Важным параметром на этой стадии является величина угла а, которая была существенно уменьшена в [22] для оптимизации процесса инициирования.

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

6. ИНИЦИИРОВАНИЕ ДЕТОНАЦИИ В ВИНТОВОЙ ТРУБЕ

Исследуется труба постоянного круглого сечения с диаметром 28 мм, состоящая из трех секций: входной секции Sb винтовой секции S2 и выходной секции S3. Винтовая секция в свою очередь состоит из трех частей. Первая часть S1 - поверхность, заметаемая окружностью при

ее движении, разбиваемом на три составляющие: равномерное движение вдоль оси X, поворот с постоянной скоростью на 90° и равномерное удаление центра окружности от оси Х. Вторая часть S22 представляет собой цилиндрическую винтовую поверхность с оборотом на 360°, намотанную на ось Х. Третья часть S| аналогична первой, только центр окружности при движении не удаляется, а равномерно возвращается к оси Х (см. рис. 5). Таким образом, поверхность винтовой секции трубы заметается окружностью в плоскости, параллельной YZ, центр которой движется по закону:

x = 28• m, mе [0;6];

28• m • cos(90°-m), mе[0;1), y = < 28 • cos (90° • m), m е[1;5),

28 •( 6 - m) cos (90°^ m), m е [5;6];

28 • m • sin (90° • m), m е [0;1), z = < 28 • sin (90° • m), m е[1;5),

28•(6-m)• sin(90°• m),mе [5;6].

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

строения ячеек в круглом сечении используется О-сетка. Полученная сетка имеет около 5 млн узлов.

В начальный момент времени в секцию входит плоская инициирующая УВ с параметрами за ее фронтом, соответствующими параметрам за фронтом УВ с числом Маха М, аналогично постановке задачи в предыдущем разделе. Исследовались случаи М = 3,4, М = 3,2, М = 3,0. На стенках установлено граничное условие не-протекания.

б

Рис. 4. Механизм инициирования детонации: а - стадия реинициирования детонации, изоповерхность плотности пропана в 0,05 кг/м3 с нанесенным на нее полем давления, а также поле давления в МПа в продольном сечении, момент времени 95 мкс; б - структура детонационной волны в выходном сечении трубы, изоповерхность плотности пропана в 0,05 кг/м3 с нанесенным на нее полем температуры в градусах Кельвина в выходном сечении трубы, момент времени 135 мкс

а

Рис. 5. Схема винтовой трубы

Результаты расчетов показали, что в трубе описанной геометрии удается инициировать детонацию во всем сечении при числе Маха инициирующей УВ, превышающем М = 3,4. В численном эксперименте с М = 3,2 наблюдалась детонация в так называемом спиновом режиме, который занимает промежуточное положение между детонацией и дефлаграцией и наблюдается вблизи детонационных пределов [1]. При М = 3,0 инициирование детонации происходило в небольшом объеме, но детонационный процесс быстро затухал, наблюдался режим медленного горения.

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

варительно нагретой инициирующей ударной волной смеси при взаимодействии со стенкой трубы приводит к возникновению второго очага самовоспламенения при числах Маха инициирующей УВ М = 3,4 и М = 3,2 (см. рис. 6). Интенсивности же инициирующей УВ с М = 3,0 оказалось недостаточно для возникновения второго очага самовоспламенения.

В связи с увеличением объема самовоспламенившейся смеси при М = 3,4 и 3,2 возникают вторичные УВ и фронт лидирующей УВ ускоряется.

При выходе ДВ в выходную цилиндрическую секцию при М = 3,4 детонация занимает все сечение, при этом детонационный фронт имеет сложную форму. При выходе же детонационной волны при М = 3,2 наблюдается спиновая детонация (см. рис. 7). Многочисленные экспериментальные исследования показывают, что траектория «головы» спина относительно трубы - спираль, наклоненная под углом примерно 45° к образующей [2]. В экспериментах такую спираль регистрируют следовым методом, покрыв перед опытом стенки трубы тонким слоем сажи; тогда ДВ оставляет на саже следы, вычерчиваемые тройными точками. В численном эксперименте подобная картина наблюдается, если отобразить на поверхности трубы поле максимальных давлений, т. е. в каждой точке вывести максимальное за все время расчета значение давления в ней. Для случая М = 3,2 угол наклона составляет примерно 50°.

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

Рис. 7. Структура спиновой детонации в выходной секции винтовой трубы. Поле максимумов давления на поверхности трубы, изоповерхность температуры 298 К (фронт лидирующей УВ), изоповерхность давления 102 атм («голова» спина). Момент времени 107,5 мкс

ЗАКЛЮЧЕНИЕ

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

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

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

ной секции распространяется в спиновом режиме. В случае дальнейшего ослабления инициирующей УВ инициирования детонации не происходит.

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

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

1. Васильев А. А., Митрофанов В. В., Топчи-

ян М. Е. Детонационные волна в газах // Физика горения и взрыва. 1987. № 5. С. 109-130.

2. Митрофанов В. В. Детонация гомогенных и гетерогенных систем. Новосиб.: ИГиЛ СО РАН, 2003. 200 с.

3. Белоцерковский О. М. Численное моделирование в механике сплошных сред. М.: Физматлит, 1994. 448 с.

4. Семенов И. В., Уткин П. С., Марков В. В.

Численное моделирование двумерных детонационных течений на многопроцессорной вычислительной технике // Вычислительные методы и программирование. 2008. Т. 9. С. 119-128.

5. Frolov S. M., Aksenov V. S., Shamshin I. O. Detonation propagation through U-bends // Nonequilibrium processes, Vol. 1: Combustion and Detonation / Eds. G. Roy et al. 2005. P. 348-364.

6. Barth T., Ohlberger M. Finite Volume Methods: Foundation and Analysis // Encyclopedia of Computational Mechanics. 2004. Vol. 1. P. 439-470.

7. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с.

8. Численное решение многомерных задач газовой динамики / С. К. Годунов [и др.]. М.: Наука, 1976. 400 с.

9. Karypis G., Kumar V. METIS 4.0: Unstructured graph partitioning and sparse matrix ordering system // Technical report, Department of Computer Science, University of Minnesota. 1998. URL: http://glaros. dtc.umn.edu/gkhome/views/metis (дата обращения 05.04.2010).

10. Walshaw C., Cross M., Everett M. G. A localized algorithm for optimising unstructured mesh partitions // International Journal of Supercomputer Applications. 1995. Vol. 17, No. 4. P. 280-295.

11. Preis R., Diekmann R. The party partitioning-library // Technical report TR-RSFB-96-024, Universi-tat-GH Paderborn. 1996.

12. Pellegrini F. PT-Scotch and LibScotch 5.1 // Technical report LaBRI. 2008. URL: http://www.labri.fr/perso/pelegrin/scotch/ (дата обращения 05.04.2010).

13. Зельдович Я. Б., Компанеец А. С. Теория детонации. М.: Госуд. изд-во техн.-теор. лит-ры, 1955. 268 с.

14. Седов Л. И., Коробейников В. П., Марков В. В. Теория распространения взрывных волн // Труды Математического института им. В. А. Стекло-ва РАН. 1986. № 175. С. 178-216.

15. Борисов А. А. Инициирование детонации в газовых и двухфазных смесях // Импульсные детонационные двигатели / Под ред. д.ф.-м.н. С. М. Фролова. М.: Торус Пресс, 2006. С. 159-186.

16. Oran E. S., Gamezo V. N. Origins of the Def-lagration-to-detonation Transition in Gas-phase Combustion // Combustion and Flame. 2007. Vol. 148. P. 4-47.

17. Фролов С. М. Быстрый переход горения в детонацию // Химическая физика. 2008. Т. 27, № 6. С.31-44

18. Frolov S. M. Initiation of Strong Reactive Shocks and Detonation by Traveling Ignition Pulses // Journal of Loss Prevention. 2005. Vol. 19, № 2-3. P. 238-244.

19. Сокращение длины и времени перехода горения в детонацию в трубе с профилированными регулярными препятствиями / С. М. Фролов [и др.] // Доклады Академии наук. 2007. Т. 415, № 4. С. 509513.

20. Семенов И. В., Уткин П. С., Марков В. В.

Численное моделирование инициирования детонации в профилированной трубе // Физика горения и взрыва. 2009. Т. 45, № 6. С. 73-81.

21. Numerical and experimental investigation of detonation initiation in profiled tubes // I. V. Semenov [et al.] // 22nd ICDERS Proceedings. 2009. CD, Paper No. 168.

22. Фролов С. М., Аксенов В. С. Инициирование газовой детонации в трубе с профилированным препятствием // Доклады Академии наук. 2009. Т. 427, № 3. С. 344-347.

23. Фролов С. М. Импульсное детонационное горение: новое поколение энергетических установок // Энергетика. 2008. Т. 3, № 41. С. 44-45.

ОБ АВТОРАХ

Семенов Илья Витальевич,

ст. науч. сотр. Ин-та автоматизации проектирования РАН. Дипл. приклад. матем. (Даль-невост. гос. ун-т, 1995). Канд. физ.-мат. наук по матем. мо-дел., числ. методам и компл. прогр. (МФТИ, 2005). Иссл. в обл. механики реагир. много-фаз. сред, паралл. вычисл.

Ахмедьянов Ильдар Фидаи-левич, мл. науч. сотр. того же ин-та. Дипл. матем. (МГУ, 2005). Иссл. в обл. параллельных вычислений.

Уткин Павел Сергеевич,

асп., мл. науч. сотр. того же ин-та, асс. каф. вычисл. матем. МФТИ. Дипл. магистр приклад. матем. и физики (МФТИ, 2008). Иссл. в обл. механики реагир. многофазных сред.

Лебедева Алина Юрьевна,

стаж.-иссл. того же ин-та. Дипл. бакалавр приклад. ма-тем. и физики (МФТИ, 2009). Иссл. в обл. механики реагир. многофазных сред, построения криволинейных сеток.