Моделирование неустановившегося движения воды в нижнем бьефе Богучанской ГЭС* Текст научной статьи по специальности «Физика»

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

Том 13, Специальный выпуск 2, 2008

Моделирование неустановившегося движения воды в нижнем бьефе Богучанской ГЭС*

Е.Д. Карепова

Институт вычислительного .моделирования СО РАН, Красноярск, Россия

e-mail: jane@icm.krasn.ru

Г. А. Федоров Сибирский федеральный университет, Красноярск, Россия

Presence of a hydropower station changes the regime of water flows. It also rises the question on an estimation of the influence of the operating conditions of a hydropower station on the tail water flow and on the modeling of an emergency flow discharge.

This paper addresses the modeling of an unsteady flow in r. Angara in the tail-water of the Boguchanskaya hydropower station. This problem is solved numerically using the Sen-Venan equations. The method of river-indicator basins developed by A.V. Ogievsky is employed for estimation of a side inflow. Results of testing of the model using the data of observations for spring and summer floods in 2004-2005 years are presented. Results of numerical modeling of various unsteady flow regimes in the tail-water of the Boguchanskaya hydropower station are provided

1. Моделирование отекания воды по руслу

При моделировании используются уравнения Сен-Венана, описывающие одномерное неустановившееся движение воды по открытому руслу под действием силы тяжести. В таких моделях считается, что центробежный эффект, связанный с извилистостью русла, пренебрежимо мал, поэтому, в частности, свободная поверхность принимается горизонтальной в каждом сечении: £ = £(t, x)(t — время, ось Ox направлена вдоль русла). Кроме того, движение предполагается медленно изменяющимся, что позволяет не учитывать местные потери напора (например, вследствие резкого сужения/расширения русла). Несмотря на одномерность, уравнения Сен-Венана учитывают параметры сечения русла в интегральных характеристиках, таких как площадь живого сечения и и осредненная по и пропускная способность русла K. Эти характеристики зависят прежде всего от уровня воды в сечении. Коэффициент шероховатости n, описывающий сопротивление подстилающей поверхности, принимается разным в русле и при выходе воды па пойму и также осредняется по и (рис. 1).

Рассмотрим следующую схематизацию русла. Зададимся декартовой системой координат таким образом, чтобы ось Oz была направлена вертикально вверх. Пусть известна геометрия рельефа русла и прилегающей поймы в некоторых сечениях реки x = Xj,

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 05-01-00579, № 08-01-00621) и Президентской программы “Ведущие научные школы РФ” (грант № НШ-3431.2008.9).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

п=0.02

Расстояние от условного нуля створа, м

Ь2(х,ф Ых,£)

Рис. 1. Пример схематизации иоиеречнш'о се- Рис. 2. Сечение русла

чвния

г = 0,..., N (рис, 2). Дно описывается функцией г = Н(х,у), а при усреднении по ширине в каком-либо смысле (см., например, [1]) г = Н(х). Пусть также свободная поверхность описывается функцией г = £(£, х), следовательно, Л.(£, х) = £(¿, х) — Н(х) — глубина реки в сечении. Дно представляется кусочно-линейной функцией с заданным для каждого линейного элемента уклоном, который считается положительным в сто-

дН

рону уменьшения отметок дна русла ?’о = —^—■

дх

Для получения дискретного аналога рассматривается следующий вид уравнений Сен-Венана:

ди д(з

Ж+& = «' (1)

оя а а? . _

lS + lh\~)+g'JllFx-gш, ~ °- (2)

Здесь £ — отметка уровеня свободной поверхности (в расчетах используется Балтийская система высот — м БС); Q — расход воды; и = и(£) — площадь живого сечения; q — поток бокового притока (расход воды через единицу длины); г/ = Q2/К2 — уклон трения,

К = Сил/Ё~н — пропускная способность русла, О = Л.1/6/и — коэффициент скоростного

напора Шези, К}г = и/\к — гидравлический радиус, Хь(^,Ь) — смоченный периметр; и

При расчетах неустановившегося течения в русле р. Ангары рассматривается только докритическое течение (Гг < 1), требующее задания по одному граничному условию в верхнем и нижнем створах. В верхнем сечении задаются наблюдаемые значения расходов Q = Q(t,xо), а в нижнем створе — зависимость Q = Q(£(xN-1)). Последняя зависимость является характеристикой замыкающего створа, и при ее получении используются натурные данные (рис. 3). Однако следует помнить, что кривая расходов

Хм

25000 0 20000 ,р

й* 15000 2

О 10000

г

га

5000

о

83 84 85 86 87

отметки уровня, м БС

89

Рис. 3. Зависимость Q(£) для участка русла в районе гидропоста Татарка, но данным наблюдений весенне-осеннего периода 2004 2005 гг.

Рис. 4. Четырехточечный шаблон

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

Для корректной постановки задачи следует также задать начальные данные по всей длине исследуемого участка русла: £ = £(0, ж), = <^(0,ж). В качестве начальных

данных используется решение уравнения Бернулли, которое для открытых русел (пьезометрическая высота совпадает с отметкой свободной поверхности) записывается как уравнение баланса для энергии в двух соседних сечениях ж*, ж*+1:

«*.) + 2^ + я„,,( = я**,) +

г = 0, - 2.

(3)

2$ ^ 2д

Здесь вторые слагаемые в левой и правой частях называются скоростным напором в сечении; скорость и берется усредненной по сечению; а* и а*+1 — средневзвешенные коэффициенты скорости, учитывающие неравномерность распределения скорости по сечению; слагаемое Иад* описывает потерю удельной энергии потока на участке. Потеря удельной энергии потока между двумя сечениями состоит из потерь па преодоление сил трения и потерь при расширении/сжатии потока (местные потери). Для учета потерь, связанных с трением, в гидравлике используется величина г;/Дх, где Дж — расстояние между сечениями. Потери при сжатии (расширении) русла учитываются как потери кинетической энергии при помощи эмпирического безразмерного коэффициента В, Таким образом, балансовые соотношения (3) можно считать системой нелинейных уравнений относительно £(0, ж*), г = 1,..., N — 1, решение которой, например методом хорд, определяет начальные отметки свободной поверхности по длине реки, соответствующие некоторым расходам воды, принимаемым постоянными па каждом участке и задаваемым по данным наблюдений. Для корректного расчета начальных отметок уровня следует задать граничное условие. При докритическом течении это отметка уровня ВОДЫ В нижнем створе £(0, жм_1).

Заданные по длине участка реки расходы и соответствующие им полученные при решении уравнений Бернулли отметки уровня будем считать начальными данными при решении уравнений Сен-Венана, Однако сами эти данные уравнениям не соответствуют. Поэтому на первых шагах по времени в верхнем створе задается постоянный расход, соответствующий начальному. Таким образом, начальные данные “пересчитываются” по модели Сен-Венана, Количество таких “подготовительных” шагов по времени должно быть достаточным для установления течения по всему руслу, поэтому оно зависит от длины рассматриваемого участка и средней скорости течения. Следует отметить, что при наличии сильно меняющегося во времени бокового притока такая процедура не даст правильных результатов, а неиспользование предварительного установления приводит к несколько большим ошибкам в первые итерации по времени. Причем ошибки тем больше, чем дальше исследуемый створ от начального и чем отчетливее картина неустановившегося течения, наблюдавшаяся в период, непосредственно предшествующий расчетам.

Поскольку на исследуемом участке реки имеются крупные боковые притоки, при моделировании следует задавать q в уравнении (1), Оценки величины притока проводятся по методу бассейнов-индикаторов А,В, Огиевского с использованием данных о ежедневных расходах воды рек-аналогов, впадающих в р, Ангару на рассматриваемых участках. Расходы аналогов Qa(t) на однородных (по ланшафтному устройству площади водосбора) участках делятся на площади водосборов аналогов Sa, и полученные отношения умножаются на площадь водосбора всего соответствующего участка Sy4, приведенную к его длине Ly4 по руслу реки:

9у,М = (4)

Sa Ly4

При дискретизации уравнений Сен-Венана (1) и (2) используется четырехточечная неявная разностная схема [2-4],

Для дискретизации рассмотрим сетку, в общем случае неравномерную по пространству (заданную вдоль фарватера русла) и равномерную по времени:

= (zf = (xi, tk) : xo < ■ ■ ■ < Xi ■■■ < xn-i, tk = kr, k = 0,1,... }. (5)

Здесь N — количество створов (сечений); х0 — входной створ; xn-1 — замыкающий створ; т — шаг по времени; hi = xi+1 — xi; i = 0,1,..., N — 2 — шаг по пространству. Рассмотрим следующую аппроксимацию функций и их производных на четырехточечном шаблоне (рис, 4):

R = \(ff + ft.). / ^ff + o (ff+1 - ff),

0/ и 9f r: l.f^-ff + 0 ((Л+? - ./?+1) - (./?+! ~ ft))

dt т ’ dx 2 hi

Здесь принято обозначение fk = f (tk,xi). Таким образом, значения функций и производных по пространству аппроксимируются в некоторой точке (xi+0.5hi5 tk+0т) каждой ячейки сетки (xi, xi+1) х (tk, tk+1), Здесь 9 — параметр схемы, В работах [2-4] доказаны безусловная устойчивость разностной схемы при 0.5 <9 < 1 и условная устойчивость при 9 = 0.5, Там же доказана сходимость для длинноволновых течений.

Дискретизация уравнений Сен-Венана (1) и (2) на сетке (5) дает систему нелинейных алгебраических уравнений, которая линеаризуется при следующих предположениях,

1, Квадрат приращения функций по времени пренебрежимо мал, т, е, отбрасываются члены второго порядка малости: (Д/,к+1)2 ~ 0 Д/,к+1 = /,к+1 — /,$■

2, Для функций, зависящих от расхода и отметки уровня, / = /(0, £) используется следующая линеаризация:

д /*+1 ~ АОк+1 + Д£*+1.

дQ Чг ас

При этих предположениях па каждом временном шаге ¿к+1 = (к + 1)т, к = 0,1,..., система (1) и (2) может быть записана в виде системы линейных алгебраических уравнений относительно приращений расхода Д0$+1 и отметки уровня Д£к+1 по времени:

А1,0 ДС^+1 + А,1,4 Дй'+1 + д2+0 ддк+! + Д;41 д## = (6)

в,1'0 Д0?+1 + В,1,4 Д£?+1 + В+ Д<гй + В2+! д Й+1 = ^;2-Л+1,

г = 0,... N — 2.

Здесь Д1,0, Д1'4, Д-й, Ан. В,1,0 В«, В,2+0, В2;4! - известные коэффициенты. Учет

граничного условия в верхнем створе 0$ = 0(£$, ж0), к = 0,1,..., изменит два первых уравнения системы (6), Задание зависимости -1 = 0(£л-1) изменит вид двух последних уравнений системы (6) и добавит еще одно. На практике известна только дискретная зависимость О = 0(£5) для некоторых значений отметок уровня {СД^. При попадании вычисленного значения -1 в интервал [£5-1, £5] значение ОЛ— аппрок-

симируется выражением

Я%-1 + @ & С^3_ 1 + ^^— {£%-1 + 1 — Сй-О •

— ?5-1

Это дополнительное тождество на неизвестные Д£л+-11 и Д0л+-11- Выражая ДОЛ— из этого тождества и подставляя его в последние два уравнения системы (6), получаем два, новых уравнения.

2. Апробация и верификация модели

Для апробации и верификации описанной модели рассмотрен участок р, Ангары (гидро-пост Богучаны — гидропост Татарка) общей протяженностью 293 км. Натурные данные были предоставлены отделом разработки и внедрения гидрометеорологических прогнозов Красноярского центра по гидрометеорологии и мониторингу окружающей среды с региональными функциями (КЦГМС-Р) Росгидромета, Для верификации модели использовались следующие исходные данные,

1, Схематизация русла и поймы, данные по 37 характерным створам, включающие батиметрию русла и высоты поймы, значения коэффициентов шероховатости Маннинга (п), различные для русла и поймы, с учетом характеристики подстилающей поверхности (рис, 1 и 5), Данные были использованы для построения вычислительной сетки и расчета параметров уравнений. Входной створ соответствовал водомерному посту Богучаны (323 км от устья р, Ангары), замыкающий створ — водомерному посту Татарка (30 км от устья р, Ангары),

расстояние от устья р. Ангары, км

Рис. 5. Схематизация русла р. Ангара на участке гидропостов Богучаны — Татарка

Рис. 6. Динамика отметок уровня и расходов на участке гидропоста Татарка, весенне-детний сезон 2005 г.

2. Зависимость во входном створе и зависимость ), построенная по данным наблюдений для открытого русла за 2004-2005 гг. в замыкающем створе (см. рис. 3), использовались в качестве граничных условий.

3. Боковой приток р. Ангары, рассчитанный по среднесуточным расходам ее притоков, выбранных в качестве рек-аналогов (4). Общая площадь водосбора промежуточного притока р. Ангары на исследуемом участке оценена в 174 000 км2. Водосбор поделен на пять участков, соответствующих рекам-аналогам: Карабулы, Иркенеева, Каменки, Тасеева, Татарки. Следует отметить, что на водосборе расположен крупный приток р. Тасеева, дающий основной вклад в общий боковой приток (до 85%).

Тестовые расчеты проводились на основе данных наблюдений в весенне-летний сезоны 2004-2005 гг. На рис. 6 отображена динамика отметок уровня и расходов на участке водомерного поста Татарка с апреля по сентябрь 2005 г. В период весеннего паводка наблюдается три ярко выраженных максимума отметок уровня — 29 апреля, 19 и 23 мая. При этом максимумы расхода наблюдаются только в мае, более того, кривая ф(£) в апреле не соответствует зависимости, изображенной на рис. 3 и используемой в расчетах для открытого русла. Дело в том, что апрельский паводок и частично майские связаны с заторными ледовыми явлениями (о чем, в частности, есть отметка в ежегоднике). Поскольку ледовые явления и связанные с ними подъемы уровня воды не описываются в рамках тестируемой модели, достоверность расчетов следует проверять только с 24.05.2005 г., когда река на всем расчетном участке освободилась ото льда. Отметим, что и значения первых расчетных суток будут иметь несколько большую погрешность, так как на момент начала расчетов сложилась картина сильно меняющегося течения, а начальные данные рассчитываются по установившейся модели.

В табл. 1 содержится анализ погрешностей расчетов отметок уровня для тех пунктов участка, для которых имеются данные наблюдений. Это водомерные посты Богучаны (316 км от устья Ангары), Каменка (207 км от устья Ангары), Рыбное (101 км от устья Ангары) и Татарка (30 км от устья Ангары). Расчеты показывают, что максимальная ошибка не превышает 56 см для участка водомерного поста Рыбное. Эта ошибка получена в начале расчета и может быть объяснена недостаточно хорошим приближением

Таблица 1. Основные характеристики тестовых расчетов отметок уровня с 24.05.2005 г. по 30.08.2005 г. в створах наблюдения

Показатель Богучаны Каменка Рыбное Татарка

Разность максимальной и минимальной наблюденных отметок, м 1.04 1.40 2.01 2.96

Максимальная ошибка расчета отметок уровня, м 0.18 0.32 0.56 0.34

Максимальная ошибка расчета отметок уровня с 28.05.2005 г., м 0.18 0.24 0.39 0.34

Максимальная ошибка в процентах от разности уровней 17.34 23.10 28.03 11.39

Максимальная ошибка в процентах от разности уровней с 28.05.2005 г. 17.08 17.37 19.52 11.39

Средняя ошибка расчета отметок уровня, м 0.10 0.16 0.13 0.10

Средняя ошибка в процентах от разности уровней 9.61 11.53 6.22 3.24

Таблица 2. Основные характеристики тестовых расчетов расходов с 24.05.2005 г. по 30.08.2005 г. в створах наблюдения

Показатель Татарка

Разность максимальной и минимальной наблюденных расходов, м3/с 11230.0

Максимальная ошибка расчета расхода, м3/с 1177.0

Максимальная ошибка в процентах от разности расходов 10.5

Средняя ошибка расчета расхода, м3/с 318.0

Средняя ошибка в процентах от разности расходов 2.8

О

О

сл

о

(Я

о

сл

о

(Я

О

сл

о

О

сл

о

О

сл

о

га

О

сл

о

га

О

сл

о

га

О

сл

о

ш

2

к

X

со

о

СВ

о

о

сл

о

(Я

о

сл

о

(Я

о

сл

о

о

сл

дата

о о сл дата

о

га

о

сл

о

га

о

сл

о

га

о

сл

Рис. 7. Рассчитанная и наблюденная динамика отметок уровня в четырех створах: Богу-чаны, Каменка, Рыбное, Татарка (316, 207, 101, 30 км от устья Ангары соответственно)

начальных данных. Средняя же ошибка по всем сравниваемым створам не превышает 16 см (Каменка), что составляет не более 11,5 % от изменения рассчитываемой величины, В табл. 2 содержатся аналогичные данные для расхода в замыкающем створе (данных наблюдений по расходам в промежуточных створах нет).

На рис, 7 показана рассчитанная и наблюденная динамика изменения отметок уровня в створах Богучаны, Каменка, Рыбное и Татарка,

Отметим, что летний период 2005 г, отличается ярко выраженными дождевыми паводками, Лето 2004 г, — более засушливо, течение в р, Ангаре близко к установившемуся, и, соответственно, ошибки расчетов по модели минимальны (3-6 % от рассчитываемых величин). Таким образом, построенную модель можно считать пригодной для прогнозных расчетов неустановившегося течения в нижнем бьефе Богучанской ГЭС,

3. Моделирование неустановившегося движения воды в нижнем бьефе Богучанской ГЭС

При численном моделировании неустановившегося движения в нижнем бьефе Богучанской ГЭС при различных режимах ее работы был рассмотрен участок р, Ангары общей протяженностью 413 км. Входной створ расположен в 750 м от ГЭС (443,25 км от устья р, Ангары), замыкающий створ соответствует водомерному посту Татарка (413,1 км от Богучанской ГЭС), Натурные данные предоставлены отделом разработки и внедрения гидрометеорологических прогнозов КЦГМС-Р Росгидромета,

Схематизация русла и поймы проводилась по данным о 52 характерных створах. Зависимость ф(£) во входном створе задавала два сценария работы БоГЭС,

Сценарий 1 включал в себя длительный рабочий режим ГЭС с суточным колебанием расхода от 2600 м3/с в ночное время до 5175 м3/с днем, а также период катастрофически высокого попуска в условиях переполнения водохранилища Богучанской ГЭС (далее — просто “высокий попуск”) продолжительностью неделю с постепенным

3

верхние графики).

Сценарий 2 использовался для анализа прохождения волны по руслу и включал в себя длительный рабочий режим, резкое, в течение трех часов, увеличение расхода с 3

часов на рабочий режим (рис, 10, верхние графики).

Расчеты производились по каждому из сценариев в двух вариантах. Во-первых, без учета бокового притока для изучения распространения суточных волн и волн попусков. Во-вторых, с наблюденным боковым притоком в период больших дождевых паводков июня 2005 г, (см, рис, 6), Даты расчета были совмещены таким образом, что максимальный сброс Богучанской ГЭС проходил на сутки раньше максимального наблюденного суммарного бокового притока по длине реки. Для расчета бокового притока водосбор поделен на семь участков, приток на которых оценивался по притоку соответствующих рек-аналогов: Чадобца, Муры, Карабулы, Иркенеева, Каменки, Тасеева, Татарки,

На рис, 8 представлена рассчитанная динамика изменения отметок уровня для некоторых створов. По оси абсцисс отложена середина, каждых суток. На одной диаграмме приведены графики для нескольких створов, расстояние дано в километрах от Богучанской ГЭС, Для сравнения на каждом рисунке сверху отображены данные по первому створу. На рис, 8 слева отображены створы первых 120 км, что позволяет отследить

в реия, сут ки в ремя, сут ки

Рис. 8. Динамика отметок уровня, рассчитанная по сценарию 1 без бокового притока

динамику распространения суточных колебаний более подробно. Вертикальные линии синхронизируют один из суточных максимумов и прохождение середины высокого попуска. Расчеты показывают, что амплитуда суточных колебаний уменьшается с 1 м (первый створ) до практически незаметных 15 см на участке гидропоста Богучаны, а

3

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

Ниже по течению суточные колебания можно считать несущественными. На рис. 8 справа показана дальнейшая динамика распространения высокого попуска, вершина которого становится все круче, а спад волны с ростом расстояния от первого створа происходит все медленнее. Следует отметить, что максимальное значение расхода с приближением к замыкающему створу меняется несущественно. Максимальная рассчитанная отметка уровня превышает принятую для межени на 2 м, а расход — более 3

3

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

2000 4-----------1---------1---------1---------1----------1---------1---------1---------1---------1---------

12 12 12 12 12 12

время, сутки

Рис. 9. Динамика распространения кратковременного аварийного попуска. Сценарий 2, без бокового притока

время, сутки время, сутки

время, сутки

Рис. 10. Динамика распространения кратковременного аварийного попуска с боковым притоком

ниже устья р, Тасеева (370 км от Богучанской ГЭС), Максимальная рассчитанная отметка уровня в этом случае превышает принятую для межени на 3,12 м, а расход — почти на 12 000 м3/с, что существенно превышает аналогичные показатели по весенним паводкам.

Вычислительный эксперимент по распространению кратковременного высокого попуска продемонстрировал характерную картину распластывания волны с одновремен-

3

3

3

попуске превышал меженный максимально на 66 см. Скорость распространения гребня волны составляет в среднем 170 км/сут.

На рис, 10 отображена динамика изменения отметок уровня (вверху) и расходов (внизу) при кратковременном высоком попуске с учетом бокового притока. Как уже обсуждалось выше, боковой приток оказывает существенное влияние только на последние 40-50 расчетных километров в связи с впадением р, Тасеева, Вычислительный эксперимент по распространению короткого попуска в этом случае дал максимальное превыше-

3

Эти показатели меньше наблюдаемых во время прохождения весенних паводков. Список литературы

[1] Грушевский М.С. Неустановившееся движение воды в реках и каналах. Л.: Гидрометео-издат, 1982. 288 с.

[2] Годунов С.К., Рябенький B.C. Разностные схемы. М.: Наука, 1977. 439 с.

[3] Fread D.L. Numerical Properties of Implicit Four-Point Finite Difference Equations of

Unsteady Flow: NOAA Technical Memorandum NWS HYDRO 18, US Department of

Commerce, NOAA, NWS, Silver Spring, MD. USA. 1974. 88 p.

[4] Liggett J.A., Cijnge J.A. Numerical methods of solution of unsteady flow equation //

Unsteady Flow in Open Channels. Chap. 4. Water Resource Pub., Fort Collins, CO, USA.

1975.

Поступила в редакцию 14 марта 2008 г.