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

2013

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика

№ 6(26)

МЕХАНИКА

УДК: 539.3, 539.4

Б.В. Бошенятов, Д.Г. Лисин

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВОЛН ТИПА ЦУНАМИ В ГИДРОДИНАМИЧЕСКОМ ЛОТКЕ

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

Ключевые слова: численное моделирование, уравнения Навье - Стокса, волны цунами, гидродинамический лоток, экспериментальные исследования.

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

В работах [2-4] приведены технические характеристики Большого гидродинамического лотка ИПРИМ РАН, предназначенного для моделирования гравитационных волн в воде, даны методики измерений и первые результаты экспериментальных исследований. Показано, что распространение волн типа цунами в гидродинамическом лотке, как и в натурных условиях, с достаточной точностью описывается линейной теорией мелкой воды в невязком приближении. Однако при взаимодействии волн с различными подводными преградами и пузырьковыми завесами в ряде случаев необходимо учитывать нелинейные и вязкие эффекты. Настоящая работа посвящена созданию и апробации (путем сравнения с экспери-

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

Рассматривается нестационарная задача о течении вязкой несжимаемой жидкости со свободной поверхностью раздела в канале переменного сечения. Течение в лотке считается двумерным, т.е. геометрия лотка имеет бесконечную длину в z-направлении, ось x направлена вдоль лотка, а ось у - вертикально вверх.

Числа Рейнольдса течений жидкости в гидродинамическом лотке могут достигать значений порядка 104, однако в наших расчетах модели турбулентности не используются. Обоснованием такого подхода служит тот факт, что эксперименты в каналах прямоугольного сечения [5] свидетельствуют о достаточно высоких числах Рейнольдса перехода в турбулентное состояние Re* = pUH/ц, где р, n, U -плотность, коэффициент динамической вязкости и скорость жидкости соответственно, Н - высота канала. При этом величина Re* увеличивается с уменьшением расстояния от входа в канал. Так, при x /Н = 60 начало перехода в турбулентное состояние соответствует значению Re1* = 8103, а конец перехода - Re2* = 1,8104. Кроме того, известно, что с уменьшением начальных возмущений в потоке число Рейнольдса перехода также увеличивается. В нашем случае (при длине волны X х 3 м и Н х 0,1 м) величина x /Н < 30, а начальные возмущения перед волной близки к нулю.

Таким образом, течение жидкости в лотке описывается уравнениями Навье -Стокса (1,3), которые решаются численным методом конечных объемов при использовании модели VOF (Volume of Fluid) [6] совместно с уравнением сохранения скалярной величины у:

трация несущей жидкости в расчетной ячейке. Значение скалярной функции у в ячейке может обозначать одно из трех состояний: у = 0 - ячейка содержит только воздух; у = 1 - ячейка содержит только воду; 0 < у < 1 - ячейка содержит границу раздела между жидкостью и газом. Таким образом, в нашем случае, у является индикатором межфазных поверхностей и свободной поверхности жидкости. Физические свойства среды рассчитываются как средневзвешенные величины в соответствии с объемными концентрациями фаз в каждой ячейке. Средневзвешенная плотность в ячейке рассчитывается как р = ур1 + (1-у)р2, где р1 - плотность несущей жидкости, р2 - плотность воздуха; соответственно вязкость п = 7Л1 + (1-у)П2. Ра - сила, обусловленная поверхностным натяжением: =а-кУу, где с = 72,8 Н/м - коэффициент поверхностного натяжения вода -

1. Основные уравнения и численный метод

VU = 0;

(1)

dY + V(U у) = 0;

dt

()

d (pU) + U V(pU) = -Vp + nV 2U + pg - pFCT,

dt

(3)

где p - давление, g = 9,81 м/с2 - ускорение силы тяжести, у - объемная концен-

Расчетная область, в соответствии с начальным состоянием уровней воды в лотке и в кессонном генераторе (см. рис. 1), включает в себя две подобласти: нижняя подобласть заполнена водой, имеющей заданную начальную конфигурацию поверхности раздела; верхняя - воздухом. В расчетной области могут устанавливаться различные (подвижные и неподвижные) преграды, макеты сооружений и пр. При t = 0 под действием силы тяжести начинается волновое движение, которое необходимо рассчитать.

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

Условия на верхней границе расчетной области устанавливаются в соответствии с технологией вычислительного алгоритма пакета OpenFOAM [7]. Открытый под лицензией GNU пакет OpenFOAM представляет собой набор написанных на C++ библиотек и макросов с открытым кодом и, по сути, является своеобразным языком программирования и удобным инструментом для численного моделирования физических процессов механики сплошных сред.

Расчеты проводились на квадратных сетках с размером стороны 2 мм. Длина расчетной области составляла 15 м; высота: 0,2 м и 0,4 м. Аппроксимационная сходимость подтверждалась совпадением (х - 0-диаграмм и профилей волны, полученных на исходной сетке и на сетке в четыре раза более грубой. Шаг по времени был нефиксированным и рассчитывался автоматически из условия, что число Куранта для уравнения переноса (2) не должно превосходить 0,6. Результаты расчетов записывались через каждые 0,1 с, а в случае необходимости - через 0,01 с.

В качестве вычислительного комплекса использовалась 8-ядерная (2 процессора Intel Xeon 2,16 ГГЦ по 4 ядра) рабочая станция с 64-разрядной архитектурой. Распараллеливание вычислительного процесса проводилось разделением расчетной области поперёк направления распространения волны на 8 приблизительно равных (по количеству ячеек) частей. Для реализации параллельного вычисления использовался интерфейс OpenMPI. Для визуализации результатов вычисления использовалась программа ParaView.

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

Большой гидродинамический лоток ИПРИМ РАН (БГЛ) имеет следующие габаритные размеры: длина - 15 м, ширина - 0,26 м, высота - 0,4 м. Эксперименты и численное моделирование волновых процессов в лотке проводились при следующих параметрах:

- Начальная глубина воды в лотке (до откачки воздуха из генератора волны) Н0 изменялась от 100 до 103 мм.

- Длина волны X и 3 м, усредненная амплитуда падающих волн А в различных экспериментах составляла от 4,5 до 15 мм.

U = 0; n ■ Vy = 0.

(4)

2. Результаты численного моделирования и сравнение с экспериментом

2.1. Условия экспериментов и численного моделирования

Для регистрации волновых процессов в лотке используются резистивные датчики уровня воды, 10-канальная измерительная аппаратура [3], четырехканальный цифровой осциллограф и двухканальный регистратор Velleman PCS 500. Установка оснащена скоростной цифровой видеокамерой Photron FASTCAM SA4 500K, скорость съемки до 3600 кадр./с при полном разрешении 1024x1024 пс.

Резистивными датчиками, которые располагались на различных расстояниях от генератора волны, измерялось смещение свободной поверхности воды в зависимости от времени t„(f). Это позволило построить волновые (х- ^-диаграммы для каждого эксперимента, определить скорости всех волн, а также, при взаимодействии с преградами, амплитудные коэффициенты отражения волн R = AR/A , где A и Ar - усредненные амплитуды падающей и отраженной волн соответственно.

2.2. Генерация волны

В отличие от широко распространенного способа генерации волн при помощи различных подвижных механизмов, например вертикальным движением дна [1] или движением наклонной стенки генератора волны [8, 9], генератор волны в гидродинамическом лотке ИПРИМ РАН не имеет подвижных элементов. Способ генерации гравитационной волны в БГЛ основан на явлении распада произвольного разрыва уровней воды в лотке и в генераторе волны, который задается в начальный момент времени. Данный способ технически реализован в генераторе волны кессонного типа (рис. 1), который представляет собой отсек лотка (длиной а = 1465 мм) с герметичной верхней крышкой (1) и передней стенкой (2), которая в рабочем состоянии погружена в воду. В верхней крышке генератора волны расположен патрубок (3) для откачки или наполнения воздухом верхнего объема генератора. Нижняя часть объема генератора, высотой 90 мм, сообщается с рабочим объемом лотка. Перед работой лоток заполняется водой до уровня Н0 > 90 мм. Затем, через патрубок (3), из верхней части генератора откачивается воздух, в результате чего устанавливается заданный перепад уровней воды п0: в генераторе (Н + п0) и уровня воды в рабочей части лотка Н. После разгерметизации (при ґ = 0) верхней части объема генератора, в рабочей части лотка формируется волна длинной X и 2а и амплитудой А х п0/2.

{ f Воздух

А L_J / /| / -

/ / /

/

П0

90 H

a

Рис. 1. Схематический чертёж генератора волны кессонного типа

Генератор волны, у которого в момент времени ґ = 0 мгновенно убирается передняя стенка, будем называть идеальным. Работа идеального и реального (кессонного) генераторов сравнивалась путем численного моделирования. На расстоянии х > а от передней стенки генератора, расчетные профили гравитационной

З

1

волны £(/) идеального и кессонного генераторов практически совпадали друг с другом.

На рис. 2 приведены результаты численного моделирования волнового процесса в БГЛ при образовании волны типа цунами идеальным генератором. При / = 0, т.е. в момент мгновенного удаления передней стенки, дано начальное распределение уровней воды в генераторе волны (х < 1,465 м) и в лотке. Внизу (под белой горизонтальной линией) дан масштаб продольной скорости и. Видно, что при / = 0 скорости в генераторе волны и лотке равны нулю. При / = 0,7 с видны две волны, которые распространяются в разные стороны от начального разрыва уровней воды. Вверху, над белой горизонтальной линией, приведены профили уровней воды в этих волнах, внизу - продольные скорости. Далее, при / = 2,0 с наблюдаем отражение волны, которая распространялась влево, от задней стенки (х = 0) генератора. В момент времени / = 3,3 с гравитационная волна, которая распространяется в сторону увеличения координаты х, почти полностью сформировалась. Скорость жидкости в лотке перед волной и за волной равна нулю.

у, м

10,12

- П 1 п ґ = 0 с

и, ли и Мадпіїисіе (т/е)

0 0.04 0.08 0.12

I

А , "~х, м

1,465

ґ = 0,7 с

Рис. 2. Численное моделирование процесса формирования волны типа цунами идеальным генератором. Величина параметра нелинейности в данном численном эксперименте составляет А/Н = 0,1

На рис. 3 дано сравнение экспериментальной осциллограммы высоты волны от времени ^ (О, измеренной датчиком уровня на расстоянии 1,5 м от передней стенки генератора волны, с расчетной зависимостью для идеального генератора (черная линия) при идентичных начальных условиях: Н = 0,102 м; По = 0,015 м. Видно, что уже на расстоянии х = а х Х/2 экспериментальная зависимость практически совпадает с расчетной зависимостью для идеального генератора, далее (при х > а) совпадение еще лучше.

Рис. 3. Сравнение экспериментальной зависимости высоты волны от времени Е, (t) (серая линия) с численным расчетом для идеального генератора (черная линия) на расстоянии 1,5 м от передней стенки генератора волны.

2.3. Распространение волн

На рис. 4 представлена характерная (х - 0-диаграмма гравитационных волн типа цунами, распространяющихся в лотке. Внизу (под диаграммой) дан схематический чертёж лотка и места (1-4) расположения датчиков в том же масштабе по координате х. Параметры эксперимента: По = 15 мм, H = 102 мм. Видно, что скорости падающей и отраженных волн равны с = 1000 мм/с = const и соответствуют скорости, вычисленной по линейной теории мелкой воды с = yjgH для глубины лотка H = 0,102 м. Параметр нелинейности, вычисленный по средней высоте падающей волны, равен A/H = 0,074. Скорости распространяющихся в лотке волн (сплошная и пунктирная линии), полученные с использованием программы численного моделирования на основе полных уравнений Навье - Стокса, практически совпали с экспериментально измеренными величинами.

Линии соответствуют результатам численного моделирования, точки - эксперименту. Черные линии - передний фронт волны, пунктирные - задняя кромка

волны. Под диаграммой представлен схематический чертеж лотка: Г - генератор волны, Б - отражающий (береговой) склон, 1—4 - места расположения датчиков уровня волны.

■ Исходная волна (фронт)

□ Исходная волна (задняя кромка) ГТТ Численное моделирование Точки - эксперимент

А Отраженная от берегового склона (фронт)

_ Отраженная от берегового склона (задняя кромка волны) ♦ Отраженная от начала генератора волна( фронт)

1 - 4 - Места расположения датчиков

Рис. 4. Диаграмма (х — Г) гравитационных волн типа цунами, распространяющихся в гидродинамическом лотке при следующих начальных параметрах: По = 15 мм, Н = 102 мм, А/Н = 0,074

На рис. 5 дано сравнение экспериментальной зависимости безразмерной скорости волны от параметра нелинейности А/Н [4] с численным моделированием на основе полных уравнений Навье - Стокса. Там же приведены теоретические зависимости: пунктирной линией - зависимость, рассчитанная на основе линейной теории мелкой воды, а штрихпунктирной линией - по нелинейной теории мелкой воды. Видно, что результаты численного моделирования (сплошная линия) лучше всего соответствуют экспериментальным данным. Нелинейная теория мелкой воды качественно правильно описывает процесс увеличения скорости волны с ростом параметра нелинейности, но дает завышенный результат. Линейная теория мелкой воды дает совпадение с экспериментом лишь при значениях А/Н < 0,1.

2.4. Взаимодействие волн с непроницаемыми подводными преградами

В данном разделе приведены результаты исследования взаимодействия волн с непроницаемыми подводными преградами двух типов (рис. 6) при малой величине параметра нелинейности А /Н < 0,1, т.е. исследовались волны малой амплиту-

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

А/Н

а Эксперимент (2011) о Эксперимент (2012)

- - - Ш Нелинейная теория мелкой воды

- ■' " Ш Линейная теория мелкой воды

—°— Численное моделирование (2Б-Навье - Стокс)

Рис. 5. Зависимость скорости длинных гравитационных волн в гидродинамическом лотке от параметра нелинейности А/Н

Рис. 6. Схематический чертеж исследуемых преград № 1 и № 2

На рис. 7 приведены экспериментальные зависимости амплитудного коэффициента отражения R волн типа цунами от параметра преграды h/H (где h - затопленная высота преграды) при значениях параметра нелинейности A/H < 0,1 (светлые точки). Там же, для сравнения, даны результаты экспериментов [4, 9] при относительно больших значениях параметра нелинейности А/Н > 0,15 (темные точки). Пунктирной линией, показан аналитический расчет по линейной одномерной теории [1].

R

0,8

0,6

0,4

0,2

0 0,2 0,4 0,6 0,8 1 ЫН

Рис. 7. Зависимость амплитудного коэффициента отражения волны типа цунами для подводных непроницаемых преград от параметра преграды Н/Н и параметра нелинейности А/Н

Видно, что линейная теория дает результаты, близкие к экспериментальным данным, лишь при малых значениях параметра преграды Н/Н. С другой стороны, численный расчет на основе полных уравнений Навье - Стокса дает хорошее совпадение с соответствующими экспериментами во всем диапазоне изменения параметра 0 < Н/Н < 1. Эксперименты при значениях параметра нелинейности А/Н в диапазоне от 0,04 до 0,1 в координатах R = _ДН/Н) обобщаются единой кривой, коэффициент отражения волн не зависит от амплитуды падающей на преграду волны.

Из рис. 7 также видно, что при значениях параметра нелинейности А/Н > 0,15 (0,15 и 0,28) коэффициент отражения волн зависит от амплитуды волны и в этих случаях необходимо учитывать нелинейные эффекты.

2.5. Трансформация сильно нелинейной волны при взаимодействии с мелководьем

Проблема описания динамики длинных гравитационных волн при распространении в мелкой воде прибрежной полосы с затоплением береговой зоны является одной из самых сложных проблем волн цунами. Это связано с необходимостью решения нестационарной задачи, а также с необходимостью учета нелинейных и вязких эффектов, т.е. в этом случае надо решать уравнения динамики в наиболее общей постановке. На рис. 8. дается сравнение результатов экспериментальных исследований трансформации сильно нелинейной волны типа цунами, при её распространении в прибрежной зоне с малым уклоном дна, с численным моделированием этих процессов в гидродинамическом лотке на основе полных уравнений Навье - Стокса. На рис. 8, а в координатах х - у (с соблюдением масштаба) показан профиль мелководной части дна (шельфа) и профиль волны в момент времени t = 8,3 с от начала генерации волны в БГЛ, полученный в результате численного моделирования. Пунктирной линией показан начальный уровень воды в лотке в момент генерации волны, при t = 0 с. А на рис. 8, с приведен соответствующий

этому времени кадр, полученный в результате скоростной киносъемки процесса. Поле зрения скоростной цифровой камеры показано белым прямоугольником. На рис. 8, Ь и 8, ё дано сравнение результатов численного расчета формы волны с экспериментом в момент времени t = 8,4 с. Из рис. 8 видно очень хорошее совпадение результатов численного моделирования с экспериментом.

Рис. 8. Обрушение волны типа цунами при взаимодействии с пологим береговым склоном. Сравнение результатов численного моделирования при ^ с: а - 8,3 и Ь - 8,4 с экспериментом - с -8,3 и ё - 8,4

Выводы

Результаты численного моделирования двумерных волновых процессов в гидродинамическом лотке ИПРИМ РАН (генерация и распространение волн типа цунами, а также их взаимодействие с непроницаемыми подводными преградами и пологим береговым склоном) на основе полных уравнений Навье - Стокса показали хорошее согласование численных расчетов с соответствующими экспериментальными данными. Таким образом, предложенная математическая модель и её программная реализация могут быть использованы для комплексного (экспериментального и численного) исследования проблем волн цунами в лабораторных условиях.

ЛИТЕРАТУРА

1. Левин Б.В., Носов М.А. Физика цунами и родственных явлений в океане. Научное издание. М.: Янус К, 2005. 360 с.

2. Бошенятов Б.В., Попов В.В., Левин Ю.К. и др. Моделирование волн типа цунами в лабораторных условиях // Сб. трудов Всероссийской конференции «Механика композиционных материалов и конструкций структурно-сложных и гетерогенных сред», ИПРИМ РАН. М.: Альянстрансатом, 2010. С. 51-60.

3. Бошенятов Б.В., Левин Ю.К., Попов В.В., Семянистый А.В. Метод измерения волн малой амплитуды на водной поверхности // ПТЭ. 2011. № 2. С. 116-117.

4. Бошенятов Б.В., Попов В.В. Экспериментальные исследования взаимодействия волн типа цунами с подводными преградами // Известия высших учебных заведений. Физика. 2012. Т. 55. № 9/3. С. 145-150.

5. Кутателадзе С.С., Миронов Б.П., Накоряков В.Е., Хабахпашева Е.М. Экспериментальные исследования пристенных турбулентных течений. Новосибирск: Наука, 1975. 166 с.

6. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free Boundaries // J. Comp. Phys. 1981. V. 39. Р. 201.

7. http://openfoam.org/

8. Shih-Chun Hsiao, Ting-Chieh Lin. Tsunami-like solitary waves impinging and overtopping an impermeable seawall: Experiment and RANS modeling // Coastal Engineering. 2010. V. 57. Is.1. Р. 1-18.

9. Фридман А.Н., Альперович Л.С., Шемер Л. и др. О подавлении волн цунами подводными барьерами // УФН. 2010. Т. 180. № 8. С. 843-850.

Статья поступила 15.09.2013 г.

Boshenyatov B.V., Lisin D.G. NUMERICAL SIMULATION OF TSUNAMI TYPE WAVES IN A HYDRODYNAMIC CHANNEL. Results of numerical simulation of generation and propagation of long gravitation waves, as well as their interaction with underwater obstacles in a hydrodynamic channel are presented. As a mathematical model, the two-dimensional Navier - Stokes equations in the incompressible liquid approximation are used. Results of the numerical simulation agree with experimental data.

Keywords: numerical simulation, Navier - Stokes equations, tsunami waves, hydrodynamic channel, experimental investigations.

BOSHENYATOVBorisVladimirovich (Institute of Applied Mechanics RAS)

E-mail: bosbosh@mail.ru

LISIN Dmitry Gennadievich (Institute of Applied Mechanics RAS)

E-mail: 4lsdima@gmail.com