46 Вестник СамГУ — Естественнонаучная серия. 2007. №2(52)
УДК 532.546+628.516
ОДНОМЕРНАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ НЕСМЕШИВАЮЩИХСЯ С ВОДОЙ ОРГАНИЧЕСКИХ ЗАГРЯЗНЕНИЙ В ЗОНЕ АЭРАЦИИ1
© 2007 М.Н. Бардина2
В работе исследуются одномерные задачи о вертикальных потоках. Построено несколько точных решений методами бегущей волны, разделения переменных, подобия и другими. Полученные аналитические решения позволяют делать выводы общего характера, прогнозировать динамику процессов, могут быть использованы в качестве основы для ’’тестовых” задач при проверке корректности и оценке точности численных методов.
Введение
Жидкие неводные загрязнители (NAPL) составляют специфический класс загрязнителей грунтовых вод. Будучи жидкими и образуя отдельную фазу, не смешивающуюся с водой и воздухом, NAPL способны к самостоятельному перемещению в почве, а не только к пассивному переносу потоками грунтовых вод как большинство химических и радиоактивных загрязнителей. Это определяет их специфику и заставляет выделить в особый класс [1-4].
В последние годы проблемы загрязнения грунтовых вод органическими жидкостями привлекают все большее внимание специалистов по подземной гидромеханике и математическому моделированию. Отечественная школа подземной гидравлики занимает ведущее место в мире с тех пор, как начало решения некоторых важных практических задач было положено в трудах академика Лейбензона Л.С. [5]. Однако, круг специалистов в области математического моделирования переноса органических жидкостей грунтовыми водами в Российской Федерации пока недостаточно широк. Из отечественных исследователей активная работа в этом направлении ведется
1 Представлена доктором физико-математических наук профессором
В.И. Астафьевым.
2Бардина Марина Николаевна (appmath@mrsu.ru), кафедра прикладной математики Мордовского государственного университета, 430000, г. Саранск, ул. Большевистская, 68.
Ентовым В.М. [1—4], Мироненко В.А. [6-10], Дерюгиным Ю.Н. [11], Косте-риным А.В. [11].
Сравнительно разработанной является методика предсказания пассивного переноса продуктов растворения первичного загрязнения потоком воды. Иначе дело обстоит с анализом гидродинамической стадии процесса, первичного проникновения органической жидкости в ненасыщенную зону и грунтовые воды и растеканием вдоль водоупора или на поверхности воды [1-4].
Имеется несколько подходов к решению указанной задачи. Универсальный метод приводит к сложной и, как правило, трехмерной задаче, подлежащей численному решению. При этом возникают трудности в связи с определением адекватных и детальных исходных данных. Другой подход основан на выделении ключевых элементов исследуемого класса явлений и расщеплении полной задачи на несколько подзадач. Многие важные выводы удается сделать на основе качественных оценок и простых аналитических решений [4].
Для построения точных решений нелинейных уравнений математической физики разработан ряд методов, основанных на переходе к новым переменным (зависимым и независимым). При этом обычно ставится цель: найти новые переменные, число которых меньше, чем число исходных переменных. Переход к таким переменным приводит к более простым уравнениям. В частности, поиск точных решений с частными производными с двумя независимыми переменными сводится к исследованию обыкновенных дифференциальных уравнений (или систем таких уравнений). Естественно, при указанной редукции решения обыкновенных дифференциальных уравнений дают не все решения исходного уравнения с частными производными, а лишь класс решений, обладающих некоторыми специальными свойствами. Наиболее распространенными классами точных решений, которые описываются обыкновенными дифференциальными уравнениями, являются решения типа бегущей волны и автомодельные решения. Существование этих решений обычно (но не всегда) обусловлено инвариантностью рассматриваемых уравнений относительно преобразований сдвига и растяжения-сжатия
[12].
1. Модельное уравнение
Для описания модели использованы следующие обозначения: m - пористость породы; к - абсолютная проницаемость породы; fw, fo - относительные фазовые проницаемости; иw, Ио - вязкости флюидов; Pcw, Peo - фазовые давления; pw, ро - плотности флюидов; g - ускорение силы тяжести; z - вертикальная координата; Sw и So - насыщенности воды и углеводорода соответственно; qw и qo - фазовые скорости.
Рассмотрена двухфазная фильтрация в зоне аэрации. В основе описания
лежат уравнения сохранения масс и обобщенный закон Дарси для каждой из фаз [1], [2], [11], [13-16].
dmS w dqw kfw (S)
= 0, qw =-----------gradPw, (1.1)
dt dz ^w
d(\-m)S0 . dqo n kf0(S)
-----7-----+ T" = 0. 4o =------------gradPo,
dt dz
Фазовые давления связаны соотношением на капиллярном скачке
PW = PCW + pwgz- (1.2)
Основные движущие силы: сила тяжести, сила трения, выраженная через обобщенный закон Дарси и капиллярные скачки.
Обычно при расчете многофазной фильтрации давление выражают
некоторой аппроксимирующей функцией, в настоящее работе была выбрана функция Ван-Генухтена
= (1.3)
n - 1 n
где В = ----- и у = -----> причем а и п — некоторые постоянные, завися-
n n - 1
щие от параметров породы. Константа n ~ 1 ^ 5 для подавляющего числа образцов грунта, а константа а определяется соотношением Леверетта
1 Гк , .
а----л (1-4)
о \ m
где о — коэффициент поверхностного натяжения.
Относительная фазовая проницаемость имеют сложную зависимость и для ее аппроксимации возьмем функцию, предложенную в работах Барен-блатта-Ентова [1]-[3]:
fw = Sw • (1.5)
С учетом введенных соотношений (1.2)-(1.5) закон изменения водонасы-щенности (1.1) запишется в виде
_ '11 Litf
dt dz \xw
Преобразуя уравнение (1.6), получаем следующую нелинейную зависимость
dmSw d kS}^
dz \ a ^ w
+ Pwg
(1.6)
dmSw _ ^ßy(X + ß-l) x+ß-2 / лї-1 (dS '2
dt a
kß2у (y - 1) c^+2ß-2 /cß Ay~2(dS '2
) o^+2ß-2/eß ,\Y-2 /dSw\
-sw {Sw~l) [-£-)
+l +
Y S,:M <s<j _ yr.
[iivct w V w ’ dz2
, kpwg^ cx-idSw + W Q '
цw w dz
В следующих разделах настоящей статьи приведены точные решения модельного уравнения (1.7) и его численный анализ.
2. Построение решения в общем случае
Для определения точных решений поставленной задачи использовался метод бегущей волны [12], а именно представление водонасыщенности в следующем виде (?), где ^ = г - Для указанных зависимостей
решениями будут являться следующие выражения
1 = Чо +
Ч
к
S^-ldSw
(і - sWv)1 7 («о + аіSw + а2SW)
где для краткости введены обозначения
«о =^£0ІМ (і - хв,,о),-і,
Бт^ а р^а
а\ = —тт;------. а2 = —^—•
£|3у (Зу
Пусть также
I =
к
S}^-ldSw
(і - sW)1 7 («0 + alSw + a2S}w)
(2.8)
Аналитические решения были найдены в нескольких частных случаях. Приведем некоторые из них.
2.1. X = 0, п = 2.
Тогда интеграл решения определяется формулой і
I =---------1п \ао + а2 + +
С1\
2
-_arctg
+
л/а і (а0 + а2)
аі
«о + «2
-sw,
1 л1$№ ~ л / «о + а2
V а\
у-«і («о + «2) л[8цг + ^ 1 а0 + а2 V а\
а\ л/Ецг
«0 + «2 _
если -------------> 0;
«і «о + «2
если
если
«і
ар + а2 «і
< о;
= о.
2.2. I = п = 2.
2’
Из зависимости (2.8) следует, что
1.-1
а\
— а\ - а2 ,
1 \¥ + -------------о-------1П
+2-
-2а0а1 - а1а2 + «2
-X
х
^4а0а1
- а"
аг^
2а1
+ а2
^4а0а1 - а\
если
4а0 а1 - «2 Аа\
> 0;
1 1п 2«1 л/5 цг + С12 — л ^-4а0а1 + а\
4а0 а1 + а2 2а\ + а2+ л ^-4а0а1 + а\
если
если
4а0 а1 - а2 4 а\
4а0 а1 - а2 4 а\
< 0;
= 0.
(2.9)
+
1
2
1
1
2
1
2.3. ао — 0.
Равенство нулю постоянной а0 означает, что заданы нулевые или единичные начальные данные по водонасыщенности. Тогда интеграл (2.8) значительно упрощается
1= Г----------гг2-----------¿Буг
(1_5,1г) Г(а1+а2^^1)
и данный случай допускает еще несколько вариантов. Так, если X — 1, то
При X = 2, п = 2 имеем
I = — л/Бцг------------
а2
аі
8 ЦТ Н-------Г- ІП («1 + «2^ цг) -
а2 а\
а\ . 1а2
—аг^ —Бцг,
аі
аі 1 . І21"
2а і 1
- Л \_a\_
V а2
Л 1 £1 V «2
аі
если — >0;
а2
аі п
если — <0;
а2
аі п
если — = 0.
а2
Положим «2 = 0, X = /|3 + 1, где I > 0 — целое. Если п — рациональное число, справедливо соотношение
' = ~2с!<-1>'т7їММ
в у+ч
ж)
і=0
Также рассмотрен вопрос построения автомодельных решений 5 & (¿, 0 = = (Ю, где = ¿‘^. В общем случае задача сводится к поиску интеграла
■у Х+в~ 1
/= | --^(2.10)
/¡7
не разрешимого в квадратурах. Придавая отдельным постоянным конкретные фиксированные значения, можем найти частные решения исходной задачи. Не смотря на то, что интеграл (2.10) является следствием (2.8) при аі = 0, решения второго не являются следствиями первого, поскольку большинство из них не определено при аі = 0.
2
а
2
3. Аналитические решения частных задач
Исследуем теперь частный случай задачи, когда относительная фазовая проницаемость является постоянной и давление — линейным
относительно насыщенности Рс№ = С& (1 - )• При таких зависимостях
получаем следующее уравнение
дтБ^ _ к/^С№ д25У /пил
—— _ Тт~- (А11,)
дt дz2
Данное выражение является классическим в уравнениях математической физики и других разделах математики. Это уравнение обладает большим запасом точных решений, разобранных многими авторами. Ниже приведены лишь два наиболее важных.
Построение решения задачи (3.11) в виде бегущей волны дает следующую явную зависимость неизвестной функции от переменной 1
тР\1№
к/№С№\ц№0 к/№С№\ц№0 ьгшгш ^ ^
=------^------+ 5>о +-----^------еКМ^№ .
тОЦ& тОЦ&
Автомодельное решение имеет вид
1 т^.№
г- С 1 4К г г2
Бцг = ¿Vо + ут —е^С№ <%' 1 =
VI
Предположим теперь, что относительная фазовая проницаемость является линейной величиной относительно насыщенности и, что
давление — линейное относительно насыщенности Ро№ = С& (1 - 5&)- При таких зависимостях получаем нелинейное уравнение относительно насыщенности воды
дтБ & кС& / дБ & ^2 кС^ д25 & кр&§ дБ &
-<Ь Ш ~ +
------------------------I ---------------- I —
дХ \ дг ) дг2 дг
Метод бегущей волны дает неявную зависимость насыщенности от переменной 1:
? = Чо +---;—— (51 ж ~ Б жо) ~
oWgk + тБццг
кСш о<
+
----------------1п № о + {?wgk + тБ\1№) 5У | +
{рwgk + тБ^)
кС&¥&0 ,1 0 / , гл \о I
-----------------т 1п |\|/^о^ ¡¥0 + \9wgk + тБ\1цг) ^ т \ •
{рwgk + тБ^)
+ь
Автомодельные решения вида (1), 1 = приводят к решениям
либо стационарного характера, либо постоянного. Обобщенные автомодельные решения вида = 1аф (1), 1 = дают следующие результаты
тц№ г2 р^ кР^2
6кС& Х 3С& 6тц&С&
4. Численное решение модельного уравнения и параметрическая зависимость от коэффициентов
Проанализируем влияние тех или иных параметров фильтрационного процесса на протекание указанного процесса. Ниже приведены рисунки, отражающие важные особенности фильтрационного процесса.
В качестве тестовой решалась задача со следующими данными: вязкость воды = 0.001 Н ■ с ■ м-2, плотность воды р& = 1000 кг ■ м-3, постоянная т = 0.3, начальные данные у&0 = 1, 5&0 = 0.1, абсолютная проницаемость
к = 10-12 м2, поверхностное натяжение а^л = 2.5 ■ 10-2 Н/м, скорость волны Б = 3.47 ■ 10-5 м/с. Расчет проводился с шагом 0.01. Сравнение функций, аппроксимирующих относительную фазовую проницаемость, проиллюстрировано на рис. 1. Расчет насыщенности воды для различной глубины зоны аэрации приведен на рис. 2. Поведение профиля волны на глубине Ь = 10 м в зависимости от параметров п, X и пористости породы т показано на рис. 3. И наконец, расчет профиля на последовательные моменты времени — рис. 4.
Рис. 1. Зависимость фазовой проницаемости от насыщенности воды в различных моделях: 1) модель Баренблатта—Ентова; 2) модель Фелми
5. Заключение
В настоящей статье рассмотрена задача движения не смешивающихся с водой органических загрязнителей как задача теории фильтрации. Было построено несколько частных решений для задачи, записанной в общем виде методами бегущей волны и подобия (автомодельные решения), а также решения для задач, записанных в двух конкретных случаях. При этом применялись методы бегущей волны, подобия и другие. Численный анализ модельного уравнения с использованием метода Рунге—Кутта 4 порядка дал результаты, хорошо соответствующие точным решениям.
Рис. 2. Расчет насыщенности воды для различной глубины зоны аэрации для т = 0.3, п = 2, \ = 2, 1) Ь = 5, 2) Ь = 10, 3) Ь = 15, 4) Ь = 20, 5) Ь = 25, 6) Ь = 30
Рис. 3. Расчет насыщенности для Ь = 10 м в зависимости от п, 1) п = 2, 2) п = 3, 3) п = 4, 4) п = 5
Рис. 4. Расчет профиля на последовательные моменты времени 1) / = 1 ч, 2) / =
= 2 ч, 3) / = 3 ч, 4) / = 4 ч, 5) / = 5 ч, 6) / = 6 ч, 7) / = 7 ч, 8) / = 8 ч,
9) / = 9 ч
Литература
[1] Баренблатт, Г.И. Движение жидкостей и газов в природных пластах / Г.И. Баренблатт, В.М. Ентов, В.М. Рыжик. - М.: Недра, 1984. - 211 с.
[2] Баренблатт, Г.И. Теория нестационарной фильтрации жидкости и газа / Г.И. Баренблатт, В.М. Ентов, В.М. Рыжик. - М.: Недра, 1972. -288 с.
[3] Ентов, В.М. К задаче о динамике тонкой оторочки в физико-химической подземной гидромеханике. ПММ
[4] Гольдештейн, Р.В. Качественные методы в механике сплошных сред / Р.В. Гольдештейн, В.М. Ентов. - М.: Наука, 1989. - 224 с.
[5] Маскет, М. Течение однородных жидкостей в пористой среде / М. Мас-кет. - М.-Ижевск: Институт компьютерных исследований, 2004. -640 с.
[6] Мироненко, В.А. Загрязнение подземными углеводородами / В.А. Ми-роненко, Н.С. Петров // Геоэкология. Инженерная геология. Гидрогеология. Геокриология. - 1995. - Вып. 1. - С. 3 - 27.
[7] Мироненко, В.А. Проблемы гидрогеоэкологии. Том 1. Теоретическое изучение и моделирование геомиграционных процессов / В.А. Мироненко, В.Г. Румынин. - М.: Издательство Московского государственного горного университета, 1998. - 611 с.
[8] Мироненко, В.А. Проблемы гидрогеоэкологии. Том 2. Опытно-миграционные исследования / В.А. Мироненко, В.Г. Румынин. - М.: Издательство Московского государственного горного университета, 2002. -394 с.
[9] Мироненко, В.А. Динамика подземных вод / В.А. Мироненко. -М.: Издательство Московского государственного горного университета, 2001. - 519 с.
[10] Мироненко, В.А. Основы гидрогеомеханики / В.А. Мироненко, В.М. Шестаков. - М.: Недра, 1974. - 296 с.
[11] Дерюгин, Ю.Н. Математическое моделирование загрязнения грунтовых и подземных вод / Ю.Н. Дерюгин, А.В. Костерин // Вопросы Атомной Науки и Техники. Сер. Математическое моделирование физических процессов. - 1998. - Вып 4. - C. 26 - 30.
[12] Полянин, А.Д. Методы решения нелинейных уравнений математической физики и механики / А.Д. Полянин, В.Ф. Зайцев, А.И. Журов. -М.: Физматлит, 2005. - 256 с.
[13] Басниев, К.С. Нефтегазовая гидромеханика / К.С. Басниев, Н.М. Дмитриев, Г.Д. Розенберг. - М.-Ижевск: Институт компьютерных исследований, 2005. - 544 с.
[14] Хасанов, М.М. Нелинейные и неравновесные эффекты в реологически сложных средах / М.М. Хасанов, Г.Т. Булгакова. - М.-Ижевск: Институт компьютерных исследований, 2003. - 288 с.
[15] Чарный, И.А. Подземная гидрогазодинамика / И.А. Чарный. - М.: Государственное научно-техническое издательство нефтяной и горно-топливной литературы, 1963. - 396 с.
[16] Чарный, И.А. Подземная гидромеханика / И.А. Чарный. М.;Л.: ОГИЗ - Государственное издательтво технико-теоретической литературы, 1948. - 196 с.
Поступила в редакцию 13/XT7/2006; Paper received 13/X///2006.
в окончательном варианте — 13/X///2006. Paper accepted 13/X///2006.
ONE-DIMENSIONAL MODEL OF ORGANIC POLLUTION DISTRIBUTION NOT MIXING UP WITH WATER IN A ZONE OF AERATION3
© 2007 M.N. Bardina4
In the paper vertical streams one-dimensional problems are studied.
Exact solutions are constructed by methods of a running wave, division of variables, similarity and others. The obtained analytical solutions permit to make conclusions of the common character, to predict dynamics of processes and can be used as a basis for ”test” tasks for checking the correctness and to estimating of accuracy of numerical methods.
3Communicated by Dr. Sci. (Phys. & Math.) Prof. V.I. Astafjev.
4Bardina Marina Nickolaevna (appmath@mrsu.ru), Dept. of Applied Mathematics, Mordovian State University, Saransk, 430000, Russia.