CC BY
91
УДК 532.5: 519.652
МЕТОД ЕСТЕСТВЕННЫХ СОСЕДЕЙ ДЛЯ РЕШЕНИЯ ЗАДАЧ ВЯЗКОЙ И ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ
К. Е. Афанасьев, Т. С. Рейн, С. Н. Карабцев (keafa@kemsu.ru)
Данная работа посвящена разработке численных алгоритмов на основе метода естественных соседей для решения задач динамики вязкой и идеальной несжимаемых жидкостей со свободными границами.
The present work is devote to development of numerical algorithms bases on natural element method for solving of hydrodynamics problems with free boundaries.
Ключевые слова: гидродинамические нагрузки, большие деформации, свободные границы, идеальная жидкость, вязкая жидкость, математическое моделирование, комплексы программ.
Введение
Моделирование течений несжимаемой жидкости со свободными границами, сопровождающихся нелинейными эффектами и большими деформациями расчетной области, невозможно без использования современных численных методов. Физические эксперименты для изучения подобных течений оказываются сложными и дорогостоящими, а быстрота протекания реальных процессов делает численные методы практически единственным источником информации о картине течения. Наибольшее распространение в задачах моделирования течений несжимаемой жидкости получили методы конечных разностей, конечных и граничных элементов, комплексный метод граничных элементов, а также метод контрольных объемов [6, 9]. Эти методы принадлежат классу сеточных и обладают одним общим недостатком. В процессе расчета сетка, на которой строится решение, сохраняет свою узловую связность, что, при больших деформациях расчетной области, может приводить к ее вырожденности.
С ростом производительности компьютеров развитие получили бессеточные методы, которые аппроксимируют уравнения в частных производных, основываясь только на наборе узлов, без знания дополнительной информации о структуре сетки. В таких методах отношение соседства частиц не фиксировано и может со временем изменяться, то есть частицы, бывшие соседями в начальный момент времени, могут со временем расходиться достаточно далеко друг от друга. Характерными представителями этой группы методов являются метод сглаженных частиц (SPH - Smoothed Particle Hydrodynamics) [16], полунеявный метод движущихся частиц (MPS -Moving Particle Semi-implicit) [15], метод Лагранже-во-Эйлеровых частиц [10]. Данные методы позволяют достаточно точно воспроизводить кинематику течений, однако получение динамических характеристик, необходимых для расчета гидродинамических нагрузок, является трудоемкой задачей. К общим недостаткам бессеточных методов можно отнести сравнительно невысокую точность и трудность введения граничных условий.
Эти обстоятельства заставили исследователей искать новые методы, сочетающие в себе идеи и возможности бессеточного подхода, но вместе с тем обладающие достоинствами сеточных методов. Первым из бессеточных методов нового поколения
появился бессеточный метод естественных соседей (Natural Element Method) [19]. Особенность метода NEM в том, что для стационарных задач он является обычным (классическими) методом Галеркина, то есть является сеточным. Для нестационарных задач применяется подход Лагранжа к описанию среды, именно: на каждом шаге по времени по найденному на предыдущем шаге положению узлов строится сетка, определяющая новую структуру соседей для каждой узловой точки области, на которой решается аппроксимированная система уравнений. В силу этого, метод NEM сохраняет некоторые преимущества классического метода Галеркина, а именно: простоту функций формы в области определения, непрерывность между элементами, легкость введения граничных условий. При этом имеет все достоинства бессеточных методов, так как функции формы метода естественных соседей зависят только от положения узловых точек.
В настоящей работе рассматривается численный алгоритм моделирования движения вязкой и идеальной жидкости со свободными границами методом естественных соседей. Приводятся результаты решения задачи о колебаниях жидкости в прямоугольном бассейне и определяются нагрузки на стенки бассейна. Представлено сравнение результатов решения для вязкой жидкости с результатами, полученными при решении задачи движения идеальной жидкости в полной постановке для уравнений Эйлера и комплексного потенциала.
§ 1. Постановка задачи
Выпишем полную нелинейную постановку задачи о движении идеальной несжимаемой жидкости со свободной границей.
В расчетной области течения D, представленной конечным набором узлов, ограниченной свободной поверхностью Г0 и твердыми границами Г, Г2 и Г3 (рис. 1), задано течение идеальной несжимаемой жидкости, описываемое системой уравнений Эйлера и уравнением неразрывности:
^ ^ + f,, ко в D, , = й.
Dt pdxt '
диі / dxj = 0 , x(t) є D .
(1)
(2)
| Вестник КемГУ № 2 2009 Математика
Рис. 1. Схема области течения и основные обозначения
Здесь:
х(г) = (х1 (г), х2 (г)) - пространственные координаты,
и(х, г) = (м1 (х, г), и2 (х, г)) - вектор скорости,
р(х, г) - давление,
р - плотность,
^ = (/1, /2) = (0, - g) - вектор массовых сил.
Движение расчетных узлов во всей области описывается уравнением вида:
ёх1 / & = и1, х(г) е Б, I = 1,2. (3)
На свободной поверхности Г0 выполняется динамическое условие р (х, г) = раШ , на твердых стенках Г1, Г2, Г3 выполняется условие непротекания и • п = 0, где п = (п1,п2) - внешняя нормаль к границе жидкости.
Ниже приводится общая постановка задачи о движении вязкой несжимаемой жидкости. Пусть в области течения Б происходит движение ньютоновской вязкой несжимаемой жидкости, описываемое системой уравнений Навье-Стокса и уравнением неразрывности (2). В постановке Эйлера уравнения Навье-Стокса будут иметь следующий вид:
^ = _Цр+Л—(^и.) + /, , = 173.1 = ¡3. (4)
Бг р дх1 р дх. дх.
В системе уравнений (2), (4) искомыми функциями являются давление р и вектор скорости и , параметрами - плотность р , вектор внешних сил f и коэффициент динамической вязкости /л, по повторяющемуся индексу 1 производится суммирование.
Запишем граничные условия для системы уравнений Навье-Стокса (2),(4). Так как жидкость вязкая, то на твердых стенках Г1, Г2 и Г3 выполняется условие прилипания: и1 = 0, и = 1,2 . На свободной поверхности Г0, как и случае идеальной жидкости, выполняется динамическое условие р = рат .
Для нестационарной задачи о движении, как вязкой, так и идеальной жидкости, необходимо задать положение расчетных узлов х(0) = х0 и распределе-
ние неизвестных функций во всей области течения: и(х,0) = и0 (х), р(х, 0) = р0 (х) .
§ 2. Алгоритм решения
1. Дискретизация по времени
Интегрирование по времени систем уравнений Эйлера и Навье-Стокса, описывающих поведение идеальной и вязкой жидкостей соответственно, представляет некоторые трудности, когда жидкость несжимаемая или слабосжимаемая. В таком случае не может быть использован явный метод интегрирования по времени, в частности, не удается устранить нефизические осцилляции функции давления. Для систем дифференциальных уравнений в частных производных в работе [12] был предложен в общем виде метод расщепления (метод дробных шагов). Суть этого метода заключается в разбиении физического процесса на два: конвекцию-диффузию и вклад давления. На первом этапе в уравнении движения учитываются только конвективные члены, в результате чего выделяется фиктивная переменная и* - предиктор скорости. Второй этап заключается в добавлении к и* члена (-AtVp), который будет обеспечивать соленоидальность и на шаге по времени t + At. Для расщепления уравнения движения Навье-Стокса используется выражение:
£)„ и"и _ + ++ _ м * + и* _ _"
Бг Аг Аг
1 д л д дии
= (-
р дх д р дх д дх.
/ /
+ 1/2
(5)
р дх . р дх . дх .
= [ф ( х:, ? ) ]
г+1/2 = 1(фи1 + ф"),
где Аг = гп+1 - ги - шаг по времени; под символом фи понимается выражение фп =ф(хп, г"), соответст-
___ п /,п п\ п+1 / ,п+1 п+1 \
венно, и * = и1 (г , х ), и1 = и (г , х ).
Разделив слагаемые в уравнении (5) в соответствии с идеей метода расщепления, получим следующие выражения для предиктора и корректора скорости:
и д Г)и*+1/2
и* = и+ /иг +------------------(—*------)Аг,
р дх д дх д
лг д
(6)
+ + Л ^ (7)
р дх1
При расщеплении уравнения движения Эйлера, вследствие отсутствия в нем слагаемого, отвечающего за вязкость, уравнение (6) перепишется в виде: и* = и+ + А/ . (8)
Для интегрирования уравнения неразрывности вводится фиктивная функция плотности р, удовлетворяющая уравне нию Б р Бг + рди1/ дх1 = 0 [13], тогда
бр р -р _р -р +р -р Бг ~
А г
А г
р"+1 - р*_ + р -р" = _р 5(и"+1 - и,. + и,. )
А г
А/
Йг.
5(и. - и.) ди.
= -р----------------------р--------
дх, дх.
(9)
где р , так же как и при расщеплении уравнения движения, является фиктивной переменной. Из выражения (9) получим:
р - р _ 8и{
Аг
дх,.
р_+1-р___рд( -А*
Аг
(10)
(11)
Используя (7) и (11), можно записать уравнение Пуассона для функции давления:
^—.и + 1 _* +2 2_
Р + = +_. (12)
Аг дх1
Подстановка функции р* из уравнения (10) в (12) приводит последнее к виду:
р ди* д2рп_
Аг2
+ -
А? дх.
дхд
(13)
Описанный выше метод расщепления относится к классу проекционных методов [4].
Самый простой способ ввести условие несжимаемости - записать равенство: рп+1 = рп = р° = р. Тогда первое слагаемое в левой части уравнения (13) обратится в ноль. Тем не менее при лагранже-вом подходе к описанию среды лучше оценить это слагаемое, чтобы избежать накопления численной ошибки на каждом временном шаге. Перепишем условие несжимаемости, полагая, что на временном шаге гп+1 плотность равна начальной, то есть рп+1 = р0 = р . Однако, из-за вычислительной ошибки, это равенство не обязательно будет выполняться, поэтому значение плотности рп должно оцениваться на каждом временном шаге с учетом упомянутой выше ошибки. Таким образом, уравнение (13) может быть переписано в виде:
р ~р Аг2
р дии д дрп
А г дхи дх д дх д
(14)
Алгоритм движения по времени состоит из следующих шагов.
I) задание на п -ом временном слое распределения узлов области в виде массива;
II) построение аппроксимирующих функций формы;
III) вычисление предиктора скорости и из системы (6), в случае моделирования движения вязкой жидкости, и из системы (8) - в случае идеальной;
IV) решение уравнения Пуассона (14) для опре-
п+1
деления давления р ;
V) вычисление нового значения скорости ип+1 из уравнения (7) с учетом найденного на шаге IV давления;
VI) вычисление нового положения частиц на (п +1)-ом временном шаге: хп+1 = хп + ип+1Аг и далее на шаг I.
2. Основные этапы метода естественных соседей
Метод естественных соседей был предложен Л. Траверсони в 1994 году для решения задач теории пластичности [9]. Этот метод представляет собой разновидность метода Галеркина. Как и в методе Галеркина, неизвестные функции аппроксимируются следующим образом: q = ЫТQ .
Для формирования дискретной системы уравнений используется метод взвешенных невязок в интегральной форме. При этом интегралы берутся по элементам расширенной триангуляции Делоне [3].
Можно выделить следующие основные шаги метода естественных соседей:
- расчетная область разбивается ячейками Вороного первого порядка [7];
- для текущего набора узлов определяется граница расчетной области [15];
- для каждого узла, используя дискретизацию Вороного, определяется множество его естественных соседей;
- генерируются элементы расширенной триангуляции Делоне;
- с помощью метода Галеркина записывается слабая форма уравнений движения и неразрывности. Для аппроксимации неизвестных функций вычисляются функции формы Сибсона и Лапласа и их производные. Для вычисления интегралов используются квадратуры Хаммера;
- полученная система линейных алгебраических уравнений решается методом сопряженных градиентов с предобусловливанием.
3. Интерполяционные функции Сибсона и Лапласа
Интерполяционные функции Сибсона, основанные на понятии естественных соседей [7], базируются на диаграммах Вороного первого и второго порядков и определяются через отношение площадей многоугольников [18]:
а, (х) = А1 (х)/А(х), А(х) = ^А (х), I = 1, Ь . (15)
3=1
Здесь Ь - число естественных соседей для точки х = (х1, х2); А(х) - площадь ячейки Вороного первого порядка, в которую входит точка х, А1 (х) - площадь пересечения ячейки Вороного второго порядка точки х и площадью А(х) (рис. 2а).
*
| Вестник КемГУ № 2 2009 Математика
Рис. 2. К построению интерполяции: а) Сибсона; б) Лапласа
Построение интерполяции Лапласа также опирается на определение соседей посредством разбиения области ячейками Вороного [13]. Пусть точка х принадлежит многоугольнику Вороного с числом сторон равным N. Обозначим длины сторон много,, I = 1, N , а высоты, опущенные через И1. Тогда интерполяционные коэффициенты Лапласа примут следующий вид:
С N л-1
угольника через из точки х на *,
X 3/ , I = 1, N .
Уз=1 У
Такой способ определения коэффициентов
(16)
а
биение области ячейками Вороного, то многие из многоугольников разбиения будут треугольными элементами, на которых функции формы Лапласа ведут себя как линейные функции. Тот факт, что точное поведение функций формы на некоторых многоугольниках известно заранее, может упростить задачу интегрирования, а также вычисления производных.
Производные функций формы Сибсона и Лапласа могут быть получены дифференцированием уравнений (15) и (16) соответственно.
На основе функций формы Лапласа можно построить еще один вид интерполяции, основанный на понятии ячеек Вороного: расширенную интерполяцию Лапласа [1]. Ее отличие от классической заключается в выборе множества соседних узлов для точки х , а именно: интерполяция проводится по едям х, удовлетворяющим критерию опи-жружности Делоне [7]. Носитель функции формы в этом случае будет совпадать с носителем интерполяционной функции Сибсона, а коэффициенты интерполяции будут представлять функции формы Лапласа.
Точность использования интерполяционных функций Сибсона и Лапласа проверялась на решении уравнения Пуассона в ограниченной области Б при заданном значении правой части по известной функции (метод пробных функций). Будем рассматривать область Б , верхняя граница Г1 которой описывается уравнением у = 0.5 • &'п(х), х е [0,2п], а боковые и нижняя границы Г 2, Г3, Г 4 являются прямыми линиями: Г2 : х = 0, у е [-1,0]; Г3 :у = -1, х е[0,2п]; Г4 :х = 2п, у е[-1,0].
В описанной выше области решалось уравнение Пуассона Аи = Ь для функции и = х • у2 с граничными условиями Дирихле на Г1 и Неймана на границах Г 2, Г3, Г 4:
и = х • у2, (х, у) еГ1; ди/дп = - у2, (х, у) еГ2; ди/дп = -2•х• у, (х,у) еГ3; ди/дп = у2, (х,у) е Г4.
В таблице 1 приводится относительная погрешность значений функции и для различного числа узлов области N.
проще и экономичнее, чем в подходе Сибсона, так как не требует вычисления площадей пересечения многоугольников. На рисунке 2б изображена схема построения интерполяции Лапласа для двумерного случая.
Несибсоновская интерполяция имеет особенность, которая выводится из ее основных свойств. Если для заданного множества точек построить раз-
Относительная погрешность функций формы (Ф.Ф.) Сибсона, Лапласа и расширенной интерполяции (Р. И.) Лапласа (%)]
Таблица 1
N 1209 3122 5552 8668 10500
Ф.Ф. Сибсона 1.1387 Ч10-3 8.0157 Ч10-4 6.4594 Ч10-4 5.1295 Ч10-4 3.7801 Ч10-4
Ф.Ф. Лапласа 1.0495 Ч10-3 5.1352 Ч10-4 6.9359 Ч10-4 7.1723 Ч10-4 7.7108 Ч10-4
Ф.Ф. Р.И.Лапласа 1.1387 Ч10-3 6.0531 Ч10-4 5.0360 Ч10-4 5.2752 Ч10-4 4.6223 Ч10-4
Из данных таблицы видно, что с увеличением числа узлов алгоритм решения задачи, использующий для вычисления неизвестных функции формы Сибсона и функции формы расширенной интерполяции Лапласа, дает более точные результаты, чем алгоритм, в основе которого лежат функции формы Лапласа. Однако при одинаковом числе узлов области несибсоновская интерполяция требует меньше времени, чем интерполяция Сибсона и расширенная интерполяция Лапласа.
4. Стабилизация условия несжимаемости в системе уравнений Навье-Стокса
Одной из главных трудностей численного моделирования нестационарных уравнений Навье-Стокса является регуляризация условия несжимаемости (2). В описанном выше методе расщепления по пространственным переменным условие несжимаемости представлено уравнением Пуассона для давления (14). Для устранения нефизических осцилляций функции давления используется метод конечных приращений (Finite Incrément Calculus [17]), основанный на понятии дифференциального приближения [11].
Устойчивость решения системы Навье-Стокса методами, основанными на методе Галеркина, обеспечивается выбором конечно-элементных пространств для скорости и давления: степени интерполяционных полиномов компонент вектора скорости и давления должны удовлетворять условию Лады-женской-Бабушки-Бреззи (ЛББ). В данной работе для аппроксимации уравнения неразрывности использовались линейные базисные функции (функции формы расширенной интерполяции Лапласа), для уравнения движения - квадратичные базисные функции (функции формы Сибсона). Построение такого обобщенного метода естественных соседей (GNEM - General Natural Element Method) приводит к удовлетворению условий ЛББ для совместной аппроксимации, что гарантирует невырожденность решения.
Движение идеальной несжимаемой жидкости моделируется методом естественных соседей (NEM), где в качестве интерполяционных функций для скорости и давления используются функции формы Сибсона.
§ 3. Тестовые расчеты
Тестирование метода естественных соседей и обобщенного метода естественных соседей проводилось на решении следующих задач: задача Л. В. Овсянникова о деформации жидкого эллипса -
со свободной границей [5] и задача о распределении давления в покоящейся жидкости конечной глубины с твердыми границами.
1. Задача Л. В. Овсянникова о деформации жидкого эллипса
В начальный момент времени область расчета представляет собой круг единичного радиуса, в который заключена несжимаемая жидкость. Деформация круга в эллипс начинается под действием начального распределения скоростей в отсутствии внешних сил.
Распределение поля скоростей следующее: и1 = а’(г)/а(г)-х, и2 = -&({)/а^)-у. С ростом г круг деформируется в эллипс
х21а2 (г) + а2 (г) у2 = 1, большая полуось которого а (г) ^ х, а малая 1/ а (г) ^ 0 при г ^ х. Коэффициент а (г) находится из решения задачи:
(а (г) +1)а (г) а (г) = 2 (а (г))2, (17)
а(0) = 1, а'(0) = 1.
Обыкновенное дифференциальное уравнение
(17) второго порядка решалось методом Рунге-Кутта-Фельдберга. Полученные значения а (г) на
различных временных шагах использовались для сравнения с численными результатами, найденными методом естественных соседей.
Общая постановка задачи, как правило, приводится для идеальной несжимаемой жидкости. Основное отличие моделирования течений вязкой и идеальной жидкостей заключается в постановке граничных условий для функции скорости на жесткой границе. Так как в задаче о деформации жидкого эллипса вся граница расчетной области является свободной, то данную задачу можно использовать для тестирования алгоритма движения по времени как для идеальной жидкости, так и для вязкой, при значении коэффициента динамической вязкости
Л = 0 .
Начальное поле скоростей с учетом (17) задавалось в виде: и1 (0) = х, и2(0) = - у. Давление на границе круга постоянно: р = 0 . Расчеты проводились
с шагом по времени Аг = 1 • 10-3 для различного числа узлов области. В таблице 2 для 1350 узлов приведены значения относительного отклонения длины главной полуоси.
Таблица 2
Относительная погрешность главной полуоси эллипса Д(а (t)), %
t, [с] 0.1 0.4 0.6 1.0 1.3 1.65
Идеальная жидкость 0.0258 0.0462 0.0805 0.1005 0.1284 0.1519
Вязкая жидкость 0.0274 0.0483 0.0840 0.1107 0.1380 0.1596
| Вестник КемГУ № 2 2009 Математика
Таблица 3
Относительная погрешность значений функции давления Д(р) %
N 300 600 900 1200 1500 1800
Вязкая жидкость 0.028 0.017 0.0094 0.0051 0.0024 0.00103
Идеальная жидкость 0.026 0.014 0.0087 0.0050 0.0022 0.00094
Из представленных данных видно, что результаты решения задачи методом естественных соседей (идеальная жидкость) и обобщенным методом естественных соседей (вязкая жидкость) хорошо согласуются с результатами метода Рунге-Кутта.
2. Задача об определении давления в покоящейся жидкости
Поиск давления в методах NEM и GNEM является неотъемлемой составной частью алгоритма движения по времени, в то время как во многих численных методах, в частности, в методе КМГЭ, поиск давления является дополнительной краевой задачей на основе известного распределения поля скоростей на каждом временном шаге. В качестве тестовой рассматривалась задача об определении поля давления в покоящейся жидкости в прямоугольной области Q , верхняя граница которой являлась свободной, а боковые и нижняя - твердыми стенками.
Задача решалась для различного числа узлов расчетной области в отсутствии внешних сил и при наличии нулевого начального распределения скоростей. На свободной границе задавалось условие p = 0 , на твердых границах - dp / dn = 0 . Расчеты проводились до момента времени t = 5 с. Во все время расчета численное давление было близко к гидростатическому. В таблице 3 представлена относительная погрешность функции давления для различного числа узлов области N для случаев вязкой и идеальной жидкостей.
§ 4. Численные расчеты
Применимость методов NEM и GNEM для моделирования течений идеальной и вязкой несжимаемой жидкостей демонстрируется на решении задачи о колебаниях жидкости в прямоугольном бассейне и определении нагрузок на стенки бассейна. Задача формулируется следующим образом. В расчетной области D : x е [0;^] , y е [0;1 + 0.25cos(x)] под действием силы тяжести происходит движение жидкости (рис. 1). В начальный момент времени распределение поля скоростей u(x,0) = 0. Выбранное возвышение жидкости в начальный момент времени не приводит к образованию нелинейных режимов движения, при которых происходит формирование и дальнейшее обрушение волновых структур. Вследствие этого, профили свободной границы, а также значения гидродинамических нагрузок [19], создаваемых идеальной жидкостью на твердых стенках, сравниваются с результатами, полученными комплексным методом граничных эле-
ментов для потенциальной модели идеальной жидкости [8].
Для расчетов методом NEM использовалось 3220 узлов области, из них на жесткой границе - 209, на свободной - 105. Число узлов на свободной границе области в методе КМГЭ также было выбрано 105. Задача является неустойчивой, если используется постоянный шаг по времени. В данной работе шаг по времени выбирался из следующего условия [9]:
Д'=min {кЬ[}- <18)
где s = min |xJ - x(| - минимальное расстояние
i * J
между узлами расчетной области.
3-t= 1 81
4-t“ 2.42
5-1-3.75
Рис. 3. Положение свободной границы идеальной жидкости в различные моменты времени
На рисунке 3 приведены профили свободной границы идеальной жидкости для метода естественных соседей (сплошная линия) и комплексного метода граничных элементов (пунктирная линия).
Из приведенного сравнения видно, что профили свободной границы идеальной жидкости, полученные методом естественных соседей, хорошо согласуются с результатами, полученными комплексным методом граничных элементов.
Также для различных значений числа Рейнольдса: Яе = 400, Яе = 10000, Яе = 50000 проводилось моделирование колебаний вязкой несжимаемой жидкости обобщенным методом естественных соседей. Хорошее совпадение с приведенными на рисунке 3 профилями свободной границы наблюдается при Яе = 50000 .
| Вестник КемГУ № 2 2009 Математика
Рис. 4. Кривые динамических нагрузок
На рисунке 4 представлено сравнение хронограмм динамических нагрузок на левой вертикальной стенке области, рассчитанных КМГЭ (идеальная жидкость, потенциальная постановка, кривая 1), NEM (идеальная жидкость, нелинейная постановка, кривая 2), GNEM (вязкая жидкость, кривые 3-5 соответствующие значения числа Рейнольдса -Re = 400; 10000; 50000). Из представленных результатов видно, что хронограммы нагрузок, создаваемые идеальной жидкостью, полученные методами NEM и КМГЭ, совпадают с графической точностью. Хронограммы динамической нагрузки вязкой жидкости на твердых стенках бассейна лежат ниже кривых 1 и 2, отражающих значения нагрузок для идеальной жидкости. Для системы уравнений Навье-Стокса на жестких границах было задано условие прилипания, что означает отсутствие скольжения вдоль границы (для идеальной жидкости на жестких границах выполняется условие непротекания). В связи с этим в вязкой жидкости скорость приграничных частиц жидкости и соответственно давление на твердые границы ниже, чем в идеальной. Однако с увеличением числа Рейнольдса до 50000 значение нагрузок вязкой жидкости приближается к кривой 2, рассчитанной методом естественных соседей
Одним из важнейших преимуществ метода естественных соседей перед классическими сеточными методами является возможность численного моделирования течений, сопровождающихся большими
деформациями расчетной области. Для получения моментов обрушения изменим расчетную область следующим образом:
Б : х є [0;^],у є [-2; 1.1cos(x)].
Решение задачи осуществлялось методами ЫБМ, GNEM и КМГЭ для числа узлов области 3690, из которых 145 узлов принадлежат свободной границе, 254 - твердой. Шаг по времени также выбирался из условия 18. Комплексный метод граничных элементов позволяет проводить моделирование лишь до момента обрушения. Дальнейший расчет становится невозможен вследствие нарушения связности области. Метод естественных соседей позволяет проводить моделирование как самого процесса обрушения, так и движения жидкости после него. На рисунке 5 приведено сравнение профилей свободной границы, полученной различными методами. Кривая 1 соответствует методу ЫБМ, 2 - 0№М, 3 (пунктирная линия) - КМГЭ. Результаты сравнения показывают достаточно точное совпадение профилей границ. Наибольшее расхождение соответствует моменту времени непосредственно перед обрушением гребня волны.
Момент соприкосновения гребня волны с подошвой происходит для разных методов в разные моменты времени: для ЫБМ - г = 4,56 , для GNEM -г = 4,5, для КМГЭ - г = 4,61. Моделирование течения после момента обрушения возможно методами ЫБМ и GNEM. На рисунке 6 приведены картины течения после моментов обрушения. Закраска области выполнена в соответствии с полем давления.
На рисунке 7 приведены хронограммы динамических нагрузок, создаваемых вязкой и идеальной жидкостями на правой и левой вертикальных стенках: кривая 1 (сплошная) - расчет методом ЫБМ, кривая 2 (штрихпунктирная) - GNEM, кривая 3 (пунктирная) - КМГЭ. После момента обрушения гребня волны приведены значения нагрузок, полученных методами ЫБМ и GNEM.
| Вестник КемГУ № 2 2009 Математика
3 2 ] 0 ГЕ3.43| 3'-
1 2 1
з 3 X
3 2 ] 0 11-4.53! 2—
■
■
■
) з 1 X
Рис. 5. Сравнение профилей свободной границы: кривая 1 - НЕМ, 2 - ОМЕМ, 3 - КМГЭ
11-6.27 ■ ■
% Р 0 0.4 0.8 1.2 1 I
¡¡1 .ш
2 &1Ш§Йь.
1
0
1 3 X
Рис. 6. Картина течения после момента обрушения: а) НЕМ, б) ОМЕМ
а)
"1
1 / Г 'к л \ /2
г \ И /
5 / \ ]> \ /
2 > V ч
- / / # х — з \\ А ч\
1
0 V 1 1
6 у 8
б)
4 Л
Л ]
Л
2 !' 1 м
V | 1/ 1/ к
0 \
\ /> У/
\ V у 1
1
-4
1 - . 1
0 2 4 Т К
Рис. 7. Динамическая нагрузка: а) правая стенка, б) левая стенка
|_^естникКемГУ
Для контроля консервативности представленных методов №М и GNEM проводилась проверка закона сохранения полной энергии системы. К концу расчета относительное отклонение энергии для метода ЫБМ составило 0,7 %, а для GNEM - 0,93 %.
Заключение
В работе предложена модификация метода естественных соседей, позволяющая проводить численное моделирование задач динамики вязкой и идеальной несжимаемой жидкостей со свободными границами на всех этапах вычислительного эксперимента, включая этапы сильно нелинейного поведения течения с образованием многосвязности расчетной области. Характерной особенностью представленных методов ЫБМ и GNEM является возможность расчета давления жидкости в области и на границе, которое в дальнейшем может быть использовано для получения значений динамических нагрузок, создаваемых жидкостью на твердых стенках области.
Развитие условно-бессеточных численных методов - метода естественных соседей и обобщенного метода естественных соседей - является важной и актуальной задачей для моделирования течений идеальной и вязкой несжимаемой жидкостей, что дало возможность получать решение новых задач со свободными границами, решение которых традиционными методами или весьма затруднительно, или вообще невозможно.
Литература
1. Афанасьев, К. Е. Сравнительное исследование алгоритмов интерполяции Сибсона и Лапласа / К. Е. Афанасьев, Т. С. Рейн // Труды VII Всероссийской научно-практической конференции «Инновационные недра Кузбасса. 1Т-технологии ». - Кемерово: ИНТ, 2008. - С. 286 - 291.
2. Беликов, В. В. Несибсоновская интерполяция - новый метод интерполяции значений функции на произвольной системе точек / В. В. Беликов,
B. Д. Иванов, В. К. Конторович, С. А. Корытник, А. Ю. Семенов // Вычислительная математика и математическая физика. - 1997. - Т. 37. - № 1. -
C. 11 - 17.
3. Карабцев, С. Н. Эффективный алгоритм генерации конечноэлементной сетки для метода естественных соседей / С. Н. Карабцев, С. В. Стуколов //Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». - Кемерово: ИНТ, 2006. -С. 401 - 409.
4. Марчук, Г. И. Введение в проекционносеточные методы / Г. И. Марчук, В. И. Агощков. -М.: Наука, 1981. - 416 с.
5. Овсянников, Л. В. Общие уравнения и примеры. Задача о неустановившемся движении жидкости со свободной границей / Л. В. Овсянников. -Новосибирск: Наука, 1967.
6. Патанкар, С. Численные методы решения задач теплообмена и динамики жидкости / С. Н. Па-танкар. - М.: Энергоатомиздат, 1984. - 152 с.
7. Скворцов, А. В. Триангуляция Делоне и ее применение / А. В. Скворцов. - Томск: ТГУ, 2002. -128 с.
8. Стуколов, С. В. Решение нелинейных волновых задач гидродинамики идеальной жидкости комплексным методом граничных элементов: авто-реф. дис. ... канд. физ.-мат. наук. - Кемерово, 1999.
- 26 с.
9. Терентьев, А. Г. Численные методы в гидродинамике: учебное пособие / А. Г. Терентьев, К. Е. Афанасьев. - Чебоксары: Чуваш. ун-т., 1987. -80 с.
10. Франк, А. М. Дискретные модели несжимаемой жидкости / А. М. Франк. - М.: Физматлит, 2001. - 208 с.
11. Шокин, Ю. И. Метод дифференциального приближения / Ю. И. Шокин. - Новосибирск: Наука, 1979. - 224 с.
12. Яненко, Н. Н. Об экономичных неявных схемах (метод дробных шагов) / Н. Н. Яненко // Докл. АН СССР. - 1960. - Т. 134. - 5 с.
13. Chorin, A. Numerical solution of the Navier-Stokes equations / A. Ctoim // Math. Comp. - 1968. -Vol. 22. - P. 745 - 762.
14. Facundo, P. The meshless finite element method applied to a lagrangian particle formulation of fluid flows: In Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy. Instituto deDesarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral / P. Facundo. - 2003. -157 p.
15. Koshizuka, S. A particle method for incompressible viscous flow withfluid fragmentation / S. Koshizuka, H. Tamako, Y. Oka // Computational Fluid Dynamics Journal. - 1995. - P. 29 - 46.
16. Monaghan, J. Smoothed particle hydrodynamics / J. Monaghan // Ann. Rev. Astronand Astrophysics.
- 1992. - № 30. - P. 543 - 574.
17. Onate, E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation / E. Onate // Comput. Meth. Appl. Mech. Eng. - 2000. - Vol. 182. - № 1 - 2. -Р. 355 - 370.
18. Sibson, R. A brief description a natural neighbor interpolation / R. Sibson, V. Barnett (ed.) // Interpret multivariate data. - Chichester: John Wiley, 1981. -P. 21 - 36.
19. Sukumar, N. The natural element method in solid mechanics / N. Sukumar, B. Moran, T. Belyt-schko // Int. J. Num. Methods Eng. - 1998. - Vol. 43. -№ 5. - P. 839 - 887.
20. Traveroni, L. Natural neighbor finite elements / L. Traveroni // In International Conference on Hydraulic Enginnering Software. Hydrosoft Proceedings
2.Computational Mechanics Publications. - 1994. -Р. 291 - 297.
Рецензент - В. П. Потапов - д-р техн. наук, профессор, директор института угля и углехимии СО РАН.
