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

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

Том 13, № 1, 2008

Численное исследование стокового механизма

генерации волн в рамках модели

_ «.» *

потенциальных течении жидкости*

Д. Б. ДАМБИЕВА, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

In the paper some results of a numerical investigation of a process of generation of the surface waves by flowing of the fluid into cracks that appear in the bottom during underwater earthquakes are presented. The numerical algorithm is based on the finite-difference scheme for the equations in a moving coordinates which describes potential ideal flows with a free boundary.

Введение

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-05-64869) и Программы интеграционных фундаментальных исследований СО РАН (проекты 2006.113 и 2006.28).

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

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

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

Рассматривается плоскопараллельное потенциальное течение идеальной несжимаемой жидкости в бассейне, ограниченном снизу дном Го, заданным функцией y = -h(x), а с боков — вертикальными непроницаемыми стенками Г1 и Г2 с абсциссами x = 0 и x = L, где x, y — координаты точек в декартовой системе координат с вертикальной осью Oy и горизонтальной осью Ox, лежащей на невозмущенной свободной поверхности жидкости (рис. 1).

Свободная поверхность жидкости Г/ подлежит определению. Будем считать, что она описывается однозначной функцией y = n(x,t), где t — время. Кроме функции п требуется найти потенциал скорости — функцию &(x, y, t), удовлетворяющую в области H(t), занятой жидкостью, уравнению Лапласа. На боковых стенках бассейна ставится условие непротекания:

dx

x=0

0, дг

dx

0.

x=L

(1)

На дне бассейна Г0 задается условие для нормальной составляющей вектора скоро^ сти и = (и, V) = У ^

д^

дп '

где V = 0 для непроницаемой части дна и V = V(х, у,£) > 0 — для проницаемой д^/дп = У^ • п, У = (д/дх, д/ду), п — единичный вектор внешней нормали к Г0. На свободной поверхности Г/ ставятся условия: кинематическое

П + и Пх - V = 0, х =(х,у) е Г/,

динамическое

&+2 iv^i2+п

о, x е Г

/ •

y=n(x,t)

' —\

X L

y=-h(x)

Рис. 1. Бассейн с частично проницаемым дном

Предполагается, что в начальный момент времени t = 0 жидкость покоится (<(x,y,0) = 0) и ограничена сверху невозмущенной свободной границей (n(x, 0) = 0).

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

x = — ~ = — i = t ^ h = — — u = u v= v ~ = < ho' h0' v/h^' h^ h^ \J~gh0 \fghh0 hov/gho'

где g = const — ускорение свободного падения, h0 — характерная глубина, символом " отмечены безразмерные величины (в приведенных ранее формулах этот символ опущен).

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

t = q0, x = x(q1,q2,q0), y = y(q\q2,q0), (2)

которое в каждый момент времени t устанавливает взаимно однозначное соответствие между физической областью Q(t) и единичным квадратом Q в плоскости q1Oq2. Будем считать, что участки границы Г1 и Г2 области П — образы при отображении (2) боковых сторон Q, дно Г0 — образ нижней стороны, а свободная поверхность Г^ — образ верхней стороны Yf квадрата Q.

Для более компактной записи уравнений в криволинейной системе координат используем следующие обозначения:

Ф = kll — + , Ф = k21^-T + k22 TJ 2,

dq1 dq2 dq1 dq2

g22 g12 g11

k11 = ~J, k12 = k21 =--J, k22 = ,

g11 = X?1 + yjl, g12 = g21 = xq1 xq2 + yq1 yq2 , g22 = x?2 + yj2 ,

где J — якобиан преобразования (2). Тогда в криволинейных координатах математическая постановка задачи будет заключаться в определении потенциала <(q1,q2,t) и функции n(q1,t), описывающей свободную границу, удовлетворяющих в Q нелинейному уравнению эллиптического типа:

дФ1 дФ2 т о

^г + -q^ = 0, q =(q1,q2) g Q, (3)

nt + vrnqi - v = 0, q G Yf (4)

<t - uxt - vyt + + П = 0, q G Yf (5)

условиям на верхней стороне квадрата Q, условию непротекания на его боковых сторонах

Ф1 I i 0= 0, Ф1 I i = 0

Iq1=0 ' Iq1 = 1

кинематическому

и динамическому

и краевому условию на нижней стороне области

ф2 = у/дЦ.

я2=0

В уравнениях (4), (5) декартовы компоненты скорости выражаются через потенциал по формулам

^я1 У Я2 - ¥я2 У я1 -<РЯ1ХЯ2 + ¥я2 хяг (п\

и = ——-—-—-——, V =--——-——, (6)

и и

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

дд1 дд1 Од1

V = 1Й + и5Х + 4 € 7/" 2. Алгоритм расчета на подвижной сетке

Область Q покрывается равномерной прямоугольной сеткой 5н с шагом На и количеством узлов Ыа в направлении оси Ода (а = 1, 2). Совокупность внутренних узлов qj сетки (л — мультииндекс, л = (]1,]2)) будем обозначать 5^..

Пусть на п-м слое по времени сетка (хп, уп) построена и на ней вычислены значения сеточных функций ], рп. Вычисление решения на (п + 1)-м слое состоит из нескольких этапов. Сначала вычисляем шаг по времени по формуле: тп = кзапАхт-т, где кзап — коэффициент запаса, 0 < кзап ^ 1, Дхт;п минимальный в горизонтальном направлении шаг сетки (хп,уп).

Далее вычисляется значение потенциала на верхней стороне области П на основании разностной аппроксимации динамического условия (5):

^п+1 - IV]2

И_1] ипх . vnу . | 1 1 + пп =о

- и] хь,з - vj Ш,] + 2 + ] = °

тп 2

где л = (^,N2), л = 1,...,^1 - 1,

хп - хп-1 уп - уп-1

] ] У] У]

хг,з = ~-—, Уц = ~-— ■

Тп— 1 тп— 1

тп-1 — шаг по времени между временными слоями ¿п—1 и ¿п, а ип и V]1 вычисляются на основе аппроксимации формул (6).

После этого итерационным методом решается система разностных уравнений

( X] акРк) . = Я]> qj е 5н \ 1/, н, (7)

к ]

аппроксимирующих уравнение (3) во всех узлах сетки, за исключением узлов 7/, н, лежащих на верхней стороне квадрата 5, где потенциал уже определен. Для внутренних узлов разностные уравнения для потенциала имеют девятиточечный шаблон [3]. Для граничных узлов, не принадлежащих 7/, н, разностные уравнения (7) имеют шеститочечный или четырехточеный (для угловых узлов) шаблон. Для узлов сетки, лежащих на прообразе дна — нижней стороне квадрата 5, — правые части разностных уравнений (7) вычисляются по формуле Я] = -к1 (^^уй) .. В остальных узлах сетки Я] = 0.

Рис. 2. Расчетная сетка

Затем находим положение свободной поверхности • Для этого кинематическое условие (4) аппроксимируется явной противопоточной схемой с учетом знака контрава-риантной компоненты скорости V1. При этом в формулах (6) для вычисления декартовых компонент скорости используются значения потенциала на новом временном слое.

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

3. Результаты вычислительных экспериментов

Приведем результаты моделирования поверхностных волн в экспериментальной установке, описанной в работе [2]. Экспериментальный лоток длиной 3.1 м с вертикальными непроницаемыми стенками имеет на горизонтальном дне пластину шириной 0.28 м с шестью прорезями, каждая шириной по 0.02 м. Пластина располагается на расстоянии 1.48 м от левой стенки лотка и закрыта подвижной заслонкой. В начальный момент времени заслонка открывается и начинается сток воды в прорези. В некоторый заданный момент времени заслонка возвращается в исходное состояние, перекрывая слив воды. Эксперимент проходил при начальной глубине воды к0 = 0.13 м.

При расчетах мы вместо шести щелей рассматривали одну, по ширине равную суммарной ширине прорезей: й = 0.12 м. Поскольку дно бассейна горизонтально, то функция V из граничного условия (1) равна абсолютному значению |^о| скорости вытекания в вертикальном направлении. В экспериментальной работе [2] значение у0 не указано, поэтому в численных расчетах оно выбиралось из некоторого диапазона значений.

На рис. 3 изображены мареограммы (зависимости полной глубины Н = ц + к0 от времени) в центре щели (х = 1.54 м) при скорости вытекания |^0| = 0.3 м/с, полученные на различных расчетных сетках. Здесь Н0 — число узлов в горизонтальном направлении, приходящихся на проницаемую часть дна. Видно, что при измельчении сетки имеет место сходимость численного решения. При этом кривые 3 и 4 почти совпадают и дальнейшее измельчение сетки не вызывает изменения графиков мареограмм, поэтому все последующие расчеты проводились на сетке с числом узлов N1 = 200, N = 15, N = 20.

В рассматриваемом численном эксперименте заслонка закрывалась в момент времени £ =1 с. Из рис. 3 видно, что мгновенное прекращение слива приводит к возникновению волны значительной амплитуды над центром щели. Причиной такого всплеска является резкая перестройка течения в момент прекращения стока: до этого момента жидкость стекала в щель, поступая к ней с двух сторон в виде пары сформировавшихся потоков, которые после закрытия заслонки трансформируются в два встречных потока, движущихся над горизонтальным непроницаемым дном.

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

Рис. 3. Мареограммы в центре щели, полученные на различных расчетных сетках: 1 — N1 = 50, N0 = 5, N2 = 5; 2 — N = 100, N0 = 10, Ы2 = 10; 3 — N = 100, Щ = 10, Ы2 = 20; 4 — N = 200, N0 = 15, Щ = 20

Рис. 4. Профили свободной поверхности в различные моменты времени: 1 — Ь = 0.1 с; 2 — г = 1 с; 3 — Ь = 1.25 с; 4 — Ь = 2 с

Рис. 5. Мареограммы в центре щели для скорости вытекания: 1 — = 0.1 м/с; 2 — = 0.2 м/с; 3 — = 0.3 м/с; 4 — = 0.4 м/с; 5 — |г0| = 0.5 м/с

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

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

На рис. 6 показаны мареограммы для одной и той же скорости вытекания |^0| = 0.3 м/с, но при разной ширине щели d. Видно, что графики имеют тот же характер, что и на рис. 5, т.е. увеличение ширины щели при фиксированной скорости |^0| эквивалентно увеличению скорости вытекания при фиксированной ширине щели. Отсюда можно сделать предположение о том, что картина течения, видимо, определяется не двумя параметрами |^0| и d, а одним — расходом жидкости

На рис. 7 показано сравнение экспериментальных данных с результатами, полученными численно при d = 0.12 м и |^0| = 0.55 м/с в рамках линейной [1] и нелинейной моделей потенциальных течений. В момент времени £ = 4.72 с — мгновенного перекры-

Рис. 6. Мареограммы в центре щели при разных значениях I: 1 — I = 0.06 м; 2 — I = 0.12 м; 3 — I = 0.18 м

Рис. 7. Свободная поверхность в окрестности центра щели: сплошная линия — эксперимент [2], штриховая — линейная модель [1], штрихпунктирная — нелинейная модель

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

Заключение

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

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

Авторы благодарят Л.Б. Чубарова, по инициативе которого было предпринято настоящее исследование.

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

[1] БАБАЙлов В.В., Дамбиева Д.Б., Хлкимзянов Г.С., Чубаров Л.Б. Численное моделирование стокового механизма генерации волн цунами // Тр. междунар. конф. "Вычислительные и информационные технологии в науке, технике и образовании". Павлодар: ТОО НПФ "ЭКО". 2006. Т. 1. С. 160-171.

[2] Левин Б.В., Королев Ю.П., Королев П.Ю. Стоковый механизм образования цунами. Южно-Сахалинск, 2006 (Препр. РАН. Дальневосточное отд-ние. ИМГиГ. 21 с.)

[3] ХАкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., ШокинА Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001. 394 с.

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