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

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.