CC BY
11
Вычислительные технологии
Том 24, № 2, 2019
Нестационарные колебания ледового покрова
в замороженном канале под действием движущегося
>к
внешнего давления*
К. А. ШИШМАРЕВ1'*, Т. И. ХАБАХПАШЕВА2 1Алтайский государственный университет, Барнаул, Россия
2Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия ^Контактный e-mail: shishmarev.k@mail.ru
Рассмотрена нестационарная задача об изгибно-гравитационных волнах на поверхности замороженного ледового канала, вызванных движением внешней нагрузки. Задача решается методом преобразования Фурье вдоль канала и разложением прогибов льда по модам бегущих изгибно-гравитационных волн в канале.
Предложен алгоритм, позволяющий привести решение, полученное с помощью обратного преобразования Фурье, к аналитическим формулам в виде суммы стационарных и нестационарных прогибов при конечных временах в системе координат, движущейся вместе с нагрузкой. Исследована эволюция прогибов во времени для разных скоростей движения нагрузки. Показано, что при увеличении времени решение для нестационарных прогибов ледовой пластины выходит на установившийся режим в системе координат, связанной с нагрузкой. Исследовано влияние коэффициента затухания на образование волн.
Ключевые слова: ледовая пластина, канал, гидроупругая волна, внешняя нагрузка.
Библиографическая ссылка: Шишмарев К.А., Хабахпашева Т.И. Нестационарные колебания ледового покрова в замороженном канале под действием движущегося внешнего давления // Вычислительные технологии. 2019. Т. 24, № 2. С. 111-128. DOI: 10.25743/ICT.2019.24.2.010.
Введение
Гидроупругие волны, вызванные внешними нагрузками, хорошо изучены для неограниченного ледового покрова (см., например, обзор Squire и др. [1]). Известно, что прогибы льда существенно зависят от скорости движения нагрузки. Для неограниченного ледового покрова, в отличие от волн на поверхности жидкости, имеются две критические скорости, обусловленные существованием локального минимума на дисперсионной кривой [2]. При движении с докритической скоростью прогиб льда локализован под нагрузкой, а при сверхкритической скорости образуется система волн, уходящих от тела. Оба эти случая хорошо описываются в рамках линейной теории. Критический случай в линейной теории приводит к бесконечным прогибам, что нефизично, поэтому для
* Title translation and abstract in English can be found on page 128.
© ИВТ СО РАН, 2019.
моделирования применяются нелинейные модели [3], модели с демпфированием и затуханием [2]. Кроме того, на практике это может означать разлом льда, что потребует переформулировку задачи.
Вязкоупругие свойства льда исследовал Табата в ходе экспериментов с ледяным прямоугольным брусом, подвергавшимся напряжениям [4]. На основе анализа полученных данных Табата заключил, что для описания свойств льда лучше всего подходит реологическая модель Кельвина — Фойгта с одним параметром — так называемым временем затухания. Эта модель использована Козиным [5] для сравнения теоретических и экспериментальных результатов и признана удовлетворительной [5, 6]. Однако, чтобы получить удовлетворительное совпадение с экспериментальными данными, время затухания варьировалось от 3 до 10 с.
Аналогичные задачи, но для ледового покрова в каналах изучены менее подробно, хотя наличие стенок и условия закрепления на них льда существенно влияют на распределение прогибов и напряжений. Эти задачи имеют большое прикладное значение, поскольку лабораторные эксперименты с движением нагрузок по ледовому покрову проводятся именно в каналах. Изучение этих задач обусловлено также необходимостью разрушения ледового покрова в замороженных северных реках ранней весной для предотвращения наводнений, а также для безопасной транспортировки грузов вдоль замороженных рек и каналов [7, 8].
Наличие стенок и ограниченность ледовых пластин приводят к появлению новых граничных условий и существенному усложнению задачи. При этом меняются все характеристики гидроупругих волн, а также распределение прогибов и напряжений в плавающей пластине (см., например [9-12]).
В случае изгибно-гравитационных волн в канале граничные условия на его стенках приводят к образованию счетного множества дисперсионных кривых и соответствующих критических скоростей [13-15]. Каждая точка дисперсионной кривой соответствует некоторой линейной волне, распространяющейся вдоль канала с фиксированным профилем поперек канала. Анализ гидроупругой реакции ледового покрова канала на движущуюся нагрузку проведен в [8] для решения в виде бегущей волны согласно модели вязкоупругого материала Кельвина — Фойгта. В рамках линейной модели без затухания асимптотическое решение задачи исследовано в [16, 17], где показано, что нагрузка, движущаяся вдоль канала с заданной скоростью при больших временах, может порождать только волны, которые распространяются с такой же скоростью, а их количество всегда зависит от скорости движения нагрузки по отношению к критическим скоростям и их всегда конечное число. Результаты аналитических исследований асимптотической модели хорошо известны, к тому же они не требуют больших численных расчетов. Ожидается, что полученные для этой модели прогибы и удлинения являются оценкой сверху возможных прогибов и удлинений для моделей с демпфированием и затуханием. Более детальный анализ асимптотического решения для параметров ледового бассейна проведен в Приамурском государственном университете [17], где получено совпадение первой критической скорости со скоростью движения нагрузки, приводящей к разрушению модельного льда. Такой метод назван в [7] резонансным методом разрушения ледяного покрова. Аналогичные оценки для неограниченного льда получены в [18].
В настоящей работе рассмотрена нестационарная задача о движении внешней нагрузки из состояния покоя. Внешняя нагрузка моделируется как локализованное пятно давления, движущееся с постоянной скоростью вдоль центральной линии канала. Задача решается методом преобразования Фурье вдоль канала и с помощью разложе-
ния образа Фурье профиля волны поперек канала на суперпозицию собственных мод колебаний ледового покрова в канале. Расчеты проведены для параметров ледового бассейна, изучавшегося в [8, 14]. Основными параметрами численного решения, влияющими на результат, являются количество собственных мод, область и шаг численного интегрирования решения, полученного в квадратурах. Построено зависящее от времени решение в виде суммы стационарных и нестационарных прогибов, затухающих в отдалении от нагрузки при конечных временах. Для построенного решения исследована эволюция прогибов во времени, показано, как решение нестационарной задачи выходит на асимптотическое решение при I ^ ж. Кроме того, исследовано влияние изменения коэффициента затухания на образование волн.
1. Постановка задачи
Рассмотрим нестационарные колебания ледового покрова в канале прямоугольного сечения. Колебания вызваны движением внешней нагрузки вдоль канала по поверхности льда. Параметры канала и введенная декартова система координат (х, у, г) показаны на рис. 1. Канал наполнен невязкой несжимаемой жидкостью с плотностью р1, движение жидкости в канале считается безвихревым. Ледовый покров моделируется тонкой вязкоупругой пластиной, примороженной к стенкам канала, у = ±Ь. Толщина пластины Н постоянна, ее жесткость Б = ЕЩ/[12(1 — и2)], где Е — модуль Юнга, и — коэффициент Пуассона. Нестационарная задача о колебаниях ледового покрова в канале рассматривается в рамках линейной теории гидроупругости [1]. Прогиб ледового покрова ы(х,у,1) описывается уравнением колебаний вязкоупругой пластины [19]
где У2 = У2 • У2 = д4/дх4 + 2д4/(дх2ду2) + д4/дх4; М = р{Н — масса ледово-
го покрова на единицу площади; — плотность льда; т — коэффициент затухания; р(х,у, 0,£) — давление жидкости, оказываемое на нижнюю поверхность ледовой пластины; Р(х, у, ¿) — давление, вызванное движением внешней нагрузки, действующее на верхнюю поверхность ледового покрова; I — время.
(1)
—ж < х < ж, —Ь < у < Ь, г = 0,
Я
Рис. 1. Схематическое изображение канала и движущегося пятна давления
Давление жидкости р(х,у, 0,£) на границе лед — жидкость определяется из линеаризованного интеграла Коши — Лагранжа
р(х,у, 0,Ь) = -р1 <рг(х,у, 0, ¿) - р1дт(х,у,г), (2)
— то < х < то, —Ь < у < Ь, г = 0,
где д — ускорение свободного падения, р\ — плотность жидкости и (р(х,у, г,£) — потенциал скорости течения жидкости. Потенциал <р(х,у, удовлетворяет уравнению Лапласа в области течения жидкости
У2^, у, г, ¿) = 0, —то < х < то, —Ь<у<Ь, —Н < г < 0, (3)
линеаризованному кинематическому условию на границе лед — жидкость, граничным условиям непротекания на стенках и дне канала
V* = (г = 0), уу = 0(у = ±Ь), ^ = 0(г = —Н), (4)
условию затухания в отдалении от нагрузки для конечных времен
р — 0, |ж| — то, Ь < то. (5)
В начальный момент времени течение в канале отсутствует и потенциал <р(х,у, удовлетворяет начальным условиям
ф = 0, ^ = 0(* = 0). (6)
Прогиб ледового покрова т(х,у,1) удовлетворяет следующим краевым и начальным условиям: условие примороженного ледового покрова к стенкам канала моделируется условиями жесткого закрепления
т = 0, = 0, —то < х < то, у = ±Ь, (7)
колебания ледового покрова затухают в отдалении от нагрузки при конечных временах
т — 0, |ж| —У то>, Ь < то, (8)
в начальный момент времени прогиб льда удовлетворяет условиям
т(х,у, 0) = ■то(х,у), т1(х,у, 0) = 0. (9)
Условие для wt следует из кинематического условия (4) и предположения, что в начальный момент времени течение в канале отсутствует. Начальные условия для прогиба ледового покрова определяются из уравнения пластины (1) при £ = 0
В^\-ы0(х,у) = —р1уш0(х,у) — Р(х,у, 0). (10)
Внешнее давление Р(х, у, ¿) = Р(х — иь,у) действует на верхнюю поверхность ледового покрова и движется вдоль центральной линии канала с постоянной скоростью и. В подвижной системе координат (X, у^), X = х — внешняя нагрузка не изменяется во времени и имеет следующую форму:
Р(х — 1Л,у) = РоРг (Х/Ь) Р2(—то < х < то, —Ь<у<Ь, (11)
Рг(Х/1) = (cos(ndX/L)+ 1)/2 (a\X \ < L), P±(X/L) = 0 (Cl\X \> L),
P2(y/L) = (cosfrc*y/L) + 1)/2 (c2\y\ < L), P2(y/L) = 0 (c2\y\ > L),
где P0 — амплитуда внешнего давления; C\ и c2 — безразмерные параметры внешней нагрузки, характеризующие размер области внешнего давления. Такое "пятно" давления может быть вызвано, например, судном на воздушной подушке [5].
Решение системы (1)—(11) зависит от размеров канала H,L, параметров жидкости и льда pi, pi, hi, D, коэффициента затухания т, параметров внешней нагрузки U, с\, с2 и Р0.
Решение задачи (1)—(11), где Р(X, у) = 0 и начальные условия и условия при \ж\ ^то отброшены, изучено в [14] в виде линейных синусоидальных волн, распространяющихся вдоль канала. Дисперсионные соотношения между частотой волны шп и волновым числом к определены численно совместно с профилями волн поперек канала. Показано, что существует бесконечное число подобных волн в замороженном канале, п > 1. Эти волны являются комбинацией волн плескания поперек канала и изгибно-грави-тационных волн, распространяющихся вдоль него. В системе координат, движущейся вдоль канала с фазовой скоростью n-й волны, сп(к) = шп(к)/к, эти волны стационарные и периодические. В используемом в настоящей статье методе применяются свойства этих волн.
Рассматривается нестационарная задача с учетом нулевого коэффициента затухания (т = 0) в уравнении (1). Решением задачи является функция w(x, у, t), вычисленная при некоторых типичных параметрах задачи и конечных временах (t < то). Асимптотическое решение системы (1)—(11) получено в [17]. Это решение позволяет получить простую интерпретацию результатов и разложить прогибы льда w(x,y,t), t ^ то, на составные части: затухающий в отдалении от внешней нагрузки прогиб льда, симметричный в системе координат (Х,у); серия из конечного числа гидроупругих волн, распространяющихся от нагрузки вдоль канала, фазовая скорость которых равна скорости нагрузки, сп(к) = U.
В зависимости от скорости нагрузки число таких волн может быть разное, но конечное. Полученные аналитические формулы показывают, что короткие волны имеют групповую скорость сП(к) = dun/dk, превышающую скорость нагрузки U, а следовательно, распространяются перед нагрузкой. Соответственно длинные волны имеют групповую скорость, меньшую скорости нагрузки, и распространяются за нагрузкой. Ожидается, что решение задачи (1)—(11) при конечных временах будет обладать аналогичными свойствами формирования профилей колебаний и прогибы w(x,y,t) в системе координат (X, у) будут стремиться к асимптотическому решению при увеличении времени (алгоритм его получения описан ниже). Представленный в настоящей работе анализ для прямоугольного канала и внешней нагрузки в форме (11) может быть перенесен на непрямоугольные каналы и другие формы внешней нагрузки.
Решение задачи (1)—(11) в виде бегущей волны, w = w(X,y), ip(X,y,z), при ненулевом коэффициенте затухания (г > 0) исследовано в [8]. Использованный метод решения не требовал задания начальных условий (6), (9), (10). Задача решена численно. Показано, что прогибы льда и удлинения в ледовом покрове быстро затухают с расстоянием от нагрузки в рамках модели вязкоупругой пластины. В большинстве вычислений в [8] использовано значение т = 0.1 с. Ниже показано, что полученное решение нестационарной задачи (1)—(11) при конечных временах с ненулевым коэффициентом затухания (г > 0) при увеличении времени будет стремиться к стационарным затухающим про-
гибам, аналогично рассмотренным в [8], а при уменьшении самого коэффициента затухания — к асимптотическому решению при больших временах. Исследование самого асимптотического решения и его свойств проведено в [16, 17].
2. Метод решения
Начально-краевая задача (1)—(11) решается с помощью преобразования Фурье вдоль канала
те те
т(Х, у,1) = —= У у, ту(С, у,1) = —= ^ т(х, у, ¿х, (12)
—те —те
где £ — параметр преобразования. Уравнение (1) после преобразования Фурье принимает вид
МьУ + И {Шруууу - + ) = - ридь? - Рр& у)е—гт, (13)
где (£, у, г, Ь) и Ру(^, у)е—г^и* — соответствующие образы Фурье потенциала ^(х, у, г, £) и внешнего давления РУР1(Х)Р2(у). Образ Фурье функции прогибов льда и^(^,у,Ь) удовлетворяет начально-краевым условиям
™У = ™УУ(С,У), ™У = 0 (г = 0), ту = 0, тУУ = 0 (у = ±Ь), (14)
где (С,у) — образ Фурье начального условия для прогиба 1УУ(х,у).
Потенциал <ру(£, у, г, Ь) удовлетворяет модифицированному уравнению Гельмгольца
'АУу + & = (-Ь<у<Ь, -Н<г< 0), (15)
и начальным и граничным условиям
^ = Шуу (г = 0), VУУ = 0 (у = ±Ь), ^ = 0 (г = -Н), <ру = ^ = 0(1 = 0). (16)
Задача (13)—(16) решается методом разделения переменных £ и у для образа Фурье прогибов льда (£,у,£):
те
*>У (£,УЛ = £ ап(С,г)МС,у). (17)
П=1
Здесь фп(£,у) — собственные моды колебаний ледового покрова в канале: ап(^,1) — координаты разложения на собственные моды, определяемые в ходе решения. Функции фп(£,у) вычислены в [14]. Они описывают гидроупругие волны, распространяющиеся в канале
Ш(х, у, г) = Ке [Апфп(£, у)ет, (18)
где шп(^) — частота волны, £ > 0 — волновое число и Ап — амплитуда волны. Подставляя (17) в кинематическое условие (16) на границе лед — жидкость (г = 0), получим уравнение, из которого следует, что потенциал <ру (£, у,г,Ь) может быть записан в виде
те
4>у (С, У, г,1) = ^ ап,г(С, Фп(£, У, г). (19)
п=1
Функции Фп(^,у,г) удовлетворяют уравнению (15) и граничным условиям (16), где кинематическое условие имеет вид (£,у, 0) = фп(£,у).
Подставляя (17) и (19) в уравнение (13), приходим к дифференциальному уравнению для определения ап(^,1):
те
^ ап,и [Мф,п + Р1 Фп(С, У, 0)] +
п=1
те
+ ^ ап [Я {Фщуууу - 2фп,уу + £4фп) + мфп] = -Рр(С, у)е-*ш. (20)
п=1
Заметим, что выражения в квадратных скобках в (20) пропорциональны друг другу [14]
^ (£) [Мфп + р1&п(£, У, 0)] = [В (фщуууу - % Чп,уу + ^Фп) + Р19фп] . (21)
Здесь коэффициент пропорциональности равен квадрату частоты ш2п (£). Собственные колебания фп(^,у), п > 1, являются ортогональными в следующем смысле [14]:
ь ь
J [Мфп + Р1 Фп(С, У, 0)] 'Фш&У = $пт&п(0, *п(€) = ! [Мфп + Р1 Фп(С, У, 0)] Фп^.у. (22) - ь - ь
Функции фп(£,у) нормированы условием [14]
тах фп(^,у) = 1. (23)
-Ь<у<Ь
Данная нормировка предполагает, что функции фп(у) являются безразмерными. С учетом (22), умножая уравнение (20) на фп(£,у) и интегрируя результат по у от — Ь до Ь, приходим к системе обыкновенных дифференциальных уравнений для коэффициентов
ь
Н2п 1 Г
+ ^ (Оъ = Нп(0е-*т, Нп(0 = — Рр (и)Фп(иУу, п =1, 2,... (24)
п - ь
Начальные условия для уравнений (24) следуют из (14)
ап & 0) = , ^ & 0) = 0. (25)
< (4) аЪ
Решение задачи (24), (25) имеет вид
(С А = $иНп($) ( ^ _ \ + (26)
ап(^]= 2шп(0 и + ^ ип — ^и) + щп — т2. (26)
Слагаемые в (26) могут быть сингулярными. Особенности наблюдаются в точках £, являющихся решением уравнения шп(£) — £и = 0, которое эквивалентно уравнению сп(^) — и = 0. Следовательно, сингулярности образуются, если скорость нагрузки и достаточно большая (и > с^щ). Заметим, что минимальные фазовые скорости соответствующих собственных мод с^щ являются упорядоченными по возрастанию,
с^ < с1т < ... Если скорость и меньше с^щ, то особенности в (26) отсутствуют. В ином случае для определения точек, в которых будут особенности, и для определения количества слагаемых с особенностями в разложении (17) необходимо вычислить фазовые скорости сп(£). Эти скорости будут зависеть от параметров канала и льда.
Пример вычисления скоростей для первых трех четных по у мод при р\ = 1024 кг/м3, рг = 917 кг/м3, модуле Юнга Е = 4.2 • 109 Н/м2, коэффициенте Пуассона у = 0.3, толщине ледового покрова 10 см, глубине канала 2 м и ширине 20 м представлен на рис. 2 (более детальное описание дисперсионных соотношений, фазовых и групповых скоростей дано в [14]). Выбор скорости движения и показан горизонтальными линиями. При и = и0 особенности отсутствуют и прогиб пластины наблюдается только в окрестности пятна давления. При с^щ < Щ < с^щ дополнительно появляется одна система коротких волн, соответствующая пересечению с кривой с1. При с^щ < и2 < с^ появляются еще одна система коротких волн и система длинных волн, соответствующая точкам пересечения с кривой с2, и т.д. Отметим, что для этих значений параметров вычисленное значение с^щ совпадает с с1 (0). В общем случае с1 (0) > с^щ, и тогда при с^щ < и < с1 (0) будут порождаться две системы волн, одна за нагрузкой, а другая — перед ней. Вычисления фазовых скоростей показывают, что слагаемые в (26) являются сингулярными
I /-(1) I Лп) I /-(п) ТТ ^ (п)
в точках ±£1 для первой моды и в точках ±£1 и ±£2 , если и > с^П, для всех п > 1. Дополнительно, а1 (£,£) может быть сингулярной в точке £ = 0, потому что (С) = ) при £ ^ 0. Однако легко показать, что а1 (£) = 0(£-2), и тогда Н1 (£) = 0(£2) при £ ^ 0 в уравнении (24). Следовательно, ап(£,£) имеет конечный предел при £ ^ 0 для всех п.
В движущейся системе координат (X = х — Ш) прогибы льда (м(Х, у,Ь)) определяются как обратное преобразование Фурье, примененное к уравнению (19): 1 ~ 7
Ь)(Х,У,1) = ™п(Х,У^), Ып (Х,у,г)= е*ша,п (^Ж а,у)е*х (27)
V п=1 ]
-7
Уравнения (15) и (21) показывают, что гфп(£, у), шп(£) и ап(^) являются функциями от £2, а следовательно, это четные функции параметра £. Функция Нп(£) в (20) тоже является четной функцией от £ для рассматриваемой нагрузки (11). Интегралы в (24) могут быть сведены к интегралам от 0 до то при использовании свойств подынтегральных функций
с, II, м/с
0 0.5 1 1.5 к М"1
Рис. 2. Фазовые скорости первых трех четных мод
относительно £. Тогда функция тп(Х,у,1) в разложении (27) примет вид
Нп(ОМ£,у)
2оов(£Х) £ие*х
+
^ (£) — 2ип (^)
х
X
(
Шп(£) + £и Шп(£) — £и
+ е
-2г £Х
е-г (ооп(0+И и )г
"п(£) + Си
— е
-2г £ Х
ег(ооп(0-а и )г
^п(0 — си
%. (28)
Первое, третье и пятое слагаемые в (28) могут быть сингулярными в зависимости от выбора скорости и и номера моды п, однако эти особенности компенсируют друг друга. Это позволяет вычислить интегралы тп(Х, у, Ь) численно для конечных £ и определить прогибы льда т(Х,у,1) в заданный момент времени.
Для численного решения число мод фп(£,у) в разложении (17) ограничивалось конечным числом Мшо(1
^mod 1 ^mod
Ыр(£,У,^) = V aп(£,t)^фп(£,У), = —= V тп(Х,у,Ь).
, у 2п ,
п=1 у п=1
(29)
Заметим, что фп(£,у) нормированы условием (23), а из уравнения (26) следует, что коэффициенты ап (£, Ь) для рассматриваемой нагрузки в форме (11) стремятся к нулю при п ^ то и фиксированных ии I. Для всех конечных X и £ подынтегральные функции в (28) затухают с ростом £ и ведут себя как 0(1/(и^а) — £2и2)) при £ ^ то, где )/£ ^ то при £ ^ то для всех п. Это гарантирует быстрое затухание и существование обратного преобразования Фурье. Интегралы (28) приближаются соответствующими интегралами с конечной областью интегрирования от 0 до ^. Отрезок интегрирования [0, ^] разбивается на малые отрезки с шагом Д£ и большим количеством шагов N. На каждом шаге (£ Е [^ ,£7+1]) необходимо вычислять следующий интеграл:
£з+1
Яп(С,у)
2ооз(£Х) (ооъ^Х + ш£Ч] ео8[£Х — и^Ч]
+
^па) — еи2 и2п&)
<
Шп
(С )
(Шп(0 — £и).
%,
(30)
где Яп(^,у) = Hп(i,y)■фп(í,y), ^ = (Шп(£) + ) и Ш(п}
В случае и < с^щ этот интеграл (28) не имеет слагаемых с особенностями. Тогда для вычисления тп(Х,у,1) из (28) этот интеграл разбивается на сумму пяти интегралов. Функции, зависящие от £ и входящие в последние интегралы, приближаются линейными функциями на отрезке [^ ,£7+^, и полученные интегралы вычисляются аналитически с использованием интегрирования по частям. Например, второе слагаемое в интеграле (30) дает два интеграла, один из которых имеет вид
£¿+1
оо8(В' г) (А 7 (у)+ ^ (у))оо8(В3 (Х,г)№,
(31)
0
(С \ Шп (0+1) — Шп( 0) с и (у у , тт. , Шп( 0+1) — Шп( 0)
^ = ип( О)--д-&, У](X, ч = х + иг +--д-.
у — ^^ — Я'] о,у)
А (У) = а - о, А (у)= 3+\, ], & (у)
дс ^ де ' ¿п(оЖ(о)+&иу
Интегралы от остальных слагаемых вычисляются аналогично. Решенное таким образом уравнение (28) приводит к аналитическим формулам для вычисления тп(Х, у, ¿). Заметим, что полученные формулы состоят из стационарной части, не зависящей от времени и симметричной по X, полученной из первого слагаемого в (30), и нестационарной части. Стационарная часть ведет себя как 0(1/X) при X ^ то, а нестационарная часть — как 0(1Д) при Ь ^ то и и < сШ1п соответственно.
В случае и > сШ1п интеграл (28) имеет слагаемые с особенностями. Интеграл с регулярным слагаемым, имеющим еов[£X + ¿] в числителе, определяется по методу, описанному выше для и < сЩШ1п. Преобразуем оставшийся интеграл с особенностями.
п
-"Ш1П
Из (28) имеем
оо
р = Щ (Р л( 2со^(ех) си оо8]ех—М^) —ШЖ\ ИР (32) I =] ^, У)[-пщу—ёи2 — Ш — Си ; ^ (32)
о
все слагаемые в этом интеграле имеют особенности в двух точках £(га) и С^^ для всех мод, кроме первой (п > 1). Для первой моды особенность будет в одной точке. Выделяя симметричные области с особенностями, представим интеграл (28) в виде суммы пяти интегралов
?1п)+й1 2 ?2п)+<ь ^
г = у ...<£+/ + / + / ...<%. (33)
0 ?(П)—1 ?(П)+Й1 ?(П)-Й2 ^+¿2
Параметры, характеризующие интервалы разбиения, зависят от времени ( 1 = 1( ) и 2 = 2( ) ) и выбираются таким образом, чтобы выполнялись неравенства
|Ы£) — £ и Ж<< 1 (34)
на интервалах (^ — ^+д1) и (С2 — д2, С2+62). Первый, третий и пятый интегралы в (33) не имеют особенностей и вычисляются численно по описанному выше алгоритму. Для вычисления второго и четвертого интегралов подынтегральное выражение преобразуем с учетом (34) следующим образом:
& ) ( 2оо8(еX) ^ ео8[еX — (ип(С) — ^Ц)*] N = ^, у)\ш1 (о—еи2 шп (о шп(о — си )
Яп(С, у)
+ ^ оо8(е X) — ^ ^(сх — Шп(^2 — ^ *)
(С)ЫС) +си) ^ ' (0 V 2
(35)
Вычисление интегралов с использованием интегрирования по частям от функции (35) по интервалу (^ — 8^, ^ + 8г), г = 1, 2, приводит к аналитическим формулам, которые в сумме с аналитическими формулами при стационарной части и затухающей нестационарной части решения дают формулу для вычисления прогибов льда ,ш(Х, у, Ь) в заданный момент времени. Заметим, что интегралы от функции (35) не затухают с ростом ¿, а при £ ^ то приближаются к установившимся волнам, параметры которых не зависят от времени в движущейся системе координат (Х,у). Цель численного анализа — вычисление прогибов (29) по предложенному численному алгоритму и полученным аналитическим формулам для т(Х,у,1) при разных временах.
Для исследования эволюции прогибов льда при увеличении времени дополнительной задачей является определение поведения решения (28) рассматриваемой нестационарной задачи при больших временах. Интеграл (28) регулярный и может быть записан в виде суммы пяти интегралов в смысле главного значения. Каждый интеграл в отдельности исследуется асимптотическими методами. Если интеграл не имеет особых точек, то его вклад в формирование нестационарных прогибов имеет порядок 0(1/1) при £ ^ то. Вклад нерегулярных интегралов при больших временах состоит из двух частей: первая часть — стационарный прогиб, локализованный в области нагрузки, четный по X и затухающий при |Х| ^ то; вторая часть — конечная система волн (18) с амплитудами Ап, распространяющихся перед и за нагрузкой. Эти волны имеют волновое число к = £ЩШ, т = 1, 2, где ^Ш являются решениями уравнения шп(£) = £и, £п < Сп, и амплитуды АШ вычисляются по формуле
А™ 2^Нп(СШ) , (36)
ш
Сши (сп^ш) — и)
где ¿п — групповая скорость п-й собственной моды колебаний ледового покрова в канале. Фазовые скорости этих волн равны скорости движения нагрузки и [16]. Длинные волны с волновым числом {Ц распространяются за нагрузкой, а короткие волны с волновым числом £п — перед ней [8].
3. Численные результаты
Задача о гидроупругих волнах в канале решается численно для канала, имеющего глубину Н, равную 2 м. Ширина канала 2Ь — 20 м, толщина ледового покрова принималась равной 10 см, плотность льда рг = 917 кг/м3, модуль Юнга Е = 4.2 • 109 Н/м2, коэффициент Пуассона V = 0.3. Нагрузка (11) действует на область 2 х 2 ми движется вдоль центральной линии канала с постоянной скоростью и.
Основными параметрами вычислений являются параметры интегрирования по £, длина интервала ^ и количество шагов, а также количество мод Хто(1 в разложении функции и)р(^,у,Ь). На выбор параметров Х^ и количества шагов влияют несколько факторов. Первый — образ Фурье функции внешней нагрузки Рр(£, у). Преобразование Фурье внешней нагрузки в форме (11) затухает при £ ^ то. Скорость затухания зависит от формы и параметров нагрузки — с1 и с2. Для определения числа отрезков по £ и интервала ^ проведена серия тестовых расчетов для функции нагрузки Р(х,у,Ь). Сначала выполнялось преобразование Фурье заданной внешней нагрузки, затем функция внешней нагрузки восстанавливалась по предложенному в предыдущем разделе численному алгоритму с использованием аппроксимации подынтегральных функций
линейными функциями. Полученные численно значения сравнивались со значениями функции Р, определенной аналитически. Для с1 = 10 и с2 = 10, соответствующих области 2 х 2 м, результаты численного метода совпадают с аналитической функцией при ^ = 20 и количестве шагов больше 1000.
Вторым фактором при вычислении интегралов (28) по формулам (32), затрудняющим вычисление, являются быстрые осцилляции функций еов^Хи еов^Х—¡¿^1] с ростом При небольших значениях Ь (Ь < 10) для хорошей точности достаточно использовать значения параметров интегрирования, полученные при тестовых расчетах для функции внешней нагрузки. При увеличении Ь необходимо уменьшать шаг интегрирования по £ с сохранением длины всего интервала ^. Например, для Ь = 200 с точные результаты будут при количестве шагов больше 40 000.
Третьей проблемой численного решения является аккуратное разбиение интегралов (33). При увеличении Ь уменьшаются интервалы интегрирования второго и четвертого интегралов в (33). При вычислении этих интегралов используются разные длины шага интегрирования: первый, третий и пятый интегралы вычисляются с количеством шагов интегрирования, определенным ранее, второй и четвертый — с уменьшенным шагом для сохранения постоянного числа шагов на интервалах (^ — 5г,^г + 8г), г = 1, 2, при увеличении Поскольку при довольно больших Ь вычисление интегралов может занять значительное время на ограниченных вычислительных системах, не приспособленных для емких вычислений (например, вычисление прогибов с высокой точностью при Ь = 500 с и учетом большого числа мод может занять несколько часов, когда вычисление при Ь < 10 с занимает несколько минут), в этом случае удобнее пользоваться асимптотическим решением (36) [16, 17].
Сходимость полученных результатов по числу мод Хшо(1 в формуле (29) проверялась численно. Принципиальное значение имеет количество мод , для которых существуют системы распространяющихся от зоны нагрузки гидроупругих волн (т. е. стГп < и < стЩ+1). Эти моды вносят качественное различие в формирование прогибов льда. Количество мод зависит от выбора скорости движения нагрузки. Например, для и = 33 м/с таких мод не существует, а для и = 7 м/с такой модой является только первая (см. рис. 2 и сопутствующее обсуждение в предыдущем разделе). Остальные моды вносят корректирующий характер в формирование прогибов льда в области движения нагрузки и затухают с ростом п, £ и X. Для проверки сходимости численного решения количество Хто(1 изменялось от 3 до 15.
В работе [8] показано, что для решения в виде бегущей волны, не зависящего от времени, число мод Хто(1, равное 10, дает результаты с высокой степенью точности. В настоящей работе количество мод, если оно более значения , влияет главным образом на точность результатов в области движения нагрузки. Количество таких мод варьировалось от 5 до 10 для получения результатов высокой степени точности в области движения нагрузки.
Влияние количества мод на точность получаемых результатов показано на рис. 3, где представлены профили прогибов ледового покрова вдоль канала т* = ,ш(Х, 0,Ь)/Р0 для скорости движения нагрузки и = 7 м/с при £ = 2 с. Прогибы масштабированы амплитудой нагрузки Р0. Количество мод Хшо^ на рисунке изменяется от 1 до 3. Для указанного случая модой, для которой существует система распространяющихся волн, является только первая, п = 1. При малых временах прогибы формируются в области движения нагрузки, учета только первой моды (сплошная тонкая линия) недостаточно для получения точных результатов. При количестве мод, равном 3 (сплошная жирная
30 -20 -10 0 10 20 30 X
Рис. 3. Профили прогибов льда вдоль центральной линии канала при Ь = 2 си и = 7 м/с
линия) и более, результаты совпадают с графической степенью точности. В отдалении от нагрузки в направлении X при малых временах влияние количества мод на получаемые результаты уменьшается.
Эволюция прогибов льда во времени показана на рис. 4 и 5. На рис. 4 приведены профили прогибов ледового покрова вдоль центральной линии канала для докрити-ческой скорости и = 3 м/с, и < с^. Здесь прогиб состоит из стационарной части, затухающей в отдалении от нагрузки, и нестационарной части, затухающей со временем. Стационарный прогиб под нагрузкой быстро формируется после начала движения и не изменяется со временем. Хотя гидроупругие волны, вызываемые нагрузкой, движущейся с докритической скоростью, не формируются, в отдалении от нагрузки возникает система возмущений небольшой амплитуды, распространяющаяся от зоны нагрузки и затухающая со времени. Амплитуда возмущений перед нагрузкой больше амплитуды возмущений за ней. Эти возмущения не влияют на стационарный прогиб в области движения нагрузки. Заметим, что в рассматриваемой модели, не учитывающей вязкость льда, для докритической скорости прогиб при стационарной нагрузке (тонкая сплошная линия) немного меньше прогиба на конечных временах и асимптотике.
На рис. 5 приведены профили прогибов ледового покрова вдоль центральной линии канала для скорости и = 7 м/с, которая больше минимальной фазовой скорости первой моды и меньше минимальной фазовой скорости второй моды (с^ < и < с^). В этом случае прогиб состоит из стационарной части, затухающей вдали от нагрузки,
м>*-10"7
-20 -10 0 10 20 30 40 50 60 X
Рис. 4. Эволюция во времени прогибов льда вдоль центральной линии канала при и = 3 м/с
И'* 2
1
0 -1 -2
-100 -50 0 50 100 X
ц;*
2
1 О -1 -2
-100 -50 0 50 100 X
Рис. 5. Эволюция во времени прогибов льда вдоль центральной линии канала при и = 7 м/с
нестационарной части, затухающей со временем, и нестационарной части, устанавливающейся при больших временах. За нагрузкой формируются возмущения малой амплитуды, отходящие от нагрузки и затухающие со временем. Перед нагрузкой формируются короткие гидроупругие волны. При увеличении времени £ волны приближаются к асимптотическому решению при больших временах (36) (жирная сплошная линия на рис. 5, б)) [16, 17]. При конечных значениях £ волны затухают, а их профиль вдоль канала имеет одинаковую структуру с профилями установившегося решения для ненулевого коэффициента затухания т (см. графики и = 7 м/с в [8]). При увеличении времени растет область возмущения ледового покрова и изменение профилей волн соответствует изменению профилей установившегося решения при уменьшении коэффициента затухания. Аналогично случаю докритической скорости максимальные амплитуды прогибов льда достигаются после начала движения.
Коэффициент затухания (т > 0) в уравнении (1) является фактором, регулирующим скорость затухания волн. В работе [8] показано, что при уменьшении значения т увеличивается длина интервала наблюдаемых колебаний. При этом уменьшение коэффициента затухания приводит к увеличению времени вычисления прогибов и распределения напряжений в ледовом покрове для сохранения точности вычислений. Малые значения т могут приводить к неточным результатам, и в этом случае удобнее пользоваться аналитическими формулами для асимптотик (36) [16, 17].
На рис. 6 показаны профили установившихся стационарных частей прогибов в отсутствие нестационарных слагаемых, затухающих со временем, при разных значениях т.
а
б
Рис. 6. Профили прогибов льда вдоль центральной линии канала при и = 25 м/с и разной скорости затухания
Рассмотрен случай и = 25 м/с. Для исследуемого канала этот случай соответствует трем распространяющимся волнам, одной длинной за нагрузкой и двум коротким перед ней. Эффект затухания сильнее влияет на короткие волны. При уменьшении коэффициента затухания прогибы стремятся к асимптотикам по времени. Длинные волны за нагрузкой выходят на установившийся режим быстрее коротких волн перед нагрузкой (рис. 6, б). Результаты для больших значений г хорошо описывают стационарную часть в системе (X, у, г) нестационарных прогибов с нулевым коэффициентом затухания при небольшом времени движения. При уменьшении г увеличивается длина колебаний и, соответственно, увеличивается время, необходимое для выхода нестационарных колебаний на стационарный режим.
Заключение
Рассмотрены нестационарные колебания ледового покрова в замороженном канале. Прогибы льда вызваны движением внешней нагрузки из состояния покоя вдоль центральной линии канала с постоянной скоростью. Рассмотрен случай ледового покрова, примороженного к стенкам канала. Методом преобразования Фурье и разложения образа Фурье прогибов льда на собственные колебания ледового покрова в канале решение задачи в подвижной системе координат представлено в виде ряда, коэффициенты которого даны в квадратурах. Предложен метод вычисления этих квадратур по анали-
тическим формулам, выделяющим стационарные и нестационарные прогибы в системе координат, движущейся вместе с нагрузкой. По этим формулам численно определяются деформации ледовой пластины для различных моментов времени и различных скоростей движения внешней нагрузки.
В работе [8] показано, что формирование прогибов льда зависит от скорости движения нагрузки и ее соотношения с фазовыми скоростями распространяющихся гидроупругих волн в канале. При решении нестационарной задачи, рассмотренной в настоящей статье, в зависимости от выбора скорости движения нагрузки выделяются два основных случая. При движении с докритической скоростью решение состоит из стационарной части, симметричной в направлении движения перед и за нагрузкой, и нестационарной части, затухающей во времени. Нестационарные затухающие возмущения распространяются от зоны нагрузки и не влияют на прогиб в области движения нагрузки. При движении со скоростью, большей минимальной фазовой скорости, прогиб льда формируется из затухающей вдали от нагрузки стационарной симметричной части, затухающей со временем системы малых возмущений, удаляющихся от нагрузки, и системы гидроупругих волн.
Исследовано влияние количества мод в разложении прогибов ледовой пластины. Важным фактором является количество мод, для которых скорость движения нагрузки больше минимальной фазовой скорости. Такие моды вносят качественное отличие в формирование прогибов. Остальные моды с более высокими номерами имеют корректирующий характер и важны только для определения поведения ледового покрова вблизи нагрузки. Сходимость решения при увеличении мод проверена численно. Установлено, что для исследуемых канала и нагрузки достаточно рассмотреть 5-10 мод, чтобы получить решение с высокой точностью. Проведен анализ выбора оптимальных параметров численного решения, и исследовано их влияние на получаемые результаты. Случаи малых времен и докритических скоростей движения нагрузки не имеют существенных сложностей для вычисления прогибов ледового покрова. При увеличении времени для нагрузки, заданной в форме (11), необходимо увеличивать количество шагов интегрирования с соответствующим уменьшением длины каждого шага интегрирования в интегралах (30).
Исследована эволюция прогибов льда во времени для разной скорости движения нагрузки. Показано, что для докритической и сверхкритической скоростей движения нагрузки при увеличении времени и нулевом коэффициенте затухания г прогибы приближаются к соответствующему асимптотическому решению, исследованному в [16, 17]. Отмечено, что максимальный прогиб во всех рассмотренных случаях устанавливался при больших временах, что подтверждает возможность использования асимптотического решения для оценки максимальных прогибов и удлинений.
Исследовано влияние коэффициента затухания, описывающего вязкоупругие свойства льда на формирование прогибов ледового покрова. При уменьшении коэффициента затухания прогибы, установившиеся в двигающейся совместно с нагрузкой системе координат, приближаются к их асимптотикам по времени. Показано, что длинные волны выходят на установившийся режим быстрее коротких. При уменьшении коэффициента затухания увеличивается время, необходимое для выхода нестационарных колебаний на стационарный режим.
Благодарности. Работа выполнена при частичной финансовой поддержке РФФИ (грант № 16-08-00291).
Список литературы / References
[1] Moving loads on ice plates / V. Squire, R. Hosking, A. Kerr, P. Langhorne. Dordrecht: Kluwer Acad. Publ., 1996. 230 p.
[2] Hosking, R.J., Sneyd, A.D., Waugh, D.W. Viscoelastic response of a floating ice plate to a steadily moving load // J. of Fluid Mechanics. 1988. Vol. 196. P. 409-430.
[3] Guyenne, P., Parau, E. Computations of fully nonlinear hydroelastic solitary waves on deep water // J. of Fluid Mechanics. 2012. Vol. 713. P. 307-329.
[4] Tabata, T. Studies on visco-elastic properties of sea ice // Proc. of the Conf. Arctic Sea Ice. US, Easton, Maryland: US National Acad. of Sci., 1958. P. 139-147.
[5] Прикладные задачи динамики ледяного покрова / В.М. Козин, В.Д. Жесткая, А.В. По-горелова и др. М.: Академия естествознания, 2008. 329 c.
Applied problems of the dynamics of ice cover / V.M. Kozin, V.D. Zhestkaya, A.V. Pogorelova et al. Moscow: Akademiya Estestvoznaniya, 2008. 329 p. (In Russ.)
[6] Жесткая В.Д., Козин В.М. Исследования возможностей разрушения ледяного покрова амфибийными судами на воздушной подушке резонансным методом. Владивосток: Дальнаука, 2003. 161 с.
Zhestkaya, V.D., Kozin, V.M. Ice breaking with air-cushion vessels using a resonant method. Vladivostok: Dal'nauka, 2003. 161 p. (In Russ.)
[7] Козин В.М. Резонансный метод разрушения ледяного покрова. Изобретения и эксперименты. М.: Академия естествознания, 2007. 355 с.
Kozin, V.M. Resonant method of ice cover breaking. Inventions and experiments. Moscow: Akademiya Estestvoznaniya, 2007. 355 p. (In Russ.)
[8] Shishmarev, K., Khabakhpasheva, T., Korobkin, A. The response of ice cover to a load moving along a frozen channel // Applied Ocean Res. 2016. Vol. 59. P. 313-326.
[9] Brocklehurst, P., Korobkin, A., Parau, E.I. Hydroelastic wave diffraction by a vertical cylinder // Philosoph. Transact. of the Royal Soc. A: Math., Phys. and Eng. Sci. 2011. Vol. 369(1947). P. 2832-2851.
[10] Sturova, I.V. Motion of an external load over a semi-infinite ice sheet in the subcritical regime // Fluid Dynamics. 2018. Vol. 53(1). P. 49-58.
[11] Sturova, I.V., Tkacheva, L.A. The motion of pressure distribution over a free surface near the edge of ice sheet // IOP Conf. Ser.: Earth and Environ. Sci. 2018. Vol. 193, No. 1. 012065.
[12] Korobkin, A.A., Malenica, S., Khabakhpasheva, T. Interaction of flexural-gravity waves in ice cover with vertical walls // Philosoph. Transact. of the Royal Soc. A: Math., Phys. and Eng. Sci. 2018. Vol. 376. 20170347. Available at: https://doi.org/10.1098/rsta.2017.0347
[13] Коробкин А.А., Хабахпашева Т.И., Папин А.А. Математические модели снежно-ледового покрова. Барнаул: Изд-во АлтГУ, 2013. 116 с.
Korobkin, A.A., Khabakhpasheva, T.I., Papin, A.A. Mathematical models of snow-ice cover. Barnaul: Izd-vo AltGU, 2013. 116 p. (In Russ.)
[14] Korobkin, A., Khabakhpasheva, T., Papin, A. Waves propagating along a channel with ice cover // Europ. J. of Mechanics — B/Fluids. 2014. Vol. 47. P. 166-175.
[15] Батяев Е.А., Хабахпашева Т.И. Гидроупругие волны в канале со свободным ледовым покровом // Изв. РАН. МЖГ. 2015. № 6. C. 71-88.
Batyaev, E.A., Khabakhpasheva, T.I. Hydroelastic waves in a channel covered with a free ice sheet // Fluid Dynamics. 2015. Vol. 50(6). P. 775-788.
[16] Shishmarev, K.A., Khabakhpasheva, T.I., Korobkin, A.A. Large time behaviour of the ice cover caused by a load moving along a frozen channel // Proc. 33rd Intern. Workshop on Water Waves and Floating Bodies / Ed. Yves-Marie Scolan. April 4-7, 2018, Guidel-Plages, France. P. 185-189.
[17] Khabakhpasheva, T., Shishmarev, K., Korobkin, A. Large-time response of ice cover to a load moving along a frozen channel // Applied Ocean Research. 2019. Vol. 86. P. 154-165.
[18] Schulkes, R., Sneyd, A. Time-dependent response of floating ice to a steadily moving load // J. of Fluid Mechanics. 1988. Vol. 186. P. 25-46.
[19] Timoshenko, S., Woinowsky-Krieger, S. Theory of plates and shells. 2nd ed. New York: McGraw-Hill, 1959. 591 p.
Поступила в 'редакцию 28 декабря 2018 г.
Unsteady deflection of ice cover in a frozen channel under a moving load
Shishmarev, Konstantin A.1'*, Khabakhpasheva, Tatyana I.2 1Altai State University, Barnaul, 656049, Russia
2Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Shishmarev, Konstantin A., e-mail: shishmarev.k@mail.ru
The non-stationary problem of a moving external load along the surface of a frozen channel is studied in this paper. The channel has a rectangular cross section and is filled with an ideal fluid. The ice cover on the surface of the liquid is modelled by a thin elastic or viscoelastic plate (Kelvin —Voigt model), which is frozen to the channel walls. The external load is modelled by a local pressure moving at a constant speed along the channel. The ice-fluid-load system initially is at rest. The flow caused by both the moving load and the ice deflection is supposed to be potential. The problem is solved using the Fourier transform along the channel and by the method of normal modes. The presence of the channel walls and boundary conditions on them lead to an infinite (countable) number of dispersion relations and critical speeds for the hydroelastic waves propagating in the channel. Depending on the speed of the load, the form of the ice cover oscillations has one of two main cases. When the load moves at a speed lower than the first critical speed, the ice deflection is localised in the vicinity of the load. When the load is moving at supercritical speed, a system of hydroelastic waves propagating from the load is formed. The number of these wave systems is finite and depends on the relation of the speed of the load with the critical speeds for the channel. It is shown that as time increases, the solution for the ice deflection converge to a steady state in the coordinate system moving with the load. The influence of the retardation time on the formation of the ice deflections is investigated.
Keywords: ice plate, channel, hydroelastic waves, external load.
Cite: Shishmarev, K.A., Khabakhpasheva, T.I. Unsteady deflection of ice cover in a frozen channel under a moving load // Computational Technologies. 2019. Vol. 24, No. 2. P. 111-128. (In Russ.) DOI: 10.25743/ICT.2019.24.2.010.
Acknowledgements. This work has been partially supported by the RFBR (grant No. 16-08-00291).
Received December 28, 2018
© ICT SB RAS, 2019
