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

Вычислительные технологии

Том 16, № 6, 2011

Плановая модель течений для Новосибирского водохранилища*

В. А. Шлычков Институт водных и экологических проблем СО РАН (Новосибирский филиал), Россия e-mail: slav@ad-sbras.nsc.ru

Численная модель, основанная на двумерных (в горизонтальной плоскости) уравнениях Сен-Венана, адаптирована к морфометрическим условиям Новосибирского водохранилища. Задача воспроизведения пространственной структуры течений решается на криволинейной сетке, учитывающей особенности геометрии водного объекта, численная схема строится на основе метода конечных объемов. Модель верифицирована по данным инструментальных измерений скорости.

Ключевые слова: численное моделирование, русловый поток, водохранилище, конечно-разностная схема.

Введение

Современные подходы к обеспечению устойчивого функционирования гидротехнических сооружений предполагают широкое использование математических моделей водных объектов для целей мониторинга и расчета экстремальных режимов течения. Такие модели могут оказаться эффективными и при оценке возможных последствий катастрофических событий — залповых сбросов загрязняющих веществ, нарушения регламентированных условий эксплуатации (подобно аварии на Саяно-Шушенекой ГЭС в 2009 г.) и др. Между тем надежной гидродинамической модели Новосибирского водохранилища с возможностью детализации пространственной структуры потока до настоящего времени не создано.

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

Водохранилище как часть русла р. Оби служит основным источником воды для Новосибирска и является важнейшей составляющей экосистемы города. Длина водохранилища по судовому ходу составляет 180 км (рис. 1), полная емкость — 8.8 км3, полезный (регулируемый) объем — 4.4 км3 [2].

* Работа выполнена в рамках Интеграционного проекта № 23 СО РАН.

1. Постановка задачи

Дня расчета гидравлических параметров речного потока используется двумерная система уравнений Сен-Вепапа, В горизонтальной плоскости введем декартову систему координат с осями х, у. Поверхность руслового ложа зададим уравнением г = гь(х,у), где гь — функция, описывающая рельеф дна. Для записи уравнений плановых течений используем форму |3|, дополненную слагаемыми турбулентного перемешивания

дНи д Кии д Нпь , д( д д ди д ди

+ — + — = -дп -г— - и и + — кБ— + —НО—, дЬ дх ду дх НС 2 дх дх ду ду

дки ^ дкии ^ дкш д ^ д ^^ди

д1 дх ду ду кС23 дх дх ду ду'

д к дик дик д I дх ду

где Ь — время; ( = К + — уровень свободной поверхности; К — глубина потока; и,у — компоненты вектора горизонтальной скорости, осредпеппой по глубине; д — ускорение силы тяжести; Св — коэффициент Шези; и = л/и'2 + V2 — модуль скорости течения, I) — коэффициент турбулентного обмена.

Сформулируем краевые условия. На входном створе, который позиционируем вблизи г, Камень-на-Оби (см, рис, 1), будем считать известным суммарный речной расход Ql.

В поперечнике выходного створа, расположенного у плотины Новосибирской ГЭС, за-

К

Постаповку гидродинамической задачи замыкают начальные условия па компоненты скорости и = V = 0 и пространственное распределен не глубин К в момент Ь = 0,

2. Пространственная дискретизация и сеточная область

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

Плотина

Рис. 1. Контуры сеточной области и ноле глубин Новоеибирекснх) водохранилища в плане. 1, 2 контрольные створы измерения скоростей течения

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

В безвихревом потоке для уравнений мелкой воды справедлив интеграл Бернулли, связывающий кинетическую и потенциальную энергии с потенциалом скорости [4]. Это обстоятельство может быть использовано при построении пространственной аппроксимации уравнений на криволинейной сетке так, чтобы при тех же условиях имел место дискретный аналог интеграла Бернулли,

Для изложения метода решения задачи запишем уравнения движения из (1) в упрощенной форме, полагая D постоянной величиной:

du du du d( dv dv dv d( .

— + u— + v— + g— = DAu, — + u— + v— + g— = DAv, (2) dt dx dy dx dt dx dy dy

d u d v

где A — оператор Лапласа, В безвихревом потоке, т.е. при —— = ——, система (2)

dy dx

допускает первый интеграл [4], который в случае установившегося течения имеет вид

Е = Un* + v") + g( - D (p + p) = const. (3)

2 \dx dy)

Заметим, что соотношение (3), содержащее вязкие члены, получено впервые. Интеграл Бернулли определяет фундаментальные закономерности перераспределения гидравлических параметров в потоке и служит основой для получения широкого класса инженерных формул, В равномерных потоках интеграл Бернулли выполняется тождественно, Для описания динамики двумерных течений численную модель желательно сформулировать так, чтобы свойство (3) решения сохранялось автоматически.

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

Рассмотрим в плоскости (x, y) криволинейную сеточную область, отражающую пространственную геометрию водного зеркала. Выделим отдельный конечно-разностный бокс в виде выпуклого четырехугольника ABCD. Нормали к граням расположим, как показано на рис, 2, а их координаты обозначим m,n — например, для нормали к грани AD запишем nAD = (mAD, nAD), Длину вектора нормали удобно принять равной длине опорной грани. Центр бокса ассоциируем с индексом i, j матрицы координат x¿,j yi,j и определим в этой точке значения дискретной глубины hij. Базовые компоненты скорости позиционируем на гранях бокса ABCD и обозначим их ¿ vij_i, вводя тем самым смещенные сетки. Помимо основного ABCD рассмотрим дополнительный бокс

с центром о Т„,КС i - образованный ссродннами СТорон основного н смежного

со стороной AD бокса, и обозначим его A'B'C'D' (см, рис, 2), Элементарная область

0 *

Рис. 2. Геометрия криволинейных основного АБС О и вспомогательного А' Б' С' V' элементарных боксов и положение нормалей

А'В'СD' служит основой для расчета искомой скорости . Путем сдвига границ

2

контура ABCD в направлении грани AB аналогично формируется v-бокс для расчета

hJ 2

Проинтегрируем первое из уравнений (2) по площади A'B'CD. Используя теорему Гаусса — Остроградского и теорему о среднем в рамках метода контрольного объема, запишем результат в виде (для упрощения положим D = 0)

дЧг-У д t

wBC'uhJ - wAD'Ul_hj + wCD'u%_ 1J+1 - wAB'u%_

-g{m Q,j~m Q-I,j + ™ Q_ij+i-m Q-ij-h

AD'

CD'

.AB',

(4)

где площадь бокса, — значения искомой скорости, подлежащие

доопределению в "нештатных" целочисленных и дробных узлах, твстАВ' аппроксимируют нормальные скорости т = п •и по-разному для разных граней бокса (способ аппроксимации будет указан ниже). На прямоугольных сетках нормали ориентированы по ортам декартовой системы координат, вследствие чего тАв = тсв = 0 и дробные значения скорости в (4) исключаются. В общем случае для величин 1 необходимо сформулировать однозначный способ их расчета.

Определим вначале целочисленные как функции от базовых Воспользу-

емся подходом [5], традиционно применяемым при конструировании монотонных схем, и запишем

w

BC

1

.bc'

w

BC

u

i+i,3

к

Выражение (5) дает конструктивную основу для определения щ^. Поток импульса на грани А'В' зададим согласно выражению

ab' W UA_I

22

(б)

откуда следует расчетная формула для дробных 1 _1.

2 2

Формулы для Ущ, могут быть получены аналогично при рассмотрении урав-

нений (2) в пределах v-бoкca.

Для вывода условий консервативности в безвихревом потоке воспользуемся соотно-д и д V

шенпем = и проинтегрируем его по боксу. После подстановки результата в (4)

ду дх левая часть уравнения примет вид

¿VI+ швс' [ивс'иг,3 + Увс'щ,3) - гпАС' {иАО'иг-1,3 + +

+ - тАВ' (иАВ'+ = ... (7)

Потребуем, чтобы каждое из выражений в круглых скобках (7) представляло собой сумму квадратов, что будет эквивалентно консервативности схемы по слагаемым адвекции. Для этого сформируем неопределенные пока величины им, Vм (М = ЛВ', ВС',СО', АО' — индексы граней бокса) так, чтобы (7) перешло в квадратичную форму вида

диг

^¿—¡н1 + тВС'Кг>] ~ тАВ'Кг-^ + тСВ'к*-уц - тАВ'кг-у-ь = • • •' (8)

где К = Представим им как линейную комбинацию с неопределенными

коэффициентами

им = амщ_1^+1 + Ьмщ_1^_1 + сИщ^ + йИЩ-(9)

где 16 коэффициентов ам, Ьм, см, дм подлежат определению. Полагая, что для Vм справедливо аналогичное представление, подставим (9) в (7), Почленное сопоставление результата подстановки с формой (8) дает 10 линейных алгебраических уравнений. Еще четыре уравнения следуют из требования аппроксимации ам + Ьм + см + ¿м = 1. Оставшиеся из 16 две искомых величины служат свободными параметрами в процедуре оптимизации, которая обеспечивает корректность решения (необращение в нуль определителя системы, положительную полуопределенность искомых коэффициентов ам, Ьм, см, ¿м) при любой разумной геометрии бокса.

Перейдем к построению схемы с сохранением интеграла Бернулли, Такая схема должна быть дивергентна относительно Е (в безвихревом потоке). Учитывая связь Е = К+д( и показанную выше консервативность схемы по К, заключаем, что для обеспечения консервативности по Е способы аппроксимации величин <7_ 1 •_ 1 должны быть идентичны методикам расчета ■ Запишем соотношения (5), (6) в виде

Щ,] = + Рг^ЧЦ^, = ъ,зи^,з-1 + (10)

где вг,л 7г,Л ^ ~ Неотрицательные коэффициенты, ¡3^ = 1, 7^' += 1, и определим дискретные значения ( как

0,3 = + = 7»,.7<г-± ¿-1 + ^Сг-^Г (И)

В формулах (11) вычисляемыми величинами являются <7_1 так как поле (

по условию определено в целочисленных узлах и базовые ^г,з находятся из решения основных уравнений системы (1), Следовательно, первое из соотношений (11) является разностным уравнением первого порядка относительно <7_ 1 • при известных значениях (г^. Задача определения <7_1 ■ сводится к обращению двухдиагональной матрицы и

в общем случае оказывается плохообуеловленной. Однако значения уровня в гидравлике естественных водоемов обычно в расчетной области меняются слабо, а на отдельно взятом конечно-разностном интервале вариации ( тем более невелики, что позволяет принять вначале С,^ ] = Сад; а далее итерационно уточнить эти величины в соответствии со значениями а^, /3^. Полученные • используются для расчета 1 •_ 1 по второй формуле в (11). Расчет параметров потока на линиях уреза воды проводился методом распада разрыва по методике [6].

Таким образом, сформулирован способ расчета базовых и вспомогательных значений вектора скорости и уровня, обеспечивающий дивергентную форму конечно-разностных уравнений движения относительно полной энергии Е, фигурирующей в интеграле Верну I.т. Реализация схемы требует решения нескольких вспомогательных задач, а в целом алгоритм оказывается сложнее традиционно применяемых при моделировании процессов в водотоках [7]. Это связано с необходимостью аппроксимации дифференциальной системы в криволинейной сеточной области и с использованием смещенных сеток для разноименных полей. Преимущества смещенных сеток в данном контексте не столь очевидны — в полной мере они проявляются при интегрировании полной системы (1). В частности, для уравнения неразрывности смещенные сетки автоматически обеспечивают второй порядок точности, что особенно важно в криволинейных областях, где затруднительно определить точность схемы.

Аппроксимация по времени проводилась на основе неявных схем без использования процедур расщепления с применением метода неполной факторизации.

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

Дискретизация уравнений (1) для Новосибирского водохранилища проводилась в узлах криволинейной конечно-разностной сетки, построенной на основе решения системы квазилинейных уравнений эллиптического типа с управляющей функцией [3]. Пунктирная кривая на рис. 1 показывает границы области, в которой формировалась сетка. В качестве основных критериев при конструировании сетки использовались требования гладкости и максимально возможной ортогональности сеточных координатных линий с учетом плановой геометрии водоема. Результирующая сетка содержала около 75 тыс. узлов и обеспечивала среднее разрешение 300 х 100 м2 в продольном и поперечном к потоку направлениях. Для повышения точности описания гидравлически важных участков русла (основной водоток, область вблизи ГЭС) размеры сеточного бокса уменьшались до 130 х 30 м2,

3. Верификация модели течений

Для сопоставления расчетных параметров течения с натурными данными использованы результаты измерений гидрологических характеристик (глубин, скоростей течения и расходов воды) в Новосибирском водохранилище, выполненных 24-25 июня 2009 г. с помощью акустического доплеровского профилографа. Программно-аппаратный комплекс позволяет в динамическом режиме при движении судна по створу наблюдений получать мгновенную картину распределения скорости и направления течения.

Контрольные створы выбраны в верхней и средней частях водоема, их положение показано на рис, 1 сечениями 1, 2. Уровень свободной поверхности в дни проведения измерений равнялся 113,5 м балтийской системы, расход воды составлял 2500 м3/с. Эти данные служили в качестве краевых условий дня уравнений (1),

Поперечное сечение 1 располагается в верховьях водохранилища (12 км ниже г, Ка-мень-на-Оби), Здесь течение имеет выраженный русловый характер, размер поперечника составляет примерно 500 м. На рис.3, а показано распределение модуля скорости но длине поперечника (маршевая координата I увеличивается от правого берега к левому). Отдельные точки на рисунке отражают значения |и|, полученные в результате измерений и осредпеппые но вертикали. Заметный разброс амплитуд фактической скорости обусловлен турбулентным характером реальных течений и наличием стохастического компонента в потоке. Если провести статистическую фильтрацию случайных факторов, то получим распределение средней упорядоченной скорости переноса. Пунктирная кривая 2 па рис. 3, а построена сглаживанием данных измерений методом скользящего среднего, кривая 1 показывает профиль модельной скорости, полученный в расчете. В силу единой статистической природы (обе функции суть математическое ожидание) кривые могут быть корректно сопоставлены. Так, из рисунка следует, что достаточно малый локальный шаг сетки па поперечнике (около 50 м) позволяет удовлетворительно детализировать структуру потока. Сравнительный анализ кривых 1, 2 показывает, что модель дает несколько заниженные скорости течения (кривая 1). Ошибка модельного профиля скорости составляет в среднем 15 %, что с учетом статистической неопределенности численного решения представляется приемлемой величиной. Максимальные скорости достигаются в глубоководной части поперечника и оказываются дня обоих случаев близкими (0.7-0.8 м/с).

И, »-'-

0 0.1 0.2 0.3 0.4 {., кл.

а

|и|, м/с 0.3

0.2

0.1

0

1 2 3 4 5 км

5

Рис. 3. Расчетные (кривые 1), измеренные (отдельные точки) значения скорости и сглаженная эпюра измеренных скоростей течения (кривые 2) на контрольных створах 1 (а) и 2 (б)

Оценим точность воспроизведения скорости на створе 2, Ширина водного зеркала составляла здесь 5700 м. Данные измерений на рис, 3, б (отдельные точки) показывают, что течение характеризуется сильной неоднородностью, а скорость испытывает резкие колебания от нуля до максимума (0,3 м/с) в различных, порой соседних, узлах зондирования, Амплитуда сглаженных скоростей сравнительно невелика (0,1-0,2 м/с, см, кривую 2 на рис, 3,6), Пульсационный характер скорости можно, по-видимому, объяснить сложной трехмерной структурой потока: представленные на рисунке данные получены как результат осреднения мгновенных вертикальных профилей по глубине. Последние формируются под влиянием инерционных механизмов и эффектов нестационарности, турбулентного перемешивания, обтекания неровностей рельефа и других факторов и могут иметь вертикальное строение [1] с мелкомасштабными вариациями. Процедура осреднения данных по глубине нивелирует вертикальные особенности течения, что выражается в отмеченной сильной изменчивости средних значений на створе. При слабых течениях, характерных для рассматриваемой озерной части водохранилища, возрастает также роль стохастического фактора, связанного с ветром, который генерирует собственные возмущения скорости.

Сравнивая кривые 1,2 на рис, 3, б, констатируем, что приведенный теоретический профиль скорости в целом довольно близок к измеренному — участки спада и возрастания скорости и положение экстремумов совпадают с приемлемой точностью. Это позволяет говорить об удовлетворительном качестве воспроизведения скоростного режима водохранилища с помощью численной модели,

4. Результаты расчетов

Рассмотрим диагноз динамического режима водохранилища, полученный интегрированием системы (1) на установление. Результаты расчета отразим в терминах проточно-сти. Общая проточноеть водоема Р. характеризует быстроту обновления вод в водохранилище и в гидрологии определяется по формуле [8]

где Qs (¿) — суммарный расход па входи ом створе, V. — объем чаши водохран илища, Т — временной период (обычно один год), Проточноеть Новосибирского водохранилища по многолетним данным равна 6,7, т.е. годовой объем водного стока составляет 6,7 объемов чаши, В водоохранной проблематике интерес представляет не только суммарный показатель проточности, но и скорость водообмена в отдельных частях водохранилища (в том числе на мелководьях и в береговых зонах), т.е. показатель локальной проточности, Обобщение величины Ps проведем по формуле

т

о

т

о

где дп(х, у, ¿) — распределение расходов вдоль линии тока, V.. — локальный объем водного участка, для которого проводится расчет, В численной модели в качестве V.

Плотшга

10 км 0 5 10 15

Рис. 4. Поле коэффициента ироточиости р3 в озериой части Новосибирского водохранилища

логично принять элементарный объем, а величина рассчитывается по искомым полям расходов.

Рисунок 4 показывает распределение коэффициента проточности по акватории при-плотинной части водохранилища. Максимальное значение величина р8 достигает на русловом участке, а в озерном сегменте, где скорость течения невелика, интенсивность водообмена падает до значений 2-4. В верхнем течении водохранилища (створ I, см. рис. 1) имеет место максимум проточности, равный согласно расчетам 85.

Заключение

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

Дальнейшее развитие модели предполагается проводить в направлении детализации ветроволнового фактора. Дрейфовые течения, генерируемые ветром, можно описать поверхностными напряжениями в уравнениях движения. Систематическое воздействие волн также обусловливает формирование упорядоченного дрейфового течения (так называемый стоксов перенос), но его формально-математическое описание значительно сложнее. Учет волновых течений может оказаться важным при оценке потоков субстанции (взвеси, загрязнители) от берегов.

Автор благодарит К.В. Марусина за предоставленные натурные данные.

Список литературы

[1] Шлычков В.А. Численная модель взвесенесущего потока для Новосибирского водохранилища // Вычисл. технологии. 2010. Т. 16, № 2. С. 111-121.

[2] Савкин В.М., Двуреченская С.Я., Квон В.И., Попов П.А. Мониторинг качества воды Новосибирского водохранилища // Окружающая среда и экологическая обстановка в Новосибирском научном центре СО РАН. Новосибирск: Изд-во СО РАН, 1995. С. 141-153.

[3] Вольцингер Н.Е., Пясковский Р.В. Теория мелкой воды. Океанологические задачи и численные методы. Л.: Гидрометеоиздат, 1977. 208 с.

[4] Кочин Н.Е., Кивель И.А., Розе Н.В. Теоретическая гидромеханика. М.: Физматгиз, 1963. Ч. 1. 583 с.

[5] Карамышев В.Б. Монотонные схемы и их приложения к газовой динамике. Новосибирск: Изд-во Новосибирского гос. ун-та, 1994. 109 с.

[6] Беликов В.В., Семенов А.Ю. Численный метод распада разрыва для решения уравнений теории мелкой воды // Журн. вычисл. матем. и математической физики 1997. Т. 37, № 8. С. 1006-1019.

[7] Васильев О.Ф. Математическое моделирование гидравлических и гидрологических процессов в водоемах и водотоках (обзор работ, выполненных в Сибирском отделении РАН) // Водные ресурсы. 1999. Т. 26, № 5. С. 600-611.

[8] Китаев А.Б. Обмен вод в искусственных водоемах. Пермь: Изд-во Пермского гос. техн. ун-та, 2005. 112 с.

Поступила в редакцию 7 февраля 2011 г., с доработки — 5 апреля 2011 г.