Использование эффективных численных методов при моделировании конвективно-диффузионных процессов в средах с преобладающей конвекцией

Любое крупное воздействие на окружающую среду связано с целой системой внешне отдаленных явлений и событий, приводит в действие процессы в различных сферах нашей жизни. Обязательным правилом при выяснении как ближайших, так и отдаленных последствий решений, стал системный подход. При таком подходе существует необходимость легкого и быстрого получения знаний о функционировании объекта исследования, прогноза его дальнейшего функционирования, предсказания и описания его воздействия на окружающую среду, экономической эффективности, количественных характеристик и рекомендаций, приводящих к заданным результатам. При теоретическом исследовании многих задач возникает необходимость исследования сложных процессов, решение которых без применения численных методов практически невозможно на современном этапе развития. В настоящее время эффективным инструментом исследования процессов и явлений, рассматриваемых при решении широкого круга научных и технических задач, служит математическое моделирование, предполагающее построение и изучение математических моделей путем решения на современных вычислительных комплексах соответствующих систем линейных алгебраических уравнений. Популярность математического моделирования обусловлена, прежде всего, значительной универсальностью математических моделей, возможностью применения вычислительного эксперимента в тех сферах исследований, где натурный эксперимент занимает много времени, дорог или невозможен. Этот метод сочетает в себе многие достоинства теории и практики. Работа с математической моделью дает возможность быстро и без существенных затрат промоделировать поведение объекта исследований в любых условиях и задавать любые требуемые свойства, а с другой стороны, вычислительные эксперименты позволяют более глубоко и подробно изучить объект во всей его полноте, чем при теоретическом исследовании. Более того, математическое моделирование позволяет осуществить с помощью одного устройства (ЭВМ) решение целого класса задач, имеющих одинаковое математическое описание, обеспечивает простоту перехода от одной задачи к другой, позволяет вводить переменные параметры, возмущения и различные начальные условия.

В математическом моделировании работу с объектом можно условно разделить на три этапа: модель — алгоритм — программа. На первом этапе строится математическая модель объекта, т. е. теоретический «эквивалент» объекта. В выбранной модели должны быть отражены основные свойства объекта, выбраны характеристики, описывающие состояние объекта, и те, которыми можно пренебречь, при этом учитываются условия, при которых будет функционировать объект. Выбор модели является одним из важнейших этапов моделирования. При выборе модели следует исходить из разумного компромисса между сложностью модели, полнотой получаемых с ее помощью характеристик объекта и точностью этих характеристик. Так, если модель недостаточно точна, то ее нужно дополнить, уточнить введением новых факторов, может также оказаться, что предложенная модель слишком сложна, и те же результаты можно получить с помощью более простой модели. Так, при моделировании процессов тепломассопереноса в жидкости обычно пренебрегают сжимаемостью среды и полагают жидкость несжимаемой, а при рассмотрении, в частности, процессов распространения подводных взрывов, необходимо рассматривать сжимаемую жидкость. При выборе математической модели обычно руководствуются иерархией математических моделей [19,49], при использовании которой можно выбрать модель, учитывающую свойства рассматриваемой среды и процессов, происходящих в ней, а также условия, при которых происходит ее исследование.

Второй этап — разработка алгоритма для реализации математической модели на ЭВМ. Этот алгоритм должен удовлетворять довольно жестким и противоречивым требованиям. С одной стороны, требуется получить «электронный эквивалент» модели за минимальное число действий, т. е. достаточно быстро и точно, а с другой стороны, объем обрабатываемой информации не должен превышать ресурсных возможностей современных вычислительных систем. При этом необходимо учитывать, что алгоритм не должен искажать основные свойства модели, и в то же время легко адаптироваться к особенностям задач, решаемых на компьютерах.

При решении задач обычно заменяют непрерывную область изменения аргумента на дискретную и строят дискретный аналог дифференциального уравнения на дискретной области. Этот подход позволяет свести решение дифференциального уравнения к системе линейных алгебраических уравнений. Вопросам разностной аппроксимации дифференциальных операторов уделяется большое внимание [1,8,41,42,43,44,46,104].

Полученную систему линейных алгебраических уравнений решают с помощью численных методов. Выбор численного метода зависит от вида полученной матрицы и ресурсных возможностей вычислительной техники.

Для решения системы линейных алгебраических уравнений с симметричной матрицей существует множество хорошо обоснованных и изученных эффективных численных методов. Например, можно отдать предпочтение методу верхней релаксации — (SOR) или его ускоренным вариантам, или методу сопряженных градиентов (CG) самостоятельно или с переобуславлива-телем. Если симметричная матрица системы не является положительно определенной, то можно применить метод минимальных невязок (MINRES).

Для несимметричных систем выбор метода решения до сих пор довольно затруднителен. Если матрица системы является слабо несимметричной (т.е. нормы симметричной и кососимметричной составляющих отличаются не сильно), или же матрицей с диагональным преобладанием, то можно воспользоваться классическими итерационными методами — методом Гаусса-Зейделя или другими методами. В противном случае мы имеем дело с несимметричными матрицами, для которых норма кососимметричной составляющей преобладает над нормой симметричной.

На третьем этапе выбранные алгоритмы реализуются и создаются программы. Полученный «электронный эквивалент» объекта уже можно испытывать, задавая различные условия функционирования исходного объекта.

Математическое моделирование широко используется при описании процессов в движущихся средах. Значительная сложность явлений вынуждает ученых не ограничиваться теоретическими исследованиями, но также использовать при изучении процессов методы математического моделирования. Математические модели в движущихся средах, которые включают в себя конвективный и диффузионный перенос, описывают самые различные газои гидродинамические процессы. В частности, большое значение приобретают экологические проблемы, связанные с описанием распространения загрязнений в атмосфере и водоемах, с моделированием загрязнения грунтовых вод.

Решение краевых задач такого типа обычно разбивается на два этапа: сначала для определения параметров движения среды отдельно решается уравнение газоили гидродинамики, а затем, используя полученные результаты, решается уравнение конвекции-диффузии, описывающее процессы переноса примеси.

При описании установившихся газои гидродинамических процессов стационарное уравнение конвекции-диффузии является эллиптическим уравнением второго порядка с младшими членами. Особые трудности данная задача приобретает, когда процесс конвективного переноса преобладает над диффузионным процессом. С математической точки зрения это означает, что в дифференциальном уравнении, описывающем конвективно-диффузионный процесс, появляется малый параметр при старшей производной. Наличие малого параметра при старшей производной при несогласованности правой части дифференциального уравнения с краевыми условиями может привести к появлению в решении тонких внутренних и пограничных слоев. Эти задачи очень тяжелы для численного анализа, так как необходимо найти сильно меняющееся в небольшой области решение.

При разностной аппроксимации дифференциального уравнения конвективно-диффузионного переноса возникает система линейных алгебраических уравнений с несимметричной матрицей. В случае преобладающей конвекции использование для разностной аппроксимации противопотоковых схем приводит к системе линейных алгебраических уравнений с монотонной матрицей (М — матрицей) и сильному сглаживанию решения за счет появления в разностных уравнениях искусственной вязкости. Поэтому при исследовании данных задач используется центрально-разностная аппроксимация, которая не сглаживает решение, но сводит дифференциальное уравнение к системе линейных алгебраических уравнений с несимметричной матрицей, не обладающей диагональным преобладанием.

Для этого класса матриц большинство классических итерационных методов либо вообще не работают, либо обладают очень медленной скоростью сходимости. Даже методы подпространства Крылова в этом случае также показывают невысокую скорость сходимости, поэтому необходимо использование переобуславливания системы. Следовательно, создание эффективных численных методов и переобуславливателей для решения СЛАУ с сильно несимметричной матрицей является актуальной задачей.

Цель и задачи работы. Целью данной работы является разработка эффективных численных методов решения систем линейных алгебраических уравнений, получающихся в результате центрально-разностной аппроксимации стационарного уравнения конвекции-диффузии в средах с преобладающей конвекцией.

Методы исследования рассмотренных итерационных методов и переобуславливателей основаны на операторном подходе, предложенном академиком А. А. Самарским и развитым его учениками, а также на теории итерационных методов, понятиях и методах матричного анализа.

Научная новизна. Научная новизна работы определяется полученными теоретическими результатами и предлагаемым новым классом попеременнотреугольных кососимметрических методов решения стационарных задач конвективно-диффузионного переноса с преобладающей конвекцией и сильно меняющимся полем скоростей, а также возможностью использования беспараметрических попеременно — треугольных операторов как переобуславли-вателей для методов подпространства Крылова.

Достоверность. Представленные в диссертации теоремы и следствия имеют строгое математическое обоснование, предложенные методы теоретически исследованы и численно проверенны.

Практическая значимость. С помощью разработанных методов можно эффективно решать стационарные задачи конвективно-диффузионного переноса с дискретным пространственным оператором, матрица которого дисси-пативна (т.е. симметричная часть матрицы положительно определена).

Апробация работы. Основные результаты диссертации докладывались на международной конференции «Математическое моделирование в экологии и численные методы» EMMNA'99 (г. Ростов-на-Дону, 1999 г.) — на третьей международной конференции по вычислительной математике ENUMATH 99 (Finland, 1999 г.) — на международном семинаре, посвященном вычислительным и аналитическим методам решения сингулярно возмущенных задач (Lozenetz, Bulgaria, 1998 г.) — на IX Всероссийской школе-семинаре молодых ученых «Современные проблемы математического моделирования» (п. Абрау-Дюрсо, 2001 г.) — на Всероссийской конференции «Математическое моделирование и проблемы экологической безопасности» (п. Абрау-Дюрсо, 2000 г.) — на VIII Всероссийском совещании по проблемам построения сеток для решения задач математической физики, посвященном памяти А. Ф. Сидорова (п. Абрау-Дюрсо, 2002 г.) — на Всероссийской молодежной научной школе — конференции «Численные методы решения линейных и нелинейных краевых задач» (г. Казань, 2001 г.) — на международной конференции «Математическое моделирование и вычислительный эксперимент в механике и физике» (г.Ростов-на-Дону, 2001 г.) — на международной конференции IMMC-2002 «Итерационные методы и матричные вычисления» (г. Ростов-на-Дону, 2002 г.) — на Международной конференции по вычислительной математике ICCM-2002 (г. Новосибирск, 2002 г.), на Международной конференции по вычислительной алгебре (MILOVY 2002) (Milovy, Czech Republic, 2002 г.).

В полном объеме диссертация докладывалась на научном семинаре «Методы решения краевых задач» лаборатории вычислительного эксперимента ЮГИНФО РГУ.

Публикации. По теме диссертации опубликовано 14 печатных работ, в том числе 9 в соавторстве. Из них 1 статья в российском и 1 статья в зарубежном реферируемых журналах, 5 статей в сборниках трудов и 7 в тезисах докладов всероссийских и международных конференций.

СОДЕРЖАНИЕ РАБОТЫ.

Первая глава посвящена математическому описанию конвективно-диффузионных процессов в движущихся средах. Приведены примеры математических моделей, сводящиеся к решению эллиптических уравнений второго порядка общего вида.

Это уравнение описывает, например, математическую модель распространения тепла и примесей в движущихся средах, моделирует распространение примесей в грунтовых водах. Во втором разделе приведены эквивалентные формы записи стационарного уравнения конвекции — диффузии, рассмотрены различные способы аппроксимации операторов конвективного и диффузионного переноса.

В третьем разделе первой главы дается описание задач, на которых были протестированы предложенные итерационные методы и переобуславливате-ли.

Во второй главе диссертации проведен анализ существующих классических и современных итерационных методов решения систем линейных алгебраических уравнений.

В первом разделе второй главы приведены некоторые определения и теоремы из теории матриц и функционального анализа.

Во втором разделе описываются два основных подхода к построению итерационных методов релаксационного типа — операторный, предложенный и развитый в трудах А. А. Самарского и его учеников, при котором исследуется норма оператора перехода, и спектральный, при котором исследуется спектральный радиус оператора перехода. В третьем параграфе второго раздела рассматривается подход к построению вариационных методов, при котором решение является линейной комбинацией векторов подпространства Крылова.

В третьем разделе дано описание классических итерационных методов и условия их сходимости, а также приведены их матрицы расщепления. Рассмотрены некоторые подходы к ускорению классических итерационных методов. Так, приведено краткое описание методов MSOR, AOR, USSOR, которые являются ускоряющими процедурами для SOR, методов ускорения релаксационного типа с помощью чебышевского набора итерационных параметров, способа ускорения двухслойных вариационных методов, предложенного А. А. Самарским.

В четвертом разделе второй главы дано описание двух подходов к построению вариационных методов. Первым был рассмотрен способ построения решения в подпространстве Крылова. При таком подходе искомое решение является линейной комбинацией некоторого базиса в подпространстве Крылова. Дано описание подходов Ритца — Галеркина, подхода минимума невязок, подхода Петрова — Галеркина. Используя каждый из этих подходов можно условно разделить методы подпространства Крылова на три большие группы:

• методы, в которых минимизируется невязка;

• методы, для которых невязка ортогональна текущему подпространству Крылова;

• методы, для которых невязка ортогональна некоторому другому подпространству Крылова.

Вторым был рассмотрен операторный подход, предложенный А. А. Самарским. При этом подходе строится решение уравнения, минимизирующее норму погрешности в энергетическом пространстве в случае самосопряженного и положительно определенного оператора. Далее в данном разделе более подробно описывается обобщенный метод минимальных невязок (GMRES).

Пятый раздел второй главы посвящен описанию ускорения методов подпространства Крылова с помощью переобуславливания (или иначе предобу-славливания). Основная идея этой методики заключается в том, что исходная система линейных алгебраических уравнений трансформируется в другую систему с матрицей, которая обладает лучшими свойствами, и итерационный метод сходится быстрее. Описаны правый, левый и двухсторонний способы переобуславливания. Приведено описание классического способа выбора явного переобуславливателя, при котором строится явное приближение к обратной матрице. Переобуславливатели такого вида называются явными пере-обуславливателями. Далее более подробно рассматривается переобуславли-вание обобщенного метода минимальных невязок, для которого будут предложены переобуславливатели специального вида.

Третья глава диссертации посвящена численным методам решения сильно несимметричных систем линейных алгебраических уравнений.

В первом разделе систематизированы современные методы решения систем линейных алгебраических уравнений с сильно несимметричной матрицей, приведены основные релаксационные методы и методы вариационного типа для решения данного класса задач.

Во втором разделе дано описание базового кососимметрического алгоритма решения сильно несимметричных задач, условно названного треугольным кососимметрическим методом. Данный алгоритм, предложенный JI.A. Крукиером и развитый в трудах его учеников, является эффективным методом решения сильно несимметричных систем линейных алгебраических уравнений. Для данного класса методов приводятся достаточные условия сходимости, описываются предложенные ранее способы ускорения треугольных кососимметрических методов.

В третьем разделе третьей главы рассматривается новый класс беспараметрических кососимметрических методов, для которых оператор метода не содержит итерационный параметр.

В первом параграфе раздела приведено полученное ранее необходимое и достаточное условие сходимости треугольных методов этого класса и достаточное условие сходимости метода. Во втором параграфе доказана теорема о выборе оптимального итерационного параметра.

Третий параграф посвящен численному исследованию предложенных ранее треугольных методов. В этом же параграфе предлагается ускоряющая процедура для треугольного кососимметрического метода.

Приводятся результаты проведенных тестов, на которых было показано, что использование специальной диагональной матрицы в треугольном операторе метода позволяет увеличить скорость сходимости треугольных кососимметрических методов.

В четвертом параграфе третьего раздела предлагается попеременно — треугольный беспараметрический метод. Для методов этого класса кососиммет-ричная составляющая оператора метода, как и для базовых беспараметрических треугольных методов, равна кососимметричной составляющей исходной матрицы.

Для данного метода получено достаточное условие сходимости. Было получено ограничение, накладываемое на диагональную матрицу оператора метода, которое позволяет создавать различные варианты попеременно-треугольных методов беспараметрических методов.

Для исследования предлагаются три выбора диагональной матрицы, и, как следствие, три способа построения попеременно-треугольных беспараметрических методов.

Приведены результаты расчетов предложенных кососимметрических беспараметрических попеременно-треугольных методов и проведено их сравнение с базовым попеременно-треугольным кососимметрическим методом и SSOR.

В пятом параграфе третьего раздела проведено сравнение результатов численных экспериментов для треугольных и попеременно-треугольных беспараметрических методов.

В четвертом разделе третьей главы описывается использование беспараметрических операторов в качестве переобуславливателей для GMRES (m).

Традиционно для ускорения методов подпространства Крылова используются явные переобуславливатели, т. е. явное построение матрицы, близкой к обратной матрице системы. Далее для ускорения обобщенного метода минимальных невязок рассматривается другой подход, называемый неявным, при котором ищется переобуславливатель, близкий к самой матрице системы. К переобуславливателям такого класса относятся, в частности SSOR. Для этого нужно оценить близость матрицы, обратной к переобуславливателю, к матрице, обратной к матрице системы.

В качестве меры близости двух матриц используется норма разности двух матриц. Доказана общая теорема, в которой оценивается близость матриц, обратных к переобуславливателю, к матрице, обратной к матрице системы.

Эта теорема была использована при изучении близости рассмотренных ранее треугольных и попеременно-треугольных беспараметрических операторов к матрице, обратной к матрице системы.

Для переобуславливания обобщенного метода минимальных невязок предложено использовать два попеременно-треугольных оператора.

Приводятся результаты численных экспериментов, проведенных на тестовых задачах для различных коэффициентов при конвективных членах и при различном количестве базисных векторов подпространства Крылова, используемых GMRES (m), т=2, 5, 10, 15, 20.

В выводах отражены результаты численного исследования предложенных методов на модельных задачах, проведено их сравнение, даны рекомендации о целесообразности и эффективности использования тех или иных итерационных методов в зависимости от особенностей решаемой задачи.

К ЗАЩИТЕ ПРЕДСТАВЛЕНЫ СЛЕДУЮЩИЕ РЕЗУЛЬТАТЫ:

1. Предложен и теоретически обоснован попеременно-треугольный беспараметрический кососимметрический метод решения систем линейных алгебраических уравнений, получаемых при центрально-разностной аппроксимации стационарного уравнения конвекции-диффузии.

2. Предложены и теоретически обоснованы ускоряющие процедуры для треугольных и попеременно-треугольных кососимметрических итерационных методов.

3. Доказана и численно подтверждена эффективность использования указанных операторов в качестве переобуславливателей для обобщенного метода минимальных невязок.

Автор глубоко благодарен своему научному руководителю проректору РГУ, директору ЮГИНФО РГУ, д.ф.-м.н., проф. Крукиеру Л. А. и благодарит коллектив ЛВЭ ЮГИНФО РГУ за внимание к работе, оказанную помощь и полезные советы.

1 МОДЕЛИРОВАНИЕ ПРОЦЕССОВ В ДВИЖУЩИХСЯ СРЕДАХ.

Выводы:

Предложенный класс беспараметрических итерационных методов показал свою эффективность при решении стационарного уравнения конвекциидиффузии с преобладающей конвекцией и сильно меняющимся полем скоростей. Это можно объяснить тем, что беспараметрические треугольные и попеременно-треугольные операторы более явно учитывают поведение кососимметричной составляющей, чем SOR и SSOR. Можно отметить схожесть процедур ускорения предложенных методов с MSOR и USSOR, но в отличие от них, в операторах рассмотренных методов учитываются коэффициенты матрицы не в целом блоке, а в каждой строке.

Наилучшие результаты были показаны в случае сильно меняющегося поля скоростей и при больших числах Пекле. Это улучшение было особенно заметно при использовании предложенных операторов в качестве переобуславливателей для GMRES (IO), т=10,15,20. При небольшом числе базисных векторов иногда наблюдалась нестабильность работы предложенных переобуславливателей, особенно в случае треугольного беспараметрического оператора.

Главным недостатком методов релаксационного типа, который часто отталкивает ученых от использования данного класса методов, является наличие итерационного параметра в операторе метода, а также очень сильная зависимость скорости сходимости методов от выбора итерационного параметра. Более того, несмотря на то, что существуют аналитические способы выбора оптимального итерационного параметра в методах релаксационного типа, на практике приходится выбирать итерационный параметр эмпирическим путем, так как часто невозможно вычислить точные границы спектра матрицы. Более того, для сильно несимметричных матриц многие методы данного класса либо сходятся очень медленно, даже при оптимальном выборе итерационного параметра, либо вообще перестают сходиться.

При проведении вычислительных экспериментов обнаружилось, что предложенные треугольные и попеременно-треугольные беспараметрические методы сходятся для СЛАУ с сильно несимметричными матрицами, возникающими при центрально-разностной аппроксимации стационарного уравнения конвекции-диф фузии с преобладающей конвекцией и сильно меняющимся полем скоростей, а также не сильно чувствительны к выбору итерационного параметра.

Вместе с тем, хотя методы подпространства Крылова хорошо зарекомендовали себя для решения СЛАУ с несимметричной матрицей, для сильно несимметричных систем наблюдается очень медленная сходимость. Как показано в литературе по численным методам [61,132], использование ILU пере-обуславливания для сильно несимметричных задач представляется затруднительным [69,119], метод SOR не рекомендуется использовать в качестве пе-реобуславливателя [61], a SSOR очень чувствителен к выбору итерационного параметра [77,110] и не показывает существенно увеличение скорости сходимости для задач этого класса.

Разработанные операторы можно порекомендовать к использованию в качестве переобуславливателей как к методу простой итерации, так и для вариационных методов.

Можно считать, что предложенные беспараметрические операторы являются достаточно хорошими переобуславливателями метода обобщенных минимальных невязок с числом базисных векторов, большим 10, для решения СЛАУ с несимметричными матрицами, получающимися при центрально-разностной аппроксимации стационарного уравнения конвекции-диффузии с преобладающей конвекцией и сильно меняющимся полем скоростей,.