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

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

Том 3, № 1, 1998

МОДЕЛИРОВАНИЕ СИЛЬНО НЕЛИНЕЙНЫХ ВОЛНОВЫХ ТЕЧЕНИЙ

К. Е. Афанасьев Кемеровский государственный ниверситет, ИВЦ, Россия e-mail: keafa@kemgu.kemerovo.su

The paper considers some results of investigations of gravitation waves behavior in the homogeneous non-condensable liquid caused by the influence of obstacles moving along the bottom.

В настоящей работе проведены начатые в [1, 2] исследования движения полуцилиндра из состояния покоя по ровному дну. Как и раньше, изучаются кинематические характеристики возникающего течения, однако основное внимание уделено исследованию характеристик волны, возникающей за телом, и движения жидкости в том диапазоне чисел Фруда, при которых не существует стационарных течений. Первая группа задач условно называется “существенно нелинейными задачами”, вторая — “существенно неустановившимися задачами”.

1. Постановка задачи

В декартовой системе координат x,y рассматривается нестационарное движение идеальной несжимаемой жидкости со свободной границей. Область D(t), занятая жидкостью, ограничена гладкими поверхностями S(t), Q(t) и W. Граница S(t) является свободной поверхностью жидкости, Q(t) — граница движущегося тела, граница W представляет собой твердую стенку. Предполагается, что в области D(t) происходит безвихревое потенциальное движение однородной жидкости, имеющей постоянную плотность р. Движение жидкости порождается наличием внутри жидкости движущегося горизонтально по прямолинейному дну со скоростью U (t) полукругового цилиндра диаметром D и наличием массовых сил, обладающих потенциалом gz (g = const). Давление Pa на свободной границе постоянно и равно нулю. В начальный момент времени t = 0 S (t) = S (0). Область

течения представляет собой канал глубиной H. Схема течения показана на рис. 1.

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

ДФ = 0, x(x,y) G D(t), (1)

dx

— = УФ, x(x, y) G S (t), (2)

© К.Е. Афанасьев, 1998.

Рис. 1.

^ - 2 I I2 +(У - !)/Гг2 = ° Х(х,у) е 5(*),

(УФ,п) = °, Х(х,у) е Ж,

(3)

(4)

(УФ,п) = —и(¿) ео8(п, х), Х(х,у) е ф(£). (5)

Задача (1) - (5) является нелинейной, так как положение свободной границы в каждый момент времени неизвестно. Кроме того, на поверхности 5(¿) должно выполняться нелинейное граничное условие, а именно, интеграл Коши — Лагранжа. В данной постановке в качестве характерных линейных величин выбирались глубина канала , скорость цилиндра V = тах |и(¿)|, п — вектор внешней нормали к границе области течения. Число Фруда

Рг2

V2

(6)

2. Метод решения и вычислительная технология

Для решения краевой задачи (1)-(5) используется алгоритм движения по времени, который описан в работах [1-3]. На каждом шаге по времени методом граничных элементов [4, 5] решается граничное интегральное уравнение (записанное по границе Г= 5(¿)и^(¿)иЖ), к которому сводится уравнение Лапласа внутри области О(£).

Используя третью формулу Грина, в плоском случае можно записать:

еФ(£)+ / ф(х)^ (£,х)^Г(х) = [ С'(£,ж)^Г(ж),

]г дп(х)

(7)

1Г дп(х)

где Г(х) — граница области О(£), е — угол между касательными к границе в точке £, Ф(х) — гармоническая функция, п — внешняя по отношению к области единичная нормаль поверхности, С;(£,ж) — фундаментальное решение уравнения Лапласа:

^/(£>х) = — 2п1п |£ — х|'

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

N N max ФТ — ФЧ max VT — V4 max VyT — Vy4

max |ФТ| max | VT | max VyT

72/30 145/60 290/120 580/240 8.4E — 03 2.2E — 03 5.4E — 04 2.2E — 04 4.1E — 03 8.4E — 04 7.2E — 04 7.0E — 04 6.6E — 02 9.9E — 03 9.3E — 03 6.4E — 03

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

Точность расчета и сходимость метода проверялись по тесту, предложенному в [8]. Требуется найти решение уравнения Лапласа в области D = {0 < x < 2п; —1 < y < sin(x)}, й дФ 0

в которой на дне и вертикальных стенках ставится условие непротекания —— = 0, а на

дп

верхней границе — условие Ф(х,у) = — cos(x)ch(y + 1), правая часть в котором является гармонической функцией. Численные значения нормальной производной ФЧ, найденные методом граничных элементов, сравниваются с точным решением: Ф^(х,у) = cos (х)

— (sin(x)ch(y +1) — sh(y + 1)). В таблице приведена относительная погрешность

у 1 + cos(x)2

точного и численного значений функции Фп от числа узлов N по границе с указанием числа узлов на свободной границе Ng. Вектор скорости (VX, Vy) в каждом узле на свободной поверхности вычислялся методом, изложенным в работе [7], и сравнивался с точным значением Vr = sin(x)ch(y + 1), V,? = — cos(x)sh(y + 1).

3. Существенно нелинейные течения

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

В этой задаче наряду с геометрическими характеристиками представляют интерес и физические особенности возникающей волны. Опрокидывающаяся волна вследствие быстроты происходящего процесса, значительных ускорений и больших искривлений свободной поверхности представляет собой наиболее трудный объект для исследования волновых явлений. В данной проблеме характерным является то, что перед обрушением элемент водной поверхности становится вертикальным, а затем гребень разрушается, образуя струю, направленную в сторону движения волны. Последнее хорошо видно из расчетов, приведенных на рис. 2. Эти расчеты проводились в области [— 12Н, 6Н] по оси ОХ при глубине = 1. Для задачи, приведенной на рис. 2, а, О = Н, Гг = 0.2, на рис. 2, б — О = 0.5Н, Гг = 0.5. Количество узлов на свободной границе было 68, всего в области — 190. При этом точки на свободной границе в начальный момент времени специальным образом сгущались в

а

у/я

1.00

-2.60__________________-2.00__________________-1.40_________________-0.80_____________Х/Н

х/н

Рис. 2. Эволюция свободной границы в моменты обрушения.

окрестности вершины цилиндра.

В обзоре [9], где рассматриваются уединенные волны, говорится о том, что в опрокидывающейся волне можно выделить следующие три зоны, которые появляются еще до того, как волна станет вертикальной (рис. 3): 1 — скорость отдельных частиц жидкости превышает волновую скорость; 2 — в передней зоне волны имеется тонкий слой, в котором ускорение частиц значительно выше ускорений в остальной волне; 3 — на задней поверхности волны появляется довольно обширная зона, в которой ускорение частиц очень мало. Проведенные в работе [10] исследования задачи об опрокидывании уединенной волны подтверждают наличие существования этих зон.

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

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

Н = 1, V = -1 Б/Н = 1.0, Х/Н = [-12,6], Рг = V/(дН)1!2 = 0.2

¿- 1.5807 2-1.5942

3 - 1.6181 4 - 1.6463

5 - 1.6666 6 - 1.6900

1

Рис. 3.

Поле скоростей ■'■"-V 1 -1 ^ \ 1 \ 1 1 / ^ % / 1 Ч/ і I 1 / 1

-4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1-0 -0.5 о.о X

Поле ускорений

1 V \ 1 Л ЗІ ИУЧ 1

\ 1 4 \ 1 1 У /

_і

4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 X

Рис. 4.

4. Существенно неустановившиеся течения

При обтекании препятствий ограниченным потоком жидкости в нелинейной постановке существуют режимы, зависящие от числа Фруда, отношения радиуса Я тела к глубине Н жидкости, при которых не существует стационарных режимов [11, 12]. В [12] показано, что при докритическом течении стационарный волновой цуг может существовать лишь при значениях чисел Фруда Рг, меньших некоторого критического числа Рг* (Я/Н). При сверхкритическом режиме обтекания также существуют числа Рг** > 1, лишь выше которых возможен выход на стационарный режим течения. В этой работе доказывается, что при Я/Н = 0.2 стационарных решений в диапазоне 0.6 <Рг< 1.2 не существует.

В [13] показано и экспериментально доказано, что при околокритическом движении препятствия из состояния покоя возникает течение, не выходящее на стационарный режим. Чтобы проверить эти утверждения и понять природу получаемых течений, был выполнен цикл расчетов (для различных чисел Фруда) по движению с единичной скоростью из состояния покоя полукругового цилиндра при отношении Я/Н = 0.2; при этом расчеты проводились в областях, границы которых показаны на рис. 5. Количество узлов по свободной границе варьировалось от 160 до 270 при общем числе узлов от 280 до 390.

Расчеты проводились в диапазоне чисел Фруда от 0.1 до 100 и были условно разбиты на

три группы: I — сверхкритические течения, при которых существует стационарный режим: 1.2 <Ег< 100; II — докритические течения, при которых также существуют стационарные решения: 0.1 <Ег< 0.6; III — околокритические течения, при которых не существует стационарных решений: 0.6 <Ег< 1.2. Наиболее изученными являются течения группы I (рис. 5, а). Здесь известно, что при установившемся течении над телом существует горб, а вверх и вниз по потоку возмущения практически не распространяются. Теоретически и экпериментально это показано в работах [1, 2, 14]. При этом существует общая для этой группы закономерность выхода на стационар: в начальный момент движения за телом появляется возмущение, которое, слегка деформируясь, сносится вниз по потоку. Над телом формируется горб, который по истечении некоторого времени не изменяет своей формы и движется вместе с телом.

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

Я г=1

где Q — контур тела, п — вектор нормали к поверхности тела, N — количество точек на теле. Давление Р определяется из интеграла Коши — Лагранжа:

Р д Ф 1. |2 у

---1------1— IV Ф12 +——

р + дь +2|VФ| + Ег2

Технология вычисления частных (временных и пространственных) производных описана в работе [7].

Хорошим критерием выхода течения на стационарный режим для данной группы является сравнение формы волны, сформировавшейся над телом, с формой волны для чисто стационарной задачи. На рис. 6 для числа Фруда Ег= 2.0 представлены для сравнения две волны: сформировавшаяся над цилиндром по истечении некоторого периода времени (пунктир) и полученная в результате расчета стационарной задачи обтекания полукру-гового цилиндрического выступа потоком жидкости конечной глубины [2, 15] (сплошная линия). Видно хорошее качественное и количественное совпадение результатов. Максимальная относительная погрешность составила 0.2 %.

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

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

5. Законы сохранения

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

[ дФ

Ш = % дХ (8)

Согласно теории, величина Ш должна равняться нулю. Значение интеграла (8) вычислялось на каждом шаге по времени и находилось в пределах от 10-6 до 10-4, а при опрокидывании волны не превышало 10-3.

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

М = (к(х,і) — Н (х))^х,

ь

а

где у = к(х,і) — уравнение свободной границы, у = H(х) — уравнение дна. Этот закон выполняется на всех режимах течения с точностью до 2 %.

б

Рг = 0.1

Эволюция свободной поверхности

-50.00 -25.23 5.56 33.2

X

^•10

Рг = 0.3

^•10

Рг= 0.5

^•10

Рис. 5. а — сверхкритические течения (группа I), б — докритические течения (группа II), в — неустановившиеся течения (группа III).

Y > Fr = 2.0

1.08-

Нестационарная задача

Стационарная задача

1.06- ' / \

1.04- - \

1.02- У V

-4.0 -3.0 -2.0 -1.0 0 1.0 2.0 3.0 4.0 X

Рис. 6.

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

[1] Афанасьев К.Е., Афанасьева М. М., Терентьев А. Г. Исследование эволюции свободных границ методами конечных и граничных элементов при нестационарном движении тел в идеальной несжимаемой жидкости. Изв. АН СССР, МЖГ, №5, 1986, 8-13.

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

[3] Лаврентьев М. А., Шабат Б. В. Проблемы гидродинамики и их математические модели. Наука, М., 1973.

[4] Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. Мир, М., 1987.

[5] Бенерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках. Мир, М., 1984.

[6] Бахвалов Н.С. Численные методы. T. 1. Наука, М., 1975.

[7] Афанасьев К. Е., Самойлова Т. И. Техника использования метода граничных элементов в задачах со свободными границами. В “Вычислительные технологии”, ИВТ СО РАН, Новосибирск, вып. 7, №11, 1995, 19-37.

[8] Петров а. Г., Смолянин В. Г. Расчет нестационарных волн на поверхности тяжелой жидкости конечной глубины. ПММ, 57, №4, 1993, 137-143.

[9] Перегрин Д. Разрушение волн на отлогих берегах. Нелинейные волновые процессы. Механика: Новое в зарубежной науке. Мир, М., 42, 1987, 37-71.

[10] VlNJE T., Brevig P. Numerical simulation of breaking waves. Adv. Water Resour, 4, 1981, 77-82.

[11] Стурова И. В. Численные расчеты в задачах генерации плоских поверхностных волн. Препр. ВЦ СО АН СССР, №5, Красноярск, 1990.

[12] Forbes L. K., Schwartz L.W. Free surface flow over a simicircular obstruction. J. Fluid Mech, 114, 1982, 229-314.

[13] Wu T. Y.-T. Generation of upstream advancing solitons by moving disturbances. Ibid, 184, 1987, 75-99.

[14] Долина И. С., Ермаков С. А., ПЕлиновский Е. Н. Смещение свободной поверхности жидкости при обтекании цилиндра. ПМТФ, №4, 1988, 48-51.

[15] Гузевский Л. Г. Обтекание препятствий потоком тяжелой жидкости конечной глубины. Динамика сплошных сред с границами раздела. Чуваш. гос. ун-т. им. И. Н. Ульянова, Чебоксары, 1982, 61-69.

Поступила в редакцию 14 октября 1995 г., в переработанном виде 3 ноября 1997 г.