автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Многосеточная технология для математического моделирования тепловых и гидродинамических процессов

005061317

На правах рукописи

Мартыненко Сергей Иванович

МНОГОСЕТОЧНАЯ ТЕХНОЛОГИЯ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ТЕПЛОВЫХ И ГИДРОДИНАМИЧЕСКИХ ПРОЦЕССОВ

Специальность 05.13.18 Математическое моделирование, численные методы и комплексы программ

о ш 2013

АВТОРЕФЕРАТ диссертации на соискание учёной степени доктора физико-математических наук

Москва - 2013

005061317

Работа выполнена в федеральном государственном унитарном предприятии «Центральный институт авиационного моторостроения имени П.И. Баранова»

Официальные оппоненты: Бадриев Ильдар Бурханович, чл.-корр. АН

Республики Татарстан, доктор физико-математических наук, профессор кафедры вычислительной математики Института вычислительной математики и информационных технологий федерального государственного автономного образовательного учреждения высшего профессионального образования «Казанский (Приволжский) федеральный университет»

Жуков Виктор Тимофеевич, доктор физико-математических наук, старший научный сотрудник, заведующий отделом прикладных задач механики сплошных сред федерального государственного бюджетного учреждения науки «Институт прикладной математики имени М.В. Келдыша РАН»

Станкевич Игорь Васильевич, доктор технических наук, старший научный сотрудник, профессор кафедры прикладной математики федерального государственного образовательного учреждения высшего профессионального образования «Московский государственный технический университет имени Н.Э. Баумана»

Ведущая организация: Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт» (государственный университет).

Защита состоится 2 июля 2013 года в 11 ч. 00 мин. на заседании диссертационного совета Д 212.141.15 при ФГБОУ ВПО «Московский государственный технический университет имени Н.Э. Баумана» по адресу: Москва, Рубцовская наб., 2/18, ауд. 1006 л.

С диссертацией можно ознакомиться в библиотеке ФГБОУ ВПО «Московский государственный технический университет имени Н.Э. Баумана».

Автореферат разослан ^ -¿¿-а^У 2013 г.

Учёный секретарь диссертационного совета, кандидат технических наук, доцент

A.B. Аттетков

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность проблемы. В настоящее время методы математического моделирования широко используются для решения крупных научно-технических и хозяйственных задач. Однако численному моделированию, как и любому методу исследования, присущи определённые недостатки. В первую очередь к ним следует отнести трудоемкость написания и отладки компьютерных программ, предназначенных для решения сложных прикладных задач. Поэтому широкое распространение получил ряд коммерческих программных продуктов (ANSYS, STAR-CD, FLUENT, CFX, PHOENICS и др.), устроенных по принципу «черного ящика».

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

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

Состояние вопроса исследования. В настоящее время многосеточные методы, предложенные выдающимися отечественными математиками Р.П. Фе-доренко и Н.С. Бахваловым, стали доминирующими алгоритмами для решения многих прикладных задач. Бурное развитие многосеточных методов, начавшееся в середине 80-х годов, обусловлено их уникальной возможностью получать численное решение с оптимальными вычислительными усилиями. Развитие многосеточных методов пошло по традиционному для середины 80-х годов пути, связанным с адаптацией алгоритма к решаемой задаче. В результате многосеточные методы быстро превратились в труднообозримое семейство алгоритмов, причем применение их к решению новых задач приводило к появлению новых вариантов.

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

Другим принципиальным недостатком классических многосеточных ме-

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

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

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

Задачами исследования являются:

1. Разработка многосеточной технологии с минимальным количеством проблемно-зависимых компонент:

- разработка конструкции многосеточной технологии с высоким уровнем формализации;

- доказательство сходимости и анализ объема вычислительной работы; оценка эффективности распараллеливания вычислений;

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

2. Разработка программного обеспечения и алгоритмизация отдельных компонент технологии.

3. Тестирование технологии на различных модельных задачах.

4. Разработка и тестирование способа ускорения сходимости итерационных методов решения уравнений Навье-Стокса, основанного на декомпозиции

давления и физических аспектах гидродинамики.

5. Разработка и тестирование высокоэффективного многосеточного метода численного решения уравнений Навье-Стокса.

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

Научная новизна:

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

Разработано программное обеспечение и выполнено тестирование метода на различных (не)линейных задачах.

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

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

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

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

Апробация работы. Результаты диссертационной работы докладывались на XIII Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов для решения задач математической физики» (Пущино, 2000), Международной конференции «Optimization of Finite Element Approximations & Splines and Wavelets» (Санкт-Петербург, 2001), V International Congress on Mathematical Modelling (Dubna, 2002, 2009), Российской национальной конференции по теплообмену (Москва, 1994, 1998, 2006, 2010), Всероссийском семинаре «Сеточные методы для краевых задач и приложения»

(Казань, 2004, 2005, 2007, 2010), International Conference «Numerical geometry, grid generation and high performance computing» (Moscow, 2006, 2008), Всероссийской конференции «Необратимые процессы в природе и технике», (Москва, 2005, 2007, 2009, 2011), Международной конференции по неравновесным процессам в соплах и струях (Санкт-Петербург, 2002; Алушта, 2008), Международной конференции по вычислительной механике и современным прикладным программным системам (Истра, 2001; Владимир, 2003), Международной конференции «Авиация и космонавтика» (Москва, 2005, 2008), Международной научно-техническая конференция «Авиадвигатели XXI века» (Москва, 2005, 2010), XXI Всероссийском семинаре «Струйные, отрывные и нестационарные течения» (Новосибирск, 2007), Международной конференции «Шестые Оку-невские чтения» (Санкт-Петербург, 2008), III конгрессе «Авиация в XXI столетии» (Киев, 2008), Восьмой Международной школе-семинаре «Модели и методы аэродинамики», (Евпатория, 2008), Всероссийской конференции «Механика и наномеханика структурно-сложных и гетерогенных сред» (Москва, 2009), Fourth International Conference «Computational Methods in Applied Mathematics» (International Banach Center at Bedlgwo/Poznan, Poland, 2010), Всероссийской конференции по вычислительной математике (Новосибирск, 2011), а также на семинарах в ФГУП «Центральный институт авиационного моторостроения имени П.И. Баранова» (2006, 2011), Институте прикладной математики имени MB. Келдыша РАН (1995), Харбинском политехническом университете (Харбин, 2007, 2009) и Северо-Западном политехническом институте (Сиань, 2009).

Диссертация является составной частью фундаментальных научных исследований, проводимых при поддержке Российского фонда фундаментальных исследований (проекты №09-01-00151 и №12-01-00109).

Публикации. Основные научные результаты диссертационной работы отражены в 42 научных работах, в том числе в 15 статьях Перечня ведущих рецензируемых научных журналов и изданий ВАК, и в 16 тезисах докладов.

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

Структура и объем работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы. Работа изложена на 268 страницах, содержит 92 иллюстраций и 29 таблиц. Библиография включает 167 наименований.

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

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

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

Первая глава диссертации, состоящая из одиннадцати разделов, посвящена описанию, теоретическому обоснованию и тестированию универсальной многосеточной технологии (УМТ). В данной работе универсальность алгоритма определяется количеством проблемно-зависимых компонент. Термин «технология» означает совокупность приемов обработки, указывая на наличие ана,-литической и вычислительных частей УМТ.

В разделе 1.1 сформулирована проблема построения универсального и эффективного алгоритма при следующих ограничениях:

- минимальное использование ресурсов компьютера;

- возможность изменения способа и/или порядка аппроксимации (в процессе решения задачи);

- минимальное количество проблемно-зависимых компонент;

- возможность эффективного распараллеливания вычислений.

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

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

Самая мелкая сетка С? состоит из двух множеств точек (У(0;1) = {х*} и (0;1) = {а;'}, которые в зависимости от решаемой задачи могут быть узлами или гранями контрольных объемов.

Построение грубых сеток в УМТ осуществляется посредством удаления точек из множеств (7^0;1) и Сг(0;1) самой мелкой сетки С", как показано на рис. 1. Таким образом из самой мелкой сетки С" могут быть получены три грубые сетки С\, С\ и (З^, которые обладают следующими свойствами:

1) грубые сетки С}, и не имеют общих точек, т.е.

2) мелкая сетка С® представима в виде объединения грубых сеток (?}, С\ и Ст\

к= 1

Х2 3*3 3*4 «Е^

С? С}

с? с;

^Ст Хо Хъ «Сд

5

Хс X,

Х7 Ха Хо

3<7 Ха Хо Х|

Самая мелкая сетка

1-1-1

1-1-1

1-1-1

-НН-I

| Первая сетка первого уровня

Самая мелкая сетка

I - I т 1 ■ I - I т I - 1 ■ 1 т I

Вторая сетка первого уровня

С" —+—+

Самая мелкая сетка

1-1-1

1-1-1

1-1-1

Третья сетка первого уровня

Рис. 1. Построение грубых сеток в универсальной многосеточной технологии

3) все сетки геометрически подобны, однако шаг грубых сеток <3|, <32 иб^ в три раза больше, чем шаг сетки (ЗУ;

4) вне зависимости от способа определения сеточных функций на самой мелкой сетке каждый контрольный объем на сетках б}, и является объединением трех контрольных объемов на сетке С?®.

Самая мелкая сетка (7° образует нулевой сеточный уровень, а три грубые сетки С\, С?2 и Сз образуют первый сеточный уровень. Далее построение еще более грубых сеток осуществляется рекуррентным образом: каждая сетка С', i = 1, ...,3' уровня I рассматривается как самая мелкая сетка для трех грубых сеток j = 1,..., 3,+1, следующего уровня 1 + 1. Девять еще более грубых сеток, полученных из трех сеток первого уровня, образуют второй сеточный уровень и так далее, как показано на рис. 2. Построение грубых сеток завершается, когда на грубых сетках останется всего несколько точек х* и х1. В дальнейшем совокупность самой мелкой и всех грубых сеток будет называться мпогосеточной структурой.

Номер сеточного уровня с самыми грубыми сетками, обозначаемый как [Л, вычисляется перед построением грубых сеток по формуле

Г^ГЗГ

ь+ =

18 3

- 1

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

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

Самая мелкая сетка G? Нулевой уровень

Рис. 2. Многосеточная структура

Рис. 3. Пилообразный цикл и многосеточная итерация

состоит из 3Nt сеток, где N = 2, 3.

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

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

В разделе 1.4 основные компоненты УМТ проиллюстрированы на примере простейшего двухуровневого алгоритма. Показано, что разностные уравнения на каждом сеточном уровне могут быть записаны в виде СЛАУ

Alct = Rf^lr, / = 0,1,2,..., L+ ,

где А1 есть блочно-диагональная матрица, причём число блоков равно числу сеток данного уровня 1\ С1 есть вектор поправок; ~ оператор, проецирующий невязку г = (Ь — Ах^)0 с самой мелкой сетки (нулевой уровень) на сетки уровня I.

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

С;_х = Р,^_1Сг, 1 = Ь+,Ь+- 1,...,1,

где Р/_г_ 1 есть оператор пролонгации. Согласно рис. 4 оператору пролонгации соответствует перестановочная матрица, т.е. оператор также не зависит

от решаемой задачи.

В разделе 1.5 приведено доказательство сходимости УМТ. Сглаживающая итерация на сетках уровня I может быть записана в виде

ИЧс,('"+1) - с,М) = - - ,

где матрица IV зависит от выбора сглаживателя.

Доказательство сходимости УМТ проведено при следующих допущениях: Допущение 1. На каждом уровне I выполняется щ сглаживающих итераций, матрица итераций удовлетворяет ограничению

т.е. данный итерационный метод сходится.

Допущение 2. На уровне с самыми грубыми сетками (Ь+) разностные уравнения решаются точно, т.е.

О [ 1 1 1 1 1 1 1 1 1 1 , | 1 1 :

1111. "1! 1 Г 1 ! |

¡и. ' 1 1 " 1

Граница области Граница области

Рис. 4. Переход от более грубых сеток к мелкой сетке

Тогда многосеточные итерации УМТ можно записать в виде

при I = 0,1,...,1+ -2,

Для оценки нормы матрицы М обычно используют следующее допущение: Допущение 3 (свойство аппроксимации)

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

Результаты вычислительного эксперимента по решению первой краевой задачи для уравнения Пуассона на шестиуровневой многосеточной структуре с самой мелкой сеткой 1001 X 1001 использованы для иллюстрации результатов теоретического анализа.

Подсчёт количества арифметических операций показал, что УМТ не обладает оптимальной скоростью сходимости, т.е. для получения численного решения необходимо выполнить ^ М) арифметических операций, где N есть число узлов самой мелкой сетки. Множитель N является следствием независимости оператора пролонгации от решаемой задачи.

В разделе 1.6 рассмотрены различные варианты распараллеливания УМТ. Поскольку грубые сетки каждого уровня не имеют общих точек, то сглаживающие итерации на этих сетках могут проводиться параллельно вне зависимости от используемой сглаживающей процедуры. Распараллеливание ЛГ-мерных задач целесообразно проводить на многопроцессорном компьютере с р = 3Л',: процессорами, где параметр к ^ Ь+ называется глубиной распараллеливания.

Распараллеливание сглаживающих итераций на уровнях с грубыми сетками при к = 1 показано на рис. 5. Поскольку на данных уровнях УМТ обладает полным параллелизмом, то в ряде случаях возможно проведение дополнительных многосеточных итераций д* (д* = 3 на рис. 5).

-1

Самая мелкая сетка

Самые грубые сетки

Первый модуль

Самые грубые сетки

Второй модуль

Самые грубые сетки

3Л-ый модуль

Рис. 5. Распараллеливание сглаживающих итераций на уровнях с грубыми сетками

В предположении, что на всех уровнях выполняется одинаковое число сглаживающих итераций, получена следующая оценка ускорения (¿У и эффективности параллелизма

= рЁх = р-4 1 х+ 1-,

4- — + ар

ьо

где р = есть число используемых процессоров, д" - число дополнительных многосеточных итераций на уровнях с грубыми сетками (I ^ 1), Ь+ - номер уровня с самыми грубыми сетками, Е0 - эффективность распараллеливания сглаживающих итераций на самой мелкой сетке (I = 0), а - коэффициент, учитывающий потери на обмены данными.

Учёт потерь на обмены данными приводит к снижению ускорения и эффективности параллелизма:

ар

Ег

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

Эффективность распараллеливания УМТ, если на всех сеточных уровнях выполняется одинаковое число одних и тех же сглаживающих итераций, лежит

в пределах:

а) распараллеливание первой глубины (9 процессоров для двухмерных задач (ЛГ = 2) или 27 процессоров для трёхмерных задач (./V = 3))

Ь+ + 1 _ д*Ь+ + 1 <Е<

о

б) распараллеливание второй глубины (81 процессор для двухмерных задач (/V = 2) или 729 процессоров для трёхмерных задач (М = 3))

ь+ + 1 <Е< я\Ь+ - 1) + 2

В разделе 1.7 УМТ протестирована на примере решения двухмерных краевых задач для уравнения Пуассона (второй и четвёртый порядок аппроксимации), уравнения с квадратичной нелинейностью вида

д2и д2и - .

_ + = -/(*, у)>

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

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

В разделе 1.9 представлены результаты математического моделирования распространения тепла в корпусе микродвигателя при нестационарном режиме работы. В микродвигателях невозможно обеспечить полное смешение горючего и окислителя, поэтому микродвигатели работают на однокомпонентном топливе. Топливо поступает на катализатор, разлагается с выделением теплоты, и продукты разложения истекают в космическое пространство, создавая тягу. Конструкция микродвигателя (со снятой крышкой) показана на рис. 6. Элементы конструкции (катализатор и сопло Лаваля) укрупненно ноказаны на рис. 7. Двигатель может работать как в непрерывном, так и в импульсном режимах.

Рис. 6. Микродвигательная установка

Одной из ключевых проблем разработки микродвигателей является расчёт прогрева корпуса во время работы. Математическая модель состоит из трехмерного уравнения теплопроводности и уравнений, описывающих зависимость теплофизических свойств материала корпуса и гидразина от температуры. Для расчетов использована сетка 401 х 301 х 101. Температура наружной поверхности корпуса принималась равной 20°С, а катализатора - 800°С. Распределение температуры в плоскости контакта корпуса и крышки показано на рис. 8 и 9.

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

Рис. 7. Элементы конструкции микродвигательной установки

задачи многосеточная итерация представима в виде следующей последовательности действий:

1) аппроксимация правой части Е-модифицированной задачи на НСС;

2) вычисление невязки на НСС;

3) интерполяция невязки с НСС на ВСС;

4) аппроксимация левой части S-модифицированной краевой задачи при помощи интегро-интерполяционного метода на ВСС;

5) решение СЛАУ на ВСС при помощи УМТ;

6) интерполяция поправки с ВСС на исходную НСС;

7) выполнение s сглаживающих итераций на НСС;

8) пересчет приближения к решению, обнуление поправки;

9) проверка критерия останова,, возврат к п. 1 (если необходимо).

Данная итерация может быть записана в виде

а>+1> = (/ - W-1Aa)s+1 (i - П ^п1 П

\ D—Д Д^П /

+ ( (7 - w-'AAy+1 п л-1 п + £ (/ - W-UA)m ■ wA К, L □—д д—□ т=0 J

где матрица А получена в результате аппроксимации исходного дифференци-

X

Рис. 8. Температура корпуса микродвигателя через 0.020 сек после включения

X

Рис. 0. Температура корпуса микродвигателя через 0.058 сек после включения

ального оператора; матрица W определяется выбором сглаживающей процедуры на HCC; f| - оператор-проектор, индексы д и □ означают принадлежность к НСС и ВСС соответственно. В случае выполнения свойств сглаживания и аппроксимации

□_Д Д—

количество итераций не будет зависеть от количества узлов НСС. Интерполяция невязки с НСС на ВСС и поправки с ВСС на НСС приводит к снижеиию универсальности алгоритма. Разработана методика построения ВСС.

В разделе 1.11 предложена методика, адаптивного определения параметров релаксации, основанная на предположении о постоянстве их оптимального значения на сетках одного уровня. Вычислительные алгоритмы для решения прикладных задач часто содержат различные параметры релаксации, определение оптимальных значений которых представляет довольно трудную зада,-чу. Оптимальное значение определялось путём анализа скорости сходимости итерационного метода на одной сетке каждого уровня. На самой мелкой сетке использовалось значение, полученное путем экстраполяции с предыдущих уровней. Выполнен вычислительный эксперимент с решением первой краевой задачи для уравнения Пуассона при помощи метода последовательной верхней релаксации. Показана работоспособность методики и дан анализ объема дополнительной вычислительной работы. Адаптивное определение параметра релаксации приводит к увеличению времени счёта на ~ 10% для двухмерных задач и на ~ 3% для трёхмерных задач по сравнению со случаем, когда заранее известно его оптимальное значение.

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

Вторая глава диссертации, состоящая из трёх разделов, посвящена описанию разработанного программного обеспечения для УМТ, исполненного на алгоритмическом языке FORTRAN9Û. Дано описание подпрограмм, реализующих основные компоненты УМТ, в одно- и многопроцессорном исполнении.

В разделе 2.1 дано описание модуля RMT_2D.f90, предназначенного для решения двухмерных задач. Представлена компактная схема хранения мно-

^ СлИЛдГ1

МНОГОСЕТОЧНЫЕ МЕТОДЫ

Представление решения в виде поправки и приближения

Да

Нет

До аппроксимации После аппроксимации

УМТ

КММ

Каскадные методы

Рис. 10. Классификация многосеточных методов

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

Аппроксимация модифицированных краевых задач интегро-интерполя-ционным методом приводит к необходимости вычисления шаблонных функционалов вида

где (<7){;} есть среднее значение функции д(х) на отрезке (контрольном объеме) •7-{;}]> ~ шаг самой мелкой сетки, {¿} - отображение индексов узлов текущей сетки, на индексы самой мелкой Сетки. Доопределим функцию д(х), полагая ее равной нулю вне области П = [0,1], т.е.

0, хфП q(x), х eil

Далее определим характеристическую функцию Т(х) как

Теперь коэффициент можно записать в виде

-0.5 0.0 0.5 1.0 1.5

Рис. И. Схема вычисления определенных интегралов

где коэффициенты (<?){;} и Л^ имеют вид

= \ ! йх и л'{<} = | I Г(х) йх.

1Г

т'

По своему смыслу Л^ есть число контрольных объемов на самой мелкой сетке, которые образуют реальную часть данного контрольного объема.

Поскольку каждый контрольный объем на грубых сетках уровня I является объединением трех контрольных объемов на более мелких сетках уровня / — 1, то коэффициенты и на всех сетках текущего уровня I* вычис-

ляются при помощи следующих рекуррентных соотношений

' Щ1} = + <<7>м + (^Тнз- / = 1

Л( _ Д1-1 , Д/-1 , Д/-1 1 ' • • ■ ' •

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

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

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

В разделе 2.2 дано описание модуля РМТ_30.£90, предназначенного для решения трёхмерных задач.

В разделе 2.3 систематизированы основные особенности программной реализации УМТ.

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

В разделе 3.1 проанализированы вычислительные алгоритмы для решения упрощённых уравнений Навье-Стокса (в приближении пограничного слоя или длинного канала). На основании анализа предложено представление давления виде

p{t, X, у, z) = px{t, х) + Pv(t, у) + p'{t, z) + pxyz{t, x, y, z), (1)

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

Основная идея метода ускорения сходимости состоит в вычислении «части» давления (a. именно, px(t, х) +py(t, у) +pz(t, z)) при помощи высокоэффективных итерационных методов, предложенных ранее для решения упрощенных уравнений Навье-Стокса.

Отличительными особенностями декомпозиции давления являются:

1) Каждое из слагаемых в (1) по отдельности лишено физического смысла; физический смысл имеет лишь их сумма. Далее слагаемые вида px{t,x), pv(t,y) или pz(t, z) будут называться «одномерными», а слагаемое pxyz(t, х, у, z) «многомерным»; кавычки «» указывают на «искусственность» каждого слагаемого.

2) Согласно (1) давление представимо в виде суммы N «одномерных» слагаемых и одного «многомерного» (N = 2,3). Поэтому необходимо сформулировать N дополнительных условий для определения «одномерных» слагаемых. В качестве таких условий используются уравнения постоянств массового расхода, которые рассматриваются как априорная информация физического характера об искомом решении. Для моделирования течения в каверне с движущейся крышкой уравнения постоянств массового расхода получены путём интегрирования уравнения неразрывности по контрольным объёмам Vi и V2, показанным на рис. 12, т.е.

1 1

J u(t, х, у) dy = 0, У v(t, X, у) dx = 0, о о

где и и v есть компоненты вектора скорости.

3) Поскольку уравнения постоянств массового расхода следуют из уравнения неразрывности, то «одномерные» слагаемые должны определяться перед «многомерным».

4) Несмотря на представление давления в виде суммы N + 1 слагаемых, в уравнения движения входят градиенты только двух из них, а именно «многомерного» и соответствующего «одномерного». Например, градиент давления в

Рис. 12. Каверна с движущейся крышкой и контрольные объемы и Уч

уравнении движения по X в этом случае принимает вид

др дх

дрх дрхуг дх дх

5) Уменьшение объема вычислений (по сравнению с традиционными алгоритмами х) = у) = рг(^ г) = 0) зависит от типа течения. Максимальный выигрыш в эффективности возникает при доминировании градиентов «одномерных» слагаемых, т.е. при моделировании течений с выделенным направлением течения среды (течения в струях, соплах, каналах, обтекание тел и т.д.). Минимальный выигрыш в эффективности наблюдается при моделировании циркуляционных течений (например, течение в каверне с движущейся крышкой).

С) В трехмерном случае предлагаемый подход более эффективен, чем в двухмерном, поскольку большая часть давления отыскивается при помощи упрощенных уравнений Навье-Стокса.

7) Компоненты скорости и соответствующие им «одномерные» слагаемые всегда вычисляются совместно, «многомерное» слагаемое может быть вычислено как совместно, так и сегрегированно.

В разделе 3.2 рассмотрены пути совершенствования явных схем (схема расщепления по физическим факторам, метод искусственной сжимаемости и др.). Показано, что компоненты скорости и соответствующие «одномерные» слагаемые в представлении (1) могут быть определены совместно, сохраняя явный характер вычислений.

Вычислительный эксперимент по моделированию нестационарного течения в каверне с движущейся крышкой (Лх = Ну = к = 1/200, = /г/5, Яе = 1000) использован для иллюстрации минимального ускорения явной схемы расщепления по физическим факторам. Закон движения крышки задан в виде

и{п)

(-И-

Л100

Результаты вычислительного эксперимента показаны на рис. 13, где Т1/"' и Тс

время счёта по модифицированной и классической схемам на п-ом временном слое соответственно. Уменьшение времени счёта составило примерно 20%.

В разделе 3.3 рассмотрены пути совершенствования неявных разностных схем. В этом случае «одномерные» слагаемые в представлении (1) не могут быть определены в явном виде, поэтому для их отыскания необходимо сформулировать вспомогательную задачу.

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

а) уравнение движения по X и уравнение постоянства массового расхода

ди д(и2) д(уи) дЬ дх ду

дрх дх

дрху

дх

+ Ре +

д2и\ ду*)'

$ и(Ь,х,у) йу ■ о

0;

б) уравнение движения по У и уравнение постоянства массового расхода

ду

т +

д{ии) д(у2)

дх

ду

дру

ду

дрх

ду

+

Ре + ду2)'

J у(Ь, х, у) йх = 0 ,

где квадратные скобки означают, что производные (рху)'х и (р1у)'у являются фиксированными, а фигурные скобки означают, что уравнения движения решаются согласованно с соответствующими уравнениями постоянств массового расхода.

1.0

0.9

с:

\ 0.8 -0.7 -

0.6

Г-20

—Г-40

Г" 80

I

100

-1-Г-

120 140

I

160

—I— 180

—I 200

номер временного слоя п Рис. 13. Уменьшение времени счёта при моделировании течения в каверне

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

Далее значения «одномерных» слагаемых фиксируются и осуществляется переход к основной задаче, в которой уравнения движения принимают вид

где квадратные скобки означают, что производные (рт)'х и (ру)'у уже вычислены во вспомогательной задаче.

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

Микрокатализатор состоит из игл, расположенных в шахматном порядке и покрытых каталитически активным веществом (рис. 7). Геометрия одного из вариантов микрокатализатора показана на рис. 14. Результаты моделирования (Ре = 350, сетка 385 х 3150) показаны на рис 15 и 16. Для расчёта использован счёт па установление с применением непредобусловленного метода Уза-вы для вычисления давления. Показано, что конструкция микрокатализатора обеспечивает безвихревое течение среды за исключением последнего ряда игл. Интенсивное вихреообразование за последним рядом игл приводит при более высоких значениях числа Рейнольдса к появлению дорожки Кармана. В этом случае возможна неустойчивая работа двигателя. На основании результатов моделирования выданы рекомендации по изменению геометрии игл последнего ряда для уменьшения вихреообразования.

Экспериментально показана высокая эффективность метода при моделировании течений с выделенным направлением движения среды.

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

Рис. 14. Геометрия микрока.тнлиза.тора,

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

В разделе 3.5 систематизированы основные особенности принципа формальной декомпозиции давления, связанные с возможностью ускорения сходимости итерационных методов решения уравнений Навье-Стокса без использования проблемно-зависимых компонент.

Четвёртая глава диссертации, состоящая из семи разделов, посвящена построению многосеточного алгоритма для решения уравнений Навье Стокса на основе УМТ и декомпозиции давления.

В разделе 4.1 получен вид Е-модифицированных уравнений Навье-Сток-са для последующего их решения при помощи УМТ. Исходные уравнения в переменных «скорость-давление» записываются в следующем виде

М(У) + УР = К УК = <3,

где Л/* есть конвективно-диффузионный оператор; V, Р - вектор скорости и давление; .Р, С? - источниковые члены. Решение уравнений Навье -Стокса пред-

Рис. 15. Изолинии безразмерной функции тока вблизи первого ряда игл микрокатализатора около оси симметрии

Рис. 16. Изолинии безразмерной функции тока вблизи последнего ряда игл микрокатализатора около оси симметрии

ставляется в виде

V = Cv + V, Р = СР + Р,

где разностные аналоги функции C'y и СР будут играть роль поправок, а аналоги функций V и Р - приближений к решению в последующих многосеточных итерациях. С учётом

M{V)=M{Cy + V)=N{V)+N\Cv), Я* ¿Я,

Е-модифицированные уравнения Навье-Стокса принимают вид N*{Cy) + ЧСР = F*, VCv = G*,

где

F' = F - Af(V) - VP, G* = G - VV.

Нелинейность конвективно-диффузионного оператора Ai приводит к появлению дополнительных членов в Е-модифицирова,нных уравнениях Навье-Стокса, т.е. А/4 ф N'. Отличие от исходных уравнений состоит в виде оператора N* и источниковых членов F*, G*. Е-модификация уравнений Навье-Стокса не требует их линаризации, а численное решение при помощи УМТ - применения схемы FAS.

В разделе 4.2 рассмотрена конфигурация контрольных объёмов, используемых для аппроксимации модифицированных уравнений Навье-Стокса на разнесённой сетке.

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

В разделе 4.4 показаны особенности Е-модификации вспомогательной задачи.

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

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

Г Ла= -BTßM + f

где матрица Q есть некоторый предобуславливатель, а, ß, f, g - сеточные аналоги Cv, СР, F* и G* соответственно. Численное решение уравнений Навье-Стокса получено на 10000 временных слоях при использовании диагонального

Рис. 17. Геометрия задачи о течении в квадратной полости

предобуславливания, т.е. С} = ВБ^1 Вт, где матрица БА образована диагональными элементами матрицы А. Подобной подход применим для решения нестационарных задач в силу сильного диагонального преобладания матрицы А, вызванного наличием нестационарного члена.

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

'(О, у) =

Гю0у(0.2-и), при у <0.2, \0, при у > 0.2,

х(1, у) = и(х, 0) = и(х, 1) = 0,

и(х,1) =

100(х - 0.8)(1 - х), при х > 0.8, 0, при х < 0.8,

«(О,») = «(!,») =в(а;,0) = 0.

Тогда, с учётом декомпозиции давления, уравнения вспомогательной задачи принимают вид:

а) уравнение движения по X и уравнение постоянства массового расхода

' д{и2) д(уи) _ _ Гдрху дх ду (1х дх

1 1 X

/и(х,у)с1у= ¡и(0,у) йу- у) ооо

(сРи + Ре + ду2

в) уравнение движения по У и уравнение постоянства массового расхода

С д{иу) д(у2) _ йри

дх ду йу

1 у

/ у{х,у)йх = /и(0, о о

дрх

ду

+

Ре ^х2 + ду2)'

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

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

Сходимость многосеточных итераций при решении стационарных уравнений Навье-Стокса показана на рис. 21, где ||г„||, ||г„||, ||ги„||, ||г*|| и ||г*|| есть || • Цдо нормы невязок линеаризованных разностных аналогов уравнений движения по X, У, неразрывности и постоянств массового расхода соответственно. Метод Зейделя с блочным упорядочиванием неизвестных и метод Банки использованы в качестве сглаживающих процедур при решении вспомогательной и основной задач. Критерий останова при решении вспомогательной задачи задан в виде

тах(||гц||, |К||) < Ю-3, тах(|К|[, |К||) < Ю~в, а при решении основной задачи определяется как

тах(М,|М,|Ы|)<10-4.

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ

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

2. Создан программный комплекс для решения двух- и трёхмерных краевых задач теплопроводности и гидродинамики на структурированных сетках.

Рис. 18. Изолинии функции тока: решение вспомогательной (слева) и основной (справа) задал при Ре = 100

Рис. 19. Изобары: решение вспомогательной (рх + ри, слева) и основной (рх + рУ +рху} справа) задач при Ре = 100

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

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

Рис.. 20. Многосеточный цикл для решения уравнений Навье-Стокса

Ке=100. 101x101 10*'

—к« 10"

—м 10*'

—111,11 10*!

—иш 10*'

II СИ 10*°

10"'

10!

ю-3

10"4

ю'

ю"

10-'

Ис-ЮО. 1001x1001

2 4 6 8

номер многосеточной итерации

1(е»500, 101x101 10*4

— I1111 10й

—I11.11 10«

-l1.il 10*'

—о— 11г"|!

И ыII 10'°

V—I11||

2 4 6

номер многоееточнои итерации

4 6 8 10 12 14 номер шюгосеточпой итерации

Рис. 21. Сходимость многосеточных итераций при решении стационарных уравнений Навье-Стокса

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ ОТРАЖЕНЫ В РАБОТАХ

1. Кунбутаев Л.М., Мартыненко С.И. Упрощенная двумерная математическая модель сопряженного теплообмена при кипении криогенных жидкостей в электрически обогреваемых трубных каналах // Вестник МЭИ. 1998. Вып. 5. С.16-21.

2. Кунбутаев Л.М., Мартыненко С.И., Матлин Г.Г. Численное моделирование вынужденной двухфазной конвекции криогенных жидкостей в трубных каналах // Вестник МЭИ. 1998. Вып. 4. С.77-88.

3. Мартыненко С. И. Универсальная многосеточная технология для численного решения краевых задач на структурированных сетках // Вычислительные методы и программирование. 2000. Т.1, раздел 1. С.85-104.

4. Мартыненко С.И. Численное решение уравнений Навье-Стокса при помощи универсальной многосеточной технологии // Optimization of Finite Element Approximations & Splines and Wavelets: материалы докладов Международной конференции. СПб., 2001. C.44-4G.

5. Мартыненко С.И. Универсальная многосеточная технология для численного решения систем дифференциальных уравнений в частных производных // Вычислительные методы и программирование. 2001. Т.1, раздел 1. С.1-11.

6. Мартыненко С.И. Программное обеспечение для универсальной многосеточной технологии: строительные блоки и диагностические инструменты // Вычислительные методы и программирование. 2001. Т.4, раздел 1. С.1-6.

7. Мартыненко С.И. Программное обеспечение для универсальной многосеточной технологии // Математическое моделирование. 2002. Т. 14, №9. С.87-90.

8. Мартыненко С.И. Распараллеливание универсальной многосеточной технологии // Вычислительные методы и программирование. 2003. Т.4., раздел 1. С.45-51.

9. Мартыненко С.И. Робастные многосеточные методы: проблемы и перспективы развития // Сеточные методы для краевых задач и приложения: материалы V Всероссийского семинара. Казань, 2004. С. 150-154.

10. Мартыненко С.И. Основные принципы построения вычислительных технологий для перспективного программного обеспечения // Сеточные методы для краевых задач и приложения: материалы VI Всероссийского семинара. Казань, 2005. С. 170-173.

11. Мартыненко С.И. Анализ погрешностей вычислений при численном

решении уравнений Навье-Стокса // Необратимые процессы в природе и технике: тез. докл. III Всероссийской конференции. М. 2005. С. 141-142.

12. Martynenko S.I. Robust multigrid technique for black box software // Comp. Meth. in Appl. Math. 2006. V.6, N 4. P.413-435.

13. Мартыненко С.И. Универсальная многосеточная технология как обобщение геометрических многосеточных методов // Численная геометрия, построение расчетных сеток и высокопроизводительные вычисления: материалы Всероссийской конференции. М., 2006. С.131-140.

14. Мартыненко С.И., Яновский JI.C. Алгоритм с неформальной сегрегацией вычислений для решения уравнений Навье-Стокса на структурированных сетках // Струйные, отрывные и нестационарные течения: тез. докл. XXI Всероссийского семинара. Новосибирск, 2007. С.155-157.

15. Мартыненко С.И. Алгоритм с неформальной сегрегацией вычислений для численного решения уравнений Навье-Стокса на структурированных сетках // Сеточные методы для краевых задач и приложения: материалы VII Всероссийского семинара. Казань, 2007. С.181-185.

16. Мартыненко С.И. Совершенствование алгоритмов для решения уравнений Навье-Стокса в переменных «скорость-давление» // Ученые записки Казанского государственного университета. Сер.: Физ-матем. науки. 2007. Т. 149, кн. 4. С.105-120.

17. Мартыненко С.И. Совершенствование методов численного решения уравнений Навье-Стокса в примитивных переменных // Необратимые процессы в природе и технике: труды IV Всероссийской конференции. М., 2007. Ч. 1. С. 279-282.

18. Martynenko S.I. Development of numerical methods for Navier-Stokes equations in primitive variables // Numerical geometry, grid generation and high performance computing: proceedings of the Int. Conf. M., 2008. P.75-78.

19. Мартыненко С.И. Совершенствование вычислительных алгоритмов для решения уравнений Навье-Стокса на структурированных сетках // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2008. №2. С.78-94.

20. Мартыненко С.И. Адаптация уравнений Навье-Стокса к универсальной многосеточной технологии // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2008. Вып. 3. С.75-80.

21. Мартыненко С.И. Формализация вычислений при численном решении краевых задач // Ученые записки Казанского госу-

дарственного университета. Физ.-матем. науки. 2008. Т. 150, kii. 1. С. 76-90.

22. Бакулин В.Н., Мартыненко С.И. Универсальная многосеточная технология для численного решения задач механики сплошной среды //VI Окунев-ские чтения: материалы докладов Международной конференции. СПб., 2008. С.47-50

23. Мартыненко С.И., Ольшанский М.А. Многосеточные методы // Энциклопедия низкотемпературной плазмы. Серия Б. Справочные приложения, базы и банки данных. Т. VII - 1, ч. 1. Математическое моделирование в низкотемпературной плазме / Под ред. Ю.П. Попова. М.: Янус-К, 2008. С.23-37.

24. Martynenko S.I. , Markov А.А., Sliarov M.S., Yanovskiy L.S. Development of numerical methods for solving Navier-Stokes equations in primitive variables formulation // Авиация в XXI столетии: труды III конгр. Киев, 2008. С.92-95.

25. Мартыненко С.И. Универсальная многосеточная технология // Математическое моделирование. 2009. Т.21, №9. С.66-79.

26. Бакулин В.Н., Мартыненко С.И. Универсальная многосеточная технология для решения задач механики сплошной среды // Вестник МАИ. 2009. Т.16, №2. С.171-175.

27. Martynenko S.I. A physical approach to development of numerical methods for solving Navier-Stokes equations in primitive variables formulation // Int. J. of Сотр. Science and Math. 2009. V.2, N 4. P.291-307.

28. Yanovskaya M.L., Martynenko S.I. A physical approach for development of computational algorithms for solving the Navier-Stokes equations and its application in jet engine analysis // Mathematical Modelling and Computational Physics: proceeding of Int. Conf. Dubna, 2009. P.38

29. Мартыненко С.И. Применение универсальной многосеточной технологии к численному решению уравнений Навье-Стокса // Необратимые процессы в природе и технике: труды V Всероссийской конференции. М., 2009. 4.2. С.140-142.

30. Мартыненко С.И. Замечания о вычислении давления при численном решении уравнений Навье-Стокса // Математическое моделирование. 2010. Т.22, №3. С.105-119.

31. Martynenko S.I. Potentialities of the robust multigrid technique // Com. Meth. in Appl. Math. 2010. V.10, N 1. P.87-94.

32. Проблемы разработки микродвигательных установок / С.И. Мартыненко [и др.] // Изв. вузов. Авиационная техника. 2010. №2. С.53-55.

33. Мартыненко С.И., Вэйсин Чжоу. Повышение вычислительной эффективности алгоритмов для решения уравнений Навье-Стокса применительно к микродвигателям // Прямоточные ВРД и химмотология: сборник научных трудов ЦИАМ им. П.И. Баранова. М., 2010. С.75-82

34. Мартыненко С.И. К вопросу о сходимости универсальной многосеточной технологии // Математическое моделирование. 2010. Т.22, №10. С.18-34.

35. Мартыненко С. И. Исследование сходимости универсальной многосеточной технологии // Сеточные методы для краевых задач и приложения: материалы VIII Всероссийской конференции, посвященной 80-летию со дня рождения А.Д. Ляшко. Казань, 2010. С.278-282.

36. Вэйсин Чжоу, Шаров М.С., Мартыненко С.И. Применение универсальной многосеточной технологии к численному решению уравнений Навье-Стокса // Сеточные методы для краевых задач и приложения: материалы VIII Всероссийской конференции, посвященной 80-летию со дня рождения А.Д. Ляшко. Казань, 2010. С.391-394.

37. Мартыненко С.И. Особенности моделирования рабочих процессов в авиационных системах охлаждения канального типа // Авиадвигатели XXI века: материалы конф. М., 2010. С.409-410.

38. Мартыненко С.И. Совершенствование методов математического моделирования процессов гидродинамики и теплообмена при помощи априорной информации физического характера // Труды V Российской национальной конференции по теплообмену. М., 2010. T.I. С.93-96.

39. Мартыненко С.И. Замечания о применении явных схем для численного решения уравнений Навье-Стокса // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2011. №2. С.107-120.

40. Мартыненко С.И. Оценки эффективности распараллеливания универсальной многосеточной технологии // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2011. №4. С.63-80.

41. Martynenko S. Convergence acceleration of iterative algorithms for solving Navier-Stokes equations on structured grids // Hydrodynamics - Optimizing Methods and Tools. Croatia, 2011. P.175-200

42. Мартыненко С.И. К вопросу о доказательстве сходимости многосеточных методов // Ученые записки Казанского государственного университета. Физ.-матем. науки. 2012. Т. 154, кн. 4. С. 5-16.

Подписано в печать:

15.04.2013

Заказ № 8355 Тираж -100 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www. autoreferat. ru

Федеральное государственное унитарное предприятие «Центральный институт авиационного моторостроения имени П.И. Баранова»

На правах рукописи

05201351103

Мартыненко Сергей Иванович

МНОГОСЕТОЧНАЯ ТЕХНОЛОГИЯ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ТЕПЛОВЫХ И ГИДРОДИНАМИЧЕСКИХ ПРОЦЕССОВ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

Диссертация на соискание учёной степени доктора физико-математических наук

Москва - 2013 г.

СОДЕРЖАНИЕ

Введение б

1 Универсальная многосеточная технология 26

1.1 Постановка задачи о формализации вычислений....................26

1.2 Аналитическая часть технологии......................................28

1.3 Вычислительная часть технологии....................................29

1.3.1 Построение многосеточной структуры........................30

1.3.2 Аппроксимация краевых задач на многосеточной структуре 36

1.3.3 Аппроксимация граничных условий..........................40

1.3.4 Многосеточный цикл............................................44

1.4 Иллюстративный пример..............................................46

1.5 Доказательство сходимости............................................54

1.6 Распараллеливание вычислений ......................................64

1.6.1 Минимальное ускорение........................................73

1.6.2 Динамический цикл............................................75

1.6.3 Распараллеливание на 27 процессорах........................77

1.6.4 Распараллеливание на 729 процессорах......................79

1.6.5 Влияние обмена данными......................................80

1.7 Двухмерные модельные задачи........................................85

1.7.1 Уравнение Пуассона............................................85

1.7.2 Уравнение с квадратичной нелинейностью..................90

1.7.3 Анизотропное уравнение........................................93

1.8 Трехмерные модельные задачи........................................95

1.8.1 Аппроксимация уравнения с переменными коэффициентами 95

1.8.2 Уравнение Пуассона............................................97

1.8.3 Уравнение с переменными коэффициентами ................99

1.8.4 Нелинейное уравнение теплопроводности..........105

1.8.5 Уравнение с разрывными коэффициентами.........110

1.8.6 Неравномерные сетки.....................111

1.8.7 Сильные сглаживатели ....................117

1.9 Моделирование передачи тепла в корпусе микродвигателя .... 120

1.10 Многосеточное предобуславливание.................124

1.11 Адаптивное определение параметров релаксации..........133

1.12 Сравнение многосеточных методов..................139

2 Комплекс программ 145

2.1 Модуль RMT_2D.f90 ..........................145

2.1.1 Построение и хранение многосеточной структуры.....147

2.1.2 Вычисление отображения индексов..............151

2.1.3 Вычисление интегралов на многосеточной структуре . . . 153

2.1.4 Многосеточная оболочка....................165

2.2 Модуль RMT_3D.f90 ..........................172

2.2.1 Многосеточная оболочка....................175

2.3 Замечания о программном обеспечении...............176

3 Совершенствование алгоритмов для численного решения уравнений Навье-Стокса 177

3.1 Декомпозиция давления........................177

3.2 Явные схемы..............................185

3.2.1 Схема расщепления по физическим факторам.......185

3.2.2 Метод искусственной сжимаемости..............189

3.2.3 Схема Мак-Кормака......................194

3.3 Неявные схемы.............................196

3.3.1 Течение в каверне с движущейся крышкой.........196

3.3.2 Обтекание ступеньки......................200

3.3.3 Течение в микрокатализаторе.................203

3.3.4 Течения между пластинами с локальным сужением .... 205

3.3.5 Течение в микросопле Лаваля.................207

3.4 Течения с определяемым массовым расходом............215

3.5 Замечания о декомпозиции давления.................218

4 Многосеточный алгоритм для численного решения уравнений Навье-Стокса 220

4.1 Модификация уравнений Навье-Стокса...............220

4.2 Конфигурация контрольных объёмов................222

4.3 Аппроксимация модифицированных уравнений...........223

4.3.1 Аппроксимация уравнения неразрывности .........223

4.3.2 Аппроксимация уравнения движения ............224

4.4 Модификация вспомогательной задачи ...............231

4.5 Многосеточные циклы.........................233

4.6 Предобусловленный метод Узавы...................236

4.7 Сглаживатель Банки..........................247

Выводы 253

Список литературы 254

ОБОЗНАЧЕНИЯ

с - поправка, вычисляемая на многосеточной структуре (/) ~ среднеинтегральное значение функции / д - номер многосеточной итерации

- количество многосеточных итераций I - номер сеточного уровня

Ь+ - номер сеточного уровня с самыми грубыми сетками N - число пространственных измерений (./V = 1,2,3) N - число узлов сетки Ке - число Рейнольдса

V - контрольный объем

и^ - оптимальное значение параметра релаксации УУ - число арифметических операций (а.о.)

V - номер сглаживающей итерации

- количество сглаживающих итераций

- эквивалентное число сглаживающих итераций Д - оператор Лапласа

СОКРАЩЕНИЯ

ВСС - вспомогательная структурированная сетка

КММ - классические многосеточные методы

МКЭ - метод конечных элементов

МКР - метод конечных разностей

СЛАУ- система линейных алгебраических уравнений

НСС - неструктурированная сетка

УМТ - универсальная многосеточная технология

РСС - предобусловленный метод сопряженных градиентов

ВВЕДЕНИЕ

Актуальность темы. В настоящее время методы математического моделирования широко используются для решения крупных научно-технических и хозяйственных задач. Однако математическому моделированию, как и любому методу исследования, присущи определённые недостатки. В первую очередь к ним следует отнести трудоемкость написания и отладки компьютерных программ, предназначенных для решения сложных прикладных задач. Поэтому широкое распространение получил ряд коммерческих программных продуктов (ANSYS, STAR-CD, FLUENT, CFX, PHOENICS и др.), устроенных по принципу «черного ящика». Инженер, использующий подобные программы, только формулирует задачу, а детали вычислительного алгоритма ему могут быть даже не известны. Применительно к техническим приложениям работу автономных программ упрощенно можно представить следующим образом: конструктор проектирует некоторый узел при помощи некоторой графической программы. Затем геометрия узла переносится в вычислительный модуль, конструктор задает граничные и начальные условия, после чего проводит тепловой, прочностной, гидродинамический или другой расчет и анализирует результаты. Как правило, после анализа полученных результатов нужно внести изменения в конструкцию и повторить расчет. Обычно конструктор выполняет несколько подобных «итераций», чтобы получить оптимальную, с его точки зрения, конструкцию. Еще бблыиую практическую ценность будет представлять возможность расчета машины в целом, а не только отдельных ее узлов, поскольку поэлементное моделирование часто связано с погрешностями постановки граничных условий.

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

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

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

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

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

сти аппроксимации. Адаптивная перестройка сетки требует достаточно точной апостериорной оценки погрешности.

3. Выбор способа и порядка аппроксимации уравнений математической модели. Метод конечных разностей (МКР) и метод конечных элементов (МКЭ) стали самыми распространенными и эффективными способами аппроксимации различных задач математической физики. Однако для придания гибкости автономным программам необходимо предусмотреть возможность изменения способа и/или порядка аппроксимации в процессе решения задачи.

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

В настоящее время основная тенденция в развитии программного обеспечения, устроенного по принципу «чёрного ящика», состоит в разделении основных этапов вычислительного эксперимента на отдельные независимые задачи. Например, в алгебраических многосеточных методах решение результирующей СЛАУ осуществляется без привлечения какой-либо информации о вычислительной сетке, исходной дифференциальной задаче и способе ее аппроксимации [81, 84, 104, 105, 156]. При этом алгоритмы для решения отдельных независимых задач содержат много компонент, оптимальный выбор которых и определяет скорость сходимости. Адаптацию алгоритма к решаемой задаче трудно выполнить в программном обеспечении, устроенном по принципу «чёрного ящика».

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

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

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

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

Актуальность работы обусловлена необходимостью разработки универсального и высокоэффективного алгоритма для численного решения широкого класса (не)линейных прикладных задач на структурированных сетках.

Состояние вопроса исследования. Подавляющее большинство задач вычислительной математики сводится к решению систем линейных алгебраической уравнений высокой размерности. В настоящее время предложено огромное количество алгоритмов решения задач линейной алгебры, однако доминирующую роль получили многосеточные методы, предложенные выдающими отечественными математиками Р.П. Федоренко [63, 64] и Н.С. Бахваловым [3]. Бурное развитие многосеточных методов, начавшееся в середине 80-х годов, обусловлено их уникальной возможностью получать численное решение с оптимальными вычислительными усилиями. Развитие многосеточных методов пошло по тра-

диционному для середины 80-х годов пути, связанным с адаптацией алгоритма к решаемой задаче. В результате многосеточные методы быстро превратились в труднообозримое семейство алгоритмов, причем применение их к решению новых задач приводило к появлению новых вариантов [67, 80, 97,100,155,162,165].

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

Первоначально классические многосеточные методы (КММ) применялись к решению линейных задач. При решении нелинейных уравнений КММ можно использовать тремя способами:

1) линеаризация по Ньютону с адаптацией числа многосеточных итераций на каждую итерацию Ньютона;

2) линеаризация по Ньютону с фиксированным числом многосеточных итераций на каждый шаг по Ньютону;

3) нелинейный многосеточный метод.

Наибольшее распространение для решения отдельных нелинейных задач получила схема FAS [62, 97, 162], к преимуществам которой следует отнести:

1) избежание глобальной линеаризации;

2) отсутствие расчёта больших Якобианов;

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

Тем не менее, трудности решения нелинейных задач при помощи КММ остаются для ряда приложений.

Другим из наиболее распространённых нелинейных многосеточных алгоритмов стали каскадные многосеточные методы [77, 149, 152], которые сочетают

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

Другим принципиальным недостатком классических многосеточных методов (КММ) является трудность их распараллеливания. КММ основаны на решении серии разностных задач различной размерности, причем наименьшая размерность может быть меньше числа используемых процессоров. Это неизбежно приводит к снижению эффективности распараллеливания. В настоящее время построение параллельных классических многосеточных алгоритмов осуществляется по следующим направлениям: совпадающие итерации (одновременное выполнение сглаживающих итераций на всех сетках [73, 74, 82, 107, 108]), коррекция на нескольких грубых сетках (использование нескольких сеток для вычисления поправки [69, 89, 90, 96, 99, 133]), полная декомпозиция области (уменьшение обменов данными в каждой многосеточной итерации [70, 131, 132]), и блочная факторизация (специальный выбор точек на грубых и мелких сетках [78, 79, 88, 122]). Поэтому иная форма алгоритмизации основополагающей идеи Р.П. Федоренко и Н.С. Бахвалова должна подразумевать возможность эффективного распараллеливания вычислений на достаточно большом количестве процессоров вне зависимости от используемой сглаживающей процедуры.

Альтернативный подход, основанный на использовании информации об исходной дифференциальной задаче, вычислительной сетке и способе аппроксимации, разработан автором и получил название «Универсальная Многосеточная Технология» (УМТ) [32, 33, 34, 35, 36, 39, 42, 126, 128]. Взаимосвязь УМТ с исходной дифференциальной задачей осуществляется предварительным представлением искомого решения в виде суммы или произведения двух функций,

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