Содержание к диссертации
Введение
ГЛАВА 1. Математические модели динамики газовых и пылегазовых сред 18
1.1 Основные предположения 18
1.2 Физико-химические процессы в газовой фазе 19
1.3 Физико-химические процессы в фазе частиц 22
1.4 Механизмы взаимодействия фаз 24
1.5 Полная система уравнений динамики для пылегазовой смеси 27
1.6 База данных для рассматриваемых физико-химических процессов 28
1.6.1 Кинетика горения угольной пыли 28
1.6.2 Кинетика горения водорода 30
1.6.3 Кинетика испарения частиц льда 31
1.6.4 Физико-химические константы 31
ГЛАВА 2. Методика численного решения и программная реализация модели 33
2.1 Принципы построения численного алгоритма 33
2.2 Расчётные сетки 35
2.3 Расчёт конвективного переноса для параметров фаз 37
2.4 Расчёт межфазного взаимодействия и химических процессов 38
2.5 Расчёт величины шага по времени 38
2.6 Задача Сода 38
2.7 Тестирование модели кинетики горения угольных частиц 40
ГЛАВА 3. Результаты вычислительных экспериментов 41
3.1 Моделирование разлёта фрагментов ледяного метеорита 41
3.1.1 Определение начального поля параметров 41
3.1.2 Результаты расчётов для цилиндрического тела 46
3.1.3 Результаты расчётов для сферического тела 58
3.1.4 Выводы 70
3.2 Моделирование движения, газификации частиц и горения пылегазовои смеси за ударной волной 72
3.2.1 ПОДЪЁМ И распыление инертной дисперсной фазы из плотного слоя 72
3.2.2 Формирование течения горючей пылегазовои смеси в угольной шахте большой протяженности 77
3.2.3 Развитие слоистой детонации в трубах 81
3.2.4 Выводы 89
3.3 Моделирование работы электрохимического пульсирующего детонационного двигателя 90
3.3.1 Мелкомасштабная модель 90
3.3.2 Запуск крупномасштабного двигателя 97
3.3.3 Выводы 105
Заключение 106
Литература
- Физико-химические процессы в фазе частиц
- База данных для рассматриваемых физико-химических процессов
- Расчёт конвективного переноса для параметров фаз
- Моделирование движения, газификации частиц и горения пылегазовои смеси за ударной волной
Введение к работе
Актуальность работы. В современном мире аппарат математического моделирования нашел широкое применение в различных областях жизнедеятельности человека. Математическое моделирование и вычислительный эксперимент позволяют значительно сократить затраты на исследование сложных динамических процессов и разработку новых технических устройств, а также дают возможность получить характеристики динамики таких явлений, для которых натурный эксперимент очень трудно провести или наблюдать. Создание вычислительных комплексов программ, ориентированных на использование СуперЭВМ, для решения широкого класса сложных многомерных задач физико-химической динамики многокомпонентной среды является очень трудным, но перспективным делом. В последнее десятилетие вычислительная мощность современных СуперЭВМ с параллельной архитектурой увеличилась настолько, что стало возможным применить ранее разработанные математические модели и вычислительные алгоритмы, усложнённые введением новых элементов моделей, для более детального исследования динамики быстропротекающих физико-химических процессов в многокомпонентных смесях. В качестве исследуемых процессов, представляющих интерес, как с научной, так и с практической точки зрения, можно отметить следующие:
Физико-химические процессы, происходящие во время взаимодействия вещества метеорита с атмосферой, являются очень сложными и представляют несомненный интерес для оценки возможных последствий вторжения космических тел в атмосферу Земли и других планет. Применение аппарата математического моделирования для исследования данных процессов имеет приоритет из-за невозможности постановки натурного эксперимента.
Физико-химические процессы движения, газификации, воспламенения, горения и детонации пылегазовой смеси за распространяющейся ударной волной. Они имеют место при взрывах в угольных шахтах и элеваторах для зерна, при извержении вулканов и в ряде промышленных и экспериментальных устройств. Детальное изучение этих явлений требуется для выработки условий безопасной работы, в процессе которой могут возникать взрывоопасные концентрации горючей пыли.
Физико-химические процессы, которые имеют место при работе электрохимического пульсирующего детонационного двигателя. Данный тип двигателя появился сравнительно недавно, и изучение особенностей его работы представляется актуальным, как с точки зрения развития теории детонации, так и с точки зрения практической эксплуатации устройств подобного типа. В качестве топлива для установки используется водород, применение которого в двигательных установках, также представляет большой практический интерес, в том числе и по экологическим причинам.
Цели и задачи диссертационной работы. Разработка комплекса алгоритмов и программ для численного моделирования пространственных нестационарных физико-химических процессов в многокомпонентных одно и двухфазных смесях. Применение данного вычислительного комплекса для следующих задач:
Изучение физико-химических процессов происходящих при движении, газификации, воспламенении, горении и детонации пылегазовой смеси за распространяющейся ударной волной в круглых трубах и плоских каналах.
Исследование движения и испарения облака фрагментов разрушенного метеорита в атмосфере Земли.
Изучение функционирования электрохимического пульсирующего детонационного двигателя.
Научная новизна работы.
1) Предложены математические модели для исследования следующих быстропротекающих процессов в многокомпонентных смесях:
Взрывного распада разрушенного ледяного метеорита, движущегося в атмосфере Земли.
Формирования детонационной и квазидетонационной волн в каналах и трубах с плотными пристеночными пылевыми слоями.
Функционирования электрохимического пульсирующего детонационного двигателя оригинальной конструкции.
Создан вычислительный комплекс для решения многомерных задач физико-химической динамики многокомпонентной среды. В нём реализованы технологии расщепления и распараллеливания составляющих его вычислительных алгоритмов.
Проведено численное исследование, указанных в пункте 1, процессов физико-химической динамики многокомпонентных смесей.
Научная и практическая ценность работы. Разработанный вычислительный комплекс может использоваться для решения широкого круга научно-практических задач. В частности, прогнозирования последствий вторжения малых космических тел в атмосферы планет; оценки взрывоопасное в угольных шахтах и помещениях с наличием горючей пыли; проектирования электрохимических пульсирующих детонационных двигателей на водородном и углеводородном топливе. Блочная структура комплекса позволяет легко вводить различные подмодели, например, для описания детальной химической кинетики произвольного многокомпонентного и многофазного состава. Вычислительный комплекс адаптирован для ЭВМ с массивной параллельной архитектурой кластерного типа, что даёт возможность его эффективного применения к решению ресурсоёмких научно-технических задач.
Структура и объём диссертации. Диссертация состоит из введения, трёх глав, заключения и списка цитированной литературы. Диссертация изложена на
120 страницах, включает 2 таблицы и 51 рисунок. Список литературы содержит 132 наименования.
Во введении обосновывается актуальность темы исследования, кратко излагается содержание диссертации, указывается ее научная новизна, формулируются основные результаты работы. Даётся обзор литературы по применяемым вычислительным методам, технологиям расщепления и распараллеливания алгоритмов, системам химической кинетики, а также по экспериментальным и численным результатам, связанным с рассматриваемыми в диссертации проблемами.
В первой главе подробно изложены математические модели для описания динамики газовой и твёрдой фаз. Движение фаз рассматривается как движение взаимопроникающих и взаимодействующих континуумов, каждый из которых имеет свою скорость и температуру. Предполагается, что частицы абсолютно твёрдые, имеют сферическую форму, одинакового диаметра и состоят из однородного многокомпонентного теплопроводного материала. Для описания движения каждой из фаз применен Эйлеров подход. Для описания движения газовой фазы применяется система уравнений Навье-Стокса. Учитываются процессы переноса: диффузия, вязкость, теплопроводность, излучение. В общем случае коэффициенты переноса могут зависеть от температуры газовой фазы. Для описания движения твёрдой фазы используются уравнения Эйлера без давления, столкновениями между частицами будем пренебрегать. Учитываются силовое, тепловое, радиационное взаимодействия и обмен массой между фазами. В многокомпонентной смеси допускается возможность протекания гомогенных и гетерогенных химических реакций, скорости которых описываются Аррениусовскими зависимостями. Система замыкается уравнениями состояния для газовой и твердой фаз. Далее в первой главе описаны, используемые в работе, модели кинетики горения угольной пыли, испарения ледяных частиц и горения водорода. Приводится набор физико-химических констант, применявшихся в расчётах.
Во второй главе даётся описание принципов построения вычислительного алгоритма для решения системы уравнений, представленной в первой главе диссертации, обсуждаются применяемые расчётные сетки, метод распараллеливания вычислительного алгоритма и приводятся результаты тестовых расчётов для отдельных программных блоков комплекса.
В третьей главе приводятся результаты вычислительных экспериментов, проведённых с помощью разработанного комплекса программ. В первой части третьей главы представлены результаты расчётов разлёта разрушенного ледяного метеорита. Предлагается подход, в рамках которого разрушенный метеорит моделируется двухфазным объёмом, имеющим ос есим метричную форму. Предполагается, что фрагменты разрушенного метеорита есть ледяные сферические частицы. Рассматривается задача о движении двухфазного объёма с заданной начальной скоростью К вдоль траектории в изотермической атмосфере с переменной плотностью. Основное внимание уделяется моделированию так называемого «взрыва в полёте», в предположении, что это явление связано с разлётом фрагментов и испарением материала метеорита. Во второй части третьей главы рассмотрены задачи движения, газификации, воспламенения, горения и детонации пылегазовой смеси за распространяющейся ударной волной в плоских каналах и осесимметричных трубах с одним закрытым концом при наличии плотных слоев пыли на стенках. У закрытого конца канала или трубы имеется детонационная камера длины Lcit1 предназначенная для генерации ударной волны. Интенсивность её воздействия на процесс, при прочих равных условиях, определяется составом горючей газовой смеси, давлением и объёмом камеры. С помощью предложенной модели, исследован вопрос о подъёме и диспергировании пыли. Приводятся результаты расчётов подъёма инертной пыли с использованием различных безразмерных коэффициентов Км и Ksaf в силах Магнуса и Сэффмана. Показано наличие задержки подъёма пыли за проходящей ударной волной, а также дальнейший интенсивный подъём при учёте силы Магнуса. Рассмотрено влияние толщины пылевого слоя на процессы подъёма и диспергирования и приведены результаты расчётов для плоских каналов и круглых труб со слоями пыли разной толщины. В третьей части третьей главы приводятся результаты численного моделирования функционирования электрохимического пульсирующего детонационного двигателя. Приведена принципиальная схема работы двигателей подобного рода. Проведено исследование работы мелкомасштабной модели упрощённого варианта электрохимического пульсирующего двигателя. Оценено влияние электрического разряда на динамику процесса горения в камере сгорания.
На защиту выносятся:
1. Математические модели, описывающие следующие быстро протекающие процессы в многокомпонентных смесях:
Взрывной распад фрагментов ледяного метеорита движущегося в атмосфере Земли.
Формирование детонационной и квазидетонационной волн в каналах и трубах с плотными пристеночными пылевыми слоями.
Функционирование электрохимического пульсирующего детонационного двигателя оригинальной конструкции.
Разработанные вычислительные алгоритмы и комплекс программ для численного моделирования быстропротекающих физико-химических процессов в многокомпонентных смесях на ЭВМ с параллельной архитектурой.
Численное исследование процесса «взрыва в полёте» метеорита в предположении, что это явление связано с разлётом фрагментов и испарением материала метеорита, включающее оценку влияния различных параметров модели: начальной формы тела, возможности испарения частиц и интенсивности потери энергии с фронта головной ударной волны за счёт излучения на динамику разлёта.
Численное моделирование течений пылегазовых смесей, формирующихся за ударной волной в результате ее взаимодействия с узким пристеночным слоем пыли, позволившее выявить механизмы интенсификации горения и формирования детонационных и квазидетонационных волн. Предложение о необходимости проведения натурных экспериментов в установках, имеющих большее отношение длины к диаметру, чем используемые в настоящее время. 5. Численное моделирование течений в электрохимическом пульсирующем двигателе. Для мелкомасштабной модели исследование влияние величины секундного расхода горючего, частоты и количества разрядов, добавки горючей пыли на локальные и интегральные характеристики течения, в частности, величину силы тяги и удельного импульса. Для полномасштабной модели С. Вуйтицкого исследование по определению критической энергии инициирования детонационного горения и анализ различных сценариев запуска двигателя с учётом раздельного поступления водорода и воздуха. Оценка влияния энергии дополнительных разрядов на изменение режима работы установки.
Мощности современных вычислительных комплексов позволяют решать непосильные ранее задачи, касающиеся широкого спектра процессов в природе и технике. В течение длительного периода времени в центре внимания многих исследователей находятся проблемы моделирования тунгусской катастрофы и движения метеоритов в атмосферах планет, взрывов горючих пыле газовых смесей в угольных шахтах и на промышленных установках. Сравнительно недавно усилия теоретиков и экспериментаторов сконцентрировались на изучении возможности использования детонационного горения в экономичных и экологически чистых реактивных двигателях. Применяются традиционные экспериментальные методы на мелкомасштабных моделях и реальных объектах. Однако их возможности ограничены из-за трудоемкости, дороговизны, опасности, недостаточной информативности, а зачастую просто из-за невозможности воспроизведения, как в случае тунгусского явления. Многие трудности позволяет преодолеть численное моделирование на современных вычислительных системах с распараллеливанием вычислительного процесса и привлечением специальных методов и программ обработки данных. Поэтому вычислительный эксперимент становится основным методом исследования и предоставляет практически неограниченные возможности познания и преобразования окружающего мира. Однако вычислительный эксперимент будет давать адекватные результаты лишь в том случае, когда он базируется на хорошей математической модели того или иного явления и надежном численном алгоритме. Математическая модель строится на основе физической модели и здесь многое зависит от ее создателя — его способностей, знаний, опыта и фантазии.
Проблеме тунгусской катастрофы посвящено огромное количество работ, различающихся как по качеству физической модели, так и по математической постановке и применяемым методам исследования. Оставляя в стороне явно фантастические теории, рассчитанные на непрофессиональную аудиторию, и следуя своим научным интересам, отметим результаты группы исследователей, работавшей под руководством и в тесном сотрудничестве с В.П. Коробейниковым. Основное их достижение состоит в том, что ими была предложена теория и проведены расчёты, которые позволили воспроизвести эффект воздействия тунгусского взрыва на поверхность Земли. С целью детального изучения процесса взаимодействия малых небесных тел, метеоров и метеороидов, с атмосферами планет методом математического моделирования в работах [1,2,3] была сформулирована комплексная модель. Она включает следующие основные разделы: 1) Определение начальной части траектории, скорости и абляции тела; 2) Расчёт напряжений, деформаций и разрушения тела; 3) Движение вещества после разрушения; 4) Расчет взрывоподобного разлета остатков вещества на заключительной стадии полета; 5) Воздействие на поверхность планеты и атмосферу; 6) Определение физической природы тела и его орбиты вне атмосферы. Все разделы этой модели уже разрабатывались применительно к проблеме разрушения космических тел, в частности, в связи с
Тунгусским явлением [2-9]. В диссертации основное внимание уделяется моделированию так называемого «взрыва в полёте» ледяного метеорита, в предположении, что это явление связано с разлётом фрагментов и испарением материала метеорита. Характер этих процессов, при прочих равных условиях, может сильно зависеть от начальной формы и размера тела, размера его фрагментов, и интенсивности излучения газа, сжатого и нагретого баллистической ударной волной.
Вторая проблема, которая рассматривается в диссертации — это распространение ударных и детонационных волн в каналах и трубах с тонкими слоями угольной пыли, находящимися на их стенках. Эта проблема важна с точки зрения взрывобезопасности в угольных шахтах. Вопросы воспламенения и горения частиц угольной пыли детально рассматривались в монографиях [10, 11]. Взрывы угольной пыли интенсивно изучались экспериментально для определения принципиальной возможности формирования самоподдерживающейся детонации и условий, при которых этот процесс реализуется [12-16]. В последние годы наиболее значительные экспериментальные результаты получены авторами работ [17-19]. Теоретические исследования проводились в основном для газовзвесей предварительно распыленных частиц и касались определения критических условий - состава смеси, ее концентрации и мощности инициирующего процесс источника энергии, при которых может реализоваться быстрое горение в форме, так называемой, квазидетонации или в чисто детонационном режиме [20-23]. Вопросы моделирования процесса диспергирования пыли из слоя теоретически рассматривались в работах [24-27]. Однако полного решения проблемы с учетом горения и реальных размеров объектов, представляющих интерес, дано не было.
Последним объектом исследования в диссертации является электрохимический пульсирующий детонационный двигатель оригинальной конструкции польского профессора С. Вуйтицкого, который тесно сотрудничал с B.IL Коробейниковым. В основе функционирования этого устройства лежит концепция его автора о энергодинамике автоцикломата [28-31]. Энергодинамика автоцикломата — это класс моделей тепловых двигателей, характеризуемых отсутствием активных подвижных частей (подобных поршню и турбине), автоматическим циклическим функционированием в стационарных условиях, саморегулированием, основанном на сильной обратной связи, и использованием физических процессов для динамического взаимодействия элементов и функционирования машины в целом. Для примера, упомянем следующие устройства: водяной пульсирующий двигатель Мак Хага, используемый для создания тяги при движении судна, и содержащий паровой котел, паровой движитель, водяное сопло и систему труб (запатентован в 1916 году в США); сверхзвуковой воздушно-реактивный двигатель, разработка которого связана с конструированием двигательной системы для аэрокосмических самолетов; пульсирующий лазерный двигатель для маленьких космических ракет; горелка Геймгольца, которая используется как нагреватель и как элемент пульсирующего двигателя; электрохимический пульсирующий двигатель для создания тяги. Указанным системам присущи следующие характерные особенности: интенсивный теплоподвод, являющийся движущей силой в цикле системы; расширение газов в камере сгорания, в результате которого возникает частичный вакуум и условия для втекания свежей горючей смеси и инициирования нового цикла; перезарядка горючей смесью; клапанная система, реагирующая на градиент давления вблизи входа в камеру сгорания; инициирование горения и детонации; отсутствие активных подвижных частей; саморегулирование при стационарных условиях; динамическое взаимодействие с окружающей средой; система соединяющих труб специальной геометрии, которые помогают организовать циклическое функционирование. Эти особенности составляют ядро признаков, которые характеризуют автоцикломаты, управляемые процессом горения. Электрохимический пульсирующий реактивный двигатель является примером пульсирующего двигателя, в котором успешно осуществляется идея сочетания электрической и тепловой энергии. Здесь электрический разряд (интенсивная искра), происходящий в камере сгорания, используется для повышения эффективности превращения тепловой энергии в механическую посредством увеличения степени сжатия двигателя. Этот электрический разряд, значительно интенсивней разряда в обычном двигателе внутреннего сгорания, и может превратить горение в детонацию и взрыв. Разряд наступает в момент времени, когда фронт пламени соприкасается со специальными электродами, подсоединенными к заряженному конденсатору. Электроды, имеющие профиль Роговского, создают электрическое поле, которое усиливается вблизи их кромок. Кольцевой разряд создает в горючей смеси сходящиеся ударные волны, благодаря которым формируется быстрый процесс горения с детонацией, высоким давлением, что приводит к увеличению степени сжатия. Следует заметить, что энергия, выделяющаяся при горении топлива в течение одного цикла работы двигателя, много больше энергии, подведенной при электрическом разряде. Следует отметить, что впервые идея использования детонации в двигателях была высказаны ЯЛ. Зельдовичем в [32].
Для описания рассматриваемых в диссертации двухфазных сред используется модель взаимопроникающих и взаимодействующих континуумов, предложенная Х.А. Рахматулиным [33], развитая и обобщенная для химически активных частиц в работах [34-36]. Вопросы теории моделей многофазных сред разрабатывались многими авторами [37-44]. Численно решались интересные с точки зрения теории и практики сложные задачи о распространении взрывных и детонационных волн [45-55]. В частности, в результате численных экспериментов был обнаружен и подтвержден в специальных натурных экспериментах эффект формирования слоя с высокой концентрацией твердых частиц пыли, так называемого ро-слоя, за ударными и детонационными волнами, распространяющимися по однородной пылегазовой смеси. Как оказалось, это явление является основной причиной формирования детонации в пылегазовых смесях, в которых концентрация горючих частиц существенно ниже общепринятых предельных значений с точки зрения взрывоопасное.
С проблемами горения пыли и моделирования работы электрохимического пульсирующего двигателя тесно связаны вопросы выбора кинетики химических реакций в многокомпонентных реагирующих смесях [56-72] и вопросы инициирования самоподдерживающейся газовой детонации, в частности, с помощью электрического разряда [73-80]. Отметим уникальные результаты, полученные в рамках простой модели электрического разряда, которые позволили объяснить природу сильного роста с увеличением продолжительности разряда величины минимальной (критической) энергии, приводящей к формированию детонации. Практическая рекомендация как следствие исследования — использование коротких и достаточно энергоемких разрядов. Минимальная критическая энергия получается при мгновенном разряде. Кроме того, была установлена возможность возникновения детонации при фокусировке даже достаточно слабых сходящихся ударных волн.
В связи с потребностями практики и прогрессом вычислительных систем интенсивно развивались численные методы расчета течений с учетом сложных физико-химических процессов. Академиком О.М. Белоцерковским введено ёмкое понятие вычислительного эксперимента как мощного современного аппарата решения задач механики сплошных сред [81]. Широкое распространение и развитие получила схема С.К. Годунова [82-91]. В работах по численным методам разрабатываются, в частности, проблемы повышения точности разностных схем [92], устойчивости и решения, так называемых, жестких систем, к которым относятся и уравнения химической кинетики. Для решения сложных задач требуются большие объемы памяти и повышаются требования к быстродействию. В настоящее время, широкое распространение получил метод объединения нескольких одно- или многопроцессорных вычислительных систем в одну с помощью высокоскоростной сети обмена информацией. Использование таких систем связано с распараллеливанием вычислительного алгоритма и одновременным исполнением его параллельных частей на разных процессорах. С параллельным программированием связано много важных как теоретических, так и чисто практических вопросов [93]. Анализу алгоритмов на предмет их возможного распараллеливания посвящены работы Воеводина [94].
Достоверность и обоснованность полученных результатов подтверждается использованием в исследованиях общих уравнений механики многофазных сред, применением классического эйлерова подхода для описания фаз, сравнением результатов с экспериментальными данными, полученными независимо от автора.
Апробация результатов работы. Материалы, отражающие содержание диссертационной работы, с достаточной полнотой опубликованы в работах [106-132], докладывались на семинарах Института автоматизации проектирования РАН и на следующих научных конференциях: Международном коллоквиуме «Достижения в экспериментальной и вычислительной детонации» (International colloquium on advances in experimentation and computations of detonations, Санкт-Петербург, Россия, 1998); 17, 18 и 19 Международных коллоквиумах по динамике взрывов и реагирующих систем (ICDERS Heidelberg, Germany, 1999; Seattle, Washington, USA, 2001; Hakone, Japan, 2003); US/European Celestial Mechanics Workshop (Poznan, Poland, 2000); International Colloquium on Control of Detonation Processes (Moscow, Russia, 2000); 28th International Symposium on Combustion (Edinburgh, UK, 2000); 3 и 4 International Symposium on Hazards, Prevention, and Mitigation of Industrial Explosions (Tsukuba, Japan, 2000; Bourge, France, 2002); Ninth international conference on numerical combustion (Sorrento, Italy, 2002); Всероссийской конференции "Актуальные проблемы прикладной математики и механики", посвященной 70-летию академика А.Ф.Сидорова (Екатеринбург, Россия, 2003); International Russian-Indian workshop high performance computing in science and engineering (HPC SE 2003, Москва, Россия, 2003); International colloquium on application of detonation processes (S.-Petersburg, Russia, 2004); Международной конференции по горению и детонации (Мемориал Зельдовича II, Москва, Россия, 2004); Russian-Japan International Workshop on Turbulence and Instabilities (Moscow, Russia, 2004); Межинститутском семинаре «Проблемы математического моделирования рабочего процесса энергопреобразующих устройств» (ИХФ РАН, Москва, 2004).
Автор выражает глубокую благодарность за постоянное внимание, поддержку и совместную научную работу своим научным руководителям: члену-корреспонденту РАН Виктору Павловичу Коробейникову и д.ф.-м.н. Владимиру Васильевичу Маркову, а также академику Олегу Михайловичу Белоцерковскому, д.ф.-м.н. Валентину Анатольевичу Гущину и к.ф.-м.н. Игорю Станиславовичу Меньшову, н.с. Дмитрию Владимировичу Воробьёву за техническое сотрудничество, с.н.с. Елене Раилевне Павлюковой за помощь в решении организационных вопросов, а также всему коллективу ИАП РАН за доброжелательное отношение и содействие в работе.
Физико-химические процессы в фазе частиц
Фаза частиц характеризуется плотностями ее компонентов, температурой и вектором скорости. Приведем систему нестационарных уравнений движения для фазы частиц в прямоугольной декартовой системе координат х (i = 1,...,3).
Уравнения сохранения массы имеют вид: +А(р „JU0 ,-1,2,3 , k = l,..,NSSP (1.3.1) at дх к где и2 обозначают компоненты вектора скорости, 32к есть скорость изменения плотности к-го компонента в результате физико-химических процессов. Уравнение для числовой плотности частиц п имеет вид: — + Д-(пи1) = 0, ( = 1,2,3 (1.3.2) dt 8хЛ 1} V Уравнение изменения количества движения имеет вид: H + T(p2ui2ui) = G 2 , / = 1,...,3 ,7 = 1,-,3 (1.3.3) dt oxJ где Gl2 - скорость изменения импульса в / направлении за счет внешних объёмных воздействий. Уравнение баланса для полной энергии имеет вид: b. + JL(E1u2l) = Q2 , у = 1,2,3 (1.3.4) dt oxJ где Е2 - полная энергия, Q2 - скорость изменения полной энергии за счет внешних объёмных воздействий будет определена в 1.4. Полная энергия для фазы частиц определяется следующим выражением: З 2 NSSP Е2 = 0.5р2(«5) + р2ке2к, где е2к - внутренняя энергия к-то компонента фазы частиц. Уравнения состояния берутся в виде: еу=4+с2/2 , k = l,...,NSSP (1.3.5)
Здесь с2к - удельная теплоёмкость к-то компонента, Т2 - температура частиц.
В работе будут рассматриваться два типа частиц: угольные и ледяные. Как известно, угольные частицы имеют сложный многокомпонентный состав. Основную массу угольной частицы составляет углеродный скелет, кроме того, угольные частицы содержат в себе так называемые летучие, которые при нагревании (процессе пиролиза) частицы переходят в газовую фазу. В состав летучих входят углеводороды, вода, водород, сера, углекислый газ и др. Весь многокомпонентный состав угольной частицы разобьем на две группы. Первая для к из диапазона 1 к NVSP NSSP будет содержать летучие компоненты. Вторая для к из диапазона NVSP k NSSP будет содержать нелетучие компоненты, а именно углеродный скелет, состоящий из горючих и негорючих компонентов. Для описания скорости изменения плотности летучих в процессе пиролиза применяется следующая формула: гл = р2явк ехр(- /ЛГ2) (1-3.6) где Bk,Evk - есть эмпирические константы, R - универсальная газовая постоянная. Скорость изменения плотности горючей составляющей углеродного скелета частицы определим следующим образом: г. -ЛОДОиРсА ехр(- /ЛГ2) (1-3.7) где AsiFs,EsJl - есть эмпирические константы, р0 ,Ла - плотность и газовая постоянная кислорода.
Процесс испарения ледяных частиц метеорного тела рассматривается с помощью следующей упрощённой схемы. Предполагаем, что тепло, поступающее к частице извне, идёт на прогрев до температуры плавления некоторой массы частицы, на фазовый переход этой массы в жидкое состояние, её прогрев до температуры испарения и фазовый переход в парообразное состояние и на нагрев образовавшегося пара до температуры газа, окружающего частицу. Считается, что температура испарения зависит от давления газа, окружающего частицу. Предполагается, что между газовой и твёрдой фазами происходит силовое и тепловое взаимодействие. Силовое взаимодействие учитывается введением в систему уравнений, прежде всего, силы сопротивления F. Её проекции на оси координат можно записать в виде: ж А2 F1 = п—Свр\и[-и[)ЩиІ-4) , 1 = 1,2,3 (1.4.1) где d - диаметр частиц, CD - коэффициент сопротивления, зависящий от относительных чисел Рейнольдса Re и Маха М. В настоящей работе для CD использовались следующие эмпирические выражения: предложенное Клячко [95], и формулы Хендерсона [96], учитывающие сжимаемость газа и температуры фаз. Для М 1.0: cn =
Обмен теплом между фазами происходит за счёт теплопроводности и радиационного потока, можно описать следующими соотношениями: QT=n;rdZl iu(T, T2) , д=п я tfcr T -Г24) . (1.4.5)
Здесь Я коэффициент теплопроводности (зависит от температуры) газовой фазы, єх (0 г, 1) - эмпирический коэффициент, ств - постоянная Стефана-Больцмана, Nu - число Нуссельта, являющееся эмпирической функцией, зависящей от Re и числа Прандтля Рг, по Кнудсену и Катцу [97] вида: Nu =2+0.6 Re0 5Ргш (1.4,6) и от числа М, для Re 10000, М 6 по Фоксу [98] в форме: зз D .55l + 0.5exp(-17jlf/Re) 2ехр(-М) (1.4.7) Nu = "" гч "; + 0.459 Pr-jj Re 1 + 17M/Re 1.5
Отмеченные выше механизмы взаимодействия фаз имеют место независимо от того, присутствует или нет обмен массами между фазами. Кроме них, при наличии этого обмена, следует учесть обмен импульсом и энергией между фазами, переносимых веществом, переходящим из одной фазы в другую.
Рассматривая задачу о подъёме и диспергировании пыли из плотного слоя, оказалось необходимым ввести в рассмотрение специальные математические модели способные описать эти процессы. В.П. Коробейниковым была предложена модель, учитывающая, кроме силы сопротивления, ещё силы Магнуса и Сэффмана, которые возникают из-за вращения частиц в завихренном потоке газа [101,102]. В дальнейшем предполагается, что частицы приводятся во вращение вихрями в газе и средняя угловая скорость вращения частиц пропорциональна вектору вихря в газе. В этих предположениях сила Магнуса может быть представлена так: рл/=п уЛ[(иі-и2)хг и,] (1.4.8) где Км - безразмерный коэффициент, и, и и2 - векторы скорости в газовой и твердой фазах соответственно. Величина этого коэффициента определяется дополнительными исследованиями или эмпирически. Сделанные оценки показывают, что эта величина находится в пределах от 0.3 до 20.
База данных для рассматриваемых физико-химических процессов
Для построения конечно-разностной схемы, рассматриваемой в предыдущей главе математической модели, использовался метод конечных объёмов. Вся вычислительная область разбивалась на множество контрольных объёмов. Перепишем систему уравнений в векторном виде:
fa + dL=d%L + n + G / = 1,2,3 (2.1.1) dt дх, дхі Здесь г 1 2 3 гч 1 2 3 п т7" Z = [ р1 ji, ",P\Mi р\Щ P[U\ р » p2,k" - )П,р2 2" Pl li Р2 2 2Л fi=[... plku l,...,piulu l+Shp,piufuti+S2 p,p1ufu i +5llp,{Ex + p)u[, yLA іHj П Wj r 2 2 2 f 2 2 2 / 2 2 2 I y o J 3 a n ЛОТР Н=Г 0 -F -F -F1 F2-F2-F2 -F3-F3-FJ F +F + Fi+F F + +F + a+eJ i=l G=[...,5lJt,...,-2Ki,-2Mj,-fi 2«2,-- 2Z("0 - X -(l-a)Gc» I=I =i (2.1.2) 1 3 .ч2 NSSP ...fc32i,...,0,c52M2,w2M2, 2H2J}-ft?2 (M2) + Y,2,ke2,k-aQcY (=1 Jbl Результирующую систему дискретных уравнений можно записать в следующем виде: (2.1.3) f,r-f rj"i » g ,=g„,-"f . / = 1,-..,3 ,7 = 1,..., ЛГ
Здесь N - число контрольных объёмов во всей вычислительной области, z - -вектор решения, осреднённый Haj-ом контрольном объёме Пу.
Численное решение проводилось с использованием метода конечных объёмов и метода расщепления по физическим и химическим процессам. Разработанная явная вычислительная схема имеет второй порядок по пространственным переменным и первый по времени. Решение исходной системы дискретных уравнений находится в несколько этапов. На первом этапе находится промежуточное решение, учитывающее только процессы конвективного переноса. На втором этапе рассчитывается тепловое и силовое межфазные взаимодействия и проводится коррекция промежуточного решения. На третьем этапе рассчитывается интенсивность химических процессов внутри и между фазами (гомогенное и гетерогенное горение, выход летучих, испарение) и определяется окончательное решение на новом временном слое.
Расчёт одного шага по времени осуществляется в следующей последовательности:
1. Вычисляются значения искомых функций в фиктивных граничных ячейках, необходимые для расчёта вязких потоков через границы расчётной области;
2. Вычисляются производные искомых функций внутри расчетной области с использованием ограничителя MINMOD;
3. Вычисляются значения искомых функций и их производных на границах расчётной области, необходимые для определения невязких потоков через границы расчётной области;
4. Вычисляются вязкие и невязкие потоки искомых функций через границы расчётных ячеек и определяются предварительные значения решения на следующем временном шаге;
5. В каждой расчётной ячейке решается система линейных уравнений, определяющая влияние межфазного силового и теплового взаимодействий на значения скоростей и температур фаз на следующем временном шаге;
6. В каждой расчётной ячейке решается система кинетических уравнений, и получаются окончательные значения плотностей компонентов, скоростей и температур фаз;
7. По критерию устойчивости определяется величина нового временного шага.
Для численного моделирования применялись различные расчётные сетки. Для моделирования разлёта фрагментов ледяного метеорита расчётная область и схема разностной сетки приведены на Рис. 2.
Использовалась подвижная неравномерная вычислительная сетка со сгущением в области больших градиентов искомых функций. Передняя и задняя границы расчётной области перемещались с одинаковыми скоростями так, чтобы головная ударная волна не выходила за переднюю границу. Нижняя граница расчётной области соответствует оси х, являющейся осью симметрии течения, которая движется вдоль траектории полёта. Продольный и радиальный размеры расчётной области определялись на тестовых расчётах.
Для моделирования движения, газификации частиц и горения пылегазовой смеси за ударной волной в каналах и трубах применялись неравномерные подвижные сетки. Головная ударная волна выделялась правой границей расчётной области с использованием алгоритма предложенного С.К. Годуновым [82]. В качестве нижней границы расчётной области бралась ось симметрии для осесимметричных расчётов или нижняя стенка канала для плоского несимметричного случая. Верхняя граница расчётной области совпадала со стенкой трубы или верхней стенкой канала. В качестве левой границы расчётной области выбиралась стенка детонационной камеры. Рис. 3 демонстрирует расчётные области и схемы разностных сеток. За фронтом головной волны в зоне длины L применялась равномерная вдоль оси х сетка, далее размер ячеек увеличивался по логарифмическому закону вплоть до левой границы расчётной области. Вдоль оси у (или г для осесимметричных случаев) применялся алгоритм изложенный выше. В качестве зоны с равномерной по вертикали сеткой брался пылевой слой (или слои в несимметричном случае).
Расчётная область и схема разностной сетки для моделирования работы электрохимического пульсирующего двигателя представлены на Рис. 4.
Основной метод распараллеливания, используемый для задач математической физики и, в частности, численной газодинамики декомпозиция расчетной области на подобласти, обрабатываемые различными процессорами. Такой метод удобен при распараллеливании программ для компьютеров с распределенной памятью, поскольку каждый процессор может производить вычисления в своей подобласти, находящейся в своей локальной памяти. В разработанном вычислительном комплексе применён метод одномерной декомпозиции вычислительной области.
Расчёт конвективного переноса для параметров фаз
Движение космического тела (КТ) в атмосфере сопровождается его интенсивным разрушением под действием поверхностных аэродинамических сил и тепловых потоков, а также массовых сил инерции. Основными механизмами разрушения являются [2,3,5,6,7,8]: испарение, сублимация, плавление и термомеханическое шелушение поверхности тела, механическое дробление объёма тела на ряд отдельных, мелких фрагментов и их горение. Относительное значение различных механизмов разрушения тела изменяется вдоль его траектории и зависит от характера взаимодействия между КТ и атмосферой. При этом в верхних слоях атмосферы доминирует разрушение поверхности тела под действием тепловых потоков, а в нижних слоях механическое дробление частей тела под действием сил инерции и высоких давлений. Напряжённо-деформированное состояние предполагается квазистатическим и определяется путём численного или аналитического решения уравнений термоупругости. Для расчёта применяются как численный метод конечных элементов [3], так и аналитические способы, основанные на представлении решения в виде рядов и конечных сумм [5,9]. Для решения краевой задачи упругости с заданными поверхностными и массовыми силами разработан аналитический подход. Используется представление решения через ряды по сферическим функциям для сферических тел и через полиномы для круговых цилиндров и других видов тел вращения. Аналитические методы позволяют проводить многопараметрические расчеты, используя типовые компьютерные программы. Были выполнены расчеты многих вариантов входа КТ. Скорость тела и его абляция в верхних слоях атмосферы определяются путём решения системы уравнений физической теории метеоров. Ускорение центра масс тела и соответствующая сила инерции также вычисляются путем решения системы уравнений физической теории метеоров с учётом абляции. Предполагается, что КТ при его полёте до разрушения имеет форму тела вращения.
Дифференциальные уравнения движения КТ берутся в виде: (3.1.1) m— = -l-CDpV2S + mgsm0 , = - pV3S dt 2 dt 1 ..d9 mV2 . a dH .. . a mV — = cosp + mexoso = —Ksinp dt H + Rp dt
Здесь Со-коэффициент аэродинамического сопротивления, 5-площадь поперечного сечения тела, Rp - радиус планеты, р - плотность атмосферы, t -время, V- скорость, т - масса тела, Н - высота полета, в - угол наклона траектории полета, g - ускорение силы тяжести, аа - коэффициент абляции. Начальное значение параметров, при t=0, отметим индексом е. Начальные положения соответствуют высотам, где достаточно малы числа Кнудсена (Кп 0.5). Плотность атмосферы зависит от высоты и определяется из условия изотермичности по известным соотношениям [5]. Задача Коши для системы обыкновенных дифференциальных уравнений решалась численно, методом Рунге-Кутта. Скорости входа брались из диапазона 12 км/с V 70 км/с. При снижении КТ на него действуют аэродинамические силы, силы инерции и гравитации. Для изучения разрушения была применена модель упругого твердого тела, находящегося под воздействием аэродинамических сил, сил инерции и гравитации.
Мы рассматриваем задачу как квазистатическую и пренебрегаем производной по времени в известном уравнении Ляме теории упругости для вектора смещения и: (l-2v)Au + Vdivu=e3vCDpV +pg (3.1.2) где є - единичный вектор вдоль траектории (1,0,0), Z,v -параметры Ляме и Пуассона, R - радиус цилиндрического тела, правая часть уравнения (3.1.2) соответствует силам инерции и гравитации.
Решается упругая задача с заданными поверхностными силами, включая силы трения и нормального давления газа, обтекающего тело. Напряженное состояние тела находится путём представления решения через полиномиальные функции от пространственных координат. Коэффициенты полиномов подбираются с учетом точного удовлетворения уравнениям упругости (ЗЛ.2) и интегральному балансу действующих сил, а также приближенного удовлетворения условиям на поверхности тела [9]. Предполагается, что течение осесимметрично и вводится цилиндрическая система координат вдоль траектории полета. Тело считается однородным и изотропным. При движении вдоль траектории тело может разрушиться. Использовалось несколько критериев локального разрушения (по касательным напряжениям, по критерию Мора и другим [108]). Тело считалось разрушенным, если зона разрушения занимает значительную часть его объема (более 80%). Разработанный метод был применен к моделированию полёта и разрушения ледяных тел, каменных и металлических КТ.
Входными данными задачи являются: высота начала расчёта, масса, размеры тела, плотность тела, скорость и угол входа в атмосферу, материал тела, его прочностные свойства, критерии разрушения. Некоторые результаты расчётов представлены на Рис. 7 и Рис. 8.
Моделирование движения, газификации частиц и горения пылегазовои смеси за ударной волной
Ниже приведены результаты расчётов крупномасштабного взрыва в шахте. Высота тоннеля шахты Н=2,5м, длина тоннеля (галереи) равна 400м. (параметры экспериментальной шахты «Барбара» в Польше). У закрытого конца тоннеля имеется детонационная камера длины Д =6,6м, предназначенная для генерации ударной волны. Интенсивность её воздействия на процесс при прочих равных условиях определяется составом горючей газовой смеси, давлением и объёмом камеры. Толщина верхнего и нижнего слоев пыли hi и h2 бралась равной 1 см. Диаметр частиц пыли d = 60 мкм, плотность пыли в слое 400 кг/м . Предполагается, что к моменту времени t-О в детонационной камере, заполненной стехиометрической метано-воздушной смесью прошла детонационная волна, инициированная у закрытого конца трубы. Эта волна рассматривается как поверхность сильного разрыва с тепловыделением. Начальные распределения параметров продуктов детонации в камере определяются по автомодельному решению Тейлора о плоской одномерной детонации Чепмена-Жуге и задаются аналитическими зависимостями:
Здесь Р - давление, и - скорость, р - плотность, а — скорость звука, Dj -скорость детонационной волны, Q - удельное тепловыделение. Индексы 0 и J относятся к параметрам перед фронтом и непосредственно на фронте волны детонации, соответственно.
Кинетика горения угольных частиц и летучих, использованная в расчётах, описана в 1.6.1.
На Рис. 29 представлено изменение давлений со временем в середине канала (у= 1.25м) в точках вдоль тоннеля с координатами х=20, 40, 60, 100, 120, 160 м (начало взрыва), длина слоев пыли равна 110 м. На Рис. 30 распределения давления вдоль всего тоннеля даны для более поздних моментов времени, для положений ударной волны, соответствующих расстояниям 100, 200, 300 и 400 м. Верхний рисунок соответствует длине слоев пыли 400 м, нижний 110 м.
Ниже приведены некоторые результаты расчётов в круглых трубах различного диаметра. Рассматривается возникновение и распространение ударной волны в круглой трубе (двумерный осесимметричный случай), на внутренней поверхности которой находится тонкий слой угольной пыли. Геометрия трубы и расположение слоя пыли показаны схематически на Рис. 31. Пыль находится сначала только в слое вдоль трубы.
Толщина слоя угольной пыли h принималась равной 2 мм, диаметр частиц -96 мкм, процент летучих в частицах - 55. На Рис. 32 и Рис. 33 даны графики, иллюстрирующие процесс формирования пульсирующей детонации в трубе с радиусом 0.065м при начальном давлении в детонационной камере Рд=1 атм и её длине Lch=\3 м. Плотность угольной пыли в слое составляла 8 кг/м3.
На этих рисунках хорошо видно, что воспламенение пылегазовои смеси, приводящее к ускорению головной ударной волны и переходу к пульсирующему детонационному режиму горения, происходит на значительном расстоянии от закрытого конца трубы и реализуется в несколько этапов.
Механизм перехода к высокоскоростному горению наглядно демонстрирует Рис. 34, на котором приведены поля температуры газа в осевом сечении трубы для нескольких последовательных моментов времени, когда головная ударная волна находится на расстоянии 116, 119, 122, 125 и 128 м от закрытого конца трубы. Из рисунков следует, что переход к детонации определяется ускорением фронта горения пылегазовой смеси и его распространением в значительной зоне за фронтом головной ударной волны. В результате чего интенсивность головной ударной волны возрастает, и создаются условия для вторичного воспламенения пылегазовой смеси в непосредственной близости от её Наблюдающиеся на Рис. 356 пики давления сразу за фронтом головной ударной волны были обнаружены ранее в экспериментах [16]. На Рис. 36 приведена скорость головной ударной волны в зависимости от расстояния от закрытого конца трубы. Видно, что здесь, в целом, процесс подобен случаю узкой трубы - наблюдаются этапы затухания и последующего ускорения головной ударной волны, но менее интенсивен. Последнее обстоятельство связано, в первую очередь, с большим диаметром трубы, при котором имеет место усиление влияния поперечных волн на процесс в целом. Другими словами, в данной ситуации существенны двумерные эффекты.
На Рис. 37 приведено поле температуры газа, а на Рис. 38 поле десятичного логарифма от плотности частиц в осевом сечении трубы для нескольких последовательных моментов времени, когда головная ударная волна находится на расстоянии 50, 100, 150, 200 и 250 м от закрытого конца трубы. На Рис. 37 видно, что фронт пламени имеет сложную зигзагообразную форму, которую он приобретает в результате взаимодействия с мощными поперечными волнами. Рис. 38 показывает, что поднимаемый головной ударной волной слой пыли приводит к образованию достаточно плотной газопылевой смеси со слоистой структурой, коррелирующей с формой фронта пламени.
Результаты расчётов хорошо согласуются с экспериментами в шахтах и ударных трубах [16,17].
Предложена математическая модель, разработан численный метод и проведены расчёты течений пылегазовых смесей, формирующихся за ударной волной в результате ее взаимодействия с узким пристеночным слоем пыли. Исследовано влияние коэффициентов в силах Магнуса и Сэффмана на динамику подъёма и диспергирования пыли из плотного слоя. Получена полная картина течения в длинных трубах большого и малого диаметра с запылёнными внутренними поверхностями. Выявлены механизмы интенсификации горения и формирования детонационных и квази-детонационных волн. Установлено, что детонация может возникать на больших расстояниях от места формирования первичной ударной волны, и поэтому при исследовании процессов в реальных экспериментах для получения достоверной информации, с точки зрения взрывобезопасности, необходимо использовать трубы и галереи с длиной, на порядок превосходящей используемые в настоящее время.
