Содержание к диссертации
Введение
Глава 1. Математическое моделирование техногенных аварий с распространением тяжелых газов и разлитием жидкостей 16
1. Математическая модель распространения облаков тяжелых газов над орографически неоднородной поверхностью 16
2. Алгоритм и методы численного решения задачи 27
3. Математическое моделирование аварии 03.06.89 под Уфой 36
4. Математическое моделирование разлитии нефти 41
Глава 2. Математическое моделирование лесных пожаров 47
1. Двумерная двухфазная модель лесных пожаров 47
2. Двумерная трехфазная модель лесных пожаров 71
3. Двухъярусная модель лесных пожаров 77
4. Алгоритм и методы численного решения задачи 80
5. Результаты математического моделирования лесных пожаров по двухфазной, трехфазной и двухъярусной моделям 87
Глава 3. Математическое моделирование некоторых задач экологии, описываемых двумерными моделями механики твердого упругого тела 106
1. Математическая постановка задачи 106
2. Разностная аппроксимация задачи 122
3. Устойчивость разностной схемы 160
4. Сходимость решения разностной задачи к обобщенному решению дифференциальной задачи 188
5. Моделирование колебаний ледяного покрова под действием техногенных динамических нагрузок 231
Заключение 236
Список литературы
- Алгоритм и методы численного решения задачи
- Математическое моделирование аварии 03.06.89 под Уфой
- Двухъярусная модель лесных пожаров
- Сходимость решения разностной задачи к обобщенному решению дифференциальной задачи
Введение к работе
Целью работы является разработка математических моделей и численных методов, позволяющих построить эффективную численную реализацию и провести численное моделирование рассматриваемых задач промышленной безопасности и экологии. Работа состоит из трех глав. Первая глава посвящена математическому моделированию техногенных аварий с распространением тяжелых газов и разлитием опасных жидкостей, описываемых моделями гидрогазодинамики, в которых среда полагается термически однородной. Вторая глава посвящена математическому моделированию лесных пожаров, которые описываются многофазными моделями с газовой, твердой и дисперсной фазами, при этом среда является химически реагирующей и многотемпературной. В третьей главе рассматриваются задачи экологии, описываемые моделями механики твердого упругого тела с нестационарным уравнением колебаний тонких упругих пластин с приложениями к задаче о колебаниях ледяного покрова под действием техногенных динамических нагрузок. Таким образом, рассматриваются нестационарные постановки задач экологии и промышленной безопасности для разных сплошных сред с различными физико-химическими свойствами. В соответствии с принятой в математическом моделировании задач экологии и промышленной безопасности традицией разделения моделей на локальные (до нескольких десятков километров), мезомасштабные (сотни километров) и глобальные (тысячи километров) все рассматриваемые модели относятся к локальным моделям приповерхностного слоя. Рассматриваемые модели — это модели оперативной оценки развития ситуации, с помощью которых можно на персональных ЭВМ рассчитать динамику рассматриваемых процессов и оценить их последствия. Исходя из этих требований, математические модели всех рассматриваемых задач построены на основе единого подхода, основанного на фундаментальных законах рассматриваемых процессов и позволяющего построить эффективную численную реализацию. Для рассматриваемых задач также характерен единый подход к построению алгоритмов и методов численного решения - задачи решаются конечно-разностными методами с применением алгоритмов расщепления, как по физическим процессам, так и по пространственным переменным.
Среди актуальных сегодня задач промышленной безопасности важное значение имеют задачи, связанные с обеспечением безопасности населения и защиты окружающей среды в случае аварий на объектах хранения и транспортировки сжиженных горючих и токсических газов, а также аварий на нефтепроводах.
По данным [1] в России действует более трех тысяч промышленных объектов, на которых находится около 700 тыс. тонн химически опасных веществ, таких как хлор (около 35% от общего количества), аммиак (около 50%) и другие. Не меньшую опасность представляют и объекты (резервуары, трубопроводы) со сжиженными горючими газами, прежде всего нефтяными (пропан, бутан и их смеси). Особую опасность представляет то обстоятельство, что промышленные объекты со сжиженными горючими и токсическими газами часто расположены в непосредственной близости или на самой территории населенных пунктов. Так, более 70% предприятий химической промышленности расположены в крупных городах с населением свыше 100 тысяч человек. Пути транспортировки сжиженных газов автомобильными или железнодорожными цистернами также проходят по территории населенных пунктов. Следовательно, техногенные аварии на подобных объектах затронут население прилегающих территорий и могут нанести большой ущерб.
Примеров таких крупных аварий в мировой практике достаточно много (см. [2]). Приведем примеры двух подобных аварий, произошедших на территории СССР во второй половине прошлого века:
-10 января 1966г. в г. Горький в результате аварии из железнодорожной цистерны на сортировочной станции вылилось 23 тонны хлора, что привело к поражению более 2000 человек;
- 3 июня 1989 г. под Уфой (м.Аша) в результате разрыва трубопровода вытекло до 10 тыс. тонн пропан-бутановой смеси. Образовалось облако топливо-воздушной смеси, которое растеклось по пересеченной местности и затем воспламенилось. Зона поражения в результате взрыва составила 2.5 км . В попавших в зону поражения двух проходивших по железной дороге поездах пострадало 623 человека, погибло 575 человек.
В качестве опасных веществ в работе рассматриваются горючие и токсические газы, хранящиеся в сжиженном состоянии при повышенном давлении и при температуре окружающей среды. Эти вещества имеют критическую температуру выше, а точку кипения ниже температуры окружающей среды. При разгерметизации (разрушении) резервуара или трубопровода сжиженный газ практически мгновенно испаряется из разлития и образует облако газовоздушной смеси. В работе рассматриваются так называемые тяжелые газы, истинная плотность которых при атмосферном давлении в несколько раз выше плотности воздуха, поэтому облако газа не поднимается вверх, а распространяется над подстилающей поверхностью под действием гравитационных сил и ветра. Это достаточно широкий класс опасных газов, к которым относятся все горючие
нефтяные газы - пропан, бутан и их смеси, и самый широко используемый в промышленном производстве токсический газ - хлор. В силу невозможности проведения натурных экспериментов с взрывоопасными и токсическими газами в условиях реальных промышленных объектов единственным инструментом для исследования процесса распространения облаков тяжелых газов является метод математического моделирования. На основе математического моделирования можно прогнозировать динамику вероятных аварий на объектах хранения и транспортировки сжиженных горючих и токсических газов, проводить оценку их последствий, а также проводить реконструкцию событий и анализ уже произошедших аварий.
Работы по математическому моделированию в задаче распространения облаков тяжелых газов велись за рубежом, начиная с 60-х годов прошлого века. Как правило, создавались упрощенные математические модели, основанные на эмпирических данных (см. обзоры [3,4].) Например, широко использовались гауссовы модели, в которых концентрация газа в облаке рассчитывается по алгебраической формуле с экспоненциальным уменьшением по мере удаления от источника выброса. Или, например, модели "DENZ" [5] и "CRUNCH" [6], в которых радиус расширения облака рассчитывался из решения простейшего обыкновенного дифференциального уравнения. Будучи важной научной поддержкой инженерных оценок риска, такие модели описывают поведение интегральных параметров физических явлений, они основаны на интерполяции данных ограниченного числа экспериментов, не учитывающих всех условий реальных промышленных объектов и окружающего пространства, и поэтому не могут дать адекватную картину рассматриваемого явления. Другим направлением в математическом моделировании рассматриваемого процесса стало создание осесимметричных моделей [3,4]. Эти модели, хотя и были построены на основе законов сохранения механики сплошной среды, не могли учесть рельефа подстилающей поверхности и наличия застройки, а потому плохо описывали развитие рассматриваемых аварий в условиях реальных объектов.
В 1 главы 1 разработана двумерная математическая модель распространения облаков тяжелых газов. Созданная автором в конце 80-х годов в Институте атомной энергии (ИАЭ) им. И.В.Курчатова [135,136] эта модель явилась первой двумерной моделью, в которой система уравнений рассматриваемого процесса была построена на основе законов сохранения газовой динамики и которая учитывала все основные физические явления рассматриваемого процесса: гравитационное растекание, турбулентность, трение о подстилающую поверхность и препятствия, разбавление газа в
облаке окружающим воздухом, а также учитывала наличие ветра. Система уравнений рассматриваемой двумерной модели была построена методом осреднения трехмерных уравнений газовой динамики по высоте течения при физически обоснованных предположениях о характере течения. Проинтегрировав осредненные по Рейнольдсу трехмерные уравнения по высоте от поверхности рельефа zo(x,y) до неизвестной свободной верхней границы облака H(x,y,t) и решая затем двумерную задачу, мы получаем, по сути, трехмерную картину течения с учетом неоднородного рельефа местности, наличия препятствий, промышленной и жилой застройки. Построенная модель является важным обобщением известной модели течения несжимаемой жидкости в приближении теории мелкой воды [7] на случай газодинамических течений с переменными по (х,у) плотностью и массовыми концентрациями компонентов газовой среды. По этим параметрам можно судить о достижении предельно допустимых концентраций (ПДК) в облаке и оценивать вероятность токсического поражения персонала промышленных объектов и населения прилегающих жилых районов. При распространении облаков горючих газов, пользуясь рассчитанными в рассматриваемой модели параметрами, можно на основе инженерных методик [8] или на основе математического моделирования процесса сгорания облака провести оценку поражающих факторов при сгорании облака и оценить возможный ущерб от таких аварий.
В 2 главы 1 разработаны алгоритм и численные методы решения рассматриваемой задачи. Алгоритм численного решения задачи основан на известном методе расщепления по физическим процессам. Все отдельные подсистемы решались разностным методом на прямоугольных сетках.
Был создан программный комплекс для численного моделирования распространения облаков тяжелых газов. Первый вариант программного комплекса был создан автором в 1989-1990 гг. в ИАЭ им И.В.Курчатова в операционной среде UNIX для ЭВМ НР-1000. Программный комплекс был тестирован путем сравнений результатов расчетов с результатами натурных экспериментов по распространению фреона, проведенных на острове Торни (Великобритания) в 1982- 1984гг. [96-98]. С помощью этого программного комплекса была промоделирована стадия гравитационного растекания облака пропан-бутановой смеси по пересеченной местности при упомянутой выше аварии 03.06.89 под Уфой, и результаты моделирования [135] совпали с данными, собранными экспертами по расследованию аварии. Результаты моделирования аварии под Уфой представлены в 3 главы 1. Впоследствии в 90-х годах на основе разработанной автором модели и метода расщепления Н.П.Савенковой и С.В.Филипповой при участии
автора был создан программный комплекс для IBM-PC совместимых компьютеров в программной среде PASKAL [144-146], который затем в 2003 году был внедрен у заказчика - в Федеральном агентстве правительственной связи и информации.
В настоящее время с появлением современных суперкомпьютеров с параллельными вычислениями стало возможным моделировать рассматриваемый процесс на основе трехмерных математических моделей, однако такое моделирование может быть осуществлено на базе научно-исследовательских центров, оснащенных суперкомпьютерной техникой. В то же время для экспертных систем прогнозирования динамики таких аварий и оценки их последствий в центрах управления и реагирования при чрезвычайных ситуациях а также в исследовательских подразделениях промышленных объединений и концернов, по-прежнему используются программные пакеты для персональных ЭВМ, созданные на основе более простых для численной реализации математических моделей, либо инженерных методик (см. [1]). Подтверждением актуальности такого подхода к математическому моделированию рассматриваемого процесса сегодня является широкое распространение в развитых странах коммерческих пакетов EFFECTS (фирмы TNO) и WHAZAN (Technica Ltd).
Другой актуальной особенно для России задачей промышленной безопасности является задача борьбы с разлитиями нефти из магистральных трубопроводов. По результатам математического моделирования процесса разлития нефти может быть определено необходимое количество сил и средств ликвидации разлива. На основе сведения модели распространения тяжелых газов к упомянутой выше модели мелкой воды и соответствующей модификации программного комплекса было проведено математическое моделирование разлитии нефти по орографически неоднородной поверхности, результаты которого представлены в 4 главы 1.
Основные результаты по первой главе опубликованы в работах [135,136,143-149].
Одной из чрезвычайно важных задач экологии является задача сохранения лесных покровов Земли, как регуляторов необходимых для жизни балансных соотношений газов и влаги в атмосфере Земли. Наряду с вырубкой лесов наиболее сильный ущерб лесным массивам наносят пожары. Кроме того, при пожарах происходит уменьшение количества кислорода в атмосфере и выброс газообразных продуктов горения, прежде всего углекислого газа, что также влияет на газовый баланс в атмосфере Земли. Как отмечено в [9], согласно имеющимся оценкам около 30% тропосферного озона, окиси углерода и углекислого газа, содержащихся в атмосфере, обусловлено вкладом лесных пожаров, а связанные с лесными пожарами выбросы аэрозоля в атмосферу могут оказывать
существенное влияние на микрофизические и оптические характеристики облачного покрова и, следовательно, на климат.
В России, в силу обширности территорий, покрытых лесом, и недостатка средств по ликвидации пожаров, проблема борьбы с пожарами стоит особенно остро. В пожароопасный сезон на территории страны ежедневно возникают сотни очагов лесных пожаров. В настоящее время в России и за рубежом интенсивно ведутся аэрокосмические исследования лесных пожаров и их влияния на окружающую среду [9-12]. В агентстве МЧС России по мониторингу и прогнозированию чрезвычайных ситуаций эта работа ведётся с помощью геоинформационной системы (ГИС) [1]. Однако, для определения эффективных сценариев реагирования недостаточно мониторинга чрезвычайной ситуации, а требуется прогноз её дальнейшего развития. Такой прогноз можно дать с помощью метода математического моделирования чрезвычайной ситуации, в рассматриваемом случае - математического моделирования лесных пожаров.
В области математического моделирования лесных пожаров к настоящему времени существуют различные подходы к описанию этого сложного явления. Одним из направлений является разработка вероятностно-статистических моделей [13-18], использующих, как правило, аппарат клеточных автоматов. Однако, в таких моделях не учитываются основные физико-химические процессы, происходящие во фронте пожара. В рамках этих моделей также нельзя адекватно описать как процесс выгорания лесных горючих материалов с учетом их неоднородного распределения на местности, так и влияние переменных внешних факторов, например, ветра. Без учета этих факторов правильно описать динамику распространения лесных пожаров практически невозможно. Существует также класс эмпирических моделей расчета распространения фронта пожара [19-25], которым присущи те же недостатки. В рамках двумерных (x,z)-моделей рассчитывается распространение лесных пожаров по заданному направлению [26,36,49], что также не может дать правильную картину распространения реального лесного пожара на местности. Широко представлен класс моделей, в рамках которых изучаются отдельные явления, происходящие в зоне лесного пожара и в атмосфере над пожаром [27-32,36-39]. Например, явления поднятия термиков и образования конвективных колонок моделируются на основе двумерных осесимметричных моделей [27-30,36-39]. В работах последнего времени [31,32] на основе трехмерных моделей численно получены трехмерные когерентные вихревые структуры в атмосфере над лесным пожаром. В работах [33-35] рассматривается трехмерная однофазная газодинамическая модель лесных
пожаров без учета химических реакций, в которой восходящие потоки тепловой энергии аппроксимируются аналитической формулой, а перенос излучения не учитывается.
Фундаментальный вклад в развитие математических моделей лесных пожаров внесен А.М.Гришиным и его учениками [36-52]. Им разработана наиболее полная трехмерная модель лесных пожаров, в которой лес является девятиярусной (по высоте) многофазной (8 фаз) реакционно-способной средой, газовая фаза которой описывается системой трехмерных уравнений газовой динамики. Однако, расчет реального лесного пожара продолжительностью несколько часов, а тем более суток, по этой трехмерной модели, даже с использованием суперкомпьютеров, на сегодняшний день представляется трудноразрешимой задачей.
В настоящей работе в главе 2 построены относительно простые для численной реализации двумерные математические модели, назначение которых - описать динамику реальных лесных пожаров для экспертной оценки развития ситуации и выработки управленческих решений по тушению пожаров, а также для оценки ущерба от пожаров. В то же время эти модели отражают основные физические законы сохранения вещества, импульса и энергии, учитывают неоднородное распределение запасов лесных горючих материалов на местности и наиболее существенные для динамики пожара физические явления, происходящие в зоне пожара. В 1 главы 2 построена двумерная двухфазная модель лесных пожаров. При построении этой модели, как и в модели распространения тяжелых газов, был использован метод осреднения исходных трехмерных уравнений газовой динамики по высоте слоя лесных горючих материалов, в котором параметры среды, например, истинная плотность и влагосодержание лесных горючих материалов, сила ветра и другие, предполагались постоянными. В модели рассматривается многокомпонентная газовая фаза, состоящая из горючих газов, негорючих газов, дисперсной сажи и окислителя, многокомпонентная твердая фаза, состоящая из лесных горючих материалов и твердых продуктов разложения лесных горючих материалов в результате химических реакций. Модель учитывает конвективный перенос и турбулентное трение, теплопроводность и перенос энергии излучения в газовой фазе, межфазное трение, межфазный массо-и теплообмен и обмен лучистой энергией, а также изменение объемных долей фаз и тепловыделение за счет химических реакций в газовой и твердой фазах. Кроме того, модель учитывает неоднородность распределения запасов лесных горючих материалов по площади. Однако, построенная одноярусная модель не учитывает процессов, происходящих над этим слоем лесных горючих материалов и под ним.
Из наблюдений известно, что в реальных лесных пожарах существует механизм распространения пожара с помощью крупных горящих частиц, которые переносятся ветром над пологом леса и, опускаясь, служат причиной образования новых очагов огня перед фронтом пожара. Для учета механизма распространения пожара горящими частицами в 2 главы 2 была построена трехфазная модель лесных пожаров, в которой горящие частицы, переносимые по ветру над слоем лесных горючих материалов, рассматриваются как отдельная фаза со скоростью и температурой, отличными от газовой фазы.
Также из наблюдений за реальными лесными пожарами известно, что верховой пожар (пожар в кронах деревьев) может самостоятельно распространяться лишь при наличии достаточно сильного ветра. При отсутствии ветра такой пожар затухает. Пожары же в нижних ярусах леса более устойчивы и могут распространяться и при отсутствии ветра. Часто наблюдается явление совместного распространения верхового и низового пожаров. При этом фронт низового пожара обгоняет фронт верхового и затем поджигает верхний слой лесных горючих материалов. Для описания этого процесса в 3 главы 2 построена двухъярусная модель лесных пожаров.
Таким образом, для моделирования процесса распространения реальных лесных пожаров автором впервые разработан комплекс относительно простых для численной реализации двумерных математических моделей, построенных на основе законов сохранения механики многофазных реагирующих сред и учитывающих неоднородность распределения лесных горючих материалов по площади, наличие препятствий для распространения огня и основные физико-химические явления, происходящие в рассматриваемом процессе.
Разработан алгоритм численного решения задачи (см. 4 гл. 2), построенный по тому же принципу, что и алгоритм численного решения задачи о распространении облаков тяжелых газов, о котором сказано выше. Создан программный комплекс и проведено численное моделирование процесса распространения лесных пожаров по двумерной двухфазной, двумерной трехфазной и двухъярусной моделям лесных пожаров в условиях неоднородного распределения запасов лесных горючих материалов по площади и наличия препятствий для распространения огня, таких как, дороги, просеки, реки, поляны и т.д. Результаты моделирования представлены в 5 главы 2. Было также проведено численное исследование чувствительности моделей по параметрам. В частности, исследовано влияние на распространение лесных пожаров таких параметров, как влагосодержание, количество дисперсной сажи, скорость ветра. Также исследовано
влияние на распространение лесных пожаров отдельных физических процессов -конвективного переноса, межфазного теплообмена, межфазного обмена лучистой энергии и т.д. Программный комплекс, построенный по двумерной двухфазной модели в 2003 году, был внедрен у заказчика - в Федеральном агентстве правительственной связи и информации.. Основные результаты по второй главе опубликованы в работах [136-138, 149-154,158].
Рассмотренные в третьей главе диссертации задачи экологии, описываемые двумерным нестационарным уравнением поперечных колебаний тонких упругих пластин, представляют большой интерес как с точки зрения развития теоретических основ математического моделирования в таких задачах (разработки численных методов и их обоснования), так и с точки зрения практических приложений. С прикладной точки зрения представляет интерес задача о распространении изгибно-гравитационных волн в ледяном покрове под воздействием внешнего давления. Такие задачи возникают, когда требуется оценить прочность ледяного покрова при воздействии на него движущегося поля давления, например, при движении по льду автотранспорта или при выборе на льду взлетно-посадочных полос для самолетов [53-57]. Также представляет практический интерес задача об определении возможности всплытия подводного аппарата на поверхность моря из под ледяного покрова [58]. Интерес представляет также задача о деформации ледяного покрова под воздействием сферической волны давления, распространяющейся в воде под ледяным покровом. Применима она также к задаче геоэкологии о поперечных колебаниях океанических литосферных плит, в которой тонкую, по сравнению с континентальной, океаническую литосферу, состоящую из относительно однородной породы (базальта), можно рассматривать в приближении тонкой упругой пластины [59,157,158].
Некоторые из перечисленных выше задач ранее решались аналитически. Так в [53] рассмотрено одномерное уравнение колебаний, а перемещение по льду плоского фронта давления задается в правой части уравнения б-функцией вида P5(x-\t), где Р -постоянная сила, действующая на поверхность ледяного покрова. Для этого случая получены критические значения скорости движения нагрузки, при которых частота колебаний ледяного покрова под действием нагрузки совпадает с частотой собственных колебаний льда на поверхности воды и в системе возникает резонанс, а прогиб пластины под действием нагрузки, а, следовательно, и внутренние напряжения в пластине неограниченно возрастают и происходит ее разрушение. Тот же подход для решения двумерного уравнения содержится в [54,55]. В работе [56] решалась стационарная задача
о прогибах ледяной поверхности под действием сосредоточенной нагрузки. Полученное аналитически решение в виде интеграла затем аппроксимировалось квадратурной формулой и рассчитывалось численно. В работе [58] рассмотрена задача о деформации ледяного покрова вьщвигающимся из воды вертикальным цилиндром. Уравнение колебаний пластин рассматривалось в осесимметричном случае, сосредоточенная сила задавалась также в виде 5-функции и было получено аналитическое решение.
Таким образом, можно сделать вывод о том, что численные методы решения таких задач практически не применялись. Если говорить о применении численных методов к решению двумерного уравнения поперечных колебаний тонких упругих пластин, то, например, в строительной механике в основном применялись методы конечных элементов, разностные же методы разработаны лишь для стационарной задачи [60,61], при этом рассматривались краевые задачи первого и второго рода для прямоугольных областей.
В 1 главы 3 приводится математическая постановка рассматриваемых экологических задач.
В 2 главы 3 разработан новый разностный метод решения поставленной задачи. С помощью замены переменных вместо исходного дифференциального уравнения, имеющего второй порядок по времени, рассматривается система уравнений с производными первого порядка по времени, для которой строится разностная схема. Двумерная область с криволинейным контуром, в которой рассматриваются исходная краевая задача, аппроксимируется ступенчатой областью, состоящей из прямоугольных ячеек, и рассматривается сетка с узлами в центрах этих ячеек. Разностная аппроксимация полученной системы уравнений строится балансным методом [62] с помощью интегрирования каждого из уравнений системы по площади ячеек. В результате получается неявная двухслойная разностная схема. Полученная разностная схема значительно проще для численной реализации, чем громоздкая трехслойная разностная схема, которая получилась бы при непосредственной разностной аппроксимации исходного дифференциального уравнения. Для решения системы разностных уравнений применялся метод расщепления по пространственным переменным.
Был создан программный комплекс на языке C++ , Эффективность разработанного численного метода была проверена путем сравнения результатов расчета с точным решением задачи о цилиндрическом изгибе прямоугольной пластины постоянной толщины под действием гармонической силы, заданной в виде плоской волны давления. Результаты расчетов показали высокое быстродействие алгоритма и хорошую точность.
В рассматриваемой задаче все коэффициенты в уравнении являются либо бесконечно дифференцируемыми функциями (толщина и цилиндрическая жесткость пластины), либо постоянными, однако, правая часть уравнения, описывающая действие сил на поверхности пластины, а также изгибающий момент и перерезывающая сила на контуре пластины могут быть разрывными функциями. В силу этого, решения исходной задачи рассматриваются как обобщенные решения из соответствующих функциональных пространств. В 3 главы 3 приведены условия гладкости входных данных, обеспечивающие, как следует из [63], существование и единственность обобщенного
решения исходной дифференциальной задачи из класса Х,2(0,Г;Я4(^))пЯ2(0,7,;12(О)).
Одной из самых трудных задач при обосновании численных методов решения задач математической физики является доказательство устойчивости построенных аппроксимаций и их сходимости к обобщенному решению исходной краевой задачи с оценкой скорости сходимости. Устойчивости разностных схем для эллиптических уравнений четвертого порядка по граничным условиям первого рода посвящены работы В.Б.Андреева [64,65]. Отметим также близкие работы Ю.И.Мокина, Р.ДЛазарова [66] и Л.С.Франка [67]. Среди работ в этой области следует отметить работы А.А.Злотника [68-70], в которых доказывается сходимость метода конечных элементов к обобщенному решению краевых задач и работы других авторов [71,72]. Отметим также работы, в которых доказывается сходимость метода Галеркина к обобщенным решениям для различных задач [73-81]. Что же касается конечно-разностного метода для уравнений четвертого порядка по пространственным переменным, то вопрос о сходимости разностных аппроксимаций к обобщенному решению исходной краевой задачи исследован лишь для стационарного бигармонического уравнения с однородными краевыми условиями [61].
В 3 главы 3 доказана теорема об устойчивости решения разностной задачи по входным данным.
В 4 главы 3 доказана теорема о сходимости разностных аппроксимаций, к обобщенному решению исходной краевой задачи для уравнения колебания тонких упругих пластин с оценкой скорости сходимости. Доказательство основано на методике разработанной автором при доказательстве сходимости разностных аппроксимаций третьей краевой задачи для эллиптического уравнения с переменными коэффициентами [140]. Эта методика, в свою очередь, основана на работах А.А.Самарского и его учеников [82-84], результаты которых были позже опубликованы в монографии [61].
В 5 главы 3 рассмотрены приложения задачи о колебаниях тонких упругих пластин к некоторым задачам экологии. С помощью созданного программного комплекса было проведено численное моделирование процесса распространения колебаний в ледяном покрове под действием динамических нагрузок, моделирующих движение одного и нескольких объектов по ледяной поверхности. При одновременном движении нескольких объектов друг за другом, в зависимости от расстояний между ними, наблюдалось либо гашение волн, либо их наложение с увеличением амплитуды. Проведено численное моделирование процесса распространения колебаний в ледяном покрове под действием динамической нагрузки, моделирующей всплытие подводного аппарата из-под ледяного покрова. Вычисляя компоненты тензора напряжений по закону Гука через изгибающие моменты, рассчитываемые в программе, можно численно исследовать прочность ледяного покрова при воздействии на него различных динамических техногенных нагрузок. Основные результаты по третьей главе опубликованы в работах [139-142,155-157].
Результаты работы по теме диссертации опубликованы в 24-х работах [135-158].
Автор глубоко благодарен В.Ф.Тишкину, Б.Н.Четверушкину и своему учителю Ф.П.Васильеву за постоянное внимание и поддержку в работе, а также Р.Н.Кузьмину за поддержку, оригинальные идеи и обсуждения физических моделей и приложений рассматриваемых задач.
Автор искренне благодарен А.А.Амосову, М.М.Потапову, В.Г.Звереву, А.П.Михайлову за обсуждения отдельных вопросов по теме диссертации, способствовавшие улучшению работы, а также Г.Г.Малинецкому за поддержку в работе.
Автор благодарен Н.П.Савенковой и С.В.Филипповой за плодотворную совместную работу на протяжении ряда лет.
Автор глубоко благодарен Е.Е.Мышецкой за совместную работу по созданию программного комплекса для моделирования лесных пожаров и проведение расчетов.
Автор благодарен Ал.А.Кулешову и В.В.Мымрину за совместную работу по созданию программного комплекса и проведению расчетов в задаче о колебаниях тонких упругих пластин.
Автор искренне благодарен Т.Г.Ермаковой и Е.Е.Мышецкой за большую помощь в подготовке текста диссертации и публикаций.
Алгоритм и методы численного решения задачи
Параболические уравнения (1.41) решаются стандартными сеточными методами, например, с помощью локально-одномерной разностной схемы или схемы переменных направлений (продольно-поперечной схемы) [66]. Находим w"+1, v"+1, которые на следующих этапах рассматриваемого временного шага уже не изменяются, а также находим С+1. Уравнения этапа 5 аппроксимируются явной разностной схемой. Уравнение этапа 6, где Jc определяется формулой (1.36), решается аналитически С«+1 _ «e-[ Polv"+,-V0N p/]A /(p5) где С" - значение, полученное после решения этапа 5. Затем на этапе 7 находим р" 1 и б 1.
Итак (л+1)-ый шаг численного решения системы (1.30)-(1.34) полностью описан. В целом численный метод имеет второй порядок аппроксимации по пространству и первый - по времени.
Как уже отмечалось во введении, программный комплекс для двумерной модели распространения тяжелых газов был создан под руководством автора в 1989-1990 годах в Институте атомной энергии им. И.В.Курчатова. Программный комплекс был реализован на языке FORTRAN в операционной среде UNIX на ЭВМ HP-1000.
С целью верификации модели было проведено сравнение результатов вычислительных экспериментов с результатами натурных экспериментов.
Наиболее полная серия опытов по распространению облака тяжелого газа была вьшолнена в Англии в 1982-1984 гг. на полигоне острова Торни [96-98]. Эксперименты финансировались группой промышленных компаний топливно-энергетического направления различных стран мира. Исследовалось распространение больших объемов предварительно подогретой до температуры окружающего воздуха смеси фреона-12 и азота. Контейнер с газом объемом 2000 м представлял собой разборный пластиковый цилиндр высотой 13 м. Его ячеистая конструкция предусматривала возможность быстрого раскрытия стенок и практически мгновенного объемного истечения смеси.
На полигоне острова Торни осуществлялся постоянный контроль состояния атмосферы (температуры и влажности воздуха, скорости ветра) аппаратурой, установленной на специальной 20-ти метровой метеорологической мачте. Мачта располагалась с наветренной стороны от испытательной площадки и была снабжена, помимо тепловых измерительных приборов, пятью чашечными и двумя звуковыми анемометрами для измерения распределения скорости ветра у поверхности земли. С наветренной стороны от источника на измерительной площадке площадью около 1 км2 были размешены более 30-ти измерительных мачт высотой 10 м, каждая из которых была оснащена четырьмя стандартными датчиками определения концентрации примеси тяжелого газа в воздухе. Все датчики предварительно перед каждым экспериментом тарировались. Несколько передвижных мачт, помимо четырех стандартных датчиков концентрации, были снабжены двумя звуковыми анемометрами и несколькими высокочастотными датчиками концентрации. В некоторых опытах использовалась термо-и видеотехника Института гидродинамики имени Кармана.
Основные результаты измерений представлены (см. [96-98]) на графиках изменения во времени концентрации примеси в различных точках площадки. Звуковые анемометры использовались для контроля динамического состояния возмущенной атмосферы. Экспериментальные данные [98] по вертикальному распределению концентрации примеси во фронте облака подтверждают гипотезу о высокой степени перемешивания газа в облаке и, следовательно, предположение о слабом изменении значений параметров течения в вертикальном направлении, сформулированное в п. 1.3 1, соответствует реальности.
На рис.1.2.1 изображены результаты численного эксперимента [135] по распространению облака фреона по ровной поверхности соответствующие натурному эксперименту. При сравнении результатов численных экспериментов с натурными (см. рис. 1.2.2) было получено хорошее совпадение (порядка 5%) по скорости распространения облака и несколько худшее совпадение (порядка 12%) по концентрации газа в облаке, что
объясняется погрешностью полуэмпирической формулы для разбавления облака окружающим воздухом.
С помощью созданного программного комплекса были также проведены тестовые расчеты по обтеканию облаком тяжелого газа препятствий различной формы и высоты (см.[135]), затеканию во впадины, распространению по наклонной поверхности, распространению при наличии и при отсутствии ветра.
Как отмечено во введении, авария под Уфой (м. Аша, Башкирия, СССР) 03.06.1989г. была связана с разрывом трубопровода, истечением 5000-10000 тонн пропан-бутаной смеси, образованием облака и его распространением по пересеченной местности с последующим воспламенением.
Впервые численное моделирование аварии под Уфой по модели (1.30)-(1.34) было проведено в 1990 году на ЭВМ НР-1000. Была промоделирована стадия гравитационного растекания облака. На рис, 1.3 представлен рельеф местности и граница зоны поражения, определенные экспертами в результате расследования аварии.
Математическое моделирование аварии 03.06.89 под Уфой
При изучении такого сложного природного явления как лесной пожар целесообразно рассматривать лес как многоярусную систему (см. [36,101,102]), состоящую из органического вещества - лесного горючего материала (ЛГМ) имеющего близкие геометрические и физико-химические параметры (истинная плотность, влагосодержание, теплотворная способность и т.д.) в каждом выделенном ярусе леса. 1) Первый ярус леса - сухая трава и мелкий кустарник высотой до 2 м. 2) Второй ярус леса - подрост - деревья высотой до 6 м. 3) Третий (верхний) ярус леса - кроны деревьев высотой до 22 м (в зависимости от типа древостоя).
Снизу под первым ярусом леса находится напочвенный слой, состоящий из мхов, опавших хвоинок, сучьев и листьев высотой до 15 см. Над верхним ярусом - приземной слой атмосферы.
Лесные горючие материалы содержат углерод, большое количество различных углеводородных соединений, воду, золу и т.д. При нагревании ЛГМ до температуры свыше примерно 370К еще до воспламенения (температура воспламенения ЛГМ, согласно [36], составляет примерно 500К) происходит пиролиз - термическое разложение ЛГМ с выделением летучих горючих (окись углерода СО, метан СН4, водород Нг и др.) и негорючих (двуокись углерода СОг, азот N2, пары воды НгО и др.) газов, дисперсных частиц сажи и с образованием твердого остатка - коксика и золы. Дисперсная сажа и коксик состоят практически из чистого углерода, который впоследствии воспламеняется.
Согласно [37] лесным пожаром называется явление неуправляемого многостадийного горения в открытом пространстве, на покрытой лесом площади, в рамках которого имеют место взаимосвязанные процессы конвективного и радиационного переноса энергии, нагревания, сушки и пиролиза ЛГМ, а также горение газообразных и конденсированных продуктов пиролиза ЛГМ.
Зона пожара, расположенная между сгоревшими и несгоревшими ЛГМ, называется фронтом пожара. Фронт пожара имеет переднюю и заднюю кромки. Во фронте пожара происходит наиболее сильное изменение параметров среды, а за задней кромкой фронта могут проходить процессы тления (беспламенного горения) твердых остатков пиролиза ЛГМ.
Механизм распространения лесного пожара в общих чертах таков: во фронте пожара происходит горение летучих (газы и дисперсная сажа) и конденсированных (коксик) продуктов пиролиза ЛГМ. В результате переноса (конвективного, диффузионного и излучением) тепла происходит нагрев находящейся перед фронтом органической массы, затем сушка, пиролиз и зажигание продуктов пиролиза и далее процесс повторяется. Приведем для наглядности схему из [36] процесса распространения лесного пожара.Однако, рассмотренная схема не учитывает важного механизма распространения пожара горящими частицами коры, шишек и древесины. Горящие частицы выносятся восходящими конвективными потоками в приземный слой атмосферы над фронтом пожара на высоту до нескольких сотен метров, затем, подхваченные ветром, переносятся на значительные расстояния от фронта пожара и, опускаясь, могут образовывать новые очаги пожара перед его передним фронтом. Часто эти пятнистые очаги пожара, сливаясь, образуют ложные фронты, которые длительное время могут двигаться впереди основного фронта.
В соответствии с приведенной выше вертикальной структурой леса будем рассматривать верховые и низовые пожары, распространяющиеся соответственно в верхнем или в нижних ярусах леса, а также сложные пожары, распространяющиеся одновременно в верхнем и в нижних ярусах леса. Напочвенные и торфяные пожары в настоящей работе рассматриваться не будут.
Ниже в 1 будет построена модель распространения лесного пожара в одном слое ЛГМ (ярусе леса), пригодная для моделирования пожара в любом из трех рассмотренных выше ярусах леса. Затем в 2 будет построена модель распространения пожара в верхнем ярусе леса с учетом переноса горящих частиц ветром в приземном слое атмосферы над лесом. И, наконец, в 3 будет построена модель совместного распространения верхового и низового лесных пожаров.
В предлагаемой ниже модели лес рассматривается как одноярусная двухфазная среда, состоящая из воздуха и летучих продуктов пиролиза и горения (газовоздушная или газовая фаза) и из лесных горючих материалов (ЛГМ) и твёрдых продуктов пиролиза и горения ЛГМ (твёрдая фаза). При построении физико-математической модели двухфазной гетерогенной смеси на основе методов механики сплошной среды такая смесь представляется (см. [103-105]) как двухкомпонентный континуум с взаимопроникающим движением фаз и с межфазным обменом массой, импульсом и энергией. Газовая фаза является многокомпонентной и состоит из горючих газов (СО, Нг, СЩ и др.), негорючих газов (СОг, N2 и др.), дисперсной сажи и окислителя (Ог). При этом предполагаем, что частицы дисперсной сажи движутся вместе с газовой фазой, и при сгорании сажи процесс теплообмена проходит быстро и можно рассматривать единую температуру газовой фазы. Твердая фаза также является многокомпонентной и состоит из ЛГМ, продукта пиролиза ЛГМ- коксика и золы.
Исходная система уравнений двухфазной модели состоит из двух подсистем -подсистемы уравнений газовой динамики в эйлеровых переменных, описывающих течение газовой фазы, и подсистемы эволюционных уравнений, описывающих изменение объемных долей компонентов твердой фазы и изменение температуры твердой фазы:
Двухъярусная модель лесных пожаров
Программный комплекс для моделирования лесных пожаров по построенным выше двумерным моделям был создан под руководством автора совместно с Е.Е.Мышецкой в Институте математического моделирования РАН на языке FORTRAN.
Поскольку натурных экспериментов по распространению реальных лесных пожаров с замером параметров среды во фронте пожара не проводилось, было проведено сравнение результатов расчетов по скорости распространения пожара для двумерной двухфазной модели с результатами расчетов по одномерной модели распространения пожаров, представленными в [50], которые согласуются с данными о распространении реальных лесных пожаров. Сравнение показало хорошее совпадение результатов (в пределах 5%).
Было проведено численное исследование чувствительности модели по параметрам. Например, в численных экспериментах определялось критическое значение влагосодержания ЛГМ - wh при котором лесной пожар не распространяется и прекращается. Этот параметр влияет на теплоемкость древесины cpi, определяемую формулой (2.14). При увеличении влагосодержания ЛГМ свыше 52% процесс горения в численных экспериментах прекращался. Это значение, очевидно, является для данного эксперимента тем самым критическим значением влагосодержания ЛГМ - Wj , при котором пожар не распространяется вследствие повышенного влагосодержания ЛГМ (см. [36,48]).
В численных экспериментах было также исследовано влияние на процесс распространения пожара количества дисперсной сажи в газовой фазе. При увеличении массовой концентрации дисперсной сажи на 15% температура во фронте пожара возрастала примерно на 60К, что согласуется с экспериментальными данными о значительном вкладе дисперсной сажи в энергетику пожара, обусловленном ее высокой теплотворной и излучательной способностью.
Применение метода расщепления по физическим процессам позволило численно исследовать влияние отдельных физических процессов на распространение пожара. Так, было установлено, что при распространении верховых пожаров определяющим процессом является конвективный перенос при наличии ветра, дующего со скоростью не менее 3 м/с (при меньшей скорости пожар не распространялся и затухал). Далее за конвективным переносом следует обмен лучистой энергии, а затем теплопроводность, что также согласуется с данными [36].
На рис.2.1.1-2.2.2 представлены результаты модельных расчетов по распространению лесного пожара в одном ярусе леса (слое ЛГМ) в продуваемом лесном массиве под действием ветра, дующего со скоростью РЬ=Зм/с в направлении противоположном оси Ох, в условиях неоднородного распределения запасов ЛГМ. Источник пожара первоначально имеет форму круга. На рис.2.1.1 а схематично изображена поляна прямоугольной формы и дорога, на которых отсутствует растительность). На рис.2.1.1b видно выгорание ЛГМ и образование температурного фронта пожара. Температура за задней кромкой фронта постепенно падает. Фронт пожара обтекает поляну и распространяется через узкую дорогу шириной 3 метра. На рис.2.1.2 изображена динамика фронта пожара рассчитанная по объемной доле ЛГМ. Фоновое зеленое поле - незатронутые пожаром ЛГМ (Фі=ФіІ=0) внутри - голубое поле полностью выгоревшие ЛГМ ф1 « 0.
На рис.2.2.1а изображена дорога шириной 20 метров. Фронт пожара не может преодолеть дорогу такой ширины и распространяется вдоль нее по направлению ветра (см. рис. 2.2.16). На рис.2.2.2 изображена динамика фронта пожара рассчитанная по объемной доле ЛГМ.
Было численно исследовано взаимодействие фронта пожара с препятствием -дорогой, развернутой под разными углами к набегающему фронту пожара. На рис.2.3, 2.4 фронт низового пожара не может преодолеть дорогу шириной 3 метра развернутую под углами 217 (рис.2.3) и 143 (рис.2.4) к фронту пожара при отсутствии ветра.
На рис.2.5 фронт низового пожара преодолевает дорогу шириной 3 метра развернутую под углом 90 к фронту пожара при отсутствии ветра, что объясняется более сильным прогревом более узкого сектора ЛГМ за дорогой.
На рис.2.7.1-2.9.2 приведены результаты модельных расчетов процесса распространения верхового пожара по трехфазной модели для модельных ситуаций, аналогичных тем, которые были рассмотрены в случае двухфазной модели. Скорость ветра направлена под углом 31 к оси Ох. На рис. 2.7.1, 2.7.2 пожар обтекает поляну. На рис.2.8.1, 2.8.2 пожар распространяется через дорогу шириной 20м при скорости ветра над лесным массивом Vw =10м/с. На рис.2.9.1, 2.9.2 пожар не может преодолеть дорогу шириной 20м при скорости ветра Vw =4м/с. Точки на рисунках 2.7.1-2.9.2 -осевшие горящие частицы. Некоторые из них выносятся ветром и оседают перед фронтом пожара и, при наличии там запасов ЛГМ, образуют новые очаги пожара, которые затем сливаются с основным фронтом. Таким образом, мы наблюдаем картину пятнистого лесного пожара, подобную картине реального лесного пожара.
Сходимость решения разностной задачи к обобщенному решению дифференциальной задачи
Для численного моделирования колебаний тонких упругих пластин был создан программный комплекс на языке C++.
Эффективность разработанного численного метода была проверена путем сравнения результатов численных расчетов с точным решением задачи о цилиндрическом изгибе прямоугольной пластины постоянной толщины под действием гармонической силы, заданной в виде плоской волны давления, распространяющейся под углом а (0 а л/2) КОСИЛ; f(x У- t) = % cos [k(x cos а + у sin а - vt)], где к - волновое число, v - фазовая скорость волны. Тогда решение исходной задачи для случая D=const будет иметь такой же вид W = W0 cos[A:(;tcosa + .у sin a - vt)]. (2.146) Подставляя это решение в уравнение (1.58) находим амплитуду Wo, затем из граничных условий (1.59), (1.60) определяем Мг и Nr , а из начальных условий (1.61), (1.62) находим ф и \\f. Результаты расчетов [142] показали высокое быстродействие алгоритма и хорошую (2-4 %) точность.
В дифференциальной задаче 1 р,а,Е,а - константы, толщина пластины h(x,y) -бесконечно дифференцируемая функция, такой же функцией является D = ЕЛ3/[12(1-а)], а функция F(x,y,t), описывающая действие сил на поверхности пластины и функции Mr,Nr, описывающие моменты и перерезывающие силы на краях пластины могут быть разрывными. В силу этого, решения задачи 1 будем рассматривать как обобщенные решения из соответствующих функциональных пространств.
Рассмотрим область gT=Qx[0,T]. Граница Г области Q, предполагается достаточно гладкой. Будем рассматривать следующие функциональные пространства [63,129,130]: CTO(Q), L2(Q), Ь2(Г), Нк(П), к = \А, с нормами ф (П), ЦфЦ , фя (П), пространства CtfO I Q)), С([0,Г];Я (О)), CtfO Z r», нормы в которых определяются как ЦфЦсгго пв)= шах \Шв» где В соответствующее функциональное пространство; пространства Lp(0,T;L2(Q)), Lp(0,T;L2(r)), р = \,2 т т L2(0,T;Hk(Q)), нормы в которых \Щ\к(0 Т;В) = ДОфЦ, , M w ) = ДфЛ , а также о о пространства (О ЩО.)), Я (0,Г;Яг(Г)).
Для оценки скорости сходимости разностных аппроксимаций необходимо, чтобы обобщенное решение исходной задачи обладало достаточно высокой гладкостью. Пусть выполнены следующие условия (см. [63]): - условия гладкости входных данных - условия согласования граничных и начальных значении Для исходной задачи 1 стандартным способом получим вспомогательную априорную оценку обобщенного решения, затем получим аналогичную оценку решения разностной задачи 2, из которой будет следовать устойчивость разностной схемы
Для простоты изложения, доказательство устойчивости, а затем сходимости разностных аппроксимаций проведем для прямоугольной сеточной области ЮА = І( і».Уі) = hNl,j = l,N2\, хотя техника доказательства полностью применима и для ступенчатых сеточных областей. Будем обозначать зависящие от времени разностные функции, определенные на a h, либо на Гн как [ф ]. Для введенных выше функциональных пространств рассмотрим их сеточные аналоги 12(соЛ),Я (соЛ), І2(ГА),С(Х2(ол)),С(Ял(а)А)),С(І2(ГА)) ;7(І2((ол)),і;7(І2(Гл)), /7=1,2, І2(ЯА(соА)) с нормами
