CC BY
40
Вычислительные технологии
Том 13, № 2, 2008
Плановая динамико-стохастическая модель ледохода
В. А. Шлычков
Институт водных и экологических проблем СО РАН, Новосибирск, Россия
e-mail: slav@ad-sbras.nsc.ru
A mathematical model of the spring ice drift is formulated. Channel base current parameters are reproduced using of a deterministic-type numerical model, while ice floes floating down the river are considered as the multi-scale elements of a stochastic ensemble. Dynamics of each ice element is calculated separately and its interaction with the coast and adjacent ice floes, including those forming ice-jam, are taken into account. Jam formation scenario was calculated for the case study of the Lena river near Yakutsk.
Введение
В ходе весенних подвижек льда происходит разлом и дробление ледяного покрова в речном русле на отдельные поля и льдины, которые транспортируются потоком вниз по течению. При высокой густоте движущегося ледяного массива на сужающихся участках реки возможны остановки ледового поля, связанные с заклиниванием ядра ледохода в стесненном водотоке. В результате формируется заторное скопление льда, мощность которого нарастает со временем за счет присоединения новых масс плывущих льдин [1].
Исследования и прогнозы заторов льда на северных реках базируются, как правило, на эмпирических связях между ледовыми условиями в зимне-весенний период, температурой воздуха, расходами воды и другими количественными характеристиками половодья [2]. Еще одним средством изучения заторообразования является гидравлическое моделирование, когда натурные условия воспроизводятся с помощью физической модели уменьшенного масштаба [3]. При несомненной практической полезности такого подхода, точность его обычно невелика, а сами методики не являются универсальными, так как разработаны применительно к конкретному водному объекту с уникальными свойствами. Современные гидродинамические модели дрейфа льда в море основаны на стилизации ледового покрова в виде сплошной вязкой пленки на поверхности водоема, для которой формулируются уравнения движения и неразрывности [4].
Ниже представлена гибридная динамико-стохастическая модель весеннего ледохода в русловом потоке, в которой параметры течения воспроизводятся на основе детерминированной численной модели, а плывущие по реке льдины рассматриваются как разномасштабные элементы стохастического ансамбля. Динамика отдельного ледового элемента рассчитывается индивидуально с учетом взаимодействия с берегом и близлежащими льдинами, в том числе и находящимися в заторной области.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
1. Постановка задачи
Для расчета гидравлических параметров речного потока используется двумерная модель плановых течений. В горизонтальной плоскости введем декартову систему координат с осями х, у. Поверхность руслового ложа зададим уравнением z = Zb(x, у), где Zb — функция, описывающая рельеф дна. Уравнения плановых течений запишем в виде [5]
dhu д huu д huv д (h + zb) g
Ht + + ^ = -gh~x--C2|u| u + Tx,
dhv д huv д hvv д (h + zb) g
Ht + ^ + ПУТ = -gh^ûT - C2|u|v +Ty'
дh дuh дvh ^ (1)
дЬ дх ду '
где t — время, h — глубина потока, u,v — компоненты средней по глубине горизонтальной скорости, g — ускорение силы тяжести, Cs — коэффициент Шези, |u| = \Ju2 + v2 — модуль скорости течения, тх,ту — напряжения ветра.
Сформулируем краевые условия. На входном створе, который позиционируем на некотором поперечнике русла, будем считать известным суммарный речной расход Q\. В поперечнике выходного створа задается уровень свободной поверхности, пересчитанный в термины глубин h.
Постановку гидродинамической задачи замыкают начальные условия на компоненты скорости u = v = 0 и пространственное распределение глубин h в момент t = 0.
Задача решалась методом конечных разностей по пространству и времени. Использовались криволинейные сетки с узлами, разнесенными по граням элементарного пространственного бокса. "Расшатанные" сетки позволяют эффективно строить консервативные схемы второго порядка точности на основе центральных разностей, а применяемые неявные методы обеспечивают нужный запас устойчивости при многошаговом интегрировании.
Пространственная дискретизация дифференциальных операторов основана на представлениях о схемах с невозрастанием полной вариации (Total Variation Diminishing, TVD), гарантирующих монотонность решения путем использования перестраивающегося шаблона и подходящего выбора аппроксимации производных в различных участках численного решения. Монотонизация схемы проводилась по методу Куранта—Изаксо-на—Риса [6]. Неявная часть TVD-операторов алгоритмизировалась в соответствии с методом [7].
Важным этапом вычислительной процедуры является алгоритм разбиения расчетной области на совокупность элементарных объемов, на которых строятся дискретные соотношения. Криволинейные сетки заданной структуры эффективно адаптируются под плановую геометрию расчетной области и учитывают особенности морфометрии водотока. Для участка русла р. Лены, по которому проводился расчет, целесообразным оказалось применение дифференциально-вариационного подхода, позволяющего регулировать свойства проектируемой сетки с помощью управляющих функций. Задача формирования системы сеточных узлов в расчетной области ставится для системы квазилинейных уравнений эллиптического типа [8]. Критерием оптимальности сетки служило требование максимальной пространственной детализации в областях основного водотока. Результирующая сетка содержит около 15000 узлов и обеспечивает среднее
разрешение 83 м в продольном и 17 м в поперечном направлениях в обоих рукавах. Конфигурацию водного зеркала и структуру сетки можно видеть на рис. 1.
Конечно-разностные аналоги системы (1) получены на основе метода конечных объемов, в котором уравнения проецируются в элементарный бокс со сторонами, близкими по направлению к криволинейным границам русла. Алгебраические системы конечно-разностных уравнений решались методом неполной факторизации.
Явление ледохода в модельном представлении имитируется совокупностью свободно плавающих ледяных обломков, которые рассматриваются как система материальных тел конечных размеров на акватории. На льдину, находящуюся в ледоходном потоке, действуют следующие основные силы:
— лобовое сопротивление воды, пропорциональное линейному размеру и толщине элемента;
— касательные влекущие силы со стороны потока, воздействующие на нижнюю поверхность льдины;
— ударные силы столкновений с соседними элементами;
— силы реакции со стороны берега и касательное трение при наползании на сушу;
— силы донного трения льдины на отмелях;
— боковое давление на льдину в сплошном заторе;
— случайные возмущения (ветер, турбулентные пульсации в потоке).
Толщина кг, линейный размер (диаметр) йг и площадь льдин вг предполагаются различными для разных элементов ансамбля и задаются случайным образом. Движение отдельной льдины при моделировании ледохода отождествляется с движением ее центра масс. Взаимосвязанная система уравнений движения ансамбля выглядит следующим образом:
тг
N
-РгСе(А1Мг + А2вг)|ц - и|(ц - и) + ^ тг тп Шгп (г г - Гп) - РгсеМ М; В Ц +
й г г
п= 1
иг
(2)
где г = 1,..., N; N — количество льдин в области; гг — радиус-вектор центра масс; иг — вектор скорости г-й льдины; ргсе — плотность льда; тг — масса элемента; — случайное динамическое воздействие на льдину.
Рис. 1. Сеточная структура расчетной области и начальный стохастический ансамбль льдин на акватории (серая заливка) р. Лены
Первое слагаемое в (2) отражает форсинг со стороны речного потока: учитываются лобовое сопротивление с коэффициентом Ai и касательные напряжения на нижней поверхности льдины с коэффициентом A2. Второе слагаемое описывает изменение импульса в результате возможных столкновений льдины i с ближайшими соседними льдинами массой mn. Элементы финитной матрицы W формируются как мера взаимного перекрытия льдин. Значения Win являются ненулевыми только в случае непосредственного контакта льдин, когда центры сближаются на расстояние не менее суммы диаметров di + dn. Между льдинами возникают силы отталкивания, так что компоненты матрицы имеют вид
а дискретное во времени ударное взаимодействие аппроксимируется "эластичным" образом. Коэффициенты эластичности ^т определяют величину сил упругого взаимодействия.
Введением эластичности в модели соударений можно решить проблему одновременного взаимодействия нескольких элементов, растягивая соударения во времени и естественным образом обеспечивая перераспределения массы по акватории. При наличии внешнего потока взаимное пересечение льдин позволяет описать процессы уплотнения ледяного массива в заторе и сформировать поле бокового давления, приводящего к разрушению части льдин и образованию эффектов торошения, а также многослойной структуры затора за счет наползания льдин друг на друга. В вычислительном плане этот прием избавляет от необходимости итерационной обработки множества попарных взаимодействий и обеспечивает логическую простоту алгоритма, что весьма важно при анализе ансамбля большой размерности.
Третье слагаемое в (2) характеризует вклад сил трения льдин об отмели и о берега. Береговая полоса, в которой развиваются процессы торможения льдин, определяется как мелководная зона с глубиной Н < Н^се, где НгСе — средняя толщина льдин. Характер влияния берега на продольный и поперечный (относительно линии берега) компоненты скорости льдин неодинаков — очевидно, что нормальная скорость должна затухать быстрее касательной вблизи уреза воды. Для учета анизотропии по скоростям служит матрица В, которая в простейшем случае совпадения линии берега с осью х содержит на главной диагонали коэффициент гидравлического трения льдин о берег Ад и коэффициент затухания поперечной скорости А^.
2. Результаты расчетов
Сформулированная плановая модель ледохода апробирована на участке р. Лены в районе Якутска. При весеннем расходе ^=10400 м3/с ширина основного русла составляет 1-2 км. На рис. 1 показаны геометрия русла в плане (серая заливка), структура конечно-разностной криволинейной сетки и начальное распределение массива льдин в области. Водное зеркало имеет сравнительно сложную конфигурацию: острова делят поток на два рукава — судоходный левый и мелководный правый (направление основного течения показано стрелками), которые в половодье соединяются протокой.
Координаты начальной совокупности элементов ледового ансамбля задавались с помощью датчика случайных чисел с равномерным распределением. Габариты льдин формировались на основе нормального закона со средним значением диаметра 150 м
и дисперсией 50 м. Плотность ледохода полагалась равной 0.4. Эти соотношения поддерживались при генерации новых элементов на входном створе взамен ушедших вниз по течению. Площадь акватории в пределах расчетной области составляет 26.7 км2, суммарная площадь, покрытая льдинами, равна 10 км2. Общее число элементов ансамбля в области не превышает 1500.
В качестве начальных гидродинамических полей задавался стационарный плановый поток, полученный интегрированием уравнений (1) на установление без учета ледовых эффектов. Профиль уровенной поверхности вдоль маршевой координаты I, следующей по динамической оси потока в левом рукаве, показывает кривая 1 на рис. 2. Расчетный кинематический режим характеризуется средней скоростью 1.5-1.8 м/с, которая местами увеличивается до 2.7 м/с. Основной водоток проходит по левому рукаву, правый пропускает не более 20 % расхода.
Решение задачи воспроизведения ледохода проводилось путем совместного интегрирования систем (1) и (2). Численные эксперименты показали, что при небольшой плотности ледохода и отсутствии крупных элементов плавающий лед транзитом продвигается по руслу и покидает расчетный участок через выходной створ. Лишь в сравнительно узком правом рукаве при прохождении крупных льдин высока вероятность образования локального затора, не влияющего на проводку льда в основном (левом) водотоке.
С увеличением дисперсии выборки и размеров элементов становится возможной закупорка основного русла одной или несколькими крупными льдинами за счет заклинивания о берега. Наблюдениями установлено, что основным фактором задержки движения колотого льда становится стеснение потока вследствие, например, сужения русла. В этом случае образуется ледяная перемычка, перекрывающая русло и служащая очагом затора. Перемычка препятствует свободному сплаву льда и обусловливает скопление льдин позади очага, что со временем приводит к увеличению протяженности заторного поля.
Заторная конфигурация ансамбля в период формирования ледяного препятствия показана на рис. 3. Русло левого рукава в точке с координатой I ^ 9 км перекрывается группой крупных льдин вблизи соединяющей протоки за счет эффектов распора и трения о берега (показано стрелкой). Заметим, что задержка движения льда возникает здесь как результат естественного развития процесса взаимодействия ансамбля с геометрией русла, т. е. имитирует природное явление в рамках модельных представлений.
В условиях формирования затора по мере аккумуляции массы льда в хвостовой части увеличиваются продольные напряжения в ледовом покрове, на сплоченные ледяные
Рис. 2. Геометрия свободной поверхности безледного потока (1) и при формировании затора (2)
Рис. 3. Конфигурация ледяного ансамбля в период закупорки русла: стрелки — очаги заторов
обломки действует избыточное давление и в результате сжатия часть льдин крошится и выдавливается со своей горизонтали. Происходит торошение затора с увеличением толщины ледяного поля и его смещением вниз под действием силы тяжести. Проседание ледяного покрова приводит к сужению поперечного сечения, а увеличением шероховатости нижней поверхности обусловливается возрастание общего гидравлического сопротивления. Эти эффекты воспроизводятся в рамках эластичной модели столкновений, изложенной выше. Уменьшение живого сечения снижает пропускную способность русла, т. е. обусловливает изменение условий движения водного потока, вследствие чего происходит его торможение, а затем и подъем уровня в хвосте затора. Кривая 2 на рис. 2 иллюстрирует влияние затора на структуру свободной поверхности £, сформировавшуюся за 1.5-часовой период с момента блокирования русла. Отчетливо выделяется зона подпора, в которой наблюдается приращение уровня воды, превышающее 1 м по отношению к безледному потоку.
Заключение
Представленные результаты отражают вероятный сценарий генезиса затора торошения при ледоходе и показывают возможности использования численной модели как инструмента предсказания динамики подпорных явлений, весьма важных при описании ледовых наводнений и формирования волн прорыва. В рабочем режиме прогностические результаты могут быть получены методом статистических испытаний на основе многовариантного прогона модельных ансамблей с различными наборами характеристик, отвечающими разбросу параметров реальной ледовой обстановки.
Развитая динамико-стохастическая модель является, по-видимому, первой (возможно, одной из первых) моделей подобного рода для описания ледохода. Ее преимущество — в отказе от традиционной концепции сплошной среды, т. е. в модели поле колотого льда представляется не в виде непрерывной субстанции, а в виде ансамбля дискретных элементов, что, несомненно, ближе к натуре и позволяет с большей надежностью воспроизводить параметры ледохода. Очевидная практическая направленность разработки связана с выявлением затороопасных участков рек и оценкой вероятности закупорки русел в зависимости от статистических характеристик ледохода. Другой важный аспект — анализ эффективности мероприятий по расчистке русел для беспрепятственной проводки льда в половодье. Актуальным является расчет максимального ударного давления льда на плотины, опоры мостов и других инженерных сооружений на ре-
ках. Без принципиальных изменений модель может быть адаптирована для прогноза внезапного появления морского льда на устьевых акваториях рек за счет нагонных эффектов.
Список литературы
[1] Берденников В.П. Динамические условия образования заторов льда на реках // Труды ГГИ. 1964. Вып. 110. С. 3-11.
[2] Кильмянинов В.В. Условия формирования наводнений при заторах льда на Средней Лене в 1998 и 1999 годах // Метеорология и гидрология. 2000. № 10. С. 93-98.
[3] Бузин В.А. Заторы льда и заторные наводнения на реках. СПб.: Гидрометеоиздат, 2004. 204 с.
[4] КолЕсов С.А. Моделирование дрейфа льда в арктическом бассейне // Труды ГГИ. 1990. Вып. 420. С. 32-38.
[5] Стокер Дж.Дж. Волны на воде. М.: Изд-во иностр. лит., 1959. 618 с.
[6] Куликовский А.Л., ПогорЕлов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 606 с.
[7] Harten A. On a class of high resolution total-variation-stable finite-difference schemes // SIAM J. of Numer. Anal. 1984. Vol. 21, N 1. P. 1-23.
[8] ВольцингЕр Н.Е., Клеванный К.А., ПЕлиновский Е.Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989. 272 с.
Поступила в редакцию 27 сентября 2007 г.
