Моделирование задач гидродинамики вязкой несжимаемой жидкости со свободными границами бессеточным методом естественных соседей Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 13, № 4, 2008

Моделирование задач гидродинамики вязкой несжимаемой жидкости со свободными границами бессеточным методом

естественных соседей

К. Е. Афанасьев, Т. С. Рейн Кемеровский государственный университет, Россия e-mail: [email protected], [email protected]

This work addresses problems of continuous mechanics for a medium with free boundaries using the Natural Element Method. A brief description of the method is given. A scheme for temporal evolution is considered. Results of numerical calculations of test problems and a numerical solution for the problem of a dam collapse in the case of a fluid layer on its base are considered.

Введение

Для моделирования неустановившихся физических явлений, описываемых различными уравнениями в частных производных, создан целый класс сеточных методов, для которых на каждом временном шаге требуется знание информации о межузловой связности. Наиболее популярные представители этой группы — метод конечных элементов (МКЭ) [1] и метод контрольных объемов (МКО) [2]. Основным недостатком МКЭ и МКО является то, что для каждого временного шага сетка, на которой строится решение, не теряет своей узловой связности, что, в свою очередь, при больших деформациях жидкости может быстро приводить к вырожденности сетки. Для корректного разбиения области течения на треугольные элементы могут использоваться различные способы. При этом точность интерполяции зависит от качества триангуляции; например, в работе [3] показано, что для МКЭ якобиан преобразования треугольных и четырехугольных элементов зависит от минимального угла элемента. Дополнительно, неединственность триангуляции приводит к структуре соседей, которая может резко меняться при малом изменении координат узлов и тем самым приводить к изменениям результата интерполяции.

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

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

бессеточные [4]. К группе таких методов, хорошо зарекомендовавших себя в задачах гидродинамики, относятся, например, метод сглаженных частиц (Smoothed Particle Hydrodynamics) [5], полунеявный метод движущихся частиц (Moving Particle Semi-implicit) [6], метод лагранжево-эйлеровых частиц [7]. К недостаткам этих бессеточных методов можно отнести сравнительно невысокую точность и трудность введения граничных условий.

Эти обстоятельства заставили исследователей искать новые методы, сочетающие в себе идеи и возможности бессеточного подхода, но вместе с тем обладающие достоинствами сеточных методов. Первым из бессеточных методов нового поколения появился бессеточный метод естественных соседей (Natural Element Method) [8, 9]. Особенность метода NEM в том, что для стационарных задач он — обычный (классический) метод Галеркина, т. е. сеточный. Для нестационарных задач применяется подход Лагранжа к описанию среды. А именно, на каждом шаге по времени по найденному на предыдущем шаге положению узлов строится сетка, определяющая новую структуру соседей для каждой узловой точки области; на вновь построенной сетке снова методом Галеркина решается система уравнений. В силу этого метод NEM сохраняет некоторые преимущества классического метода Галеркина, а именно: простоту функций формы в области определения, непрерывность между элементами, легкость введения граничных условий. При этом метод NEM имеет все достоинства бессеточных методов, так как функции формы метода естественных соседей зависят только от положения узловых точек.

Метод естественных соседей был предложен Л. Траверсони в 1994 году для решения задач теории пластичности [10]. Этот метод представляет собой разновидность метода Галеркина. Для формирования дискретной системы уравнений используется метод взвешенных невязок в интегральной форме. При этом интегралы берутся по треугольным элементам, которые образуют с текущим узлом пары его естественных соседей. Множество естественных соседей для каждого узла [11], а также узлы свободной границы на новом временном шаге определяются с помощью методов sweep line и a-shape, описанных в работах [12, 13]. Для аппроксимации неизвестных функций используются функции формы Сибсона и Лапласа [14, 15]. Полученная система линейных алгебраических уравнений решается методом сопряженных градиентов с предобусловливанием.

Авторами работ [16, 17] NEM был адаптирован к решению задач механики идеальной жидкости. В настоящей работе рассматривается модификация метода естественных соседей для течений вязкой несжимаемой жидкости. При построении интерполяционных функций используется алгоритм Уотсона [13]. Решение многомерной задачи на дробных шагах сводится к решению отдельных уравнений методом сопряженных градиентов. Приведены результаты решения задачи об обрушении плотины и определении гидродинамических нагрузок на стенки бассейна. Сравнение результатов расчетов с экспериментальными данными позволяет сделать вывод об эффективности рассматриваемого метода.

1. Постановка задачи и алгоритм решения

Пусть в некоторой области течения D происходит движение ньютоновской вязкой несжимаемой жидкости, описываемое системой уравнений Навье—Стокса [18]. В форму-

лировке Эйлера такая система будет иметь вид (при этом для описания среды применяется подход Лагранжа):

Дм* д д дш , ,

рж = — +^ах; дх;+^ (1)

_ 1(2)

Здесь х(£) = (хх, х2(£)), 2 =1, 2, ^ = 1, 2, по принятому правилу по повторяющемуся

индексу проводится суммирование. В системе уравнений (1)-(2) искомыми функциями являются давление р(х, £) и вектор скорости и(х, £) = (мх(х,£),м2(х,£)). Параметрами будут плотность р, коэффициент динамической вязкости а и вектор массовых сил / = (/1,/2).

Запишем граничные условия для системы уравнений Навье—Стокса (1)-(2). Так как жидкость вязкая, то на твердых стенках Ги выполняется условие прилипания:

м = м*, г = 1, 2. (3)

На свободной границе Гст выполняется динамическое условие:

Р = Ра^. (4)

Для нестационарной задачи необходимо задать положение расчетных узлов х(°) = х0 и распределение неизвестных функций в области течения: и(х, °) = и0(х), р(х, °) = р0(х).

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

Для систем дифференциальных уравнений в частных производных в работе [19] был предложен в общем виде метод расщепления (метод дробных шагов) для построения экономичных неявных схем. Затем в докладе [20] этот метод был сформулирован для интегро-дифференциальных уравнений и многослойных схем в дробных шагах. Для интегрирования системы уравнений Навье—Стокса метод расщепления позднее был описан в работе [21]. Суть этого метода заключается в разбиении описываемого физического процесса на два этапа: конвекцию—диффузию и вклад давления. На первом этапе в уравнении движения учитываются только конвективные члены, в результате чего выделяется фиктивная переменная и* — предиктор скорости. Второй этап заключается в добавлении к и* члена (—Д£Ур), который будет обеспечивать соленоидальность и на шаге по времени £ + Д£.

Для расщепления уравнения движения Навье—Стокса используется выражение

Дм* мп(1 — мп м™+1 — и* + и* — МП ( 1 д а д дм* „ . . ,

~ —-^ = "-"-*-- = ' о Р ТТ" + /* , (5)

Д£ Д£ \ р дх* р дх; дх;

(—Рдк р+Рд; И+/Г' 2=[ф(х.£)г+х/24 (Ф-+х+Ф-),

где Д£ = £П(1 — £П — шаг по времени; под символом фП понимается выражение фП ф(жП,£П), соответственно мП = мДх^^).

Разделив слагаемые в уравнении (5) в соответствии с идеей метода расщепления,

получим следующие выражения для предиктора и корректора скорости:

a) =<+f-At+^^At; (6)

c) <+" = AApn« (7)

Для интегрирования уравнения неразрывности вводится фиктивная функция плотности р, удовлетворяющая уравнению Dp/Dt + p(dui/öxi) = 0 [21], тогда

Dn nn+l pn pn+1 n* + n* pn pn+1 n* n* pn

DP p — p _ p — p + p — p _ p — p p — p _ - ^ - — - — - + - —

Dt At At At At

= — d«+1 — u* + u*) = — d«+1 — О — d(u*), (8)

где p*, как и при расщеплении уравнения движения, является фиктивной переменной. Из уравнения (8) получим

^ = — (9)

Pn+1 — Р* = d(un+1 — u*) (10)

At = — P dXi • (10)

Используя уравнения (7) и (10), можно записать следующее уравнение Пуассона для функции давления:

Ь) ^=4(и)

Подстановка функции р* из уравнения (9) в выражение (11) приводит последнее к

виду +1 d d д

Ь) At2 +ÄtdXi = dXj dXXjp + • (12)

Самый простой способ ввести условие несжимаемости — записать равенство

pn+1 = pn = р0 = р• (13)

Тогда первое слагаемое в левой части уравнения (12) обратится в ноль. Тем не менее при использовании лагранжева подхода лучше оценить это слагаемое, чтобы избежать накопления численной ошибки на каждом временном шаге. Перепишем условие несжимаемости, полагая, что на временном шаге tn+1 плотность равна начальной, т. е.

Pn+1 = Р0 = Р• (14)

Из-за вычислительной ошибки плотность pn не обязательно будет равна начальной

0n

плотности р0, поэтому значение плотности р' должно оцениваться на каждом временном шаге с учетом упомянутой выше ошибки. Способ оценки pn будет представлен ниже. Таким образом, уравнение (12) может быть переписано в виде

р0 — pn + = д dpn+1 ) At2 + AtdXi дх; дх; • ( )

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

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

2) построение аппроксимирующих функций формы;

3) вычисление предиктора скорости и* из системы (6);

4) вычисление нового значения плотности рп;

5) решение уравнения Пуассона (15) для определения давления рп+1;

6) вычисления нового значения скорости ип+1 из уравнения (7) с учетом найденного на шаге 5 давления;

7) вычисление нового местоположения частиц: хп+1 = хп + ип+1Д£. Пространственная дискретизация. Как и в классическом методе Галеркина,

неизвестные функции аппроксимируются следующим образом:

и = ^ N (х,г)иа

I

р ЩХ,1)Рг,

(16)

где N — глобальные базисные интерполяционные функции, определенные в каждый

момент времени на множестве расчетных узлов; Цц, Рг, ^ — узловые значения функций Щ,р,р.

Интерполяционные функции метода естественных соседей (Сибсона и Лапласа) обладают следующими свойствами:

— линейно полны: х = ^ Nг(x) хг;

представляют разложение единицы ^ N1 = 1

г

удовлетворяют условию Кронекера: Nl(xj) =

1, I = 3, 0, I = 3-

Рассмотрим схему решения системы уравнений (6), (7), (15) с граничными условиями (3) и (4). Используя метод взвешенных невязок, можно записать слабую форму уравнений (6), (7), (15), при этом к граничным условиям (3) и (4) также применяется описанная выше схема расщепления:

а) I (и* - ип) - /Д1 - = 0,

v

р дxj дxj

(17)

ь) /Ч^ + ^ = ]тдх3 дх3

v v

дд Nг(рп+1^У,

с) I N1 {«+1 - и*)Д}^ = I N1 дХ{-рп+1}вУ.

V

V

Интегрируя уравнения (17) и (18) по частям, получим:

(18) (19)

/Г £) ЛТ £) *+ 2

ЩК - < - +

v v

/ди ■ 2

N1^— Щ Г = 0, (20)

Ь) N(р0 - рп)вУ + [ ^ Ш =

v

дг2}

v

дху дху

v

рп+1дУ + у N

Ги

дхг д

к + дхг (-Рп+1)

п^ Г.

(21)

С учетом (19) граничный интеграл правой части уравнения (21) можно переписать в виде

д

N

к + дх, (-р"+1)

п4Г = у N [и"+1] щ dГ.

Ги

(22)

Главные и естественные граничные условия для системы уравнений (20)-(21) и (19) будут выглядеть следующим образом:

Го- : р = РаЬш)

Ги : и"+1 = 0.

(23)

(24)

С учетом (23), (24) выразим подынтегральное выражение граничного интеграла уравнения (20):

0.

Гди*+1/2 1 1 \дп* ] 1 Гдм"+1 1

Г щ Гс = 2 п Щ дх7 . + 2 Г2 Г с я щ

Таким образом, граничные интегралы в уравнениях (20), (21) обращаются в ноль. Используя аппроксимацию (16), дискретную форму уравнения (20) запишем:

а) I N N dVЩ = I NN dVЩ+

v v

+дг/^ ¡^ - И дг/"^ ^ dvu*+2

v

р ,} дх3 дх3 7

V

Выражение (25) в матричной форме запишется следующим образом:

ми* = ми" + дг¥ - и*+1.

р

Аналогично дискретная форма уравнения Пуассона (21) будет иметь вид

1

Ь)

дг2

NNр^0dv - N Nрjndv -

\У

V

дг I дх, j ij

т дч dvpn+l.

дхг дхг 7

v v

Выражение (27) в матричной форме запишется в виде

-м -§Г) + и* = БР"+1.

дг2 дг

(25)

(26)

(27)

Г

с

Также для уравнения (19) дискретная форма будет выглядеть следующим образом:

дг Г , ЯМ

с)

N м3 ауип+1

г]

NNаущ - - ] N ауРт+1-

(29)

V V V

Выражение (29) в матричной форме будет выглядеть следующим образом:

м ип+1 = М и* - —Вт Рп+1.

Р

Матрицы в выражениях (26), (28) и (30) имеют следующий вид:

м = у NN ау, в = у ау, в = у VN ау,

v v v

^ = ^ NтМУ,

v

в 0 0

К

0 в 0 0 0 в

(30)

(31)

(32)

(33)

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

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

J Р0ау = J рпйУ,

v (4=4°) v(t=tn)

что с учетом уравнения (16) дает:

(34)

(№)тау^0 = ¿у^1

тп\тлгч&п (35)

v (4=4°) v(t=tn)

Объем, связанный с каждой узловой точкой, может быть найден следующим образом

(^п)1

ау,

(36)

вирр(Пг)

что с учетом уравнения (35) преобразуется к виду

(П0)т Ш0 = (Пп)т

Равенство (37) представляет собой закон сохранения массы. Вектор Пп можно рассматривать как вектор, содержащий объемы, связанные с каждой частицей. В данной работе Пп вычисляются с помощью диаграммы Вороного. При этом закон сохранения массы, выраженный уравнением (37), понимается следующим образом: каждая частица сохраняет свою собственную массу. Таким образом, искомый член уравнения (28) может быть найден из выражения

Пп(р° - рп) = Ппр° - П°р° = р0(Пп - (38)

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

Тестовые расчеты. Для отладки метода МЕМ решались три задачи: Пуассона при заданном значении правой части по известной функции (метод пробных функций); задача Овсянникова о деформации жидкого эллипса [22]; задача об определении давления в покоящейся жидкости конечной глубины.

Решение уравнения Пуассона. Пусть в ограниченной области О с границей Г = Г + Г^ решалось уравнение Пуассона

Д^ = Ь, (39)

со следующими граничными условиями:

<р(х) = <£(х) для (х) € Г^; д(х) = д(х) для (х) € Гд.

Будем рассматривать область О, верхняя граница Г^ которой задается уравнением у = 0.5вт(х), х € [0, 2п], а граница Гд = Г1 + Г2 + Г3, — Г^ х = 0, у € [— 1,0]; Г2: у = -1, х € [0, 2п]; Г3: х = 2п, у € [-1, 0].

В описанной выше области решалось уравнение Пуассона (39) для функции V = ху2 с граничными условиями Дирихле на Г^ и Неймана на Г1, Г2, Г3:

V = xУ2, (х,у) €

дV 2 / Ч ^ -л О < ^ ^ -п дV 2

дп -у , (х,у) € Г1; дП = -2xУ, (х,у) € Г2; дП = y, (х,у) € Г3-

В табл. 1 приводится относительная погрешность значений функции V для различного числа узлов области N.

Из данных таблицы видно, что полученные результаты показывают высокую точность МЕМ для решения уравнения Пуассона.

зЗадача о деформации жидкого эллипса. В начальный момент времени область расчета представляет собой круг единичного радиуса, в который заключена несжимаемая жидкость. Деформация круга в эллипс начинается под действием начального распределения скоростей в отсутствие внешних сил. Распределение поля скоростей следующее: и1(£) = а'(£)/а(£)х, и2(£) = -а'(£)/а(£)у. С ростом £ круг деформируется в эллипс

Таблица 1. Относительная погрешность

N 100 1000 1500 2000 2500 3000

Д(<р), % 0.0421 0.0359 0.0176 0.0081 0.0037 0.00192

x2/a2(t) + a2(t)y2 = 1, большая полуось которого a(t) ^ то, а малая — 1/a(t) ^ 0 при t ^ то. Коэффициент a(t) находится из решения задачи:

(a4(t) + 1)a (t) a" (t) = 2 (a'(t))2 , a (0) = 1, a' (0) = 1.

(1)

Обыкновенное дифференциальное уравнение (1) второго порядка решалось методом Рунге—Кутты—Фельдберга. Полученные значения а (£) на различных временных шагах использовались для сравнения с численными результатами, найденными методом естественных соседей.

Общая постановка задачи, как правило, приводится для идеальной несжимаемой жидкости. Основное отличие моделирования течений вязкой и идеальной жидкостей заключается в постановке граничных условий для функции скорости на жесткой границе. Так как в задаче о деформации жидкого эллипса вся граница расчетной области является свободной, то данную задачу можно использовать для тестирования алгоритма движения по времени как для идеальной, так и для вязкой жидкости, при значении коэффициента динамической вязкости ^ = 0. Начальное поле скоростей с учетом (1) задавалось в виде: «1(0) = х, и2(0) = —у. Давление на границе круга постоянно: р = 0. Расчеты проводились с шагом по времени Д£ = 110-3 для различного числа узлов области.

На рис. 1 представлено распределение узлов расчетной области, полученное методами МЕМ и ИКР в моменты времени £ = 0 с, £ = 0.3, £ =1.0 и £ = 1.5 с соответственно для числа узлов 1200. Ввиду симметрии области при сравнении рассматривается только верхняя полуплоскость.

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

В табл. 2 для 1350 узлов приведены значения относительного отклонения длины главной полуоси.

Распределение поля давления в покоящейся жидкости. Поиск давления в методе естественных соседей — неотъемлемая составная часть алгоритма движения по времени, в то время как во многих численных методах, в частности, в методе КМГЭ, поиск давления — дополнительная краевая задача на основе известного распределения поля

Рис. 1. Распределение частиц в моменты времени: а — t = 0 c, б — t = 0.3 c, в — t = 1.0 c, г t = 1.5 с

Таблица 2. Относительная погрешность полуосей эллипса, %

с Д(а (*)) Д(1/а (*))

0.3 0.041454 0.0414374

0.7 0.091195 0.090494

1.0 0.121245 0.118134

1.5 0.144734 0.138825

Таблица 3. Относительная погрешность поля давления

п 300 600 900 1200 1500 1800

Д(р), % 0.028 0.017 0.0094 0.0051 0.0024 0.00103

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

Задача решалась для различного числа узлов расчетной области в отсутствие внешних сил и при наличии нулевого начального распределения поля скоростей и давления. На свободной границе задавалось условие р = 0, на твердых границах — др/дп = 0. Расчеты проводились до момента времени Ь = 5 с. Во все время расчета численное давление было близко к гидростатическому. В табл. 3 представлена относительная погрешность функции давления для различного числа узлов области п.

2. Решение задачи о разрушении плотины при наличии слоя жидкости на основании

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

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

0.2

0.1

0.2

0.4

0.6

0.8

Рис. 2. Схема расчетной области

0

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

Задача является неустойчивой, если используется постоянный шаг по времени. В работе [17] приводится сравнение расчетов задачи о распаде столба жидкости методом естественных соседей, полученных при использовании постоянного и переменного временных шагов. В данной задаче шаг по времени выбирался из следующего условия [8]:

At = min < / ., I , (1)

x \|u(x,t)|J ' V J

где s = min |xj — — минимальное расстояние между узлами расчетной области.

Таблица 4. Параметры задачи

Параметр Значение Размерность

Скорость перегородки плотины 1.5 м/с

Плотность жидкости 1000 кг/м3

Коэффициент динамической вязкости 3.2 • 10-3 кг/(м • с)

Сила тяжести 9.8 м/с2

Высота слоя жидкости при основании Н 1.8 • 10-2 м

Высота слоя жидкости слева от плотины 0.15 м

О 0.2 0.4 0.6 0.8 1 X 0 0.2 0.4 0.6 0.8 1 X

Рис. 3. Разрушение плотины. 4030 узлов (переменный шаг по времени): а — £ = 0.156 с; б — £ = 0.243 с; в — £ = 0.281 с; г — £ = 0.406 с; д — £ = 0.468 с; е — £ = 0.531 с

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

Для сравнения на рис. 4 и 5 показаны фрагменты течений для 2300 и 5800 частиц соответственно для различных моментов времени. Шаг по времени также выбирался из условия (41). Из приведенных рисунков видно качественное совпадение полученных результатов.

На рис. 6 приводится сопоставление расчетов обрушения плотины с экспериментом. Данные эксперимента были взяты из работы [23].

Вычисление энергетических характеристик. Проверка закона сохранения полной энергии Е = Ек + Ер служит для оценки консервативности метода МЕМ. Кинетическая и потенциальная энергии на каждом шаге по времени вычисляются по следующим формулам [24]:

2

ти2

Ек = — ,ЕР = тду. (42)

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

n 2

Е = £( ^ + т,дуг), (43)

г=1

где N — общее число узлов области.

Рис. 4. Разрушение плотины. 2300 узлов (переменный шаг по времени): а — Ь = 0.156 с; б — Ь = 0.243 с; в — Ь = 0.281 с; г — Ь = 0.406 с; д — Ь = 0.468 с; е — Ь = 0.531 с

0 0.2 0.4 0.6 0.8 1 X 0 0.2 0.4 0.6 0.8 1 X

0 0.2 0.4 0.6 0.8 1 X 0 0.2 0.4 0.6 0.8 1 X

Рис. 5. Разрушение плотины. 5800 узлов (переменный шаг по времени): а — Ь = 0.156 с; б — Ь = 0.243 с; в — Ь = 0.281 с; г — Ь = 0.406 с; д — Ь = 0.468 с; е — Ь = 0.531 с

Таблица 5. Отклонение полной энергии системы, %

Ы/Ь, с 0 0.156 0.219 0.281 0.343 0.406 0.468 0.531

2300 0.00 0.07 0.41 0.58 0.93 1.21 1.45 1.70

4030 0.00 0.05 0.29 0.36 0.44 0.52 0.60 0.73

5800 0.00 0.03 0.24 0.31 0.40 0.45 0.49 0.56

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

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

Способ определения нагрузки подробно описан в работе [25], согласно ему нагрузка Р* на твердой границе Г вычисляется по следующей формуле:

Р = / Рш аг, (44)

г

где Рад = Р — Р0 — волновое давление, Р — распределение давления на каждом шаге по времени, Р0 — распределение давления в начальный момент времени г = 0.

Y

0.1

о

Y 0.1

0.4

0.4

0.4

0.4

Y 0.1

0.4

0.6

0.6

0.6

0.6

0.6

0.8

0.8

0.8

0.8

0.8

0.4

0.6

0.8

1 X

1 X

1 X

1 X

1 X

1 X

Рис. 6. Сравнение результатов работы авторов с экспериментальными данными в моменты времени: а — г = 0.156 с; б — г = 0.243 с; в — г = 0.281 с; г — г = 0.406 с; д — г = 0.468 с; е — г = 0.531 с; слева — экспериментальные данные, справа — текущие расчеты

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

На рис. 7 приведены хронограммы нагрузок на вертикальные границы области для значений h = 0.018 м, h = 0.038, h = 0.09 и h = 0.12 м соответственно на правой (рис. 7, а) и левой (рис. 7, б) вертикальных стенках области. При h = 0.12 м толщина слоя практически соответствует высоте столба жидкости и выбрана эмпирически на основе проведенной серии расчетов. Такая высота слоя при основании соответствует картине течения, при которой для заданной длины бассейна не нарушается связность

а

0

в

г

д

0

е

Рис. 7. Хронограммы гидродинамической нагрузки на вертикальных стенках области: а — правая стенка области; б — левая стенка области; 1 — h = 0.018 м, 2 — h = 0.038, 3 — h = 0.09, 4 — h = 0.12 м

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

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

Заключение

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

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

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

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

[1] Коннор Дж., Бреббиа К. Метод конечных элементов в механике жидкости. Л.: Судостроение, 1979. 204 с.

[2] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

[3] Терентьев А.Г., Афанасьев К.Е. Численные методы в гидродинамике: учеб. пособие. Чебоксары: Чуваш. ун-т., 1987. 80 с.

[4] Li S., Liu W.K. Meshfree and particle methods and their applications // Appl. Mech. Rev. 2002. N 55. P. 1-34.

[5] MoNAGHAN J. Smoothed particle hydrodynamics // Ann. Rev. Astron and Astrophysics. 1992. N 30. P. 543-574.

[6] Koshizuka S., Tamako H., Oka Y. A particle method for incompressible viscous flow with fluid fragmentation // Compu. Fluid Dynamics J. 1995. Vol. 4. P. 29-46.

[7] Франк А. М. Дискретные модели несжимаемой жидкости. М.: Физматлит, 2001. 208 с.

[8] 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 de Desarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral noviembre. 2003. 157 p.

[9] Onate E., Idelsohn S.R., Zienkiewicz O.C.A finite point method in computational mechanics. Applications to convective transport and fluid flow // Intern. J. for Num. Methods in Eng. 1996. Vol. 39. P. 3839-3866.

[10] Traveroni L. Natural neighbor finite elements //In Intern. Conf. on Hydraulic Enginnering Software. Hydrosoft Proceedings 2. Comput. Mechanics Publ. 1994. Р. 291-297.

[11] Скворцов А.В. Триангуляция Делоне и ее применение. Томск: ТГУ, 2002. 128 с.

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

Sukumar N., Moran B., Belytschko T. The natural element method in solid mechanics // Intern. J. Num. Methods Eng. 1998. Vol. 43, N 5. P. 839-887.

[14] Беликов В.В., Иванов В.Д., Конторович В.К. и др. Несибсоновская интерполяция — новый метод интерполяции значений функции на произвольной системе точек // Вычисл. математика и мат. физика. 199T. Т. 3T. № 1. С. 11-1T.

[1Б] SlBSON R. A vector identity for the Dirichlet Tessellation jj Math. Proc. Cambridge Phil. Soc. 19BG. Vol. BT, N 1. P. 1Б1-1ББ.

[16] Карабцев C.H., Рейн Т.е. Применение метода естественных соседей к решению задач механики жидкости со свободной поверхностью j j Матер. III Междунар. научной летней школы "Гидродинамика больших скоростей и численное моделирование". Кемерово: ИНТ, 2GG6. С. 393-4G1.

[1T] Афанасьев К.Е., Карабцев C.H., Рейн Т.С., Стуколов C.B. Метод естественных соседей на основе интерполяции Сибсона j j Вестник Томского гос. ун-та. 2GG6. № 19. С. 21G-219.

[1B] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука. 19BT. B4G с.

[19] Яненко H.H. Об экономичных неявных схемах (метод дробных шагов) jj Докл. АН СССР. 196G. Т. 134, № Б. С. 1G34-1G36.

[2G] Марчук Г.И., Яненко H.H. Применение метода расщепления (дробных шагов) для решения задач математической физики jj Докл. на Всесоюзн. конф. по вычисл. матем. (Москва, февраль 196Б). Сб.: Некоторые вопросы прикл. и вычисл. матем. Новосибирск, 1966.

[21] Chorin A. Numerical solution of the Navier—Stokes equations jj Math. Comput. 196B. Vol. 22. P. 745-762.

[22] Овсянников Л.В. Общие уравнения и примеры. Задача о неустановившемся движении жидкости со свободной границей. Новосибирск: Наука, 196T.

[23] Alejandro Crespo J.S. Effect of wet bottom on dam break evolution jj SPH European research interest community SIG. 2GGT. Vol. 6. P. 3.

[24] Фабер Т.Е. Гидроаэродинамика. М.: Постмаркет, 2GG1. ББ9 с.

[2Б] Афанасьев К.Е., Березин E.H. Анализ динамических характеристик при взаимодействии уединенной волны с препятствием j j Вычисл. технологии. 2GG4. Т. 9, № 3. С. 22-3T.

Поступила в редакцию 15 июня 2QQ7 г., в переработанном виде — 25 марта 2QQ8 г.