CC BY
19
Онлайн-доступ к журналу:
Серия «Науки о Земле»
2008. Т. 1, № 1. С. 33-42
Иркутского государственного университета
И З В Е С Т И Я
УДК 556
Численное моделирование течений и переноса примесей в водотоках Байкальского региона
Аргучинцев В. К.(аг§@ша1;кл8и.ги), Аргучинцева А. В., Убонова Л. В.
Аннотация. Предложена математическая модель для расчета скоростного поля водотоков и распространения загрязняющих ингредиентов после сброса их в определенном участке реки. Приводятся расчеты для конкретных объектов.
Ключевые слова: модель, водоток, поле скоростей, загрязнение. Введение
Гидрологические обоснования проектных решений, оценка состояния водных объектов, расчет штрафных показателей за сбросы, прогноз характеристик на какой-то отрезок времени проводятся, как правило, по стандартным (юстированным) методикам, которые включают большое количество констант (параметров), не являющихся характерными для конкретного объекта.
Для изучения мезомасштабных процессов в водотоках и водоемах сеть гидрологических наблюдений либо слишком редка, либо вообще отсутствует для интересуемых водных объектов. Кроме того, необходимо отметить, что даже если наблюдения и имеются, то они относятся только к определенным пространственным точкам, а, следовательно, не могут в достаточной степени отражать все особенности изучаемого объекта, например поле скорости при обтекании препятствий (острова, искусственные сооружения).
Теоретическая постановка проблемы
Одним из перспективных методов изучения водных объектов является математическое моделирование, которое позволяет объяснять с теоретической точки зрения количественные и качественные закономерности, решать диагностические и прогностические задачи.
В данной работе расчеты поля скоростей и распространения бытовых и ливневых сбросов проведены по авторской модели, которая является более универсальной для описания мезомасштабных процессов в водоемах [1, 2]. Гидрологические характеристики находятся на основе решения уравнений теории мелкой воды [1, 2, 3] с параметризацией влияния силы трения о дно и учетом турбулентного обмена по горизонтали:
du du du d(h + 8) . д , du д , du ru\v\
— + u— + v— = -g—-- + lv +—kx — +—к---^,
dt дх ду дх дх дх ду dy h
dv дv дv d(h + 8) л d dv d , дv rvV
— + u— + v— = -g—-- - lu +—kx — +—к---—,
дt дх ду ду дх дх ду ду h
дh дuh дvh
— +-+-= 0.
dt дх ду
Здесь t - время; u,v - компоненты вектора скорости движения воды вдоль осей правой декартовой прямоугольной системы координат х, у ;
|v| = Vu2 + v2 ; g - ускорение свободного падения; h - глубина водоема; l = 2с sin р - параметр Кориолиса; с - угловая скорость суточного вращения Земли; р - широта; кх, ку - коэффициенты турбулентного обмена
соответственно вдоль осей координат хи y; 8(х,у) - функция, описывающая рельеф дна; r - коэффициент придонного трения.
Система дифференциальных уравнений в частных производных (1) решена при начальных условиях:
^хуО) = п)(х,у), у(х,у,0) = у0(х,у), Н(х,у,0) = И°(х,у).
Контур области интегрирования состоит из твердой части и открытых границ. На твердой части контура задается условие прилипания, на открытых границах - значения искомых функций или их производных в зависимости от направления скорости течения.
С помощью уравнения неразрывности система уравнений (1) преобразуется к следующему виду:
ди 1 ди 1 duU 1 ди 1 dvU
-+ — u-+--+ — v-+--=
dt 2 дх 2 дх 2 ду 2 ду
TTd(h + S) d , ди d , ди
- g^—-+ IV + — кх-+ — ку-
дх дх дх ду ду h
и|,
dv 1 dv 1 duV 1 dv 1 dvV
-+ — u-+--+ — v-+--
dt 2 дх 2 дх 2 ду 2 ду
TTd(h + 8) ЙТТ d , dV d , dV rV\v\
- gH—-- - lU + — кх-+ — к---L-1
ду дх дх ду ду h
(2)
dh dHU dHV n
— +-+-= 0,
dt дх ду
где H = 4h ; U = Hu, V = Hv
Интегрирование системы уравнений (2) проведено в декартовой системе координат с применением метода фиктивных областей, позволяющего учитывать произвольный рельеф дна водоемов. Численный алгоритм решения задачи построен на основе метода расщепления по физическим процессам и геометрическим переменным [4].
Решение задачи на каждом временном шаге осуществляется в два этапа. На первом этапе решается следующая система уравнений:
ди 1 ди 1 дии 1 ди 1 дуи д , ди д , ди
-+—и-+--+ — V-+--= — кх-+—к -,
дТ 2 дх 2 дх 2 ду 2 ду дх дх ду ду
д¥ 1 дУ 1 диУ 1 дУ 1 дvV д , дУ д , дУ
—+—и—+--+ - V—+--= — кх —+—ку—.
д 2 дх 2 дх 2 ду 2 ду дх дх ду ду
На этом этапе для каждой искомой функции рассматривается эволюционное уравнение:
дш 2
+ Ьш = 0, где Ь = XЬт .
ОТ т=1
Для дискретизации по времени используется схема Кранка-Николсона и двуциклический метод многокомпонентного расщепления [4], который
состоит в разложении сеточного оператора Ь > 0 на более простые операторы Гт > 0 . Операторы Гт > 0 аппроксимируем со вторым порядком точности по координатам.
Введем неравномерную сетку с основными узловыми точками х{ = /Ах
(/ = 0,1,к,I +1); у3 = ,Ау (, = 0,1,к,J +1); Тп = пА (п = 0,1,к) и шагами сетки Ах , Ау , АТ. Будем также использовать вспомогательные точки х/+1/2, у,+1/2, расположенные в серединах основных интервалов. Обозначим
Ш1],к = Ш(хг,Уj ,Тп) ; Ч+1/2,},к = (и/ +1,},к + Ч,},к)/ 2 ;
V.
/, j+1/2,к
2,к = (У,,+1,к + V,},к )/2.
Приведем разностные аналоги операторов:
ы, =
иГ+1/2,;Ш+1,; -ЧП-1/2,]Щ-1,] к"г+1/2,] ш) ^Н/ЦШ )
Ы; =
2Ах Ах2
УП, +1/2Ш+1 -У, -иМн _ ^,,+1/2 Ш,,+1 )-кП/,;-1/2 к ш) 2Ау Ау ■
Тогда алгоритм расщепления имеет вид:
(Е + )уп-1/2 =[ ЕУ-1,
(ЕУ =( Е-А-ц У"-1/2,
( Е + У+1/2 =( Е-^122 У, ( Е + = ( ЕУ+1/2.
Для повышения точности расчетов здесь использована двуциклическая перестановка этапов расщепления. Для численной реализации конечно-разностных уравнений используется немонотонная прогонка [5]. На втором этапе система уравнений имеет следующий вид:
и=- н д(И+8) + уу - гШ,
дХ дх И
ш-^, (3)
дх ду И
дИ дни дну п
— +-+-= 0.
дх дх ду
Систему уравнений (3) аппроксимируем неявной разностной схемой первого порядка точности по времени:
ш+2 -^Т"+1 И22 - И22 +8 -8 ги"+2\у"\
иг,1 иг,1 _ = 1/2,1 "1-1/2,] +01+1/2,] °1-1/2,] + у+2 ', 1 ¡ ^А
А X " 5,1 Ах 1,] ,Ч]
у22 -уп+1 и22 _И22 +8 -8 гУ+
уг,1 Уг,1 _ = _Н2 ,1+1/2 "¡,1-1/2 + 0',]+1/2 и',}А/2 -ццп+2 ',1
АХ ,1 Ау
н у+2 -нп У+2 ......"....." 0
Ип+2 - ип+1 н ТГ+2 - н ТГ+2 н у+2 - н2 у+2
'н,! ",! , '+1/ 2,1 '+1 /2,1 '-1/2 '-1/2,1 | ', 1+1/2 ¡1+1/2 пг,1-1/2уг, 1-1/2
АХ Ах Ау
Подставляя компоненты скорости и и У из первых двух уравнений в последнее, получаем уравнение для Ь, которое решаем методом матричной факторизации [6]. После решения уравнения для И рассчитываются и и У.
Построенные конечно-разностные схемы абсолютно устойчивы, имеют первый порядок аппроксимации по времени и второй - по координатам.
Для решения экологических задач, связанных с охраной рек от загрязнения, используется полуэмпирическое уравнение переноса и турбулентной диффузии примеси.
Для моделирования распространения субстанций С в водоеме рассмотрим уравнение переноса и диффузии пассивной примеси [3].
дС дС дС д , дС д , дС ^
— + и— + V— = — кх — + — ку — + Г. (4)
дТ дх ду дх дх ду ду
Из-за отсутствия детальной информации из наблюдений в качестве начальных условий принимается С, равное фоновому распределению. На границах поставлены условия второго рода. При этом потоки примеси через твердые границы водоема предполагались отсутствующими.
С использованием уравнения неразрывности уравнение (7) преобразовано к следующему виду:
дя 1 дя 1 дия 1 дя 1 дта д , дя д , дя
— +—и— +--+ — V— +--= — кх — +—ку — + /, (5)
дТ 2 дх 2 дх 2 Оу 2 ду дх х дх ду у ду
где Б = СИ, / = Г Б .
Для дискретизации по времени использована схема Кранка-Николсона и двуциклический метод многокомпонентного расщепления ([4].
Объекты исследования
Для демонстрации возможностей модели приводятся фрагменты отдельных расчетов для некоторых водных объектов: р. Селенги в районе сброса целлюлозно-картонного комбината, р. Туул (Монголия, г. Улан-Батор), р. Залари (Иркутская область).
Результаты моделирования
Эксперимент 1. На рис. 1 приводится расчет скоростного поля р. Селенги на участке сброса Селенгинского целлюлозно-картонного комбината. Расчеты велись для сетки размером 100 х 70 точек с шагом 25 м. Для получения необходимой точности расчеты проводились с шагом по времени, удовлетворяющим критерию Куранта-Фридрихса-Леви относительно скорости движения.
Эксперимент 2. На рис. 2 и 3 приведены скоростные поля потока для двух участков р. Туул в районе г. Улан-Батора. На этих участках русло реки делится островами. Поле скорости использовалось для расчета загрязнения реки. На рис. 4, 5 приведены изолинии концентраций загрязняющих веществ в процентах от их значений на втоке. Из рис. 4, 5 видно, что примесь распространяется по течению реки и наибольшие концентрации загрязняющих веществ отмечаются у берегов и вокруг островов.
Рис. 1. Схема поверхностных течений участка р. Селенги
Рис. 2. Схема поверхностных течений 1-го участка р. Туул
() 1.5 м/с О 1ЧМ1 м
Рис. 3. Схема поверхностных течений 2-го участка р. Туул
I-1
Рис. 4. Поле загрязнения воды 1-го участка р. Туул
Эксперимент 3. Модель может быть применена и для решения различных экологических задач, требующих знания детальных гидрологических характеристик. Так, по предложенной математической модели выполнены конкретные расчеты для отдельных участков р. Залари в районе пос. Ново-нукутский Иркутской области. Необходимость в расчетах возникла в связи с намечаемым строительством в указанном районе нового гипсокартонного завода и потенциальным загрязнением реки его сбросами.
Для измерения глубин (сентябрь 2007 г.) использована гидрометрическая штанга длиной 1,5 м. Промерные точки назначались через 1 м. Положение точек определялось по натянутому размеченному тросу.
Рис. 5. Поле загрязнения воды 2-го участка р. Туул
Численные эксперименты проведены на расчетной сетке 102 х 102 точек с разрешающим шагом по горизонтали 4,7 м. Граничные условия на втоке для расчета скоростей потока были заданы на основе данных натурных измерений в створах. Схема потоков по модельным расчетам представлена на рис. 6, 7. На основе найденных по модели скоростей потока были рассчитаны поля распространения концентраций ингредиентов при ливневых (рис. 6) и хозяйственно-бытовых стоках (рис. 7). Необходимо отметить, что в модели для общности расчетов сброс загрязняющих веществ (ливневых и хозяйственно-бытовых) принимался равным 1, изолинии проведены в процентах от планируемого сброса. Изолиния 1 соответствует 1 % от сброса. Далее изолинии проведены с шагом 5 %. Поэтому для конкретного вида стока и конкретной примеси легко рассчитать ее концентрацию на соответствующей изолинии.
С А
0 100 200 м
Рис. 6. Схема течений и распределение примеси в районе сброса ливневых стоков
С
4
0 100 200 м
Рис. 7. Схема течений и распределение примеси в районе сброса хозяйственно-бытовых стоков
Проведенные исследования позволили нам сделать следующие выводы:
- извилистость русла на участке расположения выпусков незначительная - 1,03;
- в верхнем бьефе дамбы происходит аккумуляция наносов с уменьшением глубины потока и увеличением скоростей течения;
- скорость размыва берега в створе выпуска хозбытовых стоков составляет 1,7 м/год;
- в результате размыва русла ниже дамбы и моста уклон потока равен 7,6 м/км, что соответствует горной, а не равнинной реке.
Выводы
Предложенная модель может быть использована для расчета скоростного поля естественных водоемов и водотоков произвольной конфигурации, а также для изучения распространения загрязняющих веществ, попадающих в результате их сброса промышленными, бытовыми объектами или за счет ската ливневых потоков. Такая информация необходима местным властям и населению для решения вопросов забора питьевой воды, воды хозяйственного назначениия, организации рыборазведения и отлова. Расчетную информацию можно использовать и для решения ряда прикладных задач, возникающих, например, с позиций безопасного судоходства на реках в условиях разветвленного русла, небольших глубин и значительных колебаний уровня, а также в плане строительства и эксплуатации береговых сооружений.
Список литературы
1. Аргучинцев В. К. Численное моделирование гидрологических характеристик и процессов распространения примесей в реках / В. К. Аргучинцев, А. В. Ар-гучинцева // Докл. РАН. - 2000. - Т. 370, № 6. - С. 803-806.
2. Аргучинцев В. К. Моделирование мезомасштабных гидротермодинамических процессов и переноса антропогенных примесей в атмосфере и гидросфере региона оз. Байкал / В. К. Аргучинцев, А. В. Аргучинцева. - Иркутск : Иркут. унт, 2007. - 255 с.
3. Вольцингер Н. Е. Теория мелкой воды / Н. Е. Вольцингер, Р. В. Пясовский. - Л. : Гидрометеоиздат, 1977. - 207 с.
4. Марчук Г. И. Методы вычислительной математики / Г. И. Марчук. - М. : Наука, 1989. - 608 с.
5. Самарский А. А. Методы решения сеточных уравнений / А. А. Самарский, Е. С. Николаев. - М. : Наука, 1978. - 592 с.
6. Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики / Н. Н. Яненко. - Новосибирск : Наука, 1967. - 195 с.
Arguchintsev V. K., Arguchintseva A. V., Ubonova L. V.
Numerical modelling of streams and pollutions transport in water-currents of the Baikal region
Annotation. The mathematical model for calculation of a velocity field of water-currents and distributions of polluting components after their dump in the certain site of the river is offered. Calculations for concrete objects are resulted. Key words: model, a water-current, a velocity field, pollution.
Аргучинцев Валерий Куприянович д-р техн. наук
Иркутский государственный университет
664003, Иркутск, ул. К. Маркса, 1 зав. каф. метеорологии и охраны атмосферы тел.: (395-2) 42-67-20
Убонова Людмила Владимировна Иркутский государственный университет,
664003, Иркутск, ул. К. Маркса, 1 ведущий инженер каф. метеорологии и охраны атмосферы тел.: (395-2) 42-67-20
Аргучинцева Алла Вячеславовна д-р техн. наук
Иркутский государственный университет
664003, Иркутск, ул. К. Маркса, 1 декан географического факультета тел.: (395-2) 42-56-84
