Применение метода естественных соседей к решению задач механики жидкости со свободной поверхностью

,

Кемеровский государственный университет

Кемерово, Россия

E-mail: *****@***ru, *****@***ru

АННОТАЦИЯ The finite element method (FEM), which was conceived in the 1950s, has now become the most widely-used numerical method for computer simulation in academic research as well as in engineering practice. Finite elements provide the possibility of handling arbitrary geometries with complex boundary conditions. A practical drawback, however, is the need for regeneration of the mesh in moving boundary and large deformation problems. This is often done by the analyst, and is considered to be one of the most time-consuming tasks in a finite element analysis.

To overcome the difficulty associated with remeshing, the past decade has seen a tremendous surge in the development of a family of Galerkin and collocation-based numerical methods—these are known as particle, gridless, meshfree, or meshless methods. Meshless methods share the common characteristic of there being no explicit connectivity information between nodes; the approximant is built in a process that is transparent to the user.

The natural element method (NEM) is a Galerkin-based meshless method that is built upon the notion of natural neighbour interpolation. This interpolation scheme has very striking properties, such as its strictly interpolating character, ability to exactly interpolate piece-wise linear boundary conditions, and a well-defined and robust approximation with no user-defined parameter on non-uniform grids.

НЕ нашли? Не то? Что вы ищете?

In this paper, we review the most notable aspects of the NEM with emphasis on the recent advances achieved by the authors in its application to fluid mechanics.

ВВЕДЕНИЕ

В настоящее время в связи с возрастающей сложностью исследуемых явлений наблюдается пересмотр подходов к численному моделированию задач механики жидкости со свободными границами. Известные и широко используемые методы подвергаются модификации, методы, использовавшиеся ранее для исследования явлений в других областях знаний, находят свое применение при исследовании задач механики жидкости. Использование новых подходов позволяет расширить класс решаемых задач и получить новые результаты. К одному из наиболее сложных для моделирования классов задач гидродинамики относятся задачи со свободными границами, сопровождающиеся сильно-нелинейной деформацией границ течений. Ряд такого рода задач (в двумерной, осесимметричной и пространственной постановке) был решен методом конечных элементов и комплексным методом граничных элементов (МКЭ и КМГЭ). Однако данные методы позволяют проводить расчет подобных задач только до определенной стадии, например, обрушение волны можно рассчитать лишь до момента соприкосновения гребня волны с ее подошвой, далее проводить расчет становится невозможным в силу изменения связности области расчета. Суть методов, подобных методам конечных и граничных элементов, состоит в следующем: расчетная область разбивается на более мелкие подобласти (элементы), на каждом из элементов неизвестные функции аппроксимируются гладкими, например, линейными функциями. Основным недостатком МКЭ и КМГЭ является потребность в перестройке сетки в задачах с подвижными границами. На четырехугольных элементах при решении задач гидродинамики с большими деформациями происходит перехлест сетки, что приводит к остановке расчета. Корректное разбиение пространства на треугольные элементы может производиться различными способами и является неединственным. При этом точность интерполяции зависит от качества триангуляции: показано, что ошибка при линейной интерполяции уменьшается при увеличении минимального угла [7]. Таким образом, для того, чтобы погрешность была наименьшей, необходимо на этапе триангуляции области максимизировать острый угол в треугольнике, при этом треугольники должны быть такими, чтобы предел отношения площадей треугольников с максимальной и минимальной площадью был равен единице. Среди множества триангуляций можно выделить триангуляцию Делоне. В триангуляции Делоне достигается максимум минимального угла по всем треугольникам, хотя при этом она остается неоднозначной (достаточно рассмотреть четырехугольник, для триангуляции которого может быть проведена как одна, так и другая диагональ). Такая неоднозначность приводит к структуре соседей, которая может резко меняться при малом изменении координат узлов и, тем самым, приводить к резким изменениям в результате интерполяции и к ее неоднозначности.

Чтобы преодолеть трудности, связанные с перестройкой сетки, в последнее время наблюдается стремительное развитие в разработке новых подходов на основе методов Галеркина и методов коллокаций. Они известны как методы частиц или бессеточные методы (БМ). Общее в этих методах то, что они не имеют никакой явной информации о связи между узлами. Но, таким образом, в стандартных бессеточных методах для того, чтобы найти интерполяционные функции, требуется обеспечить узловую связность. Так стали разрабатываться методы, сочетающие в себе сеточный и бессеточный подходы и использующие преимущества каждой из методологий. Характерным представителем этой группы методов является метод естественных соседей (Natural Element Method). В данной работе рассматривается численный алгоритм решения нестационарных задач гидродинамики со свободной поверхностью методом естественных соседей.

ПОСТАНОВКА ЗАДАЧИ

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

Рис.1. К постановке нестационарной задачи

Пусть в расчетной области течения , ограниченной свободной поверхностью и твердыми стенками (рисунок 1) решается уравнение Эйлера:

. (1)

В области выполняется закон сохранения массы:

. (2)

На твердых стенках выполняется условие непротекания:

, (3)

где – внешняя нормаль к поверхности жидкости.

На свободной поверхности выполняется условие:

. (4)

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

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

АЛГОРИТМ ДВИЖЕНИЯ ПО ВРЕМЕНИ

Интегрирование по времени уравнений (1) и (2) представляет некоторые трудности, когда жидкость несжимаемая или слабо сжимаемая. В этом случае, не может быть использован явный метод интегрирования по времени, и даже при использовании неявного временного шага, не удается устранить нефизических колебаний функции давления. Для преодоления этих трудностей используется метод расщепления, предложенный в работе [7], который состоит в разбиении каждого временного шага на два:

, (5)

, (6)

где - шаг по времени; , , а и являются фиктивными переменными, необходимыми для метода расщепления.

Из уравнений (1) и (2) следует

, (7)

(8)

Используя уравнение (8) можно записать:

, (9)

что с учетом уравнения неразрывности приводится к виду:

, (10)

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

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

I) Задание границы области в виде массива точек, построение диаграммы Вороного;

II) Построение аппроксимирующих функций формы в ячейках Вороного;

III) Нахождение скорости из системы (8):

;

IV) Решение уравнения Пуассона (10) для определения давления;

V) Вычисления нового значения скорости из уравнения (7):

;

VI) Вычисление нового местоположения частиц:

. (11)

ПРОСРАНСТВЕННАЯ ДИСКРЕТИЗАЦИЯ

Лагранжева схема расщепления, описанная в предыдущем разделе имеет ряд преимуществ:

¾  Шаг I является линейным и явным. Он позволяет определить положение различных точек без пространственной дискретизации скоростей;

¾  Рассмотренная схема расщепления устраняет стандартные колебания, появляющиеся при нахождении функции давления;

¾  Из всех шагов схемы единственный неявный шаг – решение уравнения Пуассона на давление (шаг IV). Его решение может быть получено, используя схему пространственной дискретизации классического метода конечных элементов.

Рассмотрим схему решения уравнения (10) с граничными условиями (3)-(4). Как и в классическом МКЭ, неизвестные функции аппроксимируются следующим образом:

, (12)

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

Используя метод взвешенных невязок можно записать слабую форму уравнения (10):

, (13)

где - заданное на жесткой границе значение скорости.

Интегрируя по частям уравнение (13) и учитывая условие непротекания (3), получим:

, (14)

где на - главные граничные условия;

на - естественные граничные условия.

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

Метод естественных соседей был предложен Л. Траверсони в 1994 году для решения задач теории пластичности [9]. Этот метод представляет собой разновидность бессеточного метода Галеркина. Для формирования дискретной системы уравнений используется метод взвешенных невязок в интегральной форме. При этом интегралы берутся по треугольным элементам, которые образуют с текущим узлом пары его естественных соседей. Множество естественных соседей для каждого узла, а также узлы свободной границы на новом временном шаге определяются с помощью методов sweep line и - shape, описанных в работе и «Эффективный алгоритм генерации конечноэлементной сетки для метода естественных соседей», представленной в данном сборнике [2].

Сама идея NEM заключается в использовании интерполяций Сибсона и Лапласа, базирующихся на диаграмме Вороного, основными характеристиками которой являются единственность разбиения и выпуклость всех ее элементов. В настоящей работе рассматриваются Сибсоновские функции формы. Такая интерполяция использует меры, соответствующие размерности пространства, что дает более гладкие результаты, чем интерполяция Лапласа, рассматривающая меры, на единицу меньше размерности соответствующего пространства, каждая из которых, в свою очередь, дополнительно делится на линейный размер [6]. Так можно выделить следующие основные шаги метода естественных соседей:

1.  Для каждого узла определяется множество его естественных соседей;

2.  Генерируются треугольные элементы, используя текущий узел и узлы найденного в п.1 множества;

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

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

. (15)

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

Рис. 2. Разбиение области ячейками Вороного

Рис. 3. Интерполяция Сибсона

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

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

. (16)

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

, (17)

- есть площадь ячейки Вороного первого порядка вокруг точки и - площадь пересечения ячейки Вороного первого порядка точки с площадью (рисунок 3). Если точка , то и все другие функции формы равны нулю.

Производные функций формы Сибсона получаются дифференцированием уравнения (17):

. (18)

Строгое определение интерполяции Сибсона можно сформулировать в следующем виде. Пусть - отдельные точки евклидова пространства . Пусть и , где - евклидово расстояние, - ячейка Вороного для . Тогда . Если , то , где - мера Лебега в . Таким образом определенные называются коэффициентами Сибсона. При этом вводится как (запись для точки ). В работах [3,5] доказываются свойства интерполяционных функций: функции формы Сибсона положительно определены, интерполируют узловые данные и линейно полны.

Для вычисления функций формы Сибсона необходимо искать пересечение ячеек Вороного первого и второго порядков, что является очень трудоемкой задачей вычислительной геометрии. Б. Уотсоном был предложен алгоритм поиска площадей пересечения многоугольников путем разбиения многоугольника ориентированными треугольниками (в каждом из треугольников задается порядок обхода вершин) [8]. Следующим псевдокодом опишем алгоритм вычисления Сибсоновской интерполяции:

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

Определяем множество естественных соседей () и центры треугольников () для каждой точки ;

Для {

Для {Вычисляем для каждого треугольника центры описанных окружностей , и их производные по координатам };

Для {

o  определяем - номер глобального узла в диаграмме разбиения, соответствующего локальному узлу

o  определяем - позицию узла с номером в массиве натуральных соседей точки ;

o  , ;

o 

o 

} };

Для {

}.

ТЕСТИРОВАНИЕ

В работе , «КМГЭ для решения плоских задач гидродинамики и его реализация на параллельных компьютерах» описывается численное решение уравнения Лапласа комплексным методом граничных элементов [1].

Для тестирования метода естественных соседей проводилось сравнение численных результатов, полученных при решении методами NEM и КМГЭ уравнения Лапласа в области :

, (19)

с граничными условиями

Дирихле: на участке границы ;

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

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

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

Рис. 4. Геометрия области

Рис. 5. Относительная погрешность NEM и КМГЭ

В представленной на рисунке 4 области решалось уравнение Лапласа (19) для гармонической в области функции со следующими граничными условиями:

Дирихле: , ;

Неймана: ;

.

Сравнение рассматриваемого в работе метода естественных соседей с КМГЭ показало, что NEM имеет более высокую точность решения, чем комплексный метод граничных элементов при решении уравнения Лапласа с заданными граничными условиями (рисунок 5).

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

Распределение поля скоростей следующее:

С ростом круг деформируется в эллипс , большая полуось которого , а малая при .

Коэффициент находится из решения задачи

(20)

Обыкновенное дифференциальное уравнение (20) второго порядка было решено методом Рунге-Кутта-Фельберга (RKF) [5]. Полученные значения на различных временных шагах использовались для построения свободной границы жидкого эллипса и сравнения с численными расчетами.

При решении задачи методом NEM, в виду симметрии, рассматривалась только верхняя полуплоскость. Частицы заданы по всему объему жидкости.

Начальное поле скоростей с учетом (20):

.

Давление на границе круга постоянно: .

Расчеты проводились с шагом по времени до момента времени с. Относительная погрешность отклонения не превысила 0.1%.

На рисунке 6 представлено распределение частиц в моменты времени t=0 с., t=0.3 с., t=1.0 с. и t=2.2с. соответственно.

Рис. 6. Распределение частиц в моменты времени t=0 с., t=0.3 с., t=1.0 с. и t=1.5с. соответственно

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

ЗАКЛЮЧЕНИЕ

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

ЛИТЕРАТУРА

1. , КМГЭ для решения плоских задач гидродинамики и его реализация на параллельных компьютерах: Учеб. Пособие / Кемерово 2001.

2. , Стуколов алгоритм генерации конечноэлементной сетки для метода естественных соседей. // Материалы III международной научной летней школы «Гидродинамика больших скоростей и численное моделирование», Кемерово 2006.

3. Овсянников уравнения и примеры

// Задача о неустановившемся движении жидкости со свободной границей. Новосибирск: Наука, 1967.

4. , Афанасьев методы в гидродинамике: Учеб. пособие / Чуваш. ун-т. им. . Чебоксары: ЧГУ, 1987.

5. Форсайт Дж., ашинные методы математических вычислений. – М.,Мир, 1980.

6. Belikov, V. V., Ivanov V. D., Kontorovich V. K., Korytnik S. A. and Semenov A. Yu (1997). The Non-Sibsonian Interpolation: A New Method of Interpolation of the Values of a Function on an Arbitrary Set of Points, Computational Mathematics and Mathematical Physics 37(1): 9-15.

7. Facundo Del Pin The Meshless Finite Element Method Applied to a Lagrangian Particle Formulation of Fluid Flows. Instituto de Desarrollo Tecnol_ogico para la Industria Qu__mica (INTEC) Universidad Nacional del Litoral Noviembre 2003

kumar N, Moran B, Belytschko T. The natural element method in solid mechanics. Int J Num Methods Eng 1998;43(5):839-887.

9. Traversoni, L. (1994). Natural neighbor finite elements, In International Conference on Hydraulic Enginnering Software, Hydrosoft Proceedings 2:  291-297, Computational Mechanics Publications.