CC BY
60
УДК 519.63:517.958
Б01: 10.14529/ттр140203
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ДИФФУЗИИ-АДВЕКЦИИ РАДОНА В КУСОЧНОПОСТОЯННЫХ АНИЗОТРОПНЫХ СЛОИСТЫХ СРЕДАХ С ВКЛЮЧЕНИЯМИ
В.Н. Кризский, А.Р. Нафикова
Актуальность радоновой тематики в различных областях науки и практики до сих пор продолжает расти. В аспекте радиационной безопасности интерес к радону определяется необходимостью защиты человека от патогенного воздействия ионизации, генерируемой этим элементом и дочерними продуктами его распада. Другая сторона радоновой проблемы связана с тем, что радон является одним из индикаторов сейсмогеодинамической активности структур континентальной коры. В этом плане его изучение может внести существенный вклад в понимание закономерностей развития новейшей разломной тектоники и дать значимую информацию для сейсмического прогноза. Также остаются не до конца изученными вопросы, связанные с выявлением и описанием процессов и механизмов переноса радона в различных средах, факторов, обуславливающих временную и пространственную динамику радонового поля, что представляет интерес для определения месторождений углеводородов. Все это в совокупности способствует активному развитию методов математического моделирования процессов переноса радона и его дочерних продуктов распада в различных, в том числе анизотропных средах.
В работе построена математическая модель диффузии-адвекции радона в слоистых анизотропных средах с анизотропными включениями, которая представляет собой краевую задачу математической физики параболического типа. Предложен комбинированный способ решения задачи на основе методов интегральных преобразований, интегральных представлений и граничных интегральных уравнений. Построен алгоритм расчета поля объемной активности радона.
Ключевые слова: диффузия-адвекция радона; анизотропная среда; краевая задача; метод интегральных преобразований и интегральных представлений; преобразование Лапласа.
Введение
Радон, в силу своих специфических особенностей, является индикатором при различных геологических и геотехнических исследованиях. Динамические изменения концентрации радона в приповерхностном слое почвы отражают динамические изменения напряженно-деформированного состояния горного массива, что служит основой для исследования вариаций поля радона как краткосрочного предвестника сейсмических событий [1]. В геологии изотопы радона используются для поиска урановых и ториевых руд, для экологического картирования при выборе площадок под строительство промышленных и жилых сооружений. Повышенная концентрация радона над залежами углеводородов используется для поиска и оконтуривания нефтяных и газовых месторождений.
Изучение процессов распределения радона в грунте и его стока в приземный слой атмосферы связано с решением параболических краевых задач математической физики. Разработка алгоритмов решения подобного типа задач и расчета полей объемной активности радона имеет практическое значение в таких направлениях, как сейсмология, геохимия, разведочная геофизика и т.д.
Постановка задачи и способ решения
Будем рассматривать горизонтально-слоистую модель среды с локальными включениями, отражающую типовую структуру нефтеносного района (см. рисунок).
£2 о.о, Оо.о.Уо.о
Горизонтально-слоистая среда с включениями
Пусть среда разделена гладкими параметрически заданными границами ^г.о = {Yi.o(x, y)ІYi.o ^ г при л/ж2 + у2 ^ то} (г = 0, N) на горизонтальные слои
^0.0, ^1.0, • • •, ^N.0, заполненные веществом, диффузионные свойства которого описываются симметричными тензорами
(я.0 ,п. 0 иг.0\
иXX иху ^хх \
^г.0 ,п$ лг.0 I
аху ауу аух I лг.0 лг.0 лг.0 /
^хх ^ух ^ хх /
и скоростями адвекции ^.0, VI.0, • • • , 0 соответственно.
Каждый слой Пг.0 содержит Мг локальных включений 0,г.^ (] = 1 ,Мг) с границами ^г.з, заполненных веществом, физические свойства которого описываются постоянными симметричными тензорами диффузии
1 (&3 ( хх (і-3 ( ху і-3 (хг
^г-З — (■3 (ху (■3 уу 4І
1 (і'3 \(хг (■3 (уг )
и скоростями адвекции , г = 0^^ = 1, Мг.
Математическая модель переноса радона в области исследования И = У^ и^г.з ^ Е3 может быть представлена начально-краевой задачей вида:
дАгд^) = дЛу(БцУЛг,з (Р,Ь)) + щ.з 9Агд^ - Л(А](Р,Ь) — Л4.'Х>),
Р = Р(х, у, г) € Ог.] ,г = 0, N,] = 0, Мг;
((Ог.оЧ Аг.0(Р,г) ,п)+иг.оАг.о(Р,г))\ъ.0 = ((А+1.оУ Аг+1.0 (Р,Ь),п) + +^г+1.оАг+1.о{р,Ъ))\^ц.о, г = 0^ — 1;
Аг.о(Р, ^\л.0 = Аг+1.о(Р, ^)\л.0, * = 0^ — 1; (1)
((О г.] ЧАг.э (Р,г),п)+щ.з Аг.] (Р,Ь))\Ъ^ = ((А.о V Аг.о(Р,Ь),п) +
+П.оАг.о(Р,ЩЪ^ ,г = 0^^ = 1, Мг;
Аг.з (Р,Щъ,з = Аг.о(Р,Ь)\ъ^ ,г = 0^^ = 1,Мг;
Иш Ам.о(Р,г) = Ам.ж, Нт Ао.о(Р,Ь) = 0;
11т_____ Аг.о(Р,г) = А1 (Р,Ь),г = 0^;
Р€П;.0,\/ж2+У2
Аг.з(Р, 0) = 0,г = 0, N,j = 0, Мг.
Здесь Аг.] (Р, Ь) - объемная активность радона в грунте; Л - постоянная распада радона; Аг.<х - объемная активность радона, находящегося в радиоактивном равновесии с радием (226Ка) в грунте г-го слоя, которая равна Аг.<х = Кг.етАг.варг.3(1 — Пг)), Кг.ет - коэффициент эманирования радона, Аг.яа - удельная активность 226Ка, рг.в - плотность твердых частиц, пг - пористость грунта, А1 (Р, Ь) - нормальное поле радона, описывающее диффузию-адвекцию радона в слоистой среде в предположении отсутствия включений. Переменная Ь > 0 - время.
Если область Оо.о - приземный слой атмосферы, то в задаче (1) следует положить Ао.гс> = 0. При Мо > 0 включения Оол,---, ^о.ы0 могут описывать жилые и производственные сооружения.
Представим искомую функцию объемной активности радона в грунте Аг.] (Р,Ь) в виде суммы двух вспомогательных функций нормального Аг(Р, Ь) и аномального Аг.](Р, Ь) полей, т.е.
Аг.](Р, Ь) = Аг(Р, Ь) + Аг.](Р, Ь),г = 0^^ = 0, Мг, (2)
где нормальное поле радона определяется краевой задачей:
дАг(Р Ь = ёгу(Бг.о^Аг(Р, Ь)) + щ.о дАг(Р,Ь) — Л(Аг(Р, Ь) — Аг.ж),
о . ^ г.о » ^ ■‘■г \ ч ^^ г.о о
дЬ дг
Р € Ог.о,г = 0^; ((Пг^Аг(Р,Ь),п)+иг.оАг(Р,Ь))\Ж0 = ((Ог+1^Аг+1(Р,Ь),п) + щ+1.о Аг+1(Р,Ь))\Ъо ,г = 0^ — 1; Аг(Р, Ь)\л.0 = Аг+1.о(Р, Ь)\л.0,г = 0,N — 1; 11ш AN(Р,Ь) = А^; 11ш Ао(Р,Ь) = 0;
(3)
11ш Аг(Р, Ь) = Аг(г, Ь),г = 0^; Аг(Р, 0) = 0,г = 0, N,
Р&0.1,^/х2+у2
где Аг(г,Ь) - объемная активность радона в кусочно-однородной горизонтально-слоистой среде с плоско-параллельными границами г = гг, г = 0^ — 1 и коэффициентами диффузии (1г = №хо,г = 0^. Способ определения Аг(г,Ь) описан в [2].
С учетом задачи (3) аномальное поле радона удовлетворяет следующей краевой задаче:
^(Р^ = (гу(ОЧШг.] (Р, Ь)) + V] дАгд(Р,Ь) — ЛАг.](Р, Ь),
Р € Ог.], г = 0, N, j = 0, Мг; ((Ог^Аг.о(Р,Ь),п)+иг.оАг.о(Р,Ь))\ъ.0 = ((Бг+^ Аг+1.о(Р,Ь) ,п) + +Уг+1.оАг+1.о(Р,Ь))\^ц.о, г = 0^ — 1;
Аг.о(Р, Ь)\л.0 = Аг+1.о(Р, Ь)\л.0, % = 0,N — 1; (4)
((Бц VAг.j (Р,Ь),п) + иг.] Ац (Р,Ь))\ъ^ = [(Бг^ Аг.о(Р,Ь),п) +
+иг.оАг.о(Р,Ь) + фг.о(Р,Ь)]\ъ^,г = 0^^ = 1, Мг, фг.о(Р,Ь) = ((Бг.о — Бг.] ^Аг(Р,Ь),п) + (иг.о — ^ )Аг(Р,Ь); (*)
Аг.] (Р,Щ~н,з = Аг.о(Р,Ь)\л^ ,г = = 1,Мг;
11ш Аг.о(Р, Ь) = 0,г = 0, N; Аг.](Р, 0) = 0,г = 0, N, j = 0, Мг.
Р
Сделаем в задаче (4) замену вида:
Аг.] (Р,Ь) = е—миг.] (Р',Ь), (5)
где Р' = (х, у, г'),г' = г + Ь. Получим задачу:
(6)
ди ■ ■ (Р' Ь) _ _____ ______
—^ ’ = (гу(Бг.](Р', Ь)), Р € Ог.], г = 0^^ = 0, Мг;
((Бг^щ.о(Р' ,Ь),п)+иг.оЩ.о(Р' ,Ь))\у.о = ((Бг+^щ+1.о(Р' ,Ь),п) +
+Рг+1.оиг+1.о (Р', ЬЩ 0 ,г = 0^ — 1;
((Бг.]Vuг.j(Р',Ь),п) + иг.]иг.](Р',Ь))\у. , = ((Бг^щ.о(Р',Ь),п)+
+Рг.оиг.о(Р',Ь) + /фг.о(Р',Ь))\7>. .,г = 0^^ = 1,Мг; иг.] (РЖт> . 0 = иг.о (Р',Ь))\^ о ,г = 0^^ = 1,Мг;
11ш иг.о(Р',Ь) = 0,г = 0^;
Р '^<О
иг.](Р, 0) = 0,г = 0, N^ = 0, Мг.
Применим к задаче (6) способ решения, описанный в работе [3], используя интегральное преобразование Лапласа
СО
Г{Р',,) = { (7)
о
с формулой обращения
с+го
и(Р',г) = — [ Г(Р,з)е8Ь(з. (8)
2пг ,]
с— гО
Получим следующую краевую задачу:
((гу(Бг.]VFг.j(Р, в)) — вГг.](Р', в) = 0,Р' € Ог.], г = 0^^ = 0, Мг; ((Бг.оЧРг.о(Р' ,в),п) + иг.оГг.о(Р, в))\у. о = ((Бг+1^Рг+1.о(Р', в),п)+ +иг+1.оРг+1.о(Р' , в))\^ 0,г = 0, N — 1;
Рг.о(Р> , в)\ъ.о = рг+1.о(Р , в)\ъ.о,г = 0^ — 1 (9)
((Бг.]Vрг.](Р',в),п) + Гг.](Р',в))\у, . = ((Бг^Рг.о(Р',в),п) +
+иг.о Гг.о(Р',в) + Гф..0 (Р',в))\у.. ,г = 0^^ = 1,Мг.
рФг.о(Р) = ((Бг.о — Бг](Р , в),п) + (иг.о — иг.])Pi(Р, в); Гг.](Р',в)Ь , = Рг.о(Р',в)Ь ., г = 0^^ = 1,Мг;
11ш Гг.](Р',в) = 0,г = 0^,
Р
где функции Рфго(Р') и Рг(Р',в) - есть образы функций 'фг.о(Р',Ь) и Аг(Р',Ь) при преобразовании (7) соответственно.
Для решения задачи (9) рассмотрим вспомогательную задачу для функции Грина О(Р, Я) - функции точечного источника, находящегося в произвольной точке Q(xq, уя, гд) и генерирующего диффузионное поле единичной интенсивности во вмещающем пространстве (в слоистой среде без включений):
(10)
(гу(Бг^Сг.о(Р', Я)) — вСг.о(Р', Я) = —6(Р, Я),Р € Ог.о, г = 0, N;
((Бг.оУСг.о(Р' ,Я),п)+иг.оСг.о(Р' ,ОЩо о = ((Бг+1.оУСг+1.о(Р' ,О),п) + +Щ+1.оСг+1.о(Р> , °))\^' о ,г = 0^ — 1;
Сг.о(Р', Я)\^0 = Сг+1.о(Р', Я)\^0, г = 0^ — 1;
11ш Ог.](Р',Я)=0,г =0,М-
Р
Согласно [3], интегральное представление задачи (9) будет иметь вид:
N Иг ..
Г(Р',в) = ^2^2 Гг.] (Я,в)[(иг.о — иг.] )Ог.о(Р',Я)+ (ц)
г=о]=1 ,
Ъ.з
N Иг ..
+((Бг.о — Бг.] ')VОг.о(р' ,Я),nQ)]dYг,jQ + ЕЕ/ Рфг.о (Я)Ог.о(Р,Я)(^г^ -
г=о ]=1 ,
° у. .
'г.з
Здесь nQ - вектор внешней нормали к границе включения в точке Я, а граничные значения функции Гг.] (Я, в) находятся как решение системы интегральных уравнений Фред-гольма второго рода, формируемых из (9) при Р' € 7г]:
N Иг ..
рг.](Р' , в) — ^2 ^2 / рг.](Я, в)[(иг.о — иг.])Ог.о(Р' , Я)+ (12)
А_П А_1 ^
г=0 ]=1 ] %
з
N Иг
+(Бг.о — Бг] ')VОг.о(р> ,Я),nQ]dYг.jQ Рфг.о (Я)Ог.о(Р',Я)((Ъ^ ,Р' € 'Уц
г=о ]=1 ,
] У.
Таким образом, алгоритм решения исходной задачи (1) имеет вид:
Шаг 1. Определяем нормальное поле радона Ai (z,t) в горизонтально-слоистой кусочнооднородной среде с плоско-параллельными границами z = Zi = const, i = 0,N — 1, коэффициентами диффузии di = dlz0,i = 0,N и скоростями адвекции Vi.o,i = 0,N по алгоритму, описанному в работе [2].
Шаг 2. Если границы слоев z = Ji.o(x,y) = z% = const, то есть среда имеет плоскопараллельные границы, то решение задачи (3) для нормального поля радона найдено: Ai(P, t) = Ai(z, t). Иначе следует решить задачу (2), например, методом интегральных уравнений, формируя их по участкам Ji.o(x,y) = zi.
Шаг 3. Вычисляем функции ^i.o(P', t) на границах включений Yi.j, i = 0,N,j = 1, Mi по формуле (*).
Шаг 4. Для каждого из значений параметра s множества квадратурных узлов численного обращения преобразования Лапласа (в соответствии с алгоритмом в [4]) по формуле (8):
Шаг 4.1. Находим образы F^i0(P1) функций фг.о(Р, t) при преобразованиях (5) и (7).
Шаг 4.2. Находим решение задачи (10) для функции Грина. Оно может быть получено аналитически для случая однородных слоев с плоско-параллельными границами с помощью интегрального преобразования Ханкеля-Вебера.
Шаг 4.3. Формируем систему (12) и находим ее решение - граничные значения функции Fi.j(Q,s).
Шаг 4.4. По формуле (11) определяем решение задачи (9) - функцию Fi.j(Р1, s).
Шаг 4.5. Формируем слагаемое квадратурной формулы для интеграла (8), вычисляя функции u(P',t).
Шаг 5. Находим аномальное поле Ai.j(P,t) по формуле (5).
Шаг 6. Решение исходной задачи (1) - функцию Ai.j(P,t) - получаем по формуле (2).
Заключение
Построена математическая модель диффузии-адвекции радона в слоистых анизотропных средах с анизотропными включениями, которая представляет собой краевую задачу математической физики параболического типа. Предложен комбинированный способ решения задачи на основе методов интегральных преобразований, интегральных представлений и граничных интегральных уравнений. Построен алгоритм расчета поля объемной активности радона.
Литература
1. Уткин, В.И. Газовое дыхание Земли / В.И. Уткин // Соросовский образовательный журнал. - 1997. - № 1. - С. 57-64.
2. Яковлева, В.С. Численное решение уравнения диффузии-адвекции радона в многослойных геологических средах / В.С. Яковлева, Р.И. Паровик // Вестник КРАУНЦ. Физ.-мат. науки. - 2011. - № 1(2). - С. 45-55.
3. Кризский, В.Н. О способе вычисления физических полей в кусочно-анизотропных средах. Часть II. Нестационарные поля / В.Н. Кризский // Вестник Башкирского университета. - 2009. - T. 14, № 4. - C. 1302-1306.
4. Матвеева, Т.А. Некоторые методы обращения преобразования Лапласа и их приложения: дис.... канд. физ.-мат. наук. - СПб., 2003. - 117 с.
Владимир Николаевич Кризский, доктор физико-математических наук, профессор, кафедра «Математическое моделирование:», Стерлитамакский филиал Башкирского государственного университета (г. Стерлитамак, Российская Федерация), [email protected].
Альбина Ринатовна Нафикова, ассистент, кафедра «Математическое моделирование», Стерлитамакский филиал Башкирского государственного университета (г. Стерлитамак, Российская Федерация), [email protected].
Поступила в редакцию 26 декабря 2013 г.
Bulletin of the South Ural State University. Series "Mathematical Modelling, Programming & Computer Software",
2014, vol. 7, no. 2, pp. 38-45.
MSC 65Rxx DOI: 10.1452!) mmpl40203
The Mathematical Modelling of Diffusion and Advection of Radon in Piecewise Anisotropic Layered Media with Inclusions
V.N. Krizsky, Sterlitamak Branch of the Bashkir State University, Sterlitamak, Russian Federation, [email protected],
A.R. Nafikova, Sterlitamak Branch of the Bashkir State University, Sterlitamak, Russian Federation, [email protected]
The use of radon in various areas of science and technology keeps growing. In the radiation safety aspect, the interest to radon stems from the need to protect people from the pathogenic impact of ionization produced by this element and its decay products. The other part of the problem of radon has to do with the fact that radon is an indicator of seismogeodynamic activity in the continental crust. Its study can contribute substantially to the understanding of fault tectonics and yield significant information for seismic forecasts. Some insufficiently studied questions remain related to identifying and describing the processes and mechanics of radon transfer in various media, the factors shaping the temporal and spatial dynamics of the radon field, which is of interest for locating hydrocarbon deposits. All that together promotes the active development of methods for modelling mathematically the transfer of radon and its decay products in various media, including anisotropic media.
In this article we construct a mathematical model of radon diffusion in layered anisotropic media with anisotropic inclusions, which amounts to a parabolic-type boundary value problem of mathematical physics. We propose a combined method for solving the problem based on integral transformations, integral representations, and boundary integral equations.
Keywords: diffusion-advection of radon; anisotropic media; boundary problem; method of integral transformations and integral representations; Laplace transform.
References
1. Utkin V.I. [Gas Breath of EarthJ. Sorosovskiy obrazovatel’nyy zhurnal, 1997, no. 1, pp. 57-64. (in Russian)
2. Yakovleva V.S., Parovik R.I. [The Numerical Solution of the Advection-Diffusion Equation of Radon in Multilayered Geological Media]. Vestnik KRAUNTs. Fiz.-mat. nauki, 2011, no. 1 (2), pp. 45-55. (in Russian)
3. Krizskiy V.N. [About a Way of Calculation of Physical Fields in Piecewise and Anisotropic Media. Part II. Non-stationary Fields]. Vestnik Bashkirskogo universiteta, 2009, vol. 14, no. 4, pp. 1302-1306. (in Russian)
4. Matveeva T.A. Nekotorye metody obrashcheniya preobrazovaniya Laplasa i ikh prilozheniya: dis.... kand. fiz.-mat. nauk [Laplace Transform Some Methods of the Inversion of the Laplace Transform and Their Applications: Dissertation of the Candidate of Physical and Mathematical Sciences]. Sankt-Peterburg, 2003.
Received December 26, 2013
