Содержание к диссертации
Введение
ГЛАВА 1. Обзор публикаций по численному и математическому моделированию природных катастроф 14
1.1 Современные программные комплексы, моделирующие горение 14
1.2 Программные комплексы, моделирующие цунами 17
1.3 Современные технологии распараллеливания 21
ГЛАВА 2. Методы дискретизации системы нелинейных и взаимосвязанных дифференциальных уравнений 31
2.1 Метод контрольного объема 31
2.1.1 Получение дискретного аналога для одномерного случая 31
2.1.2 Получение дискретного аналога для трехмерного случая 39
2.1.3 Об оптимизации и распараллеливании метода Патанкара 45
2.2 Итерационно-интерполяционный метод 71
2.3 Дискретные модели математической физики 75
2.3.1 Разностные схемы для уравнений гиперболического типа 75
2.3.2 Разностные схемы для уравнений параболического типа 85
2.4 Применение итерационно-интерполяционного метода к уравнениям гиперболического типа 92
2.4.1 Линейные уравнения мелкой воды 92
2.4.2 Исследование на устойчивость 98
2.4.3 Система нелинейных уравнений мелкой воды 99
2.4.4 Нелинейные уравнения мелкой воды в потоках 99
2.5 Численное моделирование жестких систем обыкновенных дифференциальных уравнении 101
2.5.1 Модификация схемы Кранка-Николсона 101
2.5.2 Решение уравнений химической кинетики 102
2.5.3 О сходимости численных методов интегрирования жестких систем обыкновенных дифференциальных уравнений и выборе шага численного интегрирования 106
ГЛАВА 3. О программном комплексе моделирования цунами 108
3.1 Физическая модель 108
3.2 Математическая модель, используемая для расчета волн цунами, вызванных смещением подводных «блоков-клавиш» 108
3.3 Результаты численного моделирования 110
3.3.1 Тестирование разностной схемы, полученной на основе итерационно- интерполяционного метода
3.3.2 Пример численного моделирования волны цунами
ГЛАВА 4. О программном комплексе моделировании горения 117
4.1 Физическая модель 117
4.2 Математическая модель 119
4.2.1 Система уравнений пограничного слоя 119
4.2.2 Начальные и граничные условия для пограничного слоя 121
4.2.3 Система уравнений тепло- и массопереноса для слоя 122
4.2.4 Начальные и граничные условия для слоя 125
4.3 Моделирование взаимодействия пеллеты с горячей средой 126
4.3.1 Физическая модель зажигания пеллеты 126
4.3.2 Математическая модель зажигания пеллет 128
4.3.3 Определение прогоночных коэффициентов 131
4.3.4 Результаты численного моделирования зажигания пеллет 132
4.3.5 Об определении времени зажигания пеллеты 137
Заключение 140
Литература 142
- Программные комплексы, моделирующие цунами
- Итерационно-интерполяционный метод
- Нелинейные уравнения мелкой воды в потоках
- Тестирование разностной схемы, полученной на основе итерационно- интерполяционного метода
Введение к работе
Актуальность темы и востребованность результатов работы связаны с тем, что в последние годы существенно возросло число природных бедствий, связанных с разрушительными свойствами цунами и пожаров. Обилие лесов и непрекращающееся освоение человеком все новых территорий постоянно повышает риск возникновения пожаров, а постройка все новых объектов промышленности и жизнедеятельности у побережья морей и океанов увеличивает риск катастрофического воздействия волн цунами.
Анализ статистических данных о характере природных катастроф за последние десятилетия беспристрастно демонстрирует тенденцию к заметному росту. И в первую очередь это утверждение относится к пожарам и волнам цунами.
Поскольку природные катастрофы относятся к неотвратимым явлениям, на первый план выходит задача минимизации нанесенного ущерба. Значимым вкладом на пути решения обозначенной проблемы является развитие экспертных систем предупреждения о цунами и распространении пожаров, учитывающих региональные особенности конкретной территории и предназначенных для информационной поддержки процедуры принятия решений в чрезвычайной ситуации во время надвигающейся катастрофы и/или в момент планирования хозяйственного освоения новых территорий, проектирования, возведения и эксплуатации новых промышленных и социальных объектов.
Рассматриваемые системы должны обладать возможностями моделирования различных сценариев развития природных катастроф с использованием реальной батиметрии и топографии, а также характеристик типа и влажности леса. Кроме того, они должны наиболее рационально расходовать ресурсы вычислительной системы, а в первую очередь - центрального процессора. Для этого разработанные алгоритмы должны быть оптимизированы и распараллелены.
С развитием информационных технологий появилась возможность создания экспертных систем подобного класса.
Цель работы:
Разработка трехмерной модели распространения горения с учетом многослойности (с заданием параметров среды для каждого слоя), двухтемпературности и нескольких очагов возникновения.
Разработка эффективных численных схем на основе итерационно-интерполяционного метода для решения уравнений мелкой воды.
Разработка, эффективных алгоритмов для решения трехмерного обобщенного уравнения теплопроводности с использованием современных вычислительных технологий.
Создание вычислительного и программного обеспечения нового поколения для решения актуальных исследовательских и прикладных задач проблем моделирования природных катастроф.
Применение созданного вычислительного инструмента для решения задачи зажигания и горения слоя пеллет.
Выполнен анализ имеющихся моделей и средств для интерпретации природных катастроф, сформулированы требования к соответствующим моделям и алгоритмам. Определены особенности разработки и реализации программных комплексов, необходимых для решения разных типов исследовательских и прикладных задач обозначенных проблем.
На основе итерационно-интерполяционного метода разработаны численные схемы для решения уравнений мелкой воды.
С использованием эффективных конечностно-разностных
Проведен анализ алгоритмических; низкоуровневых оптимизаций и эффективности распараллеливания алгоритмов.
При решении задачи о зажигании и горении; слоя:- пеллет проведено математическое и численное моделирование процессов сушки пеллет, зажигания и горения в. зависимости от типа пеллет и их характеристик.
Определено время зажигания и основные, характеристики процесса зажигания и горения.
Математических моделей волновой гидродинамики;
Математических моделей механики сплошных многофазных сред;
Эффективных вычислительных конечно-разностных алгоритмов;
Итерационно-интерполяционного метода;
Оптимизированного и расширенного до трехмерного случая метода Патанкара;
Алгоритмическая и низкоуровневая оптимизация разработанных программных алгоритмов, а также их распараллеливание;
Принципов и технологий создания проблемно-ориентированных программных комплексов, характеризующихся интегрированностью моделирующих, информационных и интерфейсных компонент, обеспечивающих, в свою очередь, возможность эксплуатации систем пользователями различного уровня квалификации.
Программный комплекс моделирования возможных цунами, вызванных землетрясением на основе клавишной модели.
Программный комплекс моделирования распространения горения в многослойной двухтемпературной среде.
Программные системы и отдельные компоненты для решения исследовательских и прикладных задач природных катастроф. А именно:
Результаты решения тестовых и модельных задач, в ходе которых были выявлены ключевые характеристики и диапазон применения созданных алгоритмов и программных систем. А именно:
Результаты решения ряда важных и актуальных прикладных задач. А именно:
Для достижения поставленной цели были решены следующие задачи:
алгоритмов и оптимизации метода Патанкара,
модернизированного под трехмерный случай, разработана
трехмерная вычислительная модель описания, распространения . пламени в многослойной двухтемпературной среде.
4. Для проверки качества вычислительных алгоритмов построен набор тестовых и модельных задач. Проведены соответствующие вычислительные эксперименты и.анализ их результатов.
.5. Для решения производственных, исследовательских и информационно-обучающих задач проблем цунами и горения созданы программные комплексы, основанные на разработанных . моделях.
9: Исследовано влияние температуры и. влажности» пеллет на их температуру горения и время зажигания.
Научная новизна работы; Впервые разработана полная трехмерная модель распространения фронта пламени с. учетом двухтемпературности и нескольких очагов; возникновения на основе расширенного и оптимизированного метода Патанкара для трехмерного случая. Впервые произведена систематизация и обобщение моделей различных типов пожаров И' создана единая модель, позволяющая^ задавать несколько слоев» с различными типами пожаров,; тогда как ранее для - каждого типа пожара строилась своя конкретная модель.
На основе современных вычислительных средств разработаны оригинальные инструменты проведения вычислительных экспериментов как для; решения прикладных . задач проблем цунами, так и .горения. Разработанные программные комплексы включают в себя системы
математических моделей, их алгоритмические реализации и необходимые базы данных. Интерфейс программных комплексов позволяет начать работу с ними даже неподготовленному пользователю.
Впервые при моделировании цунами введена возможность задания нескольких динамических источников, сбор данных со спутника, а также динамическое снятие мареограмм с мареографов в заданных пользователем точках на карте.
Особое внимание при разработке программной реализации моделей уделено внимание их алгоритмической и низкоуровневой оптимизации, что практически не встречается в подобных системах, а также проведена работа по распараллеливанию вычислений.
Проведена оценка эффективности пеллет на основе моделирования процесса горения.
Достоверность. Обоснованность и достоверность полученных результатов обеспечивается математическим обоснованием используемых« методов и алгоритмов и подтверждается согласованием результатов численных расчетов как с известными экспериментальными и натуральными данными, решением ряда тестовых и модельных задач, так и с данными, полученными при* помощи других программных систем.
Практическая значимость работы определяется успешным опытом использования ее основных результатов (программ, результатов расчетов) при решении актуальных прикладных задач природных катастроф, в том числе для создания информационно-моделирующих подсистем для автоматизации действий служб предупреждения природных катастроф.
Программные комплексы, созданные автором, успешно используются для решения задач проблем цунами и горения в рамках научной деятельности, выполнения производственных договоров и федеральных целевых программ в Нижегородском государственном техническом университете им. P.E. Алексеева. Результаты численных расчетов, полученные при помощи программного комплекса "Моделирование горения слоя пеллет" составляют основу моделирования горения слоя пеллет в ОАО
«Центральный научно-исследовательский и проектный институт лесохимической промышленности», что подтверждается актом о внедрении.
Методология исследования опирается на современные информационно-вычислительные технологии, предусматривающие использование
На защиту выносятся:
а) Оптимизированное ядро программного комплекса, реализующее функции модифицированного метода
Патанкара для трехмерного случая;
б) Ядро программного комплекса, реализующее итерационно- интерполяционный метод для решения нелинейных уравнений мелкой воды;
в) Программный комплекс по моделированию процесса зажигания и горения слоя пеллет.
а) Проведено распараллеливание программной реализации метода Патанкара для систем с общей памятью;
б) Получено ускорение оптимизированного метода Патанкара;
в) Выявлен эффект замедления программной реализации метода Патанкара при использовании элементов массивов различной длинны.
а) При решении задачи о зажигании и горении слоя пеллет проведено математическое и численное моделирование процессов сушки пеллет, зажигания и горения в зависимости от типа пеллет и их характеристик.
б) Определено время зажигания и основные характеристики процесса зажигания и горения.
в) Исследовано влияние температуры на время зажигания, а влажности и зольности — на теплотворную способность пеллет.
Общая характеристика работы.
Во введении обосновывается актуальность темы исследования,
формулируется цель и ставятся основные задачи работы. Приводится обзор
научной литературы по изучаемой проблеме. Раскрывается научная новизна
и практическая значимость работы. Формулируются основные положения, выносимые на защиту. Основное внимание уделяется численным методам, использующимся для решения соответствующих задач математического моделирования: конечно-разностные алгоритмы, итерационно- интерполяционный метод, метод контрольного объема; а также специальным вычислительным приемам, обеспечивающим быстроту и качество разрабатываемых алгоритмов.
В главе 1 изложен краткий обзор публикаций по проблемам, рассматриваемым в данной диссертации, приводится подробный анализ работ и программных комплексов, моделирующих цунами и пожары. Делается вывод об актуальности постановок и решения задач.
В параграфе 1.1 рассматриваются основные приемы и технологии распараллеливания. Делается вывод о том, что в связи со спецификой задачи целесообразно применять методы распараллеливания для вычислителей с общей памятью (многопроцессорным и/или многоядерным).
В параграфе 1.2 представлен обзор по математическому и численному моделированию процессов горения в открытых системах. Рассматриваются существующие программные комплексы.
В параграфе 1.3 представлен обзор по математическому и численному моделированию генерации и распространении волн цунами. Отмечается необходимость создания комплекса, позволяющего моделировать цунами от динамического сейсмического источника на основе клавишной модели.
Программные комплексы, моделирующие цунами
Огромное количество работ посвящено проблемам математического и численного моделирования цунами от динамического источника. В настоящее время существует ряд численных моделей и программных комплексов [60, 77, 88, 89].
Одним из наиболее современных продуктов, используемых при моделировании генерации и распространении волн цунами, является NAMI DANCE. Данный программный продукт позволяет осуществлять моделирование цунами на основе уравнений мелкой воды в потоках. Источники цунами в нем моделируются как гипотетические источники заданной формы. Пользователю дается возможность определить время, прихода волны, ее форму, максимальные заплески на побережьях, определить поле скоростей и динамику водной поверхности, а также определить предполагаемый ущерб. Данный программный продукт обладает хорошей графической оболочкой, однако, несмотря на все свои достоинства, он не позволяет осуществлять моделирование волн цунами от динамических источников. Кроме того, разностные схемы, используемые при моделировании, накладывают большие ограничения на шаг по времени, что замедляет время расчетов. В связи с этим возникла необходимость создания нового программного продукта, в котором эти недостатки были бы устранены.
Как известно, формирование цунами зависит от характера и динамики смещений дна в зоне очага землетрясения, а точнее, от начальных смещений дна. Как правило, при расчетах генерации волн цунами используют сейсмические данные для определения ориентации разрыва в очаге и энергии цунами. Затем решается статическая задача гидродинамики о пересчете распределений перемещений дна на поверхность океана. Полученные смещения водной поверхности принимаются за начальные условия, и проводится расчет распространения волны в данной акватории с учетом реальной батиметрии. В настоящее время существует ряд численных моделей и программных комплексов (см., напр., [106]), позволяющих достаточно точно проводить расчеты распространения волны цунами до побережья. После Индийского цунами 26 декабря 2004 года точность таких расчетов можно оценить, сравнив ЗБ-срезы в Индийском океане со спутниковыми данными по смещению водной поверхности при распространении цунами (см., напр., [24]). Однако, встает вопрос об адекватности модели очага, применяемой при таких расчетах. Особенности генерации цунами, начальная скорость, его параметры, характеристики на берегу (особенно в ближнеполевой зоне) напрямую зависят от выбора модели для определения сейсмического очага землетрясения.
К настоящему времени разработан механизм сильных землетрясений в зонах субдукции [34, 90, 94]. Известно, что узкие сейсмические пояса Земли связаны с условиями контакта на границах крупных литосферных плит. Взаимодействие плит в зоне поддвига является ответственным за сейсмический процесс в островных дугах и активных окраинах континентов. Наиболее сильные землетрясения происходят в зонах субдукции в окрестности пологой плоскости контакта между подошвой островодужного выступа и кровлей пододвигаемой плиты. Как показывают многочисленные геоморфологические и геолого-геофизические данные, островодужный выступ состоит из отдельных крупных сегментов, образованных поперечными разломами, проникающими вплоть до кровли пододвигаемой плиты. Наличие поперечных разломов потребовало ввести новые, более мелкие элементы взаимодействия — блоки (клавиши) с фронтального края нависающей плиты. Выяснилось, что такое минимальное усложнение традиционной- схемы субдукции вполне достаточно для успешного объяснения основных закономерностей сейсмического процесса в зонах субдукции [94]. Характерный размер блоков-клавиш составляет около 100 км. Такая «нарезка» на блоки фронтальных частей островных и материковых окраин структурно определяет размер очага сильного землетрясения (см. [98]). В основном, такие очаги связаны с деформируемыми и «стреляющими» при сбросе напряжений блоками-клавишами. Но иногда, длина очагов сильнейших землетрясений соответствует нескольким смежным блокам, в которых происходит одновременная разрядка накопленной упругой энергии. Можно предположить, что в декабре 2004 г., в Индийском океане одновременно «выстрелило» сразу 8 или 10 блоков-клавиш Зондской островной дуги, и этот мощный «аккорд» вызвал образование громадного очага землетрясения и, как следствие, появление гигантского цунами.
Прохождение глубоководного желоба с океанской стороны Южных Курил делает район островов Курильской гряды потенциально сейсмически опасной зоной, причем подводные землетрясения в зоне субдукции в большинстве случаев являются цунамигенными. Возникновение цунами при таких землетрясениях таит в себе опасность для районов побережья Охотского моря, в том числе и Сахалина, где расположены как производственные объекты, так и жилые сооружения. В истории было зафиксировано ряд случаев атак цунами на побережье Сахалина, вызванных землетрясениями в районе Курильской гряды. Тестирование комплекса можно осуществлять, сравнивая его результаты с известными данными землетрясения на Суматре [8, 26, 46]. Хотя интенсивность цунами в существенной мере определяется силой подводного землетрясения, высота волны на побережье Сахалина редко достигала рекордных величин, известных в истории цунами на тихоокеанском побережье [23, 49, 92]. Например, высота волны на побережье Сахалина от Урупского землетрясения 13 октября 1963 г. достигала только полуметра. Это связано в первую очередь с тем, что распространение волны цунами в сторону Сахалина ограничено в первую очередь наличием высоких подводных хребтов (таких как например хребет Витязя), отделяющих акваторию Охотского моря от источника цунами в зоне глубоководного желоба, так что естественными волноводами для цунами служат только два глубоководных пролива: пролив Буссоль и пролив Крузенштерна.
Итерационно-интерполяционный метод
Для чистоты эксперимента были выгружены все работающие программы, отключены антивирусы и остановлены все неиспользуемые службы операционной системы. После чего на четырех-ядерном персональном компьютере был запущен пакетный файл запуска контрольных примеров. Спустя -8 часов в файлах newl.txt, new2.txt, new5.txt, old5.txt оказались результаты использования памяти и времени работы каждого варианта программы (оптимизированного и исходного) для каждого примера (№№ 1 и 5) и для каждого варианта размера расчетной сетки (от 250 до 4000 с шагом 250).
При помощи распараллеливания вычислений по второму предложенному методу на четырех-ядерном персональном компьютере мы добились уменьшения времени счета на 60%, что соответствует 2,5 кратному увеличению производительности.
Метод решения системы уравнений, предложенный профессором С. Патанкаром основывается на очень эффективном алгоритме, известным как Рис. 2.1.19. Ускорение вычислении, полученное благодаря распараллеливанию главного счетного цикла 65 TDMA (TriDiagonal-Matrix Algorithm). Название TDMA происходит из свойства матрицы коэффициентов дискретных аналогов: в этой матрице ненулевые элементы располагаются только на трех смежных диагоналях. Рассмотрим уравнение теплопроводности: г дф X v dt , дх, дх, + S . (2.1.44) Выражение для его дискретного аналога будет иметь следующий вид: а рф р = аЕфЕ + ап,фуу + + ахф5 + Ь , (2.1.45) где - индекс расчетной точки, "Е" — индекс соседней справа точки, " УУ" - индекс соседней слева точки, "ТУ" — индекс соседней сверху точки, - индекс соседней снизу точки.
Как легко заметить, при расчете очередной точки, используются значения четырех соседних точек. Если быть точнее, используются значения коэффициентов четырех соседних точек, а также значения четырех соседних ф . Иными словами, при реализации данного алгоритма в качестве программы для вычислительной системы, мы имеем дело с одновременным обращением к четырем ячейкам памяти, выполнением нескольких простых арифметических действий и записью результата в пятую ячейку памяти (ОЗУ).
Известно [82], что в силу определенных обстоятельств (сложность изготовления, а следовательно, и высокая стоимость оперативной памяти, работающей на той же частоте, что и процессор) оперативная память работает на достаточно низкой, по сравнению с процессорной, частоте. Отсюда вытекает и большое время доступа к ячейкам оперативной памяти, составляющее несколько десятков и даже сотен тактов процессорного времени. К счастью существует память с более низким временем доступа и работающая на той же частоте, что и процессор. Объем ее достаточно скромен и составляет порядка 16-128 Кб. (Для сравнения, объем оперативной памяти компьютера составляет несколько гигабайт, что в несколько десятков тысяч раз больше объема КЭШ-памяти). Данная память получила название «сверхоперативной» или КЭШ (cache). Ее основная задача состоит в том, чтобы хранить те данные, с которыми оперирует процессор в данный момент.
Копирование данных из КЭШа в ОЗУ и обратно происходит пакетами, размером в 32-128 байт, называемым КЭШ-линией. Т.е. при чтении любой ячейки из памяти, происходит чтение и копирование в КЭШ целой КЭШ- линии. После чего, если следующая затребованная ячейка оказывается в этой же КЭШ-линии, процессор берет ее уже из КЭШа. Доступ к данным, находящимся в КЭШе, происходит на 1-2 порядка быстрее. Поэтому наиболее выгодно обращаться к ячейкам в таком порядке, чтобы они все были в одной и той же КЭШ-линии. Заметим, что при полностью заполненном КЭШе происходит замещение КЭШ-линии на новую. Существуют различные алгоритмы выбора замещаемой КЭШ-линии, однако в нашем случае это не важно. Важно то, что мы должны свести к минимуму операции замещения КЭШ-линий.
В предложенном программном комплексе CONDUCT (профессора С. Патанкара) автор предлагает размещать коэффициенты дискретного аналога в отдельные массивы. Таким образом, при решении системы уравнений мы обращаемся к ячейкам памяти, находящимся в разных массивах, и как следствие, они не могут находиться рядом друг с другом. В предложенной реализации используется восемь массивов коэффициентов. И ко всем им мы обращаемся одновременно, что, конечно же, приводит к постоянному замещению КЭШ-линий разных массивов. Это связано с тем, что для увеличения скорости работы КЭШ поделен на несколько банков (чаще всего на 4), состоящих из нескольких КЭШ-линий. Количество банков определяет ассоциативность КЭШа. И в каждом банке каждая строка закреплена за определенными адресами ОЗУ, кратными количеству КЭШ-линий в банке. Т.е. если мы имеем 4 банка по 64 КЭШ-линии размером 64 байта (полный объем равен 4 64 64 = 16 кб), то каждому блоку в 64 байта ОЗУ могут соответствовать 4 КЭШ-линии (по одной из каждого банка), а каждая N-ая КЭШ-линия соответствует каждому N-ому блоку в 64 байт ОЗУ.
Как легко убедиться, если запрашиваемые ячейки будут находиться в 64-байтных блоках ОЗУ с кратными номерами (имеется в виду остаток от деления общего объема ОЗУ на размер КЭШ-банка) и их количество превысит величину ассоциативности КЭШа, то КЭШ-линии будут постоянно перезатираться новыми, что приведет к существенному падению производительности.
Для исключения подобной ситуации нам необходимо постараться расположить данные в памяти таким образом, чтобы работать с ячейками, находящимся в одной КЭШ-линии (это в идеале, а на практике — свести к минимуму число одновременно затребованных КЭШ-линий).
В данной работе предлагается организовать один массив, элементом которого является структура из восьми коэффициентов, ранее располагавшихся в совершенно разных массивах.
С помощью такой организации данных нам понадобится одновременно запросить максимум 4 КЭШ-линии, а чаще всего только 3, т.к. коэффициенты точек Е, Р и W располагаются в ОЗУ рядом и будут лежать в одной КЭШ-линии (или в двух соседних), еще две КЭШ-линии потребуются для коэффициентов точек S и N.
Во время проведения тестирования и замеров производительности для различного размера расчетной сетки было обнаружено странное поведение первоначальной версии программного комплекса CONDUCT (рис. 2.1.20).
Нелинейные уравнения мелкой воды в потоках
Система уравнений химической кинетики является жесткой, т.к. константы скоростей для разных реакций сильно различаются, а следовательно времена выхода на равновесные режимы для различных реакций могут отличаться друг от друга на несколько порядков. (2.5.4) В общем виде уравнение химической кинетики можно записать в виде Г = а(Са)-СаУ2и(Са),сс = 1,...,//, где У1а и У2а - положительные непрерывно дифференцируемые функции своих аргументов. Как показал Я.Б.Зельдович, система уравнений (2.5.4) имеет единственное стационарное состояние в // -мерном пространстве. По определению Са О при любых / 0. Видно, что правая часть состоит из положительного и отрицательного слагаемого. Поэтому поступим, как и в методе Патанкара при реализации линеаризованного источникового члена.
Равновесные значения концентраций компонентов можно определить с1С /п1 V, (С ) из условия —— = 0, т.е. С = ) (, поэтому в окрестности равновесия & У7а\Са) можно представить в виде Са (2.5.5) Л где г„ - характерное время химической релаксации а — компонента в результате т химической реакций. Если в результате химических реакций «-компонент только исчезает, то С = 0 и уравнение (2.5.5) упрощается. Для неравновесных термохимических состояний возможна оценка характерных времен химической релаксации, следовательно, и степени жесткости системы. Такой подход существенно проще, чем определение числа обусловленности матрицы правых частей для рассматриваемой системы уравнений.
В связи с тем, что при решении рассматриваемых задач разброс собственных значений матрицы Якоби для правых частей системы приводит к большим трудностям при их численном решении, использование классических методов типа Рунге-Кутты и Адамса затруднительно из-за очень жестких требований на шаг интегрирования, которое получается из условия устойчивости. Таким образом, жесткость системы уравнений предъявляет повышенные требования к устойчивости применяемых для их решения разностных схем, т.е. требуется применять пошаговые методы специальных типов.
В анализе таких методов может быть использована значительная часть теории методов аппроксимации на бесконечных интервалах, поскольку жесткая задача Коши проявляет те свойства на коротких интервалах, какие обычно наблюдаются при увеличении длины интервала до бесконечности.
Сложность рассматриваемых задач заставляет использовать не только А -устойчивые методы, но и получать экономичные разностные схемы на основе метода расщепления. Требование экономичности диктует необходимость использовать одношаговые разностные схемы.
Следуя методу Патанкара, из (2.5.4) легко можно получить неявную разностную схему. Расписывая в конечных разностях производную по времени, при условии, что линейный множитель в правой части берется на верхнем временном слое — = У1а(Са)-С+1Ка(Са), а = Выражая А? , будем иметь С" +д/1/"+1 1 = Са+ЛГК (2 1 + А Г1
При этом и сами скорости химических реакций можно считать на (п)- ом слое по времени и уточнять в результате итераций. Легко видеть, что схема (2.5.6) является -устойчивой, т.к. при А? - о мы получаем уравнение для определения равновесной концентрации. Если правую часть (2.5.4) аппроксимировать полусуммой значений источников на нижнем и верхнем временных слоях, то сг1 - с: 1 & 2 Выражая С("+1, будем иметь СГ1 = 2 . (2.5.7) 1 +
Полученная схема, также как и схема (2.5.6), является тоже А- устойчивой.
Обычно схема Кранка-Николсона считается безусловно устойчивой, объясняя этот факт тем, что решение будет получено независимо от значения шага по времени. Устойчивость же в математическом смысле просто гарантирует, что колебания в решении, в конечном счете, будут затухать. Однако, такая устойчивость не гарантирует физически правдоподобного результата. Поэтому, наряду с -устойчивостью метода, вводят понятие Ь- устойчивости, которая обеспечивает однозначное соответствие быстроубывающих компонент разностного дифференциального уравнения.
Тестирование разностной схемы, полученной на основе итерационно- интерполяционного метода
Рассмотрим однородный бесконечный массив, в котором существует область с повышенной температурой (очаг возгорания). В массиве выделим несколько слоев, каждый со своими характеристиками (рис. 4.1.1). Кроме перечисленных слоев, к модели добавляется пограничный слой (приземный слой атмосферы) с учетом внешнего потока. Систему координат введем таким образом, чтобы оси Ох и Оу лежали в горизонтальной плоскости, а ось Ох была направлена вертикально вверх. Начало координат расположим в непосредственной близости от одного из очагов возгорания. очага горения; / - его длина, а й„ - высота факела пламени (проекция / на ось Ог); Це, Ус IV,. - компоненты скорости ветра над пологом; ух,уг,у: — углы отклонения факела пламени от осей Ох, Оу, Ог соответственно, а — угол наклона почвы.
Считается, что органическая масса слоев представляет собой реакционно-способную пористую сплошную многофазную среду, состоящую из сухого органического вещества, воды в жидко-капельном состоянии, конденсированных продуктов пиролиза, обогащенных углеродом, минеральной части (золы), газовой и дисперсных фаз. Для простоты считаем СО, СН4 и другие горючие компоненты, входящие в состав летучих продуктов пиролиза одним эффективным горючим газом, а СОг и другие инертные компоненты — эффективным продуктом реакций, получая таким образом газовую фазу, состоящую из трех компонент: окислителя (02), горючего газа (в качестве эффективного горючего газа принимаем СО, как преобладающего среди горючих компонентов продуктов пиролиза), и СОо совместно с другими инертными компонентами газовой фазы. Перенос энергии от фронта горения к еще не горящему топливу в общем случае осуществляется путем кондукции, конвекции и излучения. Излучение играет значительную роль в непосредственной близости от фронта горения, в том числе излучение от факела пламени.
С точки зрения переноса излучения слой является излучающей, поглощающей и рассеивающей средой. Излучает и поглощает лучистую энергию конденсированная фаза слоя, газовая же фаза при температурах до 1500К оптически прозрачна. Длина свободного пробега излучения в слое /д= — -Я4, где х _ интегральный коэффициент поглощения, — удельная X поверхность слоя. Поскольку в слое 1м"1, то 1К «Л, где Л = ]\п 7, с„ш - 1хп гр сюя характерный размер (высота слоя), следовательно, перенос излучением можно рассматривать в диффузионном приближении. Над слоем рассматривается факел пламени. Считается, что этот факел излучает как плоская стенка длиной /, расположенная под некоторым углом я-у. к поверхности слоя (рис. 4.1.1.), величина которого зависит от скорости ветра. объемная доля газовой фазы; Т - температура газовой фазы; Са — массовые концентрации компонентов газовой фазы (а= 1 - кислород, а=2 - горючие компоненты продукта пиролиза, сс=Ъ — инертные компоненты газовой фазы, не реагирующие компоненты продукта пиролиза и водяного пара); и, V, Ж- проекции скорости над слоем по осям х, у, г соответственно; Р — давление на внешней границе пограничного слоя; Я5 — массовая скорость реакции горения летучих продуктов пиролиза; Я5а — массовые скорости образования компонентов газовой фазы; 75 - тепловой эффект реакции окисления летучих продуктов пиролиза; г5 1 - доля теплоты газофазной реакции окисления газообразных продуктов пиролиза, усвоенная конденсированной фазой; М( Мс, М — молярные массы индивидуальных компонентов, углерода и смеси в целом; Е5, кэ - энергия активации и предэкспонента химической реакций горения летучих продуктов пиролиза; - коэффициенты динамической вязкости, турбулентной теплопроводности и турбулентной диффузии соответственно; g — ускорение свободного падения.
Здесь с ,р1,(р1 - удельные теплоемкости, истинные плотности и объемные доли /-ой фазы многофазной реагирующей среды(г=1 — сухое органическое вещество, /=2 - вода в жидко капельном состоянии, /=3 - конденсированные продукты пиролиза, г-4 - минеральная часть — зола); ср5,р5,(р5 - удельная теплоемкость, истинная плотность и объемная доля газовой фазы; Т - температура газовой фазы; Т5 — температура твердой фазы; Са - массовые концентрации компонентов газовой фазы (а=1 - кислород, а=2 — горючие компоненты продукта пиролиза, а=3 - инертные компоненты газовой фазы, не реагирующие компоненты продукта пиролиза и водяного пара); II, V, УУ- проекции скорости над слоем горючих материалов по осям х, у, г соответственно; Р - давление в потоке; С/д — плотность потока излучения; К - коэффициент ослабления; к5 — спектральный коэффициент поглощения; с — скорость света; а — постоянная Стефана-Больцмана; Я/, Яо, Я3, Я5 — массовые скорости реакции пиролиза сухого органического вещества, испарения влаги, горения конденсированных и летучих продуктов пиролиза соответственно; R5a — массовые скорости образования компонентов газовой фазы; q2, q3, qs тепловые эффекты реакций испарения, горения кокса и окисления летучих продуктов пиролиза; Q - массовая скорость образования газовой фазы; As - коэффициент межфазного теплообмена; v5 l — доля теплоты газофазной реакции окисления газообразных продуктов пиролиза, усвоенная конденсированной фазой; Ма, Мс, М - молярные массы индивидуальных компонентов, углерода и смеси в целом; S - удельная поверхность фитомассы полога, cd - эмпирический коэффициент сопротивления; E.,kn (/ = 1 ..3,5) - энергии активации, и предэкспоненты химических реакций; ac,vr - коксовое число и массовая доля горючего газа в общей массе летучих продуктов пиролиза; ut,A.t,Dt - коэффициенты динамической вязкости, турбулентной теплопроводности и турбулентной диффузии соответственно; g — ускорение свободного падения.
Выше изложенная система описывает протекающие процессы в слое layer, где layer = 1 ..(A/ /aur -1), Niayer - общее количество слоев, в рассматриваемой модели включая приземный слой атмосферы (layer=Niayer).
При записи системы уравнений (4.2.18)-(4.2.34) учитывался эффект двухтемпературности среды [105], поэтому уравнения сохранения энергии были записаны отдельно для газовой и дисперсной фаз.
