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

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

Том 4, № 3, 1999

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ДОСТАТОЧНО ДЛИННЫХ ДВУМЕРНЫХ ВОЛН НА ВОДЕ В БАССЕЙНАХ С ПОЛОГИМ ДНОМ *

А. А. Литвиненко, Г. А. Хабахпашев Институт теплофизики СО РАН, Новосибирск, Россия

The work is devoted to the numerical implementation of a nonlinear evolutionary equation for three-dimensional perturbations of the free surface of not very deep layers of the incompressible viscous fluid. The calculations are shown adequately to describe the processes under study. In particular, the obtained results are in good agreement with the experimental data about the plane wave run-up on a gentle slope.

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

1. Постановка задачи и исходные уравнения

Будем считать, что, во-первых, стационарные составляющие течения жидкости равны нулю, во-вторых, "длина волны" А существенно больше, а амплитуда возмущения Па значительно меньше равновесной глубины слоя К (К/А ~ е1/2, Па/К ~ е, где е — малый параметр), в-третьих, капиллярные эффекты не велики (число Бонда Во = рдК2/а > 1, здесь р — плотность, д — ускорение свободного падения, а — поверхностное натяжение), в-четвертых, неподвижное твердое дно является слабонаклонным (УК ~ е3/2, где оператор градиента V определен в горизонтальной плоскости), и наконец, в-пятых, появляющиеся пограничные слои остаются тонкими, т. е. время прорастания пограничного слоя на всю

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-01766, СГП ВНШ, грант №96-15-96314, и СО РАН, ИГ-43-97.

© А. А. Литвиненко, Г. А. Хабахпашев, 1999.

толщину жидкости много больше характерного времени волнового процесса (число гидродинамической гомохромности Ио^ = vtw/к2 ~ е2, здесь V — кинематическая вязкость жидкости). Следовательно, возникающее течение будет потенциальным всюду за исключением узкой придонной области. Исходя из этих предположений и обычных краевых условий в работе [5] выведена система двух уравнений для возмущения уровня воды п и осредненного по глубине вектора горизонтальной составляющей скорости жидкости (и):

д'ц

(к + п)

д (и)

дь

дь

+ V дп +

+ (к + п)^ ■ (и)) + ((и) ■ V)(h + п) = 0,

(и) 2

а,

Р

к2 2 --V2-

з дг

д (и)

gV

п ЛЬг ^Ъ—Гг

ио

(1)

(2)

где и0 — начальное возмущение скорости воды. Дополнительное предположение о том, что нелинейные волны двигаются преимущественно в одном (хотя и произвольном) направлении, позволило свести эту систему к одному эволюционному уравнению [5]:

- дк^2ц - 3 gV2п2 - к2(3

1

ВО

V2ё? - д^п ■ Vk)+

дЬ2

+/п1д

V2n 1 дц

,_ + —= ^ I = 0. (3)

/с л/Ь-й к^Ь дЬ 1 К у

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

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

г

с

г

2. Численная реализация

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

3

Лггпп - дк(Лхх + Луу)пп - 2д(Лхх + Луу)(п2)П-

-д(ЛхппЛхк + ЛуппЛук) + к2 (3 - ВО) Л«[(Лхх + Луу)пп] = Гп,

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

Дх2Ду2 4 /1 1 \ к2 д^к| /Дх2Ду + Ду2Дх\ 2< Дх2 + Ду2 + ^3 - ВОу "V Дх2 + Ду2 )

3

дк + 2дпа

Несомненным достоинством уравнения (3) является тот факт, что данная неявная разностная схема обладает большей устойчивостью, чем аналогичная разностная схема для более простого уравнения (уравнения (3) без дисперсионного члена). Если рассмотреть случай одномерного движения, то условие устойчивости примет вид т2 < а Ах2 + в , где а и в — некоторые константы, зависящие только от ускорения свободного падения и глубины бассейна. Отметим, что устойчивых явных схем для уравнения (3), а также для другой модификации этого же уравнения с дисперсионным членом вида дк3( 1/3 — 1/Бо)У4п построить не удалось.

Решение итоговой системы алгебраических уравнений, возникшей после разностной аппроксимации уравнения (3) по данной схеме, в проведенных расчетах выполнялось методом Зейделя. Скорость сходимости итерационного процесса высока, так как имеется строгое диагональное преобладание. Процесс сходится до требуемой точности за 5-10 итераций, что вполне достаточно для исследуемых задач. Отметим также, что матрица итоговой системы уравнений зависит только от топографии бассейна и не изменяется в процессе решения. Поэтому можно один раз найти обратную матрицу для данного бассейна и использовать ее для различных начальных данных. Время перехода на следующий слой при расчете будет почти таким же, как и для явной схемы.

В проведенных расчетах значения искомой функции п в краевых узлах брались равными ее значениям в ближайшем внутреннем узле. Это соответствует условию отражения от вертикальной стенки. В качестве начального возмущения (на первом и втором слоях по времени) использовалось солитонное решение уравнения (3) при УК = 0 и V = 0:

где и = ^дК( 1 + п*), L = 2Нл/(1 + п*)/3п*, П* = Па/К.

Эта разностная схема тестировалась на аналитическом решении (4). Отметим, что на расстояниях 3Ь от вершины возмущение составляет 2.5 % от амплитуды. Результаты теста показали, что решение сходится не менее чем со вторым порядком по всем переменным (и по пространству, и по времени). Производился также анализ роста погрешности решения со временем (рис. 1). При этом ошибка нормировалась на амплитуду волны (шаг по времени был 0.05 с, шаг по пространству 0.05 м, глубина воды 0.01 м, амплитуда волны 0.001 м и, следовательно, характерный горизонтальный размер солитона 6Ь = 2.58 м). Таким образом, для того, чтобы погрешность была менее 1 % относительно амплитуды одиночной волны на временах около 10 секунд необходимо 40-60 узлов сетки на ее длину.

Напомним, что в эволюционном уравнении (3) член, отвечающий за учет вязкости жидкости, имеет следующий вид:

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

П = Па/ СОй К2 [(х — и1)/ь],

(4)

О 2 4 6 8 10 12 14 16 18 20 22 24 t, с

Рис. 1. Зависимость погрешности вычислений 6 от времени исследуемого процесса I.

аппроксимацию оператора Лапласа со вторым порядком Кц:

f * VVt fkT V2ndt,

'0 \A - ti k=1J(k-1)T \Zt - ti

± Г V2nkl(k - 1)T - 2nk-1(kT - *) dt, + 0(t2) =

k=1J(k- 1)t Ф - ti

£ Г* Л |(k - 1)t - t,] + Ank-1(kT - ti) du + O(T2 + Дж2 + Ду2).

(k— 1)t yt - ti

I.......7/г=п-'"2 ' л х2

к=1

После чего интегралы вычисляются в явном виде и получается итоговая формула:

п

рП = [_2 {кт^—кТ - (к - 1)т^г - (к - 1)т | -

к=1

- — Ш - кт)3/2 - [1 - (к - 1)т]3/2) (Лпк - Лпк)-3т

-2 (к7 - дД - (к - 1)т) (кЛпк - (к - 1)Лпк) .

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

3. Численные результаты

Были выполнены расчеты, результаты которых можно сравнить с опытными данными. Один из тестовых экспериментов проводился в U.S. Army Engeneer Waterways Experimental Station [7]. Бассейн имел 23.3 м в длину и 0.45 м в ширину, глубина воды до начала подъема

дна была 0.218 м. Глубина невозмущенной жидкости изменялась вдоль оси х:

0.2180

при х < 15.04,

0.218 + (0.136 - 0.218)(х - 15.04)/(19.40 - 15.04) 0.136 + (0.116 - 0.136)(х - 19.40)/(22.33 - 19.40) 0.116 + (0.057 - 0.116)(х - 22.33)/(23.23 - 22.33)

при 15.04 < х < 19.40, при 19.40 < х < 22.33, при 122.33 < х < 23.33.

Вдоль всего бассейна располагался ряд датчиков (рис. 2). Движением волнопродуктора создавалась плоская одиночная волна. При численных расчетах в качестве начального возмущения бралось решение в виде солитона (4) с амплитудой той волны, что и в опыте. На рис. 3, а начальная амплитуда равна 0.0082 м, на рис. 3, б — 0.058 м. Мареограм-мы выведены для датчиков с четвертого по восьмой (х4 = 12.64, рис. 3, а; х4 = 14.06, рис. 3, б; х5 = 15.04, х6 = 17.22, х7 = 19.40, х8 = 20.86). Графики показывают хорошее соответствие численных результатов и экспериментальных данных. Отметим, что вклады последних членов в уравнении (3), определяемых градиентом невозмущенной глубины и вязкостью жидкости, были пренебрежимо малы, а наклон дна фактически учитывался изменением как фазовой скорости очень длинных линейных волн, так и коэффициента при дисперсионном члене.

Для сравнения с другим экспериментом, приведенным в [7], численно исследована задача о воздействии конического острова на плоскую одиночную волну. Бассейн имел длину 25 м и ширину 30 м, а глубина воды равнялась 0.32 м. В его середине находился остров, представляющий собой усеченный конус высотой 0.625 м. В связи с тем, что заплеск волны на наклонный берег пока не рассчитывался, при вычислениях остров был заменен на затопленный усеченный конус (высотой 0.2 м) с поставленным на него вертикальным цилиндром (рис. 4). Для реализации отражения волны от цилиндрической части острова значение высоты возмущения п в граничном узле цилиндра бралось равным среднему арифметическому значений п в соседних узлах, лежащих за пределами цилиндра. В узлах внутри цилиндра п задавалась равной нулю. На рис. 5 показано начало взаимодействия волны с островом (а), увеличение амплитуды волны при обходе острова и возникновение отраженной волны (б), смыкание основного возмущения за островом (в), образование "шейки" между прошедшей волной и островом (г), наличие "ласточкина хвоста" за фронтом волны (д).

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

Здесь Ш — параметр, характеризующий ширину хребта. Хорошо видно, что для очень широкого хребта (Ш = 9, рис. 6, ') над ним наблюдается лишь некоторое отставание

(К/2) ехр[-(15 - х)2/2] ехр[-(20 - у)2/Ш2] при х < 15, (К/2) ехр[-(20 - у)2/Ш2] при х> 15.

I.

4

X

I.

6

I.

8 9 1®

I_и.

Рис. 2. Схема экспериментального бассейна и расположения датчиков [7].

V, м 0.010 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 -0.001 -0.002

_ л—V.

7-/ / ч \ / 7 \\

/ т \\ // // / \\ \ '1 \ !1 \

\\ \ ч // ! Vу/ \

/ 5\ I 61 !1 1 < / 7УУ

I / \\ 1 1 /Ъ ■ / \ 1 \

\ / /' Ч Д\ № // V 11 \ V 1

: / / !/ / // // \\ / \ ^ У \ ^ \ / ' \ / \\ / / \\ \\ 1

// г / \ \ Т 4 / / \ч л* / / \\ \Ч : Л н

г _._. --■" г \ ^ -- ■—

1 V ЛЛ А/ ! Ш-Н" \ / \ ;

г"

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 г,с

Г], м б

ц

¡1 ' / М л

I4 / 6 171 8

V V 1 1 ч \ || 1

м \\ ¡1 V 1 \

Л \ 1 // \ (

J У \\ Чг - V У 0 XV, ; --. --. :

0.090 0.080 0.070 0.060 0.050 0.040 0.030 0.020 0.010 0.000 -0.010 -0.020

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 t,c

Рис. 3. Эволюция нелинейной плоской уединенной волны при накате на пологий откос. - — экспериментальные данные [7],---— результаты расчетов по уравнению (3).

3.6 м

2.8 м

2.4 м

>

Уровень воды 32 см

Рис. 4. Схема модельного острова.

возмущения, увеличение его амплитуды и образование дисперсионного "хвоста" (рис. 6, а), а для достаточно узкого хребта (Ш =1, рис. 6, ') имеют место, кроме того, слабые волны типа корабельных (рис. 6, б).

Для случаев, когда вязкость жидкости существенна, исследовался вопрос о конкуренции между влияниями уменьшения глубины бассейна и нестационарного трения о дно. Как известно, результатом первого эффекта является увеличение амплитуды волны, а второго — ее уменьшение. На рис. 7 приведены мареограммы в точках, расположенных на расстоянии 5 м друг от друга. В случае а продемонстрирована ситуация, при которой эти эффекты практически уравновешивают друг друга (амплитуда одиночной волны не меняется со временем, хотя за основным фронтом четко видны вязкий и дисперсионный "хвосты"). В случае б наклон дна становится больше и диссипативные потери не могут противодействовать нарастанию амплитуды волны. Отметим, что значение кинематической вязкости V = 100 мм2/с, использовавшееся в расчетах, соответствует вязкости масел и концентрированных водных растворов глицерина.

4. Заключение

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

Авторы выражают признательность П. И. Гешеву, О.Ю. Цвелодубу и Л. Б. Чубарову за полезные обсуждения данной работы и ценные советы.

Рис. 5. Взаимодействие волны с островом при Ь = 2.5 (), 4.2 (), 6.0 (), 7.5 () и 9.0 с ().

Рис. 6. Воздействие подводного хребта на нелинейную плоскую уединенную волну.

Г], м

0.040

0.030

0.020

0.010

0.000

: д 1 2 3 4 5

\ 1

|

V ) ' { )

-0.010

0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 г,с

Г}, м

0.040

0.030

0.020

0.010

0.000

5

: а1 2 3 4

\

\ \ \

V ) ) V 1

-0.010

0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 t,c

Рис. 7. Эволюция уединенной плоской волны при учете нестационарного трения жидкости о дно. До третьего датчика дно горизонтально (Н = 0.2 м), а далее имеет постоянный наклон по оси х: дН/дх = -1/200 (а) и дН/дх = -1/100 (б).

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

[1] УизЕм Дж. Линейные и нелинейные волны. Мир, М., 1976.

[2] Марчук Ан. Г., Чубаров Л. Б., Шокин Ю. И. Численное моделирование волн цунами. Наука, Новосибирск, 1983.

[3] Вольцинер Н. е., Клеванный К. А., Шлиновский Е. Н. Длинноволновая динамика прибрежной зоны. Гидрометеоиздат, Л., 1989.

[4] Шокин Ю. И., Чубаров Л. Б., Марчук Ан. Г., Симонов К. В. Вычислительный эксперимент в проблеме цунами. Наука, Новосибирск, 1989.

[5] Хабахпашев Г. А. Нелинейное эволюционное уравнение для достаточно длинных двумерных волн на свободной поверхности вязкой жидкости. Вычислительные технологии, 2, №2, 1997, 94-102.

[6] КовЕня В. М., ЯнЕнко Н. Н. Метод расщепления в задачах газовой динамики. Наука, Новосибирск, 1981.

[7] Proceedings of the International Workshop on Long-Wave Runup Models. Friday Harbor, San Juan Island, WA, USA, September 12-16, 1995.

Поступила в редакцию 5 октября 1998 г., в переработанном виде 17 ноября 1998 г.