Исследование процессов переноса, диффузии и трансформации радиоактивных примесей, поступающих в атмосферу при авариях на объектах энергетики

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

Том 14, № 2, 2009

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

атмосферу при авариях на объектах энергетики

Г. В. Муратова, М.В. Глушлнин Южно-Российский региональный центр информатизации ЮФУ,

Ростов-на-Дону, Россия e-mail: muratova@rsu.ru, mig29m@rambler.ru

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

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

Введение

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

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

© ивт со ран, 2009.

1. Постановка задачи и численная реализация

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

1) на высокодисперсные (размер частиц менее 0.1 мкм);

2) на среднедисперсные (размеры частиц от 0.1 до 1 мкм);

3) грубодисперсные (размер частиц более 1 мкм).

Частицы аэрозолей находятся в постоянном броуновском движении, интенсивность которого увеличивается с уменьшением размеров частиц. Загрязняющие вещества в атмосфере могут изменять свое состояние. Конденсация газов может приводить к образованию аэрозолей и, наоборот, газы могут образовываться в результате испарения частиц аэрозолей. В процессе броуновского движения частицы аэрозолей могут сталкиваться друг с другом и слипаться. В результате образуются новые, более крупные, частицы. Этот процесс называется коагуляцией. Радиоактивные вещества могут попадать в атмосферу как в форме газов, так и в виде аэрозолей. К газообразным примесям относятся радиоактивные инертные газы, тритий, оксиды и сульфиды радиоактивного углерода, пары радиоактивного йода. Другие радиоактивные вещества, как правило, попадают в атмосферу в виде аэрозолей [2].

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

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

Таблица 1. Классификация атмосферных примесей [1]

Тип примеси Характеристика

Зола (грит) Крупные твердые частицы более 76 мкм

Пыль Твердые частицы размером более 1 и менее 76 мкм

Дым Твердые частицы диаметром менее 1 мкм

Дождь Жидкие частицы диаметром более 400 мкм

Морось Жидкие частицы размером более 100 и менее 400 мкм

Фог Жидкие частицы размером более 10 и менее 100 мкм

Туман Жидкие частицы размером менее 10 мкм

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

m

m

u(x,y,z) = uo(x,y)[ h-

v{x,y,z) = V0(x,y)^ h

где u0, v0 — горизонтальные компоненты скорости ветра v на высоте анемометра, м/с; ha — высота анемометра, м; m — показатель степени, зависящий от стратификации, значения его даны в табл. 2. Для экстраполяции остальных физических величин на всю рассматриваемую область используется способ, предложенный в работе [3].

Математическая модель переноса пассивных радиоактивных примесей в атмосфере [4-6] основывается на нестационарном трехмерном уравнении турбулентной диффузии для средних значений объемных активностей Ay в прямоугольных декартовых координатах. Так как число Маха M << 1, то можно полагать, что плотность воздуха постоянна (р = const) и среда несжимаема (div v = 0). Таким образом, в области

Qt = П[0,Т], П= {x е [0,Lx], y е [0,Ly], z е [0,Lz]},

где ось x направлена на юг, ось y — на восток, а ось z — вертикально вверх, уравнение, описывающее конвективный перенос, турбулентную диффузию, эмиссию, гравитационное оседание и трансформацию (радиоактивный распад) i-й радиоактивной примеси из N (где N — количество рассмотренных примесей) в атмосфере в системе декартовых координат можно представить в следующем виде:

dAy dAy dAy dAy

+ u—--+ v—---+

dt dx dy g dz

- ¿

где Лу — объемная активность г-й компоненты радиоактивной примеси, Бк/м3; и, V, — компоненты вектора скорости ветра V = ^(и^^), м/с; /ша — скорость гравитационного оседания радиоактивной примеси, м/с; К — коэффициент горизонтальной турбулентной диффузии, м2/с; Кг — коэффициент вертикальной турбулентной диффузии, м2/с; \г — постоянная радиоактивного распада г-й компоненты примеси, с-1; Яир — функция эмиссии приподнятых источников, Бк/(м3 • с).

Таблица 2. Зависимость эмпирического показателя степени т от стратификации в формуле для вычисления вертикального профиля ветра

Класс атмосферной стабильности A B C D E F G

m 0.08 0.165 0.215 0.31 0.405 0.43 0.44

Начальные условия имеют вид

(ж, у,*, 0) =

Краевые условия на границе области дП ставятся следующим образом. На свободной границе они имеют вид

дАгу

— 0 при ж = 0, ж = Ь

дж дА*

—^ = 0 при у = 0, у = Ь дУ

1 п т

— 0 при * = Ьх.

X )

У'

ному турбулентному потоку примеси — К у потока — и>д . Случай 7 = 0 соответ-

д*

Их физический смысл заключен в беспрепятственном проникновении загрязнений сквозь границу [3].

Краевые условия на подстилающей поверхности следует подбирать исходя из анализа физических процессов [7], происходящих на этой поверхности. Достаточно общим краевым условием на подстилающей поверхности при * = 0 является условие

дА*

К*~&Г + щ= ТА, + QCc0wn при * = 0,

где 7 — коэффициент поглощения примеси поверхностью, характеризующий взаимодействие радиоактивной примеси с подстилающей поверхностью, м/с; QCcown — функция эмиссии наземных источников, Бк/(м2 • с). Для общности мы допускаем здесь наличие гравитационного оседания со скоростью шд, приводящего к добавлению к вертикаль-

дА1

д*

ствует "отражению" примеси от стенки, случай 7 = то — "поглощению" примеси, а случай 0 < 7 < то — промежуточной ситуации частичного отражения и частичного поглощения.

Для вычисления коэффициента вертикальной турбулентной диффузии К используется эмпирическая формула [8]:

г |к* (Ь), *<ь,

К(г) = ] К(Ь), Ь < * < я,

I 1 + (к(Ь) — 1) ехр (я — *), * > я,

где Ь — высота приземного слоя атмосферы, м; Я — высота верхней границы пограничного слоя атмосферы, м; К (Ь) — значение коэффициента вертикальной турбулентной диффузии на высоте Ь, м2/с.

Коэффициент горизонтальной турбулентной диффузии К вычисляется по эмпирической формуле [3]:

К = а© тах(0.5, |у|)Я,

где а© — угол горизонтальной флуктуации направления ветра, рад; V — скорость ветра, м/с.

Значения Ь, Я, К(Ь) и а© зависят от стратификации атмосферы и определяются по табл. 3 в точках расположения метеостанций, а затем интерполируются на всю расчетную область.

Таблица 3. Зависимость физических величин и эмпирических параметров от стратификации атмосферы [9]

Класс атмосферной стабильности Л В С Б Е Р С

° 25 20 15 10 5 2.5 1.7

Н, м 2000 1500 1000 750 300 250 250

Л, м 250 250 150 150 150 100 100

К(Л), м2/с 160 100 70 15 5 1.5 0.13

Важным фактором, влияющим на поведение аэрозолей в атмосфере, является гравитационное оседание. Оно приводит к очищению атмосферы и загрязнению подстилающей поверхности. Скорость оседания частиц примеси под действием силы тяжести зависит от их размеров и плотности, а также от вязкости и плотности воздуха. Характер движения частиц определяется числом Рейнольдса [10]:

Ке = РТрМ ,

ц

где р — плотность воздуха, кг/м3; гр — радиус частицы, м; V — скорость ветра, м/с; ц — динамическая вязкость воздуха, кг/(м • с).

Для мелких частиц, когда Ке < 1, скорость оседания шд определяется по формуле Стокса:

2(рР - р)# 2

"д _ '

где рр — плотность частиц, кг/м3; д — ускорение свободного падения, м/с2.

Для частиц, имеющих большие размеры (Ке > 1), скорость падения может быть вычислена с помощью следующего выражения [11]:

^ _ 2(рр - р)д / ^ ке

_ 9 ц Н 16

Изменение поверхностной активности Лг8 ¿-й радиоактивной примеси на подстилающей поверхности вследствие гравитационного оседания и радиоактивного распада определяется с помощью уравнения

дА

ог + А* Агг _ 7Ау0 ,

где Аг — поверхностная активность ¿-й радиоактивной примеси, Бк/м2; 7 — коэффициент поглощения примеси поверхностью, характеризующий взаимодействие радиоактивной примеси с подстилающей поверхностью, м/с; А^0 — объемная активность ¿-го радионуклида на высоте г _ 0, Бк/м3. Коэффициент поглощения примеси поверхностью 7 зависит от свойств подстилающей поверхности.

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

tth = \ (xi, yj, zk), xi = ihx, yj = jhy, zk = khz

Lx

L,.

L2

i = 0, Nx, j = 0, Ny, k = 0, Nz, Nx = , Ny = , Nz =

hx

hy hz

Пусть QT = {tn = птп, n = 0, Nt, to = 0, tNt = T} — произвольная сетка на отрезке 0 < t < T с шагами Tn = tn+1 — tn. В разностной схеме отнесем переменные к вершинам ячеек сетки. Приближенное решение ищется в виде сеточной функции дискретных аргументов ^"jk = ^(xi,yj, Zk,tn), которая считается приближенным значением проекции искомой функции Ay (xi,yj , Zk,tn).

При численном моделировании распространения радиоактивных примесей предъявляются дополнительные требования к конечно-разностным аппроксимациям уравнения (1) и методам их решения [12]. Так как объемная активность по физическому смыслу — неотрицательная величина, целесообразно использовать так называемые монотонные схемы, позволяющие получать неотрицательные решения. При построении вычислительного алгоритма для уравнения, которое получается после перехода к безразмерным переменным в уравнении турбулентной диффузии (1), воспользуемся методом расщепления по физическим процессам и на каждом малом интервале времени [tn; tn+1] длиной тп используем схему, состоящую из трех этапов [12].

1. Перенос примеси по траекториям:

Ж

dx

?— дУ

— + и-—+ V— + (w — w„ ) —

cu a,. g dz

0.

2. Турбулентная диффузия примесей:

Ж

d

dx

= Qup + Q^i^ ) + - ( ^ ) + - ( Kzlt

d_ дУ

дУ

A

dz

z dz

3. Радиоактивный распад:

f + » = 0.

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

— + U— = 0 St dx

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

В качестве примера нелинейных разностных схем решения уравнения переноса (5) рассмотрим монотонную схему, разработанную Ван Лиром и основанную на нелинейной схеме Фромма, которая аппроксимирует исходное дифференциальное уравнение со вторым порядком точности по пространственным переменным и по времени. Идея схемы

Фромма заключается в использовании того известного факта, что решение уравнения переноса по схеме Лакса—Вендроффа, имеющего вид

'3+1 = Ь

"отстает" по времени от решения дифференциальной задачи. Заметим, что при г = 1 схема абсолютно точна. Переход с ]-го на (] + 1)-й уровень можно осуществить, предварительно пройдя по характеристике с г = 1, а затем двигаясь по времени назад, т. е. осуществив следующее преобразование:

ф3+1 = Ь

Ясно, что в этом случае фз+1 будет опережать на шаг по времени решение дифференциальной задачи. Если в качестве решения взять полусумму решений (6) и (7), то фазовая ошибка существенно будет уменьшена. Таким образом, схема, предложенная Фроммом, имеет вид

'"+1 = \

Прежде чем построить монотонную схему на основе (8), представим в развернутом виде конечно-разностную аппроксимацию (6) и (7):

^Г1 = ^ - г&г-1/2'А - г

= ^21 - г\-1/2'А - г

г 2

где '1, '2 — решения задач (9) и (10) соответственно, и

Аг-

-1/2' = 'г - Аг-1

Рассмотрим класс монотонных схем, решения которых удовлетворяют следующему условию:

3 + 1 - 3

0 < '-' < 1.

'г - Аг-1

Если специальным образом выбирать соответствующие управляющие функции, то удается построить монотонную версию схемы для (9) и (10).

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

'+ = А - а1гАг-1/2'3 - «2т(1 - г) (Аг+1/2'3 - Аг-3/2'3) +

'4

(1 - г) ГБ(Л) (А, а при иг < 0

г

+ аз -

'+1 = ' - а1гАг+1/2''3 - а24(1 + г) (Аг+3/2'3 - Аг-1/2'3) + г

+ а3г

где r = Дг+1/2 -

-— max |ui| — число Куранта, Ai_3/2—

bx -i+1

-г

s

Aj+3/2- = -г+2 - -г+1,

S

-i_1 - -i_2, Ai-1/2^

-г

-г_1,

S(Z ) = |Ai+3/2 —| - |Ai+1/2-

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

т

Схема (11), (12) устойчива, если r = max |uj| < 1. При а1 = 1, а2 = а3 = 0 получаем схему Годунова. При а1 = а2 = 1, а3 = 0 получаем немонотонную консервативную схему Фромма (второго порядка). При а1 = а2 = а3 = 1 получаем монотонную схему Ван Лира (второго порядка).

На втором этапе полученное в ходе расщепления уравнение (3) решается неявной схемой двуциклического покомпонентного расщепления [12]:

—j

j+a/6 __ ^j+(«_1)/6 1

T„ t/2

-j+4/6 - -j+2/6 !

+ 2 (^V+a/6 + Ла^+(а_1)/6) =0, a = 1, 2;

Tnt

+ 2 (Л^+4/6 + Л V+2/6) = 0;

. j+(7_a)/6 - wj+(6_a)/6 1

где

Trai/2

д

Л1 - = - ¿1*1

Л = Л^д-2 дУ V дУ

0, a = 2,1,

д

Л3- = Kz-jf-

dz

д-

dz

Численная схема (13) имеет второй порядок точности по пространственным переменным и времени.

На третьем этапе полученное в ходе расщепления уравнение (4) представляет собой однородное линейное обыкновенное дифференциальное уравнение, решение которого имеет вид

-

га+1

exp (-Arn).

2. Вычислительные эксперименты 2.1. Программный комплекс RAD

На основе построенных моделей авторами статьи создан программный комплекс RAD (RADiation), позволяющий проводить оперативные расчеты распространения загрязнения в воздушной среде. Разработанный программный комплекс включает следующие модули: модуль ввода начальных данных, расчетный модуль, база данных с системой управления, модуль визуализации результатов расчетов. При помощи расчетного модуля вычисляется поле объемной активности радиоактивной примеси в атмосфере для

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

Программа RAD создана с помощью среды разработки Eclipse IDE (версия 3.2), объектно-ориентированного языка программирования Java и библиотек классов платформы Java Platform, Standard Edition. С помощью программного комплекса RAD выполнен ряд вычислительных экспериментов, моделирующих различные штатные и нештатные ситуации работы АЭС.

2.2. Результаты расчетов

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

Вычислительный эксперимент 1. В качестве области моделирования был выбран прямоугольный параллелепипед размером 60 х 60 х 2 км. Шаг по времени вычислялся автоматически из условия устойчивости явной схемы Ван Лира. Горизонтальные шаги регулярной сетки hx = hy = 600 м, по вертикали hz = 50 м. Размер сетки 101х101х41. Рассматривался мгновенный точечный источник загрязнения, находящийся в центре области моделирования, высота которого равна 100 м. Класс атмосферной стабильности — D (нейтральная). Направление ветра — юго-восточное, скорость ветра — 7 м/с. В качестве загрязняющих веществ были взяты благородные радиоактивные газы 85Kr, 87Kr, 41Ar, 133Xe, 135Xe, для которых проведены вычислительные эксперименты. Для газообразных примесей скорость гравитационного оседания wä и вертикальная компонента скорости ветра равны 0. Для ядерных реакторов типа ВВЭР в газообразных выбросах доля благородных радиоактивных газов 133Xe (период полураспада 5.24 сут.) и 135Xe (период полураспада 9.08 ч) составляет 72.0 и 13.2% соответственно [13]. Результаты вычислительных экспериментов, моделирующих ситуацию залпового выброса 133Xe из точечного источника общей активностью 3.7 • 1013 Бк, представлены на рис. 1.

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

Вычислительный эксперимент 2. Моделируется аварийный выброс, случившийся 6 апреля 1993 года на радиохимическом заводе (РХЗ) Сибирского химического комбината (СХК) [14, 15]. В результате повреждения технологического оборудования и взрыва газов были разрушены конструкции здания и произошел выброс загрязняющих веществ через вытяжную трубу с эффективной высотой выброса 150 м и через развал стены высотой 15 м. В разных источниках [14-16] приводятся различные значения суммарной активности выброса (от 400 до 680 Ки). Анализ проб почвы [14, 15] позволил установить следующий радионуклидный состав выброса, %: 95Zr — 20, 95Nb — 42, 103Ru — 2, 106Ru — 35, 1 % составили 94Nb, 137Cs, 51Cr, 125Sb. Выброс и последующее распространение радиоактивных веществ в атмосфере происходили при юго-западном ветре с изменениями направления от 190 до 220° и скорости 8-13 м/с, осадков в момент

Рис. 1. Графики изолиний объемной активности Лу (106Бк/м3) 133Хе на высоте Н = 100 м в моменты времени: а — Ь = 20 мин; б — Ь = 40 мин; в — Ь = 60 мин

в

выброса не наблюдалось, инсоляция и облачность были умеренные. След пересек автомобильную дорогу Томск—Самусь, а также повлиял на радиационную обстановку в деревне Георгиевка. Центральная часть следа на трассе ограничена изолинией мощности экспозиционной дозы 100 мкР/ч [14]. Уровень мощности дозы в Георгиевке составил 20-40 мкР/ч [14, 17].

Для моделирования распространения облака радионуклидов и расчета мощностей эквивалентных доз были заданы следующие условия: тип выброса мгновенный; высота выброса — 150 м; суммарная активность выброса — 600 Ки; радионуклидный состав выброса, %: — 20, 95№ — 43, 103Ии — 2, 106Ил — 35; скорость ветра — 10 м/с; направление скорости ветра в процессе расчета изменялось от 190 до 220° (табл. 4). Источник загрязнения находится в центре области моделирования. В качестве области моделирования был выбран прямоугольный параллелепипед размером 60 х 60 х 1.5 км. Шаг по времени вычисляется автоматически из условия устойчивости явной схемы Ван Лира. Горизонтальные шаги регулярной сетки кх = ку = 1000 м, по вертикали кг = 50 м.

Рис. 2. Графики изолиний суммарной поверхностной активности As (Ки/км2) в моменты времени: а — £ =10 мин; б — £ = 20 мин; в — £ = 30 мин; г — £ = 40 мин; д — £ = 50 мин

Таблица 4. Изменение скорости и направления ветра

№ Скорость ветра, м/с Направление, ° Время начала, ч:мин.

1 10 220 12:58

2 10 200 13:07

3 10 180 13:09

4 10 220 13:13

Размер сетки 61 х 61 х 31. Класс атмосферной стабильности — D (нейтральная). Результаты моделирования выпадений радионуклидов на поверхность почвы представлены на рис. 2.

Результаты проведенных вычислительных экспериментов удовлетворительно согласуются с результатами, полученными ранее другими авторами [11, 17-19] для районов, схожих с районом Волгодонской АЭС.

Заключение

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

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

[1] Страус В. Промышленная очистка газов: Пер. с англ. М.: Химия, 1981. 616 с.

[2] Нормы радиационной безопасности (НРБ-99): Гигиенические нормативы. М.: Центр сан.-эпидем. нормирования, гигиен. сертификации и экспертизы Минздрава РФ, 1999. 116 с.

[3] Самарская Е.А., Сузан Д.В., Тишкин В.Ф. Построение математической модели распространения загрязнений в атмосфере // Матем. моделирование. 1997. Т. 9, № 11. С. 59-71.

[4] Берлянд М.Б. Прогноз и регулирование загрязнения атмосферы. Л.: Гидрометеоиздат, 1985. 270 с.

[5] Бруяцкий Е.В. Теория атмосферной диффузии радиоактивных выбросов. Киев: Ин-т гидромеханики НАН Украины, 2000. 443 с.

[6] Монин А.С. Атмосферная диффузия и загрязнение воздуха. М.: ИЛ, 1962. 512 с.

[7] Монин А.С., Яглом А.М. Статистическая гидромеханика. Теория турбулентности. СПб.: Гидрометеоиздат, 1992. Т. 1. 695 с.

[8] Draxler R.R. Modeling the results of two recent mesoscale dispersion experiments // Atmospheric Environment. 1979. N 13. P. 1523-1533.

[9] КАлиткин Н.Н., Карпенко Н.В., Михайлов А.П. и др. Математические модели природы и общества. М.: Физматлит, 2005. 360 с.

[10] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1970. 904 с.

[11] Берлянд М.Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы. Л.: Гидрометеоиздат, 1975. 448 с.

[12] Алоян А.Е. Динамика и кинетика газовых примесей и аэрозолей в атмосфере: Курс лекций. М.: ИВМ РАН, 2002. 201 с.

[13] Козлов В.Ф. Справочник по радиационной безопасности. М.: Энергоатомиздат, 1991. 352 с.

[14] Лысцов В.Н., Иванов А.Б., Колышкин А.Е. Радиоэкологические аспекты аварии в Томске // Атомная энергия. 1993. Т. 74, № 4. С. 364-367.

[15] САвкин М.Н., Титов А.В. Анализ радиационной обстановки на следе аварийного выброса Сибирского химического комбината // Медицина катастроф. М.: Ин-т биофизики, 1995. С. 76-84.

[16] Булатов В.И., Чирков В.А. Томская авария: мог ли быть сибирский Чернобыль? Новосибирск: ЦЭРИС, 1994. 32 с.

[17] Белов И.В., Беспалов М.С., КлочковА Л.В. и др. Сравнение моделей распространения загрязнений в атмосфере // Матем. моделирование. 1999. Т. 11, № 8. С. 52-64.

[18] НьидстАдтА Ф.Т.М., Ван ДопА Х. Атмосферная турбулентность и моделирование распространения примесей. Л.: Гидрометеоиздат, 1985. 351 с.

[19] Белов И.В., Беспалов М.С., КлочковА Л.В. и др. Транспортная модель распространения газообразных примесей в атмосфере города // Матем. моделирование. 2000. Т. 12, № 11. С. 38-46.

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