Numerical simulation of tsunami type waves in a hydrodynamic channel.pdf Известно, что волны цунами являются одним из наиболее опасных и разрушительных бедствий, которым подвержено побережье Мирового океана. Вдали от береговой линии эти волны не представляют никакой опасности, так как их высота редко превышает 1 м. Однако длина волны в десятки раз больше глубины океана и, распространяясь в океане как в мелкой воде, волна цунами вовлекает в движение всю толщу воды от поверхности до дна. Поэтому в ней сосредоточена колоссальная энергия, которая со скоростью современного воздушного лайнера может распространяться на огромные расстояния. При входе волны цунами в зону мелководья скорость переднего фронта волны резко уменьшается, а амплитуда волны - увеличивается. В бухтах или устьях рек, из-за ограничения свободного пространства с боков, наблюдаются еще большие увеличения высоты волн, до 20 м и более. Причиной появления волн цунами являются подводные землетрясения, оползни, вулканы и др. практически не предсказуемые причины. Ясно, что проводить систематические научные исследования распространения волн цунами и их взаимодействия с различными препятствиями и сооружениями в натурных условиях не представляется возможным. Поэтому при исследованиях проблем цунами широко используют математические (аналитические и численные) и физические методы моделирования [1]. В работах [2-4] приведены технические характеристики Большого гидродинамического лотка ИПРИМ РАН, предназначенного для моделирования гравитационных волн в воде, даны методики измерений и первые результаты экспериментальных исследований. Показано, что распространение волн типа цунами в гидродинамическом лотке, как и в натурных условиях, с достаточной точностью описывается линейной теорией мелкой воды в невязком приближении. Однако при взаимодействии волн с различными подводными преградами и пузырьковыми завесами в ряде случаев необходимо учитывать нелинейные и вязкие эффекты. Настоящая работа посвящена созданию и апробации (путем сравнения с экспериментами) численного метода моделирования волновых процессов в гидродинамическом лотке в наиболее общей постановке: с учетом нестационарных, нелинейных и вязких эффектов для двухфазных сред. 1. Основные уравнения и численный метод Рассматривается нестационарная задача о течении вязкой несжимаемой жидкости со свободной поверхностью раздела в канале переменного сечения. Течение в лотке считается двумерным, т.е. геометрия лотка имеет бесконечную длину в z-направлении, ось x направлена вдоль лотка, а ось y - вертикально вверх. Числа Рейнольдса течений жидкости в гидродинамическом лотке могут достигать значений порядка 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] совместно с уравнением сохранения скалярной величины у: VU = 0; (1) —+ V(U y) = 0; (2) dt d(pU) + UV(pU) = -Vp + nV 2U + pg -pFa, (3) dt где p - давление, g = 9,81 м/с2 - ускорение силы тяжести, у - объемная концентрация несущей жидкости в расчетной ячейке. Значение скалярной функции у в ячейке может обозначать одно из трех состояний: у = 0 - ячейка содержит только воздух; у = 1 - ячейка содержит только воду; 0 < у < 1 - ячейка содержит границу раздела между жидкостью и газом. Таким образом, в нашем случае, у является индикатором межфазных поверхностей и свободной поверхности жидкости. Физические свойства среды рассчитываются как средневзвешенные величины в соответствии с объемными концентрациями фаз в каждой ячейке. Средневзвешенная плотность в ячейке рассчитывается как р = ур: + (1-у)р2, где р1 - плотность несущей жидкости, р2 - плотность воздуха; соответственно вязкость n = YHi + (1—y)n2. Fa - сила, обусловленная поверхностным натяжением: Fa =a-kVY, где с = 72,8 Н/м - коэффициент поверхностного натяжения вода -tvy^ воздух; к = V • -—Ц- I = Vn . ivyu Расчетная область, в соответствии с начальным состоянием уровней воды в лотке и в кессонном генераторе (см. рис. 1), включает в себя две подобласти: нижняя подобласть заполнена водой, имеющей заданную начальную конфигурацию поверхности раздела; верхняя - воздухом. В расчетной области могут устанавливаться различные (подвижные и неподвижные) преграды, макеты сооружений и пр. При t = о под действием силы тяжести начинается волновое движение, которое необходимо рассчитать. Граничные условия на жестких стенках канала (и стенках жестких преград) устанавливаются следующими: (4) n-Vy = о. v 7 Условия на верхней границе расчетной области устанавливаются в соответствии с технологией вычислительного алгоритма пакета OpenFOAM [7]. Открытый под лицензией GNU пакет OpenFOAM представляет собой набор написанных на С++ библиотек и макросов с открытым кодом и, по сути, является своеобразным языком программирования и удобным инструментом для численного моделирования физических процессов механики сплошных сред. Расчеты проводились на квадратных сетках с размером стороны 2 мм. Длина расчетной области составляла 15 м; высота: о,2 м и о,4 м. Аппроксимационная сходимость подтверждалась совпадением (x - О-диаграмм и профилей волны, полученных на исходной сетке и на сетке в четыре раза более грубой. Шаг по времени был нефиксированным и рассчитывался автоматически из условия, что число Куранта для уравнения переноса (2) не должно превосходить о,6. Результаты расчетов записывались через каждые о,1 с, а в случае необходимости - через о,о1 с. В качестве вычислительного комплекса использовалась 8-ядерная (2 процессора Intel Xeon 2,16 ГГЦ по 4 ядра) рабочая станция с 64-разрядной архитектурой. Распараллеливание вычислительного процесса проводилось разделением расчетной области поперёк направления распространения волны на 8 приблизительно равных (по количеству ячеек) частей. Для реализации параллельного вычисления использовался интерфейс OpenMPI. Для визуализации результатов вычисления использовалась программа ParaView. Оценка точности и достоверности численных результатов проводилась сопоставлением расчетов с экспериментальными данными, полученными в идентичных условиях. 2. Результаты численного моделирования и сравнение с экспериментом 2.1. Условия экспериментов и численного моделирования Большой гидродинамический лоток ИПРИМ РАН (БГЛ) имеет следующие габаритные размеры: длина - 15 м, ширина - о,26 м, высота - о,4 м. Эксперименты и численное моделирование волновых процессов в лотке проводились при следующих параметрах: - Начальная глубина воды в лотке (до откачки воздуха из генератора волны) Но изменялась от 1оо до 1о3 мм. - Длина волны X и 3 м, усредненная амплитуда падающих волн А в различных экспериментах составляла от 4,5 до 15 мм. Для регистрации волновых процессов в лотке используются резистивные датчики уровня воды, 10-канальная измерительная аппаратура [3], четырехканальный цифровой осциллограф и двухканальный регистратор Velleman PCS 500. Установка оснащена скоростной цифровой видеокамерой Photron FASTCAM SA4 500K, скорость съемки до 3600 кадр./с при полном разрешении 1024x1024 пс. Резистивными датчиками, которые располагались на различных расстояниях от генератора волны, измерялось смещение свободной поверхности воды в зависимости от времени £„(f). Это позволило построить волновые (x - ^-диаграммы для каждого эксперимента, определить скорости всех волн, а также, при взаимодействии с преградами, амплитудные коэффициенты отражения волн R = AR/A , где A и Ar - усредненные амплитуды падающей и отраженной волн соответственно. 2.2. Генерация волны В отличие от широко распространенного способа генерации волн при помощи различных подвижных механизмов, например вертикальным движением дна [1] или движением наклонной стенки генератора волны [8, 9], генератор волны в гидродинамическом лотке ИПРИМ РАН не имеет подвижных элементов. Способ генерации гравитационной волны в БГЛ основан на явлении распада произвольного разрыва уровней воды в лотке и в генераторе волны, который задается в начальный момент времени. Данный способ технически реализован в генераторе волны кессонного типа (рис. 1), который представляет собой отсек лотка (длиной а = 1465 мм) с герметичной верхней крышкой (1) и передней стенкой (2), которая в рабочем состоянии погружена в воду. В верхней крышке генератора волны расположен патрубок (3) для откачки или наполнения воздухом верхнего объема генератора. Нижняя часть объема генератора, высотой 90 мм, сообщается с рабочим объемом лотка. Перед работой лоток заполняется водой до уровня H0 > 90 мм. Затем, через патрубок (3), из верхней части генератора откачивается воздух, в результате чего устанавливается заданный перепад уровней воды п0: в генераторе (H + п0) и уровня воды в рабочей части лотка Н. После разгерметизации (при t = 0) верхней части объема генератора, в рабочей части лотка формируется волна длинной X и 2a и амплитудой А х п0/2. 3 H { f Воздух 1 П0 90 a Рис. 1. Схематический чертёж генератора волны кессонного типа Генератор волны, у которого в момент времени t = 0 мгновенно убирается передняя стенка, будем называть идеальным. Работа идеального и реального (кессонного) генераторов сравнивалась путем численного моделирования. На расстоянии x > a от передней стенки генератора, расчетные профили гравитационной волны £(t) идеального и кессонного генераторов практически совпадали друг с другом. На рис. 2 приведены результаты численного моделирования волнового процесса в БГЛ при образовании волны типа цунами идеальным генератором. При t = о, т.е. в момент мгновенного удаления передней стенки, дано начальное распределение уровней воды в генераторе волны (x < 1,465 м) и в лотке. Внизу (под белой горизонтальной линией) дан масштаб продольной скорости U. Видно, что при t = о скорости в генераторе волны и лотке равны нулю. При t = о,7 с видны две волны, которые распространяются в разные стороны от начального разрыва уровней воды. Вверху, над белой горизонтальной линией, приведены профили уровней воды в этих волнах, внизу - продольные скорости. Далее, при t = 2,о с наблюдаем отражение волны, которая распространялась влево, от задней стенки (x = о) генератора. В момент времени t = 3,3 с гравитационная волна, которая распространяется в сторону увеличения координаты x, почти полностью сформировалась. Скорость жидкости в лотке перед волной и за волной равна нулю. у, м |о,12 t = о с [ 0,10 Рис. 2. Численное моделирование процесса формирования волны типа цунами идеальным генератором. Величина параметра нелинейности в данном численном эксперименте составляет A/H = о,1 0,104 ilfflr -Щг -L ^ 1------ -f - A-AM/WA Рис. 3. Сравнение экспериментальной зависимости высоты волны от времени £ (t) (серая линия) с численным расчетом для идеального генератора (черная линия) на расстоянии 1,5 м от передней стенки генератора волны. На рис. 3 дано сравнение экспериментальной осциллограммы высоты волны от времени £ (t), измеренной датчиком уровня на расстоянии 1,5 м от передней стенки генератора волны, с расчетной зависимостью для идеального генератора (черная линия) при идентичных начальных условиях: H = 0,102 м; По = 0,015 м. Видно, что уже на расстоянии x = a s Х/2 экспериментальная зависимость практически совпадает с расчетной зависимостью для идеального генератора, далее (при x > a) совпадение еще лучше. 4(t), м ---1------- 0,112 -=--к-------- 0,110 0,108 0,106 0,102 0,100 0 2 4 6 8 t, c 0 2 4 6 8 t, 2.3. Распространение волн На рис. 4 представлена характерная (x - ^-диаграмма гравитационных волн типа цунами, распространяющихся в лотке. Внизу (под диаграммой) дан схематический чертёж лотка и места (1-4) расположения датчиков в том же масштабе по координате x. Параметры эксперимента: п0 = 15 мм, H = 102 мм. Видно, что скорости падающей и отраженных волн равны с = 1000 мм/с = const и соответствуют скорости, вычисленной по линейной теории мелкой воды с = yjgH для глубины лотка H = 0,102 м. Параметр нелинейности, вычисленный по средней высоте падающей волны, равен A/H = 0,074. Скорости распространяющихся в лотке волн (сплошная и пунктирная линии), полученные с использованием программы численного моделирования на основе полных уравнений Навье - Стокса, практически совпали с экспериментально измеренными величинами. Линии соответствуют результатам численного моделирования, точки - эксперименту. Черные линии - передний фронт волны, пунктирные - задняя кромка волны. Под диаграммой представлен схематический чертеж лотка: Г - генератор волны, Б - отражающий (береговой) склон, 1—4 - места расположения датчиков уровня волны. 40 а m Рис. 4. Диаграмма (x — t) гравитационных волн типа цунами, распространяющихся в гидродинамическом лотке при следующих начальных параметрах: п0 = 15 мм, Н = 102 мм, A/H = 0,074 На рис. 5 дано сравнение экспериментальной зависимости безразмерной скорости волны от параметра нелинейности А/Н [4] с численным моделированием на основе полных уравнений Навье - Стокса. Там же приведены теоретические зависимости: пунктирной линией - зависимость, рассчитанная на основе линейной теории мелкой воды, а штрихпунктирной линией - по нелинейной теории мелкой воды. Видно, что результаты численного моделирования (сплошная линия) лучше всего соответствуют экспериментальным данным. Нелинейная теория мелкой воды качественно правильно описывает процесс увеличения скорости волны с ростом параметра нелинейности, но дает завышенный результат. Линейная теория мелкой воды дает совпадение с экспериментом лишь при значениях A/H < 0,1. 2.4. Взаимодействие волн с непроницаемыми подводными преградами В данном разделе приведены результаты исследования взаимодействия волн с непроницаемыми подводными преградами двух типов (рис. 6) при малой величине параметра нелинейности A /Н < 0,1, т.е. исследовались волны малой амплитуды, когда распространение волн в канале постоянного сечения с достаточной точностью (см. рис. 5) описывается линейной теорией. Дано сравнение результатов численного моделирования с соответствующими экспериментальными исследованиями [4]. 30 20 10 0 ■ Исходная волна (фронт) □ Исходная волна (задняя кромка) ГТТ Численное моделирование Точки - эксперимент А Отраженная от берегового склона (фронт) _ Отраженная от берегового склона (задняя кромка волны) + Отраженная от начала генератора волна( фронт) 1 — 4 — Места расположения датчиков А/Н ■■■■■ Эксперимент (2011) о Эксперимент (2012) - - - 1D Нелинейная теория мелкой воды " " ' 1D Линейная теория мелкой воды —°— Численное моделирование (2D-Навье - Стокс) Рис. 5. Зависимость скорости длинных гравитационных волн в гидродинамическом лотке от параметра нелинейности А/Н Рис. 6. Схематический чертеж исследуемых преград № 1 и № 2 На рис. 7 приведены экспериментальные зависимости амплитудного коэффициента отражения R волн типа цунами от параметра преграды h/H (где h - затопленная высота преграды) при значениях параметра нелинейности A/H < 0,1 (светлые точки). Там же, для сравнения, даны результаты экспериментов [4, 9] при относительно больших значениях параметра нелинейности А/H > 0,15 (темные точки). Пунктирной линией, показан аналитический расчет по линейной одномерной теории [1]. R 0,8 0,6 0,4 0,2 0 0,2 0,4 0,6 0,8 1 h/H Рис. 7. Зависимость амплитудного коэффициента отражения волны типа цунами для подводных непроницаемых преград от параметра преграды h/H и параметра нелинейности A/H Видно, что линейная теория дает результаты, близкие к экспериментальным данным, лишь при малых значениях параметра преграды h/H. С другой стороны, численный расчет на основе полных уравнений Навье - Стокса дает хорошее совпадение с соответствующими экспериментами во всем диапазоне изменения параметра 0 < h/H < 1. Эксперименты при значениях параметра нелинейности A/H в диапазоне от 0,04 до 0,1 в координатах R = f(h/H) обобщаются единой кривой, коэффициент отражения волн не зависит от амплитуды падающей на преграду волны. Из рис. 7 также видно, что при значениях параметра нелинейности A/H > 0,15 (0,15 и 0,28) коэффициент отражения волн зависит от амплитуды волны и в этих случаях необходимо учитывать нелинейные эффекты. 2.5. Трансформация сильно нелинейной волны при взаимодействии с мелководьем Проблема описания динамики длинных гравитационных волн при распространении в мелкой воде прибрежной полосы с затоплением береговой зоны является одной из самых сложных проблем волн цунами. Это связано с необходимостью решения нестационарной задачи, а также с необходимостью учета нелинейных и вязких эффектов, т.е. в этом случае надо решать уравнения динамики в наиболее общей постановке. На рис. 8. дается сравнение результатов экспериментальных исследований трансформации сильно нелинейной волны типа цунами, при её распространении в прибрежной зоне с малым уклоном дна, с численным моделированием этих процессов в гидродинамическом лотке на основе полных уравнений Навье - Стокса. На рис. 8, а в координатах x -y (с соблюдением масштаба) показан профиль мелководной части дна (шельфа) и профиль волны в момент времени t = 8,3 с от начала генерации волны в БГЛ, полученный в результате численного моделирования. Пунктирной линией показан начальный уровень воды в лотке в момент генерации волны, при t = 0 с. А на рис. 8, с приведен соответствующий этому времени кадр, полученный в результате скоростной киносъемки процесса. Поле зрения скоростной цифровой камеры показано белым прямоугольником. На рис. 8, b и 8, d дано сравнение результатов численного расчета формы волны с экспериментом в момент времени t = 8,4 с. Из рис. 8 видно очень хорошее совпадение результатов численного моделирования с экспериментом. Рис. 8. Обрушение волны типа цунами при взаимодействии с пологим береговым склоном. Сравнение результатов численного моделирования при t, с: a - 8,3 и b - 8,4 с экспериментом - c -8,3 и d - 8,4 Выводы Результаты численного моделирования двумерных волновых процессов в гидродинамическом лотке ИПРИМ РАН (генерация и распространение волн типа цунами, а также их взаимодействие с непроницаемыми подводными преградами и пологим береговым склоном) на основе полных уравнений Навье - Стокса показали хорошее согласование численных расчетов с соответствующими экспериментальными данными. Таким образом, предложенная математическая модель и её программная реализация могут быть использованы для комплексного (экспериментального и численного) исследования проблем волн цунами в лабораторных условиях.

Скачать электронную версию публикации