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

РЕЙН ТАТЬЯНА СЕРГЕЕВНА

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ВЯЗКОЙ

НЕСЖИМАЕМОЙ ЖИДКОСТИ СО СВОБОДНЫМИ ГРАНИЦАМИ

ОБОБЩЕННЫМ МЕТОДОМ ЕСТЕСТВЕННЫХ СОСЕДЕЙ

05.13.18 – математическое моделирование, численные методы

и комплексы программ

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

Кемерово 2008

Работа выполнена в ГОУ ВПО «Кемеровский государственный университет»

на кафедре ЮНЕСКО по НИТ доктор физико-математических наук,

Научный руководитель:

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

Официальные оппоненты:

профессор Хакимзянов Гаяз Салимович доктор физико-математических наук, профессор Воеводин Анатолий Федорович Томский государственный университет

Ведущая организация:

Защита состоится 10 июня 2008 г. в 10-30 на заседании диссертационного совета ДМ 003.046.01 при Институте вычислительных технологий СО РАН по адресу: 630090, г. Новосибирск, пр. Ак. Лаврентьева, 6.

С диссертацией можно ознакомиться в специализированном читальном зале вычислительной математики и информатики ГПНТБ СО РАН.

Автореферат разослан 8 мая 2008 г.

И.о. учёного секретаря диссертационного совета доктор технических наук, А.Д. Рычков профессор

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

Актуальность темы Появление вычислительных машин в 60-х годах прошлого столетия стимулировало развитие вычислительных методов в естественных науках, инженерных дисциплинах и в управлении. Появление персональных ЭВМ на рубеже 80–90-х годов заметно ускорило процессы разработки новых алгоритмов и математических моделей.

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

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

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

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

С ростом производительности компьютеров развитие получили бессеточные методы, которые аппроксимируют уравнения в частных производных, основываясь только на наборе узлов, без знания дополнительной информации о структуре сетки. В таких методах отношение соседства частиц не фиксировано и может со временем изменяться, то есть частицы, бывшие соседями в начальный момент времени, могут со временем расходиться достаточно далеко друг от друга. Характерными представителями этой группы методов являются метод сглаженных частиц (SPH – Smoothed Particle Hydrodynamics), полунеявный метод движущихся частиц (MPS – Moving Particle Semi-implicit), метод Лагранжево-Эйлеровых частиц, метод точечной интерполяции (PIM – Point Interpolation Method). Данные методы позволяют достаточно точно воспроизводить кинематику течений, однако полученные динамические характеристики, необходимые для расчета гидродинамических нагрузок, являются неточными. К общим недостаткам бессеточных методов также можно отнести и сравнительно невысокую точность, и трудность введения граничных условий.

Эти обстоятельства заставили исследователей искать новые методы, сочетающие в себе идеи и возможности бессеточного подхода, но вместе с тем обладающие достоинствами сеточных методов. Первыми из бессеточных методов нового поколения появились бессеточный метод конечных элементов (MFEM – Meshless Finite Element Method) и метод естественных соседей (NEM – Natural Element Method). Особенность методов NEM и MFEM в том, что для стационарных задач они являются обычными (классическими) методами Галеркина, то есть являются сеточными. Для нестационарных задач, в которых применяется Лагранжев подход к описанию изучаемого процесса, на каждом шаге по времени по найденному на предыдущем шаге положению узлов строится новая сетка, определяющая новую структуру соседей для каждой узловой точки области. На вновь построенной сетке аппроксимированная система уравнений снова решается методом Галеркина. В силу этого методы NEM и MFEM сохраняют некоторые преимущества классического метода Галеркина, а именно простоту функций формы в области определения, непрерывность между элементами, легкость введения граничных условий. При этом имеют все достоинства бессеточных методов, так как функции формы метода естественных соседей зависят только от положения узловых точек.

Для формирования дискретной системы уравнений используется метод взвешенных невязок с набором весовых функций, совпадающих с базисными. Интегралы берутся по элементам расширенной триангуляции Делоне1. Множество естественных соседей для каждого узла, а также узлы свободной границы на новом временном шаге определяются с помощью методов «sweep-line» и « -shape»2. Для аппроксимации неизвестных функций используются функции формы Сибсона3 и Лапласа4. Полученная система линейных алгебраических уравнений после внедрения граничных условий решается методом сопряженных градиентов с предобусловливанием.

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

Задачи исследования 1. Разработка алгоритма построения расширенной триангуляции Делоне на основе разбиения расчетной области ячейками Вороного первого порядка.

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

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

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

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

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

Farin G. Surfaces over Dirichlet tessellations // Computer Aided Geometric Design. – 1990. – Vol. 7. – Р. 281–292.

Fortune S.J. A sweepline algorithm for Voronoi diagrams // Journal Algorithmica. – 1987. – № 2. – P. 153–174.

Sibson R. A brief description a natural neighbor interpolation / R. Sibson, V. Barnett (ed.) // Interpret multivariate data. – Chichester: John Wiley, 1981. – P. 21–36.

Несибсоновская интерполяция – новый метод интерполяции значений функции на произвольной системе точек / В.В. Беликов, В.Д. Иванов, В.К. Конторович и др. // Вычислительная математика и математическая физика. – 1997. – Т. 37, № 1. – С. 11–17.

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

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

3. Проведены в полной нелинейной постановке численные эксперименты по расчету задач об обрушении плотины в зависимости от варьируемых параметров.

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

4. Разработана параллельная реализация обобщенного метода естественных соседей.

На защиту выносятся:

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

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

3. Результаты численного моделирования задачи об обрушении плотины.

4. Параллельная реализация обобщенного метода естественных соседей.

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

Основные результаты работы были использованы при выполнении следующих проектов:

проекта № 4829 «Численное моделирование течений жидкости со свободными границами современными численными методами на многопроцессорных вычислительных системах» (2005 год) по ведомственной научной программе Федерального агентства по образованию «Развитие научного потенциала высшей интеграционного проекта фундаментальных исследований Объединенного ученого совета по механике и энергетике СО РАН (2006–2008 годы) по теме «Численное моделирование нестационарного взаимодействия сложных упругих конструкций с жидкостью или газом», блок 2: «Нестационарное взаимодействие нелинейных поверхностных волн с плавающими и закрепленными упругими конструкциями», Пункт 1. «Развитие методов расчета гидродинамических нагрузок при резко нестационарном воздействии волн с большими деформациями области течения»;

проекта № 4256 «Создание типового информационно-вычислительного портала для организации учебной и научной деятельности вуза» по ведомственной научной целевой программе Федерального агентства по образованию «Развитие научного потенциала высшей школы (2006–2008 годы)» (2006–2008 годы).

Представление результатов. Основные результаты диссертации докладывались на: III Всероссийской научно-практической конференции «Информационные технологии и математическое моделирование» (Анжеро-Судженск, 2004); региональной научно-практической конференции студентов, аспирантов и молодых ученых, посвященной десятилетию Новокузнецкого филиала-института Кемеровского государственного университета (Новокузнецк, 2005); Всероссийской научно-практической конференции «Недра Кузбасса. Инновации» (Кемерово, 2006); Международной научной конференции «Наука и образование» (Белово, 2006); IX Международной летней научной школе «Гидродинамика больших скоростей и численное моделирование»

(Кемерово, 2006); II и III Российско-Германской школе по параллельным вычислениям на высокопроизводительных вычислительных системах (Новосибирск, 2006); VII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (с участием иностранных ученых) (Красноярск, 2006);

VI Международной научно-практической конференции «Инновационные недра Кузбасса. IT-технологии» (Кемерово, 2007); Международной конференции «Сопряженные задачи механики реагирующих сред, информатики и экологии» (Томск, 2007);

VII Всероссийской научно-практической конференции «Информационные недра Кузбасса. IT-технологии» (Кемерово, 2008); а также на объединенном семинаре ИВТ СО РАН «Информационно-вычислительные технологии (численные методы механики сплошной среды)» под руководством академика РАН Шокина Ю.И., профессора Ковени В.М. (Новосибирск, декабрь 2007); на научном семинаре Института гидродинамики им. М.А. Лаврентьева СО РАН «Прикладная гидродинамика» под руководством чл.-кор. РАН Пухначева В.В. (Новосибирск, февраль 2008) и на научном семинаре «Информационные технологии и математическое моделирование» под руководством профессора Афанасьева К.Е. (Кемерово, 2004–2008).

Публикации По теме диссертации опубликовано 13 работ, в том числе (в скобках в числителе указан общий объём этого типа публикаций, в знаменателе – объём, принадлежащий лично автору) 2 статьи в изданиях, рекомендуемых ВАК для предоставления основных результатов диссертации (1,44/0,71 печ. л.), 4 публикации в трудах и материалах конференций (1,65/1,4 печ. л.), 7 публикаций в тезисах конференций (0,9/0,82 печ. л.).

Личный вклад автора. Основные научные и практические результаты диссертации получены автором лично. В работе [1] автору принадлежит методика численного исследования задач динамики вязкой несжимаемой жидкости. В работах [2, 3] автором были проведены численные расчеты и реализован алгоритм построения функций формы Сибсона. В [4] автор участвовал в постановке задачи и разработке алгоритмов разбиения расчетной области ячейками Вороного. В работах [5, 6] автору принадлежит постановка задачи и анализ результатов. В работе [8] автору принадлежит реализация параллельных алгоритмов построения функций формы и формирования матрицы СЛАУ. В работе [9] автором были проведены расчеты модельных задач, а также сравнение полученных результатов с экспериментальными данными. Автор принимал участие в подготовке и представлении докладов на конференциях.

Структура и объём работы. Диссертационная работа состоит из введения, четырех глав, списка цитируемой литературы и приложения. Общий объём работы составляет 180 страниц машинописного текста, включая приложение – 10 страниц; библиографический список состоит из 156 литературных источников.

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

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

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

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

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

Второй параграф посвящен построению интерполяционных функций формы Сибсона. Рассматриваются алгоритмы решения задачи о вычислении в заданной точке x 0 = (x 0, y 0 ) значения некоторой функции f (x) по ее значениям {fk = f (x k, yk )}, заданным на фиксированной системе точек {xk }.

Карабцев, С.Н., Стуколов С.В. Эффективный алгоритм генерации конечноэлементной сетки для метода естественных соседей уколов // Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». – Кемерово: ИНТ, 2006. – С. 401–409.

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

где A(x) – площадь ячейки Вороного первого порядка, в которую входит точка x, Ai (x) – площадь пересечения ячейки Вороного второго порядка точки x и площадью A(x). Для вычисления функций формы Сибсона необходимо искать пересечение ячеек Вороного первого и второго порядков, что является весьма трудоемкой задачей вычислительной геометрии. В настоящей работе рассматривается предложенный Д.Уотсоном6 алгоритм поиска площадей пересечения многоугольников путем разбиения многоугольника ориентированными треугольниками.

Третий параграф связан с построением интерполяции Лапласа (несибсоновской интерполяции), которое также опирается на определение соседей посредством разбиения области ячейками Вороного.

Пусть точка x = (x, y ) принадлежит многоугольнику Вороного с числом сторон, равным n. Обозначим длины сторон многоугольника через si, i = 1, n, а высоты, опущенные из точки x на si, – через hi. Тогда интерполяционные коэффициенты Лапласа будут иметь следующий вид:

Такой способ определения коэффициентов i проще и экономичнее, чем в подходе Сибсона, так как не требует вычисления площадей пересечения многоугольников. Методика построения функций формы Лапласа изложена в приложении 1.2.

Производные функций формы Сибсона и Лапласа можно получить дифференцированием выражений (1) и (2) соответственно.

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

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

Watson D.F. Natural neighbor sorting on the N-dimensional sphere // Pattern recognition. – 1988. – Vol. 21, № 1. – P.63–67.

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

Тесты проводятся в областях с различной геометрией и с различным числом расчетных узлов. Рассматривается область D, верхняя граница 1 которой задается уравнением y = 0, 5 Sin(x ), x [0,2 ], а боковые и нижняя границы 2, 3, 4 являются прямыми линиями: 2 : x = 0, y [1, 0] ; 3 : y = 1, x [0,2 ] ;

В описанной выше области решается уравнение Пуассона u = b для функции u = x y 2 с граничными условиями Дирихле на 1 и Неймана на границах В качестве численного метода решения уравнения Пуассона используется метод Галеркина в слабой форме. Численное интегрирование осуществляется по элементам расширенной триангуляции Делоне с помощью квадратур Хаммера. Полученная система линейных алгебраических уравнений решается методом сопряженных градиентов с предобусловливанием.

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

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

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

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

Пусть в некоторой области течения D происходит движение ньютоновской вязкой несжимаемой жидкости, описываемое системой уравнений Навье-Стокса:

Здесь i = 1,2, j = 1,2, x(t ) = (x 1(t ), x 2 (t )), по повторяющемуся индексу j производится суммирование. В системе уравнений (3)–(4) искомыми функциями являются давление p(x, t ) и вектор скорости u = (u1(x, t ), u2 (x, t )), параметрами – плотность, коэффициент динамической вязкости и вектор массовых сил f = ( f1, f2 ).

Так как жидкость вязкая, то на твердых стенках Гu выполняется условие прилипания: ui = ui, i = 1,2. На свободной поверхности Г выполняется динамическое условие p = patm. Для нестационарной задачи задается положение расчетных узлов x(0) = x 0 и распределение неизвестных функций во всей области течения:

u(x, 0) = u 0 (x), p(x, 0) = p 0 (x).

Во втором параграфе описывается метод расщепления для интегрирования по времени системы уравнений Навье-Стокса7. В случае, когда жидкость несжимаемая, не может быть использован явный метод интегрирования по времени, но даже при использовании неявного метода не удается устранить нефизические колебания функции давления. Суть метода расщепления заключается в разбиении описываемого физического процесса на два: конвекцию-диффузию и вклад давления. На первом этапе в уравнении движения учитываются только конвективные члены, в результате чего выделяется фиктивная переменная u (x, t ) и записываются выражения для предиктора и корректора скорости:

Здесь [u(x, t )] На втором этапе решается уравнение Пуассона для давления:

Основной алгоритм движения по времени состоит из следующих шагов:

I) определяются границы области и строятся интерполяционные функции;

II) вычисляется предиктор скорости u (x, t ) из системы (5);

III) решается уравнение Пуассона (7) для определения давления p n +1(x, t ) ;

IV) вычисляется новое значение скорости un +1(x, t ) из уравнения (6) с учетом найденного на шаге III давления;

V) вычисляется новое положение узлов на (n + 1) -м временном шаге:

В третьем параграфе описывается переход к слабой форме уравнений, полученных после расщепления основных уравнений, а именно: предиктора и корректора скорости (5)–(6), а также уравнения Пуассона (7). Для аппроксимации неизвестных функций используются функции формы Сибсона и Лапласа.

Марчук Г.И., Яненко Н.Н. Применение метода расщепления (дробных шагов) для решения задач математической физики // Докл. на Всесоюз. конф. по вычисл. матем., Москва, февраль 1965. Докл. на конгр. ИФИП, НьюЙорк, май 1965.

Одной из главных трудностей численного моделирования нестационарных уравнений Навье-Стокса является регуляризация условия несжимаемости (4). В описанном выше методе расщепления по пространственным переменным условие несжимаемости представляется уравнением Пуассона для давления (7). Для устранения нефизических осцилляций функции давления используется метод конечных приращений (Finite Increment Calculus8), основанный на понятиях дифференциального приближения9. При описании метода конечных приращений для стабилизации уравнений Навье-Стокса выписываются дифференциальные приближения уравнений движения и неразрывности, а также выражения для вычисления характеристических параметров метода.

Устойчивость решения системы уравнений Навье-Стокса методами, основанными на методе Галеркина, обеспечивается выбором конечно-элементных пространств для скорости и давления: степени интерполяционных полиномов компонент вектора скорости и давления должны удовлетворять условию Ладыженской-БабушкиБреззи (ЛББ). В данной работе для аппроксимации уравнения неразрывности используются линейные базисные функции (функции формы расширенной интерполяции Лапласа), для уравнения движения применяются квадратичные базисные функции (функции формы Сибсона). Построение так называемого обобщенного метода естественных соседей (GNEM – General Natural Element Method), гибрида методов MFEM и NEM приводит к удовлетворению условий ЛББ для совместной аппроксимации, что гарантирует невырожденность решения.

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

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

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

Рассматриваются следующие задачи:

ламинарное течение в плоском канале (течение Пуазейля) – задача с твердыми границами;

задача Л.В. Овсянникова о деформации жидкого эллипса – задача со свободной границей;

движение вязкой жидкости в квадратной каверне с движущейся крышкой;

колебания жидкости в прямоугольном бассейне;

распределение поля давления в покоящейся жидкости конечной глубины.

Onate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation // Comput. Meth. Appl. Mech. Eng. – 2000. – Vol. 182, № 1–2. – Р. 355–370.

Шокин Ю.И. Метод дифференциального приближения. – Новосибирск: Наука, 1979. – 224 с.

Задача Л.В. Овсянникова10 о деформации жидкого эллипса была выбрана для тестирования алгоритма движения по времени. Давление на границе круга задается постоянным: p = 0. Общая постановка задачи, как правило, приводится для идеальной несжимаемой жидкости. Основное отличие моделирования течений вязкой и идеальной жидкостей заключается в постановке граничных условий для функции скорости на твердой границе. Так как в задаче о деформации жидкого эллипса вся граница расчетной области является свободной, то данную задачу можно использовать для тестирования алгоритма движения по времени, как для идеальной, так и для вязкой жидкости, при значении коэффициента динамической вязкости = 0. Удовлетворительные расчеты были получены для 1200 и более частиц. Расчеты проводились с шагом по времени t = 1 103 до момента времени t = 1, 5 с.

На рис. 1 приводится сравнение результатов расчета свободной границы, полученное обобщенным методом естественных соседей с решением системы обыкновенных дифференциальных уравнений для задачи Л.В. Овсянникова методом РунгеКутта-Фельдберга. На рис. 1 деформация жидкого эллипса приводится для моментов времени t = 1, 0 с, t = 1, 5 с. В численных расчетах число узлов по области бралось равным 1200, на свободной границе – 138. Ввиду симметрии области для сравнения приводится только верхняя полуплоскость. Результаты решения задачи методом естественных соседей хорошо согласуются с результатами, полученными методом РунгеКутта: значение относительной погрешности не превысило 0,15 %.

Далее были проведены расчеты задачи о движении вязкой жидкости в плоской каверне. Пусть рассматривается двумерная область квадратной формы с длиной грани l = 1. Границы области являются твердыми стенками; нижняя и боковые грани 1 3 – неподвижные; верхняя граница 4 перемещается с постоянной скоростью.

Граничные условия для данной задачи задавались следующим образом:

u1(x, t ) = 0, u2 (x, t ) = 0 – на твердых неподвижных стенках 1, 2, 3, u1(x, t ) = 1, u2 (x, t ) = 0 – на подвижной стенке 4, p n = 0 – на всех границах 1 4, p = p0 – в среднем узле нижней границы 2. Задача решается при значениях числа Рейнольдса Re = 100, 400,1000 и 1800. При этом в качестве характерной скорости выбирается скорость подвижной границы.

Проводится сравнение полученных картин течения для различных чисел Рейнольдса с расчетами других авторов. В частности, на рис. 2 приводится сравнение полученных автором линий тока для числа Рейнольдса Re = 1000 (рис. 2, а) с результатами, представленными в одной из последних работ, посвященных численному реОвсянников Л.В. Общие уравнения и примеры // Задача о неустановившемся движении жидкости со свободной границей: cб. работ теорет. отдела / Акад. наук СССР, Сиб. отд-ние; Ин-т гидродинамики. – Новосибирск:

Наука, 1967. – С. 3–75.

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

Рис. 2. Линии тока Re = 1000 : а) результаты автора; б) результаты работы [11] Для исследования скорости сходимости обобщенного метода естественных соседей проводилось сравнение числа итераций (временных шагов), за которое задача о движении жидкости в квадратной каверне сходится к стационарному решению с заданной точностью в зависимости от числа расчетных узлов и значений числа Рейнольдса.

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

Наиболее значимой характеристикой течения жидкости является давление. Такие бессеточные методы, как метод сглаженных частиц, не позволяют получить удовлетворительную картину поля давления в расчетной области. Модификация метода SPH – ISPH и метод MPS в связи с использованием модели несжимаемой жидкости способны представить более качественное распределение поля давления, но и здесь не удалось найти надежных результатов. Метод естественных соседей при решении системы уравнений Навье-Стокса позволяет восстанавливать давление с высокой точностью.

В качестве тестовой рассматривается задача об определении давления в покоящейся жидкости в прямоугольной области, верхняя граница которой является свободной, а боковые и нижняя – твердыми стенками. Задача решается для различного числа узлов области в отсутствии внешних сил, но при наличии ненулевого начального распределения скоростей. На свободной границе задается условие p = 0, на твердых границах – p n = 0. Требуется получить распределение давления на твердых стенках, а также нулевое значение компонент вектора скорости внутри области. Расчеты проводились до момента времени t = 5 с, когда процесс устанавливался. К этому моменту численное давление было близко к гидростатическому. Относительная погрешность давления не превосходила 0, 01 %, а отклонение скорости от нулевого значения – 0,007 %.

Ковеня В.М. Об одном алгоритме решения уравнений Навье-Стокса вязкой несжимаемой жидкости // ИВТ СО РАН. – Новосибирск. – 2006. – Т. 11, № 2. – С. 39–51.

Для тестирования взаимодействия вязкой жидкости конечной глубины с твердыми стенками была выбрана задача о колебаниях жидкости в прямоугольном бассейне и определении нагрузок на стенки бассейна. В расчетной области D : x [0; ], y [0; 1 + 0,25 cos(x )] жидкость движется под действием силы тяжести. В начальный момент времени распределение поля скоростей u(x, 0) = 0. Моделирование колебаний движения жидкости проводилось для различных значений числа Рейнольдса: Re = 400, Re = 10000, Re = 50000. Выбранное в начальный момент времени возвышение жидкости не приводит к образованию нелинейных режимов движения, при которых происходит формирование и дальнейшее обрушение волновых структур. Профили свободной границы, а также значения гидродинамических нагрузок, создаваемых идеальной жидкостью на твердых стенках, сравнивались с результатами, полученными комплексным методом граничных элементов (КМГЭ) для потенциальной модели идеальной жидкости12.

Одним из важнейших преимуществ метода естественных соседей перед классическими сеточными методами является возможность численного моделирования течений, сопровождающихся большими деформациями расчетной области. Для получения моментов обрушения расчетная область модернизировалась следующим образом: D : x [0; ], y [0; 2 + 1,1cos(x )]. Решение задачи осуществлялось методами GNEM и КМГЭ. Комплексный метод граничных элементов позволяет проводить моделирование лишь до момента соприкосновения гребня волны с подошвой. Дальнейший расчет становится невозможен вследствие нарушения связности области.

Обобщенный метод естественных соседей позволяет проводить моделирование течений, сопровождающихся сильными деформациями расчетной области. На рис. 3 приведены хронограммы нагрузок, создаваемых вязкой и идеальной жидкостями на правой и левой вертикальных стенках: кривая 1 (сплошная) – расчет методом GNEM, кривая 2 (пунктирная) – КМГЭ. Значения нагрузок после момента обрушения получены методом GNEM.

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

Для проверки сходимости метода показано качественное совпадение картин течений для различного числа узлов области. Приведено сравнение текущих результатов с Афанасьев К.Е., Стуколов С.В. КМГЭ для решения плоских задач гидродинамики и его реализация на параллельных компьютерах: учеб. пособие – Кемерово, 2001. – 206 с.

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

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

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

Рис. 4. Сравнение результатов работы авторов с экспериментальными данными Также проводится анализ численных результатов распределения гидродинамических нагрузок на вертикальные стенки области при изменении высоты h слоя жидкости при основании. В результате точно определено значение h, при котором момент обрушения возникает после отката образующейся в результате разрушения плотины волны.

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

Для оценки возможной силы удара жидкости по горизонтальной поверхности была проведена серия расчетов для следующего диапазона значений высоты преграды hg = 4h, 6h, 8h, 10h, 12h, где h – высота слоя жидкости при основании.

Del Pin F. The meshless finite element method applied to a lagrangian particle formulation of fluid flows // Instituto de Desarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral noviembre, 2003. – 157 p.

Crespo J.S. Effect of wet bottom on dam break evolution // SPH European research interest community SIG. – 2007.– № 6. – 3 p.

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

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

Рис. 5. Взаимодействие волны с горизонтальной преградой:

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

В первом параграфе рассматриваются вопросы геометрической декомпозиции расчетной области и балансировки загрузки процессоров, способы распределения данных по процессорам, а также особенности параллельной реализации обобщенного метода естественных соседей. Параллельный подход к численному решению задач механики жидкости при помощи обобщенного метода естественных соседей основан на декомпозиции матрицы СЛАУ на блоки, количество которых равняется числу процессоров; расчете каждым процессором своего блока и обмене данными между ними на каждом шаге по времени15. В качестве языка реализации программного комплекса был выбран Fortran 90 с расширением библиотечными функциями MPI (Message Passing Interface).

Многопроцессорные системы. Построение, развитие, обучение / К.Е. Афанасьев, В.Г. Домрачев, И.В. Ретинская и др. // М.: КУДИЦ-ОБРАЗ, 2005. – 224 с.

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

В третьем параграфе для определения эффективности и ускорения реализованного параллельного алгоритма обобщенного метода естественных соседей была проведена серия расчетов на кластерах кафедры ЮНЕСКО по НИТ Кемеровского государственного университета и СКИФ Cyberia Томского государственного университета. В качестве тестовой решается задача об обрушении плотины.

Исследование производительности и эффективности программного кода выполняется в зависимости от размерности задачи и числа процессоров.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ, ПОЛУЧЕННЫЕ В ДИССЕРТАЦИИ

1. Предложен и реализован метод построения расширенной триангуляции Делоне на основе диаграммы Вороного.

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

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

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

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

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

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

ОСНОВНЫЕ ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

Журналы, рекомендованные ВАК для представления основных научных результатов диссертации:

1. Рейн, Т.С. Метод естественных соседей для решения задач вязкой несжимаемой жидкости [Текст] / К.Е. Афанасьев, Т.С. Рейн // Вестник НГУ (Серия «Математика, механика, информатика»). – 2008. – T. 8, № 2. С. 31–38.

2. Метод естественных соседей на основе интерполяции Сибсона [Текст] / К.Е.

Афанасьев, С.Н. Карабцев, Т.С. Рейн, С.В. Стуколов // Вестник ТГУ. Выпуск «Информационные технологии и математическое моделирование» (Серия «Математика. Кибернетика. Информатика»). – 2006. – № 19. – С. 210–219.

Труды конференций:

3. Рейн, Т. С. Применение метода естественных соседей к решению задач механики жидкости со свободными поверхностями [Текст] / Т. С. Рейн, С. Н. Карабцев // Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». – Кемерово: ИНТ, 2006. – С. 393–401.

4. Карабцев, С.Н. Реализация эффективного алгоритма построения диаграмм Вороного на плоскости [Текст] / С.Н. Карабцев, Т.С. Рейн, С.В. Стуколов // Труды V Всероссийской научно-практической конференции «Недра Кузбасса.

Инновации». – Кемерово: ИНТ, 2006. – С. 114–120.

5. Рейн, Т.С. Решение модельных задач гидродинамики методом естественных соседей [Текст] / Т.С. Рейн, С.Н. Карабцев // Труды VI Всероссийской научнопрактической конференции «Инновационные недра Кузбасса. IT-технологии».– Кемерово: ИНТ, 2007. – С. 311–317.

6. Афанасьев, К.Е. Сравнительное исследование алгоритмов интерполяции Сибсона и Лапласа [Текст] / К.Е. Афанасьев, Т.С. Рейн // труды VII Всероссийской научно-практической конференции «Инновационные недра Кузбасса. ITтехнологии». – Кемерово: ИНТ, 2008. – С. 286–291.

7. Рейн, Т.С. Стабилизация условия несжимаемости при решении задач о движении вязкой несжимаемой жидкости методом естественных соседей [Текст] / Т.С. Рейн // Программа и тезисы докладов VII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (с участием иностранных ученых). – Красноярск, 2006. – С. 27.

8. Рейн, Т.С. Параллельная реализация метода естественных соседей при решении задач механики жидкости со свободной поверхностью [Текст] / Т.С. Рейн, С.Н. Карабцев // Тезисы III Российско-Германской школы по параллельным вычислениям на высокопроизводительных вычислительных системах. – Новосибирск, 2006. – С. 15.

9. Афанасьев, К.Е. Численное моделирование задачи о разрушении плотины при наличии слоя жидкости при основании методом естественных соседей [Текст] / К.Е. Афанасьев, Т.С. Рейн // Материалы международной конференции «Сопряженные задачи механики реагирующих сред, информатики и экологии». – Томск, 2007. – С. 159–160.

10. Рейн, Т.С. Движение вертикальной стенки, закрепленной на пружинах, под действием волны, формирующейся при обрушении плотины [текст] / Т.С. Рейн // Программа и тезисы докладов VIII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям. – Новосибирск, 2007. – С. 69–70.

Подписано в печать 04.05.2008. Формат 60х841/16. Бумага офсетная № 1.

Печать офсетная. Усл. печ. л. 1,0. Тираж 100 экз. Заказ № 650043, г. Кемерово, ул. Ермака, 7. Тел. 58-34-