Содержание к диссертации
Введение
Глава 1. Основной объект. Исходная система уравнений и ее расщепление 10
1.1. Основные уравнения 11
1.2. Приведение уравнений и расщепление полученной системы 14
1.3. Замечания о движении сред, находящихся в неравновесном состоянии 18
Выводы к первой главе 20
Глава 2. Чередование явных и неявных разностных схем в процессе интегрирования уравнений математической физики 22
2.1. Аналитическое исследование устойчивости разностных схем 22
2.2. Исследование устойчивости разностных схем в численных экспериментах 27
2.3. Выводы к второй главе 28
Глава 3. Математическое моделирование кинетических процессов. О моделировании кинетических процессов 30
3.1. Постановка задачи 30
3.1.1. Постановка прямой кинетической задачи 30
3.1.2. Основные понятие и величины химической кинетики 32
3.1.3.Математическая модель 34
3.2. Численные методы решения задачи 37
3.2.1. Трудности, возникающие при численном решении прямой кинетической задачи 37
3.2.2. Метод Гира 39
3.2.3. Метод Рожкова 41
3.2.4. Численное решение задачи с помощью метода Рожкова 43
3.2.5. Полученные результаты 47
3.3. Асинхронный метод 51
3.3.1. Описание метода 51
3.3.2. Применение на модельной задаче 55
3.3.3. Численное решение задачи с помощью асинхронного метода 63
3.3.4. Полученные результаты 68
3.4. Использование многопроцессорной вычислительной техники. 71
3.4.1. Параллельные вычисления 71
3.4.2.Архитектура параллельной вычислительной техники 72
3.4.3Лополигии на множестве процессоров 73
3.4.4. Виды параллелизма 73
3.4.5. Разработка параллельного приложения 74
3.4.6. Эффективность вычислений на многопроцессорной вычислительной технике 75
3.4.7. Распараллеливание вычислительного процесса 75
Заключение по третьей главе 77
Приложение к третьей главе 78
Список литературы к третьей главе 104
Глава 4. Представление о вычисляющей среде и его применение для распараллеливания алгоритмов в механике жидкостей и газов 105
4.1 Вычисляющая среда. Описывающие ее уравнения. 105
4.2. Оптимальное распараллеливание 107
Список литературы к четвертой главе 108
Глава 5. Математическое моделирование течения реагирующих сред. Оптимальное распараллеливание 109
5.1. О моделировании процессов в реагирующих потоках 109
5.2. Постановка задачи и исходные уравнения 110
5.2.1. Постановка задачи 110
5.2.2. Начальные и граничные условия 114
5.2.3. Исходные уравнения . 116
5.3. Алгоритм решения 120
5.3.1. Решение общего уравнения 120
5.3.2. Решение уравнения для определения поля давления 132
5.3.3. Функция для вычисления значения турбулентной вязкости 133
5.3.4. Некоторые примеры 134
5.4. Алгоритм распараллеливания 142
5.4.1. Простое распараллеливание 142
5.4.2. Оптимальное распараллеливание 144
5.4.3. Оценка полученных результатов 150
Заключение по пятой главе 154
Список литературы к пятой главе 155
Результаты исследования 156
- Замечания о движении сред, находящихся в неравновесном состоянии
- Исследование устойчивости разностных схем в численных экспериментах
- Трудности, возникающие при численном решении прямой кинетической задачи
- Функция для вычисления значения турбулентной вязкости
Введение к работе
Математическое моделирование течений, осложненных физико-химическими процессами, актуально в следующих областях.
1. Определение условий сжигания топлива в энергетике и на транспорте, при которых достигается максимальная энергетическая эффективность и возникает минимальное количество экологически вредных выбросов. Разработка устройств, в которых реализуются указанные режимы горения.
2. Создание конструкций, топлив и режимов работы аэродинамических двигателей, обладающих высокими силовыми характеристиками, и минимальным воздействием на атмосферу.
3. Моделирование процессов в химических и газодинамических лазерах.
4. Оптимизация химических реакторов, использующих традиционные технологии и создание новых на основе неравновесной плазмохимии.
5. Исследование экологических процессов в реках, озерах и водохранилищах.
Ограничимся перечисленными областями, хотя указанный список легко может быть расширен.
Необходимо сразу же указать, что численное моделирование процессов в указанных областях связано с очень большими вычислительными затратами.
Положение еще усугубляется, если нужно не только промоделировать работу некоторого устройства, но и найти наилучшее инженерное решение, т.е. организовать поиск по некоторым параметрам, связанный с перебором вариантов.
Указанные вычислительные трудности являются существенным препятствием на пути создания совершенных устройств, отличающихся высокой продуктивностью, низкими энергозатратами и экологической безопасностью.
Замечания о движении сред, находящихся в неравновесном состоянии
Решение задач численного моделирования химических процессов имеет высокую практическую ценность и естественнонаучный интерес. С помощью применения вычислительных систем может быть проведено численное моделирование таких актуальных задач, как моделирование процессов горения и многих других технологических процессов.
Однако при изучении и моделировании уравнений химической кинетики возникает ряд математических трудностей из-за одновременного присутствия очень медленно и очень быстро протекающих химических реакций. При расчете этих задач многие вычислительные методы дают сбой. Для этих систем разрабатываются специальные и эффективные алгоритмы, к которым в последнее время усиливается интерес.
Моделирование кинетики химических процессов рассматривается на реакции окисления метана. Выбор такого объекта исследования обусловлен тем, что в последнее время перед человечеством все острее встают вопросы загрязнения окружающей среды и истощения запасов органического топлива, прежде всего нефти. Одним из путей решения этих вопросов является использование экологически чистого, возобновляемого в приемлемые сроки вида топлива. Таковым может быть метанол, получаемый при конверсии метана.
Рассмотрим реакцию окисления метана. Механизм химического процесса представленного в табл.1.1., основан на предложенной Н.Н. Семеновым радикально-цепной схеме. Реакция (1) есть реакция зарождения цепи, реакция (2)-(5) и (7)-(12) - продолжение цепи, реакция (6) - разветвление цепи, реакции (13)-(17) - обрыва цепи [3.1]. В этой же таблице приведены значения соответствующих предэкспонент и энергий активации.
Под прямой задачей химической кинетики понимают задачу нахождения концентраций участвующих в реакции веществ в любой момент времени, исходя из известных начальных концентраций, схемы реакции и констант скоростей отдельных стадий. Необходимо решить прямую кинетическую задачу для реакции окисления метана, схема которой представлена в табл.3.1.1., при заданном соотношении веществ метана и кислорода [СУЦ]: [02] - 0.29 : 0.71. Фундаментальное понятие химической кинетики - механизм реакции. В широком смысле слова механизм - это содержательная интерпретация опытных данных, являющаяся квинтэссенцией всех накопленных сведений о протекании сложной реакции. В таком механизме должны быть вычислены отдельные этапы и стадии реакции, приведены характеристики промежуточных продуктов, описаны переходные состояния отдельных стадий, даны энергетические уровни веществ и т.д. Применительно к каталитическим реакциям должны быть охарактеризованы свойства поверхности, исследован характер адсорбции и т.п. В практике химической кинетики употребительно и другое (узкое) понимание механизма как совокупность стадий. Каждая стадия состоит из прямой и обратной реакций. Если стадии предполагаются простыми, то они состоят из элементарных реакций. В качестве кинетического закона которых, принимается закон действующих масс (или действующих поверхностей - для гетерогенных каталитических реакций) [3.6]. Элементарная реакция - это реакция, осуществляющаяся с преодолением одного энергетического барьера. Один барьер - одна элементарная реакция. Именно таким понятием механизма оперирует формальная кинетика, исследующая кинетические модели - системы дифференциальных и алгебраических уравнений, соответствующих механизму. Химическая кинетика описывает стадии химические реакции с помощью так называемых стехиометрическихуравнений: где 1} и - это стехиометрические коэффициенты исходных веществ Л, и продуктов реакции В І соответственно; (стехиометрические коэффициенты - это коэффициенты перед молекулами веществ в химическом уравнении); Л А» NB - число исходных веществ и продуктов реакции соответственно. Это символическая запись химической реакции, в которой її молекул вещества Аи h молекул вещества А , ..., 1ш молекул вещества AN, вступив в реакцию, образуют г\ молекул вещества Ви Гг молекул вещества В2, rN молекул вещества BNb. [6] Основной закон, выражающий скорость протекания реакции - закон действующих масс гласит; скорость протекания реакции пропорциональна концентрациям реагентов. Поэтому, если обозначить буквами q, концентрации соответствующих веществ, то где К- константа скорости протекания реакции N - число веществ, участвующих в реакции. Константа скорости протекания реакции - основная кинетическая характеристика; коэффициент пропорциональности в кинетическом уравнении, связывающем скорость реакции с концентрациями реагирующих веществ и их стехиометрическими коэффициентами. Для мономолекулярных реакций константа скорости имеет размерность с"1, для бимолекулярных - л/моль-с. Скорость каждой химической реакции зависит как от природы реагирующих веществ, так и от условий, в которых реакция протекает. Важнейшими из этих условий являются: концентрация, температура и присутствие катализатора. С. Аррениус впервые показал, что влияние температуры сводится к увеличению числа активных молекул, т. е. таких молекул, которые в момент столкновения обладают энергией, не меньше определенной для данной реакции величины, называемой энергией активации химической реакции. Энергия активации - это некоторое избыточное количество энергии, необходимое для вступления молекул в реакцию [3.1].
Исследование устойчивости разностных схем в численных экспериментах
Существует особенность решения систем дифференциальных уравнений, описывающих химические реакции. Она состоит в том, что временные характеристики различных переменных существенно отличаются друг от друга. Для кинетических систем такого типа характерно наличие быстро и медленно меняющихся переменных. Быстрые переменные практически мгновенно подстраиваются под изменения медленных.
Наличие быстрой и медленной подсистем определяет трудности, возникающие при численном решении прямой кинетической задачи. Для достаточно точного вычисления решения по быстрым переменным необходимо выбирать шаг интегрирования значительно меньшим, чем полное время протекания процесса, определяемое изменением медленных переменных.
Системы дифференциальных уравнений, описывающие поведение как быстрой, так и медленной подсистемы, называются жесткими. Основные сложности, возникающие при численном решении прямой кинетической задачи, связаны именно со свойством жесткости задачи Коши: У = /(у), t&[a,b], у(а) = у0, В случае прямой кинетической задачи вектор / всегда удовлетворяет условию Липшица по всем п компонентам вектора у, поэтому можно доказать, что решение задачи Коши существует и оно единственно. Численное решение задачи Коши - это вычисление последовательности приближений ys = у (ts) в точках дискретизации { /s}, где hs - шаг дискретной сетки. Метод, по которому получается численное решение, должен обладать свойствами, обеспечивающими приближение точного решения численным в точках дискретизации. Естественным требованием, предъявляемым к методу, является возможность оценки разности (К , )-5 1 во всех точках дискретизации. При уменьшении шага сетки, на которой ищется численное решение, эта разница должна уменьшаться. Такое свойство численного метода называется сходимостью. Говорят, что метод обладает сходимостью, если при стремлении к нулю шага дискретизации к нулю стремится и разность точного и численного решений во всех точках разбиения. Второе необходимое свойство численного метода - это свойство устойчивости. Оно состоит в том, что малые возмущения в начальных условиях и в разностном уравнении, возникающие, например, в результате округления и конечной длины разрядной сетки машины, подавляются или по крайне мере не увеличиваются. Один и тот же метод при разных шагах дискретизации может быть устойчивым, так и неустойчивым. Однако устойчивость метода связана не только с шагом численного интегрирования, но и со свойствами конкретной задачи Коши. Поэтому, для того чтобы охарактеризовать устойчивость метода численного интегрирования вне зависимости от конкретной задачи, его, как правило, исследуют на устойчивость при решении модельной линейной задачи Коши. у = Лу, у(0) = yQ, где X - комплексная постоянная. Это естественно, так как локальное поведение исходной задачи в первом приближении определяется решением линейной задачи с матрицей, являющейся якобианом устойчивости метода. Метод называется абсолютно устойчивым, когда для фиксированного шага интегрирования полная погрешность метода \ys y(ts)\ остается ограниченной при S - . Говорят, что задача Коши является жесткой, если в локальной области она устойчива, т.е. собственные значения якобиана имеют отрицательные действительные части: Re(/L(.) Ot и R -A /Rei-A )»!. В задачах химической кинетики локальная жесткость может достигать порядка 10+б - 10+9. Трудности решения жестких задач может накладывать требование абсолютной устойчивости метода, связанное с малыми возмущениями, возникающими в процессе реализации метода на цифровой машине. Действительно, величина шага интегрирования должна выбираться так, чтобы AaMh принадлежало области абсолютной устойчивости метода. Таким образом, шаг интегрирования согласуется с характерным временем быстрого процесса ReC-Л ), в то время как характерное время медленного процесса Ret-Д ) много больше, и необходимое число шагов интегрирования будет сравнимо с Для того чтобы избежать этого ограничения, и, чтобы выбор шага определялся лишь соображениями точности решения, при решении задач Коши используют специальные жестко устойчивые численные методы, такие как метода Гира, возможно использование метода Розенброка. Жестко устойчивый метод Гира является представителем семейства методов типа «предикатор-корректор», в которых на первом шаге вычисляется (на основе экстраполяции) предварительное решение, которое уточняется на втором шаге. В ходе вычислений постоянно проводится внутренний контроль ошибки интегрирования, на основании которого могут изменяться порядок метода и шаг интегрирования по времени. Описание метода Гира очень громоздко. Поэтому рассмотрим его лишь в общих чертах. Сначала вычисляется предикатор С - приближенное значение концентраций для нового момента времени п. Вычисление осуществляется в соответствии с алгоритмом Нордсика, с помощью интерполяционного полинома, использующего значения вектора концентраций Cn.j и его производных С } в предыдущий момент времени п-1:
Трудности, возникающие при численном решении прямой кинетической задачи
Для написания корректного и эффективно работающего параллельного приложения необходимо соблюдения следующих условий: Возможность алгоритма разбиваться на относительно независимые, требующие примерно одинакового времени выполнения блоки которые впоследствии реализуются как отдельные процессы. Параллельно работающие части алгоритма должны составлять большую часть времени выполнения части программы. При уменьшении параллельной части приложения резко снижается эффективность. Количество и объем передаваемых сообщений должны быть сведено до минимума; кроме того, допустимая частота сообщений может лимитироваться используемыми коммуникациями. Применим механизм распараллеливания к решению задачи при асинхронном подходе (см, приложение листинг № 3.4). і, При решении систем обыкновенных дифференциальных уравнений чаще используется функциональный вид параллелизма. Воспользуемся этим видом при решении данной задачи. Возьмем п = 12 процессоров и соединим их в топологию звезда. Алгоритм распараллеливания: Первый такт по времени нулевой процессор считывает информацию из файла о количестве веществ и о количестве реакций, рассылает информацию всем рабочим процессорам. Рабочие процессоры принимают эти значения, и каждый определяет количество реакций, для которых будет вычислять значения констант скорости. Второй такт по времени нулевой процессор считывает информацию о начальной концентрации веществ, значениях предэкспонент и энергии активации, рассылает каждому процессору нужное количество этих значений. Рабочие процессоры принимают значения и вычисляют константы скорости реакций. Третий такт по времени рабочие процессоры отсылают вычисленные значения, которые принимает нулевой процессор и формирует из них массив значений констант скорости реакций. Четвертый такт по времени нулевой процессор рассылает сформированный массив на предыдущем этапе всем рабочим процессорам, которые его принимают. Пятый такт по времени нулевой процессор формирует массив Name и рассылает его всем рабочим процессорам, которые его принимают. Одновременно все рабочие процессоры транслируют реакции и заполняют массив Subst. Шестой такт по времени нулевой процессор определяет количество «медленных» веществ, которое должно быть на каждом процессоре и пересылает номера этих веществ каждому процессору, которые принимают данные значения и формируют массив с номерами, соответственно «медленных» веществ. Организуем цикл по времени от TimeBegin до TimeEnd с различными шагами для каждой подсистемы. Седьмой такт по времени нулевой процессор рассылает всем рабочим процессорам начальные значения концентрации веществ, которые его принимают. Восьмой такт по времени нулевой процессор вычисляет значения «быстрых» веществ с маленьким шагом и пересылает массив значений Subst[i\.Concentr на рабочие процессоры, которые его получают. Девятый такт по времени все рабочие процессоры вычисляют свои значения «медленных» веществ и отсылают на нулевой процессор, который их принимает, записывает в массив. Десятый шаг по времени все процессоры делают шаг по времени. Одиннадцатый шаг по времени проверяют условие окончания вычислительного процесса, если TimeBegin TimeEnd, то переходим к восьмому такту по времени, иначе нулевой процессор записывает полученные значения в файл. В данной работе реализованы основные этапы вычислительного эксперимента, который заключался в нахождение концентраций участвующих в реакции веществ в любой момент времени, исходя из начальных концентраций, схемы реакции и констант скоростей отдельных стадий. Сложная реакция окисления метана была исследуемым объектом.
В результате выполнения работы построена математическая модель химического процесса, которая представляла собой жесткую систему дифференциальных уравнений. На этапе построения вычислительного алгоритма - метода приближенного решения поставленной задачи, проведено исследование двух методов - синхронного и асинхронного для определения наиболее эффективного из них. Это сравнение сначала производилось на модельной задаче, а затем для задачи о конверсии метана.
Разработаны программы, реализующие данные методы, проведены различные эксперименты. Результаты исследования, представленные в графическом виде, отражают зависимость концентраций компонент системы от времени.
Основным итогом работы стало доказательство эффективного использование асинхронного метода для решения жестких систем дифференциальных уравнений, описывающих химическую кинетику.
Функция для вычисления значения турбулентной вязкости
Эта глава посвящена вопросам оптимизации параллельных программ для моделирования физико-химических процессов. Моделирование физико-химических процессов имеет большое прикладное значение. Оно используется в фундаментальной науке и различных областях производства, а именно: в энергетике, в ракетостроении, в экологии и т.д. Однако моделирование таких процессов связано с выполнением большого объема вычислений. Для их реализации приходится применять мощные многопроцессорные вычислительные системы. При этом возникает проблема написания параллельных программ для моделирования таких процессов.
При использовании многопроцессорной вычислительной техники расчетные поля (поля, которые нужно найти при моделировании процессов, например, поля скорости, температуры и др.) обычно разделяются между процессорами на участки согласно геометрическому или конвейерному видам параллелизма. Каждый процессор вычисляет значения полей на своих участках в следующий момент времени и передает новые вычисленные значения на своих границах «соседним» процессорам. Целью исследования являлся поиск возможных путей оптимизации таких алгоритмов. Под словом оптимизация понимается уменьшение времени выполнения программы и увеличение точности получаемых результатов. В данной работе предложены некоторые возможные способы оптимизации распараллеливания.
Эти способы оптимизации были применены к задаче о распространении реагирующих сред в сосуде. Для моделирования этих процессов на компьютере был написан программный комплекс на языке Microsoft Visual C++ 6.0 с использованием интерфейса Message Passing Interface и его реализации WMPI Beta. Программный комплекс написан для платформы Microsoft NT. В ходе изложения приводятся листинги отдельных частей программ. Содержание главы следующее: В 5.2 ставится задача о распространении реагирующих сред в сосуде: т- приводятся исходные данные и граничные условия, а также приводятся исходные уравнения для решения поставленной задачи. В 5.3 рассматривается алгоритм решения поставленной задачи на компьютере: приводятся разностные уравнения и код подпрограммы для их решения в однопроцессорном варианте. Также приводятся результаты работы программы. В 5.4 рассматриваются различные способы распараллеливания: простое распараллеливание и два варианта оптимального распараллеливания. Сравниваются результаты, полученные при применении каждого из этих способов. Для исследования поставленной проблемы, были решены три задачи о моделировании процессов для движущихся и реагирующих сред. Формулировка задач следующая. Дан сосуд. В первой и во второй задачах форма сосуда изображена на рис. 5.2.1 и 5.2.2 соответственно. В сосуд ведут две трубы, а из сосуда выходит одна труба. В первую трубу поступает водяной пар, насыщенный водородом, во вторую - водяной пар, насыщенный кислородом. В сосуде газы перемешиваются, кислород и водород реагируют друг с другом с выделением тепла. В результате реакции образуется водяной пар. Из третьей трубы выходит водяной пар с остатками водорода и кислорода. В третьем случае (рис. 5.2.3) размер сосуда по оси Оу намного больше размера сосуда по оси Ох. Для упрощения задачи он принимается равным бесконечности. В него ведет множество труб (расположены слева), поступающие вещества, в которых чередуются: в одну трубу поступает водяной пар, насыщенный водородом, в другую - водяной пар, насыщенный кислородом, и т.д. Из сосуда ведет одна труба, которая расположена справа. Далее рассматривается часть сосуда (рис.5.2.4), охватывающую половины двух труб, в которую поступают вещества и часть трубы, расположенную напротив них, из которой выходит получившаяся смесь. В последнем случае учитывается тот факт, что сосуд состоит из подобных частей, границы которых (параллельны оси Ох) являются осями симметрии сосуда. В задаче требуется рассчитать поля скорости, давления, температуры и концентраций веществ в сосуде в последующие моменты времени. Начальные скорости газов, поступающих через трубы, изменяются по параболическому закону (рис. 5.2.1), то есть в центре трубы газ имеет максимальную скорость, а по краям трубы скорости равны 0. Исключение составляет третья труба в третьем случае. Поскольку она бесконечна по ширине, скорость газа в каждой ее точке постоянна. Средняя скорость газов в первой трубе взята в 2 раза меньше, чем средняя скорость во второй трубе, а средняя скорость газов во всем сосуде взята равной 1 міс. Начальное давление принято равным 100 атм. Температура стенки сосуда и поступающих газов взята равной 1000 К, за исключением температуры сред, поступающих во вторую трубу в первой задаче. Там она взята равной 1100 К. Концентрация поступающих сред вычисляется согласно закону Авогадро по формуле где Р -давление, Т температура, R - молярная газовая постоянная.
Концентрация водорода в поступающей через первую трубу смеси газов, и концентрация кислорода в поступающей через вторую трубу смеси газов принимаются равными 0.01 от концентрации этих смесей газов. Поэтому в дальнейшем в расчетах часто будут браться характеристики водяного пара, а не всей смеси газов, ввиду малости их концентрации.
Имеются 4 типа границы в сосуде: «вход» (на границе с I и II трубами), «выход» (на границе с III трубой), «граница» (оставшаяся часть границы в первой и во второй задачах) и «симметричная граница» (оставшаяся часть границы в третьей задаче). На границе типа «граница» скорость равна нулю, то есть используется граничное условие первого рода. На границах типа «вход» и «выход» также используются граничные условия первого рода; там скорость задана.
