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

УДК 532.59

С.Г. Демышев, O.A. Дымова, Н.В. Маркова

Численное моделирование приливных волн и течений в бассейне с двумя проливами

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

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

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

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

© С.Г. Демышев, O.A. Дымова, Н.В. Маркова, 2006

66 ISSN 0233- 7584. Мор. гидрофиз. журн., 2006, N9 2

1. Рассмотрим прямоугольный бассейн (ахЬ) постоянной глубины с двумя проливами шириной 2/, который заполнен несжимаемой однородной невязкой жидкостью (рис. 1). Пусть на жидкость действуют диссипативные силы и сила тяжести. Уравнения, описывающие движение такой жидкости, в линейном приближении имеют вид [5]:

ди

■ о.

dv ' Л дС

dt ду dt

+ h

du dv дх ду

-о, (1)

где t - время, и и V - составляющие скорости по осям х и у соответственно, возвышение свободной поверхности, к - глубина бассейна, g - ускорение силы тяжести, ц - коэффициент диссипации. На твердой части границы выполняется условие непротекания, а в проливах распределение скорости и (х, у, /) выбирается следующим образом:

я{у-Ь])

и(0 ,y,t)=u(a,y,t) =

и0 cos

21

Sin СПt .

(2)

где щ - максимальное значение скорости в проливе, Ь\= Ы2, а - частота волны. Задачу будем решать с нулевыми начальными условиями:

и(х, у,0) = у(х, у,0) = С(х, у,0) = 0 . (3)

АУ

X

а

Рис. 1. Схема бассейна

2. Для данной постановки (1) - (3) аналитическое решение найти не удалось, поэтому для решения поставленной задачи используются численные методы. Расчетную область (прямоугольный бассейн) разобьем на боксы с постоянными сторонами Аде, Ду. Распределение переменных, которое указано на рис. 2, соответствует сетке С в терминологии работы [6]. Граница области проводится по граням боксов, где определены скорости. При переходе от дифференциальных уравнений (1) к дискретной постановке используется известная схема «чехарда» [7]. При этом конечно-разностный аналог системы (1) можно записать в следующем виде:

h ( „ „ \ h Ах

n-1

(5)

2At

MJ

-1/,

(4)

g (r _Cn L »

-v":1

lAt

i'j __ 8 (ftl \ , л

(6)

где индекс и - показывает временной уровень, А? - шаг по времени.

V

Ау

U

------ч f С

v.................................4 f..................................i

и

V

Ах

Р и с. 2. Пространственное распределение переменных на сетке

Анализ устойчивости системы (4) - (6) показывает [8], что при выбранных параметрах сетки схема устойчива при А/ < 15 с. Для предотвращения слабой неустойчивости, присущей схеме «чехарда», используется фильтр Ас-селина [9], который для каждой искомой функции имеет вид:

щп=у/п+ 0,5о(у/"-х - 2у/п + у/п+х),

где тильдой обозначена сглаженная функция, а и - параметр фильтра, который принимается равным 0,05.

3. Для проверки численной модели была решена задача, имеющая аналитическое решение.

Рассмотрим бассейн со сторонами а = 200 км, Ъ = 100 км и глубиной И = 200 м, покрытый равномерной сеткой. Шаг по оси х составляет 2 км, по оси у - 1 км. На боковых границах (х = 0, х = а) задается постоянная скорость щ = 0,5 м/с. Движение жидкости описывается линеаризованной системой (1). В данном случае трение не учитывается (/и = 0). Решение будем искать в виде периодических по времени функций:

u{x,y,t) = u (х)cos—у sin at, v(x, у, t) = v (x)sin — y sin at, b b

C(x9y,t) = ^(x)cos-y у cos at,

(7)

здесь и, v, £ - амплитудные функции волновых характеристик,/: = 1,2,3... Подставляя (7) в (1) и решая дифференциальные уравнения, находим:

= (Ах sin Хх + А2 cos Ajc) ,

<^(jc) = — (A¡ cos Ajc - A2 sin Xjc) , gb

v (x) =—(Ax eos he - A2 sin fac), A,

(8)

где «(о) = ЭД = Ио, t^¡,A2=u0, Р =

sin ka gh b

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

и (х) = С, sinh yjc + С2 cosh ух,

Ç(x) = ——(С, cosh ух* н- С о sinh ух), (9)

gy

V (jc) = — (Q cosh ух + С2 sinh ух), Y

где С2 = w0»

sinhya gA

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

u(x9y,0) = v(x9y90) = 0, ф,у,0) = l(x)cos^y (10)

о

и граничными условиями

и(0,t) = и(а, y9t) = u0cos^ys'mot, (11)

b

v(x,0,t)=v(x,b,t) = 0. (12)

При проведении эксперимента параметр к = 2, Д/=10с. Так как схема аппроксимации по времени в (4) - (6) трехслойная, то начальные поля задаются на первых двух шагах. С помощью (8) были посчитаны волновые характеристики для волн с периодом 1 ч, а по (9) - для периода 12 ч. Результаты расчетов, выполненных на несколько периодов волны, представлены на рис. 3, 4. Сплошной линией обозначено численное решение, пунктирной - аналитическое. В ходе численного эксперимента рассматривались девять точек бассейна. Максимальное совпадение численного и аналитического решений (до пятого знака после запятой) наблюдалось в центре проливов, т.е. в точках с координатами (0, 50) и (200, 50). Поведение исследуемых функций в целом практически не изменялось, поэтому на графиках представлена одна точка (20, 10). Из рис. 3 и 4 видно, что численная модель достаточно хорошо воспроизводит аналитическое решение. Для периода 1 ч в среднем отличие между аналитическим и численным решениями составляет не более 10 % как для возвышения свободной поверхности, так и для значения скорости. Для периода 12 ч разница не превышает 4 % для возвышения уровня и 6 % для скорости. Следует отметить, что в обоих случаях наибольшие различия возникают, когда искомые функции достигают своих максимальных значений. Таким образом, можно сказать, что отличие численного решения от аналитического незначительно.

0.44 0.2 о--0.2" -0.4-

С,М

и, м/с

0.04

-0.04

0

5

а

Р и с. 3. Изменение уровня (а) и составляющей скорости п точке

(20,10) для периода 1 ч (сплошная линия - численное решени :кое)

0.01-

-0.0 г

Л ^

V *

0.04-

-0.04

(,ч

т~

40

80

I

120

Р и с. 4. Изменение уровня (а) и составляющей скорости по оси * (б) со временем в точке (20, 10) для периода 12 ч (обозначения те же, что на рис. 3)

4. После тестирования разработанной программы были исследованы колебания жидкости в бассейне с двумя узкими проливами (рис.1). Использовалась схема (4) - (6) с граничными условиями (2) и начальными условиями (3). Ширина проливов 2/ = 6 км, максимальное значение скорости в проливе щ = 50 см/с. Численные эксперименты проводились для двух параметров диссипации щ = 0,1 а, ц2 = 0,2а и изучалось влияние изменения коэффициента щ 0 = 1,2) на время установления колебаний Г, амплитуду возвышения свободной поверхности Ао и амплитуду волновой скорости установившихся колебаний У0. При этом величина волновой скорости У(х, у, /) вычислялась по формуле

V = л]и2 +у2

Результаты численных расчетов приведены в табл. 1-3.

В табл. 1 представлены результаты расчетов Г, Ао и У0 для волны с периодом 1 ч в пяти точках бассейна. Из таблицы четко прослеживается влияние диссипативных сил на исследуемые параметры. Так, при увеличении ц время установления и амплитуда колебаний во всех точках бассейна уменьшаются. Причем при возрастании ц в 2 раза величина Т уменьшается в 1,9 раза в точках (50, 25) и (50, 75) и в 2 раза в точках (150, 25) и (150, 75) и почти в 1,9 раза в центре бассейна (100, 50). При этом амплитуда колебаний уменьшается в среднем на 5% в указанных точках. Была также изучена зависимость возвышения свободной поверхности от времени в точке (100, 50) при различных коэффициентах диссипации. Уровень жидкости в центре не поднимался выше 1 см. На исследуемые параметры существенное влияние оказывает и расположение точек. Отметим, что при любом ц наименьшее время для установления колебаний требуется для точек, находящихся в центре бассейна. Наибольшее время до установления проходит в точках (150, 25) и (150, 75). Для одного и того же ц время установления в точках с одинаковой координатой х отличается не более чем на 1 мин, в точках с одинаковой координатой^ это отличие составляет 1,5 ч для ц = и 30 мин для ц = ц2- В зависимости от координаты точки амплитуда установившихся колебаний для одинаковых ц изменяется не более чем на 7%. На амплитуду скорости параметр диссипации влияет слабо. Так, У0 изменяется не более чем на 10% в указанных точках. На величину скорости заметно влияет геометрическое положение рассматриваемой точки. Максимальная скорость 4 см/с достигается в центре, скорость 2 см/с - в точках (50, 25), (50, 75) и 1,6 см/с - в точках (150, 25) и (150, 75).

Таблица 1

Зависимость времени установления колебаний Т, максимальной амплитуды возвышения свободной поверхности А0 и скорости установившихся колебаний У0 от коэффициента диссипации \х для т = 1 ч

(х, у\ км

щ т А0, см У0, см/с

12ч 32мин 0с 18,582 2,231

Ц2 6ч ЗЗмин 30с 17,737 2,114

12ч 32мин Юс 18,326 2,233

Ц2 6ч 34мин 0с 16,886 2,104

14ч 2мин 0с 19,708 1,623

Ц2 7ч Змин 50с 18,292 1,614

14ч 2мин 20с 18,713 1,621

И2 7ч 4мин Юс 17,460 1,606

10ч ЗЗмин 0с 0,775 4,481

Ц2 5ч 35мин 40с 0,73 3,983

(50,25) (50,75) (150,25) (150,75) (100,50)

В табл. 2 сведены данные, полученные для т = 12 ч. Из таблицы видно, что увеличение коэффициента диссипации приводит к уменьшению времени, которое требуется для установления колебаний. Так, при увеличении ц величина Г уменьшается в 1,3 раза для точек (50, 25), (50, 75), в 1,6 раза для точек (150, 25), (150, 75) и в 1,4 раза для центральной точки (100, 50). Увеличение ц

в пределах \i\ < |i < ц,2 не оказывает заметного влияния на амплитуду установившихся колебаний А0, ее изменение не превышает 1% для каждой из рассматриваемых точек. Наименьшее значение А0 получено в центре бассейна. Помимо коэффициента диссипации, на величину Г также влияют координаты рассматриваемой точки. Для т = 12 ч колебания быстрее всего устанавливаются в точках (150,25) и (150,75): 125 ч 52 мин при \i = Hi и 78 ч 13 мин при (I = ц2. При этом в центральной точке Т =132 ч 1 мин и 96 ч 19 мин, в точках (50,25), (50,75) - 132 ч 1 мин и 120 ч 24 мин для |х, равных щ и ц2 соответственно. Таким образом, разница между максимальным и минимальным значениями Т составляет 6 ч для ц = щ и 42 ч для ц = ц2. На амплитуду волновой скорости V0 увеличение |х в указанных пределах существенного влияния не оказывает: изменение скорости не превышает 2%. При т = 12 ч максимальное значение V0 составляет 3,5 см/с. Положение точек на величину V0 также не влияет. Для одного и того же значения ц, скорость меняется не более чем на 2% во всех указанных точках.

Та блица 2

Зависимость времени установления колебаний Г, максимальной амплитуды возвышения свободной поверхности А0 и скорости установившихся колебаний Vo от коэффициента диссипации \i для т = 12 ч

(х, у), км И. 1 1 т А0> см F0, см/с

(50,25) Hi Н2 132ч 1мин 120ч 24мин 2,673 2,648 3,462 3,427

(50,75) Hi Н2 132ч 1мин 120ч 24мин 2,663 2,637 3,446 3,393

(150,25) Hl И2 125ч 52мин 78ч 14мин 2,534 2,567 3,475 3,403

(150,75) Hl Н2 125ч 52мин 78ч 14мин 2,524 2,556 3,459 3,387

(100,50) Hl Н2 132ч 2мин 96ч 19мин 0,053 0,053 3,395 3,462

В табл. 3 приведены данные, полученные для х = 24 ч. Из нее видно, что увеличение коэффициента диссипации в 2 раза приводит к уменьшению времени установления в 1,6 раза в точках (50,25) и (50,75), в 1,2 раза - в точках (150,25), (150,75) и в 1,15 раза в центре. На амплитуду установившихся колебаний изменение коэффициента диссипации (jii < ц < |х2) оказывает слабое влияние: увеличение А0 в этом случае не превышает 1% во всех рассматриваемых точках. Для одного и того же значения (х величина Т различна в указанных точках. При этом в точках с одинаковой координатой х отличие не превышает 1,5%, тогда как в точках с одинаковой координатой у оно составляет 10%. Для ц = \i\ быстрее всего колебания устанавливаются в центре бассейна, потом в точках (150,25), (150,75) и, наконец, в точках (50,25), (50,75). Для |i = ц2 картина выглядит несколько иначе: сначала колебания устанавливаются в точках (50,25) и (50,75), затем в точках (150,25), (150,75) и позже всего в центре. На амплитуду А0 расположение рассмотренных точек для

данных значений ц влияет слабо: А0 изменяется не более чем на 4%. Минимальное значение А0 наблюдается в центре. Изменение коэффициента диссипации и координат точек не оказывает заметного влияния на величину Уо. Так, максимальная скорость установившихся колебаний достигает 3 см/с, и ее изменение не превышает 1,5% в зависимости от исходных параметров.

Таблица 3

Зависимость времени установления колебаний 7*, максимальной амплитуды возвышения свободной поверхности А0 и скорости установившихся колебаний У0 от коэффициента диссипации ц для т = 24 ч

(х,у), км

и т А0, см Уо, см/с

380ч 14мин 1,267 3,329

Ц2 240ч 46мин 1,275 3,297

384ч 38мин 1,256 3,314

Ц2 240ч 46мин 1,269 3,282

Ц| 340ч 46мин 1,225 3,335

И2 276ч 45мин 1,227 3,314

348ч ЗОмин 1,217 3,321

Ц2 288ч 46мин 1,222 3,298

М1 335ч Збмин 0,024 3,342

Ц2 293ч 2мин 0,025 3,339

(50,25) (50,75) (150,25) (150,75) (100,50)

Сравнительный анализ результатов, полученных для трех значений т, позволил установить такие зависимости Т и А0 от периода волны: для 1 ч возвышение свободной поверхности достигает максимального значения 20 см за минимальное для трех случаев время (14ч при ]и=Ц1 и 7ч при Ц=Ц2). Для т = 12 ч величина Ао в исследуемых точках достигает максимального значения 2,6 см за 132 ч при Ц=Ц] и за 120 ч при Ц=Ц2- При т = 24 ч максимальное поднятие уровня составляет 1,2 см, и достигается это значение через 384 ч для Ц=Ц1 и 288 ч для Ц=Ц2. Для центральной точки величина А о не превышает 1 мм. Таким образом, увеличение периода волны приводит к уменьшению амплитуды установившихся колебаний и к увеличению времени их установления. Амплитуда скорости У0 достигает 4 см/с для волны с периодом 1 ч, а для 12 и 24 ч составляет около 3 см/с. Для волны с периодом х = 1 ч значение скорости находилось в интервале от 1,6 до 4 см/с в зависимости от координаты. Для колебаний с периодами 12 и 24 ч во всех точках, расположенных в центральной части бассейна, амплитуда волновой скорости одинакова. Таким образом, увеличение периода волны приводит к тому, что скорость установившихся колебаний равномерно распределяется в центральной части бассейна.

По данным расчетов были построены поверхности уровня и поля волновых скоростей в случае установившихся колебаний. Коэффициент диссипации равнялся 0,2а. Уровневые поверхности представлены на рис. 5. На рис. 5, а изображен пространственный профиль свободной поверхности ¿¡(х, у, /1) для т = 1 ч. Здесь >Т- момент времени, при котором | [х, у, /) | достигает максимальных значений в проливах (х = 0, х = 200 км, 47 км < у < 53 км).

Положительные и отрицательные отклонения уровня показаны относительно невозмущенной поверхности. Амплитуда волны в проливах при г= 1 ч достигает 54 см. На поверхности имеется три узловых линии. На рис. 5,6 представлена уровневая поверхность для г= 24 ч. Видно, что максимальное значение амплитуды 4 см также достигается в проливе, уровневая поверхность имеет одну узловую линию (х=100 км). Отметим, что при х- 12 ч амплитуда волны в проливах равна 8 см, а положение узловой линии сохраняется.

Р и с. 5. Поверхности уровня после установления колебаний для волн с периодами 1 ч (а) и 24 ч (б)

Векторные поля волновых скоростей представлены на рис. 6. На рис. 6,а изображено поле скорости V в указанный выше момент времени t\ для т= 1 ч. Размер и направление стрелки соответствуют модулю и направлению вектора V. Максимальный размер стрелки отвечает значению скорости 4 см/с. Из рисунка видно, что вектор волновой скорости распределен неравномерно для колебаний с периодом r= 1 ч. Максимальные значения скорости устанавливаются в центре (х = 100 км), а также при л; = 25кмих = 175 км. Минимальные скорости, близкие к нулю, достигаются при х = 75 км и х = 125 км. Для т= 24 ч (рис. 6, б) значение вектора скорости в центральной части бассейна (50 <х < 150 км) равно 3 см/с и его направление совпадает с j направлением скорости в проливе. Отметим, что при г= 12 ч вектор волновой

скорости в центральной части бассейна равен 3,5 см/с и имеет то же направление, что и в случае 24 ч.

км —^ 4 см/с

100 80-

60-

40 —

20-

(Г

->

->

-> ->

<— ^г ^г

<Г- ^г <— <-

к к

Р и с. 6. Структура течения после установления колебаний для волн с периодами 1 ч (а) и 24 ч(б)

В заключение сформулируем основные выводы:

- увеличение коэффициента диссипации в указанных пределах приводит к уменьшению времени установления колебаний;

- изменение коэффициента диссипации слабо влияет на амплитуду возвышения свободной поверхности и структуру течения для установившихся колебаний;

- время установления колебаний зависит от положения точек в бассейне, а на амплитуду установившихся колебаний координаты точек практически не влияют;

- при изменении периода от 1 до 12 ч и выше количество узловых линий на свободной поверхности уменьшается от трех до одной;

- увеличение периода волны вызывает уменьшение амплитуды возвышения свободной поверхности, уменьшение амплитуды волновой скорости установившихся колебаний и увеличение времени, за которое эти колебания

устанавливаются.

СПИСОК ЛИТЕРАТУРЫ

1. Алексеев Д.В., Черкесов Л.В. Влияние геометрических характеристик бассейна на параметры свободных и вынужденных волн // Допов. НАН Украши- 2001.-- № 3.™ С. 115-119.

2. Алексеев Д.В., Черкесов Л.В. Исследование влияния боковых стенок бассейна на характеристики вынужденных волн // Морской гидрофизический журнал.- 2001.- №» 4.~ С. 23-3 !,

3. Алексеев ДВ., Черкесов Л.В. Исследование влияния нелинейности на сейши в ограниченном бассейне // VII Междунар. науч.-техн. конф. «Современные методы и средства океанологических исследований»: Материалы конф.- М.: ИО РАН, 2001№ 3.- С. 52-53.

4. Коновалов А.В., Черкесов Л.В. Генерация длинных нелинейных волн в замкнутом бассейне движущимися возмущениями атмосферного давления // Изв. РАН. ФАО.- 1995,- 31, № 5.™ С. 713-718.

5. Черкесов Л.В., Иванов В.А., Хартиев СМ Введение в гидродинамику и теорию волн-С»Пб.: Гидрометеоиздат, 1992 - 264 с.

6. Arakawa A. Computational design for long-term numerical intégration of the équations of fluid motion: Two-dimensional incompressible flow // J. Compilation Physics - 1966.-1 - P. 119-143,

7. Роуч П. Вычислительная гидродинамика - M.: Мир, 1980.-616 с.

8. Марчук А.Н., Чубарое Л.Б., Шокин Ю.И. Численное моделирование волн цунами - Новосибирск: Наука, 1983.- 175 с.

9. Asselin R.A. Frequency filter for time intégrations. // Mont. Weat. Rev - 1980 - 100.- P. 487-490.

Морской гидрофизический институт HAH Украины, Материал поступил

Севастополь в редакцию 10.12.04

ABSTRACT Origin and development of tidal waves and flows in a basin with two straits are studied within the framework of the linear theory of long waves. The problem is solved using the numerical model that is based on the finite-difference scheme of the second order of precision over space and time. The field of velocities on the strait - basin boundary is assumed to be known. Analysis of the results of the numerical experiments permits to obtain the dependencies of the characteristics of tidal waves and flow structure upon the period of the forcing oscillations and the dissipation coefficient. It is shown, in particular, that the wave period increase entails the decrease of the steady-state oscillation amplitude and the increase of time necessary for settling the oscillations.

76

ISSN 0233-7584. Mop. гидрофиз. жури., 2006, № 2