Моделирование совместного загрязнения грунта и атмосферы Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

УДК 504.064.2.001.18

МОДЕЛИРОВАНИЕ СОВМЕСТНОГО ЗАГРЯЗНЕНИЯ ГРУНТА И АТМОСФЕРЫ

© 2005 г. И.Э. Наац, В.С. Кириллов, П.В. Корчагин

The mathematical modeling of the unsteady pollution substance transportation in the atmosphere, the substance sedimentation to the pond surface and the permeation to the ground is provided in the article.

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

Структуру подобной модели (решающего алгоритма) можно представить в виде следующей диаграммы (рис. 1).

Вычисление

поля

концентрации

ЗВ в воздухе

(1)

Вычисление концентрации вещества, "оседающего на поверхность водоема (2)

Вычисление концентрации

ЗВ,

проникающего в грунт (3)

Вычисление поля концентрации ЗВ в грунте под водоемом (5) Вычисление концентрации ЗВ в водоеме (4)

Рис. 1. Структура решающего алогритма в задаче совместного загрязнения атмосферы, водоема и грунта

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

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

1. Вычислительная схема для задачи распространения загрязнения в атмосфере

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

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

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

X 1

соответственно. Их размерность в последующих вы-

2

числениях принимается [м /с]. Считаем задачу двумерной по пространственным координатам х и г. Уравнение распространения примеси

д+у£±-д-дКх-дКг= Б, (1) дt дх дг дх дх д дг

где q(х, г, /) - поле концентрации аэрозоля; V - скорость ветра; w - скорость оседания аэрозоля; Б( 0 = Q (0 8 (х-х 0) 8 (г-г 0) - функция точечного источника мощности Q, расположенного в точке (х0, г0). Расчетная область принимается в форме прямоугольника D={(х,z)е[0,L]х[0,H]}. В момент времени ¡=0 концентрация считается известной: q\ =q (х,г).

Граничные условия задаем следующим образом: концентрация q равна нулю на верхней границе (г=Н) и там, где воздушные потоки втекают в область, и на поверхности водоема (г=0, хе [х1,х2]). На участках границы, через которые воздушные потоки покидают область, производная q по нормали к границе равна нулю. На границе с подстилающей поверхностью задаем соотношение: wq+v■дq/дz=аlq, (г=0, х£[ х1,х2]); а - константа, характеризующая степень поглощения

вещества поверхностью земли; V - коэффициент молекулярной диффузии воздуха. В расчетах размерность указанных величин принимается следующей:

2

М=[м /с], [а ]=[м/с] (рис. 2).

Рис. 2. Расчетная область в задаче загрязнения атмосферы

Масса аэрозоля, поглощаемая водоемом в единицу времени, вычисляется по формуле:

M^) = Ц2 V ■ (х,0, {)йх, где хь х2— границы водоема.

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

Допустим, что выброс примеси происходит в светлое время суток, когда в атмосфере наблюдается неустойчивая или безразличная температурная стратификация, определяющая высотные профили скорости ветра У(2) и коэффициентов турбулентности Кх(2), К2(2). В соответствии с принятыми моделями пограничного слоя рассмотрим три случая состояния турбулентной атмосферы до уровня 2=Н=250 м.

Во всех трех случаях в слое 2е [0,ф0] (20=0,05 м -

-5 2

уровень шероховатости) У=0, Кх=Кг^=1,3-10 м /с (рис. 3).

10

[М/с]

ПМ

Н г

Рис. 3. Профили скорости ветра (а) и коэффициентов диффузии (б,в)

Пусть I - атмосфера при неустойчивой стратификации. В слое 2е[20,Ь], где ¿=20 м, скорость ветра из-

дУ и* 1

меняется по закону: -= -

dz

к z(1 + 15z / L)

1/4 '

Здесь u =0,33 м/с - динамическая скорость;

2

к=0,38-постоянная Кармана. K =l -dV/dz, где

z

-1/4

l=Kz(1+14|Ri|) , Ri= -0,3 - число Ричардсона. Выше уровня z=L скорость V постоянна. В слое ze[L,h], где h=150 м - высота приземного слоя, Kz= z43. Выше уровня z=h величина Kz сохраняется постоянной. Коэффициент турбулентной диффузии Kx изменяется следующим образом: в слое ze [z0,h] он пропорционален ln(z/z ), причем Kx(h)=Kz(h); выше уровня z=h он

совпадает с Kz.

Пусть II - атмосфера при нейтральной стратификации. В слое ze [z0,h]: V = — ln—,

к z0

ln( z / z )

Kx = u*Kh-— + v; Kz = и*к(z - z0) + v. Здесь

ln(h / z0)

u*=0,6 м/с. Выше уровня z=h величины V, Kx, Kz остаются постоянными.

Пусть III — атмосфера при нейтральной стратификации и запирающей инверсии на верхнем уровне. До высоты z=h формулы для определения Kx, Kz сов-

а

падают с II. Выше уровня ф=И коэффициенты турбулентной диффузии Кх, К2 линейно убывают с высотой, достигая при ф=Н значений, равных 0,1 ■К (И).

X

Решение уравнения (1) осуществлялось в рамках конечноэлементного метода с использованием покоординатного расщепления и вариационного подхода. В полной мере эта методика описана в [1, 2].

В расчетах, которые ниже иллюстрируются графиками (рис. 4), выбирались следующие исходные данные: высота расчетной области Н=250 м; ширина ¿=5000 м, источник примеси находится в точке с координатами (х0гг0)=(500,50). Водоем расположен на расстоянии 4000 м от левого края области; ширина водоема 2 м. Скорость осаждения частиц примеси ^=2,5 см/с. Функция источника моделируется как источник постоянной мощности 2=360 г/с, начавший действовать в момент времени 1=0.

На рис. 4 приведены установившиеся ^ >2000) поля концентрации примеси для трех случаев I, II, III состояния атмосферы.

Через некоторое время после начала работы источника примесь начинает попадать в водоем. Количество вещества М, поглощаемого водоемом в единицу времени, быстро растет и достигает некоторого установившегося значения (рис. 5).

L х

L х

L х

Рис. 4. Установившееся поле концентрации q для разных состояний атмосферы (I, II, III)

Рис. 5. Поток вещества в водоем после начала работы источника

В случае состояния I пограничного слоя стационарное значение M равно 0,032 г/с; в случае II и III -M=0,009 г/с.

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

2. Вычислительная схема для задачи распространения загрязнения в грунте

Количество вещества, переносимого в грунт за единицу времени, можно рассчитать по формуле:

I1 ={ (- f >,

(2)

где С, - концентрация вещества в водоеме, г/л; у„, -проекция вектора скорости течения жидкости в грунте на нормаль к границе водоема, м/сут; Б - коэффициент гидродинамической дисперсии; I - граница раздела водоем-грунт.

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

С(хy,фу)д- = д\ кх(Xy,дг\ +

д дх V дх I

д

II

ду

ду

+ —I ky(x,у,z,v)—~ I + —I kz(у)

ду

д

dz

д(у + z)

(3)

В этом уравнении неизвестной величиной является \\(х,у,х,/) - давление всасывания жидкости. Опишем кратко основные величины, входящие в (3). Функция С(х,у,?,\\) в левой части представлена в виде:

дОО

С\\) =-, где О - влажность грунта в точке (х,у,2/)

дд\

области интегрирования уравнения (3). В свою очередь функцию у можно связать с давлением жидкости Р формулой

У = Р / РЯ (4)

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

W =

1+ß(_wf

W< 0,

(5)

Qs, цу> 0.

В (5) в, О» п — некоторые коэффициенты, характеризующие физические свойства почвы [4].

Коэффициенты уравнения (3) в правой части кх, ку, к2 - составляющие коэффициента фильтрации вдоль соответствующих координатных осей (их размерность м/сут). При вычислениях, связанных с фильтрацией в грунтах значения этих коэффициентов можно оценить по полуэмпирическим формулам [1] того же вида что и (5), т. е.

к' \< 0,

k (И =

1 + а(_^)Г[

(6)

ks, 0.

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

Уравнение (3) дополняется начальными и граничными условиями \(X,у,2,0)=\ (X,у,2),

на границе раздела

\(X, у, 2, /) = у(х, у, 2, /)

грунт-водоем, и -^г = р(х,у, 2,/) на остальной град п

нице области, в которой решается исходное уравнение (3).

Процесс массопереноса описывается уравнением относительно концентрации растворенного вещества С(х, у, 2, /) в пределах исследуемой области

в

дС д

(

dt

д '

дх

дС

D( х, y, z,V) — дх

Л

д

д y

г дСл

D( х, y, z,V ) — д y

' дС

+—I D(x,y,z, V)— дz l дz

_ V дС+V дС+V дС

' х ^ v ^ z у, дх dy dz

(7)

В этом уравнении О - влажность грунта; Ух= -кх(\\)-д\\/<дх, ¥у= -ку(\\)-д\\/ду, У2= —кг(\)-д\/дг -компоненты скорости фильтрации потока; кх, ку, к2 -коэффициенты фильтрации; \ - давление всасывания, о котором речь шла выше в связи с решением уравнения (3); Б(х,у,1,У) - коэффициент гидродинамической дисперсии, определяемый в модели эмпирической формулой [2] 0(У)=Ом+Ок\У\а, где Бм - коэффициент

молекулярной диффузии; Бк - коэффициент конвек-

О О 0 1 /О

тивной дисперсии; \У\=(Ух +Уу +У2) - модуль скорости фильтрации жидкости в грунте; а - параметр, связанный с физическими свойствами грунта. Начальные и граничные условия уравнения (7) аналогичны тем, которые использовались при решении уравнения (3), а именно, С(х,у,г,0)=С (х,у,г) и

ЗС/дп=Я (ху^,/); С(х,у,2,/)=С (х,у,2,() на границе области решения.

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

Рассмотрим область прямоугольной формы размером Ьх=30 м, Ь2=5 м. Естественное углубление почвы имеет ширину 4 м и глубину Н=2 м, при этом уровень жидкости совпадает с уровнем почвы. Допустим, что основной грунт, в котором происходит фильтрация, имеет следующие характеристики [3]: к* = 0,01; а= 5,25; п = 2; Ох = 0,45; в = 0,745. Будем считать, что на расстоянии Ь=4м от уровня поверхности почвы расположен горизонтальный пласт мощностью а = 0,5 м с меньшими коэффициентом фильтрации и пористостью кв = 0,005; а = 5,25; п = 2; О* = 0,3; в= 0,745. При этом грунт и пласт имеют следующие характеристики по диффузии загрязняющего вещества Бк = 0,1; Бм = 0,1; а =2 и Бк = 0,01; Бм = 0,01; а=2 соответственно. Сорбция вещества грунтом невелика. Концентрация загрязняющего вещества в жидкости равна с=0,1 (г/л). Грунт считается изотропным.

Геометрическая форма области фильтрации, для которой проводится эксперимент, и точки дискретизации области приведены на рис. 6.

В начальный момент времени грунт сухой и имеет давление всасывания, равное \0= — 0,3 м, и концентрацию вещества с0=0 г/л. Граничные условия зададим следующим образом: участки 1-2, 4-5-6-7 - непроницаемы для жидкости. Следовательно, на участке 6-7, 8-1 д\/дх=0, а на участке 1-2, 5-6, 7-8 д\\/д2 = 0, на участке 2-3-4 \ =

На участке 2-3-4 с=с , где концентрация рассчитывается по формуле (2) для каждого шага, а на 4-5-6-71-2 - дс/дп=0.

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

Х

Рис. 6. Геометрическая форма области фильтрации и точки дискретизации в задаче фильтрации раствора жидкости из неровностей грунта

+

Рис. 9. Линии концентрации вещества при времени эксперимента Т=10 сут (1,2-1,2 г/л; 3-0,5 г/л; 4- 0,1 г/л)

Рис. 7. Распространение зоны полного насыщения грунта жидкостью по времени (1- Т= 2 сут; 2- Т= 5 сут; 3- Т= 10 сут; 4- Т= 15 сут; 5- Т= 20 сут)

Рис. 8. Изменение концентрации вещества в водоеме

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

V ¿^- / 1

<4ШШ

Рис. 10. Линии концентрации вещества при времени эксперимента Т=15 сут (1,2-0,7 г/л; 3-0,5 г/л; 4- 0,3 г/л)

Работа выполнена при поддержке Министерства образования, грант Е 02-12.1-175.

Литература

1. Наац И.Э. // Научные школы и научные направления СевКавГТУ. Ставрополь, 2001.

2. Семенчин Е.А., Наац В.И., Наац И.Э. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.

3. Атмосферная турбулентность и моделирование распространения примесей / Под ред. Ф.Т.М. Ньистадта и Х. Ван Допа. Л., 1985.

4. Антонцев С.Н., Епихов Г.Н., Кашеверов А.А. Системное математическое программирование процессов водообмена. Новосибирск, 1986.

5. Бэр Я., Заславски Д., Ирмей С. Физико-математические основы фильтрации воды. М., 1974.

6. Пеньковский В.А., Рыбакова С.Т. // Природные условия Западной Сибири и проблемы переброса рек в Среднюю Азию. Новосибирск, 1975.

Северо-Кавказский государственный технический университет, г. Ставрополь

6 сентября 2004 г.