Содержание к диссертации
Введение
Глава 1, Конвективно-диффузионный перенос в задачах экологии для водной и воздушной среды 23
1.1 Математические модели конвективно-диффузионного переноса в водных и воздушных средах 24
1.11 . Физическое описание процессов конвекции и диффузии 24
1.1.2 Математическая модель температурного режима в водоемах .27
1.1.3 Модели распространения загрязнения в атмосфере 35
1.1.4 Формы записи операторов диффузионного и конвективного переноса 42
1.2 Разностная аппроксимация дифференциальной задачи конвекции -диффузии 46
1.2.1 Разностные схемы для стационарной задачи конвекции-диффузии 48
1.2.2 Эффективные способы аппроксимации нестационарного уравнения конвекции-диффузии 55
1.3 Выбор методов решения систем линейных алгебраических уравнений со специальными свойствами 62
1.3.1. Общая теория итерационных методов 63
1.3.2. Классические итерационные методы 68
1.3.3. Вариационные итерационные методы 72
1.3.4. Треугольные кососимметричные методы 76
Глава 2 Многосеточный метод для задач конвекции - диффузии 82
2.1. Этапы развития многосеточного метода 83
2.2. Описание многосеточного метода 89
2.2.1. Сглаживающая процедура 94
2.2.2. Грубо-сеточная коррекция 97
2.2.3. Функция интерполяции 103
2.2.4. Функция ограничения 105
2.2.5. Многосеточный алгоритм 107
2.3 Фурье-анализ многосеточного метода 113
2.3.1 Фурье-анализ для сеточных функций и операторов 114
2.3.2 Анализ на конечной области или анализ модельной задачи (МРА) 116
2.3.3 Локальный односеточный Фурье-анализ или анализ сглаживания 120
2.3.4 Двухсеточный локальный Фурье-анализ 137
2.3.5 Обобщенный анализ сглаживания 147
2.3.6 Упрощенный двухсеточный анализ 149
2.4. Сходимость модификаций многосеточного метода для задач конвекции - диффузии с преобладающей конвекцией 152
2.5 Численные исследования модификаций многосеточного метода для задач конвекции-диффузии с преобладающей конвекцией 157
Глава 3 Математическая модель температурного распределения в Азовском море 185
3.1 Пакет прикладных программ 185
3.2 Реализация математической модели температурного распределения Азовского моря 194
3.2.1. Гидрофизические характеристики Азовского моря 195
3.2.2 Описание модели температурного распределения Азовского моря 198
3.2.3 Численные эксперименты расчета температурного распределения в Азовском море 205
Глава 4 Математическая модель распространения радиоактивных примесей в воздушной среде в районе Волгодонской АЭС 227
4.1 Актуальность моделирования процессов распространения загрязняющих радиоактивных веществ в воздушной среде 227
4.2 Математическая постановка задачи 230
4.2.1 Обзор существующих математических моделей 231
4.2.2 Анализ входных метеорологических данных 237
4.2.3 Модель переноса радионуклидов в воздушной среде 239
4.3 Использование экономичных разностных схем с треугольным оператором для решения задачи 244
4.3.1 Исследование устойчивости треугольных кососимметричных схем 244
4.3.2 Численные эксперименты исследования свойств треугольных разностных схем 249
4.4 Численные эксперименты расчета распространения радионуклидов в воздушной среде в районе Волгодонской АЭС 252
Заключение 267
Список литературы
- Физическое описание процессов конвекции и диффузии
- Фурье-анализ для сеточных функций и операторов
- Реализация математической модели температурного распределения Азовского моря
- Использование экономичных разностных схем с треугольным оператором для решения задачи
Введение к работе
На современном этапе развития информационных технологий, включающего значительный прогресс средств переработки, передачи и хранения информации, проникновения их во все сферы жизни, математическое моделирование переживает очередную ступень своего формирования, «встраиваясь » в структуру информационного общества. Наличия информации, как таковой, зачастую недостаточно для анализа ситуации, принятия управленческих решений и контроля их исполнения. Необходимы адекватные и надежные способы обработки информации. История развития методологии математического моделирования показывает - именно она предоставляет такие способы, становясь, тем самым, ядром информационных технологий, процесса информатизации общества.
В общем перечне актуальных задач, решаемых с помощью математического моделирования, экологические проблемы занимают особое место. Увеличение антропогенного воздействия на окружающую среду, вызванное интенсивным использованием природных богатств, развитием материального производства, приводит к нарушению экологического равновесия как локально - в отдельных районах земного шара, так и глобально - в масштабах планеты в целом.
Естественным средством объективного анализа возникающих проблем являются методы, основанные на построении и совместном изучении математических моделей природных систем. Использование математического моделирования и проведение вычислительного эксперимента позволяют оценить все аспекты и последствия реализации любых проектов, связанных с воздействием на природную среду, как в перспективе, так и при возникновении всевозможных кризисных и экстремальных ситуаций. Важность и актуальность этого направления исследований усиливается тем обстоятельством, что его результаты имеют непосредственный практический выход в сферу социальных и экономических отношений современного общества.
Сущность методологии математического моделирования, предложенной в работе А.А.Самарского [92], состоит в замене исходного объекта его «образом» - математической моделью - и в дальнейшем изучении модели с помощью вычислительно-логических алгоритмов, реализуемых на современной компьютерной технике. Процесс математического моделирования можно условно разбить на три этапа «модель - алгоритм- программа». При этом следует уделять внимание всем трем составляющим триады. Необходимо отметить, что нынешнее состояние вычислительной техники, современных численных методов позволяют осуществлять моделирование объектов, поведение которых описывается весьма сложными математическими зависимостями, например, нелинейными системами дифференциальных или интегральных уравнений. Но сложные вычислительные алгоритмы обладают своими внутренними свойствами, которые далеко не всегда аналогичны, даже с точностью до ошибок аппроксимации, свойствам исходной математической модели. Это может приводить к появлению эффектов, имеющих чисто вычислительную природу. Поэтому важной задачей теории численных методов является разработка вычислительных алгоритмов, исключающих или сводящих к минимуму появление подобных ситуаций. Но пока такая теория отсутствует, большое значение имеет качественное исследование модели и ее возможного поведения, возможность найти ответы на три вопроса: что в данной модели может быть, что будет обязательно и чего не будет никогда [36]. Таким образом, проблема разработки адекватных моделей, особенно для описания процессов окружающей среды, и методов, их реализующих, остается весьма актуальной.
При решении многих экологических проблем необходимо исследовать процессы в движущихся средах, основными компонентами которых являются диффузионный перенос той или иной субстанции и конвективный перенос, обусловленный движением самой среды. При моделировании процессов полагают рассматриваемую среду сплошной, т.е. представляющую собой непрерывное распределение вещества и физических характеристик его состояния. Во многих случаях можно не оговаривать, о какой именно среде идет речь, поскольку и жидкость, и газ обладают схожими свойствами -сплошностью и текучестью. Сплошность - непрерывность распределения массы и физико-механических характеристик среды - является одним из основных свойств принятой модели жидкости или газа. Второе основное свойство - легкая подвижность или текучесть среды.
Обладая общими свойствами непрерывности и легкой подвижности, жидкости и газы отличаются друг от друга по физическим свойствам, связанным с различием их внутренней молекулярной структуры. Жидкости, в отличие от газов, можно считать малосжимаемыми, а иногда, в простейшей, достаточной для описания многих гидродинамических явлений схеме - просто несжимаемыми. В противоположность жидкостям, в газах межмолекулярные расстояния велики, а силы взаимодействия между молекулами сравнительно малы. В связи с этим, газы обладают свойством значительной по сравнению с жидкостями сжимаемостью. Однако, в случае слабых перепадов давлений, малых скоростей движения и отсутствия сколько-нибудь значительных нагревов и газ можно с достаточной степенью приближения рассматривать как несжимаемый.
В качестве базовых моделей многих процессов механики жидкости и газа выступают краевые задачи для стационарных и нестационарных уравнений конвекции - диффузии. К ним можно отнести задачи гидро- и газодинамики, распространение загрязнения и температурное распределение в водоемах и атмосфере, движение подземных вод, задачи магнитной гидродинамики и др. Среди экологических процессов, основой математических моделей которых является уравнение конвекции- диффузии, следует выделить распространение загрязняющих веществ в водной и воздушной средах, как особенно актуальные и востребованные задачи.
Несмотря на фактическую разницу в явлениях и описывающих их параметрах, коэффициентах уравнений, начальных и граничных условий каждая из указанных моделей имеет одинаковые члены, характеризующие два процесса - конвекцию (т.е. перенос субстанции за счет движения среды) и диффузию (т.е. вязкостные свойства среды). Изучая поведение уравнения конвекции-диффузии как модельной задачи, можно получить достаточно много информации о поведении решения конкретных практических задач.
При решении задач о переносе тепла с большими числами Пекле, о течениях жидкости, описываемых уравнениями Навье-Стокса с большими числами Рейнольдса или задач магнитной гидродинамики с большими числами Хартмана приходится часто сталкиваться с ситуацией, когда в уравнении конвекции - диффузии коэффициент при производной второго порядка мал по сравнению с коэффициентом при первой производной. Эти задачи ставятся на парабологиперболических или параболоэллиптических поверхностях и, таким образом, обнаруживают некоторые черты дифференциальных уравнений различных типов.
Известно [30], что когда параметр при старшей производной стремится к нулю и краевые условия не согласованы с правой частью уравнения, решения таких задач характеризуются появлением пограничных слоев, т.е. резкими изменениями решения в очень малой области расчета. Причем типы этих пограничных слоев могут быть различны (внутренний, приграничный, от начальных данных и т.д.) и зависят они как от граничных условий, так и от поля скоростей (коэффициентов при первых производных).
Трудности численного решения таких задач обусловлены их двойственной природой. Когда коэффициент при старшей производной становится достаточно малым, начальная эллиптическая задача ведет себя по существу как гиперболическая вне приграничных областей, в то время как диффузионный эффект наблюдается только в слоях. Однако, при стандартной численной аппроксимации подход к решению эллиптических и гиперболических задач различается.
Если задачи конвекции-диффузии с пограничными или внутренними слоями аппроксимируются с помощью центрально-разностной схемы, численное решение может быть «зашумлено» осцилляциями, так как соответствующий разностный оператор не является монотонным. Альтернативой является аппроксимация первых производных разностями «против потока» [90]. Тогда разностный оператор получается монотонным, и численное решение свободно от осцилляции, но порядок аппроксимации более низкий 0(h) и происходит «размазывание» пограничных слоев. Схемы с искусственной диффузией (streamline upwind scheme) [199] имеют более высокий порядок аппроксимации и меньше размазывают пограничные слоя, но недостатком их является то, что они не обеспечивают хороших оценок сходимости во всей области, и погрешность решения сильно возрастает в пограничных слоях, если сетка к ним не адаптирована [155]. Для решения проблем такого типа существуют так называемые глобально равномерно сходящиеся численные методы, т.е. методы, которые сходятся равномерно по малому параметру во всей области расчета. К таким методам относятся методы экспоненциальной подгонки [30], локального сгущения сетки [112], [195] и др.
При различных методах разностной аппроксимации дифференциального уравнения конвекции - диффузии получаем системы линейных алгебраических уравнений различного типа. В случае преобладающей конвекции использование противопотоковых схем приводит к системе линейных алгебраических уравнений с монотонной М-матрицей и сильному сглаживанию решения за счет появления в разностных уравнениях искусственной вязкости.
При использовании центрально-разностной аппроксимации и при желании сохранить монотонность разностного оператора приходиться накладывать ограничение на шаг сетки. Когда коэффициенты при первых производных достаточно велики, это ограничение становится существенным. Если эти ограничения не выполнены, то матрица полученной системы линейных алгебраических уравнений не будет иметь диагонального преобладания, и использование для решения такой системы базовых итерационных методов [97] приведет к большим трудностям, т.к. условие диагонального преобладания в исходной матрице является для этих методов условием их эффективной сходимости.
Следует также учесть, что для задач конвекции - диффузии кроме метода разностной аппроксимации весьма важной является начальная форма записи уравнения [97]. Существуют три формы записи оператора конвективного переноса, которые эквиваленты для дифференциального уровня в несжимаемых средах, но после аппроксимации приводят к различным формам разностньк уравнений, отличающихся по своим свойствам.
Таким образом, использование центрально-разностной аппроксимации при решении задач конвекции - диффузии с преобладающей конвекцией сохраняет характер поведения решения, но в результате получается система линейных алгебраических уравнений с несимметричной матрицей, не имеющей диагонального преобладания. Для этого типа задач большинство классических методов либо вообще не работают, либо обладают очень медленной скоростью сходимости. Поэтому так актуальна проблема создания эффективных численных методов для решения задач конвекции-диффузии с преобладающей конвекцией.
В настоящее время для решения задач линейной алгебры существует множество различных численных методов, которые непрерывно усовершенствуются и модифицируются. Активно разрабатываются новые методы. В результате оказывается, что значительная часть созданных методов имеет право на существование, обладая своей областью применимости. При решении конкретной задачи важно выбрать наиболее подходящий для рассматриваемого класса задач метод из множества допустимых методов решения данной задачи. Этот метод, очевидно, должен обладать наилучшими характеристиками, такими как минимум времени решения задачи на компьютере (или минимум числа арифметических и логических операций при нахождении решения), вычислительной устойчивостью, т. е. устойчивостью по отношению к ошибкам округления и др. При выборе метода решения задач конвекции-диффузии необходимо учитывать перечисленные выше особенности рассматриваемого класса задач.
Одним из критериев выбора алгоритма, используемого при численном моделировании той или иной физической задачи, является объем вычислительной работы, необходимый для его реализации. Существует правило, что этот объем должен быть пропорционален реальным физическим изменениям, происходящим в моделируемой системе. Если алгоритм требует большого количества тяжелой вычислительной работы для расчета слабого эффекта или очень медленного физического процесса, то от такого «затратного» алгоритма следует, отказаться, выбрав более эффективный.
Примером «затратных» алгоритмов являются обычные итерационные методы для решения алгебраических уравнений, возникающих при численном решении уравнений в частных производных или интегро-дифференциальных уравнений. Так, практически единственным, но наиболее существенным недостатком методов Якоби и Гаусса-Зейделя, используемых для решения эллиптических задач методом сеток, является их низкая скорость сходимости. Другим примером могут служить решения нестационарных задач, с шагом по времени (выбор которого диктуется условиями устойчивости) много меньшим масштаба реального изменения решения. То есть, в общем случае, «затратным» можно назвать такой алгоритм, который требует использования очень подробных сеток, там, где на большей части расчетной области величина шага по пространству или по времени много меньше, чем реальный масштаб изменения решения.
В этом случае эффективным решением проблемы является использование многосеточного алгоритма, который позволит преодолеть главную трудность, возникающую при решении такого рода задачи - ее «жесткость». Жесткость задачи заключается в существовании нескольких компонент решения, которые имеют разный масштаб и конфликтуют друг с другом. Например, гладкие компоненты, которые можно эффективно аппроксимировать на грубых сетках, но которые плохо сходятся на мелких сетках, конфликтуют с высокочастотными компонентами, которые необходимо аппроксимировать с помощью мелких сеток. Используя несколько уровней дискретизации, многосеточный алгоритм решает конфликты такого рода, позволяя достигать большой эффективности, путем снижения объема вычислений, необходимых для получения численного решения.
Благодаря вышеуказанным свойствам многосеточный метод (Multi-Grid Method - MGM) стал в последние годы одним из эффективных и довольно универсальных итерационных методов решения задач. Он принадлежит к классу быстро сходящихся итерационных методов, является оптимальным по числу арифметических операций для достижения точности, согласованной с порядком сходимости. Скорость сходимости многосеточного метода всегда независима от числа неизвестных в системе, полученной в результате аппроксимации дифференциального уравнения, то есть многосеточный метод обладает неулучшаемой оценкой сходимости. Другая особенность метода - то, что он является своего рода шаблоном. Не существует строго определенного многосеточного алгоритма, применимого ко всем краевым задачам. Многосеточный метод устанавливает лишь структуру алгоритма, эффективность которого во многом зависит от адаптации его компонент к конкретной задаче.
Обладая высокой эффективностью, многосеточные методы допускают наиболее естественное распараллеливание и векторизацию приложений, что позволяет отнести их к наиболее перспективному и быстро развивающемуся разделу высокопроизводительных алгоритмов.
Многосеточный метод может применяться к задачам, рассматриваемым в областях произвольной формы и с различными граничными условиями. MGM может быть использован при решении сложных, несимметричных и нелинейных систем уравнений. Он может применяться для решения нестационарных параболических уравнений. В последние несколько лет ведутся активные исследования использования многосеточного метода для гиперболических уравнений.
В настоящее время многосеточные алгоритмы эффективно применяются для решения задач динамики плазмы и гидродинамики, расчета собственных значений и собственных функций дифференциальных операторов, для расчета нейтронных полей в ядерном энергетическом реакторе, для решения задач теории упругости, а так же в задачах обтекания тел достаточно сложной формы.
Обобщая вышеизложенное, сформулируем основную цель и задачи исследования.
Целью работы является создание эффективных алгоритмов реализации математических моделей конвективно-диффузионного переноса в движущихся средах и их использование в конкретных экологических задачах водной и воздушной среды.
Для достижения поставленной цели было необходимо решить следующие задачи:
• разработать математическую модель динамики температурного распределения в Азовском море, определить вычислительные алгоритмы реализации созданной модели на высокопроизводительных вычислительных системах, разработать программный комплекс и провести численные эксперименты
• разработать математическую модель распространения радиоактивных примесей в воздушной среде в районе Волгодонской АЭС
• разработать и исследовать эффективные разностные схемы для решения нестационарного уравнения конвекции-диффузии с преобладающей конвекцией
• разработать модификации многосеточного метода для решения сильно несимметричных СЛАУ, возникающих после центрально - разностной аппроксимации стационарного уравнения конвекции - диффузии с преобладающей конвекцией
• создать программный комплекс, реализующий разработанные алгоритмы для математической модели распространения радиоактивных примесей в воздушной среде в районе Волгодонской АЭС и провести численные эксперименты
Представим краткое описание содержания диссертационного исследования.
Во введении сформулирована цель и задачи работы, обоснованы актуальность исследуемой темы, кратко описано содержание работы. Обсуждаются полученные в диссертационном исследовании результаты.
Первая глава носит вспомогательный характер и состоит из трех разделов.
Первый раздел посвящен общему описанию математических моделей конвективно-диффузионного переноса в водных и воздушных средах. Особое внимание уделено моделям температурного распределения в водоеме и распространения примесей в атмосфере. Приведены три формы записи операторов диффузионного и конвективного переноса в используемом в данных моделях уравнении конвекции - диффузии. Обсуждаются особенности рассматриваемых процессов, влияние соотношения между процессами конвекции и диффузии на математическую постановку задачи и на дальнейший выбор методов решения.
Во втором разделе рассмотрена разностная аппроксимация задач конвекции - диффузии. Приведены разностные схемы для стационарной и нестационарной задачи. Описано влияние формы записи конвективных членов уравнения конвекции - диффузии на свойства матрицы, получаемой после разностной аппроксимации дифференциальной задачи.
В третьем разделе представлены методы решения сильно несимметричных систем линейных алгебраических уравнений, возникающих после центрально-разностной аппроксимации уравнения конвекции -диффузии. Представлен обзор итерационных методов решения сильно несимметричных СЛАУ. Приведены некоторые сведения из теории матриц и функционального анализа, необходимые в дальнейшем исследовании. Рассмотрены классические итерационные и вариационные методы решения СЛАУ.
Если при конструировании методов для решения задач с сильно несимметричными СЛАУ особое внимание уделять учету структуры оператора решаемой задачи, это позволяет строить и применять специальные итерационные методы, которые обладают более высокой скоростью сходимости, чем методы из общей теории. К таким методам относятся треугольные кососимметричные методы, впервые предложенные в работах Л.А. Крукиера [40]. Эффективность их применения достигается особым выбором операторов и итерационных параметров.
Треугольные кососимметричные методы могут быть использованы для решения несимметричных систем линейных алгебраических уравнений, получаемых после центрально разностной аппроксимации уравнения конвекции - диффузии с преобладающей конвекцией, в которых отсутствует диагональное преобладание. Однако при значительном доминировании конвекции в исходном уравнении мы приходим к сильно несимметричным системам, и треугольные кососимметричные методы, хотя и сходятся для этих случаев, но скорость сходимости невысока. Причем отмечено, что эти методы быстро сходятся на первых итерациях, замедляясь в дальнейшем. Это свойство методов называется сглаживающим и объясняется быстрым подавлением высокочастотных гармоник ошибки, но гораздо более медленным воздействием на низкочастотные составляющие. Именно этим свойством должны обладать итерационные методы, используемые в качестве сглаживателей в многосеточном методе. В работе автором диссертации впервые предложено использовать треугольные кососимметричные методы в качестве одной из компонент многосеточного метода.
Вторая глава посвящена описанию и исследованию модификации многосеточного метода для решения сильно несимметричных систем линейных алгебраических уравнений, получаемых после разностной аппроксимации стационарного уравнения конвекции - диффузии с преобладающей конвекцией.
Глава состоит из пяти разделов.
В первых двух разделах приведены общее описание и этапы развития многосеточного метода, основная идея которого принадлежит Р.П. Федоренко. Многосеточный метод является оптимальным по числу арифметических операций для достижения точности, согласованной с порядком сходимости. Другая особенность метода - то, что он является своего рода шаблоном. Многосеточный метод устанавливает лишь структуру алгоритма, эффективность которого во многом зависит от адаптации его компонент к конкретной задаче. Значительный вклад в теорию многосеточного метода внесли Р.П. Федоренко, Н.С. Бахвалов, ГЛ. Астраханцев, A. Brandt, W. Hackbusch, P.O. Frederikson, P. Wesseling, P. Sonneveld, U. Trottenberg, B.B. Шайдуров, а также P. Bastian, I. Bey, D. Braess, J.H. Bramble, S.C. Brenner, M. Dryja, H.C. Elman, J.E. Pasciak, R.A. Nkolaides, A. Reusken, R. Sarazin, R. Stevenson, S. P. Vanka, J, Wang, G.Wittum, J. Xu, H. Yserentant, X. Zhang, M.A. Ольшанский и многие другие.
Многосеточный алгоритм позволяет значительно повысить эффективность базового итерационного метода, комбинируя обычный итерационный процесс с приемом, называемым грубосеточной коррекцией -последовательным использованием в вычислениях более грубых сеток. Одна итерация метода включает в себя 4 наиболее важных этапа: сглаживание (smoothing), ограничение, проекция (restriction), продолжение, интерполяция (prolongation) и грубосеточная коррекция (coarse greed correction).
Сглаживающая процедура является важной компонентой многосеточного алгоритма, наиболее зависимой от решаемой задачи. В диссертационной работе приводится описание существующих сглаживателей многосеточного метода -методы Якоби, Гаусса - Зейделя, метод Якоби с весами. В работе впервые предложена новая модификация многосеточного метода, в которой используются специальные сглаживатели для решения сильно несимметричных СЛАУ, получаемых после центрально-разностной аппроксимации уравнения конвекции - диффузии.
Третий раздел главы посвящен Фурье - анализу многосеточного метода, который является важным инструментальным средством для получения количественных оценок сходимости и оптимизации различных компонент многосеточного метода. Основная идея Фурье-анализа, изложенная в этом разделе, состоит в том, чтобы представить ошибку или невязку в виде суммы некоторых периодических функций, называемых компонентами Фурье или гармониками. При этом появляется возможность оценить воздействие составляющих многосеточного метода на каждый компонент Фурье -разложения.
В диссертации проведен односеточный локальный Фурье-анализ или анализ сглаживания и двухсеточный Фурье-анализ предложенной модификации многосеточного метода, в которой в качестве сглаживателей используются методы из класса треугольных кососимметричных методов. При проведении односеточного анализа основное внимание в многосеточном цикле уделяется процедуре сглаживания, а влиянием грубо-сеточной коррекции пренебрегают или используют "идеальный" оператор грубо-сеточной коррекции. При проведении Фурье-анализа сглаживания важным моментом является вычисление коэффициента сглаживания, который показывает, насколько эффективным сглаживателем является рассматриваемый метод.
В работе проведен двухсеточный локальный Фурье-анализ (LFA). При проведении LFA были получены коэффициенты сглаживания для треугольных кососимметричных методов. Для оценки сходимости двухсеточного метода с треугольными кососимметричными сглаживателями были вычислены коэффициенты асимптотической сходимости двухсеточного метода. Проведено сравнение результатов Фурье-анализа многосеточного метода, в котором в качестве сглаживателей выбирались треугольные кососимметричные сглаживатели с многосеточным методом со стандартными сглаживателями (методом Гаусса-Зейделя и методом Якоби).
С помощью Фурье-анализа в работе были проведены исследования треугольных кососимметричных методов, показано, что они обладают хорошим сглаживающим свойством. Для эффективной работы метода в качестве сглаживателя коэффициент сглаживания ///яс должен быть меньше единицы.
Чем меньше коэффициент сглаживания, тем быстрее метод подавляет высокочастотные компоненты ошибки.
Четвертый раздел главы посвящен исследованию предложенных модификаций многосеточного метода для сильно несимметричных СЛАУ и доказательству их сходимости. В качестве исследуемой задачи рассматривается стационарное уравнение конвекции - диффузии.
Проведено исследование трех сглаживателей из рассматриваемого класса треугольных кососимметричных методов - ТКМ, ТКМ1 и ТКМ2. В диссертационной работе эти методы использованы в качестве сглаживателей многосеточного метода, который можно рассматривать как своего рода ускоряющую процедуру треугольных кососимметричных методов. Доказаны теоремы сходимости предложенной модификации многосеточного метода.
В пятом разделе главы приведены результаты численных экспериментов использования многосеточного метода для задач конвекции - диффузии с преобладающей конвекцией. Приведены результаты исследования зависимости эффективности многосеточного метода от числа сглаживающих итераций и количества уровней. Численные исследования подтвердили полученные ранее с помощью Фурье - анализа результаты об оптимальном количестве сглаживающих итераций в многосеточном методе для решения сильно несимметричных систем.
Были проведены вычислительные эксперименты для четырех модельных задач с различными векторами скорости движения среды на сетках разной размерности - от 32x32 до 512x512. Исследовалось поведение метода в зависимости от числа Пекле.
Проведенные численные исследования модификаций многосеточного метода подтвердили полученные ранее теоретические результаты и показали, что предложенная модификация многосеточного метода со сглаживателями ТКМ, ТКМ1, ТКМ2 эффективна для решения задач конвекции-диффузии с преобладающей конвекцией.
Третья глава содержит результаты моделирования динамики температурного распределения в Азовском море.
Первый раздел главы посвящен описанию программного комплекса для реализации гидрофизических моделей процессов переноса в водной и воздушной среде. В рамках диссертационной работы был разработан ряд модулей программного комплекса, реализующих предложенные подходы для математического моделирования экологических задач, основой которых являются процессы конвективно-диффузионного переноса в водной и воздушной среде. Разрабатываемый в Южно-Российском региональном центре информатизации Ростовского госуниверситета программный продукт предназначен для решения задач экологии водной и воздушной сред и состоит, как и любой стандартный пакет прикладных программ, из функциональной и сервисной частей. Программный комплекс ориентирован на расчет движения среды (гидро-, газодинамика), решение задач конвекции-диффузии, что реализуется с помощью функциональных модулей. Сервисная часть включает модули генератора карт местности, БД с системой управления, модуль визуализации результатов расчета, генератор отчетов. В рамках диссертационного исследования созданы расчетные модули, реализующие многосеточный метод и треугольные кососимметричные схемы для решения задач конвекции - диффузии.
Второй раздел главы посвящен описанию гидрофизических характеристик Азовского моря - природного объекта исследования данной главы. Азовское море имеет важное хозяйственное значение для южных регионов России. Вследствие непродуманных действий по эксплуатации ресурсов Азовского моря за последние десятилетия, в настоящее время экосистема Азовского моря выведена из состояния равновесия. Поэтому особую актуальность приобретает работа по созданию инструмента исследования и прогнозирования состояния акватории моря. Важная гидрофизическая характеристика водоема - температурное распределение -является одной из компонент в задачах теории климата, прогноза погоды, расчета энергообмена и др. В работе представлена математическая модель температурного режима Азовского моря, в которой особое внимание уделено определению функции притоков - оттоков тепла. Проведен анализ компонентов теплового баланса, определены его составляющие, оказывающие наиболее существенно влияние на температурный режим водоема в рассматриваемый неледоставный период (апрель - октябрь). При численной реализации предложенной модели использовалась противопотоковая разностная схема, поскольку для данной задачи, в которой процесс конвекции не является преобладающим, это схема наиболее эффективна. Задача решалась на высокопроизводительных вычислительных системах с использованием пакета распараллеленных итерационных методов Aztec.
Проведен ряд вычислительных экспериментов, в которых исследовались зависимости функции F(x,y,z,t) притоков-оттоков тепла от пространственных координат и от времени, влияние поля скоростей течений моря, способа задания начального распределения температуры воды на результаты расчета динамики температурного распределения. Результаты численных экспериментов сравнивались с натурными наблюдениями, погрешность составила 5%-10% для разных вычислительных экспериментов.
На основе анализа результатов проведенных вычислительных экспериментов и сравнения их с натурными данными были сделаны выводы, которые можно использовать в качестве рекомендаций для практического использования созданной модели температурного распределения Азовского моря:
1. Способ задания начального распределения температуры воды не оказывает существенного влияния на расчет динамики температурного распределения в водоеме.
2. Суммарная солнечная радиация оказывает наиболее существенное влияние на процесс распределения температуры в водоеме в рассматриваемый период времени (апрель - октябрь).
3. Характерные ветровые ситуации оказывают определенное влияние на температурное распределение в Азовском море.
Четвертая глава содержит результаты численных экспериментов для математической модели распространения радионуклидных примесей в воздушной среде в районе Волгодонской АЭС, проведенных на основе созданного программного комплекса, реализующего предложенные подходы к решению задач конвекции - диффузии с преобладающей конвекцией.
Первый раздел главы посвящен анализу актуальности моделирования процессов распространения загрязняющих радиоактивных веществ в воздушной среде. Подчеркивается особая важность такого рода исследований для регионов юга России в связи с принятым в июле 2006 года решением о продолжении строительства второго энергоблока Волгодонской АЭС.
Второй раздел главы посвящен моделированию процесса распространения примесей в воздушной среде. Исследуется модель распространения радионуклидов в районе Волгодонской АЭС. Полученные в результате моделирования и вычислительных экспериментов данные дают возможность с их помощью анализировать экологическую безопасность штатного и нештатного режимов работы Волгодонской АЭС.
В заключении приведены основные результаты, полученные в диссертационной работе.
Научная новизна диссертационного исследования заключается в разработке математических моделей температурного распределения в Азовском море и распространения радиоактивных примесей в воздушной среде в районе Волгодонской АЭС. В работе предложен новый класс треугольных кососимметричных разностных схем для решения динамических задач конвекции-диффузии с преобладающей конвекцией, модификация многосеточного метода решения сильно несимметричных систем линейных алгебраических уравнений (СЛАУ), где в качестве сглаживателя используется итерационный метод из класса треугольных кососимметричных методов.
Достоверность полученных результатов обеспечивается строгим математическим обоснованием предложенных методов и алгоритмов, качественным совпадением результатов вычислительных экспериментов с натурными наблюдениями.
К ЗАЩИТЕ ПРЕДСТАВЛЕНЫ СЛЕДУЮЩИЕ РЕЗУЛЬТАТЫ:
1. Разработана математическая модель динамики температурного распределения в Азовском море. Определены вычислительные алгоритмы реализации модели на высокопроизводительных вычислительных системах. Создан программный комплекс, проведены численные эксперименты динамики температурного распределения в Азовском море.
2. Предложен и исследован новый класс условно устойчивых и абсолютно устойчивых треугольных разностных схем решения задач конвекции-диффузии с преобладающей конвекцией.
3. Предложена, теоретически и численно исследована модификация многосеточного метода решения сильно несимметричных СЛАУ, полученных после разностной аппроксимации уравнения конвекции-диффузии с преобладающей конвекцией. Доказаны теоремы сходимости предложенной модификации многосеточного метода. На основе Фурье-анализа рассмотренных модификаций многосеточного метода исследованы способы выбора различных сглаживателей из класса треугольных кососимметричных итерационных методов. Показана зависимость скорости сходимости многосеточного метода от количества сглаживающих итераций.
4. Разработана математическая модель переноса радиоактивных примесей в воздушной среде в районе Волгодонской АЭС. Создан программный комплекс, реализующий предложенные алгоритмы. Проведены вычислительные эксперименты на основе реализованных математических моделей переноса радионуклидов в районе Волгодонской АЭС.
Физическое описание процессов конвекции и диффузии
Диффузией называется перемещение частиц в направлении убывания их концентрации, обусловленное тепловым движением. Под частицами в данном случае понимается наименьшая структурная единица вещества в рассматриваемом процессе (молекула, ион, атом). Диффузия приводит к выравниванию концентрации частиц диффундирующего вещества и равномерному заполнению частицами объема, если только неравномерное распределение не поддерживается какими-либо внешними силами, действующими на частицы. Диффузия имеет место в газах, жидкостях и твердых телах, причем, диффундировать могут как растворенные в веществе посторонние частицы, так и частицы самого вещества. Последнее явление называется самодиффузией.
Скорость диффузии определяется величиной коэффициента диффузии, который возрастает с повышением температуры, когда тепловое движение частиц становится более быстрым. С наибольшей скоростью диффузия протекает в газах. Скорость диффузии в газах определяется как скоростью теплового движения молекул, так и длиной их свободного пробега, т.е. средней длиной тех прямолинейных отрезков пути, которые проходят молекулы газа от столкновения к столкновению.
Первый закон Фика определяет количество вещества, диффундирующего в направлении убывания концентрации. Если градиент концентрации С вдоль направлении х равен dcidx, то этот закон определяет для массы вещества dm, диффундирующего за время dt через площадку S, перпендикулярную направлению х, следующую зависимость: dm = -D-S—dl. (1.1.1) dx
Знак минус показывает, что диффузия происходит в сторону уменьшения концентрации; D - коэффициент диффузии, численно измеряемый массой вещества, диффундирующего через единичную площадку за время t = \ при градиенте концентрации, равном / (имеет размерность см2 -сек 1). Коэффициент диффузии определяет скорость процесса и зависит от природы частиц и состояния диффундирующего вещества и растворителя (в растворах).
Из выражения (1.1.1) при условии постоянства коэффициента диффузии получается соотношение, называемое вторым законом Фика: dt dx = D -. (1.1.2)
Из соотношения (1.1.2) в частном случае, когда dcldt-Q (стационарный поток) следует, что концентрация диффундирующего вещества должна линейно уменьшаться вдоль направления диффузионного потока. Формулы (І.1.1) и (1.1.2) приведены для случая, когда диффузия протекает в одном направлении, однако они обобщаются и на случай трех измерений.
В случае трехмерной диффузии изменение концентрации с течением времени при постоянной температуре и отсутствии внешних сил описывается дифференциальным уравнением диффузии в частных производных: dt дх\ дх) ду 1\ М dz Если D не зависит от концентрации С, то уравнение приводится к виду: BD C+D C+D C dt дхг ду1 dz1 А. Эйнштейном было установлена зависимость коэффициента диффузии частиц в вязкой среде: D-q-k, где к -постоянная Больцмана, Т - абсолютная температура, q - коэффициент, называемый подвижностью частиц и численно равный скорости, с которой частица движется в данной среде под действием силы, равной 1. Обратная величина f=\lq носит название коэффициента трения. По закону Стокса, в случае шарообразных частиц справедливо выражение / = бпгу, где г -радиус шарообразной частицы, а ц- коэффициент вязкости среды. Следовательно, коэффициент диффузии шарообразных частиц в вязкой среде D = kT/6mjr. Напомним, что под частицами в данном случае понимается наименьшая структурная единица вещества в рассматриваемом процессе (молекула, ион, атом). Закон Стокса применим только к частицам шарообразной формы. В случае частиц произвольной формы учитывается еще и асимметрия частиц, характеризуемая отношением /: /fl, где /, /0 - коэффициенты трения, соответственно, для частиц произвольной формы и эквивалентных им (с равной массой) шарообразных частиц.
Конвекция - это перенос теплоты, массы и других физических величин в жидкостях, газах или сыпучих средах потоками вещества.
Различают следующие виды конвекции: конвекцию в атмосфере - более или менее упорядоченный воздухообмен между верхними и нижними слоями тропосферы; конвекцию тепла - перенос теплоты, обусловленный перемещением масс жидкости или газа, под влиянием различия температур в разных частях жидкости или газа; конвекцию электрического разряда - перемещение разряда вместе с телом, внутри которого или на поверхности которого распределен этот заряд; конвекцию в море - водообмен между верхними и нижними слоями в океанах и морях, вызванный изменениями плотности морской воды вследствие изменения ее температуры или солености.
Фурье-анализ для сеточных функций и операторов
Применение Фурье-анализа к многосеточным методам было предложено Брандтом [124]. Основная идея Фурье-анализа состоит в том, что ошибка (или невязка) раскладывается в сумму некоторых периодических функций, называемых компонентами Фурье. При этом появляется возможность оценить воздействие различных многосеточных компонент на каждый компонент Фурье [129], [217], [218], [221], [235].
Существует два типа Фурье-анализа, применимого к многосеточному методу. Первый - точный анализ, известный как анализ модельной задачи (МРА), который может быть применен только к некоторым модельным ситуациям, таким как прямоугольная область и постоянные коэффициенты операторов [238]. Второй тип анализа называется локальным анализом Фурье (LFA) или способом локального анализа (LMA) [124], [221], [238]. В LFA, основные дискретные операторы с постоянными коэффициентами считаются формально расширенными на бесконечную сетку. Следовательно, граничными условиями пренебрегают. Согласно общим предположениям, любой дискретный оператор, нелинейный, с непостоянными коэффициентами может быть локально линеаризован и может быть локально заменен (замораживая коэффициенты) оператором с постоянными коэффициентами [221]. Это демонстрирует большой диапазон применимости LFA и его локальную природу. В большинстве случаев приходится переходить именно к LFA особенно для оценки несимметричных задач, исследование которых является актуальным в реальной жизни.
Стандартной стратегией получения быстрого многосеточного метода является рассмотрение модельных уравнений, которые описывают большой набор задач с качественно подобным поведением. В этом случае Фурье-анализ применяется для того, чтобы выбрать базовый итерационный метод с хорошим сглаживанием ошибки (анализ сглаживания), и грубо-сеточную коррекцию (двухсеточный анализ). Эта стратегия очень успешна, пока двухсеточный метод является хорошим приближением многосеточного метода. Анализ многосеточного метода нельзя заменить рассмотрением двухсеточнго метода при исследовании работы V- и W-циклов, пред- и постсглаживания, различных сглаживателеи на различных сетках, или различной дискретизации на различных сетках. Таким образом, важно перейти от Фурье-анализа с двумя сетками к анализу с -сетками. При этом можно получить более глубокое понимание трудностей, связанных с грубо-сеточной коррекцией.
Существует два принципа для адаптации многосеточных компонентов: оптимальность и ошибкоустойчивость. Оптимальность означает, что многосеточный метод является эффективным, но строго специализированным для решения одной специфической задачи, тогда как ошибкоустойчивость говорит о том, многосеточный метод хорошо работает для большого класса задач. Эти принципы имеют и преимущества, и недостатки. Оптимизированный многосеточный метод используется тогда, когда специализированная задача должна быть решена не один раз. Ошибкоустойчивость появляется в том случае, когда должны решаться различные задачи с подобной структурой, например, в научном программном обеспечении для решения дифференциальных уравнений в частных поризводных. Трудно выбрать универсальные многосеточные компоненты для различных классов задач.
Главной целью анализа является разложение ошибки или невязки Vh:=u;-uh, іі АХ-Ци,, (2.3.1) в ряд Фурье после /-ой итерации многосеточного метода с использованием it-сеток. Здесь, uh обозначает точное дискретное решение (2.2.1) и и\ его дискретное приближение после г -ой итерации многосеточного метода с использованием -сеток. Далее каждая компонента Фурье может быть исследована и дана оценка ее влияния на эффективность многосеточного метода. Основные сведения об Фурье-анализе для сеточных функций и операторов могут быть найдены в [124], [132], [162], [170], [218], [221], [235], [238].
Реализация математической модели температурного распределения Азовского моря
Важная гидрофизическая характеристика водоема - температурное распределение, которое является одной из компонент в задачах теории климата, прогноза погоды, расчета энергообмена и др. В качестве природного объекта исследования рассматривается Азовское море.
Азовское море расположено между 3517 и 3717 с.ш. и 3339 и 3918 в.д. Климат Азовского моря относится к континентальному климату умеренных широт. Для него характерна умеренно мягкая, короткая зима и теплое продолжительное лето.
Азовское море, несмотря на свои небольшие размеры, имеет важное хозяйственное значение для южных регионов России [16]. Оно отличается высокой рыбопродуктивностью, которая обеспечивается мелководностью моря, низкой соленостью морской воды, высокой концентрацией биогенных веществ. Несмотря на уменьшение уловов ценных видов рыб за последние десятилетия, море и в настоящее время не потеряло своего рыбохозяйственного значения.
Теплый умеренно-влажный климат, песчаные дно и пляжи, месторождения лечебных грязей и минеральных вод создают благоприятные условия для использования моря в лечебно-курортных целях. Строительство Волго-Донского канала превратило Азовское море в важное звено транспортного сообщения между портами Волжско-Камского и Азово-Черноморского бассейнов.
Мелководность и внутриконтинентальное расположение моря обусловливают большую временную изменчивость гидрологических и гидрохимических характеристик.
Гидролого-гидрохимический режим моря формируется под воздействием речного стока, водообмена с Черным морем и климатических факторов. Существенную роль играет мелководность моря. Благодаря большому количеству поступающей солнечной радиации Азовское море имеет довольно высокую среднюю годовую температуру воды +11,5.
В настоящее время в акватории моря действуют 11 гидрометеостанций и 9 гидропостов. Первоначально гидрометеорологические и гидрохимические исследования проводились с целью увеличения продуктивности рыбного хозяйства и развития мореплавания. Но с увеличением антропогенного воздействия человека изменяется и направленность гидрометеорологических и гидрохимических исследований. Интенсивное развитие народного хозяйства в бассейне Азовского моря в 60-х годах (строительство водохранилищ, орошаемое земледелие, развитие водоемных отраслей промышленности и т.д.) привело к заметному нарушению стабильности гидрологического, гидрохимического режима моря и его экосистемы [57].
Азовское море является уникальным физико-географическим объектом. Изучение гидрометеорологических и гидрохимических условий Азовского моря представляет большой практический и научный интерес. Существует ряд причин, обосновывающих особую актуальность создания инструмента исследования и прогнозирования состояния акватории моря
Безвозвратное изъятие значительной части стока рек, уменьшение стока рек Дона и Кубани из-за строительства гидросооружений привело к изменению биогенного и минерального стока, качественного состава вод, поступающих в море. Сезонное нивелирование стока рек привело к резкому снижению повторяемости, площади и продолжительности затопления пойменных нерестилищ. В последние годы произошло уменьшение ареалов размножения рыб из-за труднодоступности их привычных нерестилищ и ухудшения качества воды. Строительство судоходного канала в Таганрогском заливе и других гидротехнических сооружений изменяют основные гидрофизические показатели моря.
Несмотря на некоторое улучшение гидрологических и гидрохимических условий вследствие изменения экономической ситуации в стране в последние годы, проблема Азовского моря продолжает оставаться острой.
Использование экономичных разностных схем с треугольным оператором для решения задачи
В нашем случае- все условия следствия выполнены. Теорема доказана. Теорема 4.2. Пусть оператор А диссипативен, оператор В представим в виде(4.3.5). Если выполнено неравенство (4.3.7) и при этом Bs-0.5e M 0, М=А0-;ІК,-КиІІ=-11 (4.3.9) то схема (4.3.2)-(4.3.3) устойчива в HLt L=BS -0,5 эМ.
Доказательство. Подставим в операторное неравенство (4.3.8) значение оператора В0 из (.4.3.6). Тогда выполнение неравенства (4.3.9) эквивалентно удовлетворению условия (4.3.8) теоремы 1. Таким образом, теорема 2 доказана.
В схеме (4.3.2) - (4.3.3) с оператором В вида (4.3.5) есть некоторая свобода выбора оператора Bs и параметра &{т) в рамках полученных ограничений.
Следствие 1. Пусть оператор А диссипативен, оператор В представим в виде(4.3.5), причем В8=Е (3.3.10) Тогда разностная схема (2.3.22)-(2.3.23) устойчива по начальным данным в энергетическом пространстве HL, L=E Q.5U M, M=AQ-KL+KV. при выполнении условия т (о (4.3.11) Г где у - спектральный радиус матрицы и.
Доказательство. Если Bs = Е, то по теореме 2 для устойчивости схемы (4.3.2)-(4.3.3) в пространстве HL, 1=Е-0.5а М, достаточно выполнения условия Е-$.5а)М 0 . По свойствам собственных чисел эрмитовых матриц это требование эквивалентно 1—0.5е / 0, где у - спектральный радиус матрицы М. Отсюда получаем условия следствия 1. Следствие 1 доказано. Следствие 2. Если оператор А диссипативен, оператор В представим в виде(4.3.5), причем Bs=E+0.5aA (4.3.12) тогда если AtM, M=A0-j(Kv KL\ j=-1,1 то при выполнении условия (4.3,7) схема(4.3.2)-(4.3.3)устойчива в Н,, L=E+Q.5oy{A-M).
Доказательство. По условиям теоремы 1 схема (4.3.2)-(4.3.3) устойчива в пространстве HL при выполнении неравенств (4.3.7)- (4.3.8). Если оператор Bs имеет вид (4.3.12), то L=B0-0.5wA0=E+Q,56)(A-M) и неравенство (4.3.8) будет выполнено при условии Е+0.5в (А М) 0 или Л М. Следствие 2 доказано. Замечание. Пусть Л = } - диагональная матрица, представленная в одном из двух видов і)л,=НМІ N 2)ЛяХЫ» где т0 - элементы матрицы М, M=AV + KU KI. В случае А(.. = А=Л/, i=0JJ неравенство А М выполняется по [41], а во втором случае, N когда Ли = J] «j, , / = 0,JV , это неравенство выполняется по теореме 247 Гершгорина [26]. Таким образом, условие (2.3.28) выполняется для обоих вариантов построения оператора Л.
Приведем результаты численных экспериментов, позволяющие сравнить эффективность различных вариантов построения оператора В разностной схемы: 1. В=Е явная схема 2. В=Е+тК, треугольная кососимметричная схема 1(ТКС1) 3. B=E+2TKL треугольная кососимметричная схема 2(ТКС2) 4. B=E + 0.5TA,+TKL треугольная кососимметричная схема 3(ТКСЗ) 5. В=Е + 0.5тА2+тК, треугольная кососимметричная схема 4(ТКС4) 6. В=Е + тА] +2тК, треугольная кососимметричная схема5(ТКС5) 7. B=E + TA2+2TKL треугольная кососимметричная схема 6(ТКС6) Здесь использованы следующие обозначения: А ЯЕ, Л=\\МІ M=A0+Ka-KL,M={miJYj; Л2 = к/!и=і Л# = 5Ы = =i О, i j
Эксперименты проводились для четырех тестовых задач конвекции -диффузии, коэффициенты векторов скорости которых представлены в таблице 2.1. Это задачи, на которых, как правило, исследуются поведение численных методов решения - задача с постоянными коэффициентами, с линейными, с разделяющимися переменными и быстро меняющимися коэффициентами.
Начальные и краевые условия, а также вектор f(x x2,t) выбирались таким образом, чтобы точным решением задачи была функция p(xl X2,t) = е 8Ш(Ж,)8Ш(ЯГ2) ,
Задачи решались на временном отрезке fe[0,l] с шагом по времени, который выбирался таким образом, что его увеличение приводило к потере устойчивости схемы (2.3.22)-(2.3.23) или к росту относительной погрешности более чем на 20%. Решения получены на сетке 32x32 по пространству.
В таблицах 2-5 представлены результаты численных экспериментов решения задачи конвекции - диффузии при доминировании конвективных процессов (Ре: 10, ,103, 10s). Приведены значения максимального шага по времени, при котором схема устойчива и относительная погрешность (%) решения в норме пространства L2 на последнем шаге по времени отрезка [о,1 ] для каждого из четырех векторов скорости. В таблицах знак 8 означает любой шаг по времени для абсолютно устойчивых схем. Расчеты проведены при т=\.
