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

Вычислительные технологии Том 9, № 2, 2004

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

В. М. БЕЛОЛИПЕЦКИЙ, С. Н. ГЕНОВА Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: [email protected]

A mathematical model of bed silt dynamics in rivers and channels based on Bingham approximation for viscous-plastic medium is considered. Simplified equations of slow motion of a thin viscous-plastic layer are formulated for variables averaged over the cross-section of the river. An analytical solution is constructed using the quasi-stationary approximation. Estimations of thickness of the activity (moving) layer of the bed silt are obtained. A computational algorithm for calculation of the suspension and sediment dynamics, and the traction load, which takes into account both stirring-up and sedimentation, is suggested.

Введение

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

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

*Research described in this publication was made possible in part by Award No. KY-002-X1 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF). © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

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

Теория русловых процессов подробно описывается в основополагающих работах [1-9]. Для моделирования взвешенных наносов применяется диффузионная теория. Рассмотрение в полной постановке динамики русловых процессов представляет собой достаточно сложную проблему. Процесс перемещения наносов по дну имеет довольно сложный характер. Из-за непостоянства вертикальных и горизонтальных составляющих скоростей течения движение частиц имеет прерывистый характер: частицы перекатываются, останавливаются, передвигаются по дну скачками. Для оценки расхода наносов, перемещающихся путем влечения по дну, используются эмпирические формулы [2, 4]. Например, расход русловых наносов принимают в виде функции, зависящей от глубины и скорости течения.

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

рв — рв д д2 (\

гшд =----------а . (1)

д Рв 18и К ;

Здесь рв — плотность воды, рв — плотность частиц, д — ускорение свободного падения, V — коэффициент кинематической вязкости, а — диаметр частицы. Концентрация взвеси также влияет на гидравлическую крупность частиц [1]. Геометрические характеристики частиц наносов определяются их линейными размерами и формой. Обычно в расчет вводятся средние диаметры, зависящие от формы и линейных размеров частиц [7].

Распределение концентрации взвешенных наносов по глубине потока находится по формуле Рауза — Великанова, полученной из одномерного уравнения диффузии [9]:

(1 — г 5 \™д/ки.

3 = Ч—1-0 ■ (2)

где 3 — концентрация наносов на расстоянии г от дна; 3о — концентрация наносов у

дна (при г = 5); 1 — глубина потока; 5 = а/(301); а — диаметр взвешенных частиц;

к = 0.4 — постоянная Кармана; и. = 0.4\/д11 — динамическая скорость; I — уклон свободной поверхности. Для определения 30 используется выражение, предложенное Х. Эйнштейном [10]:

2.17и>д а

30

exp[0.39(ps - pB)d/pB/h] - 1

Расход наносов Qb, перемещающихся путем влечения по дну через единицу ширины створа реки, равен произведению скорости перемещения частиц ur на толщину перемещающегося слоя, который принимают равным диаметру частиц d, и на коэффициент т, учитывающий отсутствие сплошного движения всех частиц, лежащих на рассматриваемом участке створа [4]: Qb = mur d.

Множитель т называется динамическим коэффициентом сплошности, представляющим отношение объема движущихся частиц ко всему объему частиц в слое толщиной а. Скорость движения наносов, согласно экспериментальным данным, равна иг = ив — ин, где ив — средняя по глубине скорость воды, ин — предельная скорость потока, при которой отлагаются наносы, ин определяется по формуле Г.И. Шамова [4]:

ин = ба1/311/6. (3)

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

1. Одномерное приближение для определения концентраций взвешенных наносов

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

Для описания переноса взвешенных наносов в турбулентном потоке несжимаемой жидкости применяется упрощенное уравнение

дБ1 + дБ1 _ Яяі + Ч ді Мв дх к и

. д^г Чді . Ч о г л\

+ ив_т" + т; ^• (4)

Здесь Бг — концентрация г-й фракции, кг/м3; Бщ — концентрация примеси г-й фракции, поступающей с путевым притоком д; — расход взмыва-осаждения примеси г-й фрак-

ции; ш — площадь поперечного сечения русла.

1.1. Определение массообмена между дном и водным потоком

Для получения расчетных формул воспользуемся методиками работ [3, 12]. Пусть гранулометрический состав донных отложений задан параметрами агдн, где агдн — процентное содержание г-й фракции в донных отложениях, г = 1, 2,... , п. Из условия

'Шдг < IV* (5)

определяются фракции, которые взмучиваются. Здесь /шдг — гидравлическая крупность, определяемая по формуле Стокса (1) для а = йг, 'ш* = 0.4и*, и. = 0.4^/д11. Пусть взмучиваемые фракции соответствуют г = 1, 2,..., г*. Тогда процентное содержание г-й фракции

а0 = 100а0дн/г, г = а1дн + а1дн + ... + а0*дн, г =1, 2,..., г*. Средний диаметр взмучиваемых фракций dср определяется по формуле

0*

dср = 0.01> aidi. (6)

0=1

Транспортирующая способность потока Бтр зависит от средней по глубине скорости течения воды, глубины и гидравлической крупности [2, 3, 5]:

«в

о I 0.2—— при тд <^*, Ря- Рв д ,2

Ьтр = \ дК'Шд 1Пд = — -----т^-4р. (7)

0

^ Рв 18 V

при гшд > и>*, гв

Как отмечается в работе [5], существует нижний предел скорости начала подвижности, а именно при крупности наносов около 0.2 мм и мельче минимальная предельная скорость не уменьшается. На основании этого предлагается считать, что dср > 0.2 мм.

Массообмен с дном определяется по формулам, полученным с использованием методик работ [3, 12]:

0_яо (Б0тр *50о)до, Б0тр 0.01ао5Тр, 0_я ^ ^ 0_яо1 (8)

о

где дя — массообмен с дном; Бо0 — концентрация г-й фракции у дна. При вычислении Б0тр по

(8) следует учитывать, что Б0тр не может превысить концентрацию г-й фракции в донных отложениях (Б0дн), поэтому при Б0тр > Б0дн полагают Б0тр = Б0дн. Если концентрация г-й фракции в донных отложениях равна нулю, то и Б0тр = 0. Полученные формулы для вычисления массообмена уточняют соотношения работы [13].

Изменение концентраций примесей вдоль потока обусловлено как изменением глубины и средней скорости течения, так и взаимодействием водного потока с донными отложениями. Предполагаем, что основное изменение состава донных наносов происходит за счет взмучивания и седиментации. При дя > 0 происходит поступление наносов в поток (размыв), а при дя < 0 — заиление русла.

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

2. Упрощенная модель динамики донных наносов

Водонасыщенные илы представляют собой текучую вязкопластическую среду [1]. Активный слой донных наносов можно моделировать однородной несжимаемой вязкопластической средой. Такой подход применялся в работе [14] для решения задачи о течении вязкопластической среды в деформируемых каналах, в работах [15, 16] — для моделирования движения оползней. В [15, 16] скорости движения частиц среды предполагались однонаправленными (вдоль поверхности склона). В работе [14] решались уравнения движения без учета сил тяжести. Оптимизации процессов управления вязкопластическим течением в тонком слое с изменяемыми формами границ посвящена работа [17], в которой рассмотрена плоская задача о течении вязкопластичной среды в тонком слое между двумя

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

Упрощенная модель динамики донных наносов основана на следующих предположениях:

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

2. Активный слой донных наносов моделируется вязкопластической средой.

3. Движение донных наносов вызывается эффективным сдвиговым напряжением трения (после вычитания начального напряжения сдвига) и описывается уравнениями медленных слоистых течений в приближении тонкого слоя.

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

где т — касательное напряжение между соседними слоями движущейся среды; т0 — начальное напряжение сдвига; v0 — коэффициент эффективной вязкости вязкопластической среды; р — плотность донных отложений (плотность вязкопластической среды); и — горизонтальная составляющая скорости движения среды [15, 16]. В области, где |т| < т0, скорость деформации отсутствует и, следовательно, скорость не зависит от координаты z (du/dz = 0). Если т0 = const, то уравнения движения вязкопластической среды (при условии |т | > т0) совпадают с уравнениями для ньютоновской жидкости (v = v0 ).

2.1. Постановка задачи в приближении тонкого слоя

Вводится система координат (x, z), ось x направлена по ходу водотока, ось z — вниз, начало координат (х = 0, z = 0) находится на границе раздела вода — дно, угол а — уклон дна, z = n(t, х) — уравнение свободной поверхности вязкопластической среды (поверхность дна водотока); z = — H(t,x) — уравнение водной поверхности; z = Y(t,x) — уравнение нижней границы активного слоя донных отложений (рис. 1).

Медленные слоистые движения вязкопластической среды в приближении тонкого слоя описываются уравнениями

здесь p — давление; касательное напряжение т определяется по формуле (9) (в рассматриваемой задаче du/dz < 0).

(9)

(10)

Граничные условия: на границе раздела вода — дно (г = п) выполняются кинематическое условие

дП, , дП _ дя

д£ + 1г=п дх р ( )

и динамическое условие непрерывности тензора напряжений (Л = п — Н — глубина водотока)

р|^=п РвдЛ5

тп ди\ тГр

-------+ ^*0 ^ I =------при Тгр > то,

Р *=п Р

ди\

dz /

— ) =0 при ТГр < То

z=n

или

ди Г 0 при Тгр < т0,

*=, = Рв дЛ,Р^0 дг |г=П = ( —тгр + то > То. (12)

Касательное напряжение на границе раздела вода — дно тгр определяется по формуле

Т

гр рвдив /СШ,

где ив — средняя по глубине скорость течения воды, определяемая из решения уравнений Сен-Венана [11]; Сш — коэффициент Шези; рв — плотность воды. Если ив оценить по формуле Шези, то

Тгр рв ghi0- (13)

Значение i0 для малых углов вычисляется по формуле

дп дп

io = tg a + — « sin a + —.

дх дх

Из уравнения гидростатики и граничных условий (12) находятся давление

p = g(z — п) cos a + — gh P P

и уравнения для продольной составляющей скорости

дм 15т pBg dh

я! = + g sin a----я“>

dt pdz р dx (лл\

т то 5м (14)

_ =--------+ v0tt~ •

р р dz

Для нестационарной задачи задаются начальные условия

м = u0(x,z), п = П0(х). (15)

Возможны два случая:

i тгр — то ~ п

1. -------- = г > U, активный слой донных наносов движется как вязкая жидкость;

р

г < 0 — донные наносы покоятся.

2. г < 0. Вообще говоря, существует граница z = п*, на которой касательное напряжение совпадает с т0. Если п < П* < Y, то в области п < z < п* |т| < т0 и деформация отсутствует, скорость не зависит от координаты z, среда ведет себя как твердое тело. Слой вязкопластической среды толщиной п* < z < Y движется как вязкая жидкость, а слой п < z < п* движется как твердое тело со скоростью м = м(п*). В рассматриваемой модели динамики донных наносов предполагается, что реализуется только первый случай.

Рассмотрим стационарный случай. Пусть т0 = const, тогда получается задача

д2м рв dh

= -g sin a +------g—,

dz2 р dx

- тг при тг > 0,

v —I = J (16)

0 dz lz=h I 0 при Г < 0,

м|*=7 = 0 •

В стационарном случае при qs = 0 из кинематического условия (11) следует, что 5п/5х = 0, т.е. в этом случае получается задача о стекании вязкопластической среды по наклонной плоскости под действием напряжения трения водного потока.

Задача (16) имеет решение (при Г > 0)

A 2 , Ci С2

м = z2 +--z +-----,

2v0 V0 V0

1 \ Cif \ A t 2 2\ (17)

м (п) =-----(Y - п) — ^-(Y2 - п2), V 1

v0 2v0

Ci = —(Г + Ап), C2 = - A Y2 - CiY,

„ • рв dh

A = —g sin a +-------g—•

р dx

При Г < 0 м = 0.

В нестационарном случае при тг > 0 получается задача

5м 52м

sí =V0 si2 — A

du (18)

v> = -T'

u|z=Y = 0.

Положение границы раздела вода — дно определяется из уравнения (11). Для численного решения задачи (18) используются неявная схема и метод прогонки. Рассмотренную выше задачу назовем задачей 1 с известной нижней границей активного слоя.

2.2. Задача с неизвестной границей активного слоя (задача 2)

Так как положение границы активного слоя г = 7(¿,х) неизвестно, требуется дополнительное условие к граничным условиям задачи 1. Зададим условие проскальзывания на границе раздела вода — дно:

ди

дги=п = к1(«п - мвд) при Т > 0,

ип = м|г=п = 0 при Т < 0. (19)

Здесь скорость воды на дне ивд считаем известной; к1 — коэффициент проскальзывания (к1 > Т/ивд). При к1 ^ то ип ^ ивд (условие прилипания), при к1 ^ Т/ивд ип ^

0 (полное проскальзывание). Начальные условия (15) дополняются заданием положения нижней границы активного слоя:

7 = 7о(ж). (15а)

Для определения придонной скорости воды используется логарифмический профиль продольной скорости [9]:

и =___________1.25ив________ (20)

вд 1ё{6.15[0.35(1 + 1/Л)]Л> ’ ( )

где ив — средняя по сечению скорость воды; Л = (0.35А/Л,), А — шероховатость, для

однородных наносов А = 0.7д, д — средний размер частиц. Для Л./А > 50 формула (20) упрощается:

1.25ив 1ё(6Л5й/А)'

При |тгр| < т0 донные наносы не движутся, но могут осуществляться процессы взмучивания или осаждения.

Итак, если |тгр| < т0, то и = 0, если |тгр| > т0, то решается задача (11), (14), (15), (15а),

(19).

Решение стационарной задачи при Т > 0 описывается соотношениями (17), в которых 7(¿,ж) не определена. Из (17) и (19) следует уравнение для определения толщины активного слоя Т = 7 — П

Т2 — 2(Т/А)Т + 2^оич/А = 0, ич = ивд — Т/&ь

решение которого имеет вид

„ С Т/А + ^/(Т/А)2 — 2^оич/А при А < 0,

Т = < _________________ (22)

( Т/А — \/(Т/А)2 — 2^0ип/А при А > 0.

Физический смысл имеет решение, для которого 7 > п. Знак перед корнем выбран из условия ограниченности решения при А ^ 0.

Положение границы раздела п(^,ж) определяется из численного решения уравнения

(11). Определив толщину активного слоя, можно уточнить скорость движения донных

наносов на границе раздела и оценить коэффициент проскальзывания к1 путем решения задачи 1 с известной толщиной слоя 7.

Из соотношений (17) находится средняя по толщине активного слоя скорость движения донных наносов

-=к Т—АТ> (23)

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

д (г*3гдн^ д (г*3гдн)

+ “ср дх = —^ (24)

Так как = дя, ^ 3^дн = р, из уравнения (24) получается уравнение для определения г*: % %

дг* дг* Яя ,0-,

~т + “ср аХ = — 7. (25)

Концентрации взвешенных наносов определяются уравнением (4). Для уравнений (24), (25) ставятся начальные и граничные условия.

2.3. Основные параметры донных наносов

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

т0 2 Сш

--- = г —

рв

Отсюда следует

— ин — )—•

Рв у/д

иН

Т0 = Рв£т^.

Сш

Если ин определить по формуле (3) для среднего диаметра донных наносов д, Сш — по формуле Маннинга Сш = Л,1/6/п (п — шероховатость), то для нахождения т0 получается формула

т0 = 36рвдп2д2/3. (26)

Здесь используются единицы килограмм-сила, метр, секунда.

Плотность донных отложений р зависит от пористости грунта а и определяется по формуле

р = (1 — а )ря + а Рв.

Пористость грунтов изменяется в пределах от 0.30 до 0.55, ря = 2650 кг/м3, рв = 1000 кг/м3.

В работе [1] отмечается, что вязкость илистых грунтов ^0 = р^ изменяется в пределах 0.01... 0.32 кг/м-с, начальное напряжение сдвига т0 — в пределах 0.5 ... 5.5 кг/м-с2. В работе [15] приводятся данные для глинистой суспензии (концентрация глины по весу 62.2%): ^0 = 1.1, т0 = 5.74. Для сравнения т0, вычисленное по формуле (26) для д = 1.0 • 10-3, п = 0.025, равно 2.16. Требуются дополнительные исследования по установлению зависимостей коэффициентов вязкости и начальных напряжений сдвига от гранулометрического состава и влажности донных отложений.

3. Вычислительный алгоритм динамики взвешенных и донных наносов

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

Э т а п 1. Задание гранулометрического состава донных отложений в виде (дг, агдн), где дг — диаметр частиц г-й фракции, мм; агдн — процентное содержание г-й фракции в донных отложениях, г = 1, 2,...,п. В начальный момент времени концентрация г-й фракции в донных отложениях определяется по формуле

сю раадн

гдн = 100 ’

где р — плотность донных осадков.

Э т а п 2. Определение скорости течения воды ив и глубины Н из решения уравнений Сен-Венана.

Э т а п 3. Вычисление толщины активного слоя влекомых наносов и средней по толщине слоя скорости движения донных наносов по соотношениям (13), (21)—(24).

Э т а п 4. Определение массообмена между дном и водным потоком. Из условия (5) определяются фракции, которые взмучиваются, I = tg а + дН/дх. Обозначим взмучиваемые фракции индексом г = 1, 2,..., г*, аг — процентное содержание взвешиваемых фракций, аг = 0 при г = г* + 1, г* + 2, ...,п. Процентное содержание всех взвешиваемых фракций г = а1дн + а2дн + +аг+дн. Тогда процентное содержание взвешиваемой г-й фракции аг = (100/г)агдн, г = 1, 2,..., г*. Если г = 0 (нет взмучиваемых фракций), то все аг = 0.

Массообмен между дном и придонным слоем определяется по формулам (6)-(8). На основе формулы Рауза — Великанова (2) для распределений концентраций по глубине водного потока находится средняя по глубине концентрация г-й фракции £г:

5г = вг £¿0,

где

1 х

в = I (т+—к)'-' ^ г

0

Э т а п 5. Определение концентраций взвешенных и донных наносов и положения границы раздела вода — дно.

Численные решения соответствующих уравнений (4), (11), (24), (25) основаны на схеме бегущего счета:

С” + и _^Сп+1 і і Л^ С»

Сп+1 = ^ + иВ Лж ’-7-1 + ^ + Ц ^

/ Л А ’

1 + ив-;— '

Ax J

z" I u At zra+l Atq»

z*j + Ucp AxZ*i-i p qSj

zn+1 _

z*j / At

1 + Ucp AX

Z* Sn I u At çra+1 At

zn+1 °адн j + Ucp a х°гдн,і-1 zn+1 У^г,і

S"+1 _____ AX Z*

tflHj /_ At

’ Ax

1 + ucp

nn і un ___________nn+1 +_______qn

Чдн,і + un,j ДхПДн,і-1 + p qSj

n+1 = ________________AX_____________________p

да,і = / A A :

1 + и"> AX )

где Ai — шаг по времени; Ax — шаг по пространственной переменной; индекс j соответствует координате Xj, Xj+1 = Xj + Ax; tn+1 = tn + At, /j = / (tn,Xj).

Э т а п 6. Расчет гранулометрического состава донных отложений:

S n+1

«гдн = -^100.

ідн

Р

Э т а п 7. Вычисление расходов взвешенных и влекомых наносов. Расход взвешенных наносов в створе х = хз- определяется по формуле

Qw = QSj, Sj = £ S"/1.

ул '-'з / ; г,3

г=1

Расход влекомых наносов вычисляется по формуле

^5 7Висрр-

Здесь ^ — расход воды, В — ширина русла. Общий расход наносов : + ^ь.

На следующем временном интервале расчеты повторяются со второго этапа по седьмой.

4. Примеры расчетов

4.1. Модельное русло, схематизирующее участок реки в районе г. Енисейска

Подробные данные о мутности воды и гранулометрическом составе взвешенных и донных наносов р. Енисей получены в нескольких створах. Верификация разработанной вычислительной методики проводилась для конкретного участка — г. Енисейска. Исходные данные

В. М. Белолипецкий, С. Н. Генова Таблица 1

1 2 3 4 5 6 7 8 9

0.0008 0.003 0.0075 0.03 0.075 0.15 0.35 0.75 1.5

$ 2.2 1.3 1.2 7.2 2.4 5 0 0 0

Sгдн 20 14 8 16 2 28 160 90 20

<н 1.0 0.7 0.4 0.8 0.1 1.4 8.0 4.5 1.0

агдн 0.003 0.0001 0.0001 0.0007 0.0002 0.0005 8.3706 4.7061 1.0463

10 11 12 13 14 15

3.5 7.5 15 35 75 120

0 0 0 0 0 0

^гдн 26 216 200 400 400 400

а0 гдн 1.3 10.8 10.0 20.0 20.0 20.0

агдн 1.3579 11.2787 10.4504 20.9008 20.9008 20.9008

Таблица 2

1 2 3 4 5 6 7 8

А 0.0008 0.003 0.0075 0.03 0.075 0.15 0.35 15

Sг 2.2 1.3 1.2 7.2 2.4 5 0 0 0

Sгдн 20 14 8 16 2 28 160 1752

а0дн 1.0 0.7 0.4 0.8 0.1 1.4 8.0 96.6

агдн 0.003 0.0001 0.0001 0.0007 0.0002 0.0005 8.3706 91.4766

для расчетов взяты из работы [21]. Гранулометрический состав донных отложений и их процентное содержание в начальный момент времени а°дн и в конце расчета агдн приводятся в табл. 1 (май, 1984 г.).

В результате вычислительных экспериментов было выяснено, что можно рассчитывать перенос и осаждение примесей не для всего гранулометрического состава (15 фракций), а только для части фракций, наиболее полно характеризующих все донные отложения и участвующих в процессе взмыва-осаждения, при этом процентное содержание мелких фракций при расчете на длительный период времени для 15 и 8 фракций одинаково,

поэтому крупные фракции можно объединить ( а(ДН = ^ а(ДН, табл. 2 ). Верхние индексы

в (8) и (15) показывают отношение данной величины соответственно к расчетам с 8 и 15 фракциями.

Рассматривался участок реки с постоянными вдоль русла глубиной, шириной и скоростью течения, (данные, схематизирующие реальное русло в районе г. Енисейска). По данным [21] для расхода Q = 14 500 м3/с мутность воды в районе г. Енисейска ^ = 17 г/м3.

1

В результате расчетов была получена мутность ^ 19 г/м3.

4.2. Модельное русло, схематизирующее протоку Балчуговскую

Рассматривается упрощенная схема течений в русле реки при раздвоении, когда имеются главное русло и протока (рис. 2).

Полный расход воды Q равен сумме расходов отдельных его частей:

(27)

где Qр — расход главного русла; Qп — расход протоки.

Для определения расхода воды и скорости используется формула Шези [22]:

Q = иСШ^К/, V = ^ (28)

и

или

Q = XV7/, (29)

где

К = иС^ (30)

— модуль расхода; I — уклон свободной поверхности; коэффициент Шези определяется по формуле Маннинга. Для широких русел гидравлический радиус К приближенно равен глубине потока (К ~ Л).

Формула (29) применяется для определения уклона свободной поверхности протоки и главного русла:

Q2 ^

I = I = — (31)

/р = К2 ’ /п = К2 • (31)

рп

Здесь индекс “р” относится к параметрам главного русла, индекс “п” — к параметрам протоки. Предполагается, что

1р = 1п = I,

тогда из (27), (31) следует

^ 1 + Кр/Кп • (32)

Таким образом, для заданных ип, Лп, ип, Лп из соотношений (27), (28), (32) определяются Qр , Vр , Qп , Уп.

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

Для верификации модели были проведены расчеты для протоки со следующими исходными данными: Q = 3300 м3/с, Лр = 7 м, ир = 3500 м2. Протока состоит из двух участков длиной 3 и 2.7 км: ЛП = 3 м, иП = 510 м2, ЛП = 1.2 м, иП = 900 м2 Результаты расчетов: Ур = 0.9 м/с, У1 = 0.32 м/с, У;2 = 0.18 м/с.

По найденным скоростям в протоке при заданных гранулометрическом составе и мутности в начальном створе были рассчитаны параметры взвешенных и донных наносов. В табл. 3 приведены исходные величины (гранулометрический состав — данные Ф.В. Су-хорукова [Интеграционный проект СО РАН, № 75]) и результаты расчетов.

N

Параметр 1 2 3 4 5 6 7 8 9 10

мм 0.0008 0.003 0.0075 0.03 0.075 0.15 0.35 0.75 1.5 3.5

3 £ Си ^0 8 6.8 5.4 3.8 1 0 0 0 0 0

^гдн г/м" 105.2 73.6 42 84.2 10.4 147.4 842 474 105.2 116

а0 гдн 5.26 3.68 2.1 4.21 0.52 7.37 42.1 23.7 5.26 5.8

агдн 5.34 3.74 2.26 4.10 0.66 5.95 42.72 24.02 5.34 5.88

Мутность воды практически не изменилась. Расчетный гранулометрический состав согласуется с натурными данными.

4.3. Расчеты динамики взвешенных и влекомых наносов

на участке р. Енисей от г. Красноярска до г. Енисейска

По разработанной модели были сделаны расчеты для реального русла реки Енисей. Гидравлические характеристики для задачи транспорта наносов (глубина, скорость течения) получены из решения одномерной задачи Сен-Венана. Исходные данные по мутности и гранулометрическому составу представлены в табл. 4.

Русло реки было разделено на главное русло и пойму. В главном русле большая глубина и более высокие скорости течения воды. В этом случае мелкие частицы взмучиваются и на дне увеличивается процентное содержание крупных частиц. Мутность воды на разных участках колеблется от 25 до 100 г/м3.

Проведены расчеты для поймы реки. Скорости в прибрежной зоне в несколько раз меньше, чем в основном русле. На рис. 4 показаны скорости течения и мутность воды вдоль потока, на рис. 5 — приведен расход взмыва осаждения примесей.

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

Между поступлением наносов и транспортирующей способностью реки наблюдается соответствие. Твердые частицы собираются в реке с более или менее обширных площа-

Таблица 4

Параметр 1 2 3 4 5 6 7 8

(1г, мм 0.0008 0.003 0.0075 0.03 0.075 0.15 0.35 15

3 /м г/ ^0 7.5 4 3.3 4.5 4.6 1.2 0 0

^гдн г/м" 4 8 8 10 10 28 160 1772

а0 гдн 0.2 0.4 0.4 0.5 0.5 1.4 8.0 88.6

Рис. 3.

Таблица 3

Рис. 4.

Рис. 5.

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

В фарватере русла взвешенные наносы проносятся до конца расчетного участка (г. Енисейск). В местах осаждения (на поймах и в протоках) происходит седиментация взвешенных наносов. Численными экспериментами определены прибрежные зоны с преимущественным осаждением наносов, являющиеся вероятными участками аккумуляции загрязненных наносов в донных отложениях. Осаждение происходит в протоках и прилегающих к берегам областях, отстоящих от ГЭС на расстояния 135, 150, 159, 180, 220, 226, 229, 247, 250, 261, 328, 351, 358, 365, 389, 422 457, 460 км.

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

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

[1] Боровков В.С. Русловые процессы и динамика речных потоков на урбанизированных территориях. Л.: Гидрометеоиздат. 1989. 286 с.

[2] ГришАнин К.В. Динамика русловых потоков. Л.: Гидрометеоиздат, 1979. 211 с.

[3] Караушев А.В. Методические основы оценки и регламентирования антропогенного влияния на качество поверхностных вод. Л.: Гидрометеоиздат, 1987. 285 с.

[4] Чеботарев А.И. Общая гидрология. Л.: Гидрометеоиздат, 1975. 544 с.

[5] АбАльянц С.Х. Устойчивые и переходные режимы в искусственных руслах. Л.: Гидрометеоиздат, 1981. 238 с.

[6] Барышников Н.Б. Морфология, гидрология и гидравлика пойм. Л.: Гидрометеоиздат, 1984. 280 с.

[7] ДинАмикА русловых потоков и литодинамика прибрежной зоны моря. М.: Наука, 1994. 304 с.

[8] Мирцхулава Ц.Е. Основы физики и механики эрозии русел. Л.: Гидрометеоиздат, 1988. 303 с.

[9] Россинский К.И., ДЕБольский В.К. Речные наносы. М.: Наука, 1980. 218 с.

[10] Einstein H.A. Formulas for the transportation of bedload // Prog. ASCE. 1942. N 3.

[11] Численное моделирование задач гидроледотермики водотоков / В.М. Белолипецкий, С.Н. Генова, В.Б. Туговиков, Ю.И. Шокин. Новосибирск: СО РАН, ИВТ, ВЦ (г. Красноярск), 1994. 136 с.

[12] Адесман А.В. Основные уравнения диффузионной модели общих русловых деформаций // Динамика и термика рек, водохранилищ, внутренних и окраинных морей: Тез. докл. четвертой конф. 1994. Т. 1. М. С. 171-173.

[13] Белолипецкий В.М., Генова С.Н., Гуревич К.Ю., Дегерменджи А.Г. и др. Компьютерная система для исследования динамики гидрофизических и радиоэкологических характеристик речной системы // Вычисл. технологии. 2001. Т. 6, № 2.

С. 14-24.

[14] ЧЕсноков В.М. Пространственные тонкослойные течения вязкопластичных сред в деформируемых каналах и полостях // МЖГ. 1998. № 4. С. 201-205.

[15] САгомонян А.Я. К вопросу дождевой эрозии почв // Вестн. МГУ. Сер. Математика, механика. 1995. № 5. С. 85-94.

[16] САгомонян А.Я. Движение оползней, возникающих на склонах возвышенностей под действием дождя // Изв. РАН. Механика твердого тела. 1998. № 6. С. 143-148.

[17] Петров А. Г. Об оптимизации процессов управления вязкопластическим течением в тонком слое с изменяемыми формами границ // Изв. РАН. Механика твердого тела. 1997. № 2. С. 127-132.

[18] МилитЕЕв А.Н., Базаров Д.Р. Математическая модель для расчета двумерных (в плане) деформаций русел // Водные ресурсы. 1999. Т. 26, № 1. С. 22-26.

[19] Meakin P., Sun T., Jossang T., Schwarz K. A simulation model for meandering rivers and their associated sedimentary environments // Physica. A. 1996. Vol. 233, N 3-

4. P. 606-618.

[20] Петров А.Г., Петров П.Г. Вектор расхода наносов в турбулентном потоке над размываемым дном // ПМТФ. 2000. Т. 41, № 2. С. 102-112.

[21] Ежегодные данные о режиме и ресурсах поверхностных вод суши. 1984. Ч. 1 и 2. Т. 1. Вып. 12. Обнинск. ВНИИГМИ — МЦД. 1986. 379 с.

[22] Спицын И.П., Соколова В.А. Общая и речная гидравлика. Л.: Гидрометеоиздат, 1990. 359 с.

Поступила в редакцию 11 июля 2003 г., в переработанном виде — 9 января 2004 г.