Маргвелашвили. Численное моделирование влияния паводков на перенос чернобыльских радионуклидов в растворе и на взвесях в Киевском водохранилище Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

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

Том 2, № 2, 1997

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ ПАВОДКОВ НА ПЕРЕНОС ЧЕРНОБЫЛЬСКИХ РАДИОНУКЛИДОВ В РАСТВОРЕ И НА ВЗВЕСЯХ В КИЕВСКОМ ВОДОХРАНИЛИЩЕ^

М. И. ЖЕЛЕЗНЯК, Н. Ю. МАРГВЕЛАШВИЛИ Институт проблем математических машин и систем

НАН Украины, Киев

Построена трехмерная численная модель поля течений и распределения радиоактивной примеси в Киевском водохранилище. Проведено сравнение результатов численных расчетов с данными натурных измерений. Исследуется влияние эффектов гидравлической крупности частиц на процесс формирования общей картины загрязнения вод водоемa. Анализируется влияние паводковых явлений на процессы распределения примеси.

1. Введение

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

Задачи миграции радионуклидов в водных средах в последние десятилетия получили активное развитие и представлены рядом моделей различного уровня сложности [1, 2]. Миграция радионуклидов по Днепровскому каскаду водохранилищ после Чернобыльской аварии рассчитывалась ранее в рамках камерных, одномерных русловых и двумерных плановых моделей [3-5]. Для более детализированного анализа процессов в области больших градиентов концентраций радионуклидов у дна при разных гидродинамических условиях необходимо построение трехмерной модели.

*© М.И.Железняк, Н. Ю. Маргвелашвили, 1997.

^Работа выполнена при поддержке Гранта Международного Научного Фонда №К3У100 и исследовательского контракта МАГАТЭ №7334/Ш/КБ.

Чернобыльская АЭС расположена в 32 км вверх по течению реки Припять от места ее впадения в Киевское водохранилище. Водохранилище имеет длину 80 км и средний обьем 3,7 км3, средняя глубина водоема 4 м, максимальная глубина 15 м. После Чернобыльской катастрофы в придонных отложениях водоема накопилось около 2500 Ки цезия-137. Поскольку все годы после аварии в водохранилище проходили только паводки малой и средней водности, сохраняется потенциальная угроза вымывания цезия из донных отложений в паводки большой водности.

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

2. Описание модели

Расчеты проводились на основе трехмерной численной модели циркуляции вод и транспорта примеси в водоеме. Гидродинамический блок модели с учетом предположения о гидростатическом распределении давления по глубине имеет вид:

ди + оУ + т _ 0 (1)

дх ду дх '

ди + иди + уШ + ^сV _ дП + 0. (иди) + АДи (2)

дЬ дх ду дх д дх дх и дх '

— + и— + V— + Ш— + ¡и _ — + — (и—) + АДУ. (3)

дЬ дх ду дх д ду дх У дх '

д2 д2

Здесь Д _ —— + ——, и, V, Ш — соответственно х, у, х компоненты скорости жидко-

дх2 ду2

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

Параметризация эффектов турбулентного обмена по вертикали проводится на основе гипотезы пути смешения Прандтля с учетом функциональной зависимости масштаба турбулентности от глубины V _ /2|д[//дх|. При наличии одного придонного пограничного слоя масштаб турбулентности определяется по формуле Монтгомери и имеет вид [6]: I _ ( к — х + х^), где к, х^ — глубина водоема и уровень шероховатости донной поверхности соответственно. Согласно литературным источникам [6, 9], представленные соотношения позволяют достаточно удовлетворительно рассчитывать характеристики турбулентного обмена водных систем.

Граничные условия задачи на свободной поверхности жидкости х _ п(х,у,Ь) записываются в виде

д^ ди тх дУ Тх

т _ ^ _ _ 70-

Здесь тх, ту — компоненты касательного напряжения, создаваемого на поверхности жидкости под влиянием ветрового воздействия, р0 — плотность жидкости. На донной поверхности ставятся условия прилипания:

и _0' V _0' Ш _0.

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

Уравнение транспорта взвешенных наносов имеет вид:

дБ ггдБ т ЯБ /ттг „,чдБ д , дБ, , А „

Ж + + Ч +(и' - ^Ж = &) + 4ДБ (4)

где Б — концентрация взвеси, Жо — заданная скорость осаждения частиц примеси в поле силы тяжести. Условие отсутствия потока примеси через свободную поверхность водоема записывается в виде

дБ

(Ж - Ж0)Б = V— г = п.

дг

Вертикальный поток взвешенной примеси на донной поверхности водоема определяется процессами взмучивания и осаждения частиц примеси и описывается следующей зависимостью:

дБ

V— + Ж05 = - qb, г = к - а, дг

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

^ = Г Жо(Б* - 5о), Бо >Б,, ь = Г 0, Бо > Б,,

д 1 0, Бо < Б,, д \ ЕгЖо(Б, - Бо), Бо < Б,,

где Ег — коэффициент эрозии, характеризующий защищенность дна от размывов и изменяющийся в пределах 0.0—1.0, Бо — концентрация примеси на уровне г = к - а, Б, — придонная равновесная концентрация примеси, соответствующая несущей способности однородного стационарного потока при тех же локальных значениях параметров. Расчет значений Б, проводится по формуле Бийкера [7].

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

дС Г7дС ТДдС тт,дС д. дС. ^ _ ,г.

Ж + и Ж + % + иЖ = а!('Ж) + 4ДС - АС - Б(К-С - С')' (5)

Здесь К — коэффициент распределения при условии гидродинамического и химического равновесия системы частицы примеси—жидкость, К = Св/С при Т ^ то, а^ — характеристика скорости обмена в системе жидкость—взвешенные частицы примеси, А — постоянная затухания данного радионуклида.

Граничные условия на свободной поверхности записываются в виде

дС

V—- = Ж С, г = п-дг

Диффузионный поток на донной поверхности описывается соотношением дС

и— _ р3(1 — e)Z*alíз(KdC — Сь)' х _ к — а'

где е — пористость дна, Z* — толщина слоя донных отложений, а^з — характеристика скорости обмена в системе дно—жидкость, р3 — плотность частиц прмеси.

Транспорт радионуклидов взвешенными частицами примеси описывается уравнением

дБС3 1ТдБС3 удБС^ Ш ) дБС3 _ _д_( дБС3)

дЬ дх ду 0 дх дх дх

+АДБС3 — ХБС3 — а^Б (КаС — С3). (6)

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

дБС3

(Ш — Шо)БС3 — _ 0' х _ п

дх

ШоБС3 + _ С3д3 — Сьд°' х _ к — а.

дх

Изменение толщины слоя донных отложений описывается уравнением деформации донной поверхности

Р3(1 — е) _ д3 — дь. (7)

Динамика эволюции концентрации радионуклидов в слое донных отложений описывается уравнением

М _ — СЬ) — р^(С3д — СУ). (8)

3. Численная реализация модели

Разностная аппроксимация задачи строится интегро-интерполяционным методом на разнесенной сетке [10].

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

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

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

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

Тестирование модели проводится путем сравнения результатов численных расчетов с решением одномерной линеаризованной модели и модели Экмана [8].

4. Еезультаты верификации модели и численных расчетов

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

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

Получено удовлетворительное качественное и количественное согласование распределения полей стоковых течений в водоеме. В отсутствии ветра, при уровне входного потока в 2000 м3/с максимальные значения осредненных по глубине скоростей составляют 10 см/с в верхней мелководной и 3-5 см/с в нижней глубоководной части бассейна, что соответствует обобщенным данным измерений поля течений в Киевском водохранилище.

На рис. 1 представлены картины распределения горизонтальных компонент скорости течений в плоскости приповерхностного (а) и придонного (б) сечений водоема при условии ветрового воздействия на поверхности. Направление ветра противоположно направлению стокового течения, значение модуля скорости ветра составляет 7 м/с, расход жидкости 2000 м3/с.

Расчеты транспорта взвеси и распределения радионуклидов, сорбированных частицами взвеси, проводились для фракций частиц размерами в 2 и 20 мкм. Пересчет измеренных значений концентрации радионуклидов, сорбированных всей взвесью для фракций частиц указанных размеров, проводился в рамках предположения о равномерном распределении радионуклидов по частицам разного диаметра. При численных расчетах использовались следующие значения коэффициентов обмена и распределения: а1>2 = 1/(8, 64-104) с= (сутки)—1, а1>3 = 1/(3.15 • 107) с=(год)-1, Кл = 5000 л/кг, Клъ = 3000 л/кг.

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

Рис. 1. Поле горизонтальных компонент скорости в а) приповерхностной (|^max | =8 см/с) и б) придонной (|^max | = 2.3 см/с) областях водоема.

Рис. 2. Поле концентрации радионуклидов на частицах взвеси диаметром 2 мкм (ПКи/л) при расходе жидкости 2000 м3/с.

Рис. 3. Поле концентрации радионуклидов на частицах взвеси диаметром 20 мкм (ПКи/л) при расходе жидкости 2000 м3/с.

О 5 10 15 20 25 30

Рис. 4. Поле концентрации радионуклидов на частицах взвеси диаметром 2 мкм (ПКи/л) при расходе жидкости 16000 м3/с.

0 5 10 15 20 25 30

Рис. 5. Поле концентрации радионуклидов на частицах взвеси диаметром 20 мкм (ПКи/л) при расходе жидкости 16000 м3/с.

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

Сравнение расчетных значений концентрации растворимой формы радионуклидов с результатами наблюдений на выходе из водоема за период апрель—май 1988 г. показывает расхождение расчетных и измеренных значений концентрации до 10 % в конце первого месяца и до 30 % в конце второго месяца.

На рис. 4, 5 показаны данные расчетов распределения поля концентрации радионуклидов, сорбированных частицами взвеси при экстремально больших для условий Киевского водохранилища значениях входного потока жидкости — 16000 м3/с. Максимальные значения суммарной концентрации радионуклидов на обоих фракциях частиц достигают

16 ПКи/л, что более чем на порядок превышает аналогичные значения, полученные при расчетах с входным потоком в 2000 м3/с. Наблюдается незначительное падение с 20 до

17 ПКи/л максимальных значений концентрации радионуклидов в растворе. Получено двукратное увеличение средних значений суммарной коцентрации цезия-137 в растворе и на частицах взвеси указанного выше размера.

5. Заключение

Построена трехмерная численная модель распределения поля течений и транспорта примеси в Киевском водохранилище.

Получено удовлетворительное совпадение результатов расчетов гидродинамического блока модели с данными натурных измерений.

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

Проведены расчеты распределения полей концентрации взвеси и радионуклидов при экстремально больших для условий Киевского водохранилища значениях расхода воды (16000 м3/с). Получено двукратное увеличение суммарной концентрации радионуклидов в растворе и на частицах взвеси диаметром 2 и 20 мкм по сравнению с аналогичными данными, рассчитанными при значениях расхода в 2000 м3/с.

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

[1] CODELL R. B., KEY K.T., WHELAN G. Collection of mathematical models for radionuclide dispersion in surface water. NUREG-08668, Washington, DC, 1982.

[2] Onishi Y., Serne J., Arnold E. ET AL. Critical review: radionuclide transport, water quality, mathematical modelling and radionuclide adsorption/desorption mechanism. NUREG/GR-1322, Pacific Northwest Laboratory, Richland, 1981.

[3] Cherneva J.T., Hursin S.L., Tkalich P.V., Zheleznyak M.J. A model of the hydrodynamics processes and sediment transport over complicated bottom. In "Proc. 6th National Congress of Theoretical and Applied Mechanics", Varna, 1989, 65-1-65-6.

[4] Zheleznyak M., Demchenko R., Khursin S. et al. Mathematical modeling of radionuclide dispersion in the Pripyat-Dnieper aquatic system after the Chernobyl accident. The Science of the Total Environment, 112, 1992, 89-114.

[5] Zheleznyak M., Kuzmenko Yu., Tkalich P. et al. Modelling of radionuclides transport in the set of river reservoirs (Eds. A. Peters et al.). Computational Methods in Water Resources X, Kluwer Academic Publishers, Dordrecht, The Netherlands, vol. 2, 1994, 1189-1196.

[6] MONTGOMERY R. Generalisation for cylinders of Prandtl's linear assumption for mixing length. Ann. N.Y. Acad. Sci, 44, 1981, 89-103.

[7] Rijn L. van. Sediment transport, Part II: Suspended load Transport. J. Hydraulic Engineering, 110, 1979, 1613-1641.

[8] ВольцинГЕР Н., Клеванный К., ПЕЛИНОВСКИй Е. Длинные волны на мелкой воде. Гидрометеоиздат, Л., 1985.

[9] Demchenco R. I., Zheleznyak M. J., Koziy L.A. Modelling of sedimentation and radionuclides deposition in a bottom trap. In "Proc. Tenth International Conference on Computational Methods in Water Resourses", Germany, Jul, vol. 2, 1994, 1341-1348.

[10] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Мир, M., т. 2, 1990.

Поступила в редакцию 15 сентября 1995 г.