Математическая модель для исследования устойчивости русловых процессов Текст научной статьи по специальности «Физика»

УДК 551.465.553

В.А. Шлычков, Л.Ф. Дзюбенко

ИВЭП СО РАН, Новосибирский филиал; НГАСУ, Новосибирск

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДЛЯ ИССЛЕДОВАНИЯ УСТОЙЧИВОСТИ РУСЛОВЫХ ПРОЦЕССОВ

В настоящее время многие реки оказались подвержены влиянию естественных климатических и антропогенных факторов. Это обстоятельство наряду со значительным (до 40%) сокращением числа наблюдательных постов обусловливает необходимость проведения теоретических исследований по устойчивости русловых потоков. Использование традиционных подходов, основанных на инженерных формулах и методах, не позволяет получать надежных предсказаний русловых деформаций.

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

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

В основу численной модели плановых течений положены фундаментальные законы механики жидкости, выраженные в виде двумерных уравнений мелкой воды. Введем декартову систему координат с осями х,у^, в которой ось z направлена вверх. Пусть поверхность водосбора (включая русловые углубления) задана уравнением г = 5( х, у) + \хх + \уу, где , I - малые уклоны в направлениях х,у, характеризующие понижение оси

главного русла вдоль водотока.

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

/| h ^ h

скорости U = — J udz, v = — J vdz, где u,v - компоненты вектора

8 8

горизонтальной скорости вдоль осей x,y; h - глубина потока. Отбрасывая черту у переменных, систему уравнений плановых течений запишем в виде

д hu д huu д huv .. ,д (h + 8) д ... д и д ... д и .

+---------+---------= -ghix -gh—------ +—hKx— +—hKv — + г„ -Л\u\u

д t д x д y x д x д x x д x д y y д y x u 1

д hv д huv д hvv .. , д (h + 8) д ... д v д ... д v .

+---------+--------= -ghiv - gh—------ + — hKx — + — hKv — + ^ -Л.. \u\v

д t д x д y y д y д x x д x д y y д y y u 1

д h дuh дvh

+---------+---= Ra -1, (1)

д t д x д ya

где тх,ту - напряжения ветра, Xu - коэффициент гидравлического трения, \й\ - модуль скорости течения, Kx, Ky - коэффициенты турбулентности, Ra -

интенсивность источников (осадки), I - инфильтрационный сток и испарение.

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

un = 0, K —^ = -Cu | u\us при (x,y) e Г, (2)

д n

где Г - граница области, n - внешняя нормаль к ней, un, us - нормальная и касательная составляющие вектора горизонтальной скорости, cu -коэффициент сопротивления вдоль линии берега. При наличии открытых (жидких) границ должны быть известны расходы , Q(y) в заданных сечениях

hu = Q(x), hu = Q(y) при (x, y) еГ, (3)

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

д u д v

hK — = hK — = 0 при (x,y) еГ. (4)

д n д n

Моделирование турбулентного обмена проводится на основе двухпараметрической системы уравнений полуэмпирической теории турбулентности [Монин, Яглом, 1992], включающие осредненные по z уравнения для плотности кинетической энергии турбулентности e и скорости ее диссипации в.

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

скорости ис, которая зависит от диаметра частиц и сил внутреннего трения грунта в воде.

Уравнения переноса взвешенных фракций в речной воде формулируются

в рамках плановой модели (1), описывающей динамику двумерного

вертикально осредненного потока. Вектор горизонтальной скорости и

= (и,ч) здесь считается известной функцией переменных х,у и времени £

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

(отношение объема твердых взвешенных частиц к общему объему смеси), с±

- частная мутность отдельной фракции. При этом считается с с,.

/'

Дифференциальное уравнение переноса концентрации получим из общего уравнения турбулентной диффузии твердой примеси [Караушев, 1977] путем его осреднения по вертикали

8Ьс 8Ьис 8hvc 8 , 8с 8 иг. 8с ... . .. .

~Я4~ + ^---+ ^----= ^“+ У^~ + г^(с* -с) -(1~п) Мдс + Кс + ,

81 8 х 8 у 8 х 8 х 8 у 8 у

(5)

где Ъ.(х,у,Ь) - глубина потока, Ех, Ех, - коэффициенты дисперсии в направлениях х,у, шд - гидравлическая крупность (скорость седиментации), зависящая от рассматриваемого компонента полидисперсной примеси, Ис - интенсивность внешних поступлений примеси, ^ - поток взвеси за счет поверхностных волн, величина с характеризует равновесное значение концентрации, N - коэффициент обмена между примесью и потоком в придонном слое, п - параметр, отражающий интенсивность обмена.

Величина г < 1 имеет смысл вероятности превышения мгновенной скоростью неразмывающего значения и и рассчитывается по формуле [Россинский, Дебольский, 1980]

<£-иь Г

г = —\т= |ехр| - иЬ * , (6)

2^ ис I

и

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

согласно (5) при Ис=И^=0 будет равным

г N с,

. В перегруженном

гN+(1 - г)мд

потоке с мутностью больше критической с > с,, процессы донного взмыва блокируются, и источник в правой части (5) становится отрицательным, что соответствует выведению наносов из суспензии при любом г. Представленный здесь вероятностный подход к описанию взаимодействия взвешенных и донных наносов путем введения г по (6) является развитием постановки [Атавин и др., 1995].

Описание обмена наносами проведем согласно работе [Беликов и др.,

2003], полагая в (5) N = и., где и. =| и | - скорость трения, С - параметр

C

Шези. Придонная неразмывающая скорость определяется по формуле

где d -размер частиц рассматриваемой фракции.

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

о 5 10 15 0 5 10 15

а) 6)

Рис. 1. Конфигурация русла Новосибирского водохранилища в районе с. Новопичугово и структура течения:

а) поле функции тока (стрелки показывают направление движения частиц);

б) изолинии модуля скорости (цифры показывают скорость потока в см/с).

Символом ® обозначено местоположение карьерной разработки

Ширина русла в исследуемом районе водохранилища составляет 6-8 км. На рис.1,а показана структура функции тока на русловом участке водохранилища. На входном створе русла (верхняя граница рисунка)

3 1

задавался полный расход <2=1200 м с- , соответствующий сезону осенней межени. Линии тока в стационарном течении отражают траектории части, и, вместе с тем, образуют трубки тока, в которых значение расхода воды постоянно (на рис.1,а трубки тока выделяют струи с равными расходами 100

3 1

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

Это проявляется и в поле скоростей, показанном на рис. 1,б. Расчетное поле скоростей, полученных в эксперименте, характеризовалось сравнительно небольшими значениями до 10 см/с. При таких скоростях становится ощутимым влияние поверхностных ветровых напряжений. Так, при юго-западном ветре (вдоль течения) силой 5 м/с скорость потока увеличивается вдвое, что соответствует данным наблюдений.

Анализ поля скоростей показывает, что превышение пульсационным компонентом неразмывающих значений невелико (^ < 0.02) и эрозии русла в данных условиях не происходит. Основным механизмом русловой деформации является аккумуляция донных отложений за счет взвеси, переносимой в толще потока и выпадающей в областях, где скорость потока меньше "незаиляющего" значения. Эти области расположены на мелководьях правобережья в районе карьерных работ. Тем самым, разработка дна вносит риск заиления правобережной части русла с сопутствующими русловыми деформациями.

Работа выполнена при поддержке РФФИ, проект № 05-05-98012.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Атавин А.А., Жданов Е.П., Копылов Ю.Н. Одномерная продольная модель переноса тепла и взвешенных речных наносов в водохранилище. Метеорология и гидрология, 1995, № 2, с.93-101.

2. Беликов В.В., Прудовский А.М., Середавкин В.П. К определению оптимальных параметров отстойников ГЭС. Сб. "Безопасность энергетических сооружений", вып. 12, 2003, с.166-183.

3. Караушев А.В. Теория и методы расчета речных наносов. Л., Гидрометеоиздат, 1977. 272 с.

4. Монин А.С., Яглом А.М. Статистическая гидромеханика. Т.1. СПб, Гидрометеоиздат. 1992. 694 с.

5. Россинский К.И., Дебольский В.К. Речные наносы. М., Наука, 1980. 214 с.

© В.А. Шлычков, Л.Ф. Дзюбенко, 2006