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

Вестник КРАУНЦ. Физ.-мат. науки. 2011. № 1 (2). C. 46-56

УДК 517.958

ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ ДИФФУЗИИ - АДВЕКЦИИ РАДОНА В МНОГОСЛОЙНЫХ ГЕОЛОГИЧЕСКИХ СРЕДАХ* В.С. Яковлева1, Р.И. Паровик2, 3

1 Томский политехнический университет, 634050, г. Томск, ул. Ленина, 36

2 Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, Камчатский край, c. Паратунка, ул. Мирная, 7

3 Филиал Дальневосточного федерального университета, 683031, г. Петропавловск-Камчатский, ул. Тушканова, 11/1

E-mail: jak@interact.phtd.tpu.ru, romano84@mail.ru

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

Ключевые слова: радон, диффузия, адвекция, многослойная среда, атмосфера

© Яковлева В.С., Паровик Р.И., 2011

MSC 65N80

NUMERICAL SOLUTION OF OF DIFFUSION -ADVECTION EQUATION OF RADON TRANSPORT IN MANY-LAYERED GEOLOGICAL MEDIA V.S. Yakovleva1, R.I. Parovik2, 3

1 Tomsk Polytechnic University, 634050, c. Tomsk, Lenina st., 36, Russia

2 Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch, Russian Academy of Sciences, 684034, Kamchatskiy Kray, Paratunka, Mirnaya st., 7, Russia

3 Far Eastern Branch of Federal University, 683031, Petropavlovsk-Kamchatskiy, Tushkanova st., 11/1, Russia

E-mail: jak@interact.phtd.tpu.ru, romano84@mail.ru

The solutions of stationary and non-stationary diffusion-advection equations of radon transport in many-layered geological media by integro-interpolation method are presented.

Key words: radon, model, diffusion, advection, many-layered media, atmosphere

© Yakovleva V.S., Parovik R.I., 2011

*Работа выполнена при поддержке гранта АВЦП «РНПВШ» № 2.1.1/544.

Введение

Моделирование переноса радона в геологических средах активно используется для решения многих практических задач, а также для выявления закономерностей в поведении радона [1-7]. На основе математических моделей производят оценки параметров переноса и строят прогнозы, находящие свое приложение в различных областях, например в строительстве для введения поправок к существующим нормам и правилам, в геологоразведке при поиске урансодержащих руд [8-10], в радиоэкологии при оценках радоноопасности территорий и зданий [3], в геофизике при изучении литосферно-атмосферных связей [11-15].

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

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

Задача стационарного переноса радона в многослойной геологической среде с разрывными коэффициентами

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

Постановка задачи

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

Задача. Требуется установить распределение объемной активности радона (ОА Ян) по глубине г ^ 0 многослойной геологической среды. Согласно такой постановке уравнение переноса радона в N -слойной среде будет иметь следующий вид:

й

йг

йЛ, (г) йг

йЛ^ - Я (Л,(г) -Л,,„) = 0,

й (°” (г) (г) " А (Л” й "Л».») = 0, (1)

йг (д.+, (г) йЛ”»+г,(г)) + и»+, (г) йЛ”+г,(г) - Я (Л»+, (г) -Л»+,.„) = 0,

I (^ (г) ^ (г) - Я ^(г) - Л^») = 0.

Условия на внешних границах всей многослойной структуры принимаются как:

Л,(0) = 0,11ш ЛN (г) = Л*». (2)

г—^»

Где Л»(г) - активность радона, приходящаяся на единицу объема пористой среды, Бк/м-3, нижний индекс (п = 1,...,N) указывает на номер соответствующего слоя; N - количество слоев; I» - толщина п-го слоя, м; и» - скорость адвективного переноса радона в п-ом слое, м/с; Л» - эффективный (объемный) коэффициент диффузии радона в п-ом слое, м2/с; Я - постоянная распада радона, с-1; Л»,» -объемная активность радона, находящегося в радиоактивном равновесии с 226Яа в »-ом слое, Бк/м-3, которая равна Л»,» = р»,5 (1 - п») , где К»,еш - коэф-

фициент эманирования радона для слоя », отн. ед.; Л»,^а - удельная активность 226Яа для слоя », Бк/кг ; р»,5 - плотность твердых частиц грунта, кг/м-3 для »-го слоя; п» - пористость грунта для »-го слоя.

Методика решения

Рассмотрим расчетную область задачи (1,2). Пусть толщина геологической среды, для которой будем производить расчеты, I является достаточно большой величиной. Тогда ее можно представить в виде суммы соответствующих слоев: I = I, + 12 + ... + 1»-1 + I» + 1»+, + ... + ^ . Разобьем толщину слоя на К, точек, тогда I, = г,йь где /,=0,1,2,... ,К, - 1, й, - шаг сетки для слоя I,. Для слоя 12 аналогично будем иметь 12 = /2й2, где /2 = К2, К2+1,... ,К2 - , Продолжая процесс для »-го слоя, получим: I» = г^й», где г» = К»-,,К»-, + ,...,К» - , Наконец, для последнего слоя будем иметь ^ , где г# = Ку-:,,Ку-, + 1,... ,Ку-1, Ку Для простоты поло-

жим, что й = й, = Й2=... =й»=... =ЙN, т.е. будем рассматривать равномерную сетку по пространственной координате. Тогда для всей толщины грунта будет выполнено соотношение I=1Н, где /=0, 1, 2,... , К, - ,, К,, К,+1,... ,К2 - ,, К2, К2+1,... ,К» - ,,

К», К»+1,... , - ,, К#. Область рассматриваемой задачи приведена на рис.1.

Для решения задачи (1,2) используем интегро-интерполяционный метод. Согласно, этому методу осуществляется переход от системы дифференциальных уравнений (1) к алгебраической системе, которая представлена в разностной форме. Такой переход осуществляется с помощью некоторого интегрального соотношения (уравнения баланса), выражающего закон сохранения для малого объема [16]. Входящие в уравнение баланса интегралы и производные следует заменить приближенными разностными выражениями. Таким образом, для решения задачи переноса радона в слоистых средах используем интегро-интерполяционный метод (метод баланса) построения консервативных разностных схем [17, 18].

Рис. 1. Расчетная область для задачи (1,2)

Запишем для уравнения (1) на »-ом слое уравнение баланса на отрезке г^,^ < г < г^,^ где гг■n-1/2 = й(г» - ,/2), й - шаг разностной схемы, [18]

г,п + 1/2 ^ ^ г,п+1/2

І і (вд"чг>г + /

г;и-1/2 ^—1/2

1А (г) А

ия(г)—— іг—

аг

(3)

-Я

+! /2 гг'и +!/2

Л»(г)йг = -Я J Л»,»(г)йг.

-,/2 ^г'»-,/2

Балансное соотношение (3) отражает закон сохранения для отрезка г^,^ < г < гг»+,/2• Для получения разностного уравнения из балансного соотношения (3) необходимо использовать те или иные восполнения сеточных функций. Функцию решения будем искать в целых узлах (Л(г),г = ), а диффузионный и адвективный

потоки - в полуцелых. Представим первый интеграл разностью диффузионных потоков (д(г) = -Л»(г) йЛ»(г),г = г;»±,/2) в полуцелых узловых точках и запишем их аппроксимацию в соответствии с работой [18]

_ Dn,гn+1/2Лn,г» + 1 - (Dn,гn + ,/2 + Dn,гn-1/2)Лn,гn + Dn,гn-1/2Лn,г»-1 ,

^п,^^ - qn,гn-1/2 = й . (4)

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

гг'и+1/2

ии(г)

~ ’!— / (Ап,ги - Аи,ги-1) +------’!+ / (Ая,ги+1 - Ая,ги) .

и

2

(5)

г,

іп —1/2

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

г,

іп + 1/2

г,

і п + 1/2

Я J А(г)іг ~ ЯМг-И, Я J АИ)^(г)іг ~ ЯМл

(6)

гг'и-1/2

г іп —1/2

п

Выражения (4), (5) и (6) подставим в уравнение (3), что приведет нас к системе алгебраических уравнений следующего вида:

В итоге, решение задачи (1,2) можно записать следующим образом

А1,0 = 0, А1,г- = Р1,г+1А1,г+1 + З1,г+1, А1,К = Р1^ + 1 А1,К + 1 + 31,^ + Ъ А2,г = Р2,г+1А2,г+1 + З2,г+1, Ап,г = Ря,г+1Ая,г+1 + Зп,г+1, Ая,Ки = Ря,Ки + 1Ая,Ки+1 + Зя,Ки+1

АУД = А^, З1 = 0, Р1 = 0,

г = 1,2,..., К - 1, *1, К + 1,..., * - 1, *2, *2 + 1,..., К - 1, К, * + 1,..., Ку - 1.

Задача нестационарного переноса радона в многослойной геологической среде

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

Систему (7) для п-го слоя (1) решаем методом прогонки:

An,in — pn,in+1An,in+1 + qn,in + 1, i#n — Kn—1, Kn—1 + 1, Kn—1 + 2..^ Kn 2

(8)

An+1 ,in+1 — pn+1 ,in+1 + 1 An,in+1 + 1 + qn+1,in+1 + 1 , in+1 — Kn + 2, Kn + З^.^Kn+1 1,Kn+1

pn+1 ,in+l +1

An+1 ,i — Pn+1,i+1An+1,i+1 + 4п+1,І+Ь AN—2,i — PN—2,i+1 AN—2,i+1 + qv^i+b (10)

AN—2,KN—2 — PN—2,Ку—2 + 1 AN—2,Ку—2+1 + qN^Kv^ + b AN—1,i — PN—1,i+1AN —1,i+1 + #^—М+Ъ

Постановка задачи

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

Задача. Требуется установить распределение объемной активности радона по глубине геологической среды и во времени. В нестационарном случае система уравнений (1) для N - слойной среды может быть записана в виде:

дА1(г,г) д ( хдАЦг,г)\ , .дА^г,г) ,*/,/>. Л л А

-Т" ^1 (г,г)—- и (г,г)—+ Я (А1(г,г) - А^) = 0,

дг дг \ ’ д г / ’ дг

дА2<г,г) _ Э ҐВ2 (г,() ^ (г,(} + А _ Аа = 0,

д г дг \ дг у ’ дг

дАи(г,г) д / \ дАп(г,г)\ , дА„(г,г) . , , Л Л , ч

-V ( ПП (г г) ^ ) ип (г г) ^ + Я (Ап(г) Аи,^) — 0, (11)

дг дг V ’ д г ; пк’ ' дг

^„+1^г) д^ ^„+1^ г)^ яч ^„+1^ г) Л_п

_ ~=Г“ ПП+1 г) ----^------ _ и„+1 г) -^-+ Я (А„+1 (г) _ А„+1,^) — 0,

д г д г V дг ) дг

дА*^г) _ д (в* (г,г) ) - и* (г,г) + Я (А*(г) _ А*, „) — 0.

д г д г \ д г у ’ дг

Также на внешних границах грунта задаются краевые условия первого рода

ИшА1 (г, г) = 0, Иш Ау (г, г)= Ау^ (12)

и начальное условие

Ау (г, 0) = Ау,~(г). (13)

Методика решения

Для численного анализа задачи (11,12, 13) применим метод баланса [18]. Рассмотрим расчетную область для этой задачи. Введем равномерные сетки: — г — гй,і — 0,1,2,..,К1,...,,...,К*, й — 1/К* с шагом й на полу-

бесконечном интервале 0 < г < <*>, где К* - достаточно большое число; — {г- — -т,7 — 0,1,...,Т,т — Т/Ки} с шагом т на отрезке 0 < г < Т. Они образуют общую равномерную сетку: юйт — х ют — {(гг-, г-), гг — гй, 0 < г < К*, г- — 7 т, 0 < 7 < Т}, т. е. расчетная область представляет собой прямоугольник. Введем сеточную функцию решения А (гг,?/) — А- є юйт, которую будем рассматривать в целых узлах равномерной сетки, а поток, диффузию и адвекцию - в полуцелых узлах. Согласно уже приведенной методике, после интегрирования по отрезку гв-1/2 < г < гг-и+1/2 для п-го слоя грунта с учетом аппроксимации по времени получаем трехдиагональную систему алгебраических уравнений:

— а А-/'+1 + у А-^1 — в -А 7+1 — р ■

и,гАи,г-1 + /п,гА„,г- Ри,гАи,г+1 — гя,г,

т П

.7+1/2

ап,г —

п-1/2

ти

й2

,/+1/2 п,гп 1/2

"2й

вп,г —

т П

7+1 /2 п,гп + 1/2 й2

ти

+

,/+1/2 п,гп + 1/2

~2й ;

(14)

h2

2h

f 1 + Ят, Fn,i„ — АП,іп + ЯтАп,.І

■n,.i„.

Решение системы (14) можно получить методом прогонки. Применяя эту методику для остальных слоев, получим решение:

г = 1,2,..., К - 1, *1, К + 1,..., * - 1, *2, *2 + 1,..., К - 1, Ки, + 1,..., Ку - 1.

На основе описанной численной модели переноса радона в неоднородных слоистых геологических средах был разработан алгоритм и создана программа «SimRaTran» [19], позволяющая моделировать перенос радона в пористых средах, представленных несколькими (до 20-ти) эманирующими слоями с различными параметрами. Эта программа позволяет рассчитывать распределение объемной активности радона и плотности потока радона по глубине, плотность потока радона с поверхности земли в атмосферу.

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

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

Рассмотрим среду, состоящую из трех слоев мощностью 11 = 4 м, 12 = 2 м, 1з = 8 м, параметры которых (Ои, Ки,еш, Ри,я, Пи) для упрощения анализа, являются одинаковыми. В первом и третьем слое А#а=30 Бк/м3. Во втором слое активность радия, продукта распада урана, была выбрана для расчетов из характерного для урансодержащих пород диапазона А#а= 1000 Бк/м3.

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

На рис. 2 показаны функции зависимости объемной активности радона от глубины при положительной и отрицательной скорости адвекции |и|=10-6 м/с. Там же для сравнения показана функция распределения объемной активности по глубине в

случае однородной среды, т. е. все три слоя представлены суглинком с одинаковым содержанием радия Аяа = 30 Бк/м3, а скорость адвективного потока взята равной 10-6 м/с.

Рис. 2. Кривые распределения объемной активности радона по глубине геосреды с 4-мя слоями. 1 - кривая распределения радона с и = —4 ■ 10-6 м/с; 2 - кривая распределение радона с и = 4 ■ 10—6 м/с; 3 - кривая распределения радона в однородной среде с и = 4 ■ 10—6 м/с ; 4 - кривая распределения радона в однородной среде с и = 0 м/с

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

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

На рис. 3. представлены кривые распределения объемной активности радона в слоистой среде (п = 4) в различные моменты времени ? согласно (14). Мощность слоев составляет: 1\ = 1 м; /2 = 2 м; /3 = 3 м; /4 = 4 м. Для простоты анализа физико-геологические параметры для каждого слоя брали одинаковыми, кроме со-

держаниями радия, которое составляло для 1 слоя - 90 Бк/м3, 2-го - 4 Бк/м3, 3-го

- 30 Бк/м3, 4-го - 1000 Бк/м3.

sfc'K3

» ЇМ ’СО 3(10 -100 *00 SM SM wo

Рис. 3. Распределение объемная активности радона по глубине геосреды с 4-мя слоями в моменты времени г. а) и = 4■ 10—6 м/с при: - (1) г = 107с; - (2), г = 105 с; - (3), г = 2 ■ 105 с; - (4), г = 3 ■ 105 с; в) - и = —4■ 10—6 м/с в те же моменты времени г; б) - и = —4■ 10—6 м/с в момент времени г = 107с; г)- V = 4 ■ 10—6 м/с в момент времени г = 107с

На рис. 3а представлено семейство расчетных кривых распределений объемной активности радона в геологической среде со скоростью адвекции и = 4 ■ 10—6 м/с. На рисунке показано, что на границе 1-го и 2-го слоев объемная активность радона в разные моменты времени составляет А = 25 - 50 кБк/м3, на границе 2-го и 3-го слоев - А = 25 - 100 кБк/м3, на границе 3-го и 4-го слоев - А = 400 - 450 кБк/м3.

На рис. 3в представлено семейство расчетных кривых распределений объемной активности радона в среде с отрицательной скоростью адвекции V = —4■10—6 м/с. На границе 1-го и 2-го слоев объемная активность радона уменьшилась: А = 10 - 25 кБк/м3; на границе 2-го и 3-го слоев А = 25 кБк/м3, на границе 3-го и 4-го слоев А = 200 - 250 кБк/м3. Такой эффект подтверждает тот факт, что при отрицательной скорости адвекции радона суммарный его поток уменьшается.

На (рис. 3,б,г) представлены случаи, когда при достаточно большом моменте времени (t = 107с.) кривые распределений объемной активности радона совпадают со стационарным режимом переноса радона в геологической среде (8). Это подтверждает адекватность построенного численного решения.

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

Заключение

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

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

Литература

1. Паровик Р.И., Ильин И.А., Фирстов П.П. Обобщенная одномерная модель массопереноса радона и его эксхаляция в призеный слой атмосферы // Математическое моделирование. - 2007. -T.19. -№ 11. - С. 43-50.

2. Яковлева В.С. Диффузионно-адвективный перенос радона в многослойных геологических средах // Изв. Томского политех. унив. - 2009. - Т. 315. - № 2. - С. 67-72.

3. Yakovleva V.S. A theoretical method for estimating the characteristics of radon transport in homogeneous soil // Ann. Geophys. - 2005. - № 48(1) - P. 195-198.

4. Van der Spoel W.H., Van der Graaf E.R., De Meijer R.J. Combined diffusive and advective transport of radon in a homogenous column of dry sand // Health Sci. - 1998. - № 74(1). - P. 48-63

5. Etiope G, Martinelli G. Migration of carrier and trace gases in the geosphere: an overview // Phys. Earth Planet In. - 2002. - V. 129. - P. 185-204.

6. Kozak J.A., Reeves H.W., Lewis B.A. Modeling radium and radon transport through soil and vegetation.// J. Contam Hydrol. - 2003. - V.66. - P. 179-200.

7. Nielson K.K., Rogers V.C. Multiphase radon generation and transport in porous materials. // Health Phys. - 1991. - V.60(6) - P. 807-815.

8. Граммаков А.Г., Никонов А.И., Тарфеев Г.П. Радиометрические методы поисков и разведки урановых руд. - М.: Госгеолтехиздат, 1957. - 610 с.